Novel fractional time-stepping algorithms for natural convection problems with variable density

Applied Numerical Mathematics 151 (2020) 64–84

Contents lists available at ScienceDirect

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

Novel fractional time-stepping algorithms for natural convection problems with variable density ✩ Jilian Wu a , Leilei Wei a , Xinlong Feng b,∗ a b

College of Science, Henan University of Technology, Zhengzhou, 450001, PR China College of Mathematics and Systems Science, Xinjiang University, Urumqi 830046, PR China

a r t i c l e

i n f o

Article history: Received 5 April 2019 Received in revised form 23 August 2019 Accepted 13 December 2019 Available online 31 December 2019 Keywords: Natural convection Variable density Fractional time-stepping technique Stability analysis Error estimate

a b s t r a c t This work presents novel fractional time-stepping finite element approaches for solving incompressible natural convection problems with variable density. Compared with other projection methods, the main merit of these methods is that we only need to solve one Poisson equation per time step for the pressure, which is computationally more efficient. First-order stability analyses and error estimates are proved and stability analysis of second-order version is established. Numerical tests which confirm the theoretical analyses are provided. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Natural convection (NC) problem, driven by the density differences occurring due to temperature gradients in fluids, has always been a hot topic in the field of heat transfer which is a branch of hydrodynamics. Since NC phenomena have been found in many fields, e.g., cooling of electronic components, fire control, room ventilation, nuclear reaction systems, double glass window design, solar energy device, etc., it is significant to study the problem. When the temperature variation is smaller, it can be modeled by using a Boussinesq approximation, which treats the density as a constant in terms other than the gravitational term. To the best our knowledge, researchers mostly focus on constant density model, see [6,19,25–27, 30,31,33]. However, in most geophysical flows, the temperature difference is often large enough that leads to considerable density variations under which the Boussinesq approximation is no longer valid. In these cases, we are led to consider the following model governed by the incompressible NC problems with variable density:

⎧ ρ + ∇ · (ρ u) = 0 ⎪ ⎪ ⎨ t ρ (ut + (u · ∇)u) − μu + ∇ p = ρ g ⎪∇ ·u=0 ⎪ ⎩ ρ ( T t + (u · ∇) T ) − κ  T = 0

in  × (0, T 1 ] in  × (0, T 1 ], in  × (0, T 1 ], in  × (0, T 1 ],

(a), (b) (c) (d)

(1.1)

✩ The work is supported in part by the High-Level Personal Foundation of Henan University of Technology (No. 2017BS028), the Fundamental Research Funds for the Henan Provincial Colleges and University in Henan University of Technology (No. 2018QNJH23), the Foundation of Henan Educational Committee (No. 19A110005) and the NSF of China grant (No. 11362021, 11671345, 11801146). Corresponding author. E-mail addresses: [email protected] (J. Wu), [email protected] (L. Wei), [email protected] (X. Feng).

*

https://doi.org/10.1016/j.apnum.2019.12.012 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

65

where the unknown variables are the density ρ > 0, the velocity vector u, the pressure p and the temperature T . The constants κ , μ, g and T 1 represent the thermal conductivity parameter, the dynamic viscosity coefficient, the gravitational acceleration and the fixed time, respectively. The fluid flows in a bounded domain  in Rd (d = 2 or 3) with a sufficiently smooth boundary ∂  and the solution to the above problem is found over a time interval (0, T 1 ]. The above control equations are derived, under assumptions that the fluids are non-homogeneous incompressible Newtonian fluids, that there is no external forces to do work except gravity, that there no sources of heat, that the specific heat at constant volume is a constant and that the energy depends only on the temperature, from the mass conservation, momentum conservation and the energy conservation, see [2,5,20]. Besides, some mathematical formal asymptotics are needed in deriving the (1.1d), see the Chapter 1 of [20] for more details. The system (1.1) is supplemented with the following initial and boundary conditions for ρ , u and T :

⎧ u(x, 0) = u0 (x) in  and u(x, t )| = g1 (x, t ), ⎪ ⎪ ⎪ ⎨ ρ (x, 0) = ρ0 (x) in  and ρ (x, t )|u(x,t ) = r (x, t ), ⎪  ⎪ ⎪ ⎩ T (x, 0) = T 0 (x) in  and ∂ T (x, t )  = 0, T (x, t )|2 = g 2 (x, t ), 1 ∂n

(1.2)

where  = ∂ , 1 is a regular open subset of , 2 = ∂  \ 1 and v is the inflow boundary defined by v = {x ∈  : v(x) · n < 0} for any velocity field v with n being the outward unit normal vector. Throughout this paper we assume that the boundary  is impermeable, i.e., u · n = 0 everywhere on  and v = ∅. For the first three equations of the system (1.1)-(1.2), there are precise results of mathematical theory of existence and uniqueness, and we refer the reader to see the Chapter 2 of [20]. Then it is possible to study mathematical theory of existence and uniqueness of the system (1.1)-(1.2), see, the Chapter 3.4 (and Remark 3.9) of [20] for more details. Furthermore, approximating (1.1)-(1.2) efficiently is a challenging task, which involves all the difficulties associated with the densitydependent Navier-Stokes (N-S) equations as well as the temperature equation [16,32]. Because of these difficulties, so far, very few papers have been dedicated to the mathematical analysis of the approximation of system (1.1)-(1.2). To the best of our knowledge, some papers, e.g., [28,32] have given the stability analysis of NC problems with variable density but the rigorous error analysis has not been proven yet. We attempt to develop efficient numerical methods and give mathematical analysis for (1.1)-(1.2) referring to Guermond and Salgado [16,17] in this paper. As we all know, projection methods [15] are very popular in the computational fluid dynamics community. Over the years, a large amount of efforts have been devoted to develop projection methods to solve incompressible flows with constant density and incompressible N-S equations with variable density (for example, [1,3,4,13,22] and the references therein). However, the common feature of all the projection methods referred to above is that the pressure or some related scalar quantity must be computed by solving an elliptic equation at each time step. Pleasantly, Guermond and Salgado [16,17] give an novel fractional time-stepping algorithm for solving N-S equations with variable density. The advantage of the method is that we only need to solve an Poisson equation per time step compared with projection methods. Based on the discussion above, we will develop the novel efficient fractional time-stepping method and use the FEM in space to solve NC problems with variable density. The methods require solving only one Poisson problem per time step for the pressure or some related scalar quantity. Inspired by the analysis of the full discrete time-stepping methods for solving incompressible N-S problems with variable density in [16,17], we will prove the stable and convergent analysis of the proposed schemes for NC problems with variable density. The remainder of this paper is organized as follows. Notations, along with space and time discretizations, are introduced in Section 2. We develop two first-order fractional time-stepping schemes, one in non-incremental form and the other in incremental form, for incompressible NC problems with variable density in Section 3. In Section 4, we prove the stability of two fractional time-stepping algorithms. And the error estimates of two algorithms are proved in Section 5. In Section 6, second-order version of the method is introduced and the stability of algorithm is proved. Then numerical experiments illustrating the performance of the methods are reported in Section 7. Finally, we end with a short conclusion in Section 8. 2. Preliminaries For the mathematical setting of problem (1.1)-(1.2), we introduce standard Hilbert spaces, some notations and some necessary assumptions which will be frequently used in the following sections. Note in passing that vector-valued functions and vector-valued spaces are identified with bold fonts. 2.1. Notation and assumptions Firstly, we introduce the following standard Hilbert spaces:

X := H 01 () = { v ∈ H 1 () : v = 0 on ∂ },



Q := L 20 () = {q ∈ L 2 () :

q dx = 0}, 

66

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

W := H 1 (),

X := X d , Y := ( L 2 ())d ,

H = {v ∈ Y : ∇ · v = 0, v · n|∂  = 0}. Further, we will consider the standard definitions for Sobolev spaces W m, p () equipped with norm  · m, p and semi-norm | · |m, p (m ≥ 0, p ≥ 1). The closure with respect to the norm  · W m, p of the space of C ∞ -functions compactly supported in  m, p

m ,2

is denoted by W 0 (). To simplify the notation, we use the following notations H m () := W m,2 (), H m 0 () := W 0 (), and  · m =  · m,2 , | · |m = | · |m,2 when p = 2. And W 0, p () = L p () when m = 0. The L 2 -scalar product and L 2 -norm of L 2 () := H 0 () are denoted by ·, · and  · 0 , respectively. Then we assume the boundary  of bounded domain  to be sufficiently smooth. Besides, we assume that system (1.1)-(1.2) has a unique smooth solution for all T 1 > 0 and that all the compatibility conditions required for the solution to be smooth enough are satisfied. Let τ > 0 be a time step and set tn = nτ for 0 ≤ n ≤ N = T 1 /τ , where · is the floor function. The time-discrete approximations of (ρ (tn ), u(tn ), p (tn ), T (tn )) will be denoted by (ρ n , un , pn , T n ). Let E be a normed space equipped with the norm  ·  E . For any time-dependent function θ : [0, T 1 ] → E, we denote θ n := θ(tn ), and the sequence θ 0 , θ 1 , · · · , θ n is denoted by θτ . We define the following discrete norms:

θτ  2 ( E ) := (τ

N 

θ n 2E )1/2 ,

θτ  ∞ ( E ) := max (θ n  E ).

n =0

(2.1)

0≤n≤ N

The space of functions θ : [0, T 1 ] → E such that the map (0, T 1 )  t → θ(t ) E ∈ R is L p -integrable is indifferently denoted by L p ((0, T 1 ); E ) or L p ( E ), see [17]. In the following sections, we use c to denote a generic constant whose value may change at each occurrence for convenient and simplicity. This constant may depend on the data of the problem and its exact solution, but it does not depend on the discretization parameters or the solution of the numerical scheme. 2.2. The space discretization Let Th = { K } be a uniformly regular family of triangulation of , and define the mesh size h = max {diam( K )}. To K ∈Th

construct a Galerkin approximation of (1.1)-(1.2), we consider the following four sequences of finite dimensional spaces





¯ ) : sh | K ∈ P 2 ( K ), ∀ K ∈ Th , M h = sh ∈ W ∩ C 0 (



¯ )d : vh | K ∈ P 2 ( K )d , ∀ K ∈ Th , Xh = vh ∈ X ∩ C 0 (



¯ ) : qh | K ∈ P 1 ( K ), ∀ K ∈ Th , Q h = qh ∈ Q ∩ C 0 ( Wh =





¯ ) : ϕh | K ∈ P 2 ( K ), ∀ K ∈ Th , ϕh ∈ X ∩ C 0 (

where P i ( K ) is the set of all polynomials on K of degree less than i ∈ N . Note that the pair of spaces (Xh , Q h ) satisfy a discrete inf-sup condition (cf. [8,9]); i.e., there is c > 0 independent of h such that



inf

sup

0=qh ∈ Q h 0=vh ∈Xh

· ∇ qh d  ≥ c. ∇ vh 0 qh 0  vh

According to [8,9,17], we similarly assume that the finite element space pair ( M h , Xh , Q h , W h ) satisfies following approximation properties: Lemma 2.1. There exists l ∈ N ∗ such that for all ∈ [0, l] the density space satisfies

inf s − sh 0 ≤ ch +1 s H +1 , ∀s ∈ H +1 (),

sh ∈ M h

(2.2)

the velocity-pressure spaces satisfy for all

inf v − vh 0 + hv − vh H1 ≤ ch +1 vH +1 , ∀v ∈ H +1 () ∩ H10 (),

(2.3)

inf q − qh 0 ≤ ch q H , ∀q ∈ H () ∩ L 20 (),

(2.4)

vh ∈Xh

qh ∈ Q h

and the temperature space satisfies

inf ϕ − ϕh 0 ≤ ch +1 ϕ  H +1 , ∀ϕ ∈ H +1 ().

ϕh ∈ W h

(2.5)

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

67

Now, for any t ∈ [0, T 1 ] we define the ‘Stokes projection’ ( S hu(t ), I h p (t )) ∈ Xh × Q h of the solution (u(t ), p (t )) to (1.1)-(1.2) through solving the equation



μ ∇ S hu(t ), ∇ vh + ∇ I h p (t ), vh = μ ∇ u(t ), ∇ vh − p (t ), ∇ · vh , ∀vh ∈ Xh , S hu(t ), ∇ qh = 0,

∀qh ∈ Q h .

(2.6)

Similarly, we define the ‘elliptic projection’ J hT (t ) ∈ W h of the solution T (t ) to (1.1)-(1.2) by requiring

∇ J hT (t ), ∇ ϕh = ∇ T (t ), ∇ ϕh , ∀ϕh ∈ W h .

(2.7)

From the results of [17,21], and owing to the regularization properties of the Stokes operator and the elliptic operator, the following estimates hold: Lemma 2.2. If u ∈ L β (Hl+1 () ∩ H10 ()), p ∈ L β ( H l ()), and T ∈ L β ( H l+1 () ∩ H 01 ()) for 1 ≤ β ≤ ∞, then there exists c > 0 such that

u − S huL β (L2 ) + h(u − S huL β (H1 ) + p − I h p L β ( L 2 ) ) ≤ chl+1 (uL β (Hl+1 ) + p L β (Hl ) ),

(2.8)

 T − J h T L β ( L 2 ) + h( T − J h T L β ( H 1 ) ≤ chl+1  T L β ( H l+1 ) .

(2.9)

and

Moreover, if u ∈

L β (H2 () ∩ H10 ()), p

∈

L β ( H 1 ()), and T

∈

L β ( H 2 () ∩ H 01 ()) then

 S huL β (L∞ ∩W1,3)+ I h p L β ( H 1 )+ J h T L β (L ∞ ∩W 1,3) ≤ c (uL β (H2)+ p L β ( H 1 ) + T L β (H 2) ).

(2.10)

We end this section with the introduction of the following algebraic identities to treat time derivative terms [17,23]: Lemma 2.3 (Inner product of time derivative terms). We use δ as difference of two consecutive functions, for example, for any sequence

{ Z n }nN=0 ,

δ Z n+1 = Z n+1 − Z n , δδ Z n+1 = δ(δ Z n+1 ) = Z n+1 − 2 Z n + Z n−1 , · · · , then we have

2 Z n+1 − Z n , Z n+1 =  Z n+1 20 −  Z n 20 +  Z n+1 − Z n 20 , 2 3 Z n+1 − 4 Z n + Z n−1 , Z n+1 =  Z n+1 20 + 2 Z n+1 − Z n 20

+δ 2 Z n+1 20 −  Z n 20 − 2 Z n − Z n−1 20 ,

(2.11) (2.12)

and

2 3 Z n+1 − 4 Z n + Z n−1 , Z n+1 = 3 Z n+1 20 − 4 Z n 20 +  Z n−1 20

+2δ Z n+1 20 − 2δ Z n 20 + δ 2 Z n+1 20 .

(2.13)

3. Description of first-order fractional time-stepping method In this section, we will develop first-order fractional time-stepping schemes (which are non-incremental and incremental schemes) for incompressible NC problems with variable density. We refer the reader to [16] for the heuristics for N-S problems with variable density behind the method. It is well known that the non-incremental pressure correction method is low-order accurate for constant density N-S equations, see [11,24]. To overcome this accuracy deficiency, Guermond and Salgado [16,17] introduce an incremental version of the method to solve the N-S equations with variable density. We now develop novel fractional methods to solve NC problems with variable density. sup Firstly, we define ρ0min := minx∈¯ ρ0 (x), ρ0 := supx∈¯ ρ0 (x). Henceforth we assume that ρ0min > 0 and that the approx0 imate density field ρh satisfies

χ ≤ ρh0 ≤ , where the parameters

(3.1)

χ and  are assumed to satisfy the following property:

χ ∈ (0, ρ0min ],  ∈ [ρ0sup , +∞).

(3.2)

68

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

The role of the parameters χ and  is clarified in the next subsection. We shall be more specific later about the properties sup that we expect from ρh0 and p h0 . In the numerical computation reported in section 7 we take χ = ρ0min ,  = ρ0 . Now, we construct the semi-discrete scheme of system (1.1)-(1.2) by defining the approximate sequences {ρ n }n=0,··· , N , n {u }n=0,··· , N , { pn }n=0,··· , N , and { T n }n=0,··· , N as follows: Algorithm 3.1 (Semi-discrete scheme). We set ρ 0 = ρ0 , u0 = u0 , T 0 = T 0 , and note that p 0 := p |t =0 is not part of the initial data but this quantity can be computed from (1.1b) with the known ρ0 and u0 . We refer the reader to the Remark 3.1 of [17] for more details. Then we repeat the following steps for all time index n ranging 0 to N − 1: Step 1. Update ρ n+1 as the solution of

ρ n +1 − ρ n ρ n +1 + ∇ · (ρ n+1 un ) − ∇ · un = 0. τ 2

(3.3)

Step 2. The velocity un+1 is now updated by

⎧ ∗ n+1 n n ⎨ ρ u − ρ u + ρ n+1 (un · ∇)un+1 + 1 ∇ · (ρ n+1 un )un+1 − μun+1 +∇ p  = ρ n+1 g, ⎩

τ

2

where ρ ∗ = 12 (ρ n+1 + ρ n ), p  = pn + γ δ pn , Step 3. Find φ n+1 as the solution of

γ ∈ {0, 1}.

⎧ ⎨ φ n+1 = χ ∇ · un+1 , ⎩

(3.4)

un+1 | = 0,

∂n φ

n +1

τ

(3.5)

| = 0 ,

then we can get

p n +1 = φ n +1 + γ p n .

(3.6)

Step 4. Update T n+1 as the solution of

⎧ ∗ n +1 −ρn T n 1 ⎨ρ T + ρ n+1 (un ·∇) T n+1+ ∇ ·(ρ n+1un ) T n+1 − κ  T n+1 = 0, ⎩

τ

2

(3.7)

T n + 1 | = 0 .

The equation (3.3) is obtained by approximating (1.1a) using a first-order semi-implicit discretization. The additional term

ρ n+1 2

∇ · un is consistent since it is zero if ∇ · un = 0, and its meaning will become clear when we do the stability analysis.

Remark 3.1. Note that (3.5) is a standard Poisson equation. This is the main novelty of the new fractional time-stepping technique in this paper compared with other projection methods which need to solve an elliptic equation with variable coefficient. Remark 3.2. The Algorithm 3.1 unifies the non-incremental (γ = 0) and the incremental (γ = 1) first-order semi-discrete schemes of fractional time-stepping method for NC problem with variable density. Moreover, note that for non-incremental scheme, we can write (3.3), (3.4) and (3.7) as

ρ n +1 − ρ n ρ n +1 + un · ∇ ρ n+1 + ∇ · un = 0, τ 2 ⎧ n n+1 n n +1 ⎨ ρ (u − u ) + ρ n+1 (un · ∇)un+1 + ρ (∇ · un )un+1 − μun+1 +∇ pn = ρ n+1 g, τ 4 ⎩

(3.8)

(3.9)

un+1 | = 0,

and

⎧ n n +1 − T n) ρ n +1 ⎨ ρ (T + ρ n+1 (un ·∇) T n+1+ (∇ · un ) T n+1 − κ  T n+1 = 0, ⎩

τ

T n + 1 | = 0 ,

respectively.

4

(3.10)

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

69

Remark 3.3. For the incremental scheme (γ = 1), the pressure term in the momentum equation (3.4) is a second-order extrapolation of the pressure at time t n+1 , since using (3.6) at time step n we have

p  = pn + φ n = 2pn − pn−1 . Remark 3.4. The term

ρ ∗ un+1 −ρ n un + 12 ∇ · (ρ n+1 un )un+1 is asymptotically consistent with the equation. Notice that τ

ρ ∗ un+1−ρ n un 1 + ∇ · (ρ n+1 un )un+1 = (ρ ut )n+1 + O (τ ), τ 2 ∗ n+1 n n ρ T −ρ T 1 + ∇ · (ρ n+1 un ) T n+1 = (ρ T t )n+1 + O (τ ) τ 2 if the involved functions are sufficiently smooth, and see [16] for more details. The purpose of this particular way of discretizing the quantity ρ ut and ρ T t is to further decouple the mass conservation and the momentum conservation equations. This will become clear once we do the stability analysis. Then we will give a full-discrete scheme of the Algorithm 3.1 in the following. First of all, given the initial data (ρ0 , u0 , T 0 ), we construct the approximate data (ρh0 , uh0 , ph0 , T h0 ) ∈ M h × Xh × Q h × W h , and we assume that ρh0 , uh0 and T h0 satisfies (3.1)-(3.2) and the following estimate holds:

u0 − uh0 0 + h∇(u0 − uh0 )0 ≤ chl+1 ,  T 0 − T h0 0 + h∇( T 0 − T h0 )0 ≤ chl+1 .

(3.11)

Algorithm 3.2 (Full-discrete scheme). Given (ρhn , unh , pnh , T hn ) ∈ M h × Xh × Q h × W h , we obtain the next approximations

(ρhn+1 , unh+1 , pnh+1 , T hn+1 ) ∈ M h × Xh × Q h × W h from the following four steps:

Step 1. Finding a suitable algorithm to approximate the mass conservation returns the following stability hypothesis:

ρhn+1 ∈ M h and this algorithm satisfies

χ ≤ min ρhn+1 (x), sup ρhn+1 (x) ≤  ∀n ≥ 1. ¯ x∈

Step 2. Defining

(3.12)

¯ x∈

ρh∗ = 12 (ρhn+1 + ρhn ) and ph = pnh + γ δ pnh to find unh+1 ∈ Xh for ∀vh ∈ Xh is such that 

 1   ρh∗ unh+1 − ρhn unh , vh + ρhn+1 (unh · ∇)unh+1 , vh + ∇ · (ρhn+1 unh )unh+1 , vh τ 2 







+ ∇ ph , vh + μ

∇ unh+1 , ∇ vh





n +1 g, vh h

= ρ



(3.13)

.

Step 3. Finding φhn+1 ∈ Q h for ∀qh ∈ Q h satisfies



 χ  ∇φhn+1 , ∇ qh = unh+1 , ∇ qh ,

(3.14)

τ

then we get

pnh+1 = φhn+1 + γ pnh .

(3.15)

Step 4. Finding T hn+1 ∈ W h for ∀ϕh ∈ W h is such that



 T n +1  ρh∗ T hn+1 − ρhn T hn , ϕh + ρhn+1 (unh ·∇) T hn+1 , ϕh + h ∇ ·(ρhn+1unh ), ϕh τ 2   + κ ∇ T hn+1 , ∇ ϕh = 0.

 (3.16)

Remark 3.5. The ρhn+1 is computed using the mass conservation equation which is hyperbolic in Step 1. It is generally known that Galerkin techniques are not well suited for the solution of hyperbolic problems [8,17]. However, there are many techniques aiming at addressing this issue [7,17], e.g., Galerkin least squares, discontinuous Galerkin, characteristics method, edge stabilization, entropy viscosity, subgrid viscosity, etc. We choose the technique of characteristics method [7,28] to obtain approximate density.

70

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Remark 3.6. Note that (3.13) and (3.16) can be rewritten as

ρhn

unh+1 − unh

τ



 n+1      ρh − ρhn 1 , vh + + ∇ · (ρhn+1 unh ) unh+1 , vh + ∇ ph , vh 

2

n +1 n (uh h

+ ρ

τ

· ∇)unh+1 , vh





+μ

∇ unh+1 , ∇ vnh+1





n +1 g, vh h

= ρ

 ,

(3.17)

and





  ρhn ( T hn+1 − T hn ) , ϕh + ρhn+1 (unh ·∇) T hn+1 , ϕh + κ ∇ T hn+1 , ∇ ϕh τ

 n+1   ρh − ρhn 1 + + ∇ ·(ρhn+1unh ) T hn+1 , ϕh = 0, 2 τ

(3.18)

respectively. These forms emphasize the stabilizing term, which are also used in [17,18]. Remark 3.7. For the non-incremental scheme (γ = 0), if also be rewritten as

ρhn

unh+1 − unh

τ

ρhn+1 −ρhn ρ n +1 + ∇ · (ρhn+1 unh ) − h2 ∇ · unh = 0 holds, (3.13) and (3.16) can τ



    1 , vh + ρhn+1 (∇ · unh )unh+1 , vh + ∇ ph , vh 

4

n +1 n (uh h

+ ρ

· ∇)unh+1 , vh





+μ

∇ unh+1 , ∇ vnh+1





n +1 g, vh h

= ρ

 ,

(3.19)

and





  ρhn ( T hn+1 − T hn ) , ϕh + ρhn+1 (unh ·∇) T hn+1 , ϕh + κ ∇ T hn+1 , ∇ ϕh τ

 +

1 4

(3.20)

ρhn+1 (∇ · unh ) T hn+1 , ϕh = 0,

respectively. Remark 3.8. For the incremental scheme (γ = 1), by analogy with projection methods for constant density flows [14], we can write a rotational version of the above algorithm by replacing the pressure update (3.15) by

pnh+1 = pnh + φhn+1 − μ∇ · unh+1 .

(3.21)

4. Stability analysis of first-order fractional time-stepping method In this section, we will give some proofs for the stability of first-order fractional time-stepping schemes. We consider the non-incremental scheme and incremental scheme separately. Since the stability analyses of the semi-discrete and fulldiscrete schemes can use the same procedure, apart from condition (3.12), we will only give the full-discrete version in the following subsections. We henceforth denote by  · 0 the L 2 -norm. No notational distinction is made between the norm in L 2 () and Y. 4.1. Stability analysis of the non-incremental scheme In this subsection, we will give a proof for the stability of the non-incremental version, which is obtained by setting

γ = 0. We start with the L 2 -stability of the density [4,13,16]. Lemma 4.1. For all τ > 0, h > 0, 0 ≤ N ≤ T 1 /τ − 1 and any sequence of velocities {un }n=0,1,··· , N in L∞ () with bounded divergence and satisfying un · n| = 0, the solution to the solution of (3.3) satisfies

ρ N 20 +

N −1 

ρ n+1 − ρ n 20 = ρ 0 20 ,

n =0

and if (3.1) holds, the numerical density ρ n determined by (3.3) satisfies

χ ≤ ρ n ≤ . Then we establish the stability of step 2-step 4 in Algorithm 3.2 in the following under the assumption of (3.12).

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

71

Theorem 4.1. Assume that (3.12) holds. Then, for any τ > 0, h > 0 the solution to (3.13)-(3.16) satisfies the following inequality:

σhN uhN20 +σhN T hN20 + μτ

N  n =1

+ 2κτ where σhn :=

N 

N −1

 τ2  [ ∇ pnh20 + ∇ pnh 20 ] χ n =1 n =0 N

∇ unh 20 +

∇ T hn 20 ≤ σh0 uh0 20 + C T 1 ρh0 20 + σh0 T h0 20 ,

n =1



ρhn .

Proof. Firstly, using Lemma 4.1, the Cauchy-Schwarz inequality and v0 ≤ C ∇ v0 (∀v ∈ H10 ()), we obtain

2τ ρhn+1 g, unh+1 ≤ 2τ g0 ρhn+1 0 unh+1 0 ≤ 2C τ g0 ρhn+1 0 ∇ unh+1 0

(4.1)

≤ C τ ρhn+1 20 + μτ ∇ unh+1 20 ≤ C τ ρh0 20 + μτ ∇ unh+1 20 .

Similarly to the proving process of Theorem 2.1 of [16], we then get the following priori bound of velocity and pressure

σhn+1 unh+1 20 − σhn unh 20 + μτ ∇ unh+1 20 + Secondly, we let

τ2 [∇ pnh+1 20 + ∇ pnh 20 ] ≤ C τ ρh0 20 . χ

(4.2)

ϕh = 2τ T hn+1 in equation (3.20), integrate by parts, and use the identity (2.11), we then get

σhn T hn+1 20 − σhn T hn 20 + σhn ( T hn+1 − T hn )20 + 2κτ ∇ T hn+1 20   τ + τ ρhn+1 unh · ∇| T hn+1 |2 + ρhn+1 (∇ · unh )| T hn+1 |2 = 0.

(4.3)

2





Next, multiplying the density equation obtain

σhn+1 T hn+1 20 − σhn T hn+1 20 − τ

ρhn+1 −ρhn τ



+ ∇ · (ρhn+1 unh ) −

ρhn+1 unh · ∇| T hn+1 |2 −



τ

ρhn+1 2



2

∇ · unh = 0 by τ | T hn+1 |2 and integrating by parts, we

ρhn+1 (∇ · unh )| T hn+1 |2 = 0.

(4.4)



Adding up Eqs. (4.3) and (4.4), we get

σhn+1 T hn+1 20 − σhn T hn 20 + σhn ( T hn+1 − T hn )20 + 2κτ ∇ T hn+1 20 = 0. Lastly, adding up (4.2) and (4.5), and adding up over n = 0, · · · , N − 1 gives the desired stability result.

(4.5)

2

4.2. Stability analysis of the incremental scheme In this subsection, we will give a proof for the stability of the incremental version, which is obtained by setting γ = 1. Note that equations used to determine the density and the temperature are the same as the non-incremental case. Therefore, the L 2 -stability of the density is again consequence of Lemma 4.1. Theorem 4.2. Assume that (3.12) holds, then, for any τ > 0, h > 0 the solution to (3.13)-(3.16) satisfies the following inequality:

σhN uhN 20 +σhN T hN 20 + μτ

N  n=1

+

N−1 2 

τ χ

∇ unh 20 + 2κτ

N 

∇ T hn 20 +

n=1

∇( pnh − pnh−1 )20 ≤ σh0 uh0 20 +σh0 T h0 20 +

n =1

τ2 ∇ phN 20 χ

τ2 ∇ ph0 20 + C T 1 ρh0 20 . χ

Proof. Firstly, in addition to (4.1), using the proving process of Theorem 3.1 of [16], we then obtain the following priori bound of velocity and pressure

σhn+1 unh+1 20 −σhn unh 20+

τ2 [∇ pnh+1 20 −∇ pnh 20 +∇( pnh − pnh−1 )20 ] χ + μτ ∇ unh+1 20 ≤ C τ ρh0 20 .

(4.6)

72

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Secondly, similarly, we let

 2τ

ϕh = 2τ T hn+1 in (3.16) and integrate by parts. Notice that the boundary conditions imply





1

ρhn+1 (unh · ∇) T hn+1 + ∇ · (ρhn+1 unh ) T hn+1 · T hn+1 2



=τ  =

 

ρhn+1 unh · ∇| T hn+1 |2 + ∇ · (ρhn+1 unh )| T hn+1 |2





  ∇ · ρhn+1 | T hn+1 |2 unh = 0.



Using (2.11), we have



ρh∗ T hn+1 − ρhn T hn , 2τ T hn+1 = ρhn+1 T hn+1 , T hn+1 + 2 ρhn ( T hn+1 − T hn ), T hn+1 − ρhn T hn+1 , T hn+1

τ = σhn+1 T hn+1 20 + σhn ( T hn+1 − T hn )20 − σhn T hn 20 .

Then we obtain

σhn+1 T hn+1 20 + σhn ( T hn+1 − T hn )20 − σhn T hn 20 + 2κτ ∇ T hn+1 20 = 0. The desired result is obtained by adding up (4.6) and (4.7) for n = 0, · · · , N − 1.

(4.7)

2

Remark 4.1. The stability property from Theorem 4.2 does not explicitly require the pair of spaces (Xh , Q h ) to satisfy the LBB condition. However, the LBB condition must be invoked to prove stability on the pressure in L 2 () when going through the details. We refer the reader to, e.g. [10,12,15,16] for more details. 5. Error estimate of first-order fractional time-stepping method In this section, we prove that the Algorithm 3.2 is convergent provided that the mass conservation equation is approximated appropriately. Firstly, to simplify the notation we introduce the following functions to represent the errors:



e(t ) := u(t ) − S hu(t ),

η(t ) := p (t ) − I h p (t ),ξ(t ) := T (t ) − J hT (t ),

enh

ηhn := I h pn − pnh ,

n

:= S hu

− unh ,

ξhn := J hT n − T hn .

(5.1)

The functions e(t ), η(t ) and ξ(t ) can be regarded as the interpolation errors, whereas the functions enh , ηhn and ξhn represent the approximation errors. In addition to (3.11), we make the following regularity assumptions on the exact solution of problem (1.1):

u ∈ W2,∞ (H10 () ∩ Hl0+1 ()), p ∈ W 1,∞ ( H l ()), T ∈ W 2,∞ ( H 01 () ∩ H l0+1 ()), where we have already assumed that l ≥ 1 in (2.2)-(2.5). Let us now determine the equations that control the errors. We first recall the equations [17] that control enh and taking the difference between (2.6), (3.13) and (3.14) at time step t n+1 :

(5.2)

ηhn by



    ρh∗ enh+1 − ρhn enh  , vh + μ ∇ enh+1 , ∇ vnh+1 + ∇( I h pn+1 − ph ), vh = Rn1+1 (vh ), ∀vh ∈ Xh , τ    χ    ∇ ηh , ∇ qh = enh+1 , ∇ qh + ∇ I h p  , ∇ qh , τ

(5.3) (5.4)

where the residual Rn1+1 (vh ) is decomposed as follows:

Rn1+1 (vh ) =

4 

R n1i+1 (vh ),

i =1

R n11+1 (vh ) :=



ρhn

S hun+1 − S hun

(5.5)

 − ρ n+1 unt +1 , vh ,

τ  ρ n +1 − ρ n  1 h h S hun+1 − ρtn+1 un+1 , vh , R n12+1 (vh ) := 2 τ   R n13+1 (vh ) := g(ρ n+1 − ρhn+1 ), vh ,

(5.6) (5.7) (5.8)

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84



R n14+1 (vh ) := ρhn+1 (unh · ∇)unh+1 − ρ n+1 (un+1 · ∇)un+1 , vh

+



1 2

 ∇ · (ρhn+1 unh )unh+1 − ∇ · (ρ n+1 un+1 )un+1 , vh ,

n+1 h

and where ηh = η − γ ηhn , I h p  = I h pn+1 − γ I h pn , the version is the nonincremental if γ = 0 and incremental if Secondly, we similarly obtain the equation that control ξhn by combining (2.7) and (3.16) at time step t n+1 : 

73

(5.9)

γ = 1.



3    ρh∗ ξhn+1 − ρhn ξhn , ϕh + κ ∇ξhn+1 , ∇ ϕhn+1 = Rn2+1 (ϕh ) = R n2i+1 (ϕh ), ∀ϕh ∈ W h , τ i =1

(5.10)

where the residual Rn2+1 (ϕh ) is decomposed as follows:

R n21+1 (ϕh ) :=



ρhn

J hT n+1 − J hT n

 − ρ n+1 T tn+1 , ϕh ,

(5.11)

τ  ρ n +1 − ρ n  1 h h R n22+1 (ϕh ) := J hT n+1 − ρtn+1 T n+1 , ϕh , 2 τ   n +1 n n +1 R 23 (ϕh ) := ρh (uh · ∇) T hn+1 − ρ n+1 (un+1 · ∇) T n+1 , ϕh +

1 2

(5.12)

 ∇ · (ρhn+1 unh ) T hn+1 − ∇ · (ρ n+1 un+1 ) T n+1 , ϕh .

(5.13)

Equations (5.3)-(5.13) will be used repeatedly in the error analysis. The error analysis is based on energy arguments and the first of these arguments consists of testing (5.3) and (5.10) with vh := 2τ enh+1 and ϕh := 2τ ξhn+1 , respectively. Then using the equalities 2 (ρh∗ enh+1 − ρhn enh ), enh+1 = σhn+1 enh+1 20 +

σhn δ enh+1 20 − σhn enh 20 and 2 (ρh∗ ξhn+1 − ρhn ξhn ), ξhn+1 = σhn+1 ξhn+1 20 + σhn δξhn+1 20 − σhn ξhn 20 , we can get

   σhn+1 enh+1 20 − σhn enh 20 + σhn δ enh+1 20 + 2μτ enh+1 21 + 2τ ∇ ηh , enh+1   = 2τ ∇( I h p  − I h pn+1 ), enh+1 + 2τ Rn1+1 (enh+1 ),

(5.14)

σhn+1 ξhn+1 20 − σhn ξhn 20 + σhn δξhn+1 20 + 2μτ ξhn+1 21 = 2τ Rn2+1 (ξhn+1 ).

(5.15)

Next, we will give an estimate on the consistency residual Rn1+1 (enh+1 ) and Rn2+1 (ξhn+1 ) in the following lemma. Lemma 5.1. Assume that the solution to (1.1)-(1.2) satisfies (5.2) and that the sequence of approximate densities {ρhn } satisfies (3.12). Then

|Rn1+1 (enh+1 )| ≤ c [τ + hl+1 +ρhn − ρ n 0 +

δ ρhn+1

1

− ρtn+1  H −1 ]2 + μ∇ enh+1 20 + c σhn enh 20 , τ 4 n+1 δρ 1 |Rn2+1 (ξhn+1 )| ≤ c [τ + hl+1 +ρhn − ρ n 0 + h − ρtn+1  H −1 ]2 + μ∇ξhn+1 20 + c σhn enh 20 . τ 4

(5.16) (5.17)

Proof. Firstly, we can prove (5.16) using the method slightly different from those proposed by Guermond et al. [17] because of the additional term R n13+1 (enh+1 ) ≤ c ρhn+1 − ρ n+1 0 enh+1 1 in this paper. Analogously, together with (2.9), (2.10), (3.12), (5.2), Cauchy inequality, Hölder inequality and the identity ρ u · ∇ ϕ , ϕ + 1 ∇ · (ρ u)ϕ , ϕ = 0 with ϕ = ξh , we have the following estimates: 2

R n21+1 (ξhn+1 ) =



1

ρhn δ J hT n+1 − ρ n+1 T tn+1 , ξhn+1 τ 



  1 = ρhn ( δ J hT n+1 − T tn+1 ), ξhn+1 + (ρhn − ρ n ) T tn+1 , ξhn+1 − δ ρ n+1 T tn+1 , ξhn+1

τ   1 ≤ c ξhn+1 L 6 ρhn L ∞  δ J hT n+1 − T tn+1 0 +(ρhn − ρ n 0 +δ ρ n+1 0 ) T tn+1 L 3 τ   n +1 l +1 ≤ c ξh 1 τ + h + ρhn − ρ n 0 ,

 1 1 n +1 R n22+1 (ξhn+1 ) = δ ρh J hT n+1 − ρtn+1 T n+1 , ξhn+1 2 τ

      1 1 n +1 1 = δ ρh − ρtn+1 J hT n+1 , ξhn+1 + ρtn+1 J hT n+1 − T n+1 , ξhn+1 2

τ

2

74

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

1

≤ c  δ ρhn+1 − ρtn+1  H −1  J hT n+1 ξhn+1  H 1 + ξhn+1 L 6 ρtn+1 L 3  J hT n+1 − T n+1 0 0 τ   1 n +1 n +1 n +1 l +1 ≤ c ξh 1 h +  δ ρh − ρt  H −1 , τ 

1 R n23+1 (ξhn+1 ) = − ρhn+1 (enh · ∇) J hT n+1 + ∇ · (ρhn+1 enh ) J hT n+1 , ξhn+1 2 

  1 + (ρhn+1 − ρ n+1 )( S hun · ∇) J hT n+1 + ∇ · (ρhn+1 − ρ n+1 ) S hun J hT n+1 , ξhn+1 2  n +1 n n +1 n +1 − (u · ∇) T n+1 ] + ρ [( S hu · ∇) J hT  1 + [∇ · (ρ n+1 S hun ) J hT n+1 − ∇ · (ρ n+1 un+1 ) T n+1 ], ξhn+1 2   ≤ c σhn enh 0 + ρhn+1 − ρ n+1 0 + (τ + hl+1 ) ξhn+1 1 . Therefore, summing the above inequalities, we can easily get the estimate (5.17). Then we complete the proof of Lemma 5.1. 2 Remark 5.1. To avoid lengthy technicalities of little interest and to remain as general as possible for the mass conservation equation approximation, we are not going to do the full convergence analysis. We shall also assume instead that density sequence {ρhn } ⊂ M h such that for m > 0:

ρ − ρh 2 2 (0,tn ; L 2 ) +ρt −

δ ρh

τ

2 2 (0,tn+1 ; H −1 ) ≤ c (λ)(τ + hm )2

(5.18)

+ λehτ 2 2 (0,tn+1 ;H1 ) + c (λ)σhτ ehτ 2 2 (0,tn ;L2 ) , where λ ≥ 0 can be chosen as small as necessary [17]. And we shall use characteristic method expecting that density sequence satisfies (5.18) for m > 1 in section 7. Having the above assumption, Rn1+1 (enh+1 ) and Rn2+1 (ξhn+1 ) can be simplified as follows: Corollary 5.1. Assume that (5.18) holds, then the following estimate holds under the regularity assumptions of Lemma 5.1:

2τ

N  n =0

2τ

N  n =0

1

|Rn1+1 (enh+1 )| ≤ c (τ + hmin{l+1,m})2 + μehτ 2 2 (0,tn+1 ;H1 ) + c σhτ ehτ 2 2 (0,tn ;L2 ) .

(5.19)

2

1

|Rn2+1 (ξhn+1 )| ≤ c (τ + hmin{l+1,m})2 + μξhτ 2 2 (0,tn+1 ; H 1 ) + c σhτ ehτ 2 2 (0,tn ;L2 ) .

(5.20)

2

Proof. From Lemma 5.1 and (5.18), we can easily prove the results using the method of [17].

2

5.1. Error estimate of the non-incremental scheme We now give the main error estimate of non-incremental version (γ = 0) of Algorithm 3.2 under assumption (5.18). Now, we firstly assume that the algorithm is initialized so that p h0 satisfies the following estimate:

 ph0 0 ≤ c .

(5.21)

Theorem 5.1. Assume that the solution to (1.1)-(1.2) satisfies (5.2). Let (uh )τ , ( T h )τ be the solution of (3.13)-(3.16) with γ = 0, and assume that (3.11), (3.12), (5.18) and (5.21) hold, then 1

1

uτ − uhτ  ∞ (L2 ) ≤ c (hmin{l+1,m} + τ 2 ), uτ − uhτ  2 (H1 ) ≤ c (hmin{l,m} + τ 2 ), 1 2

1 2

 T τ − T hτ  ∞ ( L 2 ) ≤ c (hmin{l+1,m} + τ ),  T τ − T hτ  2 ( H 1 ) ≤ c (hmin{l,m} + τ ).

(5.22) (5.23)

Proof. In this case p  = pnh and φhn+1 = pnh+1 . We will prove Theorem 5.1 using the method similar to [17] in the following. We proceed in following two steps:

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

75

(i) Initialization (n = 0): We consider the initial step separately, as it involves the initial quantities. Firstly, using (5.14) for k = 0 and the hypothesis (5.21), we get

3

σh1 eh1 20 + μτ eh1 21 ≤ σh0 eh0 20 + c τ + 2τ |R11 (eh1 )|, 2

and we refer the reader to [17] for more details. Analogously, using (5.15) for k = 0, we get

σh1 ξh1 20 + 2μτ ξh1 21 ≤ σh0 ξh0 20 + 2τ |R12 (ξh1 )|. (ii) General step (n ≥ 1): Using the proving method of Theorem 4.1 in [17], we get

σhτ ehτ  ∞ (L2 ) + ehτ  2 (H1 ) ≤ c (τ 1/2 + hmin{l+1,m} ).

(5.24)

On the other hand, using (5.15), we get

σhn+1 ξhn+1 20 + 2μτ ξhn+1 21 ≤ σhn ξhn 20 + 2τ |Rn2+1 (ξhn+1 )|.

(5.25)

By summing over the time steps, using (5.20) and (5.24), we obtain N N   (σhn+1 ξhn+1 20 + μτ ξhn+1 21 ) ≤ c (τ + h2 min{l+1,m} )+ σhn ξhn 20 , n =0

(5.26)

n =0

which, by the discrete Gronwall lemma, implies

σhτ ξhτ  ∞ ( L 2 ) + ξhτ  2 ( H 1 ) ≤ c (τ 1/2 + hmin{l+1,m} ).

(5.27)

Lastly, using the triangle inequality, the definition u − =e + T − =ξ + and (in the case of the

∞ (L2 )-norm, ∞ ( L 2 )-norm) assumption (3.12), we can get the claimed error estimates. Notice that it is only at this point that the interpolation error in the H1 -norm, which is of order O (hl ), is introduced. This is a well-known superconvergence effect induced by our particular choice for the pair ( S hu, I h p ); see (2.6) and [17,29]. 2 n

unh

n

enh ,

n

T hn

n

ξhn

Conjecture 5.1. We expect that further regularity assumptions combined with a standard duality argument, e.g., multiplying the error equation by Senh+1 , where S is the solution operator to the time-independent Stokes problem, should allow us to conclude that the following estimates [4,12,14,17] hold in addition to (5.22) and (5.23):

uτ − uhτ  2 (L2 ) ≤ c (hmin{l+1,m} + τ ),

 T τ − T hτ  2 ( L 2 ) ≤ c (hmin{l+1,m} + τ ).

5.2. Error estimate of the incremental scheme This subsection will show the error estimate of the incremental scheme for Algorithm 3.2 (γ = 1) based on assumption (5.18). According to [17], we assume that either the initial pressure is properly approximated, i.e.,

 ph0 − I h p 0 0 ≤ ch

l +1 2

(5.28)

, 1

or that the approximation of the initial pressure is uniformly bounded in H with respect to h, and the initial approximation of the velocity is discretely divergence free, i.e.,

 ph0 1 ≤ c , uh0 , ∇ qh = 0, ∀qh ∈ Q h .

(5.29)

Under the above conditions, we can get the following error estimate for the incremental version of Algorithm 3.2. Theorem 5.2. Assume that the solution to (1.1)-(1.2) satisfies (5.2). Let (uh )τ , ( T h )τ be the solution of (3.13)-(3.16) with γ = 1, and assume that (3.11), (3.12) and (5.18) hold. If either (5.28) or (5.29) holds, then

uτ − uhτ  ∞ (L2 ) ≤ c (hmin{l+1,m} + τ ), uτ − uhτ  2 (H1 ) ≤ c (hmin{l,m} + τ ),  T τ − T hτ  ∞ ( L 2 ) ≤ c (h

min{l+1,m}

+ τ ),  T τ − T hτ  2 ( H 1 ) ≤ c (h

min{l,m}

+ τ ).

(5.30) (5.31)

Proof. In this case p  = 2pnh − pnh−1 and φhn+1 = δ pnh+1 and we will prove the theorem in two steps. (i) Initialization (n = 0): We consider the initial step separately, as it involves the initial quantities. From (5.14), (2.8), (5.28), the Cauchy-Schwarz inequality and the condition (5.29) and using the same technique in [17], we get

σh1 eh1 20 + μτ eh1 21 ≤ c τ 2 + σh0 eh0 20 + 2τ |R11 (eh1 )|.

76

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Besides, using (5.15) for k = 0, we get

σh1 ξh1 20 + 2μτ ξh1 21 ≤ σh0 ξh0 20 + 2τ |R12 (ξh1 )|. (ii) General step (n ≥ 1): Just as in the proof of Theorem 4.2 in [17] and Theorem 5.1 in this paper, using (3.12), (5.14), (5.19), Conjecture 5.1, Cauchy-Schwarz inequality and the discrete Gronwall lemma, we get

σhτ ehτ  ∞ (L2 ) + ehτ  2 (H1 ) ≤ c (τ + hmin{l+1,m} ).

(5.32)

What’s more, by summing over the time steps, using (5.20) and Conjecture 5.1, we obtain N N   (σhn+1 ξhn+1 20 + μτ ξhn+1 21 ) ≤ c (τ + hmin{l+1,m} )2 + σhn ξhn 20 , n =0

(5.33)

n =0

which, by the discrete Gronwall lemma, implies

  σhτ ξhτ  ∞ ( L 2 ) + ξhτ  2 ( H 1 ) ≤ c τ + hmin{l+1,m} .

(5.34)

Then we finish the proof by the same technique used in the proof of Theorem 5.1 in the following step.

2

6. Second-order fractional time-stepping method In this section, we first develop second-order semi-discrete and fully-discrete fractional time-stepping schemes for incompressible NC problems with variable density in subsection 6.1. Then, we prove stability of a second-order fully-discrete fractional time-stepping scheme based on finite element method for NC problems with variable density in subsection 6.2. Now, we recall some three term recursion inequalities [17] which will be needed to prove stability in subsection 6.2. Lemma 6.1. Assume that the characteristic polynomial of the three term recursion equation

Axn+1 + Bxn + C xn−1 = g n+1 , n ≥ 2,

(6.1)

has two (not necessarily distinct) nonzero real roots r1 and r2 . Then, the generic solution to this equation has the form

xν = c 1 r1ν + c 2 r2ν +

ν l 1  ν −l  l−s s r1 r2 g , A l =2

c 1 , c 2 ∈ R.

s =2

Lemma 6.2. Assume that the coefficients of the three term recursion inequality

Ayn+1 + B yn + C yn−1 ≤ g n+1 , n ≥ 1,

(6.2)

satisfy

A > 0, C ≥ 0, A + B + C ≤ 0. Let { y }k≥0 be a solution to (6.2) with initial data y 0 and y 1 . If {xn }k≥0 solves (6.1) with initial data x0 = y 0 and x1 = y 1 , then the following estimate holds: n

y ν ≤ xν ∀ν ≥ 0. 6.1. Description of the numerical algorithm Keeping the same notation as in the previous sections, firstly, we give the second-order semi-discrete scheme of the algorithm in the following steps. Algorithm 6.1 (Semi-discrete scheme). Initialization. First, we choose a penalty parameter χ as in (3.2). For initial data (ρ 0 , u0 , p 0 , φ 0 = 0, T 0 ), we assume that minx∈¯ ρ 0 (x) > 0 and that ρ 0 satisfies (3.1). Then we compute an approximation (ρ 1 , u1 , p 1 , φ 1 = p 1 − p 0 , T 1 ) at time t = τ by using Algorithm 3.1. Besides, we define for simplicity

⎧ 3 n +1 2 n 1 1 ⎪ ∗ ⎪ − ρ + ρ n−1 = ρ n+1 + (3ρ n+1 − 4ρ n + ρ n−1 ), (a) ⎪ ⎨ ρ = 2ρ 3 6 6 4

1

p ∗ = p n + φ n − φ n −1 , ⎪ ⎪ 3 3 ⎪ ⎩ ∗ u = 2un − un−1 .

(b) (c)

(6.3)

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

77

Then, given {ρ n , un , pn , φ n , T n }n=0,..., N −1 we advance in time in the following steps: Step 1. Updating ρ n+1 as the solution of

⎧ n +1 ⎪ − 4ρ n + ρ n −1 ρ n +1 ⎨ 3ρ + u∗ · ∇ ρ n+1 + ∇ · u∗ = 0, 2τ 2 ⎪ n +1 ⎩ ρ n +1 | . u¯ n+1 = r

(6.4)

Step 2. The velocity un+1 is now updated by

⎧ ∗ n +1 3ρ u − 4ρ n+1 un + ρ n+1 un−1 un+1 ⎪ ⎪ + ρ n+1 (u∗ · ∇)un+1 + ∇ · (ρ n+1 u∗ ) ⎪ ⎨ 2τ 2 + ∇ p ∗ − μun+1 = ρ n+1 g, ⎪ ⎪ ⎪ ⎩ n +1 u | = gn1+1 .

(6.5)

Step 3. Finding φ n+1 as the solution of

⎧ ⎨ φ n+1 = 3χ ∇ · un+1 , 2τ ⎩ ∂ n φ n + 1 | = 0 .

(6.6)

Step 4. Updating pn+1 by

p n +1 = p n + φ n +1 . Step 5. Updating T

n+1

(6.7) by

⎧ ∗ n +1 − 4ρ n +1 T n + ρ n +1 T n −1 ⎪ T n +1 ⎨ 3ρ T + ρ n+1 (u∗ ·∇) T n+1+ ∇ ·(ρ n+1u∗)− κ T n+1 = 0, 2τ

2

⎪ ⎩ T n +1 | = g n +1 .  2

(6.8)

Then, we can construct fully-discrete version of Algorithm 6.1. Algorithm 6.2 (Full-discrete scheme). Initialization. Given the initial data (ρ 0 , u0 , T 0 ), we construct the approximate data (ρh0 , uh0 , p h0 , T h0 ) ∈ M h × Xh × Q h × W h , and we assume that ρh0 satisfies (3.1) and that uh0 and T h0 satisfy the following estimates:

u0 − uh0 0 + h∇(u0 − uh0 )0 ≤ chl+1 ,  T 0 − T h0 0 + h∇( T 0 − T h0 )0 ≤ chl+1 .

(6.9)

Then (ρh1 , uh1 , p h1 , φh1 = p h1 − p h0 , T h1 ) are computed by Algorithm 3.2. Besides, (ρh∗ , p h∗ , uh∗ ) is the discrete form of (ρ ∗ , p ∗ , u∗ )

defined in (6.3). Now, given (ρhn , unh , pnh , φhn = pnh − pnh−1 , T hn ) ∈ M h × Xh × Q h × Q h × W h for 1 ≤ n ≤ N − 1, we compute the next time-step approximation in the following steps. Step 1. We are not specific about the way ρhn+1 ∈ M h is computed, but we assume that (3.12) and that there is a uniform constant M so that

max

0≤n≤ N −1



ρhn+1 − ρhn L∞ ≤ M χ . τ

(6.10)

Step 2. Updating unh+1 ∈ Xh for ∀vh ∈ Xh so that the following holds:

3ρ ∗ unh+1 − 4ρhn+1 unh + ρhn+1 unh−1 2τ



n +1 ∗ (uh h

+ ρ

· ∇)unh+1

+

unh+1 2

   , vh + μ ∇ unh+1 , ∇ vh n +1 ∗ uh ), vh h

∇ · (ρ



    + ∇ ph∗ , vh = ρhn+1 g, vh .

(6.11)

Step 3. We compute the pressure correction φhn+1 ∈ Q h for ∀qh ∈ Q h so that the following holds:



 3χ  n +1  ∇φhn+1 , ∇ qh = uh , ∇ qh . 2τ

(6.12)

Step 4. Then the pressure is updated by setting

pnh+1 = pnh + φhn+1 .

(6.13)

78

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Step 5. Finally, we update T hn+1 for ∀ϕh ∈ W h such that

3ρh∗ T hn+1 − 4ρhn+1 T hn + ρhn+1 T hn−1 2τ

  T hn+1 n+1 ∗ n +1 n+1 ∗ , ϕh + ρh (uh ·∇) T h + ∇ ·(ρh uh ), ϕh 2  n +1  + κ ∇ T h , ∇ ϕh = 0.

(6.14)

  T n+1  ∇ · (ρhn+1 uh∗ ), vh and h2 ∇ · (ρhn+1uh∗ ), ϕh have been added to the equations to obtain uncondi un+1 tional stability with respect to the advection terms. A main advantage of adding the terms is ρhn+1 (uh∗ · ∇)unh+1 + h2 ∇ ·    T n+1 (ρhn+1 uh∗ ), unh+1 = 0 and ρhn+1 (uh∗ ·∇) T hn+1 + h2 ∇ ·(ρhn+1uh∗ ), T hn+1 = 0, and we refer the reader to [17] for the proof.

Remark 6.1. Terms

 un+1 h

2

3ρ ∗ un+1−4ρ n+1 un+ρ n+1 un−1

un+1

3ρ ∗ T n+1−4ρ n+1 T n+ρ n+1 T n−1

h h h h h h h h h Remark 6.2. Terms + h2 ∇ · (ρhn+1 uh∗ ) and h h + 2τ 2τ approximation of [ρh uh,t ](t n+1 ) and [ρh T h,t ](t n+1 ), respectively (see [17] for the proof).

T hn+1 ∇ 2

·(ρhn+1uh∗ ) is a second-order

Remark 6.3. We can construct a rotational version of Algorithm 6.2 by replacing the pressure update (6.13) by the following: find pnh+1 ∈ Q h such that













pnh+1 , qh = pnh + φhn+1 , qh + μ unh+1 , ∇ qh .

(6.15)

6.2. Stability analysis The objective of this subsection is to prove the stability of the second-order fractional time-stepping scheme in standard form incompressible NC problems with variable density. The novel proof approaches and techniques are firstly used to prove stability estimates of the time-dependent Stokes equation with constant density, and then used to prove stability estimates of N-S equations with variable density [17]. The main novelty is that the proof does not use the projected velocity. Theorem 6.1. Assume that (3.12) and (6.10) hold. Then, for any inequality:

σhN uhN20 + μτ

N 

N −1

 τ2  [ ∇ pnh 20 + ∇ pnh 20 ] ≤ σh0 uh0 20 + C T 1 ρh0 20 , χ n =1 n =0 N

∇ unh 20 +

n =1

τ > 0, h > 0 the solution to (6.11)-(6.14) satisfies the following

(6.16)

and

σhN T hN20 + where σhn :=

N −1 

σhn+1 T hn+1 − σhn T hn 20 + 2κτ

∇ T hn 20 ≤ σh0 T h0 20 .

(6.17)

n =1

n =0



N 

ρhn .

Proof. Firstly, we can prove slightly different from those proposed by Guermond et al. [17] because  (6.16) using the method  of the additional term 2τ g(ρ n+1 − ρhn+1 ), unh+1 ≤ c τ ρh0 20 + μτ ∇ unh+1 20 in this paper. Next, we will give the estimate on temperature in a way similar. As already mentioned above, the time derivative can be rewritten as follows:

3ρ ∗ T hn+1 − 4ρhn+1 T hn + ρhn+1 T hn−1 2τ

= ρhn+1

3T hn+1 − 4T hn + T hn−1 2τ

Using Assumption (6.10), we have the following estimate:

  (3ρhn+1 − 4ρhn + ρhn−1 ) T hn+1 , T hn+1   = 3 (ρhn+1 − ρhn )| T hn+1 |2 − (ρhn − ρhn−1 )| T hn+1 |2 

≥ −(3

 n +1 − hn h

ρ

ρ

χ

L ∞ +

≥ −4M τ  T hn+1 20 .

n − hn−1 h

ρ

ρ χ

L ∞ )σhn+1 T hn+1 20

+

T hn+1 3ρhn+1 − 4ρhn + ρhn−1 2

2τ

.

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

79

Using (2.13), Assumption (6.10) and the inequality −(a − b)2 ≥ −2(a2 + b2 ), a similar treatment gives







2 ρhn+1 (3T hn+1 − 4T hn + T hn−1 ), T hn+1 = n+1 n+1 2 T h 0 + h

ρhn+1 2T hn+1 (3T hn+1 − 4T hn + T hn−1 )





n+1 (| T hn−1 |2 − 4| T hn |2 − 2|δ T hn |2 )+ 2 hn+1 δ T hn+1 20 + hn+1 δ 2 T hn+1 20 h

=3σ

ρ

σ

 n+1 n+1 2 T h 0 h

=3σ 

+ 2σhn+1 δ T hn+1 20 + σhn+1 δ 2 T hn+1 20 +

n n −1 2 | − h |T h

+

n −1 n −1 2 | T h | +  hn−1 T hn−1 20 h

ρ

σ









ρhn | T hn−1 |2





n +1 n 2 |T h | + 4 h

−4



ρhn | T hn |2

ρ 



ρhn+1 |δ T hn |2 + 2

− 4σhn T hn 20 − 2 n+1 n+1 2 T h 0 h

ρhn+1 | T hn−1 |2 −





ρ 



σ



ρhn |δ T hn |2 − 2σhn δ T hn 20 

+ σhn+1 δ 2 T hn+1 20 − 2σhn δ T hn 20 − (4 + 4M τ )σhn T hn 20  n −1 n −1 2 + (1 − 2M τ )σh T h 0 − 2 (ρhn+1 − ρhn )| T hn − T hn−1 |2

≥3σ

n+1 δ T hn+1 20 h

+ 2σ

 n+1 n+1 2 T h 0 h

≥3σ

+ (1 − 6M τ

+ 2 hn+1 δ T hn+1 20 ) hn−1 T hn−1 20 .

σ

+ σhn+1 δ 2 T hn+1 20 − 2σhn δ T hn 20 − (4 + 8M τ )σhn T hn 20

σ

Combining the above two inequalities, we get

2 3ρ ∗ T hn+1 − 4ρhn+1 T hn + ρhn+1 T hn−1 , T hn+1

≥ (3 − 4M τ )σhn+1 T hn+1 20 − (4 + 8M τ )σhn T hn 20 + (1 − 6M τ )σhn−1 T hn−1 20 n +1 δ T hn+1 20 h

+ 2σ

n n 2 h δ T h 0

− 2σ

(6.18)

n+1 2 n+1 2 δ T h 0 . h

+ σ

Then we proceed in two steps as above: first we investigate the time steps n = 1, 2; then we investigate the cases n ≥ 3. (i) Initialization: We consider the steps n = 1, 2, as they involve the initial quantities. For n = 1 or n = 2 we set ϕh = 4τ T hn+1 in (6.14). Using (3.12), (6.10), (6.18), Remark 3.1, the Cauchy-Schwarz inequality and Young inequality, we obtain that

σhn+1 T hn+1 20 + 4κτ ∇ unh+1 20 ≤ c (σh0 T h0 20 + σh1 T h1 20 ), if

τ is small enough. (ii) General step: For n ≥ 3 notice that, we set

ϕh = 4τ T hn+1 in (6.14), from (6.18) we get

(3 − 4M τ )σhn+1 T hn+1 20 − (4 + 8M τ )σhn T hn 20 + (1 − 6M τ )σhn−1 T hn−1 20 + 2σhn+1 δ T hn+1 20 − 2σhn δ T hn 20 + σhn+1 δ 2 T hn+1 20 + 4κτ ∇ unh+1 20 ≤ 0.

(6.19)

Introducing the notation

A :=3 − 4M τ , n

B := −(4 + 8M τ ),

a := hn unh 20 , n ≥ 0, bn :=4 ∇ unh 20 +  hn δ 2 T hn 20 , dn :=2 hn δ T hn 20 , n ≥ 2,

C = 1 − 6M τ ,

σ

κτ

σ

n ≥ 1,

σ

inequality (6.19) can be rewritten as

Aan+1 + Ban + Can−1 ≤ −(bn+1 + dn+1 − dn ),

n ≥ 3.

Define g n+1 := −(bn+1 + dn+1 − dn ). If τ is small enough, this three term recursion inequality satisfies the assumptions of Lemma 6.2. The roots of the characteristic polynomial are

r1 := r2 :=

2 + 4M τ −

√

1 + 38M τ − 8M τ 2

3 − 4M τ

2 + 4M τ +

√

1 + 38M τ − 8M τ 2

3 − 4M τ

1

41M τ

3

3

= (1 −

+ O (τ 2 )),

= 1 + 9M τ + O (τ 2 )).

80

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Table 1 Error in time for the nonincremental scheme.

τ

1/10

1/20

1/40

1/80

1/160

1/320

ρ − H1

2.9552E−1 / 4.4014E−3 / 9.7743E−2 / 4.8160E−2 / 6.4360E−1 / 1.2023E−2 / 5.9351E−3 / 79.96

1.5098E−1 0.968 2.1141E−3 1.058 7.0957E−2 0.462 3.4539E−2 0.480 4.9015E−1 0.393 5.5623E−3 1.112 2.6859E−3 1.144 154.65

8.1741E−2 0.885 1.0115E−3 1.064 4.7535E−2 0.578 2.2625E−2 0.610 3.3683E−1 0.541 2.6500E−3 1.070 1.2521E−3 1.101 303.32

4.7154E−2 0.794 5.1979E−4 0.960 2.9609E−2 0.683 1.3606E−3 0.734 2.1098E−1 0.675 1.3164E−3 1.009 6.2251E−4 1.008 710.23

2.8538E−2 0.725 2.8336E−4 0.875 1.7500E−2 0.759 7.6546E−3 0.830 1.2293E−1 0.779 6.6645E−4 0.982 3.2204E−4 0.951 1549.62

1.6033E−2 0.832 1.3772E−4 1.041 1.0026E−2 0.804 4.1104E−3 0.897 6.8605E−2 0.8415 3.3753E−4 0.982 1.6866E−4 0.933 2969.77

Order ρ − L2 Order u − H1 Order u − L2 Order p − L2 Order T − H1 Order T − L2 Order CPU(s)

Both roots are positive; the first one is strictly less than one third, and the second one is greater but close to one. Then, in a way similar to the proof of the estimate on velocity in [17], we can get

aν +

1 ν 1 b + dν ≤ K(1 + ecT 1 )(a1 + a2 ). 3 3

This inequality, combined with the estimates obtained at the initialization step, implies the result.

2

7. Numerical experiments In this section, we present two numerical experiments to illustrate numerically algorithms presented in the previous sections. Throughout this section, we use the finite element spaces ( P 2 , P 2 , P 1 , P 2 ) for the density, the velocity vector, the pressure and the temperature. For the following numerical implementations, we choose g = (0, −9.8), κ = 1, χ = 1. The first one is a flow problem with manufactured analytical solutions. And the second one is a Bénard convection problem. 7.1. Examples with analytical solution Firstly, in order to test the accuracy of algorithms proposed in this paper, we consider a known analytical solution in an unit circle:

ρ (x1 , x2 , t ) = 2 + x1 cos(sin(t )) + x2 sin(sin(t )), u 1 (x1 , x2 , t ) = −x2 cos(t ), u 2 (x1 , x2 , t ) = x1 cos(t ), p (x1 , x2 , t ) = sin(x1 ) sin(x2 ) sin(t ), T (x1 , x2 , t ) = u 1 (x1 , x2 , t ) + u 2 (x1 , x2 , t ).

μ = 1, T 1 = 1, and the initial conditions of (1.1)-(1.2) are determined by the above exact solutions. log( E i / E i +1 ) with respect to the time step τ , where E i and E i +1 log(τi /τi +1 ) are the relative errors corresponding to the time steps τi and τi +1 , respectively. Computations are made on a fixed small We choose

We carry out the accuracy tests calculated by the formula

enough mesh size with different time steps so that the spatial discretization error can be negligible compared with the time error. The results are given in Tables 1–4. From Tables 1–2, we observe that the numerical results with the first-order fractional time-stepping FEM in two versions behave similarly and consistent with the Theorem 5.1, Theorem 5.2 and the results of [32]. So we only present the numerical results of first-order scheme in one form in the following test. Then we give the results of the second-order standard and rotational version of the method in Tables 3–4, respectively. The rotational version of the algorithm is that consists of using the pressure update (6.15), introduced in Remark 6.3, instead of (6.13). We observe that the convergence order of error on the velocity in the H 1 -norm and on the pressure in the L 2 -norm of rotational version is higher than that of standard formulation. 7.2. Bénard convection In this subsection, we consider the Bénard convection which is a classical fluid dynamic phenomenon. In many applications such as geophysical flows, density varies under the influence of temperature can not be ignored, which needs to use

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

81

Table 2 Error in time for the first-order incremental scheme.

τ

1/10

1/20

1/40

1/80

1/160

1/320

ρ − H1

2.9602E−1 / 4.3787E−3 / 3.2519E−2 / 1.4455E−2 / 1.6127E−1 / 1.2968E−2 / 6.7706E−3 / 79.40

1.5166 E−1 0.965 2.3600E−3 0.892 9.9463E−3 1.709 3.9573E−3 1.869 6.4951E−2 1.312 6.6456E−3 0.965 3.5205E−3 0.943 154.06

8.2536E−2 0.888 1.4370E−3 0.716 4.5683E−3 1.123 2.0635E−3 0.939 4.1691E−2 0.640 3.4135E−3 0.961 1.8250E−3 0.948 303.56

4.7987E−2 0.782 8.9102E−4 0.690 2.6349E−3 0.794 1.2928E−3 0.675 2.6365E−2 0.661 1.7378E−3 0.974 9.3521E−4 0.965 602.44

2.9182E−2 0.718 5.0695E−4 0.814 1.4351E−3 0.877 7.1907E−4 0.846 1.4736E−2 0.839 8.7646E−4 0.987 4.7330E−4 0.983 1203.01

1.6450E−2 0.827 2.6434E−4 0.939 7.4471E−4 0.946 3.7835E−4 0.926 7.8604E−3 0.907 4.4041E−4 0.993 2.3832E−4 0.990 2497.06

Order ρ − L2 Order u − H1 Order u − L2 Order p − L2 Order T − H1 Order T − L2 Order CPU(s)

Table 3 Error in time for second-order standard scheme.

τ

1/10

1/20

1/30

1/40

1/50

1/60

ρ − H1

1.4251E−1 / 1.9872E−3 / 1.9897E−2 / 8.4647E−3 / 9.7382E−2 / 3.2131E−3 / 1.6467E−3 /

3.7859E−2 1.912 4.7475E−4 2.066 5.6374E−3 1.819 1.9429E−3 2.123 2.7515E−2 1.823 7.6282E−4 2.075 3.9129E−4 2.073

1.7343E−2 1.925 2.0751E−4 2.041 2.8217E−3 1.707 8.3768E−4 2.075 1.3843E−2 1.694 3.3422E−4 2.035 1.7159E−4 2.033

9.9478E−3 1.932 1.1577E−4 2.029 1.7576E−3 1.646 4.6464E−4 2.049 8.6629E−3 1.629 1.8682E−4 2.022 9.5955E−5 2.020

6.4544E−3 1.939 7.3733E−5 2.021 1.2279E−3 1.607 2.9506E−4 2.035 6.0741E−3 1.591 1.1915E−4 2.015 6.1215E−5 2.014

4.5293E−3 1.943 5.1042E−5 2.017 9.2054E−4 1.580 2.0390E−4 2.027 4.5675E−3 1.564 8.2566E−5 2.012 4.2426E−5 2.011

Order ρ − L2 Order u − H1 Order u − L2 Order p − L2 Order T − H1 Order T − L2 Order

Table 4 Error in time for second-order rotational scheme.

τ

1/10

1/20

1/30

1/40

1/50

1/60

ρ − H1

1.4201E−1 / 1.6811E−3 / 1.0655E−2 / 5.3087E−3 / 4.6198E−2 / 2.7936E−3 / 1.4285E−3 /

3.7692E−2 1.914 4.2041E−4 2.000 2.9269E−3 1.864 1.4382E−3 1.884 1.1684E−2 1.983 7.0029E−4 1.996 3.5884E−4 1.993

1.7245E−2 1.928 1.8820E−4 1.982 1.3738E−3 1.865 6.6667E−4 1.896 5.2557E−3 1.970 3.1384E−4 1.979 1.6098E−4 1.977

9.8792E−3 1.937 1.0652E−4 1.978 8.0239E−4 1.869 3.8506E−4 1.908 2.9959E−3 1.954 1.7762E−4 1.979 9.1150E−5 1.977

6.4021E−3 1.944 6.8512E−5 1.978 5.2820E−4 1.874 2.5094E−4 1.919 1.9398E−3 1.948 1.1418E−4 1.980 5.8613E−5 1.979

4.4872E−3 1.949 4.7768E−5 1.978 3.7514E−4 1.877 1.7660E−4 1.927 1.3617E−3 1.941 7.9554E−5 1.982 4.0847E−5 1.981

Order ρ − L2 Order u − H1 Order u − L2 Order p − L2 Order T − H1 Order T − L2 Order

Fig. 1. Bénard Convection: domain and boundary conditions.

82

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

Fig. 2. Isolines of density and temperature, and velocity of Bénard convection based on first-order nonincremental scheme with variable density T h = 1, T c = 0, τ = 0.0001, Ra1 = 105 , h = 1/120 at different times: (a) t = 0.02, (b) t = 0.045, (c) t = 0.065 (from left to right). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 3. Isolines of density and temperature, and velocity of Bénard convection based on second-order rotational scheme with variable density T h = 1, T c = 0, τ = 0.001, Ra1 = 102 , h = 1/40 at different times: (a) t = 0.202, (b) t = 0.502, (c) t = 1.002 (from left to right).

Fig. 4. Isolines of density and temperature, and velocity of Bénard convection based on first-order scheme with variable density T h = 100, T c = 0, τ = 0.01, Ra1 = 10, h = 1/40 at different times: (a) t = 2.02, (b) t = 5.02, (c) t = 7.02 (from left to right).

the model (1.1). The domain and boundary conditions of Bénard convection problem [32] are displayed in Fig. 1. The top and bottom boundaries are isothermal, maintained respectively at T c (cold) and T h (hot); both vertical walls are adiabatic. And the boundary conditions for the velocity are no-slip at all boundaries. The Rayleigh number (Ra), which describes the relationship between buoyancy and viscosity within a fluid and can also be used as a criterion to predict convectional instabilities, is the most important dimensionless number for NC problem with small density variations. Thus, we introduce a similar number (Ra1 ) in the right hand side of (1.1b) to control the magnitude of buoyancy force for NC problem with large density variations in the following tests.

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

83

Fig. 5. Isolines of density and temperature, and velocity of Bénard convection based on second-order rotational scheme with variable density T h = 100, T c = 0, τ = 0.01, Ra1 = 10, h = 1/40 at different times: (a) t = 2.02, (b) t = 5.02, (c) t = 7.02 (from left to right).

Fig. 6. Isolines of density and temperature, and velocity of Bénard convection based on second-order rotational scheme with variable density T h = 100, T c = 0, τ = 0.01, μ = 0.1, Ra1 = 1, h = 1/40 at different times: (a) t = 4.02, (b) t = 8.02, (c) t = 11.02 (from left to right).

Firstly, we consider the case of T h = 1, T c = 0. The initial conditions are T 0 = 0, u0 = 0. The density, temperature and streamline functions of first-order nonincremental scheme with ρ 0 = 0.01(x2 − 1), μ = 1, Ra1 = 105 , τ = 0.0001, h = 1/120 are shown in Fig. 2 at difference times: t = 0.02, t = 0.045, t = 0.065. Then the numerical results of second-order rotational scheme with ρ 0 = (x2 + 1), μ = 1, Ra1 = 100, τ = 0.001, h = 1/40 are presented in Fig. 3 at difference times: t = 0.202, t = 0.502, t = 1.002. Secondly, we consider the case of T h = 100, T c = 0. The initial conditions are T 0 = 0, u0 = 0 and ρ 0 = (x2 + 1). The numerical results of first-order scheme are shown Fig. 4 at difference times: t = 2.02, t = 5.02, t = 7.02. Then the numerical results of second-order scheme are shown Figs. 5–6. From Figs. 2–3, we can see that our method in this paper conforms with Wu’s [32]. And from Figs. 2–6, we can see that the bigger the Ra1 is, the earlier stable the fluid tends to be. 8. Conclusion We firstly developed two first-order novel fractional time-stepping algorithms, one in non-incremental form and the other in incremental form, for the NC problem with variable density, and proved their stability and convergence. Then, we gave the second-order fractional time-stepping algorithm and proved the stability. Finally, numerical tests illustrated that algorithms in this paper are stable and reliable to deal with NC problem with variable density. Acknowledgements The authors would like to thank the editor and referees for their valuable comments and suggestions which helped us to improve the results of this paper. References [1] A. Almgren, J. Bell, P. Colella, L. Howell, M. Welcome, A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations, J. Comput. Phys. 142 (1998) 1–46.

84

J. Wu et al. / Applied Numerical Mathematics 151 (2020) 64–84

[2] G. Batchelor, An Introduction to Fluid Dynamics, vol. 1, Cambridge University Press, 1967, pp. 4–252. [3] J. Bell, D. Marcus, A second-order projection method for variable-density flows, J. Comput. Phys. 101 (1992) 334–348. [4] H. Chen, J. Mao, J. Shen, Error estimate of Gauge-Uzawa methods for incompressible flows with variable density, J. Comput. Appl. Math. 364 (2020) 112321. [5] A. Chorin, J. Marsden, A Mathematical Introduction to Fluid Mechanics, vol. 3, Springer, 1990, pp. 1–45. [6] G. De Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, Int. J. Numer. Methods Fluids 3 (3) (1983) 249–264. [7] J. Douglas Jr., T. Russell, Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM J. Numer. Anal. 19 (1982) 871–885. [8] A. Ern, J. Guermond, Theory and Practice of Finite Elements, Appl. Math. Sci., vol. 159, Springer-Verlag, New York, 2004. [9] V. Girault, P. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer Ser. Comput. Math., Springer-Verlag, Berlin, 1986. [10] J. Guermond, Some practical implementations of projection methods for Navier-Stokes equations, M2AN Math. Model. Numer. Anal. 30 (5) (1996) 637–667. [11] J. Guermond, L. Quartapelle, On the approximation of the unsteady Navier-Stokes equations by finite element projection methods, Numer. Math. 80 (5) (1998) 207–238. [12] J. Guermond, Un résultat de convergence d’ordre deux en temps pour l’approximation des équations de Navier-Stokes par une technique de projection incrémentale, M2AN Math. Model. Numer. Anal. 33 (1) (1999) 169–189. Also in C. R. Acad. Sci., Sér. 1 Math. 325 (1997) 1329–1332. [13] J. Guermond, L. Quartapelle, A projection FEM for variable density incompressible flows, J. Comput. Phys. 165 (1) (2000) 167–188. [14] J. Guermond, J. Shen, On the error estimates for the rotational pressure-correction projection methods, Math. Comput. 73 (248) (2004) 1719–1737. [15] J. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (44–47) (2006) 6011–6045. [16] J. Guermond, A. Salgado, A splitting method for incompressible flows with variable density based on a pressure Poisson equation, J. Comput. Phys. 228 (2009) 2834–2846. [17] J. Guermond, A. Salgado, Error analysis of a frational time-stepping technique for incompressible flows with variable density, SIAM J. Numer. Anal. 49 (3) (2011) 917–944. [18] F. Guillen-Gonzalez, J. Gutierrez-Santacreu, Unconditional stability and convergence of fully discrete schemes for 2D viscous fluids models with mass diffusion, Math. Comput. 77 (2008) 1495–1524. [19] P. Huang, J. Zhao, X. Feng, Highly efficient and local projecton-based stabilized finite element method for natural convection problem, Int. J. Heat Mass Transf. 83 (2015) 357–365. [20] P. Lions, Mathematical Topics in Fluid Mechanics, vol. 1, Oxford Lecture Series in Mathematics and Its Applications, vol. 3, The Clarendon Press Oxford University Press, New York, 1996. [21] Z. Luo, Theory Bases and Applications of Finite Element Methods, Science Press, Beijing, 2006 (in Chinese). [22] J. Pyo, J. Shen, Gauge-Uzawa methods for incompressible flows with variable density, J. Comput. Phys. 221 (2007) 181–197. [23] J. Pyo, Error estimates for the second order semi-discrete stabilized Gauge-Uzawa method for the Navier-Stokes equations, Int. J. Numer. Anal. Model. 10 (1) (2013) 24–41. [24] J. Shen, On error estimates of projection methods for the Navier-Stokes equations: first-order schemes, SIAM J. Numer. Anal. 29 (1992) 57–77. [25] H. Su, L. Qian, D. Gui, X. Feng, Second order fully discrete and divergence free conserving scheme for time-dependent conduction-convection equations, Int. Commun. Heat Mass Transf. 59 (2014) 120–129. [26] H. Su, X. Feng, Y. He, Defect-correction finite element method based on Crank-Nicolson extrapolation scheme for the transient conduction-convection problem with high Reynolds number, Int. Commun. Heat Mass Transf. 81 (2017) 229–249. [27] H. Su, X. Feng, Y. He, Second order fully discrete defect-correction scheme for nonstationary conduction-convection problem at high Reynolds number, Numer. Methods Partial Differ. Equ. 33 (3) (2017) 681–703. [28] W. Wang, J. Wu, X. Feng, A novel characteristic variational multiscale FEM for incompressible natural convection problem with variable density, Int. J. Numer. Methods Heat Fluid Flow 29 (2) (2019) 580–601. [29] M. Wheeler, A priori L2 error estimates for Galerkin approximations to parabolic partial differential equations, SIAM J. Numer. Anal. 10 (1973) 723–759. [30] J. Wu, P. Huang, X. Feng, D. Liu, An efficient two-step algorithm for steady-state natural convection problem, Int. J. Heat Mass Transf. 101 (2016) 387–398. [31] J. Wu, X. Feng, F. Liu, Pressure-correction projection FEM for time-dependent natural convection problem, Commun. Comput. Phys. 21 (2017) 1090–1117. [32] J. Wu, J. Shen, X. Feng, Unconditionally stable Guage-Uzawa finite element schemes for incompressible natural convection problems with variable density, J. Comput. Phys. 348 (2017) 776–789. [33] T. Zhang, X. Feng, J. Yuan, Implicit-explicit schemes of finite element method for the non-stationary thermal convection problems with temperaturedependent coefficients, Int. Commun. Heat Mass Transf. 76 (2016) 325–336.