A Monolithic Finite Element Formulation for Magnetohydrodynamics Involving a Compressible Fluid

Abstract

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.

1. Introduction

The flow of a conducting fluid in the presence of a magnetic field is termed as magnetohydrodynamics (MHD). There are several applications of MHD, including the flow of liquid metals [1,2,3,4,5,6,7,8,9,10,11,12], nanofluids [13,14,15,16,17,18], non-Newtonian fluids [19,20,21,22,23], magneto-gas dynamics [24,25,26], astrophysical fluid dynamics [27], and many biomedical devices. References [14,15] use a control-volume-based finite element strategy (CVFEM) in their formulation. References [10,11,12,24] use edge elements for discretizing the magnetic field. Ciuca et al. [26] use hybridized discontinuous Galerkin methods for ideal and resistive MHD, while Eswaran et al. [28] use a finite volume approach.

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.

2. Mathematical Formulation

2.1. 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.

2.1.1. Magnetic Field Equations

The Maxwell equations for electromagnetism under MHD assumptions reduce to the following set of equations [29]:

$$\begin{array}{c}{\mu }_m\frac{\partial \mathit{H}}{\partial t}+\nabla \times \mathit{E}=\mathbf{0},\end{array}$$

(1a)

$$\begin{array}{c}\nabla ·\left({\mu }_m\mathit{H}\right)=0,\end{array}$$

(1b)

$$\begin{array}{c}\nabla \times \mathit{H}=\mathit{j},\end{array}$$

(1c)

$$\begin{array}{c}\nabla ·\mathit{j}=0,\end{array}$$

(1d)

where $\mathit{H}$ and $\mathit{E}$ represent the magnetic and electric fields, respectively, ${\mu }_m$ is the magnetic permeability, and $\mathit{j}$ is the current density. The above equations are supplemented by Ohm’s law for a conducting fluid, which is given by

$$\mathit{j}={\sigma }_e(\mathit{E}+{\mu }_m\mathit{u}\times \mathit{H}),$$

(2)

where ${\sigma }_e$ and $\mathit{u}$ denote the electrical conductivity and the velocity of the fluid, respectively. Combining Equations (1) and (2), we obtain the governing differential equation for $\mathit{H}$ as

$$\nabla \times \left(\nabla \times \mathit{H}\right)+{\sigma }_e{\mu }_m\frac{\partial \mathit{H}}{\partial t}={\sigma }_e{\mu }_m\nabla \times \left(\mathit{u}\times \mathit{H}\right)\mathrm{on}\mathsf{\Omega },$$

(3)

where $\mathsf{\Omega }$ represents the fluid domain.

2.1.2. Coupled Equations for Compressible MHD

Let $\mathsf{\Omega }$ and $\mathsf{\Gamma }$ denote the domain and its boundary, ${\mathsf{\Gamma }}_t$ and ${\mathsf{\Gamma }}_{\mathit{u}}$ denote the parts of boundaries where tractions and velocity are prescribed, respectively, ${\mathsf{\Gamma }}_{\theta }$ and ${\mathsf{\Gamma }}_{\mathit{q}}$ represent portions of $\mathsf{\Gamma }$ where temperature and normal heat flux are prescribed, respectively, and ${\mathsf{\Gamma }}_H$ and ${\mathsf{\Gamma }}_E$ represent portions of the boundary where $\mathit{H}\times \mathit{n}$ and $(\nabla \times \mathit{H})\times \mathit{n}$ are prescribed, respectively.

We are interested in finding an approximate solution to the following initial-boundary value problem:

Determine the velocity $\mathit{u}$, density $\rho$, pressure p, temperature $\theta$, magnetic field $\mathit{H}$, stresses $\tau$, rate of deformation $\mathit{D}$, and heat flux $\mathit{q}$ on the domain $\mathsf{\Omega }$ such that

$$\begin{array}{c}\nabla \times \left(\nabla \times \mathit{H}\right)+{\sigma }_e{\mu }_m\frac{\partial \mathit{H}}{\partial t}={\sigma }_e{\mu }_m\nabla \times \left(\mathit{u}\times \mathit{H}\right),\end{array}$$

(4a)

$$\begin{array}{c}\nabla ·\left({\mu }_m\mathit{H}\right)=0,\end{array}$$

(4b)

$$\begin{array}{c}\left[\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\rho (\nabla \mathit{u})\mathit{u}+\mathit{u}\nabla ·\left(\rho \mathit{u}\right)\right]=\nabla ·\tau +\rho \mathit{b}+{\mu }_m(\nabla \times \mathit{H})\times \mathit{H}\end{array}$$

(4c)

$$\begin{array}{c}\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=\mathbf{0},\end{array}$$

(4d)

$$\begin{array}{c}{C}_v\left[\frac{\partial \left(\rho \theta \right)}{\partial t}+\rho \theta \nabla ·\mathit{u}+\mathit{u}·\left(\rho \nabla \theta +\theta \nabla \rho \right)\right]=-p\nabla ·\mathit{u}+\sigma :\mathit{D}-\nabla ·\mathit{q}+\rho Q_h+\frac{{\left|\nabla \times \mathit{H}\right|}^2}{{\sigma }_e},\end{array}$$

(4e)

$$\begin{array}{c}p=\tilde{p}(\rho ,\theta ),\end{array}$$

(4f)

where $\mathit{b}$ is the body force per unit mass, $\sigma$ is the viscous stress, and $C_v$ is the specific heat at constant volume. For a Newtonian fluid, $\tau =-p(\rho ,\theta )\mathit{I}+\sigma$, where $\sigma =\lambda \left(\mathrm{tr}\mathit{D}\right)\mathit{I}+2\nu \mathit{D}$, with $\lambda$ and $\mu$ denoting the viscosity coefficients. Similarly, we have $\mathit{q}=-\kappa \nabla \theta$, where $\kappa$ is the thermal conductivity. The last term in Equation (4c) represents the additional body force, namely the Lorentz force ${\mu }_m(\mathit{j}\times \mathit{H})$, that acts on the fluid. A term $\mathit{j}·\mathit{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 $\mathit{u}$ from the first law of thermodynamics [32], and using the properties of the scalar triple product and Equation (2), we have

$$\begin{array}{cc}\hfill \mathit{j}·\mathit{E}-{\mu }_m[\mathit{u}·(\mathit{j}\times \mathit{H})]& =\mathit{j}·[\mathit{E}+{\mu }_m\mathit{u}\times \mathit{H}]\hfill \\ \hfill & =\frac{\mathit{j}·\mathit{j}}{{\sigma }_e},\hfill \end{array}$$

which is the last term in Equation (4e). The above governing equations are to be solved subject to appropriate initial and boundary conditions.

2.2. Variational Formulation

Denoting the variations of $\mathit{H}$, $\mathit{u}$, $\rho$, $\theta$, and p with a subscript $\delta$, the variational statements corresponding to Equation (4) can be written as

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\sigma }_e{\mu }_m{\mathit{H}}_{\delta }·\frac{\partial \mathit{H}}{\partial t}d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\left(\nabla \times {\mathit{H}}_{\delta }\right)·\left(\nabla \times \mathit{H}\right)d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\left(\nabla ·{\mathit{H}}_{\delta }\right)\left(\nabla ·\mathit{H}\right)d\mathsf{\Omega }=\\ {\int }_{\mathsf{\Omega }}{\sigma }_e{\mu }_m\left(\nabla \times {\mathit{H}}_{\delta }\right)·\left(\mathit{u}\times \mathit{H}\right)d\mathsf{\Omega }\forall {\mathit{H}}_{\delta },\end{array}$$

(5a)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T\left[\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\rho (\nabla \mathit{u})\mathit{u}+\mathit{u}\nabla ·\left(\rho \mathit{u}\right)\right]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{u}}_{\delta })pd\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\mathbb{C}}_c{\mathit{D}}_{c}d\mathsf{\Omega }=\\ {\int }_{\mathsf{\Omega }}\rho {\mathit{u}}_{\delta }^T\mathit{b}d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\mu }_m{\mathit{u}}_{\delta }^T[(\nabla \times \mathit{H})\times \mathit{H}]d\mathsf{\Omega }+{\int }_{{\mathsf{\Gamma }}_t}{\mathit{u}}_{\delta }^T\overline{\mathit{t}}d\mathsf{\Gamma }\forall {\mathit{u}}_{\delta },\end{array}$$

(5b)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\left[\frac{\partial \rho }{\partial t}+\rho \nabla ·\mathit{u}+\mathit{u}·\nabla \rho \right]d\mathsf{\Omega }=0\forall {\rho }_{\delta },\end{array}$$

(5c)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }\left[\frac{\partial \left(\rho \theta \right)}{\partial t}+\rho \theta \nabla ·\mathit{u}+\mathit{u}·\left(\rho \nabla \theta +\theta \nabla \rho \right)\right]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\theta }_{\delta }\left(p\nabla ·\mathit{u}-\sigma :\mathit{D}\right)d\mathsf{\Omega }+\\ {\int }_{\mathsf{\Omega }}k\nabla {\theta }_{\delta }·\nabla \theta d\mathsf{\Omega }={\int }_{\mathsf{\Omega }}\rho {\theta }_{\delta }{Q}_{h}d\mathsf{\Omega }-{\int }_{{\mathsf{\Gamma }}_q}{\theta }_{\delta }{\overline{q}}_{n}d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}\frac{{\theta }_{\delta }{\left|\nabla \times \mathit{H}\right|}^2}{{\sigma }_e}d\mathsf{\Omega }\forall {\theta }_{\delta },\end{array}$$

(5d)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{p}_{\delta }\left[p-\tilde{p}\left(\rho ,\theta \right)\right]d\mathsf{\Omega }=0\forall p_{\delta },\end{array}$$

(5e)

where ${\mathbb{C}}_c$ and ${\mathit{D}}_c$ denote the material constitutive tensor and rate of deformation tensor expressed in an engineering form as

$${\mathit{D}}_c=\left[\begin{array}{c}{D}_{xx}\\ D_{yy}\\ D_{zz}\\ 2D_{xy}\\ 2D_{yz}\\ 2D_{xz}\end{array}\right],{\mathbb{C}}_c=\left[\begin{array}{cccccc}\lambda +2\mu & \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2\mu & 0& 0& 0\\ 0& 0& 0& \mu & 0& 0\\ 0& 0& 0& 0& \mu & 0\\ 0& 0& 0& 0& 0& \mu \end{array}\right],$$

${\overline{q}}_n$ is the normal heat flux prescribed on ${\mathsf{\Gamma }}_q$, and $\overline{\mathit{t}}$ is the prescribed traction on ${\mathsf{\Gamma }}_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 ${(\nabla ·\mathit{H})}^2$, as shown in [29]. In addition, the term

$${\int }_{{\mathsf{\Gamma }}_H}\left({\mathit{H}}_{\delta }\times \mathit{n}\right)·\left[{\sigma }_e{\mu }_m\left(\mathit{u}\times \mathit{H}\right)-\left(\nabla \times \mathit{H}\right)\right]d\mathsf{\Gamma }=-{\int }_{{\mathsf{\Gamma }}_H}{\sigma }_e({\mathit{H}}_{\delta }\times \mathit{n})·\mathit{E}d\mathsf{\Gamma },$$

has been excluded from Equation (5a), since, in the problems that we consider here, either $\mathit{H}\times \mathit{n}$ is specified on the boundary, or the surface is purely conducting.

2.3. 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:

$$\begin{array}{c}{\sigma }_e{\mu }_m\left({\int }_{\mathsf{\Omega }}{\mathit{H}}_{\delta }·\frac{{\mathit{H}}^{n+1}-{\mathit{H}}^n}{t_{\Delta }}d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·{(\mathit{u}\times \mathit{H})}_{\alpha }d\mathsf{\Omega }\right)+\\ {\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·(\nabla \times {\mathit{H}}_{\alpha })d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{H}}_{\delta })(\nabla ·{\mathit{H}}_{\alpha })d\mathsf{\Omega }=0\forall {\mathit{H}}_{\delta },\end{array}$$

(6a)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }·\left[\frac{{\left(\rho \mathit{u}\right)}^{n+1}-{\left(\rho \mathit{u}\right)}^n}{t_{\Delta }}+{\left[\rho (\nabla \mathit{u})\mathit{u}\right]}_{\alpha }+{[\mathit{u}\nabla ·\left(\rho \mathit{u}\right)]}_{\alpha }\right]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\mathbb{C}}_c{\left({\mathit{D}}_c\right)}_{\alpha }d\mathsf{\Omega }-\\ {\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{u}}_{\delta })p_{\alpha }d\mathsf{\Omega }={\int }_{\mathsf{\Omega }}{\left(\rho {\mathit{u}}_{\delta }^T\mathit{b}\right)}_{\alpha }d\mathsf{\Omega }+{\int }_{\mathsf{\Gamma }}{\left({\mathit{u}}_{\delta }^T\overline{\mathit{t}}\right)}_{\alpha }d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}{\mu }_m{\mathit{u}}_{\delta }^T{[(\nabla \times \mathit{H})\times \mathit{H}]}_{\alpha }d\mathsf{\Omega }\forall {\mathit{u}}_{\delta },\end{array}$$

(6b)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\left[\frac{{\rho }^{n+1}-{\rho }^n}{t_{\Delta }}+{(\rho \nabla ·\mathit{u})}_{\alpha }+{[\mathit{u}·(\nabla \rho )]}_{\alpha }\right]d\mathsf{\Omega }=0\forall {\rho }_{\delta },\end{array}$$

(6c)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }\left[\frac{{\left(\rho \theta \right)}^{n+1}-{\left(\rho \theta \right)}^n}{t_{\Delta }}+{(\rho \theta \nabla ·\mathit{u})}_{\alpha }+{[\mathit{u}·(\rho \nabla \theta +\theta \nabla \rho )]}_{\alpha }\right]+{\int }_{\mathsf{\Omega }}\nabla {\theta }_{\delta }·{[k\nabla \theta ]}_{\alpha }d\mathsf{\Omega }+\\ {\int }_{\mathsf{\Omega }}{\theta }_{\delta }{p}_{\alpha }\frac{\nabla ·{\mathit{u}}^{n+1}+\nabla ·{\mathit{u}}^n}{2}d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{\theta }_{\delta }\left[{\sigma }_{\alpha }:\frac{{\mathit{D}}_c\left({\mathit{u}}^{n+1}\right)+{\mathit{D}}_c\left({\mathit{u}}^n\right)}{2}\right]d\mathsf{\Omega }=\\ {\int }_{\mathsf{\Omega }}{\theta }_{\delta }{\rho }_{\alpha }{Q}_{h}d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\theta }_{\delta }\frac{{\left|\nabla \times \mathit{H}\right|}_{\alpha }^2}{{\sigma }_e}d\mathsf{\Omega }-{\int }_{{\mathsf{\Gamma }}_q}{\theta }_{\delta }{\overline{q}}_{\alpha }d\mathsf{\Gamma }\forall {\theta }_{\delta },\end{array}$$

(6d)

where ${\beta }_{\alpha }=(1-\alpha ){\beta }^n+\alpha {\beta }^{n+1}$ for any field variable $\beta$.

2.4. Linearization

A Newton–Raphson strategy is developed by linearizing the above equations. Let $({\mathit{u}}^k,{\rho }^k,{\theta }^k,p^k,{\mathit{H}}^k)$ and $({\mathit{u}}^{k+1},{\rho }^{k+1},{\theta }^{k+1},p^{k+1},{\mathit{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 $({\mathit{u}}_{\Delta },{\rho }_{\Delta },{\theta }_{\Delta },p_{\Delta },{\mathit{H}}_{\Delta })$ 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:

$$\begin{array}{c}{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{H}}_{\delta }^T{\mathit{H}}_{\Delta }d\mathsf{\Omega }+\alpha t_{\Delta }\{{\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{H}}_{\delta })(\nabla ·{\mathit{H}}_{\Delta })+(\nabla \times {\mathit{H}}_{\delta })·(\nabla \times {\mathit{H}}_{\Delta })d\mathsf{\Omega }-\\ {\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·({\mathit{u}}^k\times {\mathit{H}}_{\Delta })d\mathsf{\Omega }+{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·({\mathit{H}}^k\times {\mathit{u}}_{\Delta })d\mathsf{\Omega }\}\\ ={\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{H}}_{\delta }^T\left({\mathit{H}}^n-{\mathit{H}}^k\right)d\mathsf{\Omega }+\alpha t_{\Delta }\{{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·({\mathit{u}}^k\times {\mathit{H}}^k)d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}[(\nabla \times {\mathit{H}}_{\delta })·(\nabla \times {\mathit{H}}^k)-\\ (\nabla ·\mathit{H})(\nabla ·{\mathit{H}}^k)\left]d\mathsf{\Omega }\right\}+(1-\alpha )t_{\Delta }\{{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}(\nabla \times {\mathit{H}}_{\delta })·({\mathit{u}}^n\times {\mathit{H}}^n)d\mathsf{\Omega }-\\ {\int }_{\mathsf{\Omega }}[(\nabla \times {\mathit{H}}_{\delta })·(\nabla \times {\mathit{H}}^n)-(\nabla ·{\mathit{H}}_{\delta })(\nabla ·{\mathit{H}}^n)]d\mathsf{\Omega }\}\forall {\mathit{H}}_{\delta },\end{array}$$

(7a)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T{\rho }^k{\mathit{u}}_{\Delta }d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T{\mathit{u}}^k{\rho }_{\Delta }d\mathsf{\Omega }+\alpha t_{\Delta }\{{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T[{\rho }^k(\nabla {\mathit{u}}^k){\mathit{u}}_{\Delta }+{\rho }^k(\nabla {\mathit{u}}_{\Delta }){\mathit{u}}^k+{\mathit{u}}^k\nabla ·\left({\rho }^k{\mathit{u}}_{\Delta }\right)+\\ \nabla ·\left({\rho }^k{\mathit{u}}^k\right){\mathit{u}}_{\Delta }]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\mathbb{C}}_c\left[{\mathit{D}}_c\left({\mathit{u}}_{\Delta }\right)\right]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\mathbb{C}}_c^{\prime }{\mathit{D}}_c\left({\mathit{u}}^k\right){\theta }_{\Delta }d\mathsf{\Omega }+\\ {\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T[{\mathit{u}}^k\nabla ·\left({\mathit{u}}^k{\rho }_{\Delta }\right)+(\nabla {\mathit{u}}^k){\mathit{u}}^k{\rho }_{\Delta }-\mathit{b}{\rho }_{\Delta }]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{\mu }_m{\mathit{u}}_{\delta }^T[(\nabla \times {\mathit{H}}^k)\times {\mathit{H}}_{\Delta }-\\ {\mathit{H}}^k\times (\nabla \times {\mathit{H}}_{\Delta })]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{u}}_{\delta })p_{\Delta }d\mathsf{\Omega }\}=-{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T[{\rho }^k{\mathit{u}}^k-{\rho }^n{\mathit{u}}^n]+\alpha t_{\Delta }\{{\int }_{{\mathsf{\Gamma }}_t}{\mathit{u}}_{\delta }^T{\overline{\mathit{t}}}^{k}d\mathsf{\Gamma }+\\ {\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T[{\rho }^k\mathit{b}-{\rho }^k(\nabla {\mathit{u}}^k){\mathit{u}}^k-{\mathit{u}}^k\nabla ·\left({\rho }^k{\mathit{u}}^k\right)+{\mu }_m(\nabla \times {\mathit{H}}^k)\times {\mathit{H}}^k]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\tau }_c^{k}d\mathsf{\Omega }\}+\\ (1-\alpha )t_{\Delta }\{{\int }_{{\mathsf{\Gamma }}_t}{\mathit{u}}_{\delta }^T{\overline{\mathit{t}}}^{n}d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}{\mathit{u}}_{\delta }^T[{\rho }^n\mathit{b}-{\rho }^n(\nabla {\mathit{u}}^n){\mathit{u}}^n-{\mathit{u}}^n\nabla ·\left({\rho }^n{\mathit{u}}^n\right)+{\mu }_m(\nabla \times {\mathit{H}}^n)\times {\mathit{H}}^n]d\mathsf{\Omega }\\ -{\int }_{\mathsf{\Omega }}{\left[{\mathit{D}}_c\left({\mathit{u}}_{\delta }\right)\right]}^T{\tau }_c^{n}d\mathsf{\Omega }\}\forall {\mathit{u}}_{\delta },\end{array}$$

(7b)

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}{\rho }_{\delta }{\rho }_{\Delta }d\mathsf{\Omega }+\alpha t_{\Delta }[{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\left[{\rho }^k(\nabla ·{\mathit{u}}_{\Delta })+(\nabla {\rho }^k)·{\mathit{u}}_{\Delta }\right]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\rho }_{\delta }[(\nabla ·{\mathit{u}}^k){\rho }_{\Delta }+\\ {\mathit{u}}^k·(\nabla {\rho }_{\Delta })\left]d\mathsf{\Omega }\right]=-\alpha t_{\Delta }{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\nabla ·\left({\rho }^k{\mathit{u}}^k\right)d\mathsf{\Omega }-(1-\alpha )t_{\Delta }{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\nabla ·\left({\rho }^n{\mathit{u}}^n\right)d\mathsf{\Omega }\\ -{\int }_{\mathsf{\Omega }}{\rho }_{\delta }\left({\rho }^k-{\rho }^n\right)d\mathsf{\Omega }\forall {\rho }_{\delta },\end{array}$$

(7c)

$$\begin{array}{c}\frac{t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }[(\alpha p^k+(1-\alpha )p^n)\nabla ·{\mathit{u}}_{\Delta }-{\left[{\mathit{D}}_c\left({\mathit{u}}^n\right)\right]}^T{\mathbb{C}}_c{\mathit{D}}_c\left({\mathit{u}}_{\Delta }\right)-2\alpha {\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)\right]}^T{\mathbb{C}}_c{\mathit{D}}_c\left({\mathit{u}}_{\Delta }\right)]d\mathsf{\Omega }+\\ {\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }{\theta }^k{\rho }_{\Delta }d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }{\rho }^k{\theta }_{\Delta }d\mathsf{\Omega }-\frac{\alpha t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }{\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)\right]}^T{\mathbb{C}}_c^{\prime }[{\mathit{D}}_c\left({\mathit{u}}^k\right)+{\mathit{D}}_c\left({\mathit{u}}^n\right)]{\theta }_{\Delta }d\mathsf{\Omega }+\\ \frac{\alpha t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }[\nabla ·{\mathit{u}}^k+\nabla ·{\mathit{u}}^n]p_{\Delta }d\mathsf{\Omega }+\alpha t_{\Delta }\{{\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }[{\rho }^k{\theta }^k(\nabla ·{\mathit{u}}_{\Delta })+({\rho }^k\nabla {\theta }^k+{\theta }^k\nabla {\rho }^k)·{\mathit{u}}_{\Delta }]d\mathsf{\Omega }+\\ {\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }[{\mathit{u}}^k.\nabla {\theta }^k{\rho }_{\Delta }+(\nabla ·{\mathit{u}}^k){\theta }^k{\rho }_{\Delta }+{\theta }^k{\left({\mathit{u}}^k\right)}^T\nabla {\rho }_{\Delta }]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{\theta }_{\delta }{\mathit{Q}}_h{\rho }_{\Delta }d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{C}_v{\theta }_{\delta }\\ [{\rho }^k(\nabla ·{\mathit{u}}^k){\theta }_{\Delta }+(\nabla {\rho }^k·{\mathit{u}}^k){\theta }_{\Delta }+{\rho }^k({\mathit{u}}^k·\nabla {\theta }_{\Delta })]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\nabla {\theta }_{\delta }·[k_f\nabla {\theta }_{\Delta }+k_f^{\prime }\nabla {\theta }^k{\theta }_{\Delta }]d\mathsf{\Omega }-\\ \frac{2}{{\sigma }_e}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }(\nabla \times {\mathit{H}}^k)·(\nabla \times {\mathit{H}}_{\Delta })d\mathsf{\Omega }\}=\alpha t_{\Delta }\{{\int }_{\mathsf{\Omega }}{\theta }_{\delta }[{\rho }^k{\mathit{Q}}_h+\frac{|\nabla \times {\mathit{H}}^k{|}^2}{{\sigma }_e}-\\ C_v[{\rho }^k{\theta }^k(\nabla ·{\mathit{u}}^k)+{\mathit{u}}^k·({\rho }^k\nabla {\theta }^k+{\theta }^k\nabla {\rho }^k)]]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{k}_f\nabla {\theta }_{\delta }·\nabla {\theta }^{k}d\mathsf{\Omega }-{\int }_{\mathsf{\Gamma }}{\theta }_{\delta }{\overline{q}}^{k}d\mathsf{\Gamma }\}\\ +(1-\alpha )t_{\Delta }\{{\int }_{\mathsf{\Omega }}{\theta }_{\delta }\left[{\rho }^n{\mathit{Q}}_h+\frac{|\nabla \times {\mathit{H}}^n{|}^2}{{\sigma }_e}-C_v[{\rho }^n{\theta }^n(\nabla ·{\mathit{u}}^n)+{\mathit{u}}^n·({\rho }^n\nabla {\theta }^n+{\theta }^n\nabla {\rho }^n)]\right]d\mathsf{\Omega }-\\ {\int }_{\mathsf{\Omega }}{k}_f\nabla {\theta }_{\delta }·\nabla {\theta }^{n}d\mathsf{\Omega }-{\int }_{\mathsf{\Gamma }}{\theta }_{\delta }{\overline{q}}^{n}d\mathsf{\Gamma }\}-\frac{\alpha t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }{\left[{\tau }_c^k\right]}^T\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)+{\mathit{D}}_c\left({\mathit{u}}^n\right)\right]d\mathsf{\Omega }+\\ \frac{(1-\alpha )t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\theta }_{\delta }{\left[{\tau }_c^n\right]}^T\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)+{\mathit{D}}_c\left({\mathit{u}}^n\right)\right]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{c}_v{\theta }_{\delta }[{\rho }^k{\theta }^k-{\rho }^n{\theta }^n]d\mathsf{\Omega }\forall {\theta }_{\delta },\end{array}$$

(7d)

$$\begin{array}{c}\alpha t_{\Delta }{\int }_{\mathsf{\Omega }}{p}_{\delta }\left[p_{\Delta }-{\left.\frac{\partial \tilde{p}}{\partial \rho }\right|}_k{\rho }_{\Delta }-{\left.\frac{\partial \tilde{p}}{\partial \theta }\right|}_k{\theta }_{\Delta }\right]d\mathsf{\Omega }=-\alpha t_{\Delta }{\int }_{\mathsf{\Omega }}{p}_{\delta }\left[p^k-\tilde{p}({\rho }^k,{\theta }^k)\right]d\mathsf{\Omega }\forall p_{\delta },\end{array}$$

(7e)

where $k_f^{\prime }$ and ${\mathbb{C}}_c^{\prime }$ denote the derivatives of $k_f$ and ${\mathbb{C}}_c$ with respect to $\theta$.

2.5. Finite Element Formulation

Let the magnetic field, velocity, density, temperature, pressure fields, and their variations and increments be interpolated as

$$\begin{array}{cccccccccc}\hfill \mathit{H}& =\mathit{N}\widehat{\mathit{H}},\hfill & \hfill \mathit{u}& =\mathit{N}\widehat{\mathit{u}},\hfill & \hfill \rho & ={\mathit{N}}_{\rho }\widehat{\rho },\hfill & \hfill \theta & ={\mathit{N}}_{\theta }\widehat{\theta },\hfill & \hfill p& ={\mathit{N}}_p\widehat{\mathit{p}},\hfill \\ \hfill {\mathit{H}}_{\delta }& =\mathit{N}{\widehat{\mathit{H}}}_{\delta },\hfill & \hfill {\mathit{u}}_{\delta }& =\mathit{N}{\widehat{\mathit{u}}}_{\delta },\hfill & \hfill {\rho }_{\delta }& ={\mathit{N}}_{\rho }{\widehat{\rho }}_{\delta },\hfill & \hfill {\theta }_{\delta }& ={\mathit{N}}_{\theta }{\widehat{\theta }}_{\delta },\hfill & \hfill p_{\delta }& ={\mathit{N}}_p{\widehat{\mathit{p}}}_{\delta },\hfill \\ \hfill {\mathit{H}}_{\Delta }& =\mathit{N}{\widehat{\mathit{H}}}_{\Delta },\hfill & \hfill {\mathit{u}}_{\Delta }& =\mathit{N}{\widehat{\mathit{u}}}_{\Delta },\hfill & \hfill {\rho }_{\Delta }& ={\mathit{N}}_{\rho }{\widehat{\rho }}_{\Delta },\hfill & \hfill {\theta }_{\Delta }& ={\mathit{N}}_{\theta }{\widehat{\theta }}_{\Delta },\hfill & \hfill p_{\Delta }& ={\mathit{N}}_p{\widehat{\mathit{p}}}_{\Delta }.\hfill \end{array}$$

In the case of the 27-noded hexahedral element, the shape functions $\mathit{N}$ for $\mathit{u}$ and $\mathit{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

$$\begin{array}{cc}\hfill {\mathit{D}}_c\left({\mathit{u}}_{\Delta }\right)& =\mathit{B}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill (\nabla {\mathit{u}}_{\Delta }){\mathit{u}}^k& =\mathit{R}{\mathit{B}}_{NL}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill \nabla ·{\mathit{u}}_{\Delta }& ={\mathit{B}}_p{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill {\rho }^k(\nabla ·{\mathit{u}}_{\Delta })+(\nabla {\rho }^k)·{\mathit{u}}_{\Delta }& ={\mathit{B}}_{\rho 1}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill (\nabla ·{\mathit{u}}^k){\rho }_{\Delta }+{\mathit{u}}^k·(\nabla {\rho }_{\Delta })& ={\mathit{B}}_{\rho 2}{\widehat{\rho }}_{\Delta },\hfill \\ \hfill \nabla {\theta }_{\Delta }& ={\mathit{B}}_{\theta }{\widehat{\theta }}_{\Delta },\hfill \\ \hfill \nabla \times {\mathit{H}}_{\Delta }& ={\mathit{B}}_H{\widehat{\mathit{H}}}_{\Delta },\hfill \\ \hfill \nabla ·{\mathit{H}}_{\Delta }& ={\mathit{B}}_p{\widehat{\mathit{H}}}_{\Delta },\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill \mathit{B}& =\left[\begin{array}{ccccccccc}{N}_{1,x}& 0& 0& N_{2,x}& 0& 0& .& .& .\\ 0& N_{1,y}& 0& 0& N_{2,y}& 0& .& .& .\\ 0& 0& N_{1,z}& 0& 0& N_{2,z}& .& .& .\\ N_{1,y}& N_{1,x}& 0& N_{2,y}& N_{2,x}& 0& .& .& .\\ 0& N_{1,z}& N_{1,y}& 0& N_{2,z}& N_{2,y}& .& .& .\\ N_{1,z}& 0& N_{1,x}& N_{2,z}& 0& N_{2,x}& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_p& =\left[\begin{array}{ccccccccc}{N}_{1,x}& N_{1,y}& N_{1,z}& N_{2,x}& N_{2,y}& N_{2,z}& .& .& .\end{array}\right],\hfill \\ \hfill \mathit{R}& =\left[\begin{array}{ccccccccc}{u}_x^k& u_y^k& u_z^k& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& u_x^k& u_y^k& u_z^k& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& u_x^k& u_y^k& u_z^k\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{NL}& =\left[\begin{array}{ccccccccc}{N}_{1,x}& 0& 0& N_{2,x}& 0& 0& .& .& .\\ N_{1,y}& 0& 0& N_{2,y}& 0& 0& .& .& .\\ N_{1,z}& 0& 0& N_{2,z}& 0& 0& .& .& .\\ 0& N_{1,x}& 0& 0& N_{2,x}& 0& .& .& .\\ 0& N_{1,y}& 0& 0& N_{2,y}& 0& .& .& .\\ 0& N_{1,z}& 0& 0& N_{2,z}& 0& .& .& .\\ 0& 0& N_{1,x}& 0& 0& N_{2,x}& .& .& .\\ 0& 0& N_{1,y}& 0& 0& N_{2,y}& .& .& .\\ 0& 0& N_{1,z}& 0& 0& N_{2,z}& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_H& =\left[\begin{array}{ccccccccc}0& -N_{1,z}& N_{1,y}& 0& -N_{2,z}& N_{2,y}& .& .& .\\ N_{1,z}& 0& -N_{1,x}& N_{2,z}& 0& -N_{2,x}& .& .& .\\ -N_{1,y}& N_{1,x}& 0& -N_{2,y}& N_{2,x}& 0& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\rho 1}& =\left[\begin{array}{cccccccc}{\rho }^k{N}_{1,x}+\frac{\partial {\rho }^k}{\partial x}{N}_1& {\rho }^k{N}_{1,y}+\frac{\partial {\rho }^k}{\partial y}{N}_1& {\rho }^k{N}_{1,z}+\frac{\partial {\rho }^k}{\partial z}{N}_1& {\rho }^k{N}_{2,x}+\frac{\partial {\rho }^k}{\partial x}{N}_2& & .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\rho 2}& =\left[\begin{array}{ccccc}(\nabla ·{\mathit{u}}^k)N_1+({\mathit{u}}^k·\nabla )N_1& (\nabla ·{\mathit{u}}^k)N_2+({\mathit{u}}^k·\nabla )N_2& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\rho }& ={\mathit{B}}_{\theta }=\left[\begin{array}{ccccc}{N}_{1,x}& N_{2,x}& .& .& .\\ N_{1,y}& N_{2,y}& .& .& .\\ N_{1,z}& N_{2,z}& .& .& .\end{array}\right].\hfill \end{array}$$

The various cross-product terms occurring in the formulation are discretized as follows:

$$\begin{array}{cc}\hfill {\mathit{u}}^k\times {\mathit{H}}_{\Delta }& ={\mathit{u}}_{\mathrm{mat}}^k\mathit{N}{\widehat{\mathit{H}}}_{\Delta },\hfill \\ \hfill {\mathit{H}}^k\times {\mathit{u}}_{\Delta }& ={\mathit{H}}_{\mathrm{mat}}^k\mathit{N}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill {\mathit{u}}^k\times {\mathit{H}}^k& =-{\mathit{H}}_{\mathrm{mat}}^k{\mathit{u}}^k,\hfill \\ \hfill {\mathit{u}}^n\times {\mathit{H}}^n& =-{\mathit{H}}_{\mathrm{mat}}^n{\mathit{u}}^n,\hfill \\ \hfill {\mathit{G}}^k& =\nabla \times {\mathit{H}}^k,\hfill \\ \hfill (\nabla \times {\mathit{H}}^k)\times {\mathit{H}}^k& =-{\mathit{H}}_{\mathrm{mat}}^k{\mathit{G}}^k,\hfill \\ \hfill (\nabla \times {\mathit{H}}^n)\times {\mathit{H}}^n& =-{\mathit{H}}_{\mathrm{mat}}^n{\mathit{G}}^n,\hfill \\ \hfill {\mathit{H}}^k\times (\nabla \times {\mathit{H}}_{\Delta })& ={\mathit{H}}_{\mathrm{mat}}^k{\mathit{B}}_H{\widehat{\mathit{H}}}_{\Delta },\hfill \end{array}$$

where

$$\begin{array}{cccc}\hfill {\mathit{u}}_{\mathrm{mat}}^k& =\left[\begin{array}{ccc}0& -u_3^k& u_2^k\\ u_3^k& 0& -u_1^k\\ -u_2^k& u_1^k& 0\end{array}\right],\hfill & \hfill {\mathit{u}}_{\mathrm{mat}}^n& =\left[\begin{array}{ccc}0& -u_3^n& u_2^n\\ u_3^n& 0& -u_1^n\\ -u_2^n& u_1^n& 0\end{array}\right],\hfill \\ \hfill {\mathit{H}}_{\mathrm{mat}}^k& =\left[\begin{array}{ccc}0& -H_3^k& H_2^k\\ H_3^k& 0& -H_1^k\\ -H_2^k& H_1^k& 0\end{array}\right],\hfill & \hfill {\mathit{H}}_{\mathrm{mat}}^n& =\left[\begin{array}{ccc}0& -H_3^n& H_2^n\\ H_3^n& 0& -H_1^n\\ -H_2^n& H_1^n& 0\end{array}\right],\hfill \\ \hfill {\mathit{G}}_{\mathrm{mat}}^k& =\left[\begin{array}{ccc}0& -G_3^k& G_2^k\\ G_3^k& 0& -G_1^k\\ -G_2^k& G_1^k& 0\end{array}\right].\hfill \end{array}$$

Using the arbitrariness of the variations, the matrix form of the incremental Equation (7) can be written as

$$(\mathit{M}+\alpha t_{\Delta }{\mathit{K}}^k){\widehat{\mathit{h}}}_{\Delta }=\alpha t_{\Delta }{\mathit{f}}_{\Delta }^k+(1-\alpha )t_{\Delta }{\mathit{f}}_{\Delta }^n+\mathit{g},$$

(8)

where

$$\begin{array}{cc}\hfill \mathit{M}& =\left[\begin{array}{ccccc}{\mathit{M}}_{HH}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathit{M}}_{uu}& {\mathit{M}}_{u\rho }& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathit{M}}_{\rho \rho }& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathit{M}}_{\theta u}& {\mathit{M}}_{\theta \rho }& {\mathit{M}}_{\theta \theta }& {\mathit{M}}_{\theta p}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right],\mathit{K}=\left[\begin{array}{ccccc}{\mathit{K}}_{HH}& {\mathit{K}}_{Hu}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ {\mathit{K}}_{uH}& {\mathit{K}}_{uu}& {\mathit{K}}_{u\rho }& {\mathit{K}}_{u\theta }& K_{up}\\ \mathbf{0}& {\mathit{K}}_{\rho u}& {\mathit{K}}_{\rho \rho }& \mathbf{0}& \mathbf{0}\\ {\mathit{K}}_{\theta H}& {\mathit{K}}_{\theta u}& {\mathit{K}}_{\theta \rho }& {\mathit{K}}_{\theta \theta }& {\mathit{K}}_{\theta p}\\ \mathbf{0}& \mathbf{0}& {\mathit{K}}_{p\rho }& {\mathit{K}}_{p\theta }& {\mathit{K}}_{pp}\end{array}\right],\hfill \\ \hfill {\widehat{\mathit{h}}}_{\Delta }& =\left[\begin{array}{c}{\widehat{\mathit{H}}}_{\Delta }\\ {\widehat{\mathit{u}}}_{\Delta }\\ {\widehat{\rho }}_{\Delta }\\ {\widehat{\theta }}_{\Delta }\\ {\widehat{\mathit{p}}}_{\Delta }\end{array}\right],{\mathit{f}}_{\Delta }^k={\left[\begin{array}{c}{\mathit{f}}_H\\ {\mathit{f}}_u\\ {\mathit{f}}_{\rho }\\ {\mathit{f}}_{\theta }\\ {\mathit{f}}_p\end{array}\right]}^k,{\mathit{f}}_{\Delta }^n={\left[\begin{array}{c}{\mathit{f}}_H\\ {\mathit{f}}_u\\ {\mathit{f}}_{\rho }\\ {\mathit{f}}_{\theta }\\ \mathbf{0}\end{array}\right]}^n,\mathit{g}=\left[\begin{array}{c}{\mathit{g}}_H\\ {\mathit{g}}_u\\ {\mathit{g}}_{\rho }\\ {\mathit{g}}_{\theta }\\ \mathbf{0}\end{array}\right],\hfill \end{array}$$

with

$$\begin{array}{cc}\hfill {\mathit{M}}_{HH}& ={\int }_{\mathsf{\Omega }}{\sigma }_e{\mu }_m{\mathit{N}}^T\mathit{N}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{uu}& ={\int }_{\mathsf{\Omega }}{\rho }^k{\mathit{N}}^T\mathit{N}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{u\rho }& ={\int }_{\mathsf{\Omega }}{\mathit{N}}^T{\mathit{u}}^k{\mathit{N}}_{\rho }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{\rho \rho }& ={\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T{\mathit{N}}_{\rho }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{\theta u}& =\frac{t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}[(1-\alpha )p^n-\alpha p^k]{\mathit{N}}_{\theta }^T{\mathit{B}}_{p}d\mathsf{\Omega }-\frac{t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T{\mathit{D}}_c^T\left({\mathit{u}}^n\right){\mathbb{C}}_c\mathit{B}d\mathsf{\Omega }+\hfill \\ & \alpha t_{\Delta }{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T{\mathit{D}}_c^T\left({\mathit{u}}^k\right){\mathbb{C}}_c\mathit{B}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{\theta \rho }& ={\int }_{\mathsf{\Omega }}{C}_v{\theta }^k{\mathit{N}}_{\theta }^T{\mathit{N}}_{\rho }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{\theta \theta }& ={\int }_{\mathsf{\Omega }}{\rho }^k{C}_v{\mathit{N}}_{\theta }^T{\mathit{N}}_{\theta }d\mathsf{\Omega }-\frac{\alpha t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T{\mathit{D}}_c^T\left({\mathit{u}}^k\right){\mathbb{C}}_c^{\prime }\left[{\mathit{D}}_c\left({\mathit{u}}^n\right)-{\mathit{D}}_c\left({\mathit{u}}^k\right)\right]{\mathit{N}}_{\theta }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{M}}_{\theta p}& =\frac{\alpha t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T\left[\nabla ·{\mathit{u}}^n-\nabla ·{\mathit{u}}^k\right]{\mathit{N}}_{p}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{HH}& ={\int }_{\mathsf{\Omega }}[{\mathit{B}}_H^T{\mathit{B}}_H+{\mathit{B}}_p^T{\mathit{B}}_p]d\mathsf{\Omega }-{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T{\mathit{u}}_{\mathrm{mat}}^k\mathit{N}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{Hu}& ={\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T{\mathit{H}}_{\mathrm{mat}}^k\mathit{N}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{uH}& ={\int }_{\mathsf{\Omega }}{\mu }_m{\mathit{N}}^T\left({\mathit{H}}_{\mathrm{mat}}^k{\mathit{B}}_H-{\mathit{G}}_{\mathrm{mat}}^k\mathit{N}\right)d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{uu}& ={\int }_{\mathsf{\Omega }}{\rho }^k{\mathit{N}}^T\left[(\nabla {\mathit{u}}^k)\mathit{N}+\mathit{R}{\mathit{B}}_{NL}\right]d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\mathit{B}}^T{\mathbb{C}}_c\mathit{B}d\mathsf{\Omega }\hfill \\ \hfill & +{\int }_{\mathsf{\Omega }}\left[{\mathit{N}}^T{\mathit{u}}^k\left[{(\nabla {\rho }^k)}^T\mathit{N}+{\rho }^k{\mathit{B}}_p\right]+\left[(\nabla {\rho }^k)·{\mathit{u}}^k+{\rho }^k(\nabla ·{\mathit{u}}^k)\right]{\mathit{N}}^T\mathit{N}\right]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{u\rho }& ={\int }_{\mathsf{\Omega }}{\mathit{N}}^T\left[(\nabla {\mathit{u}}^k){\mathit{u}}^k-\mathit{b}\right]{\mathit{N}}_{\rho }d\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}{\mathit{N}}^T\left[{\mathit{u}}^k{\left({\mathit{u}}^k\right)}^T{\mathit{B}}_{\rho }+(\nabla ·{\mathit{u}}^k){\mathit{u}}^k{\mathit{N}}_{\rho }\right]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{u\theta }& ={\int }_{\mathsf{\Omega }}{\mathit{B}}^T{\mathbb{C}}_c^{\prime }{\mathit{D}}_c\left({\mathit{u}}^k\right){\mathit{N}}_{\theta }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{up}& =-{\int }_{\mathsf{\Omega }}{\mathit{B}}_p^T{\mathit{N}}_{p}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\rho u}& ={\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T{\mathit{B}}_{\rho 1}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\rho \rho }& ={\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T{\mathit{B}}_{\rho 2}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta H}& =-\frac{2}{{\sigma }_e}{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T{\left({\mathit{G}}^k\right)}^T{\mathit{B}}_{H}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta u}& ={\int }_{\mathsf{\Omega }}{C}_v{\mathit{N}}_{\theta }^T\left[{\left[{\theta }^k\nabla {\rho }^k+{\rho }^k\nabla {\theta }^k\right]}^T\mathit{N}+{\rho }^k{\theta }^k{\mathit{B}}_p\right]d\mathsf{\Omega }\hfill \\ \hfill & +{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T[p^k{\mathit{B}}_p-2{\mathit{D}}_c^T\left({\mathit{u}}^k\right){\mathbb{C}}_c\mathit{B}]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta \rho }& ={\int }_{\mathsf{\Omega }}{C}_v{\mathit{N}}_{\theta }^T\left[{\theta }^k{\left({\mathit{u}}^k\right)}^T{\mathit{B}}_{\rho }+{\theta }^k(\nabla ·{\mathit{u}}^k){\mathit{N}}_{\rho }+{\mathit{u}}^k·(\nabla {\theta }^k){\mathit{N}}_{\rho }\right]d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}{Q}_h{\mathit{N}}_{\theta }^T{\mathit{N}}_{\rho }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta \theta }& ={\int }_{\mathsf{\Omega }}{C}_v{\mathit{N}}_{\theta }^T[{\rho }^k{\left({\mathit{u}}^k\right)}^T{\mathit{B}}_{\theta }+(\nabla {\rho }^k)·{\mathit{u}}^k+{\rho }^k(\nabla ·{\mathit{u}}^k)]{\mathit{N}}_{\theta }d\mathsf{\Omega }\hfill \\ \hfill & +{\int }_{\mathsf{\Omega }}[k_f{\mathit{B}}_{\theta }^T{\mathit{B}}_{\theta }+k_f^{\prime }{\mathit{B}}_{\theta }^T\nabla {\theta }^k{\mathit{N}}_{\theta }-{\mathit{N}}_{\theta }^T{\mathit{D}}_c^T\left({\mathit{u}}^k\right){\mathbb{C}}_c^{\prime }{\mathit{D}}_c\left({\mathit{u}}^k\right){\mathit{N}}_{\theta }]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta p}& ={\int }_{\mathsf{\Omega }}(\nabla ·{\mathit{u}}^k){\mathit{N}}_{\theta }^T{\mathit{N}}_{p}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{p\rho }& =-{\int }_{\mathsf{\Omega }}{\left.\frac{\partial \tilde{p}}{\partial \rho }\right|}_k{\mathit{N}}_p^T{\mathit{N}}_{\rho }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{p\theta }& =-{\int }_{\mathsf{\Omega }}{\left.\frac{\partial \tilde{p}}{\partial \theta }\right|}_k{\mathit{N}}_p^T{\mathit{N}}_{\theta }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{pp}& ={\int }_{\mathsf{\Omega }}{\mathit{N}}_p^T{\mathit{N}}_{p}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_H^k& =-{\int }_{\mathsf{\Omega }}[{\mathit{B}}_H^T(\nabla \times {\mathit{H}}^k)+{\mathit{B}}_p^T(\nabla ·{\mathit{H}}^k)]d\mathsf{\Omega }-{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T{\mathit{H}}_{\mathrm{mat}}^k{\mathit{u}}^{k}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_u^k& ={\int }_{\mathsf{\Omega }}\{{\mathit{N}}^T\left[{\rho }^k\mathit{b}-{\rho }^k(\nabla {\mathit{u}}^k){\mathit{u}}^k-\left[(\nabla {\rho }^k)·{\mathit{u}}^k+(\nabla ·{\mathit{u}}^k){\rho }^k\right]{\mathit{u}}^k-{\mu }_m{\mathit{H}}_{\mathrm{mat}}^k{\mathit{G}}^k\right]\hfill \\ \hfill & -{\mathit{B}}^T{\tau }_c^k\}d\mathsf{\Omega }+{\int }_{\mathsf{\Gamma }}{\mathit{N}}^T{\overline{\mathit{t}}}^{k}d\mathsf{\Gamma },\hfill \\ \hfill {\mathit{f}}_{\rho }^k& =-{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T\left[{\rho }^k(\nabla ·{\mathit{u}}^k)+{\mathit{u}}^k·(\nabla {\rho }^k)\right]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_{\theta }^k& ={\int }_{\mathsf{\Omega }}\{{\rho }^k{\mathit{N}}_{\theta }^T{Q}_h-k_f{\mathit{B}}_{\theta }^T\nabla {\theta }^k-C_v\left[({\theta }^k\nabla {\rho }^k+{\rho }^k\nabla {\theta }^k)·{\mathit{u}}^k+{\rho }^k{\theta }^k(\nabla ·{\mathit{u}}^k)\right]{\mathit{N}}_{\theta }^T\hfill \\ \hfill & +{\left({\tau }_c^k\right)}^T\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)\right]{\mathit{N}}_{\theta }^T\}d\mathsf{\Omega }-{\int }_{{\mathsf{\Gamma }}_q}{\mathit{N}}_{\theta }^T{\overline{q}}^{k}d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}\frac{{\left({\mathit{G}}^k\right)}^T\left({\mathit{G}}^k\right)}{{\sigma }_e}{\mathit{N}}_{\theta }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_p^k& =-{\int }_{\mathsf{\Omega }}\left[p^k-\tilde{p}({\rho }^k,{\theta }^k)\right]{\mathit{N}}_p^{T}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_H^n& =-{\int }_{\mathsf{\Omega }}[{\mathit{B}}_H^T(\nabla \times {\mathit{H}}^n)+{\mathit{B}}_p^T(\nabla ·{\mathit{H}}^n)]d\mathsf{\Omega }-{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T{\mathit{H}}_{\mathrm{mat}}^n{\mathit{u}}^{n}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_u^n& ={\int }_{\mathsf{\Omega }}\{{\mathit{N}}^T\left[{\rho }^n\mathit{b}-{\rho }^n(\nabla {\mathit{u}}^n){\mathit{u}}^n-\left[(\nabla {\rho }^n)·{\mathit{u}}^n+(\nabla ·{\mathit{u}}^n){\rho }^n\right]{\mathit{u}}^n-{\mu }_m{\mathit{H}}_{\mathrm{mat}}^n{\mathit{G}}^n\right]\hfill \\ \hfill & -{\mathit{B}}^T{\tau }_c^n\}d\mathsf{\Omega }+{\int }_{\mathsf{\Gamma }}{\mathit{N}}^T{\overline{\mathit{t}}}^{n}d\mathsf{\Gamma },\hfill \\ \hfill {\mathit{f}}_{\rho }^n& =-{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T\left[{\rho }^n(\nabla ·{\mathit{u}}^n)+{\mathit{u}}^n·(\nabla {\rho }^n)\right]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_{\theta }^n& ={\int }_{\mathsf{\Omega }}\{{\rho }^n{\mathit{N}}_{\theta }^T{Q}_h-k_f{\mathit{B}}_{\theta }^T\nabla {\theta }^n-C_v\left[({\theta }^n\nabla {\rho }^n+{\rho }^n\nabla {\theta }^n)·{\mathit{u}}^n+{\rho }^n{\theta }^n(\nabla ·{\mathit{u}}^n)\right]{\mathit{N}}_{\theta }^T+\hfill \\ \hfill & {\left({\tau }_c^n\right)}^T\left[{\mathit{D}}_c\left({\mathit{u}}^n\right)\right]{\mathit{N}}_{\theta }^T\}d\mathsf{\Omega }-{\int }_{{\mathsf{\Gamma }}_q}{\mathit{N}}_{\theta }^T{\overline{q}}^{n}d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}\frac{{\left({\mathit{G}}^n\right)}^T\left({\mathit{G}}^n\right)}{{\sigma }_e}{\mathit{N}}_{\theta }d\mathsf{\Omega },\hfill \\ \hfill {\mathit{g}}_H& =-{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{N}}^T({\mathit{H}}^k-{\mathit{H}}^n)d\mathsf{\Omega },\hfill \\ \hfill {\mathit{g}}_u& =-{\int }_{\mathsf{\Omega }}{\mathit{N}}^T[{\rho }^k{\mathit{u}}^k-{\rho }^n{\mathit{u}}^n]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{g}}_{\rho }& =-{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\rho }^T({\rho }^k-{\rho }^n)d\mathsf{\Omega },\hfill \\ \hfill {\mathit{g}}_{\theta }& =-{\int }_{\mathsf{\Omega }}{C}_v\left[{\rho }^k{\theta }^k-{\rho }^n{\theta }^n\right]{\mathit{N}}_{\theta }^{T}d\mathsf{\Omega }+\frac{t_{\Delta }}{2}{\int }_{\mathsf{\Omega }}{\mathit{N}}_{\theta }^T[(1-\alpha ){\tau }_c^n-\alpha {\tau }_c^k][{\mathit{D}}_c\left({\mathit{u}}^k\right)-{\mathit{D}}_c\left({\mathit{u}}^n\right)]d\mathsf{\Omega }.\hfill \end{array}$$

The corresponding steady-state solution can be obtained directly without time stepping by solving

$${\mathit{K}}^k{\widehat{\mathit{h}}}_{\Delta }={\mathit{f}}_{\Delta }^k,$$

(9)

2.6. Two-Dimensional Formulation

In the case of two-dimensional flows in the x-y plane, only ${(\nabla \times \mathit{H})}_z$ is non-zero, and is given by

$$\begin{array}{c}\hfill {(\nabla \times \mathit{H})}_z=\frac{\partial \mathit{H}}{\partial y}-\frac{\partial \mathit{H}}{\partial x}.\end{array}$$

The velocity, magnetic field, rate of deformation tensor, and constitutive tensor are given by

$$\begin{array}{cc}\hfill \mathit{u}& =\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right],\mathit{H}=\left[\begin{array}{c}{H}_x\\ H_y\end{array}\right],{\mathit{D}}_c=\left[\begin{array}{c}{D}_{xx}\\ D_{yy}\\ 2D_{xy}\end{array}\right],{\mathbb{C}}_c=\left[\begin{array}{ccc}\lambda +2\mu & \lambda & 0\\ \lambda & \lambda +2\mu & 0\\ 0& 0& \mu \end{array}\right],\hfill \end{array}$$

whereas the other matrices are as follows:

$$\begin{array}{cc}\hfill \mathit{B}& =\left[\begin{array}{ccccccc}{N}_{1,x}& 0& N_{2,x}& 0& .& .& .\\ 0& N_{1,y}& 0& N_{2,y}& .& .& .\\ N_{1,y}& N_{1,x}& N_{2,y}& N_{2,x}& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_p& =\left[\begin{array}{ccccccc}{N}_{1,x}& N_{1,y}& N_{2,x}& N_{2,y}& .& .& .\end{array}\right],\hfill \\ \hfill \mathit{R}& =\left[\begin{array}{cccc}{u}_x^k& u_y^k& 0& 0\\ 0& 0& u_x^k& u_y^k\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{NL}& =\left[\begin{array}{ccccccc}{N}_{1,x}& 0& N_{2,x}& 0& .& .& .\\ N_{1,y}& 0& N_{2,y}& 0& .& .& .\\ 0& N_{1,x}& 0& N_{2,x}& .& .& .\\ 0& N_{1,y}& 0& N_{2,y}& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\rho 1}& =\left[\begin{array}{ccccccc}{\rho }^k{N}_{1,x}+\frac{\partial {\rho }^k}{\partial x}{N}_1& {\rho }^k{N}_{1,y}+\frac{\partial {\rho }^k}{\partial y}{N}_1& {\rho }^k{N}_{2,x}+\frac{\partial {\rho }^k}{\partial x}{N}_2& {\rho }^k{N}_{2,y}+\frac{\partial {\rho }^k}{\partial y}{N}_2& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\rho 2}& =\left[\begin{array}{ccccc}(\nabla ·{\mathit{u}}^k)N_1+({\mathit{u}}^k·\nabla )N_1& (\nabla ·{\mathit{u}}^k)N_2+({\mathit{u}}^k·\nabla )N_2& .& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_{\theta }& =\left[\begin{array}{cccc}{N}_{1,x}& N_{2,x}& .& .\\ N_{1,y}& N_{2,y}& .& .\end{array}\right],\hfill \\ \hfill {\mathit{B}}_H& =\left[\begin{array}{ccccccc}-N_{1,y}& N_{1,x}& -N_{2,y}& N_{2,x}& .& .& .\end{array}\right].\hfill \end{array}$$

The cross-product terms in the two-dimensional formulation are as follows:

$$\begin{array}{cc}\hfill {\mathit{H}}_{\Delta }\times {\mathit{u}}^k& =-\left[u_x^k{u}_y^k\right]\mathit{S}{\widehat{\mathit{H}}}_{\Delta },\hfill \\ \hfill {\mathit{H}}^k\times {\mathit{u}}_{\Delta }& =\left[H_x^k{H}_y^k\right]\mathit{S}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill {\mathit{u}}^k\times {\mathit{H}}^k& =\left[u_x^k{u}_y^k\right]\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right],\hfill \\ \hfill (\nabla \times {\mathit{H}}^k)\times {\mathit{H}}^k& =-{(\nabla \times {\mathit{H}}^k)}_z\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right],\hfill \\ \hfill {\mathit{H}}^k\times (\nabla \times {\mathit{H}}_{\Delta })& =\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right]{\mathit{B}}_H{\widehat{\mathit{H}}}_{\Delta },\hfill \\ \hfill {\mathit{u}}_{\Delta }\times (\nabla \times {\mathit{H}}^k)& ={(\nabla \times {\mathit{H}}^k)}_z\mathit{S}{\widehat{\mathit{u}}}_{\Delta },\hfill \\ \hfill (\nabla \times {\mathit{H}}^k)\times {\mathit{H}}_{\Delta }& =-{(\nabla \times {\mathit{H}}^k)}_z\mathit{S}{\widehat{\mathit{H}}}_{\Delta },\hfill \end{array}$$

where

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccccc}0& N_1& 0& N_2& .& .& .& 0& N_9\\ -N_1& 0& -N_2& 0& .& .& .& -N_9& 0\end{array}\right].\end{array}$$

Transient and steady-state solutions are obtained from Equations (8) and (9), respectively, where the relevant matrices are now given by

$$\begin{array}{cc}\hfill {\mathit{K}}_{uH}& ={\int }_{\mathsf{\Omega }}{\mu }_m{\mathit{N}}^T\left\{{(\nabla \times {\mathit{H}}^k)}_z\mathit{S}+\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right]{\mathit{B}}_H\right\}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{\theta H}& =-\frac{2}{{\sigma }_e}{\int }_{\mathsf{\Omega }}{(\nabla \times {\mathit{H}}^k)}_z{\mathit{N}}_{\theta }^T{\mathit{B}}_{H}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{Hu}& ={\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T\left[\begin{array}{cc}{H}_x^k& H_y^k\end{array}\right]\mathit{S}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{K}}_{HH}& ={\int }_{\mathsf{\Omega }}[{\mathit{B}}_H^T{\mathit{B}}_H+{\mathit{B}}_p^T{\mathit{B}}_p]d\mathsf{\Omega }-{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T\left[\begin{array}{cc}{u}_x^k& u_y^k\end{array}\right]\mathit{S}d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_H^k& =-{\int }_{\mathsf{\Omega }}[{(\nabla \times {\mathit{H}}^k)}_z{\mathit{B}}_H^T+{\mathit{B}}_p^T(\nabla ·{\mathit{H}}^k)]d\mathsf{\Omega }+{\sigma }_e{\mu }_m{\int }_{\mathsf{\Omega }}{\mathit{B}}_H^T{\left[\begin{array}{c}{u}_x^k\\ u_y^k\end{array}\right]}^T\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right]d\mathsf{\Omega },\hfill \\ \hfill {\mathit{f}}_u^k& ={\int }_{\mathsf{\Omega }}\{{\mathit{N}}^T\left[{\rho }^k\mathit{b}-{\rho }^k(\nabla {\mathit{u}}^k){\mathit{u}}^k-\left[(\nabla {\rho }^k)·{\mathit{u}}^k+(\nabla ·{\mathit{u}}^k){\rho }^k\right]{\mathit{u}}^k-{\mu }_m{(\nabla \times {\mathit{H}}^k)}_z\left[\begin{array}{c}{H}_y^k\\ -H_x^k\end{array}\right]\right]-\hfill \\ & {\mathit{B}}^T{\tau }_c^k\}d\mathsf{\Omega }+{\int }_{\mathsf{\Gamma }}{\mathit{N}}^T{\overline{\mathit{t}}}^{k}d\mathsf{\Gamma },\hfill \\ \hfill {\mathit{f}}_{\theta }^k& ={\int }_{\mathsf{\Omega }}\{{\rho }^k{\mathit{N}}_{\theta }^T{Q}_h-k_f{\mathit{B}}_{\theta }^T\nabla {\theta }^k-C_v\left[({\theta }^k\nabla {\rho }^k+{\rho }^k\nabla {\theta }^k)·{\mathit{u}}^k+{\rho }^k{\theta }^k(\nabla ·{\mathit{u}}^k)\right]{\mathit{N}}_{\theta }^T+\hfill \\ \hfill & {\left({\tau }_c^k\right)}^T\left[{\mathit{D}}_c\left({\mathit{u}}^k\right)\right]{\mathit{N}}_{\theta }^T\}d\mathsf{\Omega }-{\int }_{{\mathsf{\Gamma }}_q}{\mathit{N}}_{\theta }^T{\overline{q}}^{k}d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}\frac{|{(\nabla \times {\mathit{H}}^k)}_z{|}^2}{{\sigma }_e}{\mathit{N}}_{\theta }d\mathsf{\Omega }.\hfill \end{array}$$

3. Numerical Examples

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

3.1. 2D Lid-Driven Cavity Problem in the Presence of a Magnetic Field

The domain of the problem is a square of dimension $h=1\mathrm{m}$. The schematic for this problem is shown in Figure 1.

At time $t=0$ , the top surface moves with a velocity $u_0=1\mathrm{m}/\mathrm{s}$. A uniform magnetic field is applied throughout the domain at all times, such that $\mathit{H}=2{\mathit{e}}_y\mathrm{A}/\mathrm{m}$. The magnetic field component $H_x$ is prescribed to be zero at the top and bottom surfaces at all times. The properties of the conducting fluid used for this example are ${\mu }_r=7.9577\times {10}^7$, ${\sigma }_e=0.01{\left(\mathsf{\Omega }-\mathrm{m}\right)}^{-1}$, $\mu =1.0\mathrm{kg}/\mathrm{ms}$, $\rho =1000\mathrm{kg}/{\mathrm{m}}^3$ corresponding to $\mathrm{Re}=1000$, ${\mathrm{Re}}_{\mathrm{m}}=1$ and $\mathrm{Ha}=20$, where

$$\begin{array}{c}\mathrm{Re}=\frac{\rho u_{0}h}{\mu },{\mathrm{Re}}_{\mathrm{m}}={\mu }_m{\sigma }_e{u}_{0}h,Ha=\sqrt{\frac{{\sigma }_e}{\mu }}{B}_{0}h,\\ B_0={\mu }_m{H}_y,{\mu }_m={\mu }_r{\mu }_0,{\mu }_0=4\pi \times {10}^{-7}.\end{array}$$

Three different uniform meshes of nine-noded quadrilateral elements have been used to mesh the geometry, as shown in

Table 1. Convergence of the solution with respect to mesh and timestep refinement in the 2D lid-driven cavity problem in the presence of a magnetic field.

Mesh Level	Elements	${\mathit{t}}_{\mathit{\Delta }}$	${\mathbf{E}}_{\mathbf{steady}}$	${\mathbf{E}}_{\mathbf{transient}}$
i	$20\times 20$	1 s	0.1478	0.1276
ii	$40\times 40$	0.5 s	0.0668	0.0561
iii	$80\times 80$	0.25 s	0.0224	0.0284

, 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 $\alpha$ in the transient analysis is chosen as $0.5$. The solution obtained using a $120\times 120$ mesh and a time step of $0.125\mathrm{s}$ is used as the reference solution. In Table 1, the error norm ${\mathrm{E}}_{\mathrm{steady}}$ is calculated on the basis of the computed steady-state velocity component $u_x$ along the midplane $x=0.5$ as

$$E_{\mathrm{steady}}=\sqrt{\sum |{\mathit{u}}_{\mathrm{reference}}-\mathit{u}{|}_{x=0.5}^2},$$

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.

Figure 3 presents the variation of the velocity components $u_x$ and $u_y$ as a function of y at the midplane $x=0.5\mathrm{m}$, whereas Figure 4 presents the streamlines corresponding to the steady-state solution.

3.2. 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 $-\rho G$, and a magnetic field perpendicular to the two plates $\mathit{H}=(B_o/{\mu }_m){\mathit{e}}_y$. The analytical solution for an incompressible fluid is given by [33].

$$\begin{array}{cc}\hfill u_x& =V\left(1-\frac{cosh(Ha\eta )}{cosh(Ha)}\right),u_y=0,u_z=0,p=-\rho Gx-\frac{\mu H_x^2}{2},\hfill \\ \hfill H_x& =\frac{B_o{\mathrm{Re}}_{\mathrm{m}}sinh(Ha\eta )}{{\mu }_{m}Hacosh(Ha)}-\left(1+\frac{E}{B_{o}V}\right)\frac{B_o{\mathrm{Re}}_{\mathrm{m}}\eta }{{\mu }_m},H_y=\frac{B_o}{{\mu }_m},H_z=0,\hfill \end{array}$$

(10)

where ${\mathrm{Re}}_{\mathrm{m}}={\mu }_m{\sigma }_{e}Vh$, $Ha=B_{o}h\sqrt{{\sigma }_e/\mu }$ and $\eta =y/h$.

Depending on whether the walls are insulating or conducting, the constants V and E in Equation (10) are given by

$$\begin{array}{cccc}\hfill V& =\frac{\rho GHa}{{\sigma }_e{B}_o^{2}tanh(Ha)},\hfill & \hfill E& =\frac{\rho G}{\sigma B_o}\left(1-\frac{Ha}{tanh(Ha)}\right)(\mathrm{Insulating}\mathrm{walls}),\hfill \\ \hfill V& =\frac{\rho G}{{\sigma }_e{B}_o^2},\hfill & \hfill E& =0(\mathrm{Perfectly}\mathrm{conducting}\mathrm{walls}).\hfill \end{array}$$

Since the above analytical solution is for an incompressible fluid, we use water as the fluid, with the equation of state as given in [34]. The data used are the same as in [5,29], namely ${\mu }_r=9.2878$, ${\sigma }_e=7.14\times {10}^5{\left(\mathsf{\Omega }-\mathrm{m}\right)}^{-1}$, $\mu =1.5\times {10}^{-4}\mathrm{kg}/\mathrm{ms}$, $\rho G=4.85\times {10}^{-5}\mathrm{N}/{\mathrm{m}}^3$, $B_o=1.4494\times {10}^{-5}\mathrm{Wb}/{\mathrm{m}}^2$, which correspond to $Ha=5$.

We use the following boundary conditions in our simulation, which correspond to the perfectly insulating case: At $y=\pm 0.5$, we have $\mathit{u}=0$, $H_x=0$, and $T=293.15\mathrm{K}$, whereas at $x=0,1$, the traction corresponding to a uniform pressure, and $H_y=B_0/{\mu }_m$, are imposed. We use a $2\times 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\times 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.

3.3. 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 ${\mathit{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]

$$\begin{array}{cc}\hfill U& =\frac{sinh(Hay)}{2sinh(Ha/2)}-y,\hfill \end{array}$$

(11)

where

$$\begin{array}{cc}\hfill U& =\frac{{Ha}^2}{Gr}\left(\frac{u_{x}H}{\nu }\right),\hfill \\ \hfill \mathrm{Gr}& =\frac{g\beta }{{\nu }^2}\frac{(\Delta T)H^4}{L},\hfill \\ \hfill Ha& =\sqrt{\frac{{\sigma }_e}{\mu }}{B}_{0}H.\hfill \end{array}$$

The flow properties used for this problem are $T_h=30 °\mathrm{C}$, $T_c=10 °\mathrm{C}$, $C_v=10\mathrm{J}/\mathrm{kg}-\mathrm{K}$, $\mu ={10}^{-3}\mathrm{kg}/\mathrm{m}-\mathrm{s}$, $k_f=1\mathrm{W}/\mathrm{m}-\mathrm{K}$, $\beta =207\times {10}^{-6}{\mathrm{K}}^{-1}$, $\rho =998.2\mathrm{kg}/{\mathrm{m}}^3$, ${\mu }_r=7.9577\times {10}^3$, ${\sigma }_e=10{\left(\mathsf{\Omega }-\mathrm{m}\right)}^{-1}$, and $g=1.552083478\times {10}^{-4}\mathrm{m}/{\mathrm{s}}^2$, which correspond to $Pr=0.01$, $\mathrm{Gr}=2\times {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. Mesh specifications for the two-dimensional buoyancy-driven flow problem.

Ha	No. of Elements	No. of Nodes	No. of Degrees of Freedom
1	$40\times 20$	3321	15,101
10	$40\times 20$	3321	15,101
100	$100\times 75$	30,351	142,176
1000	$200\times 200$	160,801	759,201

.

Figure 7 presents a comparison between the analytical and numerical results for a large range of Ha values (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. Figure 8 and Figure 9 present, respectively, the isotherms and streamlines for $Ha=1$, $Ha=10$, and $Ha=100$.

3.4. 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\mathrm{mm}\times 30\mathrm{mm}\times 30\mathrm{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\mathrm{K}$, $T_c=278.15\mathrm{K}$, $\mu ={10}^{-3}\mathrm{kg}/\mathrm{m}-\mathrm{s}$, $\nu ={10}^{-6}{\mathrm{m}}^2/\mathrm{s}$, $\beta =207\times {10}^{-6}{\mathrm{K}}^{-1}$, $C_v=756\mathrm{J}/\mathrm{kg}-\mathrm{K}$, ${\kappa }_f=31.5\mathrm{W}/{\mathrm{m}}^2-\mathrm{K}$, ${\mu }_r=7.9577\times {10}^3$, ${\sigma }_e=3.85\times {10}^6{\left(\mathsf{\Omega }-\mathrm{m}\right)}^{-1}$, $\mathrm{Pr}=0.024$. The values of g and $B_0$ for different values of Ha are given in

Table 3. Properties used for different Ha in the 3D natural convection problem.

Ha	${\mathit{B}}_0$ (Tesla)	$\mathit{g}(\mathbf{m}/{\mathbf{s}}^2)$
0	0	9.81
92	0.049423808	20
134	0.074135713	20
460	0.247119042	20

, where

$$\begin{array}{cc}\hfill \mathrm{Gr}& =\frac{g\beta (T_h-T_c)D^3}{{\nu }^2},{\mathrm{Gr}}^{*}=\frac{g\beta }{{\nu }^2{\kappa }_f}\frac{Q_{\mathrm{net}}{D}^4}{A},Ha=\sqrt{\frac{{\sigma }_e}{\mu }}{B}_{0}D,\hfill \\ \hfill \mathrm{Nu}& =\frac{Q_{\mathrm{net}}D}{A{\kappa }_f(T_h-T_c)},B_0={\mu }_m{H}_x,{\mu }_m={\mu }_r{\mu }_0,{\mu }_0=4\pi \times {10}^{-7}.\hfill \end{array}$$

$Q_{\mathrm{net}}$ is the net rate of heat transfer $\left(\mathrm{W}\right)$, D is the side of the cube $\left(30\mathrm{mm}\right)$, and A is the cross-sectional area $\left(900{\mathrm{mm}}^2\right)$. The values are chosen such that the modified Grashoff number ${\mathrm{Gr}}^{*}$ varies from ${10}^6$ to $4.24\times {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 ${\mathrm{Gr}}^{*}$ and Ha.

4. 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.