Journal of Applied Mathematics

Volume 2013 (2013), Article ID 372906, 21 pages

http://dx.doi.org/10.1155/2013/372906

Research Article

Fully Discrete Finite Element Approximation for the Stabilized Gauge-Uzawa Method to Solve the Boussinesq Equations

Jae-Hong Pyo

Department of Mathematics, Kangwon National University, Chuncheon 200-701, Republic of Korea

Received 9 February 2013; Accepted 11 March 2013

Academic Editor: Jinyun Yuan

Copyright © 2013 Jae-Hong Pyo. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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 the work of Pyo, 2013. 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 Bénard convection problem and the thermal driven cavity flow.

1. 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 with or ,

with initial conditions ,   in , vanishing Dirichlet boundary conditions ,   on , and pressure mean-value . The forcing functions and are given, and is the vector of gravitational acceleration. The nondimensional numbers and are reciprocal of the Reynolds and Prandtl numbers, respectively, whereas is the Grashof number. The Boussinesq system (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 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 [1] to solve Navier-Stokes equations and extended to more complicated problems, the Boussinesq equations in [2] and the nonconstant density Navier-Stokes equations in [3]. 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 [4] and proved superiority for accuracy on the normal mode space, but we could not 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 [5] which is unconditionally stable for semidiscrete level to solve the Navier-Stokes equations. The goal of this paper is to extend SGUM to the Boussinesq equations (1), which model the motion of an incompressible viscous fluid due to thermal effects [6, 7]. 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 computes first an artificial velocity and then decomposes 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 (12) can not be satisfied any more. In order to solve this difficulty, we impose the discontinuous velocity on across interelement boundaries to make fulfill discrete divergence free velocity (12) automatically. We will discuss this issue at Remark 2 below. This discontinuity makes it difficult to treat nonlinear term and to apply the integration by parts, because the discontinuous solution is not included in . 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 Bénard convection problem with the same setting in [2]. In this performance, we newly find out that the number of circulations depends on the time step size . We obtain similar simulation within [2] for a relatively larger , but the behavior is changed for the small . So we conclude that the numerical result of the Bénard convection problem in [2] is not an eventual solution and a more smaller time marching step is required to get the desired simulation.

We will denote as the time marching size. Also we will use as the difference of 2 consecutive functions for example, for any sequence function ,

In order to introduce the finite element discretization, we need further notations. Let be the Sobolev space with derivatives in , set and , where or , and denote by the subspace of of functions with vanishing mean value. We indicate with the norm in and with the inner product in . Let be a shape-regular quasiuniform partition of of mesh size into closed elements [8–10]. The vector and scalar finite element spaces are

where , , and are spaces of polynomials with degree bounded uniformly with respect to [9, 10]. We stress that the space is composed of continuous functions to ensure the crucial equality

Using the following discrete counterpart of the form

we now introduce the fully discrete SGUM to solve the evolution Boussinesq equations (1).

Algorithm 1 (the fully discrete stabilized Gauge-Uzawa method). Compute , , and via any method satisfying Assumption 4 below, and set and . Repeat for .

Step 1. Set and , and then find as the solution of

Step 2. Find as the solution of

Step 3. Update and according to

Step 4. Update pressure by

Step 5. Find as the solution of

Remark 2 (discontinuity of ). We note that is a discontinuous function across inter-element boundaries and that, in light of (7) and (8), 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 Section 2 and introduce well-known lemmas. In Section 3, we prove the stability result.

Theorem 3 (stability). The SGUM is unconditionally stable in the sense that, for all , the following a priori bound holds:

We then will carry out the following optimal error estimates through several lemmas in Section 4.

Theorem 4 (error estimates). Suppose the exact solution of (1) is smooth enough and . If Assumptions 1, 3, 4, and 5 mentioned later hold, then the errors of Algorithm 1 will be bound by

Moreover, if Assumption 2 also hold, then one has

We note that the condition in Theorem 4 can be omitted for the linearized Boussinesq equations (see Remark 16). Finally, we perform numerical tests in Section 5 to check accuracy and physically relevant numerical simulations, the Bénard convection problem and the thermal driven cavity flow.

2. 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 , as well as Poisson’s equation

with boundary condition .

We now state a basic assumption about .

Assumption 1 (regularity of and ). The unique solutions of (16) and of (17) satisfy

We remark that the validity of Assumption 1 is known if is of class [11, 12], or if   is a two-dimensional convex polygon [13] and is generally believed for convex polyhedral [12].

We impose the following properties for relations between the spaces and .

Assumption 2 (discrete inf-sup). There exists a constant such that

Assumption 3 (shape regularity and quasiuniformity [8–10]). There exists a constant such that the ratio between the diameter of an element and the diameter of the largest ball contained in is bounded uniformly by , and is comparable with the meshsize for all .

In order to launch Algorithm 1, we need to set via any first-order methods which hold the following conditions.

Assumption 4 (the setting of the first step values). Let be the exact solution of (1) at . The first step value satisfies

Assumption 5 (approximability [8–10]). For each , there exist approximations such that

The following elementary but crucial relations are derived in [14].

Lemma 5 (inverse inequality). If denotes the Clément interpolant, then

Lemma 6 (div-grad relation). If  , then

Let now indicate the finite element solution of (16); namely,

Then we can find the well-known lemmas in [8–10].

Lemma 7 (error estimates for mixed FEM). Let be the solutions of (16), and, be the Stokes projection defined by (24), respectively. If Assumptions 1, 2, 3, and 5 hold, then

Proof. Inequality (25) is standard [8–10]. To establish (26), we just deal with the -norm, since the other can be treated similarly. If denotes the Clément interpolant, then and

as a consequence of an inverse estimate. This completes the proof.

Remark 8 ( stability of ). The bound is a simple by-product of (16). To see this, we add and subtract , use the stability of in , and observe that (25) implies .

Lemma 9 (error estimates for Poisson's equation [8, 10]). Let be the solution of Poisson's equation (17), and, be its finite element approximation

If Assumptions 1, 3, and 5 hold, then there exists a positive constant satisfying

We finally state without proof several properties of the nonlinear form . In view of (5), we have the following properties of for all :

Applying Sobolev imbedding lemma yields the following useful results.

Lemma 10 (bounds on nonlinear convection [1, 12]). Let with , and let . Then

In addition,

Remark 11 (Treatment of convection term). As we mention in Remark 2, is discontinuous across inter-element boundaries and so . Thus, we can't directly apply (32) anymore to treat the convection term. To solve this difficulty, we apply (34) together with (26) and inverse inequality Lemma 5.

We will use the following algebraic identities frequently to treat time derivative terms.

Lemma 12 (inner product of time derivative terms). For any sequence , one has

3. 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 (6) by using (8) and (10) as follows:

We now choose and apply (35) to obtain

where

We give thanks to (30) for eliminating the convection term. In light of , (12) and (36) yield

Before we estimate , we evaluate an inequality via choosing in (9) to get

Lemma 6 derives , and so (9) and (37) lead us to

Clearly, we have

We now attack with which comes from (8) and (12). Then we arrive at

In conjunction with which comes from (30), if we choose in (11) and use (35), then we obtain

Inserting back into (39), adding with (46), and then summing over from to lead to (13) by help of discrete Gronwall inequality and the equality . So we finish the proof of Theorem 3.

4. Error Estimates

We prove Theorem 4 which is error estimates for SGUM of Algorithm 1. This proof is carried out through several lemmas. We start to prove this theorem with defining to be the Stokes projection of the true solution at time . It means that is the solution of

We also define as the solution of

And we denote notations ,

From Lemma 7, we can deduce

In conjunction with the definition of in (26), we can derive

We now carry out error evaluation by comparing (65) with (6)–(10) and then by comparing (82) and (11). We derive strong estimates of order , and this result is instrumental in proving weak estimates of order for the errors

Then, in conjunction with (50), we can readily get the same accuracy for the errors

Additionally, we denote

We readily obtain the following properties:

Moreover, from (12),

Whence we deduce crucial orthogonality properties:

We also point out that, owing to Lemma 6, defined in (9) satisfies

In conjunction with , Assumption 5 leads to

We now estimate the first-order accuracy for velocity and temperature in Lemma 13 and then the 2nd-order accuracy for time derivative of velocity and temperature in Lemma 15. The result of Lemma 13 is instrumental to treat the convection term in proof of Lemma 15. We will use Lemmas 13 and 15 to prove optimal error accuracy in Lemma 17. Finally, we will prove pressure error estimate in Lemma 18.

Lemma 13 (reduced rate of convergence for velocity and temperature). Suppose the exact solution of (1) is smooth enough. If Assumptions 3–5 hold, then the velocity and temperature error functions satisfy

Proof. We resort to the Taylor theorem to write (1) as follows:

where

are the truncation errors. In conjunction with the definition of the Stokes projection , we readily get a weak formulation for (62), :

We replace the pressure term in (6) with (10) and then subtract it from (65) to obtain

Choosing in (66) and using (58) and (35) yield

where

We now estimate all the terms from to , respectively.  To tackle , we first add and subtract to obtain

Because of , which comes from (30), the second term of can be replaced by

If we apply Lemma 10, then we can readily obtain

Since we have according to (52), we arrive at

In light of (51) and (58), becomes

In order to estimate and , we note that the cancellation law (12) gives . Then yields

In conjunction with the definition , can be evaluated by

If we now apply inequality , then we can get

So we arrive at

In light of and (37), (9) and (59) yield

Also we readily get

The Hölder inequality and (60) yield

Inserting the above estimates into (67) gives

On the other hand, the definition (48) of leads to a weak formula of (63) as, for all ,

We now subtract (11) from (82) to derive

Choosing yields

where

To estimate , we note which comes from (30), then and yield

We use which is (60) to attack :

Inserting the above estimates into (84) gives

Adding 2 inequalities (81) and (88) and summing over from to imply

Because of , we have , and thus, Assumption 4 yields and the first 4 terms in the right-hand side can be bound by Assumption 4 and properties and which are directly deduced from the conditions in Algorithm 1. Moreover, the remaining terms can be treated by the discrete Gronwall lemma. Finally, in conjunction with (58), we conclude the desired result and complete this proof.

Remark 14 (optimal estimation). In order to get optimal accuracy, we must get rid of the terms of and by applying duality argument in Lemma 17. To do this, we first evaluate the errors for time derivative of velocity and temperature in Lemma 15. Thus, we need to evaluate optimal initial errors for the case , and so we have to recompute again (77). We start to rewrite as

In light of Assumption 4, we arrive at

We now start to estimate errors for time derivative of velocity.

Lemma 15 (error estimate for time derivative of velocity and temperature). Suppose the exact solution of (1) is smooth enough and . If Assumptions 3–5 hold, then the time derivative velocity and temperature error functions satisfy

Remark 16 (the condition ). The assumption requires to control convection terms which are used at only (100) and (117), so we can omit this condition for the linearized Boussinesq equations. However, this condition cannot be removed for nonlinear equation case, because (99) must be bounded by .

Proof. We subtract two consecutive formulas of (66) and impose to obtain

where

We now estimate each term from to   separately. The convection term can be rewritten as follows:

In estimating convection terms, we will use Lemma 10 frequently without notice. We recall to obtain

The result in Lemma 13, , is essential to treat the next 2 convection terms. Invoking (33), we have

We note which comes from (30). Then we obtain

To attack , we first note Lemma 5 which is, for any , . If we apply

which is the result of Lemma 13, then we can conclude, in light of (34),

We now estimate using :

In light of Hölder inequality, (51) yields

Integral by parts leads to

In order to tackle , marking use of , we readily get

If we now apply inequality , then we can have

So we arrive at

Invoking (9), (37) and (60) lead to

Finally, the term becomes

For the new term , we argue as in Lemma 13 to arrive at

We now insert the above estimates into (93) to obtain

To evaluate errors of , we subtract two consecutive formulas (83) and choose by as follows:

where

We treat the term first by rewriting it as follows:

In estimating convection terms, we will use Lemma 10 frequently without notice. We recall and which is the result in Lemma 13 to obtain

In conjunction with which comes from (30), we rewrite as

Because we can readily get via Lemmas 5 and 13, we conclude, in conjunction with (60),

Before we attack , we note that the assumption is required to apply . We now imply and (51) to get

In light of Hölder inequality, becomes

Inserting the above estimates into (111) yields

Adding 2 equations (110) and (119) and then summing up from to lead up to (92) and complete the proof.

We now estimate optimal accuracy for velocity and temperature.

Lemma 17 (full rate of convergence for velocity and temperature). One denotes that and are the solutions of (16) and (24) with , respectively. And let and be solutions of (17) and (28) with . Let the exact solution of (1) be smooth enough and . If Assumptions 1 and 3–5 hold, then one has

Proof. We choose in (66) to obtain

where

We now estimate all the terms from to , respectively. The convection term can be rewritten as follows:

and we denote by , for , the four terms in the right-hand side. To estimate convection terms, we will use frequently Lemma 10 without notice. Using , we can readily get

Because and on boundary, we can use (32) and so we get

In light of , we can obtain

and so we can conclude

If we apply which can be derived from Lemmas 5 and 13, then we can get, by the help of Lemma 10 and Assumption 1,

In conjunction with (51), we can have

The definition and Assumption 1 give us

If we apply (127) again, then we arrive at

On the other hand, the truncation error term becomes

For the new term , we observe ; whence

Invoking and inserting the above estimates from and into (121) lead us to

On the other hand, we choose in (83) to obtain

where

In order to estimate , we note first and which come from Lemmas 5 and 13, respectively. Then Lemma 10 and the assumption yield

Since , we can readily get

Inserting the above estimates into (136) yields

In conjunction with the discrete Gronwall inequality, adding (135) and (140) and summing over from to lead to (120).

We now estimate the pressure error in . This hinges on the error estimates for the time derivative of velocity and temperature of Lemma 15.

Lemma 18 (pressure error estimate). Let the exact solution of (1) be smooth enough and . If Assumptions 1–5 hold, then one has

Proof. We first recall again inf-sup condition in Assumption 2. Consequently, it suffices to estimate in terms of . In conjunction with (10), we can rewrite (66) as

We now proceed to estimate each term from to separately. We readily obtain

Terms and can be dealt with thanks to the aid of Lemma 10 and as follows:

In light of from Lemma 13, we can have

Integrate by parts and Hölder inequality yield

On the other hand, we have

The new term can be bound by the Hölder inequality as

Inserting the estimates from to back into (142) and employing the discrete inf-sup condition in Assumption 2, we obtain

If we now square it, multiply it by , and sum over from to , then Lemmas 13, 15, and 17 derive (141).

5. Numerical Experiments

We finally document 3 computational performances of SGUM. The first is to check accuracy, and then the next 2 examples are physically relevant numerical simulations, the Bénard convection problem and the thermal driven cavity flow. We perform the last 2 examples under the same set within [2], but we conclude with different numerical simulation for the second test, the Bénard convection problem, from that of [2]. We impose Taylor-Hood () in all 3 experiments.

Example 19 (mesh analysis). In this first experiment, we choose square domain and impose forcing term the exact solution to become

Table 1 is error decay with and . We conclude that the numerical accuracy of SGUM is optimal and consists with the result of Theorem 4.

Table 1: Error decay for Algorithm 1 with and .

1/16
1/32
1/64
1/128
1/256

0.00198001
0.000651734
0.000188113

Order
1.603153
1.792684
1.888075
1.937437

0.00837795
0.00282504
0.000825997
0.000225008

Order
1.568326
1.774063
1.876160
1.930060

0.0279984
0.00950852
0.00280233
0.00077052
0.000205078

Order
1.558052
1.762594
1.862723
1.909660

0.0341734
0.0123334
0.00376493
0.00105695
0.000283961

Order
1.470303
1.711876
1.832716
1.896142

0.197857
0.107171
0.0453294
0.0170497
0.00610541

Order
0.884544
1.241396
1.410701
1.481586

0.000215482

Order
2.000522
2.000166
2.000081
2.000630

0.000197687

Order
1.987613
1.995466
1.997187
1.996959

0.0167834
0.00839145
0.0041957
0.00209784
0.00104892

Order
1.000043
1.000009
1.000007
1.000000

Example 20 (Bénard convection). In order to explore the applicability of the SGUM, we consider the Bénard convection on the domain with forcing and . Figure 1 displays the initial and boundary conditions for velocity and temperature , as already studied in [2]. Figures 2–4 are simulations at with the nondimensional parameters , and  . Figure 2 is the result for the case which is the same condition in [2], and so it displays similar behavior within [2] including 6 circulations in the velocity stream line. However, Figures 3 and 4, the higher resolution simulations with , and , display 8 circulations in the stream line. So we conclude that the high resolution result is correct simulation and thus Figure 2 and the result in [2] are not eventual simulation.

Figure 1: Example 20: initial and boundary values of Bénard convection problem.

Figure 2: Example 20: streamlines of velocity and isolines of temperature and pressure, at time . The nondimensional parameters are ,  ,  , and the discretization parameters are ,  .

Figure 3: Example 20: streamlines of velocity and isolines of temperature and pressure, at time . The nondimensional parameters are , and , and the discretization parameters are .

Figure 4: Example 20: streamlines of velocity and isolines of temperature and pressure, at time . The nondimensional parameters are , and , and the discretization parameters are .

Example 21 (thermal driven cavity flow). We consider the thermal driven cavity flow in an enclosed square , as already studied in several papers [2, 6, 7]. The experiment is carried out with the same setting as in the work of Gresho et al. [6], which is shown in Figure 5. Figure 6 displays the evolution from rest to steady state .

Figure 5: Example 21: initial and boundary values for thermal driven cavity flow.

(a)

(b)

(c)

(d)

Figure 6: Example 21: time sequence , , , , and for the driven cavity. The first two columns are the streamlines and vector fields for velocity, and the third and fourth ones are the contour lines for pressure and temperature, respectively. The nondimensional parameters are , and , and the discretization parameters are . Note that stands for .

References

R. H. Nochetto and J.-H. Pyo, “The gauge-Uzawa finite element method. I. The Navier-Stokes equations,” SIAM Journal on Numerical Analysis, vol. 43, no. 3, pp. 1043–1068, 2005.

R. H. Nochetto and J.-H. Pyo, “The gauge-Uzawa finite element method. II. The Boussinesq equations,” Mathematical Models & Methods in Applied Sciences, vol. 16, no. 10, pp. 1599–1626, 2006.

J.-H. Pyo and J. Shen, “Gauge-Uzawa methods for incompressible flows with variable density,” Journal of Computational Physics, vol. 221, no. 1, pp. 181–197, 2007.

J.-H. Pyo and J. Shen, “Normal mode analysis of second-order projection methods for incompressible flows,” Discrete and Continuous Dynamical Systems B, vol. 5, no. 3, pp. 817–840, 2005.

J. H. Pyo, “Error estimates for the second order semi-discrete stabilized Gauge-Uzawa methodFor the Navier-Stokes equations,” International Journal of Numerical Analysis & Modeling, vol. 10, pp. 24–41, 2013.

P. M. Gresho, R. L. Lee, S. T. Chan, and R. L. Sani, “Solution of the time-dependent incompressible Navier-Stokes and Boussinesq equations using the Galerkin finite element method,” in Approximation Methods for Navier-Stokes Problems, vol. 771 of Lecture Notes in Mathematics, pp. 203–222, Springer, Berlin, Germany, 1980.

K. Onishi, T. Kuroki, and N. Tosaka, “Further development of BEM in thermal fluid dynamics,” in Boundary Element Methods in Nonlinear Fluid Dynamics, vol. 6 of Developments in Boundary Element Methods, pp. 319–345, 1990.

S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, 1994.

F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, 1991.

V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer, 1986.

P. Constantin and C. Foias, Navier-Stokes Equations, The University of Chicago Press, Chicago, Ill, USA, 1988.

J. G. Heywood and R. Rannacher, “Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization,” SIAM Journal on Numerical Analysis, vol. 19, no. 2, pp. 275–311, 1982.

R. B. Kellogg and J. E. Osborn, “A regularity result for the Stokes problem in a convex polygon,” Journal of Functional Analysis, vol. 21, no. 4, pp. 397–431, 1976.

R. Temam, Navier-Stokes Equations, AMS Chelsea, 2001.