We show the velocity solutions to the convective, skew-symmetric, and rotational Galerkin finite element formulations of the Navier-Stokes equations are identical if Scott-Vogelius elements are used, and thus all three formulations will the same pointwise divergence free solution velocity. A connection is then established between the formulations for grad-div stabilized Taylor-Hood elements: under mild restrictions, the formulations' velocity solutions converge to each other (and to the Scott-Vogelius solution) as the stabilization parameter tends to infinity. Thus the benefits of using Scott-Vogelius elements can be obtained with the less expensive Taylor-Hood elements, and moreover the benefits of all the formulations can be retained if the rotational formulation is used. Numerical examples are provided that confirm the theory

Bowers, Abigail L.

Cousins, Benjamin R.

Linke, Alexander

Rebholz, Leo G.

English

We show the velocity solutions to the convective, skew-symmetric, and
rotational Galerkin finite element formulations of the Navier-Stokes
equations are identical if Scott-Vogelius elements are used, and thus all
three formulations will the same pointwise divergence free solution velocity.
A connection is then established between the formulations for grad-div
stabilized Taylor-Hood elements: under mild restrictions, the formulations'
velocity solutions converge to each other (and to the Scott-Vogelius
solution) as the stabilization parameter tends to infinity. Thus the benefits
of using Scott-Vogelius elements can be obtained with the less expensive
Taylor-Hood elements, and moreover the benefits of all the formulations can
be retained if the rotational formulation is used. Numerical examples are
provided that confirm the theor

Repositorium für Naturwissenschaften und Technik

Weierstraß-Institutfür Angewandte Analysis und Stochastikim Forschungsverbund Berlin e. V.Preprint ISSN 0946 – 8633New connections between finite element formulations ofthe Navier-Stokes equationsAbigail L. Bowers 1 , Benjamin R. Cousins 2 , Alexander Linke 3 ,Leo G. Rebholz 4submitted: 09 June 20101 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: abowers@clemson.edu2 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: bcousins@clemson.edu3 Weierstrass Institutefor Applied Analysis and StochasticsMohrenstr. 3910117 Berlinemail: alexander.linke@wias-berlin.de4 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: rebholz@clemson.eduNo. 1516Berlin 20102010 Mathematics Subject Classification. 76D05 65M60.Key words and phrases. Navier-Stokes equations, rotational form, Scott-Vogelius element, strong massconservation.A. Bowers, B. Cousins, and L. Rebholz partially supported by National Science Foundation grantDMS0914478.Edited byWeierstraß-Institut für Angewandte Analysis und Stochastik (WIAS)Mohrenstraße 3910117 BerlinGermanyFax: +49 30 2044975E-Mail: preprint@wias-berlin.deWorld Wide Web: http://www.wias-berlin.de/AbstractWe show the velocity solutions to the convective, skew-symmetric, and rota-tional Galerkin finite element formulations of the Navier-Stokes equations are iden-tical if Scott-Vogelius elements are used, and thus all three formulations will thesame pointwise divergence free solution velocity. A connection is then establishedbetween the formulations for grad-div stabilized Taylor-Hood elements: under mildrestrictions, the formulations’ velocity solutions converge to each other (and to theScott-Vogelius solution) as the stabilization parameter tends to infinity. Thus thebenefits of using Scott-Vogelius elements can be obtained with the less expensiveTaylor-Hood elements, and moreover the benefits of all the formulations can be re-tained if the rotational formulation is used. Numerical examples are provided thatconfirm the theory.1 IntroductionWe consider finite element formulations of the Navier-Stokes equations (NSE), andfor simplicity we restrict to the equilibrium case with homogeneous Dirichlet boundaryconditions on a polygonal or polyhedral domain Ω,u ·∇u+∇p−ν∆u = f in Ω, (1.1)∇ ·u = 0 in Ω, (1.2)u = 0 on ∂Ω, (1.3)where u is velocity, and p is the zero mean pressure. The results herein are easilyextendable to the time-dependent case and more general boundary conditions.We prove that when computing solutions to the Galerkin finite element method for (1.1)-(1.3) using Scott-Vogelius (SV) elements (see below for a detailed description), thevelocity solutions obtained by convective, skew-symmetric, and rotational (1.8)-(1.9)formulations are identical. For general LBB stable elements, these formulations cangive very different answers [10], and thus this result shows that if one formulation isattractive for a particular problem (e.g. type of flow, use of particular preconditioners,conservation properties, etc.), then using SV elements can alleviate adverse conse-quences due to the chosen formulation. This equivalence result is then used to proveour second result, which establishes a new connection between the formulations forgrad-div stabilized Taylor-Hood (TH) elements. We prove that under mild restrictions,as the stabilization parameter tends to infinity, the velocity solutions of the formulationsall converge to the SV solution, and thus to each other.1The difference in accuracy between different discrete NSE formulations is well docu-mented, and is not restricted to finite element formulations. In [5, 6], Horiuti and Itamifound using finite difference methods that discretization errors in near wall regions canbe worse with rotational form, causing overall accuracy to suffer. Zang, in [16], foundin spectral methods, rotational form usage can lead to greater aliasing errors, althoughfor (Pk,Pk−2) spectral elements, although Wilhelm and Kleiser found that instabilitiescan occur if the rotational form is not used. In the finite element context, even thoughmatrix properties of the rotational form can allow for better preconditioning of the re-sulting linear systems [1, 2], its use of Bernoulli pressure can result in larger error inboth pressure and velocity solutions [10]. The results herein show how, in a particular(although not restrictive) setting, the best features of each of the formulations can beretained, even using the less expensive TH elements.Denote the L2(Ω) inner product by (·, ·), and with LBB stable finite dimensional spaces(Xh,Qh) ⊂ (H10(Ω),L20(Ω)). Listed below are the three common finite element formula-tions of (1.1)-(1.3): Find (uh, ph) ∈ (Xh,Qh) satisfying ∀(vh,qh) ∈ (Xh,Qh),Convective Form:(uh ·∇uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.4)(∇ ·uh,qh) = 0 (1.5)Skew-symmetric Form:(uh ·∇uh,vh)+12((∇ ·uh)uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.6)(∇ ·uh,qh) = 0 (1.7)Rotational Form:((∇×uh)×uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.8)(∇ ·uh,qh) = 0 (1.9)Each of the formulations arises from a consistent variational formulation of (1.1)-(1.2),followed by the usual Galerkin finite element discretization procedure. Skew symmetryadds an additional term to the convective formulation to enforce unconditional stability,and rotational form arises from the identityu ·∇u = (∇×u)×u+12∇|u|2 (1.10)being applied at the continuous level, then forming the Bernoulli pressure P = p+ 12|u|2.Thus the pressure ph in (1.9) approximates P, but in (1.7) and (1.5) it approximates p.1.1 Scott-Vogelius elementsScott-Vogelius elements have recently become interesting for approximating NSE so-lutions due to results of Zhang [17] and Burman and Linke [3]. Zhang showed that the2((Pk)d,Pdisck−1) element pair is LBB stable with optimal approximation properties providedk ≥ d and the mesh is constructed as a barycenter refinement of a quasi-uniform mesh(a mild restriction). With SV elements, the weak enforcement of mass conservationprovides pointwise mass conservation. This is because ∇ ·Xh ⊂ Qh, and thus one canchoose the special test function qh = ∇ ·uh (in (1.5),(1.7), or (1.9), resulting in‖∇ ·uh‖ = 0 =⇒ ∇ ·uh = 0.In general, e.g for TH elements, such a choice of test functions is not possible. We notea drawback of these elements is the use of a discontinuous pressure space, whichleads to significantly larger linear systems.Work in [3, 12, 4] has shown that using the SV element pair can provide excellent re-sults for approximating NSE solutions, as well as providing pointwise mass conserva-tion. Moreover, it can be shown that with SV elements, the velocity error is independentof the pressure error, whereas for most element choices the velocity error can be scaled(at least) by Re∗ pressureError [9, 10].Since only in the case of k ≥ d and the mesh condition are SV elements known to beLBB stable with optimal approximation properties, all references to SV elements willassume these restrictions are met.2 A equivalence theorem for the formulations with SVelementsWe now prove the equivalence of the velocities of the three formulations.Theorem 2.1. Suppose SV elements are used to solve the problems (1.4)-(1.5), (1.6)-(1.7), and (1.8)-(1.9), and the data (Ω, f ,ν) is such that the formulations have uniquesolutions. Then the velocity solutions of these problems are identical.Remark 2.1. Using the theorem, it is trivial to show the pressures are the same for theconvective and skew-symmetric formulations. However, this is not true for the rotationalform pressure since it approximates the Bernoulli pressure.Proof. The equivalence of the convective and skew-symmetric forms is trivial, sincewe have that ∇ · uh = 0 and thus 12((∇ · uh)uh,vh) = 0 in (1.6). Since this makes thetwo formulations identical, solution uniqueness of each formulation guarantees thatcomputed velocity solutions of these two formulations must be identical as well.For the equivalence of the rotational form’s velocity, since (Xh,Qh) is LBB stable, wepass to the equivalent Vh formulations, where the velocity is searched for inVh := {vh ∈ Xh : (∇ · vh,qh) = 0 ∀qh ∈ Qh},and the test functions are restricted to be in Vh. Thus considerVh convective form: (uh ·∇uh,vh)+ν(∇uh,∇vh) = ( f ,vh) ∀vh ∈Vh, (2.1)Vh rotational form: ((∇×uh)×uh,vh)+ν(∇uh,∇vh) = ( f ,vh) ∀vh ∈Vh. (2.2)3Using the vector identity (1.10), these two formulations are identical since Vh containsonly pointwise div-free functions and vh ∈ H10 :((∇×uh)×uh,vh) = (uh ·∇uh,vh)−12(∇|uh|2,vh)= (uh ·∇uh,vh)+12(|uh|2,∇ · vh)−12∫∂Ω|uh|2(vh ·n)ds= (uh ·∇uh,vh). (2.3)Therefore we have that these formulations are equivalent, and since solutions existuniquely, their respective velocity solutions must be identical.2.1 Numerical verificationWe present here a numerical example to verify Theorem 2.1. The test problem weconsider has solutionu =(2x2(x−1)2y(2y−1)(y−1)−2x(x−1)(2x−1)y2(y−1)2), p = y,and was taken from [8]. A plot of the true solution in given in figure 1.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.70.80.91Figure 1: LEFT: Velocity vector field and pressure contours of true solution for thenumerical experiments; RIGHT: The mesh used for first numerical experimentUsing ν = 0.01, we calculate f from the NSE and the selected solution. Then usingthe domain Ω = (0,1)2, use it to solve each of the three formulations (1.4)-(1.5), (1.6)-(1.7), and (1.8)-(1.9), we use ((P2)2,Pdisc1 ) elements and a Newton method. The mesh,shown in Figure 1, is a barycenter refinement of a quasi-uniform mesh, which provided3,300 total degrees of freedom. Denoting the computed velocity solutions by uconh , ussh ,and uroth , respectively, we found the SV velocity solutions are identical up to machineprecision, as predicted by Theorem 2.1:‖ussh −uconh ‖ = 1.063E −16, ‖uroth −uconh ‖ = 8.403E −1743 A connection between formulations using Taylor-HoodelementsRecent work in [10, 14] has shown that when using TH elements, accuracy and physicalfidelity can be increased through the use of grad-div stabilization, that is, adding theterm γ(∇ · uh,∇ · vh) to (1.4), (1.6), (1.8). This term, which can be derived for use inthe discrete scheme by adding 0 = −γ∇(∇ ·u) to the continuous NSE, and also arisesin the SUPG formulations, penalizes for lack of mass conservation and relaxes theeffect of the pressure error on the velocity error. Since rotational form uses Bernoullipressure, which is typically much more complex than usual pressure, this stabilizationseems natural and has proven quite effective when the rotational formulation is used[10]. Although its effect is less dramatic when used with the other formulations, grad-divstabilization can still have a significant positive effect [13, 7].Related work in [4] has shown that on a barycentric mesh and for k ≥ d, the grad-div stabilized skew-symmetric form solutions corresponding to a sequence of grad-divstabilization parameters γi → ∞, converge to the convective form SV solution, when theSV solution is unique; otherwise convergence is for a subsequence to a SV solution.This leads us to the following result.Theorem 3.1. Suppose the mesh is a barycenter refinement of a quasi-uniform mesh,and grad-div stabilized TH elements with k ≥ dim are used to solve (1.6)-(1.7), and(1.8)-(1.9). If the solution to (1.4)-(1.5) using SV elements (with the same k) is unique,then for any sequence of stabilization parameters γi →∞, the corresponding sequencesof grad-div stabilized TH velocity solutions of the skew-symmetric and rotational formu-lations converge to each other and to the SV solution.Remark 3.1. For TH elements, the skew-symmetric formulation should be used in placeof the convective formulation for unconditionally stable, as the convective formulation is(only) conditionally stable, e.g. for h small enough [11]. Thus for h small enough, Theo-rem 3.1 can be extended to include the convective formulation. Moreover, following [4],in this case one can show the sequence of modified pressures pγih := ph + γi(∇ ·uh) ofthe the convective and skew symmetric formulations converge to each other and to theSV pressure of the convective formulation.Proof. For either the rotational or skew symmetric formulations with grad-div stabiliza-tion and TH elements, taking vh = uh gives the a priori boundν‖∇uh‖2 + γi‖∇ ·uh‖2 ≤C(data), (3.1)for any γi chosen from the sequence. Thus both sequences {uss,ih }i and {urot,ih }i arebounded and therefore have convergent subsequences, and moreover their limits mustbe divergence free.Denote by U roth the limit of the convergent subsequence {{uroth }i} j, and consider the Vhformulation satisfied by an element of the subsequence urot, jh ,ν(∇urot, jh ,∇vh)+γ j(∇ ·urot, jh ,∇ ·vh)+((∇×urot, jh )×urot, jh ,vh) = ( f ,vh) ∀vh ∈VT Hh . (3.2)5Since V SVh ⊂ VT Hh , we can restrict vh ∈ VSVh in (3.2). This vanishes the grad-div term,and then using the uniform boundedness and finite dimensionality of urot, jh , we can takethe limit as γ j → ∞ to getν(∇U roth ,∇vh)+(Uroth ·∇Uroth ,vh) = ( f ,vh) ∀vh ∈VSVh . (3.3)Now by (3.1), ∇ ·U roth = 0. Since Uroth ∈ XT Hh = XSVh , and is divergence free, Uroth ∈VSVh ,and by (3.3) also satisfies the Vh convective formulation of the NSE with SV elements.Thus U roth is a rotational form SV solution, and by assumption and Theorem 2.1, Urothis the SV solution. Therefore all convergent subsequences of {urot,ih }i must converge tothe SV solution, and hence so must the entire sequence.From [4] and Theorem 2.1, we have that the entire sequence of skew-symmetric solu-tions also converges to the SV solution. Thus both the skew symmetric and rotationalforms’ sequences of velocity solutions converge to the same limit, and thus the theoremis proven.3.1 Numerical verificationTo experimentally verify Theorem 3.1, we use the same test problem as in Section 2.For γi = {0, 1, 10, 100, 1,000, 10,000}, we find solutions to all three grad-div stabilizedTH formulations, then calculate the differences between their velocity solutions. Resultsare given in Table 1, and confirm the theory.γ ‖ussh,γ −uconh,γ ‖ ‖ussh,γ −uroth,γ‖ ‖uSV −ussh,γ‖ ‖∇ ·ussh,γ‖0 1.3569E-06 2.7352E-04 2.8780E-04 1.5530E-021 1.5149E-08 6.9127E-06 6.9911E-06 2.1946E-0410 1.5472E-09 7.2828E-07 7.3620E-07 2.2585E-05100 1.5506E-10 7.3227E-08 7.4022E-08 2.2654E-061,000 1.5477E-11 7.3267E-09 7.4062E-09 2.2661E-0710,000 6.3895E-12 7.3306E-10 7.4106E-10 2.2661E-08Table 1: Convergence of the grad-div stabilized TH solutions for the different formula-tions toward each other, toward the SV solution, and toward a divergence free velocityfield as γ → ∞.4 Further extensionsIt is important to note that the results herein can easily be extended to the time-dependent case. With the restriction that the time-step be small enough to ensure solu-tion uniqueness, analogous theorems can be proven in a similar manner to those hereinfor the steady case. Moreover, these results are also extendable to time-dependent6equations that include NSE-type systems where velocity-pressure is approximated byLBB stable element pairs, e.g. the α models of turbulence, and MHD.One important extension is for general TH elements, i.e. not on a barycenter refinedmesh or k = 2 in 3D. In this case, similar proofs hold, with the limit solution living inthe divergence free subspace of the velocity spaces. However, one cannot expect thislimit solution to be accurate (see [4, 13]), as except in special cases, large grad-divstabilization parameters can over-stabilize.The results herein can also be extended to other element choices. For example, thetwo-dimensional (Pbubble2 ,Pdisc1 ) element is a common element choice that has foundsuccess in two-dimensional fluid flow computations. This element pair is LBB stable,but does not exactly conserve mass; it only provides local mass conservation, althoughfor many problems this is much better than only global mass conservation [15]. If grad-div stabilization is added, then the limiting results for the velocity solution would be thesolution from using (P2,Pdisc1 ) Scott-Vogelius elements. For, the cubic volume bubblefunction can be represented on the reference element as (xy− x2y− xy2)(u0,v0)T andits divergence is in the linear pressure space only, if u0 = 0 and v0 = 0. Therefore, alarge grad-div stabilization eliminates all the degrees of freedom corresponding to theelement bubble functions. Thus the different formulations of the nonlinear convectionterm would have to converge to each other. Similar arguments can also be made forother element choices such as (P2,P0), (P3,P1), and so forth, since in the limit theyconverge, e.g., to the Scott-Vogelius elements (P2,Pdisc1 ) and (P3,Pdisc2 ).References[1] M. Benzi and J. Liu. An efficient solver for the Navier-Stokes equations in rotationalform. SIAM Journal on Scientific Computing, 29:1959–1981, 2007.[2] M. Benzi and M. Olshanksii. An augmented Lagrangian-based approach to theOseen problem. SIAM Journal of Scientific Computing, 28(6):2095–2113, 2005.[3] E. Burman and A. Linke. Stabilized finite element schemes for incompressible flowusing Scott-Vogelius elements. Applied Numerical Mathematics, 58(11):1704–1719, 2008.[4] M. Case, V. Ervin, A. Linke, and L. Rebholz. Improving mass conservation in FEapproximations of the Navier Stokes equations using C0 velocity fields: A connec-tion between grad-div stabilization and Scott-Vogelius elements. Submitted, 2010.[5] K. Horiuti. Comparison of conservative and rotation forms in large eddy simulationof turbulent channel flow. J. Comput. Phys., 71:343Ð370, 1987.[6] K. Horiuti and T. Itami. Truncation error analysis of the rotation form of convectiveterms in the Navier-Stokes equations. J. Comput. Phys., 145:671–692, 1998.7[7] V. John and A. Kindl. Numerical studies of finite element variational multiscalemethods for turbulent flow simulations. Comput. Meth. Appl. Mech. Engrg.,199:841–852, 2010.[8] S. Kaya and B. Riviere. A two-grid stabilization method for solving the steady-stateNavier-Stokes equations. Numer. Meth. Part. Diff. Equations, 3(3):728–743, 2006.[9] W. Layton. An introduction to the numerical analysis of viscous incompressibleflows. SIAM, 2008.[10] W. Layton, C. Manica, M. Neda, M.A. Olshanskii, and L. Rebholz. On the accuracyof the rotation form in simulations of the Navier-Stokes equations. J. Comput.Phys., 228(5):3433–3447, 2009.[11] W. Layton and L. Tobiska. A two-level method with backtracking for the NavierÐS-tokes equations. SIAM J. Numer. Anal., 35:2035Ð2054, 1998.[12] A. Linke. Collision in a cross-shaped domain — A steady 2d Navier-Stokes ex-ample demonstrating the importance of mass conservation in CFD. Comp. Meth.Appl. Mech. Eng., 198(41–44):3278–3286, 2009.[13] M. Olshanskii, G. Lube, T. Heister, and J. Löwe. Grad-div stabilization and subgridpressure models for the incompressible Navier-Stokes equations. Comp. Meth.Appl. Mech. Eng., 198:3975–3988, 2009.[14] M.A. Olshanskii. A low order Galerkin finite element method for the Navier-Stokesequations of steady incompressible flow: A stabilization issue and iterative meth-ods. Comp. Meth. Appl. Mech. Eng., 191:5515–5536, 2002.[15] B. Riviere. Discontinuous Galerkin Methods for Solving Elliptic and ParabolicEquations: Theory and Implementation. SIAM, 2008.[16] T. Zang. On the rotation and skew-symmetric forms for incompressible flow simu-lations. Appl. Numer. Math., pages 27–40, 1991.[17] S. Zhang. A new family of stable mixed finite elements for the 3d Stokes equations.Math. Comp., 74(250):543–554, 2005.8

New connections between finite element formulations of the Navier-Stokes equations

Publications Server of the Weierstrass Institute for Applied Analysis and Stochastics

Weierstraß-Institutfür Angewandte Analysis und Stochastikim Forschungsverbund Berlin e. V.Preprint ISSN 0946 – 8633New connections between finite element formulations ofthe Navier-Stokes equationsAbigail L. Bowers 1 , Benjamin R. Cousins 2 , Alexander Linke 3 ,Leo G. Rebholz 4submitted: 09 June 20101 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: abowers@clemson.edu2 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: bcousins@clemson.edu3 Weierstrass Institutefor Applied Analysis and StochasticsMohrenstr. 3910117 Berlinemail: alexander.linke@wias-berlin.de4 Department of Mathematical SciencesClemson UniversityClemson SC 29634email: rebholz@clemson.eduNo. 1516Berlin 20102010 Mathematics Subject Classification. 76D05 65M60.Key words and phrases. Navier-Stokes equations, rotational form, Scott-Vogelius element, strong massconservation.A. Bowers, B. Cousins, and L. Rebholz partially supported by National Science Foundation grantDMS0914478.Edited byWeierstraß-Institut für Angewandte Analysis und Stochastik (WIAS)Mohrenstraße 3910117 BerlinGermanyFax: +49 30 2044975E-Mail: preprint@wias-berlin.deWorld Wide Web: http://www.wias-berlin.de/AbstractWe show the velocity solutions to the convective, skew-symmetric, and rota-tional Galerkin finite element formulations of the Navier-Stokes equations are iden-tical if Scott-Vogelius elements are used, and thus all three formulations will thesame pointwise divergence free solution velocity. A connection is then establishedbetween the formulations for grad-div stabilized Taylor-Hood elements: under mildrestrictions, the formulations’ velocity solutions converge to each other (and to theScott-Vogelius solution) as the stabilization parameter tends to infinity. Thus thebenefits of using Scott-Vogelius elements can be obtained with the less expensiveTaylor-Hood elements, and moreover the benefits of all the formulations can be re-tained if the rotational formulation is used. Numerical examples are provided thatconfirm the theory.1 IntroductionWe consider finite element formulations of the Navier-Stokes equations (NSE), andfor simplicity we restrict to the equilibrium case with homogeneous Dirichlet boundaryconditions on a polygonal or polyhedral domain Ω,u ·∇u+∇p−ν∆u = f in Ω, (1.1)∇ ·u = 0 in Ω, (1.2)u = 0 on ∂Ω, (1.3)where u is velocity, and p is the zero mean pressure. The results herein are easilyextendable to the time-dependent case and more general boundary conditions.We prove that when computing solutions to the Galerkin finite element method for (1.1)-(1.3) using Scott-Vogelius (SV) elements (see below for a detailed description), thevelocity solutions obtained by convective, skew-symmetric, and rotational (1.8)-(1.9)formulations are identical. For general LBB stable elements, these formulations cangive very different answers [10], and thus this result shows that if one formulation isattractive for a particular problem (e.g. type of flow, use of particular preconditioners,conservation properties, etc.), then using SV elements can alleviate adverse conse-quences due to the chosen formulation. This equivalence result is then used to proveour second result, which establishes a new connection between the formulations forgrad-div stabilized Taylor-Hood (TH) elements. We prove that under mild restrictions,as the stabilization parameter tends to infinity, the velocity solutions of the formulationsall converge to the SV solution, and thus to each other.1The difference in accuracy between different discrete NSE formulations is well docu-mented, and is not restricted to finite element formulations. In [5, 6], Horiuti and Itamifound using finite difference methods that discretization errors in near wall regions canbe worse with rotational form, causing overall accuracy to suffer. Zang, in [16], foundin spectral methods, rotational form usage can lead to greater aliasing errors, althoughfor (Pk,Pk−2) spectral elements, although Wilhelm and Kleiser found that instabilitiescan occur if the rotational form is not used. In the finite element context, even thoughmatrix properties of the rotational form can allow for better preconditioning of the re-sulting linear systems [1, 2], its use of Bernoulli pressure can result in larger error inboth pressure and velocity solutions [10]. The results herein show how, in a particular(although not restrictive) setting, the best features of each of the formulations can beretained, even using the less expensive TH elements.Denote the L2(Ω) inner product by (·, ·), and with LBB stable finite dimensional spaces(Xh,Qh)⊂ (H10 (Ω),L20(Ω)). Listed below are the three common finite element formula-tions of (1.1)-(1.3): Find (uh, ph) ∈ (Xh,Qh) satisfying ∀(vh,qh) ∈ (Xh,Qh),Convective Form:(uh ·∇uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.4)(∇ ·uh,qh) = 0 (1.5)Skew-symmetric Form:(uh ·∇uh,vh)+12((∇ ·uh)uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.6)(∇ ·uh,qh) = 0 (1.7)Rotational Form:((∇×uh)×uh,vh)− (ph,∇ · vh)+ν(∇uh,∇vh) = ( f ,vh) (1.8)(∇ ·uh,qh) = 0 (1.9)Each of the formulations arises from a consistent variational formulation of (1.1)-(1.2),followed by the usual Galerkin finite element discretization procedure. Skew symmetryadds an additional term to the convective formulation to enforce unconditional stability,and rotational form arises from the identityu ·∇u = (∇×u)×u+ 12∇|u|2 (1.10)being applied at the continuous level, then forming the Bernoulli pressure P = p+ 12 |u|2.Thus the pressure ph in (1.9) approximates P, but in (1.7) and (1.5) it approximates p.1.1 Scott-Vogelius elementsScott-Vogelius elements have recently become interesting for approximating NSE so-lutions due to results of Zhang [17] and Burman and Linke [3]. Zhang showed that the2((Pk)d,Pdisck−1) element pair is LBB stable with optimal approximation properties providedk≥ d and the mesh is constructed as a barycenter refinement of a quasi-uniform mesh(a mild restriction). With SV elements, the weak enforcement of mass conservationprovides pointwise mass conservation. This is because ∇ ·Xh ⊂ Qh, and thus one canchoose the special test function qh = ∇ ·uh (in (1.5),(1.7), or (1.9), resulting in‖∇ ·uh‖= 0 =⇒ ∇ ·uh = 0.In general, e.g for TH elements, such a choice of test functions is not possible. We notea drawback of these elements is the use of a discontinuous pressure space, whichleads to significantly larger linear systems.Work in [3, 12, 4] has shown that using the SV element pair can provide excellent re-sults for approximating NSE solutions, as well as providing pointwise mass conserva-tion. Moreover, it can be shown that with SV elements, the velocity error is independentof the pressure error, whereas for most element choices the velocity error can be scaled(at least) by Re∗ pressureError [9, 10].Since only in the case of k ≥ d and the mesh condition are SV elements known to beLBB stable with optimal approximation properties, all references to SV elements willassume these restrictions are met.2 A equivalence theorem for the formulations with SVelementsWe now prove the equivalence of the velocities of the three formulations.Theorem 2.1. Suppose SV elements are used to solve the problems (1.4)-(1.5), (1.6)-(1.7), and (1.8)-(1.9), and the data (Ω, f ,ν) is such that the formulations have uniquesolutions. Then the velocity solutions of these problems are identical.Remark 2.1. Using the theorem, it is trivial to show the pressures are the same for theconvective and skew-symmetric formulations. However, this is not true for the rotationalform pressure since it approximates the Bernoulli pressure.Proof. The equivalence of the convective and skew-symmetric forms is trivial, sincewe have that ∇ · uh = 0 and thus 12((∇ · uh)uh,vh) = 0 in (1.6). Since this makes thetwo formulations identical, solution uniqueness of each formulation guarantees thatcomputed velocity solutions of these two formulations must be identical as well.For the equivalence of the rotational form’s velocity, since (Xh,Qh) is LBB stable, wepass to the equivalent Vh formulations, where the velocity is searched for inVh := {vh ∈ Xh : (∇ · vh,qh) = 0 ∀qh ∈ Qh},and the test functions are restricted to be in Vh. Thus considerVh convective form: (uh ·∇uh,vh)+ν(∇uh,∇vh) = ( f ,vh) ∀vh ∈Vh, (2.1)Vh rotational form: ((∇×uh)×uh,vh)+ν(∇uh,∇vh) = ( f ,vh) ∀vh ∈Vh. (2.2)3Using the vector identity (1.10), these two formulations are identical since Vh containsonly pointwise div-free functions and vh ∈ H10 :((∇×uh)×uh,vh) = (uh ·∇uh,vh)−12(∇|uh|2,vh)= (uh ·∇uh,vh)+12(|uh|2,∇ · vh)−12∫∂Ω|uh|2(vh ·n)ds= (uh ·∇uh,vh). (2.3)Therefore we have that these formulations are equivalent, and since solutions existuniquely, their respective velocity solutions must be identical.2.1 Numerical verificationWe present here a numerical example to verify Theorem 2.1. The test problem weconsider has solutionu =(2x2(x−1)2y(2y−1)(y−1)−2x(x−1)(2x−1)y2(y−1)2), p = y,and was taken from [8]. A plot of the true solution in given in figure 1.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.70.80.91Figure 1: LEFT: Velocity vector field and pressure contours of true solution for thenumerical experiments; RIGHT: The mesh used for first numerical experimentUsing ν = 0.01, we calculate f from the NSE and the selected solution. Then usingthe domain Ω = (0,1)2, use it to solve each of the three formulations (1.4)-(1.5), (1.6)-(1.7), and (1.8)-(1.9), we use ((P2)2,Pdisc1 ) elements and a Newton method. The mesh,shown in Figure 1, is a barycenter refinement of a quasi-uniform mesh, which provided3,300 total degrees of freedom. Denoting the computed velocity solutions by uconh , ussh ,and uroth , respectively, we found the SV velocity solutions are identical up to machineprecision, as predicted by Theorem 2.1:‖ussh −uconh ‖= 1.063E−16, ‖uroth −uconh ‖= 8.403E−1743 A connection between formulations using Taylor-HoodelementsRecent work in [10, 14] has shown that when using TH elements, accuracy and physicalfidelity can be increased through the use of grad-div stabilization, that is, adding theterm γ(∇ · uh,∇ · vh) to (1.4), (1.6), (1.8). This term, which can be derived for use inthe discrete scheme by adding 0 =−γ∇(∇ ·u) to the continuous NSE, and also arisesin the SUPG formulations, penalizes for lack of mass conservation and relaxes theeffect of the pressure error on the velocity error. Since rotational form uses Bernoullipressure, which is typically much more complex than usual pressure, this stabilizationseems natural and has proven quite effective when the rotational formulation is used[10]. Although its effect is less dramatic when used with the other formulations, grad-divstabilization can still have a significant positive effect [13, 7].Related work in [4] has shown that on a barycentric mesh and for k ≥ d, the grad-div stabilized skew-symmetric form solutions corresponding to a sequence of grad-divstabilization parameters γi →∞, converge to the convective form SV solution, when theSV solution is unique; otherwise convergence is for a subsequence to a SV solution.This leads us to the following result.Theorem 3.1. Suppose the mesh is a barycenter refinement of a quasi-uniform mesh,and grad-div stabilized TH elements with k ≥ dim are used to solve (1.6)-(1.7), and(1.8)-(1.9). If the solution to (1.4)-(1.5) using SV elements (with the same k) is unique,then for any sequence of stabilization parameters γi →∞, the corresponding sequencesof grad-div stabilized TH velocity solutions of the skew-symmetric and rotational formu-lations converge to each other and to the SV solution.Remark 3.1. For TH elements, the skew-symmetric formulation should be used in placeof the convective formulation for unconditionally stable, as the convective formulation is(only) conditionally stable, e.g. for h small enough [11]. Thus for h small enough, Theo-rem 3.1 can be extended to include the convective formulation. Moreover, following [4],in this case one can show the sequence of modified pressures pγih := ph + γi(∇ ·uh) ofthe the convective and skew symmetric formulations converge to each other and to theSV pressure of the convective formulation.Proof. For either the rotational or skew symmetric formulations with grad-div stabiliza-tion and TH elements, taking vh = uh gives the a priori boundν‖∇uh‖2 + γi‖∇ ·uh‖2 ≤C(data), (3.1)for any γi chosen from the sequence. Thus both sequences {uss,ih }i and {urot,ih }i arebounded and therefore have convergent subsequences, and moreover their limits mustbe divergence free.Denote by U roth the limit of the convergent subsequence {{uroth }i} j, and consider the Vhformulation satisfied by an element of the subsequence urot, jh ,ν(∇urot, jh ,∇vh)+γ j(∇ ·urot, jh ,∇ ·vh)+((∇×urot, jh )×urot, jh ,vh) = ( f ,vh) ∀vh ∈V T Hh . (3.2)5Since V SVh ⊂ V T Hh , we can restrict vh ∈ V SVh in (3.2). This vanishes the grad-div term,and then using the uniform boundedness and finite dimensionality of urot, jh , we can takethe limit as γ j → ∞ to getν(∇U roth ,∇vh)+(U roth ·∇U roth ,vh) = ( f ,vh) ∀vh ∈V SVh . (3.3)Now by (3.1), ∇ ·U roth = 0. Since U roth ∈ XT Hh = XSVh , and is divergence free, U roth ∈V SVh ,and by (3.3) also satisfies the Vh convective formulation of the NSE with SV elements.Thus U roth is a rotational form SV solution, and by assumption and Theorem 2.1, U rothis the SV solution. Therefore all convergent subsequences of {urot,ih }i must converge tothe SV solution, and hence so must the entire sequence.From [4] and Theorem 2.1, we have that the entire sequence of skew-symmetric solu-tions also converges to the SV solution. Thus both the skew symmetric and rotationalforms’ sequences of velocity solutions converge to the same limit, and thus the theoremis proven.3.1 Numerical verificationTo experimentally verify Theorem 3.1, we use the same test problem as in Section 2.For γi = {0, 1, 10, 100, 1,000, 10,000}, we find solutions to all three grad-div stabilizedTH formulations, then calculate the differences between their velocity solutions. Resultsare given in Table 1, and confirm the theory.γ ‖ussh,γ −uconh,γ ‖ ‖ussh,γ −uroth,γ‖ ‖uSV −ussh,γ‖ ‖∇ ·ussh,γ‖0 1.3569E-06 2.7352E-04 2.8780E-04 1.5530E-021 1.5149E-08 6.9127E-06 6.9911E-06 2.1946E-0410 1.5472E-09 7.2828E-07 7.3620E-07 2.2585E-05100 1.5506E-10 7.3227E-08 7.4022E-08 2.2654E-061,000 1.5477E-11 7.3267E-09 7.4062E-09 2.2661E-0710,000 6.3895E-12 7.3306E-10 7.4106E-10 2.2661E-08Table 1: Convergence of the grad-div stabilized TH solutions for the different formula-tions toward each other, toward the SV solution, and toward a divergence free velocityfield as γ → ∞.4 Further extensionsIt is important to note that the results herein can easily be extended to the time-dependent case. With the restriction that the time-step be small enough to ensure solu-tion uniqueness, analogous theorems can be proven in a similar manner to those hereinfor the steady case. Moreover, these results are also extendable to time-dependent6equations that include NSE-type systems where velocity-pressure is approximated byLBB stable element pairs, e.g. the α models of turbulence, and MHD.One important extension is for general TH elements, i.e. not on a barycenter refinedmesh or k = 2 in 3D. In this case, similar proofs hold, with the limit solution living inthe divergence free subspace of the velocity spaces. However, one cannot expect thislimit solution to be accurate (see [4, 13]), as except in special cases, large grad-divstabilization parameters can over-stabilize.The results herein can also be extended to other element choices. For example, thetwo-dimensional (Pbubble2 ,Pdisc1 ) element is a common element choice that has foundsuccess in two-dimensional fluid flow computations. This element pair is LBB stable,but does not exactly conserve mass; it only provides local mass conservation, althoughfor many problems this is much better than only global mass conservation [15]. If grad-div stabilization is added, then the limiting results for the velocity solution would be thesolution from using (P2,Pdisc1 ) Scott-Vogelius elements. For, the cubic volume bubblefunction can be represented on the reference element as (xy− x2y− xy2)(u0,v0)T andits divergence is in the linear pressure space only, if u0 = 0 and v0 = 0. Therefore, alarge grad-div stabilization eliminates all the degrees of freedom corresponding to theelement bubble functions. Thus the different formulations of the nonlinear convectionterm would have to converge to each other. Similar arguments can also be made forother element choices such as (P2,P0), (P3,P1), and so forth, since in the limit theyconverge, e.g., to the Scott-Vogelius elements (P2,Pdisc1 ) and (P3,Pdisc2 ).References[1] M. Benzi and J. Liu. An efficient solver for the Navier-Stokes equations in rotationalform. SIAM Journal on Scientific Computing, 29:1959–1981, 2007.[2] M. Benzi and M. Olshanksii. An augmented Lagrangian-based approach to theOseen problem. SIAM Journal of Scientific Computing, 28(6):2095–2113, 2005.[3] E. Burman and A. Linke. Stabilized finite element schemes for incompressible flowusing Scott-Vogelius elements. Applied Numerical Mathematics, 58(11):1704–1719, 2008.[4] M. Case, V. Ervin, A. Linke, and L. Rebholz. Improving mass conservation in FEapproximations of the Navier Stokes equations using C0 velocity fields: A connec-tion between grad-div stabilization and Scott-Vogelius elements. Submitted, 2010.[5] K. Horiuti. Comparison of conservative and rotation forms in large eddy simulationof turbulent channel flow. J. Comput. Phys., 71:343Ð370, 1987.[6] K. Horiuti and T. Itami. Truncation error analysis of the rotation form of convectiveterms in the Navier-Stokes equations. J. Comput. Phys., 145:671–692, 1998.7[7] V. John and A. Kindl. Numerical studies of finite element variational multiscalemethods for turbulent flow simulations. Comput. Meth. Appl. Mech. Engrg.,199:841–852, 2010.[8] S. Kaya and B. Riviere. A two-grid stabilization method for solving the steady-stateNavier-Stokes equations. Numer. Meth. Part. Diff. Equations, 3(3):728–743, 2006.[9] W. Layton. An introduction to the numerical analysis of viscous incompressibleflows. SIAM, 2008.[10] W. Layton, C. Manica, M. Neda, M.A. Olshanskii, and L. Rebholz. On the accuracyof the rotation form in simulations of the Navier-Stokes equations. J. Comput.Phys., 228(5):3433–3447, 2009.[11] W. Layton and L. Tobiska. A two-level method with backtracking for the NavierÐS-tokes equations. SIAM J. Numer. Anal., 35:2035Ð2054, 1998.[12] A. Linke. Collision in a cross-shaped domain — A steady 2d Navier-Stokes ex-ample demonstrating the importance of mass conservation in CFD. Comp. Meth.Appl. Mech. Eng., 198(41–44):3278–3286, 2009.[13] M. Olshanskii, G. Lube, T. Heister, and J. Lo¨we. Grad-div stabilization and subgridpressure models for the incompressible Navier-Stokes equations. Comp. Meth.Appl. Mech. Eng., 198:3975–3988, 2009.[14] M.A. Olshanskii. A low order Galerkin finite element method for the Navier-Stokesequations of steady incompressible flow: A stabilization issue and iterative meth-ods. Comp. Meth. Appl. Mech. Eng., 191:5515–5536, 2002.[15] B. Riviere. Discontinuous Galerkin Methods for Solving Elliptic and ParabolicEquations: Theory and Implementation. SIAM, 2008.[16] T. Zang. On the rotation and skew-symmetric forms for incompressible flow simu-lations. Appl. Numer. Math., pages 27–40, 1991.[17] S. Zhang. A new family of stable mixed finite elements for the 3d Stokes equations.Math. Comp., 74(250):543–554, 2005.8

New connections between finite element formulations of the Navier--Stokes equations

We show the velocity solutions to the convective, skew-symmetric, and rotational Galerkin finite element formulations of the Navier–Stokes equations are identical if Scott–Vogelius elements are used, and thus all three formulations will be the same pointwise divergence free solution velocity. A connection is then established between the formulations for grad-div stabilized Taylor–Hood elements: under mild restrictions, the formulations’ velocity solutions converge to each other (and to the Scott–Vogelius solution) as the stabilization parameter tends to infinity. Thus the benefits of using Scott–Vogelius elements can be obtained with the less expensive Taylor–Hood elements, and moreover the benefits of all the formulations can be retained if the rotational formulation is used. Numerical examples are provided that confirm the theory

New connections between finite element formulations of the Navier–Stokes equations

https://archive.wias-berlin.de/servlets/MCRFileNodeServlet/wias_derivate_00001626/wias_preprints_1516.pdf

New connections between finite element formulations of the Navier--Stokes equations

Abstract

Similar works

Full text

Available Versions

Repositorium für Naturwissenschaften und Technik

Publications Server of the Weierstrass Institute for Applied Analysis and Stochastics

Publications Server of the Weierstrass Institute for Applied Analysis and Stochastics