MULTIBODY SIMULATION FOR GRAVITY-DRIVEN FALLING FILM USING COMPUTATIONAL FLUID DYNAMICS SIMULAÇÃO MULTICORPO DE FILME LÍQUIDO GOVERNADO POR

The annular flow two-phase regime appears in nuclear power plants cooling by water. On accidents with lost of cooling this regime determine the heat extraction conditions from the core. The experimental facility Gerador de Película Ondulatoria (GEPELON) was constructed at Universitat Politècnica de València to study this phenomenon. In this work, we use computational fluid dynamics techniques to simulate an air-water regime, specifically a gravity-driven falling liquid film in the GEPELON experimental facility. We consider a multibody geometry where the porous sintered of the input device is simulated in a realistic way. The Volume of Fluid (VOF) model was used in the resolution of the constitutive equations. We propose three geometric approximations to estimate film thickness, and simulation results were compared with experimental values. The simulation results show a good agreement with experimental measurements. Besides, no significant differences were observed between geometrical approximations for film thickness calculations.


INTRODUCTION
Two-phase regimes appear in several engineering applications like heat exchangers, oil distillation towers, and chemical and nuclear reactors. Among these gas-liquid two-phase regimes, the annular flow is particularly interesting because of its complexity. In annular flow, we have a continuous gaseous phase at the channel center with a liquid annular film between gas and channel wall. In addition, we have waves in the gas-liquid interface and disperse drops in the gas phase (FANCHI, 2018). This regime appears during normal operation of boiling water nuclear reactors and also, in all water-cooled reactors in loss-of-coolant situations. The hydrodynamic stability and heat transfers studies in the annular regime are relevant for the energetic nuclear reactors operational safety (TODREAS; KAZIMI, 2011).
The difference between the gas core and liquid film velocities generates interfacial waves and in some cases the droplet entrainment. These phenomena have a strong influence on the hydrodynamic behavior and intensify the heat transfer, increasing the problem complexity. Several researchers have studied these regimes using different approaches such as: The GEPELON (Generador de Película Ondulatória) experimental facility at the Universitat Politècnica de València was built for the analysis of this phenomenon. The main goal of this facility is the study of annular flow regimes (air-water and steam-water) on vertical pipes, reproducing the operational conditions and the accident scenarios in nuclear reactors. The GEPELON facility measures the water film thickness using conductivity sensors allowing the characterization of the interface between the phases (MUÑOZ-COBO et al., 2017).
In this work, we use a CFD approach to simulate an air-water annular regime, associated with a gravity-driven falling liquid film in a vertical pipe. A previous work was developed in (IGLESIAS et al., 2018) where the authors simulated this regime considering a simplified geometry of the GEPELON facility. This simplification includes a one-body two-dimensional Cartesian geometry, and the annular flow entrance approximated by a fake crown. The previous results suggested an improvement in the geometry and also in the simulations models. As Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 a continuity, our research introduces a multi-body simulation, where the GEPELON facility is modeled with an axis-symmetric cylindrical geometry, and a realistic approach is used for the sintered porous entrance. In addition, were considered three interface reconstruction geometrical models selecting the best match with experimental results.
The next section offers the simulation details like the GEPELON experimental facility parameters, simulation conditions, and the techniques used for the liquid film thickness calculations. The simulation results are presented and commented in section 3. Finally, section 4 brings the research conclusions and future works.

Experimental Facility
The Generador de Película Ondulatória (GEPELON) facility was built in the Thermalhydraulic research laboratory of the Energetic Engineering Institute at the Universitat Politècnica de València with several research project resources. The facility is flexible and can be configured for various air-water and steam-water regimes, also it has dimensions and operational parameters alike to those found it in industrial plants. In this work we consider a water falling film regime, as is shown in Fig. 1.
In the experimental facility is a water container at atmospheric pressure, coupled to a set of flow controllers (valves, pumps and instrumentation). The flow controllers send the water to a smaller pressured tank located at the test pipe upper extremity. In this tank the water transversely passes through a porous steel sintered filter, entering in the test pipe with an annular distribution. The test pipe central region is filled with stagnant air at atmospheric pressure.
The test section is a methacrylate pipe with 42E-03 m radius and 3.80 m longitude. Along the test pipe are four conductivity sensors, located at 0.25, 1.10, 1.58 and 3.50 meters from the inlet. The conductivity sensors gather the water film thickness over time, allowing the average film thickness determination. In this work was considered water flows between 20 L/m (0.3319 kg/s) and 2 L/m (0.0332 kg/s) with a water temperature of 25 • C and air temperature of 23 • C.

Simulation Conditions
To simulate the gravity driven water falling film regime in the GEPELON facility, we considered an axisymmetric cylindrical geometry that represents the first 0.70 m of the test section as it's shown in Fig. 2(a). The geometric model includes the porous sintered steel filter and the test pipe initial section. We highlight that the film annular inlet is modeled in a realistic Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 way because the porous media and its parameters, was considered in the simulation. The built model includes three coupled bodies, there are: the porous media, and the pipe divided into two bodies, one of the first 0.2 m, and the other for the remaining 0.5 m to the model outlet.
The model used as boundary conditions: inlet mass flow at the filter outer surface, constant pressure of the mixture in the pipe outlet, symmetric conditions in the pipe center line, open border to the atmospheric pressure in the pipe upper extremity and non-slip boundary in the pipe wall. These conditions are consistent with the physics of the modelled problem.
Next, a regular structured mesh was built in the domain, considering the element connections between the bodies' interfaces assuring the model quantities continuity. In water falling film regime the liquid phase slides down in the pipe wall (domain right boundary) and the pipe central region is occupied by the air phase. As our main goal is to estimate the water film thickness, when the mesh is created, we have to pay attention to the region near the pipe wall where the phases interaction occurs and strong gradients appears. For that reason, we used a mesh inflation with finer elements in this region, allowing that the simulation results better describe the behavior of the problem. The mesh inflation has 12 elements, being the closest to the wall, of 1E-04 m height. From the second element the height increases with a Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 Source: Authors.
1.15 ratio. A mesh representation and the inflation details are offered in Fig. 2(b) and (c). The constructed mesh has 22000 elements, 22968 nodes and a minimum orthogonal quality of 0.99. This configuration is appropriate for the addressed problem.
In this work, the Ansys Fluent version 118.1 (ANSYS, 2017a) was used. To incorporate the pressure drop inside the porous media, the momentum conservation equations include a momentum drain. This adjustment permits to reproduce the phenomenological behavior. For the sintered filter, specifically, we have pores only in radial direction, then the flow movement happens only in this direction. The pressure drops are calculated as where ∆p is the pressure drop, µ is the fluid dynamic viscosity, ρ its density, η is the porous Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 media width, and v is the fluid velocity. In the left side of eq. (1), the first term is associated to the viscous resistance and the second term to the inertial resistance. For that reason, the α coefficient is the media permeability, and C 2 inertial resistance coefficient. These values are model input data, and were calculated using an adjustment of the experimental data provided by the filter manufacturer (AMESPORE, 2016). In our simulation the permeability was 2,1712E-13 m 2 and inertial resistance coefficient was 7,5353E+08 m. In a high-density porous media with small pores and slow fluid velocity, like our case, the main contribution in the pressure drop is associated with the viscosity term, the inertial losses have a weaker influence.
To simulate the two-phases regime we used the Volume of Fluid model (VOF) and a compressive interface scheme to estimate the volume fraction for each phase. In addition, we describe the interaction between phases using the Continuum Surface Force model (CSF) and use the Wall Adhesion model for a detailed wall treatment. We do not consider the turbulence phenomenon because our simulation domain is the beginning of the test pipe, Steady-state simulations were performed for several inlet water mass flows. The convergence criteria to stop the iterative process was that the residuals were less than 1E-03 for all the unknowns, or stable residues with inlet-outlet mass imbalance lower than 1E-05. Due to the studied phenomenon instability, we used Fluent's pseudo-transient simulation mechanism with low relaxation coefficients. We highlight that to reach the convergence criteria for the lower mass flows rate, we had to use very low relaxation coefficient and also an elevated number of iterations.

Film Thickness Approximations
The CFD simulation output variables in each element are the pressure drop, the mixture's velocity, and mass flow, and also the volume fraction of the phases. However, the film thickness (δ ), our main simulation target variable, is not an output variable, then, it has to be estimated from the known quantities. To calculate the film thickness, the volumetric fraction and geometrical criteria were used. In this work were considered three geometrical approaches for the liquid profile inside the element: 1. Constant liquid profile, in this case, the film height is constant and proportional to the Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 volumetric fraction. The main advantage of this approximation is the simplicity, though, there are small discontinuities between the elements.
2. Sawtooth liquid profile, in this case, the film has tooth shape, being the film area proportional to the volumetric fraction. The film continuity is required in the element's interfaces. The main disadvantage of this approach is the non-differentiable points that appear in the center and edges of the element.
3. Spline form liquid profile, here the film profile is a bicubic spline (HOFFMAN; FRANKEL, 2001), the film area is proportional to the volumetric fraction, and smooth and continuous functions appear in all the elements.
A representation of these approximations is offered in Fig. 3. These three profile reconstructions were implemented and compared with the experimental results. Also, the reconstructions methods' computational performances were analyzed.

RESULTS AND DISCUSSION
In this work were performed steady-state CFD simulations for the conditions described in section 2.2. The inlet mass flow varies from 20 L/m to 2 L/m, with a step of 2 L/m. For inlet mass flows over 10 L/m, adjusting the Fluent pseudo-transient relaxation coefficients, the convergence is reached. However, for lower inlet mass flows the simulate became unstable and the convergence was not reached. The modeled phenomenon is dynamic, with continuous changes in the air-water interface, this behavior is significant for thinner films. Thus, we consider the instabilities in steady-state simulations consistent with the problem physics. In these cases, non-stationary simulations must be used for the film thickness characterization. We also consider that a mesh improvement, using finer elements in the wall can contribute to the simulation stability.
The computational experiments were developed in a personal computer, with the following configuration: Intel Xeon E5-1660 V3 processor with 8 physical cores and a RAM of 16 GB. The mean time consuming for each simulation was 40 minutes.
For the simulations with water flows higher than 10 L/m was calculated the liquid film thickness profile using the approximations described in section 2.3. The results at sensor 1 position (z = 0.45 m) and the corresponding experimental values are shown in Table 1 and   Table 1 and Fig. 4 we can observe a liquid film thickness decrease when the inlet mass flow is reduced. This behavior is consistent with the problem phenomenology. Comparing the numerical and experimental results we note that the relative deviations are lower than 8% in all simulated cases showing a good agreement among experimental and simulation values. The spline reconstruction method offers the minimum relative deviation and the sawtooth the Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 maximum. Between the spline and constant approximations, the results are similar, with a mean relative deviation of approximately 3%. On the other hand, if we consider computational performance, the spline approximation is very expensive. For that reason, even with the discontinuities that appears in the element interfaces, the constant profile reconstruction is the best alternative. Source: Authors.
In Fig.4 we remark that the simulated values reproduce the experimental behavior, with small differences when the profile reconstruction approximation changes. Is important to point out that these results were generated using a fine mesh, is possible that when coarse meshes are used more significant differences between reconstruction approximations, appear. Besides, a mesh inflation improvement, reducing the size of the elements near the wall can bring more stability to the simulation.

CONCLUSIONS
In this work we presented CFD simulations results for a gravity-driven falling liquid film in a vertical pipe, considering a multi-body simulation of the GEPELON experimental facility. The main contributions involve the realistic representation of the sintered porous filter located Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 285-01, 285-12, 2020 DOI: 10.21575/25254782rmetg2020vol5n61335 at the test section entrance and the study of different geometrical approximations to calculate the film thickness.
For water mass flow between 20 and 10 L/m the simulation results are consistent and show a good agreement with the experimental data. For lower water mass flows the simulation became unstable and the convergence was not reached.
All the geometrical approximations to estimate the water film thickness generate accurate results. Considering computational performance, the constant approach is more appropriated. The future works involve dynamic simulations to study the evolution of the waves inside the liquid film over time. Also, the proposed simulation models and the used mesh must be improved to reach convergence for lower mass flows.