Stability and convergence for a complete model of mass diffusion

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 aDpto. de Ciencias Básicas, Universidad del Bío-Bío, Facultad de Ciencias, Campus Fernando May, Casilla 447, Chillán, Chile bDpto. E.D.A.N., University of Sevilla, Aptdo. 1160, 41080 Sevilla, Spain cDpto. de Matemática Aplicada I, University of Sevilla, E. T. S. I. Informática, Avda. Reina Mercedes, s/n, 41012 Sevilla, Spain article info abstract 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 We propose a fully discrete scheme for approximating a three-dimensional, strongly nonlinear model of mass diffusion, also called the complete Kazhikhov–Smagulov model. TheschemeusesaC0finite-element approximation for all unknowns (density, velocity and pressure), even though the density limit, solution of the continuous problem, belongs to H2. 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. 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 Ω⊆R3a 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: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (ρu)t+∇·(ρu−λ∇ρ)⊗u−λu⊗∇ρ−μu+λ2∇·1 ρ∇ρ⊗∇ρ+∇P=ρfin Q, ∇·u=0in Q, ρt+∇·(ρu−λ∇ρ)=0inQ, (1) *Corresponding author. E-mail addresses: [email protected] (R.C. Cabrales), [email protected] (F. Guillén-González), [email protected] (J.V. Gutiérrez-Santacreu). 1This work was made while the first author was a Postdoctoral student at Dpto. E.D.A.N. 2This author’s work was partially supported by project MTM2009-12927, Spain. 1162 R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 where ρ:Q→R+is the fluid density, u:Q→R3is the incompressible (averaged) velocity field, and P:Q→Ris the fluid pressure. Moreover, frepresents 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)n i=1,b=(bi)n i=1is denoted by a⊗bwith coefficients (a⊗b)i,j=aibj. 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: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ρut+(ρu−λ∇ρ)·∇u−∇·(μ−λρ)(∇u)t+μ∇u+λ2∇·1 ρ∇ρ⊗∇ρ+∇p=ρfin Q, ∇·u=0in Q, ρt+u·∇ρ−λρ=0in Q, (2) where p=P−λu·∇ρis a new potential function. We complete the model with the boundary conditions u(x,t)=0,∂ρ ∂n(x,t)=0for(x,t)∈Σ, (3) and the initial conditions ρ(0)=ρ0(x), u(0,x)=u0(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)⩽Mfor all x∈Ω. (5) The density-dependent Navier–Stokes problem is formally obtained from system (1) by just choosing λ=0: ⎧ ⎨ ⎩ ρut+(u·∇)u−μu+∇p=ρfin Q, ∇·u=0inQ, ρt+u·∇ρ=0inQ, (6) 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 λ2term, 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 λ2terms, is studied. This scheme consists of C0finite-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 Lp(Ω) and Hs(Ω),1⩽p⩽∞,s=0,1,...,we denote the classical Lebesgue and Sobolev spaces, respectively. The usual norm in Lp(Ω) and Hs(Ω) is denoted by ·Lp(Ω) and ·Hs(Ω). The norm and inner product in L2(Ω) will be denoted by |·|and (·,·), respectively. To define the concept of weak solution, we introduce the following space of functions: H=u∈L2(Ω):∇·u=0inΩ,u·n=0onΓ, V=u∈H1 0(Ω):∇·u=0inΩ, L2 0(Ω) =p∈L2(Ω): Ω p(x)dx=0, H2 N(Ω) =ρ∈H2(Ω):∂ρ ∂n=0on∂Ω,  Ω ρ(x)dx= Ω ρ0(x)dx. It is known by Poincaré’s inequality that uH1(Ω) and |∇u|are equivalent norms on H1 0(Ω). On the other hand, H2 N(Ω) is an affine space, and ∇ρH1(Ω) and |ρ|are equivalent semi-norms on H2 N(Ω). Definition 1. Apair(ρ,u)is said to be a weak solution of problem (1)–(3)–(4) on (0,T)if: (a) u∈L∞(0,T;H)∩L2(0,T;V),ρ∈L∞(0,T;H1(Ω)) ∩L2(0,T;H2 N(Ω)),0<m⩽ρ(x,t)⩽M,∀(x,t)∈Q. ((b) ∀φ∈C1([0,T];V)such that φ(T)=0, T  0−u,ρφt+(ρu−λ∇ρ)·∇φ+μ(∇u,∇φ) −λρ(∇u)t,∇φdt −λ2 T  01 ρ∇ρ⊗∇ρ,∇φdt = T  0 (ρf,φ)dt +ρ0u0,φ(0). (c) The equation of mass diffusion (1)cis satisfied almost everywhere in Q. Definition 2. Apair(ρ,u)is said to be a weak solution of problem (6)on (0,T)if it verifies: (a) u∈L∞(0,T;L2(Ω)) ∩L2(0,T;V),ρ∈L∞(Q)with 0 <m⩽ρ(x,t)⩽Ma.e. (x,t)∈Q. (b) For all φ∈C1([0,T];V)with φ(T)=0, T  0−ρu,φ t+(u·∇)φ+μ(∇u,∇φ)dt = T  0 (ρf,φ)dt +ρ0u0,φ(0). (c) For all ϕ∈C1([0,T];H1(Ω)) with ϕ(T)=0, 1164 R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 − T  0 (ρ,ϕt)dt − T  0 (ρu,∇ϕ)dt =ρ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 u0∈H,ρ0∈H1(Ω) and f∈L2(0,T;L6/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)⩽Ma.e. (x,t)∈Q.(7) Multiplying (1)cby −λρ, integrating over Ω, integrating by parts in the convective term, and taking into account the interpolation inequality ∇ρL4(Ω) ⩽CΩρ1/2 L∞(Ω)|ρ|1/2,(8) one arrives at λd dt |∇ρ|2+λ2|ρ|2⩽C0|∇u|2.(9) Adding the momentum system (1)aby uto the density equation (1)cby 1 2u·u,integratedoverΩ, the following energy equality holds: 1 2 d dt  Ω ρ|u|2dx+μ|∇u|2=λ Ω ρ(∇u)t:∇udx+λ21 ρ∇ρ⊗∇ρ,∇u+(ρf,u). (10) Reasoning as in the simplified model (see [13]), one has λ Ω ρ(∇u)t:∇udx=λ Ωρ−M+m 2(∇u)t:∇udx⩽λ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⩽ρ⩽Mand Young’s inequality, getting λ21 ρ∇ρ⊗∇ρ,∇u ⩽ε1μλ2|ρ|2+C2 ε1 λ2 μ|∇u|2, with ε1>0 to be chosen later on. The last term of (10) is easily bounded by (ρf,u)⩽μ 2|∇u|2+Cf2 L6/5(Ω). Compiling the above estimates into (10), we arrive at d dt  Ω ρ|u|2dx+2μ|∇u|2⩽1+C1 λ μ+C2 ε1 λ2 μ2μ|∇u|2+ε1μλ2|ρ|2+Cf2 L6/5(Ω).(11) Adding up (11) to (9) multiplied by με2with ε2>0 to be chosen later on, this gives us ε2μλd dt |∇ρ|2+d dt |√ρu|2+(ε2−ε1)μλ2|ρ|2+1−C1 λ μ−C2 ε1 λ2 μ2−ε2C0μ|∇u|2⩽Cf2 L6/5(Ω). Selecting ε1=ε2/2, and ε2and λ/μsmall enough such that C1λ μ+C2 ε1(λ μ)2+ε2C0⩽1 2,weget ε2μλd dt |∇ρ|2+d dt |√ρu|2+ε2 2λ2μ|ρ|2+μ 2|∇u|2⩽Cf2 L6/5(Ω). Integrating for t∈(0,T)and bounding from below the density, we obtain R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 1165 ε2μλ∇ρ(t) 2+mu(t) 2+ t  0ε2 2μλ2ρ(s) 2+μ∇u(s) 2ds ⩽C t  0 f(s)  2 L6/5(Ω) ds,(12) On the other hand, the following “fractional in time estimate” holds [1]: T−δ  0u(t+δ) −u(t) 2dt ⩽Cδ1/2∀δ∈(0,T). This estimate implies compactness for the velocity uin L2(0,T;L2(Ω)) [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 R3with a polyhedral boundary and that there exists a family of triangulations {Th}h>0of Ωmade up of tetrahedra or hexahedra in three dimensions, so that Ω=K∈ThK.Let Wh⊂H1(Ω),Vh,˜ Vh⊂H1 0(Ω) and Mh,˜ Mh⊂L2 0(Ω) be finite-elements spaces associated to density, velocity and pressure respectively. Throughout this work we will suppose the following hypotheses: (H0) Regularity for the data: u0∈V,ρ0∈H2 N(Ω) with 0 <m⩽ρ0⩽Min Ωand f∈L2(0,T;L6/5(Ω)). Assume λ/μsufficiently small. (H1) Assume Ωan open, bounded set of R3, whose boundary is polyhedral and such that the continuous dependencies in H2-norm of the Poisson–Neumann problem and in H2×H1-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∈Wh, |∇ ¯ ρh|L3(Ω) ⩽Ch−1/2|∇ ¯ ρh|∀ ¯ ρh∈Wh, ¯ ρhL∞(Ω) ⩽Ch−1/2¯ ρhH1(Ω) ∀¯ ρh∈Wh, ∇ ¯ ρhL4(Ω) ⩽Ch−3/4|∇ ¯ ρh|∀ ¯ ρh∈Wh, •and the interpolation errors: ¯ u−˜Jh¯ uH1(Ω) +¯ u−Jh¯ uH1(Ω) ⩽Ch|¯ u|H2(Ω) ∀¯ u∈H2(Ω) ∩H1 0(Ω), |¯ p−˜ Kh¯ p|+|¯ p−Kh¯ p|⩽Ch|¯ p|H1(Ω) ∀¯ p∈H1(Ω) ∩L2 0(Ω), |¯ ρ−Ih¯ ρ|+h¯ ρ−Ih¯ ρH1(Ω) ⩽Ch2|¯ ρ|H2(Ω) ∀¯ ρ∈H2(Ω), ¯ ρ−Ih¯ ρW1,3(Ω)∩L∞(Ω) ⩽Ch1/2|¯ ρ|H2(Ω) ∀¯ ρ∈H2(Ω), ¯ ρ−Ih¯ ρW1,4(Ω) ⩽Ch1/4|¯ ρ|H2(Ω) ∀¯ ρ∈H2(Ω), where Jh,˜Jh,Kh,˜ Khand Ihare interpolation operators from H2(Ω) ∩H1 0(Ω) into Vh,H2(Ω) ∩H1 0(Ω) into ˜ Vh, H1(Ω) ∩L2 0(Ω) into Mh,H1(Ω) ∩L2 0(Ω) into ˜ Mh, and H2(Ω) into Wh, respectively. (H3) Inf–sup conditions. There exist β>0 and ˜ β>0 (independent of h) such that, ∀¯ ph∈Mhand ∀¯ qh∈˜ Mh, ¯ phL2 0(Ω) ⩽βsup ¯ uh∈Vh\{0} (¯ ph,∇·¯ uh) |∇¯ uh|, ¯ qhL2 0(Ω) ⩽˜ βsup ¯ wh∈˜ Vh\{0} (¯ qh,∇·¯ wh) |∇ ¯ wh|. 1166 R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 (H4) Compatibility condition between ˜ Mhand Wh: (Wh·Wh)∩L2 0(Ω) ⊂˜ Mh, that is, ∀¯ ρ1 h,¯ ρ2 h∈Wh,¯ ρ1 h¯ ρ2 h−1 |Ω| Ω ¯ ρ1 h(x)¯ ρ2 h(x)dx∈˜ Mh. (H5) Compatibility condition between (Mh,˜ Mh): Mh⊂˜ Mh. (H6) Stability properties |Jhu|⩽C|u|∀u∈L2(Ω), |∇ Jhu|⩽C|∇u|∀u∈H1 0(Ω), IhρH1(Ω) ⩽CρH1(Ω) ∀ρ∈H1(Ω). For instance, a way of defining the discrete spaces (Wh,Vh,Mh,˜ Vh,˜ Mh)verifying (H2)–(H6) is the following. Let {Th}h>0 be a regular, quasi-uniform family of triangulations of Ω,withh=maxK∈ThhK(hK=diameter of K), and Xl h=xh∈C0(Ω) such that xh|K∈Pl(K), ∀K∈Th. Then we define Wh=X1 h. There are several possibilities to define (Vh,Mh)[5], by using the Taylor–Hood element (P2×P1) or the mini-element (P1+bubble ×P1), for instance. For the spaces (˜ Vh,˜ Mh)we choose ˜ Vh=X3 h∩H1 0(Ω) and ˜ Mh= X2 h∩L2 0(Ω). Note that if Vh=˜ Vhand Mh=˜ Mhare chosen, we do not need to consider the projection problem (13). 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] with time step k=T/N:(tn=nk)n=N n=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 (u0 h,ρ0 h)∈Vh×Whbe approximations of (u0,ρ0)as h→0. Time step n+1:Given (un h,pn h,ρn h)∈Vh×Mh×Wh. •Find (wn h,qn h)∈˜ Vhט Mhsuch that, for each (¯ wh,¯ qh)∈˜ Vhט Mh, ∇wn h,∇¯ wh−qn h,∇·¯ wh=∇un h,∇¯ wh, ∇·wn h,¯ qh=0.(13) •Find ρn+1 h∈Whsuch that, for each ¯ ρh∈Wh: ρn+1 h−ρn h k,¯ ρh+wn h·∇ρn+1 h,¯ ρh+λ∇ρn+1 h,∇¯ ρh=0.(14) R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 1167 •Find (un+1 h,pn+1 h)∈Vh×Mhsuch that, for each (¯ uh,¯ ph)∈Vh×Mh: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ρn h un+1 h−un h k,¯ uh+1 2ρn+1 h−ρn h kun+1 h,¯ uh+cρn+1 hun h−λ∇ρn+1 h,un+1 h,¯ uh +aρn+1 h,un+1 h,¯ uh−λ21 ρn+1 h∇ρn+1 h⊗∇ρn h,∇¯ uh =ρn+1 hfn+1,¯ uh+pn+1 h,∇·¯ uh, (15) ∇·un+1 h,qh=0,(16) where we have used the following short-hand notation: fn+1=1 k tn+1  tn f(t)dt, a(ρ,u,v)=μ(∇u,∇v)+ Ω λ˜ M+˜ m 2−ρ(∇u)t:∇vdx, with ˜ M>Mand 0 <˜ m<m, and 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|∇u|2if ˜ m⩽ρ⩽˜ M,(17) a(ρ,u,v)⩽CuH1vH1if ρL∞(Ω) ⩽C, c(w,u,u)=0,(18) c(w,u,v)⩽CwL3uH1vH1.(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 (ρn+1 h,un+1 h,pn+1). This would be started with wn has an H1orthogonal projection of un honto the discrete free-divergence space related to (˜ Vh,˜ Mh)which in turn depends on ˜ Mh(this projection will guarantee an approximate discrete maximum principle for the density), second ρn+1 has a finite-element approximation of a convection–diffusion equation with wn hbeing the convective velocity, and third (un+1 h,pn+1 h)as a mixed finite-element approximation of the Navier–Stokes-like equations. In [10] scheme (13)–(16) was proposed without the term λ2(1 ρn+1 h∇ρn+1 h⊗∇ρn h,∇¯ uh). There we follow the following 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 λ2term 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 uh,kand ρh,kwe denote the piecewise constant functions taking values uh,k=un hand ρh,k=ρn hon (tn−1,tn],respectively (which we will denote by uh,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,k)→0 h k=0,(S) there exists a convergent subsequence of (uh,k,ρh,k)(denoted in the same way)as (h,k)→0towards a weak solution (u,ρ)of problem (1),(3)–(4) in the sense of Definition 1. Theorem 5. Under the hypotheses of Theorem 4and extending (H2) by the additional approximation hypothesis |¯ ρ−Ph¯ ρ|⩽Ch2/3¯ ρW1,3/2(Ω) ∀¯ ρ∈W1,3/2(Ω) (H2) (here Phis the L2-projector on Wh)and changing constraint (S) by lim (λ,h,k)→0 1 λh k=0,(S) then there exists a convergent subsequence of (uh,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 (wl h)n l=0satisfy k n l=0|∇wl h|2⩽Cd,withC d>0independent of (h,k,λ), and n, then the discrete solution ρn+1 hof (14) satisfies the pointwise estimates 0<˜ m⩽ρn+1 h⩽˜ MinΩ 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:Wh→Whdefined as: −(hρh,¯ ρh)=(∇ρh,∇¯ ρh)∀¯ ρh∈Wh.(20) Then the discrete density equation (14) can be rewritten as: ρn+1 h−ρn h k,¯ ρh+wn h·∇ρn+1 h,¯ ρh−λhρn+1 h,¯ ρ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(Ω) > 0such that, for any ρh∈Wh,onehas: ∇ρhL3(Ω) ⩽C|∇ρh|1/2|hρh|1/2,(22) ∇ρhL4(Ω) ⩽Ch1/4|hρh|+ρh1/2 L∞(Ω)|hρh|1/2.(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 ∇ρhL4(Ω) ⩽Ch1/4|hρh|+ ρ(h)  1/2 L∞(Ω)|hρh|1/2,(24) where ρ(h)∈H2(Ω) is the solution of the following elliptic problem −ρ(h)=−hρhin Ω, ∂ρ(h) ∂n∂Ω =0, Ω ρ(h)(x)dx=0.(25) Indeed (24) is based on the inverse inequality ∇ρhL4(Ω) ⩽Ch−3/4|∇ρh|, the approximation property ∇ρ−∇IhρL4(Ω) ⩽ Ch1/4ρH2(Ω), 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H1and the interpolation error ρ−IhρL∞(Ω) ⩽Ch1/2ρH2(Ω). Now, we are in a position to prove some recursive inequalities for scheme (13)–(16). Lemma 8. Fixed n:0 ⩽n⩽N−1, assume 0<˜ m⩽ρn h,ρn+1 h⩽˜ MinΩ(26) and 1 16C2 λ2μkhρn h 2+1 2μk∇un h 2⩽Cd,(27) with Cd>0independent of (h,k,λ), and n. Then, provided that h and k are sufficiently small satisfying (S), there exists a unique solution (ρn+1 h,un+1 h,pn+1 h)of scheme (13)–(16) which satisfies: ⎧ ⎪ ⎨ ⎪ ⎩ρn+1 hun+1 h 2−ρn hun h 2+ρn hun+1 h−un h 2+3μ 4k∇un+1 h 2 ⩽kC1 fn+1  2 L6/5(Ω) +εμλ2khρn+1 h 2+hρn h 2, (28) λ∇ρn+1 h 2−λ∇ρn h 2+λ∇ρn+1 h−ρn h 2+λ2 2khρn+1 h 2⩽C2k∇un h 2,(29) where C1,C2,ε>0are constants independent of (h,k,λ),andn,withεbeing arbitrarily small. Proof. Taking ¯ uh=2kun+1 hand ¯ ph=pn+1 hin (15)–(16), and using the identity (a−b,2a)=|a|2−|b|2+|a−b|2and properties (17) and (18) gives: ⎧ ⎨ ⎩ρn+1 hun+1 h 2−ρn hun h 2+ρn hun+1 h−un h 2+μ−λ( ˜ M−˜ m)k∇un+1 h 2 ⩽kC1 fn+1  2 L6/5(Ω) +A, (30) where A=2kλ2 ˜ m∇ρn+1 hL4(Ω)∇ρn hL4(Ω)|∇un+1 h|and C1=C1(μ)>0 a constant independent of (h,k,λ), and n. In view of inequality (23) and hypothesis (26), the term Acan be bounded by A⩽Ckλ2h1/2hρn+1 hhρn h+˜ M1/2h1/4hρn+1 hhρn h 1/2+hρn hhρn+1 h 1/2 +˜ Mhρn+1 h 1/2hρn h 1/2∇un+1 h := A1+A2+A3+A4. By hypothesis (27), |hρn h|⩽4C1/2 2C1/2 d/(λμ1/2k1/2), hence the term A1is bounded as follows: A1⩽CC 1/2 d μ1/2h1/2k1/2λhρn+1 h∇un+1 h⩽kC0 ε1 h kλ2μhρn+1 h 2+ε1μk∇un+1 h 2, for any ε1>0 and C0=C0(μ). Thus, by stability constraint (S), we are allowed to take hand ksmall enough such that C0 ε1 h k⩽ε2,withε2>0 being arbitrary, to arrive at A1⩽ε2λ2μkhρn+1 h 2+ε1μk∇un+1 h 2. The term A2has a similar treatment as A1, 1176 R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ − T  0ˆ ρh,k,λ,d dt ˜ ηh,kdt +λ T  0 (∇ρh,k,λ,∇ηh,k)dt − T  0 (ˆ Hh,k,λρh,k,λ,∇ηh,k)dt =ρ0h,η0 h, (39) where vh,k,˜ vh,k(respectively, ηh,k,˜ ηh,k) are suitable approximations of a free-divergence test functions v∈C1([0,T]; C∞ c(Ω)) (respectively, C1([0,T];C∞ c(Ω)))withv(T)=0(respectively,η(T)=0) such that the sequence {vh,k,λ}h,k,λ is bounded in L∞(0,T;W1,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 J:=λ2 T  01 ρh,k,λ ∇ρh,k,λ ⊗∇ˆ ρh,k,λ,∇vh,k→0as(h,k,λ)→0. Indeed, since in particular vh,kis bounded in L2(0,T;W1,3(Ω)), one has J⩽Cλ1/2λ3 T  0 ∇ρh,k,λ4 L3(Ω)1/4λ3 T  0 ∇ ˆ ρh,k,λ4 L3(Ω)1/4T  0 ∇vh,k2 L3(Ω)1/2 ⩽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 ρn+1 h∇ρn+1 h⊗∇ρn h,∇¯ uh by the fully explicit approximation λ21 ρn h∇ρn h⊗∇ρn h,∇¯ uh, 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1 ρn+1 h∇ρn+1 h⊗∇ρn+1 h,∇¯ uh, then the term λ2(1 ρn+1 h∇ρn+1 h⊗∇ρn+1 h,∇un+1 h)is estimated using the discrete interpolation inequality (23) by CΩkλ2 ˜ mh1/2hρn+1 h 2+2˜ M1/2h1/4hρn+1 h 3/2+˜ Mhρn+1 h∇un+1 h. Then, it is not clear how to control the terms CΩkλ2 ˜ mh1/2hρn+1 h 2+2˜ M1/2h1/4hρn+1 h 3/2∇un+1 h 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 software FreeFem++ [12]. In all simulations, the approximating spaces are Vh=˜ Vh=P2and Mh=˜ Mh=P1for the velocity 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 Wh=P2or Wh=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∞(L2)and L∞(H1)norms of the density, velocity and pressure, by using P1elements 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 kof the same order than hattending to a computational cost and λ=0). This truncation consists in redefining the density ρn+1 hcalculated by the continuity equation (14) as ρn+1 h=χ(ρn+1 h), 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 Qwith homogeneous slip boundary conditions and the external force fand initial values such that (ρex,uex,pex)defined by ρex =2+xcossin(t)+ysinsin(t),(40) uex =−ycos(t), xcos(t),(41) pex =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 hand 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∞(L2)and L∞(H1)norms. Firstly, in Figs. 1 and 2, we measure the error in time for the velocity, pressure, and density (this latter is approximated by P1and P2finite elements). In both cases, we consider the mesh size h=0.0234458 and take time steps ksmaller 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∞(L2)norm and for the pressure in the L∞(H1)norm. We point out that the error in time in the L∞(L2)and L∞(H1)norms is comparable with the error measured in the L2(L2)and L2(H1)norms. Secondly, in Fig. 3, we fix the time step k=10−5and 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∞(L2) norm and of order 1 in the L∞(H1)norm. 1178 R.C. Cabrales et al. / Applied Numerical Mathematics 61 (2011) 1161–1185 Fig. 2. Errors in the L∞(L2)and L∞(H1)norms for the density, velocity, and pressure, by using P2elements for the density. We consider a fixed mesh size of h=0.0234458. Fig. 3. Errors in the L∞(L2)and L∞(H1)norms for the density, velocity and pressure, by using P2elements 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∞(L2)and L∞(H1)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∞(L2)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∞(H1)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 ρ0is given by ρ0=ρM+ρm 2+ρM−ρm 2tanh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∞(L2)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 =ρmd3/2G1/2 μ, where Gis the acceleration of the gravity. We have considered the following boundary conditions for 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 =20000 by using P1and P2 elements to approximate the density. The results are presented in Figs. 6, 7, 8, and 9. For P1elements approximating the density, our results are in good qualitative agreement with those reported in [3] and [6]. On the other hand, if we use P2for 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 P1for the density, this is not our case. However, in those calculations where we use a P2for 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 P1elements 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∞(H1)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 ρ+>ρ−>0aretwo 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Ω1and dΩ2, such that |dΩ|=|dΩ1|+|dΩ2|. Then, the volume fraction of each fluid in the point xis 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 P1elements. 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 P2elements. φ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 φ,wehave 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 P1elements. 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 P2elements. ρ=φρ++(1−φ)ρ−,μ=φρ+ν++(1−φ)ρ−ν−, where ν+,ν−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.8sec,forP1elements (left) and for P2 elements (right). ⎧ ⎨ ⎩ ρ(ut+u·∇u−λ∇log(ρ)·∇u)−∇·μ∇u+(∇u)t+∇P=ρfin Q, ∇·u=0inQ, ρt+∇·(ρu−λ∇ρ)=0inQ, (43) where λ>0 is a mass diffusion coefficient, Pis the pressure, uthe velocity, and f=g(sin θ,cos θ) is the external force (depending on the gravity gand 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.4sec,forP1elements (left) and for P2elements (right). Table 1 Values of the parameters used for numerical simulations (see [4]). Parameter Value Gravity acceleration, g9.8 m/s2 Slope, θ32 ◦C Reynolds number, Re 106 Schmidt number, Sc 0.3 Initial height, h01m Obstacle height, hs0.25h0m Obstacle thickness 0.05h0m Parameters for light and heavy fluids Heavy fluid Light fluid Density, ρ±4kg/m 31kg/m 3 Kinematic viscosity, ν±10−5m2/s10 −5m2/s Re := h0gh0 ν+,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: 1. The equation for calculating the approximate density ρn+1 his the same, 2. As we saw before, we truncate the approximate density as ρn+1 h=χ(ρn+1 h), 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 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 P1and P2approximation for the density, respectively. The results for P1are 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 2Dviscous 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 3DNavier–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 Lp(0,T;B), Ann. Mat. Pura Appl. 146 (1987) 65–97.