FULLY DISCRETE FINITE ELEMENT APPROXIMATION FOR THE STABILIZED GAUGE-UZAWA METHOD TO SOLVE THE BOUSSINESQ EQUATIONS

The stabilized Gauge-Uzawa method [SGUM], which is a 2nd order projection type algorithm used to solve Navier-Stokes equations, has been newly constructed in [11]. In this paper, we apply the SGUM to the evolution Boussinesq equations, which model the thermal driven motion of incompressible fluids. We prove that SGUM is unconditionally stable and we perform error estimations on the fully discrete finite element space via variational approach for the velocity, pressure and temperature, the three physical unknowns. We conclude with numerical tests to check accuracy and physically relevant numerical simulations, the Benard convection problem and the thermal driven cavity flow.


Introduction
The stabilized Gauge-Uzawa method [SGUM] is a 2nd order projection type method to solve the evolution Navier-Stokes equations. In this paper, we extend SGUM to the evolution Boussinesq equations: given a bounded polygon Ω in R d with d = 2 or 3, (1.1) u t + (u · ∇) u + ∇p − µ u + κµ 2 gθ = f , in Ω, ∇· u = 0, in Ω, θ t + u · ∇θ − λµ θ = b, in Ω, with initial conditions u(x, 0) = u 0 , θ(x, 0) = θ 0 in Ω, vanishing Dirichlet boundary conditions u = 0, θ = 0 on ∂Ω, and pressure mean-value ∫ Ω p = 0. The forcing functions f and b are given and g is the vector of gravitational acceleration. The nondimensional numbers µ = Re −1 and λ = P r −1 are reciprocal of the Reynolds and Prandtl numbers, respectively, whereas κ is the Grashof number. The Boussinesq system (1.1) describes fluid motion due to density differences which are in turn induced by temperature gradients: hot, and thus less dense, fluid tends to rise against gravity and cooler fluid falls in its place. The simplest governing equations are thus the Navier-Stokes equations for motion of an incompressible fluid, with forcing κµ 2 gθ due to buoyancy, and the heat equation for diffusion and transport of heat. Density differences are thus ignored altogether except for buoyancy.
The projection type methods are representative solvers for the incompressible flows and the Gauge-Uzawa method is a typical projection method. The Gauge-Uzawa method was constructed in [8] to solve Navier-Stokes equations and extended to more complicated problems, the Boussinesq equations in [9] and the non-constant density Navier-Stokes equations in [13]. However, most of studies for the Gauge-Uzawa method have been limited only for the first order accuracy backward Euler time marching algorithm. The second order Gauge-Uzawa method using BDF2 scheme was introduced in [12] and proved superiority for accuracy on the normal mode space, but we couldn't get any theoretical proof via energy estimate even stability and we suffer from weak stability performance on the numerical test. Recently, we construct SGUM in [11] which is unconditionally stable for semi-discrete level to solve the Navier-Stokes equations. The goal of this paper is to extend SGUM to the Boussinesq equations (1.1), which model the motion of an incompressible viscous fluid due to thermal effects [4,10]. We will estimate errors and stability on the fully discrete finite element space. The main difficulties in the fully discrete estimation arise from losing the cancellation law due to the failing of the divergence free condition of the discrete velocity function. The strategy of projection type methods compute first an artificial velocity and then decompose it to divergence free velocity and curl free functions. However the divergence free condition cannot be preserved in discrete finite element space, and so the cancellation law (1.9) can't be satisfied any more. In order to solve this difficulty, we impose the discontinuous velocity on across inter-element boundaries to make fulfill discrete divergence free velocity (1.9) automatically. We will discuss this issue at Remark 1.1 below. This discontinuity makes it difficult to treat non-linear term and to apply the integration by parts, because the discontinuous solution is not included in H 1 (Ω). So we need to hire technical skills in proof of this paper.
One more remarkable discovery is in the second numerical test at the last section which is the Benard convection problem with the same setting in [9]. In this performance, we newly find out that the number of circulations depends on the time step size τ . We obtain similar simulation within [9] for a relatively larger τ , but the behavior is changed for the small τ . So we conclude that the numerical result of the Benard convection problem in [9] is not an eventual solution and a more smaller time marching step is required to get the desired simulation.
In order to introduce the finite element discretization we need further notations. Let H s (Ω) be the Sobolev space with s derivatives in L 2 (Ω), set L 2 (Ω) = ( L 2 (Ω) ) d and H s (Ω) = (H s (Ω)) d , where d = 2 or 3, and denote by L 2 0 (Ω) the subspace of L 2 (Ω) of functions with vanishing meanvalue. We indicate with · s the norm in H s (Ω), and with · , · the inner product in L 2 (Ω). Let T = {K} be a shape-regular quasi-uniform partition of Ω of meshsize h into closed elements K [1,2,5]. The vector and scalar finite element spaces are: where P(K), Q(K), and R(K) are spaces of polynomials with degree bounded uniformly with respect to K ∈ T [2,5]. We stress that the space P h is composed of continuous functions to ensure the crucial equality Using the following discrete counterpart of the form N (u, v, w) : we now introduce the fully discrete SGUM to solve the evolution Boussinesq equations (1.1).
Step 2: Find ψ n+1 h ∈ P h as the solution of Step 3: Update u n+1 h and q n+1 h ∈ P h according to Step 4: Update pressure p n+1 Step 5: Find θ n+1 h ∈ T h as the solution of . We note that u n+1 h is a discontinuous function across inter-element boundaries and that, in light of (1.4) and (1.5), u n+1 h is discrete divergence free in the sense that We now summarize the results of this paper along with organization. We introduce appropriate Assumptions 1-5 in §2 and introduce well known lemmas. In §3, we prove the stability result: Theorem 1 (Stability). The SGUM is unconditionally stable in the sense that, for all τ > 0, the following a priori bound holds: (1.10) u N +1 .
We then will carry out the following optimal error estimates through several lemmas in §4.
Theorem 2 (Error estimates). Suppose the exact solution of (1.1) is smooth enough and τ ≤ Ch. If Assumptions 1 and 3-5 below hold, then the errors of Algorithm 1 will be bound by Moreover, if Assumptions 2 also hold, then we have We note that the condition τ ≤ Ch in Theorem 2 can be omitted for the linearized Boussinesq equations (see Remark 4.4). Finally, we perform numerical tests in §5 to check accuracy and physically relevant numerical simulations, the Benard convection problem and the thermal driven cavity flow.

Preliminaries
In this section, we introduce 5 assumptions and well known lemmas to use in proof of main theorems. We resort to a duality argument for which is the stationary Stokes system with vanishing boundary condition v = 0, as well as the Poisson equation with boundary condition ω = 0. We now state a basic assumption about Ω.
Assumption 1 (Regularity of (v, r) and ω). The unique solutions {v, r} of (2.1) and ω of (2.2) satisfy We remark that the validity of Assumption 1 is known if ∂Ω is of class C 2 [3,6], or if ∂Ω is a two-dimensional convex polygon [7], and is generally believed for convex polyhedral [6]. We impose the following properties for relations between the spaces V h and P h .
Assumption 3 (Shape regularity and quasiuniformity [1,2,5]). There exists a constant C > 0 such that the ratio between the diameter h K of an element K ∈ T and the diameter of the largest ball contained in K is bounded uniformly by C, and h K is comparable with the meshsize h for all K ∈ T.
In order to launch Algorithm 1, we need to set (u 1 h , p 1 h , θ 1 h ) via any first order methods which hold the following conditions. Assumption 4 (The setting of the first step values). Let (u(t 1 ), p(t 1 ), θ(t 1 )) be the exact solution of (1.1) at t = t 1 . The first step value The following elementary but crucial relations are derived in [14].
Let now (v h , r h ) ∈ V h × P h indicate the finite element solution of (2.1), namely, Then we can find the well known lemmas in [1,2,5]: Proof. Inequality (2.4) is standard [1,2,5]. To establish (2.5), we just deal with the L ∞ -norm since the other can be treated similarly.
as a consequence of an inverse estimate. This completes the proof. Remark 2.4 (H 1 stability of r h ). The bound ∇r h 0 ≤ C ( v 2 + r 1 ) is a simple by-product of (2.1). To see this, we add and subtract I h q, use the stability of I h in H 1 , and observe that (2.4) Lemma 2.5 (Error Estimates for Poisson's equation [1,5]). Let ω ∈ H 1 0 (Ω) be the solution of the Poisson's equation (2.2) and ω h ∈ T h be its finite element approximation If Assumptions 1, 3, and 5 hold, then there exists a positive constant C satisfies We finally state without proof several properties of the nonlinear form N . In view of (1.2), we have the following properties of N for all Applying Sobolev imbedding Lemma yields the following useful results.
We will use the following algebraic identities frequently to treat time derivative terms.

Proof of stability
We show that the SGUM is unconditionally stable via a standard energy method in this section. We start to prove stability with rewriting the momentum equation (1.3) by using (1.5) and (1.7) as follows: We now choose w h = 4τ u n+1 h and apply (2.11) to obtain We give thanks to (2.7) for eliminating the convection term. In light of u n+1 h , (1.9) and (2.12) yield ) .
Before we estimate A 2 , we evaluate an inequality via choosing , and so (1.6) and (2.13) lead us Clearly, we have We now attack A 4 with u n+1 which comes form (1.5) and (1.9). Then we arrive at in (1.8) and use (2.11), then we obtain Inserting A 1 -A 4 back into (3.1), adding with (3.2), and then summing over n from 1 to N lead (1.10) by help of discrete Gronwall inequality and the equality u n+1 . So we finish the proof of Theorem 1.

Error estimates
We prove Theorem 2 which is error estimates for SGUM of Algorithm 1. This proof is carried out through several lemmas. We start to prove with defining We also define Θ n+1 h ∈ T h as the solution of From Lemma 2.3, we can deduce ) , In conjunction with the definition |||·||| in (2.5), we can derive We now carry out error evaluation by comparing (4.10) below with (1.3)-(1.7) and then by comparing (4.15) and (1.8). We derive strong estimates of order 1 and this result is instrumental in proving weak estimates of order 2 for the errors Then, in conjunction with (4.2), we can readily get of the same accuracy for the errors Additionally, we denote We readily obtain the following properties as well as, from (1.9), whence we deduce crucial orthogonality properties: We now estimate the first order accuracy for velocity and temperature in Lemma 4.1, and then the 2nd order accuracy for time-derivative of velocity and temperature in Lemma 4.3. The result of Lemma 4.1 is instrumental to treat the convection term in proof of Lemma 4.3. We will use Lemmas 4.1 and 4.3 to prove optimal error accuracy in Lemma 4.5. Finally, we will prove pressure error estimate in Lemma 4.6.
Proof. We resort to the Taylor theorem to write (1.1) as follows: and (4.9) We replace the pressure p n h term in (1.3) with (1.7) and then subtract it from (4.10) to obtain in (4.11) and using (4.5) and (2.11) yield We now estimate all the terms from A 1 to A 7 respectively. To tackle A 1 , we first add and subtract 2u(t n ) − u(t n−1 ) to obtain = 0, which comes from (2.7), the second term of A 1 can be replaced by .
If we apply Lemma 2.6, then we can readily obtain Since we have u(t n+1 ) 2 + G n+1 ≤ M according to (4.4), we arrive at In light of (4.3) and (4.5), A 2 becomes In order to estimate A 3 and A 4 , we note that the cancellation law (1.9) gives In conjunction with the definition ε n+1 , A 4 can be evaluated by If we now apply inequality (a + b) 2 In light of ∇ · u(t n+1 ) = 0 and (2.13), (1.6) and (4.6) yield Also we readily get The Hölder inequality and (4.7) yield Inserting the above estimates into (4.12) gives (4.14) On the other hand, the definition (4.1) of Θ n+1 h ∈ T h leads a weak formula of (4.9) as, for all φ h ∈ T h , We now subtract (1.8) from (4.15) to derive (4.16) 2 dt which is (4.7) to attack A 9 , Inserting the above estimates into (4.17) gives Adding 2 inequalities (4.14) and (4.18) and summing over n from 1 to N imply ) dt Because of q 1 = 0, we have ε 1 h = P 1 h + 3 2τ ψ 1 h = e 1 h , and thus Assumption 4 yields ∇ε 1 h 2 0 ≤ C and the first 4 terms in the right hand side can be bound by Assumption 4 and properties E 0 = 0 and q 1 = 0 which are directly deduced from the conditions in Algorithm 1. As well as the remaining terms can be treated by the discrete Gronwall lemma. Finally, in conjunction with (4.5), we conclude the desired result and complete this proof.

Remark 4.2 (Optimal estimation).
In order to get optimal accuracy, we must get rid of the terms of A 3 and A 4 by applying duality argument in Lemma 4.5. To do this, we first evaluate the errors for time-derivative of velocity and temperature in Lemma 4.3. Thus we need to evaluate optimal initial errors for the case n = 1, and so we have to recompute again (4.13). We start to rewrite A 4 as dt.
In light of Assumption 4, we arrive at (4.19) We now start to estimate errors for time-derivative of velocity.

Remark 4.4 (The condition τ ≤ Ch).
The assumption τ ≤ Ch requires to control convection terms which are used at only (4.23) and (4.26), so we can omit this condition for the linearized Boussinesq equations. However, this condition can not be removed for non-linear equation case, because (4.22) below must be bounded by h.
Proof. We subtract two consecutive formulas of (4.11) and impose w h = 4τ δ E n+1 h to obtain (4.21) We now estimate each term A 1 to A 7 separately. The convection term A 1 can be rewritten as follows: In estimating convection terms, we will use Lemma 2.6 frequently without notice. We recall u(t) 2 ≤ C to obtain The result in Lemma 4.1, 2E n − E n−1 0 ≤ C (τ + h), is essential to treat the next 2 convection terms. Invoking (2.9), we have = 0 which comes from (2.7). Then we obtain To attack B 1 , we first note Lemma 2.1 which is, for any which is the result of Lemma 4.1, then we can conclude, in light of (2.10), (4.23) In light of Hölder inequality, (4.3) yields

Integral by parts leads
In order to tackle A 4 , marking use of δ E n+1 We readily get If we now apply inequality (a + b) 2 ≤ 4a 2 + 4 3 b 2 , then we can have 4τ 2 3 ∇δδε n+1 Invoking (1.6), (2.13) and (4.7) lead Finally, the A 6 term becomes For the new term A 7 we argue as in Lemma 4.1 to arrive at We now insert the above estimates into (4.21) to obtain (4.24) To evaluate errors of δϑ n+1 h , we subtract two consecutive formulas (4.16) and choose by φ h = 4τ δϑ n+1 h = 4τ ( δϑ n+1 − δη n+1 ) as follows: , We treat the A 8 term first by rewriting it as follows: In estimating convection terms, we will use Lemma 2.6 frequently without notice. We recall θ(t) 2 ≤ C and E n 0 ≤ C (τ + h) which is the result in Lemma 4.1 to obtain = 0 which comes from (2.7), we rewrite A 8,2 + A 8,4 as Because we can readily get E n+1 h − U n+1 h L 3 (Ω) ≤ M via Lemmas 2.1 and 4.1, we conclude, in conjunction with (4.7), Before we attack B 4 , we note that the assumption τ ≤ Ch is required to apply δE n h L 3 (Ω) ≤ Cτ − 1 2 δE n h 0 . We now imply ϑ n 2 0 + τ ϑ n 2 1 ≤ C(τ 2 + h 2 ) and (4.3) to get (4.26) In light of Hölder inequality, A 9 becomes Inserting the above estimates into (4.25) yields Adding 2 equations (4.24) and (4.27) and then summing up n from 2 to N lead up to(4.20) and complete the proof.
We now estimate optimal accuracy for velocity and temperature.
We now estimate all the terms from A 1 to A 5 respectively. The convection term A 1 can be rewritten as follows: and we denote by A 1,i , for i = 1, 2, · · · , 4 the four terms in the right hand side. To estimate convection terms, we will use frequently Lemma 2.6 without notice. Using u(t n+1 ) 2 ≤ M , we can readily get Because ∇ · ( 2u(t n ) − u(t n−1 ) ) = 0 and 2u(t n ) − u(t n−1 ) = 0 on boundary, we can use (2.8) and so we get In light of ε n+1 and so we can conclude which derives from Lemmas 2.1 and 4.1, then we can get, by the help of Lemma 2.6 and Assumption 1, In conjunction with (4.3), we can have The definition E n+1 and Assumption 1 give us If we apply (4.30) again, then we arrive at On the other hand, the truncation error term becomes Invoking v 0 = 0, inserting the above estimates from A 1 and A 5 into (4.29) lead us In the other hand, we choose φ h = 4τ ω n+1 h in (4.16) to obtain (4.32) δ ∇ω n+1 where In order to estimate A 6 , we note first E n+1 L 3 (Ω) ≤ Ch − 1 2 E n+1 0 and ϑ n+1 2 1 ≤ C τ 2 +h 2 τ which come from Lemma 2.1 and Lemma 4.1, respectively. Then Lemma 2.6 and the assumption τ ≤ Ch yield 2 dt, we can readily get Inserting above estimates into (4.32) yields dt.
In conjunction with the discrete Gronwall inequality, adding (4.31) and (4.33) and summing over n from 1 to N leads (4.28).
We now estimate the pressure error in L 2 (0, T ; L 2 (Ω)). This hinges on the error estimates for the time derivative of velocity and temperature of Lemma 4.3.
Proof. We first recall again inf-sup condition in Assumption 2. Consequently, it suffices to estimate e n+1 , ∇· w in terms of ∇w 0 . In conjunction with (1.7), we can rewrite (4.11) as We now proceed to estimate each term A 1 to A 7 separately. We readily obtain Term A 3 and A 4 can be dealt with thanks to the aid of Lemma 2.6 and u(t n+1 ) 2 ≤ M as follows: In light of u n+1 Integrate by parts and Hölder inequality yield On the other hand, we have The new term A 8 can be bound by the Hölder inequality as Inserting the estimates for A 1 to A 7 back into (4.35), and employing the discrete inf-sup condition in Assumption 2, we obtain If we now square it, multiply it by τ , and sum over n from 1 to N , then Lemmas 4.1, 4.3 and 4.5 derive (4.34).

Numerical Experiments
We finally document 3 computational performance of SGUM. The first is to check accuracy and then the next 2 examples are physically relevant numerical simulations, the Benard convection problem and thermal driven cavity flow. We perform the last 2 examples under the same set within [9], but we conclude with different numerical simulation for the second test, the Benard convection problem, to that of [9]. We impose Taylor-Hood (P2 − P1) in all 3 experiments.  Table 1 is error decay with µ = 1 and τ = h. We conclude that the numerical accuracy of SGUM is optimal and consists with the result of Theorem 2.

Example 3: Thermal Driven Cavity Flow.
We consider the thermal driven cavity flow in an enclosed square Ω = [0, 1] 2 , as already studied in several papers [4,10,9]. The experiment is carried out with the same setting as in Gresho-Lee-Chan [4], which is shown in Figure 5. Figure 6 displays the evolution from rest (t=0) to steady state (t=0.2).   Figure 6. Example 3: Time sequence t = 0.003, 0.01, 0.02, 0.04, and 0.1 for the driven cavity. The first two columns are the streamlines and vector fields for velocity, and the third and fourth are the contour lines for pressure and temperature, respectively. The nondimensional parameters are κ = 10 5 , λ = 1, µ = 1, and the discretization parameters are τ = 10 −4 , h = 2 −5 . Note that u max stands for u ∞ .