Numerical simulations of high density ratio lock-exchange flows

In this paper direct numerical simulations of exchange flows of large density ratios are presented and are compared with experiments by Gröbelbauer et al. fJ. Fluid Mech. 250, 669 s1993dg. These simulations, which make use of a dynamic mesh adaptation technique, cover the whole density ratio range of the experiments and very good agreement with the experimental front velocities and the Froude number variations is obtained. Moreover, in order to establish more definitely the Froude number dependency on density ratio, the simulations were carried up to ratios of 100 compared with 21.6 accessible in experiments. An empirical law for the dense front Froude number as a function of the density parameter is proposed. The mathematical difficulty of the problem is discussed. This difficulty arises because, when the density ratio is large, the existence of a solution is dependent on a compatibility condition between the diffusion and viscous terms model. Moreover, a specific numerical technique is required to treat the finite, nonuniform divergence of the mass-averaged velocity field described by the continuity equation. © 2005 American Institute of Physics . fDOI: 10.1063/1.1849800 g


I. INTRODUCTION
Numerical simulations of gravity driven flows are relatively rare compared with the number of experiments which considered various aspects of gravity currents and of density intrusions. 1 Recent numerical simulation 2,3 of gravity currents are limited to small density differences where the Boussinesq approximation is applicable. 4In certain geophysical flows, such as avalanches or pyroclastic flows, and in industrial applications related with heavy gases, the density change across the current fronts is, however, no longer small.Since theoretical models or experimental results which hold for small density ratios can, in general, not be extrapolated to these flows, large density ratio flows need specific attention.
Direct numerical simulations of gravity currents of large density ratios seem to be nonexistent.Most of the experiments are also limited to low density ratios because these were mostly performed with liquids where it is difficult to establish large density ratios.Gröbelbauer et al. 5 conducted lock-exchange flow experiments with gases of density ratios up to 21.6.These flows exhibit some interesting behaviors.In the Boussinesq limit the flow is symmetric and the Froude number varies 6 like Fr= U F / ͱ gh = * / ͱ 2, where h is half the channel depth, U F the front velocity, * = ͱ ͑ d − ᐉ ͒ / ͑ d + ᐉ ͒, and ᐉ and d are the densities of the light and dense fluids, respectively.For large density ratios the exchange flow is asymmetric and asymptotic theories ͑for ᐉ → 0͒ give for the light front 7 Fr ᐉ ϱ =1/ ͱ 2 and for the dense front 8 Fr d ϱ =2 ͱ 2. The experimental results of Gröbel- bauer et al. clearly show this divergence in the respective Froude number values and the results seem to approach the asymptotic limits.In the lock-exchange experiments of Keller and Chyou 9 which cover density ratios up to 10 3 ͑water/air͒ for the light front, the Froude number limit 1 / ͱ 2 is not reached.The reason for this is most likely viscous effects due to the small channel dimensions used in these experiments.
Lock-exchange flows are a good test for direct numerical simulations of flows of miscible, large density difference fluids.Numerical simulations can reach larger values of the density ratio than accessible in experiments, except for the limit case of nonmiscible liquid-gas exchange flows where density ratios of order 10 3 are reached, and can give additional information about the variation of the Froude number and the structure of the intrusion fronts.However, the existence of a solution of the Navier-Stokes equations in these conditions is subject to a condition either on the density ratio compared with Schmidt number, or on the form of the viscous and diffusion terms.Furthermore, due to the unusual condition of a finite, nonuniform divergence of the massaveraged velocity field, a specific technique is needed in order to preserve this existence result when the equations are discretized.Finally, dynamic mesh adaptation is necessary when the density ratio is large.The main purpose of this paper is to derive the appropriate equations and develop a suitable numerical algorithm for treating the non-Boussinesq lock-exchange problem.Comparisons with existing laboratory experiments for density ratios up to 21.6 validate the numerical simulations, which are carried up to ratios of 100 in order to establish more definitely the Froude number dependency on density ratio.
In Sec.II the flow conditions, corresponding to the experiments of Gröbelbauer et al. are presented.The governing equations for large density ratio flows are derived in Sec.III and the numerical algorithm is presented in Sec.IV.The initial behavior of the front in the asymptotic limit of negligible ᐉ is derived in Sec.V and compared with numerical results.The numerical results of the front velocities and the variations of the Froude number are presented in Sec.VI and compared with experiments.

II. LOCK-EXCHANGE FLOW CONDITIONS
When a horizontal channel is divided into two parts by a vertical splitter plate and each chamber is filled with a fluid of different density, an intruding, gravity-driven flow occurs when the splitter plate is removed ͑see Fig. 1͒.It consists in the spread of a dense current of the heavier fluid under the lighter fluid, and of a lighter fluid current above the heavier fluid.This is referred to as lock-exchange flow.In the experiments by Gröbelbauer et al., gases with a density ratio of up to 21.6 were released in an unevenly divided horizontal channel of half-height h = 0.15 m, as shown in Fig. 1.The lock gate could be placed at a distance 20h from the right or left wall, and 10h from the other one.The passage time of either the light or dense front was measured at fixed positions on the horizontal walls of the larger chamber and the Froude number of each front for the various gas pairs was calculated.Table I lists the pairs of gases used and the range of the numerical simulations conducted.The dynamic viscosity of these gases lies between 12.57ϫ 10 −6 Pa s ͑freon 22͒, 18.64ϫ 10 −6 Pa s ͑helium͒, and 21ϫ 10 −6 Pa s ͑argon͒, while the kinematic viscosity ranges from 3.43 ϫ 10 −6 m 2 s −1 ͑freon 22͒ to 1.10ϫ 10 −4 m 2 s −1 ͑helium͒.Thus, it is natural to keep the dynamic viscosity constant in the attempt to reproduce these experiments by numerical simulation.This might be different for liquids.The theoretical formulation below is sufficiently general to include liquids provided the physical properties are known.
In a lock-exchange flow, instabilities could develop in the wall boundary layers at the top and the bottom, at the interface between the dense and light fluids and at the intrusion fronts.Concerning the wall boundary layer, it is well known 10 that for a flow past a flat plate, the boundary layer becomes turbulent for Re x տ 3.5ϫ 10 5 .The Reynolds numbers of the two fronts based on distance x Ӎ Ut are Re x,d = d U 2 t / and Re x,ᐉ = ᐉ U 2 t / .Assuming that both fronts have a velocity U of the same order of magnitude, the Reynolds numbers differ by the density ratio, Re x,d /Re x,ᐉ = d / ᐉ .The dense front boundary layer might reach the critical value, for instance, at x Ӎ 8h for a density ratio of 9.93.For larger density ratios turbulence or at least instabilities could develop at even shorter distances.A transition to Our main concern in this section is to take into account the mutual diffusion of the fluids in the nonhomogeneous, incompressible Navier-Stokes equations.

A. Mass and constituent conservation equations
The mass conservation of constituent i across the surface S of a fixed volume V can be written as where q ˜i is the part of the mass flux which is due to diffusion.Thus ⌽ d and ⌽ ᐉ obey the equations: Because of mass conservation, the left-hand side of Eq. ͑2͒ is necessarily zero.Now, in order for the right-hand side to be zero for arbitrary distributions of the constituents, we need to where D is a reference diffusivity and F is some function of the local composition, as suggested by Joseph and Renardy. 11 nondimensional form, when using these specific fluxes in ͑2͒ and ͑1a͒, the continuity equation and the corresponding equation of the volume fractions are where ReSc= Uh / D is the product of the Reynolds and Schmidt numbers, with U = ͱ ␣gh the terminal velocity of a dense fluid parcel in the light fluid.The variables are nondimensionalized by x = x ˜/ h, u = u ˜/ U, and t = t ˜U / h.Equation ͑3͒, ١ • u 0, is unusual.It arises because of the diffusion between the two species.It is readily seen from Eqs. ͑3͒ and ͑4͒ that when Sc tends to infinity, ١ • u goes to zero.Otherwise, diffusion will result in equal and opposite mass fluxes of constituents d and ᐉ across the boundary of any small volume V͑t͒ entrained by the flow velocity.As a result, since both constituents are incompressible and of different densities, the volume V͑t͒ will vary; giving ١ • u 0. Note that diffusion effects are obviously negligible for Boussinesq conditions, ␣ Ӷ 1.

B. Momentum equation
We can assume that the mixture behaves like a Newtonian fluid, with a dynamic viscosity that may depend on the local composition of the mixture ⌽.Therefore, we write ͑⌽͒ = ͑⌽͒, where is a constant reference dynamic viscosity and a nondimensional function of the composition of the mixture.Denoting Du = ͑١u + ١ u T ͒ / 2, the momentum equation 12 is and here Re= ᐉ Uh / .For lock-exchange flows and most gravity-driven flows, the boundary condition for u is either where n is the wall normal and =2Du − 2 3 ١ • uI ͑no inflow, slip condition͒.Then, for both mechanical and mathematical reasons, the boundary condition for ⌽ will be ١⌽ • n =0.
In Sec.II we have argued that for gases, ͑⌽͒Х1.However, in this case, proofs of existence of a global weak solution 13 are subject to the condition that 2ScϾ ␣, which means that as far as we know the model may be ill-posed in other situations.There is no physical reason for the Schmidt number to behave this way when ␣ varies; indeed, its value remains of order 1 for common gases.In practice, a blow-up of the numerical solution occurred within the relevant time range for lock-exchange flows for ␣ տ 60.
Bresch et al., 14 on the contrary, show that if the relation holds, then the unconditional existence of global weak solutions can be proved.This condition is never satisfied if we choose ͑⌽͒ = 1.If we take a constant kinematic viscosity = / , that is, if ͑⌽͒ =1+␣⌽, then the relation is matched for Sc=1/2, which is close to the actual Schmidt number for gas mixtures, and a diffusivity of the form F͑⌽͒ =1/͑1+␣⌽͒.This form of the mass diffusivity is a common choice, and can be shown to correspond to the case when the molecular diffusivity of species are equal and independent of the local composition of the mixture. 15n numerical simulations, a nonconstant F is nevertheless an additional difficulty, which requires a specific and computationally expensive treatment. 16Thus the numerical simulations presented here were all performed with a constant mass diffusivity ͑F =1͒.This means that condition ͑6͒ was not satisfied in most numerical simulations, but nevertheless the solutions for constant kinematic viscosity ͓͑⌽͒ =1+␣⌽͔ and F͑⌽͒ = 1 were stable in all cases ͑␣ up to 100 was tested͒.Tests were conducted for ␣ = 20.6, with = 1 and =1+␣⌽.The results showed that the choice of F has no effect on the front velocities.
It should be kept in mind that, when ␣ is large, the meaning of the Reynolds number is very different in cases of ͑⌽͒ = 1 and ͑⌽͒ =1+␣⌽.Indeed, suppose two solutions of ͑5͒, one for each choice of .The actual Reynolds number of the light front ͑i.e., that could be calculated a posteriori from measurements of the light front velocity͒ will be the same for both solutions, while the actual Reynolds number of the dense front is ␣ + 1 times larger in the case of constant than in the case of constant kinematic viscosity.This is because the kinematic viscosity of the dense fluid is ␣ + 1 times smaller.The dilemma is that one model is not able to treat density ratios of ␣ տ 60, and the other does, strictly speaking, not conform to the conditions of the experiments considered, but remains stable.
Note also that numerical simulations can be found in literature ͑e.g., Ref. 17͒ which are based on the volumeaveraged velocity v = u + ͑␣ / ReSc͒F͑⌽͒ ٌ ⌽, because this vector field is solenoidal: ١ • v = 0. Nevertheless, this choice introduces additional inertial terms 15 of higher order in the transformed momentum equation, which cannot be neglected when ␣ is large.The problem is not simplified in doing this.

IV. NUMERICAL APPROACH
The large density difference flows considered are composed of intrusion fronts, 18 where density and velocity gradients are locally steep, and of large areas away from these fronts which have a uniform density and small velocity gradients away from the walls.This calls for a method capable of automatic and unconstrained mesh adaptation, since the location of the interface between dense and light parts of the flow is unknown.However, refining the mesh in areas of steep density gradients makes it difficult to control numerical stability conditions such as ʈuʈ⌬t Ͻ⌬x, where ⌬t and ⌬x are the time step and a local mesh-resolution indicator.Thus we use the method of characteristics for the time discretization of the convective part of the equations, which is not subject to such a condition. 19For the space discretization, we have used a finite elements method, for which mesh adaptation based on the error control is well developed and which allows to use the method of characteristics because the approximation of the velocity is a continuous function.A classical choice for solving the Stokes problem is obtained with the Taylor-Hood finite element, 20 which is a piecewise quadratic approximation of the velocity and a piecewise linear one for the pressure.The volume fraction ⌽ is also discretized in a piecewise quadratic functional space.
The discretization we have used is given in more detail in the Appendix, but one technical difficulty specific to highdensity ratio Navier-Stokes equations needs to be pointed out here.The continuity Eq. ͑3͒ is ١ • u =− for some function which is one of the unknowns of the problem.Now if there is no inflow at the boundary, 21 it is clear from the divergence theorem that −͐ ⍀ dx = ͐ ⍀ ١ • udx = 0.In general, this is not true anymore for the numerical approximation h of ͑where h denotes the diameter of the largest element in the mesh͒, and we have ͐ ⍀ h dx of the same order like the numerical error.This is not sufficiently small to guarantee that an approximation u h of u exists such that ١ • u h =− h , and thus the numerical method will break down.In the Appendix we propose an additional projection step which resolves this problem without reducing the quality of approximation, and we show 16 that this is optimal in the sense of a finite elements approximation.
The mesh adaptation is an iterative process: a first guess of the solution at time t n+1 is calculated on a uniform coarse mesh, and is used to generate a new mesh on which a better approximation of the solution can be calculated.When iterated, this procedure reaches a fixed point corresponding to the best approximation space of a given dimension for the solution. 22This process is handled by the mesh generator BAMG ͑Ref.23͒ for both ⌽ and u, using refinement ratios of order 10 3 between the coarsest triangle size and the finest one.Figure 2 shows the mesh refined around the vorticity sheets of a dense intruding front.The whole of the finite elements resolution is embedded in the open-source C + + environment RHEOLEF. 24

V. ASYMPTOTIC BEHAVIOR AT THE RELEASE
Following Stoker, 8 who obtained an asymptotic solution for the dam-break flow, we have conducted an analytical study of the onset of the lock-exchange flow in the case when ␣ tends to infinity, noting that, away from the walls, viscous effects are negligible in the limit t → 0. The boundary conditions are shown in Fig. 1͑b͒, and in addition we suppose that the left boundary is at the infinity.Also, we suppose that the side walls of the channel allow a perfect slip and thus that the solution is spanwise invariant ͑in z direction͒.Note that since we neglect ᐉ , only the left part of the domain ⍀ − is considered in the calculation.Because ␣ is then infinity, we do not use the same nondimensional form as in Sec.III, but we use UЈ = ͱ gh.Thus, in Lagrange representation with ͑a , b͒ the coordinates corresponding to the initial positions of the particles, if we denote X͑a , b ; t͒ and Y͑a , b ; t͒ the displacement of the particles and p͑a , b ; t͒ the pressure, the Euler equations can be rearranged such that X tt X a + ͑Y tt + 1͒Y a + ␤p a = 0, ͑7a͒ in which the only dimensional quantities are the pressure and The initial conditions correspond to the gate in Fig. 1 with the fluid at rest, so that a Taylor expansion of the displacements around t = 0 gives X͑a , b ; t͒ = a + ␥t 2 + o͑t 2 ͒ and Y͑a , b ; t͒ = b + ␦t 2 + o͑t 2 ͒, and keeping the O͑t 2 ͒ terms in ͑7͒: ͑8c͒ Summing Eqs.͑8a͒ and ͑8b͒, and taking the derivatives with respect to b and a, respectively, yields

͑9͒
We recognize in Eqs.͑8c͒ and ͑9͒ the Cauchy-Riemann conditions, thus, it is necessary and sufficient that the complex function ␦ + i␥ be an analytic function of a + ib in its domain so that ␥ , ␦ are solutions of the problem.Now we make use of the boundary conditions.Obviously, ␦ vanishes for b = 0 and b = 2.For a = 0, using the free surface condition p = 0, the first-order term in ͑8b͒ gives ␦ = −12, and for a → −ϱ we have ␦ → 0. From Eqs. ͑8c͒ and ͑9͒, we infer that ١␥ • n ‫ץ‬⍀ − =0.
Since the system ͑8c͒ and ͑9͒ implies that ⌬␥ = ⌬␦ =0 in ⍀ − , there cannot be more than one solution for ␦, and ␥ is unique up to a constant.This constant is easy to determine, since there must be no influx from infinity, so ͐ 0 2 ͑a , b͒db tends to 0 when a tends to infinity.
Using the mapping w ¯= cosh͓͑−a + ib͒ /2͔, Stoker exhibits an analytic function which enforces the boundary conditions: Finally we obtain the initial acceleration:

͑11͒
The acceleration is independent of d , but depends only on UЈ 2 / h = g.There is a singularity in the acceleration at the junction points between the free surface and the walls.This of course would be damped by viscous forces, nevertheless we can expect a strong boundary layer at these points.Moreover, since the viscous effects propagate as t, we can compare the velocity profile of the solution of a viscous model with the analytical results outside the boundary layer.In Fig. 3 we have plotted the velocity obtained from asymptotic theory of the a = 0 particles at time t.For comparison, we have included the velocity of the particles at x = 0 at the same instant 25 obtained from a numerical simulation with ␣ = 79.Figure 3 shows that the numerical error is small.

A. Evolution of front positions
In order to validate the numerical simulations presented in this paper, the conditions of the experiments Gröbelbauer et al. were reproduced as closely as possible.No-slip boundary conditions, which are known to be the relevant conditions for gas-solid interfaces, were used for all boundaries, except when specified otherwise ͑see Fig. 9͒.The parameters of seven of these experiments reproduced by numerical simulation are shown in Table I.The characteristic Reynolds number of these flows was calculated with the viscosity ᐉ and density ᐉ of the lighter gas, Re= ᐉ Uh / ᐉ .Since the dynamic viscosity is nearly invariant, the Reynolds number for a given density ratio ␣ is directly proportional to ᐉ and, therefore, differs by nearly an order of magnitude, depending on whether the light gas is air or helium.It was found that for the large density ratio flows presented here, the influence of the characteristic Reynolds number on the front progression remains noticeable up to large values of Re.This is seen in Fig. 4, where the Froude number variation with Reynolds number is shown for two density ratios ␣ ͑␣ = 6.23 and ␣ = 20.6͒for the light and dense fronts.Values obtained by Birman et al. 4 for ␣ = 1.5 are included for comparison.For this reason, additional numerical simulations were carried out for different values of ␣, keeping the kinematic viscosity of the light fluid unchanged ͑equal to the kinematic viscosity of air͒, so that the characteristic Reynolds number in these simulations varies like ͱ ␣.
The flow is two dimensional except for the instability at the leading edge, giving rise to the so-called lobe-and-cleft structure, and, possibly, for the boundary layer instability of the dense intrusion when ␣ is large.The Kelvin-Helmholtz instability at the interface is mainly a two-dimensional process but three-dimensional structures at a smaller scale are known to develop as well.All these three-dimensional mo-tions resulting from instabilities can be assumed to have a negligible effect on the bulk properties and geometry of the exchange flow and the simulations using the Boussinesq approximation conducted by Härtel et al. 2 support this assumption.This justifies to restrict our direct numerical simulations to two dimensions.The advantage of this restriction is its much lower computational cost, thus enabling us to use much finer meshes than in a three-dimensional simulation.
In the agreement between the arrivals of the simulated and the experimental fronts is very good.For the dense front, however, a nearly constant shift in time between the calculated and measured front arrivals is observed ͓Fig.5͑b͔͒.Possible explanations for this time shift are either a large numerical inconsistency in time around t = 0 or a time lag in the measurements.The first hypothesis is eliminated by the asymptotic study of the onset of the lock exchange carried out in Sec.V, showing that the numerical solution does fit the analytical prediction.Thus, we are left to suppose that there is either a uniform time lag in the measured arrival time of the front, which may be due to a detection problem, or that the time shift is due to the opening of the gate.Gröbelbauer et al. claim that its manual opening was fast enough and did not induce a large scatter in their measurements, but their chief concern was the established front velocity and not the initial acceleration.It was pointed out that the front detection probes were located at a distance from the floor or the ceiling amounting to 25% of the total height.The foremost front considered in the numerical results may, therefore, have a consistent lead over the front position detected by the probes.This may not explain the whole difference but would account for part of it.

B. Variation of Froude number with density ratio
In Figs.sults imply that when the Reynolds number is sufficiently large, the dense front can be considered to be nondissipative in the sense of Stoker. 8Therefore, at large Reynolds number friction in the boundary layer must remain negligibly small.In order to clarify the importance of the wall boundary layers at the top and bottom of the channel we performed calculations for the same density ratio ͑␣ = 20.6͒ with slip boundary conditions on the channel walls.The results are presented in Fig. 9 where the nondimensional dense and light front velocities are plotted as functions of nondimensional distance x d .It is seen from this figure that when there is a free slip on the wall the established dense front velocity is found to be practically the same for both viscous models = 1 and =1+␣⌽.Furthermore, a no-slip wall boundary condition has practically no effect on the dense front velocity in the case of the constant dynamic viscosity model ͑ =1͒.On the contrary, for the constant kinematic viscosity model ͑ =1+␣⌽͒, the dense front progression is reduced by friction in the wall boundary layer For the light intrusion front Fig. 9 shows that for the constant dynamic viscosity model the wall boundary condition has practically no effect on the front velocity; the velocity is nearly the same with and without wall slip.The constant kinematic viscosity model does not change the wall conditions but increases the viscosity of the displaced dense fluid.
An interesting feature of the flow is the interfacial instability behind the two fronts exhibited by the numerical re-sults.Images of the intrusions are shown in Fig. 10 for three different dense front positions and four density ratios.These images show that in the Boussinesq limit ͑␣ = 0.11͒ the flow is practically symmetric and the interfacial instability is located in the central part of the flow.The start-up rolls are also clearly visible.With increasing density ratio, the instability moves more and more to the dense side and even up to the front ͑see images for ␣ = 20.6 and 39͒, which is in agreement with the stability analysis of Benjamin; 7 the light front gets more stable with increasing density ratio.The decrease of the thickness of the unstable interface ͑decrease in size of the Kelvin-Helmholtz billows as well as of the start-up rolls͒ with increasing density ratio is most likely the reason why the limit of Fr d =2 ͱ 2 is approached in spite of dissipation at the interface; as the density ratio goes to infinity, the ratio of energy dissipation rate to the kinetic energy flux of the dense intrusion goes to zero.For this limit to be reached, the dissipation in the wall boundary layer must also be negligible, which is the case for large Reynolds number ͑Fig.9͒ and as long as the boundary layer remains laminar.

C. Stability of the interface
In order to see why the interface of the light intrusion is more stable, it is of interest to evaluate the interfacial Richardson number Ri= g⌬ i / i ␦ i / ͑⌬ i U͒ 2 , where ␦ i is the interfacial shear layer thickness, ⌬ i U and ⌬ i are, respectively, the velocity and density changes across the shear layer and i is the mean interfacial density.Behind the light front,

which gives
Ri ᐉ ϳ ␦ i / h when C 1 ϳ 2. Since ␦ i / h is of order 10 −1 , with ␦ i increasing somewhat with the density ratio, Ri ᐉ is of order 10 −1 or less.The light intrusion interface should, therefore, Ri d tends to ␦ /4h.The dense intrusion interface becomes more unstable as the density ratio increases.This is also in agreement with Benjamin's stability analysis. 7Furthermore, the coherent structures move more and more with the dense front velocity with increasing density ratios, with the tendency of the structures to move closer to the front as seen in Fig. 10.
Concerning the diffusion, it is of interest to point out that the assumption that F is a constant ͑see Sec.III͒ overestimates the diffusion of light fluid into the dense one, so that the density gradient is reduced.This has, however, little effect on the value of the Richardson number, hence the interfacial instability, because ⌬ i / i ϳ 1. Simulations with F͑⌽͒ =1/͑1+␣⌽͒ for the case ␣ = 20.6 support this conclusion.

VII. CONCLUSIONS
The direct numerical simulations presented in this paper are, to our knowledge, the first simulations of exchange flows of miscible fluids of very large density ratios.The difficulty of the numerical simulation of such flows are exposed and an appropriate numerical scheme is designed.A finite element discretization is used, allowing a dynamic mesh adaptation which is an essential feature in the simulations of this type of flow.The results concerning front velocities and the related Froude number variation with density ratio are in good agreement with the experiments by Gröbelbauer et al. 5 which covered density ratios d / ᐉ ഛ 21.6.In addition, the numerical simulations were extended to density ratios of 100 and allowed to establish more definitely the dependency of the Froude numbers Fr d and Fr ᐉ on the density parameter * .A different, empirical law for the variation of the Froude number of the dense front with the density parameter is proposed.
It is found that the two fronts have a different sensitivity with respect to the viscosity model used.While the light front requires a constant dynamic viscosity model which corresponds to the physical properties of the fluids, the dense front is also fairly well simulated with a constant kinematic viscosity model.An explanation for this behavior is proposed which relies on the wall boundary layer in the case of the dense front and on the effective viscosity of the displaced dense fluid by the light front.
Because of wall friction and interfacial instability the intrusions are strictly speaking always dissipative.Nevertheless, Fig. 11 indicates that when ␣ is small ͑␣ ഛ 0.5͒, both currents would be loss free in the sense of Benjamin 7 and of Keller and Chyou; 9 the current depth is equal to h ͑half the channel height͒.At large values of ␣, the light current continues to occupy close to half the channel depth ͓Fig.11͑b͔͒ and when the Reynolds number is sufficiently large the lossfree Benjamin limit Fr ᐉ ϱ is approached; the interfacial instability is inhibited and the friction in the boundary layer is negligible.On the other hand, the dense current decreases in height and approaches the loss free Stoker solution Fr d ϱ =2 ͱ 2. This means that when the Reynolds number is large the losses due to boundary layer friction and interfacial instability are also negligibly small in the dense current.
+ ⌬t ; t n ͒, then we can define an implicit Euler scheme between t n and t n+1 = t n + ⌬t using this equality.
We cannot apply directly ͑A1͒ since we have used the unknown velocity u͑x , t + ͒ for ͓0,⌬t͔, while we only know u͑x , t͒, but we can calculate This does not affect the order of approximation in ͑A2͒.Using this, Eq. ͑4͒ yields a Poisson-like, classical problem, and Eqs.͑3͒ and ͑5͒ a Stokes-like problem.

Semidiscrete algorithm
We will restrict ourselves here to the case of closed boundaries ͑such that ͉u͉ ‫ץ‬⍀ = 0 and ‫ץ‬⌽ / ‫ץ‬n =0͒, which is not very stringent since many variable-density problems occur in such configurations.A slip condition would be a straightforward extension of this scheme, but would make the notations superfluously complicated.Thus we will search the solution ͑⌽ , u , p͒ in V ϫ V 0 d ϫ Q, with V = H 1 ͑⍀͒, V 0 = H 0 1 ͑⍀͒, and Q = ͕q L 2 ͑⍀͒ , ͐ ⍀ qdx =0͖. is an intermediate variable in V which stands for − ١ • u.
The variational formulation is written in terms of the multilinear forms: Now we discretize the problem by choosing finite element spaces V h and Q h for the approximation of V and Q.
Step ,q h ͒ = ͑ h n+1 ,q h ͒ .∀ q h Q h ͑A6b͒ Step 1 of the algorithm is more complicated than it appears if one considers that we use unstructured meshes with strong local refinements.This means that the knowledge of the coordinates of X n ͑x͒ does not give directly the element K of the mesh it belongs to, and an efficient search algorithm is necessary to determine it.Indeed, if N denotes the number of elements in our mesh, the search algorithm will be used for each degree of freedom in the mesh, that is, O͑N͒ times per time step.We propose an algorithm which allows to keep the overall cost of a time step in O͑N ln N͒, and consists for a given mesh in sorting its elements in a localization tree of depth ln N, which allows a O͑ln N͒ localization for each degree of freedom.
Step 2 is then a classical elliptic equation to solve, a multifrontal LU-type factorization is used.
Step 3 is straightforward, but as shown in Sec.IV, it does not yield an element of Q, and thus in general the Eq.͑A6b͒ has no solution if h n = ⌫ h n .Thus Step 4 performs an orthogonal projection of ⌫ h n onto Q.If the Babuška-Brezzi inf-sup condition holds, this is enough to ensure that Eq. ͑A6b͒ admits solutions.Moreover, this projection preserves the error because ⌫ h n can be shown a good approximation of ͑•,t n ͒ which is an element of Q.
In Step 5 remains a Stokes-like problem, with the difference that the right hand side of Eq. ͑A6b͒ is not zero.We use an augmented Lagrangian technique with a Uzawa iterative algorithm for problem ͑A6͒ as done in Ref. 26.
In Ref. 16 we prove that this scheme yields optimal error bounds ʈu − u h ʈ V + ʈ⌽ − ⌽ h ʈ V ഛ C͑h 2 + ⌬t͒ and that for any ജ 0, for a sufficiently fine mesh and time step we have − ഛ⌽ h n ͑x͒ ഛ 1+ for any x ⍀ and t n ͓0,T͔.We also explain the difficulty of alternatives to the projection step 4.

Fick's law
governs the diffusive fluxes of one fluid into the other with q ˜d =−D dᐉ ͑⌽ d ͒١ ˜⌽d and q ˜ᐉ =−D ᐉd ͑⌽ d ͒١ ˜⌽ᐉ = D ᐉd ͑⌽ d ͒١ ˜⌽d , where the D ij coefficients may depend on the local composition ⌽ d of the mixture.Since ⌽ d + ⌽ ᐉ =1 we can use only one volume fraction ⌽ = ⌽ d .Now, if we sum ͑1a͒ and ͑1b͒ multiplied, respectively, by d and ᐉ , we get

FIG. 2 .
FIG. 2. Local zoom in domain ⍀showing ͑a͒ the nondimensional vorticity and ͑b͒ the mesh used for its calculation; dense intruding front for ␣ = 1.99 at nondimensional time t =6.
FIG. 6. Froude number of the light front Fr ᐉ vs * in experiments and numerical simulations for both viscosity models.+, Experimental values; ‫,ؠ‬ numerical simulations with constant dynamic viscosity model ͑ =1͒ and Re= air h ͱ ␣gh / air ; ᮀ, numerical simulations with constant kinematic viscosity model ͑ =1+␣⌽͒ and Re= air h ͱ ␣gh / air ; ᭝, numerical simulations with constant dynamic viscosity model ͑ =1͒ and Re = He h ͱ ␣gh / He .Error bars for the experimental val- ues represent the discrepancies found between Figs. 2 and 6 in the paper by Gröbelbauer et al. ͑Ref.5 ͒ −•−, joins the theoretical limits for * = 0 and * = 1 according to Fr ᐉ = * / ͱ 2.

TABLE I .
Values of the density parameter * and Reynolds numbers in the experiments of Gröbelbauer et al.͑Ref.5͒ and in the numerical simulations presented here.