Hybridizable discontinuous Galerkin methods for solving the two-fluid plasma model

TL;DR

This work addresses the numerical stiffness of the coupled ion/electron fluid and Maxwell system by employing Hybridizable Discontinuous Galerkin methods to enable stable, high-order implicit time stepping. It provides a detailed HDG formulation, including problem decomposition, face-oriented fluxes, and automatic Jacobian generation, validated through convergence studies and large-scale two-fluid plasma simulations. The key contributions are the development of a scalable HDG approach with Schur complement condensation for a 42-equation, 3D-like system, and the demonstration that implicit HDG timestepping can outperform explicit schemes by enabling significantly larger time steps while preserving accuracy. The practical impact lies in enabling robust, high-fidelity simulations of multi-scale plasma dynamics that would be prohibitively expensive with traditional DG methods or explicit time stepping.

Abstract

The two-fluid plasma model has a wide range of timescales which must all be numerically resolved regardless of the timescale on which plasma dynamics occurs. The answer to solving numerically stiff systems is generally to utilize unconditionally stable implicit time advance methods. Hybridizable discontinuous Galerkin (HDG) methods have emerged as a powerful tool for solving stiff partial differential equations. The HDG framework combines the advantages of the discontinuous Galerkin (DG) method, such as high-order accuracy and flexibility in handling mixed hyperbolic/parabolic PDEs with the advantage of classical continuous finite element methods for constructing small numerically stable global systems which can be solved implicitly. In this research we quantify the numerical stability conditions for the two-fluid equations and demonstrate how HDG can be used to avoid the strict stability requirements while maintaining high order accurate results.

Figures (9)

• Figure 1: Example RKDG two-fluid linear stability eigencurves for $h = 10^{-1}$. Eigenvalues of $\stackrel{\leftrightarrow}{\boldsymbol{J}}_{S}$ are plotted as points, and are guaranteed to be pure imaginary. The DG advection stability requirements of $\stackrel{\leftrightarrow}{\boldsymbol{J}}_{A}$ grows quadratically with polynomial order while the eigenvalues of $\stackrel{\leftrightarrow}{\boldsymbol{J}}_{S}$ are constant with respect to polynomial order leading to a transitions from $\stackrel{\leftrightarrow}{\boldsymbol{J}}_{S}$ stability limited to $\stackrel{\leftrightarrow}{\boldsymbol{J}}_{A}$ stability limited by increasing the polynomial order from zero to threewarburton2008chalmers2020.
• Figure 1: Example nodal degrees of freedom in HDG finite element spaces. Black nodes correspond to element interiors projected onto $W_{h}$ and red nodes correspond to the skeleton space projected onto $M_{h}$.
• Figure 1: L-2 error norms for the linear advection equation converge at the expected optimal $O(N+1)$ rates for degree $N$ basis. The observed $10^{-5}$ minimum error is from the temporal discretization error.
• Figure 1: Two-fluid Brio-Wu shocktube properties at $t=100$. Similar features to those found in Shumlak and Loverichshumlak2003 are found, which are notably distinct from the classical MHD results.
• Figure 2: Global and face-oriented coordinate systems