Structure-preserving Finite Element Methods for Stationary MHD Models
SSTRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHDMODELS
KAIBO HU AND JINCHAO XUA
BSTRACT . In this paper, we develop a class of mixed finite element scheme for stationary magnetohydro-dynamics (MHD) models, using magnetic field B and current density j as the discretization variables. Weshow that the Gauss’s law for the magnetic field, namely ∇· B = 0 , and the energy law for the entire systemare exactly preserved in the finite element schemes. Based on some new basic estimates for H h (div) , weshow that the new finite element scheme is well-posed. Furthermore, we show the existence of solutionsto the nonlinear problems and the convergence of Picard iterations and finite element methods under someconditions.
1. I
NTRODUCTION
In this paper, we develop structure-preserving finite element discretization for the following stationaryincompressible magnetohydrodynamics (MHD) system:( . a) ( u · ∇ ) u − R − e ∆ u − S j × B + ∇ p = f , ( . b) j − R − m ∇ × B = , ( . c) ∇ × E = , ( . d) ∇ · B = 0 , ( . e) ∇ · u = 0 , where the Ohm’s law holds:( . ) j = E + u × B . Here u is the velocity of conducting fluids, p is the pressure, B is the magnetic field, E is the electricfield and j is the volume current density. Dimensionless parameters R e , R m and S are the Reynoldsnumber of fluids, magnetic field and the coupling number respectively.In the study of magnetohydrodynamics (MHD) system, it is well-known that the Gauss’s law for themagnetic field, namely ∇ · B = 0 , is an important condition in numerical computation of MHD system[5, 9]. Nonzero divergence of B will introduce a parallel force, which breaks the energy law. In ourprevious work Hu, Ma and Xu [14], we proposed a class of structure-preserving and energy-stable finiteelement discretizations that exactly preserve the magnetic Gauss’s law on the discrete level for the timedependent MHD systems. The goal of this paper is to extend such discretizations to stationary cases.Such a discretization is however not straightforward as the time-dependent and the stationary systemshave different structures. In the time-dependent problem, the Faraday’s law reads: ∂ B ∂t + ∇ × E = . Mathematics Subject Classification.
Primary 65N30, 65N12.
Key words and phrases.
Divergence-free, Stationary, MHD equations, Finite Element.This material is based upon work supported in part by the US Department of Energy Office of Science, Office of AdvancedScientific Computing Research, Applied Mathematics program under Award Number DE-SC-0014400 and by Beijing InternationalCenter for Mathematical Research of Peking University, China. a r X i v : . [ m a t h . NA ] O c t KAIBO HU AND JINCHAO XU
In [14], we chose to keep the electric field E and use the H (curl) -conforming finite element space for E and H (div) -conforming finite element space for B to discretize the above Faraday’s law as follows: B n − B n − ∆ t + ∇ × E n = . This implies that ∇ · B n = 0 holds for all n ≥ as long as it holds for n = 0 .In the stationary case, the Faraday’s law reads: ∇ × E = . In this case, we can not directly apply the technique used in [14] for the evolutionary case to preserve theGauss’s law ∇ · B = 0 exactly on the discrete level. Instead we treat the Gauss’s law as an independentequation in the whole MHD system and we then introduce a Lagrange multiplier to appropriately enforcethis law on both the continuous and the discrete level.The idea of the use of Lagrange multiplier itself is not new (see Sch¨otzau [18] and the referencetherein) and the novelty of our approach here lies in how this technique is used in combination with thetechniques developed in [14]. In Sch¨otzau [18], a magnetic multiplier r ∈ H (Ω) / R is used to imposethe Gauss’s law in the following way: (cid:90) Ω B · ∇ s = 0 , ∀ s ∈ H (Ω) / R which does not guarantee that the Gauss’s law holds strongly (namely ∇ · B h = 0 point-wise in thedomain) in the corresponding discrete case. The main difference in our approach is that the Gauss’s lawwill indeed be preserved on the discrete level strongly by using appropriate finite element discretizationof B so that B h is H (div) -conforming. On the other hand, the charge conservation ∇· j = 0 is preservedin a weak sense. The finite element de Rham sequence as studied in [1, 13, 4] plays an important role inthe construction and analysis in our paper.MHD equations admit many different variational formulations which lead to different mathematicalproperties and numerical efficiency on the discrete level. In most existing literature, variables E and j are eliminated to reduce the size of the corresponding discretized problems. In [14], we demonstratedthat it is advantageous to keep E and use it as an independent (or intermediate) discretization variable inappropriate finite element space. Indeed, this approach may lead to larger discretized systems, but thesesystems have better mathematical structures and may be solved, as illustrated in [15], more efficientlythan the corresponding smaller systems derived from traditional schemes by eliminating both E and j .In this paper, we continue and extend this study for the stationary problem. Instead of retaining E explicitly as a variable, we choose B and j as electromagnetic variables motivated by the energy law.For simplicity of exposition, we use the following homogeneous Dirichlet boundary conditions u = , B · n = 0 , j × n = . According to the Ohm’s law that j = E + u × B , the above boundary conditions are obviously equivalentto u = , B · n = 0 , E × n = . TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 3
The extension to non-homogeneous boundary conditions is straightforward and standard and the rel-evant details will not be given in this paper.The rest of the paper is organized as follows. In §
2, we present the notation and basic finite elementspaces used in the discussion. § H h (div0) functions, including regu-larity and the discrete Poincar´e’s inequality. In §
4, a new formulation based on B and j is studied. Weprove the well-posedness based on an equivalent reduced system. In §
5, we prove the analysis of theproposed algorithms based on the key technical results established in §
3. This includes the convergenceof Picard iterations and the finite element discretizations. Concluding remarks are given in § OTATION AND BASIC FINITE ELEMENT SPACES
In this section, we introduce some basic Sobolev spaces and their corresponding finite element dis-cretizations that will be used in the rest of the paper.We assume that Ω is a bounded Lipschitz polyhedron. For the ease of exposition, we further assumethat Ω is contractable, i.e. there is no nontrivial harmonic form. For general domains (non-simply-connected domain, non-connected boundary), we can solve the problem in the orthogonal complementof (discrete) harmonic forms, as in Arnold, Falk and Winther [1] for the Hodge Laplacian. Thereforesuch an assumption on the domain is to make the presentation more clear, and the methodology is alsovalid for general topology.Using the standard notation for inner product and norm of the L space ( u, v ) := (cid:90) Ω u · v d x, (cid:107) u (cid:107) := (cid:18)(cid:90) Ω | u | d x (cid:19) / , we define the following H ( D, Ω) space with a given linear operator D : H ( D, Ω) := { v ∈ L (Ω) , Dv ∈ L (Ω) } , and H ( D, Ω) := { v ∈ H ( D, Ω) , t D v = 0 on ∂ Ω } , where t D is the trace operator: t D v := v, D = grad ,v × n, D = curl ,v · n, D = div . Here H (grad , Ω) is a scalar function space, while H (curl , Ω) and H (div , Ω) are for vector valuedfunctions. We often use the following notation: L (Ω) := (cid:26) v ∈ L (Ω) : (cid:90) Ω v = 0 (cid:27) . When D = grad , we often use the notation: H (Ω) := H (grad , Ω) , H (Ω) := H (grad , Ω) . For clarity, the corresponding norms in H ( D, Ω) are denoted by (cid:107) u (cid:107) = (cid:107) u (cid:107) + (cid:107)∇ u (cid:107) , (cid:107) F (cid:107) := (cid:107) F (cid:107) + (cid:107)∇ × F (cid:107) , (cid:107) C (cid:107) := (cid:107) C (cid:107) + (cid:107)∇ · C (cid:107) . KAIBO HU AND JINCHAO XU
We will also use the space L p with norm (cid:107)·(cid:107) ,p given by (cid:107) v (cid:107) p ,p = (cid:82) Ω | v | p . For a general Banach space Y with a norm (cid:107)·(cid:107) Y , the dual space Y ∗ is equipped with the dual norm defined as (cid:107) h (cid:107) Y ∗ := sup (cid:54) = y ∈ Y (cid:104) h , y (cid:105)(cid:107) y (cid:107) Y . For the special case that Y = H (Ω) , Y ∗ = H − (Ω) and the corresponding norm is denoted by (cid:107)·(cid:107) − ,which is defined as (cid:107) f (cid:107) − := sup (cid:54) = v ∈ H (Ω) (cid:104) f , v (cid:105)(cid:107)∇ v (cid:107) . We will use C to denote the constant in the following inequality, which is a consequence of Sobolevimbedding theorem and Poincar´e’s inequality: (cid:107) u (cid:107) , ≤ C (cid:107)∇ u (cid:107) , ∀ u ∈ H (Ω) . ( . )Since the fluid convection frequently appears in the following discussions, we introduce the trilinearform L ( w ; u , v ) := 12 [(( w · ∇ ) u , v ) − (( w · ∇ ) v , u )] . When w is a known function, L ( w ; u , v ) is a bilinear form of u and v . This will occur in the Picarditeration, where w is the velocity of the last iteration step.Let T h be a triangulation of Ω , and we assume that the mesh is regular and quasi-uniform, so that theinverse estimates hold [6]. The finite element de Rham sequence is an abstract framework to unify theabove spaces and their discretizations, see e.g. Arnold, Falk, Winther [1, 2], Hiptmair [13], Bossavit [4]for more detailed discussions. Figure 1 shows the commuting diagrams we will use. Current density j ,magnetic field B and the multiplier r will be discretized in the last three spaces respectively. Figure 2shows the finite elements of the lowest order. H (grad) grad −−−−→ H (curl) curl −−−−→ H (div) div −−−−→ L (cid:121) Π grad h (cid:121) Π curl h (cid:121) Π div h (cid:121) Π h H h (grad) grad −−−−→ H h (curl) curl −−−−→ H h (div) div −−−−→ L ,h F IGURE
1. Continuous and discrete de Rham sequenceF
IGURE
2. DOF of finite element de Rham sequence of lowest orderAs we shall see, H (div) functions with vanishing divergence will play an important role in the study.So we define on the continuous level H (div0 , Ω) := { C ∈ H (div , Ω) : ∇ · C = 0 } , TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 5 and the finite element subspace H h (div0 , Ω) := { C h ∈ H h (div , Ω) : ∇ · C h = 0 } . We use V h to denote the finite element subspace of velocity u h , and Q h for pressure p h . Thereare many existing stable pairs for V h and Q h , for example, Taylor-Hood elements [11, 3]. Spaces H h (div , Ω) and L ,h (Ω) are finite element spaces from the discrete de Rham sequence. For these spaceswe use their explicit names for clarity, and use the notation V h and Q h for the fluid part to indicate thatthey are usually different from H ,h (Ω) and L ,h (Ω) in the de Rham sequence.There is a unified theory for the discrete de Rham sequence of arbitrary order [3, 1, 2]. In the case n = 3 , the lowest order elements can be represented as: R → P Λ d −→ P Λ d −→ P Λ d −→ P Λ → R → P Λ d −→ P Λ d −→ P − Λ d −→ P Λ → R → P Λ d −→ P − Λ d −→ P Λ d −→ P Λ → R → P Λ d −→ P − Λ d −→ P − Λ d −→ P Λ → The correspondence between the language of differential forms and classical finite element methods issummarized in Table 1.To link the finite element spaces, below we will require H h (curl , Ω) , H h (div , Ω) and L ,h (Ω) to bein the same sequence. k Λ kh (Ω) Classical finite element space0 P r Λ ( T ) Lagrange elements of degree ≤ r P r Λ ( T ) Nedelec 2nd-kind H (curl) elements of degree ≤ r P r Λ ( T ) Nedelec 2nd-kind H (div) elements of degree ≤ r P r Λ ( T ) discontinuous elements of degree ≤ r P − r Λ ( T ) Lagrange elements of degree ≤ r P − r Λ ( T ) Nedelec 1st-kind H (curl) elements of order r − P − r Λ ( T ) Nedelec 1st-kind H (div) elements of order r − P − r Λ ( T ) discontinuous elements of degree ≤ r − T ABLE
1. Correspondences between finite element differential forms and the classicalfinite element spaces for n = 3 (from [1])As we shall see, it is useful to group the spaces to define X h := V h × H h (curl , Ω) × H h (curl , Ω) × H h (div , Ω) . and group Q h × L ,h (Ω) to define Y h := Q h × L ,h (Ω) . For the analysis, we also need to define a reduced space, where j h and σ h (introduced below) areeliminated: ˜ X h := V h × H h (div , Ω) . KAIBO HU AND JINCHAO XU
In order to define appropriate norms, we introduce the discrete curl operator on the discrete level. Forany C h ∈ H h (div , Ω) , define ∇ h × C h ∈ H h (curl , Ω) : ( ∇ h × C h , F h ) = ( C h , ∇ × F h ) , ∀ F h ∈ H h (curl , Ω) . For any w h ∈ H h (curl , Ω) , we define ∇ h · w h ∈ H h (grad , Ω) by ( ∇ h · w h , v h ) = − ( w h , ∇ v h ) , ∀ v h ∈ H h (grad , Ω) . We define P : L (Ω) → H h (curl , Ω) to be the L projection ( P φ, F h ) = ( φ, F h ) , ∀ F h ∈ H h (curl , Ω) , φ ∈ L (Ω) . We further define (cid:107)·(cid:107) d to be a modified norm of H h (div , Ω) by (cid:107) C h (cid:107) d := (cid:107) C h (cid:107) + (cid:107)∇ · C h (cid:107) + (cid:107)∇ h × C h (cid:107) . Moreover, (cid:107)·(cid:107) c for H h (curl , Ω) is simply the L norm: (cid:107) F h (cid:107) c := (cid:107) F h (cid:107) . There are some motivations to define such a stronger norm for H h (div , Ω) and weaker norm for H h (curl , Ω) space. One technical reason is that we want the nonlinear term ∇ × ( u h × B h ) to be bounded in someproper discretization. But generally u h × B h may not belong to H h (curl) for u h ∈ H (Ω) and B h ∈ H (div , Ω) . So we choose to move the curl operator to the H h (div) test function in the varia-tional formulation to get ( u h × B h , ∇ h × C h ) . Therefore we add the weak curl norm to H h (div , Ω) space. Another motivation can be seen in the energy estimate: on the continuous level, the energy esti-mate contains j = R − m ∇ × B , but not ∇ × j . So it is natural to use L norm for the discrete variable j h .Now we define the norms for various product spaces. For Y h space, we define (cid:107) ( q, r ) (cid:107) Y := (cid:107) q (cid:107) + (cid:107) r (cid:107) . For the other product spaces, we define (cid:107) ( u h , j h , σ h , B h ) (cid:107) X := (cid:107) u h (cid:107) + (cid:107) j h (cid:107) c + (cid:107) σ h (cid:107) c + (cid:107) B h (cid:107) d , ( u h , j h , σ h , B h ) ∈ X h , and (cid:107) ( u h , B h ) (cid:107) X := (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d , ( u h , B h ) ∈ ˜ X h .
3. E
STIMATES FOR DIVERGENCE - FREE VECTOR FIELDS
In this section, we will establish some new regularity results for the strong divergence-free space H h (div0 , Ω) which will be used for our forthcoming analysis. The main ingredients used in our analysisinclude some regularity results for the space Z := H (curl , Ω) ∩ H (div0 , Ω) (c.f. [13, 18]), and for thespace X ch := { w ∈ H h (curl , Ω) : ∇ h · w h = 0 } (c.f. [13, 18]), together with some appropriately defined “Hodge mapping” ( H d below) that connects H h (div0 , Ω) with Z .We first give a preliminary result based on Hodge decomposition: Lemma 1. ∇ × Z = H (div0 , Ω) = ∇ × H (curl , Ω) . TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 7
Proof.
From the Hodge decomposition for L (Ω) : L (Ω) = ∇ H (Ω) + ∇ × H (curl , Ω) = H (curl0 , Ω) + H (div0 , Ω) . Here H (curl0 , Ω) := { F ∈ H (curl , Ω) : ∇ × F = } . Therefore H (curl , Ω) = L (Ω) ∩ H (curl , Ω)= H (curl0 , Ω) + H (div0 , Ω) ∩ H (curl , Ω)= H (curl0 , Ω) + Z . ( . )This implies H (div0 , Ω) = ∇ × H (curl , Ω) = ∇ × Z . (cid:3) We now define the “Hodge mapping” for H h (div0) functions. Let H d : H h (div0) → Z be definedby ( ∇ × ( H d B h ) , ∇ × v ) = ( ∇ h × B h , ∇ × v ) , ∀ v ∈ Z , ∀ B h ∈ H h (div0 , Ω) . ( . )Due to the Poincar´e’s inequality of Z , (cid:107) z (cid:107) (cid:46) (cid:107)∇ × z (cid:107) holds for any z ∈ Z . Therefore ( . ) uniquelydefines H d B h .From Lemma 1, we have ∇ × Z = H (div0) . Therefore ( ∇ × ( H d B h ) , w ) = ( ∇ h × B h , w ) , ∀ w ∈ H (div0) . ( . )In particular, choosing w = ∇ × ( H d B h ) , we see (cid:107)∇ × ( H d B h ) (cid:107)≤ (cid:107)∇ h × B h (cid:107) . In the following, we will use ˜ B to denote the continuous lifting of B h : ˜ B := H d B h . Moreover, H c : X ch → H (curl , Ω) ∩ H (div , Ω) is the Hodge mapping for H h (curl , Ω) [13, 18], definedby ∇ × ( H c F h ) = ∇ × F h , ∀ F h ∈ X ch . We also use the notation ˜ F to denote H c F h when F h ∈ X ch . Lemma 2 (Approximation of H d ) . If Ω is a bounded polyhedral domain in R , there exists < δ (Ω) ≤ such that (cid:107) B h − H d B h (cid:107) (cid:46) h + δ (cid:107)∇ h × B h (cid:107) , for all B h ∈ H h (div0 , Ω) .Proof. We define Π h div to be the bounded cochain projection to H h (div , Ω) [10]. Note that ∇· (cid:16) B h − Π h div ˜ B (cid:17) =0 due to the commuting diagram. Therefore there exists φ h ∈ X ch and the corresponding lifting ˜ φ := H c φ h ∈ H (curl , Ω) ∩ H (div , Ω) such that B h − Π h div ˜ B = ∇ × φ h = ∇ × ˜ φ and thereexists a positive constant < δ (Ω) ≤ such that( . ) (cid:107) φ h − ˜ φ (cid:107) (cid:46) h + δ (cid:107)∇ × φ h (cid:107) = h + δ (cid:107) B h − Π h div ˜ B (cid:107) , where the first inequality is from the approximation property of H c .From ( . ), we have ( ∇ h × B h , ˜ φ ) = ( ∇ × ˜ B , ˜ φ ) = ( ˜ B , ∇ × ˜ φ ) , KAIBO HU AND JINCHAO XU and ( B h , ∇ × φ h ) = ( ∇ h × B h , φ h ) = ( ∇ h × B h , φ h − ˜ φ ) + ( ˜ B , ∇ × ˜ φ ) . Namely, ( B h − ˜ B , B h − Π h div ˜ B ) = ( ∇ h × B h , φ h − ˜ φ ) . Thus (cid:107) B h − ˜ B (cid:107) = ( B h − ˜ B , B h − Π h div ˜ B ) + ( B h − ˜ B , Π h div ˜ B − ˜ B )= ( ∇ h × B h , φ h − ˜ φ ) + ( B h − ˜ B , Π h div ˜ B − ˜ B ) . By ( . ) and the interpolation error estimates (cid:107) ˜ B − Π h div ˜ B (cid:107) (cid:46) h + δ (cid:107) ˜ B (cid:107) + δ (cid:46) h + δ (cid:107)∇ × ˜ B (cid:107) (cid:46) h + δ (cid:107)∇ h × B h (cid:107) , we obtain (cid:12)(cid:12)(cid:12) ( ∇ h × B h , φ h − ˜ φ ) (cid:12)(cid:12)(cid:12) (cid:46) h + δ (cid:107) B h − Π h div ˜ B (cid:107)(cid:107)∇ h × B h (cid:107)≤ h + δ (cid:16) (cid:107) B h − ˜ B (cid:107) + (cid:107) ˜ B − Π h div ˜ B (cid:107) (cid:17) (cid:107)∇ h × B h (cid:107)≤ h + δ (cid:107) B h − ˜ B (cid:107)(cid:107)∇ h × B h (cid:107) + h δ (cid:107)∇ h × B h (cid:107) ≤ (cid:107) B h − ˜ B (cid:107) + 12 h δ (cid:107)∇ h × B h (cid:107) + h δ (cid:107)∇ h × B h (cid:107) , and hence (cid:107) B h − ˜ B (cid:107) (cid:46) (cid:107) ˜ B − Π h div ˜ B (cid:107) + h δ (cid:107)∇ h × B h (cid:107) . This completes the proof. (cid:3)
For nonlinear problems and their linearizations, it is technical to prove the boundedness of variationalforms, and this often requires careful estimates of regularity. The nonlinear terms in the variational formsproposed in this paper will have the form ( u h × B h , j h ) , where u h ∈ V h ⊂ H (Ω) , B h ∈ H h (div0 , Ω) and j h ∈ H h (curl , Ω) . Lemma 3.
For u h ∈ V h and B h ∈ H h (div0 , Ω) , we have the following bound: (cid:107) u h × B h (cid:107) (cid:46) (cid:107) u h (cid:107) (cid:107)∇ h × B h (cid:107) . Proof.
From Lemma 2, we have (cid:107) B h − ˜ B (cid:107) (cid:46) h + δ (cid:107)∇ h × B h (cid:107) , where < δ ≤ is a positive constant depending on the domain.Then (cid:107) u h × B h (cid:107)≤ (cid:107) u h × ( B h − ˜ B ) (cid:107) + (cid:107) u h × ˜ B (cid:107) . For the first term, (cid:107) u h × ( B h − ˜ B ) (cid:107) ≤ (cid:107) u h (cid:107) , ∞ (cid:107) B h − ˜ B (cid:107) (cid:46) h − (cid:107) u h (cid:107) , · h (cid:107)∇ h × B h (cid:107) (cid:46) (cid:107) u h (cid:107) (cid:107)∇ h × B h (cid:107) , where the second inequality comes from the inverse estimates and the approximation results. TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 9
Due to the regularity of Z [13], we have (cid:107) u h × ˜ B (cid:107) ≤ (cid:107) u h (cid:107) , (cid:107) ˜ B (cid:107) , (cid:46) (cid:107) u h (cid:107) (cid:107)∇ × ˜ B (cid:107)≤ (cid:107) u h (cid:107) (cid:107)∇ h × B h (cid:107) . This implies (cid:107) u h × B h (cid:107) (cid:46) (cid:107) u h (cid:107) (cid:107)∇ h × B h (cid:107) . (cid:3) Below we will use a positive constant C to denote the bound: (cid:107) u h × B h (cid:107)≤ C (cid:107)∇ u h (cid:107)(cid:107)∇ h × B h (cid:107) , ( . )and therefore ( u h × B h , j h ) ≤ C (cid:107)∇ u h (cid:107)(cid:107)∇ h × B h (cid:107)(cid:107) j h (cid:107) . In the discussions below, we will need discrete Poincar´e’s inequality for H h (div0 , Ω) functions. Wenote that the two dimensional case is given in [7], and the proof can be modified to adapt to the threedimensional case. We include a different proof here. Lemma 4.
For B h ∈ H h (div0 , Ω) , we have the following discrete Poincar´e’s inequality: (cid:107) B h (cid:107) (cid:46) (cid:107)∇ h × B h (cid:107) . Proof.
Because ∇ · B h = 0 , we can choose E h ∈ H h (curl) such that ∇ × E h = B h , and ∇ h · E h = 0 . From the discrete Poincar´e inequality for X ch in [2], (cid:107) E h (cid:107) curl (cid:46) (cid:107)∇ × E h (cid:107) = (cid:107) B h (cid:107) . ( . )We have (cid:107)∇ h × B h (cid:107) = sup F h ∈ H h (curl) ( ∇ h × B h , F h ) (cid:107) F h (cid:107) = sup F h ∈ H h (curl) ( B h , ∇ × F h ) (cid:107) F h (cid:107) . ( . )Therefore combining ( . ) and ( . ), we get (cid:107)∇ h × B h (cid:107)≥ ( B h , ∇ × E h ) (cid:107) E h (cid:107) , and (cid:107)∇ h × B h (cid:107) (cid:38) (cid:107) B h (cid:107) . (cid:3) Combined with L p - L p bounded interpolations (c.f. [8]), we can further establish L p estimates of H (div0) finite element functions. Theorem 1.
For bounded Lipschitz polyhedral domain Ω , we have (cid:107) B h (cid:107) , (cid:46) (cid:107)∇ h × B h (cid:107) , B h ∈ H h (div0 , Ω) . Proof.
From triangular inequality, we have (cid:107) B h (cid:107) , ≤ (cid:107) B h − Π h div H d B (cid:107) , + (cid:107) Π h div H d B (cid:107) , . From inverse estimates, interpolation error estimates and the approximation of Hodge mapping (Lemma2), (cid:107) B h − Π h div H d B h (cid:107) , (cid:46) h − / (cid:13)(cid:13) B h − Π h div H d B h (cid:13)(cid:13) (cid:46) h − / ( (cid:107) B h − H d B h (cid:107) + (cid:13)(cid:13) H d B h − Π h div H d B h (cid:13)(cid:13) ) (cid:46) (cid:107)∇ h × B h (cid:107) . Using the L stability of the interpolation operator and regularity results of Z , we have (cid:13)(cid:13) Π h div H d B h (cid:13)(cid:13) , (cid:46) (cid:107) H d B h (cid:107) , (cid:46) (cid:107)∇ × H d B h (cid:107) ≤ (cid:107)∇ h × B h (cid:107) . Then the triangular inequality implies (cid:107) B h (cid:107) , ≤ (cid:107) B h − Π h div H d B h (cid:107) , + (cid:13)(cid:13) Π h div H d B h (cid:13)(cid:13) , (cid:46) (cid:107)∇ h × B h (cid:107) . (cid:3) In following discussions, we still use a generic constant C to denote the bound (cid:107) B h (cid:107)≤ C (cid:107)∇ h × B h (cid:107) , ∀ B h ∈ H h (div , Ω) .
4. A
NEW FINITE ELEMENT FORMULATION
In Hu, Ma and Xu [14], the authors studied a numerical scheme using B and E as variables. Astraightforward analysis by Brezzi theory leads to a stringent condition on the time step size. In thissection, we propose a new finite element scheme whose well-posedness will not depend on such as-sumptions.We note that it is the variable j that appears in the energy estimate. Therefore it seems natural to use B and j as mixed variables of the electromagnetic part of the MHD system. Discretization methods basedon B and j actually have already existed in the literature. For example, some finite volume methodsusing B and j have been developed in [17, 16] where the conservation of ∇ · j = 0 was considered (butno discussion on the condition ∇ · B = 0 ), and in [19], B and j were used as variables in the simulationof liquid metal breeder blankets.We eliminate E by Ohm’s law and consider the following model:( . a) − R − e ∆ u + ∇ p + ( u · ∇ ) u + S B × j = f , ( . b) ∇ × j − ∇ × ( u × B ) = , ( . c) j − R − m ∇ × B = , ( . d) ∇ · u = 0 , ( . e) ∇ · B = 0 . The well-posedness of the continuous formulation has been shown in [18]. The author proved thatthere exists at least one solution u ∈ H (Ω) , B ∈ H (curl , Ω) ∩ H (div0 , Ω) for the nonlinear sys-tem where j is eliminated. The variational form reads: find ( u , B , p, φ ) ∈ H (Ω) × H (curl , Ω) ∩ TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 11 H (div , Ω) × L (Ω) × H (Ω) such that for any ( v , C , q, ψ ) ∈ H (Ω) × H (curl , Ω) ∩ H (div , Ω) × L (Ω) × H (Ω) , L ( u ; u , v ) + R − e ( ∇ u , ∇ v ) − s (( ∇ × B ) × B , v ) − ( p, ∇ · v ) = (cid:104) f , v (cid:105)− ( u × B , ∇ × C ) + R − m ( ∇ × B , ∇ × C ) + ( ∇ φ, C ) = 0 , ( ∇ · u , q ) = 0 , ( B , ∇ ψ ) = 0 . ( . )Considering j = R − m ∇ × B as an intermediate variable, we conclude with the existence of solutionsto ( . ): for any f ∈ (cid:0) H (Ω) (cid:1) ∗ , there exists at least one solution u ∈ H (Ω) , B ∈ H (curl , Ω) ∩ H (div0 , Ω) and j ∈ L (Ω) .4.1. Mixed finite element discretizations.
We now present our new finite element discretization of theabove system ( . ). Problem 1.
Given f ∈ V ∗ h . Find ( u h , j h , σ h , B h , p h , r h ) ∈ X h × Y h , such that for any ( v h , k h , τ h , C h , q h , s h ) ∈ X h × Y h ,( . a) R − e ( ∇ u h , ∇ v h ) + L ( u h ; u h , v h ) + S ( j h , v h × B h ) − ( p h , ∇ · v h ) = (cid:104) f , v h (cid:105) , ( . b) SR − m ( ∇ × j h , C h ) − SR − m ( ∇ × σ h , C h ) + ( r h , ∇ · C h ) = 0 , ( . c) SR − m ( σ h , τ h ) − SR − m ( u h × B h , τ h ) = 0 , ( . d) S ( j h , k h ) − SR − m ( B h , ∇ × k h ) = 0 , ( . e) − ( ∇ · u h , q h ) = 0 , ( . f ) ( ∇ · B h , s h ) = 0 . In the above scheme, an additional variable σ h is introduced to accommodate for the evaluation ofthe discrete curl operator ∇ h × which is nonlocal. This extra work comes from the nonlinear couplingterm ( ∇ × ( u × B ) , C ) , because curl operator cannot act on u × B directly.Before further discussions, we verify basic properties of the discretization and the energy estimates,which are basic and important tools in the design and analysis of numerical methods, especially fornonlinear problems. Theorem 2.
Any solution of Problem 1 satisfies (1)
Gauss’s law of magnetic field in the strong sense: ∇ · B h = 0 . (2) the Lagrange multiplier r h = 0 , hence ( . b) reduces to ∇ × ( j h − σ h ) = . (3) the energy estimates: R − e (cid:107)∇ u h (cid:107) + S (cid:107) j h (cid:107) = (cid:104) f , u (cid:105) , and R e (cid:107)∇ u h (cid:107) + S (cid:107) j h (cid:107) ≤ R e (cid:107) f (cid:107) − . Proof. (1) It is a direct consequence of ( . f), since ∇ · H h (div; Ω) = L ,h (Ω) . (2) Take C h = ∇ × j h − ∇ × σ h . From ( . b) we see ∇ × j h − ∇ × σ h = . Hence ( r h , ∇ · C h ) = 0 , ∀ C h ∈ H h (div , Ω) . This implies r h = 0 . (3) Take v h = u h , C h = B h , τ h = ∇ h × B h , k h = j h in ( . a)-( . d). Add them together, wehave R − e (cid:107)∇ u h (cid:107) + S (cid:107) j h (cid:107) + S ( j h , u h × B h ) − SR − m ( u h × B h , ∇ h × B h ) = (cid:104) f , u h (cid:105) . Again from ( . d), the last two terms on the left hand side vanish by taking k h = P ( u h × B h ) .This implies the desired result. (cid:3) From ( . c), we see σ h = P ( u h × B h ) , and from ( . d), we get j h = R − m ∇ h × B h . To provethe existence of solution of the nonlinear scheme, we formally eliminate σ h and j h using the aboveidentities, to get a system with u h , B h , p h and s h .For this purpose, we define ˜ a ( ˜ ψ h ; ˜ ξ h , ˜ η h ) := R − e ( ∇ u h , ∇ v h ) + L ( w h ; u h , v h )+ SR − m ( ∇ h × B h , v h × G h ) − SR − m ( u h × G h , ∇ h × C h )+ SR − m ( ∇ h × B h , ∇ h × C h ) , and b ( ˜ ξ h , y h ) := − ( ∇ · u h , q h ) + ( ∇ · B h , s h ) . Hereafter, ˜ ψ h , ˜ ξ h , ˜ η h , y h are short for ( w h , G h ) , ( u h , B h ) , ( v h , C h ) ∈ ˜ X h and ( q h , s h ) ∈ Y h .Eliminating j h and σ h , Problem 1 is equivalent to the following form. Problem 2.
Given ˜ θ = ( f , ) ∈ ˜ X ∗ h , find ˜ ξ h ∈ ˜ X h , x h ∈ Y h , such that ˜ a ( ˜ ξ h ; ˜ ξ h , ˜ η h ) + b ( ˜ η h , x h ) = (cid:104) ˜ θ , ˜ η h (cid:105) , ∀ ˜ η h ∈ ˜ X h , ( . ) b ( ˜ ξ h , y h ) = 0 , ∀ y h ∈ Y h . ( . ) where (cid:104) ˜ θ , ˜ η h (cid:105) := (cid:104) f , v h (cid:105) . To see the equivalence, we note that if ( u h , j h , σ h , B h , p h , r h ) ∈ X h × Y h solves Problem 1, then ( u h , B h , p h , r h ) ∈ ˜ X h × Y h solves Problem 2 with the same data and (cid:107) ( u h , B h ) (cid:107) ˜ X ≤ (cid:107) ( u h , j h , σ h , B h ) (cid:107) X .Conversely, from a solution ( u h , B h , p h , r h ) of Problem 2, we can reconstruct ( u h , ∇ h × B h , P ( u h × B h ) , B h , p h , r h ) ∈ X h × Y h which solves Problem 1 with the same data, and (cid:107) ( u h , ∇ h × B h , P ( u h × B h ) , B h ) (cid:107) X ≤ (cid:107) ( u h , B h ) (cid:107) ˜ X . Such a variational form is closely related to the “curl-formulation”, for example, in [18]. Here curloperators are replaced by its discrete version “ ∇ h × ”.The existence of solution of the nonlinear discrete scheme ( . ) can be stated as Theorem 3.
There exists at least one solution ( u h , B h , p h , r h ) ∈ ˜ X h × Y h solving Problem 2. Thereforethere exists at least one solution ( u h , j h , σ h , B h , p h , r h ) ∈ X h × Y h solving Problem 1. TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 13
It suffices to prove the existence of solution of Problem 2 under the norm (cid:107) ( u h , B h , p h , r h ) (cid:107) A := (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d + (cid:107) p h (cid:107) + (cid:107) r h (cid:107) . ( . )Define the kernel space ˜ X h by ˜ X h := { ˜ η h ∈ ˜ X h : b ( ˜ η h , y h ) = 0 , ∀ y h ∈ Y h } . Following a general routine of Brezzi theory, we first establish boundedness and inf-sup conditions ofthe variational form.
Lemma 5. (Boundedness)
With the norms given in ( . ) , b ( · , · ) is bounded and ˜ a ( · ; · , · ) is bounded in ˜ X h . From the construction of solutions in Brezzi theory, we note that it is enough to prove the boundednessof ˜ a ( · ; · , · ) in ˜ X h . Proof.
From Cauchy inequality and imbedding theorem, (( u h · ∇ ) u h , v h ) ≤ (cid:107) u h (cid:107) , (cid:107)∇ u h (cid:107)(cid:107) v h (cid:107) , (cid:46) (cid:107) u h (cid:107) (cid:107)(cid:107) u h (cid:107) (cid:107) v h (cid:107) . Similarly, (( u h · ∇ ) v h , u h ) (cid:46) (cid:107) u h (cid:107) (cid:107)(cid:107) u h (cid:107) (cid:107) v h (cid:107) . Furthermore, from Lemma 3, ( j h , v h × B h ) (cid:46) (cid:107)∇ h × B h (cid:107)(cid:107) j h (cid:107) c (cid:107) v h (cid:107) ≤ (cid:107) B h (cid:107) d (cid:107) j h (cid:107) c (cid:107) v h (cid:107) , and ( u h × B h , ∇ h × C h ) (cid:46) (cid:107)∇ h × B h (cid:107)(cid:107) u h (cid:107) (cid:107) C h (cid:107) d ≤ (cid:107) B h (cid:107) d (cid:107) u h (cid:107) (cid:107) C h (cid:107) d . The boundedness of other linear terms are obvious. (cid:3)
Here we note again that the estimate of the boundedness of ( u h × B h , ∇ h × C h ) is a major motivationof introducing the modified (cid:107)·(cid:107) c and (cid:107)·(cid:107) d norms, because u h × B h may not be in H h (curl , Ω) , so curl is actually a discrete operator acting on the H h (div) function B h . Lemma 6. (inf-sup condition of b ( · , · ) ) There exists a positive constant α such that inf y h ∈ Y h sup ˜ η h ∈ ˜ X h b ( ˜ η h , y h ) (cid:107) ˜ η h (cid:107) ˜ X (cid:107) y h (cid:107) Y ≥ α > . Proof.
It suffices to prove the following two inf-sup conditions of the pressure and magnetic multipliers:there exists constant α > such that inf q h ∈ Q h sup v h ∈ V h ( ∇ · v h , q h ) (cid:107) v h (cid:107) (cid:107) q h (cid:107) ≥ α > , inf s h ∈ L ,h (Ω) sup C h ∈ H h (div;Ω) ( ∇ · C h , s h ) (cid:107) C h (cid:107) d (cid:107) s h (cid:107) ≥ α > . The first inequality is standard for existing Stokes pairs. Now we focus on the second. The proof is athree dimensional case of the discussion in Chen et al. [7]. We include the proof here for completeness.The major difficulty is that (cid:107)·(cid:107) d is a stronger norm than (cid:107)·(cid:107) div .It is known that for any s h ∈ L ,h (Ω) , there exists v ∈ H (Ω) , such that ∇ · v = s h , and (cid:107) v (cid:107) (cid:46) (cid:107) s h (cid:107) . Let Π div and Π be the interpolation in H h (div , Ω) and L ,h (we refer to [3] for the definition, and[10] for the local bounded cochain projection, which is bounded in H (curl) and H (div) ). We denote v h = Π div v . Then ∇ · v h = ∇ · Π div v = Π ∇ · v = Π s h = s h . Note that v ∈ H (Ω) , hence Π div is well-defined and bounded. Therefore (cid:107) v h (cid:107) = (cid:107) Π div v (cid:107)≤ (cid:107) Π div (cid:107)(cid:107) v (cid:107) (cid:46) (cid:107) Π div (cid:107)(cid:107) s h (cid:107) . Now it suffices to prove (cid:107)∇ h × v h (cid:107) (cid:46) (cid:107) s h (cid:107) .In fact, using inverse inequality and approximation results (see, for example, [6] and [3]), ( ∇ h × v h , ∇ h × v h ) = ( ∇ h × v h − ∇ × v , ∇ h × v h ) + ( ∇ × v , ∇ h × v h )= ( v h − v , ∇ × ∇ h × v h ) + ( ∇ × v , ∇ h × v h ) (cid:46) h − (cid:107) v h − v (cid:107)(cid:107)∇ h × v h (cid:107) + (cid:107) v (cid:107) (cid:107)∇ h × v h (cid:107) (cid:46) (cid:107) v (cid:107) (cid:107)∇ h × v h (cid:107) (cid:46) (cid:107) s h (cid:107)(cid:107)∇ h × v h (cid:107) . Therefore (cid:107)∇ h × v h (cid:107) (cid:46) (cid:107) s h (cid:107) . This proves the desired result. (cid:3)
Next we consider to solve the subsystem related to ˜ a ( · ; · , · ) . We introduce the existence theorem fornonlinear variational forms, which is given in, for example, [11]. Since we focus on the discrete levelhere, we only given the results for finite dimensional problems. Theorem 4.
Assume that the dimension of V is finite, and there exists a positive constant α such thatbounded trilinear form a ( · ; · , · ) on V satisfies a ( v ; v , v ) ≥ α (cid:107) v (cid:107) , ∀ v ∈ V . Then the problem: given f ∈ V ∗ , find u ∈ V , such that for all v ∈ V , a ( u ; u , v ) = f ( v ) , has at least one solution. It is easy to see that ˜ a ( ˜ ξ h ; ˜ ξ h , ˜ ξ h ) = R − e (cid:107)∇ u (cid:107) + (cid:107)∇ h × B (cid:107) . From the discrete Poincar´e inequality (Lemma 4), we have ˜ a ( ˜ ξ h ; ˜ ξ h , ˜ ξ h ) (cid:38) (cid:107) ˜ ξ h (cid:107) X on ˜ X h . Thereforethe condition in Theorem 4 is satisfied with V = ˜ X h and a ( · ; · , · ) = ˜ a ( · ; · , · ) .Combining Theorem 4 with the boundedness (Lemma 5) and the inf-sup condition of b ( · , · ) (Lemma6), we have proved the existence of solution of nonlinear discrete problem (Theorem 3). TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 15
Picard iterations.
In order to solve nonlinear Problem 1, the following Picard iteration can beused:
Algorithm 1.
For n = 1 , , , . . . , given ( u n − h , B n − h ) ∈ V h × H h (div , Ω) , f ∈ V ∗ h . Find ( u nh , j nh , σ nh , B nh , p nh , r nh ) ∈ X h × Y h , such that for any ( v h , k h , τ h , C h , q h , s h ) ∈ X h × Y h ,( . a) R − e ( ∇ u nh , ∇ v h ) + L ( u n − h ; u nh , v h ) + S ( j nh , v h × B n − h ) − ( p nh , ∇ · v h ) = (cid:104) f , v h (cid:105) , ( . b) SR − m ( ∇ × j nh , C h ) − SR − m ( ∇ × σ nh , C h ) + ( r nh , ∇ · C h ) = 0 , ( . c) SR − m ( σ nh , τ h ) − SR − m ( u nh × B n − h , τ h ) = 0 , ( . d) S ( j nh , k h ) − SR − m ( B nh , ∇ × k h ) = 0 , ( . e) − ( ∇ · u nh , q h ) = 0 , ( . f ) ( ∇ · B nh , s h ) = 0 . The following basic properties of Algorithm 1 can be also established similarly.
Theorem 5.
Any solution of Algorithm 1 satisfies (1)
Gauss’s law of magnetic field in the strong sense: ∇ · B nh = 0 . (2) the Lagrange multiplier r nh = 0 , hence ( . b) reduces to ∇ × ( j nh − σ nh ) = . (3) the energy estimates: R − e (cid:107)∇ u nh (cid:107) + S (cid:107) j nh (cid:107) = (cid:104) f , u nh (cid:105) , and R e (cid:107)∇ u nh (cid:107) + S (cid:107) j nh (cid:107) ≤ R e (cid:107) f (cid:107) − . ( . )We also recast Algorithm 1 into an abstract form of Brezzi theory for the convenience of analysis.We will use ξ h , η h to denote ( u h , j h , σ h , B h ) and ( v h , k h , τ h , C h ) respectively, and use ξ − h to denote ( u − h , j − h , σ − h , B − h ) which is the solution of last iteration step (or initial guess). We assume u − h and B − h are given as known functions. For the initial guess, we assume (cid:107) u h (cid:107) , (cid:107) B h (cid:107) and (cid:107) j h (cid:107) = (cid:107)∇ h × B h (cid:107) arebounded. From the energy estimates, we know (cid:107) u − h (cid:107) , (cid:107) B − h (cid:107) and (cid:107)∇ h × B − h (cid:107) are bounded uniformlywith the iteration step.Define a ( ξ h , η h ) = R − e ( ∇ u h , ∇ v h ) + L ( u − h ; u h , v h ) + S ( j h , v h × B − h )+ SR − m ( ∇ × j h , C h ) − SR − m ( ∇ × σ h , C h ) + SR − m ( σ h , τ h ) − SR − m ( u h × B − h , τ h ) + S ( j h , k h ) − SR − m ( B h , ∇ × k h ) . The variational form with general right hand sides can be written as:
Problem 3.
Given ξ − h ∈ X h , θ = ( f , l , g , h ) ∈ X ∗ h , ψ = ( m, z ) ∈ Y ∗ h . Find ( u h , j h , σ h , B h , p h , r h ) ∈ X h × Y h , such that a ( ξ h , η h ) + b ( η h , x h ) = (cid:104) θ , η h (cid:105) , ∀ η h ∈ X h , ( . ) b ( ξ h , y h ) = (cid:104) ψ , y h (cid:105) , ∀ y h ∈ Y h . ( . ) Here (cid:104) θ , η h (cid:105) := (cid:104) f , v h (cid:105) + (cid:104) l , k h (cid:105) + (cid:104) g , τ h (cid:105) + (cid:104) h , C h (cid:105) , and (cid:104) ψ , y h (cid:105) := (cid:104) m, q h (cid:105) + (cid:104) z, s h (cid:105) . Problem 3 is equivalent to Algorithm 1 when u − = u n − , B − = B n − and l , g , h , m, z = 0 .We give the main theorem of well-posedness of the Picard iteration scheme: Theorem 6. (Well-posedness of Picard iterations)
There exists unique ( u h , j h , σ h , B h , p h , r h ) solving Problem 3, and the solution satisfies: (cid:107) ( u h , j h , σ h , B h ) (cid:107) X + (cid:107) ( p h , r h ) (cid:107) Y ≤ C ( (cid:107) ( f , l , g , h ) (cid:107) X ∗ + (cid:107) ( m, z ) (cid:107) Y ∗ ) , where C only depends on the domain, (cid:107) u − h (cid:107) and (cid:107) B − h (cid:107) d . Remark 1. If u − h and B − h are from the iterative scheme Algorithm 1, (cid:107) u − h (cid:107) and (cid:107) B − h (cid:107) d are uniformlybounded by known data, from the energy estimate ( . ) . Next we focus on the proof of this theorem. Similar to the nonlinear problem, we first formallyeliminate the variable j h by ∇ h × B h , and formally eliminate σ h to get a system with u h , B h and p h , s h as the variables (Problem 4 below). Boundedness and inf-sup condition of the bilinear form b ( · , · ) arealso similar to the nonlinear problem. Finally, we use the coercivity of the bilinear form ˜ a ( ˜ ξ − h ; · , · ) on ˜ X h to get the well-posedness of the Picard iterations. Problem 4.
Given ˜ ξ − h ∈ ˜ X h and ˜ θ = ( ˜ f , ˜ h ) ∈ ˜ X ∗ h , ˜ ψ = ( m, z ) ∈ Y ∗ h , find ˜ ξ h ∈ ˜ X h , x h ∈ Y h , suchthat ˜ a ( ˜ ξ − h ; ˜ ξ h , ˜ η h ) + ˜ b ( ˜ η h , x h ) = (cid:104) ˜ θ , ˜ η h (cid:105) , ∀ ˜ η h ∈ ˜ X h , ( . ) ˜ b ( ˜ ξ h , y h ) = (cid:104) ˜ ψ , y h (cid:105) , ∀ y h ∈ Y h . ( . ) where (cid:104) ˜ f , v h (cid:105) := (cid:104) f , v h (cid:105) − (cid:104) l , P ( v h × B − h ) (cid:105) , (cid:104) ˜ h , C h (cid:105) := (cid:104) h , C h (cid:105) − R − m (cid:104) l , ∇ h × C h (cid:105) + (cid:104) g , ∇ h × C h (cid:105) . In what follows we use (cid:107)·(cid:107) c ∗ to denote the dual norm of H h (curl , Ω) (with norm (cid:107)·(cid:107) c ): (cid:107) l (cid:107) c ∗ := sup F h ∈ H h (curl , Ω) (cid:104) l , F h (cid:105)(cid:107) F h (cid:107) c . To see ˜ f and ˜ h are bounded linear operators, we note the basic estimates: (cid:104) l , P ( v h × B − h ) (cid:105) ≤ (cid:107) l (cid:107) c ∗ (cid:107) P ( v h × B − h ) (cid:107) c ≤ (cid:107) l (cid:107) c ∗ (cid:107) v h × B − h (cid:107) (cid:46) (cid:107) l (cid:107) c ∗ (cid:107) v h (cid:107) (cid:107)∇ h × B − h (cid:107) , and (cid:104) l , ∇ h × C h (cid:105) ≤ (cid:107) l (cid:107) c ∗ (cid:107)∇ h × C h (cid:107) c ≤ (cid:107) l (cid:107) c ∗ (cid:107) C h (cid:107) d , (cid:104) g , ∇ h × C h (cid:105) ≤ (cid:107) g (cid:107) c ∗ (cid:107)∇ h × C h (cid:107) c ≤ (cid:107) g (cid:107) c ∗ (cid:107) C h (cid:107) d . In the following discussion, we will use the Riesz representation l , g ∈ H h (curl , Ω) of l , g ∈ H h (curl , Ω) ∗ which are defined by ( g , τ h ) := (cid:104) g , τ h (cid:105) , ∀ τ h ∈ H h (curl , Ω) , and ( l , k h ) := (cid:104) l , k h (cid:105) , ∀ k h ∈ H h (curl , Ω) . Note that (cid:107) g (cid:107) c = (cid:107) g (cid:107) c ∗ and (cid:107) l (cid:107) c = (cid:107) l (cid:107) c ∗ .For the relation between Problem 3 and Problem 4, we have: TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 17
Lemma 7. If ( u h , B h , p h , r h ) solves Problem 4 and (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d + (cid:107) p h (cid:107) + (cid:107) r h (cid:107) ≤ c ( (cid:107) ˜ f (cid:107) − + (cid:107) ˜ h (cid:107) H h (div) ∗ + (cid:107) ( m, z ) (cid:107) Y ∗ ) , then ( u h , j h , σ h , B h , p h , r h ) := ( u h ,R − m ∇ h × B h + S − l , P ( u h × B − h ) + S − R m g , B h , p h , r h ) ∈ X h × Y h solves Problem 3, and there exists a positive constant c , depending on c and (cid:107) B − h (cid:107) d such that (cid:107) ( u h , j h , σ h , B h ) (cid:107) X + (cid:107) ( p h , r h ) (cid:107) Y ≤ c (cid:0) (cid:107) ( f , l , g , h ) (cid:107) X ∗ + (cid:107) ( m, z ) (cid:107) Y ∗ (cid:1) . ( . ) On the other hand, if ( u h , j h , σ h , B h , p h , r h ) solves Problem 3, then ( u h , B h , p h , r h ) solves Problem4.Proof. In Problem 3, we take v h , k h , C h , q h , s h = 0 in ( . ) to see σ h = P ( u h × B − h ) + S − R m g , ( . )and take v h , τ h , C h , q h , s h = 0 in ( . ) to get j h = R − m ∇ h × B h + S − l . ( . )If ( u h , B h , p h , r h ) solves Problem 4, and (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d + (cid:107) p h (cid:107) + (cid:107) r h (cid:107) ≤ c (cid:16) (cid:107) ˜ f (cid:107) − + (cid:107) ˜ h (cid:107) H h (div) ∗ + (cid:107) ( m, z ) (cid:107) Y ∗ (cid:17) , it is easy to see from ( . ) and ( . ) that ( u h , R − m ∇ h × B h + S − l , P ( u h × B − h )+ S − R m g , B h , p h , r h ) solves Problem 3, and R − m (cid:107)∇ h × B h (cid:107) c + (cid:107) B h (cid:107) d = (cid:107) B h (cid:107) + (cid:107)∇ · B h (cid:107) +(1 + R − m ) (cid:107)∇ h × B h (cid:107) (cid:46) (cid:107) B h (cid:107) d , (cid:107) P ( u h × B − h ) (cid:107)≤ (cid:107) u h × B − h (cid:107) (cid:46) (cid:107) u h (cid:107) (cid:107)∇ h × B − h (cid:107) . This implies ( . ).On the other hand, solution of Problem 3 also solves Problem 4 by substituting ( . ) and ( . ) into( . ). (cid:3) Once the well-posedness of Problem 4 is established, the first part of Lemma 7 will imply existenceand stability of the original Problem 3, and the second part will imply the uniqueness. Hence it sufficesto prove well-posedness of Problem 4 under the norm (cid:107)·(cid:107) A (( . )).Similar to the nonlinear case, we have Lemma 8. (Boundedness) ˜ a ( ˜ ξ − h ; · , · ) is a bounded bilinear form on ˜ X h with respect to (cid:107)·(cid:107) A ( ( . ) ) We note that the bound depends on the domain and (cid:107) u − h (cid:107) , , (cid:107)∇ h × B − h (cid:107) . By the energy estimates,we know these terms are bounded by known data.The boundedness and inf-sup condition of b ( · , · ) are the same as the nonlinear problem (Lemma 5,Lemma 6).Next we show the coercivity of ˜ a ( ˜ ξ − h ; · , · ) on ˜ X h : Lemma 9.
There exists a positive constant α such that ˜ a ( ˜ ξ − h ; ˜ ξ h , ˜ ξ h ) ≥ α ( (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d ) , ∀ ˜ ξ h ∈ ˜ X h . Proof.
Taking v h = u h , C h = B h , ˜ a ( ˜ ξ − h ; ˜ ξ h , ˜ ξ h ) = R − e (cid:107)∇ u h (cid:107) + SR − m (cid:107)∇ h × B h (cid:107) . From Poincar´e’s inequality (Lemma 4) and ∇ · B h = 0 on ˜ X h : (cid:107) B h (cid:107) (cid:46) (cid:107)∇ h × B h (cid:107) . Hence (cid:107) B h (cid:107) d (cid:46) (cid:107)∇ h × B h (cid:107) , and there exists a positive constant α which only depends on the domain and R e , R m , S such that ˜ a ( ˜ ξ − h ; ˜ ξ h , ˜ ξ h ) ≥ α ( (cid:107) u h (cid:107) + (cid:107) B h (cid:107) d ) . (cid:3) From Lemma 8, Lemma 6 and Lemma 9, we have proved the well-posedness of Problem 4. FromLemma 7, this shows the well-posedness of Problem 3, and hence Algorithm 1 as a special case.5. C
ONVERGENCE ANALYSIS
Convergence of Picard iterations.
There is a general argument to prove the convergence of Picarditerations under the condition of small data, which guarantees the uniqueness of the nonlinear scheme(c.f. Girault and Raviart [11] Chapter IV, Remark 1.3; Gunzburger et al. [12] Proposition 7.1). Sincewe have established the boundedness and coercivity of the nonlinear variational form, the convergenceof Picard iteration scheme proposed in this paper can be analyzed in the same way, and a comparableresult holds. However we note that in the condition obtained in this way, the coupling number S cannotbe arbitrarily small, which seems to be contrary to the physical intuition. For example, in Gunzburger etal. [12], when we assume that the boundary data is zero, the criterion ((4.26) of [12]) is reduced to (cid:107) f (cid:107) − < S √ γ (cid:16) min (cid:16) k SR e , k R m (cid:17)(cid:17) max (cid:16) S , √ R m (cid:17) . ( . )Here we have used the notation in ( . ), with a correspondence to the original notation in [12]: N = S , M = √ SR e , F = S − f , where F is the right hand side in [12]. Furthermore, here γ , k and k are positive constants in the Sobolev imbedding and the Poincar´e’s inequality of velocity and magneticfields. Now it is easy to see that in ( . ), S cannot be arbitrarily small for fixed R e , R m and f (cid:54) = . Thecondition (2.16) in Sch¨otzau [18] is similar.Therefore in this section, we use a different approach and directly prove the convergence of the Picarditerations by contraction. As a result, we will see that the small data condition (( . ) below) will onlycontain R e and R m , but not S . The (discrete) energy law is crucial in the argument below as an a prioriestimate.A similar argument also holds on the continuous level with minor modifications. We omit the sub-script “ h ” in this section. Theorem 7.
The Picard iteration scheme (Algorithm 1) converges when (cid:107) f (cid:107) − ≤ (cid:0) C R e + 4 C R e R m (cid:1) − , ( . ) where C and C , depending only on the domain, are positive constants related to the Sobolev imbeddingand regularity estimates of H h (div) functions given in ( . ) and ( . ) . The above conditions are satisfied when the data (cid:107) f (cid:107) − is small relative to R − e and R − m . TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 19
Proof.
By the standard theory of mixed methods, it suffices to consider the convergence in X h := { η h ∈ X h : b ( η h , y h ) = 0 , ∀ y h ∈ Y h } . The equation of the n -th step can be written as L ( u n − ; u n , v ) + R − e ( ∇ u n , ∇ v ) − S ( j n × B n − , v ) = (cid:104) f , v (cid:105) , ( . ) − ( u n × B n − , ∇ h × C ) + R − m ( ∇ h × B n , ∇ h × C ) = 0 . ( . )The ( n − -th step is similarly written as L ( u n − ; u n − , v ) + R − e ( ∇ u n − , ∇ v ) − S ( j n − × B n − , v ) = (cid:104) f , v (cid:105) , ( . ) − ( u n − × B n − , ∇ h × C ) + R − m ( ∇ h × B n − , ∇ h × C ) = 0 . ( . )Define the errors e nu := u n − u n − , e nB := B n − B n − , e nj := j n − j n − . From the equation j n = R − m ∇ h × B n , we have e nj = R − m ∇ h × e nB .Subtracting ( . )-( . ) from the n -step equation ( . )-( . ), we get the error equation: (cid:0) ( u n − · ∇ ) e nu , v (cid:1) + 12 (cid:0) ( e n − u · ∇ ) u n − , v (cid:1) − (cid:0) ( u n − · ∇ ) v , e nu (cid:1) − (cid:0) ( e n − u · ∇ ) v , u n − (cid:1) + 1 R e ( ∇ e nu , ∇ v ) + S (cid:0) B n − × e nj , v (cid:1) + S (cid:0) e n − B × j n − , v (cid:1) = 0 , − ( e nu × B n − , ∇ h × C ) − ( u n − × e n − B , ∇ h × C ) + R − m ( ∇ h × e nB , ∇ h × C ) = 0 . Multiplying the second equation by SR − m , adding the above two equations and taking v = e nu , C = e nB yield (cid:0) ( e n − u · ∇ ) u n − , e nu (cid:1) − (cid:0) ( e n − u · ∇ ) e nu , u n − (cid:1) + R − e ( ∇ e nu , ∇ e nu ) + S ( e n − B × j n − , e nu ) ( . ) − SR − m (cid:0) u n − × e n − B , ∇ h × e nB (cid:1) + SR − m ( ∇ h × e nB , ∇ h × e nB ) = 0 . From the energy estimates ( . ), we know (cid:107)∇ u n (cid:107)≤ R e (cid:107) f (cid:107) − , and (cid:107) j n (cid:107)≤ (cid:18) R e S (cid:19) (cid:107) f (cid:107) − , which hold for all n > .Then we have the estimates for the nonlinear terms: (cid:12)(cid:12)(cid:12)(cid:12) (cid:0) ( e n − u · ∇ ) u n − , e nu (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) e n − u (cid:107) , (cid:107)∇ u n − (cid:107)(cid:107) e nu (cid:107) , ≤ C R e (cid:107) f (cid:107) − (cid:107)∇ e n − u (cid:107)(cid:107)∇ e nu (cid:107)≤ R e (cid:107)∇ e n − u (cid:107) + 12 C R e (cid:107) f (cid:107) − (cid:107)∇ e nu (cid:107) , (cid:12)(cid:12)(cid:12)(cid:12) (cid:0) ( e n − u · ∇ ) e nu , u n − (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) u n − (cid:107) , (cid:107)∇ e nu (cid:107)(cid:107) e n − u (cid:107) , ≤ C R e (cid:107) f (cid:107) − (cid:107)∇ e n − u (cid:107)(cid:107)∇ e nu (cid:107)≤ R e (cid:107)∇ e n − u (cid:107) + 12 C R e (cid:107) f (cid:107) − (cid:107)∇ e nu (cid:107) , (cid:12)(cid:12) S ( e n − B × j n − , e nu ) (cid:12)(cid:12) ≤ SC (cid:107)∇ h × e n − B (cid:107)(cid:107) j n − (cid:107)(cid:107)∇ e nu (cid:107)≤ SC R m (cid:107) j n − (cid:107)(cid:107) e n − j (cid:107)(cid:107)∇ e nu (cid:107)≤ SC R m (cid:18) R e S (cid:19) (cid:107) f (cid:107) − (cid:107) e n − j (cid:107)(cid:107)∇ e nu (cid:107)≤ S (cid:107) e n − j (cid:107) +2 R e C R m (cid:107) f (cid:107) − (cid:107)∇ e nu (cid:107) , and (cid:12)(cid:12) SR − m ( u n − × e n − B , ∇ h × e nB ) (cid:12)(cid:12) ≤ SR − m C (cid:107)∇ h × e n − B (cid:107)(cid:107)∇ u n − (cid:107)(cid:107)∇ h × e nB (cid:107)≤ SC R e R m (cid:107) f (cid:107) − (cid:107) e n − j (cid:107)(cid:107) e nj (cid:107)≤ S (cid:107) e n − j (cid:107) +2 SR m C R e (cid:107) f (cid:107) − (cid:107) e nj (cid:107) . Combining the above estimates with ( . ), we have (cid:18) R e − C R e (cid:107) f (cid:107) − − R e C R m (cid:107) f (cid:107) − (cid:19) (cid:107)∇ e nu (cid:107) + (cid:0) S − R m SC R e (cid:107) f (cid:107) − (cid:1) (cid:107) e nj (cid:107) ≤ R e (cid:107)∇ e n − u (cid:107) + 14 S (cid:107) e n − j (cid:107) . We define the energy functional to be E n := 12 R e (cid:107)∇ e nu (cid:107) + 12 S (cid:107) e nj (cid:107) . Therefore when R e ≥ C R e (cid:107) f (cid:107) − +2 C R e R m (cid:107) f (cid:107) − , and S ≥ R m SC R e (cid:107) f (cid:107) − , i.e. when ( . ) holds, we have E n ≤ E n − . This implies that ( u n , B n ) converges to some ( u , B ) in the norm defined by R − e (cid:107)∇ u n (cid:107) + SR − m (cid:107)∇ h × B n (cid:107) . Combined with the continuity of the trilinear form, we can take the limit and ( u , B ) is a solution of thenonlinear Problem 1.From the inf-sup condition of the velocity-pressure pair, we also have the convergence of the pressure p n . (cid:3) TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 21
Convergence of the finite element method.
We prove the convergence of the nonlinear finiteelement scheme. In the discussions below, we deal with the reduced form of the finite element schemewith variables ( u h , B h , p h , r h ) (Problem 2), then recover j h and σ h from these variables.As a routine approach for mixed methods, the proof below consists of several steps. We first subtractthe finite element solution from the true solution to obtain certain orthogonality (( . )). Then we insertan arbitrary discrete function to the orthogonality equation to get ( . ). Combining with triangularinequalities, numerical errors can be bounded by the difference of the true solution and the discretefunctions inserted above. Such an estimate is usually called quasi-orthogonality (Theorem 8). Then thefinal estimate (( . )) follows from polynomial approximation results.The analysis below also contains some new features compared with conventional error estimates formixed methods. The finite element scheme involves the discrete adjoint operator ∇ h × , which can onlybe defined for finite element functions. Therefore it is no wonder that the consistency error (cid:107)∇ × B −∇ h × B I (cid:107) will come into our analysis. Moreover, in the analysis for the nonlinear problem, we willfrequently use the key technical results established in Section § a priori estimate for bothnumerical and true solutions. Combining these key estimates and small source assumptions, which arecommon for nonlinear problems, we obtain the desired results.We begin detailed analysis by discovering the orthogonality. Subtracting the true solution of ( . )from the variational form ( . ), we have for any ( v h , C h ) ∈ ˜ X h , ( q h , s h ) ∈ Y h , [(( u h − u ) · ∇ u h , v h ) + (( u · ∇ ) ( u h − u ) , v h ) − (( u h · ∇ ) v h , u h − u ) − (( u h − u ) · ∇ v h , u )] + R − e ( ∇ ( u h − u ) , ∇ v h ) − ( p h − p, ∇ · v h ) − SR − m (( ∇ h × B h ) × B h , v h ) + SR − m (( ∇ × B ) × B , v h ) = 0 , − SR − m ( u h × B h , ∇ h × C h ) + SR − m ( ∇ × ( u × B ) , C h )+ SR − m ( ∇ × ∇ h × B h − ∇ × ∇ × B , C h ) + ( r h − r, ∇ · C h ) = 0 , ( ∇ · ( u h − u ) , q h ) = 0 , ( ∇ · ( B h − B ) , s h ) = 0 . ( . )We assume that ( u I , B I ) ∈ ˜ X h and ( p I , r I ) ∈ Y h are arbitrary discrete functions. Inserting ( u I , B I ) , ( p I , r I ) into ( . ), we get: for any ( v h , C h ) ∈ ˜ X h , ( q h , s h ) ∈ Y h , [(( u h − u I ) · ∇ u h , v h ) + (( u · ∇ ) ( u h − u I ) , v h ) − (( u h · ∇ ) v h , u h − u I ) − (( u h − u I ) · ∇ v h , u )] + R − e ( ∇ ( u h − u I ) , ∇ v h ) − ( p h − p I , ∇ · v h ) − SR − m ( ∇ h × ( B h − B I ) × B h , v h ) − SR − m (( ∇ × B ) × ( B h − B I ) , v h )= [(( u − u I ) · ∇ ) u h , v h ) + (( u · ∇ )( u − u I ) , v h ) − (( u h · ∇ ) v h , u − u I ) − (( u − u I ) · ∇ v h , u )] + R − e ( ∇ ( u − u I ) , ∇ v h ) − ( p − p I , ∇ · v h )+ SR − m (( ∇ h × B I − ∇ × B ) × B h , v h ) + SR − m (( ∇ × B ) × ( B I − B ) , v h ) , − SR − m (( u h − u I ) × B h , ∇ h × C h ) − SR − m ( u × ( B h − B I ) , ∇ h × C h )+ SR − m ( ∇ h × ( B h − B I ) , ∇ h × C h ) + ( r h − r I , ∇ · C h )= − SR − m (( u − u I ) × B h , ∇ h × C h ) + SR − m ( u × ( B I − B ) , ∇ h × C h )+ SR − m ( ∇ × ( ∇ × B − ∇ h × B I ) , C h ) + ( r − r I , ∇ · C h ) − SR − m ( ∇ × (id − P )( u × B ) , C h ) , ( ∇ · ( u h − u I ) , q h ) = ( ∇ · ( u − u I ) , q h ) , ( ∇ · ( B h − B I ) , s h ) = ( ∇ · ( B − B I ) , s h ) . Here we have used the identity ( ∇ × ( u × B ) , C h ) = ( ∇ × (id − P )( u × B ) , C h ) + ( u × B , ∇ h × C h ) . Adding the first two equations together, we can write the above system as R − e ( ∇ ( u h − u I ) , ∇ v h ) − ( p h − p I , ∇ · v h ) − SR − m ( ∇ h × ( B h − B I ) × B h , v h ) − SR − m (( u h − u I ) × B h , ∇ h × C h ) + SR − m ( ∇ h × ( B h − B I ) , ∇ h × C h )+( r h − r I , ∇ · C h ) + G ( u h , B h , u , B ; u h − u I , B h − B I ; v h , C h )= H ( u h , B h , u , B ; u − u I , B − B I , p − p I , r − r I ; v h , C h ) + SR − m (( ∇ h × B I − ∇ × B ) × B h , v h )+ SR − m ( ∇ × ( ∇ × B − ∇ h × B I ) , C h ) − SR − m ( ∇ × (id − P )( u × B ) , C h ) , ( ∇ · ( u h − u I ) , q h ) = ( ∇ · ( u − u I ) , q h ) , ( ∇ · ( B h − B I ) , s h ) = ( ∇ · ( B − B I ) , s h ) , ( . )where G ( u h , B h , u , B ; u h − u I , B h − B I ; v h , C h ) = 12 [(( u h − u I ) · ∇ u h , v h )+ (( u · ∇ ) ( u h − u I ) , v h ) − (( u h · ∇ ) v h , u h − u I ) − (( u h − u I ) · ∇ v h , u )] − SR − m (( ∇ × B ) × ( B h − B I ) , v h ) − SR − m ( u × ( B h − B I ) , ∇ h × C h ) , and H ( u h , B h , u , B ; u − u I , B − B I , p − p I , r − r I ; v h , C h ) = 12 [(( u − u I ) · ∇ ) u h , v h )+ (( u · ∇ )( u − u I ) , v h ) − (( u h · ∇ ) v h , u − u I ) − (( u − u I ) · ∇ v h , u )]+ SR − m (( ∇ × B ) × ( B I − B ) , v h ) − SR − m (( u − u I ) × B h , ∇ h × C h )+ SR − m ( u × ( B I − B ) , ∇ h × C h ) + R − e ( ∇ ( u − u I ) , ∇ v h ) − ( p − p I , ∇ · v h ) + ( r − r I , ∇ · C h ) . Thanks to the energy law and the key estimate for the regularity of B h (Theorem 1), norms (cid:107) u h (cid:107) , (cid:107) B h (cid:107) d , (cid:107) u (cid:107) , (cid:107) B (cid:107) , can be bounded by the source (cid:107) f (cid:107) − . Therefore H and G are bounded bilinear forms with coefficientswhich can be controlled by (cid:107) f (cid:107) − . Specifically, we have the boundedness | G ( u h , B h , u , B ; u h − u I , B h − B I ; v h , C h ) |≤ Γ (cid:0) (cid:107)∇ ( u − u I ) (cid:107) + (cid:107)∇ h × ( B h − B I ) (cid:107) (cid:1) / (cid:0) (cid:107)∇ v h (cid:107) + (cid:107)∇ h × C h (cid:107) (cid:1) / , and | H ( u h , B h , u , B ; u − u I , B − B I , p − p I , r − r I ; v h , C h ) | ( . ) ≤ Γ (cid:0) (cid:107) u − u I (cid:107) + (cid:107) B − B I (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) (cid:1) / (cid:107) ( v h , C h ) (cid:107) ˜ X , where Γ = C ( (cid:107)∇ u h (cid:107) + (cid:107)∇ u (cid:107) ) + SR − m C C ( (cid:107)∇ × B (cid:107) + (cid:107)∇ u (cid:107) ) , and Γ = C ( (cid:107)∇ u h (cid:107) + (cid:107)∇ u (cid:107) ) + SR − m C C (cid:107)∇ h × B h (cid:107) + SR − m C (cid:107)∇ × B (cid:107) , + SR − m (cid:107) u (cid:107) , ∞ +2 + R − e . From the energy law, we have (cid:107)∇ u h (cid:107)≤ R e (cid:107) f (cid:107) − , (cid:107)∇ u (cid:107)≤ R e (cid:107) f (cid:107) − , and (cid:107)∇ × B (cid:107)≤ (cid:114) R e R m S (cid:107) f (cid:107) − , (cid:107)∇ h × B h (cid:107)≤ (cid:114) R e R m S (cid:107) f (cid:107) − . TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 23
Therefore Γ ≤ (cid:16) C R e + √ / C C (cid:112) R e S + SR − m C C R e (cid:17) (cid:107) f (cid:107) − . There are three remaining terms on the right hand side of ( . ), i.e. I := SR − m (( ∇ h × B I − ∇ × B ) × B h , v h ) ,I := SR − m ( ∇ × ( ∇ × B − ∇ h × B I ) , C h ) , and I := SR − m ( ∇ × (id − P )( u × B ) , C h ) . Next, we estimate these three terms. The following lemma gives an estimate for the consistency term ∇ × B − ∇ h × B I . An analogous 2D version can be found in [7]. Lemma 10.
We have the estimate for the consistency of the discrete adjoint operator (cid:107)∇ × B − ∇ h × B I (cid:107) (cid:46) (cid:107) (id − P ) ∇ × B (cid:107) + h − (cid:107) B − B I (cid:107) . Proof.
We recall that P denotes the L projection to H h (curl , Ω) . We have (cid:107)∇ × B − ∇ h × B I (cid:107) = (cid:107)∇ × B − P ( ∇ × B ) + P ( ∇ × B ) − ∇ h × B I (cid:107)≤ (cid:107) (id − P ) ∇ × B (cid:107) + (cid:107) P ( ∇ × B ) − ∇ h × B I (cid:107) . For the second term, we use a dual estimate: for any φ h ∈ H h (curl , Ω) , ( P ( ∇ × B ) − ∇ h × B I , φ h ) = ( ∇ × B − ∇ h × B I , φ h )= ( B − B I , ∇ × φ h ) ≤ (cid:107) B − B I (cid:107)(cid:107)∇ × φ h (cid:107) (cid:46) h − (cid:107) B − B I (cid:107)(cid:107) φ h (cid:107) . This implies that (cid:107) P ( ∇ × B ) − ∇ h × B I (cid:107) (cid:46) h − (cid:107) B − B I (cid:107) and the desired result follows. (cid:3) Lemma 10 implies the estimate for I : | I | (cid:46) (cid:0) (cid:107) (id − P ) ∇ × B (cid:107) + h − (cid:107) B − B I (cid:107) (cid:1) (cid:107) B h (cid:107) d (cid:107) v h (cid:107) . We turn to the estimate for I : ( ∇ × ∇ × B , C h ) − ( ∇ h × B I , ∇ h × C h )= ( ∇ × ( P + id − P ) ∇ × B , C h ) − ( ∇ h × B I , ∇ h × C h )= ( ∇ × B − ∇ h × B I , ∇ h × C h ) + ( ∇ × (id − P ) ∇ × B , C h ) . Using Lemma 10 again, we get | I | (cid:46) (cid:0) h − (cid:107) B − B I (cid:107) + (cid:107) (id − P ) ∇ × B (cid:107) + (cid:107)∇ × (id − P ) ∇ × B (cid:107) (cid:1) (cid:107) C h (cid:107) d . Moreover, we have a straightforward estimate for I : | I | ≤ (cid:107)∇ × (id − P ) ( u × B ) (cid:107)(cid:107) C h (cid:107) . For any B ∈ H (div , Ω) , we define (cid:107) B (cid:107) := (cid:107) B (cid:107) + (cid:107)∇ · B (cid:107) . Lemma 11.
Assume that (cid:107) f (cid:107) − is sufficiently small. There exists C > depending on Ω , (cid:107) u (cid:107) , ∞ and (cid:107) B (cid:107) , , such that for any ( u I , B I ) ∈ ˜ X h , ( p I , r I ) ∈ Y h , (cid:107) u h − u I (cid:107) + (cid:107) B h − B I (cid:107) d + (cid:107) p h − p I (cid:107) + (cid:107) r h − r I (cid:107) ≤ C ( (cid:107) u − u I (cid:107) + (cid:107) B − B I (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) + h − (cid:107) B − B I (cid:107) + (cid:107) (id − P ) ∇ × B (cid:107) + (cid:107)∇ × (id − P ) ∇ × B (cid:107) + (cid:107)∇ × (id − P )( u × B ) (cid:107) ) . Proof.
Given ( u , B , p, r ) and ( u I , B I , p I , r I ) , the system ( . ) can be seen as equations for ( u h − u I , B h − B I , p h − p I , r h − r I ) . Compared with the nonlinear discrete system which we have analyzed,i.e. Problem 2, a new term G appears on the left hand side and the fluid convection term has beenabsorbed into G .We assume that (cid:107) f (cid:107) − ≤ min (cid:8) / R − e , / SR − m (cid:9) (cid:16) C R e + √ / C C (cid:112) R e S + SR − m C C R e (cid:17) − . ( . )A direct consequence ( . ) is R e ≤ / − and R m ≤ (cid:0) / S Γ − (cid:1) / . Then we have | G ( u h , B h , u , B ; v h , C h ; v h , C h ) | ≤ R − e (cid:107)∇ v h (cid:107) + 12 SR − m (cid:107)∇ h × C h (cid:107) , ∀ ( v h , C h ) ∈ ˜ X h , then the left hand side A ( w h , G h ; v h , C h ) := R − e ( ∇ w h , ∇ v h ) − SR − m (( ∇ h × G h ) × B h , v h ) − SR − m ( w h × B h , ∇ h × C h ) + SR − m ( ∇ h × G h , ∇ h × C h ) + G ( u h , B h , u , B ; w h , G h ; v h , C h ) defines a bounded coercive bilinear form for fixed u h , B h , u and B . The boundedness constant dependson (cid:107) u h (cid:107) , (cid:107) u (cid:107) , (cid:107)∇ h × B h (cid:107) and (cid:107)∇ × B (cid:107) , , which further depend on (cid:107) f (cid:107) − .For the right hand sides, H ( u h , B h , u , B ; u − u I , B − B I , p − p I , r − r I ; · ) can be regarded as abounded linear functional on ˜ X h for fixed u h , B h , u , B , u I , B I , and the dual norm can be bounded by Γ (cid:0) (cid:107) u − u I (cid:107) + (cid:107) B − B I (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) (cid:1) / , due to ( . ). Moreover, given u h − u I and B h − B I , ( ∇ · ( u h − u I ) , q h ) and ( ∇ · ( B h − B I ) , s h ) are bounded linear functionals on Q h and L h respectively, with dual norms (cid:107)∇ · ( u h − u I ) (cid:107) and (cid:107)∇ · ( B h − B I ) (cid:107) . From the estimates for I , I and I , dual norms of these three terms can be bounded by max (cid:8) (cid:107)∇ × (id − P ) ∇ × B (cid:107) + h − (cid:107) B − B I (cid:107) + (cid:107) (id − P ) ∇ × B (cid:107) , (cid:107)∇ × (id − P )( u × B ) (cid:107) (cid:9) , up to a positive constant.From a general argument of the Brezzi theory, we see that the norms of the solution of ( . ), i.e., (cid:107) ( u h − u I , B h − B I ) (cid:107) X + (cid:107) ( p h − p I , r h − r I ) (cid:107) Y can be bounded by the dual norm of the right hand side. This completes the proof. (cid:3) Combining triangular inequalities and the estimate (cid:107)∇ × B − ∇ h × B h (cid:107)≤ (cid:107)∇ × B − ∇ h × B I (cid:107) + (cid:107)∇ h × ( B I − B h ) (cid:107) , we obtain the following quasi-optimal estimate. Theorem 8.
Assume that the condition ( . ) holds. There exists a generic positive constant C > depending on Ω , (cid:107) f (cid:107) − , (cid:107) u (cid:107) , ∞ and (cid:107) B (cid:107) , , such that for any ( u I , B I ) ∈ ˜ X h , ( p I , r I ) ∈ Y h , (cid:107) u − u h (cid:107) + (cid:107) B − B h (cid:107) + (cid:107)∇ × B − ∇ h × B h (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) ≤ C ( (cid:107) u − u I (cid:107) + (cid:107) B − B I (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) + h − (cid:107) B − B I (cid:107) + (cid:107) (id − P ) ∇ × B (cid:107) + (cid:107)∇ × (id − P ) ∇ × B (cid:107) + (cid:107)∇ × (id − P )( u × B ) (cid:107) ) . ( . ) TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 25
We remark that (cid:107)∇ × B − ∇ h × B h (cid:107) = (cid:107) R m ( j − j h ) (cid:107) yields an L error estimate for the currentdensity j .The last step is to estimate the convergence order based on the polynomial approximation theory. Werecall the following approximation result. Lemma 12.
Assume that H h (curl , Ω) contains piecewise polynomials of degree s . Then the L projec-tion P satisfies the approximation property (cid:107) φ − P φ (cid:107) + h (cid:107)∇ × ( φ − P φ ) (cid:107) (cid:46) h s +1 (cid:107) φ (cid:107) s +1 , ∀ φ ∈ H s +1 (Ω) . The proof is almost the same as the classical result of L projections for Lagrange elements. Forcompleteness, we include the proof here. Proof.
Let Π h curl be a bounded interpolation operator to H h (curl , Ω) , for example, defined in [10]. Thenwe have (cid:107)∇ × ( φ − P φ ) (cid:107) ≤ (cid:13)(cid:13) ∇ × (cid:0) φ − Π h curl φ (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) ∇ × Π h curl ( φ − P φ ) (cid:13)(cid:13) . For the first term on the right hand side, (cid:13)(cid:13) ∇ × (cid:0) φ − Π h curl φ (cid:1)(cid:13)(cid:13) (cid:46) h s (cid:107) φ (cid:107) s +1 . For the second, we use the inverse estimate to get (cid:13)(cid:13) ∇ × Π h curl ( φ − P φ ) (cid:13)(cid:13) (cid:46) h − (cid:13)(cid:13) Π h curl ( φ − P φ ) (cid:13)(cid:13) (cid:46) h s (cid:107) φ (cid:107) s +1 . This implies h (cid:107)∇ × ( φ − P φ ) (cid:107) (cid:46) h s +1 (cid:107) φ | s +1 .On the other hand, the approximation (cid:107) φ − P φ (cid:107) (cid:46) h s +1 (cid:107) φ (cid:107) s +1 follows directly from the property of the L projection operator. This completes the proof. (cid:3) In the following discussions, we assume that H h (curl , Ω) , H h (div , Ω) and L h (Ω) contain piecewisepolynomials of degree r , r and r respectively. From the construction of discrete de Rham complexes,we have r i = r i +1 or r i = r i +1 + 1 where i = 1 , . We assume that the approximation space V h forthe velocity contains piecewise polynomials of degree s u and the discrete pressure space Q h containspiecewise polynomials of degree s p .We estimate the projection error on the right hand side of ( . ) based on Lemma 12: (cid:107)∇ × (id − P )( u × B ) (cid:107) (cid:46) h r (cid:107) u × B (cid:107) r +1 , (cid:107)∇ × (id − P ) ∇ × B (cid:107) (cid:46) h r (cid:107)∇ × B (cid:107) r +1 , Consequently, we have (cid:107) u − u h (cid:107) + (cid:107) B − B h (cid:107) + (cid:107)∇ × B − ∇ h × B h (cid:107) + (cid:107) p − p I (cid:107) + (cid:107) r − r I (cid:107) ≤ C ( h s u (cid:107) u (cid:107) s u +1 + h s p +2 (cid:107) p (cid:107) s p +1 + h r (cid:107) B (cid:107) r +1 + h r ( (cid:107) u × B (cid:107) r +1 + (cid:107)∇ × B (cid:107) r +1 ) + h r +2 (cid:107) r (cid:107) r +1 ) . ( . )Based on the error estimate ( . ), we can get balanced errors by choosing finite elements such that r = r = r + 1 = s u = s p + 1 . One particular choice is to use BDM spaces for the magnetic field B , N´ed´elec spaces of the first kind for the electric field E . The pressure multiplier p and the magneticmultiplier r may be chosen to have the same order.The above analysis excludes the lowest order Raviart-Thomas element, but includes the case of thelowest order BDM element. We believe that this restriction is only technical but a more refined estimateis beyond the scope of this paper.
6. C
ONCLUDING REMARKS
In this paper we considered the mixed finite element discretizations of the stationary MHD system.Compared to the time-dependent system, the Gauss’s law of magnetic field is an independent equationwhich cannot be derived from the Faraday’s law. Therefore classical techniques of Lagrange multipliersare employed to impose the Gauss’s law. The structure-preserving discretization proposed in this paperfor the stationary MHD system preserves both the discrete energy law and most importantly the Gauss’slaw ∇ · B = 0 .We note that we can also use a formulation based on B and E , which is similar to the time-dependentcase studied in [14]. But the well-posedness of such a formulation can only be established when theReynolds number R e is assumed to be sufficiently small. To remove such an undesirable constraint, weproposed the new formulation using B and j as the variables. Such a formulation was partially motivatedby the fact that the energy is given in terms of (cid:107) j (cid:107) rather than (cid:107) E (cid:107) .These two formulations look similar. In the finite element discretization of both cases, we have j = E + P ( u × B ) (only one variable of E and j is explicitly used in one scheme). This is an equationin H h (curl , Ω) . The current density j and the electric field E differ by a nonlinear term, which isprojected to H h (curl , Ω) . But the resulting formulations are different due to the different treatmentsof the nonlinear term P ( u × B ) in the discretization of the Lorentz force term. We note that in theformulation proposed in [14], the Lorentz force term ( j , v × B ) is discretized as ( E + u × B , v × B ) . Whereas in the formulation proposed in this paper, the corresponding discretization is as ( E + P ( u × B ) , P ( v × B )) . It is easy to see that these two discretizations are indeed different.Similar differences can be also found at other places. A key point to get well-posedness is the cancella-tion of the symmetric nonlinear coupling terms. Under such a restriction, other parts of the schemes alsohave to be different according to the different Lorentz force terms. Indeed the energy estimates of thesetwo kinds of formulations have already shown the difference. The energy estimates of the formulationin [14] involve (cid:107) E + u × B (cid:107) , while the formulation in this paper involves (cid:107) j (cid:107) = (cid:107) E + P ( u × B ) (cid:107) .As a result of these differences, a careful analysis indicates that the well-posedness of the formulationproposed in this paper can be established without any assumption on the size of R e .A CKNOWLEDGEMENT
The authors would like to thank Mr. Juncai He, Prof. Ragnar Winther and Dr. Shuonan Wu forhelpful discussions, and the anonymous referees for valuable suggestions, which have greatly improvedthe quality of the paper. R
EFERENCES1. Douglas N Arnold, Richard S Falk, and Ragnar Winther,
Finite element exterior calculus, homological techniques, and appli-cations , Acta numerica (2006), 1.2. , Finite element exterior calculus: from Hodge theory to numerical stability , Bulletin of the American MathematicalSociety (2010), no. 2, 281–354.3. Daniele Boffi, Franco Brezzi, and Michel Fortin, Mixed Finite Element Methods and Applications , Springer, 2013.4. Alain Bossavit,
Computational electromagnetism: variational formulations, complementarity, edge elements , Academic Press,1998.
TRUCTURE-PRESERVING FINITE ELEMENT METHODS FOR STATIONARY MHD MODELS 27
5. Jeremiah U Brackbill and Daniel C Barnes,
The effect of nonzero ∇ · b on the numerical solution of the magnetohydrodynamicequations , Journal of Computational Physics (1980), no. 3, 426–430.6. Susanne C Brenner and Ridgway Scott, The mathematical theory of finite element methods , vol. 15, Springer Science &Business Media, 2008.7. Long Chen, Ming Wang, and Lin Zhong,
Convergence analysis of triangular MAC schemes for two dimensional Stokesequations , Journal of Scientific Computing (2014), 1–29.8. Snorre H Christiansen, Hans Z Munthe-Kaas, and Brynjulf Owren,
Topics in structure-preserving discretization , Acta Numer-ica (2011), 1–119.9. Wenlong Dai and Paul R Woodward, On the divergence-free condition and conservation laws in numerical simulations forsupersonic magnetohydrodynamical flows , The Astrophysical Journal (1998), no. 1, 317.10. Richard S Falk and Ragnar Winther,
Local bounded cochain projections , Mathematics of Computation (2014), no. 290,2631–2656.11. Vivette Girault and Pierre-Arnaud Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms , vol. 5,Springer Science & Business Media, 2012.12. Max D Gunzburger, Amnon J Meir, and Janet S Peterson,
On the existence, uniqueness, and finite element approximationof solutions of the equations of stationary, incompressible magnetohydrodynamics , Mathematics of Computation (1991),no. 194, 523–563.13. Ralf Hiptmair, Finite elements in computational electromagnetism , Acta Numerica (2002), no. July 2003, 237–339.14. Kaibo Hu, Yicong Ma, and Jinchao Xu, Stable finite element methods preserving ∇ · B = 0 exactly for MHD models ,Numerische Mathematik (2014), 1–26.15. Yicong Ma, Kaibo Hu, Xiaozhe Hu, and Jinchao Xu, Robust preconditioners for incompressible MHD models , Journal ofComputational Physics (2016), 721–746.16. Ming-Jiu Ni and Jun-Feng Li,
A consistent and conservative scheme for incompressible MHD flows at a low magnetic Reynoldsnumber. Part III: On a staggered mesh , Journal of Computational Physics (2012), no. 2, 281–298.17. Ming-Jiu Ni, Ramakanth Munipalli, Neil B. Morley, Peter Huang, and Mohamed a. Abdou,
A current density conservativescheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: On a rectangular collocated grid system ,Journal of Computational Physics (2007), no. 1, 174–204.18. Dominik Sch¨otzau,
Mixed finite element methods for stationary incompressible magneto–hydrodynamics , Numerische Math-ematik (2004), 771–800.19. Zhiyi Yang, Tao Zhou, Hongli Chen, and Ming-Jiu Ni,
Numerical study of MHD pressure drop in rectangular ducts withinsulating coatings , Fusion Engineering and Design (2010), no. 10-12, 2059–2064.B EIJING I NTERNATIONAL C ENTER FOR M ATHEMATICAL R ESEARCH , P
EKING U NIVERSITY , B
EIJING
HINA . E-mail address : [email protected] C ENTER FOR C OMPUTATIONAL M ATHEMATICS AND A PPLICATIONS AND D EPARTMENT OF M ATHEMATICS , T HE P ENN - SYLVANIA S TATE U NIVERSITY , U
NIVERSITY P ARK , PA 16802, USA.
E-mail address ::