## Decoupled Modular Regularized VMS-POD for Darcy-Brinkman Equations

Fatma G. Eroglu, Songul Kaya, and Leo G. Rebholz

Abstract—We extend the post-processing implementa- tion of a projection based variational multiscale (VMS) method with proper orthogonal decomposition (POD) to flows governed by double diffusive convection. In the method, the stabilization terms are added to momentum equation, heat and mass transfer equations as a com- pletely decoupled separate steps. The theoretical analyses are presented. The results are verified with numerical tests on a benchmark problem.

Index Terms—post-process, variational multiscale, proper orthogonal decomposition, double-diffusive, re- duced order models.

I. INTRODUCTION

In this study, the Darcy-Brinkman equations with double diffusive convection is considered. The dimen- sionless form is given as:

ut− 2ν∇ · Du + (u · ∇)u + Da^{−1}u
+∇p = (β_{T}T + β_{C}C) g in (0, τ ] × Ω,

∇ · u = 0 in (0, τ ] × Ω,
u = 0 in (0, τ ] × ∂Ω,
Tt+ u · ∇T = γ∆T in (0, τ ] × ∂Ω,
Ct+ u · ∇C = Dc∆C in (0, τ ] × ∂Ω,
T, C = 0 on Γ_{D},

∇T · n = ∇C · n = 0 on ΓN,
u(0, x) = u_{0}, T (0, x) = T_{0}, C(0, x) = C_{0} in Ω,
(1)
where u(t, x), p(t, x), T (t, x), C(t, x) are the fluid
velocity, the pressure, the temperature, and the con-
centration fields, respectively. Let Ω ⊂ R^{d}, d ∈ {2, 3}

be a confined porous enclosure with polygonal bound- ary ∂Ω and ΓN be a regular open subset of the boundary and ΓD = ∂Ω \ ΓN. The initial veloc- ity, temperature and concentration fields are given as u0, T0, C0. The parameters in (1) are the kinematic viscosity ν > 0, inversely proportional to Re, the thermal diffusivity γ > 0, the velocity deformation

Department of Mathematics, Middle East Technical University, 06800 Ankara, Turkey, Department of Mathematics, Faculty of Sci- ence, Bartın University, 74110 Bartın, Turkey; fguler@bartin.edu.tr Department of Mathematics, Middle East Technical University, 06800 Ankara, Turkey; smerdan@metu.edu.tr.

Department of Mathematical Sciences, Clemson University, Clemson, SC, 29634; rebholz@clemson.edu; Partially supported by NSF DMS1522191 and U.S. Army Grant 65294-MA.

tensor Du = (∇u + ∇u^{T})/2, the mass diffusivity
D_{c} > 0, the Darcy number Da, and the gravitational
acceleration vector g. The solutal and the thermal
expansion coefficients are βC, and βT, respectively.

The dimensionless parameters are the Prandtl number P r, the Darcy number Da, the buoyancy ratio N , the Lewis number Le, the Schmidt number Sc, and the thermal and solutal Grashof numbers GrT and GrC, respectively. Here H is the cavity height, k the permeability, and ∆T and ∆C are the temperature and the concentration differences, respectively.

Double diffusive convection represents a form of convection driven by two different potentials with different diffusion rates. It is very important in many applications such as oceanography, meteorology and geology. The physical model is formed by forcing of momentum with both heat and mass transfer. Darcy terms defines the porosity of domain.

The basic challenges of the Darcy-Brinkman scheme come from the Navier Stokes equations (NSE). This is due to the complex behaviour of the flow at high Re and the absence of the analytical solution of NSE.

Combining momentum equation with mass and heat transfer equations makes the problem more difficult.

Thus, solving Darcy Brinkman equations accurately and efficiently remains a challenge for the compu- tational fluid dynamics community. Furthermore, the use of full order methods lead to large degrees of freedom. This causes complex algebraic systems and high computational time. To address this issue, model order reduction techniques are used.

In this study, we use Galerkin based proper orthog- onal decomposition (POD) method. The idea is find most energetic structure in the sytem which represent the snapshots. This idea can be found in Karhunen L´oeve expansion [1] , principal component analysis [2] and singular value decomposition [3]. Optimal POD basis functions are obtained by using finite element solution. Since, POD uses only most energetic structure in the system, it decreases the computational cost, process time and complexity of system. Thanks to the significant advantages, POD has been found ef- ficient method for different multhpysics problems [4], [5], [6]. Hence the application model order reduction with POD methodology to Darcy-Brinkman scheme is significant.

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

In this system, heat transfer is expressed with Rayleigh number (Ra) which is defined as ratio of buoyancy term to viscous term. The magnitude of the Ra indicates whether the flow is laminar or turbu- lent. For high Ra, the instability occurs due to the emergence of convection cells. Thus, the behaviour of the flow becomes turbulent. For such case, the VMS method can be used to eliminate the oscillation and stabilize the convective terms. Recent works [7], [8], [9], [10] show VMS-POD increase numerical accuracy.

The basic points of VMS are separation of scales as resolved small scale and resolved large scale, and eliminate the oscillations in small scales by using projection. Separation of scales is a challenge in the method. Selection of POD basis functions in descend- ing order are a remedy to this difficulty, i.e. small and large scales are decomposed naturally in POD method.

Many complex flows are solved by legacy codes, so it may be difficult to implement a new method in these flows. For such cases, the post-processing methods are easily added to legacy codes. Hence, the main objective of our study is application the post-processing VMS-POD idea by adding a separate, uncoupled and modular stabilization step for POD solution of Darcy-Brinkman system in each time step.

This work is arranged as follows. Section 2 presents the continuous variational formulation of the dou- ble diffusive Darcy-Brinkman system (1) and its dis- cretization, and here the VMS-POD variational formu- lation is defined. Section 3 is devoted to the numerical analysis of the VMS-POD formulation. Finally, Sec- tion 4 concludes the work with a summary.

II. FULL ORDER MODEL FOR THE DOUBLE DIFFUSIVEDARCY-BRINKMAN SYSTEM

Throughout the work standard notations for Sobolev
spaces and their norms will be used. The norm in
(H^{k}(Ω))^{d} is denoted by k·kk and the norms in
Lebesgue spaces (L^{p}(Ω))^{d}, 1 ≤ p < ∞, p 6= 2 by
k·kL^{p}. The space L^{2}(Ω) is equipped with the norm
and inner product k·k and (·, ·), respectively, and for
these we drop the subscripts. The continuous velocity,
pressure, temperature and concentration spaces are
denoted by

X := (H^{1}_{0}(Ω))^{d}, Q := L^{2}_{0}(Ω),
W := {S ∈ H^{1}(Ω) : S = 0 on ΓD},

Ψ := {Φ ∈ H^{1}(Ω) : Φ = 0 on ΓD},
and the divergence free space given as

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

We denote the dual space of X by H^{−1} with norm
kf k_{−1}= sup

v∈X

|(f , v)|

k∇vk.

The following notations are utilized for discrete norms

|||w|||∞,p := max

0≤n≤Mkw^{n}kp,

|||w|||m,p :=

∆t

M

X

n=0

kw^{n}k^{m}_{p}

1/m

. The variational formulation of (1) reads as follows:

Find u : (0, τ ] → X, p : (0, τ ] → Q, T : [0, τ ] → W and C : [0, τ ] → Ψ satisfying

(ut, v) + 2ν(Du, Dv) + b^{1}(u, u, v)
+(Da^{−1}u, v) − (p, ∇ · v)

= βT(gT, v) + βC(gC, v), (2)
(T_{t}, S) + b_{2}(u, T, S) + γ(∇T, ∇S) = 0, (3)
(Ct, φ) + b3(u, C, Φ) + Dc(∇C, ∇φ) = 0, (4)
for all (v, q, S, Φ) ∈ (X, Q, W, Ψ), where

b_{1}(u, v, w) := 1

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

2(((u · ∇)T, S) − ((u · ∇)S, T )) , b3(u, C, Φ) := 1

2(((u · ∇)C, Φ) − ((u · ∇)Φ, C)) , represent the skew-symmetric forms of the convective terms.

We consider a conforming finite element method for
(2)-(4), with spaces Xh ⊂ X, Qh ⊂ Q, Wh ⊂ W
and Ψh⊂ Ψ. We also assume that the pair (Xh, Qh)
satisfies the discrete inf-sup condition. It will also be
assumed for simplicity that the finite element spaces
X_{h}, Wh, Ψh are composed of piecewise polynomials
of degree at most m and Qhis composed of piecewise
polynomials of degree at most m − 1. In addition,
we assume that the spaces satisfy the interpolation
approximation properties. The discretely divergence
free space for (Xh, Qh) pairs is given by

Vh= {vh∈ Xh: (∇ · vh, qh) = 0, ∀qh∈ Qh}. (5) The inf-sup condition implies that the space Vh is a closed subspace of Xh and the formulation above involving Xh and Qh is equivalent to the following Vh formulation: Find (uh, Th, Ch) ∈ (Vh, Wh, Ψh) satisfying

(uh,t, vh) + 2ν(Du^{h}, Dv^{h})

+b_{1}(u_{h}, u_{h}, v_{h}) + (Da^{−1}u_{h}, v_{h})

= βT(gTh, vh) + βC(gCh, vh), (6) (Th,t, Sh) + b2(uh, Th, Sh)

+γ(∇Th, ∇Sh) = 0, (7) (Ch,t, Φh) + b3(uh, Ch, Φh)

+D_{c}(∇C_{h}, ∇Φ_{h}) = 0, (8)
for all (vh, S_{h}, Φ_{h}) ∈ (V_{h}, W_{h}, ψ_{h}).

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

The goal of the POD is to find low dimensional bases for velocity, temperature, concentration by solv- ing the minimization problems. The solution of the problem is obtained by using the method of snapshots.

We note that all eigenvalues are sorted in descending
order. Thus, the basis functions {ψ_{i}}^{r}_{i=1}^{1} , {φ_{i}}^{r}_{i=1}^{2} and
{η_{i}}^{r}_{i=1}^{3} correspond to the first r_{1}, r_{2} and r_{3} largest
eigenvalues {λi}^{r}_{i=1}^{1} , {µi}^{r}_{i=1}^{2} , {ξi}^{r}_{i=1}^{3} of the velocity,
the temperature, the concentration, respectively. For
simplicity, we will denote POD spaces using just r
instead of r1, r2 and r3. However, in the analysis, we
are careful to distinguish that these parameters can be
chosen independently.

Following spaces are needed for VMS-POD formu- lation.

Xr = span{ψ_{1}, ψ_{2}, . . . , ψ_{r}_{1}}, (9)
W_{r} = span{φ_{1}, φ_{2}, . . . , φ_{r}_{2}}, (10)
Ψr = span{η1, η2, . . . , ηr_{3}}, (11)
X_{R} = span{ψ_{1}, ψ_{2}, . . . , ψ_{R}

1}, (12)

WR = span{φ1, φ2, . . . , φR_{2}}, (13)
ΨR = span{η1, η2, . . . , ηR_{3}}, (14)
and

L_{R,u}= ∇X_{R}, L_{R,T} = ∇W_{R}, L_{R,C} = ∇Ψ_{R}. (15)
Note that by construction XR ⊆ Xr ⊂ Vh ⊂ X,
W_{R} ⊆ W_{r} ⊂ W_{h} ⊂ W and Ψ_{R} ⊆ Ψ_{r} ⊂ Ψ_{h} ⊂ Ψ.

The error term is decomposed by using the L^{2}projec-
tion P_{w,r} which is defined by

(w − P_{w,r}w, v_{r}) = 0, (16)
for all test functions vr which is in POD spaces.

The following lemma will be used to bound the POD projection error. The proof can be found in [8].

Lemma II.1. For true solution w^{n}at timet^{n}, we have
1

M

M

X

n=1

kw^{n}− Pw,rw^{n}k^{2}≤ C

h^{2m+2}|||w|||^{2}_{2,m+1}

+

d

X

i=r+1

λi

, (17) 1

M

M

X

n=1

k∇(w^{n}− Pw,rw^{n})k^{2}≤ C
(h^{2m}
+kSw,rk2h^{2m+2})|||w|||^{2}_{2,m+1}+ε^{2}_{u}

, (18)

whereεw,r=

s d

P

i=r+1

kψ^{w}_{i} k^{2}_{1}λ^{w}_{i} .

Now, we state the POD-Galerkin (POD-G) formu- lation of the Darcy-Brinkman double diffusive system.

Find (ur, Tr, Cr) ∈ (Xr, Wr, Ψr) satisfying
(ur,t, vr) + 2ν(Du^{r}, Dv^{r})

+ b1(ur, ur, vr) + (Da^{−1}ur, vr)

= β_{T}(gT_{r}, v_{r}) + β_{C}(gC_{r}, v_{r}), (19)
(Tr,t, Sr) + b2(ur, Tr, Sr)

+ γ(∇Tr, ∇Sr) = 0, (20)
(C_{r,t}, Φ_{r}) + b_{3}(u_{r}, C_{r}, Φ_{r})

+ Dc(∇Cr, ∇Φr) = 0, (21)
for all (vr, S_{r}, Φ_{r}) ∈ (X_{r}, W_{r}, Ψ_{r}).

We equip this system (19)-(21) with the BDF2 tem- poral discretization. We consider adding the decoupled VMS stabilization from [10], [11], where in effect ad- ditional viscosity gets added to the smaller R1, R2, R3

velocity, temperate and concentration modes in a post- processing step.

In order to prove an error estimate for the error
between the true solution and the VMS-POD solution
of the system, we use the L^{2} projection operators
Pu,r : L^{2} → LR,u, PT ,r : L^{2}→ LR,T, PC,r : L^{2} →
LR,C. They are defined by

(u − Pu,Ru, vR) = 0, (T − PT ,RT, SR) = 0, (C − PC,RC, ζR) = 0,

(22)

for all (vR, SR, ζR) ∈ (LR,u, LR,T, LR,C) For sim- plicity we use ePw,R instead of I − Pw,R.

Specifically, we post-process (u^{n+1}_{r} , T_{r}^{n+1}, C_{r}^{n+1})
by solving the following algorithm. Let g ∈
L^{2}(0, τ ; H^{−1}(Ω)) and initial conditions

(u0, T0, C0), (u1, T1, C1) ∈ ((L^{2}(Ω))^{d}, L^{2}(Ω), L^{2}(Ω))
be given in (Xr, W_{r}, Ψ_{r}).

Algorithm II.1. The post-processing VMS-POD ap- proximation for double diffusive system (1) given as:

Step 1: Find (w_{u,r}^{n+1}, w_{T ,r}^{n+1}, w_{C,r}^{n+1}) ∈ (X_{r}, W_{r}, Ψ_{r})
satisfying

(3w^{n+1}_{u,r} − 4u^{n}_{r}+ u^{n−1}_{r}

2∆t , vr) + 2ν(Dw^{n+1}u,r , Dv^{r})
+b_{1}(w^{n+1}_{u,r} , w_{u,r}^{n+1}, v_{r}) + (Da^{−1}w_{u,r}^{n+1}, v_{r})

= βT(gw_{T ,r}^{n+1}, vr) + βC(gw^{n+1}_{C,r} , vr), (23)

(3w_{T ,r}^{n+1}− 4T_{r}^{n}+ T_{r}^{n−1}

2∆t , Sr) + b2(w^{n+1}_{u,r} , w^{n+1}_{T ,r} , Sr)
+γ(∇w^{n+1}_{T ,r} , ∇Sr) = 0, (24)

(3w_{C,r}^{n+1}− 4C_{r}^{n}+ C_{r}^{n−1}

∆t , Φr) + b3(w^{n+1}_{u,r} , w^{n+1}_{C,r} , Φr)
+Dc(∇w^{n+1}_{C,r} , ∇Φr) = 0, (25)

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

for all (vr, Sr, Φr) ∈ (Xr, Wr, Ψr).

Step 2: Find (u^{n+1}_{r} , T_{r}^{n+1}, C_{r}^{n+1}) ∈ (X_{r}, W_{r}, Ψ_{r}),
for all (v_{r}, S_{r}, Φ_{r}) ∈ (X_{r}, W_{r}, Ψ_{r}):

w^{n+1}_{u,r} − u^{n+1}_{r}

∆t , vr

=

α1Peu,R∇(u^{n+1}_{r} + w^{n+1}_{u,r} )

2 , ePu,R∇vr

, (26)

w^{n+1}_{T ,r} − T_{r}^{n+1}

∆t , Sr

=

α2PeT ,R∇(T_{r}^{n+1}+ w^{n+1}_{T ,r}

2 ), ePT ,R∇Sr

, (27)

w^{n+1}_{C,r} − C_{r}^{n+1}

∆t , Φr

=

α_{3}Pe_{C,R}∇(C_{r}^{n+1}+ w_{C,r}^{n+1}

2 ), eP_{C,R}∇Φ_{r}

, (28)
whereP_{R} is the L^{2} projection into X_{R}, which is the
subset of Xr that is the span of the first R (< r)
velocity modes.

III. NUMERICAL ANALYSIS OF DOUBLE DIFFUSIVE

DARCY-BRINKMAN SYSTEM

This section is devoted to a derivation of the priori error estimation of (23)-(28). We first give the stability of solutions of (23)-(28).

Lemma III.1. (Stability) The Algorithm II.1 is stable for α1 ≤ 2ν, α2 ≤ 8γ, α3 ≤ 8Dc in the following sense: for any∆t > 0,

ku^{M +1}_{r} k^{2}+k2u^{M +1}_{r} − u^{M}_{r} k^{2}

+2ν∆tkDw^{M +1}u,r k^{2}+2Da^{−1}∆t

M

X

n=1

kw_{u,r}^{n+1}k^{2}

+2α1∆t

Peu,R∇(u^{M +1}_{r} + w^{M +1}_{u,r} )
2

2

≤ ku1k^{2}+k2u1− u0k^{2}+α1∆t

2 k∇u1k^{2}
+C^{∗}kgk^{2}_{∞}(β_{T}^{2}γ^{−1}(kT_{1}k^{2}+k2T_{1}− T0k^{2}
+α2∆t

2 k∇T1k^{2}) + β_{C}^{2}D^{−1}_{c} (kC1k^{2}
+k2C1− C0k^{2}+α3∆t

2 k∇C1k^{2})), (29)
kT_{r}^{M +1}k^{2}+4γ∆tk∇w_{T ,r}^{M +1}k^{2}+k2T_{r}^{M +1}− T_{r}^{M}k^{2}

+2α_{2}∆tk eP_{T ,R}∇(T_{r}^{M +1}+ w^{M +1}_{T ,r}

2 )k^{2}

≤ kT1k^{2}+k2T_{1}− T0k^{2}+α_{2}∆t

2 k∇T1k^{2}, (30)
kC_{r}^{M +1}k^{2}+4D_{c}∆tk∇w_{C,r}^{M +1}k^{2}+k2C_{r}^{M +1}− C_{r}^{M}k^{2}

+2α3∆tk ePC,R∇(C_{r}^{M +1}+ w_{C,r}^{M +1}

2 )k^{2}

≤ kC1k^{2}+k2C1− C0k^{2}+α3∆t

2 k∇C1k^{2}, (31)

where C^{∗}= min{ν^{−1}, Da}.

Proof:Letting Sr= w^{n+1}_{T ,r} in (24) and using skew
symmetry property yields

(3w_{T ,r}^{n+1}− 4T_{r}^{n}+ T_{r}^{n−1}
2∆t , w_{T ,r}^{n+1})
+b_{2}(w_{u,r}^{n+1}, w_{T ,r}^{n+1}, w_{T ,r}^{n+1})

+γ(∇w^{n+1}_{T ,r} , ∇w^{n+1}_{T ,r} ) = 0 (32)
Using the skew symmetry property,
b2(w^{n+1}_{u,r} , w^{n+1}_{T ,r} , w^{n+1}_{T ,r} ) = 0 and the identity:

α(3α − 4β + θ) = 1

2((α^{2}− β^{2})
+(2α − β)^{2}− (2β − θ)^{2}

+(α − 2β + θ)^{2}), (33)
we get

1

4∆tkw^{n+1}_{T ,r} k^{2}− 1

4∆tkT_{r}^{n}k^{2}+γk∇w_{T ,r}^{n+1}k^{2}
+ 1

4∆t(k2w^{n+1}_{T ,r} − T_{r}^{n}k^{2}−k2T_{r}^{n}− T_{r}^{n−1}k^{2})
+ 1

4∆tkw_{T ,r}^{n+1}− 2T_{r}^{n}+ T_{r}^{n−1}k^{2}= 0. (34)
Setting Sr= T_{r}^{n+1}+ w^{n+1}_{T ,r}

2 in (27) gives
kw^{n+1}_{T ,r} k^{2}−kT_{r}^{n+1}k^{2}=

2α2∆tk ePT ,R∇(T_{r}^{n+1}+ w^{n+1}_{T ,r}

2 )k^{2}. (35)
Substituting (35) in (34), multiplying with 4∆t and
adding and subtracting k2T_{r}^{n+1}− T_{r}^{n}k^{2} in (34) gives

kT_{r}^{n+1}k^{2}−kT_{r}^{n}k^{2}+4γ∆tk∇w_{T ,r}^{n+1}k^{2}
+(k2w_{T ,r}^{n+1}− T_{r}^{n}k^{2}−k2T_{r}^{n+1}− T_{r}^{n}k^{2})
+(k2T_{r}^{n+1}− T_{r}^{n}k^{2}−k2T_{r}^{n}− T_{r}^{n−1}k^{2})
+kw_{T ,r}^{n+1}− 2T_{r}^{n}+ T_{r}^{n−1}k^{2}

+2α_{2}∆tk eP_{T ,R}∇(T_{r}^{n+1}+ w^{n+1}_{T ,r}

2 )k^{2}= 0.(36)
Using the properties of L^{2}inner product and (35), we
rearrange the fourth and fifth terms in the right hand
side of (36). See [10] for details of the operations.

k2w^{n+1}_{T ,r} − Tr^{n}k^{2}−k2Tr^{n+1}− Tr^{n}k^{2}

= 8α2∆tk ePT ,R∇(w^{n+1}_{T ,r} + Tr^{n+1}− Tr^{n})

2 k^{2}

+8α2∆t( ePT ,R∇(Tr^{n}

2 ), ePT ,R∇(w^{n+1}_{T ,r} + Tr^{n+1}− Tr^{n})

2 )

(37)

Inserting (37) in (36) and applying the Cauchy-

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

Schwarz and Young’s inequalities gives

kT_{r}^{n+1}k^{2}−kT_{r}^{n}k^{2}+4γ∆tk∇w^{n+1}_{T ,r} k^{2}
+k2T_{r}^{n+1}− T_{r}^{n}k^{2}−k2T_{r}^{n}− T_{r}^{n−1}k^{2}
+2α2∆tk ePT ,R∇(T_{r}^{n+1}+ w_{T ,r}^{n+1}

2 )k^{2}
+kw_{T ,r}^{n+1}− 2T_{r}^{n}+ T_{r}^{n−1}k^{2}

≤ 2α2∆tk ePT ,R∇(T_{r}^{n}

2 )k^{2}. (38)
Dropping the positive seventh term and rearranging the
right hand side term of (38) yields

kT_{r}^{n+1}k^{2}−kT_{r}^{n}k^{2}+4γ∆tk∇w^{n+1}_{T ,r} k^{2}
+k2T_{r}^{n+1}− T_{r}^{n}k^{2}−k2T_{r}^{n}− T_{r}^{n−1}k^{2}
+2α_{2}∆tk eP_{T ,R}∇(T_{r}^{n+1}+ w_{T ,r}^{n+1}

2 )k^{2}

≤ 2α2∆tk ePT ,R∇(T_{r}^{n}+ w_{T ,r}^{n}
2 )k^{2}
+α2∆t

2 k ePT ,R∇w_{T ,r}^{n} k^{2}. (39)

Using k ePT ,Rk≤ 1 and α2≤ 8γ, we get

kT_{r}^{n+1}k^{2}−kT_{r}^{n}k^{2}+4γ∆tk∇w^{n+1}_{T ,r} k^{2}
+k2T_{r}^{n+1}− T_{r}^{n}k^{2}−k2T_{r}^{n}− T_{r}^{n−1}k^{2}
+2α2∆tk ePT ,R∇(T_{r}^{n+1}+ w_{T ,r}^{n+1}

2 )k^{2}

≤ 2α2∆tk ePT ,R∇(T_{r}^{n}+ w_{T ,r}^{n}
2 )k^{2}

+4γ∆tk∇w^{n}_{T ,r}k^{2}. (40)

Finally summing over the time step n = 1, . . . , M yields

kT_{r}^{M +1}k^{2}+4γ∆tk∇w^{M +1}_{T ,r} k^{2}+k2T_{r}^{M +1}− T_{r}^{M}k^{2}
+2α_{2}∆tk eP_{T ,R}∇(T_{r}^{M +1}+ w^{M +1}_{T ,r}

2 )k^{2}

≤ kT1k^{2}+k2T1− T0k^{2}+4γ∆tk∇w^{1}_{T ,r}k^{2}
+2α_{2}∆tk eP_{T ,R}∇(T1+ w^{1}_{T ,r}

2 )k^{2}. (41)

Using k eP_{T ,R}k≤ 1 and assuming w_{T ,r}^{1} = 0 and we get
the stated result (30).

In a similar manner, setting Φr= w^{n+1}_{C,r} in (25) and
using the assumption α3 ≤ 8D_{c} yields (31). Finally,
choosing vr = w^{n+1}_{u,r} in (23), using the polarization

identity, and multiplying both sides by 4∆t yields

kw^{n+1}_{u,r} k^{2}−ku^{n}_{r}k^{2}+k2w^{n+1}_{u,r} − u^{n}_{r}k^{2}

−k2u^{n}_{r} − u^{n−1}_{r} k^{2}+8ν∆tkDw^{n+1}u,r k^{2}
+kw^{n+1}_{u,r} − 2u^{n}_{r}+ u^{n−1}_{r} k^{2}

+4Da^{−1}∆tkw^{n+1}_{u,r} k^{2}

= 4∆tβ_{T}(gw^{n+1}_{T ,r} , w^{n+1}_{u,r} )

+4∆tβC(gw^{n+1}_{C,r} , w^{n+1}_{u,r} ). (42)

Note that if we let vr= ^{(u}

n+1
r +w^{n+1}_{u,r})

2 in (26), we have
kw^{n+1}_{u,r} k^{2}−ku^{n+1}_{r} k^{2}=

2α1∆tk ePu,R∇(u^{n+1}_{r} + w^{n+1}_{u,r}

2 )k^{2}. (43)
Insert (43) in (42), and apply Cauchy-Schwarz,
Young’s inequality and Poincar´e’s inequality, which
provides

ku^{n+1}_{r} k^{2}−ku^{n}_{r}k^{2}+k2w^{n+1}_{u,r} − u^{n}_{r}k^{2}

−k2u^{n}_{r} − u^{n−1}_{r} k^{2}+4ν∆tkDw^{n+1}u,r k^{2}
+kw^{n+1}_{u,r} − 2u^{n}_{r}+ u^{n−1}_{r} k^{2}

+2Da^{−1}∆tkw^{n+1}_{u,r} k^{2}
+2α1∆t

Peu,R∇(u^{n+1}_{r} + w^{n+1}_{u,r}

2 )

2

≤ C^{∗}kgk^{2}_{∞}(β^{2}_{T}∆tkw^{n+1}_{T ,r} k^{2}
+β_{C}^{2}∆tkw_{C,r}^{n+1}k^{2})

(44)

where C^{∗} = min{ν^{−1}, Da^{−1}}. Using the similar
argument with (40) for velocity and utilizing the
assumption α1≤ 2ν, we obtain

ku^{n+1}_{r} k^{2}−ku^{n}_{r}k^{2}+k2u^{n+1}_{r} − u^{n}_{r}k^{2}

−k2u^{n}_{r}− u^{n−1}_{r} k^{2}+4ν∆tkDw^{n+1}u,r k^{2}
+2Da^{−1}∆tkw^{n+1}_{u,r} k^{2}

+2α1∆t

Peu,R∇(u^{n+1}_{r} + w^{n+1}_{u,r}

2 )

2

≤ C^{∗}kgk^{2}_{∞}(β_{T}^{2}∆tk∇w^{n+1}_{T ,r} k^{2}

+β^{2}_{C}∆tk∇w_{C,r}^{n+1}k^{2}) + 4ν∆tkDwu,r^{n} k^{2}
+2α1∆t

Peu,R∇(u^{n}_{r} + w^{n}_{u,r}

2 )

2

(45)

Summing over the time steps and inserting (30) and

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

(31) yields

ku^{M +1}_{r} k^{2}+k2u^{M +1}_{r} − u^{M}_{r} k^{2}

+2ν∆tkDw^{M +1}u,r k^{2}+2Da^{−1}∆t

M

X

n=1

kw_{u,r}^{n+1}k^{2}

+2α1∆t

Peu,R∇(u^{M +1}_{r} + w^{M +1}_{u,r}

2 )

2

≤ ku^{1}_{r}k^{2}+k2u^{1}_{r}− u^{0}_{r}k^{2}+2ν∆tkDwu,r^{1} k^{2}
+2α1∆t

Peu,R∇(u^{1}_{r}+ w^{1}_{u,r}

2 )

2

+C^{∗}kgk^{2}_{∞}(β_{T}^{2}γ^{−1}(kT_{1}k^{2}+k2T_{1}− T_{0}k^{2}
+α2∆t

2 k∇T1k^{2}) + β_{C}^{2}D^{−1}_{c} (kC1k^{2}
+k2C1− C0k^{2}+α3∆t

2 k∇C1k^{2}) (46)
Using the assumption of w^{1}_{u,r} = 0 and k ePu,Rk≤ 1,
we get stated result (29).

The optimal asymptotic error estimation requires the following regularity assumptions for the true solution:

u ∈ L^{∞}(0, k; H^{m+1}(Ω)),
utt∈ L^{2}(0, T ; H^{1}(Ω)),
T, C ∈ L^{∞}(0, k; H^{m+1}(Ω)),

T_{tt}, C_{tt}∈ L^{2}(0, T ; H^{1}(Ω)),

p ∈ L^{∞}(0, k; H^{m}(Ω)). (47)
We define the discrete norms for v^{n} ∈ H^{p}(Ω), n =
0, 1, 2, ..., M as the following:

|||v|||∞,p:= max

0≤n≤Mkv^{n}kp,

|||v|||m,p:= (∆t

M

X

n=0

kv^{n}k^{m}_{p} )^{1/m}.

Theorem III.1. (Error Estimation) Suppose regularity assumptions (47) hold. Then for the sufficiently small

∆t, the error satisfies

ku^{M} − u^{M}_{r} k^{2}+kT^{M}− T_{r}^{M}k^{2}+kC^{M} − C_{r}^{M}k^{2}

≤ K

1 + h^{2m}+ (∆t)^{2}+

1 + kSu,rk2

+kST ,rk2+kSC,rk2+kSu,Rk2

+kS_{T ,R}k2+kS_{C,R}k2

h^{2m+2}

+

d

X

i=r1+1

(kψ_{i}k^{2}_{1}+1)λi+

d

X

i=r2+1

(kφik^{2}_{1}+1)µi

+

d

X

i=r_{3}+1

(kηik^{2}_{1}+1)ξi+

d

X

i=R1+1

kψ_{i}k^{2}_{1}λi

+

d

X

i=R_{2}+1

kφ_{i}k^{2}_{1}µ_{i}+

d

X

i=R_{3}+1

kηik^{2}_{1}ξ_{i}

. (48)

Proof: We begin the proof by deriving error
equations, subtracting from (2), (3), (4) to (23), (24),
(25) at time t^{n+1}, respectively, then we have

u^{n+1}_{t} −3w^{n+1}_{u,r} − 4u^{n}_{r} + u^{n−1}_{r}

2∆t , v_{r}

!

+ 2ν(D(u^{n+1}− w_{u,r}^{n+1}), Dvr)

+ b_{1}(u^{n+1}, u^{n+1}, v_{r}) − b_{1}(w^{n+1}_{u,r} , w^{n+1}_{u,r} , v_{r})
+ (Da^{−1}(u − w_{u,r}^{n+1}), vr) − (p^{n+1}, ∇ · vr)

= βT(g(T^{n+1}− w^{n+1}_{T ,r} ), vr)
+ βC(g(C^{n+1}− w^{n+1}_{C,r}), vr),

(49)
T_{t}^{n+1}−3w^{n+1}_{T ,r} − 4T_{r}^{n}+ T_{r}^{n−1}

2∆t , Sr

!

+ b_{2}(u^{n+1}, T^{n+1}, S_{r}) − b_{2}(w^{n+1}_{u,r} , w^{n+1}_{T ,r} , S_{r})
+ γ(∇(T^{n+1}− w_{T ,r}^{n+1}), ∇S_{r}) = 0, (50)

C_{t}^{n+1}−3w^{n+1}_{C,r} − 4C_{r}^{n}+ C_{r}^{n−1}

∆t , Φ_{r}

!

+ b3(u^{n+1}, C^{n+1}, Φr) − b3(w^{n+1}_{u,r} w^{n+1}_{C,r}, Φr)
+ Dc(∇(C^{n+1}− w_{C,r}^{n+1}), ∇Φr) = 0. (51)
The notations that are used in the proof are defined as
η^{n}_{u} := u^{n}− ˜U^{n}, φ^{n}_{u,r} := w^{n}_{u,r}− ˜U^{n},
θ^{n}_{u,r} := u^{n}_{r} − ˜U^{n}, E^{n}_{u,r} := u^{n}− w^{n}_{u,r}
e^{n}_{u,r} = u^{n}− u^{n}_{r},

η_{T}^{n} := T^{n}− ˜T^{n}, φ^{n}_{T ,r} := w^{n}_{T ,r}− ˜T^{n},
θ_{T ,r}^{n} := T_{r}^{n}− ˜T^{n}, E_{T ,r}^{n} := T^{n}− w_{T ,r}^{n}
e^{n}_{T ,r} = T^{n}− T_{r}^{n},

η_{C}^{n} := C^{n}− ˜C^{n}, φ^{n}_{C,r} := w^{n}_{C,r}− ˜C^{n},
θ^{n}_{C,r} := C_{r}^{n}− ˜C^{n}, E_{C,r}^{n} := C^{n}− w^{n}_{C,r}
e^{n}_{C,r} = C^{n}− C_{r}^{n},

where ( ˜U^{n}, ˜T^{n}, ˜C^{n}) are L^{2} projections of
(u^{n}, T^{n}, C^{n}) in (Xr, Wr, Ψr) at time t^{n}.

Letting S_{r}= φ^{n+1}_{T ,r} in (50), and reorganizing it, we
get

3E_{T ,r}^{n+1}− 4e^{n}_{T ,r}+ e^{n−1}_{T ,r}
2∆t , φ^{n+1}_{T ,r}

!

+γ(∇E_{T ,r}^{n+1}, ∇φ^{n+1}_{T ,r} ) + b_{2}(u^{n+1}, T^{n+1}, φ^{n+1}_{T ,r} )

−b2(w_{u,r}^{n+1}, w_{T ,r}^{n+1}, φ^{n+1}_{T ,r} )
+

T_{t}^{n+1}−3T^{n+1}− 4T^{n}+ T^{n−1}
2∆t , φ^{n+1}_{T ,r}

= 0.

(52)
Utilizing E_{T ,r}^{n+1} = η_{T}^{n+1}− φ^{n+1}_{T ,r} , e^{n}_{T ,r} = η_{T}^{n} − θ^{n}_{T ,r},

**IAENG International Journal of Applied Mathematics, 49:2, IJAM_49_2_01**

**______________________________________________________________________________________ **

(16) and (33), we get 1

4∆tkφ^{n+1}_{T ,r} k^{2}− 1

4∆tkθ^{n}_{T ,r}k^{2}
+ 1

4∆t(k2φ^{n+1}_{T ,r} − θ^{n}_{T ,r}k^{2}−k2θ^{n}_{T ,r}− θ^{n−1}_{T ,r} k^{2})
+ 1

4∆tkφ^{n+1}_{T ,r} − 2θ^{n}_{T ,r}+ θ_{T ,r}^{n−1}k^{2}+γk∇φ^{n+1}_{T ,r} k^{2}

= γ(∇η^{n+1}_{T} , ∇φ^{n+1}_{T ,r} ) + b_{2}(w_{u,r}^{n+1}, w_{T ,r}^{n+1}, φ^{n+1}_{T ,r} )

−b_{2}(u^{n+1}, T^{n+1}, φ^{n+1}_{T ,r} )

−

T_{t}^{n+1}−3T^{n+1}− 4T^{n}+ T^{n−1}
2∆t , φ^{n+1}_{T ,r}

. (53) Adding and subtracting ˜T in (27) on both sides to get

φ^{n+1}_{T ,r} − θ^{n+1}_{T ,r}

∆t , Sr

=

α_{2}Pe_{T ,R}∇(φ^{n+1}_{T ,r} + θ_{T ,r}^{n+1}+ 2 ˜T

2 ), eP_{T ,R}∇S_{r}

.
(54)
Setting S_{r}=^{φ}

n+1
T ,r +θ_{T ,r}^{n+1}

2 in (54) produces
kφ^{n+1}_{T ,r} k^{2}= kθ_{T ,r}^{n+1}k^{2}

+α2∆t

2 k ePT ,R∇(φ^{n+1}_{T ,r} + θ_{T ,r}^{n+1})k^{2}

+ ∆t(α2PeT ,R∇ ˜T^{n+1}, ePT ,R∇(φ^{n+1}_{T ,r} + θ^{n+1}_{T ,r} )). (55)
In a similar setting, for concentration we have

kφ^{n+1}_{C,r}k^{2}= kθ_{C,r}^{n+1}k^{2}
+α3∆t

2 k ePC,R∇(φ^{n+1}_{C,r} + θ_{C,r}^{n+1})k^{2}

+ ∆t(α3PeC,R∇ ˜C^{n+1}, ePC,R∇(φ^{n+1}_{C,r} + θ^{n+1}_{C,r})).

(56) Rewriting the nonlinear terms, we have

b_{2}(w^{n+1}_{u,r} , w_{T ,r}^{n+1}, φ^{n+1}_{T ,r} ) − b_{2}(u^{n+1}, T^{n+1}, φ^{n+1}_{T ,r} )

= −b2(η^{n+1}_{u} , T^{n+1}, φ^{n+1}_{T ,r} ) + b2(φ^{n+1}_{u,r} , T^{n+1}, φ^{n+1}_{T ,r} )

− b2(w^{n+1}_{u,r} , η_{T}^{n+1}, φ^{n+1}_{T ,r} ). (57)
After substituting (57) in (53) and multiplying 4∆t,
we bound the right hand side of (53) as follows

|γ(∇η_{T}^{n+1}, ∇φ^{n+1}_{T ,r} )| ≤ Cγk∇η^{n+1}_{T} k^{2}
+γ

10k∇φ^{n+1}_{T ,r} k^{2},

|b2(η^{n+1}_{u} , T^{n+1}, φ^{n+1}_{T ,r} )| ≤ C

γk∇η^{n+1}_{u} k^{2}k∇T^{n+1}k^{2}
+γ

10k∇φ^{n+1}_{T ,r} k^{2},

|b2(φ^{n+1}_{u,r} , T^{n+1}, φ^{n+1}_{T ,r} )| ≤ C

νγ^{2}kφ^{n+1}_{u,r} k^{2}k∇T^{n+1}k^{4}
+γ

10k∇φ^{n+1}_{T ,r} k^{2}
+νkDφ^{n+1}u,r k^{2},

|b2(w^{n+1}_{u,r} , η_{T}^{n+1}, φ^{n+1}_{T ,r} )| ≤ C

γk∇w^{n+1}_{u,r} k^{2}k∇η_{T}^{n+1}k^{2}
+γ

10k∇φ^{n+1}_{T ,r} k^{2},

|(α2PeT ,R∇ ˜T^{n+1}, ePT ,R∇(φ^{n+1}_{T ,r} + θ^{n+1}_{T ,r} ))|

≤ α_{2}k eP_{T ,R}∇ ˜T^{n+1}k^{2}+α2

4 k eP_{T ,R}∇(φ^{n+1}_{T ,r} + θ_{T ,r}^{n+1})k^{2}. (58)

For the last term in the right hand side of (53), using Taylor series expansion with the remainder in integral form with Cauchy Schwarz and the triangle inequality, we obtain

|(T_{t}^{n+1}−3T^{n+1}− 4T^{n}+ T^{n−1}

2∆t , φ^{n+1}_{T ,r} )|

≤ Cγ^{−1}∆tkT_{tt}k^{2}_{L}2(0,τ ;H^{1}(Ω))

+γ

10k∇φ^{n+1}_{T ,r} k^{2}. (59)

Inserting all bounds in (53), using ˜T^{n+1} = T^{n+1}−
η^{n+1}_{T} and summing over the time steps produces

kθ^{M +1}_{T ,r} k^{2}+k2θ_{T ,r}^{M +1}− θ_{T ,r}^{M} k^{2}
+α2∆t

4

M

X

n=1

k ePT ,R∇(φ^{n+1}_{T ,r} + θ^{n+1}_{T ,r} )k^{2}

+2γ∆t

M

X

n=1

k∇φ^{n+1}_{T ,r} k^{2}

≤ kθ_{T ,r}^{0} k^{2}+k2θ_{T ,r}^{1} − θ^{0}_{T ,r}k^{2}
+K∆t

γ

M

X

n=1

k∇ηTn+1

k^{2}
+α_{2}k eP_{T ,R}∇(T^{n+1}− η_{T}^{n+1})k^{2}
+γ^{−1}

M

X

n=1

k∇η_{u}^{n+1}k^{2}k∇T^{n+1}k^{2}

+ν^{−1}γ^{−2}

M

X

n=1

kφ^{n+1}_{u,r} k^{2}k∇T^{n+1}k^{4}

+γ^{−1}

M

X

n=1

k∇w^{n+1}_{u,r} k^{2}k∇η_{T}^{n+1}k^{2}

+ν

M

X

n=1

kDφ^{n+1}u,r k^{2}

+γ^{−1}∆tkTttk^{2}_{L}2(0,τ ;H^{1}(Ω))

. (60)