# Numerical simulations of high density ratio lock-exchange flows

Jocelyn Etienne, Emil J. Hopfinger, Pierre Saramito

► To cite this version:

Jocelyn Etienne, Emil J. Hopfinger, Pierre Saramito. Numerical simulations of high density ratio lock-exchange flows. Physics of Fluids, 2005, 17 (3), pp.036601:1-12. 10.1063/1.1849800. hal-00222927

HAL Id: hal-00222927

https://hal.science/hal-00222927v1

Submitted on 5 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Numerical simulations of high density ratio lock-exchange flows

Jocelyn Étienne<sup>a)</sup>

LMC-IMAG, Boîte Postale 53, 38041 Grenoble, France

Emil J. Hopfinger

LEGI-CNRS/INPG, Boîte Postale 53, 38041 Grenoble, France

Pierre Saramito

LMC-IMAG, Boîte Postale 53, 38041 Grenoble, France

(Received 18 May 2004; accepted 30 November 2004; published online 28 January 2005)

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.* [J. Fluid Mech. **250**, 669 (1993)]. 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.

[DOI: 10.1063/1.1849800]

## 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.¹ Recent numerical simulation<sup>2,3</sup> of gravity currents are limited to small density differences where the Boussinesq approximation is applicable.⁴ In 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*.⁵ 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⁶ like $Fr=U_F/\sqrt{gh}=Q^*/\sqrt{2}$, where h is half the channel depth, $U_F$ the front velocity, $Q^*=\sqrt{(\rho_d-\rho_\ell)/(\rho_d+\rho_\ell)}$, and $\rho_\ell$ and $\rho_d$ are the densities of the light and dense fluids, respectively. For large density ratios the exchange flow is asymmetric and asymptotic theories (for $\rho_\ell \to 0$) give for the light front $Fr_\ell^\infty=1/\sqrt{2}$ and for the

dense front⁸ $Fr_d^\infty=2\sqrt{2}$. The experimental results of Gröbelbauer *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⁹ which cover density ratios up to $10^3$ (water/air) for the light front, the Froude number limit $1/\sqrt{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 mass-averaged 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

<footnote><sup>a)</sup>Author to whom correspondence should be addressed. Present address: BP Institute for Multiphase Flows, University of Cambridge, Madingley Rise, Cambridge CB3 0EZ, UK. Electronic mail: jocelyn@bpi.cam.ac.uk</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Lock-exchange flow: (a) experimental setup used by Gröbelbauer *et al.* (Ref. 5), $h=0.15$ m; (b) boundary conditions in the initial configuration used in the asymptotic study (Sec. V).</figcaption></figure>

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 $\rho_\ell$ 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 $\mu$ of these gases lies between $12.57 \times 10^{-6}$ Pa s (freon 22), $18.64 \times 10^{-6}$ Pa s (helium), and $21 \times 10^{-6}$ Pa s (argon), while the kinematic viscosity $\nu$ ranges from $3.43 \times 10^{-6}$ m$^2$ s$^{-1}$ (freon 22) to $1.10 \times 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¹⁰ that for a flow past a flat plate, the boundary layer becomes turbulent for $Re_x \ge 3.5 \times 10^5$. The Reynolds numbers of the two fronts based on distance $x \approx Ut$ are $Re_{x,d} = \rho_d U^2 t / \mu$ and $Re_{x,\ell} = \rho_\ell U^2 t / \mu$. 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,\ell} = \rho_d/\rho_\ell$. The dense front boundary layer might reach the critical value, for instance, at $x \approx 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

<table><caption>TABLE I. Values of the density parameter $\rho^*$ and Reynolds numbers in the experiments of Gröbelbauer *et al.* (Ref. 5) and in the numerical simulations presented here. $\alpha=(\rho_d-\rho_\ell)/\rho_\ell$, $\rho^*=\sqrt{(\rho_d-\rho_\ell)/(\rho_d+\rho_\ell)}$, $U=\sqrt{\alpha gh}$, $Re=\rho_\ell U h / \mu_\ell$, and $\nu_\ell=\mu_\ell/\rho_\ell$.</caption><thead><tr><th>α</th><th>ρ<sup>*</sup></th><th>U(m s<sup>-1</sup>)</th><th>Re</th><th></th><th>ν<sub>ℓ</sub> (m<sup>2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>99.0</td><td>0.990</td><td>12.1</td><td>1.20 × 10<sup>5</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>79.0</td><td>0.988</td><td>10.8</td><td>1.07 × 10<sup>5</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>59.0</td><td>0.984</td><td>9.32</td><td>9.25 × 10<sup>5</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>39.0</td><td>0.975</td><td>7.58</td><td>7.53 × 10<sup>4</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">20.6</td><td rowspan="2">0.955</td><td rowspan="2">5.51</td><td>7.49 × 10<sup>3</sup></td><td>R22 and helium</td><td>1.10 × 10<sup>-4</sup></td></tr><tr><td>5.47 × 10<sup>4</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">8.93</td><td rowspan="2">0.904</td><td rowspan="2">3.63</td><td>4.93 × 10<sup>3</sup></td><td>Argon and helium</td><td>1.10 × 10<sup>-4</sup></td></tr><tr><td>2.92 × 10<sup>4</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">6.23</td><td rowspan="2">0.870</td><td rowspan="2">3.03</td><td>4.12 × 10<sup>3</sup></td><td>Air and helium</td><td>1.10 × 10<sup>-4</sup></td></tr><tr><td>2.44 × 10<sup>4</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>1.99</td><td>0.706</td><td>1.71</td><td>1.70 × 10<sup>4</sup></td><td>R22 and air</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">1.18</td><td rowspan="2">0.609</td><td rowspan="2">1.32</td><td>1.57 × 10<sup>4</sup></td><td>R22 and argon</td><td>1.26 × 10<sup>-5</sup></td></tr><tr><td>1.31 × 10<sup>4</sup></td><td>No experiment</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>0.38</td><td>0.400</td><td>0.75</td><td>7.42 × 10<sup>3</sup></td><td>Argon and air</td><td>1.51 × 10<sup>-5</sup></td></tr><tr><td>0.11</td><td>0.228</td><td>0.40</td><td>4.80 × 10<sup>3</sup></td><td>CO<sub>2</sub> and argon</td><td>1.26 × 10<sup>-5</sup></td></tr></tbody></table>
turbulence would cause a decrease of the front progression because of an increased wall shear stress. In the experiments by Gröbelbauer *et al.* this is not clearly observed.

At the interface between the dense and light fluids, shear-layer instabilities develop which give rise to Kelvin–Helmholtz billows. The smoke visualization by Gröbelbauer *et al.* indicates some instability on the interface especially close to the dense front when the density ratio is large. As long as this instability remains two dimensional, its essential features are accurately captured by the direct numerical simulations.

### III. GOVERNING EQUATIONS

Let us consider an isothermal flow of local density $\rho$ and velocity $\tilde{u}$ in a domain $\Omega$ over a time span $[0, T]$ (the symbol $\sim$ denotes dimensional counterparts of quantities and operators otherwise used in nondimensional form). For a perfect mixture of two incompressible fluids, of density $\rho_d$ (the heavier one) and of density $\rho_\ell$ (the lighter one), the local density is $\rho = \rho_d\Phi_d + \rho_\ell\Phi_\ell$ where $\Phi_d, \Phi_\ell$ are the volumic fraction of the constituents, $\Phi_d + \Phi_\ell = 1$ and both, $\rho_d$ and $\rho_\ell$, are constants. The characteristic density ratio is $\alpha = (\rho_d - \rho_\ell)/\rho_\ell$.

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

$$-\frac{\partial}{\partial \tilde{t}} \int_V \rho_i \Phi_i dV = \int_S \rho_i \Phi_i \tilde{u}_i \cdot ndS + \int_S \rho_i \tilde{q}_i \cdot ndS,$$

where $\tilde{q}_i$ is the part of the mass flux which is due to diffusion. Thus $\Phi_d$ and $\Phi_\ell$ obey the equations:

$$\frac{D\Phi_d}{D\tilde{t}} + \Phi_d \tilde{\nabla} \cdot \tilde{u} = - \tilde{\nabla} \cdot \tilde{q}_d, \quad (1a)$$

$$\frac{D\Phi_\ell}{D\tilde{t}} + \Phi_\ell \tilde{\nabla} \cdot \tilde{u} = - \tilde{\nabla} \cdot \tilde{q}_\ell. \quad (1b)$$

Fick's law governs the diffusive fluxes of one fluid into the other with $\tilde{q}_d = -D_{d\ell}(\Phi_d)\tilde{\nabla}\Phi_d$ and $\tilde{q}_\ell = -D_{\ell d}(\Phi_d)\tilde{\nabla}\Phi_\ell = D_{d\ell}(\Phi_d)\tilde{\nabla}\Phi_d$, where the $D_{ij}$ coefficients may depend on the local composition $\Phi_d$ of the mixture. Since $\Phi_d + \Phi_\ell = 1$ we can use only one volume fraction $\Phi = \Phi_d$. Now, if we sum (1a) and (1b) multiplied, respectively, by $\rho_d$ and $\rho_\ell$, we get

$$\frac{D\rho}{D\tilde{t}} + \rho \tilde{\nabla} \cdot \tilde{u} = \tilde{\nabla} \cdot [(\rho_d D_{d\ell} - \rho_\ell D_{\ell d}) \tilde{\nabla} \Phi]. \quad (2)$$

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 have $\rho_d D_{d\ell} = \rho_\ell D_{\ell d} = \rho_d D F(\Phi)$, where $D$ is a reference diffusivity and $F$ is some function of the local composition, as suggested by Joseph and Renardy.¹¹

In nondimensional form, when using these specific fluxes in (2) and (1a), the continuity equation and the corresponding equation of the volume fractions are

$$\nabla \cdot \mathbf{u} = - \frac{\alpha}{1 + \alpha\Phi} \frac{D\Phi}{Dt}, \quad (3)$$

$$\frac{D\Phi}{Dt} + \Phi \nabla \cdot \mathbf{u} = \frac{1}{\text{ReSc}} \nabla \cdot [F(\Phi) \nabla \Phi], \quad (4)$$

where $\text{ReSc} = Uh/D$ is the product of the Reynolds and Schmidt numbers, with $U = \sqrt{\alpha gh}$ the terminal velocity of a dense fluid parcel in the light fluid. The variables are nondimensionalized by $x = \tilde{x}/h$, $u = \tilde{u}/U$, and $t = \tilde{t}U/h$.

Equation (3), $\nabla \cdot \mathbf{u} \neq 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, $\nabla \cdot \mathbf{u}$ goes to zero. Otherwise, diffusion will result in equal and opposite mass fluxes of constituents $d$ and $\ell$ 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 $\nabla \cdot \mathbf{u} \neq 0$. Note that diffusion effects are obviously negligible for Boussinesq conditions, $\alpha \ll 1$.

#### B. Momentum equation

We can assume that the mixture behaves like a Newtonian fluid, with a dynamic viscosity $\mu$ that may depend on the local composition of the mixture $\Phi$. Therefore, we write $\mu(\Phi) = \eta\lambda(\Phi)$, where $\eta$ is a constant reference dynamic viscosity and $\lambda$ a nondimensional function of the composition of the mixture. Denoting $Du = (\nabla u + \nabla u^T)/2$, the momentum equation¹² is

$$(1 + \alpha \Phi) \frac{Du}{Dt} = - \nabla p + \frac{1}{\text{Re}} \nabla \cdot \left[ \lambda(\Phi) \left( 2Du - \frac{2}{3} \nabla \cdot \mathbf{u} \right) \right] - \frac{1 + \alpha \Phi}{\alpha} e_y \quad (5)$$

and here $\text{Re} = \rho_\ell Uh/\eta$. For lock-exchange flows and most gravity-driven flows, the boundary condition for $\mathbf{u}$ is either $\mathbf{u}|_{\partial\Omega}=0$ (no inflow, no-slip condition) or $\mathbf{u} \cdot \mathbf{n}=0$ and a zero wall friction $\sigma \cdot \mathbf{n} - [(\sigma \cdot \mathbf{n}) \cdot \mathbf{n}] \mathbf{n}=0$, where $\mathbf{n}$ is the wall normal and $\sigma = 2D\mathbf{u} - \frac{2}{3}\nabla \cdot \mathbf{u} \mathbf{u}$ (no inflow, slip condition). Then, for both mechanical and mathematical reasons, the boundary condition for $\Phi$ will be $\nabla\Phi \cdot \mathbf{n}=0$.

In Sec. II we have argued that for gases, $\lambda(\Phi) \cong 1$. However, in this case, proofs of existence of a global weak solution¹³ are subject to the condition that $2\text{Sc} > \alpha$, 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 $\alpha$ 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 $\alpha \ge 60$.

Bresch *et al.*,¹⁴ on the contrary, show that if the relation
$$ \nabla\lambda(\Phi) = \frac{\alpha}{2Sc}(1 + \alpha\Phi)F(\Phi)\nabla\Phi \qquad (6) $$

holds, then the unconditional existence of global weak solutions can be proved. This condition is never satisfied if we choose $\lambda(\Phi)=1$. If we take a constant kinematic viscosity $\nu=\mu/\rho$, that is, if $\lambda(\Phi)=1+\alpha\Phi$, 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(\Phi)=1/(1+\alpha\Phi)$. 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.¹⁵

In numerical simulations, a nonconstant $F$ is nevertheless an additional difficulty, which requires a specific and computationally expensive treatment.¹⁶ Thus 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 [$\lambda(\Phi)=1+\alpha\Phi$] and $F(\Phi)=1$ were stable in all cases ($\alpha$ up to 100 was tested). Tests were conducted for $\alpha=20.6$, with $\lambda=1$ and $\lambda=1+\alpha\Phi$. The results showed that the choice of $F$ has no effect on the front velocities.

It should be kept in mind that, when $\alpha$ is large, the meaning of the Reynolds number is very different in cases of $\lambda(\Phi)=1$ and $\lambda(\Phi)=1+\alpha\Phi$. Indeed, suppose two solutions of (5), one for each choice of $\lambda$. 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 $\alpha+1$ times larger in the case of constant $\lambda$ than in the case of constant kinematic viscosity. This is because the kinematic viscosity of the dense fluid is $\alpha+1$ times smaller. The dilemma is that one model is not able to treat density ratios of $\alpha \ge 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 volume-averaged velocity $v=u+(\alpha/\text{ReSc})F(\Phi)\nabla\Phi$, because this vector field is solenoidal: $\nabla \cdot v=0$. Nevertheless, this choice introduces additional inertial terms¹⁵ of higher order in the transformed momentum equation, which cannot be neglected when $\alpha$ is large. The problem is not simplified in doing this.

## IV. NUMERICAL APPROACH

The large density difference flows considered are composed of intrusion fronts,¹⁸ 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\|_{\Delta t} < \Delta x$, where $\Delta t$ and $\Delta 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.¹⁹ For 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,²⁰ which is a piecewise quadratic approximation of the velocity and a piecewise linear one for the pressure. The volume fraction $\Phi$ 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 high-density ratio Navier-Stokes equations needs to be pointed out here. The continuity Eq. (3) is $\nabla \cdot u = -\chi$ for some function $\chi$ which is one of the unknowns of the problem. Now if there is no inflow at the boundary,²¹ it is clear from the divergence theorem that $-\int_{\Omega}\chi d\mathbf{x} = \int_{\Omega}\nabla \cdot u d\mathbf{x} = 0$. In general, this is not true anymore for the numerical approximation $\chi_h$ of $\chi$ (where $h$ denotes the diameter of the largest element in the mesh), and we have $\int_{\Omega}\chi_h d\mathbf{x}$ 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 $\nabla \cdot u_h = -\chi_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¹⁶ 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.²² This process is handled by the mesh generator BAMG (Ref. 23) for both $\Phi$ 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.²⁴

## V. ASYMPTOTIC BEHAVIOR AT THE RELEASE

Following Stoker,⁸ 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 $\alpha$ tends to infinity, noting that, away from the walls, viscous effects are negligible in the limit $t \to 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 $\rho\epsilon$, only the left part of the domain $\Omega^-$ is considered in the calculation. Because $\alpha$ is then infinity, we do not use the same nondimensional form as in Sec. III, but we use $U' = \sqrt{gh}$. Thus, in Lagrange representation with $(a,b)$ the coordinates corresponding to the initial positions of
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

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 + \beta p_a = 0, \quad (7a)$$

$$X_{tt}X_b + (Y_{tt} + 1)Y_b + \beta p_b = 0, \quad (7b)$$

$$X_a Y_b - X_b Y_a = 1, \quad (7c)$$

in which the only dimensional quantities are the pressure and $\beta=h/\rho_d U'^2$.

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+\gamma t^2+o(t^2)$ and $Y(a,b;t)=b+\delta t^2+o(t^2)$, and keeping the $O(t^2)$ terms in (7):

$$2\gamma(1 + \gamma_a t^2) + (2\delta + 1)\delta_a t^2 + \beta p_a = 0, \quad (8a)$$

$$2\gamma\gamma_b t^2 + (2\delta + 1) + (1 + \delta_b t^2) + \beta p_b = 0, \quad (8b)$$

$$\gamma_a + \delta_b = 0. \quad (8c)$$

Summing Eqs. (8a) and (8b), and taking the derivatives with respect to $b$ and $a$, respectively, yields

$$\gamma_b - \delta_a = 0. \quad (9)$$

We recognize in Eqs. (8c) and (9) the Cauchy-Riemann conditions, thus, it is necessary and sufficient that the complex

function $\delta+i\gamma$ be an analytic function of $a+ib$ in its domain so that $\gamma, \delta$ are solutions of the problem.

Now we make use of the boundary conditions. Obviously, $\delta$ vanishes for $b=0$ and $b=2$. For $a=0$, using the free surface condition $p=0$, the first-order term in (8b) gives $\delta=-12$, and for $a \to -\infty$ we have $\delta \to 0$. From Eqs. (8c) and (9), we infer that $\nabla\gamma \cdot n_d\Omega_d=0$.

Since the system (8c) and (9) implies that $\Delta\gamma=\Delta\delta=0$ in $\Omega^-$, there cannot be more than one solution for $\delta$, and $\gamma$ is unique up to a constant. This constant is easy to determine, since there must be no influx from infinity, so $\int_0^2(a,b)db$ tends to 0 when $a$ tends to infinity.

Using the mapping $\bar{w}=\cosh[\pi(-a+ib)/2]$, Stoker exhibits an analytic function which enforces the boundary conditions:

$$\delta + i\gamma = -\frac{i}{2\pi}\ln\frac{\bar{w}-1}{\bar{w}+1}.$$

Finally we obtain the initial acceleration:

$$2\gamma(a,b) = \frac{1}{\pi}\ln\left(\frac{\cos^2\frac{\pi b}{4} + \sinh^2\frac{\pi a}{4}}{\sin^2\frac{\pi b}{4} + \sinh^2\frac{\pi a}{4}}\right), \quad (10)$$

$$2\delta(a,b) = -\frac{2}{\pi}\arctan\left(\frac{\sin\frac{\pi b}{2}}{\sinh\frac{\pi a}{2}}\right). \quad (11)$$

The acceleration is independent of $\rho_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 $vt$, 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 $\alpha=79$. Figure 3 shows that the numerical error is small.

# VI. RESULTS AND COMPARISON WITH EXPERIMENTS

## 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 $\mu_\ell$ and density $\rho_\ell$ of the lighter gas, $Re=\rho_\ell Uh/\mu_\ell$. Since the
<figure><img src="image_4.png" /><figcaption>FIG. 3. Velocity profile at $x=0$ and different instants of time. —, Numerical results; ---, analytical asymptotic approximation (nondimensional values).</figcaption></figure>

dynamic viscosity is nearly invariant, the Reynolds number for a given density ratio $\alpha$ is directly proportional to $\rho_{\ell}$ 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 $\alpha$ ($\alpha=6.23$ and $\alpha=20.6$) for the light and dense fronts. Values obtained by Birman et al.⁴ for $\alpha=1.5$ are included for comparison. For this reason, additional numerical simulations were carried out for different values of $\alpha$, 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 $\sqrt{\alpha}$.

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 $\alpha$ 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.² 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 Fig. 5, we compare the numerical results obtained with the constant dynamic viscosity model ($\lambda=1$) with the experimental results of Gröbelbauer et al. For the light front

<figure><img src="image_5.png" /><figcaption>FIG. 4. Froude number dependence on the Reynolds number of the light and dense fronts. Froude number of the light front $Fr_{\ell}: -^{\circ}$, $\alpha=20.6$, and $\lambda=1$; $-\Delta^{-}$, $\alpha=6.23$, and $\lambda=1$; $-\square^{-}$, Birman et al. (Ref. 4), $\alpha=1.5$, and $\lambda=1$; $+\alpha\Phi^{-}$. Froude number of the dense front $Fr_{d}: -^{\circ}$, $\alpha=20.6$, and $\lambda=1$; $\triangle$, $\alpha=6.23$, and $\lambda=1$; $\Box$, Birman et al. (Ref. 4), $\alpha=1.5$, and $\lambda=1+\alpha\Phi^{-}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison of the numerical and experimental results (in dimensional units). Plot (a) light front, plot (b) dense front. —, Numerical results for $\alpha=20.6$ and $+$, experimental values. ---, Numerical results for $\alpha=8.93$ and $\times$, experimental values. $-^{\circ}$, Numerical results for $\alpha=1.99$ and $\square$, experimental values. The Reynolds number in the simulations is the same as in the experiments.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Froude number of the light front $Fr_ℓ$ vs $\rho^*$ in experiments and numerical simulations for both viscosity models. +, Experimental values; ◦, numerical simulations with constant dynamic viscosity model ($\lambda=1$) and $Re=\rho_{air}h\sqrt{\alpha g h/\mu_{air}}$; □, numerical simulations with constant kinematic viscosity model ($\lambda=1+\alpha\Phi$) and $Re=\rho_{air}h\sqrt{\alpha g h/\mu_{air}}$; △, numerical simulations with constant dynamic viscosity model ($\lambda=1$) and $Re=\rho_{He}h\sqrt{\alpha g h/\mu_{He}}$. Error bars for the experimental values represent the discrepancies found between Figs. 2 and 6 in the paper by Gröbelbauer et al. (Ref. 5) ---, joins the theoretical limits for $\rho^*=0$ and $\rho^*=1$ according to $Fr_ℓ=\rho^*/\sqrt{2}$.</figcaption></figure>

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. 6 and 7, we compare experimental and numerical Froude numbers, $Fr_ℓ=U_ℓ/\sqrt{gh}$ and $Fr_d=U_d/\sqrt{gh}$, of both the light and dense fronts for different density ratios. We use the same density parameter as the one introduced by Gröbelbauer et al., that is,

$$ \rho^* = \sqrt{\frac{\rho_d - \rho_\ell}{\rho_d + \rho_\ell}} $$

Since the Froude number only accounts for the established velocity of the front, the shift in front arrival between numerical and experimental results has no effect.

It is seen in Fig. 6 that, for the light front, the constant dynamic viscosity model ($\lambda=1$) is in close agreement with the experiments of Gröbelbauer et al. It is interesting to note that both the numerical and experimental results concerning $Fr_ℓ$ deviate from the straight line $\rho^*/\sqrt{2}$. Nevertheless, the numerical results carried to $\alpha=59$ show that when $Re$ is large, $Re=(\rho_{air}\sqrt{gh/\mu_{air}})\sqrt{\alpha}$ (that is keeping the same light fluid), the nondissipative limit $Fr_ℓ^c=1/\sqrt{2}$ is approached. The Froude number values obtained for large $\alpha$ and for Reynolds numbers corresponding to experiments are only slightly

<figure><img src="image_8.png" /><figcaption>FIG. 7. Froude number of the dense front $Fr_d$ vs $\rho^*$ in experiments and numerical simulations for both viscosity models. For symbols see Fig. 6; ---, $2\sqrt{2}[1-(1-\rho^*)^{0.3}]$; ---, $1.8\sqrt{2}[1-(1-\rho^*)^{0.3}]$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8. Correlation law between Froude number of the dense front $Fr_d$ and $1-\rho^*$. (a) ◦, Numerical simulations with constant dynamic viscosity model ($\lambda=1$) and $Re = \rho_{air}h\sqrt{\alpha g h / \mu_{air}}$; ---, $2\sqrt{2}[1-(1-\rho^*)^{0.3}]$. (b) +, Experimental values; □, numerical simulations with constant kinematic viscosity model ($\lambda=1+\alpha\Phi$) and $Re = \rho_{air}h\sqrt{\alpha g h / \mu_{air}}$; △, numerical simulations with constant dynamic viscosity model ($\lambda=1$) and $Re = \rho_{He}h\sqrt{\alpha g h / \mu_{He}}$; ---, $1.8\sqrt{2}[1-(1-\rho^*)^{0.3}]$.</figcaption></figure>

lower. In the same figure, we show the results obtained with the constant kinematic viscosity model ($\lambda=1+\alpha\Phi$); this increases the dense fluid viscosity by the density ratio $\alpha+1$. The results indicate that this model is clearly not consistent with the experiments by Gröbelbauer et al. and neither approaches the asymptotic limit; the velocity is strongly reduced by the increase of viscosity of the dense fluid. The results come actually closer to the experimental observations of Keller and Chyou⁹ where viscous effects are important.

Figure 7 shows that for the dense front the constant kinematic viscosity model ($\lambda=1+\alpha\Phi$) as well as the dynamic viscosity model ($\lambda=1$) for Reynolds numbers corresponding to the experiments of Gröbelbauer et al. fit closely the experiments. The simulations with larger $Re$ ($\rho_\ell = \rho_{air}$) and constant dynamic viscosity give slightly higher values of $Fr_d$ when $\alpha$ is large. The data points can be closely fitted by a power law of the form

$$Fr_d = Fr_d^\infty[1 - (1 - \rho^*)^n].$$

From the logarithmic plot shown in Fig. 8(a), it is seen that the high Reynolds number results fall on a straight line over nearly two decades when $Fr_d^\infty = 2\sqrt{2}$ and $n=0.3$. The lower Reynolds number results, including the experimental points, are better approximated by $Fr_d^\infty = 1.8\sqrt{2}$ and $n=0.3$. These results imply that when the Reynolds number is sufficiently

large, the dense front can be considered to be nondissipative in the sense of Stoker.⁸ Therefore, 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 ($\alpha=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 $\lambda=1$ and $\lambda=1+\alpha\Phi$. 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 ($\lambda=1$). On the contrary, for the constant kinematic viscosity model ($\lambda=1+\alpha\Phi$), the dense front progression is reduced by friction in the wall boundary layer. The constant kinematic viscosity model increases the effective dynamic viscosity, hence decreases the effective Reynolds number, in the dense intrusion boundary layer by $\mu_d=(\alpha+1)\mu_\ell$. Consequently, the wall shear stress increases from $\tau_\mu=1/\text{Re}\partial u_\ell/\partial y$ when $\lambda=1$ to $\tau_\nu=(\alpha+1)/\text{Re}\partial u_\ell/\partial y$ if $\lambda=1+\alpha\Phi$. Thus, it can be concluded that the experimental results of Gröbelbauer et al. are probably affected by non-negligible wall friction when $\alpha$
<figure><img src="image_10.png" /><figcaption>FIG. 9. Nondimensional velocities of the dense and light fronts, respectively, $U_d$ and $U_\ell$ as functions of the nondimensional dense front position $x_d$ for $\alpha=20.6$ and $\text{Re}=\varrho_{\text{air}}h\sqrt{\alpha g h / \mu_{\text{air}}}$. Velocity $U_d$ (upper curves): -- --, no-slip boundary condition and $\lambda=1+\alpha\Phi$ (constant kinematic viscosity); --, no-slip condition and $\lambda=1$ (constant dynamic viscosity); - -, free slip condition and $\lambda=1+\alpha\Phi$; —, free slip condition and $\lambda=1$. Velocity $U_\ell$ (lower curves): - × -, no-slip condition and $\lambda=1$; \*-, free slip condition and $\lambda=1$.</figcaption></figure>

≥ 10. In the calculations with the constant dynamic viscosity model and Reynolds number larger than the experimental values, wall friction is negligible.

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 ($\alpha=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 $\alpha=20.6$ and 39), which is in agreement with the stability analysis of Benjamin;⁷ 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\sqrt{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\Delta_i\rho/\rho_i\delta_i/(\Delta_i U)^2$, where $\delta_i$ is the interfacial shear layer thickness, $\Delta_i U$ and $\Delta_i \rho$ are, respectively, the velocity and density changes across the shear layer and $\rho_i$ is the mean interfacial density. Behind the light front, $\Delta_i U = C_1 U_\ell \approx C_1 \sqrt{gh/2} \rho^*$ and $\Delta_i \rho/\rho_i \approx 2(\rho^*)^2$, which gives $Ri_\ell \sim \delta_i/h$ when $C_1 \sim 2$. Since $\delta_i/h$ is of order $10^{-1}$, with $\delta_i$ increasing somewhat with the density ratio, $Ri_\ell$ is of order $10^{-1}$ or less. The light intrusion interface should, therefore,

<figure><img src="image_11.png" /><figcaption>FIG. 10. Nondimensional vorticity maps for different density ratios and at different stages in the flow development for the constant dynamic viscosity model $\lambda(\Phi)=1$. The dense front positions are taken the same for the three density ratios.</figcaption></figure>
be weakly unstable but gets more stable with increasing $\rho^*$.
This is consistent with Benjamin's stability analysis.⁷ At the
dense side, $\Delta_i U = C_2 U_d = C_2 2\sqrt{2gh(1-\frac{4}{\sqrt{4}}1-\rho^*)}$, with $C_2$ de-
creasing from 2 to about 1.2 as the density ratio increases,
and again $\Delta_i \rho / \rho_i \approx 2(\rho^*)^2$. The dense intrusion Richardson
number is, therefore, $Ri_d=(\rho^*)^2 \delta_i / [4hC_2^2(1-\frac{4}{\sqrt{4}}1-\rho^*)^2]$.
When $\rho^* \ll 1$, we have $C_2 \approx C_1 \approx 1$, giving $Ri_d \sim \delta_i / h$ by
Taylor expansion. As $\rho^*$ goes to 1, $C_2$ grows close to 1 and
$Ri_d$ tends to $\delta / 4h$. The dense intrusion interface becomes
more unstable as the density ratio increases. This is also in
agreement with Benjamin's stability analysis.⁷ Furthermore,
the coherent structures move more and more with the dense
front velocity with increasing density ratios, with the ten-
dency 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) overestis-
mates the diffusion of light fluid into the dense one, so that
the density gradient is reduced. This has, however, little ef-
fect on the value of the Richardson number, hence the inter-
facial instability, because $\Delta_i \rho / \rho_i \sim 1$. Simulations with
$F(\Phi) = 1 / (1 + \alpha \Phi)$ for the case $\alpha = 20.6$ support this conclu-
sion.

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 dif-
ficulty 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 ad-
aptation 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.⁵
which covered density ratios $\rho_d/\rho_l \le 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_l$ on the density parameter $\rho^*$.
A different, empirical law for the variation of the Froude
number of the dense front with the density parameter is pro-
posed.

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 cor-
responds 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. Neverthe-
less, Fig. 11 indicates that when $\alpha$ is small ($\alpha \le 0.5$), both
currents would be loss free in the sense of Benjamin⁷ and of
Keller and Chyou;⁹ the current depth is equal to $h$ (half the
channel height). At large values of $\alpha$, the light current con-
tinues to occupy close to half the channel depth [Fig. 11(b)]

<figure><img src="image_12.png" /><figcaption>FIG. 11. (a,b) Nondimensional vorticity maps for the steady flow (long times). (a) $\alpha=0.11$, Re=4.80×10³ at nondimensional time $t=8$ ($\bar{t}=3$ s); (b) $\alpha=20.6$, Re=5.47×10⁴ at nondimensional time $t=46$ ($\bar{t}=1.3$ s). See Fig. 10 for the gray scale legend. (c) Nondimensional density isolines, $\alpha=20.6$, Re=5.47×10⁴ at nondimensional time $t=46$ ($\bar{t}=1.3$ s).</figcaption></figure>

and when the Reynolds number is sufficiently large the loss-
free Benjamin limit $Fr_\ell^\infty$ is approached; the interfacial insta-
bility 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^\infty$
$=2\sqrt{2}$. This means that when the Reynolds number is large
the losses due to boundary layer friction and interfacial in-
stability are also negligibly small in the dense current.

ACKNOWLEDGMENTS

This work was partially supported by a grant (BQR) of
the INPG. The first author, J.E., acknowledges a Fellowship
from the French Ministry of Education and Research.

APPENDIX: NUMERICAL SCHEME

1. Discretization in time

The method of characteristics consists in approximating
the total derivative $(\partial/\partial t + \mathbf{u} \cdot \nabla)$ by a finite difference in time
along the pathlines of the flow. First we define the pathlines
with a mapping $X(\mathbf{x}, t; t+\tau)$ between the fluid particles lo-
cated at $\mathbf{x}$ in $\Omega$ at time $t$ and the position these reach when
advected by the fluid velocity $\mathbf{u}$ over a time span $\tau$:

$$
X(\mathbf{x}, t; t + \tau) = x + \int_{t}^{t+\tau} \mathbf{u}[X(\mathbf{x}, t; s), s]ds. \quad (A1)
$$

Then it is easily shown that

$$
\begin{align}
& \left( \frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla \right) f(\mathbf{x}, t + \Delta t) \nonumber \\
& \quad = \frac{f(\mathbf{x}, t + \Delta t) - f[\mathbf{X}(\mathbf{x}, t + \Delta t; t), t]}{\Delta t} + O(\Delta t). \tag{A2}
\end{align}
$$

Thus, if we can calculate $X^n$ an approximation of $X(\cdot, t^n)$
+Δt; tⁿ), then we can define an implicit Euler scheme between tⁿ and tⁿ+<fnref n="1" /> = tⁿ + Δ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

$$
X^n = x - \int_t^{t+\Delta t} \mathbf{u}[X(\mathbf{x},t;s),t]ds = X(\mathbf{x},t+\Delta t;t) + O(\Delta t^2).
$$

This does not affect the order of approximation in (A2). Us-
ing this, Eq. (4) yields a Poisson-like, classical problem, and
Eqs. (3) and (5) a Stokes-like problem.

**2. Semidiscrete algorithm**

We will restrict ourselves here to the case of closed
boundaries (such that **u**|<sub>∂Ω</sub>=**0** and ∂Φ/∂**n**=0), which is not
very stringent since many variable-density problems occur in
such configurations. A slip condition would be a straightfor-
ward extension of this scheme, but would make the notations
superfluously complicated. Thus we will search the solution
(Φ,**u**,p) in V×V₀ᵈ×Q, with V=H<fnref n="1" />(Ω), V₀=H₀<fnref n="1" />(Ω), and Q
={q∈L<fnref n="2" />(Ω), ∫<sub>Ω</sub>qdx=0}. χ is an intermediate variable in V
which stands for −∇·**u**.

The variational formulation is written in terms of the
multilinear forms:

$$
a(\Phi, \mathbf{u}, \mathbf{v}) = \frac{1}{\Delta t} [\mathbf{u}, (1 + \alpha \Phi) \mathbf{v}] + \frac{1}{\mathrm{Re}} \left\{ 2[\mathrm{D}\mathbf{u}, \lambda(\Phi)\mathrm{D}\mathbf{v}] - \frac{2}{3} [\nabla \cdot \mathbf{u}, \lambda(\Phi) \nabla \cdot \mathbf{v}] \right\},
$$

$$
b(\mathbf{v}, \mathbf{q}) = -(q, \nabla \cdot \mathbf{v}),
$$

$$
c(\Phi, \psi) = \frac{1}{\Delta t} (\Phi, \psi) + \frac{1}{\text{ReSc}} (\nabla \Phi, \nabla \psi).
$$

Now we discretize the problem by choosing finite ele-
ment spaces $V_h$ and $Q_h$ for the approximation of $V$ and $Q$.

**Algorithm**

Initialization. n=0. Choose Φ₀h some arbitrary function
in V_h, with Φ₀h(x) ∈ [0,1] almost every x ∈ Ω and
∇Φ₀h·n|∂Ω=0, almost every x ∈ ∂Ω, and u₀h in V₀,ₙ.

Loop. n ≥ 0, assuming (Φⁿ, **u**ⁿ) are given.

Step 1. Calculate $X^n(\cdot)$ with

$$
X^n(\mathbf{x}) = \mathbf{x} - \Delta t \mathbf{u}_h^n \left( \mathbf{x} - \frac{\Delta t}{2} \mathbf{u}_h^n(\mathbf{x}) \right). \quad (A3)
$$

Step 2. Find $\Phi_h^{n+1}$ in $V_h$ such that, for all $\psi_h \in V_h$,

$$
c(\Phi_h^{n+1}, \psi_h) = \left( \Phi_h^n \chi^n + \frac{1}{\Delta t} \Phi_h^n \circ X^n, \psi_h \right). \quad (A4)
$$

Step 3. Calculate $\Gamma_h^{n+1} \in V_h$, such that, for all $\psi_h \in V_h$,

$$
(\Gamma_h^{n+1}, \psi_h) = \left( \frac{\alpha \Phi_h^{n+1} - \Phi_h^n \circ X^n}{1 + \alpha \Phi_h^n} \right), \psi_h. \quad (A5)
$$

Step 4. Calculate $\chi_h^{n+1} = \Gamma_h^{n+1} - 1/|\Omega|(\int_{\Omega} \Gamma_h^{n+1} d\mathbf{x})|^{-1}$.

Step 5. Find ($\mathbf{u}_h^{n+1}$, $p_h$) in $V_{0,h}^d \times Q_h$ such that

$$
\begin{align}
& a(\Phi_h^{n+1}, \mathbf{u}_h^{n+1}, \mathbf{v}_h) + b(\mathbf{v}_h, p_h) \nonumber \\
& \quad = \frac{1}{\Delta t} (\mathbf{u}_h^n \circ \mathbf{X}^n, v_h) - \left( \frac{1 + \alpha \Phi_h^{n+1}}{\alpha} e_z, \mathbf{v} \right), \quad \forall \mathbf{v} \in V_{0,h}^d, \tag{A6a}
\end{align}
$$

$$
b(\mathbf{u}_h^{n+1}, q_h) = (\chi_h^{n+1}, q_h). \quad \forall q_h \in Q_h \quad (A6b)
$$

Step 1 of the algorithm is more complicated than it ap-
pears if one considers that we use unstructured meshes with
strong local refinements. This means that the knowledge of
the coordinates of Xⁿ(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 $\chi_h^n = \Gamma_h^n$. Thus Step 4 performs an orthogo-
nal projection of $\Gamma_h^n$ onto Q. If the Babuška-Brezzi inf-sup
condition holds, this is enough to ensure that Eq. (A6b) ad-
mits solutions. Moreover, this projection preserves the error
because $\Gamma_h^n$ can be shown a good approximation of $\chi(\cdot, t^n)$
which is an element of Q.

In Step 5 remains a Stokes-like problem, with the differ-
ence 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<fnref n="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 ex-
plain the difficulty of alternatives to the projection step 4.

<footnote n="1">¹J. E. Simpson, *Gravity Currents in the Environment and the Laboratory*, 2nd ed. (Cambridge University Press, Cambridge, 1997).</footnote>

<footnote n="2">²C. J. J. Härtel, E. Meiburg, and F. Necker, “Analysis and direct numerical simulation of the flow at a gravity-current head: part 1,” J. Fluid Mech. **418**, 189 (2000).</footnote>

<footnote n="3">³F. Necker, C. Härtel, L. Kleiser, and E. Meiburg, “High-resolution simulations of particle-driven gravity currents,” Int. J. Multiphase Flow **28**, 279 (2002).</footnote>

<footnote n="4">⁴In an article by V. K. Birman, J. E. Martin, and E. Meiburg, J. Fluid Mech. (submitted), about the same time as the present article was submitted for publication. Numerical simulations of lock-exchange flows are presented for density ratios up to 5, and Reynolds number of Re≈1.4×10³ (in our definition of Re). No mesh adaptation was used.</footnote>

<footnote n="5">⁵H. P. Gröbelbauer, T. K. Fanneløp, and R. E. Britter, “The propagation of intrusion fronts of high density ratios,” J. Fluid Mech. **250**, 669 (1993).</footnote>

<footnote n="6">⁶C.-S. Yih, *Dynamics of Nonhomogeneous Fluids* (Macmillan, New York, 1965).</footnote>

<footnote n="7">⁷T. B. Benjamin, “Gravity currents and related phenomena,” J. Fluid Mech. **31**, 2 (1968).</footnote>

<footnote n="8">⁸J. J. Stoker, *Water Waves* (Interscience, New York, 1957).</footnote>

<footnote n="9">⁹J. J. Keller and Y. P. Chyou, “On the hydraulic lock-exchange problem,” ZAMP **42**, 874 (1991).</footnote>

<footnote n="10">¹⁰H. Schlichting and K. Gersten, *Boundary Layer Theory* (Springer, Berlin,</footnote>
2000).

<footnote n="11">¹¹D. D. Joseph and Y. Y. Renardy, *Fundamentals of Two-Fluid Dynamics. Part II: Lubricated Transport, Drops and Miscible Liquids*, Interdisciplinary Applied Mathematics Vol. 4 (Springer, New York, 1993).</footnote>

<footnote n="12">¹²Note that the choice of the reference velocity $U = \sqrt{\alpha g/h}$ implies that the Froude number is simply $\sqrt{\alpha}$ in the equations.</footnote>

<footnote n="13">¹³S. N. Antonsev, A. V. Kazhikhov, and V. N. Monakov, *Boundary Value Problems in Mechanics of Nonhomogeneous Fluids* (North-Holland, Amsterdam, 1990).</footnote>

<footnote n="14">¹⁴D. Bresch, E. H. Essoufi, and M. Sy, “De nouveaux systèmes de type Kazhikhov–Smagulov: modèles de propagation de polluants et de combustion à faible nombre de Mach,” C. R. Acad. Sci., Ser. I: Math. **335**, 973 (2002).</footnote>

<footnote n="15">¹⁵J. Camacho and H. Brenner, “On convection induced by molecular diffusion,” Ind. Eng. Chem. Res. **34**, 3326 (1995).</footnote>

<footnote n="16">¹⁶J. Étienne, “Simulation numérique d’écoulements gravitaires à fortes différences de densité. Application aux avalanches,” Ph.D. thesis, INP Grenoble, 2004.</footnote>

<footnote n="17">¹⁷F. Boyer, “Écoulements diphasiques de type Cahn-Hilliard,” Ph.D. thesis, Université Bordeaux I, 2001.</footnote>

<footnote n="18">¹⁸Note that the sharpness of this front will only be governed by convection, and that away from the intrusion area, in shear-flow areas, the interface will be diffuse.</footnote>

<footnote n="19">¹⁹O. Pironneau, *Finite Elements Methods for Fluids* (Wiley, Chichester, 1989).</footnote>

<footnote n="20">²⁰P. Hood and C. Taylor, “A numerical solution of the Navier–Stokes equations using the finite element technique,” Comput. Fluids **1**, 73 (1973).</footnote>

<footnote n="21">²¹If there is some inflow, a similar condition arises which yields the same problem.</footnote>

<footnote n="22">²²H. Borouchaki, P. L. George, F. Hecht, P. Laug, and E. Saltel, “Delaunay mesh generation governed by metric specifications. Part I: Algorithms,” Finite Elem. Anal. Design **25**, 61 (1997).</footnote>

<footnote n="23">²³BANG is a software program by F. Hecht, INRIA, Rocquencourt, France, http://www-rocq.inria.fr/gamma/cdrom/www/bamg</footnote>

<footnote n="24">²⁴RHEOLOF is a software program by P. Saramito, N. Roquet, and J. Étienne, LMC-IMAG, Grenoble, France, http://www-lmc.imag.fr/lmc-edp/Pierre.Saramito/rheolef</footnote>

<footnote n="25">²⁵Therefore, because the Eulerian coordinates of particles at $a=0$ are $(x,y)^T=t^2[\gamma(0,b), \delta(0,b)]^T$, the comparison is only valid up to the first order in $t$.</footnote>

<footnote n="26">²⁶N. Roquet and P. Saramito, “An adaptive finite element method for Bingham fluid flows around a cylinder,” Comput. Methods Appl. Mech. Eng. **192**, 3317 (2003).</footnote>