A hybridizable discontinuous Galerkin method for Stokes ow N - Description
C Nguyen a J Peraire B Cockburn Department of Aeronautics and Astronautics Massachusetts Institute of Technology Cambridge MA 02139 USA School of Mathematics University of Minnesota Minneapolis MN 55455 USA article info Article history Received 21 ID: 23795 Download Pdf
178K - views
A hybridizable discontinuous Galerkin method for Stokes ow N
C Nguyen a J Peraire B Cockburn Department of Aeronautics and Astronautics Massachusetts Institute of Technology Cambridge MA 02139 USA School of Mathematics University of Minnesota Minneapolis MN 55455 USA article info Article history Received 21
A hybridizable discontinuous Galerkin method for Stokes ow N
Download Pdf - The PPT/PDF document "A hybridizable discontinuous Galerkin me..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Presentation on theme: "A hybridizable discontinuous Galerkin method for Stokes ow N"â€” Presentation transcript:
Page 1 A hybridizable discontinuous Galerkin method for Stokes ﬂow N.C. Nguyen a, , J. Peraire , B. Cockburn Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA School of Mathematics, University of Minnesota, Minneapolis, MN 55455, USA article info Article history: Received 21 February 2009 Received in revised form 28 July 2009 Accepted 19 October 2009 Available online 1 November 2009 Keywords: Finite element methods Discontinuous Galerkin methods Hybrid/mixed methods Augmented Lagrangian Stokes ﬂow abstract
Inthispaper,weintroduceahybridizablediscontinuousGalerkinmethodforStokesﬂow.Themethodis devisedbyusingthediscontinuousGalerkinmethodologytodiscretizeavelocity–pressure–gradientfor- mulationoftheStokessystemwithappropriatechoicesofthenumericalﬂuxesandbyapplyingahybrid- ization technique to the resulting discretization. One of the main features of this approach is that it reduces the globally coupled unknowns to the numerical trace ofthe velocity and the mean of the pres- sure on the element boundaries, thereby leading to a signiﬁcant reduction in the size of the resulting
matrix.Moreover,byusinganaugmentedlagrangianmethod,thegloballycoupledunknownsarefurther reducedtothenumericaltraceofthevelocityonly.Anotherimportantfeatureisthattheapproximations of the velocity, pressure, and gradient converge with the optimal order of 1 in the -norm, when polynomialsofdegree 0areusedtorepresenttheapproximatevariables.Basedontheoptimalcon- vergence of the HDG method, we apply an element-by-element postprocessing scheme to obtain a new approximate velocity, which converges with order 2 in the -norm for 1. The postprocessing performedatthe element levelisless expensive
thanthesolutionprocedure.Numerical results are pro- vided to assess the performance of the method. 2009 Elsevier B.V. All rights reserved. 1. Introduction Inthispaper,weintroduceahybridizablediscontinuous Galer- kin (HDG) method for the Stokes system in in on where isaboundeddomainin withLipschitzboundary isaviscosity, 2 isagivensourceterm,and 2 istheDirichletboundarydata.Weassumethat isaconstantfunc- tion on and that satisﬁes the compatibility condition where denotes the outward unit normal vector on the boundary . We thus continue the study of HDG methods for second-order, symmetric
elliptic problems [5,9,11] , convection–diffusion prob- lems  , and nonlinear convection–diffusion problems  . Our long-term goal is to devise HDG methods for the incompressible Navier–Stokes equations; the consideration of the Stokes system is thus a necessary intermediate step towards this goal. Several hybridizable methods have been developed for the velocity–pressure–vorticity formulations of the Stokes ﬂow. Hybridization, as a technique to avoid the construction of diver- gence-free velocities, was introduced in  for an LDG method. The technique was then extended to a
classical mixed method for the two-dimensional and three-dimensional Stokes problems in [6,7] . This was the ﬁrst hybridization of the Stokes problem giving rise to a global system involving the degrees of freedom of unknowns (the tangential velocity and the pressure) deﬁned solely on the faces of the elements. Recently, motivated by the fact that the HDG methods for second-order symmetric elliptic equations  can be efﬁciently implemented and are more accu- rate than other DG methods [5,11] , a new class of HDG methods was introduced for the Stokes problem  . It was
shown to be hybridizable in four different ways including a tangential veloc- ity/pressure hybridization and a velocity/average pressure hybridization. Inthispaper,weapplytheapproachproposedin  todevisea new HDG method for the velocity–pressure–gradient formulation of the Stokes equations in in in on 0045-7825/$ - see front matter 2009 Elsevier B.V. All rights reserved. doi: 10.1016/j.cma.2009.10.007 Corresponding author. E-mail address: email@example.com (N.C. Nguyen). Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 Contents lists available at ScienceDirect Comput. Methods Appl. Mech.
Engrg. journal homepage: www.else vier.com/locate/cma Page 2 Here is the second-order velocity gradient tensor and is thesecond-orderidentitytensor.Theresultingmethodiscloselyre- lated to the DG method proposed in  and to the mixed method proposed in  . Indeed, the only difference between the HDG method we propose here and the DG method considered in  is the deﬁnition of the numerical traces. However, this difference al- lows us to hybridize the method, to obtain optimally convergent approximations for the velocity gradient and the pressure, and to obtain
superconvergence properties for the velocity. The method proposed in  can be hybridized to give a globally coupled sys- tem in terms of the velocity on the faces of the elements and the pressure in the whole domain,whereas our methodcan be hybrid- ized to give a system involving the velocity on the faces and the averageofthepressureontheelements.For 0,themethodpro- posedin  isnotdeﬁned.For 1,itusessmallerspacesandis able to provide approximations converging with the same order of accuracy as ours. However, it does not have the ability of handling nonmatching meshes and
variable-degree approximations typical of DG methods. TheHDGmethodappliedtotheabovesystemisdevisedasfol- lows.First,weintroducethenumericaltraceofthevelocityandthe mean of the pressure on each element as new approximate vari- ables. On each element of the triangulation, we can now express the approximate velocity, pressure and gradient in terms of the numerical trace of the velocity and the mean of pressure, which weshallreferasthelocalsolver.Furthermore,newequationshave tobeaddedtothesystemtorenderitsolvable.Theseequationsen- force the continuityof the normal componentsof the total
ﬂux. In this way, the velocity, pressure, and gradient can be expressed in an element-by-element fashion in terms of the numerical trace of the velocity and the mean of the pressure. This allows us to ob- tain the ﬁnal system involving only the numerical trace of the velocity and the mean of the pressure, thereby reducing the glob- allycoupledunknownssigniﬁcantly.Inadditiontothereductionin the number of unknowns, the HDG method possesses various advantages related to convergence properties and postprocessing of the numerical solution.
WeshowthattheHDGmethodiswell-deﬁned,thatis,thatthe numerical solution exists and is unique. We present numerical examples to examine the accuracy and convergence properties of the method. They indicate that when polynomials of degree 0areusedtorepresentthevelocity,pressure,andvelocitygra- dient, all the variables converge optimally with the order 1in the -norm for Stokes problems with smooth solution. Moreover, we use an element-by-element postprocessing to obtain a new approximate velocity which converges with order 2 in the norm for 1. Since the local postprocessing is performed at
the element level, the postprocessed velocity is less expensive to compute than the original approximate solution. Therefore, the 1 convergent pressure and 2 convergent velocity can be computed at the cost of a DG approximation using polynomials of degree Finally, we propose an efﬁcient implementation of the HDG method via the augmented Lagrangian approach by introducing a timederivativeofthepressureintothecontinuityequation.Inthis way,we can expressthe pressure in termsof the velocity,thereby eliminatingthemeanofthepressurefromthelocalsolver.Asare-
sult,wearriveatasystemintermsoftheapproximatetraceofthe velocityonly.Themainadvantageofthisimplementationstrategy is that the new system has less degrees of freedom than the origi- nal system involving both the approximate trace of the velocity and the mean of the pressure. Thepaperisorganizedasfollows.InSection weintroducethe HDGmethod for solvingthe Stokes system and apply a local post- processingtocomputeanewapproximationofthevelocity.InSec- tion we describe the detailed implementation of the HDG method. In Section we provide numerical results to assess the convergence and accuracy of the
method. Finally, in Section we present some concluding remarks. 2. The hybridizable discontinuous Galerkin method 2.1. Notation Beforedescribing theHDGmethod,we needtointroducesome notation.Wedenoteby acollectionofdisjointregularelements that partition and set f . For an element of the collection is the boundary face if the Lebesgue measure of is nonzero. For two elements and of the collection is the interior face between and if the 1 Lebesgue measure of is nonzero. We denote by and the set of interior and boundary faces, respectively. We set Let and
betheoutwardunitnormalvectorsontwoneigh- boring elements and , respectively. We use to denotethetraces of on from theinteriorof ,where , and are second-order tensorial, vectorial, and scalar functions, respectively. Then, we deﬁne the jumps as follows. For we set For ,thesetofboundaryedgesonwhich and aresingle- valued, we set where istheunitoutwardnormalto .Here denotestheusual dot product and denotes the usual dyadic or tensor product. Now, let denote the space of polynomials of degree at most onadomain andlet bethespaceofsquareintegra- ble functions on . We set , and . We introduce
discontinu- ousﬁniteelementapproximationspacesforthe gradient,velocity, and pressure as f f f In addition, we introduce a ﬁnite element approximation space for the approximate trace of the velocity f We also set f on where denotes the -projection into the space Note that consists of functions which are continuous inside the faces (or edges) and discontinuous at their borders. We further denote by the set of functions in that are con- stant on each for all elements f The mean of our approximate pressure will belong to this
space. Here the mean is deﬁned as follows. For a given function in N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 583 Page 3 , we use to deﬁne the mean of on the element bound- aries of an element as follows: Obviously, we have for any in Finally, we deﬁne various inner products for our ﬁnite element spaces.Forfunctions in ,wedenote rq if isa domain in and rq if is a domain in . Likewise, forfunctions in ,wedenote if isado- mainin and if isadomainin .Forfunc- tions in ,we denote if isa domainin and if is a domain in ; recall the
stan- dard notation tr , where tr is the trace operator. We deﬁne the volume inner products as for , and . We also deﬁne the boundary inner products as for , and 2.2. Formulation ThepointofdeparturefordevisinganHDGmethodistodeﬁne a local solver which computes approximate solutions within each element once a discrete approximation to , say , is obtained on the boundary of every mesh element. However, we note that the local Stokes problem on each element with a Dirichlet bound- ary data is not solvable. Hence, to render our local solver well-de- ﬁned we introduce a
new variable which approximates the meanofpressure ontheelementboundaries.Speciﬁcally,ourlo- cal Stokes problem on each element is the following: in in in on for the given pair . It is obvious that this local Stokes prob- lem is well-posed. We ﬁrst seek an approximation 2 such that for all h h h forall 2 .Herethenumericaltraces and are approximations to and over respectively. Furthermore, is an approximation of the mean of the pressure on the element boundary . Note that (4) is nothing but the DG discretization of the local Stokes
problem (3) . Note also thatthepresenceof inthethirdequationof (4) isnecessarytoen- forcetheidentity 0for ,since for Next, we specify the numerical trace of the form Here isthesecond-ordertensorconsistingofstabilizationparam- eters.Ithasanimportanteffectonboththestabilityandaccuracyof theresultingscheme.Theselectionofthestabilizationtensor will bediscussedinournumericalexperimentsgiveninSection .Letus brieﬂymotivatethechoiceoftheabovenumericaltrace.Wewantit to depend only on , so that we can express intermsof .AsweshalldiscussinSection 2.3 ,thisallowsus to eliminate all other variables
to obtain a weak formulation in terms of only. By adding the contribution of (4) over all the elements, and appending three suitably chosen equations, we arrive at the fol- lowing problem: ﬁnd an approximation 2 such that hi DE hi DE for all 2 . Note that the Dirichlet boundary condition has been enforced by requiring that Let us brieﬂy comment on the three added equations. The ﬁrst enforcesthecontinuityofthenormalcomponentofthenumerical trace of the total ﬂux on the interelement boundaries. The second is needed for the consistency of the method and
en- sures that the velocity can be weakly, locally divergence-free. Fi- nally, the third equation is just the average pressure constraint and is needed for the sake of well-posedness of the method. Let us brieﬂy comment on the conservative property of our numericalﬂuxes.Weobservethat issingle-valuedover since belongsto .Furthermore,if belongsto ,then the ﬁfth equation in (6) simply states that pointwise over the interior faces ; in other words, the normal component of the numerical trace is single-valued. Hence, both and are conservative ﬂuxes according to the
deﬁnition in  LetusbrieﬂymotivatethehybridizationoftheHDGmethod.At ﬁrstsight,thesystem (6) appearsunattractivesinceitinvolvestoo manyunknowns.However,itturnsoutthatbyappealingtothelo- calsolver (4) wecaneliminateallthevariables ,and toob- tain a weak formulation in terms of only. After solving for , the approximate gradient, velocity, and pressure, , can be inexpensively computed in an element-by-ele- ment fashion. Since is deﬁned on element faces and since is deﬁned as a constant function on each element boundary, the hybridization approach reduces
signiﬁcantly the number of the globally coupled unknowns. This is precisely what we are going to accomplish in the next subsection. 2.3.Characterizationofthenumericaltrace andthepressuremean Webeginbyintroducingthespeciﬁclocalsolvers.Theﬁrstlocal solver associates to the function 2 584 N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 Page 4 satisfying (4) whenweset 0and 0.Thesecondlocalsol- verassociatesto thefunction 2 satis- fying (4) whenweset 0,and .Andthethirdlocal solver associates to the function 2 satisfying
(4) when we set 0, and . More pre- cisely, we have where denotes the local solver (4) that maps into We next introduce the numerical traces We obtain the following auxiliary result whose proof is given in Appendix A Lemma 2.1. For any and , we have DE DE DE We are ready to state the characterization result. Theorem 2.1. Let be the solution of (6) We have that where 2 satisﬁes 10 and Here the forms are given by 11 for all , and Proof. Weﬁrstnotefrom (6)and(7) that satisﬁes where 2 is such that DE h 12 The desired result then follows from Lemma 2.1
and (10)–(12) 2.4. Existence and uniqueness of the numerical solution Amulti-valuedtensor issaidtobe positive-deﬁnite onaface if both branches, and ,of are positive-deﬁnite, namely, 13 When ispositive-deﬁniteon all facesof wesaythat isstrictly positive-deﬁnite and indicate this by 0on With a strictly positive-deﬁnite stabilization tensor we can provetheexistenceanduniquenessoftheHDGsolutionasfollows. We ﬁrst need to prove the well-posedness of the local solvers. Lemma 2.2. If the stabilization tensor satisﬁes the condition 0on 14 we have that
both and exist and are unique. Proof. Substituting (7) into (49) we obtain 15a DE 15b 15c 15d for all 2 . Due to the linearity, ﬁnite dimensionality,andtothefactthatthisisasquaresystem,itissuf- ﬁcientto showthat the only solutionof theabove systemfor is . Indeed, taking , and and adding the ﬁrst three equations, we get h Thisequationimpliesthat overthesimplex andthat on since we assume 0. It thus follows from (15a) that Since theaboveequationimplies isconstantover As a consequence, 0 over since 0on . Hence, from (15b) we have Since the above equation
implies is constant over We thus obtain that 0 since 0. The existence and unique- nessof canbeshowninthesamemanner.Thiscompletes the proof. Theorem 2.2. If the stabilization parameter satisﬁes the condition (14) we have that the solution of the variational formulation (10) given in Theorem 2.1 exists and is unique. Proof. The existence and uniqueness of follows if we show that the only solution of the problem (10) for and is and 0. In (10) we choose and and subtract the ﬁrst equation from the second one to obtain h i As a consequence, we can conclude that on and
on since we assume 0. N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 585 Page 5 After a simple integration by parts, the Eq. (50a) , with replaced by , now reads Thisimpliesthat isaconstanton .Asaconsequence, isalso a constant on and, since , we must have that on . Finally, we insert into the ﬁrst equation of (10) to obtain which implies 0. This completes the proof of Theorem 2.2 As a consequence of Lemma 2.2 and Theorem 2.2 , the approxi- matesolution ofthemixedformulation (6) exists and is unique. 2.5. The general form of the numerical traces We
have shown how to eliminate to obtain a weak formulationintermsof .Thekeyelementsarethelocalsol- ver (4) , the deﬁnition of the numerical trace (5) , and the jump condition 16 HereweshowthattheHDGmethodcanbeformallystatedbyafor- mulation in which the unknown variables are , and only. We proceed as follows. First, we derive an explicit expression for the numerical traces intermsof .Bythechoiceofthespace ,andassuming that the stabilization tensor is constant on each face in , the jump condition (16) implies that on Inserting (5) into the above equation, we obtain on Since both and are
positive deﬁnite, and is invertible we obtain that, on 17 Substituting this expression into (5) yields on 18 Recall that since , on the boundary faces , we have 19 Thus, we can view the HDG method as: ﬁnd an approximate solu- tion 2 such that hi DE hi 20 for all 2 , where the numerical traces, and , are given by Eqs. (17)–(19) Note that this is nothing but the weak formulation of the DG methodproposedin  .Thespacesarealsoidenticaltooursbut, aswepointedoutintheIntroduction,thenumericaltracesarediffer-
ent.Thedifferenceliesinthedeﬁnitionofthenumericaltraceforthe velocity.Indeed,incontrastwithourchoice,thenumericaltraceof thevelocityusedin  has two components:onefortheﬁrstequa- tionandanotherforthethirdequation.Theonefortheﬁrstequation lackstheterminvolvingthejumpofthetotalﬂuxandtheoneforthe third equation lacks the term involving the velocity gradient. This subtledifferenceisresponsibleforthehugedifferenceintheapprox- imationpropertiesofthemethods,aswearegoingtoseeintheSec- tionofnumericalexperiments. 2.6. Local postprocessing of the velocity We use the local
postprocessing proposed in  to obtain a new approximate velocity, of , which may converge at faster rate than the original approximation . We deﬁne the postpro- cessed approximate velocity on as the element of such that 21 To compute we need only to invert a matrix of size equal to the dimension of for each element of the triangulation Therefore, the postprocessed velocity is less expensive to compute than the original approximate velocity. 2.7. Neumann boundary condition Let us end this section by extending the method to the case when the Neumann boundary condition is en- forced in
part of the boundary ;@ . First, we require that the approximate trace belongs to f on 22 where is the Dirichlet boundary. We then replace the jump condition with h h 23 Asaresult,thebilinearforms and oftheweakformulation (10) remain unchanged, while the linear functional is now given by h 24 Hence, in order to incorporate the Neumann condition on we need only to redeﬁne the space according to (22) and modify according to (24) 3. Implementation aspects of the HDG method In this section, we describe in detail how to efﬁciently imple- ment
the HDG method via an augmented Lagrangian approach; see  andthereferencestherein.Towardsthisend,weintroduce a time derivative of the pressure into the continuity equation. In this way, we can express the pressure in term of the velocity, therebyeliminatingthemeanofthepressurefromthelocalsolver. Thus,wearriveatasystemintermsoftheapproximatetraceofthe velocityonly.Theefﬁciencyofthisimplementationstrategyliesin the fact that the new system has less degrees of freedom than the original system which involves both the approximate trace of the velocity and the mean of the pressure.
Although the HDG method can also be implemented by using the Uzawa method, we choose the augmented Lagrangian method because it is more efﬁcient thanUzawamethodforsolvingthesaddlepointsystemassociated withtheStokesproblem.Wereferto  foradetaileddiscussion. 3.1. Motivation of the method The idea of the method is to introduce an evolution problem whose limit, as time goes to inﬁnity, is nothing but the solution of the original problem. Let us show that for the continuous prob- lem. For a given initial pressure f , the evolution problem is 586 N.C. Nguyen et al./Comput.
Methods Appl. Mech. Engrg. 199 (2010) 582–597 Page 6 in 1 on 25 where is a function of and deﬁned as the solution of in 1 in 1 on 1 26 The system (25) and (26) is the time-dependent version of the ori- ginal problem (1) We proceed as follows. We ﬁrst set and , where we recall that is the solution of the steady-state original problem (1) . It then fol- lows from (25)-(26), and (1) that in 1 on 27 where is a function of and deﬁned as the solution of in 1 in 1 on 1 28 Multiplyingtheﬁrstequation (27) by andintegratingon ,we get dt where kk is the -norm
and . It follows from the equations (28) and integration by parts that Now, by the equations (28) , we have that where . Since is a self-adjoint, strongly elliptic operator, its smallest eigenvalue is strictly positive and we can write that This implies that dt and, as a consequence, that k k This shows that as time goes to inﬁnity, converges exponen- tiallyfastintimeto .Itfollowsfromthisresultandfromtheequa- tions (28) that and alsoconvergeexponentiallyfastintime to and , respectively. 3.2. Augmented Lagrangian approach The augmented Lagrangian method we use is obtained
by dis- cretizing equations (25) and (26) in time by using the backward- EulermethodandinspacebyusingtheHDGformulationdescribed Section for the spatial discretization. The main difference, how- ever, is that the local problem is now given by in 1 on in 1 in 1 on 1 29 Note that we no longer need to use the mean of the pressure since this is now a time-dependent problem. Thus, we begin deﬁning the iterative method by providing the following initial guess for the approximate pressure 30 forall .Next,givenaconstanttimestep andapressure for 1, we deﬁne the iterate as an approximation
to such that 31 Here the functions and are components of the function 2 determined by the equations DE DE 32 for all 2 , where on 33 Here and are approximations to and , respec- tively. The system of Eqs. ()()()(31)–(33) is nothing but the back- ward-Euler method for the temporal discretization and the HDG method described in Section for the spatial discretization of the continuous equations (25) and (26) 3.3. Convergence to the original HDG approximation In a similar manner as the exact solution, we show that as goes to inﬁnity, converges exponentially in time to the
original HDG approximation introduced in Section The proof of the following result can be found in Appendix B Lemma 3.1. Let and . Then the error satisﬁes k 34 where h i 35 It now follows from Lemma 3.1 and Section 3.1 that 36 where isanapproximationtothesmallesteigenvalue .Wethus obtain k k 37 which yields 38 where . This result implies that as goes to inﬁnity, converges to . The same conclusion holds for the remaining components of the approximation. Inpractice,westoptheiterationswhentherelativeerrorofthe pressure is less than a
prescribed tolerance tol , that is, when iter such that iter iter iter tol 39 N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 587 Page 7 Wenotefrom (38) that iter dependsontheartiﬁcialtimestep in such a way that iter decreases as increases. 3.4. Characterization of the approximate solution Havingestablishedthatthe iteratesoftheaugmentedLagrang- ian method converge to the approximation solution of the HDG method of Section , we now show how to efﬁciently implement the method. Towards this end, we begin by characterizing each of the iterates by
means of a hybridization technique similar to that used to hybridize the original HDG method. Weﬁrstintroducethreelocalsolvers.Theﬁrstlocalsolvermaps to 2 satisfying DE 40 for all 2 . The second local solver maps to 2 satisfying DE 41 for all 2 . The last local solver maps to 2 satisfying hi DE hi hi 42 for all 2 . Note that the local solvers here are different from the ones introduced in Section The following theorem characterizes the numerical trace as thesolutionofavariationalformulation.Theproofofthistheorem follows almost
directly from those of Theorems 2.1 and 2.2 and is thus omitted here. Theorem 3.1. Let be the solution of (31)–(33) We have that 43 where is the only function in satisfying 44 Here the forms are given by DE 45 for all Wenotethatthebilinearform issimilarto deﬁnedin (11) of the HDG method of Section except that has an additional term due to the pressure. Similarly, the functional has an addi- tional term due to the pressure. These pressure terms result from introducing the time derivative into the continuity equation. 3.5. Implementation considerations The characterization result in
Theorem 3.1 allows for an efﬁ- cient implementation of the HDG method, which we shall articu- late as follows. First, for each element , we compute the function by solvingthe local solver (40) ; thefunc- tion by solving the local solver (41) for all ele- ments of a basis of ; and the function by solving the local solver (42) for all elements of a basis of .Wethenneedtocompute foreachtimelevel To compute , we note that the matrix equation of the weak formulation (44) is of the form 46 where represents the degrees of freedom for . The matrix and vector can be formed by the usual
ﬁnite element assembly procedure once the elemental matrices and vectors are computed as follows. Let f g where dim be the set of basis functions on the faces of the boundary ofanelement .(Notethatthesebasisfunctionsarecon- structed from polynomials of degree which are deﬁned on the faces of .) The elemental matrix and vector are then given by ij DE 47 where 2 is the solution of the second local solver (42) on the element for Toverifythestoppingcriterion (39) weneedtoupdatethepres- sure as 48 where , and are calculated from (40), (41) for , and (42) for ,
respectively. If (39) does not hold, we increase 1, assemble the right-hand side of (46) , solve the system (46) for , and update the pressure according to (48) . When (39) holds, we terminate the process and compute Table 1 Implementation of the HDG method via the augmented Lagrangian approach. Implementation steps Step 1. Given tol , pick , and set 0 and Step 2. Forallof of ,compute bysolving (40) by solving (41) for forallshape functions ,1 and by solving (42) for for all shape functions f g Step 3. Calculate the elemental matrix according to (47) and form the stiffness matrix by
applying the ﬁnite element assembly procedure Step 4. Calculate the elemental vector according to (47) and form the vector by applying the ﬁnite element assembly procedure Step 5. Solve , where represents the degrees of freedom of Step 6. Compute according to (48) Step 7. If (39) does not hold, set 1 and go to Step 4 Step 8. If it does, compute according to (43) and stop 588 N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 Page 8 from (43) . The implementation of the HDG method is sum- marized in Table 1
Finally,wepointoutthedegreesoffreedomandsparsitystruc- ture of the discrete system (46) , restricting our attention to the case of a conforming triangulation (no hanging nodes). It is clear that the matrix has a block structure with square blocks of order equal to the dimension of for each face Thenumberofblockrowsandblockcolumnsisequaltothenum- ber of interior faces of the triangulation . Furthermore, on each block row, there are at most blocks that are not equal to Table 2 History of convergence of the HDG method for Degree Mesh Error Order Error Order Error Order Error Order 3.40e 0 1.16e 0
Numerical experiments WeconsidertheStokesproblemwhoseexactsolutioncoincides with the analytical solution of the incompressible Navier–Stokes equations obtained by Kovasznay in  , namely, exp cos exp sin exp where Re Re and Re is the Reynolds number. The Kovasznay ﬂow is also a solution of the Stokes problem (1) with the source term . We take Dirichlet boundary condi- tions for the velocity as the restriction of the exact solution to the domain boundary. Here the computational domain is and 1 so that the Reynolds number is Re 10.
Inourexperiments,weconsidermeshesthatarereﬁnementsof a uniform mesh of 32 congruent triangles. Each reﬁne- ment is obtained by subdividing each triangle into four congruent triangles. We say that the mesh has level if it is ob- tained from the original mesh by of these reﬁnements. On these meshes, we consider polynomials of degree to represent all the approximate variables using a nodal basis within each element, with the nodes uniformly distributed. The numerical example Fig.2.
Verticalcomponentoftheoriginalvelocity(left)andthepostprocessedvelocity(right)forpolynomialdegree 1(top), 2(middle)and 3(bottom)overthe original mesh 0. N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 591 Page 11 and meshes are taken from  to permit comparison. In all test cases, the stabilization tensor is chosen as where isapositiveconstantfunctiondeﬁnedon .Inourimple- mentation, we set tol 10 for the error tolerance. Below we present numerical results to assess the convergence and accuracy of the HDG method. We ﬁrst explore the effect of
thestabilizationparameter ontheconvergenceoftheHDGmeth- od and compare the results of the HDG method with those of the hybridized globally divergence-free LDG method  . We then demonstrate the effectiveness of the local postprocessing in improving the approximate solution. Finally, we study the effect of the artiﬁcial time step on the condition number of the dis- crete matrix and the required number of iterations. 4.1. Convergence of the HDG method We ﬁrst presentthe convergence results of the HDGmethodin Table2 for Table3 for 1,and Table4 for .We can clearly see the effect of
the stabilization parameter on the accuracy and convergence of the numerical solution. In particular, when the stabilization parameter is chosen as the approxi- mate velocity converges with the suboptimal order , while both the approximate pressure and gradient converge with the optimal order 1. On the other hand, when we set the approx- imatevelocityconvergeswiththeoptimalorder 1;however,in Fig.3. Streamlineoftheoriginalvelocity(left)and thepostprocessedvelocity(right)for polynomialdegree 1(top), 2(middle) and 3(bottom)overtheoriginal mesh 0. 592 N.C. Nguyen et al./Comput. Methods Appl. Mech.
Engrg. 199 (2010) 582–597 Page 12 thiscase,boththeapproximatepressureandgradientseemtocon- vergewithorder .Settingthestabilizationparameterto 1,all the variables converge at the optimal rate of 1. These results indicate that the optimal value of the stabilization parameter is in order of unity.It is interestingto note that the HDGmethodre- sult in approximations which converge optimally for 0 when 1, while some other DG methods such as the LDG method may not produce optimal convergent approximations in this case. Moreover,theapproximaterotationaltensor con- verges with the same order
1as for and 1. We next plot the approximate solution for different meshes and polynomial degrees. Figs. 2–6 show the two components of the approximate velocity and the streamline over the original mesh 0 and the mesh 1 for different values of and 3. We see that the approximate solution can be signiﬁcantly improved by increasing the polynomial degree or reﬁning the mesh. We now compare our results with those obtained the hybrid- ized globally divergence-free LDG method  . We note that for otherDGmethodssuchastheLDGmethod [4,12] andBassi–Rebay method [2,3] the approximate
pressure, gradient, and vorticity converge suboptimally with order . To wit, we display in Table 5 the error and order of convergence of the hybridized globally divergence-free LDG method and the HDG method with forthesameStokesproblemonthesamemeshes.Theresultsfrom thistablearetakenfrom  .Weseethattheapproximatevelocity oftheHDGmethodhasconsiderablysmallererrorsthanthatofthe hybridized globally divergence-free LDG method although they both converge with the same order 1. The approximate pres- sure of the HDG method converges with order 1 (one order Fig. 4. Horizontal component of the
original velocity (left) and the postprocessed velocity (right) for polynomial degree 1 (top), 2(middle) and 3(bottom) over the mesh 1. N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 593 Page 13 higher) and thus has signiﬁcantly smaller errors than that of the hybridized globally divergence-free LDG method. The postpro- cessedvelocityoftheHDGmethodfor 1hasroughlythesame magnitude errors as the approximate velocity of the hybridized globallydivergence-free LDGmethod for 2. These results indi- catethat theHDGmethodcanprovidethe sameaccuracywithfar less
computational cost. 4.2. Effectiveness of the local postprocessing The errorand order of convergence of the postprocessed veloc- ityarealsodisplayedin Tables2–4 .Forthecase thepost- processed velocity appears to have the same order of convergence as the originalvelocity . For the case 1 weob- serve that superconverges with order 2 for 1, which is one order higher than . It is interesting to note that for the case thepostprocessedvelocitysuperconvergeswithorder for 1, while the original velocity converges with order only. For 0, however, converges with the same order as in all cases.
Furthermore, we present in Fig. 1 the horizontal component,in Fig. 2 the vertical component, and in Fig. 3 the streamline of the original velocity and the postprocessed velocity over the original mesh 0. We also plot in Figs. 4–6 the same quantities over themesh 1.In alltheseof ﬁgures,wesee aclearimprovement of the approximations as we increase from 0to 2. Most Fig.5. Verticalcomponentoftheoriginalvelocity(left)andthepostprocessedvelocity(right)forpolynomialdegree 1(top), 2(middle)and 3(bottom)overthe mesh 1. 594 N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597
Page 14 notably, the local postprocessing improves the approximation of the velocity signiﬁcantly for 1 and 2, since the postpro- cessed velocity is clearly superior to the original velocity 4.3. Effect of the artiﬁcial time step Finally, we examine how the artiﬁcial time step affects the conditionnumberofthematrix andthenumberofiterationsre- quiredtoreachtheerrortolerance tol 10 .(Recallthatthesolu- tion procedure is started with zero pressure and terminated successfully when the relative error of the pressure is less than tol .) We deﬁne the condition
number ratio as where istheconditionnumberofthematrix .Herethecondition number is deﬁned as the ratio of the largest singular value of to thesmallestsingularvalue,whicharecomputedbyasingularvalue decomposition of We report in Tables 6 and 7 the condition number ratio and the number of iterations for convergence, respectively, as a func- tion of and for several values of .In Table 6 , we see that the ratio remains remarkably close to 1/2. As a consequence, we have that the condition number of is close to In Table 7 , we see that the number of iterations for convergence is
relativelysmallandindependentofthemeshsize andpolynomial degree . Hence, the augmented Lagrangian approach is attractive Fig. 6. Streamline of the original velocity (left) and the postprocessed velocity (right) for polynomial degree 1 (top), 2 (middle) and 3 (bottom) over the mesh 1. N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 595 Page 15 for solving the discrete Stokes system arising from the HDG discretization. 5. Conclusions In this paper, we present a hybridizable DG (HDG) method for Stokes ﬂow. We also use a local postprocessing to improve the
numerical approximations. The main features of our approach and the results of our numerical experiments can be summarized as follows. All the approximate variables converge with the optimal order 1 for an appropriate choice of the stabilization parameter, which in our particular example is of order unity. The approxi- matesolutioncanbepostprocessedtoyieldanewapproximate velocity which converges with an additional order 2 for 1. Note that these results are only observed for smooth problems. The approximate pressure and postprocessed velocity of the HDG method using polynomials of degree 1 have
accuracy comparable to the approximate pressure and velocity of the hybridizedgloballydivergence-freeLDGmethod  usingpoly- nomials of degree 1. Although the global coupled unknowns are the approximate trace of the velocity and the mean of pressure, the method can be better implemented by using the augmented Lagrangian approach since the mean of pressure is eliminated. Numerical results indicate that the number of iterations required for con- vergence is independent of both the mesh size and the polyno- mial degree. The extension of this work to the incompressible Navier–Stokes equations
constitute the subject of ongoing research. We end this paper by pointing out that the a priori error analysis of the HDG method presented here is provided in  Appendix A. Proof of Lemma 2.1 Proof. We integrate by parts the local solver (4) to obtain 49a DE 49b 49c 49d for all 2 when we consider the data . Similarly, we obtain hi 50a DE 50b hi 50c 50d for all 2 when we set the data to be . Moreover, we note that for the data . The ﬁrst identity in (8) can be derived as follows: DE DE DE DE DE DE by (50a) with and ,and (50c) with and Then DE DE by (49b) with and (49d)
.By (49a) with and (50b) with DE DE DE by (49c) with and (7) . The second identity in (8) can be de- rived as follows: DE DE DE Table 6 The condition number ratio as a function of , and Degree Mesh Artiﬁcial time step 16 .54 .52 .50 .50 .49 .51 .50 .50 .50 .50 8 .48 .47 .47 .47 .47 16 .47 .46 .46 .46 .46 32 .46 .46 .46 .46 .46 2 2 .57 .54 .53 .52 .52 4 .52 .50 .49 .49 .48 8 .49 .48 .47 .47 .47 16 .49 .47 .47 .46 .46 32 .48 .47 .47 .46 .46 Table 7 The number of iterations required for convergence as a function of , and Degree Mesh Artiﬁcial time step 16 16 12 9 7 6 4 1612976 8
1612976 16 16 12 9 8 6 32 17 12 9 8 6 2 2 16 12 9 7 6 4 1612976 8 1612986 16 17 12 9 8 6 32 17 12 9 8 6 596 N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 Page 16 by (50a) with and ,and (50c) with and Then DE DE by (50b) with (50d) ,and (7) . Finally,thelastidentityin (8) follows from the identity . This completes the proof. Appendix B. Proof of Lemma 3.1 Proof. We ﬁrst note that forall .Giventhepressureerror for 1,wehavefrom (31) and (6) that satisﬁes 51 for all . Here is the element of such that DE DE 52 for all 2 , where on Taking in
(51) , we obtain that where h . Then application of the Cauchy–Schwartz inequality yields k k Moreover, by choosing , and in (52) and summing the three equations up, we obtain h i The desired result follows from the last two equations. This com- pletes the proof. References  D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Uniﬁed analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2001) 1749–1779.  F. Bassi, S. Rebay, A high-order accurate discontinuous ﬁnite element method for the numerical
solution of the compressible Navier–Stokes equations, J. Comput. Phys. 131 (1997) 267–279.  F. Bassi, A. Crivellini, D.A. Di Pietro, S. Rebay, An artiﬁcial compressibility ﬂux for the discontinuous Galerkin solution of the incompressible Navier–Stokes equations, J. Comput. Phys. 218 (2) (2006) 794–815.  J. Carrero, B. Cockburn, D. Schtzau, Hybridized globally divergence-free LDG methods. I. The Stokes problem, Math. Comput. 75 (2006) 533–563.  B.Cockburn,B.Dong,J.Guzmn,AsuperconvergentLDG-hybridizableGalerkin methodfor second-order ellipticproblems,Math.
Comput.32(2)(2007) 233 262.  B. Cockburn, J. Gopalakrishnan, Incompressible ﬁnite elements via hybridization. Part I: the Stokes system in two space dimensions, SIAM J. Numer. Anal. 43 (4) (2005) 1627–1650.  B. Cockburn, J. Gopalakrishnan, Incompressible ﬁnite elements via hybridization. Part II: the Stokes system in three space dimensions, SIAM J. Numer. Anal. 43 (4) (2005) 1651–1672.  B. Cockburn, J. Gopalakrishnan, The derivation of hybridizable discontinuous Galerkin methodsfor Stokesﬂow,SIAMJ.Numer.Anal.47(2009) 1092–1125.  B. Cockburn, J. Gopalakrishnan,
R. Lazarov, Uniﬁed hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second- order elliptic problems, SIAM J. Numer. Anal. 47 (2009) 1319–1365.  B. Cockburn, J. Gopalakrishnan, N.C. Nguyen, J. Peraire, F.-J. Sayas, Analysis of an HDG method for Stokes ﬂow, Math. Comput., submitted for publication.  B. Cockburn, J. Guzmn, H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comput. 78 (2009) 1–24.  B.Cockburn,G.Kanschat,D.Schtzau,C.Schwab,LocaldiscontinuousGalerkin methods for
the Stokes system, SIAM J. Numer. Anal. 40 (1) (2002) 319–343.  M. Fortin, R. Glowinski, Augmented Lagrangian methods, Studies in Mathematics and Its Applications, vol. 15, North-Holland Publishing Co., Amsterdam, 1983. Applications to the numerical solution of boundary value problems, Translated from the French by B. Hunt and D.C. Spicer.  L.I.G. Kovasznay, Laminar ﬂow behind two-dimensional grid, Proc. Cambridge Philos. Soc. 44 (1948) 58–62.  N.C. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous Galerkin method for linear
convection–diffusion equations, J. Comput. Phys. 228 (2009) 3232–3254.  N.C. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous Galerkin method for nonlinear convection–diffusion equations, J. Comput. Phys. 228 (2009) 8841–8855.  Rolf Stenberg, Some new families of ﬁnite elements for the Stokes equations, Numer. Math. 56 (1990) 827–838. N.C. Nguyen et al./Comput. Methods Appl. Mech. Engrg. 199 (2010) 582–597 597