Modelling of turbulence-interface interactions in stratified two-phase flows

In this work we describe interactions between turbulence and water-air surface in the ensemble-averaged picture where, instead of a sharp interface between the phases we deal with a surface layer where the probability of the surface position is nonzero. Changes of the turbulent kinetic energy and the characteristic size of the eddy influences the width of the surface layer. We present a numerical solution of the equation for the intermittency function which describes the probability of finding the water phase at the given point and time.


Introduction
In the present work we investigate stratified turbulent air-water flows where the phases are separated by a deformable, but non-broken interface. We consider the statistical approach where the ensemble averaging of physical quantities, including these connected with the fluctuating surface, is performed. After the averaging we do not deal with a sharp boundary between the two phases but receive a layer where the probability of the surface position is non-zero. We consider the intermittency function denoting the probability of finding the water phase at a given point of the flow and at the given time. The function α is defined as the ensemble average of the Heaviside function α(x,t) = χ(x,t) where χ(x,t) = 1 if the water phase is present at x and t 0 otherwise (1) The region were 0 < α(x,t) < 1 is called the intermittency region or the surface layer. The width of the surface layer depends on the disrupting action of turbulent eddies on one hand and stabilizing role of gravity and the surface tension on the other hand, cf. Ref. [1]. In our approach these two competing mechanisms are modelled, respectively, by the diffusion and contraction terms in the evolution equation for the function α, cf. [2,3]. Hence, we are not only interested in finding the mean surface position as it is done in the so-called macroscopic approaches, most common up to date, but we stay at the mesoscopic level of description by modelling the evolution of the probability function for the surface position α. In this work we present results of simulations of a surface layer changing in time due to the changing kinetic energy and the characteristic eddy size and compare them with results of apriori tests of two eddies reaching and deforming the surface. The kinetic energy and the eddy size form the apriori simulations are used to determine coefficients in the equation for and calculate its profile within the intermittency region.

Model for the turbulence-interface interactions
If the ensemble averaging is applied to the advection equation of the Heaviside function χ we obtain where u is the ensemble-averaged one-fluid velocity u = u w χ + u a (1 − χ), u w and u a are velocities in, respectively, water and air phase, u is the velocity fluctuation. As it was proposed in Ref. [3] and [4], the unclosed term u · ∇χ can be modelled by two counter-acting terms: diffusion and contraction where coefficients d t and C t depend on turbulence statistics andn = −∇α/|∇α| is a vector normal to the isosurface of α. In the 1D stationary case and for constant coefficients D t and C t Eq.
(2) has an analytical solution whereŷ s denotes a mean position of the interface. In such a case the difference between top and bottom limits of the intermittency layer defined as the region where the probability of finding the surface equals 99.6% reads In the case of the Gaussian error function such limits are obtained for t − b =ŷ s ± 3σ where σ is the standard deviation of the corresponding Gaussian distribution. We proposed in Ref. [3] to model the coefficient D t as the product of the velocity scale q and the characteristic length scale of the eddy L D t = c D qL. The coefficient C t = c C q and the ratio c D /c C can be estimated from the formula (4) such that the proper width of the intermittency region could be obtained.

Test case vortex-interface interactions
In order to perform apriori tests we used the second order finite-volume code FASTEST with the interface capturing volume of fluid method based on the compressive high resolution scheme methodology, see Ref. [4,5]. The aim of the test case presented in figure 1 is to reconstruct the temporal evolution of a two-dimensional gravitational wave generated by the velocity field of two counter rotating vortices initially localized beneath interface. To satisfy flow conditions in a regime of scarified flow the size of the computational domain is chosen to be 20 × 20m. The air-water interface (ρ a = 1.205kg/m 3 , µ a = 1.8 · 10 −5 kg/(ms), ρ w = 998kg/m 3 , µ w = 8.8 · 10 −4 kg/(ms), g = −9.81m/s 2 ) is initially flat and located at y = 15m. The initial intensities equal Ω = 0.511/s for the first and Ω = −0.511/s for the second vortex, (x 0 , y 0 ) is the position of the center of the vortex, x 0 = 9.5m for the first and x 0 = 10.5m for the second vortex, y 0 = 12.5m. The length scale of the eddy is estimated as L = r 0 = 1m. These two parameters locate the present simulation in the so-called scarified regime, cf. Ref. [1], i.e. such flow where the surface is deformed, but not broken by the turbulent eddies. Further details of the simulations are given in Ref. [3]. The computational domain is discretized with about 6 · 10 6 CVs. The CFL condition is kept below 0.5 during whole simulation what requires the times step size ∆t = 10 −4 s. Such mesh resolution assures up to 32 control points in the intermittency region where averaging is performed and the interface-turbulence interaction model is verified. The vorticity magnitude of the eddies approaching the surface and the corresponding deformation of the surface are presented in figures 1 and 2.   The statistics like the intermittency function or the kinetic energy were calculated by averaging in x direction. As the eddies approach the interface they deform it and in the averaged picture the profile of α function is smoothed creating a surface region between the top and the bottom position of the surface, cf. figure 3. In this stage of the flow the width of the intermittency region increases, hence in the following we call it the diffusion stage. As the vortex structure reaches the surface it splits into two eddies travelling in opposite directions. This is connected with the decrease of the intermittency region width caused on one hand by the decrease of the characteristic length scale of the eddy structure and, on the other hand by the stabilizing role of gravity. We will further refer to this phase of the flow calling it contraction stage.

Numerical solution of the α-equation
Modeling of the evolution of the intermittency region in space and time introduce additional difficulties to the numerical solution of the model Eq. (2). Beside typical problems with the advection of the function α = χ known from the interface capturing VOF methods (e.g. restrictions in time step size due to the CFL condition or using flux limiter to avoid Gibbs phenomenon) the new length and time scales are introduced to the solution. The new length scale representing the characteristic thickness of the intermittency region is given by Eq. (4). The new time scale τ ∼ D t /C 2 t is the characteristic time of the intermittency region evolution. In our numerical procedure both new flow field features must be accounted for. The first one by providing the appropriate mesh resolution in the vicinity of the interface, the second one by introduction of the suitable time integration procedure that can handle variation of the intermittency region thickness in time. The solution to Eq. (2) is implemented into the second- order accurate finite volume flow solver FASTEST and is related to the previous implementation of the interface capturing VOF method see Ref. [4] and [5]. The interface capturing VOF method in FASTEST is based on the compressive high resolution scheme M-CICSAM that keeps the interface sharp due to the downwind differencing scheme present in the discretization of the convective term. The numerical method used in the present work to solve Eq. (2) is based on a different idea. It is built around the numerical techniques that are employed in the diffuse interface methods for example presented in Ref. [6]. In such an approach the thickness of the interface between two fluids, cf. Eq. (4) is proportional to where ∆x is the mesh size. We notice that Eq. (5) provides the condition that allows to determine the number of grid points in the vicinity of the interface, when model coefficients D t , C t are given. To keep the interface thickness constant during advection, the additional equation is introduced to reinitialize function α. We notice that reinitialization of the mass conserving level-set function α in Ref. [6] has only numerical meaning i.e. it should keep thickness of the diffuse interface constant and equal to ε during its advection. In the present work we will introduce a physical interpretation to the reinitialization procedure.
Let us split the solution of the Eq. (2) into two steps. In the first step we carry out the advection of the α function ∂ α ∂t in the averaged velocity field u obtained from the RANS solver. In this step we use implicit time and space discretization, where the second order TTL scheme is used for discretization in time and the second order MUSCL flux limiter is used within a deferred correction to keep α bounded function. In the second step, the evolution of the intermittency region must be reconstructed. Here we introduce (reinitialization) equation where the inner artificial time τ ∼ D t /C 2 t is associated with the characteristic time scale of the intermittency region evolution, V D is the diffusion and V C is compression velocity of the intermittency region. Additionally, we make an assumption that coefficients D t , C t remain constant on a single time step ∆t as they are determined from the turbulence model integrated in time along with the Navier-Stokes equations. Eq. (7) is solved in time τ using explicit low-storage, third order accurate TVD Runge-Kutta scheme. The spatial discretization is carried out in framework of the finite volume method and therefore the discretized Eq. (7) in the control volume P becomes where central differencing interpolation is used to obtain values at the faces f of the control volume. If the numerical procedure presented above is used only to keep the interface compact D t = ∆x/2V D and C t = V C , see Eq. (7). Moreover, it is assumed that V D = V C = const where constant value is selected to be u max to increase the rate of convergence in inner time τ. Since in our model in general V D V C in what follows we will discuss different regimes of the solution of Eq. (7) dependent on the sign of the velocity V = V D −V C . Let us first define second order norm used to measure distance between solutions in time τ where subscript i denotes inner time iterations within time τ. Three situations must be considered: In this case thickness of the intermittency region remains constant since the RHS in Eq. (7) is zero. We choose number of the inner time steps in the reinitialization procedure N τ = 4 with inner time step ∆τ = D max t /(2C min t ) 2 where D max t is maximal and C min t is minimal value of the coefficients in the whole computational domain. We checked the 4 inner time iterations are sufficient to achieve constant interface thickness during its advection by Eq. (6) see also Ref. [6].  7) can occur as is discussed in Ref. [6]. Let us assume that thickness of the interface on the n-th time step ∆t and i-th inner time step ∆τ is d n i = 24/ √ 2πD n t /C n t . If before next inner iteration of Eq. (8) we find that d n i+1 = d n i + V ∆τ < ε than we set V D = V C = const and N τ = 4. If it is not the case Eq. (8) is solved until the prescribed tolerance is achieved i.e. L 1 ≤ 10 −3 ∆τ where ∆τ = D t /(4C 2 t ). To assess the performance of the reinitialization algorithm with variable ε that should mimic the modification of the intermittency region thickness due to the presence of the turbulent structures we set thickness of the interface as a function of spatial variable. Let us introduce a test case where we set the distribution of the thickness to be a continuous function The initial interface profile is set to ε ini = ∆x in the whole computational domain. Results are presented in figure 4. It can be noticed that the distribution of the function α is well reconstructed. We found out that about two hundred inner time iterations are required in order to obtain convergence of the results up to tolerance L 1 ≤ 10 −5 ∆τ.  Width of the intermittency layer obtained from simulations with D t = q(t)L(t) and C t = const (line with symbols) in comparison to apriori results.

Validation tests
In this section a numerical solution of Eq. (2) using the method described in section 4 is presented. The diffusion and contraction coefficients were provided every 0.1s from the apriori simulations described in section 3 and Ref. [3]. The model transport equation (2) was solved in artificial time τ until a convergent solution was obtained. To get reasonable results for the intermittency layer, the condition D t /C t ≥ ∆x/2 should be satisfied. On the considered finest grid of 64 cells this condition was first fulfilled at t 0 = 5.6s and we took the profile of α from the apriori test at t 0 = 5.6s as the initial condition. Top and bottom heights of the intermittency layer were calculated at α top = 0.998, α bottom = 0.002, which corresponds to the 99, 6% probability of finding the surface within the intermittency zone. As stated above, values for D t and C t are prescribed by turbulent statistics obtained in the apriori test. Results for three different possibilities in choosing these statistics are presented in the following. In general, we use a turbulent intensity q = 2k/3, where k denotes the turbulent kinetic energy, and a turbulent length scale L to define D t . The turbulent kinetic energy is calculated in two ways. In the first approach we use the one-fluid-velocity fluctuations, leading to In the second approach we calculate k from the phase-weighted velocity fluctuations of the water phase, which yieldsq = 2k/3, wherek Weighted variables have also been proposed for modelling turbulent intermittent flows, cf. Ref.
[7], which show mathematical similarities with the present model for the evolution of α. The turbulent length scale L is considered to be the eddy size, i.e. L(t < 12.5s) = 2m as long as both eddies are rising to the surface and L(t > 19.9s) = 1m once the vortex structure has reached the surface and split into two eddies travelling parallel to the surface. Between 12.5s < L < 19.9s, L(t) is defined by a linear interpolation. We first averaged the kinetic energy from Eq. (11) in y direction over the intermittency layer and calculated the corresponding intensity q(t) and the coefficient D t = c D q(t)L(t). The model constant c D = 4.125 · 10 −5 was evaluated from the apriori test, cf. Ref. [3]. The antidiffusion coefficient C t = 1.915 · 10 −3 m/s was kept constant thorough the simulations. Its value was determined from Eq. (4) considering the width of the intermittency layer at t = 25.5s, which, from the apriori tests, is known to be t − b ≈ 0.04m. Simulations with these coefficients yield reasonable width of the intermittency layer, cf. figure 5. Figure 6 presents the development of top and bottom limits of the surface layer. It can be seen that for D t being a function of time only, a solution symmetric with respect to the mean surface position is obtained. In fact, top and bottom limits are non-symmetric especially during the contraction stage. We attempted to improve the results by considering y-dependence, i.e., calculating the intensity of fluctuations from Eq. (11), keeping the same value of C t . Unfortunately, the agreement between the results presented in figure 6 is not satisfactory. For the third approach presented here, the diffusion coefficient is defined by D t = c Dq (y,t)L(t), calculatingq(y,t) from Eq. (12), while again C t is not changed. In figure 7 the evolution of top and bottom limits in time is presented. It can be seen that the model accounts very well for the position of the intermittency regions upper limit, i.e. the limit of the air-phase side. At the diffusion stage, the width of the intermittency layer is overpredicted, while at contraction stage it is underpredicted. Profiles of α at several times are shown in figure 8. The results generally mimic the apriori α-profiles, however, the concave sections especially at t = 20s, 25s, 14.96 < y < 15m, are not reproduced. Obtaining such reasonable results from defining D t by the phase-weighted turbulent kinetic energy shows that the influence of the water phase on the evolution of the intermittency layer is stronger than the influence of the air phase. The results may be improved by choosing C t = C t (α) as stated in Ref. [3], which agrees with our results presented here.

Conclusions
In the present work evolution of the intermittency region where the probability of the surface position is not zero has been modelled by solving numerically α equation (2). Equation (2)     terms in comparison to the standard formulas used to calculate the position of the interface in VOF methods, namely, the diffusion term due to the disturbing action of the turbulent eddies and contraction due to the stabilizing role of the gravity and surface tension. Numerical difficulties connected with the solution of the α equation were discussed and an appropriate numerical method to solve this equation has been proposed. Next, the numerical scheme has been used to model evolution of the surface layer changing its width due to the interaction with vortices. Satisfactory results of the surface layer width t − b changing in time have been obtained. In order to improve results for the non-symmetric behaviour of the top and bottom limits of the surface layer phase-filtered kinetic energy has been used to determine the model coefficient D t Results of numerical solution of Eq. (2) are presented for the first time here. An on-going further work is connected with the implementation of the model for turbulence-interface interactions into the Navier-Stokes solver in the two-way coupling approach and proposing a closure for the Reynolds-stress terms in the vicinity of the surface.