Numerical simulation of evaporating two-phase flow : application to concentrated solar plants with direct steam generation

Numerical simulations using CFD are conducted on a boiling two-phase flow in order to study the changes in flow patterns during evaporation. A model for heat and mass transfer at the tube inner wall and at the liquid-gas interface is presented. Transport of two custom scalars is solved: one stands for the enthalpy fields in the flow, the other represents a new dispersed vapor phase in the liquid. A correlation is used to model heat and mass transfer at the tube inner wall. The dispersed phase is created at the surface in the liquid and flows up to the liquid-vapor interface. There, it is transformed into actual vapor phase. The multiphase VOF model is validated for the creation of slugs in an horizontal tube for an adiabatic flow. Results are presented for a subcooled boiling flow in a bend.


Introduction
The present study deals with numerical simulation of boiling two-phase flows using Computational Fluid Dynamics (CFD).This complex phenomenon occurs in a large range of applications, including nuclear, petrochemical, chemical, refrigeration and thermal energy industries, and more particularly Concentrated Solar Plants (CSP).The use of solar power to produce energy by means of thermodynamics is less known than by using the photovoltaic effect.The concentrated solar power technologies are based on mirrors concentrating solar radiation on a pipe, heating a working fluid up to high temperature (about 600 K).This fluid is evaporated to produce electricity thanks to a steam turbine coupled to a generator.
The aim of this paper is to present a CFD model, coupling heat and mass transfer at the tube wall and flow pattern simulations.Here, a detailed mechanistic model for heat and mass transfer at the wall will not be used, because focus will be on the liquid-vapor interface and the two-phase flow patterns.Instead, a correlation will be implemented in order to retrieve the quantities of heat and mass transferred at the interface.
The technological and scientific context will be presented, followed by a description of the custom thermal model that is implemented.Then an adiabatic validation of the CFD commercial code will be shown, followed by some preliminary results of our thermal model.

Context
Concentrated Solar Power has been investigated since the 1970s, and is an attractive way to produce electricity, in particular when associated with a heat storage system.The sunlight can be focused on a single point (parabolic and tower systems) or linearly (parabolic-trough and linear Fresnel collectors), which is the situation in this study.Most of the commercial plants use synthetic oil as the working fluid, but the aim of the next generation is to directly boil water as the working fluid in order to reduce operating costs.However two-phase flow occurring in the boiler tube could lead to uncontrolled flow and heat transfer due to its inherent instabilities.
Two-phase flow in tubes can be studied using flow pattern maps, where visually-defined flow patterns are placed on a graph with gas superficial velocity in abscissa and liquid superficial velocity in ordinate.The superficial velocity of the gas phase is defined using the gas mass flow, the total cross-sectional area of the tube and the gas density by the equation ( 1): Most of the existing maps were developed for adiabatic flows like the Taitel & Dukler map [1] for horizontal flows.One of the first non-adiabatic map, based on refrigerants data is that of Kattan et al. [2], which is continuously updated by the addition of new fluids, leading to the Wojtan's map [3].Taitel & Dukler and Wojtan maps [1;3] are plotted in figure 1, along with nominal path of a CSP taken as a continuous evaporation in a single tube subjected to a constant heat flux.Those two maps are not really adapted to high pressure water-steam flow, but it can be seen that boundaries belong to the same orders of magnitude.At low steam quality (low steam superficial velocity) the flow is expected to be intermittent, with large vapor slugs moving along the tube.At a higher quality, the flow is supposed to be annular, with vapor flowing at a high velocity in the center of the tube and liquid flowing in an annulus or being stratified at the bottom of the tube.The aim of the present work is to simulate evaporating two-phase flows in tubes using CFD commercial code Ansys Fluent 14.5 and to compare the results to existing flow pattern maps.The intermittent adiabatic flow regime will be validated, and then a heat and mass transfer model will be developed and implemented.

Problem description
A horizontal gas-liquid flow in a circular tube is modeled.Velocity inlet and pressure outlet boundary conditions are applied, combined to a constant heat flux imposed on the outer wall of the tube, heating the fluid up to saturation, leading to phase change.

Basic transport equations
The objective of the present model is to clearly detect the interface between liquid and vapor, in order to accurately predict the two-phase flow patterns.
In order to track the interface, and to save computational resources, the VOF multiphase model is used.A single set of conservation equations is solved in transient regime, for the mixture (m) of the two phases [4]: In these equations, ρ is the density, v the velocity vector, p the pressure, μ the dynamic viscosity, g the gravity acceleration, E the total energy, k the thermal conductivity and T the temperature.An additional transport equation of the volume fraction is solved for the vapor phase only (α v ): The volume fraction α v represents here the volumetric fraction of vapor in the total cell volume.The liquid phase volume fraction is computed using the following condition: The surface tension at the interface is the result of a force balance between molecules of the two fluids.In the VOF model, the surface tension is modeled by adding a momentum source term to the momentum equation, in the cells that contain the liquid-vapor interface.The interface is captured using the Geo-reconstruct scheme, which represents the interface using a piecewise-linear approach.The turbulence is taken into account using a RNG k-ε model.

Custom model presentation
Along with the basic VOF model, a custom model is implemented using User-Defined Functions (UDF).It consists in different modules playing a role at all the stages of the numerical resolution.

New scalars
In the case of phase change, the enthalpy field is not well solved by the solver, because the latent heat of phase change is not taken into account.To fix this problem, a custom User-Defined Scalar (UDS) is set-up and solved by the solver at the end of the iterative process.It stands for the enthalpy of the mixture h m , and is defined as: The fluid temperature is then calculated from the enthalpy field.It should be noted that for a mixture-like model, the enthalpy of the mixture increases continuously with the addition of energy to the system.
In order not to compute each vapor bubble created at the wall, a virtual third phase is defined, called "dispersed vapor phase" (figure 2).It is represented as a volume fraction of dispersed vapor in each liquid cell, written α d , leaving the wall where nucleation takes place, up to the liquid-vapor interface which is calculated by the VOF model.

Figure 2. Problem definition
This volume fraction is transported by a conservation equation using another UDS, the following way: S d is a volumetric source term including all the creation and suppression terms of the dispersed vapor.The velocity of this phase is the addition of a relative velocity to the mixture velocity computed by the VOF model [5]: The relative velocity is computed at the beginning of each iteration by a simple force equilibrium on a bubble [6].The basic frame of the transport of a scalar φ using UDS in Fluent is the following: In order to adapt equation (8) to equation (10), a new variable is built and the advective flux term of equation ( 10) is customized by a UDF.The dispersed phase transport has been validated by a comparison with the multiphase mixture model with slip velocity already implemented in Ansys Fluent.

Heat transfer at the wall
In order to take into account the creation of vapor due to evaporation at the tube inner wall, a coupled calculation is conducted at the wall: the heat transfer coefficient and source terms of the two UDS are computed simultaneously at the beginning of each iteration.Firstly, the heat transferred by the surface boiling is computed at the beginning using the well-known Rohsenow correlation: where q b is the transferred heat, μ l is the liquid viscosity, C pl is the liquid specific heat capacity, h lv is the latent heat, A S is the surface on which the flux is calculated, Pr l is the liquid Prandtl number, ρ is the density, T w and T sat are respectively the wall and the saturation temperatures, σ is the surface tension and g is the gravity acceleration.This correlation gives the transferred heat as a function of fluid properties and wall superheat.It has been chosen a semi-empirical approach rather than a mechanistic partitioning approach in order to reduce computational costs and to avoid closure models.
From equation (11), a heat transfer coefficient is then computed, and applied at the tube inner wall.The heat flux transmitted to the fluid is introduced as a volumetric source term in the custom enthalpy transport equation (7).A volumetric mass source term is introduced into the "dispersed vapor" volume fraction equation (8).It can be calculated in term of mass flux by dividing the transferred heat by the latent heat of vaporization of the fluid: The created "dispersed vapor" is then transported up to the liquid-vapor interface, and has to be transformed into VOF vapor phase.

Liquid-vapor interface
The mass transfer takes place at the interface, transforming liquid into vapor while deleting "dispersed vapor" phase.Two steps are needed: first, the interface has to be detected, and the source terms are calculated.The interface is defined by two cells: one having a vapor volume fraction below 1% (donor), adjacent to one with volume fraction above 1% (acceptor).In figure 3, some cells on the interface are represented assuming a planar problem.Vapor cells are represented in white and liquid cells in grey.The source terms are represented in the cells, the arrows showing the fluxes trough the faces.A donor (cell number 4 for example) has a sink term, distributed to the acceptors (cells number 1,3 and 7) which have a source term.

Figure 3. Interface description
A first loop is made on every cell of the fluid domain to detect which cells are near the interface and how many neighbors they have.A second loop computes the source terms according to the previous information.All the "dispersed vapor" is taken from the donor cell (sink term in equation ( 8)), and the vapor is distributed into the acceptor cells (source term in equation ( 5)).

Subcooled case
The model presented above works well in the case of an existing VOF interface, which is not the case when subcooled liquid enters the domain.In this case, instead of using a nucleation site on the tube wall to create a VOF vapor bubble, dispersed vapor is converted into actual VOF vapor when it is concentrated in a region.It reflects dispersed bubbles coalescence, and transition from dispersed phase to separated phases.

Source terms balance
Due to the number of interphase terms used in the model, custom mass and energy equilibrium are verified at each time step, including the new scalars equations.Cells and faces loops are run at the end of each time step in order to verify the balance between the different source terms in the system.

Adiabatic validation 4.1 Stability analysis
Although two-phase slug flows were extensively studied experimentally, their complexity and the computational costs made that only few studies were conducted to simulate 3D slug flows in horizontal tubes, and even fewer have compared it to experimental data [7][8][9].It can be shown using an instability analysis based on the Kelvin-Helmholtz theory [1], that waves are formed when the velocity difference between liquid and vapor exceeds a calculated value.
Kelvin-Helmholtz theory is based on a stability analysis of an interface between two fluids (liquid and vapor for example).To simplify the problem, a linear, inviscid and incompressible analysis will be conducted on an infinitesimal perturbation of the interface [10].It consists in solving a force balance (see figure 4) for the wave velocity [1;10;11].

Figure 4. Forces acting on a wave at the gas-liquid interface
The forces acting are gravity, surface tension and inertia.The last tends to increase the wave amplitude due to vapor acceleration lowering the pressure above it, while the two others stabilize it.A perturbation flowing at the surface of the liquid in the z direction can be descripted by a wave equation: where c is the complex velocity: The amplitude of the wave varies in time with the rate , so it grows or lowers depending on the sign of the complex part of the velocity.After some simplifications [11], it can be shown that an instability grows if the following relation between the difference in velocities and the height of the liquid in the tube is satisfied: In figure 5 is depicted with a dark grey line the value of this limit as a function of the liquid fraction in the tube.It has been numerically verified that below a liquid fraction of 50%, the flow is stratified as predicted by Taitel and Dukler [1].This is due to the fact that in this case the vapor acceleration is no longer sufficient to grow a perturbation.Different working points are studied in a 10 meterlong circular horizontal tube of 0.054 meter intern diameter by varying superficial velocities of liquid and vapor while keeping a constant void fraction of 50% at the inlet.The difference in velocities and the liquid fraction are monitored at different distances from the inlet and plotted along with the stability limit (equation ( 15)).As the axial distance increases, the points keep away from the inlet conditions and reach an equilibrium around the calculated stability limit due to a force balance at the interface and mass flow conservation.The majority of them stay below the stability limit, but when a point goes over the limit, a wave is observed as seen in figure 6.  Stability limits are therefore well simulated by the code as verified with a Kelvin-Helmholtz analysis.Waves and then slugs are created and amplified by hydrodynamic instabilities without imposing any oscillating boundary condition at the inlet.

Frequency validation
In order to validate quantitatively the formation of waves growing into slugs, simulations are compared with data from the literature [12].Their experimental installation is a 37 meter-long tube, with an internal diameter of 0.078 meter.A two-phase flow of air and water at atmospheric pressure with different superficial velocities at the inlet results in slugs into the tube.The two inlets are perpendicular to the main axis and separated by a blade.Slug frequency measurements are made at different positions from the inlet.
Numerical simulations were conducted with 3D VOF model under the geometrical and boundary conditions of the literature [12], during a physical time of 60 seconds.In order to save computational resources, half the domain was simulated.The computational mesh is composed of 450 000 cells and 150 for a cross-section area.The results are presented in figure 7. The calculated slug frequencies converge to stable values corresponding to the experiment.The first calculated waves are detected further downstream in the tube, which could be attributed to very stable numerical boundary conditions compared to experimental conditions: hydraulic oscillations due to the pump, flow meters and other technical equipment lead to oscillation at the inlet, contributing to wave formation.The design of the injection, with the flows arriving perpendicularly to the main tube axis also creates large instabilities downstream to the inlet and is not simulated.The effect of the increase of slug frequencies at the beginning of the tube is still present in simulations, but the effect is less pronounced than in the experiment.However the results show that the code is able to reproduce quantitatively the slug formation process.

First thermal model results
The main objective of this study is to simulate evaporating two-phase flows and the evolution of flow patterns during the evaporation process in macrochannels.Very few experimental studies about visualization of such kind of flows currently exist in the literature, due to the difficulties to simultaneously bring heat to the system and observe it through a transparent wall.Yang et al. [13] provides heat by the mean of Joule effect in a transparent conductive oxide coated tube and observe at different flow rates and heat fluxes the flow pattern evolution in a coiled tube.They also compared their experiment with CFD simulation using a temperature imbalance model combined to mass transfer driven by a non-physical coefficient.More recently, experimental work was conducted to describe the evolution of flow patterns for a broader range of operating conditions but without photos of the visualizations [14].
First results of the model descripted above are presented in figures 8-11, which are vertical median planes and front views of the 3D simulated domain.The configuration of Yang [13] is chosen because of its potential validation data.Subcooled liquid enters a 6 mm intern and 8 mm extern diameter tube made of quartz, heated by a constant heat flux.enthalpy in the first cell exceeds saturation enthalpy of the liquid.Dispersed vapor is created at this point as seen in figure 8, and coalesces into vapor calculated by the VOF model (figure 9).At first, small bubbles are created, then coalescing into slugs.The bubbles are stuck at the top of the tube due to gravity.Temperature field is presented in figure 12.The subcooled liquid enters the domain and is heated up to saturation temperature.The temperature at the outlet is close to the saturation temperature.The wall is superheated and is hotter at the top of the tube than at the bottom because of the poor heat exchange with the vapor.

Conclusion
Numerical simulation using CFD was conducted on an evaporating two-phase flow.A new model for heat and mass transfer at the tube inner wall and at the liquid-gas interface was presented.A correlation is used to model heat and mass transfer at the tube inner wall.A dispersed phase is created at the surface in the liquid and flows up to the liquid-vapor interface.There, it is transformed into actual vapor phase.
After having validated the 3D CFD multiphase VOF model for the creation of slugs in an horizontal tube, results including heat transfer are presented for a subcooled boiling flow in a bend.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences,

Figure 5 .
Figure 5. Stability limit and comparison with working points.

Figure 6 .
Figure 6.Numerically obtained interface evolution at a distance of 3 meters from the inlet, as a function of time.

Figure 7 .
Figure 7.Comparison between experimental and numerical slug frequencies.

Figure 9 .
Figure 9. VOF vapor void fraction, side view.Liquid is in white color and vapor in black.
For calculation time reason, only two straight parts and the first bend are simulated yet.The working fluid is refrigerant R141b, boiling at 307.65 K at atmospheric pressure.Subcooled liquid enters the domain at a temperature of 299.15 K.It is heated up to saturated temperature by the heat flux which passes through the wall.Vapor is created when HEAT 2014 03003-p.5

Figure 10 .
Figure 10.Left: volume fraction of VOF vapor.Right: volume fraction of dispersed vapor.

Figure 10
Figure10represents the front view of the system, and figure11shows a cross-section of the tube right in the middle of the bend.Due to the pressure gradient created by the liquid, dispersed and VOF vapor are shifted in the inner part of the bend, which explain that they are not visible in the bend in figure8.

Figure 11 .
Figure 11.Left: Dispersed phase volume fraction.Right: mixture pressure gradient and velocity plane components.