A Nonconforming Finite Element Method for Fourth Order Curl Equations in R^3

In this paper we present a nonconforming finite element method for solving fourth order curl equations in three dimensions arising from magnetohydrodynamics models. We show that the method has an optimal error estimate for a model problem involving both curl^2 and curl^4 operators. The element has a very small number of degrees of freedom and it imposes the inter-element continuity along the tangential direction which is appropriate for the approximation of magnetic fields. We also provide explicit formulae of basis functions for this element.


Introduction
The magnetohydrodynamics (MHD) equations describe macroscopic dynamics of electrically conducting fluid that moves in a magnetic field. MHD model is governed by Navier-Stokes equations coupled with Maxwell equations through Ohm's law and Lorentz force. As an example, a resistive MHD system is described by the following equations: where ρ is the mass density, u is the velocity, p is the pressure, B is the magnetic induction field, η is the resistivity, η 2 is the hyper-resistivity, µ 0 is the magnetic permeability of free space, and µ is the viscosity. The primary variables in MHD equations are fluid velocity u and magnetic field B.
MHD model has widespread applications in thermonuclear fusion, magnetospheric and solar physics, plasma physics, geophysics, and astrophysics. Mathematical modeling and numerical simulations of MHD have attracted much research effort in the past few decades. Various numerical algorithms have been used in MHD simulations; examples include finite difference methods, finite volume methods, finite element methods, and Fourier-based spectral and pseudo-spectral methods [27]. In [12,13,14,15,21], two-dimensional, incompressible MHD problems are studied in terms of finite element approximations of the stream function-vorticity advection formulation. Since MHD flow often develop sharp interfaces, adaptive h-refinement techniques have been applied in MHD simulations [16,25,33]. Finite element computations of MHD problems in three-dimensions have been reported in [5,6,17,23,24,32].
In the existing finite element discreitzations for the above MHD model, a standard pair of stable or stabilized finite element spaces are often used to discretize the velocity and pressure variables in the fluid equations. For the magnetic field variable B, however, at least two approaches are possible when the fourth order term (∇×) 4 B is not presented in the model, namely when electron viscosity η 2 = 0. One approach is to use the standard edge element ( [24]) and the other approach is to use the Lagrange element after replacing (∇×) 2 B by −∆B ( [6,23,32]). Both these approaches will become more difficult when the fourth order term (∇×) 4 B is presented. We may still replace (∇×) 4 B by a biharmonic operator ∆ 2 B. But discretizing a biharmonic operator in three dimensions is challenging. It requires 220 degrees of freedom per element if a conforming finite element is used. One possible way to reduce the number of degrees of freedom is to use nonconforming discretizations which allow weaker inter-element smoothness constraints but still provid convergent approximations. Among the class of nonconforming finite elements for fourth order problems, Morley-type elements are special in the sense that they provide approximations with polynomials of minimal degree [18,28]. In [29], a systematic construction of Morley-type elements is provided for solving 2m-th order partial differential equations in R n . In particular, we may apply the element in [29] with n = 3 and m = 2 consisting of piecewise quadratic elements to our system of biharmonic equations. This amounts to 30 degrees of freedom on each element. This element provides a reasonable discretization of MHD equations when it is appropriate to replace (curl) 4 operator by the biharmonic operator. This approach, however, may lead to difficulty for certain boundary conditions in practical applications. Indeed the treatment of boundary conditions is also an issue for the second order problem if (∇×) 2 B is replaced by −∆B [10].
One more natural approach is to discretize the fourth order curl operator by some generalized higher order edge elements. But such type of edge elements are not available in the literature. The construction of such type of edge element is the main goal of this paper.
Another possible approach to deal with the fourth order term (∇×) 4 B is to use operator splitting technique. Namely, one can introduce an intermediate variable σ = (∇×) 2 B and then reduce the original problem to a system of second order equations. However, it is known that for some problems, such a technique cannot be applied. For example, when modeling the bending of simply supported plate on non-convex polygonal domains, the original biharmonic problem is not equivalent to the lower order system of two Poisson equations [2,22]. In view of this, we consider discretizing the fourth order problem directly.
In this paper, we investigate MHD equations that contain both fourth-order term and second-order term. In the literature, the major tool used for performing MHD simulations involving a fourth order equation has been the pseudo-spectral method [1]. By choosing an appropriate formulation, we are able to construct a finite element approximation for this problem. This is a nonconforming finite element that involves only a small number of degrees of freedom.
The rest of this paper is organized as follows. In Section 2, we describe a simplified model problem and the corresponding variational formulation. In Section 3, we construct basis functions and provide the convergence analysis. Finally, in Section 4, we give some concluding remarks.

Model Problem
In the following, we introduce model problems for the fourth-order magnetic induction equations described above. Assume that Ω ⊂ R 3 is a bounded polyhedron. By considering a semi-discretization in time and then ignoring the nonlinear terms, we obtain the following equations: where ∇ · f = 0, and the parameters α, β, γ > 0. We consider homogeneous boundary conditions, The above choice of boundary conditions arise naturally in the variational formulation given below. On the other hand, in the numerical simulations of the problem with pseudo-spectral method, one often uses periodic boundary conditions, e.g., [1,7].
It is worth pointing out that the parameter α is usually much smaller than either β or γ. This fact imposes some difficulties in designing robust numerical methods, as have been studied in the context of biharmonic problems, e.g., [20,30].
The above fourth-order curl equations also arise from an interior transmission problem in the study of inverse scattering problems for inhomogeneous medium, e.g., [3].
In order to provide an appropriate framework for our analysis, we define the following function spaces: V is a Hilbert space with scalar product and norm given by The following lemma gives a sufficient condition for a piecewisely defined function to be an element in V .
Using the following identity: and ∇ · u = 0, the first equation in (2.1) can be rewritten in the following form: Multiplying Equation (2.3) by the test function v and using integration by parts, we obtain the following variational formulation: where the bilinear form a(·, ·) defined on V × V is given by The well-posedness of the above variational problem follows from the Lax-Milgram lemma.
The next lemma indicates that the weak solution satisfies the divergence-free constraint.

Lemma 2.2.
Assume ∇ · f = 0, and let u be the solution of problem (2.4). Then

A Nonconforming Finite Element
In this section, we construct a nonconforming finite element to solve the fourthorder equation. One of the advantages for using a nonconforming element is that the number of degrees of freedom is small compared to that for conforming elements. The following construction is based on Nédélec elements of the first family that consist of incomplete vector polynomials [19]. The advantage of using incomplete vector polynomial space is that it provides the same order of convergence in terms of energy norms as the one given by corresponding complete polynomial space. In the following, we define the degrees of freedom in a special way to ensure that the consistency error estimate holds.
where P 2 is the space of homogeneous multivariate polynomials of degree 2; • Σ K is the set of degrees of freedom, see Figure 1, edge degrees of freedom: where τ is the unit tangential vector along the edge e, face degrees of freedom: where n is the unit normal vector to the face f , In the above finite element triple, the space P K is the same as the second order Nédélec element of the first family for H(curl) problem. The difference is the definition of the second set of degrees of freedom. It is designed specifically to ensure consistency for the fourth-order problems. The total number of the degrees of freedom for this element is 20, which is the same as the dimension of the polynomial space R 2 (K).
It should be pointed out that the scaling factor 1/|f | 2 in the definition of the second set of degrees of freedom is associated with the construction of basis functions to be given later.
The next lemma given in [19] describes a relation between edge integrals and face integrals which will be useful in the error analysis.
Lemma 3.2. If u ∈ R 2 (K) is such that the edge degrees of freedom (3.1) vanish, then Proof. Given u ∈ R 2 (K) satisfies e u · τ q ds = 0, ∀ edge e ⊂ K.
where u T is the tangential part of u, and ∇ f × and ∇ f × are surface vector curl and surface scalar curl, respectively. Let q be a constant. Notice that we conclude, As a direct consequence of Lemma 3.2, if both the edge degrees of freedom (3.1) and face degrees of freedom (3.2) vanish, then The polynomial space R 2 (K) has the following property [8].
We recall that the finite element (K, P K , Σ K ) is said to be unisolvent if a function in P K can be uniquely determined by specifying values for degrees of freedom in Σ K . Proof. It is sufficient to prove that, given u ∈ R 2 (K), Obviously, ∇(∇ × u) is a constant vector. Then using (3.3) and integration by parts, we obtain This implies that Using again (3.3), we have By Lemma 3.3, we have u = ∇p, with p ∈ P 2 (K).
This implies ∂p/∂τ = 0 on each edge e. Hence, p is constant and u = 0.
In the following, we construct the basis functions. The explicit form of these basis functions not only is useful for implementation, but also instrumental for the interpolation error estimate. 3.1. Basis functions. The main idea of the construction is to consider linear combinations of basis functions of a related Nédélec element. Let K be an arbitrary tetrahedron with four vertices a i , a j , a k and a l , see Figure 2. The corresponding barycentric coordinates are given by λ i , λ j , λ k , and λ l , respectively. On each of the four faces, say face l (with vertices a i , a j , a k ), we choose the following two tangential direction vectors: . The edge degrees of freedom on edge e ij (with vertices a i and a j ) are defined explicitly by: where τ is the unit direction vector of edge e ij , s is an arc length parameter. The face degrees of freedom are defined as: where n l is the unit outward normal vector of the face f l .
We recall that the basis functions of the second order Nédélec element of the first family in barycentric coordinates are (see, e.g., [9], [26], [31]): (1) Two basis functions on each edge e ij : Two basis functions on each face f l : In the following, we list a few useful facts about the geometry of a tetrahedron.
(1) The unit outward normal vector of face f l is given by − ∇λ l ∇λ l .
(2) The two tangential vectors of face f l are given by q ij and q ik .
(3) Let h l be the height of the tetrahedron corresponding to the face f l , then (4) Let |K| be the volume of the tetrahedron K, then Next, we construct basis functions in barycentric coordinates. They provide a set of dual basis functions with respect to the prescribed degrees of freedom.
Step 1. Construct eight basis functions {φ lij } corresponding to the face degrees of freedom such that We use the basis functions of the second order Nédélec element as building blocks as they automatically satisfy the first condition (3.4). Using the facts listed above, we find that the basis functions corresponding to the facial degrees of freedom on face f l are given by the following: Similarly, M lik (φ lij ) = 0, and M f l (φ lij ) = 0 where l = l.
Step 2. Construct twelve basis functions {ψ (t) ij |1 ≤ i < j ≤ 4, t = 1, 2} corresponding to the edge degrees of freedom such that Here, we use the edge basis functions of the second order Nédélec element as building blocks since they satisfy condition (3.6). Since ∇ × L ji = 0, L ji automatically satisfy condition (3.7).
For functions L ij , we need to subtract from them a linear combination of face basis functions so that (3.6) and (3.7) hold. This can be done because by construction, our face basis functions have no edge moments. This strategy for constructing basis functions can be found in [9,26].
Finally, we can write the basis functions of the new element as the following: (1) Two basis functions on each face l (1 ≤ l ≤ 4): where L lij = λ i (λ j ∇λ k − λ k ∇λ j ).

Convergence analysis. Let
be a triangulation of the domain Ω. On this triangulation we introduce the finite element space V h and define the discrete norm · h by Consider the following discrete bilinear form: It is straightforward to verify that the bilinear form a h satisfies The nonconforming finite element discretization of problem (3.10) is: The convergence of the above finite element approximation can be analyzed through the following second Strang lemma [4].
where the first term on the right-hand side is called the interpolation error and the second term is called the consistency error.
In order to estimate the consistency error we first define an average operator P f on a face f by Since for any v h ∈ V h , the quantity f ∇ × v h dA is continuous, we know that P f is well-defined for ∇ × v h . The following two lemmas are standard results.
Next, we estimate the interpolation error and consistency error separately.
3.2.1. Interpolation error estimate. Let K and K f be the two tetrahedra sharing a common face f , r K be the local interpolation operator for the second order Nédélec element of the first family, namely, given u ∈ V , define r K u such that e r K u · τ ds = e u · τ ds, ∀ edge e ⊂ K, Lemma 3.8. Given u ∈ V , let u I be defined as above, then Proof. Let r h u be the global interpolation operator defined by r h u| K = r K u. By triangle inequality, By the interpolation error estimate of Nédélec element, we have Notice that on each tetrahedron K, where {φ mnp } are basis functions on face f , and {M mnp (·)} are degrees of freedom on face f . Using Lemma 3.7 and q np L 2 (f ) = O(h 2 ), we get Notice ∇λ i = O(h −1 ), and φ mnp 2 0,K = O(h 7 ), by Cauchy-Schwarz inequality, we have Hence, By inverse inequality, we have Combining these estimates, we get r h u − u I h h|∇ × u| 2,Ω , and the desired estimate follows.
Remark: We note that the error estimate (3.9) indicates that r h u and u I are super-close. Such type of estimate can not usually be obtained by the standard scaling argument (using Bramble-Hilbert lemma). In our proof, we made use of the detailed information of the basis functions constructed in the previous section.
3.2.2. Consistency error estimate. Given a tetrahedron K, in addition to the local interpolation operator r K , we introduce another local interpolation operatorr K corresponding to the first order Nédélec element of the second family, namely,r K u ∈ (P 1 (K)) 3 ⊂ R 2 (K), and e ((r K u) · τ ) q ds = e (u · τ ) q ds, ∀ q ∈ P 1 (e), ∀ edge e ⊂ K.
Consider two tetrahedra K and K f that share a common face f .
where n is the unit outward normal vector of ∂K.
Consider the decomposition (see [8]): where div w K = 0, w K · n| ∂K = 0, and p K ∈ P 2 (K). The following Lemma 3.9 can be found in [11]: As a consequence of Lemma 3.9, we have the following estimate.
Proof. Using the interpolation operators defined above, we obtaiñ Hence,r Now, we can show the following lemma, which is critical for the consistency error estimate.
Lemma 3.11. For ϕ ∈ H(curl; Ω), Proof. By the interpolation error estimates of the Nédélec elements Lemma 3.10, and Equation (3.10), we have Next, we show the consistency error estimate for the nonconforming finite element approximation defined above. Proof. By applying integration by parts, we get
Finally, we have the following convergence result. Proof. Using the second Strang lemma, and previous lemmas, the desired inequality follows.