Consistent scalar transport with front capturing methods: application to two-phase heat transfer

Accurate numerical simulations of heat transfers in 3D liquid-gas flows are of first importance in multiple industrial applications such as droplet evaporation in combustion chambers, spray cooling or propellant behavior in cryogenic tanks. Liquid-gas flows are characterized by the discontinuity of properties (e.g. viscosity, thermal diffusivities, ...) and flow variables (e.g. pressure, energy, chemical composition) across the interface. Robust and accurate algorithms are then necessary to transport the flow variables consistently with the interface. This work presents an algebraic interface capturing method which enables a consistent transport of scalars with the interface. This two-fluid approach ensures conservation of the transported scalars while controlling accurately the flux at the interface. The interface represented by a hyperbolic tangent profile is transported and reinitialized without geometric reconstruction. The scalars are transported separately in each phase and a reinitialization step is performed in order to impose the correct flux at the interface. The method is implemented in the YALES2 low-Mach number flow solver and takes advantage of adaptive unstructured grids to handle complex geometries. It has been assessed on various test cases and compared to analytical solutions.


Introduction
Liquid-gas flows are ubiquitous in nature and industrial systems. Numerous numerical methods have been designed to track and capture liquid/gas interfaces. Modeling scalar transport and diffusion in such flows such as temperature or species mass fractions is a challenging task because the numerical methods used for the interface capturing and scalar transport are different. However, accurate multi-physics liquid/gas flow simulations require a consistent transport of the interface and of the scalars while ensuring the proper jump conditions on the scalars and their gradients at the interface. The latter condition is particularly difficult to satisfy, especially when the gradients at the interface are sharp. In one-fluid methods, two main families of methods have been designed to deal with the jump conditions: smooth and sharp methods. In smoothed interface methods, the jump conditions are smoothed out on several control volumes in order to use standard numerical schemes. In sharp methods, such as the Ghost-Fluid Method (GFM) [1], the interface is considered to have zero thickness and the jump conditions are directly built into modified numerical schemes at the interface. The GFM enables for instance to impose a jump condition on the scalars at the interface or a jump on the scalar gradients or both. While this approach is suited for certain configurations [2,3], this method is not conservative and it is not satisfactory when the interface flux must be controlled accurately. One way to improve this scalar transport is to rely on a two-fluid method in which two different scalars for each phase are transported and coupled together at the interface. In this type of method, minimizing numerical errors and imposing the correct flux at the interface are challenging. Recently, promising results have been obtained by Jain et al. [4] following the two-fluid approach and using the phase-field method to represent the interface. However, their method imposes a no-flux boundary condition at the interface. The original method presented in this paper allows to alleviate these issues. The interface is modelled by a hyperbolic tangent profile, which represents the liquid volume fraction and thus avoiding geometrical reconstructions. The scalars are weighted by their respective volume fraction and transported over the whole domain with a source term which represents the flux between the two phases. The transport of the scalars is supplemented by a reinitialization equation which imposes the correct scalar profile along the interface normal. The paper is organized as follows: first, the framework, the scalar transport and reinitialization equations are presented. Second, the interface flux model is exposed. Then, the various validation cases are detailed and discussed. Finally, the method is applied to a 3D head-on collision case with dynamic mesh adaptation.

ACLS framework
The proposed framework builds upon the ACLS [5] in which a conservative level set method is used to capture the interface transport. The liquid-gas interface Γ is represented by the liquid volume fraction ψ(x, t) and assumed to have a hyperbolic tangent profile along the normal: where 2ε is the thickness of the profile and φ(x, t) = ± min is the signed distance function to the interface. In the present framework, which builds upon the ACLS method, the interface position is not reconstructed from the liquid volume fraction ψ(x, t) as in volumeof-fluid methods with geometric interface reconstruction. Instead, the interface corresponds to the iso-surface: Γ(t) = {x ∈ R 3 | ψ(x, t) = 0.5}. Making the assumption that the flow is incompressible, the interface is advected and reinitialized [5] in conservative form: In this work, the alternative form of the reinitialization from Chiodi et al. [9] is preferred: where n = ∇φ/ ∇φ 2 is the interface normal and φ map = ε ln (ψ/(1 − ψ)) is a mapped signeddistance function for ψ ∈]0; 1[. The signed-distance function φ is reconstructed at nodes in a narrow-band around the interface based on the Geometric Projection Marker Method [6].

Scalar transport
A general framework is now derived for scalar transport and diffusion in liquid-gas flows. The starting point is the heat equation with homogeneous properties which is valid in each phase: where the subscript l and g denote the liquid and gas phases respectively. T is the temperature, u is the velocity of the flow and ρ, c p and λ are the density, the isobaric specific heat per unit mass and the thermal conductivity, respectively. Assuming no phase change, the heat flux and the temperature are continuous across the interface: where the subscript Γ indicates that the terms are evaluated at the interface. To ensure proper conservation of the liquid and gas temperatures, Eq. 4 is multiplied by the volume fraction to obtain a two-fluid model: where ψ l = ψ and ψ g = 1 − ψ are the volume fractions of the liquid and gas phases and β · ∇ψ i = λ i ∇T i · ∇ψ i is related to the heat flux at the interface. Eq. 6 is a transport equation for the scalar ψ i T i that must be solved in each phase. Away from the interface, Eq. 6 reduces to Eq. 4. Therefore, the proposed method has no impact away from the interface. Transporting ψ i T i instead of T i within a finite-volume framework guarantees to have the heat conservation in the domain without external heat fluxes. Moreover, heat flux at the interface (Eq. 5) which is a boundary condition for Eq. 4 appears now as a source term through β · ∇ψ i . A model that ensures that the correct heat transfer from one phase to the other will be derived hereafter.

Scalar reinitialization
The numerical solving of Eq. 6 generates errors in the transported scalars, which may be inconsistent with the errors from the interface transport and reinitialization of ψ. The proposed solution is based on the reinitialization from Chiodi et al. [9]: a reinitialization step is performed in order to reshape the transported scalars and to ensure that the correct Neumann condition is imposed in a narrow band around the interface. The scalar reinitialization is written: where τ is a pseudo-time and α = 4ψ(1 − ψ). The solving of Eq. 7 imposes the same interface flux as the one in the source term of Eq. 6 for consistency. In the divergence operator, the temperature gradient and the term with β are those that would be found in a classical Hamilton-Jacobi equation, and εn and ψ i are present for consistency with the interface reinitialization. The function α = 4ψ(1 − ψ) is a mask function equal to one at the interface position, and which controls where the reinitialization is applied. It is important to note that Eq. 7 is conservative and the proper Neumann condition is retrieved when converged. Imposing a Dirichlet boundary condition on the interface temperature has to be performed indirectly in the interface flux model.

Interface flux model
In order to solve Eq. 6, it is necessary to model β = λ i ∇T i or more precisely, its mean value over a time step j. This model has to take into account that the two temperatures on each side of the interface are not necessarily the same. A first model can be obtained from the analytical solution of the 1D two-phase heat diffusion equation on semi-infinite domains with two initially constant temperatures [11]: i=l,g, where ∆t is the time step. T Γ,l and T Γ,g are the temperatures of the liquid and gas phases in the vicinity of the interface. e l and e g are the thermal effusivities of the liquid and gas phases, given by e i = λ i ρ i c p i . This expression is injected in the source term of the transport (Eq. 6) and reinitialization (Eq. 7) equations. The model indicates that the heat flux between the fluids is proportional to the temperature difference and that the fluid with the lowest effusivity drives the temperature gradient at the interface. The model is based on two semi-infinite domains, which is not the case in the finite-volume framework where the flux is imposed locally. In order to ensure the boundedness of the liquid and gaseous temperatures, this flux model (Eq. 8) has to be limited. A physically based criterion can be derived using the interface temperature [11]: where T n Γ represents the interface temperature and the exponent n corresponds to the n-th time step. Eq. 9 indicates that the temperature at the interface is the barycenter of the liquid and gaseous temperatures in the vicinity of the interface weighted by their mutual effusivities. Assuming that the gas temperature is higher than the liquid temperature, the interface flux j · n (Eq. 8) is negative and the following relations must always be satisfied: In the finite-volume framework, the previous conditions are written: where η∆x is a characteristic length of the control volume. Replacing T n+1 Γ by its expression (Eq. 9) leads to: where (j · n) n can be either positive or negative. The previous equation imposes a limit on the value of the interface flux. If the modeled flux satisfies the boundedness criterion (Eq. 12), the interface flux is given by Eq. 8. Else, the interface flux is limited to satisfy Eq. 12. The control volume characteristic size η∆x is set to half of the mean length of the edges connected to a node. β has to be defined in the narrow-band of the level-set to solve Eq. 7. The Fast-Marching Method is used to extend β, it ensures the propagation of the interface flux β in the whole narrow-band.

Verification & Validation
The method presented in this paper has been assessed on different test cases. 1D results are first presented in which the numerical solutions are compared with analytical solutions. The method is then applied to a 2D static and moving droplet.

1D cases with imposed flux
Consider a 1D domain Ω = [−0.5; 0.5] of length L = 1. Initially water at 300 K fills the subdomain x ≥ 0 and air at 350 K fills the sub-domain x ≤ 0. Both fluids are static and a heat flux β ·n = 12 W/m 2 is imposed at the interface. Dirichlet boundary conditions given by T (−0.5, t) = 350 K and T (0.5, t) = 300 K are imposed. Fig. 1 shows the analytical solution on semi-infinite domains [11] and the numerical solution given by T (x, t) = ψ(x)T l (x, t)+(1−ψ(x))T g (x, t) after t=400 s. The agreement between the analytical solution and the numerical one is very good.
To assess the accuracy of the method, the errors based on the L 2 and L ∞ norms are computed as follows: where L is the length of the domain, N is the number of nodes and L i is the length of the i-th control volume. T ex,i and T i are the analytical [11] and numerical solutions respectively. Figure 2 shows the errors based on the L 2 and L ∞ norms for different values of ε. L 2 and L ∞ norms converge and saturate a high resolution. The finite domain size is suspected here.

1D cases with modeled flux
The method is now assessed using the interface flux model derived above. The configuration is similar to the previous section. Fig. 3 compares the analytical solution on semi-infinite domains [11] to the numerical solution after t = 400 s. Our method correctly reproduces the behaviour The mesh resolution is L/∆x = 901. of the analytical solution and the gradient at the interface is well predicted. The errors based on the L 2 and L ∞ norms depicted in Fig. 4 show a convergence order between 1 and 2 for the different interface thicknesses investigated. The method can also be applied to a water-oil case. This case is particularly challenging since the interface flux is very high. Fig. 5 compares the analytical solution to the numerical one after t = 10000s. Our method accurately predicts the interface gradients in each phase.  The mesh resolution is L/∆x = 901.

2D static and convected droplet
This section presents the numerical results obtained for static and moving water droplet at 300 K in a hot air environment at 350 K. Simulations are performed on a 2D triangle-based unstructured grid. Figure 6, depicting the temperature field for the static droplet case, shows that the liquid temperature in the droplet is almost not affected by heat transfers, while air temperature is strongly modified in the droplet vicinity and diffuses towards the far field. As mentioned previously, this observation may be accounted for by the high effusivity ratio between liquid and gas leading to temperature gradient driven by the gas phase.     3D head-on droplet collision using adaptative mesh refinement A challenging test case for the validation of two-phase modeling approaches is the interface merging scenario [7]. The case investigated here corresponds to the experimental test of [12] where two water droplets of equal size collide with Weber number We=23 and Ohnesorge number Oh=0.0047. Numerically, the water droplets are initialized with two different temperatures: T=300 K (upper droplet) and T=325 K (lower droplet) in a hot air environment at temperature T=350 K. The dynamic and parallel mesh adaptation strategy introduced by [14] and adapted to two-phase flows in [15] is used. The targeted metric in the interface vicinity is ∆x = 6 µm. The proposed approach is robust and accurate except at the droplet merging, a configuration which is not accounted for in the flux model. To prevent to have unrealistic fluxes, a flux limit at 1 × 10 5 W/m 2 is imposed and activates only at merging. Instantaneous mesh, interface droplet and temperature field of the colliding droplets are shown in Fig. 9.

Conclusions
An original two-fluid approach to handle transport of scalars in two-phase flows is presented here and applied to heat transfer problems. The volume fraction weighted phase scalars are transported consistently with the thickened interface represented by an hyperbolic tangent profile. Contrarily to classic one-fluid model, the interfacial flux is imposed using a reinitialization