Divergence-free H(div)-conforming hierarchical bases for magnetohydrodynamics (MHD)

0.17 MB·

by  Wei Cai

Divergence-free H(div)-conforming hierarchical bases
for Magnetohydrodynamics (MHD)
Wei Cai†∗, Jian Wu and Jianguo Xin
†Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223, USA.
2
1
0 Abstract
2
In order to solve the magnetohydrodynamics (MHD) equations with a H(div)-conforming
t
c element, a novel approach is proposed to ensure the exact divergence-free condition on the
O magnetic field. The idea is to add on each element an extra interior bubble function from
higherorderhierarchicalH(div)-conformingbasis. FoursuchhierarchicalbasesfortheH(div)-
0
conforming quadrilateral, triangular, hexahedral and tetrahedral elements are either proposed
2
(inthecaseoftetrahedral)orreviewed. Numericalresultshavebeenpresentedtoshowthelinear

] independence of the basis functions for the two simplicial elements. Good matrix conditioning
A
has been confirmed numerically up to the fourth order for the triangularelement and up to the
N
third order for the tetrahedral element.
.
h Keywords: Hierarchical bases, H(div)-conforming elements, Divergence-free condition
t
a
m
1 Introduction
[

1 The magnetohydrodynamics (MHD) equations describe the dynamics of a charged system under
v
the interaction with a magnetic field and the conservation of the mass, momentum and energy for
5
7 the plasma system. Such a dynamics is considered constrained as the magnetic field of the system
5
isevolved withtheconstraint ofzerodivergence, namely, ∇·B= 0. Numericalmodelingofplasmas
5
. has shown that the observance of the zero divergence of the magnetic field plays an important role
0
in reproducing the correct physics in the plasma fluid [1]. Various numerical techniques have been
1
2 devised to ensure the computed magnetic field to maintain divergence-free [2]. In the original work
1 of [1] a projection approach was used to correct the magnetic field to have a zero divergence.
:
v A more natural way to satisfy this constraint is through a class of the so-called constrained
i
X transport (CT) numerical methods based on the ideas in [3]. As noted in [4], a piecewise H(div)
r vector field on a finite element triangulation of a spatial domain can be a global H(div) field if
a
and only if the normal components on the interface of adjacent elements are continuous. Thus, in
most of the CT algorithms for the MHD, the surface averaged magnetic flux over the surfaces of a
3-D element is used to represent the magnetic field while the volume averaged conserved quantities
(mass, momentum, and energy) are used.
In the two seminal papers [5, 6], N´ed´elec proposed to use quantities (moments of tangential
components of vector fields) on edges and faces to define the finite dimensional space in H(div)
andH(curl). Thespecificconstruction of thebasis functionsinbothspaces, specifically in H(div),
can be done in various ways such as the hierarchical type of basis proposed in [7] for both H(div)
and H(curl). Unfortunately, the proposed hierarchical basis for H(div) in [7] turns out to be
erroneous as it can be easily checked that for quadratic polynomial approximation the proposed
edge-based basis functions happens to be linearly dependent.
∗Corresponding author. Tel.: [phone]; fax: [phone]. E-mail address: [email protected] (W. Cai).
1
In this paper, we will first present a new hierarchical basis functions in H(div) in tetrahedral,
for completeness, together with a review of the hierarchical basis of H(div) for rectangles in 2-D
and for cube in 3-D. An important common feature of the hierarchical basis functions is the fact
that for order p ≥ 3 in the case of simplexes, the basis functions will include interior bubble basis
functions which have zero normal components on the whole boundary of the element boundaries.
Therefore, a simple way to enforce the divergence-free on each element can be easily accomplished
by adding one single (p+1)-th order (any q-th order, q ≥ max(m,p +1),m = 4 for tetrahedral
element, m = 3 for triangular element, and m = 2 for quadrilateral or hexahedral element) interior
bubble function to a p-th order H(div) basis. This extra bubble basis will be able to satisfy the
local divergence-free condition. Due to the fact of normal continuity of the p-th order basis across
the element interface and zero normal component of the added-in p+1-th order interior bubble
functions, the augmented function space satisfies divergence-free condition globally.
The rest of the paper is organized as follows. The constructions of the locally divergence-
free bases are given in Section 2-5 for rectangular and triangular elements in 2-D, and cubic and
tetrahedral elements in 3-D. The divergence-free condition is discussed in Section 6. Numerical
results on the matrix conditioning are given in Section 7. Concluding remarks are given in Section
8.
2 Basis functions for the quadrilateral element
In [8] Zaglmayr gave a hierarchical basis for quadrilateral H(div)-conforming element. In this
section we summarize the result in [8].
The basis functions are constructed on the reference element - a unit square Q := [0,1]2 [8]
with vertexes of V (0,0), V (1,0), V (1,1) and V (0,1). The coordinate system for the reference
1 2 3 4
element is in terms of the variables (ξ,η). A bilinear function λ which is associated with a specific
i
vertex V has been utilized for the construction. The multiplicative factors of the bilinear function
i
λ have been used to form the linear function σ , viz.
i i
λ := (1−ξ)(1−η), σ := (1−ξ)+(1−η),
1 1
λ := ξ(1−η), σ := ξ+(1−η),
2 2
λ := ξη, σ := ξ+η,
3 3
λ := (1−ξ)η, σ := (1−ξ)+η. (1)
4 4
The bilinear function has this favorable property
λ | = δ , (2)
i Vj ij
where δ is the Kronecker delta. The edge e := [V ,V ] which points from vertex V to vertex V is
ij i j i j
parameterized by
ζ := σ −σ ∈ [−1,1]. (3)
e j i
For convenience of basis construction, the linear edge-extension parameter is also defined, viz.
λ := λ +λ ∈ [0,1], (4)
e i j
2
which is one on edge e := [V ,V ] and zero on the opposite edge. Note that the unit tangential
i j
vector τ and the outward unit normal vector n can be deduced as
e e
1
~τ = ∇ζ , n = ∇λ . (5)
e e e e
2
2.1 Edge-based functions
These functions are further grouped into two categories: the lowest-order and higher-order func-
tions.
Lowest-order functions:
These functions are associated with the four edges. By construction each function is perpendic-
ular to the associated edge and has unit normal component on the associated edge. Furthermore
the divergence of each function is unit. The shape function is given by
1
ψRT0 = λ (∇×ζ ), i= {1,2,3,4}, (6)
ei 2 ei ei
which has the property
~τ ·ψRT0 = 0, n ·ψRT0| = 1, ∇·ψRT0 = 1. (7)
ei ei ei ei ei ei
Higher-order functions:
The function is taken to be a curl field of a scalar function in order to be free of divergence.
The basis function is
ψj+1 = ∇×(λ L (ζ )), i = {1,2,3,4}, 0 ≤ j ≤ p−1, (8)
ei ei j+2 ei
where the function L (•) is the so-called integrated Legendre polynomial of degree n [8]. The
n
obvious property is
∇·ψj+1 = 0, i = {1,2,3,4}, 0 ≤ j ≤ p−1. (9)
ei
2.2 Interior functions
The interior functions are classified into three categories.
Type 1: (curl field)
The shape functions are given as
ψQ1 = ∇×(L (2ξ −1)L (2η−1)), 0 ≤ i,j ≤ p−1. (10)
ij i+2 j+2
These functions are divergence free, viz.
∇·ψQ1 = 0, 0 ≤ i,j ≤ p−1. (11)
ij
3
Type 2:
The formula of these functions is given as
ψQ2 = L (2ξ −1)L′ (2η−1)ˆi +L′ (2ξ −1)L (2η−1)ˆj , 0 ≤ i,j ≤ p−1. (12)
ij i+2 j+2 ξ i+2 j+2 η
On the boundary of the reference element the normal component of these functions vanishes, viz.
n ·ψQ2| = 0, 0 ≤ i,j ≤ p−1, (13)
e ij e
which justifies to be one type of interior functions.
Type 3:
The formula of these functions is
ψQξ3 = L (2ξ −1)ˆi , ψQη3 = L (2η−1)ˆj , 0 ≤ i≤ p−1. (14)
i i+2 ξ i i+2 η
Again on the boundary of the reference element the normal component of these functions vanishes,
viz.
n ·ψQξ3| = 0, n ·ψQη3| = 0, 0 ≤i ≤ p−1. (15)
e i e e i e
3 Basis functions for the hexahedral element
The reference element is defined for a unit cube H := [0,1]3 in [8]. The vertexes of the cube are
V (0,0,0), V (1,0,0), V (1,1,0), V (0,1,0), V (0,0,1), V (1,0,1), V (1,1,1), and V (0,1,1). The
1 2 3 4 5 6 7 8
basis functions are expressed in terms of the trilinear function λ , which is one at the vertex V
j j
and zero at all other vertexes. Along with the linear function σ and in terms of the coordinates
j
variables (ξ,η,ζ) they are given as
λ := (1−ξ)(1−η)(1−ζ), σ := (1−ξ)+(1−η)+(1−ζ),
1 1
λ := ξ(1−η)(1−ζ), σ := ξ+(1−η)+(1−ζ),
2 2
λ := ξη(1−ζ), σ := ξ+η+(1−ζ),
3 3
λ := (1−ξ)η(1−ζ), σ := (1−ξ)+η+(1−ζ).
4 4
λ := (1−ξ)(1−η)ζ, σ := (1−ξ)+(1−η)+ζ,
5 5
λ := ξ(1−η)ζ, σ := ξ+(1−η)+ζ,
6 6
λ := ξηζ, σ := ξ+η+ζ,
7 7
λ := (1−ξ)ηζ, σ := (1−ξ)+η+ζ. (16)
8 8
The edge e := [V ,V ] which points from vertex V to vertex V is parameterized by
i j i j
µ := σ −σ ∈ [−1,1]. (17)
e j i
The tangential vector associated with the edge e is given by ~τ = 1∇µ . The edge extension
e 2 e
parameter λ := λ +λ ∈ [0,1] is one on the edge e and zero on all other edges that are parallel to
e i j
4
the edge e. The face f = [V ,V ,V ,V ] where the vertexes V and V are not connected by an edge
i j k ℓ i k
can be parameterized by
(ξ ,η ) := (σ −σ ,σ −σ ) ∈ [−1,1]×[−1,1]. (18)
f f i j i ℓ
The linear face extension parameter λ = λ +λ +λ +λ is equal to one on the face f and zero
f i j k ℓ
on the opposite face. The outward unit normal vector of face f can be obtained by n = ∇λ .
f f
3.1 Face-based functions
In this subsection we record the results in [8]. We have also fixed one error in [8]. These functions
are associated with the six faces whose formulas are classified into two groups.
Lowest-order Raviart-Thomas functions:
ψRT0 = λ n , i = 1,2,··· ,6. (19)
fi f f
Higher-order functions: (divergence-free)
These functions are constructed as curl fields of certain vectors in order to be divergence-free.
The formulas are given as [8]
f
ψ k = ∇×(λ (L (η )∇L (ξ )−L (ξ )∇L (η ))), 0 ≤ i,j ≤ p−1, k = 1,2,··· ,6.
i,j f j+2 f i+2 f i+2 f j+2 f
(20)
f
ψ k = ∇×(λ L (ξ )∇η ), 0≤ i ≤ p−1, k = 1,2,··· ,6. (21)
i f i+2 f f
f
ψ k = ∇×(λ L (η )∇ξ ), 0 ≤ j ≤ p−1, k = 1,2,··· ,6. (22)
j f j+2 f f
3.2 Interior functions
The interior functions are further classified into three categories. The triplet (ξ ,η ,ζ ) := (2ξ −
1 2 3
1,2η−1,2ζ−1) is used in the formulas. The function ℓ (•) is the classical un-normalized Legendre
n
polynomialof degreen. While wehererecord theresultsin [8], wehave implemented thecorrection
of a number of mistakes in [8] as well.
Type 1: (divergence-free)
These functions are taken to be the curl fields of certain vector functions that are associated
with the H(curl)-conforming element.
ψH11 = 4L (ξ )ℓ (η )ℓ (ζ )ˆi −4ℓ (ξ )ℓ (η )L (ζ )kˆ , 0 ≤ i,j,k ≤ p−1.
i,j,k i+2 1 j+1 2 k+1 3 ξ i+1 1 j+1 2 k+2 3 ζ
ψH21 = 4ℓ (ξ )L (η )ℓ (ζ )ˆj −4ℓ (ξ )ℓ (η )L (ζ )kˆ , 0 ≤ i,j,k ≤ p−1.
i,j,k i+1 1 j+2 2 k+1 3 η i+1 1 j+1 2 k+2 3 ζ
ψH31 = 2ℓ (η )L (ζ )kˆ −2L (η )ℓ (ζ )ˆj , 0 ≤ j,k ≤ p−1.
j,k j+1 2 k+2 3 ζ j+2 2 k+1 3 η
5
ψH41 = 2L (ξ )ℓ (ζ )ˆi −2ℓ (ξ )L (ζ )kˆ , 0 ≤ i,k ≤ p−1.
i,k i+2 1 k+1 3 ξ i+1 1 k+2 3 ζ
ψH51 = 2L (ξ )ℓ (η )ˆi −2ℓ (ξ )L (η )ˆj , 0 ≤ i,j ≤ p−1. (23)
i,j i+2 1 j+1 2 ξ i+1 1 j+2 2 η
Type 2:
These functions are linear combinations of certain components in the above type.
ψH12 = L (ξ )ℓ (η )ℓ (ζ )ˆi +ℓ (ξ )L (η )ℓ (ζ )ˆj , 0 ≤ i,j,k ≤ p−1.
i,j,k i+2 1 j+1 2 k+1 3 ξ i+1 1 j+2 2 k+1 3 η
ψH22 = L (η )ℓ (ζ )ˆj +ℓ (η )L (ζ )kˆ , 0≤ j,k ≤ p−1.
j,k j+2 2 k+1 3 η j+1 2 k+2 3 ζ
ψH32 = ℓ (ξ )L (ζ )kˆ +L (ξ )ℓ (ζ )ˆi , 0 ≤ i,k ≤ p−1.
i,k i+1 1 k+2 3 ζ i+2 1 k+1 3 ξ
ψH42 = L (ξ )ℓ (η )ˆi +ℓ (ξ )L (η )ˆj , 0 ≤ i,j ≤ p−1. (24)
i,j i+2 1 j+1 2 ξ i+1 1 j+2 2 η
Type 3:
These functions are taken as certain components in Type 2.
ψH13 = L (ξ )ˆi , 0 ≤ i≤ p−1.
i i+2 1 ξ
ψH23 = L (η )ˆj , 0≤ j ≤ p−1.
j j+2 2 η
ψH33 = L (ζ )kˆ , 0 ≤ k ≤ p−1. (25)
k k+2 3 ζ
Two remarks are in place.
• All the interior basis functions are linearly independent, which can be verified easily.
• The normal traces of these interior functions vanish on the boundary ∂H of the reference
hexahedral element. This is dueto the fact that on a certain face f either one of the standard
unit vectors is perpendicular to the normal vector of the face n or to the fact that the
f
integrated Legendre polynomials are evaluated at 1 or −1, which is zero.
4 Basis functions for the triangular element
The result on the basis construction has been reported in [9]. For the completeness of the current
study, we record the basis functions in this section.
Any point inthe 2-simplex K2 is uniquelylocated in terms of the local coordinate system (ξ,η).
The vertexes are numbered as v (0,0),v (1,0),v (0,1). The barycentric coordinates are given as
0 1 2
λ := 1−ξ−η, λ := ξ, λ := η. (26)
0 1 2
The directed tangent on a generic edge e = [j ,j ] is similarly defined as in (40) for the three-
j 1 2
dimensional case. In the same manner the edge is also parameterized as in (41). A generic edge
can be uniquely identified with
e := [j ,j ], j = {0,1}, j < j ≤ 2, j = j +j . (27)
j 1 2 1 1 2 1 2
6
The two-dimensional vectorial curl operator of a scalar quantity, which is used in our construction,
needs a proper definition. We use the one given in the book [11], viz.
τ
∂u ∂u
curl(u) := ∇×u:= , − (28)
∂η ∂ξ
(cid:20) (cid:21)
Based upon the newly created shape functions for the three-dimensional H(div)-conforming tetra-
hedral elements and using the technique of dimension reduction we construct the basis for the
H(div)-conforming triangular elements in two dimensions. However, it is easy to see that the two
groups for the face functions cannot be appropriately modified for our purpose. Instead we borrow
the idea of Zaglmayr in the dissertation [8], viz., we combine the edge-based shape functions in [8]
with our newly constructed edge-based and bubble interior functions. In [8] Zaglmayr had applied
the so-called scaled integrated Legendre polynomials in the construction, viz.
x
ξ
Ls(x,t) := tn−1 ℓ dξ, n≥ 2, t ∈ (0,1], (29)
n n−1 t
Z (cid:18) (cid:19)
−t
where ℓ (x) is the n-th order Legender polynomial.
n
4.1 Edge functions
For the completeness of our basis construction, in this subsection we record the results in [8].
Associated with each edge the formulas for these functions are given as
ΦN0 = λ ∇×λ −λ ∇×λ (30)
e[k1,k2] k2 k1 k1 k2
for the lowest-order approximation and
Φje[k1,k2] = ∇× Lsj+2(γek,λk2 +λk1) , j = 0,··· ,p−1 (31)
(cid:0) (cid:1)
for higher-order approximations.
4.2 Interior functions
The interior functions are further classified into two categories: edge-based and bubble interior
functions. By construction the normal component of each interior function vanishes on either edge
of the reference 2-simplex K2, viz.
e t
n j ·Φ = 0, j = {1,2,3}, (32)
e
where n j is the unit outward normal vector to edge ej.
Edge-based interior functions:
The tangential component of each edge-based function does not vanish on the associated only
edge e := [k ,k ] but vanishes the other two edges, viz.
k 1 2
τej ·Φt,i = 0, e 6= e , (33)
e[k1,k2] j k
7
e
where τ j is the directed tangent along the edge ej := [j1,j2]. The following basis functions are
proposed here as
e
Φt,i = C λ λ (1−λ )iP(0,2) 2λk2 −1 τ k , (34)
e[k1,k2] i k1 k2 k1 i (cid:18)1−λk1 (cid:19) |τek|
(0,2)
where the function P (•) is the classical un-normalized Jacobi polynomial of degree i with a
i
single variable [10], and the scaling coefficient is given by
C = 2(i+2)(i+3)(2i+3)(2i+5), i = 0,1,··· ,p−2. (35)
i
p
The following orthonormal property of edge-based interior functions can be proved
t,m t,n
< Φe[k1,k2],Φe[k1,k2] > |K2 = δmn, {m,n} = 0,1,··· ,p−2. (36)
Interior bubble functions:
Theinterior bubblefunctions vanishon theentire boundary∂K2 of the reference 2-simplex K2.
The formulas of these functions are given as
λ −λ
Φt,~ei = C λ λ λ (1−λ )mP(2,2) 1 2 P(2m+5,2)(2λ −1)~e , i = 1,2, (37)
m,n m,n 0 1 2 0 m 1−λ n 0 i
(cid:18) 0 (cid:19)
where
(m+3)(m+4)(2m+5)(2m+n+6)(2m+n+7)(2m+2n+8)
C = ,
m,n
s (m+1)(m+2)(n+1)(n+2)
and
0 ≤ {m,n},m+n ≤ p−3.
One can again prove the orthonormal property of the interior bubble functions
< Φtm,~e1i,n1,Φtm,~e2j,n2 > |K2 = δm1m2δn1n2, (38)
where
0≤ {m ,m ,n ,n },m +n ,m +n ≤ p−3, {i,j} = 1,2.
1 2 1 2 1 1 2 2
The following table shows the decomposition of the space (P (K))2 for the H(div)-conforming
n
triangular element.
Decomposition Dimension
Edge functions 3(n+1)
Edge-based interior functions 3(n−1)
Interior bubble functions (n−2)(n−1)
Total (n+1)(n+2)
Table 1: Decomposition of the space (P (K))2 for the H(div)-conforming triangular element.
n
8
5 Basis functions for the tetrahedral element
Our constructions are motivated by the work on the construction of H(div)-conforming hierarchi-
cal bases for tetrahedral elements [7]. We construct shape functions for the H(div)-conforming
tetrahedral element on the canonical reference 3-simplex. The shape functions are grouped into
several categories based upon their geometrical entities on the reference 3-simplex [7]. The basis
functions in each category are constructed so that they are orthonormal on the reference element.
Any point in the 3-simplex K3 is uniquely located in terms of the local coordinate system
(ξ,η,ζ). Thevertexes are numbered as v (0,0,0),v (1,0,0),v (0,1,0),v (0,0,1). The barycentric
0 1 2 3
coordinates are given as
λ := 1−ξ−η−ζ, λ := ξ, λ := η, λ := ζ. (39)
0 1 2 3
The directed tangent on a generic edge e =[j ,j ] is defined as
j 1 2
τej := τ[j1,j2] = v −v , j < j . (40)
j2 j1 1 2
The edge is parameterized as
γej := λj2 −λj1, j1 < j2. (41)
A generic edge can be uniquely identified with
e := [j ,j ], j = 0,1,2, j <j ≤ 3, j = j +j +sign(j ), (42)
j 1 2 1 1 2 1 2 1
where sign(0) = 0. Each face on the 3-simplex can be identified by the associated three vertexes,
and is uniquely defined as
f := [j ,j ,j ], 0 ≤ {j ,j ,j ,j } ≤ 3, j < j < j . (43)
j1 2 3 4 1 2 3 4 2 3 4
The standard bases in Rn are noted as ~e , i= 1,··· ,n, and n = {2,3}.
i
5.1 Face functions
The face functions are further grouped into two categories: edge-based face functions and face
bubble functions.
Edge-based face functions:
These functions are associated with the three edges of a certain face f , and by construction
j1
all have non-zero normal components only on the associated face f , viz.
j1
nfjk ·Φefj[1k,1i,k2] = 0, jk 6= j1, (44)
f
where njk is the unit outward normal vector to face fj .
k
The edge-based face functions for higher order have been proposed in [7] as follows:
Φefj[1k,1i,k2] = li(γek)λk1∇λk2 ×∇λk3, i= 0,··· ,p−1. (45)
e
9
For instance, for the face opposite to the vertex v (0,0,0), f := [1,2,3], the face functions related
0 0
to edge e[1,2] are given by
Φf0,i = l (λ −λ )λ ∇λ ×∇λ , i = 0,··· ,p−1. (46)
e[1,2] i 3 2 1 2 3
However, itcanbeeasilycheckedthatthebasisfunctiongivenin(45)infactarenotindependent
e
for p = 2 (for example, the sum of all the basis on a given face in fact equals to zero) and thus the
proposed basis function is not complete. To remedy this degeneracy, two types of constructions of
new hierarchical high-order independent edge-based face functions will be presented here.
• Type One: high-order independent edge-based face functions
In [9], the following orthonormal basis functions are given as
Φfj1,i = C λ (1−λ )iP(3,0) 2λk2 −1 ∇λk1 ×∇λk2 , (47)
e[k1,k2] i k3 k1 i 1−λ |∇λ ×∇λ |
(cid:18) k1 (cid:19) k1 k2
where
C = 3(2i+4)(2i+5), i= 0,1,··· ,p−1,
i
and p
k ={j ,j }, k = {j ,j }, k < k , k = {j ,j ,j }\{k ,k }.
1 2 3 2 3 4 1 2 3 2 3 4 1 2
One can prove the orthonormal property of these edge-based face functions
< Φefj[1k,1m,k2],Φefj[1k,1n,k2] > |K3 = δmn, {m,n} = 0,1,··· ,p−1, (48)
where δ is the Kronecker delta. Note that with our new construction, the edge-based face
mn
functions are all linearly independent, which is also verified by the fact that in the spectrum of the
mass (Gram) matrix, none of the eigenvalues is zero.
• Type Two: high-order independent edge-based face functions
An alternative approach using the idea of recursion from [7] can also be used to construct
independent edge-based face functions as follows.
For p = 1, for each edge we have one face function for this edge as proposed in [7]
Φfj1,0 = λ ∇λ ×∇λ , (49)
e[k1,k2] k1 k2 k3
and for p = 2, one additional new basis function can be constructed as
e
Φfj1,1 = λ λ ∇λ ×∇λ , (50)
e[k1,k2] k1 k2 k3 k1
which can be shown to satisfy the condition (44), and and for p ≥ 3, the basis functions are given
e
by
Φefj[1k,1i,+k21] ≡ ℓi(γek)Φefj[1k,11,k2]+ℓi−1(γek)Φefj[1k,10,k2] (51)
e = ℓi(γek)e[λk1λk2∇λk3 ×∇λk1e]+ℓi−1(γek)[λk1∇λk2 ×∇λk3], i = 1,··· ,p−2.
10