# Simulations numériques à haute résolution d'écoulements gravitaires à fortes variations de densité

Jocelyn Etienne, Emil J. Hopfinger, Pierre Saramito

► To cite this version:

Jocelyn Etienne, Emil J. Hopfinger, Pierre Saramito. Simulations numériques à haute résolution d'écoulements gravitaires à fortes variations de densité. 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-00137087

HAL Id: hal-00137087

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

Submitted on 25 Apr 2023

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
# Simulations numériques à haute résolution d'écoulements gravitaires à fortes variations de densité

Jocelyn Étienne<fnref n="1" />,<fnref n="3" />, Emil Hopfinger² & Pierre Saramito<fnref n="1" />

## Résumé :

Des écoulements gravitaires de mélanges binaires de fluides miscibles sont simulés par éléments finis couplés à la méthode des caractéristiques (Lagrange-Galerkin), avec adaptation dynamique de maillage. Les simulations d'écoulements d'échange sont en excellent accord avec les experiences (rapports de densité inférieurs a 20) et sont conduites jusqu'à des rapports de 100.

## Abstract :

Gravity currents of binary mixtures of miscible fluids are simulated by a finite elements method, coupled to the method of characteristics (Lagrange-Galerkin), with dynamic mesh adaptation. Simulations of lock-exchange flows are in excellent agreement with experiments (density ratios lower than 20) and are conducted up to density ratios of 100.

## Mots-clefs :

### Éléments finis, maillage adaptatif, mélanges de fluides, écoulements géophysiques.

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]. Numerical simulation [2, 3] of gravity currents are limited to small density differences where the Boussinesq approximation is applicable, and very recent ones [4] are limited to density ratio of 6. 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 can be, however, much larger. 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.

Lock-exchange flows (see Fig. 1) 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 non-miscible 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, non-uniform 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

<footnote n="1">¹ Laboratoire de Modélisation et Calcul (IMAG) et ² Laboratoire des Écoulements Géophysiques et Industriels, Boîte Postale 53, 38041 Grenoble Cedex</footnote>

<footnote n="3">³ Adresse actuelle: BP Institute for Multiphase Flows, University of Cambridge, Cambridge CB3 0EZ, Grande-Bretagne
jocelyn@bpi.cam.ac.uk</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1 – Lock-exchange fbw set-up used by Gröbelbauer et al. [5], h = 0.15 m.</figcaption></figure>

lock-exchange problem, which allows density ratios up to 100 in order to establish more defini-
tely the Froude number dependency on density ratio.

1 Governing equations

For a perfect mixture of two incompressible fluids, of density ρ<sub>d</sub> (the heavier one) and of
density ρ<sub>ℓ</sub> (the lighter one), the local density is ρ = ρ<sub>d</sub>Φ + ρ<sub>ℓ</sub>(1 − Φ) where Φ and 1 − Φ are
the volumic fraction of respectively the dense and light constituents, and both, ρ<sub>d</sub> and ρ<sub>ℓ</sub>, are
constants. The characteristic density ratio is α = (ρ<sub>d</sub> − ρ<sub>ℓ</sub>)/ρ<sub>ℓ</sub>.

The mutual diffusion between the constituents introduces a diffusive flux, which is governed
by Fick's law and gives the constituent evolution equation :

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

which, combined with the mass conservation equation

$$
\frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{u} = 0,
$$

yields a nonzero divergence condition for the velocity:

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

Here 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 and D a reference diffusivity of the pair of fluids. The variables are non-dimensionalized by x = ˜x/h, u = ˜u/U and t = t̃U/h. Equation (2), ∇·u ≠ 0, is unusual. It arises because of the diffusion between the two species. It is readily seen from equations (2, 1) 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.

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 non-dimensional function of the
composition of the mixture. Denoting D**u** = (∇**u** + ∇**u**^T)/2, the momentum equation is:

$$
(1 + \alpha \Phi) \frac{D\mathbf{u}}{Dt} = -\nabla p + \frac{1}{\text{Re}} \nabla \cdot \left[ \lambda(\Phi) \left( 2\mathcal{D}\mathbf{u} - \frac{2}{3} \nabla \cdot \mathbf{u} \mathbf{I} \right) \right] - \frac{1 + \alpha \Phi}{\alpha} e_y \quad (3)
$$
and here $Re = \rho_l U h / \eta$. For lock-exchange flows and most gravity-driven flows, the boundary condition for $u$ is either $u|_{\partial\Omega} = 0$ (no inflow, no slip condition) or $u \cdot n = 0$ and a zero wall friction $\sigma \cdot n - [(\sigma \cdot n) \cdot n]n = 0$, where $n$ is the wall normal and $\sigma = 2D u - \frac{2}{3} \nabla \cdot u \mathbf{I}$ (no inflow, slip condition). Then, for both mechanical and mathematical reasons, the boundary condition for $\Phi$ will be $\nabla\Phi \cdot n = 0$.

For gas pairs, we know that $\lambda(\Phi) \cong 1$. However, in this case, proofs of existence of a global weak solution [6] are subject to the condition that $2Sc > \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 occured within the relevant time-range for lock-exchange flows for $\alpha \gtrsim 60$.

Bresch et al. [7], on the contrary, show that if the relation

$$ \nabla\lambda(\Phi) = \frac{\alpha}{2Sc} (1 + \alpha\Phi) F(\Phi) \nabla\Phi \quad (4) $$

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[8].

Note also that numerical simulations can be found in literature (e.g. [9]) which are based on the volume-averaged velocity $v = u + (\alpha/ReSc)F(\Phi)\nabla\Phi$, because this vector field is solenoidal: $\nabla \cdot v = 0$. Nevertheless, this choice introduces additional inertial terms[8] of higher order in the transformed momentum equation, which cannot be neglected when $\alpha$ is large. The problem is not simplified in doing this.

## 2 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

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Local zoom in domain $\Omega$ showing (a) the non-dimensional vorticity and (b) the mesh used for its calculation; dense intruding front for $\alpha = 1.99$ at non-dimensional time $t = 6$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Froude number of, (a) the light front $Fr_\ell$ and, (b) the dense front $Fr_d$ versus $\varrho^*$ in experiments and numerical simulations for two viscosity models. +, experimental values; ◦, numerical simulations with constant dynamic viscosity model ($\lambda = 1$) and $Re = \varrho_{air}h\sqrt{\alpha gh}/\mu_{air}$; □, numerical simulations with constant kinematic viscosity model ($\lambda = 1 + \alpha\Phi$) and $Re = \varrho_{air}h\sqrt{\alpha gh}/\mu_{air}$; △, numerical simulations with constant dynamic viscosity model ($\lambda = 1$) and $Re = \varrho_{He}h\sqrt{\alpha gh}/\mu_{He}$. Error bars for the experimental values represent the discrepancies found between Figures 2 and 6 in the article by Gröbelbauer et al. [5] In (a) – -, joins the theoretical limits for $\varrho^* = 0$ and $\varrho^* = 1$ according to $Fr_\ell = \varrho^*/\sqrt{2}$; in (b) – -, $2\sqrt{2}(1 - (1 - \varrho^*)^{0.3})$; - - -, $1.8\sqrt{2}(1 - (1 - \varrho^*)^{0.3})$.</figcaption></figure>

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-discretisation of the convective part of the equations, which is not subject to such a condition [10]. For the space discretisation, 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 [11], which is a piecewise quadratic approximation of the velocity and a piecewise linear one for the pressure. The volume fraction $\Phi$ is also discretised in a piecewise quadratic functional space.

The discretisation we have used is given in more detail in [12], but one technical difficulty specific to high-density ratio Navier-Stokes equations needs to be pointed out here. The continuity equation (2) 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
<figure><img src="image_5.png" /><figcaption>FIG. 4 – (a, b and c), non-dimensional vorticity maps: (a), $\alpha = 0.11$, Re = $4.80 \cdot 10^3$ at non-dimensional time $t = 8$; (b,c), $\alpha = 20.6$, Re = $5.47 \cdot 10^4$ at non-dimensional time (b) $t = 8$, and (c) $t = 46$; (d), non-dimensional density iso-lines, $\alpha = 20.6$, Re = $5.47 \cdot 10^4$ at non-dimensional time $t = 46$.</figcaption></figure>

will break down. In [12] we propose an additional projection step which resolves this problem without reducing the quality of approximation, and we show [13] 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 [14]. This process is handled by the mesh generator BAMG [15] 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 [16].

### 3 Results and conclusions

The results concerning front velocities and the related Froude number variation with density ratio are shown in Fig. 3 in terms of the density parameter $\rho^* = \sqrt{(\rho_d - \rho_\ell)/(\rho_d + \rho_\ell)}$, and for Reynolds numbers taken accordingly to the density ratio, $Re = (\rho_\ell\sqrt{gh}/\mu_\ell)\sqrt{\alpha}$, ranging from $5 \cdot 10^3$ ($\alpha = 0.11$) to $1.2 \cdot 10^5$ ($\alpha = 100$). These results are in good agreement with the experiments by Gröbelbauer et al. [5] which covered density ratios $\rho_d/\rho_\ell \le 21.6$ (for a more detailed discussion, see [12]). 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_\ell$ on the density parameter $\rho^*$. A new, empirical law for the variation of the Froude number of the dense front with the density parameter is proposed.

Because of wall friction and interfacial instability the intrusions are strictly speaking always dissipative. Nevertheless, Figure 4 indicates that when $\alpha$ is small ($\alpha \le 0.5$), both currents would be loss free in the sense of Benjamin[17] and of Keller and Chyou[18]; the current depth is equal to $h$ (half the channel height). At large values of $\alpha$, the light current continues to occupy close to half the channel depth (Figure 4c) and when the Reynolds number is sufficiently large the loss-free Benjamin limit $Fr_\ell^\infty$ 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\sqrt{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.

Références

[1] J.E. Simpson. *Gravity currents in the Environment and the Laboratory*. Cambridge University Press, 1997. 2nd ed.

[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. flow topology and front speed for slip and no-slip boundaries. *J. Fluid Mech.*, 418:189–212, 2000.

[3] F. Necker, C. Härtel, L. Kleiser, and E. Meiburg. High-resolution simulations of particle-driven gravity currents. *Int. J. Multiphase Flow*, 28(2):279–300, 2002.

[4] V. K. Birman, J. E. Martin, and E. Meiburg. The non-Boussinesq lock-exchange problem. *J. Fluid Mech.*, accepted in 2004.

[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–687, 1993.

[6] S. N. Antonsev, A. V. Kazhikhov, and V. N. Monakov. *Boundary value problems in mechanics of nonhomogeneous fluids*. North-Holland, Amsterdam, 1990.

[7] 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. Paris, Ser. I*, 335:973–978, 2002.

[8] J. Camacho and H. Brenner. On convection induced by molecular diffusion. *Ind. Eng. Chem. Res.*, 34:3326–3335, 1995.

[9] F. Boyer. *Écoulements diphasiques de type Cahn-Hilliard*. PhD thesis, Université Bordeaux I, 2001.

[10] O. Pironneau. *Finite Elements Methods for Fluids*. John Wiley & Sons, Chichester, 1989.

[11] P. Hood and C. Taylor. A numerical solution of the Navier-Stokes equations using the finite element technique. *Comput. & Fluids*, 1(1):73–100, 1973.

[12] J. Étienne, E. J. Hopfinger, and P. Saramito. Numerical simulations of high density ratio lock-exchange flows. *Phys. Fluids*, 17(2), 2005.

[13] J. Étienne. *Simulation numérique d'écoulements gravitaires à fortes différences de densité*. *Application aux avalanches*. PhD thesis, INP Grenoble, September 2004.

[14] 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. Des.*, 25:61–83, 1997.

[15] F. Hecht. bamg: Bidimensional anisotropic mesh generator. Technical report, INRIA, Rocquencourt, France, 1997. http://www-rocq1.inria.fr/gamma/cdrom/www/bamg.

[16] P. Saramito, N. Roquet, and J. Étienne. rheolef user's manual. Technical report, LMC-IMAG, Grenoble, France, 2003. http://www-lmc.imag.fr/lmc-edp/Pierre.Saramito/rheolef.

[17] T. B. Benjamin. Gravity currents and related phenomena. *J. Fluid Mech.*, 31(2):209–248, 1968.

[18] J. J. Keller and Y. P. Chyou. On the hydraulic lock-exchange problem. *Z. Angew. Math. Phys.*, 42:874–910, 1991.