A Monolithic Finite Element Formulation for Magnetohydrodynamics Involving a Compressible Fluid

This work develops a new monolithic finite-element-based strategy for magnetohydrodynamics (MHD) involving a compressible fluid based on a continuous velocity–pressure formulation. The entire formulation is within a nodal finite element framework, and is directly in terms of physical variables. The exact linearization of the variational formulation ensures a quadratic rate of convergence in the vicinity of the solution. Both steady-state and transient formulations are presented for two- and three-dimensional flows. Several benchmark problems are presented, and comparisons are carried out against analytical solutions, experimental data, or against other numerical schemes for MHD. We show a good coarse-mesh accuracy and robustness of the proposed strategy, even at high Hartmann numbers.

The salient features of this work are as follows: • Most of the MHD strategies treat the fluid as incompressible [28,29]. Since compressibility effects play a key role in applications such as magneto-gas dynamics, in this work, we focus on developing a MHD strategy for compressible fluids. Moreover, we show that the strategy works very well, even when the fluid is almost incompressible (with a stiff equation of state); thus, the same strategy can be used even when the fluid is almost incompressible; • We develop a monolithic strategy based on the fully coupled equations for the fluid and the magnetic fields. The inclusion of the coupling terms yields a faster rate of convergence compared to a staggered strategy, where the fluid and the magnetic field variables are solved for in a sequential manner; • The formulation is based on primitive flow variables, such as velocity, density, temperature, pressure, and the magnetic field, which makes the implementation simple. The treatment of the boundary conditions is also very straightforward, as these are generally specified in terms of the primitive variables or their duals; • In contrast to the work in references [1,2,[4][5][6][7][8] that use a stabilized formulation, we use a stable formulation based on an appropriate choice of interpolations for the various field variables. We use higher order interpolation functions for the velocity and magnetic fields, as compared to the pressure, density, and temperature field variables. This ensures the satisfaction of the inf-sup stability conditions. No stabilizing terms need to be used, and no parameters need to be adjusted in the proposed formulation; • The governing equations for the fluid, namely the continuity, Navier-Stokes, energy, and state equations are considered in their entirety without any approximations, such as the Boussinesq approximation; • Joule heating effects have been included in the formulation.

Governing Equations for MHD Involving a Compressible Fluid
We first present the equations for the magnetic fields, and then the coupled equations to be solved on the entire domain.

Magnetic Field Equations
The Maxwell equations for electromagnetism under MHD assumptions reduce to the following set of equations [29]: where H and E represent the magnetic and electric fields, respectively, µ m is the magnetic permeability, and j is the current density. The above equations are supplemented by Ohm's law for a conducting fluid, which is given by where σ e and u denote the electrical conductivity and the velocity of the fluid, respectively. Combining Equations (1) and (2), we obtain the governing differential equation for H as where Ω represents the fluid domain.

Coupled Equations for Compressible MHD
Let Ω and Γ denote the domain and its boundary, Γ t and Γ u denote the parts of boundaries where tractions and velocity are prescribed, respectively, Γ θ and Γ q represent portions of Γ where temperature and normal heat flux are prescribed, respectively, and Γ H and Γ E represent portions of the boundary where H × n and (∇ × H) × n are prescribed, respectively.
We are interested in finding an approximate solution to the following initial-boundary value problem: Determine the velocity u, density ρ, pressure p, temperature θ, magnetic field H, stresses τ, rate of deformation D, and heat flux q on the domain Ω such that where b is the body force per unit mass, σ is the viscous stress, and C v is the specific heat at constant volume. For a Newtonian fluid, τ = −p(ρ, θ)I + σ, where σ = λ(tr D)I + 2νD, with λ and µ denoting the viscosity coefficients. Similarly, we have q = −κ∇θ, where κ is the thermal conductivity. The last term in Equation (4c) represents the additional body force, namely the Lorentz force µ m (j × H), that acts on the fluid. A term j · E, which represents Joule heating [30], gets added to the first law of thermodynamics [31] for a fluid due to magnetic effects. After subtracting the dot product of the momentum equation given by Equation (4c) with the velocity field u from the first law of thermodynamics [32], and using the properties of the scalar triple product and Equation (2), we have which is the last term in Equation (4e). The above governing equations are to be solved subject to appropriate initial and boundary conditions.

Variational Formulation
Denoting the variations of H, u, ρ, θ, and p with a subscript δ, the variational statements corresponding to Equations (4) can be written as where C c and D c denote the material constitutive tensor and rate of deformation tensor expressed in an engineering form as , q n is the normal heat flux prescribed on Γ q , andt is the prescribed traction on Γ t . Note that the weak implementation of the state equation given by Equation (5e) is critical in ensuring that a lower-order interpolation for the pressure can be used in order to satisfy the inf-sup conditions, which would otherwise be dictated by the interpolations being used for the density and temperature. Equation (5a) has been derived using a penalty term of the type (∇ · H) 2 , as shown in [29]. In addition, the term has been excluded from Equation (5a), since, in the problems that we consider here, either H × n is specified on the boundary, or the surface is purely conducting.

Time Stepping Strategy
The time discretization on the time interval [t n , t n+1 ] is carried out by using the following general trapezoidal rule: where β α = (1 − α)β n + αβ n+1 for any field variable β.

Linearization
A Newton-Raphson strategy is developed by linearizing the above equations. Let (u k , ρ k , θ k , p k , H k ) and (u k+1 , ρ k+1 , θ k+1 , p k+1 , H k+1 ) represent the values of the velocity, density, temperature, pressure, and magnetic field, respectively, at the k'th and k + 1'th iterations at time step n + 1. Let (u ∆ , ρ ∆ , θ ∆ , p ∆ , H ∆ ) denote the increments in the respective field variables. The chain and product rules are applied to obtain the linearized form of Equation (6) as follows: where k f and C c denote the derivatives of k f and C c with respect to θ.

Finite Element Formulation
Let the magnetic field, velocity, density, temperature, pressure fields, and their variations and increments be interpolated as In the case of the 27-noded hexahedral element, the shape functions N for u and H are the standard triquadratic Lagrange interpolation functions, while the density, temperature, and pressure fields are interpolated using a trilinear continuous interpolation in order to obtain a stable discretization. We thus have The various cross-product terms occurring in the formulation are discretized as follows: Using the arbitrariness of the variations, the matrix form of the incremental Equation (7) can be written as where The corresponding steady-state solution can be obtained directly without time stepping by solving

Two-Dimensional Formulation
In the case of two-dimensional flows in the x-y plane, only (∇ × H) z is non-zero, and is given by The velocity, magnetic field, rate of deformation tensor, and constitutive tensor are given by whereas the other matrices are as follows: The cross-product terms in the two-dimensional formulation are as follows: Transient and steady-state solutions are obtained from Equations (8) and (9), respectively, where the relevant matrices are now given by

Numerical Examples
We now demonstrate the robust performance of the developed strategy by means of several two and three-dimensional examples.

2D Lid-Driven Cavity Problem in the Presence of a Magnetic Field
The domain of the problem is a square of dimension h = 1 m. The schematic for this problem is shown in Figure 1.
Three different uniform meshes of nine-noded quadrilateral elements have been used to mesh the geometry, as shown in Table 1, to study the convergence of the solution with respect to mesh refinement. Ten loadsteps have been used for conducting the steady-state analyses. At every loadstep, convergence is achieved within a maximum of six iterations for all of the meshes considered. In the transient case, convergence is achieved within five iterations at each time step. The trapezoidal rule parameter α in the transient analysis is chosen as 0.5. The solution obtained using a 120 × 120 mesh and a time step of 0.125 s is used as the reference solution. In Table 1, the error norm E steady is calculated on the basis of the computed steady-state velocity component u x along the midplane x = 0.5 as where the summation is over the number of nodal points along the midplane. For the transient case, we have calculated the error as the average of the above-mentioned error measure over the first 5 seconds. From Figure 2, the rate of convergence of the error is found to be approximately equal to 1.360 and 1.084 for the steady and transient cases, respectively.

Hartman-Poiseuille Flow
This is one of the standard benchmark problems used in the literature. A conducting fluid flows between two parallel plates under the influence of a pressure gradient given by −ρG, and a magnetic field perpendicular to the two plates H = (B o /µ m )e y . The analytical solution for an incompressible fluid is given by [33].
where Re m = µ m σ e Vh, Ha = B o h σ e /µ and η = y/h. Depending on whether the walls are insulating or conducting, the constants V and E in Equation (10) , E = 0 (Perfectly conducting walls).
We use the following boundary conditions in our simulation, which correspond to the perfectly insulating case: At y = ±0.5, we have u = 0, H x = 0, and T = 293.15 K, whereas at x = 0, 1, the traction corresponding to a uniform pressure, and H y = B 0 /µ m , are imposed. We use a 2 × 64 element mesh for Ha = 1, 2, 5, 10, 20 and 100. Figure 5 shows almost an exact match of the numerical results obtained with the analytical solution. Convergence is achieved within the first four iterations for all of the Hartmann numbers. The same problem has been solved by Nizar [6] using a mesh of 1600 nodes, by Shadid et al. [1] using a 200 × 200 element mesh for Ha = 20, and by Eswaran et al. [28] using a mesh of 5151 nodes, whereas, in the present formulation, only a 645 node mesh is used.

Two-Dimensional Buoyancy-Driven Flow Problem
We consider the case of buoyancy-driven flow inside a rectangular cavity with an aspect ratio of H/L = 0.25. The schematic for this problem is shown in Figure 6. The vertical walls are at constant temperature T h and T c , whereas the horizontal walls are thermally insulated. The magnetic field B 0 acts along the y-direction. The analytical solution for the non-dimensional velocity along the x-direction along the center vertical plane is [35,36] where U = Ha 2 Gr The flow properties used for this problem are T h = 30 • C, T c = 10 • C, C v = 10 J/kg-K, µ = 10 −3 kg/m-s, k f = 1 W/m-K, β = 207 × 10 −6 K −1 , ρ = 998.2 kg/m 3 , µ r = 7.9577 × 10 3 , σ e = 10 (Ω-m) −1 , and g = 1.552083478 × 10 −4 m/s 2 , which correspond to Pr = 0.01, Gr = 2 × 10 4 , and Ha = 1000 (B 0 = 20 Tesla). Lower Hartmann numbers are simulated using a corresponding lower value for B 0 .
The mesh specifications for the different Hartmann numbers considered are given in Table 2.   (1-1000). Eswaran et al. [28] used a graded mesh with 25,351 nodes for Ha = 10, whereas, in the present case, we use a coarse uniform mesh with 3321 nodes to obtain a solution that is in very good agreement with the analytical solution. Figures 8 and 9 present, respectively, the isotherms and streamlines for Ha = 1, Ha = 10, and Ha = 100.

3D Natural Convection Problem in the Presence of a Magnetic Field
This problem is based on experiments performed by Okada et al. [37]. Gallium as a working substance is kept in a 30 mm × 30 mm × 30 mm cubical cavity. A temperature difference (T h − T c ) is applied across the fluid, and a magnetic field B 0 is applied in a direction parallel to the temperature difference as shown in Figure 10. The flow parameters used are T h = 308.15 K, T c = 278.15 K, µ = 10 −3 kg/m-s, ν = 10 −6 m 2 /s, β = 207 × 10 −6 K −1 , C v = 756 J/kg-K, κ f = 31.5 W/m 2 -K, µ r = 7.9577 × 10 3 , σ e = 3.85 × 10 6 (Ω-m) −1 , Pr = 0.024. The values of g and B 0 for different values of Ha are given in Table 3, where Q net is the net rate of heat transfer (W), D is the side of the cube (30 mm), and A is the cross-sectional area (900 mm 2 ). The values are chosen such that the modified Grashoff number Gr * varies from 10 6 to 4.24 × 10 6 for different Ha as in [37]. Meng et al. [38] have solved this problem numerically for buoyant MHD flows at high Ha. We use a graded mesh of 23,273 nodes and 2592 27-noded hexahedral elements for all Ha values to demonstrate the good coarse-mesh accuracy of the present formulation. Figure 11 shows a very good agreement for the Nusselt number at various Gr * and Ha.  Figure 11. Comparison of the numerical results with the experimental results of Okada et al. [37] and computational results of Meng et al. [38] for Nusselt number variation with Gr * at different Ha for the 3D natural convection problem.

Discussion
An appropriate choice of the interpolation functions for the various field variables is critical in ensuring the stability of the resulting numerical scheme. The weak implementation of the state equation used in this work ensures that a lower order interpolation can be used for the pressure, density, and temperature fields compared to the interpolation for the velocity and magnetic fields. The resulting strategy yields stable solutions for all of the field variables without the need for any stabilizing terms in the formulation, or without the need to adjust any parameters, even with coarse meshes. The monolithic nature of the strategy and the exact linearization of the variational formulation ensure a faster rate of convergence compared to a staggered approach.
Based on the numerical solutions presented, it can be said that the present formulation shows a very good coarse-mesh accuracy, and is robust and efficient even at very high Hartmann numbers.