Mixed Discontinuous Galerkin Finite Element Methods for Incompressible Magnetohydrodynamics by Xiaoxi Wei B.Sc., Inner Mongolia University, China, 2002 M.Sc., Inner Mongolia University, China, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mathematics) The University of British Columbia (Vancouver) May 2011 c© Xiaoxi Wei, 2011 Abstract We develop and analyze mixed discontinuous Galerkin finite element methods for the numerical approximation of incompressible magnetohydrodynamics problems. Incompressible magnetohydrodynamics is the area of physics that is concerned with the behaviour of electrically conducting, resistive, incompressible and viscous fluids in the presence of electromagnetic fields. It is modelled by a system of nonlinear partial differential equations, which couples the Navier-Stokes equations with the Maxwell equations. In the first part of this thesis, we introduce an interior penalty discontinuous Galerkin method for the numerical approximation of a linearized incompressible magnetohydrodynamics problem. The fluid unknowns are discretized with the dis- continuous Pk-Pk−1 element pair, whereas the magnetic variables are approxi- mated by discontinuous Pk-Pk+1 elements. Under minimal regularity assump- tions, we carry out a complete a priori error analysis and prove that the energy norm error is optimally convergent in the mesh size in general polyhedral domains, thus guaranteeing the numerical resolution of the strongest magnetic singularities in non-convex domains. In the second part of this thesis, we propose and analyze a new mixed dis- continuous Galerkin finite element method for the approximation of a fully non- linear incompressible magnetohydrodynamics model. The velocity field is now discretized by divergence-conforming Brezzi-Douglas-Marini elements, and the magnetic field by curl-conforming Nédélec elements. In addition to correctly cap- turing magnetic singularities, the method yields exactly divergence-free velocity approximations, and is thus energy-stable. We show that the energy norm error is convergent in the mesh size in possibly non-convex polyhedra, and derive slightly ii suboptimal a priori error estimates under minimal regularity and small data as- sumptions. Finally, in the third part of this thesis, we present two extensions of our dis- cretization techniques to time-dependent incompressible magnetohydrodynamics problems and to Stokes problems with nonstandard boundary conditions. All our discretizations and theoretical results are computationally validated through comprehensive sets of numerical experiments. iii Preface Most of the results presented in this thesis have been published in the following two papers: (a) P. Houston, D. Schötzau, and X. Wei. A mixed DG method for linearized incompressible magnetohydrodynamics. Journal of Scientific Computing, 40:281–314, 2009. (b) C. Greif, D. Li, D. Schötzau, and X.Wei. A mixed finite element method with exactly divergence-free velocities for incompressible magnetohydrodynam- ics. Computer Methods in Applied Mechanics and Engineering, 199:2840– 2855, 2010. Chapter 2 is based on paper (a), while Chapter 3 describes the research published in paper (b). Chapter 4 of this thesis presents extensions that are to be submitted. I am the sole author of the results there. I have taken a central role in all aspects of these publications. I was the lead author in paper (a), which is joint with my research supervisor Prof. Dominik Schötzau and Prof. Paul Houston (University of Nottingham). My main contri- bution has been the construction and analysis of the method, and the theoretical derivation of the a priori error estimates, while the numerical tests have been im- plemented and carried out by Prof. Houston. Paper (b) is joint with Prof. Chen Greif (University of British Columbia), Prof. Schötzau, and Dan Li who was a Ph.D. student in the UBC Department of Computer Science supervised jointly by Prof. Greif and Prof. Schötzau. My focus in this paper has been on the numerical analysis of the method and the derivation of iv the a priori error estimates, while my colleague Li has contributed to the numerical experiments. v Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Governing equations of incompressible magnetohydrodynamics . 1 1.2 Background and motivation . . . . . . . . . . . . . . . . . . . . . 6 1.3 Overview of thesis . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.1 A mixed discontinuous Galerkin method for linearized mag- netohydrodynamics . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 A mixed finite element method with exactly divergence- free velocities for nonlinear magnetohydrodynamics . . . 12 1.3.3 Extensions to time-dependent magnetohydrodynamics and Stokes problem . . . . . . . . . . . . . . . . . . . . . . . 16 2 A mixed discontinuous Galerkin method for linearized incompress- ible magnetohydrodynamics . . . . . . . . . . . . . . . . . . . . . . 19 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Discretization of a model problem . . . . . . . . . . . . . . . . . 21 vi 2.2.1 An MHD model problem . . . . . . . . . . . . . . . . . . 21 2.2.2 Meshes and trace operators . . . . . . . . . . . . . . . . . 24 2.2.3 Interior penalty formulation . . . . . . . . . . . . . . . . 25 2.2.4 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 A priori error estimates . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Proof of the error estimates . . . . . . . . . . . . . . . . . . . . . 32 2.4.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.2 Error bounds . . . . . . . . . . . . . . . . . . . . . . . . 36 2.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.5.1 Smooth solutions . . . . . . . . . . . . . . . . . . . . . . 51 2.5.2 Example 3: 2D problem with a singular solution . . . . . 55 2.5.3 Hartmann channel flow . . . . . . . . . . . . . . . . . . . 58 2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3 A mixed finite element method with exactly divergence-free veloci- ties for nonlinear incompressible magnetohydrodynamics . . . . . . 66 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.2 Variational formulation of an MHD problem . . . . . . . . . . . . 70 3.3 Mixed finite element discretization . . . . . . . . . . . . . . . . . 71 3.3.1 Mixed discretization . . . . . . . . . . . . . . . . . . . . 71 3.3.2 Stability properties . . . . . . . . . . . . . . . . . . . . . 75 3.3.3 Existence and uniqueness of discrete solutions . . . . . . 77 3.4 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.4.1 Main results . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.4.2 Continuity . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.4.3 Preliminary error estimates . . . . . . . . . . . . . . . . . 86 3.4.4 Proof of Theorem 3.4.2 . . . . . . . . . . . . . . . . . . . 91 3.4.5 Proof of Theorem 3.4.3 . . . . . . . . . . . . . . . . . . . 93 3.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.5.1 Example 1: two-dimensional problem with a smooth solution 94 3.5.2 Example 2: two-dimensional problem with a singular so- lution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 3.5.3 Hartmann channel flow . . . . . . . . . . . . . . . . . . . 97 vii 3.5.4 Driven cavity flow . . . . . . . . . . . . . . . . . . . . . 104 3.5.5 Example 7: two-dimensional MHD flow over a step . . . . 107 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.1 Time-dependent incompressible magnetohydrodynamics . . . . . 112 4.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.1.2 Weak formulation . . . . . . . . . . . . . . . . . . . . . . 114 4.1.3 Space discretization . . . . . . . . . . . . . . . . . . . . 115 4.1.4 Time discretization . . . . . . . . . . . . . . . . . . . . . 116 4.1.5 Numerical tests . . . . . . . . . . . . . . . . . . . . . . . 117 4.1.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.2 An exactly divergence-free method for the Stokes equations with nonstandard boundary conditions . . . . . . . . . . . . . . . . . . 121 4.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.2.2 Weak formulation and well-posedness . . . . . . . . . . . 125 4.2.3 Finite element approximation . . . . . . . . . . . . . . . 126 4.2.4 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . 128 4.2.5 Numerical examples . . . . . . . . . . . . . . . . . . . . 133 4.2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 134 5 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . 135 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 viii List of Tables Table 1.1 Smooth solution. Convergence of uh and ph in the energy norm. 11 Table 1.2 Smooth solution. Convergence of bh and rh in the energy norm. 11 Table 1.3 Magnetic singularity in L-shaped domain. Convergence of bh and rh in the energy norm. . . . . . . . . . . . . . . . . . . . . 12 Table 1.4 Magnetic singularity in L-shaped domain. Energy norm errors of bh and rh. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Table 1.5 3D channel flow. Convergence of uh and ph in the energy norm. 15 Table 1.6 3D channel flow. Convergence of bh and rh in the energy norm. 16 Table 2.1 Example 2. Convergence of ‖eu‖L2(Ω), ‖eu‖V , and ‖ep‖L2(Ω). . 55 Table 2.2 Example 2. Convergence of ‖eb‖L2(Ω) and ‖eb‖C. . . . . . . . 56 Table 2.3 Example 2. Convergence of ‖er‖L2(Ω) and ‖er‖S. . . . . . . . . 56 Table 2.4 Example 3. Convergence of ‖eu‖L2(Ω), ‖eu‖V , and ‖ep‖L2(Ω). . 58 Table 2.5 Example 3. Convergence of ‖eb‖L2(Ω) and ‖eb‖C. . . . . . . . 59 Table 2.6 Example 3. Convergence of ‖er‖L2(Ω) and ‖er‖S. . . . . . . . . 60 Table 3.1 Discrete inf-sup constants for Vh×Qh. . . . . . . . . . . . . . 77 Table 3.2 Example 1. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). . 95 Table 3.3 Example 1. Convergence of ‖eb‖L2(Ω) and ‖eb‖H(curl;Ω). . . . . 95 Table 3.4 Example 1. Convergence of ‖er‖L2(Ω) and ‖∇er‖L2(Ω). . . . . . 95 Table 3.5 Example 2. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). . 98 Table 3.6 Example 2. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω). 98 Table 3.7 Example 3. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). . 101 Table 3.8 Example 3. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω).101 Table 3.9 Example 4. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). . 104 ix Table 3.10 Example 4. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω).104 Table 4.1 Example 1. Convergence in time of max 1≤n≤N ‖u(tn)− unh‖L2(Ω) and max 1≤n≤N ‖b(tn)−bnh‖L2(Ω). . . . . . . . . . . . . . . . . . . . 118 Table 4.2 Example 1. Convergence of |eu|1,h, ‖eu‖L2(Ω), and ‖ep‖L2(Ω). . 133 Table 4.3 Example 2. Convergence of |eu|1,h, ‖eu‖L2(Ω), and ‖ep‖L2(Ω). . 134 x List of Figures Figure 1.1 MHD flow meter. . . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.2 Electromagnetic pump. . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.3 Contours of (a) b1; (b) b2; (c) nodal approximation of b1; (d) nodal approximation of b2. . . . . . . . . . . . . . . . . . . . 7 Figure 1.4 Elemental degrees of freedom for (a) BDM1(Th); (b) Ned1(Th). 13 Figure 1.5 2D channel flow. Numerical approximations of (a) velocity; (b) normalized magnetic field. . . . . . . . . . . . . . . . . . 15 Figure 1.6 Transient MHD flow. Numerical approximations of velocity at time (a) t = 0.01; (b) t = 1. . . . . . . . . . . . . . . . . . . 17 Figure 1.7 Transient MHD flow. Numerical approximations of normal- ized magnetic field at time (a) t = 0.01; (b) t = 1. . . . . . . 17 Figure 2.1 Example 1. (a) Problem domain; (b) Initial unstructured trian- gular mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Figure 2.2 Example 1. Convergence with h-refinement: (a) ‖eu‖V ; (b) ‖eb‖C; (c) ‖er‖S; (d) ‖eu‖V +‖eb‖C +‖er‖S. . . . . . . . . . 53 Figure 2.3 Example 1. Convergence with h-refinement: (a) ‖eu‖L2(Ω); (b) ‖ep‖L2(Ω); (c) ‖eb‖L2(Ω); (d) ‖er‖L2(Ω). . . . . . . . . . . . . . 54 Figure 2.4 Example 4. Initial unstructured triangular mesh. . . . . . . . . 58 Figure 2.5 Example 4. Convergence with h-refinement: (a) ‖eu‖V ; (b) ‖eb‖C; (c) ‖er‖S; (d) ‖ep‖L2(Ω). . . . . . . . . . . . . . . . . . 61 Figure 2.6 Example 4. DG solution computed on the finest mesh with k = 1: (a) Velocity field; (b) Magnetic field. . . . . . . . . . . 62 xi Figure 2.7 Example 4. DG solution computed on the finest mesh with k = 1. Slices along x = 5, −1 ≤ y ≤ 1 of the solution: (a) First component of the velocity field; (b) First component of the magnetic field. . . . . . . . . . . . . . . . . . . . . . . . 62 Figure 2.8 Example 5. DG solution computed on a uniform tetrahedral mesh with k = 1: (a) Velocity field; (b) Magnetic field. . . . . 63 Figure 2.9 Example 5. DG solution computed on a uniform tetrahedral mesh with k = 1. Slices along x = 5, −1 ≤ y ≤ 1, z = 0 of the solution: (a) First component of the velocity field; (b) First component of the magnetic field. . . . . . . . . . . . . . . . . 64 Figure 3.1 Example 1. Convergence history of the Picard iteration for the grid sequence defined in Tables 3.2–3.4. . . . . . . . . . . . . 96 Figure 3.2 Example 2. Numerical approximations of (a) velocity; (b) pressure contours. . . . . . . . . . . . . . . . . . . . . . . . . 99 Figure 3.3 Example 2. Numerical approximations of (a) magnetic field; (b) contours of the first component of the magnetic field; (c) contours of the second component of the magnetic field. . . . 99 Figure 3.4 Example 3. Slices along x = 5, −1≤ y≤ 1: (a) Velocity com- ponent u(y); (b) Magnetic component b(y). . . . . . . . . . . 101 Figure 3.5 Example 3. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. . . . . . . . . . . . . . . . . . . . . . 102 Figure 3.6 Example 4. Slices along x = 5, −2 ≤ y ≤ 2, and z = 0: (a) Velocity component u(y,0); (b) Magnetic component b(y,0). 105 Figure 3.7 Example 4. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. . . . . . . . . . . . . . . . . . . . . . 106 Figure 3.8 Example 5. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. . . . . . . . . . . . . . . . . . . . . . 107 Figure 3.9 Example 5. Numerical approximations of (a) contours of the first velocity components; (b) streamlines of velocity. . . . . . 108 Figure 3.10 Example 6. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. . . . . . . . . . . . . . . . . . . . . . 109 xii Figure 3.11 Example 7. Numerical approximations of (a) velocity; (b) nor- malized magnetic field; (c) pressure contours. . . . . . . . . . 109 Figure 3.12 Example 7. Velocity flow vectors and streamlines zoomed in behind the step for (a) κ=2.5e4; (b) κ=1e5. . . . . . . . . . . 110 Figure 4.1 Example 2. Numerical approximations of velocity at time (a) t = 0.01; (b) t = 0.1; (c) t = 1. . . . . . . . . . . . . . . . . 120 Figure 4.2 Example 2. Numerical approximations of normalized mag- netic field at time (a) t = 0.01; (b) t = 0.1; (c) t = 1. . . . . . 121 xiii Acknowledgments First and foremost I would like to express my sincerest gratitude to my supervi- sor Prof. Dominik Schötzau for all he has done for me; for his constant support throughout the program, for the abundant ideas and advices he has provided in all academic matters, and for his patience. Without his help, this thesis could not have been possible. I would also like to thank my teachers and committee members Prof. Ian Frigaard, Prof. Chen Greif, Prof. Anthony Peirce, and Prof. Brian Wetton for all I have learned from them. I would like to thank my coauthors Prof. Greif, Prof. Paul Houston, and Dan Li, and my fellow graduate student Liang Zhu for the valuable discussions and ideas. Special thanks to Prof. Greif for his great help and support. Last but not least, I thank the mathematics community and department at UBC for providing an excellent learning and working environment. xiv Chapter 1 Introduction This thesis is concerned with the development, the analysis and the implementation of finite element methods (FEMs) for the numerical approximation of incompress- ible magnetohydrodynamics (MHD) problems. In this introductory chapter, we first present the governing equations of incompressible MHD. We then provide the background and the motivation for our work, followed by a summary of our main results. 1.1 Governing equations of incompressible magnetohydrodynamics Magnetohydrodynamics is the area of physics that describes the behaviour of elec- trically conducting fluids (such as liquid metals, plasmas, salt water, etc.) in the presence of electromagnetic fields; cf. [27, 35, 65, 67, 76]. It combines electro- magnetism and fluid dynamics through two fundamental coupling effects: first, the motion of a conducting material in a transverse magnetic field induces an elec- tric current that modifies the existing electromagnetic field. Second, the electric current and the magnetic field generate a mechanical force on the fluid (Lorentz force). This force accelerates the fluid particles in the direction normal to both the magnetic field and the electric current. We are mainly interested in incompressible magnetohydrodynamics. Here, the electrically conducting viscous fluid is incompressible, and the fluid’s electric re- 1 Figure 1.1: MHD flow meter. sistivity cannot be neglected. Incompressible MHD has a number of technologi- cal and industrial applications such as metallurgical engineering, electromagnetic pumping, stirring of liquid metals, aluminum electrolysis, and measuring flow quantities based on induction; cf. [25, 35]. Figure 1.1 shows the sketch of an MHD flow meter. There, the flow in direction of u transverse to the applied magnetic field b induces the current I, which can then be measured and used to determine the flow speed. Figure 1.2 illustrates the principle of an electromagnetic pump, where the applied magnetic field b and the current density j move the fluid along the pipe in direction of the induced Lorentz force fL. For more applications of incompressible MHD, we refer the reader to, e.g., [67]. Figure 1.2: Electromagnetic pump. In this thesis, we study a standard formulation of the incompressible MHD equations as derived in [3]; see also [33, 35, 48]. That is, we neglect phenom- ena involving high frequency as well as the convection current, and consider an 2 isotropic and homogeneous fluid with constant material properties. We eliminate the electric field and formulate the problem in terms of the velocity, the pressure, and the magnetic field. The corresponding mathematical model forms a system of nonlinear partial differential equations (PDEs) where the incompressible Navier- Stokes equations are coupled with the Maxwell equations. They are derived as follows. The fluid motion is governed by the incompressible Navier-Stokes equations: ∂tu− ηρ ∆u +(u ·∇)u + 1 ρ ∇p = fL + f, (1.1a) ∇ ·u = 0. (1.1b) Here u and p are the velocity and pressure of the fluid, fL is the Lorentz force and f an external body force (both per unit mass). The positive parameters η and ρ are the viscosity and density of the fluid, respectively. The incompressibility constraint (1.1b) corresponds to conservation of mass. Electromagnetic effects are modelled by Maxwell’s equations: ∂tb + ∇× e = 0, Faraday’s law, (1.2a) −∂t(δe)+ ∇× ( 1µ b) = j, Maxwell-Ampère’s law, (1.2b) ∇ · (δe) =ρe, Gauss’ law, (1.2c) ∇ ·b = 0, Gauss’ law for magnetism. (1.2d) The fields appearing in system (1.2) are: the magnetic field b, the electric field e, the electric current density j (per unit area) and the electric charge density ρe (per unit volume). The positive parameters δ and µ are the electric permittivity and magnetic permeability, respectively. The above system has to be supplemented by Ohm’s law. If the densities of positive and negative charges are equal in any sizable region (quasi-neutrality as- 3 sumption as in [3]), the convection current is omitted and Ohm’s law becomes j = θ(e+ u×b). (1.3) Here θ (positive) denotes the electric conductivity of the fluid, and the term u×b represents the electric field induced by the flow. Furthermore we assume that phenomena involving high frequency are not con- sidered. The displacement current ∂t(δe) in (1.2b) can then be neglected, leading to the simplified Ampère’s law ∇× ( 1µ b) = j. (1.4) Combining (1.3)–(1.4) and solving for e gives e = 1 θ j−u×b = 1 θ µ ∇×b−u×b. Now we substitute this expression into (1.2a) to obtain ∂tb + 1 θ µ ∇× (∇×b)−∇× (u×b) = 0, where the term ∇× (u×b) accounts for the first coupling effect. Based on (1.4), the Lorentz force fL in (1.1a) is given by fL = 1 ρ j×b = 1 ρµ (∇×b)×b. Inserting this form of fL into the momentum equation (1.1a) yields ∂tu− ηρ ∆u +(u ·∇)u + 1 ρ ∇p− 1 ρµ (∇×b)×b = f. Here, the term 1ρµ (∇×b)×b incorporates the second coupling effect. After non-dimensionalization and without changing notation, we obtain the 4 following incompressible MHD system: ∂tu−ν∆u +(u ·∇)u + ∇p−κ(∇×b)×b = f, (1.5a) ∇ ·u = 0, (1.5b) ∂tb + νm∇× (∇×b)−∇× (u×b) = 0, (1.5c) ∇ ·b = 0. (1.5d) It has to be supplemented with suitable initial and boundary conditions. The different regimes of (1.5) are characterized by the three non-dimensional numbers ν = Re−1, νm = Rm−1 and κ . The parameter Re is the hydrodynamic Reynolds number. It is defined as Re = ρU 2/L ηU/L2 , and represents the ratio of iner- tial to viscous forces. Here, L and U are the characteristic length and velocity of the problem, respectively. The second parameter Rm = θ µLU is the magnetic Reynolds number. It measures how much the magnetic field is influenced by the flow motion. The non-dimensional coefficient κ is called the coupling number. It is sometimes expressed as a function of the Hartmann number Ha as κ = Ha2 ReRm . The Hartmann number Ha = √ θ B2U ηU/L2 is the root of the ratio of magnetic to viscous forces, where B is the characteristic magnetic induction. It is a measure of the effect of the magnetic field on the flow. For liquid metal flows, the ratio between Rm and Re is typically small. Mer- cury, for instance, has a ratio of order 10−7. In many applications, the mag- nitudes of Rm and κ are of order one, whereas Re can be substantially larger. For example, in aluminum electrolysis, typical values are Rm = 10−1, κ = 1, and Re = 105; cf. [3, 35]. For further discussion of these parameters, we refer the reader also to [71]. From a numerical point of view, the main difficulties relevant to space dis- cretization are already present for the stationary version of (1.5). Thus, we shall 5 mostly consider the following stationary incompressible MHD system: −ν∆u +(u ·∇)u + ∇p−κ(∇×b)×b = f, (1.6a) ∇ ·u = 0, (1.6b) νm∇× (∇×b)−∇× (u×b) = 0, (1.6c) ∇ ·b = 0. (1.6d) In addition, if implicit time-stepping is employed, problems of the form (1.6) (with additional zero order terms) will have to be solved in each time step. 1.2 Background and motivation The numerical approximation of incompressible MHD problems as in (1.5) or (1.6) requires discretizing systems of coupled PDEs. In this thesis, we focus on finite element methods for doing so, since they have become the methods of choice for many large-scale applications. One of the main computational challenges is dealing with the large null-space of the curl-curl operator in (1.5c) or (1.6c). Based on the vector identity −∆b = ∇× (∇×b)−∇(∇ ·b) (1.7) and since ∇ · b = 0, it seems feasible to apply an augmentation technique to re- place the curl-curl operator in (1.5c) or (1.6c) by the vector Laplacian. This would then allow one to use the standard nodal (i.e., H1-conforming) finite elements for the approximation of the magnetic field, which are continuous in all compo- nents over inter-elemental faces. Indeed, various FEM discretizations based on this approach have been proposed for both linear and nonlinear MHD systems, see, e.g., [3, 33, 48]. However, it is well-known that in non-convex polyhedra, the magnetic field may have regularity below H1. A straightforwardly applied nodal FEM discretization, albeit stable, will then converge to a magnetic field that misses certain singular solution components induced by reentrant corners; see [23, 24] and the references therein. This is illustrated in Figures 1.3 for a two-dimensional MHD 6 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 b1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 b2 (c) (d) Figure 1.3: Contours of (a) b1; (b) b2; (c) nodal approximation of b1; (d) nodal approximation of b2. problem on an L-shaped domain. The magnetic field b = (b1,b2), whose contours are shown in Figures 1.3(a) and 1.3(b), represents the strongest magnetic corner singularities at the reentrant corner. It belongs to H2/3. However, it is evident that the nodal approximations obtained on a sufficiently fine mesh and depicted in Figures 1.3(c) and 1.3(d) do not correctly resolve this singularity. A number of remedies have been proposed for electromagnetic problems in isolation, for example the weighted regularization approach in [24] or the ap- proach in [7], whereby the divergence of the electric field is stabilized in H−α with 12 < α < 1. In [50], a weighted regularization nodal finite element method has been introduced and analyzed for the full incompressible MHD system. In the recent work [73], a new variational setting for the formulation of incom- pressible MHD problems has been proposed. It is based on a mixed approach for the discretization of the Maxwell operator. In the stationary case, it amounts to 7 solving a system of PDEs of the following form: −ν ∆u +(u ·∇)u + ∇p−κ (∇×b)×b = f, (1.8a) ∇ ·u = 0, (1.8b) κνm ∇× (∇×b)+ ∇r−κ ∇× (u×b) = g, (1.8c) ∇ ·b = 0. (1.8d) Here, the additional variable r is the Lagrange multiplier related to the divergence constraint of the magnetic field, cf. [52, 66], and g is an additional source term. Indeed, by taking the divergence of (1.8c), we see that ∆r = ∇ ·g. In particular, if g is divergence-free, we obtain r ≡ 0, (1.9) and in this case the sole purpose of the multiplier is to ensure stability. In this setting, the above mentioned difficulties associated with nodal elements are seam- lessly avoided, because the magnetic field can now be sought in a curl-conforming Sobolev space, which is the natural choice especially in the presence of reentrant corners [52, 66]. As a consequence, it is now possible to relax the continuity re- quirements for the approximate magnetic field across inter-elemental faces of the underlying mesh. In particular, one is allowed to enforce only tangential continu- ity (weakly or strongly), thus making it possible to design finite element methods that are able to correctly resolve the strongest magnetic singularities. One option is to use discontinuous Galerkin (DG) methods, which are based on completely discontinuous approximating finite element spaces. Tangential con- tinuity of the approximation to the magnetic field is ensured in a weak sense by introducing in the discrete bilinear forms suitable flux terms over elemental bound- aries. Over the last two decades, DG methods have become an integral part of com- putational fluid mechanics and electromagnetics, see [14, 15, 16, 28, 51] and the 8 references therein. Central features of discontinuous Galerkin methods are their ro- bustness in convection-dominated regimes, their flexibility in the mesh design, their natural way of handling high order and adaptivity, and the fact that the approxima- tion of magnetic field can be based on standard polynomial shape functions. DG methods have been successfully applied to both ideal and viscous compressible MHD problems [64, 77]. In our work [56], the first interior penalty discontinuous Galerkin finite element method has been developed and analyzed for a linearized version of (1.8). In [46], a similar interior penalty technique is applied to enforce the tangential continuity of the magnetic variable across domains with different electromagnetic properties, while nodal elements are employed in the interior. Another option is to employ curl-conforming Nédélec elements for approxi- mating the magnetic field [68]. For these elements, tangential continuity across inter-elemental faces is enforced strongly through appropriate elemental degrees of freedom. In our work [43], a method of this type has been proposed and analyzed for the fully nonlinear system (1.8). The corresponding fluid discretization is based on the divergence-conforming discontinuous Galerkin approach of [20]; it yields exactly divergence-free velocity approximations, and ensures the energy-stability of the resulting method. In [30], a similar discretization employing Nédélec ele- ments for the magnetic field combined with conforming elements for the fluid vari- ables has been presented for a nonlinear MHD problem involving five unknowns (velocity, pressure, magnetic field, electric current and potential). In the recent work [70], various fully discrete schemes based on Nédélec elements for b have been theoretically studied for the time-dependent MHD system (1.5). Other numerical methods for the discretization of the equations of incompress- ible MHD can be found in the literature. We mention here the MHD problem analyzed in [44, 45], which stems from the dynamo effect and deals with domains consisting of both insulating and conducting regions. In [3], the long-term dissi- pativity of time-stepping algorithms for transient incompressible MHD problems has been studied. Convergence results for time-stepping methods involving nodal discretizations in convex domains can be found in [70]. 9 1.3 Overview of thesis In this section we outline the contents and summarize the main results presented in this thesis. Chapter 2 is devoted to the development and analysis of a DG method for a linearized version of (1.8). In Chapter 3, we focus on the fully nonlinear problem (1.8). In Chapter 4, we present two extensions of our discretization tech- niques to (a) the time-dependent problem (1.5), and (b) the Stokes equations with nonstandard boundary conditions. Conclusions and open problems related to our work are presented in Chapter 5. 1.3.1 A mixed discontinuous Galerkin method for linearized magnetohydrodynamics In Chapter 2, we propose and analyze the first interior penalty DG finite element method for a linearized variant of (1.8), whereby all the variables are approximated in discontinuous finite element spaces. More specifically, for k ≥ 1, our method yields approximations uh ∈Pk(Th), ph ∈Pk−1(Th), bh ∈Pk(Th), rh ∈Pk+1(Th), where Pk(Th) is the standard (vector- or scalar-valued) discontinuous finite el- ement space of order k over a tetrahedral triangulation Th of mesh size h of the computational domain. The components of the proposed discretization are bor- rowed from the DG methods available for incompressible fluid flow problems [17, 18, 19, 49, 74], and from the discontinuous element pair proposed for the mixed Maxwell operator in [54, 55], see also [53], combined with a discontinuous discretization of the coupling terms. Our main result is an a priori error estimate for the proposed method under the following smoothness assumption on the analytical solution: u ∈Hσ+1, p ∈ Hσ , (1.10a) b ∈Hτ , ∇×b ∈ Hτ , r ∈ Hτ+1, (1.10b) 10 for regularity parameters σ ,τ > 12 . These regularity assumptions are minimal in the sense that they are satisfied by the strongest hydrostatic and magnetic singular- ities in possibly non-convex polyhedral domains; see [2, 26] and [23, 24]. More precisely, we prove that |||(u−uh,b−bh, p− ph,r− rh)|||E = O(hmin{σ ,τ ,k}), (1.11) where ||| · |||E is a suitably defined energy norm and h the mesh size. In particu- lar, this shows convergence of order O(hmin{σ ,τ}) for non-smooth solutions, and order O(hk) for smooth ones. This rate is optimal in the approximation of the velocity, the pressure and the magnetic field. k DOFs uh/ph |||u−uh|||E rate |||p− ph|||E rate 4,608/384 5.823e-1 – 4.564e-1 – 1 36,864/3,072 2.896e-1 1.01 2.606e-1 0.81 294,912/24,576 1.442e-1 1.01 1.400e-1 0.90 1,440/192 1.262e-1 – 3.424e-1 – 2 11,520/1,536 3.162e-2 2.00 8.488e-2 2.01 92,160/12,288 7.822e-3 2.02 2.127e-2 2.00 Table 1.1: Smooth solution. Convergence of uh and ph in the energy norm. k DOFs bh/rh |||b−bh|||E rate |||r− rh|||E rate 4,608/3,840 3.445e-1 – 5.098e-1 – 1 36,864/30,720 1.668e-1 1.05 1.363e-1 1.90 294,912/245,760 8.184e-2 1.03 3.405e-2 2.00 1,440/960 4.920e-2 – 2.559e-1 – 2 11,520/7,680 1.146e-2 2.10 3.430e-2 2.90 92,160/61,440 2.767e-3 2.05 4.210e-3 3.03 Table 1.2: Smooth solution. Convergence of bh and rh in the energy norm. In Tables 1.1 and 1.2, we show representative convergence results for a three- dimensional problem with a smooth solution, employing k = 1 and k = 2. The experimental orders of convergence for the velocity, the pressure, and the magnetic 11 field are of order O(hk), confirming our theoretical results in (1.11). They are of order O(hk+1) for the error in the multiplier r. While optimal and expected, this rate is not reflected in our error estimates. We also note that in Section 2.5.1.2 the L2-norm errors in the velocity and the magnetic field are observed to be of the optimal order O(hk+1). For the L-shaped domain problem with the singular magnetic field shown in Figure 1.3, we obtain the numbers listed in Table 1.3. They show that our method correctly captures the singular behaviour of b at the reentrant corner, at convergence rates that are actually better than expected. On the other hand, the rates for r are of order O(h2/3), which is in agreement with the regularity of the magnetic field. A version of this chapter has been published in [56]. k DOFs bh/rh |||b−bh|||E rate |||r− rh|||E rate 2,304/2,304 1.112e-1 – 1.265 – 1 9,216/9,216 5.280e-2 1.07 8.384e-1 0.59 36,864/36,864 2.657e-2 0.99 5.387e-1 0.64 4,608/3,840 8.065e-2 1.15 1.412 – 2 18,432/15,360 3.654e-2 1.14 9.193e-1 0.62 73,728/61,440 1.766e-2 1.05 5.868e-1 0.65 7,680/5,760 6.363e-2 – 1.559 – 3 30,720/23,040 2.812e-2 1.18 1.008 0.63 122,880/92,160 1.320e-2 1.09 6.415e-1 0.65 Table 1.3: Magnetic singularity in L-shaped domain. Convergence of bh and rh in the energy norm. 1.3.2 A mixed finite element method with exactly divergence-free velocities for nonlinear magnetohydrodynamics While it is in principle possible to extend the fully discontinuous approach of Chap- ter 2 to the nonlinear setting, however in Chapter 3, for the reasons presented below, we have chosen to suitably modify the DG discretization by proceeding as follows. First, we replace the discontinuous Pk elements for the velocity by divergence- conforming Brezzi-Douglas-Marini (BDM) elements of degree k. These elements are continuous in normal direction over inter-elemental faces, and hence the DG 12 approach is now only employed to enforce the tangential continuity of the veloc- ity field. This choice of elements has been introduced and analyzed in [19, 20] in the context of the incompressible Navier-Stokes equations. It has the attractive property that it yields exactly divergence-free velocity approximations, which nat- urally guarantees the energy-stability of the discretization. In [20], a detailed error analysis of exactly divergence-free discretizations can be found for a variety of DG methods; they have been shown to be inf-sup stable and optimally convergent in natural norms. We also refer the reader to [21, 59] for further aspects. Second, we replace the completely discontinuous Pk-Pk+1 elements for the discretization of the magnetic variables b and r by the conforming first family Nédélec pair of order k. That is, we seek the approximate magnetic field in a curl- conforming finite element space of piecewise polynomials of degree k, and the discrete Lagrange multiplier in the H1-conforming nodal space of order k. This choice of elements has the advantage that it reduces the total number of coupled degrees of freedom. Based on the discrete Helmholtz decomposition [66], it also yields rh ≡ 0 (1.12) for a divergence-free right-hand g, thus mimicking the continuous scenario in (1.9). (a) (b) Figure 1.4: Elemental degrees of freedom for (a) BDM1(Th); (b) Ned1(Th). In summary, we now determine approximations to (u,b, p,r) in the finite ele- 13 ment spaces uh ∈ BDMk(Th), ph ∈Pk−1(Th), (1.13a) bh ∈Nedk(Th), rh ∈Pck (Th), (1.13b) where Pck (Th) denotes the space of continuous piecewise polynomials of degree at most k. The lowest-order BDM and Nédélec elements (k = 1) are illustrated in Figures 1.4(a) and 1.4(b), respectively. Our main results are then proofs of the following theoretical properties of the resulting finite element discretization. First, we show the existence and uniqueness of discrete solutions under a standard small data assumption. Second, we show convergence in general (possibly non-convex) polyhedral domains under the min- imal regularity assumption (1.10). That is, as the mesh size h tends to zero, we have |||(u−uh,b−bh, p− ph,r− rh)|||E → 0. Finally, we carry out an error analysis and prove the following error estimate |||(u−uh,b−bh)|||E = O(h min{σ ,τ ,k}−ε ), in 2D, O(hmin{σ ,τ ,k}− 12 ), in 3D, for ε > 0 arbitrarily small. This estimate is nearly optimal for two-dimensional problems, but falls short by half a power of h in three dimensions. This loss in optimality stems from the use of Sobolev embeddings and inverse estimates to establish the continuity of the variational forms associated with the nonlinear cou- pling terms. However, all our numerical experiments indicate optimal convergence rates, in both two and three dimensions. Let us again present some representative numerical results for the proposed method. In Table 1.4 we show the numbers obtained for the same singular prob- lem on the L-shaped domain as in Figure 1.3. The approximate magnetic field is convergent of order O(h2/3), as predicted by our theory, whereas the L2-norm of the approximate multiplier rh is now zero up to machine accuracy, as expected from (1.12). As the mesh is refined, the actual values of the L2-norm of rh slightly 14 increase, which is likely due to the increased condition numbers of the correspond- ing linear systems. DOFs bh/rh |||b−bh|||E rate |||rh|||L2 2,368/833 7.473e-2 – 4.260e-11 9,344/3,201 4.754e-2 0.65 1.406e-10 37,120/12,545 3.013e-2 0.66 3.018e-10 Table 1.4: Magnetic singularity in L-shaped domain. Energy norm errors of bh and rh. (a) (b) Figure 1.5: 2D channel flow. Numerical approximations of (a) velocity; (b) normalized magnetic field. DOFs uh/ph |||u−uh|||E rate |||p− ph|||E rate 2,592/384 0.9561 – 8.194 – 19,584/3,072 0.4903 0.96 2.837 1.53 152,064/24,576 0.2484 0.98 1.091 1.38 Table 1.5: 3D channel flow. Convergence of uh and ph in the energy norm. In Figure 1.5, we show the approximate velocity and the (normalized) mag- netic field of a unidirectional channel flow driven by a pressure gradient under the transverse magnetic field (0,1). The solution of this problem is smooth. The 15 DOFs bh/rh |||b−bh|||E rate |||rh|||L2 604/125 2.579e-5 – 1.013e-10 4,184/729 1.464e-5 0.82 4.098e-10 31,024/4,913 7.543e-6 0.96 1.795e-9 Table 1.6: 3D channel flow. Convergence of bh and rh in the energy norm. velocity field only has a component in x-direction, with a profile that can be ex- pressed in terms of a hyperbolic cosine (rather than the parabolic profile obtained for the Navier-Stokes problem in isolation). The magnetic field is of the form b = ((b(y),1) where b(y) measures the deviation from (0,1). The approximations shown in Figure 1.5 are in excellent agreement with the exact solutions. A three- dimensional version of this problem has also been tested. In Tables 1.5 and 1.6, we list the computational results. We see convergence rates of order one in u and b, and of order slightly better than one in p. We observe again that rh is zero up to machine accuracy; the slight increase of the actual errors as the mesh is refined is again due to the increased condition numbers. A version of this chapter has been published in [43]. 1.3.3 Extensions to time-dependent magnetohydrodynamics and Stokes problem In Chapter 4, we present two extensions of the techniques introduced in Chapters 2 and 3. 1.3.3.1 Time-dependent magnetohydrodynamics computations In [70], theoretical aspects of fully discrete schemes for the time-dependent MHD problem (1.5) have been established. Our goal here is to implement one of these methods and show a number of computational tests. We employ the elements proposed in (1.13), and discretize in time using the implicit Euler scheme. As a result, a linearized (but still coupled) version of prob- lem (1.8) has to be solved in each time step. In Figures 1.6 and 1.7, we show the evolution of the approximate velocity and the (normalized) magnetic field of a channel flow driven by a moving top wall 16 (a) (b) Figure 1.6: Transient MHD flow. Numerical approximations of velocity at time (a) t = 0.01; (b) t = 1. under a transverse magnetic field. The analytical solutions of this problem can be expressed by Fourier series. The approximations shown in Figures 1.6 and 1.7 are in excellent agreement with the exact solutions. (a) (b) Figure 1.7: Transient MHD flow. Numerical approximations of normalized magnetic field at time (a) t = 0.01; (b) t = 1. 1.3.3.2 An optimal error estimate of a curl-curl formulation of the Stokes equations with nonstandard boundary conditions In the spirit of the method proposed in [20], we investigate an exactly divergence- free DG method for the numerical approximation of the Stokes equations with the nonstandard boundary conditions u ·n = 0 and (∇×u)×n = 0. These boundary conditions are of practical interest in computational fluid dynam- ics; cf. [8, 37]; they correspond to normal velocity (no-penetration) and tangential vorticity conditions. They lend themselves naturally to numerical methods that are based on reducing the vector Laplacian to the curl-curl operator, see also (1.7), and 17 rewriting the Stokes equations in the curl-curl formulation ν∇× (∇×u)+ ∇p = f, ∇ ·u = 0. Again, we use the H(div)-conforming BDM finite element space of degree k for the velocity approximation [11], along with a discontinuous pressure space of degree k−1. The tangential continuity of the approximate velocity field is enforced through an interior penalty approach. We establish a crucial norm-equivalence property and show convergence of order O(hk) in the mesh size for the broken H1- norm of the velocity error, despite the fact that the method does not give any explicit control on gradients. We confirm the theoretical results on a series of numerical tests. 18 Chapter 2 A mixed discontinuous Galerkin method for linearized incompressible magnetohydrodynamics In this chapter, we introduce and analyze the first interior penalty discontinuous Galerkin method for the numerical discretization of a stationary incompressible magnetohydrodynamics model problem. The fluid unknowns are discretized with inf-sup stable discontinuous Pdk -Pk−1 elements whereas the magnetic part of the equations is approximated by discontinuous Pdk -Pk+1 elements. We carry out a complete a priori error analysis of the method and prove that the energy norm error is optimally convergent in the mesh size. These results are verified in a series of numerical experiments. 2.1 Introduction The numerical simulation of incompressible MHD problems requires discretizing a system of partial differential equations that couples the incompressible Navier- Stokes equations with Maxwell’s equations. Various finite element methods (FEM) can be found in the literature where the magnetic field is approximated by standard 19 H1-conforming finite elements, see, e.g., [3, 33, 45, 48] and the references therein. However, in non-convex polyhedra of engineering interest, the magnetic field may have regularity below H1 and a nodal FEM discretization, albeit stable, can con- verge to a magnetic field that misses certain singular solution components induced by reentrant vertices or edges; see [24]. In the recent work [73], this drawback of nodal elements was overcome by the use of Nédélec elements for the approxima- tion of the magnetic field. Thereby, a new variational setting for the formulation of incompressible MHD problems was proposed. It is based on a mixed approach for the discretization of the Maxwell operator and introduces a Lagrange multiplier related to the divergence constraint of the magnetic field, cf. [66]. Over the last two decades, discontinuous Galerkin methods have become an in- tegral part of computational fluid mechanics and computational electromagnetics, see [14, 16, 51] and the references therein. DG methods are extremely versatile and flexible; they can deal robustly with partial differential equations of almost any kind, as well as with equations whose type changes within the computational domain. Their intrinsic stability properties make them naturally suited for prob- lems where convection is dominant. Moreover, discontinuous Galerkin methods can easily handle irregularly refined meshes and variable approximation degrees (hp-adaptivity). The DG approximations of magnetic or electric fields can be based on standard polynomial shape functions, in contrast to curl-conforming or divergence-conforming elements commonly used in computational electromagnet- ics. DG methods have already been successfully applied to both ideal and viscous compressible MHD problems [64, 77]. In this chapter, we propose and analyze an interior penalty DG method for a linearized incompressible MHD model problem based on the mixed formulation introduced in [73]. Our method combines the DG discretizations that have been developed recently for incompressible flow problems and Maxwell’s equations. More specifically, the fluid unknowns are approximated using mixed discontinu- ous Pdk -Pk−1 elements [18, 19, 74] while the magnetic variables are discretized with the Pdk -Pk+1 element pair proposed and analyzed in [54, 55]. We carry out a complete a priori error analysis for the proposed DG method, and show that the energy error in all variables is convergent of order O(hk) in the mesh size h. Our results further show that the proposed DG method is able to correctly resolve the 20 strongest magnetic singularities in non-convex polyhedral domains, in contrast to H1-conforming elements. The rest of the chapter is organized as follows. In Section 2.2, we introduce an interior penalty DG method for a linearized incompressible MHD model prob- lem. In Section 2.3, we state and discuss a priori error estimates for the method. Section 2.4 is devoted to the detailed proof of these estimates. In Section 2.5, we present a series of numerical experiments validating our theoretical results. Finally, we present some concluding remarks in Section 2.6. 2.2 Discretization of a model problem 2.2.1 An MHD model problem We consider the following linear and stationary MHD system based on the mixed formulation proposed in [73]: find the velocity field u, the pressure p, the magnetic field b, and the scalar potential r such that −ν ∆u +(w ·∇)u + γu + ∇p−κ (∇×b)×d = f in Ω, (2.1a) ∇ ·u = 0 in Ω, (2.1b) κ νm ∇× (∇×b)+ ∇r−κ ∇× (u×d) = g in Ω, (2.1c) ∇ ·b = 0 in Ω. (2.1d) Here, Ω is a bounded simply-connected Lipschitz polytope in Rd (d = 2 or 3). In the two-dimensional case, the curl operator ∇× applied to a vector b = (b1,b2) is determined by ∇× (b1,b2,0), while the curl of a scalar function r is given by ∇ × (0, 0, r). That is, ∇×b = ∂b2∂x − ∂b1∂y , and ∇× r = ( ∂ r∂y ,− ∂ r∂x). The cross prod- uct × is defined similarly. The function w ∈W 1,∞(Ω)d is a prescribed convective field, and d ∈ L∞(Ω)d a given magnetic field. Typically these fields come from a linearization process. The right-hand sides f and g are vector-valued source terms in L2(Ω)d . The scalar function γ belongs to L∞(Ω). We further assume that there 21 is a positive constant γ⋆ such that γ0(x) := γ(x)− 12∇ ·w(x)≥ γ⋆ > 0, x ∈Ω. (2.2) Remark 2.2.1 The positivity of γ⋆ in (2.2) is a purely technical (but standard) assumption that facilitates dealing with the convection term in the error analysis. However, in the absence of a reaction term (γ ≡ 0), the parameter γ⋆ must be allowed to be zero. While the proposed DG method is stable and well-defined in this case as well, the error analysis becomes more involved and requires additional (duality) arguments, see [18] for the Oseen operator. Without loss of generality, we may assume that the length scale of Ω, and the L∞-norms of w and d are one. In many engineering applications (such as aluminum electrolysis), the magnitudes of Rm and κ are of order one, whereas Re can be substantially larger; cf. [3]. We will focus on this case and not make explicit our error estimates with respect to νm and κ . We suppose that the boundary Γ of Ω is connected, and can be partitioned into two disjoint parts. That is, we have Γ = ΓD∪ΓN with ΓD∩ΓN = /0. Throughout, we assume that ΓD satisfies ∫ ΓD ds > 0. Denoting by n the unit out- ward normal on the boundary, we then supplement the MHD system (2.1) with the following boundary conditions: u = uD on ΓD, (2.3a) (pI−ν∇u)n = pNn on ΓN , (2.3b) n×b = n×bD on Γ, (2.3c) r = 0 on Γ. (2.3d) Here, I is the identity matrix in Rd×d. We assume that pN ∈ L2(ΓN) and that uD and bD are restrictions to the boundary of sufficiently smooth divergence-free func- tions in Ω. Finally, we notice that, if ΓN = /0, the datum uD must satisfy the com- 22 patibility condition ∫ Γ uD ·nds = 0. Remark 2.2.2 The magnetic boundary conditions (2.3) are geared towards Hart- mann flow problems with insulating wall conditions in the setting of [35, Sec- tion 3.7.1]. The method can be readily extended to other types of electromagnetic boundary conditions. For example, with only minimal changes, it is possible to specify both b ·n and (∇×b)×n on Γ, corresponding to perfectly conducting wall conditions; cf. [33, 35, 73]. Let Γ− = {x ∈ Γ : w(x) ·n(x) < 0} be the inflow boundary of Γ. We adopt the (physically reasonable) hypothesis that w(x) ·n(x)≥ 0 for all x ∈ ΓN . Obviously, we then have Γ− ⊆ ΓD. Remark 2.2.3 As in [73], the scalar potential r is the Lagrange multiplier associ- ated with the magnetic divergence constraint. By taking the divergence of (2.1c), we see that ∆r = ∇ ·g in Ω, r = 0 on Γ. In particular, we have r = 0 provided that the function g is divergence-free. In this case, the MHD problem (2.1) is the same as the one studied in [33] or the linearized version of the one considered in [48]. To define the weak formulation of the MHD system, we introduce the Sobolev spaces V = { v ∈ H1(Ω)d : v = 0 on ΓD } , C = H0(curl;Ω) = { c ∈ L2(Ω)d : ∇× c ∈ L2(Ω)d , n× c = 0 on Γ } , S = H10 (Ω) = {s ∈ H1(Ω) : s = 0 on Γ}, and Q = L2(Ω). In the case where ΓN = /0, we also need to enforce the mean val- ues of functions in Q to be zero. We denote by (·, ·)Ω the inner product in L2(Ω) or L2(Ω)d , and by 〈·, ·〉Γ′ the inner product in L2(Γ′) or L2(Γ′)d for Γ′ ⊆ Γ. The weak formulation of the incompressible MHD system (2.1) then consists in de- termining u ∈ H1(Ω)d , b ∈ H(curl;Ω), p ∈ Q and r ∈ S, with u = uD on ΓD 23 and n×b = n×bD on Γ, such that A(u,v)+ O(u,v)+C(v,b)+ B(v, p) = (f,v)Ω−〈pNn,v〉ΓN , B(u,q) = 0, M(b,c)−C(u,c)+ D(c,r) = (g,c)Ω, D(b,s) = 0 for all (v,c,q,s) ∈ V×C×Q×S. Here, the bilinear forms are given by A(u,v) = ∫ Ω ν ∇u : ∇vdx, O(u,v) = ∫ Ω ((w ·∇)u + γ u) ·vdx, M(b,c) = ∫ Ω κ νm(∇×b) · (∇× c)dx, C(v,b) = ∫ Ω κ (v×d) · (∇×b)dx, B(u,q) =− ∫ Ω (∇ ·u)qdx, D(b,s) = ∫ Ω b ·∇sdx. Under the above assumptions, the well-posedness of this problem follows from standard stability properties and the theory of mixed finite elements; see also [73] and the references therein. 2.2.2 Meshes and trace operators We consider a family of regular and shape-regular triangulations Th that partition the domain Ω into simplices {K} (i.e., triangles for d = 2 and tetrahedra for d = 3). The index h is indicative of the mesh size h which is defined as h = maxK∈Th hK , where hK is the diameter of K. We denote by Fh the set of all edges (d = 2) or faces (d = 3) of Th. In the following, we generically refer to elements in Fh as faces. We also denote by F Ih the set of all interior faces of Th, and by FBh the set of all boundary faces. We always assume that FBh can be divided into two disjoint sets FDh and FNh of Dirichlet and Neumann faces, respectively. That is, we assume that FBh = FDh ∪FNh , where ΓD = ∪F∈FDh F and ΓN = ∪F∈FNh F . As usual, hF denotes the diameter of the face F . Finally, we write nK to denote the unit outward normal on the boundary ∂K of an element K. Next, we introduce the average and jump operators. To do so, let F = ∂K∩∂K′ 24 be an interior face shared by K and K′ and let x ∈ F . Let ϕ be a generic piecewise smooth function (scalar-, vector- or tensor-valued) and denote by ϕ and ϕ ′ the traces of ϕ on F taken from within the interior of K and K′, respectively. Then, we define the mean value of ϕ at x ∈ F as {{ϕ}}= (ϕ + ϕ ′)/2. Furthermore, let ψ be a piecewise smooth function and φ a piecewise smooth vector-valued field. Analogously, we define the following jumps at x ∈ F : [[ψ ]] = ψ nK + ψ ′nK′ , [[φ ]] = φ ⊗nK + φ ′⊗nK′ , [[φ ]]T = nK ×φ + nK′×φ ′, [[φ ]]N = φ ·nK + φ ′ ·nK′ , where φ ⊗n = (φin j)1≤i, j≤d . On a boundary face F = ∂K∩Γ, we set accordingly {{ϕ}}= ϕ , [[ψ ]] = ψ n, [[φ ]] = u⊗n, [[φ ]]T = n×φ , [[φ ]]N = φ ·n. 2.2.3 Interior penalty formulation For k ≥ 1, let Pk(K) denote the space of polynomials of total degree at most k on K. We then define Pk(Th) = { p ∈ L2(Ω) : p|K ∈Pk(K), K ∈Th }. The corresponding vector-valued function space is denoted by Pk(Th)d with d = 2 or 3. Now we introduce the finite element spaces Vh = Pk(Th)d , Ch = Pk(Th)d , Qh = Pk−1(Th), Sh = Pk+1(Th), where we also impose zero mean value of the functions in Qh in the case ΓN = /0. We consider the following discontinuous Galerkin method: find (uh,bh, ph,rh) 25 in Vh×Ch×Qh×Sh such that Ah(uh,v)+ Oh(uh,v)+Ch(v,bh)+ Bh(v, ph) = Fh(v), (2.4a) Bh(uh,q) = 〈uD ·n,q〉ΓD , (2.4b) Mh(bh,c)−Ch(uh,c)+ Dh(c,rh) = Gh(c), (2.4c) Dh(bh,s)− Jh(rh,s) = 0 (2.4d) for all (v,c,q,s) ∈ Vh×Ch×Qh× Sh. Here, the forms Ah, Oh and Bh are related to the discretization of the Oseen operator. We take the ones proposed and studied in [18, 19, 49, 74]. The forms Mh, Dh and Jh are related to the discretization of the Maxwell operator. We choose the ones corresponding to the non-stabilized Pdk −Pk+1 interior penalty methods proposed and analyzed in [54, 55]. Finally, the form Ch couples the Maxwell equations to the Oseen problem. These forms are defined next. First, the form Ah is chosen as the standard interior penalty form Ah(u,v) = ∑ K∈Th ∫ K ν ∇u : ∇vdx− ∑ F∈F Ih∪FDh ∫ F {{ν∇u}} : [[v]]ds − ∑ F∈F Ih∪FDh ∫ F {{ν∇v}} : [[u]]ds+ ∑ F∈F Ih∪FDh νa0 hF ∫ F [[u]] : [[v]]ds. The parameter a0 > 0 is a sufficiently large stabilization parameter; see Proposi- tion 2.2.4 below. For the convective form, we take the usual upwind form defined by Oh(u,v) = ∑ K∈Th ∫ K ((w ·∇)u + γ u) ·vdx + ∑ K∈Th ∫ ∂K−\Γ− w ·nK(ue−u) ·vds− ∫ Γ− w ·nu ·vds. Here, ue is the value of the trace of u taken from the exterior of K and ∂K− = {x ∈ ∂K : w(x) ·nK(x) < 0} is the inflow boundary of K. The form Bh related to the 26 divergence constraint on u is defined by Bh(u,q) =− ∑ K∈Th ∫ K (∇ ·u)qdx+ ∑ F∈F Ih∪FDh ∫ F {{q}}[[u]]N ds. Next, we define the forms for the discretization of the Maxwell operator. The form Mh for the curl-curl operator is given by Mh(b,c) = ∑ K∈Th ∫ K κ νm(∇×b) · (∇× c)dx− ∑ F∈Fh ∫ F {{κνm∇×b}} · [[c]]T ds − ∑ F∈Fh ∫ F {{κνm∇× c}} · [[b]]T ds+ ∑ F∈Fh κνmm0 hF ∫ F [[b]]T · [[c]]T ds. As for the diffusion form, to ensure stability, the stabilization parameter m0 > 0 must be chosen large enough, see Proposition 2.2.4 below. The form Dh for the divergence constraint on b is given by Dh(b,s) = ∑ K∈Th ∫ K b ·∇sdx− ∑ F∈Fh ∫ F {{b}} · [[s]]ds. The form Jh is a stabilization term that ensures the H1-conformity of the multi- plier rh. It is given by Jh(r,s) = ∑ F∈Fh s0 κνmhF ∫ F [[r]] · [[s]]ds, with s0 > 0 denoting a positive stabilization parameter. The dependence on νm and κ is chosen so as to suitably balance the multiplier terms in our error analysis. Finally, for the coupling form Ch in our DG formulation, we take a discontinu- ous Galerkin version of the bilinear form C, namely: Ch(v,b) = ∑ K∈Th κ ∫ K (v×d) · (∇×b)dx− ∑ F∈F Ih∪FNh κ ∫ F {{v×d}} · [[b]]T ds. 27 With these forms, the source terms Fh(v) and Gh(c) must be chosen as Fh(v) = ∫ Ω f ·vdx− ∑ F∈FDh ∫ F ν∇v : (uD⊗n)ds + ∑ F∈FDh νa0 hF ∫ F uD ·vds− ∑ F∈FNh ∫ F κ(v×d) · (n×bD)ds − ∫ Γ− w ·nuD ·vds− ∑ F∈FNh ∫ F pNn ·vds, and Gh(c) = ∫ Ω g · cdx− ∑ F∈FBh ∫ F κνm(∇× c) · (n×bD)ds + ∑ F∈FBh κνmm0 hF ∫ F (n×bD) · (n× c)ds − ∑ F∈FDh ∫ F κ(uD×d) · (n× c)ds, respectively. 2.2.4 Stability The stability properties of the above DG forms have been well established in the recent literature on DG methods. To review them, we introduce the norms ‖u‖21,h = ∑ K∈Th ‖∇u‖2L2(K) + ∑ F∈F Ih∪FDh h−1F ‖[[u]]‖2L2(F), ‖u‖2V = ν‖u‖21,h +‖γ 1 2 0 u‖2L2(Ω) + 1 2 ∑F∈Fh ‖|w ·n| 1 2 [[u]]‖2L2(F). 28 In the last term, n denotes any of the outward normals on F . For the magnetic field, we define |b|2C = κνm ∑ K∈Th ‖∇×b‖2L2(K) + κνm ∑ F∈Fh h−1F ‖[[b]]T‖2L2(F), ‖b‖2C = κνm‖b‖2L2(Ω) + |b|2C. Finally, on the magnetic multiplier space, we introduce ‖r‖2S = κ−1ν−1m ∑ K∈Th ‖∇r‖2L2(K) + κ−1ν−1m ∑ F∈Fh h−1F ‖[[r]]‖2L2(F). First, we recall the following coercivity properties of Ah, Oh and Mh, see [5, 18, 54] and the references therein. Proposition 2.2.4 Under the assumption (2.2), there is a threshold value a0 > 0, independent of the mesh size, ν , νm and κ , such that for every a0 ≥ a0 there is a constant C > 0 independent of the mesh size, ν , νm and κ such that Ah(u,u)+ Oh(u,u)≥C‖u‖2V , u ∈ Vh. Moreover, there is a threshold value m0 > 0, independent of the mesh size, ν , νm, and κ , such that for m0 ≥ m0 there is a constant C > 0 independent of the mesh size, ν , νm and κ such that Mh(b,b)≥C|b|2C, b ∈ Ch. Next, we recall that the velocity/pressure pair Vh×Qh is inf-sup stable; cf. [11, Remark II.2.10] and [49, Proposition 10]: Proposition 2.2.5 There is a stability constant C > 0 independent of the mesh size, ν , νm and κ such that inf p∈Qh\{0} sup v∈Vh\{0} Bh(v, p) ‖v‖1,h‖p‖L2(Ω) ≥C > 0. 29 There is no inf-sup condition available for the pair Ch × Sh. However, the underlying conforming spaces are stable, see [66]. To discuss this, we introduce the conforming spaces Cch = Ch∩C and Sch = Sh∩S. The space Cch is the Nédélec finite element space of the second family of order k [66, 68], with zero tangential trace on Γ. The space Sch is the space of continuous piecewise polynomials of degree at most k + 1, with zero trace on Γ. Thus, we may decompose Ch and Sh into Ch = Cch⊕C⊥h , Sh = Sch⊕S⊥h , (2.5) respectively. Obviously, the norms of the jumps |b|2C⊥ = κνm ∑ F∈Fh h−1F ‖[[b]]T‖2L2(F), |r|2S⊥ = κ−1ν−1m ∑ F∈Fh h−1F ‖[[r]]‖2L2(F), define norms on C⊥h and S⊥h , respectively. The following norm-equivalence results from [54, Theorem 4.1] are essential to our error analysis. Proposition 2.2.6 There is a constant C > 0, independent of the mesh size, ν , νm and κ , such that C‖b‖C ≤ |b|C⊥ ≤ ‖b‖C, C‖r‖S ≤ |r|S⊥ ≤ ‖r‖S, for any b ∈ C⊥h and r ∈ S⊥h . For the conforming pair Cch×Sch, the following properties of the forms M and D hold true, see [54, Lemma 5.3] for a proof. Proposition 2.2.7 There exists a constant C > 0 independent of the mesh size, ν , νm and κ such that M(b,b)≥C‖b‖2C, for any b in Xch, where Xch = {b ∈ Cch : D(b,s) = 0 ∀s ∈ Sch }. (2.6) Furthermore, there exists a second constant C > 0 independent of the mesh size, ν , 30 νm and κ such that inf r∈Sch\{0} sup c∈Cch\{0} D(c,r) ‖c‖C‖r‖S ≥C > 0. Employing the above properties and applying arguments as in [54, Proposition 3.3], existence and uniqueness of discrete solutions can be readily shown provided that a0 ≥ a0, m0 ≥ m0 and s0 > 0. 2.3 A priori error estimates In this section, we present and discuss the main result of this chapter: a priori error estimates for the proposed DG method. We shall suppose that the solution (u,b, p,r) of the MHD problem satisfies the regularity properties (u, p) ∈ Hσ+1(Ω)d ×Hσ(Ω), for σ > 12 , (2.7a) (b,∇×b,r) ∈Hτ(Ω)d ×Hτ(Ω)d ×Hτ+1(Ω), for τ > 12 . (2.7b) These regularity assumptions are realistic. This can be seen from the regularity properties of the Maxwell operator and the linearized Navier-Stokes operator in polyhedral domains, respectively, see [2, 26]. In particular, the strongest magnetic singularities satisfy assumption (2.7b). The following theorem represents the main result of this chapter. Theorem 2.3.1 Let the solution (u,b, p,r) of the MHD problem satisfy the regu- larity assumptions stated in (2.7). Further, let (uh,bh, ph,rh) denote the DG ap- proximation defined in (2.4). Assuming (2.2) holds, the errors can be bounded by ‖u−uh‖V +‖b−bh‖C +‖p− ph‖L2(Ω) +‖r− rh‖S ≤Chmin{σ ,k} ( ‖u‖Hσ+1(Ω) + ν− 1 2 ‖p‖Hσ (Ω) ) +Chmin{τ ,k} ( ‖b‖Hτ (Ω) +‖∇×b‖Hτ(Ω) +‖r‖Hτ+1(Ω) ) , 31 where C is a positive constant, independent of the mesh size and ν . Remark 2.3.2 For smooth solutions, the estimate in Theorem 2.3.1 ensures con- vergence rates of order O(hk) in the mesh size h. This rate is optimal in the ap- proximation of the velocity, the pressure and the magnetic field in the respective norms, but suboptimal by one order in the approximation of the multiplier r with respect to the norm ‖ · ‖S. This is due to the fact that we are using polynomials of degree k + 1 to approximate r. The same suboptimal result is observed for the conforming Nédélec family of the second type [66, 68]. On the other hand, the use of polynomials of degree k + 1 for the magnetic multiplier leads to optimal convergence rates in the L2-error in the magnetic field b, in contrast to the use of polynomials of degree k, cf. the discussion in [53, 54]. Remark 2.3.3 Our error estimates also hold in the case where the conforming Nédélec pair Cch×Sch is used for the approximation of b and r. While these spaces have less degrees of freedom than their discontinuous counterparts Ch and Sh, the use of discontinuous approximations for the magnetic field has several advantages. For example, DG approximations can be based on standard polynomial shape func- tions which greatly facilitates the implementation of higher-order elements and magnetic boundary conditions. Moreover, they are naturally suited to deal with irregularly refined meshes and variable approximation degrees (hp-adaptivity). 2.4 Proof of the error estimates 2.4.1 Preliminaries For the purpose of our analysis, we set V(h) = V+Vh, C(h) = C+Ch and S(h) = S + Sh. Using the lifting operators constructed in [5, 74] and [53, 54], it is then possible to extend the discrete bilinear forms Ah, Bh, Mh, Dh to bilinear forms Ãh : V(h)×V(h)→R, B̃h : V(h)×Q→R, M̃h : C(h)×C(h)→R and D̃h : C(h)× 32 S(h)→ R, respectively. The extended forms are continuous: |Ãh(u,v)| ≤C ν‖u‖1,h ‖v‖1,h ∀u,v ∈V(h), (2.8) |M̃h(b,c)| ≤C |b|C |c|C ∀b,c ∈C(h), (2.9) |B̃h(u, p)| ≤C‖u‖1,h ‖p‖L2(Ω) ∀u ∈ V(h), p ∈ Q, (2.10) |D̃h(b,r)| ≤C‖b‖C ‖r‖S ∀b ∈ C(h), r ∈ S(h), (2.11) with constants C > 0 that are independent of the mesh size, ν , νm and κ . Moreover, the extended forms are constructed in such a way that Ãh(u,v) = Ah(u,v), B̃h(u, p) = Bh(u, p), M̃h(b,c) = Mh(b,c), D̃h(b,r) = Dh(b,r), (2.12) for all discrete functions u,v ∈Vh, b,c ∈ Ch, p ∈ Qh and r ∈ Sh, as well as Ãh(u,v) = A(u,v), B̃h(u, p) = B(u, p), M̃h(b,c) = M(b,c), D̃h(b,r) = D(b,r), (2.13) for all u,v ∈ V, b,c ∈ C, p ∈Q and r ∈ S. Suppose now that (u,b, p,r) is the solution of the MHD equations. We define, for any v ∈ Vh, c ∈ Ch and s ∈ Sh, the following functionals RA(v) = Ãh(u,v)+ Oh(u,v)+Ch(v,b)+ B̃h(v, p)−Fh(v), RM(c) = M̃h(b,c)−Ch(u,c)+ D̃h(c,r)−Gh(c), RD(s) = D̃h(b,s)− Jh(r,s). The terms RA, RM and RD measure how well the analytical solution satisfies the DG formulation when it is rewritten in terms of the extended bilinear forms. In- deed, if now (uh,bh, ph,rh) is the DG approximation and eu = u−uh, eb = b−bh, ep = p− ph and er = r− rh are the errors, then the following error equations can 33 be shown to hold: RA(v) = Ãh(eu,v)+ Oh(eu,v)+Ch(v,eb)+ B̃h(v,ep), RM(c) = M̃h(eb,c)−Ch(eu,c)+ D̃h(c,er), RD(s) = D̃h(eb,s)− Jh(er,s), for any v ∈ Vh, c ∈ Ch, and s ∈ Sh. We remark that the third equation of our DG method is consistent when it is rewritten in terms of the form B̃h. That is, from the definition of B̃h in [74] we see that B̃h(u−uh,q) = 0, q ∈ Qh. Proceeding as in [53, 54, 74], we readily obtain the following bound. Proposition 2.4.1 Let the solution (u,b, p,r) of the MHD problem satisfy the smooth- ness assumptions in (2.7). Then, we have |RA(v)| ≤ 2ν 1 2 ‖v‖1,hE (u,b, p), |RM(c)| ≤ |c|C⊥E (u,b, p), |RD(s)| ≤ |s|S⊥E (u,b, p) for all v ∈Vh, c ∈ Ch and s ∈ Sh, where E (u,b, p) can be bounded by E (u,b, p) ≤ Chmin{σ ,k} ( ν 1 2 ‖u‖Hσ+1(Ω) + ν− 1 2 ‖p‖Hσ (Ω) ) + Chmin{τ ,k} (‖b‖Hτ (Ω) +‖∇×b‖Hτ(Ω)) , with a constant C > 0 that is independent of the mesh size and ν . Let us also establish some continuity properties of the coupling and convection forms Ch and Oh. To that end, we need to introduce the following trace semi-norms: |u|2−1/2,Fh = ∑ K∈Th hK‖u‖2L2(∂K), |u|21/2,Fh = ∑ K∈Th h−1K ‖u‖2L2(∂K). 34 Proposition 2.4.2 The coupling form Ch satisfies |Ch(u,b)| ≤C ( ‖u‖L2(Ω) + |u|−1/2,Fh ) |b|C ∀u ∈ V(h), b ∈ C(h), |Ch(u,b)| ≤C‖u‖L2(Ω) |b|C ∀u ∈ Vh, b ∈C(h). Moreover, the convection form Oh can be bounded by |Oh(u,v)| ≤C ( ‖u‖1,h +‖u‖L2(Ω) + |u|1/2,Fh ) ‖v‖L2(Ω) for all u ∈ V(h), v ∈ Vh. The constants C > 0 are independent of the mesh size and ν . Proof: Applying the Cauchy-Schwarz inequality and taking into account the shape- regularity of the meshes, we immediately obtain the first estimate |Ch(u,b)| ≤C‖u‖L2(Ω)‖∇×b‖L2(Ω) +C ( ∑ K∈Th hK‖u‖2L2(∂K) ) 1 2 |b|C. To prove the second continuity estimate, we use the following discrete trace inequality: for any polynomial u ∈Pk(K), K ∈Th, we have ‖u‖L2(∂K) ≤Ch− 1 2 K ‖u‖L2(K). (2.14) The constant C > 0 only depends on the polynomial degree k and the shape- regularity constants of the meshes. We then readily obtain |u|−1/2,Fh ≤C‖u‖L2(Ω) ∀u ∈ Vh, from which the second estimate follows. To prove the estimate for Oh, we apply the Cauchy-Schwarz inequality to the form Oh, take into account the shape-regularity of the meshes and use the 35 bound (2.14). This results in |Oh(u,v)| ≤ ( ‖w‖L∞(Ω)‖u‖1,h +‖γ‖L∞(Ω)‖u‖L2(Ω) ) ‖v‖L2(Ω) +C‖w‖L∞(Ω) ( ∑ K∈Th h−1K ‖u‖2L2(∂K) ) 1 2 ( ∑ K∈Th hK‖v‖2L2(∂K) ) 1 2 ≤C ( ‖u‖1,h +‖u‖L2(Ω) + |u|1/2,Fh ) ‖v‖L2(Ω). This proves the estimate for Oh. 2.4.2 Error bounds Let us now denote by (u,b, p,r) the solution of the MHD problem and by (uh,bh, ph,rh) its DG approximation. We split the velocity error as follows: eu = u−uh = (u−ΠBu)+ (ΠBu−uh) = η u + ξ u, (2.15) where we use the Brezzi-Douglas-Marini (BDM) projection ΠB onto Vh∩H(div;Ω) of degree k for the approximation of the velocity, see [11, Proposition III.3.6]. This allows us to use an exactly divergence-free approximation of the velocity and to de- couple the velocity error from the pressure error which is crucial for bounding the convection and coupling terms. For the other fields, specific approximations will be chosen at a later point. For notational convenience, we introduce the product norm |||(u,b, p,r)|||2 = ‖u‖21,h +‖u‖2L2(Ω) + |u|21/2,Fh + ∑ F∈Fh ‖|w ·n| 12 [[u]]‖2L2(F) +‖b‖2C + ν−1‖p‖2L2(Ω) +‖r‖2S. Finally, we decompose bh and rh into bh = bch + b⊥h , rh = rch + r⊥h , with bch ∈ Cch, b⊥h ∈ C⊥h , rch ∈ Sch and r⊥h ∈ S⊥h , in accordance to (2.5). 36 2.4.2.1 Error in u and b We first prove two technical lemmas. Lemma 2.4.3 There are constants C > 0 and Cε > 0 independent of the mesh size and ν such that ‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ ≤ ε‖b−bh‖2C +CE (u,b, p)2 +Cε |||(η u,b− c, p−q,r− s)|||2, for any ε > 0, c ∈ Cch, q ∈ Qh, and s ∈ Sch. The constant Cε depends on ε . Proof: Fix c ∈ Cch, q ∈ Qh, s ∈ Sch and ε > 0. As in (2.15), we write eb = b−bh = (b− c)+ (c−bh)= η b + ξ b, ep = p− ph = (p−q)+ (q− ph) = ηp + ξp, (2.16) er = r− rh = (r− s)+ (s− rh) = ηr + ξr. We now proceed in the following steps. Step 1: We first observe that, since the functions bch and c are conforming in Cch, we have |b⊥h |C⊥ = |c−bch−b⊥h |C⊥ = |c−bh|C⊥ ≤ |ξ b|C. Similarly, from the conformity of rch and s in Sch, |r⊥h |S⊥ = |s− rch− r⊥h |S⊥ = |s− rh|S⊥ = |ξr|S⊥ . (2.17) Taking into account these two bounds, we have ‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ ≤ ‖ξ u‖2V + |ξ b|2C + |ξr|2S⊥ . (2.18) To bound the right-hand side above, we observe (2.12), use the stability results for Ãh + Oh, M̃h in Proposition 2.2.4, the fact that Jh(ξr,ξr) = s0|ξr|2S⊥ , and add and 37 subtract the coupling and multiplier terms. Thereby, we obtain C1 (‖ξ u‖2V + |ξ b|2C + |ξr|2S⊥) ≤ Ãh(ξ u,ξ u)+ Oh(ξ u,ξ u)+Ch(ξ u,ξ b)+ B̃h(ξ u,ξp) + M̃h(ξ b,ξ b)−Ch(ξ u,ξ b)+ D̃h(ξ b,ξr) − B̃h(ξ u,ξp)− D̃h(ξ b,ξr)+ Jh(ξr,ξr). From this estimate and the error equations in Section 2.4.1, we now readily con- clude that C1 (‖ξ u‖2V + |ξ b|2C + |ξr|2S⊥)≤ T1 + T2 + T3 + T4, (2.19) where T1 = RA(ξ u)− Ã(η u,ξ u)−Oh(η u,ξ u)−Ch(ξ u,η b)− B̃h(ξ u,ηp), T2 = RM(ξ b)− M̃h(η b,ξ b)+Ch(η u,ξ b)− D̃h(ξ b,ηr), T3 = B̃h(η u,ξp), T4 = −RD(ξr)+ D̃h(η b,ξr)− Jh(ηr,ξr). Step 2: We now bound the terms T1−T4 under the additional assumption that c belongs to the kernel Xch defined in (2.6). To bound T1, we use the estimate of RA in Proposition 2.4.1, the continuity of Ãh and B̃h in (2.8) and (2.10), respectively, and the bounds for Ch and Oh in Proposition 2.4.2. Upon application of the arithmetic-geometric mean inequality, we readily obtain that |T1| ≤C‖ξ u‖V ( E (u,b, p)+ |||(η u,η b,ηp,0)||| ) ≤ C1 2 ‖ξ u‖2V +CE (u,b, p)2 +C|||(η u,η b,ηp,0)|||2. 38 Similarly, from Proposition 2.4.1, (2.9), (2.11) and Proposition 2.4.2, we have |T2| ≤C|ξ b|C ( E (u,b, p)+ |||(η u,η b,0,0)||| ) +C‖ξ b‖C‖ηr‖S ≤C|ξ b|C ( E (u,b, p)+ |||(η u,η b,0,0)||| ) +C‖eb‖C‖ηr‖S +C‖η b‖C‖ηr‖S. Using the arithmetic-geometric mean inequality again, we have that, for all ε > 0, |T2| ≤ C12 |ξ b| 2 C + C1 2 ε‖b−bh‖2C +CE (u,b, p)2 +Cε |||(η u,η b,0,ηr)|||2. Next, we claim that T3 = 0. To see this, we note that η u = u−ΠBu belongs to H(div;Ω). It follows that [[η u]]N = 0 on interior faces. In addition, by virtue of [11, Proposition III.3.7] and since ∇ ·u = 0, we have that ∇ ·η u = 0 in Ω. Then, using the definition of B̃h in [74] and the defining properties of the BDM projection (cf. [11, Proposition III.3.6]), we conclude that T3 = B̃h(η u,ξp) = ∑ F∈FDh ∫ F η u ·nξp ds = 0; thereby, proving that T3 = 0. For the term T4, we first note, since s∈ Sch, we have Jh(ηr,ξr) = 0. Furthermore, D̃h(η b,ξr) = D̃h(η b,s− rh) = D̃h(η b,s− rch)− D̃h(η b,r⊥h ). From property (2.13), we conclude that D̃h(η b,s− rch) = D(b,s− rch)−D(c,s− rch). Both terms on the right-hand side are zero: the first one due to the fourth equa- tion in the weak formulation of the MHD problem and the second one due to the assumption that c ∈ Xch. As a consequence, we obtain T4 =−RD(ξr)− D̃h(η b,r⊥h ). 39 From Proposition 2.4.1 and the continuity of D̃h in (2.11), |T4| ≤C|ξr|S⊥E (u,b, p)+C‖η b‖C‖r⊥h ‖S. The norm-equivalence in Proposition 2.2.6 and the identity (2.17) yield ‖r⊥h ‖S ≤C|r⊥h |S⊥ = C|ξr|S⊥ . These results and the arithmetic-geometric mean inequality readily show that |T4| ≤ C12 |ξr| 2 S⊥ +C ( E (u,b, p)2 +‖η b‖2C ) . Combining (2.19) and the bounds for T1 through T4 implies that C1 2 (‖ξ u‖2V + |ξ b|2C + |ξr|2S⊥) ≤ C1 2 ε‖b−bh‖2C +CE (u,b, p)2 +Cε |||(η u,η b,ηp,ηr)|||2. Dividing the previous estimate by C12 and using (2.18) yield ‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ ≤ ε‖b−bh‖2C +CE (u,b, p)2 +Cε |||(η u,b− c, p−q,r− s)|||2, (2.20) provided that c ∈Xch. Step 3: We show that, in estimate (2.20), the approximation c ∈ Xch can be replaced by any c ∈ Cch. To that end, take c ∈Cch and look for a ∈ Cch such that D̃h(a,s) = D̃h(b− c,s) ∀s ∈ Sch. By (2.11), the right-hand side is a continuous functional on Sch. Since Xch is non- empty and D̃h = D on Cch× Sch, cf. (2.13), the inf-sup condition for D in Proposi- tion 2.2.7 implies that there exists at least one non-trivial solution a∈Cch satisfying ‖a‖C ≤C‖b− c‖C, 40 with a constant C > 0 only depending on the continuity constant of D̃h and the discrete inf-sup constant of D on Cch × Sch in Proposition 2.2.7, see [11, Equa- tion II.2.20]. By construction, we have a + c ∈ Xch, since, due to (2.13) and the fourth equation of the weak formulation of the MHD system, there holds D̃h(b,s) = D(b,s) = 0 ∀s ∈ Sch. Consequently, (a+ c) can be used as an approximation in (2.20). In addition, ‖b− (a+ c)‖C ≤ ‖b− c‖C +‖a‖C ≤C‖b− c‖C, and inequality (2.20) holds for any c ∈Cch, which completes the proof. Lemma 2.4.4 There exists a constant C > 0 independent of the mesh size and ν such that ‖b−bh‖2C ≤C (‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ + |||(η u,b− c,0,r− s)|||2), for any c ∈Cch and s ∈ Sch. Proof: Let c ∈ Cc and s ∈ Sch. Again, we split the errors in b and r into two parts and adopt the same notation as in (2.16). We now proceed in two steps. Step 1: We first consider the case where the approximation c ∈ Cch to the mag- netic field b is such that c−bch ∈ Xch. (2.21) It can be readily shown that non-trivial approximations of this type exist. To show this, consider the problem: find c ∈ Cch such that D̃h(c,s) = D̃h(bch,s) ∀s ∈ Sch. (2.22) As before, the right-hand side is a continuous functional on Sch, cf. (2.11), and the discrete inf-sup condition for D in Proposition 2.2.7, cf. (2.13), ensures that problem (2.22) admits at least one non-trivial solution c ∈ Cch which then satisfies property (2.21). 41 Now let c ∈ Cch be such that (2.21) holds. We decompose the function ξ b = c−bh into ξ b = ξ cb + ξ⊥b , ξ cb ∈ Cch, ξ ⊥b ∈C⊥h , (2.23) according to (2.5). Since the approximation c belongs to the conforming space Cch, we have ξ cb = (c−bch), ξ ⊥b =−b⊥h . (2.24) Next, we bound ‖ξ cb‖C. Due to the coercivity of M̃h on Ch, see Proposi- tion 2.2.4 and (2.12), and the fact that ξ b = ξ cb−b⊥h , we have C1‖ξ cb‖2C ≤ M̃h(ξ cb,ξ cb) = M̃h(η b + ξ b,ξ cb)− M̃h(η b,ξ cb)+ M̃h(b⊥h ,ξ cb). Using the error equation from Section 2.4.1, we obtain that M̃h(η b + ξ b,ξ cb) = RM(ξ cb)+Ch(η u + ξ u,ξ cb)− D̃h(ξ cb,ηr + ξr). The term RM(ξ cb) is zero because ξ cb ∈Cch. The term−D̃h(ξ cb,ηr +ξr) can be sim- plified as follows: since ξ cb = c−bch and s ∈ Sch, we deduce from (2.21) and (2.12) that −D̃h(ξ cb,ηr + ξr) = −D̃h(ξ cb,ηr)− D̃h(ξ cb,s− rch− r⊥h ) = −D̃h(ξ cb,ηr)+ D̃h(ξ cb,r⊥h ). From the previous discussion, we conclude that C1‖ξ cb‖2C ≤ S1 + S2, (2.25) where S1 = −M̃h(η b,ξ cb)+Ch(η u,ξ cb)− D̃h(ξ cb,ηr), S2 = M̃h(b⊥h ,ξ cb)+Ch(ξ u,ξ cb)+ D̃h(ξ cb,r⊥h ). 42 The continuity properties of M̃h, D̃h and Ch in (2.9), (2.11) and Proposition 2.4.2, respectively, and the arithmetic-geometric mean inequality yield |S1| ≤ C14 ‖ξ c b‖2C +C (|η b|2C +‖η u‖2L2(Ω) + |η u|2−1/2,Fh +‖ηr‖2S). Similarly, |S2| ≤ C14 ‖ξ c b‖2C +C (‖b⊥h ‖2C +‖ξ u‖2L2(Ω) +‖r⊥h ‖2S) ≤ C1 4 ‖ξ cb‖2C +C (|b⊥h |2C⊥ +‖ξ u‖2V + |r⊥h |2S⊥), where we have also used the norm-equivalence results in Proposition 2.2.6. Combining (2.25) with the estimates for S1 and S3, we conclude that ‖ξ cb‖2C ≤C (‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ + |||(η u,η b,0,ηr)|||2). Therefore, the previous estimate, the decomposition (2.23)–(2.24), the triangle in- equality and the norm equivalence in Proposition 2.2.6 yield that ‖b−bh‖2C ≤C(‖η b‖2C +‖ξ b‖2C) ≤C(‖η b‖2C +‖ξ cb‖2C + |b⊥h |2C⊥) ≤C(‖ξ u‖2V + |b⊥h |2C⊥ + |r⊥h |2S⊥ + |||(η u,b− c,0,r− s)|||2 ) , for any s ∈ Sch and c ∈ Cch satisfying (2.21). Step 2: We now show that the last bound of Step 1 holds for c ∈ Cch arbitrary. Proceeding as before, we can find a non-trivial function a ∈Cch such that D̃h(a,s) = D̃h(b− c−b ⊥ h ,s) ∀s ∈ Sch, ‖a‖C ≤C (‖b− c‖C +‖b⊥h ‖C) . (2.26) 43 Then, due to the properties in (2.12), (2.13) and the weak formulation of the equa- tions and the DG discretization, we have D((a+ c)−bch,s) = D̃h((a+ c)−bch,s) = D̃h(b−bh,s) = D(b,s)−Dh(bh,s) = 0, for any s ∈ Sch. Hence, (a + c)−bch ∈ Xch and a + c satisfies (2.21). It can then be used as an approximation in the last inequality of Step 1. In view of (2.26) and the norm-equivalence in Proposition 2.2.6, we obtain ‖b− (a+ c)‖C ≤ ‖b− c‖C +‖a‖C ≤C‖b− c‖C + |b⊥h |C⊥ . It follows that the last inequality of Step 1 holds for any approximation c ∈ Cch, which completes the proof of the lemma. We are now ready to bound the errors in u and b. Theorem 2.4.5 There exists a constant C > 0, independent of the mesh size and ν such that ‖u−uh‖V +‖b−bh‖C + |b⊥h |C⊥ + |r⊥h |S⊥ ≤CE (u,b, p)+C|||(η u,b− c, p−q,r− s)|||, for any c ∈Cch, q ∈Qh and s ∈ Sch. Proof: Fix c ∈ Cch, q ∈ Qh and s ∈ Sch. Decomposing the errors as in (2.15) and (2.16), we obtain from the triangle inequality, Lemma 2.4.4 and Lemma 2.4.3: ‖u−uh‖2V +‖b−bh‖2C + |b⊥h |2C⊥ + |r⊥h |2S⊥ ≤C ( ‖η u‖2V +‖ξ u‖2V +‖b−bh‖2C + |b⊥h |2C⊥ + |r⊥h |2S⊥ ) ≤C ( ‖ξ u‖2V + |||(η u,η b,0,ηr)|||2 + |b⊥h |2C⊥ + |r⊥h |2S⊥ ) ≤Cε‖b−bh‖2C +C ( E (u, p,b)2 +Cε |||(η u,b− c, p−q,r− s)|||2 ) . 44 Choosing ε = 12C and bringing the term 1 2‖b−bh‖2C to the left-hand side now read- ily implies the assertion. 2.4.2.2 Error in p and r Next, we bound the errors in the pressure p and the multiplier r. Proposition 2.4.6 There is a constant C > 0 independent of the mesh size and ν such that ‖p− ph‖L2(Ω) ≤C (E (u,b, p)+ |||(η u,b− c, p−q,r− s)|||) , for any c ∈Cch, q ∈Qh and s ∈ Sch. Proof: We begin by recalling the Poincaré inequality for piecewise smooth func- tions cf. [9, Remark 1.1]: ‖v‖L2(Ω) ≤C‖v‖1,h, v ∈ V(h). (2.27) Let now c ∈Cch, q ∈Qh and s∈ Sch. As before, we split the errors into two parts and adopt the same notation as in (2.16). Obviously, by the triangle inequality ‖p− ph‖L2(Ω) ≤ ‖ηp‖L2(Ω) +‖ξp‖L2(Ω). (2.28) We must then further estimate ‖ξp‖L2(Ω). To this end, we make use of the continu- ous inf-sup condition over V×Q; see, e.g., [11]. Therefore, we conclude that there is v ∈ V such that C‖ξp‖L2(Ω) ≤ B̃h(v,ξp) and ‖v‖H1(Ω) ≤ 1, (2.29) where we have also used (2.13). We now set vh = ΠBv with ΠB denoting the BDM projection into Vh ∩H(div;Ω). By using the definition of the extended form B̃h in [74], (2.12), (2.13) and the properties of the BDM projection, we readily obtain B̃h(vh,ξp) = B̃h(v,ξp). 45 In addition, the approximation property of the BDM projection and (2.29) guaran- tee that ‖vh‖1,h ≤ ‖v−vh‖1,h +‖v‖1,h ≤C‖v‖H1(Ω) ≤C. (2.30) We use the error equations from Section 2.4.1 to obtain B̃h(vh,ξp) = T1 + T2 + T3 + T4, (2.31) where T1 = RA(vh)− Ãh(η u,vh)−Oh(η u,vh)− B̃h(vh,ηp), T2 = −Ch(vh,b−bh), T3 = −Ãh(ξ u,vh), T4 = −Oh(ξ u,vh). Let us bound T1, T2, T3 and T4. For T1, we use Proposition 2.4.1, the continuity results for Ãh, B̃h and Oh in (2.8), (2.10) and Proposition 2.4.2, respectively, combined with the Poincaré in- equality (2.27). We obtain |T1| ≤C‖vh‖1,h ( ν 1 2 E (u,b, p)+ |||(η u,0,ηp,0)||| ) ≤C ( ν 1 2 E (u,b, p)+ |||(η u,0,ηp,0)||| ) , where we have also used (2.30). To estimate T2, we use the continuity of Ch in Proposition 2.4.2 and the Poincaré inequality (2.27): |T2| ≤C‖vh‖L2(Ω)|b−bh|C ≤C‖vh‖1,h‖b−bh‖C. From the bound for b−bh in Theorem 2.4.5 and (2.30), we thus conclude |T2| ≤C ( E (u,b, p)+ |||(η u,η b,ηp,ηr)||| ) . 46 To bound T3, we proceed similarly: we use the continuity of Ãh, the bound for ξ u in Lemma 2.4.3 (with ε = 1), the error bound b−bh in Theorem 2.4.5 and (2.30). We readily conclude that |T3| ≤Cν 1 2 (E (u,b, p)+ |||(η u,η b,ηp,ηr)|||) . The term T4 is the reason for introducing the continuous velocity field v in (2.28). To bound this term, we proceed as follows. We integrate by parts the form Oh and write Oh(ξ u,vh) = T4,1 + T4,2 + T4,3, with T4,1 = − ∑ K∈Th ∫ K (w ·∇)vh ·ξ u dx+ ∑ K∈Th ∫ K (γ−∇ ·w)ξ u ·vh dx, T4,2 = ∑ K∈Th ∫ ∂K+\Γ+ w ·nKξ u · (vh−veh)ds, T4,3 = ∫ Γ+ w ·nξ u ·vh ds, where Γ+ = {x ∈ Γ : w(x) ·n(x) ≥ 0} and ∂K+ = {x ∈ ∂K : w(x) ·nK ≥ 0}. Em- ploying the Poincaré inequality (2.27), the term T4,1 can be readily bounded by |T4,1| ≤C‖ξ u‖L2(Ω) ( ‖vh‖1,h +‖vh‖L2(Ω) ) ≤C‖ξ u‖V‖vh‖1,h. For the term T4,2, we use arguments as in the proof of Proposition 2.4.2 and the discrete trace inequality (2.14) to obtain |T4,2| ≤ C‖w‖L∞(Ω) ( ∑ K∈Th hK‖ξ u‖2L2(∂K) ) 1 2 ( ∑ F∈F Ih h−1F ‖[[vh]]‖2L2(F) ) 1 2 ≤ C‖w‖L∞(Ω)‖ξ u‖L2(Ω)‖vh‖1,h ≤C‖ξ u‖V‖vh‖1,h. Finally, the term T4,3 can be written as T4,3 = ∫ Γ+ w ·nξ u · (vh−v)ds+ ∫ Γ+ w ·nξ u ·vds≡ T4,3,1 + T4,3,2, 47 with the continuous velocity field v from (2.29). For the first integral above, we use the approximation properties of the BDM projection in [11, Proposition III.3.6] and obtain T4,3,1 ≤ (∫ Γ |w ·n||ξ u|2 ds ) 1 2 (∫ Γ |w ·n||v−vh|2 ds ) 1 2 ≤ Ch 12 ‖ξ u‖V‖v‖H1(Ω). To estimate the second integral, we use the trace theorem for functions in H1(Ω). This yields T4,3,2 ≤ (∫ Γ |w ·n||ξ u|2 ds ) 1 2 (∫ Γ |w ·n||v|2 ds ) 1 2 ≤C‖ξ u‖V‖v‖H1(Ω). As a consequence, we see that T4,3 ≤C‖ξ u‖V‖v‖H1(Ω). Hence, from the above estimates, (2.29) and (2.30), we conclude that |T4| ≤C‖ξ u‖V . From Lemma 2.4.3 (with ε = 1) and Theorem 2.4.5, we then obtain: |T4| ≤C (E (u,b, p)+ |||(b− c, p−q,r− s)|||) . The desired estimate for the pressure now follows from (2.28), (2.31) and the above estimates for T1 through T4. Finally, we bound the error in r. Proposition 2.4.7 There is a constant C > 0 independent of the mesh size and ν such that ‖r− rh‖S ≤C ( E (u,b, p)+ |||(η u,b− c, p−q,r− s)||| ) , for any c ∈Cch, q ∈Qh and s ∈ Sch. 48 Proof: Let c∈Cch, q∈Qh and s∈ Sch. As before, we adopt the notation from (2.16). By the triangle inequality, we have ‖r− rh‖S ≤ ‖ηr‖S +‖ξr‖S. To bound the term ‖ξr‖S, we decompose ξr into ξr = ξ cr + ξ⊥r , ξ cr ∈ Sch, ξ⊥r ∈ S⊥h , (2.32) according to (2.5). Since s belongs to the conforming space Sch, we have ξ cr = (s− rch), ξ⊥r =−r⊥h . (2.33) By the triangle inequality and the norm-equivalence in Proposition 2.2.6, we have ‖ξr‖S ≤ ‖ξ cr ‖S +‖r⊥h ‖S ≤ ‖ξ cr ‖S +C|r⊥h |S⊥ . The latter term can be bounded by Theorem 2.4.5: |r⊥h |S⊥ ≤C ( E (u,b, p)+ |||(η u,b− c, p−q,r− s)||| ) . To bound the former term, we use (2.13) and the inf-sup condition for D in Propo- sition 2.2.7. Thereby, we obtain C‖ξ cr ‖S ≤ sup c∈Cch\{0} D̃h(c,ξ cr ) ‖c‖C . Using (2.32) and (2.33), we write D̃h(c,ξ cr ) = D̃h(c,ηr + ξr)− D̃h(c,ηr)+ D̃h(c,r⊥h ), and use the error equation in Section 2.4.1 to conclude that D̃h(c,ξ cr ) = RM(c)− M̃h(η b,c)+Ch(η u,c) −M̃h(ξ b,c)+Ch(ξ u,c)− D̃h(c,ηr)+ D̃h(c,r⊥h ) 49 for all c ∈ Cch. We note that RM(c) = 0 for c ∈ Cch. The continuity properties and the norm-equivalence in Proposition 2.2.6 then yield |D̃h(c,ξ cr )| ≤C‖c‖C (|||(η u,η b,0,ηr)|||+‖ξ b‖C +‖ξ u‖V +‖r⊥h ‖S) ≤C‖c‖C (|||(η u,η b,0,ηr)|||+‖b−bh‖C +‖ξ u‖V + |r⊥h |S⊥). Combining the above estimates with Lemma 2.4.3 (with ε = 1) and Theorem 2.4.5 readily gives the assertion. 2.4.2.3 Proof of Theorem 2.3.1 We are now ready to complete the proof of Theorem 2.3.1. In (2.15), the approxi- mation for the velocity u has already be chosen to be the BDM projection ΠBu of degree k. We now approximate the remaining fields as follows: c = ΠNb, q = Πk−1 p, s = ΠSr, (2.34) where ΠN is the H(curl;Ω)-conforming Nédélec projection of the second kind of degree k onto Cch; see [66, 68], Πk−1 the L2-projection of degree k− 1 onto Qh and ΠS the H1-conforming nodal interpolation operator of degree k + 1 into Sh. The approximation properties of these operators immediately yield the following result. Proposition 2.4.8 Choosing the interpolants as in (2.34), there holds |||(η u,b− c, p−q,r− s)||| ≤Chmin{σ ,k} ( ‖u‖Hσ+1(Ω) + ν− 1 2 ‖p‖Hσ (Ω) ) +Chmin{τ ,k} ( ‖b‖Hτ (Ω) +‖∇×b‖Hτ(Ω) +‖r‖Hτ+1(Ω) ) , with a constant C > 0 that is independent of the mesh size and ν . The error estimate in Theorem 2.3.1 now follows directly from the error es- timates in Theorem 2.4.5, Proposition 2.4.6 and Proposition 2.4.7, in conjunction with the approximation results in Proposition 2.4.1 and Proposition 2.4.8. 50 2.5 Numerical results In this section we present a series of numerical experiments to highlight the practi- cal performance of the mixed DG method introduced in this chapter for the numer- ical approximation of incompressible MHD problems. Throughout this section, we select the stabilization parameters as follows: a0 = α k2, m0 = µ k2 and s0 = 1, α ,µ > 0, cf. [54], for example. To ensure stability of the underlying DG method we set α = µ = 10 in 2D; for 3D simulations, it is necessary to increase α and µ to α = µ = 20. All computations have been performed using the AptoFEM finite element soft- ware package; see [36] for details. In order to solve the resulting system of lin- ear equations, we have employed a variety of open-source software: for relatively small numbers of degrees of freedom, we exploit the MUltifrontal Massively Par- allel Solver (MUMPS), see [1] for details; for larger problems, we have used both the out of core version of PARDISO [72], as well as the (parallel) additive Schwarz preconditioned GMRES solver available in PETSc [6]. 2.5.1 Smooth solutions First, we verify the theoretical error bound stated in Theorem 2.3.1 for problems with smooth analytical solutions. 2.5.1.1 Example 1: 2D problem in an L-shaped domain The first example we consider is a two-dimensional version of the MHD problem (2.1). While the Navier-Stokes operator has the same form in two dimensions, some care is required for the curl-curl operator and the coupling terms in the equa- tions; see [66, Page 51] and [53] for details. We consider the L-shaped domain Ω = (−1,1)2 \ ([0,1)× (−1,0]) with ΓN = {(1,y) : y∈ (0,1)} and ΓD = Γ\ΓN , cf. Figure 2.1(a). We set ν = νm = κ = 1, w = (2,1), γ = 0, d = (x,−y), and choose the forcing functions f and g and boundary conditions so that the analytical solution of the two-dimensional variant of (2.1) is 51 y 0 −1 x −1 1 1 (a) (b) Figure 2.1: Example 1. (a) Problem domain; (b) Initial unstructured triangu- lar mesh. of the form u(x,y) = (−(ycos y+ siny)ex,ysin yex), p(x,y) = 2ex siny, b(x,y) = (−(ycos y+ siny)ex,ysin yex), r(x,y) =−sinpix sinpiy. Here, we investigate the asymptotic convergence of the interior penalty DG method on a sequence of successively finer quasi-uniform unstructured triangular meshes for k = 1,2,3,4. In each case the meshes are constructed by uniformly refining the initial mesh depicted in Figure 2.1(b). In Figure 2.2 we plot the norms ‖·‖V , ‖·‖C, and ‖·‖S of the errors eu = u−uh, eb = b− bh, and er = r− rh, respectively, against the square root of the number of degrees of freedom in the finite element space Vh ×Ch ×Qh × Sh. Here, we observe that both ‖eu‖V and ‖eb‖C converge to zero, for each fixed k, at the optimal rate O(hk), as the mesh is refined, in accordance with Theorem 2.3.1. On the other hand, for this mixed-order method, ‖er‖S converges at the rate O(hk+1), for each k, as h tends to zero; this rate is indeed optimal, though this is not reflected by Theorem 2.3.1, cf. also [54]. Additionally, in Figure 2.2(d), we plot the sum of the three error contributions with respect to the square root of the number of degrees of freedom in the finite element space. Clearly, as above, this converges to zero at the optimal rate predicted by Theorem 2.3.1. 52 102 10−6 10−4 10−2 100 1 1 1 2 1 3 1 4 k=1 k=2 k=3 k=4 ‖e u ‖ V √ Degrees of Freedom 102 10−8 10−6 10−4 10−2 100 1 1 1 2 1 3 1 4 k=1 k=2 k=3 k=4 ‖e b‖ C √ Degrees of Freedom (a) (b) 102 10−8 10−6 10−4 10−2 100 1 2 1 3 1 4 1 5 k=1 k=2 k=3 k=4 ‖e r ‖ S √ Degrees of Freedom 102 10−6 10−4 10−2 100 1 1 1 2 1 3 1 4 k=1 k=2 k=3 k=4‖e u ‖ V + ‖e b‖ C + ‖e r ‖ S √ Degrees of Freedom (c) (d) Figure 2.2: Example 1. Convergence with h-refinement: (a) ‖eu‖V ; (b) ‖eb‖C; (c) ‖er‖S; (d) ‖eu‖V +‖eb‖C +‖er‖S. Second, we highlight the optimality of the proposed mixed method when the components of the error are measured in terms of the L2-norm. From Figure 2.3 we observe that the L2-norm of the error in both the approximation to the velocity field u and the magnetic field b tend to zero at the expected optimal rate O(hk+1), for each k, as h tends to zero. In agreement with Theorem 2.3.1, for each fixed k, the L2-norm of the error in the pressure p, denoted by ep = p− ph, tends to zero at the optimal rate O(hk) as the mesh is enriched, while ‖er‖L2(Ω) is of order O(hk+2) as h tends to zero. 53 102 10−10 10−8 10−6 10−4 10−2 100 1 2 1 3 1 4 1 5 k=1 k=2 k=3 k=4 ‖e u ‖ L 2 ( Ω ) √ Degrees of Freedom 102 10−6 10−4 10−2 100 1 1 1 2 1 3 1 4 k=1 k=2 k=3 k=4 ‖e p‖ L2 (Ω ) √ Degrees of Freedom (a) (b) 102 10−10 10−8 10−6 10−4 10−2 100 1 2 1 3 1 4 1 5 k=1 k=2 k=3 k=4 ‖e b‖ L2 (Ω ) √ Degrees of Freedom 102 10−10 10−8 10−6 10−4 10−2 100 1 3 1 4 1 5 1 6 k=1 k=2 k=3 k=4 ‖e r ‖ L 2 ( Ω ) √ Degrees of Freedom (c) (d) Figure 2.3: Example 1. Convergence with h-refinement: (a) ‖eu‖L2(Ω); (b) ‖ep‖L2(Ω); (c) ‖eb‖L2(Ω); (d) ‖er‖L2(Ω). 2.5.1.2 Example 2: 3D problem in the unit cube The second example is a 3D problem with a smooth analytical solution. Here, we set Ω = (0,1)3 ⊂R3 with ΓD = Γ and ΓN = /0, ν = νm = κ = 1, w = (2,1,1), γ = 0, and d = (x,−y,1), and select f and g, together with appropriate inhomogeneous boundary conditions, so that the solution of the incompressible MHD system (2.1) is given by u = (−(ycos y+ siny+ zcos z)ex,ysin yex,zsin zex), b = (−(ycos y+ siny+ zcos(z))ex,ysin yex,zsin zex), p = 2ex(sin y+ sinz)− p0, r = sinpix sinpiy sinpiz, 54 k DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖V l ‖ep‖L2(Ω) l 576/48 6.400e-2 – 1.168 – 7.321e-1 – 1 4,608/384 1.647e-2 1.96 5.823e-1 1.00 4.564e-1 0.68 36,864/3,072 4.213e-3 1.97 2.896e-1 1.01 2.606e-1 0.81 294,912/24,576 1.072e-3 1.97 1.442e-1 1.01 1.400e-1 0.90 1,440/192 5.082e-3 – 1.262e-1 – 3.424e-1 – 2 11,520/1,536 6.802e-4 2.90 3.162e-2 2.00 8.488e-2 2.01 92,160/12,288 8.417e-5 3.01 7.822e-3 2.02 2.127e-2 2.00 Table 2.1: Example 2. Convergence of ‖eu‖L2(Ω), ‖eu‖V , and ‖ep‖L2(Ω). where p0 = 4(−1+ e+ cos1− ecos1). In Table 2.1 we investigate the asymptotic rate of convergence of the error in the approximation of the hydrostatic variables; here, l denotes the computed rate of convergence. To this end, we show ‖eu‖L2(Ω), ‖eu‖V , and ‖ep‖L2(Ω) computed on a sequence of uniformly refined tetrahedral meshes for k = 1,2. As in the previous example, we again observe optimal rates of convergence for all three measures of the error. Indeed, in accordance with Theorem 2.3.1, both ‖eu‖V and ‖ep‖L2(Ω) tend to zero at the optimal rate O(hk), for each fixed k, as the mesh is refined. Additionally, we observe that ‖eu‖L2(Ω) is of optimal order O(hk+1) as h tends to zero. The corresponding errors for the magnetic variables are shown in Tables 2.2 & 2.3. Here, we clearly observe the optimality of the approximation to the magnetic field b. Indeed, from Table 2.2 we observe that ‖eb‖L2(Ω) and ‖eb‖C converge to zero at the optimal rates O(hk+1) and O(hk), respectively, for each fixed k, as the mesh is refined. As in the previous example, we again observe that ‖er‖L2(Ω) and ‖er‖S are of order O(hk+2) and O(hk+1), respectively, as the mesh is uniformly refined. 2.5.2 Example 3: 2D problem with a singular solution To verify the ability of the proposed interior penalty DG method to capture the strongest magnetic (and hydrostatic) singularities, we consider a problem in which the precise regularity of the analytical solution is known. To this end, we again 55 k DOFs bh ‖eb‖L2(Ω) l ‖eb‖C l 576 7.289e-2 – 7.324e-1 – 1 4,608 2.076e-2 1.81 3.445e-1 1.09 36,864 5.486e-3 1.92 1.668e-1 1.05 294,912 1.399e-3 1.97 8.184e-2 1.03 1,440 6.082e-3 – 4.920e-2 – 2 11,520 7.953e-4 2.94 1.146e-2 2.10 92,160 1.006e-4 2.98 2.767e-3 2.05 Table 2.2: Example 2. Convergence of ‖eb‖L2(Ω) and ‖eb‖C. k DOFs rh ‖er‖L2(Ω) l ‖er‖S l 480 1.038e-1 – 1.546 – 1 3,840 1.766e-2 2.56 5.098e-1 1.60 30,720 2.350e-3 2.91 1.363e-1 1.90 245,760 2.924e-4 3.01 3.405e-2 2.00 960 1.327e-2 – 2.559e-1 – 2 7,680 8.567e-4 3.95 3.430e-2 2.90 61,440 5.135e-5 4.06 4.210e-3 3.03 Table 2.3: Example 2. Convergence of ‖er‖L2(Ω) and ‖er‖S. let Ω be the L-shaped domain employed in Example 1 above with ΓN = {(1,y) : y ∈ (0,1)} and ΓD = Γ \ΓN . We choose ν = νm = κ = 1, and set w = 0, γ = 0 and d = (−1,1). Hence, the Navier-Stokes operator coincides with the Stokes equations. We further choose f and g, and appropriate inhomogeneous boundary conditions so that the solution to this problem is given by the strongest corner singularities for the underlying elliptic operators. That is, in polar coordinates (ρ ,φ) around the origin, the hydrostatic solution components u and p are taken to be u(ρ ,φ) = ρλ ((1+ λ )sin(φ)ψ(φ)+ cos(φ)ψ ′(φ)) ρλ (−(1+ λ )cos(φ)ψ(φ)+ sin(φ)ψ ′(φ)) , p(ρ ,φ) =−ρλ−1((1+ λ )2ψ ′(φ)+ ψ ′′′(φ))/(1−λ ), 56 where ψ(φ) =sin((1+ λ )φ)cos(λω)/(1+ λ )− cos((1+ λ )φ) − sin((1−λ )φ)cos(λω)/(1−λ )+ cos((1−λ )φ), with λ ≈ 0.54448373678246 and ω = 3pi/2. The magnetic pair (b,r) is taken as b(x) = ∇(ρ2/3 sin(2/3φ)), r(x) ≡ 0. We point out that the magnetic field b does not belong to H1(Ω)2 and thus cannot be correctly captured by nodal elements. In fact, for this example, we have that (u, p) ∈ H1+λ (Ω)2×Hλ (Ω) and b ∈ H2/3(Ω)2. Thus, the limiting regularity ex- ponent, cf. (2.7a) and (2.7b), appearing in Theorem 2.3.1 is λ which stems from the regularity of the hydrostatic variables. In this example we study the asymptotic convergence of the interior penalty DG method on the sequence of successively finer quasi-uniform unstructured triangular meshes employed in Example 1, cf. Figure 2.1(b), with k = 1,2,3. Table 2.4 presents the L2-norm of the error in both the computed velocity uh and pressure ph, as well as the ‖ · ‖V -norm of the error in uh. In agreement with Theorem 2.3.1 we see that both ‖eu‖V and ‖ep‖L2(Ω) tend to zero at the optimal rate O(hλ ) as h tends to zero. The rate of convergence of ‖eu‖L2(Ω) is observed to be between O(h1.2) and O(h1.5) approximately as the mesh is uniformly refined. From Table 2.5 we observe that both ‖eb‖L2(Ω) and ‖eb‖C are approximately O(h) as h tends to zero. For this latter error, this rate is higher than what we would expect from Theorem 2.3.1. However, this same behaviour of the error has also been observed in the case of simply approximating the time-harmonic Maxwell operator in isolation, cf. [54]. In contrast, from Table 2.6, we observe that ‖er‖S converges to zero at the rate O(h2/3) as the mesh is refined. In terms of the numerical approximation of the time-harmonic Maxwell operator in isolation, this rate is indeed optimal, cf. [54], though this is not reflected in Theorem 2.3.1. Finally, we note that the L2-norm of the error in the approximation to the variable r tends to zero at the rate O(h4/3) as h tends to zero. 57 k DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖V l ‖ep‖L2(Ω) l 144/24 1.311e-1 – 1.910 – 1.443 – 576/96 4.638e-2 1.50 1.352 0.50 1.064 0.44 1 2,304/384 1.632e-2 1.51 9.419e-1 0.52 7.690e-1 0.47 9,216/1,536 5.837e-3 1.48 6.510e-1 0.53 5.436e-1 0.50 36,864/6,144 2.120e-3 1.46 4.482e-1 0.54 3.789e-1 0.52 288/72 6.089e-2 – 1.075 – 1.520 – 1,152/288 2.405e-2 1.34 7.382e-1 0.54 9.010e-1 0.75 2 4,608/1,152 9.434e-3 1.35 5.065e-1 0.54 5.852e-1 0.62 18,432/4,608 3.837e-3 1.30 3.474e-1 0.54 3.910e-1 0.58 73,728/18,432 1.630e-3 1.24 2.383e-1 0.54 2.649e-1 0.56 480/144 3.094e-2 – 7.498e-1 – 9.219e-1 – 1,920/576 1.198e-2 1.37 5.151e-1 0.54 5.809e-1 0.67 3 7,680/2,304 4.844e-3 1.31 3.532e-1 0.54 3.779e-1 0.62 30,720/9,216 2.054e-3 1.24 2.422e-1 0.54 2.527e-1 0.58 122,880/36,864 9.046e-4 1.18 1.661e-1 0.54 1.713e-1 0.56 Table 2.4: Example 3. Convergence of ‖eu‖L2(Ω), ‖eu‖V , and ‖ep‖L2(Ω). Figure 2.4: Example 4. Initial unstructured triangular mesh. 2.5.3 Hartmann channel flow Finally, we consider 2D and 3D Hartmann channel flow problems; cf. [35]. Note that assumption (2.2) is not satisfied in these examples. However, the particular structure of the solutions implies that (w ·∇)u = 0 and assumption (2.2) is not relevant in this context. 58 k DOFs bh ‖eb‖L2(Ω) l ‖eb‖C l 144 2.601e-1 – 4.091e-1 – 576 1.492e-1 0.80 2.291e-1 0.84 1 2,304 7.699e-2 0.96 1.112e-1 1.04 9,216 4.038e-2 0.93 5.280e-2 1.07 36,864 2.265e-2 0.84 2.657e-2 0.99 288 2.244e-1 – 3.649e-1 – 1,152 1.124e-1 1.00 1.785e-1 1.03 2 4,608 5.371e-2 1.07 8.065e-2 1.15 18,432 2.679e-2 1.00 3.654e-2 1.14 73,728 1.452e-2 0.88 1.766e-2 1.05 480 1.888e-1 – 3.090e-1 – 1,920 9.013e-2 1.07 1.448e-2 1.09 3 7,680 4.156e-2 1.12 6.363e-2 1.19 30,720 2.004e-2 1.05 2.812e-2 1.18 122,880 1.055e-2 0.93 1.320e-2 1.09 Table 2.5: Example 3. Convergence of ‖eb‖L2(Ω) and ‖eb‖C. 2.5.3.1 Example 4: 2D Hartmann flow In the domain Ω = (0,L)× (−1,1), L ≫ 1, we consider the steady 2D unidirec- tional flow under a constant pressure gradient −G in the x-direction. We set w= ( G νHatanh(Ha) ( 1− cosh(yHa) cosh(Ha) ) ,0 ) , d= ( G κ ( sinh(yHa) sinh(Ha) −y ) ,1 ) , γ = 0, f = 0, and g = 0. Additionally, we impose the boundary conditions u = 0 on y =±1, (pI−ν∇u)n = pNn on x = 0 and x = L, n×b = n×bD on Γ, where bD = (0,1), pN =−Gx− G 2 2κ ( sinh(yHa) sinh(Ha) − y )2 + p0, 59 k DOFs rh ‖er‖L2(Ω) l ‖er‖S l 144 2.397e-1 – 2.107 – 576 1.150e-1 1.06 1.768 0.25 1 2,304 4.860e-2 1.24 1.265 0.48 9,216 1.946e-2 1.32 8.384e-1 0.59 36,864 7.664e-3 1.34 5.387e-1 0.64 240 1.944e-1 – 2.728 – 960 8.588e-2 1.18 2.066 0.40 2 3,840 3.498e-2 1.30 1.412 0.55 15,360 1.382e-2 1.34 9.193e-1 0.62 61,440 5.419e-3 1.35 5.868e-1 0.65 360 1.621e-1 – 3.188 – 1,440 6.932e-2 1.23 2.323 0.46 3 5,760 2.784e-2 1.32 1.559 0.58 23,040 1.095e-2 1.35 1.008 0.63 92,160 4.290e-3 1.35 6.415e-1 0.65 Table 2.6: Example 3. Convergence of ‖er‖L2(Ω) and ‖er‖S. and p0 is any constant. The analytical solution to the incompressible MHD equa- tions is given by u = w, b = d, p = pN , r ≡ 0, where κ = ννmHa2. We note that the fluid always moves in the direction in which the pressure decreases. We set L = 10, ν = νm = 0.1, Ha = 10, G = 0.5, and p0 = 10. First, in Figure 2.5 we investigate the asymptotic convergence of the interior penalty DG method on a sequence of successively finer quasi-uniform unstructured triangular meshes for k = 1,2,3. In each case the meshes are constructed by uni- formly refining the initial unstructured mesh depicted in Figure 2.4. Here, we plot the norms ‖ · ‖V , ‖ · ‖C, ‖ · ‖S, and ‖ · ‖L2(Ω) of the errors eu, eb, er, and ep, respec- tively, with respect to the square root of the number of degrees of freedom in the finite element space Vh×Ch×Qh×Sh. As in the previous examples presented in Section 2.5.1, we observe that ‖eu‖V , ‖eb‖C and ‖ep‖L2(Ω) converge to zero, for each fixed k, at the optimal rate O(hk) as the mesh is refined, in accordance with Theorem 2.3.1, while ‖er‖S converges at the rate O(hk+1), for each k, as h tends to zero. Moreover, we note that the L2-norms of the error in the approximation to u, b and r tend to zero optimally, cf. Section 2.5.1; for brevity, these results have 60 102 103 10−3 10−2 10−1 100 1 1 1 2 1 3 k=1 k=2 k=3 ‖e u ‖ V √ Degrees of Freedom 102 103 10−3 10−2 10−1 100 1 1 1 2 1 3 k=1 k=2 k=3 ‖e b‖ C √ Degrees of Freedom (a) (b) 102 103 10−4 10−3 10−2 10−1 100 1 2 1 31 4 k=1 k=2 k=3 ‖e r ‖ S √ Degrees of Freedom 102 103 10−3 10−2 10−1 100 1 1 1 21 3 k=1 k=2 k=3 ‖e p‖ L2 (Ω ) √ Degrees of Freedom (c) (d) Figure 2.5: Example 4. Convergence with h-refinement: (a) ‖eu‖V ; (b) ‖eb‖C; (c) ‖er‖S; (d) ‖ep‖L2(Ω). been omitted. Finally, in Figures 2.6 & 2.7 we show the DG solution computed on the finest mesh with 41216 elements, employing k = 1; thereby, the total number of degrees of freedom employed in the finite element space Vh ×Ch ×Qh × Sh is 783104. In particular, from Figure 2.7, we observe extremely good agreement between the computed and analytical solutions of the first components in the velocity and mag- netic fields. 2.5.3.2 Example 5: 3D Hartmann flow In this final example, we consider the steady 3D unidirectional flow in the rectan- gular duct given by Ω = [0,L]× [−y0,y0]× [−z0,z0] with y0,z0 ≪ L. We take w = u 61 (a) (b) Figure 2.6: Example 4. DG solution computed on the finest mesh with k = 1: (a) Velocity field; (b) Magnetic field. −1 −0.5 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 DG Solution Analytical Solution −1 −0.5 0 0.5 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 DG Solution Analytical Solution (a) (b) Figure 2.7: Example 4. DG solution computed on the finest mesh with k = 1. Slices along x = 5, −1 ≤ y ≤ 1 of the solution: (a) First component of the velocity field; (b) First component of the magnetic field. (cf. below), f = g = 0, γ = 0, d = (0,1,0), and consider solutions of the form u = (u(y,z),0,0), b = (b(y,z),1,0), p =−Gx+ p0, r ≡ 0. 62 (a) (b) Figure 2.8: Example 5. DG solution computed on a uniform tetrahedral mesh with k = 1: (a) Velocity field; (b) Magnetic field. We enforce the boundary conditions u = 0 for y =±y0 and z =±z0, (pI−ν∇u)n = pNn for x = 0 and x = L, n×b = n×bD on Γ, with pN = −Gx + p0 and bD = (0,1,0). As before, G and p0 are arbitrary con- stants. For this channel problem, the analytical solution can be expressed by Fourier series; for details, we refer to [34]. Here, we set L = 10, y0 = 2, z0 = 1, 63 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 DG Solution Analytical Solution −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 DG Solution Analytical Solution (a) (b) Figure 2.9: Example 5. DG solution computed on a uniform tetrahedral mesh with k = 1. Slices along x = 5,−1≤ y≤ 1, z = 0 of the solution: (a) First component of the velocity field; (b) First component of the magnetic field. ν = νm = 0.1, Ha = 5, G = 0.5, and p0 = 10. In Figures 2.8 & 2.9 we show the DG solution computed on a uniform tetra- hedral mesh comprising of 30720 elements with the polynomial degree k = 1; this results in a total of 1075200 degrees of freedom in the finite element space Vh×Ch×Qh×Sh. In particular, from Figure 2.9, we observe that there is reason- ably good agreement between the computed and analytical solutions of the first components in the velocity and magnetic fields on this relatively coarse mesh. However, here we do observe some over-shoots in the computed solution, which are particularly evident in the approximation to the magnetic field. 2.6 Conclusions In this chapter, we have proposed and analyzed a mixed DG method for a linear incompressible magnetohydrodynamics problem. We have derived a priori error estimates and computationally verified them with a set of numerical examples. We have further tested the methods for channel flow problems in both two- and three- dimensions. As mentioned in Remark 2.3.3, our analysis can be easily adapted to the case where conforming Nédélec elements are used for the approximation of the magnetic unknowns. It can also be extended to the divergence-conforming 64 Pdk −Pk−1 elements proposed in [20] for the approximation of the fluid variables. The systematic development of efficient solution methods and preconditioners for the proposed discretization is the subject of ongoing research. Indeed, a wide variety of efficient and robust saddle point solvers are available in the literature for both the Oseen operator and the Maxwell operator. However, these approaches need to be extended to MHD problems of the form (2.1). We also refer the reader to [48] and [73, Section 3.4] for discussions of iterative strategies that amount to the solution of a sequence of decoupled Navier-Stokes and Maxwell problems. Ongo- ing work also includes extensions of the method to fully nonlinear incompressible MHD systems, which will be presented in Chapter 3. 65 Chapter 3 A mixed finite element method with exactly divergence-free velocities for nonlinear incompressible magnetohydrodynamics We introduce and analyze a mixed finite element method for the numerical dis- cretization of a stationary incompressible magnetohydrodynamics problem, in two and three dimensions. The velocity field is discretized using divergence-conforming Brezzi-Douglas-Marini (BDM) elements and the magnetic field is approximated by curl-conforming Nédélec elements. The H1-continuity of the velocity field is enforced by a DG approach. Central features of the method are that it produces ex- actly divergence-free velocity approximations, and captures the strongest magnetic singularities. We prove that the energy norm error is convergent in the mesh size in general Lipschitz polyhedra under minimal regularity assumptions, and derive nearly optimal a priori error estimates for the two-dimensional case. We present a comprehensive set of numerical experiments, which indicate optimal convergence of the proposed method for two-dimensional as well as three-dimensional prob- 66 lems. 3.1 Introduction We consider a standard form of the incompressible MHD equations as derived in [3, Section 2]; see also [33, 35, 48]. That is, we neglect phenomena involving high frequency as well as the convection current, and consider a non-polarizable, non-magnetizable and homogeneous medium. In addition, to make the curl-curl operator arising in the Maxwell equations amenable to discretization with Nédélec elements, we use the mixed formulation proposed in [73]. The governing equations are then of the form −ν ∆u +(u ·∇)u + ∇p−κ (∇×b)×b = f in Ω, (3.1a) ∇ ·u = 0 in Ω, (3.1b) κνm ∇× (∇×b)+ ∇r−κ ∇× (u×b) = g in Ω, (3.1c) ∇ ·b = 0 in Ω. (3.1d) Here, u is the velocity, b the magnetic field, p the hydrostatic pressure, and r is a Lagrange multiplier associated with the divergence constraint on the magnetic field b. The functions f and g represent external force terms. For simplicity, we consider the following pure Dirichlet (i.e., Γ = ΓD) homogeneous boundary condi- tions: u = 0 on Γ, (3.2a) n×b = 0 on Γ, (3.2b) r = 0 on Γ, (3.2c) with n being the unit outward normal on Γ. By taking the divergence of the mag- netostatic equation (3.1c), we obtain the Poisson problem ∆r = ∇ ·g in Ω, r = 0 on Γ. (3.3) 67 Since g is typically divergence-free in physical applications, the multiplier r is typically zero and its primary purpose is to ensure stability; see [29, Section 3]. Various finite element methods for discretizing linear and nonlinear MHD sys- tems can be found in the literature. The magnetic field is often approximated by standard nodal (i.e., H1-conforming) finite elements [3, 33, 45, 46, 48]. However, since the strongest magnetic singularities have regularity below H1, straightfor- wardly applied nodal elements may fail to resolve them in non-convex polyhedral domains; see [23] and the references therein. A number of remedies have been proposed for electromagnetic problems, for example the weighted regularization approach in [24] or the approach in [7], whereby the divergence of the electric field is stabilized in H−α with 12 < α < 1. In [50], weighted regularization has been applied to a full incompressible MHD system. In the mixed formulation of [73] the above mentioned difficulties associated with nodal elements are seamlessly avoided without the need for stabilizing the divergence. This approach amounts to introducing the Lagrange multiplier r, and yields the PDE system (3.1). As a result, it is possible to use curl-conforming Nédélec elements for approximating the magnetic field. For these elements, only tangential continuity is enforced across inter-elemental faces. This makes this ap- proach feasible in situations of highly singular magnetic fields [52, 66, 68]. In the context of incompressible magnetohydrodynamics, a related mixed approach for the discretization of the magnetic unknowns has been presented in [30]. We are interested in discretizations for incompressible MHD problems that are based on discontinuous Galerkin methods; see, e.g., the surveys [14, 16, 28] and the references therein. In [46], an interior penalty technique is applied to enforce continuity of the magnetic variable across domains with different electromagnetic properties, while nodal elements are employed in the interior. A full DG method is proposed in [56] for a linearized variant of the system (3.1), whereby all the vari- ables are approximated in discontinuous finite element spaces, based on existing discretizations for the Oseen and Maxwell equations [18, 19, 54]. However, this approach requires a large number of degrees of freedom. Furthermore, a straight- forward extension to the nonlinear setting in a locally conservative fashion would require a post-processing procedure for smoothing the DG velocity approximations throughout the nonlinear iteration [19]. 68 In this chapter we design a new finite element discretization, in an attempt to overcome the above mentioned difficulties. Instead of discontinuous elements for all unknowns, we use divergence-conforming Brezzi-Douglas-Marini (BDM) el- ements [11, 20] for the approximation of the velocity field, and curl-conforming Nédélec elements [68] for the magnetic field, thereby substantially reducing the total number of the coupled degrees of freedom. The H1-continuity of the velocity field is again enforced by a DG technique. A central feature of this discretiza- tion is that it yields exactly divergence-free velocity approximations, guaranteeing stability of the linearized system within each Picard iteration, without any other modifications. We note that divergence-conforming discretizations have been an- alyzed for the incompressible Navier-Stokes equations in [20]. For the magnetic approximation we have a discrete version of the desirable property (3.3), in contrast to the method presented in [56]. We prove well-posedness of our discretization, and show convergence under minimal regularity assumptions. Thus, our method captures the strongest mag- netic singularities in non-convex polyhedra. Our numerical results clearly indicate optimal convergence rates in two and three dimensions, but we manage to show (nearly) optimal estimates only for the two-dimensional case. Specific details on this are given in Section 3.4 and are summarized in the conclusions in Section 3.6. We note that our method converges optimally for the linearized version of (3.1), as follows from the arguments in [56, Remark 3.3]. The rest of the chapter is structured as follows. In Section 3.2 we state the well- posedness of the variational formulation of (3.1). Section 3.3 is devoted to the finite element discretization; the existence and uniqueness of approximate solutions are proved. In Section 3.4 we present and prove the main results–convergence and a priori error estimates. In Section 3.5 we present a series of numerical experiments validating the theoretical results. In Section 3.6 we end with some concluding remarks. 69 3.2 Variational formulation of an MHD problem Upon setting V = H10 (Ω)d = {v ∈ H1(Ω)d : v = 0 on Γ}, Q = L20(Ω) = {q ∈ L2(Ω) : (q ,1)Ω = 0}, C = H0(curl;Ω) and S = H10 (Ω), the variational formulation of the incompressible MHD system (3.1)–(3.2) amounts to finding (u,b, p,r) ∈V×C×Q×S such that A(u,v)+ O(u,u,v)+C(b,v,b)+ B(v, p) = (f,v)Ω, (3.4a) B(u,q) = 0, (3.4b) M(b,c)−C(b,u,c)+ D(c,r) = (g,c)Ω, (3.4c) D(b,s) = 0, (3.4d) for all (v,c,q,s) ∈ V×C×Q×S. The variational forms are given by A(u,v) = ∫ Ω ν ∇u : ∇vdx, O(w,u,v) = ∫ Ω (w ·∇)u ·vdx, M(b,c) = ∫ Ω κνm(∇×b) · (∇× c)dx, C(d,v,b) = ∫ Ω κ (v×d) · (∇×b)dx, B(u,q) =− ∫ Ω (∇ ·u)qdx, D(b,s) = ∫ Ω b ·∇sdx. To discuss the well-posedness of the mixed formulation (3.4), we introduce the product norms ‖(u,b)‖V×C = ( ν‖u‖2H1(Ω) + κνm‖b‖2H(curl;Ω) ) 1 2 , (u,b) ∈ V×C, ‖(p,r)‖Q×S = ( 1 ν ‖p‖2L2(Ω) + 1 κνm ‖r‖2H1(Ω) ) 1 2 , (p,r) ∈ Q×S. 70 Here, the curl-norm is defined by ‖b‖H(curl;Ω) = ( ‖b‖2L2(Ω) +‖∇×b‖2L2(Ω) ) 1 2 . Furthermore, we define the norm of the source terms by |||(f,g)|||= ( ‖f‖2L2(Ω) +‖g‖2L2(Ω) ) 1 2 . Finally, we introduce the parameters ¯ν = min{ν ,κνm}, ¯κ = max{1,κ}. The following result can be found in [73, Corollary 2.18 and Remark 2.14]. Theorem 3.2.1 There is a constant c1 > 0 only depending on Ω such that for small data with c1 ¯κ ¯ν−2|||(f,g)||| < 1, the MHD problem (3.4) has a unique so- lution (u,b, p,r) ∈V×C×Q×S. Moreover, we have the stability bound ‖(u,b)‖V×C ≤ c2 |||(f,g)||| ¯ν 1 2 , for a constant c2 > 0 only depending on Ω. 3.3 Mixed finite element discretization In this section, we introduce a mixed finite element method that employs divergence- conforming elements for the approximation of the velocity field and curl-conforming elements for the magnetic field. The H1-continuity of the velocity is enforced by a DG technique. 3.3.1 Mixed discretization Recall the notation on meshes and traces in Section 2.2.2. We assume in addition that the triangulation Th is quasi-uniform. For k ≥ 1, we wish to approximate the 71 solution of (3.1)–(3.2) by finite element functions (uh,bh, ph,rh)∈Vh×Ch×Qh× Sh, where Vh = {v ∈ H0(div;Ω) : v|K ∈Pk(K)d , K ∈Th }, Ch = {c ∈ H0(curl;Ω) : c|K ∈Pk−1(K)d ⊕Rk(K), K ∈Th }, Qh = {q ∈ L20(Ω) : q|K ∈Pk−1(K), K ∈Th }, Sh = {s ∈ H10 (Ω) : s|K ∈Pk(K), K ∈Th }. (3.5) Here, we denote by H0(div;Ω) the space H0(div;Ω) = { v ∈ L2(Ω)d : ∇ ·v ∈ L2(Ω), v ·n = 0 on Γ } , and by Rk(K) the space of homogeneous vector polynomials of total degree k that are orthogonal to x. The space Vh is the divergence-conforming Brezzi-Douglas-Marini (BDM) space (see [11, Section III.3] for details); it has degrees of freedom specified for the normal components of functions along faces. The space Ch represents the first family of curl-conforming Nédélec elements (cf. [68, Chapter 5]); its degrees of freedom are defined for the tangential components of functions along faces. We notice that the finite element spaces Ch, Qh and Sh are conforming in C, Q and S, respectively, while Vh is non-conforming in V. Now we consider the following finite element method: find (uh,bh, ph,rh) ∈ Vh×Ch×Qh×Sh such that Ah(uh,v)+ Oh(uh,uh,v)+C(bh,v,bh)+ B(v, ph) = (f,v)Ω, (3.6a) B(uh,q) = 0, (3.6b) M(bh,c)−C(bh,uh,c)+ D(c,rh) = (g,c)Ω, (3.6c) D(bh,s) = 0, (3.6d) for all (v,c,q,s) ∈ Vh×Ch×Qh×Sh. The form Ah associated with the Laplacian is chosen as the standard interior 72 penalty form [4, 5]: Ah(u,v) = ∑ K∈Th ∫ K ν∇u : ∇vdx− ∑ F∈Fh ∫ F {{ν∇u}} : [[v]]ds − ∑ F∈Fh ∫ F {{ν∇v}} : [[u]]ds+ ∑ F∈Fh a0ν hF ∫ F [[u]] : [[v]]ds. Here, a0 > 0 is the interior penalty stabilization parameter; it has to be chosen larger than a threshold value which is independent of h, ν , κ and νm. For the convection form, we take the standard upwind form [62]: Oh(w,u,v) = ∑ K∈Th ∫ K (w ·∇)u ·vdx + ∑ K∈Th ∫ ∂K\Γ 1 2 (w ·nK−|w ·nK|)(ue−u) ·vds − ∫ Γ 1 2 (w ·n−|w ·n|)u ·vds. Here, ue is the trace of u taken from the exterior of K. The remaining forms are the same as in the continuous case. Notice that due to the presence of the upwind terms the form Oh(w,u,v) is not linear in the first argument; see also Lemma 3.4.6 and (3.8). By choosing the divergence-conforming BDM elements as the approximating space for the velocity, the method gives exactly divergence-free velocity approx- imations; cf. [20]. Moreover, the Lagrange multiplier rh vanishes identically for divergence-free source terms, thereby mimicking the continuous property in (3.3). Proposition 3.3.1 Let (uh, bh, ph, rh) solve (3.6). Then we have: (i) ∇ ·uh = 0 in Ω. (ii) the Lagrange multiplier rh is the solution of (∇rh,∇s)Ω = (g,∇s)Ω ∀ s ∈ Sh. In particular, if g is solenoidal, then rh ≡ 0. 73 Proof: To prove item (i), we proceed as in [20]. We note that ∇ ·uh has vanishing mean value on Ω, and is a discontinuous polynomial of degree k− 1. Thus, we have ∇ · uh ∈ Qh. Equation (3.6b) then implies that ∇ · uh is orthogonal to all functions q ∈ Qh. Therefore, it is equal to zero. To prove item (ii), we take c = ∇s in equation (3.6c) (noting that ∇Sh ⊂ Ch) and obtain (g,∇s)Ω = M(bh,∇s)−C(bh,uh,∇s)+ D(∇s,rh) = D(∇s,rh). Here, we have used the fact that ∇×∇s = 0. Therefore, rh satisfies (∇rh,∇s)Ω = (g,∇s)Ω ∀ s ∈ Sh. Since (g,∇s)Ω = (∇ ·g,s)Ω, we have rh ≡ 0 provided that ∇ ·g = 0. For our analysis, it will be convenient to introduce the following product forms: Ah(u,b;v,c) = Ah(u,v)+ M(b,c), Oh(w,d;u,b;v,c) = Oh(w,u,v)+C(d,v,b)−C(d,u,c), B(u,b;q,s) = B(u,q)+ D(b,s), L (v,c) = (f,v)Ω +(g,c)Ω. Then, the mixed discretization (3.6) is equivalent to the following saddle-point system: find (uh,bh, ph,rh) ∈ Vh×Ch×Qh×Sh such that Ah(uh,bh;v,c)+Oh(uh,bh;uh,bh;v,c)+B(v,c; ph,rh) = L (v,c), B(uh,bh;q,s) = 0 for all (v,c,q,s) ∈ Vh×Ch×Qh×Sh. 74 3.3.2 Stability properties To discuss the stability properties of the finite element formulation (3.6), we intro- duce the discrete H1-norm for the hydrostatic velocity: ‖u‖1,h = ( ∑ K∈Th ‖∇u‖2L2(K) + ∑ F∈Fh h−1F ‖[[u]]‖2L2(F) ) 1 2 . (3.7) We further define ‖(u,b)‖Vh×Ch = ( ν‖u‖21,h + κνm‖b‖2H(curl;Ω) ) 1 2 . First, we note that the forms Ah and B are continuous over the finite element spaces: |Ah(u,b;v,c)| ≤CA ‖(u,b)‖Vh×Ch ‖(v,c)‖Vh×Ch , |B(u,b; p,r)| ≤CB‖(u,b)‖Vh×Ch‖(p,r)‖Q×S, for all u,v∈Vh, b,c∈Ch, p ∈Qh, r ∈ Sh, with constants CA ,CB > 0 independent of h, ν , κ and νm. Next, we introduce the following spaces of (discretely) divergence-free func- tions: Jh = {u ∈ Vh : B(u,q) = 0 ∀q ∈ Qh }, Xh = {b ∈ Ch : D(b,s) = 0 ∀ s ∈ Sh }. For the form Oh, we then have the following continuity result: there exists a constant CO > 0 independent of h, ν , κ and νm such that, for any w1,w2 ∈ Vh, u,v ∈ Vh, d1,d2 ∈ Xh, and b,c ∈ Ch, we have |Oh(w1,d1;u,b;v,c)−Oh(w2,d2;u,b;v,c)| ≤ CO ¯κ ¯ν 3 2 ‖(w1−w2,d1−d2)‖Vh×Ch ‖(u,b)‖Vh×Ch ‖(v,c)‖Vh×Ch ; (3.8) 75 see also Proposition 3.4.8 for a more detailed discussion. Furthermore, the following stability properties of Ah and Oh hold; cf. [5, 18, 52, Theorem 4.7] and the references therein: Ah(u,b;u,b)≥CC ‖(u,b)‖2Vh×Ch ∀ (u,b) ∈ Vh×Xh, (3.9) Oh(w,d;u,b;u,b) = Oh(w,u,u)≥ 0 ∀w ∈ Jh, u ∈ Vh, b,d ∈ Ch, (3.10) with a constant CC > 0 independent of h, ν , κ and νm. Finally, let us address the inf-sup stability of the forms B and D. For the form B we have the following result [49, Proposition 10]: inf p∈Qh\{0} sup v∈Vh\{0} B(v, p) ‖v‖1,h‖p‖L2(Ω) = λh ≥C > 0, (3.11) where C is independent of h, ν , κ and νm. Moreover, since ∇Sh ⊂ Ch, there holds [54, Lemma 5.3]: inf r∈Sh\{0} sup c∈Ch\{0} D(c,r) ‖c‖H(curl;Ω)‖r‖H1(Ω) = µh ≥C > 0, (3.12) for a constant C independent of h, ν , κ and νm. An immediate consequence of (3.11) and (3.12) is the following inf-sup condi- tion for the product form B: inf (p,r)∈Qh×Sh\{(0,0)} sup (v,c)∈Vh×Ch\{(0,0)} B(v,c; p,r) ‖(v,c)‖Vh×Ch ‖(p,r)‖Q×S ≥CS > 0, (3.13) where the stability constant CS is independent of h, ν , κ and νm. In Table 3.1, we show the discrete inf-sup constants λh in (3.11) for the velocity- pressure pair Vh×Qh defined in (3.5). We use the lowest-order BDM elements on Ω = (−1,1)2 and compute the discrete inf-sup constants λh for a sequence of suc- cessively refined uniform triangular meshes. The inf-sup constants are obtained by solving a generalized eigenvalue problem related to the matrix representation of the bilinear form B and the norms in (3.11); cf. [11, page 75]. Table 3.1 illustrates that the discrete inf-sup constants are approaching a positive lower bound as the 76 mesh is refined. DOFs uh/ph 112/32 416/128 1,600/512 6,272/2,048 24,832/8,192 λh 1.273e-1 1.251e-1 1.241e-1 1.236e-1 1.233e-1 Table 3.1: Discrete inf-sup constants for Vh×Qh. 3.3.3 Existence and uniqueness of discrete solutions In the following theorem, we state the unique solvability of the method (3.6) under a discrete version of the smallness assumption in Theorem 3.2.1. The proof of this result follows along the same lines as [73, Theorem 2.12], using the stability properties outlined in Section 3.3.2. Theorem 3.3.2 There is a constant C1 > 0 independent of h, ν , κ and νm such that for small data with C1 ¯κ ¯ν−2|||(f,g)||| < 1, the mixed finite element discretiza- tion (3.6) has a unique solution (uh,bh, ph,rh) ∈ Vh ×Ch ×Qh × Sh. Moreover, there is a constant C2 > 0 independent of h, ν , κ and νm such that ‖(uh,bh)‖Vh×Ch ≤C2 |||(f,g)||| ¯ν 1 2 . The solution of (3.6) can be found by employing the following Picard-type iteration: given (un−1h ,b n−1 h ) ∈Vh×Ch, let (unh,bnh, pnh,rnh) in Vh×Ch×Qh×Sh be the solution of the linearized Oseen-type problem Ah(unh,v)+ Oh(un−1h ,u n h,v)+C(bn−1h ,v,b n h)+ B(v, p n h) = (f,v)Ω,(3.14a) B(unh,q) = 0, (3.14b) M(bnh,c)−C(bn−1h ,unh,c)+ D(c,rnh) = (g,c)Ω,(3.14c) D(bnh,s) = 0, (3.14d) for all (v,c,q,s) ∈ Vh×Ch×Qh×Sh. Theorem 3.3.2 guarantees the convergence of the iterates {(unh,bnh, pnh,rnh)}n≥0 to the solution (uh,bh, ph,rh) of (3.6) for any initial guess (u0h,b0h) ∈ Vh ×Ch 77 with exactly divergence-free u0h, provided that the small data assumption in The- orem 3.3.2 is satisfied. However, the scheme is only linearly convergent, as we illustrate in Section 3.5. Remark 3.3.3 A more efficient nonlinear solver such as Newton’s method can also be used for solving (3.6); see, e.g., [34, 35, 48]. When upwinding is not incorporated, Newton’s method can be straightforwardly applied. However, when upwind terms are included, adapting the nonlinear iteration to our discretization is more delicate, since it requires additional linearization of the convection form Oh(w,u,v) in the first argument. This remains an item for future investigation. 3.4 Error analysis In this section, we present the main results of this chapter, namely the conver- gence of finite element approximations and a priori error estimates for the two- dimensional version of our MHD problem. We provide detailed proofs in Sec- tions 3.4.2 through 3.4.5. 3.4.1 Main results Our first result is a convergence result. To state it, we suppose the solution (u,b, p,r) of (3.1)–(3.2) possesses the smoothness (u, p) ∈ Hσ+1(Ω)d ×Hσ(Ω), (3.15a) (b,∇×b,r) ∈ Hτ(Ω)d ×Hτ(Ω)d ×Hτ+1(Ω), (3.15b) for σ ,τ > 12 . Remark 3.4.1 The regularity assumption (3.15b) is minimal in the sense that it is satisfied by the strongest singularities of the Maxwell operator in polyhedral domains; cf. [23, 24]. Similarly, the regularity (3.15a) holds true for the strongest singularities of the Stokes operator in polyhedral domains; see [2, 26]. In view of these results, we expect (3.15) to be the minimal smoothness of solutions to the MHD system (3.1)-(3.2) in general Lipschitz polyhedra. However, we do not have a full proof of this conjecture. 78 Theorem 3.4.2 Let (u,b, p,r) and (uh,bh, ph,rh) be the solutions of (3.1)–(3.2) and (3.6), respectively, obtained on a sequence of quasi-uniform meshes {Th}h>0 of mesh size h. Assume (3.15) and that ¯κ ¯ν−2|||(f,g)||| is sufficiently small. Then we have lim h→0 ‖(u−uh,b−bh)‖Vh×Ch = 0, limh→0‖(p− ph,r− rh)‖Q×S = 0. Theorem 3.4.2 guarantees that the method (3.6) gives correct solutions pro- vided that the (minimal) smoothness assumption (3.15) is satisfied and the data is sufficiently small. In particular, it ensures convergence in situations where straight- forwardly applied nodal elements for the approximation of b are not capable of correctly capturing the singular solution components. Next, we present a priori error estimates for the two-dimensional version of the MHD problem (3.6). Theorem 3.4.3 Let Ω ⊂ R2 be a simply-connected Lipschitz polygon with a con- nected boundary Γ. Under the same assumption as in Theorem 3.4.2, we have the following error estimates for any ε > 0: ‖(u−uh,b−bh)‖Vh×Ch ≤Cε hmin{σ ,τ ,k}−ε ( ν 1 2 ‖u‖Hσ+1(Ω) +(κνm) 1 2 ‖b‖Hτ (Ω) +(κνm) 1 2 ‖∇×b‖Hτ (Ω) ) +Chmin{σ ,τ ,k} ( ν− 1 2 ‖p‖Hσ (Ω) +(κνm)− 1 2 ‖r‖Hτ+1(Ω) ) , and ‖(p− ph,r− rh)‖Q×S ≤Cε hmin{σ ,τ ,k}−2ε ( ν 1 2 ‖u‖Hσ+1(Ω) +(κνm) 1 2 ‖b‖Hτ (Ω) +(κνm) 1 2 ‖∇×b‖Hτ (Ω) ) +Cεhmin{σ ,τ ,k}−ε ( ν− 1 2 ‖p‖Hσ (Ω) +(κνm)− 1 2 ‖r‖Hτ+1(Ω) ) . Here, the constants C and Cε are independent of h, ν , κ and νm. While Cε depends on ε , the constant C does not. 79 The convergence rates in Theorem 3.4.3 are optimal in the mesh size, up to a loss of O(hε ) for ε arbitrarily small. This loss stems from the use of the Sobolev embedding of H1(Ω) into Lp(Ω), for all p ≥ 1, but not into L∞(Ω); cf. [39]. To bridge this gap, we use inverse estimates to establish the continuity of the nonlinear coupling form; see the proof of Lemma 3.4.7. In addition, the constant Cε might become unbounded as ε tends to zero. However, in our numerical experiments this constant is observed to stay bounded. In fact, we observe optimal rates of convergence in all our tests, for both smooth and non-smooth solutions. Full details are given in Section 3.5. Remark 3.4.4 Our technique of proof is applicable to three-dimensional prob- lems. However, since in three dimensions the Sobolev embeddings are more re- strictive, the use of the inverse estimates leads to convergence rates that fall short half a power of h for the error in u and b, and a full power of h for the error in p and r (i.e., Theorem 3.4.3 holds with ε = 12 ). To see this, we carry out the proof of Theorem 3.4.3 simultaneously for d = 2 and d = 3. We emphasize, however, that in our numerical tests, optimal convergence rates are observed for three-dimensional problems with smooth solutions. Remark 3.4.5 For the linearized variant of the MHD system (3.1), our method converges optimally in the mesh size h, as follows from [56, Remark 3.3]. That is, the estimates of Theorem 3.4.3 hold true without any loss, both in two and three dimensions. However, there we make stronger smoothness assumptions on the linearized magnetic field. Therefore, this optimality cannot be straightforwardly carried over to the nonlinear setting. The proofs of Theorems 3.4.2 and 3.4.3 are presented in the next four subsec- tions. 3.4.2 Continuity We begin by revisiting the continuity properties of the forms in a more general setting. To that end, we introduce the space V(h) = V + Vh, 80 and endow it with the norm ‖·‖1,h. We then make use of an auxiliary form Ãh(u,v) constructed as in [19, 74] via the use of suitable lifting operators. It is defined as Ãh(u,v) = ∫ Ω ν ( ∇hu : ∇hv−L (v) : ∇hu−L (u) : ∇hv ) dx + ∑ F∈Fh a0ν hF ∫ F [[u]] : [[v]]ds. Here, L : V(h)→ Σh = {σ ∈ L2(Ω)d×d : σ |K ∈Pk(K)d×d , K ∈Th } is the lifting operator given by∫ Ω L (u) : σ dx = ∑ F∈Fh ∫ F [[u]] : {{σ}}ds ∀σ ∈ Σh, and ∇h is the elementwise gradient operator. By construction, the form Ãh(u,v) satisfies Ãh(u,v) = A(u,v) ∀u,v ∈ V, Ãh(u,v) = Ah(u,v) ∀u,v ∈ Vh. (3.16) Furthermore, using arguments similar to those in [19, 74], the form Ãh(u,v) can be shown to be bounded on V(h)×V(h). Then, by setting Ãh(u,b;v,c) = Ãh(u,v)+ M(b,c), we readily obtain |Ãh(u,b;v,c)| ≤C‖(u,b)‖Vh×Ch ‖(v,c)‖Vh×Ch (3.17) for u,v ∈V(h) and b,c ∈ C. Moreover, we have |B(v,c;q,s)| ≤C‖(v,c)‖Vh×Ch‖(q,s)‖Q×S (3.18) for (v,c,q,s) ∈ V(h)×C×Q×S. In (3.17)–(3.18) and in the following, we denote by C a generic (positive) con- stant that is independent of the mesh size h and the parameters ν , κ and νm. 81 Next, we state the continuity of the convection term. The proof of this result follows similarly to the ones in [60, Proposition 4.15] and [19, Proposition 4.2]. Lemma 3.4.6 There holds: |Oh(w1,u,v)−Oh(w2,u,v)| ≤C‖w1−w2‖1,h ‖u‖1,h ‖v‖1,h for all w1,w2 ∈ V(h), and u,v ∈V(h). In the sequel, we shall analyze the two- and three-dimensional cases simulta- neously (see also Remark 3.4.4). To do so, we introduce the function ℓ(d) given by ℓ(d) = h −ε , d = 2, h− 12 , d = 3. Here, ε > 0 is a fixed number. The function ℓ(d) will indicate the loss of con- vergence rates for both the two-dimensional and three-dimensional cases. We also denote by Cd > 0 a generic constant independent of h, ν , κ and νm, but depen- dent on the dimension d. In particular, for d = 2 it depends on ε and might be unbounded as ε → 0. By introducing the kernel X ={b ∈ C : D(b,s) = 0 ∀ s ∈ S} , we state and prove the continuity of the coupling form C(d,u,c) for several cases. Lemma 3.4.7 There holds: (i) for d ∈ X∪Xh,u ∈ Vh and c ∈ C: |C(d,u,c)| ≤Cκ ‖d‖H(curl;Ω) ‖u‖1,h ‖c‖H(curl;Ω). (ii) for d ∈ X∪Xh,u ∈ V(h) and c ∈ Ch: |C(d,u,c)| ≤Cκ ‖d‖H(curl;Ω) ‖u‖1,h ‖c‖H(curl;Ω). 82 (iii) for d ∈C,u ∈ V and c ∈Ch: |C(d,u,c)| ≤Cd ℓ(d)κ ‖d‖L2(Ω)‖u‖H1(Ω) ‖c‖H(curl;Ω) . (iv) for d ∈ C,u ∈Vh and c ∈ C: |C(d,u,c)| ≤Cd ℓ(d)κ ‖d‖L2(Ω) ‖u‖1,h ‖c‖H(curl;Ω). Proof: We proceed in two steps. Step 1. We first discuss preliminary results that will be used in the proof. From the Poincaré inequality in [52, Corollary 4.4], there holds ‖b‖L2(Ω) ≤C‖∇×b‖L2(Ω) ∀b ∈ X. (3.19) Next, we recall the inverse inequality (cf. [10, Lemma 4.5.3]): for any u ∈ Pk(K), there holds ‖u‖Ln1 (K) ≤Ch d( 1n1− 1 n2 ) K ‖u‖Ln2 (K) ∀1≤ n1,n2 ≤ ∞. (3.20) Further, let H : Xh →X be the Hodge operator that maps discretely divergence- free functions into exactly divergence-free functions in such a way that ∇×Hd = ∇×d. (3.21) It satisfies the following approximation property (cf. [52, Lemma 4.5]): there ex- ists τ > 12 such that ‖d−Hd‖L2(Ω) ≤Chτ‖∇×d‖L2(Ω) ∀d ∈ Xh. (3.22) 83 Finally, we present the following Sobolev embeddings: ‖u‖Lm(d)(Ω) ≤Cd ‖u‖H1(Ω) ∀u ∈ H1(Ω)d , (3.23a) ‖u‖Lm(d)(Ω) ≤Cd ‖u‖1,h ∀u ∈ V(h), (3.23b) ‖d‖L3(Ω) ≤C‖d‖H(curl;Ω) ∀d ∈ X. (3.23c) Here, m(2) = 2/ε and m(3) = 6. The embedding (3.23a) is a standard result, while the embedding (3.23b) follows similarly to [40, 61]. Inequality (3.23c) follows from [2, Proposition 3.7]. Step 2. We are now ready to prove the bounds in the lemma. The proof of inequality (i) can be found in [73, Proposition 3.2]. To establish the second inequality, we follow [73, Lemma 2.6] and first show it for d ∈ X,u ∈ V(h) and c ∈ Ch. This is done by applying Hölder’s inequality and the Sobolev embeddings (3.23b) and (3.23c). We obtain |C(d,u,c)| ≤ κ‖d‖L3(Ω)‖u‖L6(Ω)‖∇× c‖L2(Ω) ≤Cκ‖d‖H(curl;Ω)‖u‖1,h‖c‖H(curl;Ω) . (3.24) Second, if d ∈ Xh, we decompose it into d = (d−Hd)+ Hd, where H is the Hodge operator in (3.21). We then rewrite C(d,u,c) as C(d,u,c) = C(d−Hd,u,c)+C(Hd,u,c). (3.25) Because Hd ∈ X, we can apply the previous argument, (3.24), and bound the last term of (3.25) by |C(Hd,u,c)| ≤Cκ‖Hd‖H(curl;Ω)‖u‖1,h‖c‖H(curl;Ω) ≤Cκ‖d‖H(curl;Ω)‖u‖1,h‖c‖H(curl;Ω). (3.26) In (3.26), we have used the Poincaré inequality (3.19) and property (3.21) of the 84 Hodge operator. For the first term on the right-hand side of (3.25), we obtain from Hölder’s inequality, the Sobolev embedding (3.23b) and the approximation property (3.22) that |C(d−Hd,u,c)| ≤ κ‖d−Hd‖L2(Ω)‖u‖L6(Ω)‖∇× c‖L3(Ω) ≤Cκhτ‖∇×d‖L2(Ω)‖u‖1,h‖∇× c‖L3(Ω), for τ > 12 . Finally, we apply the inverse estimate (3.20) to achieve |C(d−Hd,u,c)| ≤Cκ‖d‖H(curl;Ω)‖u‖1,h‖∇× c‖L2(Ω), (3.27) for both d = 2 and d = 3. Referring to (3.25), (3.26) and (3.27) proves the assertion of item (ii). To verify item (iii), we define m∗(d) such that 1 m(d) + 1 m∗(d) = 1 2 . Then we apply Hölder’s inequality, the Sobolev embedding (3.23a) and the inverse estimate (3.20) to conclude that |C(d,u,c)| ≤ κ‖d‖L2(Ω)‖u‖Lm(d)(Ω)‖∇× c‖Lm∗(d)(Ω) ≤Cd κ‖d‖L2(Ω)‖u‖H1(Ω)ℓ(d)‖∇× c‖L2(Ω). The proof of item (iv) is similar to that of item (iii): |C(d,u,c)| ≤ κ‖d‖L2(Ω)‖u‖L∞(Ω)‖∇× c‖L2(Ω) ≤Cd κ‖d‖L2(Ω)ℓ(d)‖u‖Lm(d)(Ω)‖c‖H(curl;Ω) . Applying (3.23b) finishes the proof. As an immediate consequence of Lemmas 3.4.6 and 3.4.7, the form Oh satisfies the following continuity properties. 85 Proposition 3.4.8 There is a constant CO > 0 independent of h, ν , κ and νm such that there holds: (i) for w1,w2 ∈ V(h), d1,d2 ∈ X∪Xh, (u,b) ∈ V(h)×C, and (v,c) ∈ Vh× Ch: |Oh(w1,d1;u,b;v,c)−Oh(w2,d2;u,b;v,c)| ≤ CO ¯κ ¯ν 3 2 ‖(w1−w2,d1−d2)‖Vh×Ch ‖(u,b)‖Vh×Ch ‖(v,c)‖Vh×Ch . (ii) for w1,w2 ∈ V(h), d1,d2 ∈ C, (u,b) ∈ V×C, and (v,c) ∈ Vh×Ch: |Oh(w1,d1;u,b;v,c)−Oh(w2,d2;u,b;v,c)| ≤Cd ℓ(d)CO ¯κ ¯ν 3 2 ‖(w1−w2,d1−d2)‖Vh×Ch ‖(u,b)‖V×C ‖(v,c)‖Vh×Ch . 3.4.3 Preliminary error estimates In this subsection, we present two lemmas for estimating the errors. Let (u,b, p,r) and (uh,bh, ph,rh) be the solutions of (3.1)–(3.2) and (3.6), respectively. We begin by defining the residual RA(v) = Ãh(u,v)+ Oh(u,u,v)+C(b,v,b)+ B(v, p)− (f,v)Ω (3.28) for any v ∈ Vh. It measures how well the exact solution satisfies the finite element formulation expressed in terms of the auxiliary form Ãh in (3.16). We have the following upper bound for the residual (cf. [74]): RA(v)≤ ν 1 2 ‖v‖1,hE (u) with E (u)≤Chmin{σ ,k}ν 1 2 ‖u‖Hσ+1(Ω). (3.29) In the following, we shall denote the errors by eu = u−uh, eb = b−bh, ep = p− ph, er = r− rh. 86 We shall also decompose the errors into eu = η u + ξ u = (u−v)+ (v−uh), eb = η b + ξ b = (b− c)+ (c−bh), ep = ηp + ξp = (p−q)+ (q− ph), er = ηr + ξr = (r− s)+ (s− rh), (3.30) for a discrete function (v,c,q,s) ∈Vh×Ch×Qh×Sh to be specified later. Lemma 3.4.9 Assume that max{c2,C2}COCC ¯κ|||(f,g)||| ¯ν2 < 1 2 . (3.31) Then there holds ‖(u−uh,b−bh)‖Vh×Ch ≤Cd ℓ(d) inf (v,c)∈Vh×Ch ‖(u−v,b− c)‖Vh×Ch +C ( E (u)+ inf (q,s)∈Qh×Sh ‖(p−q,r− s)‖Q×S ) . Proof: We proceed in two steps. Step 1. In the error decomposition (3.30), we first consider (v,c) ∈ Jh ×Xh. Clearly, we also have (ξ u,ξ b) ∈ Jh×Xh. In view of the residual equation (3.28), we obtain RA(ξ u) =Ãh(eu,eb;ξ u,ξ b)+Oh(u,b;u,b;ξ u,ξ b) −Oh(uh,bh;uh,bh;ξ u,ξ b)+B(ξ u,ξ b;ep,er) =Ãh(eu,eb;ξ u,ξ b)+Oh(u,b;u,b;ξ u,ξ b)−Oh(uh,bh;u,b;ξ u,ξ b) +Oh(uh,bh;eu,eb;ξ u,ξ b)+B(ξ u,ξ b;ep,er). Because uh ∈ Jh (see Proposition 3.3.1), the stability of Oh in (3.10) guarantees that Oh(uh,bh;ξ u,ξ b;ξ u,ξ b)≥ 0. 87 Therefore, we have Ãh(ξ u,ξ b;ξ u,ξ b)+Oh(v,c;u,b;ξ u,ξ b)−Oh(uh,bh;u,b;ξ u,ξ b) ≤RA(ξ u)− Ãh(η u,η b;ξ u,ξ b)−B(ξ u,ξ b;ep,er) +Oh(v,c;u,b;ξ u,ξ b)−Oh(u,b;u,b;ξ u,ξ b) −Oh(uh,bh;η u,η b;ξ u,ξ b). (3.32) From the coercivity of Ah in (3.9) and the continuity of Oh in Proposition 3.4.8 (i), the left-hand side of equation (3.32) can be bounded by l.h.s. of (3.32)≥CC ‖(ξ u,ξ b)‖2Vh×Ch − CO ¯κ ¯ν 3 2 ‖(u,b)‖Vh×Ch ‖(ξ u,ξ b)‖2Vh×Ch . Next, we estimate ‖(u,b)‖Vh×Ch using the stability bound in Theorem 3.2.1 (noting that ‖(u,b)‖Vh×Ch ≤ ‖(u,b)‖V×C). We obtain l.h.s. of (3.32)≥ ( CC− c2CO ¯κ|||(f,g)||| ¯ν2 ) ‖(ξ u,ξ b)‖2Vh×Ch . In view of assumption (3.31), we then have l.h.s. of (3.32)≥ 1 2 CC ‖(ξ u,ξ b)‖2Vh×Ch . For the right-hand side of (3.32), we note that (since ξ u and ξ b are in the kernels Jh and Xh, respectively) B(ξ u,ξ b;ep,er) = B(ξ u,ξ b;ηp,ηr). Then, to bound the right-hand side of (3.32), we use the continuity properties of Ãh, B and Oh in (3.17), (3.18) and Proposition 3.4.8, respectively, as well as 88 the estimate for RA(ξ u) in (3.29). We readily obtain r.h.s. of (3.32)≤‖(ξ u,ξ b)‖Vh×Ch ( E (u)+C‖(η u,η b)‖Vh×Ch +‖(ηp,ηr)‖Q×S +Cd ℓ(d) CO ¯κ ¯ν 3 2 ‖(η u,η b)‖Vh×Ch ‖(u,b)‖V×C + CO ¯κ ¯ν 3 2 ‖(uh,bh)‖Vh×Ch ‖(η u,η b)‖Vh×Ch ) . Next, we employ the stability bounds in Theorems 3.2.1 and 3.3.2 for ‖(u,b)‖V×C and ‖(uh,bh)‖Vh×Ch , respectively, apply the small data assumption (3.31), and com- bine the lower and upper bounds of (3.32) into the estimate ‖(ξ u,ξ b)‖Vh×Ch ≤Cd ℓ(d)‖(η u,η b)‖Vh×Ch +C ( E (u)+‖(ηp,ηr)‖Q×S ) . From the triangle inequality, we thus obtain the error bound ‖(u−uh, b−bh)‖Vh×Ch ≤Cd ℓ(d)‖(u−v,b− c)‖Vh×Ch +C ( E (u)+‖(p−q,r− s)‖Q×S ) , (3.33) for any (v,c) ∈ Jh×Xh and (q,s) ∈Qh×Sh. Step 2. Next, we replace (v,c) ∈ Jh×Xh in (3.30) by (v,c) ∈Vh×Ch. To that end, let (v,c) ∈Vh×Ch, and we look for (w,d) ∈Vh×Ch such that B(w,d;q,s) = B(u−v,b− c;q,s) ∀ (q,s) ∈ Qh×Sh. Since the right-hand side is a continuous functional on Qh×Sh, we conclude from the inf-sup condition of B in (3.13) that there exists at least one non-trivial solu- tion (w,d) of this problem satisfying the bound ‖(w,d)‖Vh×Ch ≤C‖(u−v,b− c)‖Vh×Ch . By construction, (w + v,d + c) ∈ Jh×Xh. Therefore, it can be inserted into (3.33). 89 With the help of the triangle inequality, we readily see that ‖(u−uh,b−bh)‖Vh×Ch ≤ ‖(u−v,b− c)‖Vh×Ch +‖(w + v−uh,d + c−bh)‖Vh×Ch +‖(w,d)‖Vh×Ch ≤Cd ℓ(d)‖(u−v,b− c)‖Vh×Ch +C ( E (u)+‖(p−q,r− s)‖Q×S ) . This completes the proof. Next, we present the following result for the multipliers. Lemma 3.4.10 Assume (3.31). Then there holds ‖(p− ph,r− rh)‖Q×S ≤C ( E (u)+ inf (q,s)∈Qh×Sh ‖(p−q,r− s)‖Q×S +‖(u−uh,b−bh)‖Vh×Ch + sup (v,c)∈Vh×Ch\{(0,0)} |C(b−bh,v,b)−C(b−bh,u,c)| ‖(v,c)‖Vh×Ch ) . Proof: For any (q,s) ∈ Qh×Sh, we recall from (3.30) that ep = ξp + ηp, er = ξr + ηr. The inf-sup condition for B in (3.13) and the triangle inequality guarantee that ‖(ξp,ξr)‖Q×S ≤C sup (v,c)∈Vh×Ch\{(0,0)} B(v,c;ξp,ξr) ‖(v,c)‖Vh×Ch ≤C(T1 + T2), where T1 = sup (v,c)∈Vh×Ch\{(0,0)} B(v,c;ηp,ηr) ‖(v,c)‖Vh×Ch , T2 = sup (v,c)∈Vh×Ch\{(0,0)} B(v,c;ep,er) ‖(v,c)‖Vh×Ch . 90 Using the continuity of B in (3.18), T1 can be easily bounded by T1 ≤C‖(ηp,ηr)‖Q×S. For T2, we make use of the weak formulation and the residual equation (3.28) and write out the form Oh into its individual parts. We obtain B(v,c;ep,er) = RA(v)− Ãh(eu,eb;v,c)−Oh(uh,bh;eu,eb;v,c) −Oh(u,b;u,b;v,c)+Oh(uh,bh;u,b;v,c) = RA(v)− Ãh(eu,eb;v,c)−Oh(uh,bh;eu,eb;v,c) −Oh(u,u,v)+ Oh(uh,u,v)−C(eb,v,b)+C(eb,u,c). Applying the bound (3.29) and the continuity properties of Ãh, Oh and Oh in (3.17), Proposition 3.4.8 (i) and Lemma 3.4.6, respectively, we conclude that T2 ≤ E (u)+C‖(eu,eb)‖Vh×Ch + CO ¯κ ¯ν 3 2 ‖(eu,eb)‖Vh×Ch‖(uh,bh)‖Vh×Ch + CO ¯κ ¯ν 3 2 ‖(eu,0)‖Vh×Ch‖(u,0)‖Vh×Ch + sup (v,c)∈Vh×Ch\{(0,0)} |C(eb,v,b)−C(eb,u,c)| ‖(v,c)‖Vh×Ch . Using the small data assumption (3.31), the assertion follows. 3.4.4 Proof of Theorem 3.4.2 In this subsection, we prove the convergence result stated in Theorem 3.4.2. In view of Lemma 3.4.9, the convergence of uh and bh is obtained under the smoothness assumption (3.15) by using the standard approximation properties of the finite element spaces Vh, Ch, Qh and Sh, respectively. This proves the first statement of Theorem 3.4.2. Next, we show the convergence of the multipliers in the energy norm ‖ · ‖Q×S. From Lemmas 3.4.10 and 3.4.9, it only remains to show that sup (v,c)∈Vh×Ch\{(0,0)} |C(b−bh,v,b)−C(b−bh,u,c)| ‖(v,c)‖Vh×Ch → 0 as h→ 0. (3.34) 91 Recalling the Hodge operator H from (3.21), we write C(b−bh,v,b) = C(b−Hbh,v,b)+C(Hbh−bh,v,b). (3.35) The first term on the right-hand side of (3.35) tends to zero due to Lemma 3.4.7 (i) and the fact that ‖b−Hbh‖H(curl;Ω) ≤ ‖b−bh‖H(curl;Ω) +‖bh−Hbh‖H(curl;Ω) → 0 as h→ 0. Here, we have applied the triangle inequality, the properties of the Hodge operator in (3.21) and (3.22), and the stability bound in Theorem 3.3.2. For the last term of (3.35), we first utilize item (iv) of Lemma 3.4.7 and then the approximation result (3.22), to get |C(Hbh−bh,v,b)| ≤Ch− 1 2 κ‖Hbh−bh‖L2(Ω)‖v‖1,h‖b‖H(curl;Ω) ≤Chτ− 12 κ ‖∇×bh‖L2(Ω)‖v‖1,h‖b‖H(curl;Ω) ≤Chτ− 12 κ ¯ν 3 2 ‖(0,bh)‖Vh×Ch‖(v,0)‖Vh×Ch‖(0,b)‖V×C ≤Chτ− 12 ‖(v,0)‖Vh×Ch |||(f,g)||| ¯ν 1 2 . In the last step, we have applied the stability bounds in Theorems 3.2.1 and 3.3.2, as well as the small data assumption (3.31). Since ¯ν− 12 |||(f,g)||| ≤ ¯κ ¯ν−2|||(f,g)||| and τ > 12 , we obtain sup (v,c)∈Vh×Ch\{(0,0)} |C(b−bh,v,b)| ‖(v,c)‖Vh×Ch → 0 as h→ 0. A similar argument shows that sup (v,c)∈Vh×Ch\{(0,0)} |C(b−bh,u,c)| ‖(v,c)‖Vh×Ch → 0 as h→ 0. Therefore (3.34) holds true, and the convergence of the multipliers is obtained. 92 3.4.5 Proof of Theorem 3.4.3 In this subsection we prove the a priori error estimates in Theorem 3.4.3. As before, we consider the cases d = 2 and d = 3 simultaneously. Based on Lemma 3.4.9, we choose v as the BDM projection of u, c the Nédélec projection of b, q and s the L2-projections of p and r, respectively. We then apply the approximation properties of these projections in [11, Proposition III.3.6], [66, Theorem 5.41] and [13], and the estimates for the errors in the velocity and mag- netic fields are readily obtained. To prove the error estimate for the multipliers, we first apply Proposition 3.4.8 (ii) to bound the supremum in the estimate of Lemma 3.4.10: sup (v,c)∈Vh×Ch\{(0,0)} |C(b−bh,v,b)−C(b−bh,u,c)| ‖(v,c)‖Vh×Ch ≤Cd ℓ(d)CO ¯κ ¯ν 3 2 ‖(0,b−bh)‖Vh×Ch‖(u,b)‖V×C. Utilizing the stability bound in Theorem 3.2.1, we obtain ‖(p− ph,r− rh)‖Q×S ≤C ( E (u)+ inf (q,s)∈Qh×Sh ‖(p−q,r− s)‖Q×S ) +Cd ℓ(d)‖(u−uh,b−bh)‖Vh×Ch ≤Cd ℓ(d)2 inf (v,c)∈Vh×Ch ‖(u−v,b− c)‖Vh×Ch +Cd ℓ(d) ( E (u)+ inf (q,s)∈Qh×Sh ‖(p−q,r− s)‖Q×S ) . Again, we choose v as the BDM projection of u, c the Nédélec projection of b, q and s the L2-projections of p and r, respectively. As before, the approximation properties of these projections finish the proof. 3.5 Numerical results In this section we present a series of numerical experiments. Our computations have been carried out using MATLAB, with direct linear solvers. The primary 93 purpose of our experiments is to confirm optimal convergence rates of our method. We start by considering one problem with a smooth solution and a second one with a singular solution. Then, we consider the numerical approximations of two- and three-dimensional Hartmann channel flow and driven cavity flow problems. Finally, we present results for another benchmark problem: MHD flow over a step in two dimensions. Throughout this section, the lowest-order BDM and Nédélec elements are em- ployed and the interior penalty stabilization parameter is a0 = 10. The Picard it- eration described in Section 3.4.5 is used to solve the nonlinear systems. For all the examples, we solve a Stokes problem and the Maxwell equations, decoupled, to obtain an initial guess. The tolerance for the Picard iterations is chosen as 1e-5. We test our method on problems with mixed Dirichlet and Neumann boundary conditions in the hydrostatic variables, even though the analysis has been carried out solely for the Dirichlet case. 3.5.1 Example 1: two-dimensional problem with a smooth solution First, we verify the theoretical results stated in Theorems 3.4.2 and 3.4.3 for a problem with a smooth analytical solution. We consider the following two-dimensional problem. We set Ω = (−1,1)2 with ΓN = {(1,y) : y ∈ (−1,1)}, ΓD = Γ\ΓN , ν = κ = 1, νm = 1e4, and choose the source terms f, g and the boundary conditions so that the analytical solution is of the form u(x,y) = (y2,x2), p(x,y) = x, b(x,y) = (1− y2,1− x2), r(x,y) = (1− x2)(1− y2). We construct this example with r 6= 0 to show the convergence rate in rh; later examples will feature a divergence-free g and a vanishing r; cf. Proposition 3.3.1. In Tables 3.2–3.4, we investigate the asymptotic rates of convergence of the errors in the approximations of the hydrostatic and magnetic variables; here, l de- notes the experimental rate of convergence. We observe that ‖u− uh‖1,h, ‖p− ph‖L2(Ω), ‖b−bh‖H(curl;Ω) and ‖∇(r− rh)‖L2(Ω) converge to zero as the mesh is re- 94 DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖1,h l ‖ep‖L2(Ω) l 112/32 3.893e-2 – 8.297e-1 – 1.297 – 416/128 1.016e-2 1.94 4.105e-1 1.01 3.734e-1 1.78 1,600/512 2.707e-3 1.91 2.045e-1 1.01 1.293e-1 1.53 6,272/2,048 7.087e-4 1.93 1.021e-1 1.00 5.475e-2 1.24 24,832/8,192 1.813e-4 1.97 5.104e-2 1.00 2.597e-2 1.08 98,816/32,768 4.578e-5 1.99 2.552e-2 1.00 1.281e-2 1.02 Table 3.2: Example 1. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). DOFs bh/rh ‖eb‖L2(Ω) l ‖eb‖H(curl;Ω) l 56/25 4.720e-1 – 9.431e-1 – 208/81 2.358e-1 1.00 4.714e-1 1.00 800/289 1.179e-1 1.00 2.357e-1 1.00 3,136/1,089 5.893e-2 1.00 1.179e-1 1.00 12,416/4,225 2.946e-2 1.00 5.893e-2 1.00 49,408/16,641 1.473e-2 1.00 2.946e-2 1.00 Table 3.3: Example 1. Convergence of ‖eb‖L2(Ω) and ‖eb‖H(curl;Ω) . DOFs bh/rh ‖er‖L2(Ω) l ‖∇er‖L2(Ω) l 56/25 1.673e-1 – 9.391e-1 – 208/81 4.433e-2 1.92 4.824e-1 0.96 800/289 1.125e-2 1.98 2.429e-1 0.99 3,136/1,089 2.822e-3 1.99 1.216e-1 1.00 12,416/4,225 7.062e-4 2.00 6.085e-2 1.00 49,408/16,641 1.766e-4 2.00 3.043e-2 1.00 Table 3.4: Example 1. Convergence of ‖er‖L2(Ω) and ‖∇er‖L2(Ω). fined, in accordance with Theorem 3.4.2. The rate of convergence is O(h). Notice that we obtain the optimal rate in this numerical experiment, even though Theorem 3.4.3 predicts a sub-optimal rate with a loss of O(hε ). Additionally, ‖u−uh‖L2(Ω) and ‖r− rh‖L2(Ω) converge at rate O(h2) as h tends to zero, which is also optimal. In Figure 3.1 we show the convergence history of the Picard iterations for the grid sequence considered in this example. The plot depicts the number of iterations 95 against the differences between consecutive iterates corresponding to the approxi- mated vector coefficients, measured in a normalized discrete 2-norm and labelled as ‘Tolerance’ in the plot. As expected, convergence is linear and the iteration count is fairly insensitive to the size of the grid. A very similar behaviour has been observed in all of our other experiments, in 2D as well as in 3D. 1 2 3 4 5 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 Iteration To le ra nc e Grid 1 Grid 2 Grid 3 Grid 4 Grid 5 Grid 6 Figure 3.1: Example 1. Convergence history of the Picard iteration for the grid sequence defined in Tables 3.2–3.4. 3.5.2 Example 2: two-dimensional problem with a singular solution In order to verify the capability of the proposed method to capture singularities in two dimensions, we consider a problem in the L-shaped domain Ω = (−1,1)2 \ ([0,1)× (−1,0]) with ΓN = {(1,y) : y ∈ (0,1)}, ΓD = Γ\ΓN , and set ν = κ = 1, νm = 1e4. We choose the forcing terms and the boundary conditions such that the analytic solution is given by the strongest corner singularities for the underlying el- liptic operators. In polar coordinates (ρ ,φ), the hydrostatic solution components u and p are then given by u(ρ ,φ) = ρλ ((1+ λ )sin(φ)ψ(φ)+ cos(φ)ψ ′(φ)) ρλ (−(1+ λ )cos(φ)ψ(φ)+ sin(φ)ψ ′(φ)) , p(ρ ,φ) =−ρλ−1((1+ λ )2ψ ′(φ)+ ψ ′′′(φ))/(1−λ ), 96 where ψ(φ) =sin((1+ λ )φ)cos(λω)/(1+ λ )− cos((1+ λ )φ) − sin((1−λ )φ)cos(λω)/(1−λ )+ cos((1−λ )φ), λ ≈ 0.54448373678246 and ω = 32pi . The magnetic pair (b,r) is given by b(ρ ,φ) = ∇(ρ2/3 sin(2/3φ)), r(ρ ,φ) ≡ 0. For this example, we have that (u, p) ∈ H1+λ(Ω)2 × Hλ (Ω) and b∈H2/3(Ω)2. Note that straightforwardly applied nodal elements cannot correctly resolve the magnetic field. In Tables 3.5–3.6, we investigate the asymptotic rates of conver- gence of the errors in the approximations of the hydrostatic and magnetic variables. Again, we observe that the discrete solution converges to the exact one as the mesh size h approaches zero, in accordance with Theorem 3.4.2. The results show full agreement with the optimal rates for ‖u− uh‖1,h and ‖b− bh‖H(curl;Ω). For the pressure, we also see that the rate for ‖p− ph‖L2(Ω) is approaching the optimal rate, albeit more slowly. Additionally, we observe the L2-norm of rh is zero be- cause g is divergence-free, in accordance with Proposition 3.3.1. As the mesh is refined, the actual values of the L2-norm of rh slightly increase, which is likely due to the increased condition numbers of the corresponding linear systems. In Figures 3.2–3.3, we show the solution computed on the finest mesh with 24,576 elements; the total number of degrees of freedom employed in the finite element space Vh×Ch×Qh× Sh is 148,481. The results show that our solution captures the strongest corner singularities and are comparable to the results in [56]. 3.5.3 Hartmann channel flow Next, we consider Hartmann channel flow problems in two and three dimensions; cf. [35]. In these examples, we denote by Ha the Hartmann number, which is defined as Ha = √ κ ννm . 97 DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖1,h l ‖ep‖L2(Ω) l 88/24 2.159e-1 – 2.468 – 15.91 – 320/96 1.781e-1 0.28 1.880 0.39 9.328 0.77 1,216/384 1.204e-1 0.56 1.368 0.46 5.387 0.79 4,736/1,536 6.816e-1 0.82 9.588e-1 0.51 3.301 0.71 18,688/6,144 3.490e-2 0.97 6.627e-1 0.53 2.124 0.64 74,240/24,576 1.705e-2 1.03 4.559e-1 0.54 1.408 0.59 Table 3.5: Example 2. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). DOFs bh/rh ‖eb‖L2(Ω) l ‖eb‖H(curl;Ω) l ‖rh‖L2(Ω) 44/21 2.796e-1 – 2.796e-1 – 2.162e-12 160/65 1.814e-1 0.62 1.814e-1 0.62 6.188e-12 608/225 1.169e-1 0.63 1.169e-1 0.63 2.289e-11 2,368/833 7.473e-2 0.65 7.473e-2 0.65 4.260e-11 9,344/3,201 4.754e-2 0.65 4.754e-2 0.65 1.406e-10 37,120/12,545 3.013e-2 0.66 3.013e-2 0.66 3.018e-10 Table 3.6: Example 2. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω). 3.5.3.1 Example 3: two-dimensional Hartmann flow Consider the two-dimensional Hartmann flow problem, which involves a steady unidirectional flow in the channel Ω = (0,10)× (−1,1) under the influence of the constant transverse magnetic field bD = (0,1). The MHD solution then takes the form: u(x,y) = (u(y),0), p(x,y) =−Gx+ p0(y), b(x,y) = (b(y),1), r(x,y) ≡ 0. (3.36) 98 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) Figure 3.2: Example 2. Numerical approximations of (a) velocity; (b) pres- sure contours. −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) (c) Figure 3.3: Example 2. Numerical approximations of (a) magnetic field; (b) contours of the first component of the magnetic field; (c) contours of the second component of the magnetic field. We impose the following boundary conditions: u = 0 on y =±1, (pI−ν∇u)n = pNn on x = 0 and x = 10, n×b = n×bD on Γ, r = 0 on Γ, 99 where pN(x,y) = p(x,y) =−Gx− G 2 2κ ( sinh(yHa) sinh(Ha) − y )2 . The exact solution is given by (3.36) with u(y) = G νHatanh(Ha) ( 1− cosh(yHa) cosh(Ha) ) , b(y) = G κ ( sinh(yHa) sinh(Ha) − y ) , p0(y) =−G 2 2κ ( sinh(yHa) sinh(Ha) − y )2 . We note that p0(y) and −κb(y) 2 2 are the same up to an additive constant. In Tables 3.7–3.8 and Figures 3.4–3.5, we set ν = κ = 1, νm = 1e4, and G = 10. We observe that rh ≡ 0, as predicted in Proposition 3.3.1, and ‖u− uh‖1,h, ‖p− ph‖L2(Ω) and ‖b− bh‖H(curl;Ω) converge to zero at the optimal rate O(h) as the mesh is refined. Moreover, we note that the L2-norms of the errors in the approximations of u, b and p tend to zero optimally as well. In Figures 3.4–3.5 we show the solution computed on the mesh with 32,768 elements; the total number of degrees of freedom employed in the finite element space Vh×Ch×Qh×Sh is 197,633. In order to show the directions of vectors, in Figure 3.5(b) and later figures, b is normalized such that the largest magnitude of each component is 1 in the computational domain. The computed and analytical solutions of the first components in the velocity and magnetic fields are virtually indistinguishable; see Figure 3.4. 3.5.3.2 Example 4: three-dimensional Hartmann flow In this example, we consider the three-dimensional unidirectional flow in the rect- angular duct given by Ω = (0,L)× (−y0,y0)× (−z0,z0) with y0,z0 ≪ L under the influence of the constant transverse magnetic field bD = (0,1,0). We take f = g = 0 100 DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖1,h l ‖ep‖L2(Ω) l 416/128 2.028e-1 – 3.215 – 13.97 – 1,600/512 5.169e-2 1.97 1.611 1.00 6.986 1.00 6,272/2,048 1.306e-2 1.99 0.8061 1.00 3.493 1.00 24,832/8,192 3.282e-3 1.99 0.4033 1.00 1.747 1.00 98,816/32,768 8.227e-4 2.00 0.2017 1.00 0.8734 1.00 Table 3.7: Example 3. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). DOFs bh/rh ‖eb‖L2(Ω) l ‖eb‖H(curl;Ω) l ‖rh‖L2(Ω) 208/81 1.679e-4 – 2.259e-4 – 3.868e-12 800/289 8.605e-5 0.96 1.148e-4 0.98 1.746e-11 3,136/1,089 4.328e-5 0.99 5.761e-4 0.99 3.627e-11 12,416/4,225 2.167e-5 1.00 2.883e-5 1.00 9.424e-11 49,408/16,641 1.084e-5 1.00 1.442e-5 1.00 2.401e-10 Table 3.8: Example 3. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω). −1 −0.5 0 0.5 1 0 1 2 3 4 5 6 y u (y) FE Solution Analytical Solution −1 −0.5 0 0.5 1 −8 −6 −4 −2 0 2 4 6 8 x 10−5 y b(y ) FE Solution Analytical Solution (a) (b) Figure 3.4: Example 3. Slices along x = 5, −1≤ y≤ 1: (a) Velocity compo- nent u(y); (b) Magnetic component b(y). and consider solutions of the form u(x,y,z) = (u(y,z),0,0), p(x,y,z) =−Gx+ p0(y,z), b(x,y,z) = (b(y,z),1,0), r(x,y,z) ≡ 0. 101 (a) (b) Figure 3.5: Example 3. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. We enforce the boundary conditions u = 0 for y =±y0 and z =±z0, (pI−ν∇u)n = pNn for x = 0 and x = L, n×b = n×bD on Γ, r = 0 on Γ, with pN(x,y,z) = −Gx− κb(y,z) 2 2 + 10. The function b(y,z) is given by the Fourier series b(y,z) = ∞ ∑ n=0 bn(y)cos(λnz), 102 where λn = (2n+ 1)pi 2z0 , bn(y) = ν κ ( An λ 2n − p21 p1 sinh(p1y)+ Bn λ 2n − p22 p2 sinh(p2y) ) , p21,2 = λ 2n + Ha2/2±Ha √ λ 2n + Ha2/4, An = −p1(λ 2n − p22) ∆n un(y0)sinh(p2y0), Bn = p2(λ 2n − p21) ∆n un(y0)sinh(p1y0), ∆n = p2(λ 2n − p21)sinh(p1y0)cosh(p2y0)− p1(λ 2n − p22)sinh(p2y0)cosh(p1y0), un(y0) = −2G νλ 3n z0 sin(λnz0). The functions u(y,z) and p0(y,z) can be also expressed by Fourier series; for de- tails, see [34]. In fact, p0(y,z) and−κb(y,z) 2 2 are identical up to an additive constant. Note also that p(x,y,z) = pN(x,y,z). In our tests, we set L = 10, y0 = 2, z0 = 1, ν = κ = 1, νm = 1e4 and G = 0.5. In Tables 3.9–3.10, we investigate the asymptotic rates of convergence of the errors in the approximations of the hydrostatic and magnetic variables. Again, we observe that the finite element solution converges to the exact solution as the mesh size h approaches zero, in accordance with Theorem 3.4.2. We observe the results show good agreement with the optimal rates for ‖u−uh‖1,h and ‖b−bh‖H(curl;Ω). For the pressure, we also see that the rate for ‖p− ph‖L2(Ω) is approaching the optimal rate, although more slowly. Additionally, we observe the L2-norm of rh is zero because g is divergence-free, in accordance with Proposition 3.3.1. In Figures 3.6–3.7 we show the solution computed on a uniform tetrahedral mesh of 24,576 elements; this results in a total of 212,577 degrees of freedom in the finite element space Vh ×Ch ×Qh × Sh. We observe that the computed and analytical solutions are in good agreement on this relatively coarse mesh; see Figure 3.6. 103 DOFs uh/ph ‖eu‖L2(Ω) l ‖eu‖1,h l ‖ep‖L2(Ω) l 360/48 3.959e-1 – 1.829 – 30.89 – 2,592/384 1.320e-1 1.58 0.9561 0.94 8.194 1.91 19,584/3,072 3.609e-2 1.87 0.4903 0.96 2.837 1.53 152,064/24,576 9.590e-3 1.91 0.2484 0.98 1.091 1.38 Table 3.9: Example 4. Convergence of ‖eu‖L2(Ω), ‖eu‖1,h, and ‖ep‖L2(Ω). DOFs bh/rh ‖eb‖L2(Ω) l ‖eb‖H(curl;Ω) l ‖rh‖L2(Ω) 98/27 1.850e-5 – 3.219e-5 – 9.855e-12 604/125 1.565e-5 0.24 2.579e-5 0.32 1.013e-10 4,184/729 8.592e-6 0.86 1.464e-5 0.82 4.098e-10 31,024/4,913 4.411e-6 0.96 7.543e-6 0.96 1.795e-9 Table 3.10: Example 4. Convergence of ‖eb‖L2(Ω), ‖eb‖H(curl;Ω), and ‖rh‖L2(Ω). 3.5.4 Driven cavity flow Let us consider a classic test problem used in fluid dynamics, known as driven- cavity flow. It is a model of the flow in a cavity with the lid moving in one direction; cf. [31, Chapter 5.1.3] and [69]. 3.5.4.1 Example 5: two-dimensional driven cavity flow In this example, we consider the two-dimensional domain Ω = (−1,1)2 with ΓD = Γ, and set the source terms to be zero. The boundary conditions are prescribed as follows: u = 0 on x =±1 and y =−1, u = (1,0) on y = 1, n×b = n×bD on Γ, r = 0 on Γ, 104 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.05 0.1 0.15 0.2 0.25 0.3 y u (y, 0) FE Solution Analytical Solution −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −3 −2 −1 0 1 2 3 x 10−6 y b(y ,0) FE Solution Analytical Solution (a) (b) Figure 3.6: Example 4. Slices along x = 5, −2 ≤ y ≤ 2, and z = 0: (a) Velocity component u(y,0); (b) Magnetic component b(y,0). where bD = (1,0). We set ν = 1e-2, νm = 1e5, κ = 1e5, which simulate liquid metal type flows. Figures 3.8–3.9 show the solution computed on a mesh with 8,192 elements and 49,665 degrees of freedom. Figure 3.8(a) shows that the circulation created by the moving lid; Figure 3.8(b) shows the magnetic field changes direction due to the coupling effect. Figure 3.9(a) demonstrates the boundary layer formation in terms of the first component of the velocity. Streamlines for the velocity field are displayed in Figure 3.9(b). The computed solution agrees with the solution in the literature [69]. 3.5.4.2 Example 6: three-dimensional driven cavity flow The problem we consider is the three-dimensional driven cavity flow in the domain Ω = (−1,1)3 with ΓD = ∂Ω. The source terms are set to be zero. The boundary conditions are prescribed as follows: u = 0 on x =±1, y =±1 and z =−1, u = (1,0,0) on z = 1, n×b = n×bD on Γ, r = 0 on Γ, 105 (a) (b) Figure 3.7: Example 4. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. where bD = (1,0,0). We set ν = 1e-2, νm = 1e5, κ = 1e5 and obtain Figure 3.10 on a uniform tetrahedral mesh comprising 24,576 elements; this results in a total of 212,577 degrees of freedom. The flow vectors on slices demonstrate a similar behaviour to 106 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) Figure 3.8: Example 5. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. the two-dimensional scenario in Section 3.5.4.1; see Figure 3.8. 3.5.5 Example 7: two-dimensional MHD flow over a step The example we present here is another classical problem of a flow over a step under a transverse magnetic field; cf. [33]. The magnetic field tends to damp the vortex of the fluid after the step. The domain is Ω = (−0.25,0.75)×(−0.125,0.125)\(−0.25,0]×(−0.125,0], with ΓN = {(0.75,y) : y ∈ (−0.125,0.125)} and ΓD = Γ \ΓN. We set f = g = 0, 107 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) Figure 3.9: Example 5. Numerical approximations of (a) contours of the first velocity components; (b) streamlines of velocity. and choose ν = 1e-2, νm = 1e5, κ = 2.5e4. The boundary data are given by u = 0 on y =±0.125, {(x,0) : x ∈ (−0.25,0)}, u = 0 on {(0,y): y ∈ (−0.125,0)}, u = (−25.6y(y−0.125),0) on x =−0.25, (pI−ν∇u)n = pNn on x = 0.75, n×b = n×bD on Γ, r = 0 on Γ, where pN = 0 and bD = (0,1). Figures 3.11–3.12 show the solution computed on a mesh with 7,168 elements and 43,649 degrees of freedom. It is evident from Figure 3.11 that the flow field is correctly captured; the magnetic field changes directions due to the coupling effect; the pressure drops behind the step. Figure 3.12 shows the velocity field in terms of stream lines. The recirculation after the step decreases as the coupling coefficient κ increases. We observe that our numerical method reproduces this damping effect 108 (a) (b) Figure 3.10: Example 6. Numerical approximations of (a) velocity; (b) nor- malized magnetic field. without any oscillation in the numerical solution. The computed solutions agree with the solutions in the literature [22, 33]. (a) (b) (c) Figure 3.11: Example 7. Numerical approximations of (a) velocity; (b) nor- malized magnetic field; (c) pressure contours. 109 −0.05 0 0.05 −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 0 −0.05 0 0.05 −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 0 (a) (b) Figure 3.12: Example 7. Velocity flow vectors and streamlines zoomed in behind the step for (a) κ=2.5e4; (b) κ=1e5. 3.6 Conclusions We have introduced a new mixed finite element method for the numerical dis- cretization of a stationary incompressible magnetohydrodynamics problem, with divergence-conforming BDM elements and curl-conforming Nédélec elements for the velocity and magnetic fields, respectively. The approximation of the velocity field is exactly mass conservative. We have shown the well-posedness of the dis- crete formulation under a standard small data assumption, and convergence of the approximations under minimal regularity assumptions. We have proved that the energy norm error is convergent in the mesh size in general Lipschitz polyhedra, and have derived a priori error estimates. As shown in detail in Section 3.4, in the two-dimensional case there is a loss of O(hε ) in the theoretical error estimates. In the three-dimensional case our error estimates end up falling short by half a power of h for the errors in u and b, and by a full power in p and r. Nevertheless, the numerical experiments of Section 3.5 show optimal convergence in all cases. This probably indicates that the sub-optimality is a mere artifact of our technique of proof, which relies on inverse estimates to establish the continuity of the nonlinear coupling form. Furthermore, the numerical experiments indicate that the constant Cε in Theorem 3.4.3 stays bounded, even 110 though this is not guaranteed by the analysis. Altogether, the computed results are in excellent agreement with results in the literature, and the method correctly resolves the strongest magnetic singularities in non-convex domains. But there is a need to further pursue the theoretical issue of sub-optimal convergence rates. Based on the theoretical results in [70], we expect the same good performance of our discretization and solution techniques to carry over to the dynamic prob- lem, provided that the nonlinear terms are treated (semi)implicitly. This will be discussed further in detail in Chapter 4. We also mention the issue of higher order elements. Here, too, we do not expect any deviation from our current computa- tional results. In particular, we expect to see optimal convergence rates for smooth solutions. The scope of our work can be broadened in a number of additional directions. A very important issue is the investigation of efficient linear solvers for large-scale problems. In such settings iterative solvers are necessary, and this brings up the need for deriving effective and scalable preconditioners. While there are efficient solution techniques for the Navier-Stokes equations as well as for the curl-curl operator, the primary challenge is how to deal with the coupling term, especially when coupling is strong. Preliminary work on this is currently underway. Another item for future work is the derivation of a nonlinear solver that con- verges more rapidly than the Picard iteration used in our experiments. As we have pointed out in Remark 3.3.3, developing the Newton iteration for our discretization is somewhat delicate and is subject of ongoing investigation. 111 Chapter 4 Extensions In this chapter, we describe two extensions of the disretization techniques presented in Chapters 2 and 3. First, we consider the fully nonlinear time-dependent MHD problem introduced in Section 1.1. We employ the implicit Euler method for the temporal discretization, and our DG approach with exactly divergence-free veloci- ties for the spatial discretization, which as before ensures the energy-stability of the fully discrete scheme. We present a preliminary set of computational experiments. Second, we investigate an exactly divergence-free DG method for the approxi- mation of the Stokes problem with nonstandard boundary conditions. These bound- ary conditions are naturally suited for approximating the Stokes equations using a curl-curl formulation. We establish a crucial norm-equivalence property for the proposed discretization, and prove optimal convergence of our method in the bro- ken H1-norm. The theoretical results are confirmed in a series of numerical tests. 4.1 Time-dependent incompressible magnetohydrodynamics 4.1.1 Introduction A number of papers on time-dependent MHD systems can be found in the litera- ture. In [3], a class of time-stepping algorithms for the transient incompressible MHD equations has been analyzed, with focus on the long-term dissipative struc- 112 ture of the underlying dynamical system. In [46], the second-order BDF2 scheme has been applied for the time discretization of the MHD system. The Navier-Stokes equations have been decoupled from the Maxwell equations by taking the nonlin- ear terms fully explicitly. Theoretical aspects of other fully discrete schemes for time-dependent MHD problems have been presented in [70]. We also refer the reader to the monograph [35]. In this section, we implement a fully discrete scheme for the numerical approx- imation of the time-dependent incompressible MHD equations, where the implicit Euler method is employed for time stepping, and the divergence/curl-conforming finite elements presented in Chapter 3 are used for the spatial discretization. The nonlinear terms are discretized semi-implicitly. Hence, in each time step, a lin- earized (but coupled) MHD system similar to the one presented in (3.14) needs to be solved. The goal of this section is to present a preliminary set of numerical computations and to demonstrate that our DG approach is in principle applicable to transient problems. In a simply-connected Lipschitz domain Ω ∈ Rd (d = 2 or 3) and for T > 0, we consider the time-dependent MHD equations ∂tu−ν ∆u +(u ·∇)u + ∇p−κ (∇×b)×b = f in Ω× (0,T ), (4.1a) ∇ ·u = 0 in Ω× (0,T ), (4.1b) κ∂tb + κνm ∇× (∇×b)−κ ∇× (u×b) = 0 in Ω× (0,T ), (4.1c) ∇ ·b = 0 in Ω× (0,T ), (4.1d) subject to the following boundary and initial conditions u = 0 on Γ× (0,T ), (4.2a) n×b = 0 on Γ× (0,T ), (4.2b) u(0) = u0 on Ω, (4.2c) b(0) = b0 on Ω. (4.2d) 113 Note that the incompressibility constraint of the magnetic field in (4.1d) is implic- itly implied by the magnetic equation (4.1c). To see this, we take the divergence of (4.1c) and obtain ∂t(∇ ·b) = 0 ∀ t ∈ (0,T ). Therefore, the constraint (4.1d) is satisfied as long as the initial magnetic field b0 is solenoidal. 4.1.2 Weak formulation Upon setting V = H10 (Ω)d , C = H0(curl;Ω), Q = L20(Ω), the variational formulation of problem (4.1)–(4.2) is: for every t ∈ (0,T ), find u(t) ∈ V, b(t) ∈ C and p(t) ∈ Q such that (∂tu(t),v)Ω + A(u(t),v)+ O(u(t),u(t),v) (4.3a) +C(b(t),v,b(t))+B(v, p(t)) = (f(t),v)Ω, B(u(t),q) = 0, (4.3b) κ(∂tb(t),c)Ω + M(b(t),c)−C(b(t),u(t),c) = 0, (4.3c) u(0) = u0, b(0) = b0, (4.3d) for all (v,c,q) ∈ V×C×Q. The weak forms are given by A(u,v) = ∫ Ω ν ∇u : ∇vdx, O(w,u,v) = ∫ Ω (w ·∇)u ·vdx, M(b,c) = ∫ Ω κνm(∇×b) · (∇× c)dx, C(d,v,b) = ∫ Ω κ (v×d) · (∇×b)dx, B(u,q) =− ∫ Ω (∇ ·u)qdx. 114 Since ∇H10 (Ω)⊂C, we may choose c = ∇s with s∈H10 (Ω) in (4.3c) and obtain ∂t(b,∇s)Ω = 0 ∀s ∈H10 (Ω), t ∈ (0,T ). Therefore the incompressibility of the weak solution b is satisfied automatically for divergence-free b0, and there is no need to introduce a Lagrange multiplier as in the stationary case considered in Chapters 2 and 3. In addition, we have the following energy identity, cf. [70]: 1 2 ∂t(‖u‖2L2(Ω) + κ‖b‖2L2(Ω))+ ν‖∇u‖2L2(Ω) + κνm‖∇×b‖2L2(Ω) = (f,u)Ω, (4.4) which is an immediate consequence of setting v = u in (4.3a) and c = b in (4.3c), making the sum of both equations, and using the divergence constraint (4.1b). 4.1.3 Space discretization Under the same assumptions on meshes and traces as in Section 3.3.1, we look for approximations to u(t), b(t) and p(t) in the same finite element spaces as in (3.5). That is, uh(t) ∈ Vh = {v ∈ H0(div;Ω) : v|K ∈Pk(K)d , K ∈Th }, bh(t) ∈ Ch = {c ∈ H0(curl;Ω) : c|K ∈Pk−1(K)d ⊕Rk(K), K ∈Th }, ph(t) ∈ Qh = {q ∈ L20(Ω) : q|K ∈Pk−1(K), K ∈Th }, with k ≥ 1. The semi-discrete scheme then reads: for every t ∈ (0,T ), find uh(t), bh(t) and ph(t) such that (∂tuh(t),v)Ω + Ah(uh(t),v)+ Oh(uh(t),uh(t),v) (4.5a) +C(bh(t),v,bh(t))+B(v, ph(t)) = (f(t),v)Ω, B(uh(t),q) = 0, (4.5b) κ(∂tbh(t),c)Ω + M(bh(t),c)−C(bh(t),uh(t),c) = 0, (4.5c) uh(0) = u0h, bh(0) = b0h, (4.5d) 115 for all (v,c,q) ∈Vh×Ch×Qh. Here, u0h and b0h are suitable approximations of u0 and b0, respectively, with ∇ ·u0h = 0. The form Ah corresponding to the Laplacian is again the interior penalty form [4, 5]: Ah(u,v) = ∑ K∈Th ∫ K ν∇u : ∇vdx− ∑ F∈Fh ∫ F {{ν∇u}} : [[v]]ds − ∑ F∈Fh ∫ F {{ν∇v}} : [[u]]ds+ ∑ F∈Fh a0ν hF ∫ F [[u]] : [[v]]ds with a0 > 0 the interior penalty stabilization parameter. It has to be chosen larger than a threshold value which is independent of h, ν , κ and νm. For the convection term, we take the standard upwind form [62]: Oh(w,u,v) = ∑ K∈Th ∫ K (w ·∇)u ·vdx + ∑ K∈Th ∫ ∂K\Γ 1 2 (w ·nK−|w ·nK|)(ue−u) ·vds − ∫ Γ 1 2 (w ·n−|w ·n|)u ·vds. Here, ue is the trace of u taken from the exterior of K. The remaining forms are the same as in the continuous case. 4.1.4 Time discretization Let us now introduce the time step ∆t = T/N, where N is the number of time levels, and set tn = n∆t for 0 ≤ n ≤ N. We discretize in time using the implicit Euler scheme, and denote by unh, bnh, pnh the approximations to uh(tn), bh(tn) and ph(tn) in (4.5). By setting u0h = u0h, b0h = b0h, and taking the nonlinear and coupling terms semi-implicitly, the fully discrete solutions are then found by solving the following linearized (but still coupled) system on each time level 1 ≤ n ≤ N: find unh, bnh 116 and pnh such that 1 ∆t (unh,v)Ω + Ah(u n h,v)+ Oh(un−1h ,u n h,v) (4.6a) +C(bn−1h ,v,b n h)+ B(v, p n h) = (fn,v)Ω + 1 ∆t (u n−1 h ,v)Ω, B(unh,q) = 0, (4.6b) κ ∆t (b n h,c)Ω + M(bnh,c)−C(bn−1h ,unh,c) = κ ∆t (b n−1 h ,c)Ω, (4.6c) for all (v,c,q) ∈ Vh×Ch×Qh, where fn = f(tn). Closely related systems of this form have been analyzed in Chapter 3. As shown in Proposition 3.3.1, our discretization guarantees that ∇ ·unh = 0 ∀1≤ n≤ N. Therefore by choosing (v,c) = (unh,bnh), and adding equations (4.6a) and (4.6c), we obtain the discrete energy property which mimics (4.4); see also [70]: 1 2 D ( ‖unh‖2L2(Ω) + κ‖bnh‖2L2(Ω) ) + ∆t 2 ( ‖Dunh‖2L2(Ω) + κ‖Dbnh‖2L2(Ω) ) +CCν‖unh‖21,h + κνm‖∇×bnh‖2L2(Ω) ≤ (fn,unh)Ω, where ‖ · ‖1,h is the discrete H1-norm defined in (3.7) and CC is the coercivity constant in (3.9). Furthermore, D is the difference quotient operator given by Dun = un−un−1 ∆t ∀1≤ n≤ N. 4.1.5 Numerical tests We test our method on a series of computational experiments using the lowest-order BDM and Nédélec elements (k = 1) for two-dimensional systems. The computa- tion has been carried out using MATLAB with direct linear solvers. We will specify (pI− ν∇u)n on the Neumann part of the boundary denoted by ΓN , where I is the identity matrix and n is the unit outward normal on Γ. Our discretization can easily 117 be extended to Neumann boundary conditions of this type. The Dirichlet bound- ary is then denoted by ΓD. We choose the stabilization parameter a0 = 10 in all cases. The initial values u0h and b0h are set to be the BDM projection of u0 and the Nédélec projection of b0, respectively. Notice that ∇ ·u0h = 0 provided that ∇ ·u0 = 0. 4.1.5.1 Example 1: Convergence test for time discretization For T = 1, we consider a two-dimensional problem on Ω = (−1,1)2 with ΓN = {(1,y) : y ∈ (−1,1)} and ΓD = Γ\ΓN . The parameters ν , νm and κ are all set to one. We choose the right-hand side source terms and boundary conditions so that the analytical solution is given by u(x,y, t) = (x+ t2, −y+ cos t), p(x,y, t) = 10+ sin t, b(x,y, t) = (1− y+ sint, x+ et). For this example, the spatial discretization is exact; therefore only temporal er- rors are introduced by our scheme. This problem is thus nicely suited to test the accuracy of the temporal discretization. ∆t max 1≤n≤N ‖u(tn)−unh‖L2(Ω) l max1≤n≤N ‖b(tn)−b n h‖L2(Ω) l 1/4 6.060e-2 – 4.358e-1 – 1/8 3.221e-2 0.91 2.229e-1 0.97 1/16 1.657e-2 0.96 1.127e-1 0.98 1/32 8.403e-3 0.98 5.667e-2 0.99 Table 4.1: Example 1. Convergence in time of max 1≤n≤N ‖u(tn) − unh‖L2(Ω) and max 1≤n≤N ‖b(tn)−bnh‖L2(Ω). In Table 4.1, we present the asymptotic rates of convergence in time of the er- rors in the approximations of the velocity and the magnetic field. Here, we denote by u(tn) and b(tn) the exact solutions of u and b at time tn, respectively, and by l the experimental rate of convergence. The mesh size is set to h = 1/4 (since spatial errors are zero up to machine accuracy). We evaluate the max 1≤n≤N ‖·‖L2(Ω)-norm (dis- 118 crete maximum norm in time and L2-norm in space) errors of the approximations. We clearly see first-order convergence in time as expected. 4.1.5.2 Example 2: Couette channel flow Next, we consider the two-dimensional Couette channel flow problem, cf. [75], where the conducting fluid is contained in the channel Ω = (0,5)× (0,1) with the top wall y = 1 moving with velocity uD = (1,0). A uniform magnetic field bD = (0,1) is applied across the channel. The MHD solution takes the form u(x,y, t) = (u(y, t),0), b(x,y, t) = (b(y, t),1), p(x,y, t) = p(y, t). (4.7) We set f = 0 and impose the following boundary conditions u = 0 on y = 0, u = uD on y = 1, (pI−ν∇u)n = pNn on x = 0 and x = 5, n×b = n×bD on Γ, and initial conditions u(0) = b(0) = 0 on Ω. Here, pN = 1− κ2 ( −1 νmHa sinh (y−1)Ha 2 ϕ(y, t) )2 , with ϕ(y, t) = sinh yHa 2 sinh Ha2 + 2pi ∞ ∑ n=1 (−1)nn λn sin(npiy)e−λnνt , λn = Ha2 4 +(npi)2. 119 In the case ν = νm, the solution of this problem can be obtained in exact and is given by (4.7) with u(y, t) = cosh (y−1)Ha 2 ϕ(y, t), b(y, t) = −1 νmHa sinh (y−1)Ha 2 ϕ(y, t), p(y, t) = 1− κ 2 b(y, t)2. We note that p(y, t) = pN , as in Hartmann flow problems. (a) (b) (c) Figure 4.1: Example 2. Numerical approximations of velocity at time (a) t = 0.01; (b) t = 0.1; (c) t = 1. In Figures 4.1 and 4.2, we show the evolution of the computed velocity and the magnetic field, respectively. In all these tests, the parameters are chosen as ν = νm = 1, Ha = 5 and κ = ννmHa2. We observe that at an early stage t = 0.01, only the local velocity profile is changed by the effect of the moving top wall as illustrated in Figure 4.1(a), whereas at time t = 0.1, the entire field is affected, see Figure 4.1(b). The approximate velocity in Figure 4.1(c) is close to the steady state solution. In Figure 4.2, in order to show the change of the induced magnetic field b(y, t), b = (b(y, t),1) is normalized such that the largest magnitude of each component is 1 in the computational domain. The computed magnetic field at t = 1 in Figure 4.2(c) is again an almost steady profile. 120 (a) (b) (c) Figure 4.2: Example 2. Numerical approximations of normalized magnetic field at time (a) t = 0.01; (b) t = 0.1; (c) t = 1. 4.1.6 Conclusions In this section, we have presented and implemented a discrete scheme for the time- dependent MHD problem, where the time stepping is realized via the implicit Euler scheme, and the space discretization is based on the H(div)-conforming BDM el- ements and the H(curl)-conforming Nédélec elements of the first kind proposed in Chapter 3. The approximate velocities are automatically divergence-free, thus guaranteeing the stability of our discretization. We have validated the accuracy of our approach through a series of numerical tests. In each time step, a fully coupled linearized MHD problem has to be solved, which is computationally expensive. To make our approach feasible in practice, it is therefore mandatory to develop effi- cient solvers/preconditioners or further decoupling strategies as in [46, 70]. While a wide variety of efficient and robust solvers is available for both the Navier-Stokes and the Maxwell subproblems in (4.6), see [31, 41, 42, 57], the effective treatment of the coupling terms remains an open question. 4.2 An exactly divergence-free method for the Stokes equations with nonstandard boundary conditions We analyze an exactly divergence-free DG method for the Stokes equations with nonstandard boundary conditions. To incorporate the boundary conditions, a curl- curl formulation of the Stokes problem is employed. Our main result is an optimal 121 a priori error estimate in the broken H1-norm for the velocity and the L2-norm for the pressure. The theoretical results are verified in a set of numerical experiments. 4.2.1 Introduction In [19, 20], a new class of exactly divergence-free finite element methods for in- compressible fluid flow problems have been introduced. The methods are based on approximating the velocity field in a divergence-conforming finite element space and the pressure in a properly matched discontinuous space, combined with a dis- continuous Galerkin discretization of the Laplace operator appearing in the Navier- Stokes equations to enforce the full H1-continuity of the discrete velocity. The resulting numerical schemes then have the desirable property that the approximate velocity field is exactly divergence-free over the computational domain. In [20], a detailed error analysis of such methods can be found for a variety of DG methods; they have been all shown to be inf-sup stable and optimally convergent in natural norms. We also refer the reader to [21, 43, 59] for further aspects and applications. In this section, we adopt an exactly divergence-free DG method for the numer- ical solution of the Stokes equations −ν∆u + ∇p = f in Ω, (4.8a) ∇ ·u = 0 in Ω, (4.8b) subject to the nonstandard boundary conditions u ·n = 0 and (∇×u)×n = 0 on Γ. (4.9) Here, u is the velocity field of the fluid, p the pressure, f ∈ L2(Ω)d represents an external body force. The boundary conditions are normal velocity (no-penetration) and tangential vorticity conditions. We shall make the additional (smoothness) assumptions on Ω (more details can be found below). Assumption 4.2.1 The domain Ω is convex or has a boundary of class C 1,1 in the sense of [2, Notation 2.1]. The boundary conditions (4.9) lend themselves naturally to numerical methods 122 that are based on rewriting the vector Laplacian in terms of the curl-curl operator. Indeed, by using the well-known vector identity −∆u = ∇× (∇×u)−∇(∇ ·u), the Stokes equations (4.8)–(4.9) can be written in the form ν∇× (∇×u)+ ∇p = f in Ω, (4.10a) ∇ ·u = 0 in Ω, (4.10b) u ·n = 0 on Γ, (4.10c) (∇×u)×n = 0 on Γ. (4.10d) Various finite element methods incorporating nonstandard boundary conditions for the Stokes (or Navier-Stokes) equations can be found in the literature. In [58], a theoretical analysis of a class of such methods has been carried out. The finite element methods therein have been designed to handle either the mixed velocity, vorticity and pressure boundary conditions, or the velocity and traction boundary conditions. In [8], the same boundary conditions as in (4.9) have been consid- ered, and formulations which decouple the pressure from the Stokes system have been investigated. In [37], a mixed method has been proposed for the Navier- Stokes equations using a vector potential-vorticity formulation. We mention here also [38] where curl-conforming elements and nodal elements are employed for the approximations of velocity and pressure, respectively. For a detailed discussion of alternate boundary conditions and formulations of the viscous term, we refer the reader to [47]. In the spirit of [20], we develop and analyze an exactly divergence-free DG method for the numerical approximation of (4.10) on triangular/tetrahedral meshes. In particular, we use a Brezzi-Douglas-Marini (BDM) finite element space of de- gree k for the velocity approximation [11], along with a discontinuous pressure space of degree k−1. The curl-curl operator is then discretized using a standard in- terior penalty approach. The proposed DG method has several important features. First, the resulting velocity approximation is exactly divergence-free. Second, the 123 boundary conditions (4.10c)–(4.10d) can be enforced in a seemingly straightfor- ward manner, (4.10c) essentially and (4.10d) naturally. Third, in contrast to some of the aforementioned papers, the pressure is sought in L20(Ω). As usual for curl-curl formulations, our analysis requires the relatively strong regularity properties stated in Assumption 4.2.1. The velocity field u in (4.10) can then be shown to belong to H1(Ω)d . One of the main contributions of our work here is to prove an estimate for the error measured in a broken H1-norm for the ve- locity and in the L2-norm for the pressure. The critical issue arising here is proving the H1-stability of the DG discretization based on the curl-curl formulation (4.10). We do this by employing an averaging operator as introduced in [61], and by split- ting the solution into an H1-conforming part and a remainder. A similar idea can also be found in [12] for the simpler no-slip boundary condition u = 0 on Γ, in the context of hybridized LDG methods. We point out that problems of the form (4.10) also appear in mixed formula- tions of Maxwell’s equations [53]. For smooth or convex domains, we expect our DG approach to work in the Maxwell context as well. However, when Assump- tion 4.2.1 is not satisfied, it is well known that u in (4.10) might have singular components that are not in H1(Ω)d any longer. The analysis of related methods for such problems remains an open question. We mention that, for Maxwell’s equa- tions in mixed form, DG methods based on completely discontinuous spaces can be found in [53, 54]. The outline of Section 4.2 is as follows. In Section 4.2.2, we review a vari- ational formulation suitable for (4.10). In Section 4.2.3, we discretize it using a mixed Pdk -Pk−1 element pair that yields exactly solenoidal velocity approxima- tions. Section 4.2.4 is devoted to the error analysis of the proposed method. A set of numerical examples is presented in Section 4.2.5. 124 4.2.2 Weak formulation and well-posedness In this section, we present a mixed formulation of (4.10) that is based on the fol- lowing two spaces V = {v ∈ L2(Ω)d : ∇×v ∈ L2(Ω)d , ∇ ·v ∈ L2(Ω), v ·n = 0 on Γ}, Q = L20(Ω) = {q ∈ L2(Ω) : (q ,1)Ω = 0}. These spaces are endowed with the norms ‖u‖2V = ‖∇×u‖2L2(Ω) +‖∇ ·u‖2L2(Ω) and ‖p‖L2(Ω). From the Poincaré inequality in [32, Proposition 7.4], there holds ‖∇×u‖L2(Ω) +‖∇ ·u‖L2(Ω) ≥C ( ‖u‖L2(Ω) +‖∇×u‖L2(Ω) +‖∇ ·u‖L2(Ω) ) , for all u ∈ V, with a constant C > 0 only depending on Ω. It follows that ‖ · ‖V is indeed a norm on V. It further follows from [2, Theorems 2.9 and 2.17] that the space V is continuously embedded in H1(Ω)d . Hence, there are constants c1 and c2 only depending on Ω such that c1‖u‖H1(Ω) ≤ ‖u‖V ≤ c2‖u‖H1(Ω), u ∈ V. (4.11) Now, we consider the variational problem: find (u, p) ∈ V×Q such that A(u,v)+ B(v, p) = (f,v)Ω, (4.12a) B(u,q) = 0, (4.12b) for all (v,q) ∈V×Q. Here, the forms A(u,v) and B(u,q) are given by A(u,v) = ∫ Ω ν(∇×u) · (∇×v)dx, B(u,q) =− ∫ Ω (∇ ·u)qdx. To discuss the well-posedness of (4.12), we first note that the forms A and B 125 are clearly continuous: |A(u,v)| ≤ ν‖u‖V‖v‖V, u,v ∈V, (4.13a) |B(u, p)| ≤ ‖u‖V‖p‖L2(Ω), u ∈ V, p ∈ Q. (4.13b) Next, we observe that the form A is coercive on the kernel of B. That is, by defining J = {u ∈ V : B(u,q) = 0 ∀q ∈ Q}, we have A(u,u)≥ ν‖u‖2V, u ∈ J. (4.14) Finally, let us prove the following inf-sup condition for the form B. Lemma 4.2.2 There holds inf p∈Q\{0} sup v∈V\{0} B(v, p) ‖v‖V‖p‖L2(Ω) ≥C > 0, with a constant C only depending on Ω. Proof: Given p ∈Q, let φ ∈H1(Ω) be the solution of the Neumann problem −∆φ = p in Ω, ∇φ ·n = 0 on Γ. Setting v = ∇φ , we have v ∈ V, B(v, p) = ‖p‖2L2(Ω) and ‖v‖V = ‖p‖L2(Ω). These properties now readily yield the desired inf-sup condition. The theory of mixed finite element methods, see for example [11, 39], along with the stability properties in (4.13), (4.14) and Lemma 4.2.2 imply that the weak formulation (4.12) has a unique solution (u, p) ∈ V×Q and is well-posed. More- over, in view of the embedding in (4.11), the velocity field u actually belongs to H1(Ω)d . 4.2.3 Finite element approximation In this section, we propose a finite element method for (4.12) that is based on Brezzi-Douglas-Marini spaces for the approximation of the velocities, and on stan- 126 dard discontinuous finite element spaces for the pressure. 4.2.3.1 Mixed formulation For a polynomial degree k ≥ 1, we now wish to approximate the solution of (4.12) by finite element functions (uh, ph) ∈ Vh×Qh, where Vh = {u ∈ H0(div;Ω) : u|K ∈Pk(K)d , K ∈Th }, Qh = { p ∈ L20(Ω) : p|K ∈Pk−1(K), K ∈Th }. The space Vh is the divergence-conforming Brezzi-Douglas-Marini (BDM) space (see [11, Section III.3] for details); it has degrees of freedom specified for the normal components of functions along faces. Then, we introduce the finite element method: find (uh, ph) ∈ Vh ×Qh such that Ah(uh,v)+ Bh(v, ph) = (f,v)Ω, (4.15a) Bh(uh,q) = 0, (4.15b) for all (v,q) ∈Vh×Qh. The forms Ah and Bh are given by: Ah(u,v) = ∑ K∈Th ∫ K ν(∇×u) · (∇×v)dx− ∑ F∈F Ih ∫ F {{ν∇×u}} · [[v]]T ds − ∑ F∈F Ih ∫ F {{ν∇×v}} · [[u]]T ds+ ∑ F∈F Ih ∫ F a0ν hF [[u]]T · [[v]]T ds Bh(u,q) = B(u,q). The parameter a0 > 0 is an interior penalty stabilization parameter; it has to be chosen larger than a threshold value that is independent of h and ν ; see also [4]. Remark 4.2.3 Notice that the boundary conditions (4.10c) are enforced essen- tially, whereas (4.10d) is incorporated as natural boundary conditions. As in [19, 20], the approximate velocity is exactly divergence-free. 127 Proposition 4.2.4 The approximate velocity field uh ∈ Vh obtained from (4.15) is exactly divergence-free. 4.2.4 Error analysis In this section, we present the error analysis of the mixed finite element method (4.15). In particular, we derive an optimal estimate for the error measured in a broken H1- norm in the velocity, and L2-norm in the pressure. 4.2.4.1 Averaging operator As in [40, 60, 61], our main technical tool is an averaging operator. To introduce it, we define Vch = Vh∩H1(Ω)d . Proposition 4.2.5 There is an averaging operator Ih : Vh → Vch such that ∑ K∈Th ‖u− Ihu‖2L2(K) ≤C ∑ F∈F Ih hF‖[[u]]T‖2L2(F), ∑ K∈Th ‖∇(u− Ihu)‖2L2(K) ≤C ∑ F∈F Ih h−1F ‖[[u]]T‖2L2(F), with a constant C > 0 independent of the mesh size. Proof: We first denote by Wh the fully discontinuous scalar space Wh = Pk(Th), and set W ch = Wh ∩H1(Ω). Let N (Th) be a set of nodes of Th that is unisolvent in the sense that it allows us to define a Lagrange basis of W ch with respect to these nodes. For a node N ∈N (Th), the set of adjacent elements is given by δN = {K ∈Th : N is a vertex of K }. 128 For a DG function u ∈ Wh, we now define the averaging operator Ihu ∈ W ch by setting its value at each node to (Ihu)(N) = 1 card(δN) ( ∑ K∈δN u|K(N) ) , N ∈N (Th). In [61], averaging operators of this type have been introduced and analyzed which map DG functions into Wh ∩H10 (Ω), i.e., into conforming piecewise polynomial functions with zero Dirichlet boundary conditions. A slight modification of the (scaling) arguments there show that the following inequalities hold ∑ K∈Th ‖u− Ihu‖2L2(K) ≤C ∑ F∈F Ih hF‖[[u]]‖2L2(F), ∑ K∈Th ‖∇(u− Ihu)‖2L2(K) ≤C ∑ F∈F Ih h−1F ‖[[u]]‖2L2(F), where [[u]] is the usual jump of a scalar function over a face F . For more details on these averaging operators, we also refer the reader to [40, 60] and the references therein. For the vector-valued discontinuous space Wh = Pk(Th)d and its H1-conforming part Wch = Wh∩H1(Ω)d , we define Ih componentwise as Ih = (Ih, Ih, Ih). Clearly, for any u ∈Wh, we have ∑ K∈Th ‖u− Ihu‖2L2(K) ≤C ∑ F∈F Ih hF‖[[u]]‖2L2(F), ∑ K∈Th ‖∇(u− Ihu)‖2L2(K) ≤C ∑ F∈F Ih h−1F ‖[[u]]‖2L2(F). Then, let K be an element and F a face on K. For any vector u, we have the following orthogonal decomposition on F: u = (u ·nK)nK +(nK ×u)×nK, where nK is the unit outward normal on F . It follows from this identity that if u ∈ Vh is divergence-conforming, then the jump of u over F reduces to the tangential 129 jump, i.e., [[u]] = [[u]]T . This finishes the proof. 4.2.4.2 Stability We introduce the following discrete semi-norms for the velocity |u|21,h = ∑ K∈Th ‖∇u‖2L2(K) + ∑ F∈F Ih h−1F ‖[[u]]‖2L2(F), |u|2curl,h = ∑ K∈Th ‖∇×u‖2L2(K) + ∑ F∈F Ih h−1F ‖[[u]]T‖2L2(F), and define the norm |||u|||21,h = ‖u‖2L2(Ω) + |u|21,h, The bilinear forms Ah and Bh are now continuous over the finite element spaces: |Ah(u,v)| ≤C ν |u|curl,h|v|curl,h, (4.16) |Bh(u, p)| ≤C |u|1,h‖p‖L2(Ω), (4.17) for all u,v ∈ Vh, p ∈ Qh, with a constant C > 0 independent of h and ν . Next, we define the discrete kernel Jh = {u ∈Vh : Bh(u,q) = 0 ∀q ∈Qh }, and prove the following crucial coercivity result. Proposition 4.2.6 If the interior penalty parameter a0 is sufficiently large inde- pendently of h and ν , then we have Ah(u,u)≥C ν |||u|||21,h, u ∈ Jh, with a constant C > 0 independent of the mesh size h and the parameter ν . Proof: It is well-known that there is a threshold value a∗0 > 0 independent of h and 130 ν such that for a0 ≥ a∗0 there holds Ah(u,u)≥Cν |u|2curl,h, u ∈ Vh, (4.18) with C > 0 independent of h and ν ; cf. [4, 5]. Let now u ∈ Jh. As in the proof of Proposition 4.2.4, we conclude that u is exactly divergence-free. Then, by the triangle inequality, |||u|||21,h ≤C (|||Ihu|||21,h + |||u− Ihu|||21,h)=C(‖Ihu‖2H1(Ω) + |||u− Ihu|||21,h) . (4.19) Then, we claim that the H1-norm of Ihu on the right-hand side of (4.19) can be bounded by ‖Ihu‖2H1(Ω) ≤C|u|2curl,h. (4.20) To prove (4.20), we use the embedding (4.11), the triangle inequality and the fact that u is divergence-free. We obtain ‖Ihu‖2H1(Ω) ≤C ( ‖∇× (Ihu)‖2L2(Ω) +‖∇ · (Ihu)‖2L2(Ω) ) ≤C( ∑ K∈Th ( ‖∇×u‖2L2(K) + T ) . where T = ∑ K∈Th ( ‖∇× (u− Ihu)‖2L2(K) +‖∇ · (u− Ihu)‖2L2(K) ) . According to Proposition 4.2.5, we readily have T ≤C ∑ F∈F Ih h−1F ‖[[u]]T‖2L2(F), and therefore (4.20) holds. Next, we claim that |||u− Ihu|||21,h ≤C|u|2curl,h. (4.21) 131 To prove this bound, we start by noticing that ‖[[u− Ihu]]‖2L2(F) = ‖[[u]]‖2L2(F) = ‖[[u]]T‖2L2(F). Now we apply again Proposition 4.2.5 and obtain |||u− Ihu|||21,h ≤C ∑ F∈F Ih h−1F ‖[[u]]T‖2L2(F). The assertion now follows by referring to (4.18), (4.19), (4.20) and (4.21). Finally, let us address the inf-sup stability of the form Bh. We have the follow- ing result [49, Proposition 10]: Lemma 4.2.7 inf p∈Qh\{0} sup v∈Vh\{0} Bh(v, p) |||v|||1,h‖p‖L2(Ω) ≥C > 0, for a constant C > 0 independent of h and ν . Proof: The proof is analogous to that of [49, Proposition 10]. Proposition 4.2.6, Lemma 4.2.7, together with equations (4.16) and (4.17) guarantee the existence and uniqueness of the solution of (4.15). 4.2.4.3 A priori error estimates It is now straighforward to prove error estimates for the Galerkin method (4.15). In- deed, by proceeding as in [11] and using the stability properties of Section 4.2.4.2, we readily obtain the following result. Theorem 4.2.8 Suppose the solution (u, p) of (4.10) possesses the smoothness (u, p) ∈Hσ+1(Ω)d ×Hσ(Ω), (4.22) for σ > 12 . Then the following error estimates hold: ν 1 2 |||u−uh|||1,h +ν− 12 ‖p− ph‖L2(Ω)≤Chmin{σ ,k} ( ν 1 2 ‖u‖Hσ+1(Ω) +ν− 1 2 ‖p‖Hσ (Ω) ) , 132 where the constant C > 0 is independent of h and ν . 4.2.5 Numerical examples In this section we present a series of numerical experiments to confirm the opti- mal convergence rates of Theorem 4.2.8. For all the examples, the lowest-order BDM elements are employed, i.e., we take k = 1. The interior penalty stabilization parameter a0 is set to 10 for two-dimensional test problems and to 20 for three di- mensions. Whenever necessary, we enforce inhomogeneous boundary conditions u ·n = uN and (∇×u)×n = ω T in a standard fashion. In this case, the load vector on the right-hand side has to be adjusted accordingly. The resulting linear systems are solved using the preconditioned MINRES iterative method, with tolerance set to 1e-6 (in the vector 2-norm of the relative residual). As a preconditioner, we have chosen the approach of [41] for saddle-point systems with highly singular (1,1)-blocks. 4.2.5.1 Example 1: a two-dimensional problem We consider the following problem on Ω = (−1,1)2. We set ν = 1, and choose the source term f and the boundary conditions uN and ω T so that the analytical solution is given by u(x,y) = (−ex(ycos y+ siny),exysin y), p(x,y) = 2ex siny. DOFs uh/ph |eu|1,h l ‖eu‖L2(Ω) l ‖ep‖L2(Ω) l 416/128 8.5664e-1 – 8.7865e-2 – 2.1706 – 1,600/512 4.2240e-1 1.02 2.4900e-2 1.82 1.1633 0.90 6,272/2,048 2.0988e-1 1.01 6.4988e-3 1.94 5.9879e-1 0.96 24,832/8,192 1.0476e-1 1.00 1.6488e-3 1.98 3.0300e-1 0.98 98,816/32,768 5.2369e-2 1.00 4.1444e-4 1.99 1.5229e-1 0.99 Table 4.2: Example 1. Convergence of |eu|1,h, ‖eu‖L2(Ω), and ‖ep‖L2(Ω). In Table 4.2, we investigate the asymptotic rates of convergence of the errors in the approximations of the velocity and pressure. Here, l denotes the experimental 133 convergence rate. We observe that |u−uh|1,h and ‖p− ph‖L2(Ω) converge to zero with first order in the mesh size, while ‖u−uh‖L2(Ω) converges with second order. 4.2.5.2 Example 2: a three-dimensional problem Next, we consider a three-dimensional problem in the cube Ω = (−1,1)3. We set ν = 1, and choose the source term f and the boundary conditions uN and ω T so that the analytical solution is given by u(x,y,z) = (y2,x2− y,z+ y2), p(x,y,z) = z. DOFs uh/ph |eu|1,h l ‖eu‖L2(Ω) l ‖ep‖L2(Ω) l 360/48 2.7139 – 2.0165e-1 – 1.5716 – 2,592/384 1.3621 0.99 5.4302e-2 1.89 6.0147e-1 1.39 19,584/3,072 6.8314e-1 1.00 1.4425e-2 1.91 2.4423e-1 1.30 152,064/24,576 3.4218e-1 1.00 3.7167e-3 1.96 1.0525e-1 1.21 Table 4.3: Example 2. Convergence of |eu|1,h, ‖eu‖L2(Ω), and ‖ep‖L2(Ω). Table 4.3 shows the asymptotic rates of convergence of the errors measured in the appropriate norms. Again, we see that |u−uh|1,h and ‖u−uh‖L2(Ω) converge with order O(h) and O(h2), respectively, as h tends to zero. The L2-norm of the error in the pressure converges at a rate slightly higher than O(h). 4.2.6 Conclusions We have proposed and analyzed a H(div)-conforming discretization for a curl-curl formulation of the Stokes equations with nonstandard boundary conditions. The tangential continuity of the approximate velocity field is enforced through a DG approach. We have proved optimal convergence rates for problems with smooth solutions. The theoretical results have been verified in two- and three-dimensional test problems. 134 Chapter 5 Conclusions and future work 5.1 Conclusions In this thesis, we have developed, analyzed and numerically tested mixed discon- tinuous Galerkin finite element methods for the numerical approximation of in- compressible magnetohydrodynamics problems. In Chapter 2, we have presented the first interior penalty DG method for a linearized stationary incompressible MHD problem, whereby all the variables are approximated in discontinuous finite element spaces. We have derived a priori er- ror estimates for the energy norm error in general (possibly non-convex) polyhedral domains, and have computationally verified them in a set of numerical examples. The theoretical convergence rates are optimal in the approximation of the veloc- ity u, the pressure p and the magnetic field b, but suboptimal by one order in the approximation of the multiplier r related to the divergence constraint of the mag- netic field. This is due to the fact that we have used polynomials of one degree higher to approximate r, as in the second family of Nédélec elements [68]. On the other hand, this choice leads to optimal L2-approximations of the magnetic field. Numerically, however, optimal asymptotic rates of convergence have been observed for all variables. In principle, it is possible to extend the fully discontinuous approach of Chap- ter 2 to the nonlinear setting, following the ideas presented in [19]. However, in Chapter 3, we have chosen to modify the fully discontinuous discretization and 135 have introduced a new mixed DG method for the numerical discretization of a fully nonlinear stationary incompressible magnetohydrodynamics problem. This proposed approach is based on employing divergence-conforming BDM elements and curl-conforming Nédélec elements for the velocity and magnetic fields, re- spectively, while DG techniques are used to enforce the tangential continuity of the discrete velocity fields. There have been two main reasons for these modifi- cations. First, the approximation of the velocity field is now exactly mass conser- vative, which ensures the energy-stability of the method in a straightforward man- ner. Second, the discrete magnetic field can be split naturally into a (discretely) divergence-free part and a (discrete) gradient, thus exactly mimicking desirable properties of the underlying PDE system. We have shown convergence in general Lipschitz polyhedra of the approximations for small data under a minimal regular- ity assumption. In addition, we have derived a priori error estimates. As shown in detail in Section 3.4, in the two-dimensional case there is a loss of order O(hε ) in the theoretical error estimates. In the three-dimensional case our error estimates end up falling short by half a power of h. We have successfully tested the method for a number of benchmark problems available in the literature: Hartmann flow, driven cavity, and flow over a step. The numerical experiments show optimal con- vergence in all cases, and the method correctly resolves the strongest magnetic singularities in non-convex domains. In the first part of Chapter 4 (Section 4.1), we have extended our discretization techniques to the time-dependent MHD problem, and have implemented a fully dis- crete scheme. We have computationally tested the method in a series of numerical experiments, and have verified convergence of the finite element approximations in time. In the second part of Chapter 4 (Section 4.2), we have analyzed an exactly divergence-free DG method for a curl-curl formulation of the Stokes equations with nonstandard boundary conditions. By establishing a crucial norm-equivalence property, we have shown optimal convergence in the broken H1-norm error of the velocity. Our theoretical results have been computationally validated in a set of numerical tests for problems with smooth solutions in two and three dimensions. 136 5.2 Future work In this section, we identify a number of items that remain open for future research. • One of the most important challenges of the approaches presented in this thesis is that solving the linear systems of equations is very computationally expensive. For large-scale problems, iterative linear solvers become manda- tory, which brings up the need for developing effective and scalable precon- ditioners. Indeed, a wide variety of efficient and robust saddle point solvers are available in the literature for both the Navier-Stokes and the Maxwell op- erators; see [31, 41, 42, 57]. However, the primary challenge remains how to deal with the coupling term. A few preliminary ideas for preconditioning the fully coupled linear systems have been numerically explored in [63, Section 4.3], with mixed results, and there is room for a lot more research devoted to this aspect. • For fully nonlinear systems, solvers which converge faster than the Picard it- eration are preferred. For example, it would be desirable to develop Newton- type methods which are super-linearly convergent, but do not require the ex- plicit calculation of the Jacobian matrix. However, as we have pointed out in Remark 3.3.3, developing the Newton iteration for our discretization is somewhat delicate due to the presence of upwind convection and remains an open issue. • Another open question is the optimality of the error estimates derived in Sec- tion 3.4. Although only suboptimal error estimates have been proved so far, numerical experiments show optimal convergence in all cases. This probably indicates that the sub-optimality is a mere artifact of our technique of proof, which relies on inverse estimates to establish the continuity of the nonlinear coupling form. In order to overcome this difficulty, new Sobolev-type em- bedding results and approximation properties of (discretely) divergence-free function spaces have to be established. • For the time-dependent MHD problem in Section 4.1, we have chosen to solve a fully coupled system within each time step. Future work could in- 137 volve comparing and analyzing different decoupling mechanisms for time- stepping. Higher-order time discretization schemes should also be consid- ered. • In this thesis, we have designed methods with exactly divergence-free veloc- ity approximations. We believe it should also be possible to design methods with exactly divergence-free magnetic approximations. Indeed, for prob- lems with smooth magnetic fields, this could be done straightforwardly in the spirit of the method presented in Section 4.2. However, the approach there cannot be adopted to the more interesting case of problems with singu- lar magnetic fields. In fact, the design and analysis of exactly divergence-free methods in this case remains an issue open for future research. 138 Bibliography [1] P. Amestoy, I. Duff, and J.-Y. L’Excellent. Multifrontal parallel distributed symmetric and unsymmetric solvers. Comput. Methods Appl. Mech. Engrg., 184:501–520, 2000. → pages 51 [2] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault. Vector potentials in three-dimensional non-smooth domains. Math. Meth. Appl. Sci., 21: 823–864, 1998. → pages 11, 31, 78, 84, 122, 125 [3] F. Armero and J. Simo. Long-term dissipativity of time-stepping algorithms for an abstract evolution equation with applications to the incompressible MHD and Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 131:41–90, 1996. → pages 2, 4, 5, 6, 9, 20, 22, 67, 68, 112 [4] D. Arnold. An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal., 19:742–760, 1982. → pages 73, 116, 127, 131 [5] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39:1749–1779, 2001. → pages 29, 32, 73, 76, 116, 131 [6] S. Balay, K. Buschelman, W. Gropp, D. Kaushik, M. Knepley, L. McInnes, B. Smith, and H. Zhang. PETSc Web page, 2001. http://www.mcs.anl.gov/petsc. → pages 51 [7] A. Bonito and J.-L. Guermond. Approximation of the eigenvalue problem for the time harmonic Maxwell system by continuous Lagrange finite elements. Math. Comp., posted on February 4, 2011, PII S 0025-5718(2011)02464-6 (to appear in print). → pages 7, 68 [8] J. Bramble and P. Lee. On variational formulations for the Stokes equations with nonstandard boundary conditions. RAIRO Modél. Math. Anal. Numér., 28(7):903–919, 1994. → pages 17, 123 139 [9] S. Brenner. Poincaré-Friedrichs inequalities for piecewise H1-functions. SIAM J. Numer. Anal., 41:306–324, 2003. → pages 45 [10] S. Brenner and L. Scott. The Mathematical Theory of Finite Element Methods. Springer-Verlag, 1994. → pages 83 [11] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. In Springer Series in Computational Mathematics, volume 15. Springer-Verlag, New York, 1991. → pages 18, 29, 36, 39, 41, 45, 48, 69, 72, 76, 93, 123, 126, 127, 132 [12] J. Carrero, B. Cockburn, and D. Schötzau. Hybridized globally divergence-free LDG methods. Part I: The Stokes problem. Math. Comp., 75:533–563, 2006. → pages 124 [13] P. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, Amsterdam, 1978. → pages 93 [14] B. Cockburn and C.-W. Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. J. Sci. Comput., 16:173–261, 2001. → pages 8, 20, 68 [15] B. Cockburn, G. Karniadakis, and C.-W. Shu. The development of discontinuous Galerkin methods. In B. Cockburn, G. Karniadakis, and C.-W. Shu, editors, Discontinuous Galerkin Methods: Theory, Computation and Applications, volume 11 of Lect. Notes Comput. Sci. Engrg., pages 3–50. Springer-Verlag, 2000. → pages 8 [16] B. Cockburn, G. Karniadakis, and C.-W. Shu, editors. Discontinuous Galerkin Methods. Theory, Computation and Applications, volume 11 of Lect. Notes Comput. Sci. Eng. Springer-Verlag, 2000. → pages 8, 20, 68 [17] B. Cockburn, G. Kanschat, D. Schötzau, and C. Schwab. Local discontinuous Galerkin methods for the Stokes system. SIAM J. Numer. Anal., 40:319–343, 2002. → pages 10 [18] B. Cockburn, G. Kanschat, and D. Schötzau. Local discontinuous Galerkin methods for the Oseen equations. Math. Comp., 73:569–593, 2004. → pages 10, 20, 22, 26, 29, 68, 76 [19] B. Cockburn, G. Kanschat, and D. Schötzau. A locally conservative LDG method for the incompressible Navier-Stokes equations. Math. Comp., 74: 1067–1095, 2005. → pages 10, 13, 20, 26, 68, 81, 82, 122, 127, 135 140 [20] B. Cockburn, G. Kanschat, and D. Schötzau. A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations. J. Sci. Comput., 31:61–73, 2007. → pages 9, 13, 17, 65, 69, 73, 74, 122, 123, 127 [21] B. Cockburn, G. Kanschat, and D. Schötzau. An equal-order DG method for the incompressible Navier-Stokes equations. J. Sci. Comput., 40:188–210, 2009. → pages 13, 122 [22] R. Codina and N. Silva. Stabilized finite element approximation of the stationary magneto-hydrodynamics equations. Comp. Mech., 38:344–355, 2006. → pages 109 [23] M. Costabel and M. Dauge. Singularities of electromagnetic fields in polyhedral domains. Arch. Rational Mech. Anal., 151:221–276, 2000. → pages 6, 11, 68, 78 [24] M. Costabel and M. Dauge. Weighted regularization of Maxwell equations in polyhedral domains. Numer. Math., 93:239–277, 2002. → pages 6, 7, 11, 20, 68, 78 [25] T. G. Cowling. Magnetohydrodynamics. Adam Hilger, England, 1976. → pages 2 [26] M. Dauge. Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. Part I: Linearized equations. SIAM J. Math. Anal., 20:74–97, 1989. → pages 11, 31, 78 [27] P. A. Davidson. An Introduction to Magnetohydrodynamics. Cambridge University Press, 2001. → pages 1 [28] C. Dawson. Foreword for the special issue on discontinuous Galerkin methods. Comput. Methods Appl. Mech. Engrg., 195:3183, 2006. → pages 8, 68 [29] L. Demkowicz and L. Vardapetyan. Modeling of electromagnetic absorption/scattering problems using hp–adaptive finite elements. Comput. Methods Appl. Mech. Engrg., 152:103–124, 1998. → pages 68 [30] M. Discacciati. Numerical approximation of a steady MHD problem. In U. Langer, M. Discacciati, D. E. Keys, O. B. Widlund, and W. Zulehner, editors, Domain Decomposition Methods in Science and Engineering XVII, volume 60 of Lect. Notes Comput. Sci. Eng., pages 313–320. Springer-Verlag, 2008. → pages 9, 68 141 [31] H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics. Numerical Mathematics and Scientific Computation. Oxford University Press, 2006. → pages 104, 121, 137 [32] P. Fernandes and G. Gilardi. Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions. Math. Models Methods Appl. Sci., 7:957–991, 1997. → pages 125 [33] J.-F. Gerbeau. A stabilized finite element method for the incompressible magnetohydrodynamic equations. Numer. Math., 87:83–111, 2000. → pages 2, 6, 20, 23, 67, 68, 107, 109 [34] J.-F. Gerbeau. Problèmes mathématiques et numériques posés par la modélisation de l’électrolyse de l’aluminium. PhD thesis, ´Ecole Nationale des Ponts et Chaussées, 1998. → pages 63, 78, 103 [35] J.-F. Gerbeau, C. L. Bris, and T. Lelièvre. Mathematical Methods for the Magnetohydrodynamics of Liquid Metals. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2006. → pages 1, 2, 5, 23, 58, 67, 78, 97, 113 [36] S. Giani, E. Hall, and P. Houston. AptoFEM Users Manual. Technical report, University of Nottingham, (in preparation). → pages 51 [37] V. Girault. Incompressible finite element methods for Navier-Stokes equations with nonstandard boundary conditions in R3. Math. Comp., 51: 55–74, 1988. → pages 17, 123 [38] V. Girault. Curl-conforming finite element methods for Navier-Stokes equations with non-standard boundary conditions in R3. In J. G. Heywood, K. Masuda, R. Rautmann, and V. A. Solonnikov, editors, The Navier-Stokes Equations: Theory and Numerical Methods, volume 1431 of Lect. Notes Math., pages 201–218. Springer-Verlag, 1990. → pages 123 [39] V. Girault and P. Raviart. Finite Element Methods for Navier-Stokes Equations, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, New York, 1986. → pages 80, 126 [40] V. Girault, B. Rivière, and M. Wheeler. A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier-Stokes problems. Math. Comp., 74:53–84, 2005. → pages 84, 128, 129 142 [41] C. Greif and D. Schötzau. Preconditioners for saddle point linear systems with highly singular (1,1) blocks. ETNA, 22:114–121, 2006. → pages 121, 133, 137 [42] C. Greif and D. Schötzau. Preconditioners for the discretized time-harmonic maxwell equations in mixed form. Numer. Linear Algebra Appl., 14: 281–297, 2007. → pages 121, 137 [43] C. Greif, D. Li, D. Schötzau, and X. Wei. A mixed finite element method with exactly divergence-free velocities for incompressible magnetohydrodynamics. Comput. Methods Appl. Mech. Engrg., 199: 2840–2855, 2010. → pages 9, 16, 122 [44] J.-L. Guermond and P. Minev. Mixed finite element approximation of an MHD problem involving conducting and insulating regions: the 2D case. Math. Model. Numer. Anal., 36:517–536, 2002. → pages 9 [45] J.-L. Guermond and P. Minev. Mixed finite element approximation of an MHD problem involving conducting and insulating regions: the 3D case. Num. Meth. Part. Diff. Eqs., 19:709–731, 2003. → pages 9, 20, 68 [46] J.-L. Guermond, R. Laguerre, J. Léorat, and C. Nore. Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method. J. Comput. Phys., 228:2739–2757, 2009. → pages 9, 68, 113, 121 [47] M. Gunzburger. Finite Element Methods for Viscous Incompressible Flows: A Guide to Theory, Practice, and Algorithms. Academic Press, San Diego, 1989. → pages 123 [48] M. Gunzburger, A. Meir, and J. Peterson. On the existence and uniqueness and finite element approximation of solutions of the equations of stationary incompressible magnetohydrodynamics. Math. Comp., 56:523–563, 1991. → pages 2, 6, 20, 23, 65, 67, 68, 78 [49] P. Hansbo and M. Larson. Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’s method. Comput. Methods Appl. Mech. Engrg., 191:1895–1908, 2002. → pages 10, 26, 29, 76, 132 [50] U. Hasler, A. Schneebeli, and D. Schötzau. Mixed finite element approximation of incompressible MHD problems based on weighted regularization. Appl. Numer. Math., 51:19–45, 2004. → pages 7, 68 143 [51] J. Hesthaven and T. Warburton. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, volume 54 of Texts in Applied Mathematics. Springer, 2008. → pages 8, 20 [52] R. Hiptmair. Finite elements in computational electromagnetism. Acta Numerica, 11:237–339, 2002. → pages 8, 68, 76, 83 [53] P. Houston, I. Perugia, and D. Schötzau. Mixed discontinuous Galerkin approximation of the Maxwell operator. SIAM J. Numer. Anal., 42:434–459, 2004. → pages 10, 32, 34, 51, 124 [54] P. Houston, I. Perugia, and D. Schötzau. Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation. J. Sci. Comput., 22:315–346, 2005. → pages 10, 20, 26, 29, 30, 31, 32, 34, 51, 52, 57, 68, 76, 124 [55] P. Houston, I. Perugia, and D. Schötzau. Mixed discontinuous Galerkin approximation of the Maxwell operator: The indefinite case. Math. Model. Numer. Anal., 39:727–754, 2005. → pages 10, 20, 26 [56] P. Houston, D. Schötzau, and X. Wei. A mixed DG method for linearized incompressible magnetohydrodynamics. J. Sci. Comput., 40:281–314, 2009. → pages 9, 12, 68, 69, 80, 97 [57] Q. Hu and J. Zou. Substructuring preconditioners for saddle-point problems arising from Maxwells equations in three dimensions. Math. Comp., 73: 35–61, 2004. → pages 121, 137 [58] T. Hughes and L. Franca. A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: symmetric formulations that converge for all velocity/pressure spaces. Comput. Methods Appl. Mech. Engrg., 65: 85–96, 1987. → pages 123 [59] G. Kanschat and D. Schötzau. Energy norm a posteriori error estimation for divergence-free discontinuous Galerkin approximations of the Navier-Stokes equations. Int. J. Numer. Meth. Fluids, 57:1093–1113, 2008. → pages 13, 122 [60] O. Karakashian and W. Jureidini. A nonconforming finite element method for the stationary Navier-Stokes equations. SIAM J. Numer. Anal., 35: 93–120, 1998. → pages 82, 128, 129 144 [61] O. Karakashian and F. Pascal. A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems. SIAM J. Numer. Anal., 41:2374–2399, 2003. → pages 84, 124, 128, 129 [62] P. LeSaint and P. Raviart. On a finite element method for solving the neutron transport equation. In C. de Boor, editor, Mathematical Aspects of Finite Elements in Partial Differential Equations, pages 89–123. Academic Press, New York, 1974. → pages 73, 116 [63] D. Li. Numerical Solution of the Time-harmonic Maxwell Equations and Incompressible Magnetohydrodynamics Problems. PhD thesis, The University of British Columbia, 2010. → pages 137 [64] F. Li and C.-W. Shu. Locally divergence-free discontinuous Galerkin methods for MHD equations. J. Sci. Comput., 22-23:413–442, 2005. → pages 9, 20 [65] H. K. Messerle. Magnetohydrodynamic Electrical Power Generation. John Wiley & Sons, Chichester, 1995. → pages 1 [66] P. Monk. Finite Element Methods for Maxwell’s Equations. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2003. → pages 8, 13, 20, 30, 32, 50, 51, 68, 93 [67] U. Müller and L. Bühler. Magnetofluiddynamics in Channels and Containers. Springer-Verlag, Berlin, 2001. → pages 1, 2 [68] J. C. Nédélec. Mixed finite elements in R3. Numer. Math., 50:57–81, 1980. → pages 9, 30, 32, 50, 68, 69, 72, 135 [69] A. I. Nesliturk, S. H. Aydin, and M. Tezer-Sezgin. Two-level finite element method with a stabilizing subgrid for the incompressible MHD equations. Int. J. Numer. Meth. Fluids, 58:551–572, 2008. → pages 104, 105 [70] A. Prohl. Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamic system. Math. Model. Numer. Anal., 42:1065–1087, 2008. → pages 9, 16, 111, 113, 115, 117, 121 [71] P. H. Roberts. An Introduction to Magnetohydrodynamics. Longmans, London, 1967. → pages 5 [72] O. Schenk and K. Gärtner. Solving unsymmetric sparse systems of linear equations with PARDISO. J. Fut. Gen. Comp. Sys., 20(3):475–487, 2004. → pages 51 145 [73] D. Schötzau. Mixed finite element methods for incompressible magneto-hydrodynamics. Numer. Math., 96:771–800, 2004. → pages 7, 20, 21, 23, 24, 65, 67, 68, 71, 77, 84 [74] D. Schötzau, C. Schwab, and A. Toselli. Mixed hp-DGFEM for incompressible flows. SIAM J. Numer. Anal., 40:2171–2194, 2003. → pages 10, 20, 26, 32, 34, 39, 45, 81, 86 [75] G. W. Schwab and A. Sherman. Engineering Magnetohydrodynamics. McGraw-Hill, New York, 1965. → pages 119 [76] A. Soward, editor. Magnetohydrodynamics and the Earth’s Core, Selected Works of Paul Roberts. The Fluid Mechanics of Astrophysics and Geophysics. Taylor & Francis, London and New York, 2003. → pages 1 [77] T. Warburton and G. Karniadakis. A discontinuous Galerkin method for the viscous MHD equations. J. Comput. Phys., 152:608–641, 1999. → pages 9, 20 146
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mixed discontinuous Galerkin finite element methods...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mixed discontinuous Galerkin finite element methods for incompressible magnetohydrodynamics Wei, Xiaoxi 2011
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Mixed discontinuous Galerkin finite element methods for incompressible magnetohydrodynamics |
Creator |
Wei, Xiaoxi |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | We develop and analyze mixed discontinuous Galerkin finite element methods for the numerical approximation of incompressible magnetohydrodynamics problems. Incompressible magnetohydrodynamics is the area of physics that is concerned with the behaviour of electrically conducting, resistive, incompressible and viscous fluids in the presence of electromagnetic fields. It is modelled by a system of nonlinear partial differential equations, which couples the Navier-Stokes equations with the Maxwell equations. In the first part of this thesis, we introduce an interior penalty discontinuous Galerkin method for the numerical approximation of a linearized incompressible magnetohydrodynamics problem. The fluid unknowns are discretized with the discontinuous ℘k-℘k-1 element pair, whereas the magnetic variables are approximated by discontinuous ℘k-℘k+1 elements. Under minimal regularity assumptions, we carry out a complete a priori error analysis and prove that the energy norm error is optimally convergent in the mesh size in general polyhedral domains, thus guaranteeing the numerical resolution of the strongest magnetic singularities in non-convex domains. In the second part of this thesis, we propose and analyze a new mixed discontinuous Galerkin finite element method for the approximation of a fully nonlinear incompressible magnetohydrodynamics model. The velocity field is now discretized by divergence-conforming Brezzi-Douglas-Marini elements, and the magnetic field by curl-conforming Nédélec elements. In addition to correctly capturing magnetic singularities, the method yields exactly divergence-free velocity approximations, and is thus energy-stable. We show that the energy norm error is convergent in the mesh size in possibly non-convex polyhedra, and derive slightly suboptimal a priori error estimates under minimal regularity and small data assumptions. Finally, in the third part of this thesis, we present two extensions of our discretization techniques to time-dependent incompressible magnetohydrodynamics problems and to Stokes problems with nonstandard boundary conditions. All our discretizations and theoretical results are computationally validated through comprehensive sets of numerical experiments. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-05-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0080512 |
URI | http://hdl.handle.net/2429/34789 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_wei_xiaoxi.pdf [ 1.81MB ]
- Metadata
- JSON: 24-1.0080512.json
- JSON-LD: 24-1.0080512-ld.json
- RDF/XML (Pretty): 24-1.0080512-rdf.xml
- RDF/JSON: 24-1.0080512-rdf.json
- Turtle: 24-1.0080512-turtle.txt
- N-Triples: 24-1.0080512-rdf-ntriples.txt
- Original Record: 24-1.0080512-source.json
- Full Text
- 24-1.0080512-fulltext.txt
- Citation
- 24-1.0080512.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0080512/manifest