# Long Time Stability Of A Linearly Extrapolated

## Transcript Of Long Time Stability Of A Linearly Extrapolated

INTERNATIONAL JOURNAL OF NUMERICAL ANALYSIS AND MODELING Volume 17, Number 1, Pages 24–41

c 2020 Institute for Scientiﬁc Computing and Information

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLENDED BDF SCHEME FOR MULTIPHYSICS FLOWS

AYTEKIN C¸ IBIK, FATMA G. EROGLU, AND SONGUL KAYA∗

Abstract. This paper investigates the long time stability behavior of multiphysics ﬂow problems, namely the Navier-Stokes equations, natural convection and double-diﬀusive convection equations with an extrapolated blended BDF time-stepping scheme. This scheme combines the two-step BDF and three-step BDF time stepping schemes. We prove unconditional long time stability theorems for each of these ﬂow systems. Various numerical tests are given for large time step sizes in long time intervals in order to support theoretical results.

Key words. Blended BDF, long time stability, Navier-Stokes, natural convection, doublediﬀusive.

1. Introduction

Most of the engineering and applied science problem involves the combination of some diﬀerent physical problems such as ﬂuid ﬂow, heat transfer, mass transfer, and electromagnetic eﬀect. These kinds of problems are mostly referred as multiphysics problems. From a mathematical point of view, these problems yield systems of coupled single physics equations. In our case, many important applications require the accurate solution of multiphysics coupling with Navier-Stokes equations. Since the simulation of Navier-Stokes equations has its own diﬃculties, the coupling among involved equations yield more complex problems. One possibility of improving numerical simulations is to develop algorithms which are reliable and robust. In addition, designed numerical schemes should capture the long time dynamics of the system in a right way. Thus, it is of practical interest to have an algorithm which is stable over the required long time intervals.

In recent years, considerable amount of eﬀort has been spent to understand the long time behavior of the numerical schemes for multiphysics problems. For such works, we refer to [3, 6, 19, 20, 22, 24]. In particular, for Navier-Stokes equations, the Crank-Nicholson in [11, 12, 19], the implicit Euler in [22], two-step Backward Diﬀerentiation (BDF2) in [2] and fractional step/projection methods in [18] are chosen as temporal discretization in order to show the long time stability. In this respect, an extrapolated two step BDF scheme for a velocity-vorticity form of Navier-Stokes equations has been investigated in [10]. The long time stability of partitioned methods for the fully evolutionary Stokes-Darcy problem in [13], for the time dependent MHD system in [14, 20] and for the double-diﬀusive convection in [21] were also established based on implicit-explicit and backward diﬀerentiation schemes.

This study concerns the behavior of the solutions of multiphysics problems for longer time simulations and time step restrictions on these solutions. In order to solve the systems numerically, a ﬁnite element in space discretization is employed along with a rather new time stepping approach so called an extrapolated

Received by the editors March 11, 2016, and in revised form, March 17, 2019. 2000 Mathematics Subject Classiﬁcation. 65M12,65M60, 65N30. *Corresponding Author.

24

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

25

blended Backward Diﬀerentiation (BLEBDF) temporal discretization. Such kind of a scheme combines BDF2 and a three-step BDF method in order to not only preserve A-stability and second order accuracy but also have a smaller constant in truncation error terms, [23]. Along with the mentioned time-stepping strategy, the three-step extrapolation is used in order to linearize the convective terms in the system of PDE’s. Thus, the solution of only one linear system of equations is encountered at each time step which reduces the compilation time and memory cost in simulations.

In [17], Ravindran considers the stability and convergence of a double diﬀusive convection system with the blended BDF scheme in short time intervals. The current study attempts to extend the works above to study the notion of long time stability by combining the BLEBDF idea and its eﬀects on several multiphysics ﬂows such as Navier-Stokes, natural convection and double-diﬀusive convection. In this work, we will provide the unconditonal long time L2 stability property of BLEBDF method for each of ﬂow systems, when they are discretized spatially with ﬁnite element method. To the best of authors’ knowledge, this is the ﬁrst study that investigates the long time stability of the ﬁnite element solutions of multiphysics ﬂows involving a blended BDF time-stepping approach along with a third order linear extrapolation idea.

The plan of the paper is as follows. In Section 2, we state the notations with some mathematical preliminaries. Section 3 is reserved for proving the unconditional long time stability of Navier-Stokes equations under the employment of extrapolated blended BDF temporal discretization along with numerical experiments. Similar results are established for natural convection and double-diﬀusive convection equations in Section 4 and Section 5, respectively. Finally, we state some conclusion remarks in the last section.

2. Notations and Preliminaries

Let Ω ⊂ Rd, d ∈ {2, 3}, be open, connected domain bounded by Lipschitz boundary ∂Ω. Throughout the paper standard notations for Sobolev spaces and their norms will be used, c.f. Adams [1]. The norm in (Hk(Ω))d is denoted by · k and the norm in Lebesgue spaces (Lp(Ω))d, 1 ≤ p < ∞, p = 2 by · Lp and p = ∞ by

· ∞ . The space L2(Ω) is equipped with the norm and inner product · and (·, ·), respectively, and for these we drop the subscripts. Vector-valued functions will be identiﬁed by bold face. The norm in dual space H−1 of H01(Ω) is denoted by · −1. The continuous velocity, pressure, temperature and concentration spaces are denoted by

X := (H10(Ω))d, Q := L20(Ω), W := H01(Ω), Ψ := H01(Ω),

and the divergence free space

V = {v ∈ X : (∇ · v, q) = 0, ∀q ∈ Q}.

We recall also the Poincar´e-Friedrichs inequality as

v ≤ CP ∇v , ∀v ∈ X.

For each multiphysics ﬂow problems, we consider a regular, conforming family Πh of triangulations of domain with maximum diameter h for spatial discretization. Assume Xh ⊂ X, Qh ⊂ Q, Wh ⊂ W and Ψh ⊂ Ψ be ﬁnite element spaces such that the spaces (Xh, Qh) satisfy the discrete inf-sup condition needed for stability of the discrete pressure, [5]. The discretely divergence free space for (Xh, Qh) pair

26

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

is given by

(1)

Vh = {vh ∈ Xh : (∇ · vh, qh) = 0, ∀qh ∈ Qh}.

The dual spaces of Vh, Wh and Ψh are given by Vh∗ , Wh∗ and Ψ∗h, and their norms

are denoted by · V∗ , · W ∗ and · Ψ∗ respectively. We also need the following

h

h

h

space in the analysis

(2) L∞(R+, Vh∗ ) := {f : Ωd × R+ → Rd, ∃ K < ∞, a.e. t > 0, f (t) Vh∗ < K}.

Similar spaces for Wh∗ and Ψ∗h will be used throughout the analysis. To simplify the analysis, we utilize the G-stability framework as in [8]. For third order backward

diﬀerentiation, the positive deﬁnite matrix G-matrix and the associated norm can

be obtained as

1

19

−12

3

G = −12 10 −3 ,

12 3 −3 1

W

2 G

=

(W ,

GW ).

W ∈ R3n.

The following relations are well-known and required for the analysis.

Lemma 2.1. G-norm and L2 norm are equivalent in the sense that there exist 0 < Cl < Cu positive constants such that

(3)

Cl

W

2 G

≤

W 2 ≤ Cu W 2G.

Proof. See reference [8].

Lemma 2.2. The following estimates hold:

u 2

1 u 2 − 5 v 2 − 5 w 2,

(4)

v ≥

w

3

12

12

G

u 2

19 u 2 + 10 v 2 + 1 w 2.

(5)

v ≤

w

12

12

12

G

Proof. Using the properties of inner product and norm, we have

u 2 v =

wG

(6)

u u

v,Gv

w

w

= 1 u 2− 5 v 2− 5 w 2+ u−v 2

3

12

12

+ 3 u + w 2 + 3 v − w 2.

12

12

Dropping last three non-negative terms yields (4). In a similar manner, deleting non-positive terms in (6) and using the triangle inequality gives (5).

Lemma 2.3. The symmetric matrix G ∈ R3n×3n which is given above satisﬁes the following equality ∀{whj }Nj=0 ⊂ L2(Ω):

53 whn+1 − 25 whn + whn−1 − 61 whn−2, whn+1

(7)

= Wn+1 2 − Wn 2 + 1 wn+1 − 3wn + 3wn−1 − wn−2 2,

G

G 12 h

h

h

h

where Wn+1 = [whn+1 whn whn−1] and Wn = [whn whn−1 whn−2] .

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

27

Proof. To prove this, we use a well-known identity (see [17]) to obtain

35 whn+1 − 52 whn + whn−1 − 61 whn−2, whn+1

=

whn+1 2 + 3whn+1 − whn 2 + 3whn+1 − 3whn + whn−1 2

−( whn 2 + 3whn − whn−1 2 + 3whn − 3whn−1 + whn−2 2) (8) + 112 whn+1 − 3whn + 3whn−1 − whn−2 2.

Using the properties of G-norm and inner product, one gets

whn+1 2 + 3whn+1 − whn 2 + 3whn+1 − 3whn + whn−1 2

(9) Obtaining

whn+1 whn+1

whn+1 2

=

whn

,

G

w

n h

= whn

= Wn+1 2G.

whn−1

whn−1

whn−1 G

Wn

2 G

similar with (9) and inserting them in (8) yields (7).

3. Long time stability of Navier-Stokes equations with BLEBDF

The goal of this section is to show that our scheme is unconditionally long-time stable. That is, the solutions remain bounded without any time step restriction. We ﬁrst study an extrapolated blended BDF method for discretizing the incompressible NSE:

ut − ν∆u + (u · ∇)u + ∇p = f in Ω,

∇ · u = 0 in Ω,

(10)

u = 0 on ∂Ω,

u(0, x) = u0 in Ω,

p dx = 0.

Ω

where u is the velocity ﬁeld, p is the ﬂuid pressure, ν is the kinematic viscosity and f denotes the body force. The variational formulation of (10) reads as follows: Find (u, p) ∈ (X, Q) satisfying

(11) (ut, v) + ν(∇u, ∇v) + b1(u, u, v) − ((p∇, ∇· u·,vq)) == (0f , v) ∀∀ vq ∈∈ QX,,

where

(12)

b1(u, v, w) := 1 (((u · ∇)v, w) − ((u · ∇)w, v)) ,

2

represents the skew-symmetric form of the convective term. Note that the convec-

tive term has the well known property

(13)

b1(u, v, v) = 0

for all u, v ∈ X, which simpliﬁes the analysis. The studied time discretization method uses diﬀerent time discretizations for

diﬀerent terms. The scheme discretizes in time via a BDF2 and a BDF3 scheme whereas the nonlinear terms are treated via a third-order extrapolation formula. One consequence of lagging the nonlinear term to previous time levels is to avoid solving nonlinear equations. A class of this type of blending BDF schemes which are also known as a class of optimized second-order BDF schemes are proposed in

28

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

[15, 23]. We now consider one of the special case of this family proposed in [23],

with an error constant half as large as BDF2 scheme.

Based on the weak formulation (10), the fully discrete approximation of it reads as follows. Given f , u0h = u−h 1 = u−h 2 = u0, for any time step ∆t > 0, ﬁnd (unh+1, pnh+1) ∈ (Xh, Qh) such that for n ≥ 0

53 unh+1 − 52 unh ∆+tuhn−1 − 16 uhn−2 , vh + ν(∇unh+1, ∇vh)

(14) (15)

+b1(3unh − 3uhn−1 + uhn−2, unh+1, vh) − (pnh+1, ∇ · vh) = (f n+1, vh), (∇ · unh+1, qh) = 0

for all (vh, qh) ∈ (Xh, Qh). We now analyze the long time stability over 0 ≤ tn < ∞ and show that the

scheme is uniformly bounded in time in L2 norm. Our proof is based on a G-norm

and the associated estimation (7).

Theorem 3.1. (Unconditional long time stability of NSE in L2) Assume f ∈

L∞(R+, Vh∗), then the approximation (14)-(15) is long time stable in the following sense: for any ∆t > 0,

(16)

Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2 ≤ (1 + α)−(n+1)( U0 2

G4

h

16

h

G

+ ν∆t ∇u0 2 + ν∆t ∇u−1 2) + max{ 8Cp2 , 2ν−1∆t } f 2

.

4

h

16

h

Clν2 3

L∞(R+;Vh∗ )

where U0 = [u0 u−1 u−2] and α = min{ Clν∆t , 3 }.

hh

h

16Cp2 4

Remark 3.1. Theorem implies that the result holds for any large n, since the constants are independent of n.

Proof. Setting vh = unh+1 in (14) and qh = pnh+1 in (15) the BLEBDF scheme and using the skew-symmetry property (13), we obtain

(17)

35 unh+1

−

25 unh

+

uhn−1

−

1 6

uhn−2

,

un+1

∆t

h

+ ν(∇unh+1, ∇unh+1) = (f n+1, unh+1).

From (7), we get

Un+1

2 G

−

(18)

Un 2 + 1 G 12

unh+1 − 3unh + 3uhn−1 − uhn−2 2 + ν∆t ∇unh+1 2 = ∆t(f n+1, unh+1).

where Un+1 = [unh+1 unh uhn−1] and Un = [unh uhn−1 uhn−2] . Applying the Cauchy-Schwarz and Young’s inequalities for the right hand side of (18) gives

Un+1

2 G

−

(19)

Un 2 + 1 G 12

unh+1 − 3unh + 3uhn−1 − uhn−2 2 + ν∆2 t ≤ ν−1∆t 2

∇unh+1 2

f n+1

2 V

∗

.

h

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

29

Adding ν∆4 t ∇unh 2 + ν1∆6t ∇uhn−1 2 to the both sides and discarding the third positive term in the left hand side of (19) yield

Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2 + ν∆t

G4

h

16

h

16

∇unh+1 2

(20)

+ ∇unh 2 + ∇uhn−1 2 + 3ν1∆6 t ∇unh+1 2 + ν∆8 t ∇unh 2

≤ U 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2 + ν−1∆t f n+1 2 .

nG

4

h

16

h

2

V∗

h

The terms in the left hand side of (20) can be rearranged by applying the Poincar´e-

Friedrichs and the equivalent norm property (3):

ν∆t 16

∇unh+1 2 + ∇unh 2 + ∇uhn−1 2 + 3ν1∆6 t ∇unh+1 2

+ ν∆8 t ∇unh 2

≥ Clν∆t Un+1 2 + 3ν∆t ∇un+1 2 + ν∆t ∇un 2

16Cp2

G 16

h

8

h

(21)

≥ α( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2),

G4

h

16

h

where α = min{ Clν∆t , 3 }. Inserting (21) in (20) and using induction along with 16Cp2 4

f i 2Vh∗ ≤ f 2L∞(R+;Vh∗ ), ∀i = 1, . . . , n + 1 gives

( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ (1 + α)−(n+1)( U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2)

G4

h

16

h

+ (1 + α)−21ν−1∆t f 2L∞(R+;Vh∗) (1 + α)−n

(22)

+(1 + α)−(n−1) + · · · + 1 .

Since | 1 | < 1, then we get 1+α

(23)

( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ (1 + α)−(n+1)( U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2)

G4

h

16

h

+ max{ 8Cp2 , 2ν−1∆t } f 2

Clν2 3

L∞(R+;Vh∗ )

which is the required result.

3.1. Numerical experiment for Navier-Stokes equations. Throughout this paper, we carry out tests for each ﬂow problem separately. We expose the evolution of the L2 norm of solution variables in long time intervals for each variable according to related problems. After each numerical test, we provide a table showing the CPU-times of relevant simulations according to varying ∆t and compare these CPU-times with classical BDF2 methods’ CPU-times in order to reveal the computational advantage of the method. We choose the inf-sup stable Scott-Vogelious

30

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

ﬁnite element pair for velocity-pressure couple, which is known to be discretely divergence free. We refer [4], for the details of this selection. Throughout all our simulations, we use public license ﬁnite element software package FreeFem++, [9].

We now perform a numerical test in order to verify theoretical long time stability result obtained in Theorem 3.1. We set Ω = (0, 1)2 with a coarse mesh resolution of 16 × 16. We calculate the approximate solution in time interval [0, 400] and calculate the L2 norms for diﬀerent ∆t and varying ν. We pick the initial velocity and the forcing term for this test problem as :

(24) u = csoins((ππxx)) scions((ππyy)) , f = 2yx2yccooss((xxyy22))++scions((xx))scino(sy(y)) .

In Figure 1, we present the evolution of un+1 in long time interval [0, 400]. We give the stability results for two large time step instances for each case. As could be easily deduced from the ﬁgure, the solutions we obtain from the proposed scheme are long time stable even by using large time step sizes. Although slight deviations seem to occur due to small viscosity for the case ν = 0.001, still the solution is within an interval of stability. Notice that, these solutions are obtained for a very coarse mesh resolution of 16 × 16. Hence, small oscillations inside a narrow interval does not degrade the stability of the solution.

|| un+1|| h

0.02 0.018 0.016 0.014 0.012

0.01 0.008 0.006 0.004 0.002

t =1 t=0.25

0

50 100 150 200 250 300 350 400

Time

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

|| un+1|| h

|| un+1|| h

0.55 0.5

0.45 0.4

0.35 0.3

0.25 0.2

0.15 0.1

0.05 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

1

0.8

0.6

0.4

0.2

0 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

Figure 1. Evolution of the L2 norm of the velocity solution for varying ν. ν = 1, ν = 0.004 (upper left to right) and ν = 0.002, ν = 0.001 (lower left to right).

|| un+1|| h

We also compare the CPU times of the scheme (14)-(15) and run the same problem in a shorter time interval [0, 150]. As could be observed from Table 1, the BLEBDF scheme has an advantage in terms of simulation CPU times when compared with a classical BDF2 scheme. When ∆t becomes smaller, the diﬀerence

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

31

between the CPU times are increasing. This shows the promise of the method, especially when smaller time step sizes are used.

Table 1. Comparison of the CPU-times (seconds) of classical BDF2 scheme and BLEBDF for the Navier-Stokes equations with ν = 0.001 on a time interval [0, 150].

∆t BDF2 BLEBDF

1 1.37

0.87

0.1 12

8

0.01 131.5 78.7

4. Long time stability of natural convection equations with BLEBDF

In this section, we prove that the BLEBDF scheme for natural convection equations is also long time stable. The unsteady natural convection system in Ω with partitioned boundary ∂Ω = ΓT ∪ ΓB with ΓT ∩ ΓB = ∅, the so called Boussinesq equations are given by

ut − ν∆u + (u · ∇)u + ∇p = Ri T e2 + f in Ω,

∇·u = 0

in Ω,

(25)

Tt − ∇ · (κ∇T ) + (u · ∇)T = γ

u(0, x) = u0, T (0, x) = T0

in Ω, in Ω,

u = 0 on ∂Ω, T = 0 on ΓT , ∂T = 0 ∂n

on ΓB,

where u, p, T are the ﬂuid velocity, the pressure and the temperature, respectively. The parameters in (25) are the kinematic viscosity ν, the thermal conductivity parameter κ > 0 and the Richardson number Ri and the unit vector is given by e2. The prescribed body forces are f and γ. The initial velocity and temperature are u0 and T0, respectively.

By using similar notations as in Section 3, given f , γ, u0h = u−h 1 = u−h 2 = u0 and Th0 = Th−1 = Th−2 = T0 ﬁnd (unh+1, pnh+1, Thn+1) ∈ (Xh, Qh, Wh) for n ≥ 1 satisfying

(26) (27)

35 unh+1 − 52 unh ∆+tuhn−1 − 16 uhn−2 , vh + ν(∇unh+1, ∇vh)

+b1(u∗, unh+1, vh) − (ph, ∇ · unh+1) = Ri (T ∗e2, vh) + (f n+1, vh), (qh, ∇ · vhn+1) = 0

35 Thn+1 − 25 Thn ∆+tThn−1 − 16 Thn−2 , Sh + κ(∇Thn+1, ∇Sh)

(28)

+b2(u∗, Thn+1, Sh) = (γn+1, Sh),

for all (vh, qh, Sh) ∈ (Xh, Qh, Wh) where u∗ = 3unh − 3uhn−1 + uhn−2 and T ∗ = 3Thn − 3Thn−1 + Thn−2. Herein the related skew-symmetric form is given by

(29)

b2(u, T, S) := 1 (((u · ∇)T, S) − ((u · ∇)S, T ))

2

for all u ∈ X, T, S ∈ W.

32

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

Theorem 4.1. (Unconditional long time stability of natural convection in L2) As-

sume f ∈ L∞(R+, Vh∗ ) and γ ∈ L∞(R+, Wh∗), then the approximation (26)-(28) is long time stable in the following sense: for any ∆t > 0,

Un+1 2 + Tn+1 2 + ν∆t ∇un+1 2 + κ∆t ∇T n+1 2

G

G4

h

4

h

≤ (1 + α)−(n+1) U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2

G4

h

16

h

+(Kα + (1 + β)−1) (1 + β)−n( T0 2 + κ∆t ∇T 0 2 + κ∆t ∇T −1 2)

G4

h

16

h

+(K + 1) max{ 8Cp2 , 2κ−1∆t } γ 2

α

Clκ2 3

L∞(R+ ,Wh∗ )

(30)

16Cp2 4ν−1∆t

+ max{ Clν2 ,

3

} f L∞(R+,V∗ ) h

where K = 49CuCp2ν−1Ri2∆t , U = [u0 u−1 u−2], T = [T 0 T −1 T −2],

α

α

0

hh

h

0

hh

h

α = min{ Clν∆t , 3 }, β = min{ Clκ∆t , 3 },

16Cp2 4

16Cp2 4

Proof. The proof starts with a stability bound for temperature. Letting Sh = Thn+1

in

(28)

and

using

the

skew

symmetry

property

b

2

(u∗

,

T

n+1 h

,

T

n+1 h

)

=

0.

Along

with

(7) and multiplying with ∆t yields

Tn+1

2 G

−

(31)

Tn 2 + 1 G 12

Thn+1 − 3Thn + 3Thn−1 − Thn−2 2 + κ∆t ∇Thn+1 2 = ∆t(γn+1, Thn+1).

where Tn+1 = [Thn+1 Thn Thn−1] and Tn = [Thn Thn−1 Thn−2] . Applying the Cauchy-Schwarz and Young’s inequalities lead to

Tn+1 2 − Tn 2 + 1 T n+1 − 3T n + 3T n−1 − T n−2 2 + κ∆t ∇T n+1 2

G

G 12 h

h

h

h

2

h

(32)

≤ κ−1∆t γn+1 2 .

2

Wh∗

Adding κ∆4 t ∇Thn 2 + κ1∆6t ∇Thn−1 2 to the both sides and dropping the nonnegative terms produce

(33)

( Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2) + κ∆t

G4

h

16

h

16

∇Thn+1 2

+ ∇Thn 2 + ∇Thn−1 2 + 3κ1∆6 t ∇Thn+1 2 + κ∆8 t ∇Thn 2

≤ T 2 + κ∆t ∇T n 2 + κ∆t ∇T n−1 2 + κ−1∆t γn+1 2 .

nG

4

h

16

h

2

Wh∗

The estimation of (33) follows closely that of (21). Again, the last ﬁve terms of the left hand side of (33) can be rearranged by the applying the Poincar´e-Friedrichs and the equivalent norm property (3) as

(34)

κ1∆6t ( ∇Thn+1 2 + ∇Thn 2 + ∇Thn−1 2)

+ 3κ1∆6 t ∇Thn+1 2 + κ∆8 t ∇Thn 2

≥ β( Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2),

G4

h

16

h

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

33

where β = min{ Clκ∆t , 3 }. Arguing as before and inserting (34) in (33) and using 16Cp2 4

induction leads to

(35)

Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2

G4

h

16

h

≤ (1 + β)−(n+1)( T0 2 + κ∆t ∇T 0 2 + κ∆t ∇T −1 2)

G4

h

16

h

+ κ−21β∆t γ 2L∞(R+,Wh∗),

which is the result for long the time stability for the temperature. Letting vh = unh+1 in (26), and qh = pnh+1 in (27), one obtains

35 unh+1

−

25 unh

+

uhn−1

−

1 6

uhn−2

,

un+1

∆t

h

+ ν ∇unh+1 2

(36)

= Ri (T ∗e2, unh+1) + (f n+1, unh+1).

Repeating the arguments used for obtaining (19) yields

(37)

Un+1 2 − Un 2 + 1 un+1 − 3un + 3un−1 − un−2 2

G

G 12 h

h

h

h

+ ν∆2 t ∇unh+1 2

≤ Cp2ν−1Ri2∆t T ∗ 2 e2 2 + ν−1∆t f n+1 Vh∗ .

Note that e2 is a unit vec√tor, i.e., e2 = 1 and using the deﬁnition of T ∗ and (3), we get T ∗ ≤ 7 Tn ≤ 7 Cu Tn G. Then, one has

(38)

Un+1

2 G

−

Un 2 + 1 un+1 − 3un + 3un−1 − un−2 2 + ν∆t

G 12 h

h

h

h

2

≤ 49CuCp2ν−1Ri2∆t Tn 2G + ν−1∆t f n+1 Vh∗ .

∇unh+1 2

Arguing as in (20), one gets the following estimation for (38)

(39)

(1 + α)( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ ( Un 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2)

G4

h

16

h

+49CuCp2ν−1Ri2∆t Tn 2G + ν−1∆t f n+1 Vh∗

where α = min{ Clν∆t , 3 }. Putting n instead of (n + 1) in (35) and inserting it in 16Cp2 4

(39) yields

(40)

(1 + α)( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ ( Un 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2)

G4

h

16

h

+49CuC2ν−1Ri2∆t (1 + β)−n( T0 2 + κ∆t ∇T 0 2

p

G4

h

+ κ1∆6t ∇Th−1 2) + κ−21β∆t γ 2L∞(R+,Wh∗)

+ ν−1∆t f n+1 V∗ . h

Applying induction and adding (35) to (40) completes the proof of the theorem.

c 2020 Institute for Scientiﬁc Computing and Information

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLENDED BDF SCHEME FOR MULTIPHYSICS FLOWS

AYTEKIN C¸ IBIK, FATMA G. EROGLU, AND SONGUL KAYA∗

Abstract. This paper investigates the long time stability behavior of multiphysics ﬂow problems, namely the Navier-Stokes equations, natural convection and double-diﬀusive convection equations with an extrapolated blended BDF time-stepping scheme. This scheme combines the two-step BDF and three-step BDF time stepping schemes. We prove unconditional long time stability theorems for each of these ﬂow systems. Various numerical tests are given for large time step sizes in long time intervals in order to support theoretical results.

Key words. Blended BDF, long time stability, Navier-Stokes, natural convection, doublediﬀusive.

1. Introduction

Most of the engineering and applied science problem involves the combination of some diﬀerent physical problems such as ﬂuid ﬂow, heat transfer, mass transfer, and electromagnetic eﬀect. These kinds of problems are mostly referred as multiphysics problems. From a mathematical point of view, these problems yield systems of coupled single physics equations. In our case, many important applications require the accurate solution of multiphysics coupling with Navier-Stokes equations. Since the simulation of Navier-Stokes equations has its own diﬃculties, the coupling among involved equations yield more complex problems. One possibility of improving numerical simulations is to develop algorithms which are reliable and robust. In addition, designed numerical schemes should capture the long time dynamics of the system in a right way. Thus, it is of practical interest to have an algorithm which is stable over the required long time intervals.

In recent years, considerable amount of eﬀort has been spent to understand the long time behavior of the numerical schemes for multiphysics problems. For such works, we refer to [3, 6, 19, 20, 22, 24]. In particular, for Navier-Stokes equations, the Crank-Nicholson in [11, 12, 19], the implicit Euler in [22], two-step Backward Diﬀerentiation (BDF2) in [2] and fractional step/projection methods in [18] are chosen as temporal discretization in order to show the long time stability. In this respect, an extrapolated two step BDF scheme for a velocity-vorticity form of Navier-Stokes equations has been investigated in [10]. The long time stability of partitioned methods for the fully evolutionary Stokes-Darcy problem in [13], for the time dependent MHD system in [14, 20] and for the double-diﬀusive convection in [21] were also established based on implicit-explicit and backward diﬀerentiation schemes.

This study concerns the behavior of the solutions of multiphysics problems for longer time simulations and time step restrictions on these solutions. In order to solve the systems numerically, a ﬁnite element in space discretization is employed along with a rather new time stepping approach so called an extrapolated

Received by the editors March 11, 2016, and in revised form, March 17, 2019. 2000 Mathematics Subject Classiﬁcation. 65M12,65M60, 65N30. *Corresponding Author.

24

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

25

blended Backward Diﬀerentiation (BLEBDF) temporal discretization. Such kind of a scheme combines BDF2 and a three-step BDF method in order to not only preserve A-stability and second order accuracy but also have a smaller constant in truncation error terms, [23]. Along with the mentioned time-stepping strategy, the three-step extrapolation is used in order to linearize the convective terms in the system of PDE’s. Thus, the solution of only one linear system of equations is encountered at each time step which reduces the compilation time and memory cost in simulations.

In [17], Ravindran considers the stability and convergence of a double diﬀusive convection system with the blended BDF scheme in short time intervals. The current study attempts to extend the works above to study the notion of long time stability by combining the BLEBDF idea and its eﬀects on several multiphysics ﬂows such as Navier-Stokes, natural convection and double-diﬀusive convection. In this work, we will provide the unconditonal long time L2 stability property of BLEBDF method for each of ﬂow systems, when they are discretized spatially with ﬁnite element method. To the best of authors’ knowledge, this is the ﬁrst study that investigates the long time stability of the ﬁnite element solutions of multiphysics ﬂows involving a blended BDF time-stepping approach along with a third order linear extrapolation idea.

The plan of the paper is as follows. In Section 2, we state the notations with some mathematical preliminaries. Section 3 is reserved for proving the unconditional long time stability of Navier-Stokes equations under the employment of extrapolated blended BDF temporal discretization along with numerical experiments. Similar results are established for natural convection and double-diﬀusive convection equations in Section 4 and Section 5, respectively. Finally, we state some conclusion remarks in the last section.

2. Notations and Preliminaries

Let Ω ⊂ Rd, d ∈ {2, 3}, be open, connected domain bounded by Lipschitz boundary ∂Ω. Throughout the paper standard notations for Sobolev spaces and their norms will be used, c.f. Adams [1]. The norm in (Hk(Ω))d is denoted by · k and the norm in Lebesgue spaces (Lp(Ω))d, 1 ≤ p < ∞, p = 2 by · Lp and p = ∞ by

· ∞ . The space L2(Ω) is equipped with the norm and inner product · and (·, ·), respectively, and for these we drop the subscripts. Vector-valued functions will be identiﬁed by bold face. The norm in dual space H−1 of H01(Ω) is denoted by · −1. The continuous velocity, pressure, temperature and concentration spaces are denoted by

X := (H10(Ω))d, Q := L20(Ω), W := H01(Ω), Ψ := H01(Ω),

and the divergence free space

V = {v ∈ X : (∇ · v, q) = 0, ∀q ∈ Q}.

We recall also the Poincar´e-Friedrichs inequality as

v ≤ CP ∇v , ∀v ∈ X.

For each multiphysics ﬂow problems, we consider a regular, conforming family Πh of triangulations of domain with maximum diameter h for spatial discretization. Assume Xh ⊂ X, Qh ⊂ Q, Wh ⊂ W and Ψh ⊂ Ψ be ﬁnite element spaces such that the spaces (Xh, Qh) satisfy the discrete inf-sup condition needed for stability of the discrete pressure, [5]. The discretely divergence free space for (Xh, Qh) pair

26

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

is given by

(1)

Vh = {vh ∈ Xh : (∇ · vh, qh) = 0, ∀qh ∈ Qh}.

The dual spaces of Vh, Wh and Ψh are given by Vh∗ , Wh∗ and Ψ∗h, and their norms

are denoted by · V∗ , · W ∗ and · Ψ∗ respectively. We also need the following

h

h

h

space in the analysis

(2) L∞(R+, Vh∗ ) := {f : Ωd × R+ → Rd, ∃ K < ∞, a.e. t > 0, f (t) Vh∗ < K}.

Similar spaces for Wh∗ and Ψ∗h will be used throughout the analysis. To simplify the analysis, we utilize the G-stability framework as in [8]. For third order backward

diﬀerentiation, the positive deﬁnite matrix G-matrix and the associated norm can

be obtained as

1

19

−12

3

G = −12 10 −3 ,

12 3 −3 1

W

2 G

=

(W ,

GW ).

W ∈ R3n.

The following relations are well-known and required for the analysis.

Lemma 2.1. G-norm and L2 norm are equivalent in the sense that there exist 0 < Cl < Cu positive constants such that

(3)

Cl

W

2 G

≤

W 2 ≤ Cu W 2G.

Proof. See reference [8].

Lemma 2.2. The following estimates hold:

u 2

1 u 2 − 5 v 2 − 5 w 2,

(4)

v ≥

w

3

12

12

G

u 2

19 u 2 + 10 v 2 + 1 w 2.

(5)

v ≤

w

12

12

12

G

Proof. Using the properties of inner product and norm, we have

u 2 v =

wG

(6)

u u

v,Gv

w

w

= 1 u 2− 5 v 2− 5 w 2+ u−v 2

3

12

12

+ 3 u + w 2 + 3 v − w 2.

12

12

Dropping last three non-negative terms yields (4). In a similar manner, deleting non-positive terms in (6) and using the triangle inequality gives (5).

Lemma 2.3. The symmetric matrix G ∈ R3n×3n which is given above satisﬁes the following equality ∀{whj }Nj=0 ⊂ L2(Ω):

53 whn+1 − 25 whn + whn−1 − 61 whn−2, whn+1

(7)

= Wn+1 2 − Wn 2 + 1 wn+1 − 3wn + 3wn−1 − wn−2 2,

G

G 12 h

h

h

h

where Wn+1 = [whn+1 whn whn−1] and Wn = [whn whn−1 whn−2] .

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

27

Proof. To prove this, we use a well-known identity (see [17]) to obtain

35 whn+1 − 52 whn + whn−1 − 61 whn−2, whn+1

=

whn+1 2 + 3whn+1 − whn 2 + 3whn+1 − 3whn + whn−1 2

−( whn 2 + 3whn − whn−1 2 + 3whn − 3whn−1 + whn−2 2) (8) + 112 whn+1 − 3whn + 3whn−1 − whn−2 2.

Using the properties of G-norm and inner product, one gets

whn+1 2 + 3whn+1 − whn 2 + 3whn+1 − 3whn + whn−1 2

(9) Obtaining

whn+1 whn+1

whn+1 2

=

whn

,

G

w

n h

= whn

= Wn+1 2G.

whn−1

whn−1

whn−1 G

Wn

2 G

similar with (9) and inserting them in (8) yields (7).

3. Long time stability of Navier-Stokes equations with BLEBDF

The goal of this section is to show that our scheme is unconditionally long-time stable. That is, the solutions remain bounded without any time step restriction. We ﬁrst study an extrapolated blended BDF method for discretizing the incompressible NSE:

ut − ν∆u + (u · ∇)u + ∇p = f in Ω,

∇ · u = 0 in Ω,

(10)

u = 0 on ∂Ω,

u(0, x) = u0 in Ω,

p dx = 0.

Ω

where u is the velocity ﬁeld, p is the ﬂuid pressure, ν is the kinematic viscosity and f denotes the body force. The variational formulation of (10) reads as follows: Find (u, p) ∈ (X, Q) satisfying

(11) (ut, v) + ν(∇u, ∇v) + b1(u, u, v) − ((p∇, ∇· u·,vq)) == (0f , v) ∀∀ vq ∈∈ QX,,

where

(12)

b1(u, v, w) := 1 (((u · ∇)v, w) − ((u · ∇)w, v)) ,

2

represents the skew-symmetric form of the convective term. Note that the convec-

tive term has the well known property

(13)

b1(u, v, v) = 0

for all u, v ∈ X, which simpliﬁes the analysis. The studied time discretization method uses diﬀerent time discretizations for

diﬀerent terms. The scheme discretizes in time via a BDF2 and a BDF3 scheme whereas the nonlinear terms are treated via a third-order extrapolation formula. One consequence of lagging the nonlinear term to previous time levels is to avoid solving nonlinear equations. A class of this type of blending BDF schemes which are also known as a class of optimized second-order BDF schemes are proposed in

28

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

[15, 23]. We now consider one of the special case of this family proposed in [23],

with an error constant half as large as BDF2 scheme.

Based on the weak formulation (10), the fully discrete approximation of it reads as follows. Given f , u0h = u−h 1 = u−h 2 = u0, for any time step ∆t > 0, ﬁnd (unh+1, pnh+1) ∈ (Xh, Qh) such that for n ≥ 0

53 unh+1 − 52 unh ∆+tuhn−1 − 16 uhn−2 , vh + ν(∇unh+1, ∇vh)

(14) (15)

+b1(3unh − 3uhn−1 + uhn−2, unh+1, vh) − (pnh+1, ∇ · vh) = (f n+1, vh), (∇ · unh+1, qh) = 0

for all (vh, qh) ∈ (Xh, Qh). We now analyze the long time stability over 0 ≤ tn < ∞ and show that the

scheme is uniformly bounded in time in L2 norm. Our proof is based on a G-norm

and the associated estimation (7).

Theorem 3.1. (Unconditional long time stability of NSE in L2) Assume f ∈

L∞(R+, Vh∗), then the approximation (14)-(15) is long time stable in the following sense: for any ∆t > 0,

(16)

Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2 ≤ (1 + α)−(n+1)( U0 2

G4

h

16

h

G

+ ν∆t ∇u0 2 + ν∆t ∇u−1 2) + max{ 8Cp2 , 2ν−1∆t } f 2

.

4

h

16

h

Clν2 3

L∞(R+;Vh∗ )

where U0 = [u0 u−1 u−2] and α = min{ Clν∆t , 3 }.

hh

h

16Cp2 4

Remark 3.1. Theorem implies that the result holds for any large n, since the constants are independent of n.

Proof. Setting vh = unh+1 in (14) and qh = pnh+1 in (15) the BLEBDF scheme and using the skew-symmetry property (13), we obtain

(17)

35 unh+1

−

25 unh

+

uhn−1

−

1 6

uhn−2

,

un+1

∆t

h

+ ν(∇unh+1, ∇unh+1) = (f n+1, unh+1).

From (7), we get

Un+1

2 G

−

(18)

Un 2 + 1 G 12

unh+1 − 3unh + 3uhn−1 − uhn−2 2 + ν∆t ∇unh+1 2 = ∆t(f n+1, unh+1).

where Un+1 = [unh+1 unh uhn−1] and Un = [unh uhn−1 uhn−2] . Applying the Cauchy-Schwarz and Young’s inequalities for the right hand side of (18) gives

Un+1

2 G

−

(19)

Un 2 + 1 G 12

unh+1 − 3unh + 3uhn−1 − uhn−2 2 + ν∆2 t ≤ ν−1∆t 2

∇unh+1 2

f n+1

2 V

∗

.

h

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

29

Adding ν∆4 t ∇unh 2 + ν1∆6t ∇uhn−1 2 to the both sides and discarding the third positive term in the left hand side of (19) yield

Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2 + ν∆t

G4

h

16

h

16

∇unh+1 2

(20)

+ ∇unh 2 + ∇uhn−1 2 + 3ν1∆6 t ∇unh+1 2 + ν∆8 t ∇unh 2

≤ U 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2 + ν−1∆t f n+1 2 .

nG

4

h

16

h

2

V∗

h

The terms in the left hand side of (20) can be rearranged by applying the Poincar´e-

Friedrichs and the equivalent norm property (3):

ν∆t 16

∇unh+1 2 + ∇unh 2 + ∇uhn−1 2 + 3ν1∆6 t ∇unh+1 2

+ ν∆8 t ∇unh 2

≥ Clν∆t Un+1 2 + 3ν∆t ∇un+1 2 + ν∆t ∇un 2

16Cp2

G 16

h

8

h

(21)

≥ α( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2),

G4

h

16

h

where α = min{ Clν∆t , 3 }. Inserting (21) in (20) and using induction along with 16Cp2 4

f i 2Vh∗ ≤ f 2L∞(R+;Vh∗ ), ∀i = 1, . . . , n + 1 gives

( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ (1 + α)−(n+1)( U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2)

G4

h

16

h

+ (1 + α)−21ν−1∆t f 2L∞(R+;Vh∗) (1 + α)−n

(22)

+(1 + α)−(n−1) + · · · + 1 .

Since | 1 | < 1, then we get 1+α

(23)

( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ (1 + α)−(n+1)( U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2)

G4

h

16

h

+ max{ 8Cp2 , 2ν−1∆t } f 2

Clν2 3

L∞(R+;Vh∗ )

which is the required result.

3.1. Numerical experiment for Navier-Stokes equations. Throughout this paper, we carry out tests for each ﬂow problem separately. We expose the evolution of the L2 norm of solution variables in long time intervals for each variable according to related problems. After each numerical test, we provide a table showing the CPU-times of relevant simulations according to varying ∆t and compare these CPU-times with classical BDF2 methods’ CPU-times in order to reveal the computational advantage of the method. We choose the inf-sup stable Scott-Vogelious

30

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

ﬁnite element pair for velocity-pressure couple, which is known to be discretely divergence free. We refer [4], for the details of this selection. Throughout all our simulations, we use public license ﬁnite element software package FreeFem++, [9].

We now perform a numerical test in order to verify theoretical long time stability result obtained in Theorem 3.1. We set Ω = (0, 1)2 with a coarse mesh resolution of 16 × 16. We calculate the approximate solution in time interval [0, 400] and calculate the L2 norms for diﬀerent ∆t and varying ν. We pick the initial velocity and the forcing term for this test problem as :

(24) u = csoins((ππxx)) scions((ππyy)) , f = 2yx2yccooss((xxyy22))++scions((xx))scino(sy(y)) .

In Figure 1, we present the evolution of un+1 in long time interval [0, 400]. We give the stability results for two large time step instances for each case. As could be easily deduced from the ﬁgure, the solutions we obtain from the proposed scheme are long time stable even by using large time step sizes. Although slight deviations seem to occur due to small viscosity for the case ν = 0.001, still the solution is within an interval of stability. Notice that, these solutions are obtained for a very coarse mesh resolution of 16 × 16. Hence, small oscillations inside a narrow interval does not degrade the stability of the solution.

|| un+1|| h

0.02 0.018 0.016 0.014 0.012

0.01 0.008 0.006 0.004 0.002

t =1 t=0.25

0

50 100 150 200 250 300 350 400

Time

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

|| un+1|| h

|| un+1|| h

0.55 0.5

0.45 0.4

0.35 0.3

0.25 0.2

0.15 0.1

0.05 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

1

0.8

0.6

0.4

0.2

0 0

t=1 t=0.25

50 100 150 200 250 300 350 400 Time

Figure 1. Evolution of the L2 norm of the velocity solution for varying ν. ν = 1, ν = 0.004 (upper left to right) and ν = 0.002, ν = 0.001 (lower left to right).

|| un+1|| h

We also compare the CPU times of the scheme (14)-(15) and run the same problem in a shorter time interval [0, 150]. As could be observed from Table 1, the BLEBDF scheme has an advantage in terms of simulation CPU times when compared with a classical BDF2 scheme. When ∆t becomes smaller, the diﬀerence

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

31

between the CPU times are increasing. This shows the promise of the method, especially when smaller time step sizes are used.

Table 1. Comparison of the CPU-times (seconds) of classical BDF2 scheme and BLEBDF for the Navier-Stokes equations with ν = 0.001 on a time interval [0, 150].

∆t BDF2 BLEBDF

1 1.37

0.87

0.1 12

8

0.01 131.5 78.7

4. Long time stability of natural convection equations with BLEBDF

In this section, we prove that the BLEBDF scheme for natural convection equations is also long time stable. The unsteady natural convection system in Ω with partitioned boundary ∂Ω = ΓT ∪ ΓB with ΓT ∩ ΓB = ∅, the so called Boussinesq equations are given by

ut − ν∆u + (u · ∇)u + ∇p = Ri T e2 + f in Ω,

∇·u = 0

in Ω,

(25)

Tt − ∇ · (κ∇T ) + (u · ∇)T = γ

u(0, x) = u0, T (0, x) = T0

in Ω, in Ω,

u = 0 on ∂Ω, T = 0 on ΓT , ∂T = 0 ∂n

on ΓB,

where u, p, T are the ﬂuid velocity, the pressure and the temperature, respectively. The parameters in (25) are the kinematic viscosity ν, the thermal conductivity parameter κ > 0 and the Richardson number Ri and the unit vector is given by e2. The prescribed body forces are f and γ. The initial velocity and temperature are u0 and T0, respectively.

By using similar notations as in Section 3, given f , γ, u0h = u−h 1 = u−h 2 = u0 and Th0 = Th−1 = Th−2 = T0 ﬁnd (unh+1, pnh+1, Thn+1) ∈ (Xh, Qh, Wh) for n ≥ 1 satisfying

(26) (27)

35 unh+1 − 52 unh ∆+tuhn−1 − 16 uhn−2 , vh + ν(∇unh+1, ∇vh)

+b1(u∗, unh+1, vh) − (ph, ∇ · unh+1) = Ri (T ∗e2, vh) + (f n+1, vh), (qh, ∇ · vhn+1) = 0

35 Thn+1 − 25 Thn ∆+tThn−1 − 16 Thn−2 , Sh + κ(∇Thn+1, ∇Sh)

(28)

+b2(u∗, Thn+1, Sh) = (γn+1, Sh),

for all (vh, qh, Sh) ∈ (Xh, Qh, Wh) where u∗ = 3unh − 3uhn−1 + uhn−2 and T ∗ = 3Thn − 3Thn−1 + Thn−2. Herein the related skew-symmetric form is given by

(29)

b2(u, T, S) := 1 (((u · ∇)T, S) − ((u · ∇)S, T ))

2

for all u ∈ X, T, S ∈ W.

32

A. C¸ IBIK, F. G. EROGLU, AND S. KAYA

Theorem 4.1. (Unconditional long time stability of natural convection in L2) As-

sume f ∈ L∞(R+, Vh∗ ) and γ ∈ L∞(R+, Wh∗), then the approximation (26)-(28) is long time stable in the following sense: for any ∆t > 0,

Un+1 2 + Tn+1 2 + ν∆t ∇un+1 2 + κ∆t ∇T n+1 2

G

G4

h

4

h

≤ (1 + α)−(n+1) U0 2 + ν∆t ∇u0 2 + ν∆t ∇u−1 2

G4

h

16

h

+(Kα + (1 + β)−1) (1 + β)−n( T0 2 + κ∆t ∇T 0 2 + κ∆t ∇T −1 2)

G4

h

16

h

+(K + 1) max{ 8Cp2 , 2κ−1∆t } γ 2

α

Clκ2 3

L∞(R+ ,Wh∗ )

(30)

16Cp2 4ν−1∆t

+ max{ Clν2 ,

3

} f L∞(R+,V∗ ) h

where K = 49CuCp2ν−1Ri2∆t , U = [u0 u−1 u−2], T = [T 0 T −1 T −2],

α

α

0

hh

h

0

hh

h

α = min{ Clν∆t , 3 }, β = min{ Clκ∆t , 3 },

16Cp2 4

16Cp2 4

Proof. The proof starts with a stability bound for temperature. Letting Sh = Thn+1

in

(28)

and

using

the

skew

symmetry

property

b

2

(u∗

,

T

n+1 h

,

T

n+1 h

)

=

0.

Along

with

(7) and multiplying with ∆t yields

Tn+1

2 G

−

(31)

Tn 2 + 1 G 12

Thn+1 − 3Thn + 3Thn−1 − Thn−2 2 + κ∆t ∇Thn+1 2 = ∆t(γn+1, Thn+1).

where Tn+1 = [Thn+1 Thn Thn−1] and Tn = [Thn Thn−1 Thn−2] . Applying the Cauchy-Schwarz and Young’s inequalities lead to

Tn+1 2 − Tn 2 + 1 T n+1 − 3T n + 3T n−1 − T n−2 2 + κ∆t ∇T n+1 2

G

G 12 h

h

h

h

2

h

(32)

≤ κ−1∆t γn+1 2 .

2

Wh∗

Adding κ∆4 t ∇Thn 2 + κ1∆6t ∇Thn−1 2 to the both sides and dropping the nonnegative terms produce

(33)

( Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2) + κ∆t

G4

h

16

h

16

∇Thn+1 2

+ ∇Thn 2 + ∇Thn−1 2 + 3κ1∆6 t ∇Thn+1 2 + κ∆8 t ∇Thn 2

≤ T 2 + κ∆t ∇T n 2 + κ∆t ∇T n−1 2 + κ−1∆t γn+1 2 .

nG

4

h

16

h

2

Wh∗

The estimation of (33) follows closely that of (21). Again, the last ﬁve terms of the left hand side of (33) can be rearranged by the applying the Poincar´e-Friedrichs and the equivalent norm property (3) as

(34)

κ1∆6t ( ∇Thn+1 2 + ∇Thn 2 + ∇Thn−1 2)

+ 3κ1∆6 t ∇Thn+1 2 + κ∆8 t ∇Thn 2

≥ β( Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2),

G4

h

16

h

LONG TIME STABILITY OF A LINEARLY EXTRAPOLATED BLEBDF SCHEME

33

where β = min{ Clκ∆t , 3 }. Arguing as before and inserting (34) in (33) and using 16Cp2 4

induction leads to

(35)

Tn+1 2 + κ∆t ∇T n+1 2 + κ∆t ∇T n 2

G4

h

16

h

≤ (1 + β)−(n+1)( T0 2 + κ∆t ∇T 0 2 + κ∆t ∇T −1 2)

G4

h

16

h

+ κ−21β∆t γ 2L∞(R+,Wh∗),

which is the result for long the time stability for the temperature. Letting vh = unh+1 in (26), and qh = pnh+1 in (27), one obtains

35 unh+1

−

25 unh

+

uhn−1

−

1 6

uhn−2

,

un+1

∆t

h

+ ν ∇unh+1 2

(36)

= Ri (T ∗e2, unh+1) + (f n+1, unh+1).

Repeating the arguments used for obtaining (19) yields

(37)

Un+1 2 − Un 2 + 1 un+1 − 3un + 3un−1 − un−2 2

G

G 12 h

h

h

h

+ ν∆2 t ∇unh+1 2

≤ Cp2ν−1Ri2∆t T ∗ 2 e2 2 + ν−1∆t f n+1 Vh∗ .

Note that e2 is a unit vec√tor, i.e., e2 = 1 and using the deﬁnition of T ∗ and (3), we get T ∗ ≤ 7 Tn ≤ 7 Cu Tn G. Then, one has

(38)

Un+1

2 G

−

Un 2 + 1 un+1 − 3un + 3un−1 − un−2 2 + ν∆t

G 12 h

h

h

h

2

≤ 49CuCp2ν−1Ri2∆t Tn 2G + ν−1∆t f n+1 Vh∗ .

∇unh+1 2

Arguing as in (20), one gets the following estimation for (38)

(39)

(1 + α)( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ ( Un 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2)

G4

h

16

h

+49CuCp2ν−1Ri2∆t Tn 2G + ν−1∆t f n+1 Vh∗

where α = min{ Clν∆t , 3 }. Putting n instead of (n + 1) in (35) and inserting it in 16Cp2 4

(39) yields

(40)

(1 + α)( Un+1 2 + ν∆t ∇un+1 2 + ν∆t ∇un 2)

G4

h

16

h

≤ ( Un 2 + ν∆t ∇un 2 + ν∆t ∇un−1 2)

G4

h

16

h

+49CuC2ν−1Ri2∆t (1 + β)−n( T0 2 + κ∆t ∇T 0 2

p

G4

h

+ κ1∆6t ∇Th−1 2) + κ−21β∆t γ 2L∞(R+,Wh∗)

+ ν−1∆t f n+1 V∗ . h

Applying induction and adding (35) to (40) completes the proof of the theorem.