Science, Faculty of
Mathematics, Department of
DSpace
UBCV
Qin, Tong
Master of Science - MSc
University of British Columbia
In this thesis the exactly divergence-free finite element method developed by Cockburn, Scheotzau and Kanschat in [13] and [14] is studied. This method is first reviewed in the context of Stokes problem. An interior penalty discontinuous Galerkin approach is formulated and analysed in the unified framework established in [13], [14] and [37]. Then we extend the method to non-isothermal flow problems, in particular, to a generalised Boussinesq equation. Following the work by Ricardo, Scheotzau and Qin in [34], the method is formulated and the numerical analysis is reviewed. Numerical examples are implemented and presented, which verify the theoretical error estimates and the exactly divergence-free property.
An Exactly Divergence-free Finite Element Method forNon-isothermal Flow ProblemsbyTong QinB.S. University of Science and Technology of China, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Mathematics)The University of British Columbia(Vancouver)August 2013c? Tong Qin, 2013AbstractIn this thesis the exactly divergence-free finite element method developed in [13] and [14]is studied. This method is first reviewed in the context of Stokes problem. An interiorpenalty discontinuous Galerkin approach is formulated and analysed in the framework of[13], [14] and [37]. Then we extend the method to non-isothermal flow problems, in partic-ular, to a generalised Boussinesq equation. Following [34], the method is formulated andthe numerical analysis is reviewed. Numerical examples are implemented and presented,which verify the theoretical error estimates and the exactly divergence-free property.iiPrefaceIn this dissertation, the theoretical analysis in Chapter 1 is a review and synthesis ofthe work in [13], [14] and [37]. And the numerical experiments in Section 2.5 are carriedout independently by the author T. Qin. The Chapter 2 is a summary and review of thepaper [34] which is a collaboration with R. Oyarzu?a and D. Scho?tzau. The formulationof the method is an extension of the work in [14]. The theoretical analysis is attributedto Oyarzu?a and Scho?tzau. The author formulated the fixed-point iteration in Section 3.7and extended the code to the generalized Boussinesq equation in Section 3.8. All thenumerical experiments in this thesis are conducted with the deal.II finite element library[5]. A journal paper related with the work in Chapter 2 has been submitted to andaccepted by the IMA Journal of Applied Mathematics.iiiContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Divergence-Free Elements for the Stokes Problem . . . . . . . . . . . . 32.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Preliminaries and Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Divergence-Free Elements for the Stokes Problem . . . . . . . . . . . . . . . 62.3.1 The Stokes Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3.2 Discrete Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3.3 Formulation of the interior penalty DG Method . . . . . . . . . . . . 102.4 Theoretical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4.1 Stability Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.4.2 Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.5 Numerical Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Non-isothermal Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Physical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.2.1 Model I: Non-Isothermal Stokes-Oldroyd System . . . . . . . . . . . 233.2.2 Model II: Buoyancy Driven Flows . . . . . . . . . . . . . . . . . . . 25iv3.3 Model Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3.1 Stability Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3.2 Unique Solvability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.4 Discrete Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.5 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.5.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.5.2 Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5.3 Coercivity and inf-sup condition . . . . . . . . . . . . . . . . . . . . 333.6 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.6.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.6.2 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.7 Fixed-point Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.8 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.8.1 Smooth Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.9 Conclusions and future directions . . . . . . . . . . . . . . . . . . . . . . . . 41Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42vList of Tables2.1 Numerical results for Kovasznay flow (? = 1) . . . . . . . . . . . . . . . . . 193.1 Numerical results I: IPDG for the non-isothermal problem (3.7) with fixed-point iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 Numerical results II: IPDG for the non-isothermal problem (3.7) with fixed-point iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40viList of Figures2.1 The illustration for the distribution of dofs of the RT0 (left) and RT1(right) elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 The convergent plots for the numerical errors . . . . . . . . . . . . . . . . . 203.1 The convergent plots for the numerical errors . . . . . . . . . . . . . . . . . 393.2 The convergent plots for the numerical errors . . . . . . . . . . . . . . . . . 393.3 The residual plot for k = 1, l = 1, ..., 4 . . . . . . . . . . . . . . . . . . . . . 40viiAcknowledgementsFirst, I would like to express my sincere gratitude to my supervisor, Prof. DominikScho?tzau, for his thoughtful guidance, consistent encouragement and patient discussionthroughout my graduate study and the writing of this thesis. I could not have imaginedhaving a better advisor and mentor for my master?s study.I would also like to express my gratitude to my friends, office mates and room mates.They not only bring me help in both work and life but also make my life colourful. Spe-cial thanks to my girlfriend Yang Chen, without whose love and insistence, I could notget this far.Last but not least, I want to express my sincere appreciation to all the helpful staff of theMathematics Department and the Institute of Applied Mathematics at UBC, for provid-ing excellent seminars and talks to enjoy and for providing a comfortable environment towork in. They really make me feel like a member of a big family!Tong QinThe University of British ColumbiaAugust 2013viiiChapter 1IntroductionWhen simulating incompressible flow with finite element methods, the exact satisfactionof the incompressibility constraint ? ?u = 0 is highly desirable. It yields provable energy-stable methods without the need to modify nonlinear terms in the underlying differentialequations as is done in the classical methods for the incompressible Navier-Stokes equa-tions in [40] and [41]. In contrast, failure to satisfy the constraint exactly can causeundesired instabilities especially in nonlinear problems, for example in collision flows incross-shaped domains [32], and in the context of rotation-free forces acting on elasticmaterials [4]. In recent decades, many efforts have been made in designing efficient meth-ods that are exactly divergence-free for incompressible flow problems in both conformingand the nonconforming setting. We mention here only [39, 42, 9, 22] and the referencestherein. In particular, in 2007, Cockburn, Kanschat and Scho?tzau developed and anal-ysed a new exactly divergence-free element method for the incompressible Navier-Stokesproblem, which is based on using divergence-conforming elements for the velocities anddiscontinuous elements for the pressure. The local approximation space pair V ? Q ischosen to satisfy the inclusion ? ? V ? Q. As shown in [14] and motivated by [13],this method is proved to be energy-stable, optimally convergent, locally conservative andyields exactly divergence-free velocity approximations. In [21], this method was furtherextended to an incompressible magnetohydrodynamics problem. In the first part of thisthesis, this method is reviewed in the context of the analysis of the Stokes problem forquadrilateral meshes. Specifically, we follow [14] and employ the interior penalty discon-tinuous Galerkin (IPDG) method [1] to enforce the H1-conformity of the velocity field.The local approximation space pair for the velocity and pressure is chosen to be theRaviart-Thomas element and the tensor product polynomial spaces (these are defined inSection 2.2) respectively, with equal approximation order k. The a priori error estimateis reviewed under the framework of [37] and is tested by numerical experiments.The second part is devoted to extending this exactly divergence-free method to non-1isothermal flow problems. The simulation of non-isothermal flows is receiving more andmore attention in many scientific and engineering branches such as the simulation of theconvection process of the mantle in geophysical studies [29], the simulation of coolingprocesses in the device design and that of the polymer melting processes in industrialpractice [26, 30], to name a few. Due to the importance of the model, numerical meth-ods for this type of problem have been actively developed, see [15, 19, 17, 18] and thereferences therein. But none of them share the property of being exactly divergence-freein contrast to the method under consideration. In the second part of this thesis, we fol-low the work in [34] and consider the generalised Boussinesq equation, which couples thesteady incompressible Navier-Stokes equation and a convection-diffusion equation througha nonlinear temperature-dependent viscosity and through a buoyancy term acting in theopposite direction of the gravity. The nonlinearity is addressed with a fixed-point ap-proach. We extend the IPDG method from the first part as well as follow the work in[14] to discretize the linearized incompressible Navier-Stokes equation in each iteration.For the convection-diffusion equation, the standard conforming method is employed. Thecrucial aspect in the theoretical analysis is the construction of a suitable lifting of the tem-perature boundary data into the computational domain. A more restrictive assumptionhas to be made on the size of the solution. In the end of both parts, numerical experimentsare conducted with the deal.II finite element library with the goal of verifying theoreticalerror estimates and the exactly divergence-free property.2Chapter 2Divergence-Free Elements for theStokes Problem2.1 IntroductionIn this chapter we recall and review the exactly divergence-free elements developed in [13]and [14] for the Stokes problem of incompressible fluid flow:???u+?p = f , in ? (2.1a)? ? u = 0, in ? (2.1b)u = g, on ? = ??, (2.1c)where ? is the computational domain in R2 and ? > 0 denotes the kinematic viscosity,u the velocity, p the pressure and f ? (L2(?))2 the unit external volumetric force actingon the fluid. Furthermore the boundary datum g ? H1/2(?)2 satisfies the compatibilitycondition:??g ? n ds = 0. (2.2)The exact satisfaction of the divergence constraint (2.1b) is a highly desirable propertyin the simulation of incompressible flows. It yields energy-stable methods without theneed to modify nonlinear terms in the underlying differential equations as is done in thecontext of the incompressible Navier-Stokes equations in [40] and [41]. Not enforcing thisconstraint exactly may result in instabilities when simulating the incompressible steadyflow problems, especially when a rotation-free part occurs in the external force vectoror in a nonlinear convection term; see [32]. However, most of the stable finite elementpairs, such as the Taylor-Hood element [24] and the MINI element [2] do not satisfy thedivergence constraint exactly.3In the context of conforming finite element methods, in 1985, Scott and Vogelius [39]constructed inf-sup stable and exact divergence-free pairs by using continuous Pk ele-ments, which are polynomials of maximal degree k, for the approximation of the velocityand the discontinuous Pk?1 elements for the pressure. This pair is also employed in [32]as a comparison with the Taylor-Hood element to illustrate the importance of mass con-servation in incompressible fluid simulation. Since ? ?Vh ? Qh, it is obvious to see theexact divergence-free property. But this pair is proved to be stable on triangular gridswithout singular vertices with the polynomial degree satisfying k ? 4. For smaller k, themeshes need to be Hsieh-Clough-Tocher or Powell-Sabin triangulations; see the review pa-per by Zhang [42]. Recently, Guzma?n and Neilan overcame this difficulty and constructedfinite element pairs for the Stokes problem on general triangular meshes that are inf-supstable, optimally convergent and divergence-free; see [22] for 2D and [23] for 3D. Theidea is to enrich the H(div; ?)-conforming elements locally with divergence-free rationalshape-functions. These additional functions will help to enforce the continuity. For ap-proaches with exactly divergence-free velocity approximation in the isogeometric analysiswe refer to [9] which is based on a B-spline generalization of the Raviart-Thomas elementsand to [4] for the application in the context of elasticity of nearly incompressible materials.In this chapter, we will use a discontinuous Galerkin approach to exactly enforce the con-straint (2.1b). This approach for Navier-Stokes equations was first introduced in [13] and[14] in the context of local discontinuous Galerkin (LDG) methods. The resulting methodsare locally conservative, energy-stable, optimally convergent and exactly divergence-free.As pointed in [13], we adopt here the interior penalty approach for the viscous term. Inthis chapter, we review the formulation and analysis of this method by solving the modelStokes problem with Raviart-Thomas elements of order k for the velocity approximationand tensor product polynomials of degree k for the pressure.2.2 Preliminaries and NotationsIn order to formulate and analyse the DG methods, we introduce the following notationsthat will be used throughout this thesis.Let ? denote a polygonal domain in R2 with Lipschitz boundary ? = ??. Supposewe have a quasi-uniform and regular (see [11]) triangulation Th of the domain ? into ele-ments {K}. We denote the elemental diameter by hK and, as usual, define the meshsizeh = maxK?Th hK . We further denote by EIh the set of all interior edges of Th and by EBhthe set of all boundary edges. We set Eh = EIh ? EBh .4The standard definitions are used for the Lebesgue and Sobolev spaces Lp(?) andW r,p(?),where r ? 0 and p ? [1,?], with the norms ? ? ?Lp(?) and ? ? ?W r,p(?), respectively. Wewrite Hr(?) instead of W r,2(?) and denote the associated seminorms by | ? |r,?. Forvector-valued function, we denote Hr(?) = Hr(?)2 and Lr(?) = Lr(?)2. Moreover, wewill use the vector-valued Hilbert spacesH(div; ?) = {v ? (L2(?))2,? ? v ? L2(?)}with the norm?v?2div,? = ?v?20,? + ?? ? v?20,?On the reference element K?, for nonnegative integers k and l, the spaces Qk,l(K?) consistof polynomials of degree at most k and l in the first and second variables, respectively.When k = l, the notation Qk(K?) is adopted instead of Qk,k(K?). And we use Pk(K) andPk(e) to denote the space of polynomials of maximal degree k defined in any element Kand any edge e, respectively. We assume that each element K ? Th is the image of thereference element K? = (?1, 1)2 under the invertible affine mapping FK : K? ? K, withFK(x?) = BK x?+ bK , BK ? R2?2, bK ? R2.To define the DG method, we introduce notations associated with traces. Let K+ andK? be two adjacent elements of Th. Let x be an arbitrary point of the interior edgee = ?K+ ? ?K? ? EIh . Let ? be a piecewise smooth scalar-, vector-, or matrix-valuedfunction and let us denote by ?? the traces of ? on e taken from within the interior ofK?. Then, we define the mean value {{?}} at x ? e as{{?}} := 12(?+ + ??).where ? can be any piecewise smooth scalar-, vector- or tensor-valued function. Further,we define the jumps [[?]] at x ? e asfor e ? EIh :q := (q+ ? q?)nK+, v := (v+ ? v?)? nK+, v := (v+ ? v?) ? nK+,for e ? EBh :q := qn, v := v ? n, v := v ? n,Here, q is any scalar-valued function, v is a vector-valued one and nK denotes an outwardunit normal vector on the boundary ?K of element K and (v ? n)ij = vinj .5We end this section by a summation identity from [3], which will be used to transformthe local formulation the DG methods to the global one:?K?Th??Kqv ? n ds =?Eh{{v}} ? q ds+?EIh{{q}} v ds (2.3)2.3 Divergence-Free Elements for the Stokes Problem2.3.1 The Stokes ProblemFor simplicity and without loss of generality, from this section on, we take ? = 1. Weintroduce the solution spaces V and Q which are defined to beV = H1(?), Q = L20(?) ={q ? L2(?) :??q dx = 0}(2.4)The weak formulation of (2.1) is to find (u, p) ? V ?Q with u|? = g such that{A(u,v) +B(v, q) = (f ,v)B(u, q) = 0(2.5)for all (v, q) ? H10(?)? L20(?) where the bilinear forms A(?, ?) and B(?, ?) are defined asA(u,v) =???u : ?v dx, B(u, p) = ???p? ? udx (2.6)The well-posedness of (2.5) is ensured by the continuity of A(?, ?) and B(?, ?), the coercivityof A(?, ?), and the following inf-sup conditioninf0=q?L20(?)sup0 =v?H10 (?)d??? q? ? v|v|1?q?0? ? > 0 (2.7)with an inf-sup constant ? only depending on ?; see [8].2.3.2 Discrete SpacesWe will use the following discrete spaces when formulating the DG method.Uh = {v ? L2(?)2 : v|K ? Qk+1(K), K ? Th} (2.8)Vh = {v ? H(div; ?) : v|K ? RTk(K), K ? Th; v ? n = gh on ?} (2.9)V0h = {v ? H(div; ?) : v|K ? RTk(K), K ? Th; v ? n = 0 on ?} (2.10)Qh = {q ? Q : q ? Qk(K), K ? Th} (2.11)6+4Figure 2.1: The illustration for the distribution of dofs of the RT0 (left) and RT1 (right)elementswhere gh is the L2 projection of g ? n into the space Pk(e), or equivalently,?eg ? n?ds =?egh?ds, ?? ? Pk(e), ?e ? EBh (2.12)We will see the reason for this definition in following sections, especially in the proofsof Proposition 2.3.3 (the divergence-free property) and Theorem 2.4.7 (the a priori errorestimate).On the reference element K? = (?1, 1)2, the Raviart-Thomas space of order k (RTk)is defined (see [8, Section 3.2]) byRTk(K?) = Qk+1,k(K?)?Qk,k+1(K?)with the corresponding interpolation operator ??k satisfying ??ku ? RTk and??????K?(??ku?? u?) ? v? = 0, ?v? ? Qk?1,k(K?)?Qk,k?1(K?)??m( ??ku?? u?) ? n? q? ds? = 0, ?q? ? Pk(?m), m = 1, ..., 4(2.13)where ?m, m = 1, ..., 4 are the boundary edges of K?. In Figure 2.1, we take the RT0(K?)and RT1(K?) elements as examples to illustrate the distribution of the degrees of freedom(dofs) of the Raviart-Thomas element. The arrows on the boundary represent the dofsrelated with the normal component. And the number in the cell is the number of dofslocated inside the cell.On a general physical element K, it?s easy to check that the affine transformation FK7does not preserve the degrees of freedom??mu ? n q ds, for q ? Pk(?m). The Piola trans-formation, which is defined asPK u? =1detBKBK u? ? F?1K , (2.14)overcomes this difficulty. It was first introduced in the continuum mechanics to map vec-tors between Eulerian and Lagrangian coordinates while preserving the total flux overthe boundary, i.e.,??G? v? ? n? ds? =??G v ? n ds, where G? is a reference volume (Lagrangiancoordinate) and G is the deformed one (Eulerian coordinate); see [12].In general, we have following properties of the Piola transformation.Lemma 2.3.1 For q ? H1(K?), v? ? H(?div, K?) and q := q? ? F?1K , v = P(v?), the Piolatransformation preserves the following quantities:?Kv ? ?q dx =?K?v? ? ??q? dx?, (2.15)?K? ? v q dx =?K??? ? v?q? dx? (2.16)??Kv ? n q ds =??K?v? ? n?q? ds? (2.17)Proof: For the affine transformation, the proof is very simple. We add it here forcompleteness. An application of the chain rule gives the following result? = B?TK ?? (2.18)Then we have?Kv ? ?q dx =?K?(1detBKBK u?)? (B?TK ??q?) detBK dx?=?K?(BK?1BK) : u?? ??q? dx?=?K?u? ? ??q? dx?.For the equation (2.16), first, we notice that?v = B?TK ??(1detBKBK v?)=1detBKB?TK ??v?BTK .8Since the trace is invariant under a similarity transformation, we further have? ? v = 1detBK?? ? v?. (2.19)Then equation (2.16) follows by a simple substitution.At last, it is easy to see that equation (2.17) can be obtained by divergence theorem andequations (2.15) and (2.16).Remark From the proof above we see that in Piola transformation, the factor 1detBK isintroduced to cancel out the scaling factor. And the matrix BK is for balancing the factorB?TK appearing in the gradient operator in (2.18).From the last identity we see that the Piola transformation preserves the normal trancesin H?1/2(?K?), which makes it a suitable transformation for the H(div) spaces. Thereforewe can define the RTk(K) space as PK(RTk(K?)). For more properties of both the RTkelement and the Piola transformation, one can refer to [8].For nonnegative integer k, l, the tensor product polynomial spaces Qk,l(K) on the physicalelement K are defined via the affine mapping FK asQk,l(K) := {q : q(x, y) = q? ? F?1K (x, y), q? ? Qk,l(K)} (2.20)and when k = l, we denote Qk(K) = Qk,k(K). Then it is obvious to get Vh ? Uh andV0h ? Uh.We end this section by an important relationship between RTk(K) and Qk(K), which,we will see, plays a key role in ensuring the exact divergence-free property.Proposition 2.3.2 For affine elements K ? Th, the local approximation space pair RTk(K)and Qk(K) satisfies the following inclusion:? ?RTk(K) ? Qk(K) (2.21)Proof: First on the reference element K, by the definition of RTk(K?) it is obvious tosee that?? ?RTk(K?) = Qk(K?) (2.22)In the proof of the last proposition we have obtained? ? v = 1detBK?? ? v?, (2.23)9for v ? RTk(K). Due to (2.22), we know that there exists a polynomial q? ? Qk(K?) suchthat ?? ? v? = q?. Then by the definition of Qk(K), we have? ? v = 1detBKq? ? F?1K ? Qk(K)which implies the conclusion.2.3.3 Formulation of the interior penalty DG MethodWe denote (uh, ph) ? Vh ? Qh to be the numerical approximation. First, we follow theprocedure in [3] and discretize each component of the viscous term with the symmetricinterior penalty method [1]. Then we get the elliptic formAh(uh,v) =???huh : ?hvdx+?0h?EIhuh : v ds??EIhuh : {{?hv}} ds ??EIhv : {{?huh}} ds+?0h?EBhuh ? v ds ??EBh(uh ? n) : ?hv ds??EBh(v ? n) : ?huh ds (2.24)for v in V0h. The penalty parameter ?0 is chosen sufficiently large but independent of h.The broken gradient operator ?h is defined as ?hv|K = ?v.Pressure GradientFor the pressure gradient, integration by parts over a single element K gives?Kv ? ?phdx =??K(nK ? v)p?h ds??K(? ? v)ph ds, v ? V0h (2.25)The function p?h is the numerical trace, which, as in [13], is taken to bep?h = {{ph}}, for e ? EIh ; p?h = ph, for e ? EBh . (2.26)By summing (2.25) over K ? Th and employing the key summation formula (2.3), we canget??v ? ?hph =?Eh{{v}} ? p?h +?EIhv {{p?h}} ???(?h ? v)ph dx (2.27)Since v belongs to V0h which is H(div) conforming with a zero normal component on theboundary and the numerical trace p?h is single-valued on edges, we can simplify (2.27) into??v ? ?hph = ???ph?h ? v dx (2.28)10and if we defineBh(v, q) := ???(? ? v)q dx, (2.29)we have ??v ? ?ph dx = Bh(v, ph) (2.30)Incompressibility ConstraintBy testing the incompressibility constraint (2.1b) against a smooth function v ? V0h overany element K ? Th and by integrating by parts, we have the flux form as below0 =?Kq? ? uh dx = ??Kuh ? ?q dx+??Ku?ph ? nKq ds (2.31)where u?ph is the numerical trace related with the incompressibility constraint. The fol-lowing proposition states that if we pick the simplest form for u?ph and with the inclusioncondition (2.21), the approximate velocity is exactly divergence-free with the incompress-ibility constraint only weakly enforced.Proposition 2.3.3 Take the numerical traces u?ph to beu?ph = {{uh}} on EIh, u?ph = g on EBh (2.32)Then1. For uh ? Vh and q ? Qh, (2.31) can be rewritten into?Bh(uh, q) = 0 (2.33)2. If uh ? Vh satisfies (2.33) for all q ? Qh, then uh is exactly divergence-free, i.e.? ? uh ? 0, in ? (2.34)Proof: 1) By the summation identity (2.3) we can sum (2.31) over K ? Th and get0 = ???uh ? ?hq dx+?Eh{{u?ph}} ? q +?EIhu?ph{{q}}= ???u ? ?hq dx+?Ehu?ph ? q ds (2.35)11Then a second integration by parts of the first term gives0 = ?(???q? ? uh dx+?Eh{{uh}} ? q +?EIhuh {{q}})+??u?p ? q ds (2.36)=??q? ? uh ??Ehu?p ? q??EIh0 ? {{q}}+?Ehu?ph ? q= ?Bh(uh, q) (2.37)2) The inclusion condition (2.21), the divergence theorem, the essential boundary conditionin Vh and the compatibility condition (2.2) give us? ? uh|K ? Qk(K) (2.38)??? ? uh dx =??uh ? n ds =??g ? n ds = 0. (2.39)These imply that ?h ? uh ? Qh, which combined with the fact that?Bh(uh, q) =??q?h ? uh dx = 0, ?q ? Qhimplies? ? uh ? 0The interior penalty DG method for the Stokes problem now consists in finding (uh, ph) ?Vh ?Qh which satisfiesAh(uh,v) +Bh(v, ph) = Fh(v) (2.40a)?Bh(uh, q) = 0 (2.40b)for all (v, q) ? V0h ?Qh and the linear forms are defined asAh(u,v) =???hu : ?hvdx+?0h?EIhu : v??EIhu : {{?hv}} ??EIhv : {{?hu}}+?0h?EBhu ? v ??EBh(u? n) : ?hv ??EBh(v ? n) : ?hu (2.41)Bh(v, q) = ???q? ? v dx (2.42)Fh(uh) =??f ? vdx ??EBh(g ? n) : ?hv ds+?0h?EBhg ? vds (2.43)122.4 Theoretical AnalysisWe first introduce the full spaceU(h) := H1(?)2 +Uh (2.44)endowed with the DG norm?v?21,h :=?K?Th??v?20,K +?e?Eh?e?0h?1|v|2 ds. (2.45)where ?0 is the penalty parameter independent of h. We also introduce the followingsubsets of U(h):V0(h) := H1(?)2 +V0h (2.46)V(h) := H1(?)2 +Vh (2.47)To facilitate the analysis of the method we need to introduce an auxiliary space ?h definedby?K := {? ? L2(?)2?2 : ? ? Qk+1(K)2?2, K ? Th} (2.48)Obviously, we have ?Uh ? ?h.In [3], lifting operators, which lift the edge terms into the cell, are introduced in or-der to cast all the existing DG methods for elliptic problems into a unified framework. In[37], these operators are further applied to build a unified framework for analysing severalDG methods for the saddle-point problems. The idea is to recast the DG methods into aperturbed form with the help of the lifting operators and then derive the abstract errorestimates which is a variant of the Strang?s lemma. The estimates of jump terms over theedges are converted to those of lifting operators.Here, we follow [37], for e ? Eh we introduce the lifting operator Le : V(h) ? ?hdefined by??Le(v) : ? dx =?ev : {{?}} ds, ?? ? ?h. (2.49)Globally, we define L : V(h) ? ?h asL :=?e?EhLe13Then, we define a perturbed form A?h(?, ?) basing on the bilinear form Ah(?, ?) in (2.41):for u, v ? V(h), V0(h)A?h(u,v) :=??[?hu : ?hv ? L(u) : ?hv ? L(v) : ?hu] dx+ ??Ehu : v ds. (2.50)We notice that when restricted to the discrete space Vh, we haveA?h(u,v) = Ah(u,v), ?u, v ? Vh, V0h.Hence, using the perturbed form does not change the discrete variational problem (2.40).2.4.1 Stability ResultsLet?s first review some useful stability results of the bilinear forms, upon which the proofof the well-posedness and the analysis of the a priori error is based. The first result isrelated to the boundedness and coercivity of the bilinear form Ah(?, ?).Proposition 2.4.1 For the perturbed form A?h(?, ?), there exist constant c1 and c2 suchthat|A?h(u,v)| ? c1?u?1,h?v?1,h, ?u ? V(h), ?v ? V0hA?h(v,v) ? c2?v?21,h, ?u ? V0hProof: Since V0h, Vh ? Uh = {v ? L2(?)2 : v|K ? Qk+1(K),K ? Th}, this propositionis a direct consequence of [37, Lemma 7.5, Lemma 7.6].The next proposition is about the boundedness of Bh(v, q).Proposition 2.4.2 The form Bh is continuous|Bh(v, q)| ??2?v?1,h?q?0 ?(v, q) ? Vh0 ? L2(?) (2.51)and the velocity-pressure pair V0h ?Qh is inf-sup stableinfq?Qhsupv?V0hBh(v, q)?v?1,h?q?0? ? > 0 (2.52)for a constant ? independent of the mesh size.14Proof: For the continuity, by the definition of Bh(?, ?) and ?h and the Cauchy-Schwarzinequality, we have|Bh(v, q)| =???????h ? v q dx????=???????K?Th?K? ? v q????????K?Th?? ? v?0,K?q?0,K????K?Th?? ? v?20,K??12???K?Th?q?20,K??12????K?Th?K2d?i=1(?ivi)2??12?q?0??2???K?Th|?v|21,K??12?q?0??2?v?1,h?q?0The proof of the inf-sup condition can be found in [37, Theorem 6.12].It is easy to verify that the stability results in this section imply the unique solvabil-ity of the approximation problem (2.40).2.4.2 Error EstimatesThe a priori error estimate is based on an abstract error estimate and some standardapproximation results. We make the regularity assumption that the exact solution (u, p)of the Stokes system (2.1) satisfies(u, p) ? Hs+1(Th)2 ?Hs(Th), s ? 1 (2.53)and u|? = g where Hr(Th), for r > 0, is the broken Sobolev space with the definitionHr(Th) = {v ? L2(?) : v|K ? Hr(K), K ? Th} (2.54)Abstract Error EstimateWe follow [37] and [13] and first present the abstract error estimate in terms of the ex-pressionRh(u, p;v) := sup0 =v?V0h|Rh(u, p;v)|?v?1,h(2.55)15where the residual Rh(u, p;v) is defined asRh(u, p;v) = A?h(u,v) +Bh(v, p) ? Fh(v) (2.56)for v ? V0h.We need to emphasize that, due to the introduction of the perturbed form A?h, we can notconclude that Rh(u, p;v) = 0 in general. The reason for this is that the term L(v) : ?huis inconsistent since ?hu ? ?h. In the next proposition, we will see that the residual isactually optimally convergent.Proposition 2.4.3 Assume that the exact solution (u, p) of the Stokes system (2.1) sat-isfy (2.53). Let Q and Q be the L2 projections onto ?h and Qh, respectively. Then theresidual in Rh(u, p;v) is given byRh(u, p;v) =?Eh{{?u?Q(?u)}} : v ds, ?v ? V0h (2.57)Moreover, we have thatRh(u, p)2 ? Ch2min(s,k)?u?2s+1 (2.58)with a constant C > 0 independent of h.Proof: The proof can be found in [37, Proposition 8.1].In the next proposition we present the abstract error estimate which is a variant of Strang?slemma.Proposition 2.4.4 Let (u, p) ? V ? Q be the exact solution of the Stokes problem and(uh, ph) ? Vh ?Qh the approximate solution of the DG problem (2.40). We have?u? uh?1,h ? Cu[infv?Vh?u? v?1,h + infq?Qh?p? q?0 +Rh(u, p)](2.59)?p ? ph?0 ? Cp[infq?Qh?p? q?0 + ?u? uh?1,h +Rh(u, p)](2.60)where Cu, Cp > 0 only depend on the stability parameters.Proof: The proof can be found in [37, Proposition 4.1] and is based on the stabilityresults of the bilinear forms Ah(?, ?) and Bh(?, ?).To prove the a priori estimate we need the following standard results:16Lemma 2.4.5 [8, Proposition 3.6] Let K be an affine element and W (K) = {v ?Lr(K))|? ?v ? L2(?)}, r > 2. The interpolant operator mapping ?K : W (K) ? RTk(K)is defined to be ?Kv = ??kv? ?P?1K where v? = v ?PK . There exists a constant c dependingonly on k and on the shape of K, such that, for any v ? Hs+1(K), we have?v ? ?Kv?0,K ? chmin{s,k}+1K |v|s+1,K (2.61)?v ? ?Kv?1,K ? chmin{s,k}K |v|s+1,K (2.62)The last lemma will help us deal with jump terms in the definition of the ? ? ?1,h, as wewill see in the proof of Theorem 2.4.7.Lemma 2.4.6 [36, Lemma 1.2.14] For any function u ? H1(Th), we have?u?0,?K ? C(h1/2K |u|1,K + h?1/2K ?u?0,K) (2.63)where K is any element in Th and constant C > 0 independent of meshsize.Then based on this we are ready to state the apriori error estimate resultTheorem 2.4.7 (a priori error estimate) Let (u, p) be the exact solution of the Stokesproblem satisfying the regularity assumption (2.53) and (uh, ph) ? Vh ? Qh the finiteelement approximation. We have the following error estimates?u? uh?h ? Chmin{s,k}?u?s+1 (2.64)?p ? ph?0 ? Chmin{s,k}?u?s+1 (2.65)where the constant C > 0 independent of h.Proof: By the standard approximation estimates for the L2-projection (see [20, Ap-pendix A]) and Proposition 2.4.3 we haveinfq?Qh?p? q?0 ? Chmin{s,k+1}?q?s (2.66)Rh(u, p) ? Chmin{s,k}?u?s+1 (2.67)According to the abstract error estimates (2.59) and (2.60), it then remains to bound theinfimum infv?Vh ?u?v?1,h. We take v = ?u, where? is the global interpolation operatormapping from H(div,?) ? (Lr(?))d r > 2, into RTk(?,Th) := {v ? H(div,?), v|K ?RTk(K)} by setting?u|K = ?Ku17where ?K is the local interpolant operator defined in Lemma 2.4.5. By (2.12), the defi-nition of the Raviart-Thomas interpolation operator ?? in (2.13) and the property of thePiola transformation (2.17), it is not hard to see that ?u ? Vh.If we set ?u = u ? ?u, then by Lemma 2.4.5, Cauchy-Schwarz inequality and by ap-plying Lemma 2.4.6 to each component we can get??u?21,h =?K?Th|?u|21,K +?E?| ?u ? n |2 ds??K?Th??u?21,K +?K?Th??K?0h?1K |?u ? n|2 ds? Ch2min{s,k}?u?2s+1 +?K?Th?0h?1K ??u?20,?K? Ch2min{s,k}?u?2s+1 +?K?ThC(hK |?u|21,K + h?1K ??u?20,K)? Ch2min{s,k}?u?2s+1 +?K?ThCh2min{s,k}+1K ?u?2s+1,K? Ch2min{s,k}?u?2s+1which impliesinfv?Vh?u? v?h ? ?e??h ? Chs?u?s+1 (2.68)Then by combining the abstract error estimates and the approximation results (2.67),(2.66) we can get the a priori error estimates (2.64),(2.65).2.5 Numerical ExperimentIn this section numerical experiments are conducted with the deal.II finite element library[5]. We test the method with the Kovasznay Flow [28]; see also the numerical experimentsin [13] and [14]. The computational domain is taken to be? =[?12,32]? [0, 2]and is covered by a uniform rectangular mesh. At level l, the meshsize is h = 2?l. Weprescribe the analytical solutionu1(x, y) = 1? e?x cos(2?y),u2(x, y) =?2?e?x sin(2?y).p(x, y) = ?12e2?x ? p?18where? =?8?2??1 +???2 + 16?2The viscosity ? is taken to be 1 and the corresponding right hand side isf =???(?2 ? 4?2)e?x cos(2?y)? ?e2?x(2??? ?32? )e?x sin(2?y)???In this experiment the penalty parameter is chosen to be?0 = 2k(k + 3)As emphasized in [14], the normal component of the inhomogeneous Dirichlet boundarycondition is interpolated in a strong way on a set of quadrature points on which the RTkmoments are integrated exactly. The L2-norm error for the pressure and the velocity ask cycle ?ep?0 ?eu?0 ?eu?1,h ?? ? uh??k=14 4.543e-01 1.08 8.733e-03 1.99 1.631e+00 1.64 1.705e-135 1.928e-01 1.24 2.227e-03 1.97 6.964e-01 1.23 5.116e-136 7.071e-02 1.45 5.660e-04 1.98 3.254e-01 1.10 1.020e-127 2.339e-02 1.60 1.431e-04 1.98 1.571e-01 1.05 3.197e-12k=21 1.945e+00 - 6.433e-02 - 6.954e+00 - 5.826e-132 3.307e-01 2.56 7.889e-03 3.03 1.225e+00 2.50 1.070e-123 4.675e-02 2.82 9.562e-04 3.04 3.099e-01 1.98 2.727e-124 6.343e-03 2.88 1.166e-04 3.04 7.919e-02 1.97 5.177e-125 8.816e-04 2.85 1.436e-05 3.02 2.005e-02 1.98 1.013e-116 1.290e-04 2.77 1.780e-06 3.01 5.047e-03 1.99 3.980e-11k=31 3.627e-01 - 5.738e-03 - 4.281e-01 - 1.495e-122 3.243e-02 3.48 3.689e-04 3.96 4.260e-02 3.33 4.182e-123 2.424e-03 3.74 2.318e-05 3.99 4.698e-03 3.18 7.731e-124 1.799e-04 3.75 1.450e-06 4.00 5.471e-04 3.10 2.009e-115 1.421e-05 3.66 9.070e-08 4.00 6.481e-05 3.08 5.423e-11k=41 5.201e-02 - 3.950e-04 - 3.961e-02 - 2.087e-122 2.312e-03 4.49 1.219e-05 5.02 2.846e-03 3.80 5.684e-123 8.393e-05 4.78 3.782e-07 5.01 1.824e-04 3.96 1.695e-114 2.943e-06 4.83 1.177e-08 5.01 1.143e-05 4.00 2.846e-115 1.068e-07 4.78 3.670e-10 5.00 7.144e-07 4.00 6.732e-11Table 2.1: Numerical results for Kovasznay flow (? = 1)well as the error in DG-norm for the velocity are listed in Table 2.5.1 for polynomials of19100 101 102 10310?1010?810?610?410?2100log Nlog(error) k=1k=2k=3k=4(a) convergent plot for ?eu?0100 101 102 10310?710?610?510?410?310?210?1100101log Nlog (error) k=1k=2k=3k=4(b) convergent plot for?ep?0100 101 102 10310?710?610?510?410?310?210?1100101log Nlog(error) k=1k=2k=3k=4(c) convergent plot for?eu?1,hFigure 2.2: The convergent plots for the numerical errors20degree k = 1, ..., 4. The convergence plots are presented in Figure 2.2. We observe optimalorders of convergence for the error of the velocity in both norms, while for the pressure theorders of convergence in the L2-norm are much better than expected in Theorem 2.4.7.This is also observed in [14]. The reason for this is that the abstract error estimate for thepressure (2.60) is bounded by the error of the velocity in ? ? ?1,h norm, which is of orderO(hk). However, in L2-norm, the optimal approximation order of the space Qk shouldbe O(hk+1). Therefore, an order between O(hk) and O(hk+1) is observed. Moreover, thedivergence of the approximate velocity in the L?-norm is evaluated at a set of (k + 1)2quadrature points in each element K. We see from the last column in the Table 2.5.1 thatthe approximate velocity is exactly divergence-free.21Chapter 3Non-isothermal Problem3.1 IntroductionWhen a fluid is subjected to a temperature change, its material properties, such as densityand viscosity, changes accordingly. If these changes are large enough to have a substantialinfluence on the flow field, the flow is considered to be non-isothermal. The numericalsimulation of the non-isothermal incompressible flow has received much attention recentlyin both natural sciences and engineering branches. Relevant industrial applications in-clude polymer processing, cooling processes, the design of heat exchangers and chemicalreactors; see e.g. [30], [26] and [35]. Moreover, in geophysics, the convection processes inthe Earth?s mantle are well described by incompressible fluid flow driven by temperature-induced small density differences; see [29]. The numerical simulation of these processes isa key piece in understanding the dynamics, composition, history and interaction betweenthe mantle and the lithosphere.In the last decade many finite element methods have been developed and analysed forsuch problems. In [15], a mixed finite element method is developed for the non-isothermalStokes-Oldroyd equations. For the same problem, a dual mixed finite element is con-structed and analysed in [19], which has properties analogous to finite volume methods,namely, local conservation of momentum and mass. In [18] a conforming finite elementmethod is presented and analysed for approximating time-dependent non-isothermal in-compressible flow problems. The time is disretized by the backward Euler method. Butnone of these methods have the exactly divergence-free property as our approach does.In this chapter, we first illustrate how the non-isothermal effects are introduced intothe equations by reviewing two classical physical models. Then we review the work of[34] which extends the exactly divergence-free element developed in [13] and [14] to ageneralized Boussinesq system. The model problem couples the stationary incompressible22Navier-Stokes equation and a convection-diffusion equation. The coupling is through atemperature-dependent viscosity, which is also nonlinear, and through a buoyancy termacting in the opposite direction of gravity.We extend the IPDG method described in Chapter 1 to this coupled system and bydoing so, the exactly divergence-free constraint is preserved on the discrete level. Theformulation of the numerical method and the a priori error estimate are first reviewed.Then numerical examples are implemented to verify the theoretical results.3.2 Physical ModelsIn this section, we briefly discuss the physical approximations that introduce the non-isothermal effects into incompressible flow models. The first one is the constitute relationof the Oldroyd-B type for modelling the viscoelastic problems. The second one is theBoussinesq approximation for linearizing buoyancy effect when considering the buoyancydriven flows. These two modelling techniques are introduced separately in the context oftwo different models.3.2.1 Model I: Non-Isothermal Stokes-Oldroyd SystemThe viscoelastic problem occurs in a variety of applications, including polymer process-ing. The complexity of the governing equations makes both the mathematical analysisand the associated numerical methods especially difficult. The current effort is to modelviscoelastic flows as the solution of a modified Stokes problem; see [6, 10]. In [15], thenon-isothermal Stokes problem is modified in the same way as in [6].We consider a fluid flowing in a bounded and connected domain ? ? R2, whose boundaryis denoted by ? and we assume the stationary and creeping flow hypothesis. Let thevelocity be denoted by u, pressure by p, extra stress tensor by ?, rate of strain tensor byD(u) = 12 (?u + ?uT ) and temperature by T . In the Oldroyd-B model the extra stresstensor is split into the Newtonian solvent part ?s = ?sD(u) and the non-Newtonianpolymer part ?p, i.e.? = ?s + ?p.?p is given by?p +We?a?p?t= 2?D(u); a ? [?1, 1], ? ? (0, 1) (3.1)where We is the Weissenberg number and the derivative?a?p?t is defined by???t= ? t + (u ? ?)? + ga(? ,?u), a ? [?1, 1] (3.2)23andga(? ,?u) = ??(u)? ?(u)? ? a(D(u)? + ?D(u)).This is a variation of the upper-convected time derivative (see [27]), which is the materialderivative of the tensor property in the coordinate system rotating and stretching withthe fluid.In the limit case We = 0 we can get the Oldroyd-Stokes problem?? ? [(2?p + 2?s)D(u)] +?p = f1 in ?? ? u = 0 in ?u = g1 on ?(3.3)where ?p and ?s are the polymetric viscosity and the solvent one respectively and > 0is a small parameter that makes the effect of ?s much smaller than that of ?p.Many flows of interest in polymetric melt processing are non-isothermal. We followthe derivation in [15] for the non-isothermal Oldroyd-Stokes system. By introducingtemperature-dependent coefficients ?p(T ) and ?s(T ) and by combining the Oldroyd-Stokessystem (3.3) and a convection-diffusion equation for the temperature, we can obtain?? ? [(2?p(T ) + 2?s(T )D(u)] +?p = f1 in ?? ? u = 0 in ????T + u ? ?T = f2 in ?u = g1 on ?T = g2 on ?(3.4)where ? is the thermal conductivity coefficient, f1 is the body force, f2 the heat source.The dependence of the polymer and solvent viscosity on the temperature is given by theArrhenius equations (see [25]):?p(T ) = a1 exp(b1T), ?s(T ) = a2 exp(b2T), b1 = 0. (3.5)The constants ai, bi, i = 1, 2 are chosen such that0 < ?s ? 1, 0 < ?s < 1243.2.2 Model II: Buoyancy Driven FlowsIf the fluid density varies with temperature, a flow can be induced in the presence of gravity.This is known as buoyancy-driven flow. Many flow phenomena are driven by buoyancy,and such flows are especially important in the modelling of the convection processes in theEarth?s mantle; see [38] . When the viscous friction forces in the fluids are large comparedto buoyancy forces, for example when simulating the mantle, the motion of the fluid isslow and inertial terms can be neglected; see [29] and the book [38]. Combining this withthe steady-state assumption, we obtain governing equations in the following form?? ? (2?D(u)) +?p = g ?(T ) in ?? ? u = 0 in ??? ? ??T + u ? ?T = f2 in ?u = g1 on ?T = g2 on ?(3.6)where ? is the viscosity, g the gravity acceleration, ?(T ) the temperature-dependent den-sity. Other notations are the same as those in last section.The standard assumption made to simplify the analysis of the system above is the Boussi-nesq approximation, in which the density differences are neglected, except in the buoyancyterm, where they appear in terms multiplied by g. If all accelerations involved in the floware small compared to g, the dependence of the density on T in the buoyancy term canbe assumed to be linear:?(T )? ?0 = ??0?(T ? T0).Here, T0 is a reference temperature, e.g., the temperature on a boundary, and ? is thecoefficient of expansion for the fluid. The corresponding buoyancy force is then given asfB = ??0g?(T ? T0).Then the nonlinear equation of the conservation of momentum can be linearized into?? ? (2?D(u)) +?p = g?0(1? ?(T ? T0))? ?0g?TIf we express the gravity in the potential form g = ??? and introduce the effectivepressure peff = p+ ??0(1? ?(T ? T0)), the equation above can be further rewritten as?? ? (2?D(u)) +?peff = ??0g?T.253.3 Model ProblemIn this section, to ease the mathematical analysis, we consider the following generalizedmodel problem, which captures all the non-isothermal effects, with zero heat source andhomogeneous velocity field on the boundary:?? ? (?(T )?u) + u ? ?u+?p = jT in ? (3.7a)? ? u = 0 in ? (3.7b)?? ? (?(T )?T ) + u ? ?T = 0 in ? (3.7c)u = 0 on ? (3.7d)T = gT on ? (3.7e)where ? is the computational domain in R2 and T is the temperature, u the velocity, pthe pressure, j ? L2(?) is a constant vector in the opposite direction of the gravity accel-eration. The boundary datum gT is assumed to be in C(??). For the rest of this chapter,we will use the same notations for the mesh, function spaces and numerical traces as thosedefined in Chapter 1.The temperature-dependent parameters ?(T ) and ?(T ) denote the effective viscosity andheat conductivity, respectively. We assume they are Lipschitz continuous and bounded,that is|?(T1)? ?(T2)| ? L?|T1 ? T2|, |?(T1)? ?(T2)| ? L?|T1 ? T2| (3.8)for all values of T1 and T2, with the Lipschitz constants L?, L? > 0. And there existconstants ?1, ?2, ?1, ?2 > 0, such that0 < ?1 ? ?(T ) ? ?2, 0 < ?1 ? ?(T ) ? ?2 (3.9)for all T .The variational formulation of the problem (3.7) is to find (u, p, T ) ? H1(?)2 ? L20(?) ?H1(?) with T |? = gT , such thatA(u,v;T ) +B(v, p) +ON (u,v;u) ?D(T,v) = 0 (3.10a)?B(u, q) = 0 (3.10b)C(T, S;T ) +OH(T, S;u) = 0 (3.10c)26for all (v, q, S) ? H10 (?)2 ? L20(?)?H10 (?), and where the bilinear and linear forms areA(u,v;O) =???(O)?u : vdx B(v, q) =??(?h ? u)q dxON (u,v;w) =??(w ? ?)u ? v dx OH(T, S;u) =??(u ? ?T )S dxD(T,v) =??T j ? v dx C(T, S;O) =???(O)?T ? ?S dx3.3.1 Stability ResultsFirst, let us review some stability results of the forms appearing in (3.10) from [34].Boundedness and ContinuityDue to the assumption (3.9), we have the following continuity properties:|A(u,v; ?) ? ?2?u?1,??v?1,?, u,v ? H1(?)2, (3.11)|C(T, S; ?)| ? ?2?T?1,??S?1,?, T, S ? H1(?), (3.12)|B(v, q)| ? CB?v?1,??q?0,?, v ? H1(?)2, q ? L2(?). (3.13)Moreover, from the Lipschitz continuity of ? and ? in (3.8) and Ho?lder?s inequality, thefollowing Lipschitz continuity properties hold for the form A and C: for T1, T2 ? H1(?),u ? W 1,?(?), T ? W 1,?(?), we have|A(u,v;T1)?A(u,v;T2)| ? LA?u?W1,?(?)?T1 ? T2?1,??v?1,?, v ? H1(?) (3.14)|C(T, S;T1)? C(T, S;T2)| ? LT ?T?W 1,?(?)?T1 ? T2?1,??S?1,?, S ? H1(?) (3.15)Due to Ho?lder?s inequality and the standard Sobolev embeddings, the trilinear forms ONand OH have the following bounds:|ON (u,v;w)| ? CN?w?1,??u?1,??v?1,?, w,u,v ? H1(?) (3.16)|OH(T, S;w)| ? CH?w?1,??T?1,??S?1,?, w ? H1(?), T, S ? H1(?) (3.17)Also, we have the bound on D|D(T,v)| ? CD?j?0,??T?1,??v?1,?, ? ? H1(?), v ? H1(?) (3.18)Coercivity and Inf-sup conditionNext, we review the coercivity of the linear forms. By the Poincare? inequality and thebounds in (3.9), we haveA(u,u; ?) ? ?A?u?21,?, v ? H10(?), (3.19)C(T, T ; ?) ? ?C?T?1,?, T ? H10 (?) (3.20)27To discuss the coercivity of the forms OH and ON , we need to introduce the kernelX = {v ? H1(?) : B(v, q) = 0, ?q ? L20(?)} = {v ? H1(?) : ? ? v ? 0} (3.21)Then integration by parts gives us,ON (u,u;w) = 0, w ? X, v ? H1(?) (3.22)OH(T, T ;w) = 0, w ? X, T ? H1(?) (3.23)Finally the bilinear form B satisfies the continuous inf-sup conditionsupv?H10(?)\{0}B(v, q)?v?1,?? ??q?0,?, ?q ? L20(?) (3.24)with an inf-sup constant ? > 0 only depending on ?; see [20].3.3.2 Unique SolvabilityIn this section, we review the results related to the existence and uniqueness of the solutionto the system (3.7). To that end, it?s enough to study the reduced problem on the kernelX which is defined in (3.21). The reduced problem consists of finding (u, p) ? V?H1(?)such that T |? = gT andA(u,v;T ) +ON (u,v;u) ?D(T,v) = 0, (3.25a)C(T, S;T ) +OH(T, S;u) = 0 (3.25b)for all (v, T ) ? X?H10 (?).First, we have the following equivalence property; see [20], which justifies our consid-eration of the reduced problem.Lemma 3.3.1 If (u, p, T ) ? H1(?) ? L20? ?H1(?) is a solution of (3.10), then u ? Xand (u, T ) is also a solution of (3.25). Conversely, if (u, T ) ? X?H1(?) is a solution of(3.25), then there exists a unique pressure p ? L20(?) such that (u, p, T ) is a solution of(3.28).The next theorem ensures the solvability of the reduced problem (3.25).Theorem 3.3.2 (Solvability) Assume (3.9) and (3.8). Then for any j ? L2(?), thereexists a lifting T1 ? H1(?) of g1 ? H1/2(?) satisfying T1|? = g1, such that the prob-lem (3.10) has a solution (u, T ) ? H10(?) ? H1(?), where T = T1 + T0 and T0 ? H10 .Furthermore, there exist constants Cu and CH depending only on ?j?0,? and the stabilityconstants ?A, ?C , ?2, ?2, CD, such that?u?1,? ? Cu?T1?1,? ?T?1,? ? CH?T1?1,? (3.26)28Under additional smoothness and smallness assumptions on the solutions, the followinguniqueness theorem holds.Theorem 3.3.3 Let (u, T ) ? [X ? W1,?] be a solution to problem (3.25), and assumethere exists a sufficiently small constant 0 < M < min{ ?ACSC?+K ,?CL?+K }, where K =(CTC? + L?)/2, such thatmax{?j?0,?, ?u?W1,??T?W 1,?} < M. (3.27)Then, the solution is unique.Proof: The proof can be found in [34, Theorem 2.3], which is motivated by a similarargument in [6] for Stokes-Oldroyd problems.3.4 Discrete ProblemTo formulate the discontinuous Galerkin method, we first introduce the following discretespacesVh := {v ? H(div; ?) : v|K ? RTk(K), ?K ? Th, v ? n = 0, on ?}Qh := {q ? L20(?) : q|K ? Qk(K), ?K ? Th}Eh := {S ? C0(??) : S|K ? Qk(K), ?K ? Th; S = gTh , on ?}E0h := {S ? C0(??) : S|K ? Qk(K), ?K ? Th; S = 0, on ?}?h := {? ? C(??) : ?|e ? Qk(e), e ? EBh }where gTh ? ?h is the approximate boundary datum which will be defined in Section ??.We have used the same element pair for the approximation of the velocity and the pres-sure. And the temperature is approximated with the standard conforming method.By following the discretization of the viscous term, pressure gradient and the incom-pressibility constraint in Chapter 1 and by using the standard upwind form in [31] for thenonlinear convection term defined byONh (u,v;w) =?K?Th?K(w ? ?)u ? v dx+?K?Th??K\??12(w ? nK ? |w ? nK |)(ue ? u) ? v dsand conforming discretization for the convection-diffusion equation (3.10c), we define theDG method as: find (uh, ph, Th) ? Vh ?Qh ?Eh, such thatAh(uh,v;Th) +ONh (uh,v;uh) +Bh(v, ph)?Dh(Th,v) = 0 (3.28a)?Bh(uh, q) = 0 (3.28b)C(Th, S;Th) +OH(Th, S;uh) = 0 (3.28c)29for all (v, q, S) ? V0h ?Qh ? E0h, whereAh(u,v;O) =???(O)?hu : ?hv dx???(?h ? v)p dx+?0h?EIh?(T )u : v ds??EIh?(O)u : {{?hv}} ds ??EIh?(O)v : {{?hu}} ds+?0h?EBh?(O)u ? v ds??EBh?(O)(u ? n) : ?hv ds??EB?(O)(v ? n) : ?hu ds (3.29)ONh (u,v;w) =?K?Th?K(w ? ?)u ? v dx+?K?Th??K\??12(w ? nK ? |w ? nK |)(ue ? u) ? v ds????12(w ? nK ? |w ? nK |)(u ? v) ds (3.30)Bh(v, q) =??(?h ? u)q dx (3.31)Dh(v, T ) =??T j ? v dx (3.32)C(T, S;O) =???(O)?T ? ?S dx (3.33)OH(T, S;w) =??(w ? ?T )S dx (3.34)For the boundary datum gT , we assumegT ? C(??). (3.35)For continuous function in C(??), we can define the Larange interpolation operator I :C(??) ? Eh. Its restriction to the boundary nodes is denoted by I? : C(??) ? ?h. Wetake the discrete boundary datum to be gTh = I?gT .Since we have employed the same approximation spaces for the velocity and pressureand the same incompressibility form Bh(?, ?) as in Chapter 1, we should expect the exactsatisfaction of the incompressibility constraint ? ? uh ? 0 to be preserved. Actually, wehave the following proposition:30Proposition 3.4.1 The approximate velocity uh ? Vh obtained by (3.28) is exactlydivergence-free, i.e., it satisfies? ? uh ? 0 (3.36)Proof: The proof is exactly the same as that of Proposition 3.3 in Chapter 1.As discussed in [13], an important consequence of Proposition 3.4.1 is the provable energy-stability of the numerical method (3.28), without modifying the convective terms in theincompressible Navier-Stokes equation. Moreover, the exact conservation of mass is ofgreat importance in the simulation of the incompressible flow and nearly incompressiblematerials; see [32] and [4].3.5 StabilityIn this section, we review the stability results in [34, Section 3.3] for the linear formsappearing in the discrete problem (3.28) and adapt them to the RTk ?Qk space pair.3.5.1 PreliminariesWe restrict the discussion in the broken Sobolev spacesHr(Th) = {v ? L2(?) : v|K ? Hr(K), K ? Th}, r = 1 or 2 (3.37)endowed with the norms?v?21,h =?K?Th??hv?20,K +?e?Eh?0h?1e ?v?20,e, v ? H1(Th) (3.38)?v?22,h = ?v?21,h +?K?Thh2k|v|22,K , v ? H2(Th) (3.39)For K ? Th, q ? Qk(K), by the standard inverse estimate (see e.g. [7, Lemma 4.5.3]), wehave |q|2,K ? Ch?1|q|1,K . Then it?s easy to see that?v?2,h ? C?v?1,h, v ? Vh (3.40)In the analysis of the stability, we need the following embedding result for d = 2 (see [20,Lemma 6.2]): For v ? H1(Th) and ?p ? [1,?), there exists a constant C > 0 such that?v?Lp(?) ? C?v?1,h (3.41)Moreover, we need the broken C1(Th)-space given byC1(Th) = {u ? H1(Th) : u|K ? C1(K), K ? Th} (3.42)31equipped with the broken W1,?-norm?u?W1,?(Th) = maxK?Th?u?W1,?(K) (3.43)At last, we assume the convective velocity w ? H0(div0,?), where the space H0(div0,?)is defined to beH0(div0,?) := {w ? H(div; ?) : ? ?w = 0, w ? n = 0} (3.44)3.5.2 ContinuitySince we employ conforming elements for the convection-diffusion equation (3.62c), thecontinuity of C(Th, S; ?) and the Lipschitz continuity of C(?, ?;T ) are direct consequencesof (3.12) and (3.15). For the elliptic form Ah, since we have the Qk(K)-version traceinequality [37, Lemma 7.1]?q?0,?m ? Ch?1/2K ?q?0,K , ?q ? Qk(K). (3.45)as well as the inverse estimate (3.40), we proceed in the standard argument in [3] to obtainthe following result.Lemma 3.5.1 (Boundedness of Ah) For the form Ah defined in (3.29), there holds|Ah(u,v; ?)| ? C?u?2,h?v?1,h, u ? H2(Th), v ? Vh, (3.46)|Ah(u,v; ?)| ? C?A?u?1,h?v?1,h, u, v ? Vh (3.47)The Lipschitz contiuity of Ah(?, ?;T ) is presented in the next lemmaLemma 3.5.2 (Lipschitz continuity of Ah) Let T1, T2 ? H1(?), u ? H10(?) and v ?Vh, then there holds|Ah(u,v;T1)?Ah(u,v;T2)| ? LAh L??T1 ? T2?1,??u?1,?(Th)?v?1,h (3.48)Proof: The proof can be found in [34, Lemma 3.3].Moreover, we have the boundedness of Bh and Dh:Lemma 3.5.3 We have the following bounds for the forms Bh and Dh.|Bh(v, q)| ? C?B?v?1,h?q?0,?, v ? H1(Th), q ? L20(?), (3.49)|Dh(S,v)| ? C?D?j?0,??S?1,??v?1,h, v ? H1(Th), S ? H1(?) (3.50)32Proof: Equation (3.49) is a straightforward application of the Cauchy-Schwarz inequal-ity. For the estimate of Dh, by Ho?lder?s inequality and the embedding inequality (3.41),we have|Dh(S,v)| ? ?j?0,??Sv?0? ? ?j?0,??S?L4,??v?L4,? ? C?D?j?0,??S?1,??v?1,hAt last, we review the continuity results for the convection forms ONh and OH .Lemma 3.5.4 Let w1, w2 ? H0(div0; ?), u ? H2(?) and v ? Vh. Then we have|ONh (u,v;w1)?ONh (u,v;w2)| ? LONh ?w1 ?w2?1,h?u?1,h?v?1,h (3.51)Proof: Please refer to [13, Proposition 4.2].For the conforming convective form OH(T, S;w), in the case that w ? H0(div0; ?), wehave the following variant of the result in (3.17) .Lemma 3.5.5 Let w ? H0(div0; ?) and T, S ? H1(?), we have|OH(T, S;w)| ? C?H?T?1,??w?1,h?S?1,? (3.52)Proof: The proof can be found in [34, Lemma 3.5], which involves the application ofintegration by parts, Ho?lder?s inequality and the embedding result in (3.41).3.5.3 Coercivity and inf-sup conditionIn this section we review the coercivity results of the forms Ah, C, ONh and OH as well asthe inf-sup condition for Bh.Lemma 3.5.6 For S ? E0h and v ? Vh we have the following resultsC(S, S; ?) ? ?C?S?1,? (3.53)Ah(v,v; ?) ? ??A?v?1,h (3.54)Proof: Since E0h ? H10 (?), the first inequality is a direct consequence of (3.20). Bynoticing that the perturbed form A?h = Ah when restricted to the discrete space Vh, thesecond one follows the covercivity result, Proposition 2.4.1 in Chapter 1.For the convective terms ONh and OH we have the following positivity results.33Lemma 3.5.7 Let w ? H0(div0; ?). Then we haveONh (u,u;w) =12?e?EIh?e|w ? n|| u? n |2 ds ? 0, u ? Vh, (3.55)OH(T, T ;w) = 0, T ? H1(?) (3.56)Proof: Since w ? H0(div0; ?) and if we replace v with u in (3.30), we can integrate thecell terms by parts and obtain?K?Th?K(w ? ?)u ? u dx =?K?Th?K12w ? ?|u|2 dx = 12?K?Th??K|u|2w ? n ds=12(?EIh{{|u|2}} w ds+?Eh{{w}} ?u|2ds)=12?EIh{{w}} ?u|2ds=12?EIhw ? n1(|u1|2 ? |u2|2) dsHere we have used ? ? w = 0 in ?, the summation identity, w ? n = 0 on ? and thecontinuity of w ? n across the inner edges. And we choose the unit normal vector n1 oneach inner edge such that w ? n1 ? 0.Moreover, the other edge terms in the ONh (u,u;w) can be rewritten into?K?Th??K\??12(w ? nK ? |w ? nK |)(ue ? u) ? u dx =?e?EIh?e(w ? n)(ue ? u) ? u dx.It is obvious that w ?n ? 0. Then in order to be consistent with our choice of the normaldirection above, we add the subindex as below?K?Th??K\??12(w ? nK ? |w ? nK |)(ue ? u) ? u dx=?e?EIh?e(w ? n2)(u1 ? u2) ? u2 dx=??EIh(w ? n1)(u1 ? u2 ? |u2|2) dx34Then we can sum the cell terms and the edge terms together and arrives atONh (u,u;w) =12?EIhw ? n1(|u1|2 ? |u2|2) ds ??EIh(w ? n1)(u1 ? u2 ? |u2|2) dx=12?EIh(w ? n1)|u1 ? u2|2 ds=12?EIh|w ? n|| u? n |2 ds ? 0The positivity of ON can be very easily obtained by integration by parts and the proper-ties of the space H0(div0; ?).We end this section with the inf-sup condition of the bilinear form B(?, ?):infq?Qhsupv?Vh\{0}B(v, q)?v?1,h?q?0,?? ? (3.57)where the inf-sup constant ? > 0 is independent of h. This condition is same as the onein Chapter 1. The proof can be found in [37].3.6 Error AnalysisIn this section, we review the a priori error estimate for the DG method (3.28). Undersuitable smoothness and smallness assumptions, we derive the estimate in a straightfor-ward way without introducing the lifting operators; see [3].3.6.1 AssumptionsSmoothness assumption:Let (u, p, T ) be the exact solution of (3.7), we assumefor k = 1 : u ? C1(??) ?H2(?)2 ?X, p ? H1(?), T ? W 1,?(?) ?H2(?)for k = 2 : u ? Hk+1(?)2 ?X, p ? Hk(?), T ? Hk+1(?)(3.58)Smallness assumption:max{?j?0,?, ?u?W 1,?(?), ?T?W 1,?(?)} ? M (3.59)for a sufficiently small number M .353.6.2 Main resultWe are ready to present the error estimate resultTheorem 3.6.1 Let (u, p, T ) be a solution of (3.7) and (uh, ph, Th) be the approximatesolution to (3.28). Under the smoothness assumption (3.58), smallness assumption (3.59)and the assumption on the boundary datum, we have the following estimate:?u? uh?2,h + ?? ? ?h?1,? ? Chk(?u?k+1,? + ?T?k+1,?), (3.60)?p? ph?0,? ? Chk(?p?k,? + ?u?k+1,? + ?T?k+1,?) (3.61)where C > 0 is a generic constant which is independent of h.Proof: In section 3.5, we have obtained the same set of stability results for our choice ofthe velocity-pressure element pair as those in [34, Section 5]. We can follow the same linein the proof of [34, Theorem 5.1] and with the standard approximation results to arriveat the conclusion. We omit the details here.3.7 Fixed-point IterationThe nonlinear nature of the problem requires us to invoke iterative methods and solvea linearized problem in each iteration. Here we follow [14] and [21] and employ thefixed-point method. We start from a initial guess T0 and u0h which are the solutions ofa Stokes problem and a Laplace problem, respectively. Provided we have obtained theapproximate velocity, pressure and temperature xnh = (unh, pnh, Tnh ) from the nth iteration,then the following linearized problem will be solved for xn+1h = (uh, ph, Th)Ah(uh,v;Tnh ) +ONh (uh,v;unh) +Bh(v, ph)?Dh(T nh ,v) = 0 (3.62a)?Bh(uh, q) = 0 (3.62b)C(Th, S;Tnh ) +OH(Th, S;unh) = 0 (3.62c)for all test functions (v, q, S), which is a decoupled system consisting of a Oseen equationand a convection-diffusion equation. The iteration is terminated once the difference of theentire coefficient vector coefficients between two consecutive iterates is sufficiently small,i.e.?coeffn+1 ? coeffn?l2 ? tol (3.63)where ? ? ?l2 is the standard l2 norm over RNdofs and ?tol? is the fixed tolerance which istaken to be 10?8 in the following experiments.363.8 Numerical ExperimentsIn this section we present the numerical results for the coupled system (3.7). The codeis developed with the deal.II finite element library [5]. Since the goal is to confirm theoptimal convergence rates for our method and the exactly divergence-free approximatevelocity, we employ the direct linear solver provided by the UMFPACK [16].3.8.1 Smooth SolutionThe computational domain is taken to be? = (?1, 1)2and is decomposed uniformly into rectangular mesh. At the lth level, the meshsize ish = 2?l. The test solution is taken to beu1 = sin(y)u2 = sin(x)p = 1 + sin(xy)T = 1 + cos(xy)And we take the temperature-dependent parameters to be?(T ) = exp(?T ), ?(T ) = exp(T ) (3.64)and the vector j = (0, 1)T . Appropriate right hand side functions are assigned to balancethe equation and discretized in the discrete problem. In the numerical methods, thepenalty parameter is chosen as?0 = (k + 1)2where k is the approximation order.Initial GuessFor the fixed-point iteration described in Section 3.7, we take the initial guess to bex0h = 0. In this way, in the first iteration, we actually solve a decoupled Stokes problemand a Laplace problem.Since the boundary values of u = gu is inhomogeneous, the essential boundary condi-tion uh ? n = guh ? n in Vh must be enforced strongly in each iteration. To this end,we follow the discussion in Chapter 1 for the enforcement of inhomogeneous boundaryconditions for the Stokes problem.37Numerical ResultsThe numerical results of are presented in Table 3.8.1 and Table 3.8.2 as well as in theconvergent plots Figure 3.1, 3.2. The expected optimal convergence rates for the error ofthe velocity as well as the temperature in the L2-norm are observed. While for the pressure,just like the results for the Stoke?s problem, the convergence order is still between O(hk)and O(hk+1). In Table 3.8.1, the errors for the velocity in the DG norm and for thetemperature in the H1 norm are presented. The convergence order O(hk) is observed.Moreover, the exactly divergence-free property is verified by evaluating ?? ? uh?? over aset of quadrature points. At last, the residual plot for the case of k = 1 is presented inFigure 1. We see that the iteration numbers increase moderately as we refine the meshand the residuals decrease in a linear fashion. For k > 1, in our experiment, the samephenomenon is observed.k l ?ep?0 ?eT ?0 ?eu?011 9.86e-02 - 7.16e-02 - 6.75e-02 -2 3.50e-02 1.49 1.92e-02 1.90 2.12e-02 1.673 1.32e-02 1.41 4.96e-03 1.95 5.02e-03 2.084 4.62e-03 1.52 1.25e-03 1.98 1.21e-03 2.065 1.53e-03 1.60 3.14e-04 2.00 3.09e-04 1.9721 1.044e-02 - 6.328e-03 - 1.101e-02 -2 1.611e-03 2.70 6.998e-04 3.18 1.306e-03 3.083 2.389e-04 2.75 8.498e-05 3.04 1.600e-04 3.034 3.735e-05 2.68 1.055e-05 3.01 1.975e-05 3.0231 6.910e-04 - 1.520e-04 - 4.713e-04 -2 7.407e-05 3.22 1.009e-05 3.91 3.295e-05 3.843 7.033e-06 3.40 6.455e-07 3.97 2.049e-06 4.014 6.318e-07 3.48 4.083e-08 3.98 1.242e-07 4.0441 6.456e-05 - 1.201e-05 - 4.212e-05 -2 2.375e-06 4.76 3.402e-07 5.14 1.121e-06 5.233 9.623e-08 4.63 1.057e-08 5.01 3.270e-08 5.10Table 3.1: Numerical results I: IPDG for the non-isothermal problem (3.7) with fixed-pointiteration38100 101 10210?810?710?610?510?410?310?210?1log Nlog error k=1k=2k=3k=4(a) convergent plot for ?eu?0100 101 10210?810?710?610?510?410?310?210?1log Nlog error k=1k=2k=3k=4(b) convergent plot for?ep?0100 101 10210?810?710?610?510?410?310?210?1log Nlog error k=1k=2k=3k=4(c) convergent plot for?eT?0Figure 3.1: The convergent plots for the numerical errors100 101 10210?610?510?410?310?210?1100log Nlog (error) k=1k=2k=3k=4(a) convergent plot for ?eu?1,h100 101 10210?710?610?510?410?310?210?1100log Nlog error k=1k=2k=3k=4(b) convergent plot for?eT?1,?Figure 3.2: The convergent plots for the numerical errors39k l ?eT?1,? ?eu?1,h ?? ? uh??11 2.201e-01 - 5.982e-01 - 1.443e-152 9.332e-02 1.24 2.831e-01 1.08 5.329e-153 4.267e-02 1.13 1.315e-01 1.11 1.199e-144 2.079e-02 1.04 6.195e-02 1.09 2.931e-145 1.030e-02 1.01 2.927e-02 1.08 6.584e-1421 4.108e-02 - 1.591e-01 - 4.732e-152 9.169e-03 2.16 4.546e-02 1.81 1.521e-143 2.223e-03 2.04 1.199e-02 1.92 2.429e-144 5.513e-04 2.01 3.060e-03 1.97 5.556e-1431 1.451e-03 - 7.545e-03 - 1.484e-142 1.805e-04 3.01 9.239e-04 3.03 3.040e-143 2.240e-05 3.01 1.077e-04 3.10 6.105e-144 2.791e-06 3.00 1.206e-05 3.16 1.131e-1341 1.703e-04 - 9.295e-04 - 4.174e-142 9.917e-06 4.10 6.241e-05 3.90 7.412e-143 5.981e-07 4.05 4.010e-06 3.96 1.645e-13Table 3.2: Numerical results II: IPDG for the non-isothermal problem (3.7) with fixed-point iterationFigure 3.3: The residual plot for k = 1, l = 1, ..., 41 2 3 4 5 6 7 8 9 1010?1010?810?610?410?2100102iterationresidual grid 1grid 2grid 3grid 4403.9 Conclusions and future directionsIn this chapter we first review the formulation and the error analysis in [34] for the exactlydivergence-free finite element methods [14], with an H(div) conforming element for theapproximate velocity, for a generalised Boussinesq equation [33]. Numerical experimentsbased on a fixed-point iteration are implemented which verify the theoretical results, espe-cially the error estimates and the exactly divergence-free property. Future work involves:? More numerical experiments. In this project, only smooth functions in the 2D unit-square domain with Dirichlet boundary conditions is considered. We need to furthertest the code on singular solutions and more realistic examples in general domainswith mixed boundary conditions such as the staircase flows, contraction flows andflows in 3D.? Efficient nonlinear iteration methods and linear solvers. The first-order fixed pointiteration can be replaced with the more efficient Newton?s method, while in doingso, the convection term need to be discretized wisely. Even though the direct solversare super stable, as the mesh is refined or the dimension increased to 3D, the storagerequirement for the increasingly large system matrix will be a bottle-neck for theefficiency of the method. Iterative methods need to be considered and efficientpreconditioners need to be designed in order to accelerate the convergence.? More effort in theoretical analysis. As emphasized in [34], the uniqueness of theapproximate solution is still a open problem. Even for the existence and the errorestimate, very restrictive regularity assumptions (3.58), (3.35) and smallness as-sumption (3.59) have to be made. It would be another challenge to get rid of theseassumptions, while preserving the same optimal error estimates.41Bibliography[1] D. N. Arnold. An interior penalty finite element method with discontinuous elements.SIAM J. Numer. Anal., 19:742?760, 1982.[2] D. N. Arnold, F. Brezzi, and M. Fortin. A stable finite element for the Stokesequations. Calcolo, 21:337?344, 1984.[3] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, and L. Donatella Marini.Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J.Numer. Anal., 39(5):1749?1779, 2002.[4] F. Auricchio, L. Beirao Da Veiga, C. Lovadina, and A. Reali. The importance of theexact satisfaction of the incompressibility constraint in nonlinear elasticity: mixedFEMs verses NURBS-based approximations. Comp. Math. Appl. Mech. and Eng.,199:314?323, 2010.[5] W. Bangerth and G. Kanschat. Concepts for object-oriented finite element software-the deal.II library. Preprint 43, SFB 359, 1999.[6] J. Baranger and D. Sandri. A formulation of the Stoke?s problem and the linear elas-ticity equations suggested by the Oldroyd model for viscoelastic flow. MathematicalModelling and Numerical Analysis, 1:331?345, 1992.[7] S. C. Brenner and L. R. Scott. The Mathematical Theory of Finite Element Methods.Springer, 2002.[8] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer Verlag,1991.[9] A. Buffa, C.de Falco, and G. Sangalli. IsoGeometric analysis: Stable elements forthe 2D Stokes equation. International Journal for Numerical Methods in Fluids, 65:1407?1422, 2011.[10] A. E. Caola, Y. Joo, and R. C. Armstrong. Highly parallel time integration ofviscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 100:191?216, 2001.42[11] P. G. Ciarlet. The Finite Element Methods for Elliptic Problems. NORTH-HOLLAND, 1978.[12] P. G. Ciarlet. Mathematics Elasticity. NORTH-HOLLAND, 1988.[13] B. Cockburn, G. Kanschat, and D. Scho?tzau. A locally conservative LDG methodfor the incompressible Navier-Stokes equations. Math. Comp., 74:1067?1095, 2005.[14] B. Cockburn, G. Kanschat, and D. Scho?tzau. A note on discontinuous Galerkindivergence-free solutions of the Navier-Stokes equations. Journal of Scientific Com-puting, 31:61?73, 2007.[15] C. Cox, H. Lee, and D. Szurley. Finite element approximation of the non-isothermalStokes-Oldroyd equations. International Journal of Numerical Analysis and Mod-elling, 1:1?18, 2007.[16] T. A. Davis. Algorithm 832: UMFPACK V4.3?an unsymmetric-pattern multifrontalmethod. ACM Transactions on Mathematical Software, 30:196?199, 2004.[17] C. e. P., J. Thomas, S. Blancher, and R. Creff. The steady Navier-Stokes/energy sys-tem with temperature-dependent viscosity-Part 1: Analysis of the continuous prob-lem. International Journal for Numerical Methods in Fluids, 56:63?89, 2008.[18] C. e. Pe?rez, J. Thomas, S. Blancher, and R. Creff. The steady Navier-Stokes/energysystem with temperature-dependent viscosity-Part 2: The discrete problem and nu-merical experiments. International Journal for Numerical Methods in Fluids, 56:91?114, 2008.[19] M. Farhloul and A. Zine. A dual mixed formulation for non-isothermal Oldroyd-Stokesproblem. Mathematical Modelling of Natural Phenomena, 6 (05):130?156, 2011.[20] V. Girault and P. A. Raviart. Finite element methods for Navier-Stokes equations.Springer, New York, 1986.[21] C. Greif, D. Li, D. Scho?tzau, and X. Wei. A mixed finite element method with ex-actly divergence-free velocities for incompressible magnetohydrodynamics. ComputerMethods in Applied Mechanics and Engineering, 199:2840?2855, 2010.[22] J. Guzma?n and M. Neilan. Conforming and divergence free Stokes elements on generaltriangular meshes. Mathematics of Computation, 2011.[23] J. Guzma?n and M. Neilan. Conforming and divergence-free Stokes elements in threedimensions . IMA Journal of Numerical Analysis, 2013.43[24] P. Hood and C. Taylor. Numerical solution of the Navier-Stokes equations using thefinite element technique. Compt. Fluids, 1:1?28, 1973.[25] J. E. House. Principles of Chemical Kinetics. Academic Press, Elsevier, 2004.[26] Y. Joo, J. Sun, M. Smith, R. Armstrong, R. Brown, and R. Ross. Two-dimensionalnumerical analysis of non-isothermal melt spinning with and without phase transition.Journal of Non-Newtonian Fluid Mechanics, 102:37?70, 2002.[27] R. E. Khayat and M. O. Starzewski. On the objective rate of heat and stress fluxes.Connection with micro/nano-scale heat convection. Discrete and Continuous Dy-namical Systems-Series B, 15 (04):991?998, 2011.[28] L. Kovasznay. Laminar flow behind a two-dimensional grid. Proc. Camb. Philos. Soc,44:58?62, 1948.[29] M. Kronbichler, T. Heister, and W. Bangerth. High accuracy mantle convectionsimulation through modern numerical methods. Geophys. J. Int., 191:12?29, 2012.[30] K. Kunisch and X. Marduel. Optimal control of non-isothermal viscoelastic fluid flow.Journal of Non-Newtonian Fluid Mechanics, 88:261?301, 2000.[31] P. LeSaint and P. A. Raviart. On a finite element method for solving the neutrontransport equation. 1974.[32] A. Linke. Collision in a cross-shaped domain A steady 2D Navier-Stokes exampledemonstrating the importance of mass conservation in CFD. Comp. Math. Appl.Mech. and Eng., 198:3278?3286, 2009.[33] S. A. Lorca and J. L. Boldrini. Stationary solutions for generalized boussinesq models.Journal of Differential Equations, 124:389?406, 1996.[34] R. Oyarzu?a, T. Qin, and D. Scho?tzau. An exactly divergence-free finite elementmethod for a generalized Boussinesq problem. IMA Journal of Applied Mathematics,accepted.[35] G. W. M. Peters and F. T. O. Baaijens. Modelling of non-isothermal viscoelasticflow. Journal of Non-Newtonian Fluid Mechanics, 68:205?254, 1997.[36] D. A. D. Pietro and A. Ern. Mathematical Aspects of the Discontinuous GalerkinMethods. Springer, 2012.[37] D. Scho?tzau, C. Schwab, and A. Toselli. hp-DGFEM for incompressible flows. SIAMJ. Numer. Anal., 40:561?571, 2003.44[38] G. Schubert, D. L. Turcotte, and P. Olson. Mantle Convection in the Earth andPlanets, Part 1. Cambridge University Press, Cambridge, 2001.[39] L. R. Scott and M. Vogelius. Conforming finite element methods for incompressibleand nearly incompressible continua. Lectures in Applied Mathematics, 22:221?244,1985.[40] R. Te?mam. Sur l?approximation des solutions des e?quations de Navier-Stokes. C. R.Acad. Sci. Paris Se?r., A 216:219?221, 1966.[41] R. Te?mam. Une me?thode d?approximation des solutions de la solutions des e?quationsde Navier-Stokes. Bull. Soc. Math. France, 98:115?152, 1968.[42] S. Zhang. Some divergence-free finite elements for the stationary Stokes equations.In Y. Y. Li, C. W. Shu, R. G. Ye, and K. Zuo, editors, Advances in Mathematics andIts Applications, pages 261?282. the University of Science and Technology of ChinaPress, Hefei, 2008.45
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
An exactly divergence-free finite element method for non-isothermal flow problems
Text