Numerical simulation of the condensation of R-113 in an inclined tube using the VOF method

. This paper presents the results of the numerical simulation of R-113 condensation in an inclined 32 mm inner diameter tube 4 m long, using the VOF method. The dimensions of the tube correspond to the dimensions of the test section in the CJSC NPVP « Turbokon » test facility, where a series of experimental studies are to be carried out. The R-113 condensation regime is examined at a pressure of 0.33 MPa, a total mass flowrate of 180 kg/(m 2 s), and a tube inclination angle of 20°. Data are presented on the distribution of local heat transfer characteristics along the length and over the inner surface of the channel.


Introduction
Designing heat exchange equipment for advanced heat recovery complexes requires reliable data on heat transfer coefficients in condensation of refrigerants in tubes.The Volume of Fluid (VOF) method [1], supplemented with models and algorithms to calculate phase transitions provides detailed information on condensation processes in tubes.Freon condensation in vertical and horizontal tubes was simulated in previous studies [2] and [3], respectively.Comparison with experimental data and predictions by existing methods showed the effectiveness of the VOF method in predicting the distribution of the basic local characteristics of the flow during condensation of refrigerants.This paper presents data on the condensation of refrigerant R 113 in an inclined tube.The dimensions of the tube correspond to the dimensions of the test section of the CJSC NPVP «Turbokon» test facility [4], where a series of experimental studies are to be performed.

Problem Formulation
Figure 1 shows a schematic of the model.Condensation of refrigerant R-113 is studied in an inclined 32 mm inner diameter tube of the test section with a length l0 = 4 m.There was an adiabatic entrance section with a length of lin = 0.1 m upstream of the test section and an adiabatic section with a length of lout = 0.4 m downstream of the test section.The R113 condensation characteristics were investigated at a pressure of 0.33 MPa, a mass flowrate of up to 180 kg/(m 2 s), and an inclination angle of 20°.The flow rate of R-113 and cooling conditions were selected to attain complete condensation of the vapour within the tube.Boundary conditions of the third kind were used on the inner wall *Corresponding author: minkokb@gmail.com of the channel.Therefore, the heat transfer coefficient from the inner wall of the tube to the cooling water and the cooling water temperature were set.The heat transfer coefficient to water was assumed to be 2500 W/(m 2 K), and the water temperature was considered constant along the length of the external cooling channel and equal to 28.5 °C.The thickness of the tube wall required to determine the thermal resistance of the wall was considered equal to 3 mm.According to [3], this approach to setting the boundary conditions adequately reproduce the features of the wall temperature and the heat flux distributions on the inner surface of the tube found in the experimental study [5] into condensation of various refrigerants in horizontal channels.

Mathematical Model
The VOF algorithm implemented in the in-house CFD code ANES was used [6].The mathematical model included the conservation equation for the volume fraction of the liquid phase and mass, momentum, and energy conservation equations: ∂γ/∂τ+∇•(Uγ)=Mgl/ρl ∂(ρcpT)/∂τ+∇•(ρcpUT-λeff∇T)=Qgl (4) Here, γvolume fraction of the liquid phase, τtime, U -velocity vector, Mgl -mass source due to the phase transition, ρdensity, μeff=μ+μteffective dynamic viscosity, μ -molecular dynamic viscosity, μt -turbulent dynamic viscosity, Prgh -static pressure minus the local hydrostatic column pressure, Fb,kprojection of the vector Fb = [-(ρl-ρg)(g•x) + σK]∇γ on the k axis, g -gravity acceleration vector, x -radiusvector of a point in the computational domain, cpspecific heat capacity, Ttemperature, λeff = λ + λteffective thermal conductivity, λ -molecular thermal conductivity, λt -turbulent thermal conductivity, Qglheat sources due to the phase transition.The indices "l" and "g" refer to the liquid and vapor phases, the properties without indices were determined based on the equation: The heat sources Qgl were calculated using the Lee model [7] modified in [8]: Here, h lg is the latent heat of the phase change, ΔTgl is the difference between the saturation temperature Tsat and the temperature of the mixture in a "two-phase" cell, Δxcv is the minimum size along coordinates axes of the current control volume (CV).
CICSAM algorithm [9] was used to solve equation (1).The hydrodynamic equations were solved using the PIMPLE algorithm [6].Unlike the PISO algorithm [10], the PIMPLE algorithm includes iterations at a time step, but their number is much less than in the SIMPLE [11] algorithm since under-relaxation may be abandoned in this case for small time step values.The numerical experiments suggest that the density on the cell faces should be calculated using the values of γ on the faces obtained by solving equation ( 1) by the CICSAM method to satisfy the continuity equation written for mass fluxes.The generally accepted calculation of the value of γ on a face based on the linear interpolation of values in its neighboring cells gives appreciable errors in problems with mass transfer at the interface and highdensity ratio ρl/ρg.
Turbulent viscosity and thermal diffusivity were determined from the Menter k-ω SST models [12] and Prt = 0.9.This model was selected because of the results presented in [2], where the ability of different turbulence models to predict the distribution of local heat transfer along the length of the vertical channel during condensation of R-113 vapour was tested.Comparison with experimental data [13,14] showed that different versions of the k-ω SST models [12] properly predict the change in flow regimes in a film and well reproduce the behaviour of the heat transfer coefficient along the channel recorded in the experiment.
At the inlet, parameters were set, which correspond to a fully developed turbulent flow of saturated vapour.The no-slip wall condition was used.The boundary conditions for the pressure at the outlet of the tube were established taking into account the hydrostatic pressure distribution along the vertical diameter of the tube.If during the simulation a reverse flow appeared at the outlet boundary, it was assumed that the condensate at saturation temperature was ingested into the tube.Preliminary testing has shown that only this approach for determining the outlet boundary conditions can reproduce regimes with complete condensation of the vapour within the tube.These conditions are likely to be the closest to the experimental conditions.If, on the other hand, the ingestion of saturated vapour was simulated at the outlet boundary, then near the outlet region (z > 3.3 m) a stable vapour flow countercurrent with the condensate film flow was formed.In this case, the complete condensation of all incoming vapour occurred within a distance of 3.3 m from the inlet boundary.

Mesh parameters
The computational grid was selected as recommended in [3].Two meshes were used to grid independence test.For the coarse mesh, the number of CVs along the angular coordinate in the computational domain was 16, the characteristic streamwise length of the CV corresponded to 20 capillary constants, and a thin film of condensate on the upper generatrix of the tube near the test section inlet contained about 5 CVs.The number of CVs along each axis has been doubled for the fine mesh.
Figures 2 shows distributions of the vapour void fraction φ obtained on coarse and fine meshes.The distributions obtained on different meshes practically coincide.Also, it can be seen that complete vapour condensation occurs within a length of about 3.3m.The thermal equilibrium vapour quality at the outlet xb = -0.17. Figure 4 shows the dependence of the vapour void fraction on the vapour quality in the steam condensation section.Additionally, Fig. 4 presents the dependence of the vapour void fraction from [12].  Figure 5 shows the distribution of the heat transfer coefficient along the channel length.The heat transfer coefficient in the region of single-phase convection after complete condensation of the vapour was calculated considering the supercooling of the condensate relative to the saturation temperature.Figure 5 also reports the heat transfer coefficient corresponding to single-phase fluid convection.It can be seen that the heat transfer coefficient changes slightly along the length of the channel and remains at a level of 1.5 kW/(m 2 K) within almost the entire length of the test section.The inlet peak is caused by the effect of a short section, where the annular flow pattern occurs.In this case, the film thickness is basically determined by the friction at the interface.Further, a stratified flow pattern arises in the tube, in which the heat transfer coefficient slightly due to an increase in the area occupied by the rivulet in the vicinity of the lower generatrix of the channel.

Figures 3 show distributions
Figure 6 shows the distribution of the local heat flux over the inner surface of the tube.It can be seen that the distribution of the local heat flux along the tube perimeter is quite uniform.Even in the zone occupied by the rivulet, the heat flux is comparable to that on the upper generatrix of the channel.It is interesting to note that the heat flux distribution has a local maximum qw on the lower generatrix of the tube just upstream of the cross-section of complete condensation of the vapour.The reason for the appearance of the maximum is an increase in the flow of condensate in the rivulet.Figures 7 and 8 show the distribution of the velocity Uz normalized to the mass-average velocity at the inlet U0 and the relative turbulent viscosity μt/μ along the vertical diameter of the tube at various distances from the tube inlet.Such high values of turbulent viscosity are due to the fact that the Reynolds number at the inlet is 470•10 3 .The appearance of a complex swirling flow in the channel should also be noted.Figure 9 shows the 2D cross-section streamlines in a section 1.5 m away from the inlet.

Conclusion
The VOF method was used to simulate the condensation of R-113 refrigerant in an inclined tube.The influence of the boundary conditions at the outlet section of the tube in the case of a reverse vapour flow is studied.It is shown that complete condensation of vapour under examined conditions can be achieved only with the ingestion of the condensate, but not saturated vapour, from the outlet boundary.
Under the selected operating conditions, a stratified flow pattern occupies most of the tube length.In this case, the heat flux in the vicinity of the rivulet is quite high, and it should be considered in determining the heat transfer coefficient averaged over the perimeter.

Fig. 2 .
Fig. 2. Distributions of the vapour void fraction φ along the length of the channel during R-113 condensation at 0.33 MPa and a total mass flowrate of 180 kg/(m 2 s).
Figures 3 show distributions of the vapour quality x, and the thermal equilibrium vapour quality xb along the length of the tube.The thermal equilibrium vapour quality at the outlet xb = -0.17.Figure4shows the dependence of the vapour void fraction on the vapour quality in the steam condensation section.Additionally, Fig.4presents the dependence of the vapour void fraction from[12].

Fig. 3 .
Fig. 3. Distributions of the vapour quality (x) and the thermal equilibrium vapour quality (xb) along the length of the channel during R-113 condensation at 0.33 MPa and a total mass flowrate of 180 kg/(m 2 s).

Fig. 4 .
Fig. 4. Dependences of the vapour void fraction (φ) on the vapour quality x for R-113 condensation at 0.33 MPa and total mass flowrate of 180 kg/(m 2 s).

Fig. 5 .
Fig. 5. Distributions of the heat transfer coefficient along the length of the channel during R-113 condensation at 0.33 MPa and a total mass flowrate of 180 kg/(m 2 s).

Fig. 6 .
Fig. 6.Distribution of the heat flux over the inner surface of the tube during R-113 condensation at 0.33 MPa and a mass flowrate of 180 kg/(m 2 s).

Fig. 7 .
Fig. 7. Distributions of the velocity Uz/U0 along vertical tube diameter for R-113 condensation at 0.33 MPa and total mass flowrate of 180 kg/(m 2 s).

Fig. 8 .
Fig. 8. Distributions of the relative turbulent viscosity along the vertical diameter of the tube for R-113 condensation at 0.33 MPa and total mass flowrate of 180 kg/(m 2 s).

Fig. 9 .
Fig. 9. 2D cross-section streamlines in a section 1.5 m away from the inlet.The interface is shown as a black line.R-113 condensation at 0.33 MPa and total mass flowrate of 180 kg/(m 2 s).