Analysis of the hydrostatic Stokes problem and finite-element approximation in unstructured meshes

Analysis of the hydrostatic Stokes problem and finite-element approximation in unstructured meshes∗ F. Guill´en-Gonz´alez†and J.R. Rodr´ıguez-Galv´an‡ September 10, 2014 Abstract The stability of velocity and pressure mixed finite-element approximations in general meshes of the hydrostatic Stokes problem is studied, where two “inf-sup” conditions appear associated to the two constraints of the problem; namely incompressibility and hydrostatic pressure. Since these two constraints have different properties, it is not easy to choose finite element spaces satisfying both. From the analytical point of view, two main results are established; the stability of an anisotropic approximation of the velocity (using different spaces for horizontal and vertical velocities) with piecewise constant pressures, and the unstability of standard (isotropic) approximations which are stable for the Stokes problem, like the mini-element or the Taylor-Hood element. Moreover, we give some numerical simulations, which agree with the previous analytical results and allow us to conjecture the stability of some anisotropic approximations of the velocity with continuous piecewise linear pressure in unstructured meshes. Keywords. Inf-sup condition, incompressible fluids, hydrostatic pressure, primitive equations, finite elements, unstructured meshes, anisotropic Stokes equations. AMS classification (2010). 35Q35, 65N12, 65N30, 76D07. Contents 1 Introduction 2 2 Preliminaries and problem setting 4 3 Well-posedness of the hydrostatic Stokes problem: a first proof 7 4 A new proof using the saddle point framework 11 ∗The first author has been partially supported by projects MTM2009-12927 and MTM2012-32325 (DGIMEC, Spain) and the second one by the research group FQM-315 (Junta de Andaluc´ıa). †Departamento EDAN and IMUS. Universidad de Sevilla. Aptdo. 1160, 41080 Sevilla (Spain). [email protected] ‡Departamento de Matem´aticas. Universidad de C´adiz. Pol´ıgono R´ıo San Pedro S/N, 11510 Puerto Real, C´adiz (Spain). [email protected] 1 5 Conforming finite elements 13 5.1 Necessary and sufficient conditions for stability ................ 13 5.2 Sufficient conditions for “incompressible” stability .............. 14 5.3 Necessary conditions for stability ........................ 14 6 Applications and first numerical results 15 6.1 Unstable approximations of the Stokes problem with vanishing vertical viscosity ....................................... 15 6.2 Unstable approximations of the hydrostatic Stokes problem ......... 18 6.3 Stability of (P2,P1) – P0for the hydrostatic Stokes problem ......... 20 6.4 About the stability of (P1,b,P1) – P1and (P2,P1) – P1for the hydrostatic Stokes problem .................................. 23 7 Numerical simulations with different domains 27 1 Introduction The equations of geophysical fluid dynamics, modeling the motion of the ocean and the atmosphere, are derived from the conservation laws from physics. In the case of large scale ocean, it is considered as made up of a slightly compressible fluid endowed with Coriolis force and the following set of laws (see e.g. [CB09,LTW92a,LTW92b,Ped87, Tem03,TZ04]): conservation of momentum, conservation of mass (continuity equation), thermodynamics equation (conservation of energy), convection-diffusion equation for the salinity and a state equation for the density (dependent on temperature and salinity). The resulting system is too complicated for the purpose of fluid dynamics and, from a practical point of view, unaffordable. In consequence, some simplifications must be introduced. In this paper (see section 2) Cartesian coordinates in a shallow domain with the “rigid-lid” hypothesis (flat surface) is assumed. Moreover, for simplicity constant density (independent on temperature and salinity effects) is also assumed, although the results of this paper can be extended to a density-dependent case under the Boussinesq approximation. Finally, if the vertical momentum equation is reduced to the hydrostatic balance ∂zp= 0 (including the linear in depth hydrostatic pressure in the potential p), the well known primitive (or hydrostatic Navier-Stokes) equations (see (11)–(13) in Section 2) can be derived, in an oceanic domain Ω ⊂R3. This paper focus on the theoretical and numerical study of the linear steady case, arriving at the so-called hydrostatic Stokes equations in Ω: −∆νu+∇xp=f,(1) ∂zp= 0,(2) ∇ · (u, v)=0,(3) enclosed with the boundary conditions given in (8)–(9). Here (u, v):Ω−→ R3is the velocity field (u= (u1, u2)t), p: Ω −→ Ris the potential (depending on the pressure), and f= (f1, f2)tis an external force. We are denoting ∇x= (∂x, ∂y)t,∇x·u=∂xu1+∂yu2, ∇ · (u, v) = ∇x·u+∂zvand ∆ν=νx∂2 xx +νy∂2 yy +νz∂2 zz, where ν= (νx, νy, νz) is the 2 (adimensional kinematic) viscosity. For simplicity, isotropic diffusion with νx=νy=νz= 1 will be considered. The main difference between the hydrostatic Stokes equations (1)–(3) and the “classical” Stokes equations, is the disappearance of vin the vertical momentum equation, which is reduced to the hydrostatic equation ∂zp= 0. In fact, owing to this reason, vlosses its role as a coercive variable and the formulation of (1)–(3) in a saddle-point framework is not straightforward (see Section 4). Also that reason makes more difficult the numerical approximation of (1)–(3), because an additional “hydrostatic” inf-sup condition (see Section 5) will be introduced, implying in particular instability of usual finite element spaces (like Taylor-Hood P2–P1), see section 6below. In short, the purpose of this paper is to delve into the analysis of the hydrostatic Stokes (or primitive) equations of the ocean and to apply the knowledge acquired to obtain some conclusions about their approximation by finite element techniques. Specifically, the stability of some (customary or new) finite element spaces, based in usual linear or quadratic polynomial approximations in conventional (structured or not) meshes, is confirmed (and denied in some cases, notably for some standard Stokes stable finite element combinations). Hence the interest of this work: opening up the possibility of approximating primitive equations with standard P1,P1-bubble or P2polynomial approximations in usual not structured meshes and by using standard software tools. This paper is organized as follows. In section 2, we give an introduction to the hydrostatic Stokes (or primitive) equations of the ocean, we fix notation and summarize previous approaches in literature. In Section 3, a first proof of well-posedness of the hydrostatic Stokes problem (1)–(3) is provided (similar to the one presented by P. Az´erad in [Az´e94,Az´e96])). The key is the introduction of two inf-sup conditions depending on the two constraints (2) and (3), respectively. One of them plays the role of the Stokes inf-sup condition (bounding the pressure as a function of the velocity), while the other one bounds the vertical velocity as a function of the pressure. This latter “hydrostatic” inf-sup condition has no parallelism in the Stokes case and, in this sense, the hydrostatic Stokes problem needs more constraints to get stability than the Stokes one. In Section 4, we present a second proof of the well-posedness of (1)-(3), based on the saddle-point theory for mixed finite elements (see for instance [BF91]). This proof is sharper in the sense that it illustrates the role of the hydrostatic inf-sup condition to enforce the coercivity of the bilinear form underlying to the saddle-point formulation of the hydrostatic Stokes problem. In Section 5, the previous ideas are applied to the analysis of finite element approximations of the hydrostatic Stokes problem (1)-(3) and, in Section 6, different combinations of velocity/pressure finite elements are studied, both from the numerical analysis point of view and through numerical simulations. In fact, we prove that when u, v, p are approximated respectively by continuous P2,P2and P1finite elements (which will be denoted by (P2,P2) – P1), the resulting scheme is not stable (although it is stable for the Stokes problem). The same applies to (P1,b,P1,b) – P1. On the other hand, we show that the approximation of u, v, p by (P2,P1) – P0(which is Stokes-stable, see [Ste90]) is also hydrostatic-stable and we approach numerically its accuracy order in space. Moreover, we conjecture the hydrostatic stability of (P1,b,P1) – P1and (P2,P1) – P1combinations in most unstructured meshes (see section 6.4), based on the Stokes-stable results provided by the authors in [GR12]. We also check numerically its accuracy order. 3 Finally, in Section 7we exploit the benefits of using the previous stable finite elements for the hydrostatic Stokes equations (1)-(3) in unstructured meshes, developing computational tests in different domains (with or without sidewall talus) and we apply mesh adaptivity techniques. Note that these adaptive techniques would not have been easy to implement for most schemes given in literature (where a vertical structured mesh is necessary). For the purposes of scientific openness, the source code of numerical tests presented below and the graphics contained in current paper have been published, under a free license, in https://bitbucket.org/rrgalvan/hydstokes-unstructured. 2 Preliminaries and problem setting Let us start delving into the primitive equations of the ocean, fixing notation and summarizing previous approaches in literature. First of all, we fix a Cartesian system of coordinates, supposing that the domain under consideration has dimensions much shorter than the Earth’s radius (a length scale not exceeding 1000 km is acceptable [CB09]). Moreover, taking into account the disparity between horizontal and vertical scales of large scale geophysical domains, the resulting domain, denoted by Ωε⊂R3, is supposed to satisfy that the so-called aspect ratio ε=vertical scale horizontal scale is very small, for example a few kms over some thousand kms, that is ε≃10−3,10−4. We denote the spatial points as (x, zε)∈Ωε, with x= (x, y). The description of the fluid system is not complete until a relation between density and pressure is also provided. Although in the oceanic water the density is a complicated function of pressure, temperature and salinity, in this work we focus in the case of an incompressible fluid, such as pure water at ordinary pressures and temperatures, where constant density can be assumed, leaving the density-dependent case for a future work. In this case, the set of laws corresponding to the conservation of mass and momentum are completed and can be studied separately. They constitute the Navier-Stokes equations for large scale ocean, where due to the difference between horizontal and vertical scales, anisotropic eddy viscosities (see [CB09,Ped87]) are considered, denoting them by ν= (νx, νy, ε2νz), where νx,νyand νzhave O(1) scale. The unknowns are the 3D velocity field, (u, vε):Ωε×(0, T)−→ R3and the pressure, p: Ωε×(0, T)−→ R. It is then natural to perform the geometrical scaling z=zε/ε (see [AG01,BL92]), obtaining the adimensional domain (figure 1): Ω = (x, z)∈R3/x= (x, y)∈S, −D(x)< z < 0,(4) where S⊂R2is a bounded domain (the surface) and D:S→R+is a function describing the bottom depth. Here, the rigid lid hypothesis has been assumed: the true free surface of the ocean has been approximated by means of a flat averaged surface. We will denote ΓS=S× {0}the surface boundary, Γb={(x,−D(x)) /x= (x, y)∈S}the bottom boundary and Γl={(x, z)/x∈∂S, −D(x)< z < 0}the lateral walls (if it exists). 4 Figure 1: Adimensional rigid-lid domain Ω Thus, if the kinematic scaling v=vε/ε is performed, one has the next anisotropic Navier-Stokes equations in the time-space domain (0, T)×Ω: ∂tu+ (u· ∇x)u+v∂zu−∆νu+∇xp=f,(5) ε2{∂tv+ (u· ∇x)v+v∂zv−∆νv}+∂zp=g, (6) ∇x·u+∂zv= 0.(7) The term f= (f1, f2)tmodels other phenomena, like the Coriolis acceleration (in this case f=f(u)). The term ginvolves both the force due to the gravity and other effects like explicit Coriolis (of order ε, cf. [AG01]). In this paper, Coriolis effect will not be considered (because this term is linear and orthogonal to the flow, hence results presented below are not affected). Then g, represents only the force due to gravity that, as is usual, will be written in the potential form g=∂zϕ, for a given function ϕ, which can be incorporated to the pressure term (redefining p=p−ϕ). Consequently, we can take g= 0 in (6). The system is endowed with initial conditions (u, v)|t=0 = (u0, v0) and adequate boundary conditions. On the surface Γs, we will consider the rigid lid hypothesis and a winddriven condition (by the stress of a function gs), while adherence conditions are considered for uand von the bottom Γband for uon the sidewall Γl(if it exists): νz∂zu|ΓS=gs, v|ΓS= 0,(8) u|Γb∪Γl= 0, v|Γb= 0,(9) ∇v·nx|Γl= 0.(10) Here nxis the horizontal component of the outwards normal vector. The choice of the slip condition (10) on Γl, instead of a Dirichlet one, prevents boundary layers on Γl, which may appear for small εdue to the loss of regularity of the vertical velocity which follows from the ε2-terms in (6). In this work, instead of supposing εsmall, we will center directly in the limit case ε= 0, arriving at the so-called hydrostatic Navier-Stokes equations or primitive equations (of the Ocean) in (0, T)×Ω: ∂tu+ (u· ∇x)u+v∂zu−∆νu+∇xp=f,(11) ∂zp= 0,(12) ∇x·u+∂zv= 0,(13) 5 endowed with boundary conditions (8) and (9). Note that the vertical velocity vloses its regularity in the x-direction and then boundary conditions on sidewalls (10) disappear. Equations (11)-(13) can be obtained as limit of the anisotropic Navier-Stokes equations (5)–(7) when εtends to zero. This fact is justified on rigorous mathematical grounds in [BL92,AG01]. Many works about existence and regularity results of (11)–(13) are based on replacing this problem by the next (equivalent) reduced problem: find u: Ω ×(0, T)→R2, the horizontal velocity and ps:S×(0, T)→Ran (artificial) surface pressure, satisfying ∂tu+ (u· ∇x)u+v∂zu−∆νu+∇xps=f,in (0, T)×Ω,(14) ∇x· hui= 0,in (0, T)×S, (15) u|Γb∪Γl= 0, νz∂zu|ΓS=gs,(16) where vand huiare defined by the vertical integrations: v(x, z, t) = Z0 z ∇x·u(x, z0, t)dz0,hui(x, t) = Z0 −D(x) u(x, z, t)dz. The existence of weak solutions of (14)-(16) was obtained in [LTW92b]. Further extensions can be found, for instance, in [CG00,Ort04] and references therein. With respect to regularity results, the existence of strong solutions (with H2(Ω)-regularity for the horizontal velocity) was treated in [Zia95] for the linear stationary case. For the nonlinear time-dependent problem, the existence and uniqueness of local in time strong solutions is given in [GMR01]. Finally, in [CT07] it is shown that this strong solution is global in time, imposing interesting physical boundary conditions which include wind driven stress on the surface and slip friction on the bottom. In practice reduced problem (14)-(16) presents some advantages over the original hydrostatic one (11)-(13); the vertical velocity can be computed as a function of ∇x·u, hence the number of unknowns is reduced, and the dimension of the pressure is also reduced. Therefore, equations (14)-(16) has been the preferred framework for developing numerical schemes, see e.g. [CG00,CR,CR05,GR05,CGS12]. But, in order to approximate by finite elements, the reduced formulation (14)-(16) is non-local (requiring the calculus of global z-integrals) which imposes a strong vertical structure to the mesh and, in practice, makes difficult to use “standard” finite element software tools. From this point of view, the original (local) hydrostatic Navier-Stokes problem (11)–(13) has the advantage of allowing one to approximate by standard unstructured finite element meshes and to use a mesh adaptive tool. The drawback is that (as occurs for the anisotropic Navier-Stokes system (5)–(7) when ε→0) the strong anisotropy of the hydrostatic equations (11)–(13) imply that many finite element pairs which are stable for the Navier-Stokes problem, for instance P2–P1or P1,b –P1, are not stable in the hydrostatic case (see Section 6). This paper focus on the theoretic and numerical study of the linear steady case (1)–(3), which arises in linearized time-discrete schemes of (11)–(13). For example, considering a semi-implicit time scheme, implicit for the diffusion and explicit for the convection, we arrive at equations like αu−∆u+∇xp=Fin Ω, 6 where α= 1/dt, with dt > 0 the time step, and the explicit approximations are collected in F. If α= 0 this equation coincides with (1). The case α6= 0 does not introduce any additional difficulty with respect to the results of this paper. 3 Well-posedness of the hydrostatic Stokes problem: a first proof In this section, well-posedness of problem (1)-(3), (8)-(9) is shown, by means of an argument introduced by P. Azerad in [Az´e94,Az´e96]. Let us define spaces H1 b,l(Ω) = {u∈H1(Ω)2/u|Γb∪Γl= 0}, H1 z,0(Ω) = {v∈L2(Ω) / ∂zv∈L2(Ω), v|ΓS∪Γb= 0}, L2 0(Ω) = {p∈L2(Ω) /ZΩ p= 0}, and denote U=H1 b,l(Ω), V =H1 z,0(Ω), P =L2 0(Ω).(17) H1 b,l(Ω) is endowed with the L2(Ω)-norm of the gradient k∇uk=qk∂xu1k2+k∂yu1k2+k∂xu2k2+k∂yu2k2 (hereafter k · k denotes the L2(Ω)-norm) which is equivalent to the standard H1(Ω)-norm, owing to the homogeneous Dirichlet conditions for uon Γl∪Γb. Furthermore H1 z,0(Ω) is endowed with k∂zvk, which is a norm owing to the homogeneous Dirichlet condition of v on Γs∪Γband the vertical Poincar´e inequality: kvk ≤ Ck∂zvk ∀v∈H1 z,0(Ω). Definition 1. It will be said that (u, v, p)∈U×V×Pis a weak solution of the hydrostatic Stokes problem (1)-(3),(8)-(9)if (∇u,∇¯ u)−(p, ∇x¯ u) = hf,¯ ui,∀u∈U,(18) (p, ∂z¯v) = 0,∀v∈V, (19) (∇ · (u, v),¯p) = 0,∀p∈P. (20) Hereafter, (·,·) denotes the L2(Ω) inner product and h·,·i the H−1 b,l (Ω)×H1 b,l(Ω) duality. Note that, in the variational formulation (18)-(20), u∈Uis the only “coercive” variable, hence v∈Vand p∈Pplay the role of two Lagrange multipliers related to the hydrostatic and incompressibility constraints (2) and (3), respectively. It is then natural to consider the following (inf-sup) conditions: sup 06=(u,v)∈U×V (∇ · (u, v), p) k∇u, ∂zvk≥βPkpk ∀p∈P, (IS)P sup 06=p∈P (∂zv, p) kpk≥βVk∂zvk ∀v∈V, (IS)V where βP, βV>0 and k∇u, ∂zvk=pk∇uk2+k∂zvk2is the euclidean norm related to the vectorial space H1 b,l(Ω) ×H1 z,0(Ω). 7 Remark 1.Condition (IS)Pis equivalent to (see for instance [BF91]): for each p∈P, there exists (e u,ev)∈U×Vsuch that (∇ · (e u,ev), p)≥βPkpk2and k∇e u, ∂zevk≤kpk. Condition (IS)Vis equivalent to: for each v∈V, there exists ep∈Psuch that (∂zv, ep)≥βVk∂zvk2and kepk≤k∂zvk. Remark 2.Inf-sup condition (IS)Pis rather similar to the well-known inf-sup Stokes condition: there exits β > 0 such that sup 06=(u,v)∈H1 0(Ω)3 (∇ · (u, v), p) k∇u,∇vk≥βkpk ∀ p∈P. (IS) The difference is that in (IS)Pthe H1-norm of vis replaced by the H1 z-norm k∂zvk. Lemma 1. If U,Vand Pare defined by (17), then there exist two positive constants, βPand βV>0, such that (IS)Pand (IS)Vhold. In fact, βPcan be taken as the Stokes inf-sup constant and βV= 1. Proof. (IS)Pis consequence of the Stokes inf-sup (IS). Indeed, given p∈P=L2 0(Ω), sup 06=(u,v)∈U×V (∇ · (u, v), p) k∇u, ∂zvk≥sup 06=(u,v)∈H1 0(Ω)3 (∇ · (u, v), p) k∇u,∇vk≥βkpk, hence (IS)Pholds with βP=β(the Stokes inf-sup constant). On the other hand, (IS)Vcan be obtained directly using that ∂zv∈Pfor each v∈V. Indeed, taking ep=∂zv, then ep∈L2(Ω) and ZΩep=ZSZ0 −D(x) ∂zv(x, z)dz dx=ZSv(x,0) −v(x,−D(x)) dx= 0, where the last equality is a consequence of the Dirichlet conditions for vgiven in (8) and (9). Then, we have found ep∈Psuch that (ep, ∂zv) = k∂zvk2and kepk=k∂zvk, what implies (IS)V(with βV= 1). Lemma 2. The following inequality holds k∇x·uk ≤ k∇xuk,∀u∈H1 b,l(Ω). The proof of this inequality is based on the integration by parts formula (recall that x= (x, y)t)ZΩ ∂xu1∂yu2=ZΩ ∂yu1∂xu2, which is valid for regular functions vanishing on Γb∪Γl. Now, we are in position to prove existence, uniqueness and energy estimates for the hydrostatic Stokes problem (1)-(3), (8)-(9). Theorem 3. The following statements are equivalent: 8 1. U,Vand Psatisfy (IS)Pand (IS)V. 2. Problem (1)-(3),(8)-(9)is well-posed in U×V×P. In this case, there exists an unique weak solution (u, v, p)∈U×V×Pof (1)-(3),(8)-(9) and the following estimates hold: k∇uk≤kfkU0,k∂zvk ≤ 1 βVkfkU0,kpk ≤ 2 βPkfkU0,(21) where kfkU0denotes the dual norm of U=H1 b,l(Ω). Proof. We split the proof into four steps. Step 1. Rewriting the variational formulation (18)-(20). We consider the Hilbert space X=H1 b,l(Ω) ×H1 z,0(Ω) ×L2 0(Ω) endowed with the norm kχk2 X=k∇uk2+k∂zvk2+kpk2for each χ= (u, v, p)∈X. Let us define the following bilinear and linear continuous forms, for each χ= (u, v, p)∈Xand χ= (u, v, p)∈X: A(χ,χ)=(∇u,∇u)−(p, ∇ · (u, v)) −(∇ · (u, v), p),(22) L(χ) = hf,uiU0,U(23) Then, weak formulation (18)–(20) can be written as: find χ∈Xsuch that A(χ,χ) = L(χ)∀χ∈X.(24) Since A(·,·) is symmetric, for the existence and uniqueness of (18)–(20) (i.e. of (24)) it is sufficient (and necessary) to show the next generalized elliptic hypothesis of the socalled Banach-Necas-Babuska’s Theorem (see for instance [EG04, Theorem 2.6]): there exists γ > 0 such that sup 06=χ∈X A(χ,χ) kχkX ≥γkχkX∀χ∈X.(25) Step 2. Proof of the abstract inf-sup condition (25). Let χ= (u, v, p)∈X. If we consider χ1= (u, v, −p), then A(χ,χ1) = (∇u,∇u)−(∇ · (u, v), p)+(∇ · (u, v), p) = k∇uk2.(26) Owing to (IS)P, fixed p∈P, there exists (e u,ev)∈U×Vsuch that k∇e u, ∂zevk ≤ kpkand (∇ · (e u,ev), p)≥βPkpk2. Then, if we define χ2= (−e u,−ev, 0) ∈X, A(χ,χ2) = −(∇u,∇e u)+(∇ · (e u,ev), p)≥ −k∇uk k∇e uk+βPkpk2 ≥ −k∇uk kpk+βPkpk2≥βP 2kpk2−1 2βPk∇uk2.(27) Now, owing to (IS)V, fixed v∈V, there exists ep∈Psuch that kepk≤k∂zvkand (∂zv, ep)≥βVk∂zvk2. 9 equations related to (5)-(7), −∆νu+∇xp=fin Ω,(44) −ε2∆νv+∂zp=gin Ω,(45) ∇x·u+∂zv= 0 in Ω,(46) endowed with boundary conditions (8)–(10). But when ε→0, this system tends to the hydrostatic Stokes problem (1)–(3) (see [AG01] for the nonlinear unsteady case and [BL92] for the steady case). According to Theorem 6, when a finite element approximation is made for the limit problem (1)–(3), its stability requires a choice of finite element spaces verifying not only the condition (IS)P h, which is related to (IS)h, but also the constraint (IS)V h, which is a new restriction, not required in the standard Stokes problem. Hence, one can expect ill-behavior of (44)–(46) when ε→0 for some of the standard Stokes-stable FE approximations, like Taylor-Hood (which is denoted in this work by (P2,P2) – P1) and P1-bubble elements (denoted by (P1,b,P1,b) – P1). This expectation will be confirmed in Section 6.2 for the limit case ε= 0, both numerically and analytically. We present before some numerical tests showing the ill-behavior of standard Stokes-stable finite elements approximations applied to the anisotropic Stokes problem when ε > 0 is small enough. Test 1 (Anisotropic Stokes cavity test when ε→0).We have considered an unstructured mesh with 3.792 triangles in the square (0,1) ×(−1,0) ⊂R2and finite element approximations of the anisotropic Stokes problem (44)–(46) with boundary conditions (8)–(10), f= 0,g= 0 and gs= 1. The scheme has been slightly altered by introducing a standard pressure penalization technique, transforming divergence equation (46) in (∂xu+∂yv, p) + δ(p, p) = 0 ∀p∈Ph,(47) where δis a small penalization parameter (in this example, δ= 10−10). With this approximation, the condition RΩp= 0 is satisfied implicitly and it is no longer necessary to include it in the definition of Ph, taking Ph={ph∈C0(Ω) / p|T∈ P1,∀T∈ Th}. Note that this penalization technique is widely used in FE implementations of Stokes mixed problem (see e.g. examples in [PHLHM], where δ= 10−10 and δ= 10−8are chosen). We have just translated it to the “quasi-hydrostatic” case (ε≤10−3). Note that δmust not be either “too high” (for avoiding the introduction of numerical viscosity) or “too small” (for vanishing the mean of pin Ω). After doing numerical experiments for different values of δ∈(10−20,10−1), where the mean of pin Ω(i.e. RΩp/ RΩ1) and other residual indicators are tested (namely incompressibility k∇·(u, v)kL2and hydrostatic k∂zp−gkL2restrictions) we confirm that δ≃10−10 is adequate, obtaining a small mean of p(for instance 2.56815 ·10−7for ε= 10−3and similar values for other choices of ε) while the rest of indicators are not significantly affected with respect to the not penalized case (δ= 0). The experiments have been programmed using the FreeFem++ software [PHLHM]. Velocity streamlines obtained for both (P2,P2)–P1and (P1,b,P1,b)–P1finite elements are shown, respectively, in Figures 2and 3, each one for ε= 10−2,10−3and 10−4. In both cases one can see that the results are correct for ε= 10−2, but the velocity field start 16 (a) ε= 10−2(b) ε= 10−3(c) ε= 10−4 Figure 2: Velocity streamlines for (P2,P2) – P1when ε→0 (a) ε= 10−2(b) ε= 10−3(c) ε= 10−4 Figure 3: Velocity streamlines for (P1,b,P1,b) – P1when ε→0 17 (a) uh(b) vh(c) ph Figure 4: Vertical velocity instability for (P2,P2) – P1and ε= 10−4. presenting small instabilities for ε= 10−3and big oscillations for ε= 10−4. Recall that values ε= 10−3and ε= 10−4are standard scales in Oceanography (Section 2). Figure 4shows isolines of uh,vhand phfor (P2,P2)–P1and ε= 10−4, where oscillations are only appreciated for the vertical velocity (Figure 4b). This suggest that the reason for instabilities of the velocity field (in figure 2c) is the fact that (P2,P2)–P1does not verify an inf-sup condition bounding vhlike (IS)V h. The results for (P1,b,P1,b)–P1are analogue. 6.2 Unstable approximations of the hydrostatic Stokes problem Here we will show that some standard Stokes-stable finite element spaces lose their stability when ε= 0 is taken in the vertical momentum equation (45). Let us consider the hydrostatic Stokes equations (1)–(3) with boundary conditions (8)– (9). First, uh,vh,phwill be approximated by (P2,P2) – P1elements and, using the analytic results obtained in Section 5, we will prove that this combination is not hydrostaticstable. A numerical experiment confirming this fact will be shown. Finally, this theoretical analysis will be extended to (P2,P2) – P0and (P1,b,P1,b) – P1. The (P2,P2)–P1case. Let Uh={uh∈H1 b,l ∩C0(Ω) / u|T∈ P2,∀T∈ Th},(48) Vh={vh∈H1 z,0(Ω) ∩C0(Ω) / v|T∈ P2,∀T∈ Th},(49) Ph={ph∈L2 0(Ω) ∩C0(Ω) / p|T∈ P1,∀T∈ Th}.(50) We will show that the second necessary condition given in Lemma 9does not hold and then (IS)V his not verified. Note that Theorem 6states that stability of the discrete hydrostatic Stokes problem is equivalent to both (IS)P hand (IS)V hare verified. Let us consider, for example, the structured mesh shown in Figure 5. The DOF for uh are q0, x1, . . . , x8, y3, y4and q2(because uhhas not got restrictions on Γs), then Nu= 12. On the other hand, since vhhas not got restrictions on Γl, then Nv= 15 (the DOF for vhare q0, x1, ..., x8,y1, y2, y5, y6and q4, q8). Finally, the DOF for phare only x0,q1, ..., q8, 18 q0 q1 q2 q3 q4 q5q6q7 q8 x1 x2 x3 x4 x5x6x7 x8 y1 y2 y3 y4 y5 y6 Figure 5: Degrees of freedom (P2,P2) – P1in a 2 ×2 structured mesh then Np= 9 (in fact, Np= 8 considering the constraint RΩph= 0). Then Np< Nvand Lemma 9says that (IS)V hdoes not hold. This fact can be generalized to an analogous structured mesh Thdefined by n×nsubintervals, with n∈N,n > 2. Then the number of vertices is (n+ 1)2, of which (n−1)2 are in the interior of Ω. The number of midpoints is n(3n+ 2), of which 2nare in ΓS∪Γb and 2nin Γl. As consequence, Nu= (n−1)2+ (n−1) + n(3n+ 2) −3n= 4n2−2n, Nv= (n−1)2+ 2(n−1) + n(3n+ 2) −2n= 4n2, and, taking into account the null mean restriction of functions in Ph,Np= (n+ 1)2−1 = n2+ 2n. Therefore Np< Nvand (IS)V h does not hold. Then, from Theorem 6, approximation (P2,P2) – P1for uh, vh, phis not hydrostatic Stokes-stable. On the other hand, defining US h=Uh∩H1 0(Ω) and VS h=Vh∩H1 0(Ω), since combination (Uh, V S h) – Phis Stokes-stable, then, by using Lemma 8,(IS)P hholds. In order to check numerically the predictions of the numerical analysis, let us show some computational simulations. Test 2 (Hydrostatic Stokes, (P2,P2) – P1and (P1,b,P1,b) – P1, cavity test).We discretize the domain of Test 1, with h≃1/40, by both structured and unstructured meshes, approaching hydrostatic Stokes scheme (1)–(2)–(3)–(8)–(9) using finite element spaces Uh, Vh,Phgiven in (48)-(50). As in Test 1(Section 6.1) we set f= 0,g= 0,gs= 1 and pressure penalization in the divergence equation (47), with parameter δ= 10−10. The graphics obtained for velocity and pressure in the hydrostatic case, ε= 0, are indistinguishable from the ones obtained in Test 1 for ε= 10−4(Figure 4). In particular, oscillations for the vertical velocity vhare evident, but not for uhor ph. These results agree with the analysis made in this paper, since we have proven that (P2,P2)–P1does not verify (IS)V h, and then we cannot hope stability for vh. On the other hand, the pressure does not present oscillations, which it can be expected because (IS)P hholds. The cases (P2,P2)–P0and (P1,b,P1,b)–P1are analogous to (P2,P2) – P1: they both verify (IS)P h, because they are stable-Stokes, but none of them verify (IS)V h, because the space for vertical velocities is “too rich” with respect to the pressure space. For example, for (P2,P2) – P0, since Ph={ph∈L2 0(Ω) / p|T∈ P0,∀T∈ Th},(51) the number of DOF for pressure in a structured mesh defined by n×nsub-intervals 19 (similar to Figure 5) is Np= 2n2−1 . But we have already seen that Nv= 4n2, therefore Np< Nvand (IS)V hdoes not hold. With regard to (P1,b,P1,b) – P1, if we take Uh={uh∈H1 b,l ∩C0(Ω) / uh|T∈ P1,b,∀T∈ Th},(52) Vh={vh∈H1 z,0(Ω) ∩C0(Ω) / vh|T∈ P1,b,∀T∈ Th},(53) then Nv= (n−1)2+ 2(n−1) + 2n2= 3n2. But we have already seen that Np=n2+ 2n, consequently Np< Nvtherefore (P1,b,P1,b) – P1is not hydrostatic Stokes–stable. The numerical tests using (P1,b,P1,b) – P1and (P2,P2) – P0are analogous to the case (P2,P2) – P1, presenting spurious oscillations for the vertical velocity, which agree with the numerical analysis. 6.3 Stability of (P2,P1)–P0for the hydrostatic Stokes problem Let Xh=Uh×Vh×Ph, with Uhdefined in (48), phin (51) and Vh={vh∈H1 z,0(Ω) ∩C0(Ω) / v|T∈ P1,∀T∈ Th},(54) Theorem 10. (P2,P1)–P0satisfies (IS)P hand (IS)V h. Proof. (IS)P hfollows from Lemma 8because (P2,P1) – P0is Stokes-stable, see [Ste90] for a proof based on the Stenberg’s macroelement technique. (IS)V hcan be easily verified, by means of a similar argument to the second part of the proof of Lemma 1. Indeed, for each vh∈Vh, we consider eph=∂zvh∈L2 0(Ω). Then eph|T∈ P0for each T∈ Thand, ZΩeph=ZSZ0 −D(x) ∂zvh(x, z)dz dx =ZS (vh(x, 0) −vh(x, −D(x))) dx = 0, hence, eph∈Phand of course (eph, ∂zvh) = k∂zvhk2and kephk=k∂zvhk. Consequently (IS)V hholds with γV= 1 and then Xhis hydrostatic-stable. Corollary 11 (Error estimates and (P2,P1) – P0).If the exact solution (u, v, p)∈H2(Ω)× H2(Ω) ×H1(Ω), then there is a constant C > 0such that k∇(u−uh), ∂z(v−vh), p −phk ≤ C hku, v, pkH2×H2×H1. Proof. It is a direct consequence of Corollary 7. Remark 7.It is clear (see the proof of Theorem 10) that (P1,P1) – P0verifies (IS)V h. But, on the other hand, it is well known that (P1,P1) – P0is not Stokes-stable in general meshes. Using that, it can be proved that (P1,b,P1) – P0is not Stokes-stable [GR12]. Although it is not clear (see Remark 6)if(IS)P hholds or not, our numerical experiments suggest that (IS)P his not verified (see Test 3). Therefore, it seems rather probable that both (P1,P1) – P0and (P1,b,P1) – P0are not hydrostatic-stable. In this sense, (P2,P1) – P0can be qualified as a “minimal” stable finite element combination with P0pressures, both for Stokes and hydrostatic Stokes problems. 20 (a) Velocity magnitude and streamlines (b) Pressure surface plot and isolines Figure 6: Velocity and pressure for (P1,P1) – P0hydrostatic Stokes (a) Velocity streamlines (b) Pressure Figure 7: Velocity and pressure for (P2,P1) – P0hydrostatic Stokes Numerical tests Test 3 ((P1,P1) – P0cavity test).We adapt Test 2to (P1,P1)–P0. As commented in Remark 7,(P1,P1)–P0verifies (IS)V hbut not the discrete Stokes inf-sup condition (IS)h. According to our numerical tests, we conjecture that (IS)P his not verified too. For example, Figure 6shows the results for (P1,P1)–P0, where spurious oscillations for the pressure appears. By the contrary, velocity streamlines (Figure 6a) seem correct, what is coherent with the fact that (IS)V hholds. Test 4 ((P2,P1) – P0cavity test).We adapt Test 2to (P2,P1)–P0. Velocity streamlines and pressure contours have been plotted in Figure 7. Results are acceptable, presenting no spurious oscillations, what is coherent with previous theoretical analysis. Note that, although a relatively coarse mesh (h≃1/40) was employed, the qualitative behavior of flow circulation and the hydrostatic pressure agree with results obtained in literature (using integro-differential schemes in structured meshes, see e.g. [CGS12,GR05,Kim13]). Test 5 (Approximation of accuraty orders for (P2,P1) – P0).We have considered the 21 Figure 8: Velocity and pressure errors for (P2,P1) – P0hydrostatic Stokes h22−32−42−52−62−72−8 uku−uhkL21.099 1.709 1.952 1.876 1.988 2.019 ku−uhkH1 00.580 0.838 1.018 0.926 0.989 1.093 vkv−vhkL21.156 1.500 1.924 1.476 1.457 2.072 k∂zv−∂zvhkL20.770 0.894 1.032 0.941 1.004 1.081 pkp−phkL20.479 0.790 1.022 0.905 0.977 1.110 Table 1: Error orders for velocities and pressure with (P2,P1) – P0 following exact solution of the hydrostatic problem (1)–(3) in the unit squared: u(x, z) = cos(2πx) sin(2πz)−sin(2πz), v(x, z) = −u(z, x),(55) p(x, z)=2πcos(2πx).(56) Note that condition RΩp= 0 is satisfied and (u, v)|∂Ω= 0. Using (P2,P1)–P0, the approximated solution has been calculated and the absolute error for different norms has been computed, using different structured mesh sizes, obtaining the Figure 8(in logarithmic scale). On the other hand, Table 1shows concrete values representing the different error orders obtained, which are calculated from the operation log(eh2/eh1)/log(h2/h1)when h1< h2travel through the mesh sizes. Both the slope of the segments in Figure 8and the values in Table 1suggest order O(h)for (uh, vh, ph)in H1×H1 z×L2, in agreement with Corollary 11. On the other hand, order O(h2)for velocity (uh, vh)in L2is approximated numerically. Note that, though uh is taken in P2, optimal order for uh(order O(h3)in L2and O(h2)in H1) is not obtained, due to the influence of the lower approximation of vin P1. These simulations, have been carried out in structured meshes, but the results in unstructured ones are similar. 22 6.4 About the stability of (P1,b,P1)–P1and (P2,P1)–P1for the hydrostatic Stokes problem In this section we will use P1,b finite elements to approximate horizontal velocity (choosing Uhas in (48) and (52), respectively) and P1-continuous elements for the vertical velocity and for the pressure (choosing Vhas in (54) and Phas in (50)). This choice denoted as (P1,b,P1) – P1results from the elimination of some degrees of freedom from the “classical” mini-element P1,b –P1(specifically, one bubble in each triangle for the vertical velocity). Similarly (P2,P1) – P1results from the elimination of some degrees of freedom from the Taylor-Hood element P2–P1, for vertical velocity. The inf-sup condition (IS)P h It can be shown that (P1,b,P1) – P1and (P2,P1) – P1are Stokes-stable, in which we have relied so far for obtaining (IS)P h, for most usual unstructured meshes. Specifically, for those meshes that “are unstructured in the vertical direction” or, more precisely, mesh families that “do not tend to a vertical structure”, as shown in [GR12] by using the macroelement technique [Ste90]. More in detail: we say that a macroelement is x–structured (y–structured) if it can be split by an hyperplane orthogonal to the OX axis (OY axis), otherwise it is said x– unstructured (y–unstructured). For instance, the mesh shown in figure 12a is x–structured and y–unstructured. Under some geometrical assumptions for the mesh family Th, a sufficient condition for (P1,b,P1) – P1to be Stokes-stable is [GR12]: Thcan be covered by a family of macroelements, Mh, verifying the following conditions: 1. The macroelements in Mhcontain one only interior element. 2. The macroelements in Mhare y–unstructured and furthermore they “do not tend to y–structured macroelements” when h→0. Also (P2,P1) – P1is Stokes-stable in this case, as long as another condition is also verified, which is related to the fact that some kind of “weak structure” is not presented in the mesh (and thus this condition is valid in common unstructured meshes). See [GR12] for more details. In practice, in unstructured meshes built by usual mesh generators previous conditions can be assumed, as our numerical tests confirm for both (P1,b,P1) – P1and (P2,P1) – P1 elements (see tests below). The inf-sup condition (IS)V h Up to our knowledge, there is not any theoretical result showing or denying condition (IS)V h for (P1b,P1) – P1and (P2,P1) – P1. But, taking into account the following numerical tests, we can formulate some conjectures. Indeed, according to Tests 6and 7,(IS)V hcould be satisfied in unstructured meshes but with γV=γV(h). On the other hand, according to Test 8, we conjecture that (IS)V h is satisfied in meshes which are x–structured but not y–structured, with γVindependent of h. 23 (a) Velocity streamlines (b) Vertical velocity, vh(c) Pressure, ph Figure 9: (P2,P1) – P1approximation in a vertical structured mesh. Oscillations can be observed in pressure (but not in velocity). (a) Velocity streamlines (b) Vertical velocity, vh(c) Pressure, ph Figure 10: (P2,P1) – P1approximation in an unstructured mesh. No oscillations in vhand ph. Numerical tests Test 6 ((P1,b,P1) – P1and (P2,P1) – P1cavity tests).We consider the cavity test of Test 2 for either a structured 40 ×40 mesh or an unstructured one with h≃1/40. Both (P1b,P1)–P1and (P2,P1)–P1have been used with similar results. For example, in (P2,P1)–P1case, Figures 9and 10 show that the qualitative behavior of horizontal and vertical velocity are correct both for structured and unstructured meshes. If the mesh is refined, similar pictures are obtained. Also, the behavior of pressure is correct in the unstructured mesh (Figure 10c), but in the structured mesh the pressure present distinguishable oscillations (Figure 9c), according to the fact of (IS)h(and probably (IS)P h) is not verified in these types of meshes. Test 7 (Error orders for (P1,b,P1) – P1and (P2,P1) – P1in unstructured meshes).We have repeated test 5, in which an exact solution was provided. The approximated solution and the absolute error for different norms have been calculated using unstructured meshes (as they are computed by the FreeFem++ environment), with mesh sizes in the range h≃2−2, ..., 2−8. The results for (P1,b,P1)–P1and 24 (a) Errors for (P1,b,P1) – P1(b) Errors for (P2,P1) – P1 Figure 11: Velocity and pressure errors for (P1,b,P1) – P1and (P2,P1) – P1hydrostatic Stokes hmax 2−32−42−52−62−72−8 uku−uhkL22.442 1.998 2.187 1.945 2.044 2.059 ku−uhkH1 01.229 1.001 1.102 0.976 1.020 1.027 vkv−vhkL21.752 1.688 1.629 1.510 1.578 1.623 k∂zv−∂zvhkL20.463 1.117 0.490 0.327 0.508 0.715 pkp−phkL21.702 1.423 1.767 1.903 1.643 1.893 Table 2: Error orders for velocities and pressure ((P1,b,P1) – P1approximation) (P2,P1)–P1are shown in Figures 11a and 11b, respectively. They allow us to deduce some conclusions, that are supported by the error orders given in Tables 2and 3. In the case of (P1,b,P1)–P1, optimal order for uhis obtained (order O(h2)in L2and O(h)in H1). In the case of vh, optimal order is not clear, it seems that the orders obtained are slightly greater than O(h1.5)in L2and O(h0.5)in H1 z. Also the order for phin L2is in the interval (3/2,2), confirming that, as expected, (IS)P hholds for unstructured meshes. On the other hand, for (P2,P1)–P1, orders for uhare around O(h2.5)in L2and O(h1.5) in H1, then optimal orders (O(h3)in L2and O(h2)in H1) are not obtained in this case. The order for pressure seems, like in the previous case, to be greater than O(h3/2). But the behavior of vhis clearly worse than in the previous (P1,b,P1)–P1case (around order O(h)in L2and without order in H1 z). Test 8 (Error orders for (P1,b,P1) – P1and (P2,P1) – P1in a x–structured/y–unstructured mesh).We have repeated Test 7but now for the x–structured/y–unstructured given in Figure 12a. In these meshes, for (P1,b,P1)–P1, we can assure that (IS)P his verified (in usual unstructured meshes, it was enounced above), but we do not know whether (IS)V h holds or not. The results for (P1,b,P1)–P1(see Figure 12b and Table 5) suggest optimal order O(h)for vhin H1 z(Ω) and then allow us to conjecture that (IS)V hholds. 25 [Ped87] J. Pedlosky. Geophysical Fluid Dynamics. Springer-Verlag, New York, USA, 1987. [PHLHM] O. Pironneau, F. Hecht, A. Le Hyaric, and J. Morice. FreeFEM++, http: //www.freefem.org/. [Ste90] R. Stenberg. A technique for analysing finite elements methods for viscuous incompressible flow. International Journal for Numerical Methods in Fluids, 11:835–948, 1990. [Tem77] R. Temam. Navier-Stokes equations: theory and numerical analysis. Amer. Mathematical Society, 1977. [Tem03] R. Temam. Some mathematical aspects of geophysical fluid dynamic equations. Milan J. Math., 71:175–198, 2003. [TZ04] R. Temam and M. Ziane. Some mathematical problems in geophysical fluid dynamics. In Handbook of Mathematical Fluid Dynamics, volume 3, pages 535–658. Friedlander and D. Serre Editors, Elsevier, 2004. [Zia95] M. Ziane. Regularity results for stokes type systems. Applicable Analysis, 58:263–292, 1995. 32 View publication statsView publication stats