Stability and convergence for a complete model of mass diffusion

Applied Numerical Mathematics 61 (2011) 1161–1185

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Stability and convergence for a complete model of mass diffusion R.C. Cabrales a,1 , F. Guillén-González b,∗,2 , J.V. Gutiérrez-Santacreu c,2 a b c

Dpto. de Ciencias Básicas, Universidad del Bío-Bío, Facultad de Ciencias, Campus Fernando May, Casilla 447, Chillán, Chile Dpto. E.D.A.N., University of Sevilla, Aptdo. 1160, 41080 Sevilla, Spain Dpto. de Matemática Aplicada I, University of Sevilla, E. T. S. I. Informática, Avda. Reina Mercedes, s/n, 41012 Sevilla, Spain

a r t i c l e

i n f o

Article history: Received 15 November 2007 Received in revised form 1 February 2011 Accepted 6 June 2011 Available online 24 August 2011 Keywords: Three-dimensional Kazhikhov–Smagulov model Density-dependent Navier–Stokes problem Finite elements Stability Convergence

a b s t r a c t We propose a fully discrete scheme for approximating a three-dimensional, strongly nonlinear model of mass diffusion, also called the complete Kazhikhov–Smagulov model. The scheme uses a C 0 finite-element approximation for all unknowns (density, velocity and pressure), even though the density limit, solution of the continuous problem, belongs to H 2 . A first-order time discretization is used such that, at each time step, one only needs to solve two decoupled linear problems for the discrete density and the velocity–pressure, separately. We extend to the complete model, some stability and convergence results already obtained by the last two authors for a simplified model where λ2 -terms are not considered, λ being the mass diffusion coefficient. Now, different arguments must be introduced, based mainly on an induction process with respect to the time step, obtaining at the same time the three main properties of the scheme: an approximate discrete maximum principle for the density, weak estimates for the velocity and strong ones for the density. Furthermore, the convergence towards a weak solution of the density-dependent Navier–Stokes problem is also obtained as λ → 0 (jointly with the space and time parameters). Finally, some numerical computations prove the practical usefulness of the scheme. © 2011 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction 1.1. The model The Kazhikhov–Smagulov equations describe the motion of a viscous, incompressible fluid having two different densities to be subject to a diffusion effect which is modeled by Fick’s law. Assume that the fluid under consideration fills Ω ⊆ R3 a bounded domain with (sufficiently regular) boundary Γ . The evolution of such a fluid is followed during the time interval [0, T ] for 0 < T < ∞. We adopt the convection that bold-face letters denote vectorial elements and use the notation Q = Ω × (0, T ) and Σ = Γ × (0, T ). Then the complete Kazhikhov–Smagulov is given in conservative form by:

 ⎧   1 ⎪ 2 ⎪ ∇ρ ⊗ ∇ρ + ∇ P = ρ f ⎨ (ρ u )t + ∇ · (ρ u − λ∇ ρ ) ⊗ u − λu ⊗ ∇ ρ − μu + λ ∇ ·

ρ

⎪ ∇·u=0 ⎪ ⎩ ρt + ∇ · (ρ u − λ∇ ρ ) = 0

* 1 2

in Q , in Q , in Q ,

Corresponding author. E-mail addresses: [email protected] (R.C. Cabrales), [email protected] (F. Guillén-González), [email protected] (J.V. Gutiérrez-Santacreu). This work was made while the first author was a Postdoctoral student at Dpto. E.D.A.N. This author’s work was partially supported by project MTM2009-12927, Spain.

0168-9274/$30.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2011.06.017

(1)

1162

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

where ρ : Q → R+ is the fluid density, u : Q → R3 is the incompressible (averaged) velocity field, and P : Q → R is the fluid pressure. Moreover, f represents volume external forces that are applied to the fluid, and μ > 0 and λ > 0 stand for the kinematic viscosity and mass diffusion coefficients, respectively. The tensorial product matrix of two vectors a = (ai )ni=1 , b = (b i )ni=1 is denoted by a ⊗ b with coefficients (a ⊗ b)i , j = ai b j . Using in (1) the equalities

    (ρ u )t + ∇ · (ρ u − λ∇ ρ ) ⊗ u = ρ ut + (ρ u − λ∇ ρ ) · ∇ u , and

  −λ∇ · (u ⊗ ∇ ρ ) = −λ(u · ∇)∇ ρ = −λ∇(u · ∇ ρ ) + λ∇ · ρ (∇ u )t , one arrives at the following (non-conservative) formulation:

 ⎧     1 ⎪ t 2 ⎪ ∇ρ ⊗ ∇ρ + ∇ p = ρ f ⎨ ρ ut + (ρ u − λ∇ ρ ) · ∇ u − ∇ · (μ − λρ )(∇ u ) + μ∇ u + λ ∇ ·

ρ

⎪ ∇·u=0 ⎪ ⎩ ρt + u · ∇ ρ − λρ = 0

in Q , in Q ,

(2)

in Q ,

where p = P − λu · ∇ ρ is a new potential function. We complete the model with the boundary conditions

u ( x , t ) = 0,

∂ρ (x, t ) = 0 for (x, t ) ∈ Σ, ∂n

(3)

and the initial conditions

ρ (0) = ρ0 (x),

u (0, x) = u 0 (x) for x ∈ Ω,

(4)

where n(x) the outwards unit normal vector to Ω at the point x ∈ Γ . Throughout this paper, the initial density will be assumed to satisfy:

0 < m  ρ0 (x)  M

for all x ∈ Ω.

(5)

The density-dependent Navier–Stokes problem is formally obtained from system (1) by just choosing λ = 0:

 ⎧  ⎨ ρ ut + (u · ∇)u − μu + ∇ p = ρ f ∇·u=0 ⎩ ρt + u · ∇ ρ = 0

in Q , in Q ,

(6)

in Q ,

jointly with initial conditions (4) and only Dirichlet boundary conditions for the velocity. It is important to observe that the Neumann boundary condition for the density does not now make sense. The Kazhikhov–Smagulov equations (1) can be seen as a regularization of the density-dependent Navier–Stokes equation (6). 1.2. Known results The first authors who dealt with the mathematical analysis for problem (1) in its simplified version, i.e., dropping the λ2 term, were Kazhikhov and Smagulov [13]. They proved, via a semi-Galerkin method, the existence of global-in-time weak solutions, under the hypothesis on the coefficients: λ < 2μ/( M − m), and the existence of local-in-time strong solutions (which is global-in-time for two dimensions). The complete model (1) were first studied by Beirão da Veiga [2] who established the existence of local-in-time strong solutions using linearization and a fixed point argument. Later, under the same method, Secchi [15] proved the existence of global-in-time strong solutions in the two-dimensional domains if λ/μ is small enough. Moreover, he established the asymptotic behavior, as λ → 0, towards global-in-time strong solutions for the density-dependent Navier–Stokes problem. Assuming nonnegative initial density, Guillén-González [7] proved the existence of global-in-time weak solutions as well as the asymptotic behavior, as λ → 0, towards global-in-time weak solutions of the density-dependent Navier–Stokes problem. In [8], an iterative method is used to prove the existence and regularity of strong solutions for (1), obtaining moreover, some convergence rates. With regard to problem (6), existence of global-in-time weak solutions was proven by Lions in [14]. The regularity of solutions was obtained by Antontsev, Kazhikhov, and Monakhov [1]. There are not many numerical schemes to approximate problems (1) and (6). In [9], an unconditionally stable, convergent, linear numerical scheme for the two-dimensional simplified model (1), without λ2 terms, is studied. This scheme consists of C 0 finite-element spatial approximation combined with the linearized backward-in-time Euler method. It is particularly interesting that this time discretization decouples the computation of the discrete density from the velocity–pressure pair. For the three-dimensional case, a conditionally stable, convergent fully discrete scheme is studied in [10]. The main differences

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1163

given in [10] with respect to [9] are to prove an approximate discrete maximum principle for the discrete density and the asymptotic behavior as λ → 0 towards a weak solution of the density-dependent Navier–Stokes problem (6). In [11] a linear fully discrete scheme is analyzed from the point of view of error estimates. First-order error estimates are proven assuming a compatibility condition between the velocity and density spaces related to hypothesis (H4) below. Moreover, the use of regularity which requires an extra compatibility condition at t = 0 for solutions to be approximated is avoided. 1.3. Outline The rest of the paper is divided as follows. In Section 2, we give the main ideas for the mathematical analysis of problem (1). In Section 3, we describe the numerical scheme and announce the two main results. In Section 4, by means of an induction argument, we first prove pointwise estimates for the density, and then energy estimates for the velocity and strong estimates for the density are obtained, using the discrete Laplacian of the density. Afterwards, we show the compactness for the density and velocity in Section 5 and the passage to the limit in Section 6, concluding the proof of Theorem 4. In Section 7, we study the asymptotic behavior as the diffusion parameter goes to zero, proving Theorem 5. Finally, Section 8 is devoted to presenting some numerical experiences. 2. Mathematical analysis of the complete Kazhikhov–Smagulov model We list here some standard notation used throughout the paper. By L p (Ω) and H s (Ω), 1  p  ∞, s = 0, 1, . . . , we denote the classical Lebesgue and Sobolev spaces, respectively. The usual norm in L p (Ω) and H s (Ω) is denoted by · L p (Ω) and · H s (Ω) . The norm and inner product in L 2 (Ω) will be denoted by | · | and (·, ·), respectively. To define the concept of weak solution, we introduce the following space of functions:





H = u ∈ L 2 (Ω): ∇ · u = 0 in Ω, u · n = 0 on Γ ,





V = u ∈ H 10 (Ω): ∇ · u = 0 in Ω ,



L 20 (Ω) = p ∈ L 2 (Ω): Ω



H 2N (Ω)



p (x) dx = 0 ,

∂ρ = ρ ∈ H (Ω): = 0 on ∂Ω, ∂n

2

ρ (x) dx = Ω



ρ0 (x) dx . Ω

It is known by Poincaré’s inequality that u H 1 (Ω) and |∇ u | are equivalent norms on H 10 (Ω). On the other hand, H 2N (Ω)

is an affine space, and ∇ ρ H 1 (Ω) and |ρ | are equivalent semi-norms on H 2N (Ω).

Definition 1. A pair (ρ , u ) is said to be a weak solution of problem (1)–(3)–(4) on (0, T ) if: (a) u ∈ L ∞ (0, T ; H ) ∩ L 2 (0, T ; V ), ρ ∈ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 2N (Ω)), 0 < m  ρ (x, t )  M, ∀(x, t ) ∈ Q . ((b) ∀φ ∈ C 1 ([0, T ]; V ) such that φ( T ) = 0,

T

    − u , ρ φt + (ρ u − λ∇ ρ ) · ∇φ + μ(∇ u , ∇φ) − λ ρ (∇ u )t , ∇φ dt − λ2

0

T 

1

ρ 0

T =

  (ρ f , φ) dt + ρ0 u 0 , φ(0) .

0

(c) The equation of mass diffusion (1)c is satisfied almost everywhere in Q . Definition 2. A pair (ρ , u ) is said to be a weak solution of problem (6) on (0, T ) if it verifies: (a) u ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; V ), ρ ∈ L ∞ ( Q ) with 0 < m  ρ (x, t )  M a.e. (x, t ) ∈ Q . (b) For all φ ∈ C 1 ([0, T ]; V ) with φ( T ) = 0,

T

   − ρ u , φt + (u · ∇)φ + μ(∇ u , ∇φ) dt =

0

(c) For all

T 0

ϕ ∈ C ([0, T ]; H (Ω)) with ϕ ( T ) = 0, 1

1

  (ρ f , φ) dt + ρ0 u 0 , φ(0) .

∇ ρ ⊗ ∇ ρ , ∇φ dt

1164

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

T −

T (ρ , ϕt ) dt −

0

  (ρ u , ∇ ϕ ) dt = ρ0 , ϕ (0) .

0

Here, we only give an outline of the proof of the existence of weak solutions to (1) for the reader’s convenience [15]. Some ideas will be used later. Theorem 3. Let u 0 ∈ H , ρ0 ∈ H 1 (Ω) and f ∈ L 2 (0, T ; L 6/5 (Ω)). If λ/μ is sufficiently small, then there exists a weak solution of problem (1), (3)–(4) on (0, T ). Proof. We proceed formally, assuming (ρ , u ) to be a sufficiently regular solution of (1), (3)–(4). Then the maximum principle applied to equation (1)c , together with (5), leads to

0 < m  ρ ( x, t )  M

a.e. (x, t ) ∈ Q .

(7)

Multiplying (1)c by −λρ , integrating over Ω , integrating by parts in the convective term, and taking into account the interpolation inequality 1/2

∇ ρ L 4 (Ω)  C Ω ρ L ∞ (Ω) |ρ |1/2 ,

(8)

one arrives at

λ

d dt

|∇ ρ |2 + λ2 |ρ |2  C 0 |∇ u |2 .

(9)

Adding the momentum system (1)a by u to the density equation (1)c by equality holds:

1 d

ρ |u |2 dx + μ|∇ u |2 = λ

2 dt





Ω

ρ (∇ u )t : ∇ u dx + λ2

1

ρ

1 u 2

· u, integrated over Ω , the following energy

∇ ρ ⊗ ∇ ρ , ∇ u + (ρ f , u ).

(10)

Ω

Reasoning as in the simplified model (see [13]), one has



ρ (∇ u )t : ∇ u dx = λ

λ Ω

ρ−

M +m



2

(∇ u )t : ∇ u dx  λ

M −m 2

|∇ u |2 .

Ω

To estimate the second term on the right-hand side of (10), we use the interpolation inequality (8), the pointwise estimate for the density m  ρ  M and Young’s inequality, getting

    1 C 2 λ2 λ2  ∇ ρ ⊗ ∇ ρ , ∇ u   ε1 μ λ2 |ρ |2 + |∇ u |2 ,

ρ

ε1 μ

with ε1 > 0 to be chosen later on. The last term of (10) is easily bounded by

(ρ f , u ) 

μ 2

|∇ u |2 + C f 2L 6/5 (Ω) .

Compiling the above estimates into (10), we arrive at

d



ρ |u |2 dx + 2μ|∇ u |2  1 + C 1

dt

λ

μ

+

C 2 λ2

ε1 μ2



μ|∇ u |2 + ε1 μλ2 |ρ |2 + C f 2L 6/5 (Ω) .

Ω

με2 with ε2 > 0 to be chosen later on, this gives us  d d √ λ C 2 λ2 2 2 ε2 μλ |∇ ρ |2 + | ρ u |2 + (ε2 − ε1 )μλ2 |ρ |2 + 1 − C 1 − − ε C 2 0 μ|∇ u |  C f L 6/5 (Ω) . 2

Adding up (11) to (9) multiplied by

dt

Selecting

μ

dt

ε1 = ε2 /2, and ε2 and λ/μ small enough such that C 1 μ +

ε2 μλ

λ

d dt

|∇ ρ |2 +

d √

ε2

dt

2

| ρ u |2 +

λ2 μ|ρ |2 +

μ 2

C2

ε1 μ

( λ )2

ε1 μ

+ ε2 C 0  12 , we get

|∇ u |2  C f 2L 6/5 (Ω) .

Integrating for t ∈ (0, T ) and bounding from below the density, we obtain

(11)

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185





2



2 

ε2 μλ∇ ρ (t ) + mu (t ) +

t 

ε2 2



2



2

μλ2 ρ (s) + μ∇ u (s)



t ds  C

0

1165

   f (s)26/5 L

(Ω)

ds,

(12)

0

On the other hand, the following “fractional in time estimate” holds [1]: T −δ

  u (t + δ) − u (t )2 dt  C δ 1/2 ∀δ ∈ (0, T ). 0

This estimate implies compactness for the velocity u in L 2 (0, T ; L 2 (Ω)) [16]. Then the existence of weak solutions can be deduced in a standard form [1] by using the Faedo–Galerkin method. 2 3. The finite-element approximation In this section we set out our assumptions on the discretization of the Kazhikhov–Smagulov problem. Then we state our main results, Theorems 4 and 5. 3.1. Hypotheses From now on, we assume that Ω is a bounded domain of R3 with a polyhedral boundary and that there  exists a family of triangulations {T h }h>0 of Ω made up of tetrahedra or hexahedra in three dimensions, so that Ω = K ∈Th K . Let

˜ h ⊂ L 2 (Ω) be finite-elements spaces associated to density, velocity and pressure W h ⊂ H 1 (Ω), V h , V˜ h ⊂ H 10 (Ω) and M h , M 0 respectively. Throughout this work we will suppose the following hypotheses:

(H0) Regularity for the data: u 0 ∈ V , ρ0 ∈ H 2N (Ω) with 0 < m  ρ0  M in Ω and f ∈ L 2 (0, T ; L 6/5 (Ω)). Assume λ/μ sufficiently small. (H1) Assume Ω an open, bounded set of R3 , whose boundary is polyhedral and such that the continuous dependencies in H 2 -norm of the Poisson–Neumann problem and in H 2 × H 1 -norm of the Stokes hold (see (25) for a Poisson–Neumann problem). This is verified for example if Ω is convex [5]. (H2) The triangulation of Ω and the discrete spaces verify • the inverse inequalities:

|∇ ρ¯h |  Ch−1 |ρ¯h | ∀ρ¯h ∈ W h , |∇ ρ¯h | L 3 (Ω)  Ch−1/2 |∇ ρ¯h | ∀ρ¯h ∈ W h , ρ¯h L ∞ (Ω)  Ch−1/2 ρ¯h H 1 (Ω) ∀ρ¯h ∈ W h , ∇ ρ¯h L 4 (Ω)  Ch−3/4 |∇ ρ¯h | ∀ρ¯h ∈ W h , • and the interpolation errors:

u¯ − ˜J h u¯ H 1 (Ω) + u¯ − J h u¯ H 1 (Ω)  Ch|u¯ | H 2 (Ω) ∀u¯ ∈ H 2 (Ω) ∩ H 10 (Ω), | p¯ − K˜ h p¯ | + | p¯ − K h p¯ |  Ch| p¯ | H 1 (Ω) ∀ p¯ ∈ H 1 (Ω) ∩ L 20 (Ω), |ρ¯ − I h ρ¯ | + h ρ¯ − I h ρ¯ H 1 (Ω)  Ch2 |ρ¯ | H 2 (Ω) ∀ρ¯ ∈ H 2 (Ω), ρ¯ − I h ρ¯ W 1,3 (Ω)∩L ∞ (Ω)  Ch1/2 |ρ¯ | H 2 (Ω) ∀ρ¯ ∈ H 2 (Ω), ρ¯ − I h ρ¯ W 1,4 (Ω)  Ch1/4 |ρ¯ | H 2 (Ω) ∀ρ¯ ∈ H 2 (Ω), where J h , ˜J h , K h , K˜ h and I h are interpolation operators from H 2 (Ω) ∩ H 10 (Ω) into V h , H 2 (Ω) ∩ H 10 (Ω) into V˜ h ,

˜ h , and H 2 (Ω) into W h , respectively. H 1 (Ω) ∩ L 20 (Ω) into M h , H 1 (Ω) ∩ L 20 (Ω) into M

˜ h, (H3) Inf–sup conditions. There exist β > 0 and β˜ > 0 (independent of h) such that, ∀ p¯ h ∈ M h and ∀¯qh ∈ M

p¯ h L 2 (Ω)  β 0

¯qh L 2 (Ω)  β˜ 0

( p¯ h , ∇ · u¯ h ) , |∇ u¯ h | u¯ h ∈ V h \{0} sup

¯ h) (¯qh , ∇ · w . ¯ |∇ w h| \{0}

sup ¯ h ∈ V˜ h w

1166

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

˜ h and W h : (H4) Compatibility condition between M

˜ h, ( W h · W h ) ∩ L 20 (Ω) ⊂ M that is,

¯h1 ,

¯h2

∀ρ ρ ∈ W h ,

¯h1 ¯h2

ρ ρ −

1

|Ω|

˜ h. ρ¯h1 (x)ρ¯h2 (x) dx ∈ M

Ω

˜ h ): (H5) Compatibility condition between ( M h , M

˜ h. Mh ⊂ M (H6) Stability properties

| J h u |  C |u | ∀u ∈ L 2 (Ω), |∇ J h u |  C |∇ u | ∀u ∈ H 10 (Ω), I h ρ H 1 (Ω)  C ρ H 1 (Ω) ∀ρ ∈ H 1 (Ω). ˜ h ) verifying (H2)–(H6) is the following. Let {Th }h>0 For instance, a way of defining the discrete spaces ( W h , V h , M h , V˜ h , M be a regular, quasi-uniform family of triangulations of Ω , with h = max K ∈Th h K (h K = diameter of K ), and





X hl = xh ∈ C 0 (Ω) such that xh | K ∈ Pl ( K ), ∀ K ∈ Th . Then we define W h = X h1 . There are several possibilities to define ( V h , M h ) [5], by using the Taylor–Hood element (P2 × P1 )

˜h= ˜ h ) we choose V˜ h = X 3 ∩ H 1 (Ω) and M or the mini-element (P1 + bubble × P1 ), for instance. For the spaces ( V˜ h , M 0 h X h2 ∩ L 20 (Ω).

˜ h are chosen, we do not need to consider the projection problem (13). Note that if V h = V˜ h and M h = M

3.2. Main results The aim of this work is to prove the existence of a weak solution to problem (1) in a fully discrete setting by using finite elements. The essential part of such a proof lies in obtaining energy estimates for the scheme independent of the discrete parameters from which one can infer their weak convergence. Then a compactness argument provides the strong convergence. The first works in that line for the simplified problem (1) were developed in [9] and [10]; unfortunately, arguments there do not carry over to the present case due to the troublesome term λ2 ∇ · ( ρ1 ∇ ρ ⊗ ∇ ρ ) in (1)1 . To construct finite-element approximations of weak solutions to (1), we must face two difficulties. The first thing is that the maximum principle (7) has to be satisfied by the finite-element approximation for the density. The second thing is how to establish the discrete version of the energy estimate (12). To overcome these difficulties, we propose the following scheme. For simplicity, we choose a uniform partition of [0, T ] N with time step k = T / N: (tn = nk)nn= =0 . We thus consider a backward Euler type scheme which is implicit with respect to the diffusion terms and semi-implicit with respect to the convective terms. The λ2 -term requires special attention since depending on how we integrate it in time we will be able or not to prove the stability of the scheme. Then the approximations of (1) remain as follows. Initialization: Let (uh0 , ρh0 ) ∈ V h × W h be approximations of (u 0 , ρ0 ) as h → 0. Time step n + 1: Given (unh , pnh , ρhn ) ∈ V h × M h × W h .

˜ h such that, for each ( w ˜ h, ¯ h , q¯ h ) ∈ V˜ h × M • Find ( w nh , qnh ) ∈ V˜ h × M

 

     ¯ h − qnh , ∇ · w ¯ h = ∇ unh , ∇ w ¯h , ∇ w nh , ∇ w  ∇ · w nh , q¯ h = 0.

(13)

• Find ρhn+1 ∈ W h such that, for each ρ¯h ∈ W h :



ρhn+1 − ρhn k

    ¯ , ρh + w nh · ∇ ρhn+1 , ρ¯h + λ ∇ ρhn+1 , ∇ ρ¯h = 0.

(14)

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1167

• Find (unh+1 , pnh+1 ) ∈ V h × M h such that, for each (u¯ h , p¯ h ) ∈ V h × M h :

⎧  n +1   un+1 − unh ⎪ 1 ρh − ρhn n+1 ⎪ n h ⎪ ρh , u¯ h + u h , u¯ h + c ρhn+1 unh − λ∇ ρhn+1 , unh+1 , u¯ h ⎪ ⎪ k 2 k ⎪ ⎨   n +1 n +1  1 n +1 2 n ¯ ¯ + a ρ , u , u ∇ ρ ⊗ ∇ ρ , ∇ u − λ ⎪ h h h h h ⎪ ⎪ ρhn+1 h ⎪ ⎪ ⎪     ⎩ = ρhn+1 f n+1 , u¯ h + pnh+1 , ∇ · u¯ h ,   ∇ · unh+1 , qh = 0,

(15)

(16)

where we have used the following short-hand notation:

f

n +1

=

1

t n +1

f (t ) dt ,

k tn

a(ρ , u , v ) = μ(∇ u , ∇ v ) +

λ

 ˜ ˜ M +m 2

− ρ (∇ u )t : ∇ v dx,

Ω

˜ > M and 0 < m ˜ < m, and with M c( w , u, v ) =

1  2

   ( w · ∇)u , v − ( w · ∇) v , u .

We list here some results related to coercivity and continuity properties for the trilinear form defined above which will be used later:

 a(ρ , u , u ) 

μ−λ

˜ −m ˜ M 2

a(ρ , u , v )  C u H 1 v H 1

|∇ u |2

˜, ˜ ρ M if m

(17)

if ρ L ∞ (Ω)  C ,

c ( w , u , u ) = 0,

(18)

c ( w , u , v )  C w L 3 u H 1 v H 1 .

(19)

Here and below, we denote by C , with or without subscript, different positive constants, always independent of the discrete parameters (k, h) and, eventually, depending on the diffusion parameter λ. In this last case, we will denote the constants by C λ . Scheme (13)–(16) has the following main features. At each time step, three linear systems which need to be solved separately to compute (ρhn+1 , unh+1 , pn+1 ). This would be started with w nh as an H 1 orthogonal projection of unh onto the

˜ h (this projection will guarantee an approxi˜ h ) which in turn depends on M discrete free-divergence space related to ( V˜ h , M mate discrete maximum principle for the density), second ρhn+1 as a finite-element approximation of a convection–diffusion equation with w nh being the convective velocity, and third (unh+1 , pnh+1 ) as a mixed finite-element approximation of the Navier–Stokes-like equations. In [10] scheme (13)–(16) was proposed without the term λ2 ( n1+1 ∇ ρhn+1 ⊗ ∇ ρhn , ∇ u¯ h ). There we follow the following ρh

strategy to establish stability and convergence. Firstly, starting from a slight variant of the truncated scheme developed in [9], we obtain weak estimates for the velocity; afterwards an approximate discrete maximum principle is established for the density; as a consequence of this, the truncation operator is unnecessary. Secondly, strong estimates for the density are attained based on the discrete interpolation (23) instead of the discrete version of the Gagliardo–Nirenberg interpolation for two dimensions used in [9]. Finally, the convergence is established by compactness arguments. The incorporation of this λ2 term introduces new difficulties, obligating to change the strategy of the main proofs. Now, the relation with a truncated scheme does not work due to the weak estimates for the velocity and the strong ones for the density are not obtained in an independent way as seen in the proof of Theorem 3. We now make an induction process with respect to the time step to obtain, the three main properties of the scheme: an approximate discrete maximum principle, and weak and strong estimates for velocity and density, respectively. Namely, fixed a time step, we obtain firstly pointwise estimates for the density, and then weak estimates for the velocity and strong estimates for the density at the same time (see Lemma 8 below). We will see that scheme (13)–(16) is conditionally stable and convergent by imposing the constraint (S) below (which was already used in [10]). We will also obtain, as in [10], the convergence towards a weak solution of the density-dependent Navier–Stokes equations, as the diffusion parameter goes to zero jointly with the space and time parameters according to constraint (S ) below. Since this argument is rather similar to that of [10], we only explicit here the pass to the limit of the λ2 -terms.

1168

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

By u h,k and ρh,k we denote the piecewise constant functions taking values u h,k = unh and ρh,k = ρhn on (tn−1 , tn ], respectively (which we will denote by u h,k,λ , ρh,k,λ when studied the asymptotic behavior with respect to λ). Thus the following two main results will be proved in this paper. Theorem 4. Under the hypotheses (H0)–(H6), and the constraint

lim

h

(h,k)→0 k

= 0,

(S)

there exists a convergent subsequence of (u h,k , ρh,k ) (denoted in the same way) as (h, k) → 0 towards a weak solution (u , ρ ) of problem (1), (3)–(4) in the sense of Definition 1. Theorem 5. Under the hypotheses of Theorem 4 and extending (H2) by the additional approximation hypothesis

(H2)

|ρ¯ − P h ρ¯ |  Ch2/3 ρ¯ W 1,3/2 (Ω) ∀ρ¯ ∈ W 1,3/2 (Ω) (here P h is the L 2 -projector on W h ) and changing constraint (S) by

1

lim

(λ,h,k)→0



λ

h k

(S )

= 0,

then there exists a convergent subsequence of (u h,k,λ , ρh,k,λ ), as (h, k, λ) → 0, towards a weak solution (u , ρ ) of the densitydependent Navier–Stokes problem in the sense of Definition 2. Finally, we show some numerical computations in order to prove the practical usefulness of scheme (13)–(16). In fact, we use scheme (13)–(16), with λ = 0, in order to compute some rates of convergence and the instability of Rayleigh–Taylor to the density-dependent Navier–Stokes problem (6), and a slight adaptation of scheme (13)–(16) in order to simulate a powder-snow avalanche. 4. A priori estimates Since (13)–(16) are a sequence of three square linear systems, to prove existence and uniqueness it suffices to prove only uniqueness, which will be a consequence of the stability of scheme (14)–(16) given in this section. 4.1. Approximate discrete maximum principle The proof of the following lemma can be found in [10]. Lemma 6. Fixed n: 0  n  N − 1, if the discrete velocities ( w lh )ln=0 satisfy k n+1 h

and n, then the discrete solution ρ

˜ ˜  ρhn+1  M 0
n

l 2 l=0 |∇ w h |

 C d , with C d > 0 independent of (h, k, λ),

of (14) satisfies the pointwise estimates

in Ω

provided k and h sufficiently small satisfying (S). 4.2. Weak estimates for the velocity and strong ones for the density Consider the linear operator h : W h → W h defined as:

−(h ρh , ρ¯h ) = (∇ ρh , ∇ ρ¯h ) ∀ρ¯h ∈ W h .

(20)

Then the discrete density equation (14) can be rewritten as:



ρhn+1 − ρhn k

    , ρ¯h + w nh · ∇ ρhn+1 , ρ¯h − λ h ρhn+1 , ρ¯h = 0.

(21)

Before we proceed any further, we need to establish the discrete version of the Gagliardo–Nirenberg interpolation and (8). Lemma 7. There exists C = C (Ω) > 0 such that, for any ρh ∈ W h , one has:

∇ ρh L 3 (Ω)  C |∇ ρh |1/2 |h ρh |1/2 ,   1/2 ∇ ρh L 4 (Ω)  C h1/4 |h ρh | + ρh L ∞ (Ω) |h ρh |1/2 .

(22) (23)

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1169

Proof. Estimate (22) is obtained in [10]. To prove (23), one firstly obtains

 1/2   ∇ ρh L 4 (Ω)  C h1/4 |h ρh | + ρ (h) L ∞ (Ω) |h ρh |1/2 , where

(24)

ρ (h) ∈ H 2 (Ω) is the solution of the following elliptic problem 

∂ ρ (h)  −ρ (h) = −h ρh in Ω, = 0 , ρ (h)(x) dx = 0. ∂ n ∂Ω

(25)

Ω

Indeed (24) is based on the inverse inequality ∇ ρh L 4 (Ω)  Ch−3/4 |∇ ρh |, the approximation property ∇ ρ − ∇ I h ρ L 4 (Ω)  Ch1/4 ρ H 2 (Ω) , and the interpolation inequality (8). Next, to obtain (23) from (24), it is necessary to reason as in [10] by comparing ρ (h) with ρh , and using the inverse inequality ρh L ∞ (Ω)  Ch−1/2 ρh H 1 and the interpolation error ρ − I h ρ L ∞ (Ω)  Ch1/2 ρ H 2 (Ω) . Now, we are in a position to prove some recursive inequalities for scheme (13)–(16). Lemma 8. Fixed n: 0  n  N − 1, assume

˜ ˜  ρhn , ρhn+1  M 0
in Ω

(26)

and

1 16C 2

 2 1  2 λ2 μkh ρhn  + μk∇ unh   C d ,

(27)

2

with C d > 0 independent of (h, k, λ), and n. Then, provided that h and k are sufficiently small satisfying (S), there exists a unique solution (ρhn+1 , unh+1 , pnh+1 ) of scheme (13)–(16) which satisfies:

⎧           ⎪ ⎨  ρ n+1 un+1 2 −  ρ n un 2 +  ρ n un+1 − un 2 + 3μ k∇ un+1 2 h h h h h h h h 4        ⎪ n+1 2 2  n 2 ⎩  kC 1  f n+1 26/5   , + εμ λ k ρ +  ρ h h h h L (Ω)

(28)

2 2 2  2     2 λ2  λ∇ ρhn+1  − λ∇ ρhn  + λ∇ ρhn+1 − ρhn  + kh ρhn+1   C 2k∇ unh  ,

(29)

2

where C 1 , C 2 , ε > 0 are constants independent of (h, k, λ), and n, with ε being arbitrarily small. Proof. Taking u¯ h = 2kunh+1 and p¯ h = pnh+1 in (15)–(16), and using the identity (a − b, 2a) = |a|2 − |b|2 + |a − b|2 and properties (17) and (18) gives:

⎧    2         ⎨  ρ n+1 un+1 2 −  ρ n un 2 +  ρ n un+1 − un 2 + μ − λ( M ˜ −m ˜ ) k∇ unh+1  h h h h h h h  2 ⎩  kC 1  f n+1  L 6/5 (Ω) + A ,

(30)

where A = 2k λm˜ ∇ ρhn+1 L 4 (Ω) ∇ ρhn L 4 (Ω) |∇ unh+1 | and C 1 = C 1 (μ) > 0 a constant independent of (h, k, λ), and n. In view of inequality (23) and hypothesis (26), the term A can be bounded by 2





















˜ 1/2 h1/4 h ρ n+1 h ρ n 1/2 + h ρ n h ρ n+1 1/2 A  Ckλ2 h1/2 h ρhn+1 h ρhn  + M h h h h



      ˜ h ρ n+1 1/2 h ρ n 1/2 ∇ un+1  +M h h h

:= A 1 + A 2 + A 3 + A 4 . 1/ 2

1/ 2

By hypothesis (27), |h ρhn |  4C 2 C d /(λμ1/2 k1/2 ), hence the term A 1 is bounded as follows: 1/2

A1 

C Cd

1/2

μ







h1/2k1/2 λh ρhn+1 ∇ unh+1   k

ε1 k

 2  2 λ2 μh ρhn+1  + ε1 μk∇ unh+1  ,

ε1 > 0 and C 0 = C 0 (μ). Thus, by stability constraint (S), we are allowed to take h and k small enough such that  ε2 , with ε2 > 0 being arbitrary, to arrive at

for any C0 h ε1 k

C0 h



2



2

A 1  ε2 λ2 μkh ρhn+1  + ε1 μk∇ unh+1  . The term A 2 has a similar treatment as A 1 ,

1170

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1/4

A 2  Ckλ2 h1/4

Cd

λ1/2

1/4 k1/4

μ

       h ρ n+1 ∇ un+1   ε2 λ2 μkh ρ n+1 2 + ε1 μk∇ un+1 2 , h h h h

˜ , Ω, μ). ˜,M provided that h and k satisfy εC λ( hk )1/2  ε2 , where C = C (m 1 Analogously, we estimate the term A 3 as 1/4

Cd

A 3  2Ckλ2 h1/4

λ1/2

1/4 k1/4

μ

      h ρ n 1/2 h ρ n+1 1/2 ∇ un+1  h h h

2  2   2   ε2 μλ k h ρhn  + h ρhn+1  + ε1 μk∇ unh+1  , 2

˜ , Ω, μ). ˜,M provided that h and k satisfy εC λ( hk )1/2  ε2 , where C = C (m 1 For A 4 , we bound as follows:







 

2 2 A 4  ε2 μλ2 k h ρhn+1  + h ρhn  + kC

λ2 

μ

2 ∇ unh+1  ,

˜ , Ω, ε3 ). ˜,M where C = C (m Using the previous estimates for the A i terms in (30), we get

  n+1 n+1 2  n n 2  n  n+1  2  λ ˜ λ2 n 2  ρ  − ρ u  + ρ u ˜ u − u + 1 − ( M − m ) − C − 4 ε μk∇ unh+1  1 h h h h h h h 2  2  kC 1  f n+1  6/5 L

μ

μ

2  2    + ε2 μλ k 4h ρhn+1  + 2h ρhn  . (Ω) 2

(31)

Accordingly, choosing ε1 , ε2 and λ/μ sufficiently small, one arrives at (28). To obtain (29), write scheme (14) as (21). Then take ρ¯h = λkh ρhn+1 as test functions in (21), ρ¯h = λ(ρhn+1 − ρhn ) in (20) for

ρh = ρhn+1 and use the identity (a − b, a) = 12 (|a|2 − |b|2 + |a − b|2 ) to obtain:

2  2    2    λ  n+1 2 − λ∇ ρhn  + λ∇ ρhn+1 − ρhn  + λ2kh ρhn+1  = λk w nh · ∇ ρhn+1 , h ρhn+1 := I . ∇ ρh 2

(32)

As −ρ n+1 (h) = −h ρhn+1 , where

ρ n+1 (h) is defined in (25) for ρh = ρhn+1 , we may write the term I of (32) as:       I = λk w nh · ∇ ρ n+1 (h), ρ n+1 (h) + λk w nh · ∇ ρhn+1 − ρ n+1 (h) , h ρhn+1 .

We estimate I as in [10]:

























˜ h ρ n+1   C λk∇ un h1/2 h ρ n+1 2 + C 2k∇ un 2 + I  C λk∇ unh  h1/2 ρhn+1  + M h h h h where we have used the stability property |∇ w nh |  |∇ unh | of scheme (13). Making use of the hypothesis

1 2

μk|∇ unh |2  C d and selecting h and k such that

C

λ



h k

1 4

 2 λ2kh ρhn+1  ,

 14 , we get (29). 2

To prove the global stability of scheme (13)–(16), we consider the approximations of the data (ρ0 , u 0 , f ) to verify the following properties:

˜, ˜  ρh0  M 0


2

μλ∇ ρh0   K 1 ,

C2 2 1  μk∇ uh0   K 3 , 2 C 1k

N −1  n =0

 n +1  2 f  6/5 L

 0 0 2  ρ u   K2, h

1 16C 2 (Ω)

h

 2 λ2 μkh ρh0   K 4 ,

 K5,

where K i > 0 are constants independent of (h, k, λ). These properties can be guaranteed by hypotheses (H0), (H2), and (H6), considering

ρh0 = I h ρ0 and uh0 = J h u 0 . Indeed, the pointwise estimates for ρ are deduced from 0 < m  ρ0  M < ∞ and the error estimate ρ0 − I h ρ0 L ∞ (Ω)  Ch1/2 ρ0 H 2 (Ω) . In view of the stability (H6), there exist K 1 > 0, K 2 > 0, and K 3 > 0 such that C1 μλ|∇ ρh0 |2  K 1 , 2 0 h

 | ρh0 uh0 |2  K 2 , and μk|∇ uh0 |2  K 3 . To prove λ2 μk|h ρh0 |2  K 4 , we take ρh = ρh0 and ρ¯h = −h ρh0 into (20), arriving

at

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1171

        h ρ 0 2  − ∇ ρ 0 , ∇h ρ 0 = ∇ ρ0 − ∇ ρ 0 , ∇h ρ 0 − ∇ ρ0 , ∇h ρ 0 h h h h h h     C  |∇ ρ0 − ∇ I h ρ0 |h ρh0  + |ρ0 |h ρh0  h    C ρ0 H 2 (Ω) h ρh0 , where we have used the inverse inequality ρ¯h H 1 (Ω)  Ch−1 |ρ¯h | and the interpolation error |∇(ρ¯ − I h ρ¯ )|  Ch ρ¯ H 2 (Ω) . Therefore, we may find a constant K 4 > 0 such that λ2 μk|h ρ0h |2  K 4 . Finally, the bound C 1 k 2

easily deduced from the regularity L (0, T , L

6/ 5

 N −1 n=0

f n+1 2L 6/5 (Ω)  K 5 is

(Ω)) for f .

Theorem 9. Assume (S) and (H0)–(H6). Then, for h and k small enough (depending on λ), and for each n = 0, . . . , N − 1, the discrete problem (14)–(16) has a unique solution (ρhn+1 , unh+1 , pnh+1 ) which satisfies the following estimates:

˜ ˜  ρhn+1  M 0
1 4C 2

in Ω

(33)

    n+1 2  n+1 n+1 2  n 2  n n 2  2 1         +  ρhn unh+1 − unh  μλ ∇ ρh + ρh uh − μλ ∇ ρh + ρh uh 4C 2

+

1 4C 2

 2   2 3μ  n+1 2 1 μλ∇ ρhn+1 − ρhn  + μλ2kh ρhn+1  k ∇ u h  +

 2  C 1k f n+1  6/5 L

4

8C 2

 2 1  2 + εμλ kh ρhn  + μk∇ unh  , (Ω) 2

(34)

4

where C 1 , C 2 , ε > 0 are constants independent of (h, k, λ), and n, with ε > 0 being sufficiently small (appearing in Lemma 8). Proof. We proceed by induction on n. Let us define C d = K 1 + K 2 + K 3 + K 4 + K 5 . With such a constant C d , the hypotheses ˜ holds. In virtue of Lemma 8 (since μλ2 k|h ρ 0 |2  C d and the ˜  ρh1  M of Lemma 6 are verified for n = 0, then 0 < m h

˜ we get (28) and (29) for n = 0. Next, multiplying ˜  ρh0 , ρh1  M), discrete maximum principle for ρh0 and ρh1 hold, that is, m μ (29) by 4C for n = 0, and adding it to (28) for n = 0, we arrive at (34) for n = 0. 2 = 0, . . . , n − 1. Adding (34) for l = 0, . . . , n − 1, we obtain We that the induction hypothesis is true for l n now assume μk l=1 |∇ ulh |2  C d ; hence, due to |∇ w nh |  |∇ unh |, one has μk ln=1 |∇ w lh |2  C d . Thus, we are in the situation of Lemma 6 ˜ holds. Accordingly, by Lemma 8, we get estimates (28) and (29). Balancing both estimates as in ˜  ρhn+1  M for n, then m case n = 0, we arrive at (34). 2 From Theorem 9, it is easy to obtain the following estimates:

Corollary 10. Under the assumptions of Theorem 9, the solution (ρhn+1 , unh+1 ) of the discrete problem (14)–(16) satisfies the following estimates:

  (i) max unh   C , 0n N

(ii) k

N   n 2 ∇ u   C , h

N −1 

 n +1 2 u − un   C ,

(iii)

h

n =0

h

n =0

2    2 (iv) max ρhn  L ∞ + λ∇ ρhn   C , 0n N

(v) kλ2

N −1 

  h ρ n+1 2  C , h

n =0

(vi) λ

N −1 

  n +1 2 ∇ ρ − ρn   C , h

h

n =0

where C > 0 is a constant independent of (k, h, λ). The proof of the following corollary is based on the properties of the discrete Stokes operator (13); for complete details, see [10]. Corollary 11. The following estimates hold:

  (vii) max  w nh   C , 0n N

(viii) k

N    ∇ w n 2  C h n =0

where C > 0 is independent of (k, h, λ). To study the convergence of scheme (13)–(16) towards a solution of (1), (3)–(4), let us define the following auxiliary functions:

1172

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

ˆ h,k , ρh,k , ρˆh,k , and ph,k as the piecewise constant functions in time, taking values Definition 12. One defines u h,k , uˆ h,k , w unh+1 , unh , w nh , ρhn+1 , ρhn , and pnh+1 on (tn , tn+1 ], respectively. In addition, we define ρ˜h,k ∈ C 0 ([0, T ]; V h ) as the continuous piecewise linear function in time such that ρ˜h,k (tn ) = ρhn . Then, by the estimates of Lemma 6, Corollaries 10 and 11, and inequality (22), the following result holds (see [10]): Lemma 13. Under the hypotheses of Theorem 9, the following estimates (independent of h and k) hold:

    ˆ h,k }h,k are bounded in L ∞ 0, T ; L 2 (Ω) ∩ L 2 0, T ; H 10 (Ω) , {uh,k }h,k , {uˆ h,k }h,k , { w   {ρh,k }h,k , {ρˆh,k }h,k , {ρ˜h,k }h,k are bounded in L ∞ 0, T ; H 1 (Ω) ∩ L ∞ (Ω) ,   {ρh,k }h,k is bounded in L 4 0, T ; W 1,3 (Ω) . ˆ h,k }h,k {ρh,k }h,k , {ρˆh,k }h,k , {ρ˜h,k }h,k (denoted in the same way), and Furthermore, there exist a subsequence of {u h,k }h,k , {uˆ h,k }h,k , { w limit functions u, ρ verifying the following weak convergence, as (h, k) → 0:



uh,k → u ,

uˆ h,k → u ,

ˆ h,k → u w

ρh,k → ρ ,

ρˆh,k → ρ ,

ρ˜h,k → ρ in

in









L 2 0, T ; H 10 (Ω) -weak, L

 ∞



0, T ; L 2 (Ω) -weak∗,

L ∞ ( Q )-weak∗,





L ∞ 0, T ; H 1 (Ω) -weak∗,

ρˆh,k → ρ in L 4 0, T ; W 1,3 (Ω) -weak. 5. Compactness To establish the convergence of scheme (13)–(16), we need the compactness for the discrete velocity and density. As in [10], one can deduce the following compactness for the discrete density,





ρh,k , ρˆh,k → ρ in L 2 0, T ; L 2 (Ω) -strong as (h, k) → 0,

(35)

based on the time derivative estimate

 N  n +1   ρh − ρhn 4/3    Cλ, k   k n =0

where C λ > 0 depends only on the data (ρ0 , u 0 , f , μ, λ), but is independent of h and k. Moreover, since the discrete density is bounded in L ∞ ( Q ), one also gets the strong convergence in L p ( Q ) for any p < ∞. On the other hand, comparing Eq. (20) for ρh = ρhn+1 and ρ¯h = ρhn+1 with its limiting equation (see [9]), one can obtain the convergence of the norm L 2 (0, T ; L 2 (Ω)) of ∇ ρh,k towards the same norm of ∇ ρ from (35). Consequently, one has

ρh,k − ρ L 2 (0,T ; H 1 (Ω)) → 0 as (h, k) → 0. Finally, the strong convergences for the velocity

ˆ h,k → u uh,k , uˆ h,k , w





in L p 0, T ; L 2 (Ω) -strong as (h, k) → 0,

∀p < ∞

can be deduced, as in [10], by means of the following result: Proposition 14. In the hypotheses of Theorem 9, the following estimate holds:

T −δ     ρh,k (t + δ) uh,k (t + δ) − uh,k (t ) 2 dt  C λ δ 1/4 ∀δ : 0 < δ < T ,

(36)

0

where C > 0 is independent of (h, k, δ), but depends on λ. Proof. Since ρh,k and u h,k are piecewise constant functions, it suffices to suppose that δ is a multiple of time step k, that is, δ = rk, for r = 0, . . . , N, and to demonstrate

k

N −r    m+r  m+r 2   C λ (rk)1/4 ∀r: 0  r  N .  ρ u h − um h h m =0

(37)

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1173

A proof of (37), without considering the λ2 -term, has been done in [10]. Then we only focus on the λ2 -term. ρ n+1 −ρ n

Adding to both sides of (15) the term 12 ( h k h , unh+1 · u¯ h ), multiplying by k, summing for n = m, . . . , m − 1 + r, taking +r u¯ h = um − um as test function, and using the equality h h









+r m+r +r m+r ρhm+r um − ρhm um um − um − ρhm um h = ρh h + ρh h, h h

we get (see [10]):

⎧      ⎪  ρ m+r (um+r − um )2 = − ρ m+r − ρ m , um · um+r − um ⎪ h h h h ⎪ h h h h ⎪ ⎪ ⎪ m −1+r ⎪ ⎪

 n+1 n+1 m+r   n +1 n  ⎪ +r ⎪ ⎪ −k a ρh , u h , u h − um uh − λ∇ ρhn+1 , unh+1 , um − um ⎪ h + c ρh h h ⎪ ⎪ n=m ⎨   m −1+r  n+1 n+1 m+r  1 ρhn+1 − ρhn n+1  m+r  ⎪ m m ⎪ ρh f , uh − uh + , uh · uh − uh ⎪ ⎪ +k 2 k ⎪ ⎪ n=m ⎪ ⎪ ⎪ m −1+r  ⎪ ⎪  m+r  1 ⎪ n +1 2 n m ⎪ ⎪ ∇ ρh ⊗ ∇ ρh , ∇ uh − uh . ⎩ −λ k n +1

(38)

ρh

n=m

Multiplying by k, summing (38) for m = 0, . . . , N − r, and bounding adequately on the right-hand side, we are going to get the desired estimate (37). Indeed,

  N −r m−1+r      m+r 1  2 2  n +1 n m I := λ k ∇ρ ⊗ ∇ ρh , ∇ uh − uh    ρ n +1 h n=m m =0



=

h

N −r m −1+r 2 2  

λ k ˜ m

  ∇ ρ n +1  h

m=0 n=m

N −1 λ2 k 2  ˜ m

n 

L 4 (Ω)

 n ∇ ρ  h

L 4 (Ω)

  m+r   (Fubini) ∇ u − um h h

 n +1       ∇ ρ  4 ∇ ρ n  4 ∇ um+r − um , h L (Ω) h h h L (Ω)

n=0 m=n−r +1

where

n=

⎧ ⎨0 ⎩

if n < 0,

n

if 0  n  N − r ,

N −r

if n > N − r .

Using the discrete interpolation inequality (23), we obtain

I  C λ2 k 2

N −1 

n 

         h ρ n+1 h ρ n  + h ρ n+1 h ρ n 1/2 ∇ um+r − um  h h h h h h

n=0 m=n−r +1

+ C λ2 k 2

N −1 

n 

          h ρ n h ρ n+1 1/2 + h ρ n+1 1/2 h ρ n 1/2 ∇ um+r − um  h h h h h h

n=0 m=n−r +1

:= I 1 + I 2 + I 3 + I 4 . We estimate I 1 as follows:

 N −1   λ2k   n+1  n h ρh h ρh k I1  C ˜ m n =0

 1/2

 C (rk)

2

λ k

N −1  n =0

  h ρ n+1 2 h

n 

  m+r 2  ∇ u − um h h

m=n−r +1

1/2 

2

λ k

N −1 

  h ρ n 2 h

1/2

1/2 

n 

1/2 k

m=n−r +1

 C (rk)1/2 .

n =0

The rest of the I i terms (i = 2, 3, 4) can be bounded in a similar way to conclude the proof.

2

1174

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

6. Passing to the limit in the momentum system Let C c∞ (Ω) be the space of infinitely differentiable functions with compact support in Ω . Consider v ∈ C 1 ([0, T ]; C c∞ (Ω)) with ∇ · v = 0 and v ( T ) = 0. We define v nh as the Stokes projection of v (t n ) onto V h , which provides the following result (a proof of this lemma can be seen in [9]): Lemma 15. Let u¯ ∈ C c∞ (Ω). Then there exists u¯ h ∈ V h such that

u¯ h → u¯

in H 10 (Ω) and

(∇ · u¯ h , qh ) = (∇ · u¯ , qh ) ∀qh ∈ M h .

Let v h,k ∈ L ∞ (0, T ; V h ) be the piecewise constant function taking values v nh+1 on (tn , tn+1 ] and let v˜ h,k ∈ C 0 ([0, T ]; V h ) be the continuous piecewise linear function such that v˜ h,k (tn ) = v nh . It holds that, as (h, k) → 0,

v h,k → v





in L ∞ 0, T ; H 10 (Ω)

and



v˜ h,k → v



in W 1,∞ 0, T ; H 10 (Ω) .

Taking u¯ h = v nh+1 as a test function in (15) and using Definition 12, we arrive at the variational formulation [10]:

⎧

T 

T ⎪ ⎪   ∂ ⎪ 0 ⎪ ρˆh,k uˆ h,k , v˜ h,k − ρ0h u 0h , v h + a(ρh,k , uh,k , v h,k ) − ⎪ ⎪ ⎪ ∂t ⎪ ⎪ ⎪ 0 0 ⎪ ⎪ ⎪ ⎪

T

T  ⎨ 1 2 ∇ ρh,k ⊗ ∇ ρˆh,k , ∇ v h,k + c (ρh,k uˆ h,k − λ∇ ρh,k , u h,k , v h,k ) + λ ⎪ ρh,k ⎪ ⎪ ⎪ 0 0 ⎪ ⎪ ⎪ ⎪

T

T  ⎪ ⎪ ∂ 1 ⎪ ⎪ ρ˜h,k , uh,k · v h,k . = (ρh,k f k , v h,k ) + ⎪ ⎪ ⎩ 2 ∂t 0

0

We just analyze the passage to the limit of the term λ2 for

T 0

( ρ1 ∇ ρh,k ⊗ ∇ ρˆh,k , ∇ v h,k ). As ρ0 ∈ H 2N (Ω) the same statements h,k

ρh,k can be assured to ρˆh,k . On one hand, taking ρh = ρhn+1 and ρ¯h = −h ρhn+1 into (20), we get         h ρ n+1 2  ∇ ρ n+1 ∇h ρ n+1   C ∇ ρ n+1 h ρ n+1 . h h h h h h

|h hn+1 |

ρ  n +1  ∇ ρ 

Therefore,

h



C |∇ hn+1 |. h

ρ In view of the discrete interpolation inequality (23), we infer  n+1 1/4  3/4  1/2  C ∇ ρh  h ρhn+1  + C ρhn+1  , L 4 (Ω)

which shows that {∇ ρh,k }h,k is bounded in L 8/3 (0, T ; L 4 (Ω)). On the other hand, ∇ ρh,k is compact in L 2 (0, T ; L 2 (Ω)) and bounded in L ∞ (0, T ; L 2 (Ω)), then ∇ ρh,k is compact in L p (0, T ; L 2 (Ω)) for any p < ∞. Therefore, the convergence ∇ ρh,k ⊗ ∇ ρˆh,k → ∇ ρ ⊗ ∇ ρ in L 4/3 (0, T ; L 2 (Ω))-weak holds as (h, k) → 0. Now, as 1/ρh,k → 1/ρ weakly- in L ∞ ( Q ) and strongly in L 2 ( Q ), then 1/ρh,k → 1/ρ strongly in L p ( Q ) for any p < ∞. Thus, we have 1/ρh,k ∇ ρh,k ⊗ ∇ ρˆh,k → 1/ρ ∇ ρ ⊗ ∇ ρ weakly in L 4/3 (0, T ; L 2 (Ω)). Finally, using that v h,k → v strongly in L ∞ (0, T ; H 10 (Ω)), one gets the desired convergence:

T  0

1

ρh,k



T  1 ∇ ρh,k ⊗ ∇ ρˆh,k , ∇ v h,k → ∇ρ ⊗ ∇ρ, ∇ v .

ρ

0

The convergence of the rest of terms follows as in [9] and [10]. Then the proof of Theorem 4 is finished. 7. Asymptotic behavior λ → 0 This section is devoted to the proof of Theorem 5. 7.1. Uniform estimates with respect to (h, k, λ) Following the derivation of Lemma 8, one observes that to obtain estimates independent of (h, k, λ) requires to assume  1/λ h/k to be sufficiently small (that is constraint (S )). This restriction appears in proving estimate (29) and the pointwise estimates of Lemma 6 (see [10]). That is the reason for what was said that k and h depend on λ in Theorem 9 for a fixed λ. Thus, from Section 4, considering the piecewise functions associated to the scheme denoted explicitly by the subscripts h, k and λ, one deduces the following estimates (independent of h, k, and λ):

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1175

    ˆ h,k,λ }h,k,λ are bounded in L ∞ 0, T ; L 2 (Ω) ∩ L 2 0, T ; H 10 (Ω) , {w

{uh,k,λ }h,k,λ ,

{uˆ h,k,λ }h,k,λ ,

{ρh,k,λ }h,k,λ ,

{ρˆh,k,λ }h,k,λ are bounded in L ∞ ( Q ),

  λ1/2 {ρˆh,k,λ }h,k,λ are bounded in L ∞ 0, T ; H 1 (Ω) ,   λ3/4 {ρh,k,λ }h,k,λ is bounded in L 4 0, T ; W 1,3 (Ω) ,    ∂ λ3/4 ρ˜h,k,λ is bounded in L 4/3 0, T ; L 2 (Ω) , ∂t h,k,λ

λ1/2 {ρh,k,λ }h,k,λ ,

whereas h, k and λ are small enough satisfying (S ). In addition,

˜, ˜  ρh,k,λ  M 0
√ ρ˜h,k,λ − ρh,k,λ L 2 (0,T ; L 2 (Ω))  ρh,k,λ − ρˆh,k,λ L 2 (0,T ; L 2 (Ω))  C k, √ uh,k,λ − uˆ h,k,λ L 2 (0,T ; H 1 (Ω))  C k. 0

Lemma 16. Assume (S ), (H0), (H1), (H2) + (H2) , (H3)–(H5). Then, provided that h, k and λ are small enough, there exist a sequence of {u h,k,λ }h,k , {uˆ h,k,λ }h,k, λ , {ρh,k,λ }h,k,λ , {ρˆh,k,λ }h,k (denoted in the same manner) and limit functions u and ρ verifying the following convergence, as (h, k, λ) → 0:







L 2 0, T ; H 10 (Ω) -weak,

uh,k,λ → u ,

uˆ h,k,λ → u ,

ˆ h,k,λ → u in H

ρh,k,λ → ρ ,

ρˆh,k,λ → ρ ,

ρ˜h,k,λ → ρ in L ∞ ( Q )-weak∗.





L ∞ 0, T ; L 2 (Ω) -weak∗,

7.2. Compactness The following estimate holds: T −δ

    ρh,k,λ (t + δ) uh,k,λ (t + δ) − uh,k,λ (t ) 2 dt  C δ 1/4 ,

∀δ : 0 < δ < T ,

0

with C > 0 independent of h, k, δ and λ. Indeed, a proof of this inequality has been done in [10] (by using in particular the additional hypothesis (H2) ), where the λ2 -term is not considered, but this λ2 -term has been estimated independent of λ in the proof of Proposition 14. ˜ one deduces the compactness result: In particular, since ρh,k,λ  m,

uh,k,λ → u

and

ˆ h,k,λ → u w





in L 2 0, T ; L 2 (Ω)

as (h, k, λ) → 0.

7.3. Passing to the limit We write the scheme in the following conservative form [10]:

⎧

T 

T ⎪ ⎪   0 0 ∂ ⎪ ⎪ ρˆh,k,λ uˆ h,k,λ , v˜ h,k − ρh uh , v˜ h,k (0) + a(ρh,k,λ , uh,k,λ , v h,k ) − ⎪ ⎪ ⎪ ∂t ⎪ ⎪ ⎪ 0 0 ⎪ ⎪ ⎪ T ⎪

T  ⎪ ⎪ 1 ⎪ 2 ⎪ ∇ρ ⊗ ∇ ρˆh,k,λ , ∇ v h,k − (ρh,k,λ uˆ h,k − λ∇ ρh,k,λ ⊗ uh,k,λ , ∇ v h,k ) + λ ⎪ ⎪ ⎪ ρh,k,λ h,k,λ ⎪ ⎪ ⎪ 0 0 ⎪ ⎪ ⎪ ⎪

T

T  ⎨ ∂ 1 ρ˜h,k , uh,k,λ · v h,k − Q h (uh,k,λ · v h,k ) = (ρh,k,λ f k , v h,k ) + ⎪ 2 ∂t ⎪ ⎪ ⎪ 0 0 ⎪ ⎪ ⎪ ⎪

T ⎪ ⎪   1  ⎪ ⎪ ρh,k,λ uˆ h,k , ∇ uh,k,λ · v h,k − Q h (uh,k,λ · v h,k ) − ⎪ ⎪ ⎪ 2 ⎪ ⎪ 0 ⎪ ⎪ ⎪ ⎪ ⎪

T ⎪ ⎪   λ ⎪ ⎪ − k ρh,k,λ , uh,k,λ · v h,k − Q h (uh,k,λ · v h,k ) , ⎪ ⎪ ⎩ 2 0

1176

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

⎧

T 

T ⎪ ⎪ d ⎪ ⎪ ˆ ˜ ρ , η dt + λ (∇ ρh,k,λ , ∇ ηh,k ) dt − ⎪ h,k,λ h,k ⎪ ⎪ dt ⎨ 0

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

0

T −

(39)

 0



( Hˆ h,k,λ ρh,k,λ , ∇ ηh,k ) dt = ρ0h , ηh ,

0

where v h,k , v˜ h,k (respectively, ηh,k , η˜ h,k ) are suitable approximations of a free-divergence test functions v ∈ C 1 ([0, T ]; C c∞ (Ω)) (respectively, C 1 ([0, T ]; C c∞ (Ω))) with v ( T ) = 0 (respectively, η( T ) = 0) such that the sequence { v h,k,λ }h,k,λ is bounded in L ∞ (0, T ; W 1,3 (Ω) ∩ L ∞ (Ω)). Here, we only pass to the limit in the λ2 -term, because the convergence of the rest of terms has been made in [10]. We must prove that

T  J := λ



1

2

ρh,k,λ

0

∇ ρh,k,λ ⊗ ∇ ρˆh,k,λ , ∇ v h,k → 0 as (h, k, λ) → 0.

Indeed, since in particular v h,k is bounded in L 2 (0, T ; W 1,3 (Ω)), one has

 J  Cλ

1/2

1/4 

T λ

3

4 h,k,λ L 3 (Ω)

∇ ρ

1/4  T

T λ

ˆh,k,λ 4L 3 (Ω)

3

0

∇ ρ 0

1/2 ∇ v h,k 2L 3 (Ω)

0

 C λ1/2 → 0. Therefore, the proof of Theorem 5 is concluded. Remark 17. Replacing the semi-implicit approximation of the λ2 -term



λ

2

1

ρhn+1

n +1 h

∇ρ



n ¯ h , ∇ uh

⊗ ∇ρ

by the fully explicit approximation



λ

2

1

ρhn



n ¯ h , ∇ uh

n h

∇ρ ⊗ ∇ρ

,

one may establish the same stability, compactness and convergence results obtained previously in this work. On the contrary, if we consider the fully implicit approximation



λ2

n +1 n +1 ¯ ∇ ρ ⊗ ∇ ρ , ∇ u , h h h n +1 1

ρh

then the term λ2 (

CΩ k

1

ρhn+1

∇ ρhn+1 ⊗ ∇ ρhn+1 , ∇ unh+1 ) is estimated using the discrete interpolation inequality (23) by

2      λ2  1/2  ˜ 1/2 h1/4 h ρ n+1 3/2 + M ˜ h ρ n+1  ∇ un+1 . h ρhn+1  + 2 M h h h h ˜ m

Then, it is not clear how to control the terms

CΩ k

2     λ2  1/2  ˜ 1/2 h1/4 h ρ n+1 3/2 ∇ un+1  h ρhn+1  + 2 M h h h ˜ m

in order to obtain the stability estimates. 8. Numerical results In this section we present three type of numerical results that show that the scheme presented before does provide good approximations for fluids with variable density. In fact, we are going to approximate both the mass diffusion problem and the limit case λ = 0 i.e. the density-dependent Navier–Stokes problem. The scheme was implemented by using the free ˜ h = P1 for the velocity software FreeFem++ [12]. In all simulations, the approximating spaces are V h = V˜ h = P2 and M h = M and pressure, respectively, hence the projection step (i.e. the first step of the scheme) is not necessary. We consider two cases for the density: either W h = P2 or W h = P1 . Note that, with this choice, hypothesis (H4) is not verified.

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1177

Fig. 1. Errors in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms of the density, velocity and pressure, by using P1 elements for the density. We consider a fixed mesh size h = 0.0234458.

The first numerical result is related to a convergence test that illustrates that the scheme for λ = 0 leads to optimal rates of convergence for smooth solutions, and numerical rates of convergence with respect to λ → 0+ are also presented. The second numerical result is concerned to the so-called Rayleigh–Taylor instability for the density-dependent Navier– Stokes problem, with Reynolds number 5000 and 20 000, where has been needed to introduce a truncation procedure for the density in order to prevent numerical oscillations and to preserve stability in calculations (taking into account that stability constraints imposed in the previous numerical analysis ( S  ) is not verified in our numerical computations, because we choose k of the same order than h attending to a computational cost and λ = 0). This truncation consists in redefining the density ρhn+1 calculated by the continuity equation (14) as ρhn+1 = χ (ρhn+1 ), where χ represents the characteristic function on the interval [ρmin , ρmax ] defined by the physics of the problem. The last numerical computations are related to a powder snow avalanche modeled by a mass diffusion fluids with density-dependent viscosity (see (43) below). In this case, where Reynolds number is taken of order 106 , we also truncate the density and moreover, an adaptive mesh procedure is implemented, based on the residual of the discrete density equation. 8.1. Rates of convergence Let Ω be the unit disk. We consider the density-dependent Navier–Stokes equations (6) in Q with homogeneous slip boundary conditions and the external force f and initial values such that (ρex , u ex , p ex ) defined by









ρex = 2 + x cos sin(t ) + y sin sin(t ) , 

(40)

 u ex = − y cos(t ), x cos(t ) ,

(41)

p ex = sin(x) sin( y ) sin(t ).

(42)

is an analytical exact solution (see [6]). We will use it as a test for studying the error behavior of the scheme for the problem with λ = 0, with respect to the mesh size h and time step k, and when λ → 0+ . 8.1.1. Study for λ = 0 Our aim is to verify the capability of the scheme to solve numerically the density-dependent Navier–Stokes equations. We show the error between the numerical solutions obtained by the scheme, for λ = 0, and the exact solution given in (40)–(42) in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms. Firstly, in Figs. 1 and 2, we measure the error in time for the velocity, pressure, and density (this latter is approximated by P1 and P2 finite elements). In both cases, we consider the mesh size h = 0.0234458 and take time steps k smaller and smaller. The error in time for all unknowns is of order 1. Quantitatively, the best approximation is for the velocity and the worst is for the density in the L ∞ ( L 2 ) norm and for the pressure in the L ∞ ( H 1 ) norm. We point out that the error in time in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms is comparable with the error measured in the L 2 ( L 2 ) and L 2 ( H 1 ) norms. Secondly, in Fig. 3, we fix the time step k = 10−5 and we consider the mesh sizes h1 = 0.4484541, h2 = 0.3068574, and h3 = 0.2410898. The velocity error in space is of order 2 in both norms, and the pressure error is of order 2 in the L ∞ ( L 2 ) norm and of order 1 in the L ∞ ( H 1 ) norm.

1178

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

Fig. 2. Errors in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms for the density, velocity, and pressure, by using P2 elements for the density. We consider a fixed mesh size of h = 0.0234458.

Fig. 3. Errors in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms for the density, velocity and pressure, by using P2 elements for the density with a fixed time step k = 10−5 . The considered mesh sizes h1 = 0.4484541, h2 = 0.3068574, and h3 = 0.2410898.

8.1.2. Study for λ → 0+ Now we investigate the convergence behavior of the numerical solutions to the complete Kazhikhov–Smagulov system when λ → 0+ . In this case, we show the error between the numerical solution given by the scheme, for λ > 0, and the exact solution (40)–(42) in the L ∞ ( L 2 ) and L ∞ ( H 1 ) norms. These results are presented in Figs. 4 and 5. The error with respect to λ is of order 1 for the velocity and pressure in the L ∞ ( L 2 ) norm as shown in Fig. 4. As can also be seen in Fig. 4, this error is quantitatively best for the velocity followed by the pressure, and worst for the density. In Fig. 5, we depict the error with respect to λ in the L ∞ ( H 1 ) norm being of order 1/2 for the density and of order 1 for the velocity and pressure. 8.2. The Rayleigh–Taylor instability As a complement to the study for λ = 0 of our scheme, we use it to solve the problem of the viscous Rayleigh–Taylor instability. Numerical results for this problem has been reported by several authors (see [6] and [3], and the references therein). Since the solution has symmetries, we show the results for Ω = (0, d/2) × (−2d, 2d), the half of the original domain given by (−d/2, d/2) × (−2d, 2d). The fluid, subject to gravity, is initially at rest and its initial density ρ0 is given by

ρ0 =

ρM + ρm 2

+

ρM − ρm 2

 tanh

y + η cos(2π x/d) 0.01d

,

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1179

Fig. 4. Errors in the L ∞ ( L 2 ) norm for the density, velocity and pressure.

with

ρM > ρm > 0 and η > 0. The problem depends on the Atwood number At defined as At =

ρM − ρm , ρM + ρm

and the Reynolds number Re defined as

Re =

ρm d3/2 G 1/2 , μ

where G is the acceleration of the gravity. We have considered the following boundary conditionsfor the velocity: no slip on horizontal boundaries and slip on vertical boundaries. As in [3], the time is scaled as t = t  d/ AtG. We select d = 1, η = 0.1, ρm = 1 and ρ M = 3. Thus, At = 0.5. We compute solutions for mesh size h = 0.01, time step k = 0.01, Re = 5000 and Re = 20 000 by using P1 and P2 elements to approximate the density. The results are presented in Figs. 6, 7, 8, and 9. For P1 elements approximating the density, our results are in good qualitative agreement with those reported in [3] and [6]. On the other hand, if we use P2 for the density, our results differ for t  2 in the upper right part of the interface where the density changes. In our case, we observe extra rolls up forming. In [3] it is pointed out that for preventing perturbations in the velocity, it is recommended to use a triangular mesh with alternate directions. But, as we see in plots of Figs. 6, where we use P1 for the density, this is not our case. However, in those calculations where we use a P2 for the density (see Fig. 7), we can see it. We are not sure what is the origin of these extra rolls up; we conject that the truncation procedure and the mesh adaptation (in the case of the avalanches) for P1 elements introduce more numerical dissipation than for P2 .

1180

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

Fig. 5. Errors in the L ∞ ( H 1 ) norm for the density, velocity and pressure.

8.3. Powder snow avalanche We use our scheme to simulate powder snow avalanches. Some numerical simulations can be seen for example in [4] and the references therein. In particular, we adopt the test problem introduced in [4]. In this case, the domain Ω is the √ rectangle (− 3, 7) × (0, 3/2) minus an obstacle O = (5, 5 + 0.05h0 ) × (0, 0.25h0 ) (h0 > 0). Initially we consider a two-fluid at rest with density given by

ρ = ρ + 1ω+ + ρ − 1ω− , √

where ω+ = (− 3, 0) × (0, 1/2), ω− = Ω \ ω+ , 1ω is the characteristic function of the set ω , and ρ + > ρ − > 0 are two different constants (the density of the heavy √ and light fluid, respectively). The boundary conditions for √ the velocity are: no slip condition on the lateral surfaces Γl = {− 3, 7} × (0, 3/2) and on the top of the domain Γt = (− 3, 7) × {3/2} and slip on the bottom ∂Ω \ (Γl ∪ Γt ). The viscosity μ depends on the density. In fact, by following [4], if we have a two-fluid flow, we can consider the volume fraction of each fluid, denoted φi , i = 1, 2. These variables are defined in the following way. Let dΩ be an elementary fluid volume surrounding an interior point x ∈ dΩ , filled with two fluids of volume dΩ1 and dΩ2 , such that

|dΩ| = |dΩ1 | + |dΩ2 |. Then, the volume fraction of each fluid in the point x is defined as

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1181

Fig. 6. Time caption for the Rayleigh–Taylor instability for At = 0.5 (density ratio 3) and Re = 5000. The density is approximated by using P1 elements.

Fig. 7. Time caption for the Rayleigh–Taylor instability for At = 0.5 (density ratio 3) and Re = 5000. The density is approximated by using P2 elements.

φ i ( x) =

lim

|dΩ|→0 x∈dΩ

|dΩi | . |dΩ|

We retain only the heavy fluid volume fraction, denoted as φ (and 1 − φ will be the light fluid volume fraction). By using φ , we have

1182

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

Fig. 8. Time caption for the Rayleigh–Taylor instability for At = 0.5 (density ratio 3) and Re = 20 000. The density is approximated by using P1 elements.

Fig. 9. Time caption for the Rayleigh–Taylor instability for At = 0.5 (density ratio 3) and Re = 20 000. The density is approximated by using P2 elements.

ρ = φ ρ + + (1 − φ)ρ − , where

μ = φ ρ + ν + + (1 − φ)ρ − ν − ,

ν + , ν − are the corresponding kinematic viscosities of each fluid. In this case, the model, in non-conservative form, is

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1183

Fig. 10. Time caption for the density in powder snow avalanche for t = 0.2, 0.4, 0.8, 1.0, 1.2, 1.4, 1.6, 1.7 and 1.8 sec, for P1 elements (left) and for P2 elements (right).

   ⎧ t ⎨ ρ (ut + u · ∇ u − λ∇ log(ρ ) · ∇ u ) − ∇ · μ ∇ u + (∇ u ) + ∇ P = ρ f ∇ · u = 0 in Q , ⎩ ρt + ∇ · (ρ u − λ∇ ρ ) = 0 in Q ,

in Q , (43)

where λ > 0 is a mass diffusion coefficient, P is the pressure, u the velocity, and f = g (sin θ, cos θ) is the external force (depending on the gravity g and the angle θ of the slope of the domain with respect to the horizontal direction). Note that the λ2 -term in (15) is omitted to solve the corresponding momentum equations (43) with density-dependent viscosity. We consider the governing equations in their dimensionless form as in [4]. This procedure introduces two scaling parameters: the Reynolds number Re and the Schmidt number Sc defined as

1184

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

Fig. 11. Time caption for the density in powder snow avalanche for t = 2.0, 2.2 and 2.4 sec, for P1 elements (left) and for P2 elements (right). Table 1 Values of the parameters used for numerical simulations (see [4]).

Re :=

h0



gh0

ν+

,

Parameter

Value

Gravity acceleration, g Slope, θ Reynolds number, Re Schmidt number, Sc Initial height, h0 Obstacle height, h s Obstacle thickness

9.8 m/s2 32 ◦ C 106 0.3 1m 0.25h0 m 0.05h0 m

Parameters for light and heavy fluids

Heavy fluid

Light fluid

Density, ρ ± Kinematic viscosity,

4 kg/m3 10−5 m2 /s

1 kg/m3 10−5 m2 /s

Sc :=

ν+ λ

ν±

.

For this system we consider a little variant of the scheme presented in Section 3.2 to follow the complex dynamics of the phenomenon. The main ideas of this modified scheme are:

ρhn+1 is the same, 2. As we saw before, we truncate the approximate density as ρhn+1 = χ (ρhn+1 ), where χ is the characteristic function on the interval [ρ − , ρ + ] (ρ − = 1 and ρ + = 4, see Table 1), and we use an adaptive mesh procedure based on the residual 1. The equation for calculating the approximate density

for the density equation (see [12] for details). We put in Table 1 the values considered for the simulations. The time step is fixed to k = 0.005. We present the results in Figs. 10 and 11 for problem (43) by using P1 and P2 approximation for the density, respectively. The results for P1 are qualitatively similar to those reported in [4]. These differences could be again caused by different numerical dissipation of both approximations with respect to the density truncation and the adaptive procedure. References [1] S.N. Antontsev, A.V. Kazhikhov, V.N. Monakhov, Boundary Value Problems in Mechanics of Nonhomogeneous Fluids, Studies in Mathematical and Its Applications, vol. 22, North-Holland Publishing Co., Amsterdam, 1990. [2] H. Beirão da Veiga, Diffusion on viscous fluids, existence and asymptotic properties of solutions, Ann. Sc. Norm. Sup. Pisa 10 (1983) 341–355. [3] C. Calgaro, E. Crausé, Th. Goudon, An hybrid finite volume-finite element method for variable density incompressible flow, J. Comput. Phys. 227 (2008) 4671–4696. [4] D. Dutykh, C. Acary-Robert, D. Bresch, Numerical simulation of powder-snow avalanche interaction with an obstacle, Applied Mathematical Modelling, electronically available at http://hal.archives-ouvertes.fr/docs/00/35/88/80/PDF/AvalSim_Dutykh-AR-Bresch.pdf. [5] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, 1986. [6] J.-L. Guermond, L. Quartapelle, A projection FEM for variable density incompressible flows, J. Comput. Phys. 165 (1) (2000) 167–188. [7] F. Guillén-González, Sobre un modelo asintótico de difusión de masa para fluidos incompresibles, viscoso y no homogéneos, in: Proceedings of the Third Catalan Days on Applied Mathematics, ISBN 84-87029-87-6, 1996, pp. 103–114. [8] F. Guillén-González, P. Damázio, M.A. Rojas-Medar, Approach of regular solutions for incompressible fluids with mass diffusion by an interative method, J. Math. Anal. Appl. 326 (1) (2007) 468–487. [9] F. Guillén-González, J.V. Gutiérrez-Santacreu, Unconditional stability and convergence of a fully discrete scheme for 2D viscous fluids models with mass diffusion, Math. Comp. 77 (263) (2008) 1495–1524.

R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185

1185

[10] F. Guillén-González, J.V. Gutiérrez-Santacreu, Conditional stability and convergence of a fully discrete scheme for 3D Navier–Stokes equations with mass diffusion, SIAM J. Num. Anal. 46 (5) (2008) 2276–2308. [11] F. Guillén-González, J.V. Gutiérrez-Santacreu, Error estimates of a linear decoupled Euler–FEM scheme for a mass diffusion model, Numer. Math. 117 (2011) 333–371. [12] F. Hecht, FreeFem++ manual, available at http://www.freefem.org/ff++/. [13] A. Kazhikhov, Sh. Smagulov, The correctness of boundary value problems in a diffusion model of an inhomogeneous fluid, Sov. Phys. Dokl. 22 (1) (1977) 249–252. [14] P.-L. Lions, Mathematical Topics in Fluid Dynamics, vol. 2, Incompressible Models, Oxford University Press, United Kingdom, 1996. [15] P. Secchi, On the motion of viscous fluids in the presence of diffusion, SIAM J. Math. Anal. 19 (1988) 22–31. [16] J. Simon, Compact sets in the space L p (0, T ; B ), Ann. Mat. Pura Appl. 146 (1987) 65–97.