Mixed finite element approximation of incompressible MHD problems based on weighted regularization.
We introduce and analyze a new mixed finite element method for the numerical approximation of stationary incompressible magneto-hydrodynamics (MHD) problems in polygonal and polyhedral domains. The method is based on standard inf-sup stable elements for the discretization of the hydrodynamic unknowns and on nodal elements for the discretization of the magnetic variables. In order to achieve convergence in non-convex domains, the magnetic bilinear form is suitably modified using the weighted regularization technique recently developed in [Numer. Math. 93 (2002) 239]. We first discuss the well-posedness of this approach and establish a novel existence and uniqueness result for non-linear MHD problems with small data. We then derive quasi-optimal error bounds for the proposed finite element method and show the convergence of the approximate solutions in non-convex domains. The theoretical results are confirmed in a series of numerical experiments for a linear two-dimensional Oseen-type MHD problem, demonstrating that weighted regularization is indispensable for the resolution of the strongest magnetic singularities.
Introduction
Incompressible magneto-hydrodynamics (MHD) describes the flow of a vis-
cous, incompressible and electrically conducting fluid. The governing partial
di#erential equations are obtained by coupling the incompressible Navier-Stokes
equations with Maxwell's equations and arise in several engineering
applications such as, for example, liquid metals in magnetic pumps or aluminum
electrolysis, see, e.g., [20]. In the stationary case, the resulting multi-field
problem is of the form: find the velocity field u(x), the hydrodynamic
pressure and the magnetic field
in# ,
in# ,
div
in# ,
in# ,
supplemented with suitable boundary conditions on #
Here,# is a bounded
domain in R 3 , and the functions f and g are given source terms, with g being
solenoidal. Furthermore, R e is the hydrodynamic Reynolds number, Rm the
magnetic Reynolds number, and S the coupling number. These numbers are
defined by
denoting the characteristic values of the magnetic field and the
velocity, respectively. The parameter L is the characteristic length scale of the
problem. The constants # and # represent the density and the viscosity of the
fluid, and - and # are the magnetic permeability and the electric conductivity,
respectively. In industrial applications, one typically has R e
Over the last few years, several finite element methods (FEM) to numerically
solve the incompressible MHD equations and linearizations thereof have been
proposed that are based on nodal (i.e., H
1(#79260428# finite elements for
the magnetic field b, combined with standard discretizations of the hydrodynamic
unknowns u and p. We mention here [1,13,16-18] and the references
cited therein. However, it has been known for some time that in non-convex
polyhedra of engineering practice the magnetic field components may have
regularity below H
1(# and that nodal FEM discretizations of the magnetic
operator, albeit stable, can converge quasi-optimally to a magnetic field that
certain singular solution components induced by reentrant vertices or
edges (for more details, see, e.g., [7] and the references cited therein). Con-
sequently, in non-convex domains, setting the magnetic components of the
incompressible MHD equations in H
1(# leads to a well-posed problem where
the magnetic field lacks certain singular (but physical) solution components.
A possible way to overcome these di#culties was recently proposed in [24,25]
by the use of N-ed-elec's elements for the magnetic field b and by the introduction
of an additional Lagrange multiplier related to the constraint div
In this paper, we propose a new mixed finite element approximation for incompressible
MHD problems. Our method is also based on nodal elements for
the magnetic field b, and employs standard inf-sup stable elements for the unknowns
u and p. However, as opposed to the approaches mentioned above, we
modify the magnetic bilinear form using the weighted regularization technique
recently developed by Costabel and Dauge in [7]. This allows us to account
for the possible low regularity of the magnetic field in non-convex domains.
We first discuss the well-posedness of this approach and show the existence
and uniqueness of weak solutions for small data. We then carry out an error
analysis for the proposed finite element method and show that it leads to
quasi-optimal error bounds in the mesh-size. Finally, we show the convergence
of the approximate solutions in non-convex domains where the components of
the magnetic fields may have regularity below H
1(#4 Our theoretical results
are confirmed in a series of numerical experiments for a linear Oseen-type
MHD problem in two dimensions.
The outline of the paper is as follows. In Section 2, we introduce a weighted
regularization approach for incompressible MHD problems and show the well-posedness
of the underlying weak formulation. Our finite element approximation
is proposed and analyzed in Section 3. A series of numerical results
for a two-dimensional MHD problem is presented in Section 4. We end our
presentation with concluding remarks in Section 5.
Throughout the paper, we use the following notation: For a Lipschitz domain
we denote by L p (D), 1 # p #, the Lebesgue space of
p-integrable functions, endowed with the norm # L p (D) . We write L p
loc(# to
denote the space of functions that are locally p-integrable. We further make
use of the subspace L 2
0 (D) of L 2 (D) defined by
For s # 0, we denote by H s (D) the standard L 2 -based Sobolev space of order s
and write # H s (D) for its norm. The closure of
(smooth functions with
compact support) in H s (D) is denoted by H s
(D). We write H -s (D) for the
dual space of H s
(D), equipped with the dual norm # H -s (D) . For a generic
function space X(D) we write X(D) n , to denote vector fields whose
components belong to X(D). Without further specification, these spaces are
equipped with the usual product norms which we simply denote in the same
way as the norms in X(D).
Weighted Regularization of Incompressible MHD Problems
In this section, we introduce the governing equations of stationary incompressible
magneto-hydrodynamics, derive a weak formulation using the weighted
regularization technique of [7], and establish the existence and uniqueness of
solutions for small data.
2.1 Incompressible MHD Equations
R 3 be a bounded Lipschitz polyhedron. We assume throughout
is simply-connected and that its boundary
# is connected. We consider stationary
incompressible MHD problems of the following form: Given forcing
terms f and g in L
find the velocity field
field and the pressure p such that
in# , (2.2)
div
# . (2.6)
Here, n is the outward normal unit vector on #
# For simplicity, we have
imposed no-slip boundary conditions on the velocity field u and perfectly insulating
magnetic boundary conditions on b. We comment on inhomogeneous
boundary conditions in Remark 2.19 below.
By taking the divergence of equation (2.2), we see that the datum g has to be
solenoidal. Thus, we assume throughout that g # H(div; # and
div
in#
Here,
| div g # L
endowed with the norm
.
Due to [14, Theorem I.3.4], there is a stream function # H
3 such that
Remark 2.1 We also consider the two-dimensional analogue of the MHD
problem (2.1)-(2.6). However, in two dimensions, the definition of the curl-
operators requires some care. For vector fields
we define the vector product w- and the scalar-valued curl-
operator curl
Furthermore, for a scalar function c, we set
vector-valued curl-operator is given by curl
c , -# x 1
c ). The two-dimensional analogue of (2.1)-(2.6) then reads as
follows: Find the velocity field
the pressure p such that
in# , (2.8)
div
# . (2.13)
Here, we denote by t the counterclockwise oriented unit tangent vector on
# .
Note that, by identifying the two-dimensional vector field
extension
is easy to see that curl curl
b.
Similarly,
b) and b - curl
b.
2.2 Weighted Spaces
To derive a weak formulation for (2.1)-(2.6) based on weighted regularization,
we need to introduce the weighted Sobolev spaces from [7].
To this end, we denote by C the set of all corners of the
polyhedron# For
At each corner c there is a ball B(c, R c )
of radius R c such that V c
=# B(c, R c ) is a cone which, in local spherical
coordinates, is of the form V
use the function r c as the radial coordinate while # c is the angular coordinate
with values in G c , a spherical polygonal domain in the unit sphere S 2 .
Furthermore, let E denote the set of all (open) edges of # Then, for each
point x of an edge e, there is a ball B(x, R x ) of radius R x such that V e
# B(x, R x ) is di#eomorphic to a wedge # e (x) -R with # e (x) being a plane
sector with opening angle # e # (0, 2#). This angle is intrinsic and is called the
opening angle
of# at the edge e. We set r e is equivalent to
the radial coordinate in # e (x). Let c and c # be the two endpoints of an edge e.
Then we define # e by r e (x). Note that # e is equivalent to
r e /r c in V c , to r e /r c # in V c # , and to r e outside of V c # V c # .
Next, we introduce weight vectors # which are collections of real numbers
{# c } c#C # e } e#E . We set For two weight vectors
# and #, we use the notation # to mean that # c # c and # e # e
for all c # C and e # E . Similarly, # is the weight vector given by the
components # c -# c and # e -# e . For constants # 1 and # 2 , we
to mean that # 1
With a weight vector # we associate the weight function
r #c
c
r #e
e
# . (2.14)
Moreover, we need the distance function E). For s # N 0 , we
define the weighted space
loc(# d |#|-s # v # L
0 with |# s}, (2.15)
equipped with the norm
|#s
.
A special role is played by the weight vector # dir given by
# dir
# dir
Here, following [7, Section 4.3], we have set
# dir
with - dir
denoting the smallest Dirichlet eigenvalue of the Laplace-Beltrami
operator in the cone G c . Note that we always have # dir
there holds
|# dir
| <2 - # dir
for a parameter # dir
We remark that # dir
# approaches zero if one of the
opening angles
of# approaches 2#.
The following result holds; see [7, Theorem 4.1].
Theorem 2.2 Let # dir be the weight vector in (2.16)-(2.17). Then for any
weight vector # with
the Laplacian with Dirichlet boundary conditions is an isomorphism from
Remark 2.3 In a Lipschitz
polygon# R 2 the definitions of # and the
corresponding weighted spaces are easier as the edges do not need to be taken
into account. To define # , let C be the set of all corners
of# . The opening
angle at the corner c # C is denoted by # c , with # c # (0, 2#). Then, for a
weight vector # c } c#C , the weight function # is given by
r c (x) #c , r c
Introducing the distance function d(x) by C), the spaces V s
are defined as in (2.15). In two dimensions, the critical weight vector # dir is
# dir
The result of Theorem 2.2 then holds true, with # dir given in (2.21).
Next, we introduce the Sobolev spaces
| curl b # L
as well as
H(curl;# | n -
and endow them with the norm
For a weight vector #, we further define the space
(curl;# | div b # V 0
and equip it with the norm
X# curl b# 2
Our finite element discretization will be based on the subspace H
#givenby
equipped with the norm # X# .
The following result is crucial; see [7, Theorem 5.1].
Theorem 2.4 Let # be a weight vector satisfying (2.19). Then the space
# is dense in X
Remark 2.5 The same result holds true in
defined as in (2.20) and (2.21), respectively.
Finally, let us show a Poincar-e-type inequality that will be needed in our
analysis. We set
X# curl b# 2
and have the following result.
Proposition 2.6 Let # be a weight vector satisfying (2.19). There holds:
(i) | - | X# is a norm on X
# .
(ii) There exists a constant C > 0 only depending
on# and # such that
|b| X# C#b#X#b # X
# .
To prove the first assertion, it is su#cient to show that |b|
we conclude that curl
in # we also have div using that n -
and the Poincar-e-type inequality from [12, Proposition 7.4], we obtain
with a constant C > 0 only depending on # Thus, which shows the
first assertion.
For the second assertion let # H 1
0(# be the solution of
# and the Laplacian is an isomorphism from
#1 see Theorem 2.2, we have # V 2
# and
for a constant C > 0 solely depending
on# and #.
By setting b #, we have div b in # as
well as n - b
As before, the inequality in [12, Proposition 7.4]
yields
for a constant C > 0 only depending on #
Referring to (2.25) and (2.26) gives
for any b # X
#7 with a constant C > 0 solely depending
on# and #.
Equation (2.27) implies the second assertion. #
We point out that the results in Proposition 2.6 can be easily adapted to the
two-dimensional case.
2.3 Weak Formulation
For a weight vector 0 # 1, we define the following weak form for the MHD
problem
# and p # L 2
0(# such that
a s (u, v)
am (b, c) - c 2 (b; u, c) = (g, c), (2.28)
for all v # H 1
# and q # L 2
0(#2 Here, we use the following
forms:
a s (u,
am (b, c)
Here, we have incorporated the regularization term
the magnetic form am , following [7]. The parameter D is a positive constant
that can be used to dimensionalize the regularization term and to balance
it with the curl-curl term. Furthermore, we use the standard anti-symmetric
for the discretization of the non-linear convection term in
the Navier-Stokes operator; see, e.g., [27, Chapter II] for details. The trilinear
forms c 1 and c 2 arise due to the coupling terms in (2.1) and (2.2); we show in
Section 2.4 below that these forms are well-defined for suitable choices of #.
The regularization term ensures that the magnetic field is solenoidal.
Proposition 2.7 Assume (2.19) and let (u, b, p) be a solution of (2.28). Then
we have div
#1 the problem
# is
well-posed and, as before, has a unique weak solution #
0(#4 By
Choosing as a test function in the second
equation of formulation (2.28) yields
0,# .
Here, we have used that curl Furthermore, by integration by parts,
since
# and div assumed in (2.7). Thus, since # (x) > 0
in# and D > 0, we have div
For the purpose of our analysis, we rewrite the formulation (2.28) in the compact
Find (u, b, p) # H 1
0(# such that
for all (v, c, q) # H 1
0(#9 Here, we use the forms
O(w, d; u, b; v, c)
The adaptation of the forms in (2.29)-(2.34) and the weak formulation (2.28)
to two-dimensional MHD problems of the form (2.8)-(2.13) is straightforward.
2.4 Well-Posedness
We show that the variational formulation (2.28) is well-posed and uniquely
solvable for small data. We begin by establishing the continuity of the forms.
Lemma 2.8 For any weight vector #, the forms a s , am , and b s satisfy the
following continuity properties:
3 .
Furthermore, there exists a constant C depending
on# such that
3 .
Proof : The continuity properties of the forms a s , am and b s follow straight-forwardly
using Cauchy-Schwarz inequalities and the fact that # div v# L
. The continuity property for follows from the continuous embedding
4(# and H-older's inequality; see, e.g., [14, Chapter IV].#
To show the continuity of the forms c 1 and c 2 , we make use of the following
embedding result, which is valid for suitable values of # (see also Remark 2.11
below). For each corner c # C, we for the set of all edges that contain
the corner c.
Lemma 2.9 Let the weight vector # satisfy
Then we have H
-# and #v# 0,-# C#v# H |#| for a constant
only depending
on# and #.
As described in [7, Section 4.1], we can
e
e#Ec
Here, V 0 is a subregion
of# away from corners and edges, and V 0
e is a subregion
of# such that V 0
e does not contain any corners or parts of any other edge
than e. The subregion V 0
c is such that c # V 0
c and e # V 0
edge
Finally, for any edge e # E c the subregion V e (c) is such that V e (c) only
contains c and parts of e. Note that this decomposition is not unique and that
the di#erent subregions may be overlapping.
As in [7, Equation (4.9)], we define the distance functions
d C E).
Furthermore, we can choose exponents # C and # E such that |# C # 0,
c#e
The weight w-# is then equivalent to
Let # be a weight vector satisfying (2.36) and let v be in H
|#1 We may
assume that 1/2 > |#| > 0, the case being trivial. We obtain
Here, we have used that # e # c for all c # E c , dist(x, #
the continuous embedding in [15, Theorem 1.4.4.3] ensures
that
which completes the proof. #
Lemma 2.10 Let the weight vector # satisfy (2.36). Then there is a constant
depending
on# and # such that
3 .
We start by noting that, from [7, Theorem 2.2] any field d # X
#canbedecomposedas
with d 0 # H
0(#7 Furthermore, there exists a constant
depending
on# and #, such that
Let us now establish the assertion for c 1 ; the proof for c 2 is completely analo-
gous. We write c 1 (d; b, v) as
We first bound c 1 (d inequality, the continuous embedding
of H
into L
4(#7 and the estimate (2.39), we have
with a constant C > 0 only depending
on# and #.
Next, we bound c 1 (#; b, v). To do so, we first note that, since |#| < 1/2,
we have H
0(#0 see [15, Corollary 1.4.4.5]. Thus, from Lemma 2.9
and duality we have that
with a constant only depending
on# and #. Hence, since # V 0
In view of |#| < 1/2, the elliptic shift theorem for polyhedral domains implies
that
1/2+# and # H
for a parameter # > 0; see [8]. Further, we have that H
1/2+# is continuously
embedded into
q(# for an exponent q > 3; see [14, Theorem I.1.3 and
Definition I.1.2]. We can then find a second exponent p < 6 such that
Using H-older's inequality with these exponents and the continuous
embedding of H
into L
p(#3 we obtain
.
This, together with (2.39)-(2.41), proves the assertion for c 1 . #
Remark 2.11 The assumptions in (2.36) restrict the choice of # to quite a
small range. In view of (2.18), this is particularly evident when we simultaneously
seek to fulfill (2.36) and (2.19), that is,
see also Theorem 2.17 below. The upper bound 1/2 shows up because of the use
of the embedding result in Lemma 2.9. Whether or not this upper bound can
be improved with a di#erent analysis technique remains an open question. We
also point out that for linear MHD problems no restrictions on # are necessary.
Next, let us show that for the two-dimensional analogues of the forms c 1 and c 2
in (2.33)-(2.34) it is possible to obtain a result with less restrictions on #.
Lemma 2.12
Let# be a polygon in R 2 and let the two-dimensional weight
vector # in (2.20) satisfy 0 # < 1. Then there is a constant C c depending
on# and # such that the two-dimensional analogues of the forms c 1 and c 2
.
As in the proof of Lemma 2.10, we can write
0(#5 The contribution c 1 (d
as in Lemma 2.10. To bound c 1 (#; b, v), we proceed as follows. First note
that
2 and # 1,# C# 2,#. From [26, Proposition 25],
we obtain
for a constant C > 0 solely depending
on# and #. Let p and q be parameters
inequality, Rellich's
embedding theorem and the above estimate, we obtain
|
with a constant C > 0 solely depending on #, and the exponents p and q.
It remains to show that q can be chosen so that # -1 # v# L
To this end, let s and s # be two other parameters with 1 < s # s < # and
We have
c be a small neighborhood of the corner c # C. In local polar coordinates
(r c , #) at the point c, there holds
|# | -qs # dx # C # Vc
r -qs #c+1
c dr c d#,
provided that qs # c < 2. The constant C only depends on # Since we have
there is a parameter # > 0 (depending on #) such that the
condition qs # c < 2 is fulfilled for #. With this choice,
we obtain
with a constant C > 0 depending
on# and #. Combining (2.42) and (2.43),
using Rellich's embedding theorem and the estimate in (2.39), results in
for a constant C > 0 depending
on# and #.
This yields the result for c 1 , the proof for c 2 is analogous. #
Remark 2.13 We point out that the proof in Lemma 2.12 is based on the
continuous embedding of H
into L
q(# for all q # 1. Since in the three-dimensional
1(# is continuously embedded into
q(# only for q # [1, 6], a
similar argument in three dimensions shows the continuity of c 1 and c 2 only
for polyhedral domains whose maximal opening angle is smaller than 2#/3.
Next, we address the coercivity of the forms a s and am .
Lemma 2.14 Let # be a weight vector satisfying (2.19). Then:
a s (u, u) # C
with a constant C 1 > 0 only depending
on# , and a constant C 2 > 0 only
depending
on# and #.
Proof : The coercivity of a s is a standard property and the coercivity of am
follows from the inequality in Proposition 2.6. #
Finally, we recall the following inf-sup condition for the form b s ; see, e.g., [14,
Section I.5.1].
Lemma 2.15 There is a constant # > 0, only depending
on# , such that
q#L 2sup
#.
For notational convenience, we introduce the space
and endow it with the norm
X# . (2.44)
We then have the following stability results for the forms in (2.35).
Proposition 2.16 Let # be a weight vector satisfying (2.19). There holds:
(i) There are continuity constants CA , CB solely depending on the data and
on D such that
# .
If we additionally assume (2.36) to hold, then there is a constant C C
depending on the data, the domain, and the weight # such that
|O(w, d; u, b; v, c)| # C C #(w, d)# W#(u, b)# W#(v, c)# W# ,
for any (w, d) # W
# , and (v, c) # W
# .
(ii) There is a coercivity constant # > 0, depending on the data, the parameter
D, the domain, and #, such that
# .
We have
# .
(iv) We have the skew-symmetry property
O(w, d; u, b; u,
# .
(v) There holds, for the same constant # as in Lemma 2.15 and independently
of #,
sup
(v,c)#W#
B(q; v, c)
#(v, c)# W#q# L
#.
Proof : The continuity and coercivity in (i) and (ii) are immediate consequences
of Lemma 2.8, Lemma 2.10 and Lemma 2.14, respectively. The continuity
property in (iii) holds since
To see the skew-symmetry property of O in (iv), it is enough to note that
b). The latter identity follows since
d) curl b.
It only remains to establish the inf-sup condition in (v) for the form B. To see
this, fix q # L 2
0(#6 From Lemma 2.15, there is an element v # H 1
3 such
that
.
We obtain
2(#(v, 0)# X#q# L
and the inf-sup condition for B follows. #
Proceeding as in the proof of [23, Theorem 10.1.1], we obtain from Proposition
2.16 the following existence and uniqueness result for small data.
Theorem 2.17 Let # be a weight vector satisfying (2.19) and (2.36). Assume
further that
< 1. (2.45)
Then the weak formulation in (2.28) has a unique solution (u, b, p) # H 1
0(# and we have the stability bounds
#p# L
with #, CA , CB , C C , and L denoting the stability constants from Proposition
2.16.
Remark 2.18 As has been pointed out in Remark 2.11, the restrictions on #
in (2.36) are most likely suboptimal. For the two-dimensional MHD problem
in (2.8)-(2.13), on the other hand, the continuity result in Lemma 2.12 can be
invoked and the result of Theorem 2.17 is obtained for any weight # satisfying
# dir < # < 1, with # dir defined in (2.21).
Remark 2.19 The extension of the result in Theorem 2.17 to MHD problems
with inhomogeneous boundary conditions is not straightforward. While
it is easily possible to lift inhomogeneous boundary data in a divergence-free
fashion into the domain, see [14, Lemma IV.2.3] for the velocity field and [22,
Proposition A.1] for the magnetic field, these liftings a#ect the size of the
data for which existence and uniqueness of solutions can be proved; see [14,
Section IV.2.1]. For velocity boundary data, this e#ect can be minimized and
controlled by using the so-called Hopf construction which yields a divergence-free
lifting
3 such that |
arbitrarily small
relative to #v# 2
analogous Hopf-type liftings for magnetic boundary data seem not to be available
in the literature and remain to be constructed.
3 Finite Element Approximation
In this section, we introduce and analyze the finite element approximation
of the mixed formulation in (2.28). We derive quasi-optimal error bounds in
the energy norm and show that the weighted regularization technique ensures
convergence of the approximation in possibly non-convex domains.
3.1 Galerkin Approximation
We choose conforming finite element spaces V h
#4 and
0(#3 and endow them with the norms # H
respectively. Here, we use the index h to denote the discretization parameter.
We generically refer to it as the mesh-size.
Throughout, we assume that the pair V h
rise to an inf-sup stable
Stokes discretization, that is, we assume that there is a constant # h > 0
independent of the mesh-size h, such that
sup
# h . (3.1)
A wide variety of spaces V h and L h fulfilling (3.1) have been proposed in
the literature; we refer to [5, Chapter IV], [14, Chapter II] and the references
cited therein. A specific choice of finite element spaces based on Hood-Taylor
elements will be discussed in Section 3.2 below.
Given a weight vector #, the finite element approximation of (2.28) is: Find
h such that
a s (u h , v)
am (b h , c) - c 2 (b
for all
# and q # L h . As before, problem (3.2) is equivalent to:
Find
h such that
for all (v, c,
Introducing the space W h
endowed with the norm # W# , we
have the following discrete inf-sup condition for the form B
sup
(v,c)#W h
B(q; v, c)
#(v, c)# W#q# L 2
with the same inf-sup constant # h > 0 as in (3.1). Condition (3.4) can be
proved by using arguments that are completely analogous to those on the
continuous level.
The discrete version of Theorem 2.17 is then an immediate consequence.
Corollary 3.1 Let # be a weight vector satisfying (2.19) and (2.36). Let the
smallness assumption (2.45) be satisfied. Then the finite element formulation
in (3.2) has a unique solution
we have the
stability bounds
with #, CA , CB , C C , and L denoting the stability constants from Proposition
2.16, and with # h denoting the discrete inf-sup constant from (3.1).
As before, in the two-dimensional case, this result holds for weight vectors #
with # dir < # < 1.
Remark 3.2 The solution
of the finite element
formulation (3.2) can be found by the following Picard iteration: Given
h be the
solution of the linearized problem
a s (u n+1
am (b n+1
for all
# and q # L h . Under the smallness assumption in (2.45),
the sequence {(u n
converges to the solution (u h , b h , p h ) of (3.2).
Other procedures based on Newton's method are possible as well; cf. [18].
We point out that, if the linearized problems above are strongly convection-
dominated, it might be necessary for their e#cient solution to include additional
stabilization terms along the lines of [13]. As our analysis is mainly
concerned with the incorporation of the divergence constraint div via the
weighted regularization approach, this point is not further investigated in this
paper.
We derive quasi-optimal error bounds for the proposed finite element approx-
imation. To this end, we introduce on W
0(# the norm ||| (-) ||| #
given by
.
The following theorem holds.
Theorem 3.3 Let # be a weight vector satisfying (2.19) and (2.36). Assume
further that
. (3.5)
Let (u, b, p) be the (unique) solution of (2.28), and let (u h , b h , p h ) its finite
element approximition obtained by (3.2). Then we have the quasi-optimal error
bound
(v,c,q)#V h
with a constant C > 0 independent of the mesh-size h.
We proceed in several steps.
Step 1: We first note that we have the error equation
for any (v, c) # W h
# .
Step 2: Set ker B
# | B(q; v,
We claim that
for any (v, c) # ker B h and q # L h , with a constant C > 0 that is independent
of the mesh-size.
To see (3.6), fix (v, c) # ker B h and q # L h . Clearly, v - u h # ker B h . Using
the error equation in Step 1, it can be easily seen that
We first estimate the left-hand side of (3.7) from below. To this end, we use the
coercivity and continuity properties in Proposition 2.16, the stability bound
in Theorem 2.17, and the smallness assumption in (3.5), and obtain
l.h.s. of
W#
W#
W# .
To bound the right-hand side of (3.7) from above, we first note that, because
Using the continuity properties in Proposition 2.16 and the bounds in Theorem
2.17 and Corollary 3.1 we get
r.h.s. of (3.7) #
Combining (3.8) and (3.9) results in
.
Since
the assertion (3.6) follows.
Step 3: It is well-known that we can use the discrete inf-sup condition in (3.4)
in order to establish the approximation result (3.6) for any (v, c) # W h
e.g., [5,14]. This proves the quasi-optimality of the error #(u-u h , b-b h )# W# .
Step 4: It remains to bound the error in the pressure. To do so, let q # L h .
The inf-sup condition (3.4) yields
(v,c)#W h
#(v, c)# W#
(v,c)#W h
#(v, c)# W#
(v,c)#W h
#(v, c)# W#
(v,c)#W h
#v, c# W#
Using the error equation from Step 1, the continuity properties in Property
2.16, and the stability bounds in Theorem 2.17 and Corollary 3.1, we
have
#(v, c)# W#
Thus,
#(v, c)# W#
Moreover, employing the continuity of B,
#(v, c)# W#
.
We obtain
# .
The estimate and the assertion follows then with the bounds of Step 3. #
Let now {V h
} h be sequences of finite element spaces,
chosen such that
Note that due to Theorem 2.4 and the density of C # functions with vanishing
trace on
0(#5 see [7, Proposition 3.2 and Theorem 4.1],
the density assumption for b is justified, provided that the weight vector #
satisfies (2.19).
Corollary 3.4 Assume (3.5) and that the weight vector # satisfies
and (2.36). Then, we have for the above sequence of spaces
lim
In two dimensions, the same result holds for weights # satisfying # dir < # < 1,
with # dir given in (2.21).
Theorem 3.3 ensures convergence of the finite element approximation in non-convex
polyhedra as h # 0, provided that the weight vector is properly chosen.
The choice # weighted regulariza-
tion), for example, does not lead to convergent FEM solutions in non-convex
polygons. This is due to the fact that, without weighted regularization, the
space H
# is known to be a closed subspace of X
# and the strongest
magnetic singularities lie in the complement X
#7 see [7]. Hence, it
is impossible to correctly capture the magnetic fields. This behavior is clearly
confirmed in our numerical results in Section 4, demonstrating that weighted
regularization is indispensable in non-convex domains.
3.2 Convergence Rates for Hood-Taylor Elements
In this section, we present a specific finite element family based on Hood-
Taylor elements for the unknowns u and p, and discuss the corresponding
convergence rates.
To this end, let T {K} be a regular and shape-regular partition
of#
into hexahedral elements {K}. We assume that each element K is a#nely
equivalent to the reference cube -
. We denote by hK the diameter
of element K and set For an approximation order k # 2,
we introduce the following finite element spaces
denotes the space of polynomials of degree # k in each variable
on K. The velocity-pressure pair V h
is referred to as "Hood-Taylor"
elements. It is well-known that the spaces V h and L h satisfy the discrete inf-sup
condition in (3.1); see [4].
For this family, let us discuss the convergence rates that can be expected
from Theorem 3.3. We first consider the case of a smooth solution (u, b, p).
Standard approximation properties then give straightforwardly the following
optimal convergence rates.
Corollary 3.5 Let the exact solution (u, b, p) of (2.1)-(2.6) satisfy
-H
.
Under the assumptions of Theorem 3.3, there holds
with a constant C > 0 independent of the mesh-size h.
Next, we show that positive convergence rates are still possible for solutions
that exhibit singularities at the corners of the
domain# To this end, we
consider a model situation where we assume that the exact solution (u, b, p)
can be decomposed into a regular and a singular part, according to
sing , b sing , p sing ).
Such decompositions can be found in, e.g., [9,21,19,6,7] in the context of
the Navier-Stokes equations, linearizations thereof, and Maxwell's equations.
Analogous results then hold for linear MHD problems, as the one considered
in Section 4 below. However, for the nonlinear MHD problems under consider-
ation, decompositions of the above type do not seem to be available in detail.
Thus, here we assume that the regular part is smooth and satisfies
.
It is then clear that there is an interpolant (v reg , c reg , q reg
that
The part (u sing , b sing , p sing ) consists of the singular functions. These functions
have a low global regularity, but are typically smooth in the interior of the
domain. This behavior can be described best in terms of the following limits
of weighted spaces
m#N
# .
We then assume that the velocity-pressure singularities belong to
sing
sing
# dir-1 .
Note that u sing
sing
arguments similar to those
in, e. g., [26] and the references therein, an interpolant (v sing , q sing
can be constructed such that
sing
sing
sing
sing
for an exponent # 1 > 0 depending on # dir .
Concerning the magnetic field b sing , we assume, in agreement with [6,7], that
it consists of Neumann singularities of the Laplace operator and singularities
that are gradients of Dirichlet singularities of the Laplacian. That is,
.
Here, similarly to (2.16)-(2.17), # neu are the minimal singularity exponents
for the Laplacian with Neumann boundary conditions; see [7, Section 6]. Note
that b sing
As in [7, Section 7 and Section 8], it is possible to construct an approximation
c sing
# such that
#b sing
sing
for a parameter # 2 # (0, 1), provided that we have k # k 0 for a su#ciently large
threshold value k 0 . The restrictions on the polynomial degree k are due to the
proof in [7, Section 7] where it is necessary to construct H 2 -conforming interpolants
of the gradient components of the magnetic singularities. However, our
numerical results in two dimensions indicate that positive convergence rates
are already achieved for k # 2 on rectangular grids.
In the model situation described above, we have the following result.
Corollary 3.6 Under the above assumptions and those in Theorem 3.3, there
holds
The constant C > 0 is independent of the mesh-size h.
Similar convergence results are obtained in the two-dimensional case.
In this section, we present several numerical experiments for the linear two-dimensional
Oseen-type MHD problem
in# , (4.2)
div
in# , (4.3)
in# , (4.4)
where d is a given smooth magnetic field and w a smooth flow field. Problems
of this type arise in each step of the Picard iteration in Remark 3.2. On non-convex
domains, problem (4.1)-(4.4) already exhibits magnetic singularities
with regularity below H
Hence, it is well suited to test the performance
of the proposed finite element method.
We approximate (4.1)-(4.4) using the two-dimensional analogue of the Taylor-
Hood family (3.10) on square meshes. Our implementation is based on the
finite element library deal.II; see [3,2]. It provides powerful C ++ classes for
handling meshes and degrees of freedom, and for solving the resulting linear
systems of equations. In our experiments, we have solved these systems by
BICGSTAB, using a simple Jacobi preconditioner. While this worked well in
our examples, we point out that the systematic design and analysis of e#-
cient solvers for the weighted regularization approach proposed in this paper
remain open issues. For comprehensive discussions of e#cient preconditioning
and solution techniques for incompressible Navier-Stokes discretizations in the
absence of electro-magnetic e#ects, we refer the reader to, e.g., [11,10,28] and
the references therein.
x
1Fig. 1. L-shaped
domain#
Throughout, we consider the L-shaped
domain# with opening angle 3#/2
shown in Figure 1. We always set R
prescribe the right-hand sides f , g, as well as the field d and w. Furthermore,
we allow for non-homogeneous Dirichlet boundary conditions for u and b - t on
the boundary
# of # These conditions are taken into account in the usual
fashion by interpolating the boundary data at the corresponding nodal degrees
of freedom. We always take choose the weight function # in the
bilinear form am as # that we are
varying in our experiments. We point out that, for the linear MHD problem
in (4.1)-(4.4), the theoretical results of the previous sections hold without any
restrictions on #.
4.1 Smooth Solution
In our first experiment, we validate the a priori error bounds in Corollary 3.5
for a smooth solution. We solve the problem (4.1)-(4.4) with
and with f , g and the boundary data chosen so that the exact solution
Note that u, b, and the corresponding right-hand side g are solenoidal.
We compute finite element approximations to this MHD solution using two-dimensional
Hood-Taylor elements on a sequence of successively
refined square meshes {T i } i#1 , referring to the index i as cycle i. The number
of elements in the mesh T i is proportional to 2 2i ; the mesh-size h i of T i is
thus proportional to 2 -i . If e i denotes the error in some component of the
approximation on cycle i (in a suitable norm), the corresponding numerical
rate of convergence is given by
log(e
log(h
In

Table

1 and Table 2, we show the errors in the indicated norms for the
hydrodynamic variables (u, p) and the magnetic field b, respectively, obtained
with exponents We also list the number of degrees
of freedom (dofs) for each of the solution components. For all choices of #,
the rates in the H 1 -error in u, the L 2 -error in p, and in the X # -error in b are
of order two, in full agreement with the results of Corollary 3.5. As can be
expected, the convergence rates in the L 2 -errors of u and b are of one order
higher and of optimal third order. The di#erence in the results with respect
to the di#erent values of # is minimal and almost negligible, indicating that
the weighted regularization term has no influence on the performance of the
proposed method if the solution is smooth.
cycle dofs in u/p L 2 -error in u H 1 -error in u L 2 -error in p

Table
Smooth solution: Errors and convergence rates in (u, p).
4.2 Non-Smooth Solution
Next, we consider the MHD problem (4.1)-(4.4) where all the solution components
have corner singularities at the origin. We set again
and choose the data so that the magnetic field b is given by the strongest
singularity of the curl-curl operator for the L-shaped domain in Figure 1,
namely
denoting the standard polar coordinates. Evidently, curl
We also point out that b # H
. The hydrodynamic pair (u, p)
is taken to be the strongest corner singularity of the Stokes operator for the
cycle dofs in b L 2 -error in b X # -error in b

Table
Smooth solution: Errors and convergence rates in b.
L-shaped domain in Figure 1. This singularity is given by
with
The exponent # is the smallest positive solution of sin(#3#/2)+#
0, which is # 0.54448373678246. Note that u is solenoidal, and that (u, p) #
-H
#6 We consider the same sequence of meshes as in Section 4.1.
The performance of the proposed method is shown in Table 3 and Table 4;
it now strongly depends on the choice of the exponent #. The choice
(no regularization) does not lead to convergent FEM solutions and the errors
in the field b do not decrease at all. This is due to the fact that for
the space H
# is known to be a closed subspace of X
# and that the
singular solution (4.6) lies in the complement X
#8 see [7]. Hence,
is impossible to correctly capture the singular magnetic field. The
tables clearly show convergence of the method for
the weighted regularization is switched on, hereby confirming the results of
Corollary 3.4 and Corollary 3.6. The convergence rates in the H 1 -norm for u
and the L 2 -norm for p are of the expected order #. The L 2 -norm in u converges
with twice that order. In contrast to the restrictions on the polynomial degree
in the theoretical results, good convergence rates for b are already obtained
for quadratic elements. For # close to one, this order is close to the order 2/3;
cf. the discussion in [7].
We point out that the errors are only slightly better with increased approximation
order k, and that the rates remain comparable. In all our tests we
started to observe convergence as soon as # > 1/3, the lower critical bound
in (2.21) that is required for convergence. However, the best results in the
were always obtained with the upper bound # = 1.
It is evident that the weighted regularization term is needed to correctly capture
the singular behavior of the solution.
cycle dofs in u/p L 2 -error in u H 1 -error in u L 2 -error in p

Table
Non-smooth solution: Errors and convergence rates in (u, p).
cycle dofs in b L 2 -error in b X # -error in b

Table
Non-smooth solution: Errors and convergence rates in b.
Conclusions
In this paper, we have introduced and analyzed a new finite element method
for incompressible MHD problems in polygonal and polyhedral domains. The
method employs nodal elements for the discretization of the magnetic fields
and standard inf-sup elements for the hydrodynamic variables. In order to account
for singular solution behavior, the magnetic bilinear form has been modified
using the weighted regularization technique recently developed in [7]. The
analysis in this paper shows that this approach leads to convergent schemes in
non-convex domains. Our two-dimensional numerical results on an L-shaped
domain confirm that the weighted regularization approach is indispensable
for the numerical resolution of singular solution components in the magnetic
fields.
We finally point out that this paper is mainly concerned with the discretization
of the elliptic operator underlying the MHD problems under consideration. For
strongly convection-dominated problems, additional stabilization techniques
might be necessary to ensure the robustness of the schemes; see, e.g., [13] and
the references therein.



--R



Concepts for object-oriented finite element software - the deal
Stability of higher-order Hood-Taylor methods
Mixed and hybrid finite element methods
Singularities of electromagnetic fields in polyhedral domains
Weighted regularization of Maxwell equations in polyhedral domains
Elliptic boundary values problems in corner domains
Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners
Iterative methods for problems in computational fluid dynamics
Performance and analysis of saddle point preconditioners for the discrete steady-state
Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions
A stabilized finite element method for the incompressible magnetohydrodynamic equations
Finite element methods for
Ellipitc problems in non-smooth domains
Mixed finite element approximation of an MHD problem involving conducting and insulating regions: the 2D case
Mixed finite element approximation of an MHD problem involving conducting and insulating regions: the 3D case
On the existence and uniqueness and finite element approximation of solutions of the equations of stationary incompressible magnetohydrodynamics
Analytic regularity of Stokes flow in polygonal domains
The electromagnetics of fluids
Regularity and FEM-error estimates for general boundary value problems of the

Numerical approximation of partial di
Mixed finite elements for incompressible magneto-hydrodynamics
Mixed finite element methods for stationary incompressible magneto-hydrodynamics
Exponential convergence in a Galerkin least squares hp-FEM for Stokes flow
Theory and numerical analysis

--TR
Stationary stokes and Navier-Stokes systems on two-or three-dimensional domains with corners. Part I: linearized equations
Stability of higher-order Hood-Taylor methods
Mixed and hybrid finite element methods
The<i>hp</i>-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations
Mixed finite element methods for stationary incompressible magnetoMYAMPERSAND#x2013;hydrodynamics
weighted regularization;mixed methods;incompressible magneto-hydrodynamics