Decoupled structure-preserving discretization of incompressible MHD equations with general boundary conditions

TL;DR

This work develops a mixed finite element, structure-preserving discretization for the incompressible MHD equations with general boundary conditions and a leapfrog-type temporal decoupling between the fluid and Maxwell parts, achieving second-order time accuracy and optimal spatial convergence. The semi-discrete formulation preserves mass and charge exactly and weakly enforces Gauss's law for magnetism, while the energy balance is maintained only in ideal or fully coupled limits; the decoupled scheme partially linearizes the Maxwell equations and confines nonlinear iterations to the fluid system. Numerical tests, including manufactured solutions, the Orszag–Tang vortex, and a lid-driven cavity, validate convergence rates and the preservation of key invariants, demonstrating computational efficiency and physical fidelity. The approach provides a practical framework for robust MHD simulations with general boundary conditions and offers avenues to extend invariants preservation (e.g., cross-helicity, magnetic helicity) in future work.

Abstract

In the framework of a mixed finite element method, a structure-preserving formulation for incompressible magnetohydrodynamic (MHD) equations with general boundary conditions is proposed. A leapfrog-type temporal scheme fully decouples the fluid part from the Maxwell part by means of staggered discrete time sequences and, in doing so, partially linearizes the system. Conservation and dissipation properties of the formulation before and after the decoupling are analyzed. We demonstrate optimal spatial and second-order temporal accuracy, as well as conservation and dissipation properties, of the proposed method using manufactured solutions, and apply it to the benchmark Orszag-Tang and lid-driven cavity cases.

Figures (8)

• Figure 1: An illustration of the leapfrog temporal scheme for the decoupled formulation. The time-step sequence is $\hat{s}^{0}\to S^{1}\to\hat{S}^1\to S^2\to \hat{S}^2\to\cdots$. Steps $S^{1}, S^{2},S^3, \cdots$ refer to the step 1, see \ref{['Eq: wf step 1']}, of the decoupled formulation, and steps $\hat{S}^{1}, \hat{S}^{2},\hat{S}^3, \cdots$ refer to the step 2, see \ref{['Eq: wf step 2']}, of the decoupled formulation. The pre-step $\hat{s}^{0}$ computes $\boldsymbol{H}_{h}^{\frac{1}{2}}$. A similar scheme has been used in a dual-field discretization of incompressible Navier-Stokes equations, see ZHANG2022110868.
• Figure 2: Results of $\left\| \boldsymbol{u}^k_h \right\|_{H(\mathrm{div})\text{-error}}$ and $\left\| \boldsymbol{H} ^{k+\frac{1}{2}}_h \right\|_{H(\mathrm{curl})\text{-error}}$ for the temporal convergence tests at $N=3$, $K=16$, $\Delta t \in \left\lbrace\frac{1}{5},\frac{1}{6},\cdots,\frac{1}{10},\frac{1}{12},\cdots,\frac{1}{36},\frac{1}{40},\cdots,\frac{1}{56}\right\rbrace$.
• Figure 3: Results of $\left\| \boldsymbol{u}^k_h \right\|_{H(\mathrm{div})\text{-error}}$, $\left\| P^{k-\frac{1}{2}}_h\right\|_{L^2\text{-error}}$, $\left\| \boldsymbol{\omega}^k_h \right\|_{H(\mathrm{curl})\text{-error}}$ and $\left\| \boldsymbol{H} ^{k+\frac{1}{2}}_h \right\|_{H(\mathrm{curl})\text{-error}}$ for the spatial convergence tests. We use $K\in\left\lbrace12,14,\cdots,32\right\rbrace$ for $N=1$, $K\in\left\lbrace10,12,\cdots,20\right\rbrace$ for $N=2$ and $K\in\left\lbrace8,10,\cdots,14\right\rbrace$ for $N=3$. And $\Delta t = \frac{1}{100}$ is employed to avoid pollution of the temporal discretization error, see Fig. \ref{['fig:temporal convergence tests']}.
• Figure 4: Results of $\left\|\nabla\cdot\boldsymbol{H}^{k+\frac{1}{2}}_h\right\|_{L^2}$ for the spatial convergence tests. We use $K\in\left\lbrace10,12,\cdots,20\right\rbrace$ for $N=2$, $K\in\left\lbrace8,10,\cdots,14\right\rbrace$ for $N=3$, and $\Delta t = \frac{1}{100}$. In this work, the employed finite-dimensional space for $\boldsymbol{H}_h \in C(\Omega)$ is the mimetic spectral space which is a polynomial space $\mathcal{P}^{N-1,N,N}\times\mathcal{P}^{N,N-1,N}\times\mathcal{P}^{N,N,N-1}$ on orthogonal meshes. So we can compute exact divergence of $\boldsymbol{H}_h$ element-wise, and, when $N=1$, $\nabla\cdot\boldsymbol{H}_h\equiv 0$. For the present results, $\left\|\nabla\cdot\boldsymbol{H}^{k+\frac{1}{2}}_h\right\|_{L^2}:=\sqrt{\sum_i H_i^2}$ where $H_i$ is the $L^2$-norm of $\nabla\cdot\boldsymbol{H}_h^{k+\frac{1}{2}}$ in the $i$th element.
• Figure 5: Some results of $\left\|\nabla\cdot \boldsymbol{u}_{h}^{k}\right\|_{L^2}$, $\tilde{\mathcal{E}}_{h}^k$, $\left|\dfrac{\tilde{\mathcal{E}}^k_h-\tilde{\mathcal{E}}^{k-1}_h}{\Delta t} - \mathcal{D}^{k+\frac{1}{2}}_h\right|$ over time for the conservation and dissipation tests at $\mathsf{c} = 1$, $N=2$, $h=\frac{1}{8}$ and $\Delta t=\frac{1}{50}$, where $\mathcal{D}_{h}^{k-\frac{1}{2}} := - \mathrm{R}_f^{-1}\mathcal{S}_{h}^{k-\frac{1}{2}} - \mathsf{c}\mathrm{R}_m^{-1}\tilde{\mathcal{J}}_{h}^{k-\frac{1}{2}} + \mathsf{c}\left( \mathcal{A}_{h}^{k-\frac{1}{2}} -\tilde{\mathcal{A}}^{k-\frac{1}{2}}_{h}\right)$, cf. \ref{['Eq: discrete energy dissipation']}. The discrete energy $\mathcal{E}^k_h + (\mathcal{E}^k_h-\mathcal{E}^0_h)\times10^{10}$ ($\color{blue}---$) for the ideal case of a Crank-Nicolson discretization of the coupled formulation (see \ref{['App: coupled discretization']}) is also shown in the top part of the bottom-left diagram.