Local Projection and Grad-Div Stabilization

The Ritz Galerkin scheme (2.13)-(2.14) is prone to suffering from instabilities occurring in the numerical solution. This might be due to dominating advection or by a violation of the discrete inf-sup condition (2.16). Many stabilization techniques have been proposed in order to cope with these spurious oscillations. Besides, instabilities in the discrete velocity can occur due to a poor mass conservation of the velocity-pressure ansatz spaces at high Reynolds numbers, see [Lin09]. For instance, this becomes relevant for conforming mixed finite element methods.

For the steady Navier-Stokes equations or related auxiliary problems as the Oseen model, the widely used residual-based stabilization (RBS) methods add consistent stabilization terms to the variational formulation in the sense that the additional terms vanish for the exact strong solution. RBS methods penalize the residual of the differential equation. The non-symmetric form of the stabilization terms and the occurrence of second order deriva-tives in the residual are drawbacks regarding the efficiency of this method. An overview about RBS methods and other stabilization techniques for can be found in [RST08].

A common way is a combination of pressure-stabilizing / Petrov-Galerkin (PSPG) and Streamline-Upwind Petrov–Galerkin (SUPG) for advection together with a stabilization of the divergence constraint (grad-div). This technique is studied, for example, in [LR06a].

The SUPG method relies on testing the residual with the streamline derivative and was introduced in [BH82], PSPG was considered in [JS86,HFB86].

The so-called grad-div stabilization is an additional element-wise stabilization of the di-vergence constraint. It enhances the discrete mass conservation and reduces the effect of the pressure error on the velocity error (cf. [GLOS05,CELR11]). In case of advection dom-inated flow, it plays an important role for robustness.

Due to the mentioned drawbacks of RBS methods, other stabilization techniques are con-sidered in the literature. [BBJL07] gives an overview over different stabilization techniques, discusses the residual-based SUPG/PSPG method and presents a symmetric stabilization technique that is related to variational multiscale (VMS) methods introduced by [HMJ00]

(see also [BB06]). The key idea of VMS methods is a separation of scales into large scales, small resolved scales and small unresolved scales. VMS methods model the influence of the unresolved scales on the resolved scales, where it is often assumed that the unresolved scales only influence the small resolved scales. Only parts of the residual are used. There-fore, the consistency of the method is not guaranteed. Instead, an approximate Galerkin orthogonality holds. So convergence rates of (quasi-)optimal order can be shown. Note that many stabilization techniques can be interpreted as VMS methods.

Similar to VMS methods, local projection based stabilization (LPS) methods rely on the idea to separate the discrete function spaces into small resolved and large resolved scales and to add stabilization terms only on the small scales. The stabilization terms can be

interpreted as models for the influence of the unresolved scales that act on the smallest resolved scales. This method has the interesting features of adding solely symmetric terms to the formulation and avoiding the computation of second derivatives of the basis func-tions. The work of [MST07] provides a special interpolation operator, that is important for the numerical analysis of the LPS method. In [LRL08], a unified numerical analysis for finite element discretizations of the Oseen problem using the LPS method is given. Equal-order and inf-sup stable velocity-pressure ansatz spaces are taken into account. The LPS method considered in [BBJL07] stabilizes the pressure as well as the convective terms and is thus applicable to equal-order elements for velocity and pressure (for these elements, the discrete inf-sup condition does not hold). In [BL09], several analytical results for finite element methods for incompressible flow problems with local projection stabilization are discussed.

Since local projection and grad-div stabilization have proven useful for a large variety of flow problems, we want to apply them to the Oberbeck-Boussinesq model (2.7). It is a common procedure to transfer models introduced for e.g. the Navier-Stokes problem to non-isothermal flow. For example, in [LL12], a projection-based variational multiscale method is applied to large-eddy simulation of the Oberbeck-Boussinesq model.

Let us formulate the Oberbeck-Boussinesq model with streamline-upwind local projection stabilization (LPS SU) and grad-div stabilization. Since both velocity and temperature are considered in the advection-dominated regime, we want to add LPS in order to sta-bilize both quantities. We assume inf-sup stable discrete velocity and pressure, hence, no stabilization for the pressure is applied.

From now on, {T_h} is an admissible and shape-regular family of triangulations into d-simplices, quadrilaterals (d= 2) or hexahedra (d= 3). Let {M_h}and {L_h} be families of shape-regular macro decompositions of Ω for velocity and temperature. They represent the coarse scales in velocity and temperature. In [MST07] and later in [KL09], two approaches are mentioned for choosing the space of large scales. In the so-called two-level approach, the large scales are defined by using a coarse mesh. The coarse mesh M_h is constructed such that each macro-elementM ∈ M_h is the union of one or more neighboring elements T ∈ T_h. So M_h arises by coarsening of the original mesh T_h or, equivalently, T_h is de-rived from M_h by barycentric refinement of d-simplices or regular (dyadic) refinement of quadrilaterals and hexahedra. In the one-level LPS-approach, the coarse scales can be represented via a lower order finite elements space onT_h. Another way is to enrich the fine spaces. We can use the same abstract framework by setting M_h =T_h.L_h is constructed analogously for the temperature.

Throughout this thesis, we suppose that the following assumption holds true.

Assumption 2.2.9 (Fine and coarse triangulations).

Let {T_h}, {M_h}, {L_h} be admissible and shape-regular families of non-overlapping tri-angulations into d-simplices, quadrilaterals (d= 2) or hexahedra (d = 3). The so-called macro elements M ∈ M_h, L ∈ L_h denote the union of one or more neighboring cells T ∈ T_h: There isnT_h<∞ such that allM and L are formed as a conjunction of at most nT_h cellsT ∈ T_h. Denote by h_T, h_M and h_L the diameters of cells T ∈ T_h, M ∈ M_h and L∈ L_h, respectively. In addition, we require that there are constants C₁,C₂>0such that

hT ≤hM ≤C1hT, hT ≤hL≤C2hT ∀ T ⊂M, T ⊂L, M ∈ M_h, L∈ L_h. Denote by F_T: ˆT → T the reference mapping. We require that F_T is bijective and its Jacobian is bounded for {T_h} according to

∃c1, c2 >0 : c1h^d_T ≤ |detDFT(ˆx)| ≤c2h^d_T ∀xˆ ∈Tˆ (2.20) with constantsc₁, c₂ >0 independent of the cell diameter h_T.

Obviously, this is true for one-level methods. In case of two-level methods, it holds for barycentric or regular refinement (cf. [MT14]).

Definition 2.2.10 (Fine and coarse finite element spaces).

We denote by Y_h^u, Y_h^θ ⊂H¹(Ω)∩L^∞(Ω) finite element spaces of functions that are con-tinuous onT_h. We consider the conforming finite element spaces

V_h = [Y_h^u]^d∩V, Q_h ⊂Y_h^p∩Q, Θ_h=Y_h^θ∩Θ

for velocity, pressure and temperature, where Y_h^p is a finite element space of continuous or discontinuous functions on T_h. Moreover, letD^u_M

h ⊂[L^∞(Ω)]^d, D^θ_L

h ⊂L^∞(Ω) denote discontinuous finite element spaces onM_h foru_h and onL_h forθ_h, respectively. We set

D^u_M ={v_h|_M: vh∈D^u_Mh}, D^θ_L={ψ_h|_L: ψh∈D_L^θ_h}.

In the following, we often write for combinations of finite element spaces (Vh/D^u_Mh)∧Qh∧(Θ_h/D^θ_Lh), or (Vh/D^u_M)∧Qh∧(Θ_h/D_L^θ).

If no LPS is applied, we omit the respective coarse space in the above notation.

Definition 2.2.11 (Fluctuation operators).

ForM ∈ M_h and L∈ L_h, let π_M^u : [L²(M)]^d→D^u_M, π^θ_L:L²(L)→D^θ_L be the orthogonal

L²-projections onto the respective macro spaces. The so-called fluctuation operators are defined by

κ^u_M: [L²(M)]^d→[L²(M)]^d, κ^θ_L:L²(L)→L²(L), κ^u_M :=Id−π^u_M, κ^θ_L:=Id−π^θ_L.

E_h is defined as the set of inner element faces E 6∈ ∂Ω of T_h. We denote by h_E the diameter of the face E ∈ E_h. For two cells TE and T_E⁰ shared by E, let nE be the unit normal vector pointing from TE into T_E⁰. For piecewise smooth functions wh, we denote by [w_h]_E := (w_h|_TE)|_E−(w_h|_T⁰

E)|_E the jump over the face E. Note that this is unique up to a sign.

Let uh ∈Vh. For all macro elements M ∈ M_h and L ∈ L_h, we denote the element-wise averaged streamline directions by u_M ∈R^d,u_L∈R^d. For instance, we can choose

uM := 1

|M|

Z

M

u_h(x) dx, uL:= 1

|L|

Z

L

u_h(x) dx.

The semi-discrete stabilized Oberbeck-Boussinesq model reads:

Find (u_h, p_h, θ_h) : (0, T)→V_h×Q_h×Θ_h such that for all (v_h, q_h, ψ_h)∈V_h×Q_h×Θ_h: (∂tuh,vh) + (ν∇u_h,∇v_h) +cu(uh;uh,vh)−(ph,∇ ·vh) + (∇·u_h, qh)

+(βgθ_h,v_h) +su(u_h;u_h,v_h) +t_h(u_h;u_h,v_h) +i_h(p_h, q_h) = (fu,v_h), (2.21) (∂tθ_h, ψ_h) + (α∇θ_h,∇ψ_h) +c_θ(u_h;θ_h, ψ_h) +s_θ(u_h;θ_h, ψ_h) = (f_θ, ψ_h) (2.22) with the streamline-upwind (SUPG)-type stabilizations su,sθ, the grad-div stabilization t_h and the pressure jump stabilizationi_h according to

su(w_h;u,v) := ^X

M∈M_h

τ_M^u(wM)(κ^u_M((wM · ∇)u), κ^u_M((wM· ∇)v))_M, s_θ(w_h;θ, ψ) := ^X

L∈L_h

τ_L^θ(w_L)(κ^θ_L((w_L· ∇)θ), κ^θ_L((w_L· ∇)ψ))_L, th(wh;u,v) := ^X

M∈M_h

γM(wM)(∇ ·u,∇ ·v)M, i_h(p, q) := ^X

E∈E_h

φE([p]E,[q]E)_E

with non-negative stabilization parameters τ_M^u,τ_L^θ,γM,φE. Note that the pressure stabi-lization takes care of pressure jumps in case a discontinuous ansatz space Q_h is chosen.

The set of stabilization parametersτ_M^u(u_h),τ_L^θ(u_h),γ_M(u_h), andφ_E has to be determined later on.

Taking the discrete inf-sup stability from Assumption 2.2.1into account, we can look for weakly solenoidal solutionsu_h ∈V^div_h and reformulate the stabilized Oberbeck-Boussinesq model (2.21)-(2.22) as follows:

Find (uh, ph, θh) : (0, T)→V^div_h ×Qh×Θ_h such that for all (vh, qh, ψh)∈Vh×Qh×Θ_h: (∂_tu_h,v_h) + (ν∇u_h,∇v_h) +c_u(u_h;u_h,v_h)−(p_h,∇ ·v_h) + (βgθ_h,v_h)

+s_u(u_h;u_h,v_h) +t_h(u_h;u_h,v_h) +i_h(p_h, q_h) = (f_u,v_h), (∂tθh, ψh) + (α∇θ_h,∇ψ_h) +cθ(uh;θh, ψh) +sθ(uh;θh, ψh) = (fθ, ψh).

We utilize this formulation later for the analysis.