1. technology and computing
2. software

Institut für Numerische und Angewandte Mathematik
Nodal-based nite element methods with local projection
stabilization for linearized incompressible magnetohydrodynamics
Wacker, B., Arndt, D., Lube, G.
Nr. 3
Preprint-Serie des
Instituts für Numerische und Angewandte Mathematik
Lotzestr. 16-18
D - 37083 Göttingen
Article Submitted to ***Journal title***
Nodal-based finite element methods with local
projection stabilization for linearized
incompressible magnetohydrodynamics
Benjamin Wacker, Daniel Arndt, Gert Lube
Institute for Numerical and Applied Mathematics,
Georg-August University of Göttingen, Lotzestr. 16-18, D-37083 Göttingen,
b.wacker/d.arndt/[email protected]
Abstract
We consider the linearized resistive magnetohydrodynamics (MHD) model
in incompressible media and its numerical solution. Conforming nodal-based
finite element approximations of the MHD model where inf-sup stable as well
as equal-order finite elements for the velocity, pressure, the magnetic field and
the magnetic pseudo-pressure are investigated. As opposed to a residual-based
stabilization method by Badia et al. in J. Comp. Phys., Vol. 234, (2013), 399–
416, we consider a local projection stabilization for the numerical solution. A
detailed stability and error analysis for the arising discrete problem is given.
Some numerical experiments like Hartmann’s MHD problem and singular
solutions are examined.
1. Introduction
The time-dependent resistive incompressible MHD problem is given by
∂t u − ν∆u + (u · ∇) u + ∇p − (∇ × b) × (%b) = fu ,
%∂t b + λ∇ × ∇ × b + ∇r − ∇ × (u × (%b)) = fb ,
∇ · u = 0,
∇·b=0
(1.1)
(1.2)
in a bounded Lipschitz polyhedral domain Ω ⊂ Rd with d ∈ {2, 3}. For a given
density % and certain force terms fu and fb , we seek the velocity field u, the magnetic
field b, the pressure p and the magnetic pseudo-pressure r.
The standard approach to the numerical solution of problem (1.1)-(1.2) consists of
curl-conforming finite element methods (FEM), see [11] for a review. There is a
revival of nodal-based FEM for the Maxwell problem starting from the paper [9]
by Costabel/Dauge. Stabilization techniques of residual type based on nodal-based
1
FEM were considered by Codina/Badia in [2] and Bonito et al. in [7, 8]. Extensions of
nodal-based stabilized FEM to the resistive MHD model can be found in [3–5]. The
focus of the latter papers is on equal-order interpolation of all unknown variables.
Extending our approach in [19], the goal here is to analyze nodal-based FEM with
conforming inf-sup stable element pairs for the approximation of (u, p) and (b, r).
These discrete models are augmented by standard global stabilization techniques
for the divergence-free constraints and local projection stabilization (LPS) of the
other first order terms in the momentum and induction equations. For comparison
purposes, we consider the case of equal-order approximation of all unknowns as well.
The paper is organized as follows: After the introduction in Sec. 1, we introduce
the continuous resistive MHD model and its linearization in Sec. 2. Then, Sec. 3 is
concerned with the discretized and LPS stabilized MHD model. Stability estimates
for the linearized discrete MHD problem are given in Sec. 4. Error estimates for infsup stable and equal-order interpolaton are given in Sec. 5. In Sec. 6, we describe the
application to the time-dependent MHD model and present some standard numerical
examples. Finally, a summary and conclusions are given in Sec. 7.
We will denote by (·, ·)G the inner product in L2 (G) and by k · km,p,G with m ∈ N0
and p ≥ 1 the norm on the Sobolev space W m,p (G) for subdomains G ⊆ Ω. In case
of G = Ω we will omit index Ω.
2. Continuous resistive MHD model
Starting with a linearization of the MHD model (1.1)-(1.2), we derive the weak
form of the linearized stationary problem. We follow the time discretization and
linearization approach in [3]. Let θ ∈ [0, 1] and set
un+1 − un
,
∆t
bn+1 − bn
=
,
∆t
u˙ n+1 =
un+θ = θun+1 + (1 − θ) un ,
b˙ n+1
bn+θ = θbn+1 + (1 − θ) bn .
ˆ n time extrapolations to un+θ and bn+θ , we receive
ˆ n and b
Denoting by u
ˆ n = fu ,
u˙ n+1 − ν∆un+θ + (ˆ
un · ∇)un+θ + ∇pn+θ − ∇ × bn+θ × %b
(2.1)
∇ · un+θ = 0,
ˆ n ) = fb ,
%b˙ n+1 + λ∇ × ∇ × bn+θ + ∇rn+θ − ∇ × (un+θ × %b
(2.2)
∇ · bn+θ = 0.
(2.3)
(2.4)
By examining
theequilibrium point via u˙ n+1 = 0 and b˙ n+1 = 0, we linearize around
ˆ n . Setting u = un+θ , b = bn+θ , p = pn+θ and r = rn+θ , this yields
ˆ n , %b
(a, d) = u
−ν∆u + (a · ∇)u + ∇p − (∇ × b) × d = fu ,
λ∇ × (∇ × b) + ∇r − ∇ × (u × d) = fb ,
2
∇ · u = 0,
∇ · b = 0.
(2.5)
(2.6)
For clarity, we use homogeneous boundary conditions. We then introduce the spaces
h
1
id
V = v ∈ H (Ω)
: v = 0 on ∂Ω ,
C = {c ∈ H (curl ; Ω) : n × c = 0 on ∂Ω} ,
Q = L20 (Ω) ,
(2.7)
S = H01 (Ω)
(2.8)
and the bilinear and linear forms
AG (U, V) = ν (∇u, ∇v) + ha · ∇u, vi − (p, ∇ · v) − h(∇ × b) × d, vi
+ (∇ · u, q) − (b, ∇s)
(2.9)
+λ (∇ × b, ∇ × c) + (∇r, c) − h∇ × (u × d) , ci
FG (V) = hfu , vi + hfb , ci.
(2.10)
together with appropriate inner and dual products (·, ·) and h·, ·i. The variational
problem reads now: Find U := (u, b, p, r) ∈ V × C × Q × S such that
AG (U, V) = FG (V)
(2.11)
for all V := (v, c, q, s) ∈ V × C × Q × S. The well-posedness of problem (2.11) is
shown in [4, Thm. 1].
3. Local projection stabilization method
Here we introduce conforming FE subspaces
Vh × Qh × Ch × Sh ⊆ V × Q × C × S.
{Th }h is a family of exact shape-regular decompositions of Ω. For a simplex T ∈ Th
or a quadrilateral/hexahedron T in Rd , let Tˆ be the reference unit simplex or unit
cube (−1, 1)d . The bijective reference mapping FT : Tˆ → T is affine for simplices
and multi-linear for quadrilaterals/hexahedra. Pl resp. Ql with l ∈ N0 are the sets
of polynomials of degree ≤ l resp. of polynomials of degree ≤ l in each variable and
let
(
Pl (Tˆ) on simplices Tˆ
Rl (Tˆ) :=
Ql (Tˆ) on quadrilaterals/hexahedra Tˆ.
:
We define Yh,−l
= {vh ∈ L2 (Ω) : vh |T ◦ FT ∈ Rl (Tˆ) ∀T ∈ Th } and Yh,l := Yh,−l ∩
W 1,2 (Ω) and assume the following inverse and approximation properties.
(A.1): Let the FE space Yh,k satisfy the local inverse inequality
k∇vh k0,2,T ≤ Ch−1
T kvh k0,2,T
∀vh ∈ Yh,k , ∀T ∈ Th .
(A.2): There are interpolation operators j1 : W 1,2 (Ω) → Yh,k and j2 : W 1,2 (Ω) →
Yh,−k such that for all T ∈ Th , for all w ∈ W 1,2 (Ω) with 1 ≤ l ≤ k + 1:
kw − j1 wk0,2,T + hT k∇(w − j1 w)k0,2,T ≤ ChlT kwkl,2,T
3
and for all q ∈ W l,2 (Ω) with 1 ≤ l ≤ k:
kq − j2 qk0,2,T + hT k∇(q − j2 q)k0,2,T ≤ ChlT kqkl,2,T .
Our main goal is to apply inf-sup stable spaces Vh × Qh , Ch × Sh with
(A.3): (Discrete inf-sup condition)
∃β 6= β(h) > 0 s.t. :
(∇ · vh , qh )
≥ βkqh k0,2 ∀qh ∈ Qh .
vh ∈Vh \{0} k∇vh k0,2
sup
(3.1)
Then we write for convenience Vh = Rk instead of Vh = [Yh,k ]d ∩ V and Qh = Rk−1
instead of Qh = Yh,k−1 ∩ Q, k ≥ 2. A similar notation will be used for Ch and Sh .
For comparison purposes, we consider as well the case of equal-order interpolation
with Vh and Ch as above (but k ≥ 1) and Qh = Yh,k ∩ Q and Sh = Yh,k ∩ S, k ≥ 1.
For the FEM stabilization, we introduce a macro grid Mh = {M } with Mh =
Th (one-level method) or Mh = T2h (two-level method). In the latter case, the
decomposition Th is derived from Mh by barycentric refinement of d-simplices or
regular (dyadic) refinement of quadrilaterals or hexahedra.
On M ∈ Mh we define local FE spaces
s
DM
:= {v ∈ [L2 (M )]d : v|M ∈ Rs },
s ∈ {0, . . . , k}.
(3.2)
s
s
. Moreover,
: [L2 (M )]d → DM
The local orthogonal L2 -projection are given by πM
s
s
we denote by κM := id − πM the so-called fluctuation operator.
(A.4): The fluctuation operators κsM provide the approximation property (depending
s
on DM
and s ∈ {0, . . . , k}):
kκSM wk0,2,M ≤ ChlM kwkl,2,M ,
∀w ∈ [W 1,2 (Ω)]d , M ∈ Mh , l = 0, . . . , s.
(3.3)
S
.
A sufficient condition for (A.4) is Ps−1 ⊂ DM
Let Uh := (uh , bh , ph , rh ) and Vh := (vh , ch , qh , sh ) with Uh , Vh ∈ Vh ×Ch ×Qh ×Sh .
Then the local projection stabilization (LPS) term in the inf-sup stable case reads
ISS
Slps
(Uh , Vh ) =
X
{τ1 (κsM ((aM · ∇)uh ), κsM ((aM · ∇)vh ))M
M
+ τ2 (∇ · uh , ∇ · vh )M
+ τ3 (κsM ((∇ × bh ) × dM ), κsM ((∇ × ch ) × dM ))M
+ τ4 (κsM (∇ × (uh × dM )), κsM (∇ × (vh × dM )))M
+ τ5 (∇rh , ∇sh )M + τ6 (∇ · bh , ∇ · ch )M
(3.4)
whereas the equal-order LPS terms are given by
EO
ISS
Slps
(Uh , Vh ) = Slps
(Uh , Vh ) +
X
τ7 (κsM (∇ph ), κsM (∇qh ))M .
(3.5)
M
Here we suppose that aM and dM are piecewise constant approximations of the fields
4
a and d on M ∈ Mh with |aM | ≤ kak0,∞,M and |dM | ≤ kdk0,∞,M . The parameters
τj for j ∈ {1, 2, 3, 4, 5, 6, 7} are piecewise-constants on the cells M of Mh . Then the
stabilized problem reads: Find Uh ∈ Vh × Ch × Qh × Sh such that
Astab (Uh , Vh ) := AG (Uh , Vh ) + Slps (Uh , Vh ) = FG (Vh )
(3.6)
holds for all Vh ∈ Vh × Ch × Qh × Sh .
4. Stability analysis for the LPS method
4.1. Some auxiliary results
We first examine some properties of the LPS terms Slps in (3.6).
Lemma 4.1. For U := (u, b, p, r) , V := (v, c, q, s) ∈ V × C × Q × S we obtain
Slps (U, U) ≥ 0,
(4.1)
1
2
1
2
|Slps (U, V)| ≤ (Slps (U, U)) (Slps (V, V)) .
(4.2)
Proof: The first property follows from the definition (3.4) or (3.5).
The Cauchy-Schwarz inequality for the L2 -scalar product implies
EO
Slps (U, V)
1 EO
≤ Slps
(U, U)
2
1
EO
Slps
(V, V)
2
for (3.5) and similarly for (3.4) by setting τ7 = 0.
2
For the solutions U ∈ V × C × Q × S resp. Uh ∈ Vh × Ch × Qh × Sh of the continuous
resp. discrete problems, we receive the approximate Galerkin orthogonality.
Lemma 4.2. Let U and Uh be the given solutions as described above. For an arbitrary test function Vh ∈ Vh × Ch × Qh × Sh ⊂ V × C × Q × S holds
(AG + Slps ) (U − Uh , Vh ) = Slps (U, Vh ).
(4.3)
Proof: Subtracting problems (2.11) and (3.6) we obtain
AG (U − Uh , Vh ) − Slps (Uh , Vh ) = 0.
Adding Slps (U, Vh ), we find the approximate Galerkin orthogonality (4.3).
2
4.2. Stability estimates
We now conclude the unique existence of the stationary problem (3.6) by a detailed
stability analysis. Let V := (v, c, q, s) ∈ V × C × Q × S. Assuming ∇ · a = 0 in the
5
linearization at the equilibrium point (or otherwise introduce the skew-symmetric
form of the nonlinear term in the momentum equation), integration by parts gives
1
ha · ∇v, vi = − h(∇ · a) v, vi = 0.
2
(4.4)
h(∇ × c) × d, vi = −h∇ × (v × d) , ci.
(4.5)
Moreover, we have
Hence symmetric testing yields
AG (V, V) = ν (∇v, ∇v) − (q, ∇ · v) + ha · ∇v, vi − h(∇ × c) × d, vi
+ (∇ · v, q) − (c, ∇s)
+ λ (∇ × c, ∇ × c) − h∇ × (v × d) , ci + (∇s, c)
= νk∇vk20,2 + λk∇ × ck20,2 .
We therefore define the semi-norms
kVkG :=
νk∇vk20,2 + λk∇ × ck20,2
1
2
,
1
2
kVkStab := Slps (V, V) , .
(4.6)
(4.7)
Our next task is to show stability of the semi-discretized stabilized problem for each
solution variable separately. Let Uh ∈ Vh × Ch × Qh × Sh be the solution of the
LPS-problem (3.6). Let a, d ∈ [L∞ (Ω)]d in what follows.
(i) At first, we want to give an estimate for the magnetic field which depends on the
velocity field at first. As proposed in [4], we define the norm
√
√
λ
kck0,2 + λk∇ × ck0,2
kckC :=
(4.8)
L0
on the space C of the magnetic fields and the norm
1
L0
kskS := √ ksk0,2 + √ k∇sk0,2
λ
λ
(4.9)
on the space S of the magnetic pseudo-pressures. Here, L0 defines a length scale of
the problem which has to be chosen. We define the Maxwell bilinear form by
CM axw ((b, r) , (c, s)) := λ (∇ × b, ∇ × c) + (∇r, c) − (b, ∇s)
= λ (∇ × b, ∇ × c) + (∇r, c) + (∇ · b, s)
= AG ((0, b, 0, r) , (0, c, 0, s))
(4.10)
and obtain the following inf-sup condition.
Lemma 4.3. For the continuous Maxwell problem, the Babuska-Brezzi-condition
inf
sup
(b,r)∈C×S (c,s)∈C×S
CM axw ((b, r) , (c, s))
≥ βm
(kbkC + krkS ) (kckC + kskS )
holds with a constant βm > 0.
6
(4.11)
Proof: The proof is given in [2, Thm. 2.1] which uses [10, Thm. 2.34] in combination
with [10, Propos. 2.36].
2
By Lemma 4.3 we can show stability of the discretized magnetic field.
Lemma 4.4. The discrete magnetic field solution bh ∈ Ch fulfills the estimation
√ !
!
λ
1
kbh kC ≤
kUh kG + max C √
kUh kStab .
(4.12)
M
βm
L0 τ6
Proof: Let (bh , 0) ∈ Ch × Sh ⊂ C × S. By (4.11) there exists a unique (c, s) ∈ C × S
with kckC + kskS = 1, hence
βm kbh kC ≤ CM axw ((bh , 0) , (c, s))
= AG ((0, bh , 0, 0) , (0, c, 0, s))
= λ (∇ × bh , ∇ × c) + (∇ · bh , s)
= CM axw ((bh , 0) , (c, 0)) + CM axw ((bh , 0) , (0, s)) .
|
{z
}
=:I
|
{z
}
=:II
For the first term, we get the estimate
√
√
√
|I| ≤ λk∇ × bh k0,2 λk∇ × ck0,2 ≤ λk∇ × bh k0,2 ≤ kUh kG .
|
{z
≤1
(4.13)
(4.14)
}
We introduce the Scott-Zhang interpolation operator ΠSZ with its stability and
approximation properties to obtain for the second term
|II| = |(∇ · bh , s − ΠSZ (s) + ΠSZ (s))|
!1
2
≤
X
τ6 k∇ ·
bh k20,2,M
X
M
M
1
ks − ΠSZ (s)k20,2,M
τ6
!1
1
+
τ6 k∇ · bh k20,2,M
kΠSZ (s)k20,2,M
τ
6
M
M
! √
λ L0 ksk0,2
hM + C
≤ kUh kStab · max
·
· √
√
M
τ6
L0
λ
2
X
!1
2
!1
2
X
|
√ !
C λ
≤ max
kUh kStab .
√
M
L0 · τ6
{z
≤1
}
Putting (4.14) and (4.15) into (4.13) leads to the assertion
√ !
!
1
C λ
kbh kC ≤
kUh kG + max
kUh kStab .
√
M
βm
L0 τ6
7
(4.15)
2
(4.16)
Squaring (4.12) and setting A := 1 + maxM
√ 2
C √λ
,
L0 τ6
Young’s inequality leads to
2A 2
2
+
kU
k
kU
k
h
h
Stab
G
2
βm
2A
≤ 2 (hfu , uh i + hfb , bh i)
βm
2A
≤ 2 (kfu k0,2 kuh k0,2 + kfb k0,2 kbh k0,2 )
βm
!
2
ν
AL
A L20 kfu k20,2
1
0
+ 2 kuh k20,2 +
≤ 2
kfb k20,2 + kbh k2C ,
2
βm
ν
L0
λβm
2
kbh k2C ≤
(4.17)
hence
kbh k2C
2A
≤ 2
βm
L20 kfu k20,2
ν
AL20
+ 2 kuh k20,2 +
kfb k20,2 .
2
ν
L0
λβm
!
(4.18)
(ii) Now, we want to give an estimate for the fluid-velocity field. We use
√
√
ν
kvkV :=
kvk0,2 + νk∇vk0,2
L0
as a norm on the velocity space V .
Lemma 4.5. The discrete fluid-velocity solution uh ∈ Vh fulfills
kuh k2V
1 L20
64A2 B 2 1
kfu k20,2 +
+
≤ 8B +
2
2 ν
βm
2
2
where the constant A is defined as in (4.17) and B := 1 +
!
L20
kfb k20,2
λ
CP
L0
2
.
Proof: Testing symmetrically Vh = Uh ∈ Vh × Ch × Qh × Sh , we get
(AG + Slps ) (Uh , Uh ) = kUh k2G + kUh k2Stab .
Poincare’s inequality, (4.18) and the definition of B imply
kuh k2V
≤ νBk∇uh k20,2
n
≤ B kUh k2G + kUh k2Stab
o
= B {hfu , uh i + hfb , bh i}
≤ B {kfu k0,2 kuh k0,2 + kfb k0,2 kbh k0,2 }
1
BL0
4B 2 L20
≤
kfu k20,2 + kuh k2V + √ kfb k0,2 kbh kC
ν
4
λ
√
2 2
4B L0
1
4 2ABL0
βm
√
=
kfu k20,2 + kuh k2V +
kfb k0,2 √ kbh kC
ν
4
4 2A
βm λ
8
(4.19)
2
1
32A2 B 2 L20
4B 2 L20
βm
2
kfu k20,2 + kuh k2V +
+
kf
k
kbh k2C
b
0,2
2
2
ν
4
λβm
32A
2 2
4B L0
1
kfu k20,2 + kuh k2V
≤
ν
2
1 L20
32A2 B 2 L20
1 L20
2
2
kf
k
+
kfb k20,2 .
+
kf
k
+
u 0,2
b 0,2
2
λβm
4 ν
4 λ
≤
This leads to
1 L20
64A2 B 2 1
kuh k2V ≤ 8B 2 +
kfu k20,2 +
+
2
2 ν
βm
2
!
L20
kfb k20,2 .
λ
(4.20)
This, together with (4.18), shows the stability, hence uniqueness and existence of
the discrete magnetic and velocity fields.
2
(iii) Now we derive an estimate for the kinematic pressure in terms of the velocity
field.
Lemma 4.6. There exists a unique discrete pressure solution ph ∈ Qh for (3.6)
such that
kph k0,2
√
!
Cp kak0,∞ Cp kdk0,∞
√
√
kfu k0,2 +
ν+
+
kUh kG
ν
λ
)!
!
(
q
√
Cκ hM
kUh kStab
+ max (CA τ1 kak0,∞,M ) + max CA τ2 d + max √
M
M
M
τ7
(4.21)
1
≤
βf
where constant CA stems from the assumptions on the interpolator juh : V → Vh .
Proof: For inf-sup stable spaces Vh × Qh , the discrete Babuska-Brezzi-condition for
the fluid part implies for all ph ∈ Qh the existence of a discrete vh ∈ Vh such that
∇ · vh = −ph ,
|vh |1,2 ≤ CAT H kph k0,2 ,
CAT H :=
1
.
βf
(4.22)
In case of equal-order interpolation, the continuous Babuska-Brezzi-condition for
the fluid part ensures for all ph ∈ Qh the unique existence of v ∈ V such that
∇ · v = −ph ,
|v|1,2 ≤
1
kph k0,2 .
βf
Then we set vh := juh v. The assumptions on the interpolation operator imply
|vh |1,2 ≤ |v|1,2 + C · |v|1,2 ≤
1
· (1 + C) · kph k0,2 =: CAEO · kph k0,2 .
βf
Considering
(AG + Slps ) (Uh , (vh , 0, 0, 0)) = hfu , vh i,
9
(4.23)
we end up with
kph k20,2 ≤ Cp kfu k0,2 |vh |1,2 + νk∇uh k0,2 |vh |1,2 + |Slps ((uh , 0, 0, 0) , (vh , 0, 0, 0))|
+ |(a · ∇uh , vh )| + |((∇ × bh ) × d, vh )| + |(ph , ∇ · (vh − v))| ,
(4.24)
where the last term disappears for inf-sup stable spaces Vh × Qh . We obtain
− (a · ∇uh , vh ) = (a · ∇vh , uh )
≤ Cp kak0,∞ |uh |1,2 |vh |1,2 ,
((∇ × bh ) × d, vh ) ≤ Cp kdk0,∞ k∇ × bh k0,2 |vh |1,2 ,
k(vh , 0, 0, 0)k2Stab =
(4.25)
(4.26)
Xn
τ1 kκ (a · ∇vh )k20,2,M + τ2 k∇ · vh k20,2,M
o
(4.27)
M
≤
max τ1 kak20,∞,M
M
|vh |21,2 + max (τ2 d) |vh |21,2 .(4.28)
M
Integrating (ph , ∇ · (vh − v)) by parts for (continuous) pressure ph , we see
(∇ph , vh − v) ≤
X
M
!1
Cκ2 h2M
τ7
2
· kUh kStab |vh |1,2
(
Cκ · hM
≤ CAEO max
√
M
τ7
(4.29)
)
kUh kStab · kph k0,2 ,
Summarizing (4.24)-(4.29) and dividing by kph k0,2 , we obtain
kph k0,2 ≤ CA kfu k0,2 +
(
+CA
√
!
Cp kak0,∞ Cp kdk0,∞
√
√
ν+
+
kUh kG
ν
λ
!
(
q
√
Cκ hM
max ( τ1 kak0,∞,M ) + max τ2 d + max √
M
M
M
τ7
))
kUh kStab .
Let us remark that the last term in the bracket disappears for inf-sup stable spaces
Vh × Qh . Moreover, the terms kUh kG and kUh kStab are bounded due to (i) and (ii).
This shows existence and uniqueness of the discrete fluid pressure.
2
(iv) Finally, we prove that the discrete magnetic pseudo-pressure is unique.
Lemma 4.7. There exists a unique discrete pseudo-pressure solution rh ∈ Sh for
(3.6) such that
(L0 + CP )L0
krh kS ≤
kbh kC .
(4.30)
λ min τ5
M
Proof: We see that
0 = (AG + Slps ) (Uh , (0, 0, 0, rh )) = − (bh , ∇rh ) +
X
M
10
τ5 k∇rh k20,2,M ,
hence
min τ5 k∇rh k20,2 ≤
M
X
τ5 k∇rh k20,2,M = (bh , ∇rh ) ≤ kbh k0,2 k∇rh k0,2 .
M
By using (4.8), we conclude that
k∇rh k0,2 ≤
1
L0
kbh k0,2 ≤ √
kbh kC .
min τ5
λ
min
τ
5
M
(4.31)
M
This and (4.9) imply the assertion and thus the unique existence of rh .
2
We summarize the results of steps (i)-(iv) in the following theorem.
Theorem 4.8. Let a, d ∈ L∞ (Ω), ∇ · a = 0, fu ∈ [L2 (Ω)]d , fb ∈ [L2 (Ω)]d . Then
there exists a unique discrete solution Uh ∈ Vh × Ch × Qh × Sh which fulfills (3.6).
Proof: The stability results of the Lemmata 4.4-4.7 imply uniqueness of (uh , bh , ph , rh )
and therefore existence of the discrete solution.
2
5. Quasi-optimal error estimates
Here our goal is to derive quasi-optimal error estimates for the LPS-method (3.6).
5.1. Preliminary results
Let U ∈ V × C × Q × S resp. Uh ∈ Vh × Ch × Qh × Sh be the solutions of the
continuous resp. of the discrete problems. Lemma 4.2 implies for the consistency
error that
(AG + Slps ) (U − Uh , Vh ) = Slps (U, Vh )
(5.1)
for all Vh ∈ Vh × Ch × Qh × Sh . Let J := ju , jb , j p , j r be appropriate interpolation
operators. Then we can decompose the error as follows
U − Uh = (U − JU) + (JU − Uh ) = E − Eh
(5.2)
with E := (εu , εb , εp , εr ) and Eh := (eu , eb , ep , er ). By setting Vh = Eh in equation
(5.1) and using (5.2), we get
(AG + Slps ) (U − Uh , Eh ) = (AG + Slps ) (E − Eh , Eh ) = Slps (U, Eh ) ,
which leads to
AG (Eh , Eh ) + Slps (Eh , Eh ) = −Slps (U, Eh ) − AG (E, Eh ) + Slps (E, Eh ) ,
hence
kEh k2G + kEh k2Stab ≤
Slps (U, Eh ) −AG (E, Eh ) − Slps (E, Eh ) |
{z
} |
{z
}|
{z
}
=:II
=:I
follows. We start with the following estimate.
11
=:−III
(5.3)
Lemma 5.1. With the above introduced notation, we get
kEh k2G + kEh k2Stab ≤ (kEkStab + kUkStab ) kEh kStab + kEkG kEh kG + |IV |
(5.4)
with
IV
= (a · ∇εu , eu ) − ((∇ × εb ) × d, eu ) − (∇ × (εu × d) , eb )
− (εp , ∇ · eu ) + (∇ · εu , ep ) + (∇εr , ep ) − (εb , ∇er ) .
(5.5)
Proof: Lemma 4.1 provides
1
1
|I| ≤ (Slps (U, U)) 2 (Slps (Eh , Eh )) 2 ≤ kUkStab kEh kStab ,
(5.6)
|III| = |Slps (E, Eh )| ≤ kEkStab kEh kStab .
(5.7)
and similarly
We also receive the estimate
−II ≤ kEkG kEh kG + (a · ∇εu , eu ) − ((∇ × εb ) × d, eu ) − (∇ × (εu × d) , eb )
− (εp , ∇ · eu ) + (∇ · εu , ep ) + (∇εr , ep ) − (εb , ∇er ) .
(5.8)
Using the definition of IV in (5.5), we can summarize the estimation as
kEh k2G + kEh k2Stab ≤ (kEkStab + kUkStab ) kEh kStab + kEkG kEh kG + |IV | .
(5.9)
Now we have to estimate the non-symmetric terms in IV . Integration by parts and
Cauchy-Schwarz inequality gives:
(a · ∇εu , eu ) = − a · ∇eu , εu ≤
X
M
− εp , ∇ · e u
≤
X
M
− (εb , ∇er ) ≤
X
M
1
kεp k20,2,M
τ2
1
1
kεb k20,2,M
τ5
1
2
2
1
kak20,∞,M
2
2
kεu k0,2,M kEh kG ,
ν
kEh kStab ,
kEh kStab ,
− (∇ × (εu × d) , eb ) = (εu , (∇ × eb ) × d) ≤
X
M
− (ep , ∇ · εu ) = (∇ep , εu ) ≤
X
M
1
kdk20,∞,M
2
2
kεu k0,2,M kEh kG ,
λ
1
kεu k20,2,M
τ7
!1
2
kEh kStab
h
In the inf-sup stable case, the latter term − (ep , ∇ · εu ) vanishes provided ju ∈ Vdiv
.
r
The term (∇εr , eb ) vanishes since r = j r ≡ 0.
h
id
Let d ∈ W1,∞ (Ω) . By formula ∇ × (e × f ) = f · ∇e − f (∇ · e) − e · ∇f + e (∇ · f ),
the inequalities of Cauchy, Schwarz and Poincare, it follows
−((∇ × εb ) × d, eu ) =
X
(εb , ∇ × (eu × d))M
M
12
(5.10)
√
1
(1 + d)2 (kdk0,∞,M + k∇dk0,∞,M )2 keb k20,2,M kEh kG .
ν
!
X
≤
M
We then summarize the above equations. Using Young’s inequality, we obtain
kEh k2G + kEh k2Stab ≤ S12 + S22 ,
with
S1
1
X1
2
1
:= kEkG +
kak20,∞,M kεu k20,2,M +
kdk20,∞,M kεu k20,2,M
ν
λ
M
M
1
X
2 2
√
2
2
−1
+
ν
1 + d kdk0,∞,M + k∇dk0,∞,M kεb k0,2,M ,
X
1
2
M
S2 := kEkStab + kUkStab +
X
M
X
+
M
1
kεu k20,2,M
τ7
1
kεp k20,2,M
τ2
1
2
+
X
M
1
kεb k20,2,M
τ5
1
2
!1
2
.
5.2. Error estimates for smooth solutions without LPS-compatibility
The approximation properties of the FE spaces, see [13], and the local L2 -projector
yield for U = (u, b, p, r) ∈ [H k+1 (Ω)]d × [H k+1 (Ω)]d × H k+1 (Ω) × H k+1 (Ω) that in
case of equal-order interpolation of all unknowns
S12
≤ C
X
h2k
M
M
kdk20,∞,M h2M 2
kak20,∞,M h2M +λ
|u|k+1,2,M
ν 1+
ν2
λ2
2
h2M |b|2k+1,2,M ,
(5.11)
kdk0,∞,M + k∇dk0,∞,M
ν
X
h2 2(s−k)
2k
τ1 |aM |2 + τ4 |dM |2 |u|2s+1,2,M
≤ C
hM τ2 + M |u|2k+1,2,M + hM
τ7
M
+ λ+
S22
2(s−k)
+hM
τ7 |p|2s+1,2,M +
h2M 2
2(s−k)
|p|k+1,2,M + hM
τ3 |dM |2 |b|2s+1,2,M
τ2
(5.12)
h2M 2
+ τ6 +
|b|k+1,2,M .
τ5
In case of inf-sup stable elements with U ∈ [H k+1 (Ω)]d ×[H k+1 (Ω)]d ×H k (Ω)×H k (Ω),
we obtain
S22
≤ C
X
h2k
M
2(s−k)
τ2 |u|2k+1,2,M + hM
τ1 |aM |2 + τ4 |dM |2 |u|2s+1,2,M
(5.13)
M
1 2
h2M 2
2(s−k)
2
2
+ |p|k,2,M + hM
τ3 |dM | |b|s+1,2,M + τ6 +
|b|k+1,2,M .
τ2
τ5
Denote the local fluid and magnetic Reynolds numbers by
Ref,M := kak0,∞,M hM /ν,
Rem,M := kdk0,∞,M hM /λ.
13
respectively. We will call an error estimate to be semi-robust of order m if the
coefficients multiplying the corresponding Sobolev norms of the solutions are of
order hm uniformly w.r.t. the problem data.
5.2.1. Inf-sup stable case
An equilibration of the terms of S12 in (5.11) according to
kak20,∞,M h2M kdk20,∞,M h2M
+
λ
∼ 1
ν2
λ2
2
h2 ∼ 1
λ + M kdk0,∞,M + k∇dk0,∞,M
ν
ν 1+
leads to the following (mild) restrictions on the local mesh width hM :
√
√
√
νRef,M ≤ C,
λRem,M ≤ C, hM (kdk0,∞,M + k∇dk0,∞,M ) ≤ C ν.
(5.14)
As we are in the case of inf-sup stable elements, it holds τ7 = 0. Equilibration of
terms in S22 in (5.13) and comparison to S12 yields
2(k−s)
τ1 |aM |2 + τ4 |dM |2 ∼ hM
2(k−s)
, τ2 ∼ 1, τ3 |dM |2 ∼ hM
, τ6 +
h2M
∼ 1.
τ5
(5.15)
This leads to the following bounds on the stabilization parameters:
2(k−s)
0 ≤ τ1 ≤ ChM
2(k−s)
/|aM |2 , τ2 ∼ 1, 0 ≤ τ3 , τ4 ≤ ChM
/|dM |2 .
(5.16)
Moreover, (5.15) suggests the balance τ5 τ6 ∼ h2M . Assuming L0 = CP , Lemma 4.7
delivers the reasonable choice τ5 ≥ CL20 λ−1 . Then we recover the proposal by [3], [4]
with
τ5 ∼ L20 /λ,
τ6 ∼ h2M λ/L20 .
(5.17)
Theorem 5.2. Assume that the solution (u, b, p) of (2.11) belongs to [H k+1 (Ω)]d ×
[H k+1 (Ω)]d × H k (Ω) and that ju u ∈ Vhdiv . Further, let the LPS parameters be chosen
according to condition (5.16)-(5.17) and that the local mesh width hM is chosen such
that (5.14) is valid for the inf-sup stable case. Then we obtain (using r ≡ 0)
kUh − JUk2G + kUh − JUk2Stab ≤ C
X
2
2
2
h2k
M |u|k+1,2,M + |b|k+1,2,M + |p|k,2,M .
M
Remark 5.3. To recover the full norms of the errors at hand, we have to adapt
the proofs of the stability analysis. Starting with the velocity field error u − uh ,
Poincare’s inequality implies
CP
ku − uh kV ≤ 1 +
L0
14
kU − Uh kG .
(5.18)
Proceeding as in Lemma 4.4, we receive for the error of the magnetic field:
1 √
kb − bh kC ≤
λk∇ × (b − bh )k0,2
βm
( √
)
!1
2
X
Cλ
1
+
max
τ6 k∇ · (b − bh )k20,2,M
√
βm M
L0 τ6
M
( √
)
1
1
Cλ
≤
kU − Uh kG +
max
kU − Uh kStab .
√
βm
βm M
L0 τ6
(5.19)
Similarly to Lemma 4.6, we get for the kinematic pressure
kp − ph k0,2
√
!
Cp kak∞ Cp kdk∞
≤ CA
+ √
kU − Uh kG
ν+ √
ν
λ
q
√
+ max ( τ1 kak0,∞,M ) + τ2 d kU − Uh kStab .
(5.20)
M
Using the interpolation properties, one also gets the desired results for the full
norms recovered from the Galerkin and the stabilized norm. This is also valid in
the following cases.
2
5.2.2. Equal-order case
(i) Error estimates of order k:
A first variant is to derive error estimates of order k as in the previous case. Then
the same (mild) mesh width restrictions follow from equilibration of terms in S12 :
√
√
√
νRef,M ≤ C,
λRem,M ≤ C, hM (kdk0,∞,M + k∇dk0,∞,M ) ≤ C ν. (5.21)
Equilibration of terms in S22 in (5.12) yields
h2
h2
h2M
2(k−s)
∼ 1, M ∼ 1, τ6 + M ∼ 1, τ1 |aM |2 + τ4 |dM |2 + τ3 |dM |2 + τ7 ∼ hM
.
τ7
τ2
τ5
(5.22)
This leads to the following bounds on the stabilization parameters:
τ2 +
2(k−s)
0 ≤ τ1 ≤ ChM
2(k−s)
/|aM |2 , 0 ≤ τ3 , τ4 ≤ ChM
/|dM |2 ,
2(k−s)
Ch2M ≤ τ7 ≤ ChM
Ch2M ≤ τ2 ≤ C,
(5.23)
.
Moreover, (5.22) suggests the balance τ5 τ6 ∼ h2M . Arguing as in the previous case,
we recover the proposal by [3], [4] with
τ6 ∼ h2M λ/L20 .
τ5 ∼ L20 /λ,
(5.24)
Theorem 5.4. Assume that the solution (u, b, p) of (2.11) belongs to [H k+1 (Ω)]d ×
[H k+1 (Ω)]d × H k+1 (Ω) and that ju u ∈ Vhdiv . Further, let the LPS parameters be
chosen according to condition (5.23)-(5.24) and that the local mesh width hM is as
in (5.25). Then, for the equal-order approximations, we obtain (using r ≡ 0)
kUh −
JUk2G
+ kUh −
JUk2Stab
≤C
X
h2k
M
M
15
|u|2k+1,2,M
+
|b|2k+1,2,M
+
|p|2k+1,2,M
.
(ii) Error estimates of order k + 12 :
Another variant is, in case of ν, λ ≤ hM , to derive error estimates of order k + 1/2.
Equilibration of terms in S12 leads to the standard (strong) mesh width restrictions:
Ref,M ≤ C, Rem,M ≤ C, hM (kdk0,∞,M + k∇dk0,∞,M ) ≤ Cν.
(5.25)
Equilibrating terms of S22 in (5.12), i.e. of
h2M
2(s−k)
, hM
τ1 |aM |2 + τ4 |dM |2 ∼ hM |aM | + ν,
τ7
h2
h2
2(s−k)
2(s−k)
hM
τ3 |dM |2 , τ6 + M , M , hM
τ7 ∼ hM |dM | + λ.
τ5
τ2
τ2 +
leads to
2(k−s)
C(hM |aM | + ν)hM
0 ≤ τ1 ≤
|aM |2
2(k−s)
C(hM |dM | + λ)hM
0 ≤ τ3 ≤
|dM |2
,
,
2(k−s)
0 ≤ τ4 ≤
Ch2M
C(hM |aM | + ν)hM
,
≤ τ2 ≤ C (hM |aM | + ν) ,
|dM |2
hM |dM | + λ
Ch2M
2(k−s)
≤ τ7 ≤ C(hM |d|M + λ)hM
.
hM |aM | + ν
(5.26)
Regarding the parameters τ5 and τ6 we repeat the argument of the previous cases
to obtain (see also [3],[4])
τ5 ∼ L20 /λ,
τ6 ∼ h2M λ/L20 .
(5.27)
Theorem 5.5. Assume that the solution (u, b, p) of (2.11) belongs to [H k+1 (Ω)]d ×
[H k+1 (Ω)]d × H k+1 (Ω) and that ju u ∈ Vhdiv . Further, let the LPS parameters be
chosen according to condition (5.26)-(5.27) and that the local mesh width hM is as
in (5.25). Then, for the equal-order approximations, we obtain (using r ≡ 0)
kUh −JUk2G +kUh −JUk2Stab
≤C
X
h2k
M (ν+λ+hM )
|u|2k+1,2,M +|b|2k+1,2,M +|p|2k+1,2,M
M
5.3. Error estimates for smooth solutions with LPS-compatibility
The restrictions (5.14) on the mesh width are not convincing. For simplicity, we
assume here elementwise constant fields a|M = aM and d|M = bM . This is a technical assumption for the stationary case in order to control certain L2 -terms. In the
time-dependent case one can omit this assumption as the L2 -terms can be controlled
via application of the Gronwall lemma, see [1] for the fluid part.
Let us assume the following LPS-compatibility condition between projection spaces
s
DM
and FE space Yh,k . Set Yh,k (M ) := {wh ∈ [Yh,k ]d : wh = 0 in Ω \ M }.
(A.5): There exist constants γ > 0 such that
(wh , vh )M
≥γ
vh ∈DM w ∈Y (M ) kvh k0,2,M kwh k0,2,M
h
h,k
inf s
sup
16
∀h > 0, ∀M ∈ Mh .
(5.28)
.
Following Thm. 2.2 in [17], one can show the existence of interpolation operators
ju : V → Vh , jb : C → Ch , j p : Q → Qh satisfying the orthogonality conditions
s
(v − ju v, ζh ) = 0 ∀v ∈ V and ∀ζh ∈ DM
,
(5.29)
s
= 0 ∀c ∈ C and ∀ηh ∈ DM
,
(5.30)
s
.
DM
(5.31)
c − jb c, ηh
p
(p − j p, ρh ) = 0 ∀p ∈ Q and ∀ρh ∈
Moreover, the interpolation properties of (A.2) for the operator j1 , j2 are maintained
in case of smooth solutions, but the right hand side Sobolev norms k · kl,2,T have to
be replaced by k · kl,2,ωT with a suitable patch ωT ⊃ T .
Sufficient conditions on the meshes Th , Mh , the FE and projection spaces for (5.29)(5.31) can be found in [17] or [1]. In particular, for the one-level approach with
Th = Mh , one has to enrich the velocity space by local bubble functions [17].
Another implication is that ju u 6∈ Vhdiv , hence the mixed term (ep , ∇ · εu ) has to
be considered. Moreover, a careful selection of the pressure spaces Qh is required.
The critical mixed term vanishes for continuous pressure space Qh = Pk−1 . (In case
of discontinuous space Qh = P−(k−1) , one can introduce additional pressure jump
terms across interior edges to handle it, see [1].)
5.3.1. Inf-sup stable case
(5.29)-(5.30) allow modified estimates of the skew-symmetric terms
(a · ∇εu , eu ) = −
X
κuh
(aM · ∇eu ) , εu ≤
X
M
− ((∇ × εb ) × d, eu ) =
M
X
εb , κbh
∇ × eu × dM
≤
1
kεu k20,2,M
τ1
X
M
− (∇ × (εu × d) , eb ) =
X
M
εu , κbh
((∇ × eb ) × dM ) ≤
M
X
M
1
2
kEh kStab ,
1
kεb k20,2,M
τ4
1
kεu k20,2,M
τ3
1
2
1
2
kEh kStab ,
kEh kStab .
Then, a modification of (5.11) leads to
S12 ≤ C
X
h2k
M
h2M
h2
h2
+ M |u|2k+1,2,M + λ + M |b|2k+1,2,M
τ1
τ3
τ4
ν+
M
(5.32)
and it still holds
S22 ≤ C
X
2(s−k)
2
h2k
M τ2 |u|k+1,2,ωM + hM
τ1 |aM |2 + τ4 |dM |2 |u|2s+1,2,ωM
(5.33)
M
1
h2 2(s−k)
+ |p|2k,2,ωM + τ6 + M |b|2k+1,2,ωM + hM
τ3 |dM |2 |b|2s+1,2,ωM .
τ2
τ5
A calibration of the parameters in (5.32) and (5.33) gives
2(k−s)
ch2M ≤ τ1 ≤ ChM
/|aM |2 ,
2(k−s)
Ch2M ≤ τ3 , τ4 ≤ ChM
τ2 ∼ 1,
/|dM |2
(5.34)
and as before
τ5 ∼ L20 /λ,
τ6 ∼ h2M λ/L20 .
17
(5.35)
Theorem 5.6. Let the orthogonality conditions (5.29)-(5.30) be valid. Assume that
the solution (u, b, p) of (2.11) belongs to [H k+1 (Ω)]d × [H k+1 (Ω)]d × H k (Ω). Further, let the LPS parameters be chosen according to conditions (5.34)-(5.35). Then
we obtain the quasi-optimal error estimate in Theorem 5.2 without the mesh-width
restrictions (5.14).
5.3.2. Equal-order case
We obtain
S12 ≤ C
X
h2k
M
h2
h2
h2M
+ M |u|2k+1,2,ωM + λ + M |b|2k+1,2,ωM
τ1
τ3
τ4
ν+
M
and
S22
≤ C
X
h2k
M
M
2(s−k)
+ τ7 hM
τ2 +
h2M 2
2(s−k)
|u|k+1,2,ωM + hM
τ1 |aM |2 + τ4 |dM |2 |u|2s+1,2,ωM
τ7
|p|2s,2,ωM +
h2 h2M 2
2(s−k)
|p|k+1,2,ωM + τ6 + M |b|2k+1,2,ωM + hM
τ3 |dM |2 |b|2s+1,2,ωM .
τ2
τ5
Similarly as in Subsec. 5.2.2 (ii), and here with s = k, we want to equilibrate
ν+
h2M
h2
h2
+ τ1 |aM |2 + τ4 |dM |2 + M + τ2 + M ∼ ν + hM |aM | + λ|d|M ,
τ1
τ3
τ7
h2
τ7 + M ∼ ν + hM |aM |,
τ2
2
2
h
h
λ + M + τ6 + M + τ3 |dM |2 ∼ λ + hM |dM |.
τ4
τ5
This leads to
hM
hM
hM
τ1 ∼
min {1, Ref,M } , τ3 ∼
min {1, Ref,M } , τ4 ∼
min {1, Rem,M } .
|aM |
|dM |
|dM |
(5.36)
The balance of τ2 and τ7 gives
τ2 ∼ τ7 ∼ max {hM |aM |, ν} .
(5.37)
As in the previous cases it holds
τ5 ∼ L20 /λ,
τ6 ∼ h2M λ/L20 .
(5.38)
Theorem 5.7. Let the orthogonality conditions (5.29)-(5.30) be valid. Assume that
the solution (u, b, p) of (2.11) belongs to [H k+1 (Ω)]d ×[H k+1 (Ω)]d ×H k+1 (Ω). Further,
let the LPS parameters be chosen according to conditions (5.36)-(5.38). Then we
obtain the quasi-optimal error estimate
kUh − JUk2G
+ kUh − JUk2Stab
≤C
X
h2k
Mσ
|u|2k+1,2,ωM
+ |b|2k+1,2,ωM
+ |p|2k+1,2,ωM
.
M
with σ := ν + λ + hM (|aM | + |dM |) and without the mesh-width restrictions (5.14).
18
5.4. Error estimates for singular solutions
A relevant problem of MHD is the occurence of nonsingular solutions b ∈ C with
b 6∈ [W 1,2 (Ω)]d , e.g. in non-convex domains. For simplicity, we restrict ourselves to
the stationary Maxwell problem with u ≡ 0, p ≡ 0, ∂t u ≡ 0. Following Badia/Codina
in [2], we introduce the Maxwell norm
k|(b, r)k| := kbkC + krkS .
A different estimate of the term M τ6 (∇ · (b − jb b, ∇ · ch )M with an appropriate
interpolation operator jb is given in [2], leading to
P

k|(bh −b, rh −r)k| ≤ C
inf
(wh ,sh )∈Ch ×Sh
k|(b − wh , r − sh )k| +
X
M
λhM
kb − wh k20,2,∂M
2
L0
!1 
According to Lemma 3.11 in [2], a solution b ∈ H0 (curl; Ω) ∩ H(div; Ω) can be
written as b = b0 + ∇φ with b0 ∈ [W 1+r,2 (Ω)]d ∩ H0 (curl; Ω) and φ ∈ W01,2 (Ω) ∩
W r,2 (Ω) for some r > 21 . Moreover, the following condition is assumed.
(A.6): There exists a FE space Gh defined on Th such that ∇φh ∈ Ch for any
φh ∈ Gh and that
inf kφ − φh ks,2,M ≤ Cht−s
M kφkt,2,M
φh ∈GH
for all M ∈ Th , for φ ∈ W t,2 (M ) and 0 ≤ s ≤ t ≤ k + 1.
This leads to
k|(b − bh , r − rh )k| ≤ C
X
√ t
λhM kb0 k1,2,M +
M ∈Th
√
λ
L1−
0
!
ht−
M kφk1+t,2,M .
Some variants of simplicial and quadrilateral/hexahedral elements fulfilling (A.6)
are discussed in [2] and [9].
The analysis in the present paper shows that the approach in [2] for the pure Maxwell
problem is compatible with the analysis of the LPS method for the stationary linearized MHD problem. In particular, the so-called cross-box elements can be handled
as a two-level LPS method where the LPS-compatibility condition (A.5) is valid.
6. Numerical examples
For the simulations, we apply the C++-FEM-packages deal.II [6] and FreeFEM++
[14]. The BDF(2)-scheme is used for the time discretization.
In this paper, we restrict ourselves to problems with analytical solutions in order
to test the error analysis. The main goal is to show the relevance of the parameter
design for τ2 , τ5 and τ6 . We investigate a two-dimensional time-dependent analytical
problem as in [5], the two-dimensional Hartmann channel flow [11], and the singular
Maxwell solution in a L-shaped domain [2].
19
2
.
6.1. Time discretization and Implementation
The rotational incremental version of a pressure-correction scheme [12] is applied
for the kinematic pressure, i.e., we introduce an auxiliary velocity vector and solve
n
e ht
the velocity equation with explicit treatment of the pressure: Find u
∈ Vh from
n−1
n−2
n
e ht
3u
− 4uht
+ uht
n
n
e ht
e ht
, vh + ν(∇u
, ∇vh ) + ((an · ∇)u
, vh )
2∆t
!
+
X
n
n
n
n
e ht
e ht
(τ1 (κsM ((aM
· ∇)u
), κsM ((aM
· ∇)vh ))M + τ2 (∇ · u
, ∇ · vh ) M
M
(6.1)
n
e ht
+τ4 (κsM (∇ × (u
× dnM )), κsM (∇ × (vh × dnM )))M )
∗
∗
= (fun , vh ) − (∇pn−1
ht , vh ) − (∇ × bht ) × bht , vh )
n−2
for all vh ∈ Vh where b∗ht := 2bn−1
ht − bht is an extrapolation.
Afterwards the solution is projected into the space of weakly solenoidal solutions:
n
, φnht ) ∈ (Vfh × Qh ) such that for all qh ∈ Qh
Find (uht
n
n
e ht
3uht
− 3u
+ ∇φnht , ∇qh = 0,
2∆t
n
(∇ · uht
, qh ) = 0
!
(6.2)
where Vfh differs from Vh by the boundary condition v · n|∂Ω = 0 for all v ∈ Vfh . Then
set the new pressure according to
n−1
n
e ht
pnht = φnht + pht
− ν∇ · u
.
For the magnetic part we apply the usual incremental implementation. We first solve
e n ∈ C from
for the magnetic field and treat the pseudopressure explicitly: Find b
h
ht
e n − 4bn−1 + bn−2
3b
ht
ht
ht
en , ∇ × c
, ch + λ ∇ × b
h
ht
2∆t
!
+
X
e n ) × dn ), κs ((∇ × c ) × dn )
τ3 κsM ((∇ × b
h
ht
M
M
M
M
en , ∇ · c
+τ6 ∇ · b
h
ht
=
(fbn , ch )
−
n−1
∇rht
, ch
M
(6.3)
M
n
e ht
+ (∇ × (u
× d) , ch )
for all ch ∈ Ch . Finally, we project this magnetic field in the space of weakly
n
solenoidal solutions: Find (bnht , rht
) ∈ (Ch × Sh ) such that for all sh ∈ Sh
n
X
3bnht − 3beht
n−1
n
+ ∇(rht
− rht
), ∇sh +
τ5 (∇rh , ∇sh )M = 0,
2∆t
M
!
(∇ ·
20
bnht , sh )
= 0.
(6.4)
6.2. Analytical two-dimensional time-dependent problem
The problem was solved on the domain Ω = (−1, 1) × (−1, 1) where the solution
πt
t
πt
t
u (t, x, y) = b (t, x, y) := x sin
exp
, −4x3 y sin
exp
10
25
10
25
2
2
:=
p (x, y)
x +y ,
r (x, y) := 0
4
T
,
does not depend on the relevant parameters λ and ν. The force terms are calculated
by evaluating the analytical solution. The boundary conditions for u are prescribed
by the given analytical solution for ∂Ω whereas the boundary conditions for b were
set by n × b = n × bD for ∂Ω where bD is the analytical solution for b. The initial
conditions for all unknowns are given by the analytical solution as well.
We have used the time interval [0, 0.1] with time-step size ∆t = 10−4 in order to
avoid order-reduction effects for the spatial convergence rates. The coefficients of
the physical problem are set as ν = 10−6 and λ = 10−6 in Fig. 1. The FE spaces
[Q)2 ]2 /Q1 are used for Vh /Qh and Ch /Sh as the classical Taylor-Hood elements. The
2
L2
only stabilization parameters in usage are τ2 = 1, τ5 = hL2λ and τ6 = λ0 where we
0
set L0 = 2 as length scale.
We see in Fig. 1 the quasi-optimal error estimates as predicted by the error analysis.
This holds also true for even smaller viscosity ν and magnetic coefficient λ which
demonstrates the robustness of our proposed stabilization method.
6.3. Two-dimensional Hartmann channel flow
The two-dimensional Hartmann channel flow problem was solved on the domain
Ω = (0, 10) × (−1, 1) where the stationary solution given by


10Re

u (x, y) := 
 0.5
Re
λ



tanh
 sinh

b (x, y) := 
10 
sinh
Re
λ
Re
λ
0.5 Re
λ
cosh

0.5 
1.0 −
y
0.5 cosh

T


Re
λ
0.5 
Re
λ
y 

T

0.5  , 0

,

− y
 , 1 ,
and p (x, y) := −10x − 21 [b1 (x, y)]2 and r (x, y) := 0 depends on λ and ν. The righthand-side force terms are given by fu = 0 and fb = 0. The physical coefficients
are set as ν = 0.5 and λ = 0.5, i.e. Re = 2. The boundary conditions for u are
prescribed by the given analytical solution for ∂Ω whereas the boundary conditions
for b were set by n × b = n × bD on ∂Ω where bD is the analytical solution for b
as well. The initial conditions are given by u0 = 0, b0 = 0, p0 = 0 and r0 = 0.
We have used the time interval [0, 10] with time-step ∆t = 0.1. The FE spaces
21
10 2
10 1
L2 (u)
L2 (b)
h3
10 0
H 1 (u)
L2 (curl(b))
h2
10 0
10 -1
10 -2
10 -2
10 -4
10 -3
10 -6
10 -8
10 -2
10 -4
10 -1
10 0
10 -5
10 -2
10 1
10 1
10 -1
10 0
10 1
10 1
L2 (div(u))
2
L (div(b))
2
h
10 0
10
10 -1
L2 (p)
2
h
0
10 -1
10 -2
10 -2
10 -3
10 -3
10 -4
10 -4
10 -5
10 -6
10 -2
10 -1
10 0
10 -5
10 -2
10 1
10 -1
10 0
10 1
Figure 1: Errors for the 2D-time-dependent analytical problem
[Q2 ]2 /Q1 are used for Vh /Qh and Ch /Sh as classical Taylor-Hood elements. The
2
L2
only stabilization parameters in usage are τ2 = 1, τ5 = hL2λ and τ6 = λ0 where we
0
set L0 = 1 as the characteristic length scale.
We see in Fig. 2 the desired quasi-optimal error estimates as predicted by the error
analysis. The worse results on the finest meshes result from the coarse time-step.
6.4. Two-dimensional singular solution in the L-shaped domain
Here we investigate the singular solution in the L-shaped domain Ω = (−1, 1)2 \
{[0, 1] × [−1, 0]} where the exact solutions are given in polar coordinates (d, ϕ) by
2ϕ
b (d, ϕ) := d sin
3
2
3
, r (d, ϕ) = 0.
The physical parameter is given by λ = 1 at first for Maxwell’s
equationi where
h 2
2
u = 0 and p = 0. It is also fb = 0. It is well-known that b ∈ W 3 ,2 (Ω) holds
h
i2
and we don’t have b ∈ W1,2 (Ω) . The Dirichlet conditions for r is described by
r = 0 on ∂Ω and the boundary conditions for the magnetic field are n × b = n × bD
where bD is the exact solution on ∂Ω. The above presented problem is solved by a
stationary solver written in FreeFEM++. The FE spaces in usage are P1 /P1 in
22
10 0
10 1
L2 (u)
L2 (b)
h3
10 -1
H 1 (u)
L2 (curl(b))
h2
10 0
10 -2
10 -1
10 -3
10 -2
10 -4
10 -5
10 -2
10 -1
10 0
10
10 -3
10 -2
10 0
10 -1
10 0
10 0
L2 (div(u))
2
L (div(b))
2
h
-1
L2 (p)
2
h
10 -2
10 -1
10 -3
10 -4
10 -2
10 -5
10 -6
10 -2
10 -1
10 0
10 -3
10 -2
10 -1
10 0
Figure 2: Errors for the 2D stationary Hartmann flow problem
the case of equal-order interpolation P2 /P1 in the Taylor-Hood case for b and r. We
first use different meshes for the calculations as seen in the Fig. 3 below.
We receive the following results with the rate of convergence (CRC) in brackets.
−4
h
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
2
1.71E − 1
1.93E − 2
7.94E − 3
3.91E − 2
3.65
h
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
2−4
1.60E − 1
1.55E − 2
6.25E − 3
3.04E − 2
3.52
Cross-box, Equal-order
2
2−6
1.09E − 1 (0.65) 6.91E − 2 (0.66)
7.66E − 3 (1.33) 3.04E − 3 (1.33)
3.46E − 3 (1.20) 1.45E − 3 (1.25)
2.49E − 2 (0.65) 1.57E − 2 (0.66)
4.57 (−0.32)
5.75 (−0.33)
Cross-box, Taylor-Hood
2−5
2−6
1.02E − 1 (0.65) 6.45E − 2 (0.66)
6.13E − 3 (1.34) 2.42E − 3 (1.34)
2.74E − 3 (1.19) 1.15E − 3 (1.25)
1.94E − 2 (0.65) 1.23E − 2 (0.66)
4.41 (−0.33)
5.54 (−0.33)
−5
23
2−7
4.36E − 2 (0.66)
1.23E − 3 (1.31)
5.96E − 4 (1.28)
9.92E − 3 (0.66)
7.23 (−0.33)
2−7
4.07E − 2 (0.66)
9.61E − 4 (1.33)
4.71E − 4 (1.29)
7.74E − 3 (0.67)
6.98 (−0.33)
Figure 3: Different meshes for the calculations: (left top) Cross-box, (right top) P1 standard mesh,
(left bottom) P1 right modified mesh, (right bottom) P1 left modified mesh
h
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
2−4
3.62E − 1
1.19E − 1
2.33E − 2
7.98E − 2
3.06
h
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
2−4
2.19E − 1
2.87E − 2
1.00E − 2
3.81E − 2
3.18
P1 standard mesh, Equal-order
2−5
2−6
2.86E − 1 (0.34) 2.23E − 1 (0.36)
7.98E − 2 (0.58) 5.68E − 2 (0.49)
1.46E − 2 (0.67) 8.83E − 3 (0.73)
6.39E − 2 (0.32) 5.07E − 2 (0.33)
3.94 (−0.36)
5.24 (−0.41)
P1 standard mesh, Taylor-Hood
2−5
2−6
1.41E − 1 (0.64) 8.98E − 2 (0.65)
1.16E − 2 (1.31) 4.72E − 3 (1.30)
4.58E − 3 (1.13) 1.96E − 3 (1.22)
2.46E − 2 (0.63) 1.56E − 2 (0.66)
3.99 (−0.33)
5.02 (−0.33)
2−7
1.74E − 1 (0.36)
4.19E − 2 (0.44)
5.30E − 3 (0.74)
4.00E − 2 (0.34)
7.12 (−0.44)
2−7
5.68E − 2 (0.66)
1.92E − 3 (1.30)
8.13E − 4 (1.27)
9.89E − 3 (0.66)
6.32 (−0.33)
For the right- and left-modified P1 -meshes, we obtain similar convergence rates (not
shown) as for the standard P1 -mesh. The same optimal convergence rates can also
be obtained for the 3D-case in the case of Taylor-Hood elements.
From our calculation, we can observe that the cross-box element works fine in accordance to [2]. We surprisingly also see that the Taylor-Hood approximation in P2 /P1
24
reaches the asymptotically optimal convergence rates as well in contrast to the equal
order approximation in P1 /P1 on the other given meshes.
Now we look at the choice of the length scale L0 . We use the standard P1 mesh
to check different length scales. We further take the Taylor-Hood pair P2 /P1 into
account for our calculations. Let λ = 1 be the physical parameter.
L0 = 100
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
h = 2−4
2.19E − 1
2.87E − 2
1.00E − 2
3.81E − 2
3.18
h = 2−5
1.41E − 1 (0.64)
1.16E − 2 (1.31)
4.58E − 3 (1.13)
2.46E − 2 (0.63)
3.99 (−0.33)
h = 2−6
8.98E − 2 (0.65)
4.72E − 3 (1.30)
1.96E − 3 (1.22)
1.56E − 2 (0.66)
5.02 (−0.33)
L0 = 101
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
h = 2−4
2.25E − 1
1.30E − 3
4.15E − 4
1.58E − 3
3.28
h = 2−5
1.43E − 1 (0.65)
5.29E − 4 (1.30)
1.86E − 4 (1.16)
9.97E − 4 (0.66)
4.04 (−0.30)
h = 2−6
9.02E − 2 (0.66)
2.18E − 4 (1.28)
7.89E − 5 (1.24)
6.29E − 4 (0.66)
5.05 (−0.32)
L0 = 105
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
h = 2−4
h = 2−5
h = 2−6
2.25E − 1 1.43E − 1 (0.65) 9.08E − 2 (0.66)
1.33E − 11 4.97E − 12 (1.42) 1.96E − 12 (1.34)
4.14E − 12 1.87E − 12 (1.15) 8.10E − 13 (1.21)
1.58E − 11 1.00E − 11 (0.66) 6.39E − 12 (0.65)
3.28
4.03 (−0.30)
5.11 (−0.34)
L0 = 2 · 105
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
h = 2−4
h = 2−5
h = 2−6
2.26E − 1 1.43E − 1 (0.65) 1.04E − 1 (0.46)
3.19E − 12 1.27E − 12 (1.33) 1.25E − 12 (0.02)
1.04E − 12 4.62E − 13 (1.17) 2.81E − 13 (0.72)
3.95E − 12 2.49E − 12 (0.67) 1.99E − 12 (0.32)
3.28
4.07 (−0.31)
5.41 (−0.41)
L0 = 106
kb − bh k0
k∇ × (b − bh )k0
kr − rh k0
k∇ (r − rh )k0
k∇ · (b − bh )k0
h = 2−4
2.60E − 1
1.66E − 13
4.53E − 14
1.96E − 13
3.33
h = 2−5
h = 2−6
3.65E − 1 (−0.49)
2.96E − 13 (−0.83)
4.84E − 14 (−0.10)
3.38E − 13 (−0.79)
10.04 (−1.59)
-
We observe that the results are reasonable for a wide range of L0 ∈ [100 , 105 ] with
25
dramatic improvement for the curl-error k∇ × (b − bh )k0,2 and for the pseudopressure with the tendency towards to machine accuracy for increasing values of L0 .
For higher values of L0 , we found a breakdown of the results for the magnetic field
and the pseudo-pressure.
Finally, we considered the choice of length scale L0 in dependence of parameter λ
(not shown). Exemplarily, for λ = 10−6 , we obtained similar results as for λ = 1 in
the range of L0 ∈ [100 , 103 ]. For L0 = 104 we found the mentioned breakdown of
kb − bh k0,2 on the finest mesh and for L0 ≥ 1.5 · 104 for all unknowns with exception
of the pseudo-pressure.
This demonstrates that the model is in some way sensitive in accordance to the
length scale in relation to the mesh refinement.
7. Summary, conclusions and outlook
We considered nodal-based FE solutions to the linearized MHD problem in case of
inf-sup stable approximation of kinematic velocity and pressure and of magnetic field
and pseudo-pressure. Contrary to previous residual-based stabilization, the focus
was here on symmetric local projection stabilization. For comparison purposes, we
considered the case of equal-order interpolation of all unknowns as well.
We obtained the following results and conclusions:
i) For the stabilization of the divergence-free constraint of the Maxwell part and
of the magnetic pseudo-pressure, one can basically take the parameter design
proposed by Badia/Codina in [2]. The selection of a characteristic length scale
still allows some freedom.
ii) For the important case of singular magnetic solutions, we followed the approach proposed in [2]. An interesting observation is that – contrary to Q2 /Q1
interpolation of the magnetic field and of the magnetic pseudo-pressure – the
application of P2 /P1 gives the theoretical error orders without using cross-box
elements. Moreover, cross-box elements fulfill the LPS compatibility condition
(A.5) in the framework of two-level LPS methods.
iii) For methods without LPS compatibility condition (A.5), we derived semirobust error estimates of order k for inf-sup-stable approximations provided
some (relatively mild) mesh width restrictions are fulfilled. Given the parameter design for the Maxwell part, see i), a relatively wide range for the stabilization parameters of the skew-symmetric terms (convective and Lorentz terms)
is possible. In particular, vanishing stabilization of these terms is possible. We
observed a reasonable behavior of the method for the latter situation for the
separated cases of the momentum equation, see [1], and of the induction equation with Lorentz term, see [15].
iv) The mentioned mesh width restrictions can be completely removed in case of
LPS compatibility condition (A.5) thus leading to semi-robust error estimates
of order k. In case of one-level LPS methods, one has to enrich the FE spaces.
26
The allowed range of the stabilization parameters for the skew-symmtric terms
is shrinked due to lower bounds on the parameters. Further research is required
in order to fix the parameter sets, in particular for high fluid and/or magnetic
Reynolds numbers.
v) In case of equal-order approximations of all unknowns one obtains a very clear
design of the parameter sets provided methods with LPS compatibility condition (A.5) are considered. The parameter design is simplified for LPS stabilization as compared to the residual-based approach, see, e.g., [3–5].
We plan to extend the approach to 3D examples and to more complex problems.
In particular, we will consider the necessity of LPS terms for the skew-symmetric
terms, see [1] for the pure fluid part.
Acknowledgements: The authors would like to thank Utku Kaya for providing
the numerical results with FreeFEM++ and helpful discussions.
The research of D. Arndt is supported by CRC 963 via German Research Foundation
(DFG). The research of B. Wacker was supported by RTG 1023 via German Research
Foundation (DFG).
References
1. D. Arndt, H. Dallmann, and G. Lube, Local projection FEM stabilization for the time-dependent incompressible Navier-Stokes problem, accepted for:
Numer. Meths. Part. Diff. Equat., (2014).
2. S. Badia, R. Codina, A nodal-based finite element approximation of the
Maxwell problem suitable for singular solutions, SIAM J. on Numer. Anal.,
Vol. 50 (2012), 398–417.
3. S. Badia, R. Codina and R. Planas, On an unconditionally convergent
stabilized finite element approximation of resistive magnetohydrodynamics, J.
Comp. Phys., Vol. 234, (2013), 399–416.
4. S. Badia, R. Codina and R. Planas, Analysis of an unconditionally convergent stabilized finite element formulation for incompressible magnetohydrodynamics, submitted.
5. S. Badia, R. Planas and J. V. Gutierrez-Santacreu, Unconditionally stable operator splitting algorithms for the incompressible magnetohydrodynamics (MHD) system discretized by a stabilized finite element formulation
based on projections, Intern. J. Numer. Meth. Engrg., 93 (2013), 302-328.
6. W. Bangerth, T. Heister and G. Kanschat, deal.II - Differential Equations Analysis Library, Technical Reference, http://www.dealii.org.
7. A. Bonito, J.-L. Guermond, Approximation of the eigenvalue problem for
time harmonic Maxwell systems by continuous Lagrange finite elements, Math.
Comp. 80 (2011) 276, 1887-1910.
27
8. A. Bonito, J.-L. Guermond, F. Luddens, Regularity of the Maxwell equations in heterogeneous media and Lipschitz domains, J. Math. Anal. Appl. 408
(2013) 2, 498-512.
9. M. Costabel, M. Dauge, Weighted regularization of Maxwell equations in
polyhedral domains, Numer. Math. 93 (2002) 2, 239-278.
10. A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements,
Springer-Verlag, New York, (2004).
11. C. Greif, D. Li, D. Schötzau and X. Wei, A mixed finite element method
with exactly divergence-free velocities for incompressible magnetohydrodynamics,
Comput. Methods Appl. Mech. Engrg., Vol. 199 (2010), 2840-2855.
12. J.L. Guermond and J. Shen, On the error estimates for the rotational
pressure-correction projection methods, Math. Comp. 73 (2004) 248, 1719-1737.
13. V. Girault, R. Scott, A quasi-local interpolation operator preserving the
discrete divergence, Calcolo, Vol. 40 (2003), 1-19.
14. F. Hecht, New development in FreeFem++, J. Numer. Math., 20 (2012), No.
3-4, 251-265. 65Y15.
15. U. Kaya, Numerical simulation of the induction equation using Lagrangian
finite elements, Master Thesis, Georg-August University Göttingen, NAM, 2014.
16. G. Lube, G. Rapin and J. Löwe ,Local projection stabilization for incompressible flows: Equal-order vs. inf-sup stable interpolation, ETNA 32 (2008),
106-122.
17. G. Matthies, P. Skrzypacz and L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, ESAIM:M2AN
41 (2007) 4, 713-742.
18. G. Matthies and L. Tobiska, Local projection type stabilisation applied to
inf-sup stable discretisations of the Oseen problem, IMA J. Numer. Anal. (2014)
doi:10.1093/imanum/drt064.
19. B. Wacker, G. Lube, A local projection stabilization FEM for the linearized
stationary MHD problem, in: Numerical Mathematics and Advanced Applications, Springer Lect. Notes Comput. Sci. Engin. 103 (2015), 765-774.
28
Institut für Numerische und Angewandte Mathematik
Universität Göttingen
Lotzestr. 16-18
D - 37083 Göttingen
Telefon:
Telefax:
0551/394512
0551/393944
Email: [email protected] URL: http://www.num.math.uni-goettingen.de
Verzeichnis der erschienenen Preprints 2015
Number
Authors
Title
2015-1
Behrends, S., Hübner, R., Schö- Norm Bounds and Underestimators for Unconsbel, A.
trained Polynomial Integer Minimization
2015-2
J. Harbering, A. Ranade, M. Single Track Train Scheduling
Schmidt
2015-3
Wacker, B., Arndt, D., Lube, Nodal-based nite element methods with local
G.
projection stabilization for linearized incompressible magnetohydrodynamics