Equal-Order Stabilized Finite Element Approximation of the p-Stokes Equations on Anisotropic Cartesian Meshes

Josefin Ahlkrona; Malte Braack

doi:10.1515/cmam-2018-0260

Open Access Published by De Gruyter February 15, 2019

Equal-Order Stabilized Finite Element Approximation of the p-Stokes Equations on Anisotropic Cartesian Meshes

Josefin Ahlkrona and Malte Braack

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2018-0260

Abstract

The p-Stokes equations with power-law exponent p ∈ ( 1 , 2 ) describes non-Newtonian, shear-thinning, incompressible flow. In many industrial applications and natural settings, shear-thinning flow takes place in very thin domains. To account for such anisotropic domains in simulations, we here study an equal-order bi-linear anisotropic finite element discretization of the p-Stokes equations, and extend a non-linear Local Projection Stabilization to anisotropic meshes. We prove an a priori estimate and illustrate the results with two numerical examples, one confirming the rate of convergence predicted by the a-priori analysis, and one showing the advantages of an anisotropic stabilization compared to an isotropic one.

Keywords: Anisotropic Meshes; Finite Element Method; Shear-Thinning; Local Projection Stabilization

MSC 2010: 65N30; 65M12; 76A05; 76D07

1 Introduction

Shear-thinning flow in thin domains is common in industry and nature, and numerical simulations is an important tool to understand such problems. The viscosity of shear-thinning (or pseudo-plastic) fluids increases when shear rates decrease, which renders numerical simulations both more computationally demanding, and harder to analyze. Thin domains further increase the complexity of the problem, and tends to reduce the accuracy and stability of numerical methods. Examples of shear thinning flows taking place on thin domains in industry are paint drying on a wall, hot (metal) rolling, and polymer coating processes. An interesting example from nature is the flow of ice in ice sheets, such as the Greenland Ice Sheet or the Antarctic Ice Sheet. There are also non-Newtonian fluids for which the viscosity decreases for decreasing shear, so called shear-thickening (or dilatant) fluids. As such those are much less common than shear-thinning fluids, we here focus on the latter.

The equations that model shear-dependent (steady) flow are the p-Stokes equations,

(1.1) - ∇ ⋅ 𝐒 ⁢ ( 𝐃𝐯 ) + ∇ ⁡ π = 𝐟 , ∇ ⋅ 𝐯 = 0 , } in Ω ,

where 𝐯 = ( v x , v y ) is the velocity field in a Cartesian coordinate system ( x , y ) , π is the pressure, 𝐒 is the shear-dependent extra stress tensor, and 𝐟 a given forcing. Here, the two-dimensional domain Ω ∈ ℝ 2 is open and bounded with polygonal boundary ∂ ⁡ Ω (for simplicity). We will consider Dirichlet condition on the entire boundary ∂ ⁡ Ω and enforce a pressure of zero-mean value,

𝐯 = 𝐠 on ⁢ ∂ ⁡ Ω , ∫ Ω π = 0 .

In order for the equations to earn the name p-Stokes, the stress tensor 𝐒 must exhibit a so-called ( p , ϵ ) -structure (defined in Assumption 1 in Section 2.5). A common example of a stress-tensor with ( p , ϵ ) structure is

(1.2) 𝐒 ⁢ ( ∇ ⁡ 𝐯 ) = μ 0 ⁢ ( ϵ 2 + | 𝐃𝐯 | 2 ) p - 2 2 ⁢ 𝐃𝐯 ,

where 𝐃𝐯 = ( ∇ ⁡ 𝐯 ) sym = 1 2 ⁢ ( ∇ ⁡ 𝐯 + ( ∇ ⁡ 𝐯 ) T ) is the strain rate tensor. Note that S ⁢ ( ∇ ⁡ 𝐯 ) = S ⁢ ( 𝐃𝐯 ) by definition. If ϵ = 0 , equation (1.2) describes a so-called power-law fluid, and for ϵ > 0 it represents a Carreau-fluid. The viscosity-parameter μ 0 scales the magnitude of the viscosity μ := μ 0 ⁢ ( ϵ 2 + | 𝐃𝐯 | 2 ) p - 2 p , we will set μ 0 = 1.0 for brevity of the presentation. The exponent p ∈ ( 1 , ∞ ) determines the non-linearity of the material. If p = 2 , the viscosity μ is constant ( μ = μ 0 ), and the system of (1.1) is simply the common steady Stokes problem describing Newtonian fluids. If p ∈ ( 1 , 2 ) , the fluid is shear-thinning, and for p > 2 the fluid is shear-thickening. Hence we consider the range p ∈ ( 1 , 2 ) in this work. Parts of the analysis is identical for p ∈ ( 2 , ∞ ) , while other parts would require extensions.

The finite element discretization of the p-Stokes equations on isotropic meshes has been studied since the early 1990s, see [3] and [23]. The analysis is closely related to the analysis of the p-Laplacian, and involves interpolation estimates in quasi-norms of the type in [3]. Such estimates can be found if the problem is transferred from classical Sobolev spaces to Orlicz–Sobolev spaces (see [14]). Using this approach, optimal a priori estimates for the shear-thinning case were derived independently by [6] and [18, 19]. In [6], inf-sup stable elements are used, while in [18, 19] stabilized equal-order finite elements are employed. The stabilization there is a Local Projection Stabilization (LPS) [5] modified to suit the p-Stokes equations. Further, we like to refer to the work of Diening, Kreuzer and Süli [13], where even more general models of rheology are discretized with inf-sup stable finite elements.

Since we are studying thin domains, the x-direction of Ω is assumed to be much wider than the y-direction. In order to resolve thin domains in a computationally efficient manner, a finite element discretization should employ anisotropic meshes. Anisotropic meshes are also useful for capturing a type of boundary layer that may occur in shear-thinning free surface flow due to the varying viscosity. These boundary layers are in general not as sharp as boundary layers in convective problems, but nevertheless need to be resolved [1, 20].

The introduction of anisotropic meshes complicates both analysis and computations. Anisotropic meshes are quite well-studied for Newtonian fluids, see e.g. [2, 22, 7]. In [10] the LPS stabilization was extended to anisotropic meshes for Newtonian fluids. However, in the non-Newtonian case, there is to our knowledge no study on anisotropic meshes. In this work we extend the stabilization of [19] to anisotropic meshes by combining the anisotropic stabilization term in [10] with the idea behind the non-linear one in [18]. We focus on (bi-)linear elements, as the regularity of the solutions are limited, which may render high-order elements suboptimal [18]. In order to derive a priori estimates, we extend the interpolation estimates in the Orlicz–Sobolev spaces of [14] to anisotropic meshes.

An example of an application where a stabilization term especially adapted to both a non-linear rheology and anisotropic meshes would be of benefit is ice sheet modelling. Ice sheets are very thin compared to the vast areas they cover, and under long time scales ice sheet dynamics can be described as a shear-thinning thin film flow with p ≈ 1.3 . Stabilized low-order equal finite elements on anisotropic meshes is a common approach in this field, but the lack of specialized stabilization terms is problematic [17].

The paper is structured as follows. In Section 2 the finite element formulation of the p-Stokes equations on anisotropic meshes is introduced. Some fundamental tools for working with the p-Stokes system are also described here, such as N-functions that are needed to transfer to Orlicz–Sobolev spaces. Section 3 is dedicated to the error analysis. First, a modified interpolation operator is constructed and interpolation estimates are shown for anisotropic meshes, including the estimates in quasi-norms that are needed for the p-Stokes system. Then the well-posedness of the problem is shown and a priori estimates for the velocity and pressure are derived. Numerical experiments in Section 4 verify the error analysis and compares the performance of isotropic stabilization to that of anisotropic stabilization. The paper concludes in Section 5 with a summary.

2 Finite Element Formulation

2.1 Notation

The spatial dimension is denoted by d, and as we work in two dimensions d = 2 , although a generalization to d = 3 with almost uniform spatial horizontal mesh sizes is possible (see details below). The scalar product between two vectors 𝐪 and 𝐩 is 𝐪 ⋅ 𝐩 . Between two tensors 𝐐 : 𝐏 ∈ ℝ d × d it is defined as 𝐐 : 𝐏 = ∑ i , j = 1 d Q i ⁢ j ⁢ P i ⁢ j , and we use | 𝐐 | = ( 𝐐 : 𝐐 ) 1 2 . The symmetric part of a tensor 𝐐 is 𝐐 sym := 1 2 ⁢ ( 𝐐 + 𝐐 T ) . We will sometimes use a multi-index notation with a multi-index α = ( α 1 , α 2 ) to write partial derivatives as ∂ α := ∂ x α 1 ⁡ ∂ y α 2 .

For ν ∈ [ 1 , ∞ ] , L ν ⁢ ( Ω ) denotes the Lebesgue space and W m , ν ⁢ ( Ω ) is a Sobolev space of order m. The space L 0 ν ⁢ ( Ω ) is defined as { q ∈ L ν ⁢ ( Ω ) : ∫ Ω q = 0 } while W 0 m , ν ⁢ ( Ω ) is the standard Sobolev space of order m with vanishing traces on ∂ ⁡ Ω . The Lebesgue measure of a domain ω ⊂ Ω is written | ω | . The L ν ⁢ ( ω ) -norm is ∥ ⋅ ∥ ν ; ω and the notation for the W m , ν ⁢ ( ω ) -norm is ∥ ⋅ ∥ m , ν ; ω . The integral ∫ ω u ⁢ v ⁢ 𝑑 ω is denoted ( u , v ) ω . The norm in the vector space 𝐖 m , ν ⁢ ( Ω ) := [ W m , ν ⁢ ( Ω ) ] d is ∥ 𝐪 ∥ m , ν , Ω = ( ∑ 1 ≤ i ≤ d ∑ 0 ≤ | α | ≤ m ∥ ∂ α ⁡ q i ∥ ν , Ω ν ) 1 ν . If there exists C > 0 so that f ≤ C ⁢ g , we will write f ≲ g , and we use the notation f ∼ g if there exist constants c , C > 0 so that c ⁢ f ≤ g ≤ C ⁢ f . For brevity of presentation we sometimes write integrals of some function g as ∫ g instead of ∫ g ⁢ 𝑑 x .

2.2 Weak Formulation and Discretization

The weak formulation of the p-Stokes system is: Find velocity and pressure ( 𝐯 , π ) ∈ 𝓧 × 𝒬 so that

(2.1) ( 𝐒 ⁢ ( 𝐃𝐯 ) , 𝐃 ⁢ 𝝎 ) Ω - ( π , ∇ ⋅ 𝝎 ) Ω = ( 𝐟 , 𝝎 ) Ω for all ⁢ 𝝎 ∈ 𝓧 ,

(2.2) ( ∇ ⋅ 𝐯 , q ) Ω = 0 for all ⁢ q ∈ 𝒬 .

The natural spaces for the solution are 𝓧 := 𝐖 0 1 , p ⁢ ( Ω ) and 𝒬 := L 0 p ′ ⁢ ( Ω ) , where p ′ = p p - 1 . We only state the case of homogeneous Dirichlet condition, because the case of inhomogeneous boundary data can be reformulated in the form (2.1)–(2.2) with modified right-hand side 𝐟 .

Now let 𝕋 h be a decomposition of Ω into rectangles. For simplicity and brevity of presentation, we assume the mesh to be globally uniform. On the decomposition 𝕋 h we introduce a finite element space of bi-linear, continuous functions

X h := X h ⁢ ( 𝕋 h ) := { ω ∈ C ⁢ ( Ω ¯ ) : ω | K ∘ F K ∈ ℚ 1 ⁢ ( K ^ ) ⁢ for all ⁢ K ∈ 𝕋 h } ,

where the space ℚ 1 ⁢ ( K ^ ) := span ⁡ { 1 ^ , x ^ , y ^ , x ^ ⁢ y ^ : ( x ^ , y ^ ) ∈ K ^ } is defined on the reference quadrilateral K ^ := ( 0 , 1 ) d , and F K maps from the reference element to an element K

F K ⁢ ( 𝐱 ^ ) = ( x 0 y 0 ) + ( h x 0 0 h y ) ⁢ ( x ^ y ^ ) ,

where the width of a quadrilateral in the x-direction is h x , and in the y-direction h y . We introduce the equal-order finite element spaces

𝓧 h := 𝓧 ∩ X h d and 𝒬 h := 𝒬 ∩ X h

and formulate a (unstabilized) discretized problem as: Find velocity and pressure ( 𝐯 h , π h ) ∈ 𝓧 h × 𝒬 h , so that

(2.3) ( 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ 𝝎 h ) Ω - ( π h , ∇ ⋅ 𝝎 h ) Ω = ( 𝐟 , 𝝎 h ) Ω for all ⁢ 𝝎 h ∈ 𝓧 h ,

(2.4) ( ∇ ⋅ 𝐯 h , q h ) Ω = 0 for all ⁢ q h ∈ 𝒬 h .

The equal-order finite element pair 𝓧 h × 𝒬 h violates the inf-sup condition [12], and to avoid spurious oscillations in the pressure, stabilization is needed. The most common argument for using stabilized equal-order elements instead of inf-sup stable elements is that they are easier to implement. Additionally, an inf-sup stable implementation on highly anisotropic meshes requires Q k × Q k - 2 -elements in general (rather than Q k × Q k - 1 -elements). If a continuous pressure is required, third-order polynomials should thus be used to approximate the velocity space. Such a high polynomial order in the velocity is usually sup-optimal because the low-order pressure approximation reduces the overall accuracy, and because of the low regularity of the solution [19]. Generalization to d = 3 with horizontal mesh sizes h x ∼ h y and smaller vertical mesh size h z ≪ h x is straightforward. However, we restrict here to the two-dimensional case for ease of presentation.

$Figure 1 A mesh 𝕋 h {\mathbb{T}_{h}} consisting of elements K. The element is anisotropic since it is longer in the x-direction than in the y-direction. Each element belongs to a patch M.$

Figure 1

A mesh 𝕋 h consisting of elements K. The element is anisotropic since it is longer in the x-direction than in the y-direction. Each element belongs to a patch M.

2.3 Stabilization

In this work we stabilize (2.3) with a non-linear stabilization term s ⁢ ( π h , q h ) of local projection type

(2.5) ( 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ 𝝎 h ) Ω - ( π h , ∇ ⋅ 𝝎 h ) Ω = ( 𝐟 , 𝝎 h ) Ω for all ⁢ 𝝎 h ∈ 𝓧 h ,

(2.6) ( ∇ ⋅ 𝐯 h , q h ) Ω + s ⁢ ( π h , q h ) = 0 for all ⁢ q h ∈ 𝒬 h .

To define the stabilization term we introduce some further notations. The mesh is arranged in a patch-wise manner, see Figure 1. Hence, we have a second shape-regular decomposition 𝕄 h of Ω where each element (or patch) M ∈ 𝕄 h consists of four rectangles of 𝕋 h so that 𝕄 h = 𝕋 2 ⁢ h . As in [5], we also define a space Y h consisting of piecewise constants on each patch M,

Y h := { ω ∈ L 2 ⁢ ( Ω ) : ω M ∘ F M ∈ ℚ 0 ⁢ ( M ^ ) } .

The restriction of the space Y h to a patch M is denoted by Y h ⁢ ( M ) := { q h | M : q h ∈ Y h } and we denote

X h 0 ⁢ ( M ) := { ω h | M : ω h ∈ X h , ω h = 0 ⁢ on ⁢ Ω ∖ M } .

We define a fluctuation operator as θ h := id - P h , where id is the identity and P h : L ν → Y h is projecting to a smaller space Y h so that ( P h ⁢ q h ) | M := P M ⁢ ( q h | M ) for all M ∈ 𝕄 h , where P M is the local projection P M : L ν ⁢ ( M ) → Y h ⁢ ( M ) . The projection in this setting is simply

P h ⁢ q h = 1 | M | ⁢ ∫ M q h .

The stabilization term is given by

(2.7) s ⁢ ( π h , q h ) := α 0 ⁢ ∑ M ∈ 𝕄 ( 𝐌 s ⁢ 𝜽 h ⁢ ∇ ⁡ π h , 𝜽 h ⁢ ∇ ⁡ q h ) M ,

(2.8) 𝐌 s ⁢ ( ∂ x ⁡ π h , ∂ y ⁡ π h ) := [ h x 2 ⁢ ( 1 τ ⁢ ( τ + | θ h ⁢ ∂ x ⁡ π h | ) ) p ′ - 2 0 0 h y 2 ⁢ ( 1 τ ⁢ ( τ + h y h x ⁢ | θ h ⁢ ∂ y ⁡ π h | ) ) p ′ - 2 ] ,

where τ ≥ 0 and α 0 ≥ 0 are user-defined parameters and 𝜽 denotes the fluctuation operator applied to vector-valued functions in a component-wise manner.

2.4 Anisotropic Meshes

The important point in this work is that the mesh is allowed to be anisotropic, so that

𝒜 := h x h y ≫ 1 .

That is why the factor h y h x and the factor h y 2 is introduced in (2.8). Without the former one, we will see that the discretization errors are proportional to ( h x h y ) p ′ - 2 2 . Such a factor can become extremely large as p → 1 (i.e. p ′ → ∞ and h x h y ≫ 1 . Exchanging the latter h y 2 factor to the more standard h x 2 factor increases errors even more. Note that for h x = h y = h the stabilization term is identical to that of [19], and for p = p ′ = 2 it is identical to the stabilization of [5] (if h x = h y = h ) or [10] (if h x > h y ).

On these anisotropic meshes, it is important that the following local inf-sup condition holds, with a constant that is independent of the aspect ratio 𝒜 .

Lemma 1.

For ν ≥ 1 there exists β loc > 0 independent of h x and h y such that

inf q h ∈ Y h ⁢ ( M ) ⁡ sup ω h ∈ X h 0 ⁢ ( M ) ⁡ ( ω h , q h ) M ∥ ω h ∥ ν , M ⁢ ∥ q h ∥ ν ′ , M ≥ β loc > 0

for all h x > 0 , h y > 0 and all M ∈ M , where ν ′ = ν ν - 1 .

Here it is important that β loc is independent of the aspect ratio 𝒜 . The parameter ν will in practice equal to p.

Proof.

The proof is similar to that of [21, Lemma 3.2], which considers isotropic elements and ν = ν ′ = 2 . The extension from ν = 2 to ν ∈ ( 1 , ∞ ) is straightforward, see [18]. In order to extend the proof to anisotropic meshes, we consider that the proof in [21] is based on transferring back and forth between reference elements and physical elements, using an estimation of the determinant of the Jacobian of the transformation from the reference elements. On our simple anisotropic meshes we can use that we clearly have | det D F K | = h x h y . Following the approach of [21], we have the desired estimate. Details are given in the Appendix. ∎

The last property that is required for the LPS-stabilization is for the fluctuation operator θ h to be bounded, so that

∥ θ h ⁢ q ∥ 0 , ν ; M ≲ ∥ q ∥ 0 , ν ; M for all M ∈ 𝕄 and all ⁢ q ∈ L ν ⁢ ( Ω ) .

This property follows directly from [18], also on anisotropic meshes.

2.5 Some Useful Help-Functions and Properties of 𝐒

2.5.1 The Help-Functions φ, 𝓕 and Their Relation to the Stress Tensor 𝐒

Since the stress-tensor is non-linear, some extra functions and properties are required for the error analysis. The extra help functions that will be used in the analysis are the N-function φ , φ * ∈ C 1 ⁢ ( ℝ 0 + , ℝ 0 + ) , the function 𝓕 : ℝ d × d → ℝ sym d × d , and 𝓖 1 , 𝓖 2 : ℝ → ℝ . The functions φ, φ * and 𝓕 are related to the continuous formulation, and their properties and their relations to 𝐒 are therefore identical to those in [6, 18, 19]. The N-function φ and the complementary N-function φ * are defined as

φ ⁢ ( t ) = 1 p ⁢ t p and φ * ⁢ ( t ) = 1 p ′ ⁢ t p ′ .

The shifted N-functions φ a ∈ C 1 ⁢ ( ℝ 0 + , ℝ 0 + ) are defined as

φ a ⁢ ( t ) = ∫ 0 t φ a ′ ⁢ ( s ) ⁢ 𝑑 s , φ a ′ ⁢ ( t ) := φ ′ ⁢ ( a + t ) ⁢ t a + t ,

and similarly for φ * . In particular, it is useful to apply a shift of ϵ to φ to get the φ ϵ -function

φ ϵ ⁢ ( t ) = ∫ 0 t ( ϵ + s ) p - 2 ⁢ s ⁢ 𝑑 s ,

which is of interest since it has a structure that is related to the stress-tensor. The function 𝓕 also has a structure reminiscent of the stress-tensor. It is defined as

(2.9) 𝓕 ⁢ ( 𝐏 ) = ( ϵ + | 𝐏 sym | ) p - 2 2 ⁢ 𝐏 sym .

In the analysis, we will exploit some relations between 𝐒 , 𝐃 , φ and 𝓕 as given in e.g. [19], summarized in Lemmas 13 and 14 in the Appendix. These relations allow for working with the help-functions in the error-analysis instead of with the stress-tensor directly. An example of how this can be advantageous is that since φ is convex one can use Jensen’s inequality, in cases where one in a linear analysis would have used the triangle inequality. Also, φ satisfies the Δ 2 -condition, i.e. there exists a c > 0 such that for all t ≥ 0 , φ ⁢ ( 2 ⁢ t ) ≤ c ⁢ φ ⁢ ( t ) . The smallest constant is called Δ 2 . We have that Δ 2 < ∞ , implying that φ ⁢ ( t ) ∼ φ ⁢ ( c ⁢ t ) (c being any constant), which clearly is a desired property in analysis. Even though we will consider a Carreau-type fluid with a stress tensor 𝐒 of the form (1.2) in our numerical experiments, Lemmas 13 and 14 are valid for any stress-tensor that fulfills the following assumption, with p ≤ 2 .

Assumption 1 (On the Stress Tensor 𝐒 ).

The extra stress tensor

𝐒 : ℝ d × d → ℝ sym d × d

belongs to C 0 ⁢ ( ℝ d × d , ℝ sym d × d ) ∩ C 1 ⁢ ( ℝ d × d ∖ { 0 } , ℝ sym d × d ) and satisfies S ⁢ ( 𝐐 ) = S ⁢ ( 𝐐 sym ) and 𝐒 ⁢ ( 𝟎 ) = 𝟎 . The stress tensor is of ( p , ϵ ) -structure, meaning that there exist p ∈ ( 1 , ∞ ) , ϵ ∈ [ 0 , ∞ ) , and constants σ 0 , σ 1 > 0 such that for all 𝐐 , 𝐏 ∈ ℝ d × d with 𝐐 ≠ 0 ,

∑ i , j , k , l = 1 d ∂ k ⁢ l ⁡ S i ⁢ j ⁢ ( 𝐐 ) ⁢ P i ⁢ j ⁢ P k ⁢ l ≥ σ 0 ⁢ ( ϵ + | 𝐐 sym | ) p - 2 ⁢ | 𝐏 sym | 2 ,

∂ k ⁢ l ⁡ S i ⁢ j ⁢ ( 𝐐 ) ≤ σ 1 ⁢ ( ϵ + | 𝐐 sym | ) p - 2 for all i , j , k , l ∈ { 1 , … , d } .

The constants in our error estimates will often depend on σ 0 and σ 1 . For brevity of presentation we will, except for in this subsection, not state this in the remainder of this paper.

2.5.2 The Help-Functions 𝓖 1 , 𝓖 2 and Their Relation to the Stabilization Term and φ *

To analyze the effect of the stabilization, some relations between the non-linear stabilization term s h ⁢ ( p , q ) , the help functions 𝓖 1 and 𝓖 2 and the complementary N-function φ * ⁢ ( t ) = t p ′ p ′ are required. We here adapt the analysis of [18, 19] for the isotropic problem to anisotropic meshes. In particular, we use two functions 𝓖 1 and 𝓖 2

𝓖 1 ⁢ ( 𝐪 ) = ( τ + | 𝐪 | ) p ′ - 2 2 ⁢ 𝐪 , 𝓖 2 ⁢ ( 𝐪 ) = ( τ 2 + | 𝐪 | ) p ′ - 2 2 ⁢ 𝐪 ,

instead of one ( 𝓖 in [18, 19]), and we revisit [19, Lemmas 3.1–3.3] to state them in an anisotropic setting. Here τ 2 = h x h y ⁢ τ . The function 𝓖 1 will be used to analyze the x-components of the error, and 𝓖 2 will be used for the y-components.

Lemma 2.

Let U ⊆ Ω be a measurable subset and let g , q ∈ L p ′ ⁢ ( Ω ) , with p ∈ ( 1 , 2 ] . Then

∥ 𝐠 - 𝐪 ∥ p ′ , U p ′ ≲ ∥ 𝓖 i ⁢ ( 𝐠 ) - 𝓖 i ⁢ ( 𝐪 ) ∥ 2 , U 2 , i = 1 , 2 ,

where the constant c only depends on p.

Proof.

The proof is based on a vector-valued version of Lemma 13, see [18] for details. ∎

Lemma 3.

For all π , q ∈ W 1 , p ′ ⁢ ( Ω ) there holds

s ⁢ ( π , π - q ) - s ⁢ ( q , π - q )

∼ α 0 ⁢ ∑ M ∈ 𝕄 h ( h x 2 ⁢ ∥ 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ π ) - 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ q ) ∥ 2 ; M 2 + ( h y h x ) p ′ - 2 ⁢ h y 2 ⁢ ∥ 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ π ) - 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ q ) ∥ 2 ; M 2 )

where φ is the complementary N-function. The constants depend on p and τ.

Proof.

We have

s ⁢ ( π , π - q ) - s ⁢ ( q , π - q ) = α 0 ⁢ ∑ M ∈ 𝕄 h ( M s , 11 ⁢ ( π ) ⁢ θ h ⁢ ∂ x ⁡ π - M s , 11 ⁢ ( q ) ⁢ θ h ⁢ ∂ x ⁡ q , θ h ⁢ ∂ x ⁡ ( π - q ) )

+ α 0 ⁢ ∑ M ∈ 𝕄 h ( M s , 22 ⁢ ( π ) ⁢ θ h ⁢ ∂ y ⁡ π - M s , 22 ⁢ ( q ) ⁢ θ h ⁢ ∂ y ⁡ q , θ h ⁢ ∂ y ⁡ ( π - q ) ) .

Rearranging terms, applying a scalar-valued version of Lemma 13 with p, ϵ, 𝓕 and φ replaced by p ′ , τ, 𝓖 1 , and φ * and remembering the definition of 𝓖 1 yields for the x-component, for each M ∈ 𝕄 h ,

( M s , 11 ⁢ ( π ) ⁢ θ h ⁢ ∂ x ⁡ π - M s , 11 ⁢ ( q ) ⁢ θ h ⁢ ∂ x ⁡ q , θ h ⁢ ∂ x ⁡ ( π - q ) ) M

= ( 1 τ ) p ′ - 2 ⁢ h x 2 ⁢ ( ( τ + | θ h ⁢ ∂ x ⁡ π | ) p ′ - 2 ⁢ θ h ⁢ ∂ x ⁡ π - ( τ + | θ h ⁢ ∂ x ⁡ q | ) p ′ - 2 ⁢ θ h ⁢ ∂ x ⁡ q , θ h ⁢ ∂ x ⁡ ( π - q ) ) M

∼ h x 2 ⁢ ∥ 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ π ) - 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ q ) ∥ 2 ; M 2

∼ h x 2 ⁢ ∫ M φ τ + | θ h ⁢ ∂ x ⁡ π | * ⁢ ( | θ h ⁢ ∂ x ⁡ π - θ h ⁢ ∂ x ⁡ q | ) ⁢ 𝑑 x .

We use the same approach to bound the y-component, but first we rewrite the factor M s , 22 ⁢ ( π ) so that

( τ + h y h x ⁢ | θ h ⁢ ∂ y ⁡ π | ) p ′ - 2 = ( h y h x ) p ′ - 2 ⁢ ( τ 2 + | θ h ⁢ ∂ y ⁡ π | ) p ′ - 2 .

The y-component is then

( M s , 22 ⁢ ( π ) ⁢ θ h ⁢ ∂ y ⁡ π - M s , 22 ⁢ ( q ) ⁢ θ h ⁢ ∂ y ⁡ q , θ h ⁢ ∂ y ⁡ ( π - q ) ) M

= ( 1 τ ) p ′ - 2 ⁢ h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ( ( τ 2 + | θ h ⁢ ∂ y ⁡ π | ) p ′ - 2 ⁢ θ h ⁢ ∂ y ⁡ π - ( τ 2 + | θ h ⁢ ∂ y ⁡ q | ) p ′ - 2 ⁢ θ h ⁢ ∂ y ⁡ q , θ h ⁢ ∂ y ⁡ ( π - q ) ) M

∼ h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ π ) - 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ q ) ∥ 2 ; M 2

∼ h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ∫ M φ τ 2 + | θ h ⁢ ∂ y ⁡ π | * ⁢ ( | θ h ⁢ ∂ y ⁡ π - θ h ⁢ ∂ y ⁡ q | ) ⁢ 𝑑 x . ∎

Lemma 4.

For all δ > 0 there exists a constant c = c ⁢ ( δ , p , α 0 , τ ) such that for all π , q ∈ W 1 , p ′ ⁢ ( Ω ) ,

| s ⁢ ( π , π - q ) | ≲ α 0 ⁢ c δ ⁢ h x 2 ⁢ ∥ τ + | ∂ x ⁡ π | ∥ p ′ ; Ω p ′ + α 0 ⁢ c δ ⁢ h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ τ 2 + | ∂ y ⁡ π | ∥ p ′ ; Ω p ′

+ α 0 δ ∑ M ∈ 𝕄 h ( h x 2 ∥ 𝓖 1 ( θ h ∂ x π ) - 𝓖 1 ( θ h ∂ x q ) ∥ 2 ; M 2

+ h y 2 ( h y h x ) p ′ - 2 ∥ 𝓖 2 ( θ h ∂ y π ) - 𝓖 2 ( θ h ∂ y q ) ∥ 2 ; M 2 ) .

Proof.

Again rewriting

( τ + h y h x ⁢ | θ h ⁢ ∂ y ⁡ π | ) p ′ - 2 = ( h y h x ) p ′ - 2 ⁢ ( τ 2 + | θ h ⁢ ∂ y ⁡ π | ) p ′ - 2 ,

working component-wise and applying a scalar-valued version of Lemma 13 with p,φ, and ϵ replaced by p ′ , φ * , and τ, we have

+ ( h y h x ) p ′ - 2 ⁢ h y 2 ⁢ ∫ M ( φ * ) τ 2 + θ h ⁢ ∂ y ⁡ π ′ ⁢ ( | θ h ⁢ ∂ y ⁡ π | ) ⁢ | θ h ⁢ ∂ y ⁡ π - θ h ⁢ ∂ y ⁡ q | ⁢ 𝑑 x .

Now we apply Young’s inequality for N-functions from Lemma 14. For the x-component this gives

∑ M ∈ 𝕄 h x 2 ⁢ ∫ M ( φ * ) τ + θ h ⁢ ∂ x ⁡ π ′ ⁢ ( | θ h ⁢ ∂ x ⁡ π | ) ⁢ | θ h ⁢ ∂ x ⁡ π - θ h ⁢ ∂ x ⁡ q | ⁢ 𝑑 x

≲ c δ ⁢ h x 2 ⁢ ∫ Ω φ τ + θ h ⁢ ∂ x ⁡ π * ⁢ ( | θ h ⁢ ∂ x ⁡ π | ) ⁢ 𝑑 x + δ ⁢ ∑ M ∈ 𝕄 h x 2 ⁢ ∥ 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ π ) - 𝓖 1 ⁢ ( θ h ⁢ ∂ x ⁡ q ) ∥ 2 ; M 2 .

The integral on the left-hand side can be bounded with Lemma 3, and Lemma 13 (change of shift), and using that φ * ⁢ ( t ) := t p ′ p ′ , and lastly the stability of the fluctuation operator,

∫ Ω φ τ + θ h ⁢ ∂ x ⁡ π * ⁢ ( | θ h ⁢ ∂ x ⁡ π | ) ⁢ 𝑑 x ≲ ∫ Ω φ * ⁢ ( τ + | θ h ⁢ ∂ x ⁡ π | ) ⁢ 𝑑 x ≲ ∥ τ + | ∂ x ⁡ π | ∥ p ′ ; Ω p ′ .

Similarly, for the y-component,

∑ M ∈ 𝕄 ( h y h x ) p ′ - 2 ⁢ h y 2 ⁢ ∫ M ( φ * ) τ 2 + θ h ⁢ ∂ y ⁡ π ′ ⁢ ( | θ h ⁢ ∂ y ⁡ π | ) ⁢ | θ h ⁢ ∂ y ⁡ π - θ h ⁢ ∂ y ⁡ q | ⁢ 𝑑 x

≲ c δ ⁢ ( h y h x ) p ′ - 2 ⁢ h y 2 ⁢ ∫ Ω ( φ * ) τ 2 + θ h ⁢ ∂ y ⁡ π ⁢ ( | θ h ⁢ ∂ y ⁡ π | ) ⁢ 𝑑 x + δ ⁢ ∑ M ∈ 𝕄 ( h y h x ) p ′ - 2 ⁢ h y 2 ⁢ ∥ 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ π ) - 𝓖 2 ⁢ ( θ h ⁢ ∂ y ⁡ q ) ∥ 2 ; M 2

and

∫ Ω ( φ * ) τ 2 + θ h ⁢ ∂ y ⁡ π ⁢ ( | θ h ⁢ ∂ y ⁡ π | ) ⁢ 𝑑 x ≲ ∫ Ω ( φ * ) ⁢ ( τ 2 + | θ h ⁢ ∂ y ⁡ π | ) ⁢ 𝑑 x ≲ ∥ τ 2 + | ∂ y ⁡ π | ∥ p ′ ; Ω p ′ . ∎

3 Error Analysis

3.1 Modified Interpolation Operator

3.1.1 Standard Properties in an Anisotropic p-Stokes Setting

The idea behind the error analysis of LPS-methods is to construct an interpolation operator j h : W 1 , p ⁢ ( Ω ) → X h (and 𝐣 h : 𝐖 0 1 , ν ⁢ ( Ω ) → 𝓧 h ) so that the interpolation error is orthogonal to Y h (or 𝐘 h := [ Y h ] d ). This interpolation operator should also fulfill some useful approximation and stability properties. The following lemma states these standard LPS properties in an anisotropic p-Stokes setting. The parameter ν will in practice be equal either p or p ′ . For brevity of presentation, we only present the properties for the scalar j h but it also extends to 𝐣 h and the vector-valued case.

Lemma 5.

Let ν ≥ 1 and let the two spaces X h and Y h satisfy Lemma 1. Then there exist an interpolation operator j h : W 1 , ν ⁢ ( Ω ) → X h satisfying the following properties:

Orthogonality with respect to Y h : For all ω ∈ W 1 , ν ⁢ ( Ω ) ,

( ω - j h ⁢ ω , q h ) Ω = 0 for all ⁢ q h ∈ Y h .
Approximation properties: For all M ∈ 𝕄 h and all ω ∈ W 1 , ν ⁢ ( Ω ) ,

(3.1) ∥ ω - j h ⁢ ω ∥ 0 , ν ; M ≲ h x ⁢ ∥ ∂ x ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) ,

(3.2) ∥ ∇ ⁡ ( ω - j h ⁢ ω ) ∥ ν , M ≲ ∥ ω ∥ 1 , ν ; Λ ⁢ ( M ) + h x h y ⁢ ∥ ∂ x ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + ∥ ∂ y ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) ,

where Λ ⁢ ( M ) = int ⁢ ⋃ { M ′ ¯ : M ′ ∈ 𝕄 h , M ¯ ∩ M ′ ¯ ≠ ∅ } is a neighborhood of the patch M . For all ω ∈ W 2 , ν ⁢ ( Ω ) ,

(3.3) ∥ ω - j h ⁢ ω ∥ 0 , ν ; M ≲ h x 2 ⁢ ∥ ∂ x 2 ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h y 2 ⁢ ∥ ∂ y 2 ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h x ⁢ h y ⁢ ∥ ∂ x ⁢ y ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) ,

(3.4) ∥ ∇ ⁡ ( ω - j h ⁢ ω ) ∥ ν , M ≲ h x ⁢ ∥ ∂ x ⁡ ω ∥ 1 , ν ; Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ ω ∥ 1 , ν ; Λ ⁢ ( M ) + h x 2 h y ⁢ ∥ ∂ x 2 ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) .

The constants depend on the inf-sup constant β loc of Lemma 1.
Stability properties: For all ω ∈ W 1 , ν ⁢ ( Ω ) ,

(3.5) ∥ ∇ ⁡ j h ⁢ ω ∥ ν , M ≲ ( 1 + h x h y ) ⁢ ∥ ∇ ⁡ ω ∥ ν , M ,

(3.6) ∥ ∂ x ⁡ j h ⁢ ω ∥ ν , M ≲ ∥ ∂ x ⁡ ω ∥ ν ; Λ ⁢ ( M ) + h x - 1 ⁢ h y ⁢ ∥ ∂ y ⁡ ω ∥ ν , Λ ⁢ ( M ) ,

(3.7) ∥ ∂ y ⁡ j h ⁢ ω ∥ ν , M ≲ ∥ ∂ y ⁡ ω ∥ ν ; Λ ⁢ ( M ) + h y - 1 ⁢ h x ⁢ ∥ ∂ x ⁡ ω ∥ ν , Λ ⁢ ( M ) .

Proof.

The interpolation operator j h (as well as 𝐣 h ) is constructed from two parts, so that j h ⁢ ω | M = i h ⁢ ω + z h ⁢ ( ω ) . The base interpolator i h ensures approximation properties, and z h ⁢ ( ω ) is added to ensure the orthogonality property. Hence, Λ ⁢ ( M ) is a patch of patches similarly defined as for the anisotropic Scott–Zhang interpolator in [2].

On isotropic meshes, the Scott–Zhang interpolation is a possible choice for i h , see e.g. [18] and [21]. On anisotropic elements, we instead use the anisotropic Scott–Zhang-type interpolators described in [2, 4, 8] which all coincide for the non-sheared elements that we study here. The anisotropic Scott–Zhang interpolation operator i h is constructed by integrating only along the long edges in the x-direction, instead of along any edge as is the case for the original Scott–Zhang interpolation operator. Various approximation and stability properties of i h have been shown in [2, 4, 8] and we choose the, for our purpose, most useful ones (see Appendix).

For z h we know that due to Lemma 1 there exists a unique z h ⁢ ( ω ) ∈ Z h ⁢ ( M ) ⊥ which ensures that

( ω - i h ⁢ ω - z h ⁢ ( ω ) , q h ) M = 0

(where Z h ⁢ ( M ) := { ω h ∈ X h 0 ⁢ ( M ) : ( ω h , q h ) M = 0 ⁢ for all ⁢ q h ∈ Y h ⁢ ( M ) } ), see [16, 18, 21]. This z h ⁢ ( ω ) also fulfills

(3.8) ∥ z h ⁢ ( ω ) ∥ ν ; M ≤ 1 β loc ⁢ ∥ ω - i h ⁢ ω ∥ ν ; M .

Here β loc is the inf-sup constant from Lemma 1. Using the triangle inequality, the approximation property of z h (3.8) and the approximation properties of i h , we arrive at the approximation properties of j h in the L ν -norm, i.e. (3.1) and (3.3). To bound the gradient of the interpolation error we need to bound the gradient of z h . We do this by combining (3.8) with the inverse inequality from [9], which yields

(3.9) ∥ ∂ x ⁡ z h ⁢ ( ω ) ∥ ν , M ≲ 1 h x ⁢ ∥ z h ⁢ ( ω ) ∥ ν , M ≲ 1 β loc ⁢ h x ⁢ ∥ ω - i h ⁢ ω ∥ ν , M ,

(3.10) ∥ ∂ y ⁡ z h ⁢ ( ω ) ∥ ν , M ≲ 1 h y ⁢ ∥ z h ⁢ ( ω ) ∥ ν , M ≲ 1 β loc ⁢ h y ⁢ ∥ ω - i h ⁢ ω ∥ ν , M .

Again combining with the triangle inequality

∥ ∇ ⁡ ( ω - j h ⁢ ω ) ∥ ν , M ≤ ∥ ∇ ⁡ ( ω - i h ⁢ ω ) ∥ ν , M + ∥ ∇ ⁡ z h ⁢ ( ω ) ∥ ν , M

and the approximation properties of i h , we have the desired estimates (3.2) and (3.4). Using that

∥ j h ⁢ ω ∥ ν , M = ∥ j h ⁢ ω + ω - ω ∥ ν , M ≤ ∥ ω ∥ ν , M + ∥ ω - j h ⁢ ω ∥ ν , M

together with the approximation property (3.1), we have also L ν -stability.

For W 1 , ν -stability we use the stability of the base interpolator i h of equation (C.2) and (C.3), or alternatively (C.4) and (C.5) together with the inverse estimates (3.9) and (3.10). ∎

Given j h , it is easy to construct a corresponding vector valued interpolation operator 𝐣 h : 𝐖 0 1 , ν ⁢ ( Ω ) → 𝓧 h which renders the interpolation error orthogonal to 𝐘 h .

3.1.2 Approximation Property in Orlicz Spaces on Anisotropic Meshes

The approximation and stability properties of j h described above are not fundamentally different from those needed in the linear Stokes case. However, for the p-Stokes case, we need an additional approximation property in Orlicz spaces, more specifically in terms of the help function 𝓕 . On isotropic meshes such an interpolation operator was shown to exist in [18]. Here we construct operators j h , 𝐣 h suitable for anisotropic meshes.

Lemma 6.

Let ν ≥ 1 , let X h and Y h satisfy Lemma 1 and let F be defined by (2.9). If F ⁢ ( D ⁢ 𝛚 ) ∈ [ W 1 , 2 ⁢ ( Ω ) ] d × d , then there exists an interpolation operator j : W 0 1 , ν ⁢ ( Ω ) → X h such that for all M ∈ M ,

(3.11) ∥ 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝝎 h ) ∥ 2 ; M ≲ h x h y ⁢ ( h x ⁢ ∥ ∂ x ⁡ 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) ∥ 2 , Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) ∥ 2 , Λ ⁢ ( M ) ) .

The remainder of this subsection is dedicated to prove Lemma 6. For this purpose we transfer from working with 𝓕 , to working in Orlicz–Sobolev spaces with N-functions. We denote classical Orlicz spaces L φ ⁢ ( Ω ) and Orlicz–Sobolev spaces with W 1 , φ ⁢ ( Ω ) . We have that f ∈ L φ ⁢ ( Ω ) if and only if ∫ Ω φ ⁢ ( | f | ) ⁢ 𝑑 x < ∞ and f ∈ W 1 , φ ⁢ ( Ω ) if and only if ∇ i ⁡ f ∈ L φ ⁢ ( Ω ) for 0 ≤ i ≤ 1 . We then need interpolation estimates in terms of N-functions, which requires Orlicz-stability, Orlicz-approximability and Orlicz-continuity. For the isotropic case, these properties are shown in [14], and we here extend them to anisotropic meshes (for low-order approximations).

Orlicz-Stability.

We start by proving Orlicz-stability. We first need to show that j h fulfills an auxiliary stability property corresponding to [14, Assumption 4.1].

Lemma 7.

The interpolation operator j h satisfies the following W 1 , 1 -stability:

⨍ M | 𝐣 h ⁢ 𝝎 | + ⨍ M h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ≲ ⨍ Λ ⁢ ( M ) | 𝝎 | + h x h y ⁢ h x ⁢ ⨍ Λ ⁢ ( M ) | ∇ ⁡ 𝝎 | ,

where ⨍ U ω := 1 | U | ⁢ ∫ U ω denotes the mean value over U.

Proof.

The interpolation estimate (3.1) with ν = 1 bounds the first term

(3.12) ⨍ M | 𝐣 h ⁢ 𝝎 | ≤ ⨍ M | 𝐣 h ⁢ 𝝎 - 𝝎 | + ⨍ M | 𝝎 | ≲ ( 1 + 1 β loc ) ⁢ ⨍ Λ ⁢ ( M ) ∑ | α | = 1 h α ⁢ | D α ⁢ ω | + ⨍ M | 𝝎 | ≲ ∑ | α | ≤ 1 ⨍ Λ ⁢ ( M ) h α ⁢ | D α ⁢ ω | .

Here, α is a multi-index. For the gradient terms we use the directional stability estimates (3.6) and (3.7) to get

⨍ M | ∂ x ⁡ 𝐣 h ⁢ 𝝎 ⁢ | + ⨍ M | ⁢ ∂ y ⁡ 𝐣 h ⁢ 𝝎 | ≲ ⨍ Λ ⁢ ( M ) | ∂ x ⁡ 𝝎 ⁢ | + h y h x ⁢ ⨍ Λ ⁢ ( M ) | ⁢ ∂ y ⁡ 𝝎 ⁢ | + h x h y ⁢ ⨍ Λ ⁢ ( M ) | ⁢ ∂ x ⁡ 𝝎 ⁢ | + ⨍ Λ ⁢ ( M ) | ⁢ ∂ y ⁡ 𝝎 |

(3.13) = ( 1 + h x h y ) ⁢ ⨍ Λ ⁢ ( M ) | ∂ x ⁡ 𝝎 ⁢ | + ( 1 + h y h x ) ⁢ ⨍ Λ ⁢ ( M ) | ⁢ ∂ y ⁡ 𝝎 | .

Since M is as anisotropic as Λ ⁢ ( M ) , combining (3.12) and (3.13) we have

⨍ M | 𝐣 h ⁢ 𝝎 | + ⨍ M h x ⁢ | ∂ x ⁡ 𝐣 h ⁢ 𝝎 ⁢ | + ⨍ M h y | ⁢ ∂ y ⁡ 𝐣 h ⁢ 𝝎 | ≲ ⨍ Λ ⁢ ( M ) | 𝝎 | + ( h x + h x 2 h y ) ⁢ ⨍ Λ ⁢ ( M ) | ∂ x ⁡ 𝝎 ⁢ | + ( h y + h y 2 h x ) ⁢ ⨍ Λ ⁢ ( M ) | ⁢ ∂ y ⁡ 𝝎 | .

Rearranging and using that h x h y ≥ h x ≥ h y completes the proof. ∎

Lemma 8 (First-Order Anisotropic Orlicz-Stability).

For a j h that satisfies Lemma 7 it holds

⨍ M φ ⁢ ( | 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ) ≲ ⨍ Λ ⁢ ( M ) φ ⁢ ( | 𝝎 | ) + ⨍ Λ ⁢ ( M ) φ ⁢ ( h x 2 h y ⁢ | ∇ ⁡ 𝝎 | ) .

Proof.

Since it holds that

| 𝝎 h ⁢ ( x ) | ≲ ⨍ K | 𝝎 h ⁢ ( y ) | , | ∇ ⁡ 𝝎 h ⁢ ( x ) | ≲ ⨍ K | ∇ ⁡ 𝝎 h ⁢ ( y ) | ,

we have

⨍ M φ ⁢ ( | 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ) ≲ ⨍ M φ ⁢ ( ⨍ M | 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( ⨍ M h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ) .

Using Lemma 7, the convexity of φ to exchange the order of summation and integration, the Δ 2 -condition to move constants, and Jensen’s inequality for integrals, we get

⨍ M φ ⁢ ( | 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ) ≲ ⨍ M φ ⁢ ( ⨍ Λ ⁢ ( M ) | 𝝎 | + ⨍ Λ ⁢ ( M ) h x 2 h y ⁢ | ∇ ⁡ 𝝎 | ) ≲ c ⁢ ⨍ M φ ⁢ ( ⨍ Λ ⁢ ( M ) | 𝝎 | ) + ⨍ M φ ⁢ ( ⨍ Λ ⁢ ( M ) h x 2 h y ⁢ | ∇ ⁡ 𝝎 | )

≲ c ⁢ ⨍ Λ ⁢ ( M ) φ ⁢ ( | 𝝎 | ) + ⨍ Λ ⁢ ( M ) φ ⁢ ( h x 2 h y ⁢ | ∇ ⁡ 𝝎 | ) . ∎

Orlicz-Approximability.

To prove Orlicz-approximability, we need the Orlicz-stability proven above, but also first-order, anisotropic versions of [14, Lemma 3.1, Theorem 3.2 and Corollary 3.3], collected in the following lemma.

Lemma 9.

For all 𝛚 ∈ ( W 1 , φ ⁢ ( K ) ) d there exists Q ⁢ 𝛚 ∈ ( ( P 0 ) ⁢ ( K ) ) d such that

⨍ K φ ⁢ ( | 𝝎 - 𝐐 ⁢ 𝝎 | ) + ⨍ K φ ⁢ ( h x ⁢ | ∇ ⁡ ( 𝝎 - 𝐐 ⁢ 𝝎 ) | ) ≲ ⨍ K φ ⁢ ( h x ⁢ | ∇ ⁡ 𝝎 | ) .

It also holds on patches M and neighborhoods of elements and patches Λ ⁢ ( K ) , Λ ⁢ ( M ) .

Proof.

The bound of the second term is trivial since 𝐐𝐯 ∈ ( ( 𝒫 0 ) ⁢ ( K ) ) N so that ∇ ⁡ ( 𝐐𝐯 ) = 0 . In the isotropic setting in [14] the first term is bounded by setting 𝐐𝐯 as an averaged Taylor polynomial, but on anisotropic elements this will result in a constant which depends heavily on the aspect ratio, see [11]. As we do not have higher-order derivatives (as in [14]) we can instead simply use that 𝐐𝐯 is a constant, and apply an anisotropic version of the Poincaré-Wirtinger inequality (see e.g. [8, Lemma 2] ), which states that there exists a constant c so that

(3.14) ∥ 𝐯 - 𝐜 ∥ ν ≲ h x ⁢ ∥ ∂ x ⁡ 𝐯 ∥ ν + h y ⁢ ∥ ∂ y ⁡ 𝐯 ∥ ν ,

for ν ≥ 1 . The first term is thus bounded by

⨍ K φ ⁢ ( | 𝐯 - 𝐐𝐯 | ) ≲ ⨍ K φ ⁢ ( h x ⁢ | ∂ x ⁡ 𝐯 ⁢ | + h y | ⁢ ∂ y ⁡ 𝐯 | ) ≲ ⨍ K φ ⁢ ( h x ⁢ | ∇ ⁡ 𝐯 | ) . ∎

Given the above lemma and Orlicz-stability, Orlicz-approximability follows.

Lemma 10 (Anisotropic Orlicz-Approximability).

Proof.

Adding and subtracting a 𝔭 ∈ V h , using the convexity of φ, the Δ 2 -condition and that 𝐣 h is a projection, we have

LHS := ⨍ M φ ⁢ ( | 𝝎 - 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ ( 𝝎 - 𝐣 h ⁢ 𝝎 ) | )

= ⨍ M φ ⁢ ( | 𝝎 - 𝔭 + 𝔭 - 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ ( 𝝎 - 𝔭 + 𝔭 - 𝐣 h ⁢ 𝝎 ) | )

≲ ⨍ M φ ⁢ ( | 𝝎 - 𝔭 | ) + ⨍ M φ ⁢ ( | 𝔭 - 𝐣 h ⁢ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ ( 𝝎 - 𝔭 ) | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ ( 𝔭 - 𝐣 h ⁢ 𝝎 ) | ) .

Applying Lemma 8 (Orlicz-stability) on the second and fourth term with 𝝎 replaced with 𝝎 - 𝔭 , and then using Lemma 9 ( 𝔭 is arbitrary), we get

LHS ≲ ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝝎 | ) + ⨍ Λ ⁢ ( M ) φ ⁢ ( h x ⁢ | ∇ ⁡ 𝝎 | ) + ⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝝎 | ) + ⨍ Λ ⁢ ( M ) φ ⁢ ( h x 2 h y ⁢ | ∇ ⁡ 𝝎 | ) .

Rearranging, using that h x 2 h y > h x and the Δ 2 -condition, and shifting all integration areas from M to Λ ⁢ ( M ) , we get the desired result. ∎

Corollary 1 (Orlicz-Continuity).

For all 𝛚 ∈ ( W 1 , φ ⁢ ( M ) ) N ,

⨍ M φ ⁢ ( h x ⁢ | ∇ ⁡ 𝐣 h ⁢ 𝝎 | ) ≲ ⨍ Λ ⁢ ( M ) φ ⁢ ( h x 2 h y ⁢ | ∇ ⁡ 𝝎 | ) .

Proof.

By ignoring the low-order term on the left-hand side in Lemma 10 and using the triangle equality and convexity of φ, we have the result. ∎

Corollary 2.

For all v ∈ ( W 1 , φ ⁢ ( M ) ) N ,

⨍ M φ a ⁢ ( | ∇ ⁡ 𝐣 h ⁢ 𝐯 | ) ≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) φ a ⁢ ( | ∇ ⁡ 𝐯 | ) .

Proof.

The result can be derived from Corollary 1 by choosing 𝐯 = h x ⁢ 𝝎 and using the Δ 2 -condition. ∎

Orlicz-Continuity and Orlicz-Approximability in Terms of 𝓕

We derive the approximation property of j h in terms of 𝓕 , i.e. we prove Lemma 6. The proof mostly follows the (lengthy) isotropic proof of [18, Lemma 4.4]. Here we only give details on the parts that differs, i.e. when the anisotropic Orlicz-continuity and Poincaré’s inequality is used.

Proof of Lemma 6 .

We start by transferring from 𝓕 to working with φ (see [18, p. 66]) for details on this step,

⨍ M | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝝎 ) | 2 ⁢ 𝑑 𝐱 ≲ ⨍ M φ ϵ + | 𝐃𝐪 | ⁢ ( | 𝐃𝐣 h ⁢ ( 𝐪 - 𝝎 ) | ) ⏟ = ⁣ : T 1 + ⨍ M φ ϵ + | 𝐃𝐪 | ( | 𝐃 𝝎 - 𝐃𝐪 ) | ) ⏟ = ⁣ : T 2 ,

where 𝐪 is arbitrary linear polynomial. The first term is bounded using Corollary 2

T 1 ≲ ⨍ M φ ϵ + | 𝐃𝐪 | ⁢ ( | 𝐃𝐣 h ⁢ ( 𝐪 - 𝝎 ) | ) ≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | ∇ ⁡ ( 𝐪 - 𝝎 ) | ) .

Choosing ∇ ⁡ q := ⨍ Λ ⁢ ( M ) ∇ ⁡ 𝝎 , and using a version of Korn’s inequality for N-functions (see [15, Theorem 6.13]), we have

h x h y ⁢ ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | ∇ ⁡ ( 𝐪 - 𝝎 ) | ) = h x h y ⁢ ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | ∇ ⁡ 𝝎 - ⨍ Λ ⁢ ( M ) ∇ ⁡ 𝝎 | )

≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | 𝐃 ⁢ 𝝎 - ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 | ) .

The second term T 2 is also bounded by ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | 𝐃 ⁢ 𝝎 - ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 | ) since D ⁢ 𝐪 = ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 . Adding T 1 + T 2 and using the relation between φ and 𝓕 (Lemma 13), we have

⨍ M | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝝎 ) | 2 ≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) φ ϵ + | 𝐃𝐪 | ⁢ ( | 𝐃 ⁢ 𝝎 - ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 | )

(3.15) ≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 ) | 2 .

Furthermore, it is possible to show that

(3.16) ⨍ Λ ⁢ ( M ) | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( ⨍ Λ ⁢ ( M ) 𝐃 ⁢ 𝝎 ) | 2 ≲ ⨍ Λ ⁢ ( M ) | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - ⨍ Λ ⁢ ( M ) 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) | 2 ,

see details on this step on [18, pp. 66–67] (the anisotropy of the elements does not alter these steps). Inserting (3.16) into (3.15), we get

⨍ M | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝝎 ) | 2 ≲ h x h y ⁢ ⨍ Λ ⁢ ( M ) | 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - ⨍ Λ ⁢ ( M ) 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) | 2 ,

or equivalently since Λ ⁢ ( M ) and M are equally aniostropic,

∥ 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝝎 ) ∥ 2 2 ≲ h x h y ⁢ ∥ 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) - ⨍ Λ ⁢ ( M ) 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) ∥ 2 2 .

Note that ⨍ Λ ⁢ ( M ) 𝓕 ⁢ ( 𝐃 ⁢ 𝝎 ) is a patch-wise constant. Thus the anisotropic Poincaré–Wirtinger inequality (3.14) applies and we have the desired result. ∎

3.2 Well-Posedness

In this subsection we show that there exists a (unique) solution to the discrete problem (2.3), which is bounded by the data of the problem. In order to do so, we first show that a discrete inf-sup condition holds.

Lemma 11 (Discrete inf-sup condition).

For all q h ∈ Q h and all ν ′ ∈ ] 2 , ∞ [ it holds that

(3.17) β c ⁢ ∥ q h ∥ ν ′ ; Ω ≤ h x h y ⁢ ( sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( ∇ ⋅ 𝝎 h , q h ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ ν ; Ω + ∑ M ∈ 𝕄 ∑ | α | = 1 h α ⁢ ∥ θ h ⁢ ( D α ⁢ q h ) ∥ ν ′ , M ) ,

where β c is the constant of a continuous inf-sup condition satisfied by the spaces X and Q .

Proof.

We know that the pair 𝓧 × 𝒬 satisfies the continuous inf-sup condition (see e.g. [18]), i.e. there exists a positive constant β c so that

β c ⁢ ∥ q ∥ ν ′ ; Ω ≤ sup 𝝎 ∈ 𝐖 0 1 , ν ⁢ ( Ω ) ⁢ ( ∇ ⋅ 𝝎 , q ) Ω ∥ ∇ ⁡ 𝝎 ∥ ν ; Ω for all ⁢ q ∈ L 0 ν ′ ⁢ ( Ω ) ,

and since 𝒬 h ⊂ L 0 ν ′ ⁢ ( Ω ) , this holds also for all q h ∈ 𝒬 h . Now we split the right-hand side into two terms by adding and subtracting 𝐣 h ⁢ 𝝎 . The first term is then multiplied by ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ = 1 to get

β c ⁢ ∥ q h ∥ ν ′ ≤ sup 𝝎 ∈ 𝐖 0 1 , ν ⁢ ( Ω ) ⁡ ( ∇ ⋅ 𝐣 𝒉 ⁢ 𝝎 , q h ) Ω ⁢ ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ ν ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ ν ⁢ ∥ ∇ ⁡ 𝝎 ∥ ν + sup 𝝎 𝒉 ∈ 𝐖 0 1 , p ⁢ ( Ω ) ⁢ ( ∇ ⋅ ( 𝝎 - 𝐣 𝒉 ⁢ 𝝎 ) , q h ) Ω ∥ ∇ ⁡ 𝝎 ∥ ν for all ⁢ q h ∈ L 0 ν ′ ⁢ ( Ω ) .

The first term can be bounded by using the stability property of the interpolation operator (3.5):

sup 𝝎 ∈ 𝐖 0 1 , ν ⁢ ( Ω ) ⁢ ( ∇ ⋅ 𝐣 𝒉 ⁢ 𝝎 , q h ) Ω ⁢ ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ ν ∥ ∇ ⁡ 𝐣 h ⁢ 𝝎 ∥ ν ⁢ ∥ ∇ ⁡ 𝝎 ∥ ν ≲ h x h y ⁢ sup 𝝎 𝒉 ∈ 𝐗 h ⁢ ( ∇ ⋅ 𝝎 𝒉 , q h ) Ω ∥ ∇ ⁡ 𝝎 h ∥ ν .

Let us now estimate the nominator of the second term. Integrating by parts (remembering the Dirichlet boundary conditions, and that 𝐣 h preserves boundary conditions) followed by subtracting P h ⁢ ∇ ⁡ q h from ∇ ⁡ q h (exploiting the orthogonality property of 𝐣 h ) yields

| ( ∇ ⋅ ( 𝝎 - 𝐣 𝒉 𝝎 ) , q h ) Ω | = | ( 𝝎 - 𝐣 𝒉 𝝎 ) , ∇ q h ) Ω | = | ( ( 𝝎 - 𝐣 𝒉 𝝎 ) , 𝜽 h ( ∇ q h ) ) Ω | .

Separation of the components gives

| ( ∇ ⋅ ( 𝝎 - 𝐣 𝒉 ⁢ 𝝎 ) , q h ) Ω | ≲ ∑ M ∈ 𝕄 ( 1 h x ⁢ ∥ 𝝎 𝒙 - 𝐣 𝒉 ⁢ 𝝎 𝒙 ∥ ν , M ⁢ h x ⁢ ∥ θ h ⁢ ( ∂ x ⁡ q h ) ∥ ν ′ , M + 1 h y ⁢ ∥ 𝝎 𝒚 - 𝐣 𝒉 ⁢ 𝝎 𝒚 ∥ ν , M ⁢ h y ⁢ ∥ θ h ⁢ ( ∂ y ⁡ q h ) ∥ ν ′ , M ) .

Applying Hölder’s inequality and adding extra 1 h x ν ⁢ ∥ ω x - j h ⁢ ω x ∥ ν , M ν and 1 h y ν ⁢ ∥ ω y - j h ⁢ ω y ∥ ν , M ν , we get

| ( 𝝎 - 𝐣 𝒉 𝝎 ) , ∇ q h ) Ω | ≲ ( ∑ M ∈ 𝕄 1 h x ν ∥ ω x - j h ω x ∥ ν , M ν ) 1 ν ( ∑ M ∈ 𝕄 h x ν ′ ∥ θ h ( ∂ x q h ) ∥ ν ′ , M ν ′ ) 1 ν ′

+ ( ∑ M ∈ 𝕄 1 h y ν ⁢ ∥ ω y - j h ⁢ ω y ∥ ν , M ν ) 1 ν ⁢ ( ∑ M ∈ 𝕄 h y ν ′ ⁢ ∥ θ h ⁢ ( ∂ y ⁡ q h ) ∥ ν ′ , M ν ′ ) 1 ν ′

≲ ( ∑ M ∈ 𝕄 1 h x ν ⁢ ∥ ω x - j h ⁢ ω x ∥ ν , M ν + 1 h y ν ⁢ ∥ ω y - j h ⁢ ω y ∥ ν , M ν ) 1 ν ⁢ ( ∑ M ∈ 𝕄 ∑ | α | = 1 h α ⁢ ∥ θ h ⁢ ( D α ⁢ q h ) ∥ ν ′ , M ) .

Using that h x > h y , and the approximation properties of the 𝐣 h we can bound the terms in the sum of the first factor on the right-hand side as

1 h x ν ⁢ ∥ ω x - j h ⁢ ω x ∥ ν , M ν + 1 h y ν ⁢ ∥ ω y - j h ⁢ ω y ∥ ν , M ν ≲ h x ν h y ν ⁢ ( ∥ ∂ x ⁡ 𝝎 ∥ ν , M + ∥ ∂ y ⁡ 𝝎 ∥ ν , M + ∥ ∂ x ⁡ ω y ∥ ν , M + ∥ ∂ y ⁡ ω y ∥ ν , M ) ν

≤ ( h x h y ⁢ ∥ ∇ ⁡ 𝝎 ∥ ν , M ) ν .

Dividing by ∥ ∇ ⁡ 𝝎 ∥ ν , M completes the proof. ∎

Lemma 12 (Existence and Stability).

The solution of (2.5) exists, is unique and satisfies

(3.18) ∥ 𝐯 h ∥ 1 , p p + ∥ 𝐒 ⁢ ( 𝐃𝐯 h ) ∥ p ′ p ′ + s h ⁢ ( π h , π h ) ≤ C 1 ,

(3.19) β c ⁢ ∥ π h ∥ p ′ ≲ h x h y ⁢ C 2 ,

where the constant C 1 depend on Ω , f , p , ϵ , and C 2 additionally on α 0 and τ.

Proof.

The proof of existence and uniqueness of the solution is identical to the proof on isotropic meshes. The proof of the velocity estimate (3.18) is not altered although the mesh is anisotropic, and we refer the interested reader to the proof of [19, Lemma 3.6]. To bound the pressure, we start by inserting the discrete equation into the discrete inf-sup condition of Lemma 11:

(3.20) β c ⁢ ∥ π h ∥ p ′ ≲ h x h y ⁢ sup 𝝎 𝒉 ∈ 𝐗 h ⁢ ( 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ 𝝎 h ) Ω - ( 𝐟 , 𝝎 h ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ p + h x h y ⁢ ∑ M ∈ 𝕄 ∑ | α | = 1 h α ⁢ ∥ θ h ⁢ ( D α ⁢ π h ) ∥ p ′ , M .

We now treat the three terms of (3.20), i.e. the first term and the two components of the second term, separately. Using Hölder’s and Poincaré’s inequality, that ∥ D ⁢ u ∥ < C ⁢ ∥ ∇ ⁡ u ∥ and ∥ 𝐃𝐯 h ∥ p p ≳ ∥ 𝐒 ( 𝐃𝐯 h ∥ p ′ p ′ , we bound the first term by

sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ 𝝎 h ) Ω - ( 𝐟 , 𝝎 h ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ p ≲ sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( ∥ 𝐒 ⁢ ( 𝐃𝐯 h ) ∥ p ′ ⁢ ∥ ∇ ⁡ 𝝎 h ∥ p ∥ ∇ ⁡ 𝝎 𝒉 ∥ p + C ⁢ ( Ω , p ) ⁢ ∥ 𝐟 ∥ p ′ ⁢ ∥ ∇ ⁡ 𝝎 h ∥ p ∥ ∇ ⁡ 𝝎 𝒉 ∥ p )

≲ ∥ 𝐃𝐯 h ∥ p p - 1 + C ⁢ ( Ω , p ) ⁢ ∥ 𝐟 ∥ p ′ .

For the x-component of the second term of (3.20) we have considering that p ′ > 2 ,

∑ M ∈ 𝕄 h x ⁢ ∥ θ h ⁢ ( ∂ x ⁡ π h ) ∥ p ′ , M ≲ ( ∑ M ∈ 𝕄 h x p ′ ⁢ ∥ θ h ⁢ ( ∂ x ⁡ π h ) ∥ p ′ , M p ′ ) 1 p ′

≲ h x 1 - 2 p ′ ⁢ ( ∑ M ∈ 𝕄 ∫ M h x 2 ⁢ ( τ + | θ h ⁢ ( ∂ x ⁡ π h ) | ) p ′ - 2 ⁢ | θ h ⁢ ( ∂ x ⁡ π h ) | 2 ⏟ x ⁢ -component of ⁢ s ⁢ ( π h , π h ) ) 1 p ′ .

For the y-component we have

( ∑ M ∈ 𝕄 h y p ′ ⁢ ∥ θ h ⁢ ( ∂ y ⁡ π h ) ∥ p ′ , M p ′ ) 1 p ′ ≲ ( h y p ′ - 2 ⁢ ∑ M ∈ 𝕄 h y 2 ⁢ ∫ M ( τ 2 + | θ h ⁢ ( ∂ y ⁡ π h ) | p ′ - 2 ) ⁢ | θ h ⁢ ( ∂ y ⁡ π h ) | 2 ) 1 p ′

≲ h y 1 - 2 p ′ ⁢ ( h x h y ) 1 - 2 p ′ ⁢ ( ∫ M h y 2 ⁢ ( τ + h x h y ⁢ | θ h ⁢ ( ∂ y ⁡ π h ) | p ′ - 2 ) ⁢ | θ h ⁢ ( ∂ y ⁡ π h ) | p ′ , M 2 ⏟ y ⁢ -component of ⁢ s ⁢ ( π h , π h ) ) 1 p ′ .

By combining the three terms while considering that p ′ > 2 and inserting (3.18), we finally arrive at

h y h x ⁢ β c ⁢ ∥ π h ∥ p ′ ≲ ∥ 𝐃𝐯 h ∥ p p - 1 + ∥ 𝐟 ∥ p ′ + s ⁢ ( π h , π h ) 1 p ′ .

Estimate (3.19) now follows as a consequence of (3.18). ∎

3.3 A Priori Error Estimate

Theorem 1 (A Priori Error Estimates).

For p ∈ ] 1 , 2 ] , ϵ ∈ [ 0 , ϵ 0 ] , F ⁢ ( Dv ) ∈ W 1 , 2 ⁢ ( Ω ) d × d and π ∈ W 1 , p ′ ⁢ ( Ω ) the following a priori error estimates hold:

(3.21) ∥ 𝐯 - 𝐯 h ∥ 1 , p ≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 ≤ C v ⁢ h x ⁢ 𝒜 p 2 ,

(3.22) ∥ π - π h ∥ p ′ ≤ C p ⁢ h x 2 p ′ ⁢ 𝒜 p ,

with the aspect ratio A = h x h y ≫ 1 and p ′ = p p - 1 . The constants depend on α 0 , p, ϵ 0 , Ω, f, τ, ∥ ∂ x ⁡ π ∥ p ′ , ∥ ∂ y ⁡ ( π ) ∥ p ′ , ∥ ∂ x ⁡ F ⁢ ( Dv ) ∥ 2 , ∥ ∂ y ⁡ F ⁢ ( Dv ) ∥ 2 , and ∥ ∂ α ⁡ v x ∥ p , ∥ ∂ α ⁡ v y ∥ p with | α | = 2 . In the pressure estimate (3.22) the constant also depends on β c .

Proof.

We start with proving the velocity estimate (3.21). By Poincaré’s inequality, generalized Korn’s inequality (see e.g. [18, equation (2.1)]) and Lemma 13 in the Appendix the velocity error is can be bounded by

∥ 𝐯 - 𝐯 h ∥ 1 , p ≲ ∥ 𝐃𝐯 - 𝐃𝐯 h ∥ p

≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2

≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) ∥ 2 + ∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 .

The first term of the right-hand side is bounded by the interpolation estimate of Lemma 6. In order to bound the second term involving the projection error, we introduce the following quantity:

𝔈 lps 2 := ∥ 𝓕 ( 𝐃𝐣 h 𝐯 ) - 𝓕 ( 𝐃𝐯 h ) ∥ 2 2 + ∑ M ∈ 𝕄 ( h x 2 ∥ 𝓖 1 ( θ h ∂ x j h π ) - 𝓖 1 ( θ h ∂ x π h ) ∥ 2 ; M 2

+ ( h y h x ) p ′ - 2 h y 2 ∥ 𝓖 2 ( θ h ∂ y j h π ) - 𝓖 2 ( θ h ∂ y π h ) ∥ 2 ; M 2 ) ,

which due to Lemma 4 is equivalent to

(3.23) 𝔈 lps 2 ∼ ∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 2 + s ⁢ ( j h ⁢ π , j h ⁢ π - π h ) ⏟ = ⁣ : I 1 - s ⁢ ( π h , j h ⁢ π - π h ) .

Using Lemma 13 and that

𝐒 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 h ) = 𝐒 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 ) + 𝐒 ⁢ ( 𝐃𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 h )

together with Galerkin-orthogonality, we get

∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 2 ≲ ( 𝐒 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) )

= ( 𝐒 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 ) , 𝐃 ⁢ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) ) ⏟ = ⁣ : I 2 + ( π - π h , ∇ ⋅ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) )

+ ( ∇ ⋅ ( 𝐯 - 𝐯 h ) , j h ⁢ π - π h ) + s ⁢ ( π h , j h ⁢ π - π h ) .

Some algebraic manipulations on the second and third term (exploiting that ∇ ⋅ 𝐯 = 0 ) and inserting into (3.23) yields

𝔈 lps 2 ≲ I 1 + I 2 + ( π - j h ⁢ π , ∇ ⋅ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) ) ⏟ = ⁣ : I 3 + ( j h ⁢ π - π h , ∇ ⋅ ( 𝐣 h ⁢ 𝐯 - 𝐯 ) ) ⏟ = ⁣ : I 4 .

We now proceed to bound each of the four terms I 1 , I 2 , I 3 , and I 4 .

For the term I 1 we use Lemma 4 and the definition of the 𝔈 lps (considering that h x ≥ h y and p ′ ≥ 2 ):

I 1 ≤ δ 1 ⁢ 𝔈 lps 2 + c δ 1 ⁢ ( h x 2 ⁢ ∥ τ + | ∂ x ⁡ ( j h ⁢ π ) | ∥ p ′ ; Ω p ′ + h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ τ 2 + | ∂ y ⁡ ( j h ⁢ π ) | ∥ p ′ ; Ω p ′ ) .

Note that, in order for Lemma 4 to be valid, the additional regularity of the pressure, π ∈ W 1 , p ′ , is required. The stability estimates (3.6)–(3.7) and τ 2 = h x h y ⁢ τ yields

I 1 ≤ δ 1 ⁢ 𝔈 lps 2 + c δ 1 ⁢ ( h x 2 ⁢ τ p ′ ⁢ | Ω | + h x 2 ⁢ ∥ ∂ x ⁡ ( j h ⁢ π ) ∥ p ′ ; Ω p ′ + h y 2 ¯ ⁢ ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ ∂ y ⁡ ( j h ⁢ π ) ∥ p ′ ; Ω p ′ )

≤ δ 1 ⁢ 𝔈 lps 2 + c δ 1 ⁢ ( h x 2 ⁢ τ p ′ ⁢ | Ω | + ( h x 2 + h y 2 ¯ ⁢ ⁢ ( h y h x ) p ′ - 2 ⁢ ( h x h y ) p ′ ) ⁢ ∥ ∂ x ⁡ π ∥ p ′ p ′ ) + c δ 1 ⁢ ( h x 2 ⁢ ( h y h x ) p ′ + h y 2 ⁢ ⁢ ( h y h x ) p ′ - 2 ) ⁢ ∥ ∂ y ⁡ π ∥ p ′ p ′

≤ δ 1 ⁢ 𝔈 lps 2 + 2 ⁢ c δ 1 ⁢ ( h x 2 ⁢ τ p ′ ⁢ | Ω | + h x 2 ⁢ ∥ ∂ x ⁡ π ∥ p ′ p ′ + h x 2 ⁢ ( h y h x ) p ′ ⁢ ∥ ∂ y ⁡ π ∥ p ′ p ′ )

≲ δ 1 ⁢ 𝔈 lps 2 + c δ 1 ⁢ h x 2 ⁢ ( 1 + 𝒜 - p ′ ) .

The factors underlined with a wiggly line stems from the inverse aspect ratio introduced in the non-linear part of the stabilization, see (2.7). These balances the aspect ratio factor arising from the stability estimates of j h , to prevent the error from being proportional to 𝒜 p ′ - 2 . The terms underlined by a straight line stems from the h y 2 factor in (2.7), which clearly prevents terms containing h x 2 ⁢ ∥ ∂ y ⁡ π ∥ p ′ to enter.

To bound I 2 , we follow the steps in the proof of [19, Theorem 3.1], namely applying Lemma 13 and using Young’s inequality – these steps are identical on anisotropic meshes. We then apply the interpolation property (3.11) for 𝓕 , assuming that 𝓕 ⁢ ( 𝐃𝐯 ) ∈ W 1 , 2 ⁢ ( Ω ) d × d to get

I 2 ≲ c δ 2 ⁢ ∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 ) ∥ 2 2 + δ 2 ⁢ ∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 2

≲ δ 2 ⁢ 𝔈 lps 2 + c δ 2 ⁢ ( h x h y ⁢ h x 2 ⁢ ∥ ∂ x ⁡ 𝓕 ⁢ ( 𝐃𝐯 ) ∥ 2 2 + h x h y ⁢ h y 2 ⁢ ∥ ∂ y ⁡ 𝓕 ⁢ ( 𝐃𝐯 ) ∥ 2 2 )

≲ δ 2 ⁢ 𝔈 lps 2 + c δ 2 ⁢ h x 2 ⁢ ( 𝒜 + 𝒜 - 1 ) .

For bounding the third term I 3 , we first use generalized Korn’s inequality, Lemma 13 and the definition of 𝔈 lps to obtain the following estimate:

(3.24) ∥ ∇ ⁡ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) ∥ p ≤ ∥ 𝓕 ⁢ ( 𝐃𝐣 h ⁢ 𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 ⁢ ∥ ϵ + | 𝐃𝐣 h ⁢ 𝐯 | + | 𝐃𝐯 h | ∥ p 1 - p 2 ≤ 𝔈 lps ⁢ ∥ ϵ + | 𝐃𝐣 h ⁢ 𝐯 | + | 𝐃𝐯 h | ∥ p 1 - p 2 .

Applying Hölder’s inequality and Young’s inequality to I 3 , and then inserting (3.24) and the interpolation property on the pressure term (assuming π ∈ W 1 , p ′ ), we have

I 3 ≤ ∥ j h ⁢ π - π ∥ p ′ ⁢ ∥ ∇ ⁡ ( 𝐣 h ⁢ 𝐯 - 𝐯 h ) ∥ p

≤ δ 3 ⁢ 𝔈 lps 2 + c δ 3 ⁢ ∥ ϵ + | 𝐃𝐣 h ⁢ 𝐯 | + | 𝐃𝐯 h | ∥ p 2 - p ⁢ ∥ j h ⁢ π - π ∥ p ′ 2

≤ δ 3 ⁢ 𝔈 lps 2 + c δ 3 ⁢ ( ϵ ⁢ | Ω | 1 p + ∥ 𝐃𝐣 h ⁢ 𝐯 ∥ p + ∥ 𝐃𝐯 h ∥ p ) 2 - p ⁢ ( h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ ) 2

≤ δ 3 ⁢ 𝔈 lps 2 + c δ 3 ⁢ ( ϵ ⁢ | Ω | 1 p + ( 1 + h x h y ) ⁢ ∥ ∇ ⁡ 𝐯 ∥ p + ∥ ∇ ⁡ 𝐯 h ∥ p ) 2 - p ⁢ ( h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ ) 2 .

The stability of the interpolation operator, the stability of the continuous solution, and Lemma 12 give ∥ ∇ ⁡ 𝐣 h ⁢ 𝐯 ∥ p ≲ h x h y ⁢ C 1 ⁢ ( 𝐟 , ϵ , Ω , p ) and ∥ ∇ ⁡ 𝐯 h ∥ p ≤ C 2 ⁢ ( 𝐟 , ϵ , Ω , p ) . Using these bounds, we get

I 3 ≤ δ 3 ⁢ 𝔈 lps 2 + c δ 3 ⁢ ( h x h y ) 2 - p ⁢ ( h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ ) 2

≲ δ 3 ⁢ 𝔈 lps 2 + c δ 3 ⁢ h x 2 ⁢ ( 𝒜 2 - p + 𝒜 - p ) .

For the last term, I 4 , we use integration by parts, the orthogonality of 𝐣 h , and the triangle inequality to split into components

I 4 = | ( 𝐣 h ⁢ 𝐯 - 𝐯 , ∇ ⁡ ( j h ⁢ π - π h ) ) |

= | ( 𝐣 h ⁢ 𝐯 - 𝐯 , 𝜽 h ⁢ ( ∇ ⁡ ( j h ⁢ π - π h ) ) ) |

≲ | ( j h ⁢ v x - v x , θ h ⁢ ( ∂ x ⁡ ( j h ⁢ π - π h ) ) ) | + | ( j h ⁢ v y - v y , θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ) | .

The x-component is bounded by applying Hölder’s inequality followed by Young’s inequality

I 4 ⁢ x := | ( j h ⁢ v x - v x , θ h ⁢ ( ∂ x ⁡ ( j h ⁢ π - π h ) ) ) |

≲ ∑ M ∈ 𝕄 h h x - 2 p ′ ∥ j h v x - v x ∥ p ; M h x 2 p ′ ∥ θ h ( ∂ x ( j h π - π h ) ) ) ∥ p ′ ; M

≲ ∑ M ∈ 𝕄 h c δ 4 h x - 2 ⁢ ( p - 1 ) ∥ j h v x - v x ∥ p ; M p + δ 4 ∑ M ∈ 𝕄 h h x 2 ∥ θ h ( ∂ x ( j h π - π h ) ) ) ∥ p ′ p ′ ; M .

Using the interpolation properties of j h , the regularity 𝐯 ∈ W 2 , p , and the fact that

∥ θ h ( ∂ x ( j h π - π h ) ) ) ∥ p ′ p ′ ; M ≤ ∥ 𝓖 1 ( θ h ( ∂ x j h π ) ) - 𝓖 1 ( θ h ( ∂ x π h ) ) ∥ 2 2 ; M ,

which is in turn bounded by 𝔈 lps , we arrive at

I 4 ⁢ x ≲ δ 4 ⁢ 𝔈 lps 2 + c δ 4 ⁢ h x 2 - 2 ⁢ p ⁢ ( h x 2 ⁢ ∥ ∂ x ⁢ x 2 ⁡ v x ∥ p ; Ω + h y 2 ⁢ ∥ ∂ y ⁢ y 2 ⁡ v x ∥ p ; Ω + h x ⁢ h y ⁢ ∥ ∂ x ⁢ y 2 ⁡ v x ∥ p ; Ω ) p

≲ δ 4 ⁢ 𝔈 lps 2 + c δ 4 ⁢ h x 2 ⁢ ( 1 + 𝒜 - 2 ⁢ p + 𝒜 - p ) .

Similarly for the y-component, I 4 ⁢ y := | ( j h ⁢ v y - v y , θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ) | , we have

I 4 ⁢ y ≲ ∑ M ∈ 𝕄 h h y - 2 p ′ ⁢ ( h y h x ) 2 - p ′ p ′ ⁢ ∥ j h ⁢ v y - v y ∥ p ; M ⁢ h y 2 p ′ ⁢ ( h y h x ) p ′ - 2 p ′ ⁢ ∥ θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ ; M

≲ c δ 4 ⁢ ∑ M ∈ 𝕄 h h y - 2 ⁢ ( p - 1 ) ¯ ⁢ ⁢ ( h y h x ) p - 2 ⁢ ∥ j h ⁢ v y - v y ∥ p ; M p + δ 4 ⁢ ∑ M ∈ 𝕄 h h y 2 ¯ ⁢ ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ ; M p ′

≲ δ 4 ⁢ 𝔈 lps 2 + c δ 4 ⁢ ⁢ h y - 2 ⁢ ( p - 1 ) ⁢ ⁢ ( h y h x ) p - 2 ⁢ ( h x 2 ⁢ ∥ ∂ x ⁢ x 2 ⁡ v y ∥ p ; Ω + h y 2 ⁢ ∥ ∂ y ⁢ y 2 ⁡ v y ∥ p ; Ω + h x ⁢ h y ⁢ ∥ ∂ x ⁢ y 2 ⁡ v y ∥ p ; Ω ) p

≲ δ 4 ⁢ 𝔈 lps 2 + c δ 4 ⁢ h y - p ⁢ h x 2 - p ⁢ ( h x 2 ⁢ ∥ ∂ x ⁢ x 2 ⁡ v y ∥ p ; Ω + h y 2 ⁢ ∥ ∂ y ⁢ y 2 ⁡ v y ∥ p ; Ω + h x ⁢ h y ⁢ ∥ ∂ x ⁢ y 2 ⁡ v y ∥ p ; Ω ) p

≲ δ 4 ⁢ 𝔈 lps 2 + c δ 4 ⁢ h x 2 ⁢ ( 𝒜 p + 𝒜 - p + 1 ) .

Again the factors underlined with a wiggly line stems from the inverse aspect ratio introduced in the non-linear part of the stabilization and the terms with a straight lines stems from choosing h y 2 instead of h x 2 in the y-component of (2.7).

Collecting I 1 - I 4 and combining constants independent of the element aspect ratio and stability parameter, and choosing δ 1 , δ 2 , δ 3 , δ 4 small enough, we get

𝔈 lps 2 ≲ h x 2 ⁢ ( 1 + 𝒜 + 𝒜 - 1 + 𝒜 - p ′ + 𝒜 2 - p + 𝒜 - p + 𝒜 - 2 ⁢ p + 𝒜 - p + 𝒜 p ) .

Since we assumed in this work an aspect ratio 𝒜 > 1 , the dominant term is 𝒜 p . Hence, we arrived at

(3.25) 𝔈 lps ≲ h x ⁢ 𝒜 p 2 ,

where the constant depends on α 0 , p, ϵ 0 , Ω, f and τ. This is exactly the assertion (3.21).

To arrive at the pressure error estimate (3.22), we first notice that the projection error j h ⁢ π - π h is in 𝒬 h , so that we can use the discrete inf-sup condition of (3.17) to get

𝒜 - 1 ⁢ β c ⁢ ∥ j h ⁢ π - π h ∥ p ′ ≲ sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( ∇ ⋅ 𝝎 h , ( j h ⁢ π - π h ) ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ p + ( ∑ M ∈ 𝕄 h x p ′ ⁢ ∥ θ h ⁢ ( ∂ x ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ , M p ′ ) 1 p ′

+ ( ∑ M ∈ 𝕄 h y p ′ ⁢ ∥ θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ , M p ′ ) 1 p ′

= : J 1 + J 2 + J 3 .

By inserting the continuous and discrete equation, using Hölder’s inequality, followed by (B.1) and the approximation properties of j h (3.1), we get

J 1 ≲ sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( 𝐒 ⁢ ( 𝐃𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 h ) , 𝐃 ⁢ ω h ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ p + sup 𝝎 𝒉 ∈ 𝐗 h ⁡ ( j h ⁢ π - π , ∇ ⋅ 𝝎 h ) Ω ∥ ∇ ⁡ 𝝎 𝒉 ∥ p

≲ ∥ 𝐒 ⁢ ( 𝐃𝐯 ) - 𝐒 ⁢ ( 𝐃𝐯 h ) ∥ p ′ + ∥ j h ⁢ π - π ∥ p ′

≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐯 h ) ∥ 2 2 p ′ + h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′

≲ 𝔈 lps 2 p ′ + h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ .

The second and third term J 2 and J 3 are estimated by considering that h x ≥ h y , and p ′ ≥ 2 , using Lemma 2, and the definition of 𝔈 lps 2 / p ′ :

J 2 ≲ h x 1 - 2 p ′ ⁢ ( ∑ M ∈ 𝕄 h x 2 ⁢ ∥ θ h ⁢ ( ∂ x ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ , M p ′ ) 1 p ′

≲ h x 1 - 2 p ′ ⁢ ( ∑ M ∈ 𝕄 h x 2 ⁢ ∥ 𝓖 1 ⁢ ( θ h ⁢ ( ∂ x ⁡ j h ⁢ π ) ) - 𝓖 1 ⁢ ( θ h ⁢ ( ∂ x ⁡ π h ) ) ∥ 2 ; M 2 ) 1 p ′

≲ h x 1 - 2 p ′ ⁢ 𝔈 lps 2 p ′

and

J 3 ≲ h y 1 - 2 p ′ ⁢ ( ∑ M ∈ 𝕄 h y 2 ⁢ ∥ θ h ⁢ ( ∂ y ⁡ ( j h ⁢ π - π h ) ) ∥ p ′ , M p ′ ) 1 p ′

≲ h x 1 - 2 p ′ ⁢ ( ∑ M ∈ 𝕄 h y 2 ⁢ ( h y h x ) p ′ - 2 ⁢ ∥ 𝓖 2 ⁢ ( θ h ⁢ ( ∂ y ⁡ j h ⁢ π ) ) - 𝓖 2 ⁢ ( θ h ⁢ ( ∂ y ⁡ π h ) ) ∥ 2 ; M 2 ) 1 p ′

≲ h x 1 - 2 p ′ ⁢ 𝔈 lps 2 p ′ .

Combining the estimates of J 1 and J 2 with the interpolation error, we get the following estimate of the discretization-error in the pressure:

∥ π - π h ∥ p ′ ≤ ∥ j h ⁢ π - π ∥ p ′ + ∥ j h ⁢ π - π h ∥ p ′

≲ ( 1 + 𝒜 β c ) ⁢ ( h x ⁢ ∥ ∂ x ⁡ π ∥ p ′ + h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ ) + 𝒜 β c ⁢ ( 1 + h x 1 - 2 p ′ ) ⁢ 𝔈 lps 2 p ′ .

Inserting (3.25) gives (3.22):

∥ π - π h ∥ p ′ ≲ 𝒜 ⁢ h x + 𝒜 ⁢ ( 1 + h x 1 - 2 p ′ ) ⁢ ( h x ⁢ 𝒜 p 2 ) 2 p ′

≲ 𝒜 ⁢ h x + 𝒜 1 + p / p ′ ⁢ ( h x 2 p ′ + h x )

≲ 𝒜 p ⁢ h x 2 p ′ . ∎

Remark.

For a stabilization without the h y h x -factor in 𝐌 s , 22 in (2.8), i.e. for a stabilization on the form

s ( π h , q h ) := α 0 ∑ M ∈ 𝕄 ( 𝐌 s 𝜽 h ∇ π h , 𝜽 h ∇ q h ) M ) ,

𝐌 s ⁢ ( ∂ x ⁡ π h , ∂ x ⁡ π h ) := [ h x 2 ⁢ ( 1 τ ⁢ ( τ + | θ h ⁢ ∂ x ⁡ π h | ) ) p ′ - 2 0 0 h y 2 ⁢ ( 1 τ ⁢ ( τ + | θ h ⁢ ∂ y ⁡ π h | ) ) p ′ - 2 , ] ,

we get, by setting the factors marked by awiggly line equal to one

𝔈 lps ≲ h x ⁢ ( 𝒜 p ′ - 2 2 + 𝒜 p - 1 ) .

The first term dominates for p < 3 2 , and as p → 1 we have that 1 2 ⁢ ( p ′ - 2 ) → ∞ . In this sense, the inaccuracies related to the anisotropy of the mesh is amplified by the non-linearity of the material. For an isotropic stabilization, i.e. with

(3.26) 𝐌 s ⁢ ( ∂ x ⁡ π h , ∂ x ⁡ π h ) = [ h x 2 ⁢ ( 1 τ ⁢ ( τ + | θ h ⁢ ∂ x ⁡ π h | ) ) p ′ - 2 0 0 h x 2 ⁢ ( 1 τ ⁢ ( τ + | θ h ⁢ ∂ y ⁡ π h | ) ) p ′ - 2 ] ,

we replace h y by h x in the terms underlined with a straight line (as well as setting the factors underlined by wiggly lines equal to one) in the estimates of I 1 and I 4 to see that

𝔈 lps ≲ h x ⁢ ( 𝒜 p ′ 2 ) .

Also terms with h x ⁢ ∥ ∂ y ⁡ π ∥ p ′ ⁣ p ′ / 2 will now arise from I 1 , where ∥ ∂ y ⁡ π ∥ p ′ is typically large in a thin domain. The anisotropic stabilization instead involves more well-balanced terms h y ⁢ ∥ ∂ y ⁡ π ∥ p ′ ⁣ p ′ / 2 .

4 Numerical Experiments

We present numerical results for two different configurations. The first experiment is designed to show how the convergence rate depends on p. For this purpose, we have chosen a problem for which the projection error, not the interpolation error, dominates the pressure error. The second experiment is designed to compare the robustness and accuracy of the isotropic and anisotropic stabilization. It is constructed to qualitatively mimic the solution of a typical thin film flow scenario seen in e.g. glaciology but still be regular enough to allow for a wide spectrum of stabilization parameters α 0 to be tested. For both experiments, we have chosen the highest aspect ratio 𝒜 that the (non-specialized) linear solver allows for. The first configuration has an aspect ratio of the cell sizes of 𝒜 = 100 , the second one of 𝒜 = 1000 . The regularization parameter in (1.2) is chosen as ϵ = 10 - 5 .

4.1 Experiment 1 – Convergence

We study how the numerical error decreases as the resolution increases. The experiment is similar to the convergence experiment in [18], but adapted to suit anisotropic problems. The domain is Ω := ( - L 2 , L 2 ) × ( - H 2 , H 2 ) , with L = 1 and H = 0.01 , and the parameter scaling the viscosity is μ 0 = 1.0 . The problem is chosen so that the exact solution is

𝐯 = ( L ⁢ ( ( x L ) 2 + ( y H ) 2 ) a - 1 2 ⁢ y H - H ⁢ ( ( x L ) 2 + ( y H ) 2 ) a - 1 2 ⁢ x L ) , π = - ( ( x L ) 2 + ( y H ) 2 ) b 2 ⁢ x L ⁢ y H .

In order for the velocity to fulfill the regularity assumption 𝓕 ⁢ ( 𝐃𝐯 ) ∈ L 2 , it is required that a > 1 . Following [18], we set a = 1.01 . For the pressure to be regular enough, we need b > - 1 - p ′ 2 . We set b = 0.1 , which fulfills the regularity requirement but is low enough to ensure that the pressure error is dominated by the projection error ( ∼ h 2 p ′ ), not the interpolation error ( ∼ h 2 ) for the mesh resolutions tested here. The right-hand side is given as 𝐟 = - ∇ ⋅ 𝐒 ⁢ ( 𝐃𝐯 ) + ∇ ⁡ π and the Dirichlet boundary condition is 𝐠 = 𝐯 | ∂ ⁡ Ω . We set α 0 = 0.01 , and τ = 1.0 .

Table 1

Results of Experiment 1, verification of Theorem 1.

	N × N	∥ p - p h ∥ p ′	conv.	∥ v x - v x , h ∥ 1 , p	conv.	∥ v y - v y , h ∥ 1 , p	conv.
p = 1.1	16	2.94 ⁢ e - 03	5.06	3.57 ⁢ e - 03	0.98	1.50 ⁢ e - 05	1.10
	64	3.03 ⁢ e - 04	3.28	1.78 ⁢ e - 03	1.01	7.07 ⁢ e - 06	1.09
	256	7.32 ⁢ e - 05	2.05	8.75 ⁢ e - 04	1.02	3.58 ⁢ e - 06	0.98
	1024	1.77 ⁢ e - 05	2.05	4.33 ⁢ e - 04	1.02	1.79 ⁢ e - 06	1.00
	4096	1.14 ⁢ e - 05	0.63	2.15 ⁢ e - 04	1.01	8.47 ⁢ e - 07	1.08
	16384	1.01 ⁢ e - 05	0.17	1.07 ⁢ e - 04	1.01	4.29 ⁢ e - 07	0.98
	65536	8.96 ⁢ e - 06	0.18	5.33 ⁢ e - 05	1.00	2.32 ⁢ e - 07	0.89
	theory		0.18		1		1
p = 1.5	16	2.74 ⁢ e - 03	2.88	1.31 ⁢ e - 02	0.90	5.32 ⁢ e - 05	0.95
	64	1.22 ⁢ e - 04	4.49	6.80 ⁢ e - 03	0.94	2.73 ⁢ e - 05	0.96
	256	9.03 ⁢ e - 05	0.43	3.48 ⁢ e - 03	0.97	1.45 ⁢ e - 05	0.92
	1024	1.67 ⁢ e - 05	2.43	1.77 ⁢ e - 03	0.98	7.67 ⁢ e - 06	0.91
	4096	8.48 ⁢ e - 06	0.98	8.94 ⁢ e - 04	0.98	4.00 ⁢ e - 06	0.94
	16384	5.11 ⁢ e - 06	0.73	4.50 ⁢ e - 04	0.99	1.98 ⁢ e - 06	1.01
	65536	3.14 ⁢ e - 06	0.70	2.26 ⁢ e - 04	0.99	9.92 ⁢ e - 07	1.00
	theory		0.67		1		1
p = 1.9	16	6.32 ⁢ e - 03	0.52	2.89 ⁢ e - 02	0.83	1.16 ⁢ e - 04	0.83
	64	3.77 ⁢ e - 04	4.07	1.58 ⁢ e - 02	0.87	6.42 ⁢ e - 05	0.85
	256	2.62 ⁢ e - 04	0.53	8.49 ⁢ e - 03	0.90	3.57 ⁢ e - 05	0.85
	1024	4.70 ⁢ e - 05	2.48	4.50 ⁢ e - 03	0.92	1.97 ⁢ e - 05	0.86
	4096	1.79 ⁢ e - 05	1.39	2.36 ⁢ e - 03	0.93	1.08 ⁢ e - 05	0.86
	16384	8.81 ⁢ e - 06	1.02	1.23 ⁢ e - 03	0.94	5.80 ⁢ e - 06	0.90
	65536	4.45 ⁢ e - 06	0.98	6.38 ⁢ e - 04	0.95	2.98 ⁢ e - 06	0.96
	theory		0.95		1		1

The number of elements in each dimension is N, so that the element aspect ratio 𝒜 = h x h y = 100 is independent of N and equals the domain aspect ratio. As the convergence rate of the pressure is expected to depend on the parameter p, we test three different cases, namely p = 1.1 , p = 1.5 , and p = 1.9 on seven different meshes, with N = { 16 , 64 , 256 , 1034 , 4096 , 16384 , 65536 } . The errors ∥ p - p h ∥ p ′ , ∥ v x - v x , h ∥ p , ∥ v y - v y , h ∥ p and their respective order of convergence are shown in Table 1. For the velocity error we observe convergence order close to 𝒪 ⁢ ( h ) for all choices of p and all meshes which corresponds quite well with the theoretical order. The error estimate in Theorem 1 with respect to the pressure π gives ∥ p - p h ∥ p ′ = 𝒪 ⁢ ( h 2 p ′ ) = 𝒪 ⁢ ( h 2 ⁢ ( 1 p - 1 ) ) . In the worst case, p = 1.1 , the numerical order decays on finer meshes, but stabilizes at O ⁢ ( h 0.18 ) which exactly corresponds with the theoretical order. For moderate p = 1.9 , the observed order of convergence in the pressure is on the finest mesh O ⁢ ( h 0.98 ) which is slightly above the theoretical order O ⁢ ( h 0.95 ) . Hence, we state that, as the mesh is refined, the convergence order approaches the expected, theoretical convergence order.

4.2 Experiment 2 – Robustness

The anisotropic stabilization was designed to balance the order of magnitude of the x- and y-component of the stabilization term. We thus expect that the anisotropic stabilization increases the range of stabilization parameters α 0 that does not lead to under-nor over-stabilization, compared to an isotropic stabilization. To verify this, we run simulations with the anisotropic stabilization (2.8), a fully isotropic stabilization (3.26), and a isotropic stabilization which still includes the factor 𝒜 2 - p ′ that was introduced to suppress error terms proportional to 𝒜 p ′ or 𝒜 p ′ - 2 , i.e. a stabilization on the form

(4.1) s ( π h , q h ) := α 0 ∑ M ∈ 𝕄 ( 𝐌 s 𝜽 h ∇ π h , 𝜽 h ∇ q h ) M ) ,

(4.2) 𝐌 s ⁢ ( ∂ x ⁡ π h , ∂ x ⁡ π h ) := [ h x 2 ⁢ ( 1 τ ⁢ ( τ + h y h x ⁢ | θ h ⁢ ∂ x ⁡ π h | ) ) p ′ - 2 0 0 h x 2 ⁢ ( 1 τ ⁢ ( τ + h y h x ⁢ | θ h ⁢ ∂ y ⁡ π h | ) ) p ′ - 2 ] .

We will call this the semi-isotropic stabilization. For each method, we vary α 0 and measure the error in velocity and pressure.

The domain is Ω := ( 0 , L ) × ( 0 , H ) , with L = 1 and two different choices H = 0.001 and H = 0.01 . We choose a right-hand side 𝐟 = - ∇ ⋅ 𝐒 ⁢ ( 𝐃𝐯 ) + ∇ ⁡ π and Dirichlet boundary conditions 𝐠 = 𝐯 | ∂ ⁡ Ω so that the exact solution ( 𝐯 , π ) is

𝐯 = ( sin ( f x L ) cos ( f y H ) - H L ⁢ cos ⁡ ( f ⁢ x L ) ⁢ sin ⁡ ( f ⁢ y H ) ) , π = c p ⁢ sin ⁡ ( π ⁢ x L ) ⁢ cos ⁡ ( π ⁢ π ⁢ y H ) ,

where c p = 100 and f = 0.01 ⁢ π . The large pressure component and small vertical component is typical of many non-Newtonian thin-film problems, found for instance in glaciology. We set p = 1.4 , μ 0 = 0.1 , and τ = c p . The number of intervals in each direction is N = 32 . The mesh and problem is chosen so that the solver converges for a large range of α 0 for both, isotropic and anisotropic stabilization.

The fully isotropic stabilization (3.26) does not converge. This is indeed the reason why we implemented the semi-isotropic stabilization (4.1). The results for both the semi-isotropic and the anisotropic stabilizations are shown in Figures 2 and 3 for H = 10 - 3 and H = 10 - 2 ( 𝒜 = 1000 and 𝒜 = 100 ), respectively. For both methods a too small stabilization parameter α 0 leads to under-stabilization of the pressure, while a too large α 0 leads to errors in the velocity components, especially in the vertical component. The errors in the pressure are due to the spurious oscillations associated to element choices that do not fulfill the inf-sup condition. The errors in the velocity components are due to the consistency error caused by the extra stabilization term. This consistency error is present to some degree for all α 0 > 0 , but we only call the problem over-stabilized when it is the dominant error. The vertical velocity component is evidently more sensitive to over-stabilization than the horizontal velocity component. This sensitivity of the vertical component is typical of flows in thin domains, and are observed in applications such as glaciology [17].

For the semi-isotropic stabilization, there is no choice of α 0 that does not lead to neither under-stabilization of the pressure nor over-stabilization of the velocity, meaning that there is no α 0 which is high enough to suppress spurious oscillations in the pressure at the same time as it is low enough to avoid significant consistency errors in the velocity components. In contrast to this, the anisotropic stabilization is less sensitive to the choice of α 0 and it is possible to choose an α 0 which yields a stable pressure solution while at the same time avoiding large over-stabilization errors in the velocity. The optimal stabilization parameter for the anisotropic stabilization parameter is rather small (around 10 - 4 ), while for the semi-isotropic stabilization it is even smaller (around 10 - 6 ). For anisotropic stabilization, the optimal value of α 0 seems to be fairly independent of the aspect ratio. Note that the computational domains differ between the two cases 𝒜 = 1000 and 𝒜 = 100 . In practice, the optimal stabilization parameter is of course not known, but the moderate value for the anisotropic stabilization is more reasonable. For common choices of α 0 = 0.1 to α 0 = 1 , the error in the vertical velocity is reduced by two to three order of magnitudes by using the anisotropic stabilization instead of the semi-isotropic one. The problem of choosing optimal values of α 0 with the (semi-)isotropic stabilization methods (e.g. LPS, Galerkin-Least-Squares (GLS), Interior Penalty (IP)) has also been reported in [17], where the accuracy in pressure is quite sensitive with respect to α 0 even for moderate aspect ratios ( ≈ 80 ) for a different configuration. Hence, anisotropic stabilization is beneficial to be used.

$Figure 2 Results of Experiment 2 with aspect ratio 1000. The error as a function of the stabilization parameter α 0 {\alpha_{0}} for both semi-isotropic (left) as in (3.26) and anisotropic stabilization (right) as in (2.8). The fully isotropic stabilization did not converge. The range of α 0 {\alpha_{0}} for which the problem is under-stabilized or over-stabilized is indicated by ↔ {\leftrightarrow} .$

Figure 2

Results of Experiment 2 with aspect ratio 1000. The error as a function of the stabilization parameter α 0 for both semi-isotropic (left) as in (3.26) and anisotropic stabilization (right) as in (2.8). The fully isotropic stabilization did not converge. The range of α 0 for which the problem is under-stabilized or over-stabilized is indicated by ↔ .

$Figure 3 Results of Experiment 2 with aspect ratio 100. The error as a function of the stabilization parameter α 0 {\alpha_{0}} for both semi-isotropic (left) as in (3.26) and anisotropic stabilization (right) as in (2.8). The fully isotropic stabilization did not converge. The range of α 0 {\alpha_{0}} for which the problem is under-stabilized or over-stabilized is indicated by ↔ {\leftrightarrow} .$

Figure 3

Results of Experiment 2 with aspect ratio 100. The error as a function of the stabilization parameter α 0 for both semi-isotropic (left) as in (3.26) and anisotropic stabilization (right) as in (2.8). The fully isotropic stabilization did not converge. The range of α 0 for which the problem is under-stabilized or over-stabilized is indicated by ↔ .

5 Summary and Conclusion

We have derived a priori estimates for an equal-order bi-linear finite element discretization of the p-Stokes equations on anisotropic meshes, with power-law exponent p ∈ ( 1 , 2 ) . To stabilize the inf-sup unstable equal-order discretization, we designed an anisotropic local projection stabilization for the p-Stokes equations. Without a specially designed stabilization term errors related to the anisotropy of the problem is amplified by the non-linearity of the p-Stokes equations.

Our error estimates where confirmed in numerical experiments. They show in particular that our anisotropic stabilization is more robust with respect to user defined parameters than an semi-isotropic stabilization, while a fully isotropic standard stabilization did not converge. For typical parameter choices, an anisotropic stabilization reduces errors the velocity by several orders of magnitudes. Such results are relevant to applications in e.g. glaciology and industrial material processing.

Funding statement: This project was funded by the Cluster of Excellence 80 “The Future Ocean”. The “Future Ocean” is funded within the framework of the Excellence Initiative by the Deutsche Forschungsgemeinschaft (DFG) on behalf of the German federal and state governments.

A Proof of Lemma 1

Proof.

The proof is an extension of the proof of [21, Lemma 3.2]. For the meshes we consider | det D F K | = h x h y and | det D F M | ∼ h x h y . We then have

∥ q h ∥ ν ′ , M ν ′ = ∫ M q ( x ) ν ′ d x = ∫ M ^ q ^ h ( x ^ ) ν ′ | det D F K | d x ^ ∼ h x h y ∥ q h ^ ∥ ν ′ , M ν ′ .

To estimate ( ω h , q h ) M , let b ^ be the hat function associated with the node in the middle of a patch M and choose ω h ⁢ ( x ) := ( q h ^ ⋅ b ^ ) ∘ F M - 1 ⁢ ( x ) , where q h ∈ Y h ⁢ ( M ) and q h ^ ∈ Q 0 ⁢ ( M ^ ) . Then ω ^ h ⁢ ( x ) := ( q h ^ ⋅ b ^ ) ∈ X ^ h 0 ⁢ ( M ^ ) since ( q ^ h ⋅ ω ^ h ) is continuous on the closure of M ^ , ( q ^ h ⋅ ω ^ h ) | K ^ ∈ Q 1 ⁢ ( K ^ ) when K ^ belongs to the patch M ^ , and b ^ | ∂ ⁡ M ^ = 0 . Therefore ( ω ^ h , q ^ h ) M ^ = ( q ^ h ⁢ ( x ^ ) , b ^ ⁢ ( x ^ ) ⁢ q ^ h ⁢ ( x ^ ) ) M ^ . Since b is positive ( q ^ h ⁢ ( x ^ ) , q ^ h ⁢ ( x ^ ) ⁢ b ^ ⁢ ( x ^ ) ) M ^ 1 2 is a norm. Because all norms on finite-dimensional spaces are equivalent ∥ q h ^ ( b ^ ) ∥ M ^ ≳ ∥ q h ^ ∥ ν ′ , M ^ and we are working on a reference patch the constant does not depend on the cell sizes. Therefore

( ω h , q h ) M ≳ h x ⁢ h y ⁢ ∫ M ^ q ^ h ⁢ ( x ^ ) 2 ⁢ b ^ ⁢ ( x ^ ) ⁢ 𝑑 x ^ ≳ h x ⁢ h y ⁢ ∥ q ^ h ∥ ν ′ , M ^ 2 .

Now we estimate ∥ ω h ∥ M . We have | b ^ | ≤ 1 for all x ^ ∈ M ^ since it is a hat function on the reference element. Then

∥ ω h ∥ ν , M ν = ∫ M ^ b ^ ν ( x ^ ) q h ^ ν | det D F K | d x ^ ≤ 1 ∫ M ^ q h ^ ν | det D F K | d x ^ ≲ h x h y ∥ q ^ h ∥ ν , M ^ ν .

Again using the argument that all norms on finite-dimensional spaces are equivalent and we are working on a reference patch, we have that ∥ q h ^ ∥ ν , M ^ ∼ ∥ q h ^ ∥ ν ′ , M ^ , where the constant is independent of the cell sizes. Combining everything, we have the following: for all q h ∈ Y h ⁢ ( M ) there exists ω h ∈ X h 0 ⁢ ( M ) such that

( q h , ω h ) M ∥ ω h ∥ ν , M ⁢ ∥ q h ∥ ν ′ , M ≥ β loc ⁢ h x ⁢ h y ⁢ ∥ q ^ h ∥ ν ′ , M ^ 2 ( h x ⁢ h y ⁢ ∥ q ^ h ∥ ν ′ , M ^ ν ) 1 ν ⁢ ( h x ⁢ h y ⁢ ∥ q h ^ ∥ ν ′ , M ν ′ ) 1 ν ′ = β loc ⁢ h x ⁢ h y ( h x ⁢ h y ) 1 ν ⁢ ( h x ⁢ h y ) 1 ν ′ = β loc ,

where the constant β loc is independent of the aspect ratio. ∎

B Relations Between 𝐒 , 𝐃 and the Help Functions φ , 𝓕

Lemma 13 (Relations Between 𝐒 , 𝐃 , φ and ℱ ).

Given that S is of ( p , ϵ ) -structure, i.e. that Assumption 1 is fulfilled, the following relations are valid:

It holds uniformly for all 𝐏 , 𝐐 ∈ ℝ d × d that

( 𝐒 ⁢ ( 𝐏 ) - 𝐒 ⁢ ( 𝐐 ) ) : ( 𝐏 - 𝐐 ) ∼ ( ϵ + | 𝐏 sym | + | 𝐐 sym | ) p - 2 ⁢ | 𝐏 sym - 𝐐 sym | 2

∼ φ ϵ + | 𝐏 sym | ⁢ ( | 𝐏 sym - 𝐐 sym | )

∼ | 𝓕 ⁢ ( 𝐏 ) - 𝓕 ⁢ ( 𝐐 ) | 2 ,

| 𝐒 ⁢ ( 𝐏 ) - 𝐒 ⁢ ( 𝐐 ) | ∼ ( ϵ + | 𝐏 sym | + | 𝐐 sym | ) p - 2 ⁢ | 𝐏 sym - 𝐐 sym |

∼ φ ϵ + | 𝐏 sym | ′ ⁢ ( | 𝐏 sym - 𝐐 sym | ) .
For 𝐮 , 𝐯 ∈ 𝐖 1 , p the stress-tensor and help-function 𝓕 relates to the quasi-norm, or so-called “natural distance”, of [ 3 ] given by | 𝐯 | ( p , 𝐮 ) 2 := ∫ Ω ( ϵ + | 𝐃𝐮 | + | 𝐃𝐯 | ) p - 2 ⁢ | 𝐃𝐯 | 2 as

( 𝐒 ⁢ ( 𝐃𝐮 ) - 𝐒 ⁢ ( 𝐃𝐯 ) , 𝐃𝐮 - 𝐃𝐯 ) Ω ∼ ∥ 𝐅 ⁢ ( 𝐃𝐮 ) - 𝐅 ⁢ ( 𝐃𝐯 ) ∥ 2 2

∼ ∫ Ω φ ϵ + | 𝐃𝐮 | ⁢ ( | 𝐃𝐮 - 𝐃𝐯 | )

∼ | 𝐮 - 𝐯 | ( p , 𝐮 ) 2 .
For p ∈ ( 1 , 2 ] and for all 𝐯 , 𝐮 ∈ 𝐖 1 , p ,

∥ 𝐃 ⁢ ( 𝐯 - 𝐮 ) ∥ p 2 ≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐮 ) ∥ 2 2 ⁢ ∥ ϵ + | 𝐃𝐯 | + | 𝐃𝐮 | ∥ p 2 - p ,

∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐮 ) ∥ 2 2 ≲ ∥ 𝐃 ⁢ ( 𝐯 - 𝐮 ) ∥ p p ,

(B.1) ∥ 𝐒 ⁢ ( 𝐃𝐯 ) - 𝐒 ⁢ ( 𝐃𝐮 ) ∥ p ′ ≲ ∥ 𝓕 ⁢ ( 𝐃𝐯 ) - 𝓕 ⁢ ( 𝐃𝐮 ) ∥ 2 2 p ′ .

The constants depend on p, σ 0 and σ 1 . In particular, it does not depend on ϵ.

Proof.

See [14] and [19]. ∎

Lemma 14 (Properties of φ).

The following statements hold:

Due to the convexity of φ we have uniformly in t , s , a ∈ ℝ 0 + that

φ a ⁢ ( t + s ) ≤ φ a ⁢ ( s ) + φ a ⁢ ( t ) .
Young-type inequality: For all δ > 0 there exists c ⁢ ( δ ) > 0 that only depends on p and δ so that for all s , t , a ≥ 0 ,

s ⁢ φ a ′ ⁢ ( t ) + φ a ′ ⁢ ( s ) ⁢ t ≤ δ ⁢ φ a ⁢ ( s ) + c ⁢ ( δ ) ⁢ φ a ⁢ ( t ) .
Change of shift: for all 𝐏 , 𝐐 ∈ ℝ d × d , ϵ ≥ 0 we have

φ ϵ + | 𝐏 | ⁢ ( | 𝐏 - 𝐐 | ) ∼ φ ϵ + | 𝐐 | ⁢ ( | 𝐏 - 𝐐 | ) ,

where the constant depends only on p . The same holds if φ is replaced by φ ′ . Furthermore, it holds that for t ≥ 0 ,

φ ϵ + | 𝐏 | ⁢ ( t ) ∼ c ⁢ ( δ ) ⁢ φ ϵ + | 𝐐 | ⁢ ( t ) + δ ⁢ φ ϵ + | 𝐐 | ⁢ ( | 𝐏 - 𝐐 | ) ,

where δ > 0 and c ⁢ ( δ ) > 0 only depend on p and δ.

Proof.

These properties are shown in [14] and [19]. ∎

Further relations between 𝐒 , φ, and 𝓕 are described in [18].

C Approximation and Stability Properties of i h

From [2, Theorem 3.3] we have that for ω ∈ W 1 , ν ,

∥ ω - i h ⁢ ω ∥ 0 , ν ; M ≲ h x ⁢ ∥ ∂ x ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) ,

∥ ω - i h ⁢ ω ∥ 1 , ν ; M ≲ ∥ ω ∥ 1 , ν ; Λ ⁢ ( M ) ,

and for ω ∈ W 2 , ν ,

∥ ω - i h ⁢ ω ∥ 0 , ν ; M ≲ h x 2 ⁢ ∥ ∂ x 2 ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h y 2 ⁢ ∥ ∂ y 2 ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) + h x ⁢ h y ⁢ ∥ ∂ x ⁢ y ⁡ ω ∥ 0 , ν ; Λ ⁢ ( M ) ,

∥ ω - i h ⁢ ω ∥ 1 , ν ; M ≲ h x ⁢ ∥ ∂ x ⁡ ω ∥ 1 , ν ; Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ ω ∥ 1 , ν ; Λ ⁢ ( M ) .

Using ∥ i h ⁢ ω ∥ = ∥ ω - i h - ω ∥ ≤ ∥ ω - i h ⁢ ω ∥ + ∥ ω ∥ and the approximation properties, we get the stability properties

(C.1) ∥ i h ⁢ ω ∥ M ≲ ∥ ω ∥ Λ ⁢ ( M ) + h x ⁢ ∥ ∂ x ⁡ ω ∥ Λ ⁢ ( M ) + h y ⁢ ∥ ∂ y ⁡ ω ∥ Λ ⁢ ( M ) ,

(C.2) ∥ ∂ x ⁡ i h ⁢ ω ∥ ν ; M ≲ ∥ ω ∥ 1 , ν ; Λ ⁢ ( M ) ,

(C.3) ∥ ∂ y ⁡ i h ⁢ ω ∥ ν ; M ≲ ∥ ω ∥ 1 , ν ; Λ ⁢ ( M ) .

In [10] there is another useful stability estimate

(C.4) ∥ ∂ x ⁡ i h ⁢ ω ∥ ν ; M ≲ ∥ ∂ x ⁡ ω ∥ ν ; Λ ⁢ ( M ) + h y h x ⁢ ∥ ∂ y ⁡ ω ∥ ν ; Λ ⁢ ( M ) ,

(C.5) ∥ ∂ y ⁡ i h ⁢ ω ∥ ν ; M ≲ ∥ ∂ y ⁡ ω ∥ ν ; Λ ⁢ ( M ) .

These properties also extend to vector-valued functions 𝝎 ∈ 𝐖 0 1 , ν .

Acknowledgements

We would like to express our gratitude to the thorough reviewers who helped us improve the manuscript.

References

[1] J. Ahlkrona, N. Kirchner and P. Lötstedt, A numerical study of scaling relations for non-Newtonian thin-film flows with applications in ice sheet modelling, Quart. J. Mech. Appl. Math. 66 (2013), no. 4, 417–435. 10.1093/qjmam/hbt009Search in Google Scholar

[2] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., B. G. Teubner, Stuttgart, 1999. Search in Google Scholar

[3] J. W. Barrett and W. B. Liu, Quasi-norm error bounds for the finite element approximation of a non-Newtonian flow, Numer. Math. 68 (1994), no. 4, 437–456. 10.1007/s002110050071Search in Google Scholar

[4] R. Becker, An adaptive finite element method for the incompressible Navier–Stokes equations on time-dependent domains, PhD thesis, Heidelberg University, 1995. Search in Google Scholar

[5] R. Becker and M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (2001), no. 4, 173–199. 10.1007/s10092-001-8180-4Search in Google Scholar

[6] L. Belenki, L. C. Berselli, L. Diening and M. Růžička, On the finite element approximation of p-Stokes systems, SIAM J. Numer. Anal. 50 (2012), no. 2, 373–397. 10.1137/10080436XSearch in Google Scholar

[7] J. Blasco, An anisotropic GLS-stabilized finite element method for incompressible flow problems, Comput. Methods Appl. Mech. Engrg. 197 (2008), no. 45–48, 3712–3723. 10.1016/j.cma.2008.02.031Search in Google Scholar

[8] M. Braack, Anisotropic H 1 -stable projections on quadrilateral meshes, Numerical Mathematics and Advanced Applications, Springer, Berlin (2006), 495–503. 10.1007/978-3-540-34288-5_45Search in Google Scholar

[9] M. Braack, G. Lube and L. Röhe, Divergence preserving interpolation on anisotropic quadrilateral meshes, Comput. Methods Appl. Math. 12 (2012), no. 2, 123–138. 10.2478/cmam-2012-0016Search in Google Scholar

[10] M. Braack and T. Richter, Local projection stabilization for the Stokes system on anisotropic quadrilateral meshes, Numerical Mathematics and Advanced Applications, Springer, Berlin (2006), 770–778. 10.1007/978-3-540-34288-5_75Search in Google Scholar

[11] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Texts Appl. Math. 15, Springer, New York, 2008. 10.1007/978-0-387-75934-0Search in Google Scholar

[12] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Ser. Comput. Math. 15, Springer, New York, 1991. 10.1007/978-1-4612-3172-1Search in Google Scholar

[13] L. Diening, C. Kreuzer and E. Süli, Finite element approximation of steady flows of incompressible fluids with implicit power-law-like rheology, SIAM J. Numer. Anal. 51 (2013), no. 2, 984–1015. 10.1137/120873133Search in Google Scholar

[14] L. Diening and M. Růžička, Interpolation operators in Orlicz–Sobolev spaces, Numer. Math. 107 (2007), no. 1, 107–129. 10.1007/s00211-007-0079-9Search in Google Scholar

[15] L. Diening, M. Růžička and K. Schumacher, A decomposition technique for John domains, Ann. Acad. Sci. Fenn. Math. 35 (2010), no. 1, 87–114. 10.5186/aasfm.2010.3506Search in Google Scholar

[16] V. Girault and P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer Ser. Comput. Math. 5, Springer, Berlin, 1986. 10.1007/978-3-642-61623-5Search in Google Scholar

[17] C. Helanow and J. Ahlkrona, Stabilized equal low-order finite elements in ice sheet modeling—accuracy and robustness, Comput. Geosci. 22 (2018), no. 4, 951–974. 10.1007/s10596-017-9713-5Search in Google Scholar

[18] A. Hirn, Finite element approximation of problems in non-Newtonian fluid mechanics, PhD thesis, Heidelberg Univeryity and Charles University, 2011. Search in Google Scholar

[19] A. Hirn, Approximation of the p-Stokes equations with equal-order finite elements, J. Math. Fluid Mech. 15 (2013), no. 1, 65–88. 10.1007/s00021-012-0095-0Search in Google Scholar

[20] R. E. Johnson and R. M. McMeeking, Near-surface flow in glaciers obeying Glen’s law, Quart. J. Mech. Appl. Math. 37 (1984), no. 2, 273–291. 10.1093/qjmam/37.2.273Search in Google Scholar

[21] G. Matthies, P. Skrzypacz and L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, M2AN Math. Model. Numer. Anal. 41 (2007), no. 4, 713–742. 10.1051/m2an:2007038Search in Google Scholar

[22] S. Micheletti, S. Perotto and M. Picasso, Stabilized finite elements on anisotropic meshes: A priori error estimates for the advection-diffusion and the Stokes problems, SIAM J. Numer. Anal. 41 (2003), no. 3, 1131–1162. 10.1137/S0036142902403759Search in Google Scholar

[23] D. Sandri, Sur l’approximation numérique des écoulements quasi-newtoniens dont la viscosité suit la loi puissance ou la loi de Carreau, RAIRO Modél. Math. Anal. Numér. 27 (1993), no. 2, 131–155. 10.1051/m2an/1993270201311Search in Google Scholar

Received: 2018-10-08

Accepted: 2019-01-28

Published Online: 2019-02-15

Published in Print: 2020-01-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

Equal-Order Stabilized Finite Element Approximation of the p-Stokes Equations on Anisotropic Cartesian Meshes

Abstract

1 Introduction

2 Finite Element Formulation

2.1 Notation

2.2 Weak Formulation and Discretization

2.3 Stabilization

2.4 Anisotropic Meshes

Lemma 1.

Proof.

2.5 Some Useful Help-Functions and Properties of 𝐒

2.5.1 The Help-Functions φ, 𝓕 and Their Relation to the Stress Tensor 𝐒

Assumption 1 (On the Stress Tensor 𝐒 ).

2.5.2 The Help-Functions 𝓖 1 , 𝓖 2 and Their Relation to the Stabilization Term and φ *

Lemma 2.

Proof.

Lemma 3.

Proof.

Lemma 4.

Proof.

3 Error Analysis

3.1 Modified Interpolation Operator

3.1.1 Standard Properties in an Anisotropic p-Stokes Setting

Lemma 5.

Proof.

3.1.2 Approximation Property in Orlicz Spaces on Anisotropic Meshes

Lemma 6.

Orlicz-Stability.

Lemma 7.

Proof.

Lemma 8 (First-Order Anisotropic Orlicz-Stability).

Proof.

Orlicz-Approximability.

Lemma 9.

Proof.

Lemma 10 (Anisotropic Orlicz-Approximability).

Proof.

Corollary 1 (Orlicz-Continuity).

Proof.

Corollary 2.

Proof.

Orlicz-Continuity and Orlicz-Approximability in Terms of 𝓕

Proof of Lemma 6 .

3.2 Well-Posedness

Lemma 11 (Discrete inf-sup condition).

Proof.

Lemma 12 (Existence and Stability).

Proof.

3.3 A Priori Error Estimate

Theorem 1 (A Priori Error Estimates).

Proof.

Remark.

4 Numerical Experiments

4.1 Experiment 1 – Convergence

4.2 Experiment 2 – Robustness

5 Summary and Conclusion

A Proof of Lemma 1

Proof.

B Relations Between 𝐒 , 𝐃 and the Help Functions φ , 𝓕

Lemma 13 (Relations Between 𝐒 , 𝐃 , φ and ℱ ).

Proof.

Lemma 14 (Properties of φ).

Proof.

C Approximation and Stability Properties of i h

Acknowledgements

References

Journal and Issue

Articles in the same Issue