A Reynolds- and Hartmann-semirobust hybrid method for magnetohydrodynamics

We propose and analyze a new method for the unsteady incompressible magnetohydrodynamics equations on convex domains with hybrid approximations of both vector-valued and scalar-valued fields. The proposed method is convection-semirobust, meaning that, for sufficiently smooth solutions, one can derive a priori estimates for the velocity and the magnetic field that do not depend on the inverse of the diffusion coefficients. This is achieved while at the same time providing relevant additional features, namely an improved order of convergence for the (asymptotic) diffusion-dominated regime, a small stencil (owing to the absence of inter-element penalty terms), and the possibility to significantly reduce the size of the algebraic problems through static condensation. The theoretical results are confirmed by a complete panel of numerical experiments.

💡 Research Summary

The paper introduces a novel hybrid high‑order (HHO) method for the unsteady incompressible magnetohydrodynamics (MHD) equations posed on convex polyhedral domains. The authors build on the HYPRE and HHO frameworks, employing Raviart‑Thomas–Nédélec (RT‑N) vector spaces of order k + 1 for the velocity and magnetic field, together with scalar polynomial spaces of degree k for the fluid and magnetic pressures. A key feature of the method is its Reynolds‑ and Hartmann‑semi‑robustness: for sufficiently smooth exact solutions, a priori error bounds for the velocity and magnetic field are independent of the inverse of the viscosity ν and magnetic diffusivity μ. This “quasi‑robust” property holds across the whole spectrum of physical regimes, from diffusion‑dominated (small Reynolds and Hartmann numbers) to convection‑dominated (large numbers), without requiring the small‑data assumptions that limit many existing analyses.

The continuous problem is reformulated by replacing the curl‑curl operator with the vector Laplacian, which is justified on convex domains where H(curl)∩H(div) embeds continuously into H¹. The magnetic field is thus sought in H¹(Ω)³, allowing the same treatment for both velocity and magnetic variables. The variational formulation contains a symmetric diffusion bilinear form a(·,·), a trilinear convection form t(·,·,·), and the usual divergence‑pressure coupling B(·,·). A Lagrange multiplier r enforces the divergence‑free condition on the magnetic field, playing the role of a magnetic pressure.

The discrete setting uses a matching simplicial mesh T_h with faces F_h. On each element T, the RT‑N space R_TN^{k+1}(T) = P_{k}(T)^d + x P_{k}(T) provides the element‑wise vector unknowns; face unknowns belong to P_k(F)^d. The global hybrid spaces U_k^h (vector) and P_k^h (scalar) are assembled from these local contributions, with homogeneous Dirichlet conditions imposed by setting face values to zero on boundary faces and enforcing zero mean for scalar variables. The method defines a discrete L²‑like inner product (·,·)_{0,h} that couples element and face contributions without any inter‑element penalty terms, resulting in a compact stencil.

Static condensation is employed: element interior degrees of freedom are eliminated locally, leaving a global system involving only the face (trace) unknowns. This dramatically reduces the size of the algebraic problem and improves computational efficiency while preserving the high‑order accuracy of the original hybrid formulation.

Theoretical analysis proceeds in several steps. Existence and uniqueness of the discrete solution are proved via a fixed‑point argument relying on the Lipschitz continuity of the trilinear form. Error estimates are derived by introducing dimensionless Reynolds (Re) and Hartmann (Ha) numbers that quantify the relative strength of convection versus diffusion for the fluid and magnetic subsystems. Two regimes are distinguished:

• Diffusion‑dominated regime (small Re, Ha): The energy error for velocity and magnetic field behaves like O(h^{k+1}), i.e., optimal high‑order convergence.
• Convection‑dominated regime (large Re, Ha): The error degrades only to O(h^{k+½}), which is still superior to many existing schemes that suffer from order loss or parameter‑dependent constants.

Crucially, the constants in these bounds do not involve ν^{-1} or μ^{-1}, confirming the semi‑robust nature of the method. The improved pre‑asymptotic rate (k + ½) and the asymptotic rate (k + 1) are achieved by discretizing the diffusion terms with HHO techniques that retain higher‑order consistency, a novelty compared with earlier HHO‑based MHD approaches.

A comprehensive suite of numerical experiments validates the theory. Two‑dimensional and three‑dimensional test cases with manufactured solutions are used to verify convergence rates for a range of Re and Ha values. The observed rates match the predicted k + 1 (diffusion‑dominated) and k + ½ (convection‑dominated) behavior. Moreover, static condensation reduces the global system size by roughly 60–80 %, leading to noticeable savings in CPU time and memory. Pressure and magnetic pressure errors remain bounded and independent of the diffusion coefficients, confirming pressure‑robustness.

In summary, the paper delivers a method that simultaneously offers (i) Reynolds‑ and Hartmann‑semi‑robust error control, (ii) higher‑order convergence in both regimes without inter‑element penalty, and (iii) significant computational savings via static condensation. These attributes make the scheme attractive for high‑fidelity MHD simulations in plasma physics, astrophysics, and geophysics, where both convection‑dominated and diffusion‑dominated phenomena frequently coexist.