Energy-consistent finite difference schemes for compressible hydrodynamics and magnetohydrodynamics using nonlinear filtering
aa r X i v : . [ phy s i c s . c o m p - ph ] F e b Energy-consistent finite difference schemesfor compressible hydrodynamics and magnetohydrodynamicsusing nonlinear filtering
Haruhisa Iijima
Institute for Space-Earth Environmental Research, Nagoya University, Furocho, Chikusa-ku, Nagoya, Aichi464-8601, JapanNational Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Abstract
In this paper, an energy-consistent finite difference scheme for the compressible hydrody-namic and magnetohydrodynamic (MHD) equations is introduced. For the compressiblemagnetohydrodynamics, an energy-consistent finite difference formulation is derived usingthe product rule for the spatial difference. The conservation properties of the internal, ki-netic, and magnetic energy equations can be satisfied in the discrete level without explicitlysolving the total energy equation. The shock waves and discontinuities in the numerical solu-tion are stabilized by nonlinear filtering schemes. An energy-consistent discretization of thefiltering schemes is also derived by introducing the viscous and resistive heating rates. Theresulting energy-consistent formulation can be implemented with the various kinds of centraldifference, nonlinear filtering, and time integration schemes. The second- and fifth-orderschemes are implemented based on the proposed formulation. The conservation propertiesand the robustness of the present schemes are demonstrated via one- and two-dimensionalnumerical tests. The proposed schemes successfully handle the most stringent problems inextremely high Mach number and low beta conditions.
Keywords:
Hydrodynamics, Magnetohydrodynamics, Computational plasma physics,Finite differences, Shock capturing, Spatial filtering, Skew-symmetric form, Secondaryconservative
1. Introduction
The fully conservative finite difference scheme has been used as a tool for the numericalsimulation of the incompressible flow due to its very small numerical dissipations and thestability in the long-term simulations [1, 2, 3]. In such schemes, the quadratic split form (alsoknown as the skew-symmetric split form) of the convective term is used to achieve the kinetic
Email address: [email protected] ( Haruhisa Iijima )
Preprint submitted to Journal of Computational Physics February 26, 2021 nergy conservation in the discrete sense while maintaining the momentum conservationwithout directly solving the kinetic energy equation. The quadratic and cubic split operatorscan also reduce the aliasing error produced in the product of the multiple variables [4, 5].Several researchers have attempted to extend the split form of the convective term tothe low Mach number flows (see a review by Pirozzoli [6]). The Jameson-Schmidt-Turkelscheme [7] implicitly used the quadratic split operator in the finite volume formulation[8]. The kinetic energy preserving (KEP) scheme is constructed using the quadratic splitoperator [9, 10, 11]. The quadratic split operator can even achieve the consistency betweenthe internal energy and kinetic energy in the discrete sense [12, 13]. The concept of thequadratic split operator can be extended into the nonuniform cylindrical coordinates [14].The temporal derivatives can be also discretized using the quadratic split operator to obtainthe spatiotemporal conservation properties [13]. The conservation of the entropy can alsobe achieved in addition to the energy-consistency [15].Fewer studies investigated the application of the energy-consistent schemes to the highMach number flows. In high-speed flow simulations, additional numerical diffusion schemesare required to stabilize the shock waves and discontinuities. Jameson et al. [7] used anartificial viscosity to stabilize their skew-symmetric finite volume scheme. Ducros et al. [8]suggested the hybridization between the central and upwind schemes. Yee [16] employeda nonlinear filtering approach for the entropy conserving scheme. However, the energyconsistency of these diffusion schemes was rarely discussed explicitly, except in the work byShiroto et al. [17].In recent years, the application of the quadratic split operator has increasingly beenapplied to the plasma physics, to achieve the secondary conservation properties. Ni etal. [18, 19, 20] constructed the current density conservative scheme for the incompressiblemagnetohydrodynamic (MHD) simulations. Shiroto et al. [17] applied the approach tothe two-temperature plasma flows. There are several applications of the quadratic splitoperator on the kinetic plasma equations [21, 22]. However, an energy-consistent scheme forapplication to compressible MHD simulations has not yet been constructed.In the compressible MHD simulations, the consistency between the internal, kinetic, andmagnetic energies is essential to accurately reproduce the plasma dynamics and energetics.Many numerical schemes solve the total (internal + kinetic + magnetic) energy equation tosatisfy the jump condition near the shocks and discontinuities. However, in such schemes, thediscretization errors of the kinetic and magnetic energy equations are imposed on the internalenergy. As a result, the evolution of the internal energy becomes inevitably erroneous whenthe internal energy is much smaller than the kinetic or magnetic energy. Such numericalsimulations yield energetically inaccurate solutions, such as the negative pressure.This study aims to construct an energy-consistent formulation for the compressible hy-drodynamic and MHD equations while taking account the nonlinear filtering. Hence, an fullyenergy-consistent finite difference scheme is devised. Our implementation of this scheme withthe appropriate nonlinear filtering scheme provides superior numerical robustness, especiallyin the extremely high Mach number and low plasma beta (internal over magnetic energies)conditions.This paper is organized as follows. Section 2 presents the energy-consistent formulations2f the governing equations as well as the nonlinear filtering flux. Note that the proposedformulation is independent from the details of the finite difference operators, the filteringflux, and the temporal integration method. Section 3 describes our implementation of theproposed scheme. In Section 4, several numerical tests are performed, especially focusing onthe energy-consistency and numerical robustness of the proposed scheme. Finally, Section5 presents the conclusions of this article.
2. Energy-consistent scheme for compressible hydrodynamic and MHD equa-tions
The basic equations of the ideal, fully compressible MHD system in the Cartesian coor-dinates x i ( i = 1 , ,
3) can be written in a conservative form as ∂ρ∂t + ∂ρV j ∂x j = 0 , (1) ∂ρV i ∂t + ∂ρV j V i ∂x j + ∂P∂x i + 12 ∂B j B j ∂x i − ∂B i B j ∂x j = 0 , (2) ∂∂t (cid:18) e + 12 ρV i V i + 12 B i B i (cid:19) + ∂∂x j (cid:20)(cid:18) e + P + 12 ρV i V i + B i B i (cid:19) V j − B j B i V i (cid:21) = 0 , (3) ∂B i ∂t + ∂ ( V j B i − V i B j ) ∂x j = 0 . (4)Here, the summation rule is assumed. The variables ρ , e , and P represent the mass density,the internal energy density, and the gas pressure, respectively, and V i and B i represent the x i -component of the velocity field and the magnetic field, respectively. The absence of themagnetic monopole provides a constraint on the magnetic field: ∂B j ∂x j = 0 . (5)The above equations are closed by the equation of states for the gas pressure P ( ρ, e ). Thehydrodynamic equations can be derived by neglecting all components of the magnetic field(setting B i = 0 for any i ).The conservation law of the momentum expressed in Eqs. (2) can be rewritten as ∂ρV i ∂t + ∂ρV j V i ∂x j = F Pi + F Li , (6)where F Pi is the pressure gradient force in the i -th direction F Pi = − ∂P∂x i (7)3nd F Li is the Lorentz force in the i -th direction F L,ci = − B j (cid:18) ∂B j ∂x i − ∂B i ∂x j (cid:19) + B i ∂B j ∂x j , (8)respectively. The last term of Eq. (8) disappears when the solenoidal condition of Eq. (5)is satisfied. The resultant Lorentz force can be written as F Li = − B j (cid:18) ∂B j ∂x i − ∂B i ∂x j (cid:19) . (9)In this form, the Lorentz force is normal to the direction of the magnetic field ( F Li B i = 0).In numerical simulations, the solenoidal constraint of the magnetic field is sometimesviolated due to discretization error. Violation of the ∇· B = 0 condition produces theinconsistency between the conservative and non-conservative Lorentz force, which yieldsartificial field-aligned forces B i ∇· B .From the basic equations (Eqs. (1)-(4)), three independent energy equations, namely,internal, kinetic, and magnetic energy equations, can be derived: ∂e∂t + ∂ ( e + P ) V j ∂x j = − W P , (10) ∂ ( ρV i V i / ∂t + ∂ ( ρV i V i V j / ∂x j = W P + W L , (11) ∂ ( B i B i / ∂t + ∂ ( B i B i V j − B j B i V i ) ∂x j = − W L , (12)respectively. Here, we defined the work done by the pressure gradient force W P = V i F Pi (13)and the work done by the Lorentz force W L = V i F Li . (14)It is clear that the equations for the internal, kinetic, and magnetic energies (Eqs. (10)–(12)) consist of terms pertaining to the energy transport (i.e., the divergence of the enthalpy,kinetic energy, and Poynting fluxes) and the works done by the pressure gradient and Lorentzforces ( W Pi and W Li ). Summation of the three energy equations retrieves the law of the totalenergy conservation expressed in Eq. (3).The momentum and magnetic flux conservation laws are solved in many numericalschemes that employ the divergence formulation of the MHD equations. However, the in-ternal, kinetic, and magnetic energy equations are not solved as independent equations and,thus, the consistency among them is sometimes violated. In stringent problems, violationof the energy-consistency sometimes yields a negative pressure.4 .2. Discrete operators This work essentially follows the notation for the finite difference operators given inMorinishi et al. [2]. The finite difference operator with stencil n acting on a variable Φ withrespect to x is defined by δ n Φ δ n x ≡ Φ ( x + nh / , x , x ) − Φ ( x − nh / , x , x ) nh , (15)where h is the grid spacing in the x -direction. The interpolation operator with stencil n with respect to x is defined asΦ nx ≡ Φ ( x + nh / , x , x ) + Φ ( x − nh / , x , x )2 . (16)In addition, the special interpolation operator of the product of two variables Φ and Ψ withstencil n with respect to x is defined asΦΨ : nx ≡
12 Φ ( x + nh / , x , x ) Ψ ( x − nh / , x , x )+ 12 Φ ( x − nh / , x , x ) Ψ ( x + nh / , x , x ) . (17)These operators are similarly defined in the x - and x -directions.The finite difference operator defined by Eq. (15) is second-order accurate in space. Thehigher-order finite difference operator can be represented by the weighted average of the n -stencil operator as δ Φ δx ≡ N c X n =1 c n δ n Φ δ n x , (18)where the coefficients c n satisfy P n c n = 1. For example, the standard five-stencil fourth-order accurate central difference operator ( N c = 2, c = 4 /
3, and c = − /
3) can beexpressed as δ Φ δx = 43 δ Φ δ x − δ Φ δ x , (19)as shown in Morinishi et al. [2]. These operators are similarly defined in other directions.In Sec. 3.1, two specific implementation of the operator δ/δx j are provided. We use theoperator δ/δx j as a building block of our formulation to simplify the notation.The n -stencil finite difference operator satisfies the following three product rules:Φ δ n Ψ δ n x j + Ψ δ n Φ δ n x j = δ n ΦΨ : nx j δ n x j , (20)Φ nx j δ n Ψ δ n x j + Ψ nx j δ n Φ δ n x j = δ n ΦΨ δ n x j , (21)5nd Ψ δ n Φ δ n x j nx j + Φ δ n Ψ δ n x j = δ n Φ nx j Ψ δ n x j . (22)The right-hand sides of these equations are conservative. In this study, these product rulesare frequently used to analyze the conservation properties of the finite difference scheme.The conservation property of the higher-order operator also satisfies similar relations, suchas Φ δ Ψ δx j + Ψ δ Φ δx j = N c X n =1 c n δ n ΦΨ : nx j δ n x j . (23)It is clear that the right-hand side of the equation is conservative.For future reference, we note the commutation rule of the finite difference operatorexpressed as δ n δ n x i (cid:18) δ n Φ δ n x j (cid:19) = δ n δ n x j (cid:18) δ n Φ δ n x i (cid:19) . (24)The higher-order operator also satisfies this rule. We propose an energy-consistent formulation of the fully compressible magnetohydrody-namic equations expressed as: ∂ρ∂t + δρV j δx j = 0 , (25) ∂ρV i ∂t + 12 δρV j V i δx j + 12 ρV j δV i δx j + 12 V i δρV j δx j = − δPδx i − B j (cid:18) δB j δx i − δB i δx j (cid:19) , (26) ∂e∂t + δ ( e + P ) V j δx j = V j δPδx j , (27) ∂B i ∂t + δ ( V j B i − V i B j ) δx j = 0 . (28)Except the energy equation, the hydrodynamic component of this formulation is similar tothe second scheme proposed by Kok [12]. In Kok’s scheme, the divergence of the enthalpyflux is expressed using the quadratic split operator. We use the simple divergence form forthe enthalpy flux in the internal energy equation for the consistency with the divergence ofthe mass flux in Eq. (25) and the curl of the electric field in Eqs. (28). One advantage ofour formulation is that the pressure perturbation is not excited at the contact discontinuitywhere P , e , and V j are constant in space. The energy conservation property is not affectedby this change. Note that the entropy conservation can also be achieved by changing thediscretization of the equation of the continuity and the internal energy as shown by Kuyaet al. [15]. Given that the purpose of this paper is to demonstrate the robustness of theenergy-consistent scheme in stringent problems, the details of various formulations for thehydrodynamic equations are not discussed. 6 .4. Momentum conservation in the proposed formulation The convective term in the equations of motion (Eqs. (26))conserves the volume-averagedmomentum as shown in the previous studies [12, 13]. Using the product rule (Eq. (23)), theconvective term can be rewritten as12 δρV j V i δx j + 12 ρV j δV i δx j + 12 V i δρV j δx j = 12 δρV j V i δx j + N c X n =1 c n δ n ( ρV j ) V i : nx j δ n x j . (29)The right-hand side is conservative, and the convective term conserves momentum.The Lorentz force in the equations of motion can be also rewritten in the conservativeform as F L i = − B j (cid:18) δB j δx i − δB i δx j (cid:19) = − N c X n =1 c n δ n B j B j : nx i δ n x i − δ n B i B j : nx j δ n x j ! − B i δB j δx j . (30)The last term disappears when the numerical solenoidal condition δB j δx j = 0 (31)is satisfied, i.e., the magnetic field divergence is zero.By removing the − B i ∇· B term in the right-hand side of the momentum equations (seeEq. (8)), an energy-consistent formulation of the MHD equations can be constructed, whichconserves momentum even when the numerical solenoidal condition is violated. However,this modification produces artificial forces parallel to the magnetic field lines as well asan additional heating term V i B i ∇· B in the internal energy equation to maintain the totalenergy conservation. Therefore, this field-aligned artificial force may yield an unphysicalsolution. To avoid this problem, the above form of the Lorentz force is selected in thisstudy. This subsection demonstrates the total energy conservation in the proposed formulation.As discussed in previous studies, the skew-symmetric convective term conserves the kineticenergy with the appropriate discretization of the mass conservation law. The convectiveterm multiplied by V i can be written as V i (cid:20) δρV j V i δx j + 12 ρV j δV i δx j + 12 V i δρV j δx j (cid:21) = 12 V i δρV j V i δx j + 12 ρV j V i δV i δx j + V i V i δρV j δx j = N c X n =1 c n δ n ( ρV j V i ) V i : nx j δ n x j + V i V i δρV j δx j . (32)7he kinetic energy equation is derived from the equations of motion Eqs. (26) and theequation of continuity Eq. (25) as ∂ρV i V i / ∂t = V i ∂ρV i ∂t − V i V i ∂ρ∂t = − N c X n =1 c n δ n ( ρV j V i ) V i : nx j δ n x j + W P + W L . (33)The first term of the last equation represents the divergence of the kinetic energy flux. Thesecond term of Eq. (32) is canceled with the equation of continuity.The work done by the pressure gradient force W P = − V j δPδx j (34)represents the interaction between the internal and kinetic energy equations. Apparently,this is exactly equivalent to the right-hand side of the internal energy equation Eq. (27),apart from the difference in sign.The magnetic energy equation is derived by multiplying the induction equation Eq. (28)by B i , such that ∂B i B i / ∂t = B i ∂B i ∂t = − B i δ ( V j B i − V i B j ) δx j = − N c X n =1 c n δ n B i ( V j B i − V i B j ) : nx i δ n x i − W L , (35)where the first term in the last equation is (minus of) the divergence of the Poynting fluxand the second term represents the work done by the Lorentz force W L = − V i B j (cid:18) δB j δx i − δB i δx j (cid:19) . (36)Comparison with the discretized Lorentz force of Eq. (30) reveals that the interactionbetween the kinetic and magnetic energy equations is consistent in the discretization level.Note that we did not use the solenoidal condition of Eq. (31) in the derivation of themagnetic energy equation. Thus, the proposed formulation conserves the total energy evenwhen the ∇· B = 0 condition is violated. Shock waves and discontinuities are often involved in the solutions of compressible hydro-dynamics and magnetohydrodynamics. In the approach adopted in this work, the numericalsimulation is stabilized by diffusing the sharp gradient through introducing of the nonlin-ear filtering flux in the semi-discrete formulation of Eqs. (25)–(28). An energy-consistentformulation is obtained even when this nonlinear filtering flux is introduced.8e design an energy-consistent formulation for the filtering flux defined at the center ofthe cell face. Note that many linear and nonlinear filtering schemes can be reduced to thisformulation. The filtering process of the variable U can be expressed as ∂U∂t = δ D j ( U ) δ x j , (37)where D j ( U ) is the filtering flux in the j -th direction defined at the cell edge. The typicalfirst-order filtering flux can be expressed as D j ( U ) = η ( U ) δ Uδ x j , (38)where η ( U ) is the diffusion coefficient of U . The construction of the higher-order filteringflux is described in Section 3. The energy-consistent property is not affected by the specificimplementation of the filtering flux D j ( U ) as long as the filtering flux is defined at the cellface. In this study, we apply this filtering scheme on the mass density ρ , the momentumdensity ρV i , the magnetic flux density B i , and the internal energy density e .When the filtering scheme is applied on the equations of motion or the induction equa-tions, the viscous heating Q vis or the resistive heating Q res should be included as heatingterms in the internal energy equation to maintain the total energy conservation. We proposethe energy-consistent discretization of the viscous and resistive heating rates correspondingto the filtering flux as Q vis = D j ( ρV i ) δ V i δ x j x j − D j ( ρ )2 δ V i V i δ x j x j (39)and Q res = D j ( B i ) δ B i δ x j x j , (40)respectively. These expressions can be rewritten as Q vis = δ δ x j (cid:20) V i x j D j ( ρV i ) − V i V i x j D j ( ρ ) (cid:21) − V i δ D j ( ρV i ) δ x j + V i V i δ D j ( ρ ) δ x j (41)and Q res = δ B i x j D j ( B i ) δ x j − B i δ D j ( B i ) δ x j . (42)Here, the first terms in these two equations are the kinetic and magnetic energy transportby the filtering fluxes, respectively. The other terms correspond to the time variation of thekinetic and magnetic energies produced by the filtering flux. Clearly, these heating ratesare energy-consistent in terms of the exchange among the internal, kinetic, and magneticenergies. 9he viscous and resistive heating rates defined in Eqs. (39) and (40) become positiveunder the specific condition. We assume the filtering flux of the momentum density ρV i evaluated from those of the mass density and velocity field as D j ( ρV i ) = ρ x j D j ( V i ) + V i x j D j ( ρ ) . (43)Then, the viscous heating rate Eq. (39) can be rewritten as Q vis = h ρ x j D j ( V i ) + V i x j D j ( ρ ) i δ V i δ x j x j − D j ( ρ )2 δ V i V i δ x j x j = ρ x j D j ( V i ) δ V i δ x j x j , (44)where the product rule of Eq. (21) for Φ = Ψ = V i is used. The above viscous heatingrate becomes positive if the sign of D j ( V i ) is the same as that of δ V i /δ x j . This conditionis often satisfied considering that the D j ( V i ) is the filtering flux acting to diffuse V i in the x j -direction. Similarly, the resistive heating rate Eq. (40) becomes positive when D j ( B i )and δ B i /δ x j have the same sign. In the multi-dimensional simulations of the magnetohydrodynamics, the solenoidal condi-tion ∇· B = 0 should be preserved throughout the time integration. Violating the solenoidalrule may produce artificial inconsistencies in the numerical solution.The central difference scheme is known to produce no numerical magnetic monopole [23]as δδx i (cid:18) ∂B i ∂t (cid:19) = − δδx i (cid:20) δ ( V j B i − V i B j ) δx j (cid:21) = − δδx i (cid:18) δV j B i δx j (cid:19) + δδx i (cid:18) δV i B j δx j (cid:19) = 0 . (45)Here, we use the commutation rule of the central difference operators (Eq. (24)). On theother hand, the filtering flux violates the solenoidal condition ∇· B = 0 as δδx k (cid:20) δ D j ( B i ) δ x j (cid:21) = 0 . (46)We need some treatment to control the magnetic monopole on the numerical solution.In this study, the hyperbolic/parabolic divergence cleaning method [24] is used to controlthe divergence of the magnetic field. The numerical magnetic monopole is transported anddiffused by introducing an additional equation of the generalized Lagrange multiplier (GLM) ψ as ∂ψ∂t = − c ψ δB i δx i − ψτ ψ , (47)10here the propagation speed c ψ and the damping rate τ ψ are the tunable parameters thatcontrol the efficiency of the divergence cleaning. Correspondingly, the induction equationsare corrected as ∂B i ∂t = − δψδx i , (48)where, for the sake of clarity, we dropped the terms containing the curl of the electric fieldand the filtering flux on the magnetic flux density. It is clear that the divergence cleaningmethod conserves the magnetic flux in a volume.Although the hyperbolic/parabolic divergence cleaning method can efficiently reduce thedivergence of the magnetic field, ∇· B still remains in the numerical solution. The numericaldiscretization should be energy-consistent even when the small numerical error of ∇· B exists.This method produces the additional term in the magnetic energy equation as ∂B i / ∂t = − B i δψδx i = − N c X n =1 c n δ n B i ψ : nx i δ n x i − Q ψ , (49)where we define the heating term due to the divergence cleaning method Q ψ = − ψ δB i δx i . (50)In the steady state, this term becomes positive, as Q ψ → τ ψ c ψ (cid:18) δB i δx i (cid:19) ≥ . (51)Introduction of Q ψ in the internal energy equation provides the energy-consistent formulationwith the hyperbolic/parabolic divergence cleaning. The energy-consistent formulation of the compressible MHD equations is summarized asfollows: ∂ρ∂t + δρV j δx j = δ D j ( ρ ) δ x j , (52) ∂ρV i ∂t + 12 δρV j V i δx j + 12 ρV j δV i δx j + 12 V i δρV j δx j + δPδx i + B j δB j δx i − B j δB i δx j = δ D j ( ρV i ) δ x j , (53) ∂e∂t + δ ( e + P ) V j δx j = − V j δPδx j + δ D j ( e ) δ x j + Q vis + Q res + Q ψ , (54) ∂B i ∂t + δ ( V j B i − V i B j ) δx j + δψδx i = δ D j ( B i ) δ x j , (55) ∂ψ∂t + c ψ δB j δx j = − ψτ ψ + δ D j ( ψ ) δ x j . (56)11he viscous and resistive heating rates ( Q vis and Q res ) are defined in Eqs. (39) and (40),respectively. The heating caused by the numerical magnetic monopole Q ψ is defined inEq. (50). In this formulation, the internal, kinetic, and magnetic energy equations areconsistent even if the nonlinear filtering flux and the divergence cleaning are introduced. Theenergy-consistent formulation of the compressible hydrodynamic equations can be derivedby setting B i , ψ , Q res , and Q ψ to be zero in Eqs. (52)–(54). The proposed formulationis independent of the specific implementations of the temporal discretization, the centraldifference operator δ/δx j , and the nonlinear filtering flux D j ( U ) for each variable U . Thedetails of our implementation are described in Section 3.
3. Numerical implementation
This section describes the details of our implementation of the proposed formulation. Inthe proposed formulation Eqs. (52)-(56), the time derivative is not discretized.
In this subsection, we describe our implementation and performance of the base schemesusing the one-dimensional scalar advection equation: ∂U∂t + a ∂U∂x = 0 , (57)where U is the scalar variable and a is the advection speed. The semi-discrete form of theabove equation can be written as ∂U∂t + a δUδx = δ D ( U ) δ x . (58)We implemented the second- and fifth-order schemes to show the general performance ofthe proposed formulation. The important feature of these schemes is the small numericalovershoot, which seems to contribute the high robustness of our formulation shown in Sec. 4.We speculate that the similar high robustness can be achieved with other combination of thecentral difference operator and filtering flux, as long as the produced numerical overshoot issufficiently small.In the second-order scheme, we used the second-order central difference operator δ/δx j defined in Eq. (18) by setting N c = 1 and c = 1. The second-order filtering flux was takenfrom the Jameson-Schmidt-Turkel (JST) scheme [7, 25]. The JST scheme is total variationdiminishing if the limiter function is chosen appropriately [26]. Among the variant of theJST scheme, we used the symmetric limited positive (SLIP) scheme [27], which is provento be positivity preserving for the scalar conservation law. The filtering flux of the SLIPscheme for the scalar advection equation can be written as D JST1 ( U ; x + h /
2) = | a | (cid:2) ∆ U x + h / − L (cid:0) ∆ U x +3 h / , ∆ U x − h / (cid:1)(cid:3) , (59)12here ∆ U x = U ( x + h / − U ( x − h /
2) is the difference of U between two adjacent cellsand L ( a, b ) is the limiter function. We used van Leer’s limiter (originaly introduced in [28],implemented as described by [27]) for the limiter function L .To demonstrate the applicability of the proposed scheme in the higher order, we alsoimplemented the fifth-order scheme using the sixth-order central difference operator δ/δx j for N c = 3, c = 3 / c = − /
5, and c = 1 /
10. The filtering flux is based on the fifth-orderweighted essentially non-oscillatory (WENO) reconstruction [29]. The WENO scheme wasoriginally developed in the context of the upwind scheme and is not intended to be usedas the filtering flux for the central difference method. The filtering flux in the x -directionbased on the WENO reconstruction can be written as D WENO1 ( U ; x + h /
2) = c diff | a | (cid:2) U R ( x + h / − U L ( x + h / (cid:3) , (60)where U R and U L are the upwind value of U reconstructed from the positive and negative x -direction, respectively. A positive diffusion parameter c diff of order unity is introduced tocontrol the numerical diffusion by the filtering flux.The one-dimensional scalar advection equation Eq. (57) conserves its energy U / ∂U / ∂t + a ∂U / ∂x = 0 . (61)The corresponding semi-discrete equation can be obtained by multiplying Eq. (58) by U as U ∂U∂t + a N c X n =1 c n δ n U / : nx δ n x − δ U x D ( U ) δ x = − D ( U ) δ Uδ x x , (62)where the right-hand side represents the energy loss by the filtering flux D ( U ). Notethat the similar terms of the energy loss in the kinetic and magnetic energy equations in ourMHD formulation are treated as the viscous and resistive heating rates in the internal energyequation to achieve the total energy conservation even if the filtering flux is introduced (seeSec. 2.6). The semi-discrete scheme Eq. (58) conserves the energy U / U ∂U/∂t = ∂ ( U / ∂t cannot be satisfied bya time integration scheme in general, like most of the Runge-Kutta methods. One simpleexception is the implicit midpoint method, which is symmetric in temporal direction. Similarto the advection equation, the total energy conservation in the semi-discrete formulation ofthe MHD equations in Eqs. (52)–(56) produces the conservation error of the total energy bythe time integration scheme. More detailed descriptions on the relation between the energyconservation and the temporal integration method can be found in [12, 13].In all test problems described in this paper, we used the third-order optimal stronglystability preserving (SSP) Runge-Kutta method [30] for the SSP property in the nonlinearequations [31] and its frequent use with the fifth-order WENO scheme [29]. We also foundthat the four-step second-order Runge-Kutta method by Jameson & Baker [32] also performs13ell even in the stringent problems in Secs. 4.5 and 4.6. Therefore, we speculate that therobustness of the proposed formulation may be weakly affected by the choice of the timeintegration scheme. Figure 1: Results of one-dimensional scalar advection problem using 200 grid points. Shown are the resultsat t = 1 for (a) JST scheme, (b) sixth-order central difference method with the filtering flux by the WENO5scheme with c diff = 2, and (c) sixth-order central difference method with the filtering flux by the WENO5scheme with c diff = 1. The solid lines indicate the analytical solution. The diamond symbols indicate thenumerical results. The Courant number of 0.3 is used. Figure 1 shows the solution of the one-dimensional scalar advection problem after oneperiod. The initial profile is identical to that of the test problem in Suresh & Huynh [33].The results show that both the second-order scheme with the JST filtering flux and the sixth-order scheme with the fifth-order WENO filtering flux efficiently preserve the monotonicityof the solution.The dependence on the diffusion parameter c diff of the WENO filter becomes clear nearthe discontinuity. Figure 2 shows solutions for the same problem as Fig. 1, but showing14 igure 2: One-dimensional scalar advection problem using 200 grid points with sixth-order central differencemethod with the filtering flux by the WENO5 scheme. Shown are the results one step after the integration(at t = 0 . c diff = 2 (blue dotted lines with plus) and c diff = 1 (red dashed lines with cross). Thesolid lines indicate the analytical solution. The Courant number of 0.3 is used. one step after the time integration (at t = 0 . c diff = 1. We found that the overshoot can be efficiently damped bydoubling the numerical diffusion in the WENO filter ( c diff = 2). However, small overshootof about 10 − still remains even if c diff = 2. We note that this value of c diff = 2 isnot guaranteed theoretically nor optimized for general problems. However, this choice wasfound to show high numerical stability in all test problems presented in Sec. 4. Becausethe formulation proposed in this study is independent of the specific implementation of thefiltering flux, better methods of evaluating the filtering flux, especially with less numericalovershoot and higher-order spatial convergence, should be investigated in future studies. For the primitive variables in the MHD equations with the hyperbolic/parabolic diver-gence cleaning method, W i = ( ρ, e, V , V , V , B , B , B , ψ ), the one-dimensional form of theevolution equations can be written as ∂W l ∂t + X n A ln ∂W n ∂x = 0 , (63)where A is the coefficient matrix for the primitive variables. The matrix A can be diagonal-ized as A ln = X m R lm λ m L mn , (64)where L and R are the left and right eigenmatrices [34, 35], respectively, and λ m indicatesthe eigenvalues of the matrix A .There are several ways to extend the filtering flux for the scalar advection equation tothe system equations. Here, we implemented two methods, namely, the LLF-type and the15oe-type filtering schemes. The first approach to evaluate the filtering flux is to diffuse allvariables with the same characteristic speed, as in the local Lax-Friedrichs (LLF) scheme[36, 37]. The filtering flux of the primitive variables in the x -direction at the location x + h / D ( W i ; x + h /
2) = a max D [ W ( x + h / ± h / , W ( x + h / ± h / , · · · ] , (65)where a max = max ( | λ | , ..., | λ | ) is the maximum phase speed of the waves in the systemand the discrete operator D represents the limiting process in the JST and WENO5 filteringfluxes. The filtering flux for the momentum density was computed from that of the massdensity and the velocity field using Eq. (43).The second filtering flux implemented in this study is based on Roe’s approximate Rie-mann solver [38, 26, 16], which is given by D ( W i ; x + h /
2) = 12 X m R lm | λ m | D [ Q m ( x + h / ± h / , Q m ( x + h / ± h / , · · · ] , (66)where Q m ( x ) = X n L mn W n ( x ) (67)is the characteristic variables with the phase speed of λ m . The left and right eigenmatrices R ij and L jk and the eigenvalues λ j are evaluated at the fixed location x + h /
2. In ourimplementation, these eigenmatrices and eigenvalues are derived from the arithmetic averageof the primitive variables at the adjacent grids, i.e., W i x ( x + h /
2) = ( W i ( x ) + W i ( x + h )) /
2. Although the original Roe’s scheme [38] employs the Roe-average, the arithmeticaverage has often been utilized in practical problems, and it seems to work well [39]. Wealso used the first entropy fix suggested by Harten & Hyman [40] to evaluate the eigenvalues.The filtering flux of the momentum density was derived using Eq. (43).In the Roe-type method, the characteristic variables are used for the filtered variablesinstead of the primitive variables to improve the monotonicity of the numerical solution[41]. The positivity of the viscous and resistive heating rates discussed in Sec. 2.6 maybe violated when the Roe-type method is applied. We actually found that the Roe-typemethod provides less robustness than the LLF-type method when the internal energy ismuch smaller than the kinetic or magnetic energy.
The time step size ∆t is determined from the Courant-Friedrichs-Lewy (CFL) condition ∆t = min x ,x ,x σ CFL " N D X i =1 max( | λ i, ( x , x , x ) | , | λ i, ( x , x , x ) | ) h i − , (68)where λ i, and λ i, are the eigenvalues for the fast magneto-acoustic waves calculated for the x i -direction, σ CFL is the Courant number of order unity, and N D = 1 , , σ CFL = 0 .
3. Although most problems can be solved with larger σ CFL like 0 . .
5, thevalue of 0 . c ψ is computed from theCFL condition at each time step as c ψ = σ CFL " ∆t N D X i =1 h i − . (69)The damping time τ ψ was set to 10 ∆t . More detailed discussion on the choice of the dampingtime τ ψ can be found in [24, 42].
4. Test problems
Based on the newly proposed formulation described in Sec. 2.8, four schemes wereimplemented, namely, JST-LLF, JST-Roe, WENO5-LLF, and WENO5-Roe schemes. Inthe JST-LLF and JST-Roe schemes, the second-order central difference with the filteringflux based on the SLIP version of the Jameson-Schmidt-Turkel scheme was used. In the twoWENO5 schemes, the sixth-order central difference with the filtering flux based on the fifth-order WENO reconstruction is used. The JST-LLF and WENO5-LLF schemes are based onthe LLF-type filtering scheme, and the JST-Roe and WENO5-Roe schemes are implementedwith the Roe-type filtering scheme. See Section 3 for the details of our implementation. Wefound that the LLF-type filtering scheme becomes more numerically robust than the Roe-type scheme. Most of the test problems presented in this study can be solved with both theLLF- and Roe-type filtering schemes, unless otherwise noted. The exceptions are the moststringent, non-standard test problems described in Secs. 4.5 and 4.6.All four schemes were integrated using the three-step third-order optimal SSP Runge–Kutta method [43]. The Courant number σ CFL of 0 . We analyze the spatial accuracy of the four proposed schemes (JST-LLF, JST-Roe,WENO5-LLF, and WENO5-Roe) using the magnetized iso-density vortex problem proposedby Balsara [44]. As argued in [45], the Gaussian taper of variables in the original problemproduces a small error near the boundary, which may reduce the convergence rate of higher-order schemes. To avoid this effect, we chose a version of this problem described in Mignoneet al. [42]. The initial condition is described by ρ = 1, ( V , V ) = (1 ,
1) + ( − x , x ) κ e q (1 − r ) ,( B , B ) = ( − x , x ) µ e q (1 − r ) , and p = 1 + [ µ (1 − qr ) − κ ρ ]e q (1 − r ) / (4 q ), where r = p x + x , κ = 1 / (2 π ), and µ = 1 / (2 π ). In the actual implementation, we use the mag-netic vector potential A = µ e q (1 − r ) / (2 q ) to compute the magnetic field as ( B , B ) =( δA /δx , − δA /δx ) so that the discretization error of ∇· B becomes zero in the initial17 able 1: L norm errors and corresponding convergence rates for the two-dimensional magnetized iso-density vortex problem. The errors are measured in the x -component of the magnetic field B and thedivergence of the magnetic field ∇· B . Note that we set q = 0 . q = 1 . Method Number of mesh L error of B L order of B L error of ∇· B L order of ∇· B JST-LLF 32 ×
32 1 . × − . × − ×
64 5 . × − .
05 4 . × − − . ×
128 1 . × − .
72 1 . × − . ×
256 4 . × − .
81 4 . × − . ×
32 7 . × − . × − ×
64 2 . × − .
39 4 . × − . ×
128 9 . × − .
69 1 . × − . ×
256 2 . × − .
74 3 . × − . ×
32 7 . × − . × − ×
64 2 . × − .
97 1 . × − ×
128 5 . × − .
98 3 . × − ×
256 1 . × − .
00 7 . × − WENO5-LLF 32 ×
32 4 . × − . × − ×
64 5 . × − .
06 2 . × − . ×
128 3 . × − .
08 8 . × − . ×
256 1 . × − .
37 3 . × − . ×
32 2 . × − . × − ×
64 2 . × − .
45 1 . × − . ×
128 1 . × − .
13 6 . × − . ×
256 6 . × − .
49 2 . × − . ×
32 6 . × − . × − ×
64 1 . × − .
90 2 . × − ×
128 1 . × − .
90 5 . × − ×
256 4 . × − .
18 2 . × − condition. The simulation evolves over a time of 10 units in the two-dimensional domainof [ − , × [ − ,
5] with periodic boundary conditions. Following [42], we set q = 0 . q = 1 . σ CFL = 0 . L errors and corresponding convergence rates measured in B and ∇· B . The L error of a quantity U is computed as an average of the absolute valueof the difference between the initial and final solutions over the whole domain. All schemesconverge to the designed order of accuracy in both B and ∇· B . The errors and convergencerates of the magnetic field B can be directly compared to Table 3 of [42]. As expected, oursecond-order schemes provide larger errors than the third-order schemes described in [42].The errors and convergence rates by the second-order total variation diminishing scheme inTable 5 of [44] seems to be comparable or slightly worse than our JST-LLF scheme. Ourfifth-order schemes also provides larger L errors than the fifth-order schemes in [42]. Thisdifference is caused by the larger numerical diffusion by the enhancement of the filteringflux (by setting c diff = 2) and the less diffusive reconstruction by WENO-Z [46] or MP5 [33]18chemes used in [42].For comparison, we also show the results of two additional central difference schemesCD2 and CD6 in Table 1. The CD2 and CD6 schemes are the second- and sixth-ordercentral difference schemes derived by setting the filtering flux to be zero in the JST andWENO5 schemes. The convergence rates of these central schemes are independent from thefiltering scheme and the numerical divergence of the magnetic field. Both schemes are lessdiffusive and rapidly converge to the designed order of accuracy. As the amount of ∇· B produced by the central difference schemes (CD2 and CD6) is limited by the machine epsilonof the floating-point operation, the order of convergence for ∇· B is not shown. The slightreduction in the convergence rate between 128 ×
128 and 256 ×
256 grid points in the CD6scheme is caused by the third-order error of the time integration scheme. If we compute thecase of 256 ×
256 grid points of the CD6 scheme with smaller CFL number of 0 . L error and convergence rate of B is 2 . × − and 6 .
00, respectively.
Figure 3: One-dimensional Sod’s shock tube problem using 100 grid points. Shown are the mass density(left panel) and the viscous heating rate (right panel) at t = 0 . The hydrodynamic Sod’s shock tube problem [47] is used to measure the performanceof the proposed schemes near the shock wave and discontinuity. Figure 3 shows the one-dimensional simulation of the Sod’s problem calculated with 100 grid points. The proposedschemes (JST-Roe and WENO5-Roe) successfully captured the shock front, contact discon-tinuity, and rarefaction wave. One advantage of the proposed scheme is that the numericalviscous heating can be derived explicitly. The result indicates that the viscous heating wasconcentrated near the shock front near x ∼ .
9, with a positive sign indicating the heatingprocess. The viscous heating was negligible at the contact discontinuity or rarefaction wave.19 igure 4: One-dimensional Shu–Osher’s shock tube problem using 300 grid points. Shown are the massdensity (left panel) and the viscous heating rate (right panel) at t = 0 .
18 for the JST-Roe scheme (bluediamond) and the WENO5-Roe scheme (red cross). The solid lines indicate the reference solution calculatedwith the WENO5-Roe scheme using 4000 grid points.
The advantage of the higher-order scheme can be apparent in the Shu–Osher’s shocktube problem [43]. Figure 4 shows the results of the one-dimensional Shu-Osher’s test using300 grid points. The higher-order WENO5-Roe scheme could resolve the wavy pattern afterthe shock front better than the JST-Roe scheme. The viscous heating is concentrated nearthe shock front at x ∼ . In the one-dimensional MHD test problems, the discretization error of the magnetic fielddivergence becomes zero. Thus, we can assess the performance of the proposed schemeswithout the influence of the ∇· B error.The Dai–Woodward’s MHD shock tube problem [48] involves various shock waves anddiscontinuities (i.e., the fast shocks, rotational discontinuities, and slow shocks propagat-ing from each side of the contact discontinuity). Figure 5 shows the results of the one-dimensional Dai-Woodward shock tube problem with 200 grid points. The JST-Roe andWENO5-Roe schemes were sufficiently robust and captured all shocks and discontinuitiesin this problem. No visible numerical overshoots were observed. The viscous and resistiveheating rates are prominently located near the two fast shock fronts.The second MHD shock tube problem considered in this study is the MHD analog of theSod’s shock tube problem originally introduced by Brio & Wu [49], with the specific heatratio γ = 2. Our test problem is a variant of this problem with γ = 5 / x ∼ .
5. The numericalsolutions calculated with both the JST-Roe and WENO5-Roe schemes agree well with thepreviously reported results. Most of the viscous and resistive heating rates appeared near20 igure 5: One-dimensional Dai–Woodward’s shock tube problem using 200 grid points. Shown are themass density (top left panel), the y-component of magnetic field (top right panel), the viscous heating rate(bottom left panel), and the resistive heating rate (bottom right panel) at t = 0 . the compound shock and slow shock at x ∼ .
65. Note that the ratio between the viscousand resistive heating rates depends on the dissipation methods (e.g., LLF or Roe).
Multi-dimensional MHD problems suffer from the numerical error in the divergence ofthe magnetic field. The robustness and energy-conservation of the proposed schemes in thetwo-dimensional domain are tested using the Orszag–Tang vortex problem [50]. There areminor variations of the numerical settings in the previous literatures. The numerical setupused in this paper is identical to that of Ryu et al. [51, 52]. The numerical solution isshown in Fig. 7. The spatial profiles of gas and magnetic pressures agreed well with theprevious reports (e.g., Fig. 3 of Ryu et al. [52]). The viscous and resistive heating rateswere concentrated near the shock fronts.When the total energy equation is not directly solved in the MHD schemes, the totalenergy conservation may be numerically violated. The proposed scheme is designed to satisfythe total energy conservation in terms of the spatial discretization. However, the totalenergy conservation may be violated by the numerical error of the temporal discretization21 igure 6: One-dimensional Brio–Wu’s shock tube problem using 200 grid points. Shown are the massdensity (top left panel), the y-component of magnetic field (top right panel), the viscous heating rate(bottom left panel), and the resistive heating rate (bottom right panel) at t = 0 . as discussed in Sec. 3.1. In our formulation, the conservation error of the total energydepends only on the discretization error of the time integration scheme. Because the third-order Runge-Kutta method was used for time integration, the conservation error of the totalenergy should converge by third order against the size of time step used in the simulations.The convergence of the total energy conservation in terms of the Courant number of thetime step criterion is shown in Fig. 8. Here, we employed the lower spatial resolutionwith 128 ×
128 grid points to reduce the computational cost. The observed third-orderconvergence of the conservation error of the total energy is consistent with our expectation.We also checked the conservation error of the momentum. Due to the spatial symmetryof this test problem, the conservation error of the momentum was kept constant in theround-off error. Instead, we evaluated the physical part of the Lorentz force Eq. 9 andthe artificial field-aligned force B i ∇· B . The standard deviation (in time and space) of the x -component of the field-aligned force was several tens of percent of the physical part of theLorentz force. The amount of this field-aligned force is a direct consequence of the numerical ∇· B error, which depends on the filtering scheme, the spatial resolution, and the parametersof the divergence cleaning method ( c ψ and τ ψ in this study). As we discussed in Sec. 2.7,22 igure 7: Two-dimensional Orszag–Tang vortex problem calculated with the WENO5-Roe scheme using256 ×
256 grid points. Shown are (a) the gas pressure, (b) the magnetic pressure, (c) the viscous heatingrate, and (d) the resistive heating rate at t = 0 . the present scheme can be modified so that the momentum conservation is satisfied evenwhen ∇· B is nonzero. This alternative version may be used in the problems where the strictmomentum conservation by the Lorentz force is required.One of the advantages of the proposed scheme is that the analysis of the energy variationsin the numerical solution can be accurate and straightforward. Figure 9 shows the temporalvariation of the internal, kinetic, and magnetic energies and the contributions to their vari-ation. The initial kinetic energy is converted into the internal and magnetic energies duringthe time evolution. The prominent energy exchange is through the works by the pressuregradient and Lorentz forces. Both the pressure gradient force the viscous heating convertthe kinetic energy into the internal energy. The temporal change of the magnetic energyis dominated by the work done by the Lorentz force. The resistive heating has only minoreffect in the energy conversion. The heating from the ∇· B error, i.e., Q ψ , is negligible (lessthan 4 × − as the volume-averaged value). This small contribution of Q ψ can be explainedthat this value is roughly proportional to the square of the divergence of the magnetic field.For each energy equation, all the contributions are fully balanced, indicating the accuracy23 igure 8: Conservation error of the total energy for Orszag–Tang vortex problem calculated with theWENO5-Roe scheme using 128 ×
128 grid points. Shown is difference between the initial ( t = 0) and finalstep ( t = 0 .
48) of the volume-averaged total energy density. of the energy analysis presented here.
The two-dimensional MHD blast wave problem proposed by Balsara & Spicer [53] is oneof the most stringent among the standard test problems. In the original version, the plasmabeta (i.e., the ratio between the gas pressure and magnetic pressure) in the ambient is setto approximately 2 . × − . This problem is a severe numerical benchmark for the MHDsolvers because the positivity of the gas pressure can be easily violated due to the shockpropagation in the extremely strong magnetic field. We found that the JST-Roe, JST-LLF,and WENO5-LLF schemes can handle this problem successfully, whereas the WENO5-Roescheme failed due to the negative pressure during the time integration.To assess the further robustness of the proposed schemes, we calculated the MHD blastproblem in the non-standard setting with a lower plasma beta. The initially uniform massdensity, velocity, and magnetic field are imposed as ( ρ, V , V , B , B ) = (1 , , , / √ π, − except in a central circle of radius 0 .
1, where thepressure is set to 1000. The only difference from the original problem [53] is the much lowergas pressure in the ambient region. The resultant plasma beta in the ambient region is2 . × − , which is the most stringent setup among the family of the MHD blast waveproblems tested in previous literatures [53, 35, 54, 55]. The numerical solution obtainedby 256 ×
256 grid points with JST-Roe scheme is shown in Fig. 10. The JST-Roe and JST-LLF schemes can handle this problem, but WENO5-LLF scheme causes a negative pressureduring the time integration. We found that the two JST schemes are extremely robust inthis problem and stable even when we set the ambient gas pressure to 10 − (i.e., the ambient24 igure 9: Contributions in energy equations for Orszag–Tang vortex problem calculated with the WENO5-Roe scheme using 256 ×
256 grid points. Shown are (a) the temporal variation of volume-averaged energydensities, (b) the gain and loss of internal energy, (c) the gain and loss of kinetic energy, and (d) the gainand loss of magnetic energy. plasma beta of 2 . × − ). The two-dimensional MHD rotor problem originally proposed by Balsara & Spicer [53]involves a rotating disk with an initial magnetic field perpendicular to the rotation axis.The strong rotational discontinuities along with the shocks and rare factions are generatedby the shearing and expansion motion of the rotating disk. Among the several variations ofthis problem, the second rotor test described by T´oth [23] was selected for this study. All ofthe proposed schemes (JST-Roe, JST-LLF, WENO5-Roe, and WENO5-LLF) are capableof solving this problem.A non-standard alternative of the rotor problem with more stringent parameter was alsocarried out. The numerical setup is identical to that of the T´oth’s second rotor problem [23]except the initial gas pressure and linear taper. The initial gas pressure is uniformly set to5 × − whereas it is 0 . . × − .The linear taper is applied between a radius of 0 . .
13 so that the rotor’svelocity and mass density linearly approach to the background over six grid points when the25 igure 10: Two-dimensional MHD blast problem calculated with the JST-Roe scheme using 256 ×
256 gridpoints. Shown are (a) the mass density, (b) the gas pressure, (c) the Mach number, and (d) the plasmabeta at t = 0 .
01. The range of the color map is determined by the maximum and minimum values at thesnapshot (except the lower range of 0 . problem is solved with 200 ×
200 grid points. Figure 11 shows the numerical solution solvedby the WENO5-LLF scheme. The expansion of the rotating disk caused the region with verylow pressure near the center of the numerical domain. At the final snapshot ( t = 0 . × − and the local Mach number exceeds 150 in thelow-pressure region. These values of the Mach number and plasma beta are more severethan those used in previous studies [53, 35, 54, 55], making the problem a severe benchmarkfor the numerical robustness.We found that two JST schemes and the WENO5-Roe scheme failed to handle thisproblem with the initial gas pressure of 5 × − . By changing the initial gas pressure from0 . × − multiplying a factor of 0 . − , 10 − , and 10 − , respectively. The poor robustness of the JST-Roe and WENO5-Roeschemes can be explained that the viscous and resistive heating can be negative for the Roe-type model. The higher robustness of the WENO5-LLF scheme over the JST-LLF schememay have been caused by the higher gradient of the solution (see panel a in Fig. 11). We26 igure 11: Two-dimensional MHD rotor problem calculated with the WENO5-LLF scheme using 200 × t = 0 . . found that the JST-LLF and JST-Roe schemes can solve the problem with the initial gaspressure of 5 × − and 5 × − , respectively, when the finer numerical mesh of 400 × ×
400 grid points.
5. Conclusions
We have presented an energy-consistent formulation of the compressible magnetohydro-dynamic equations. The transport and interaction of the internal, kinetic, and magneticenergies are satisfied in the discrete sense, while maintaining the standard conservationproperty of the mass, momentum, magnetic flux, and total energy. These characteristicshave been accomplished through the application of a simple constructive strategy of usingthe discrete versions of the product rule. The shock capturing was achieved by introducingthe nonlinear filtering flux in the discretized equations. The energy-consistent formulation of27he nonlinear filtering flux for both hydrodynamic and MHD equations was also suggested.The viscous and resistive heating rates become positive when the filtering flux satisfies thespecific conditions.The proposed formulations were implemented with the spatially second-order and thesixth-order central difference operators. The filtering flux is developed based the second-order JST scheme and fifth-order WENO reconstruction with the LLF-type and Roe-typefiltering schemes. All of these schemes are integrated by the third-order SSP Runge-Kuttamethod. The combination of the energy-consistent finite difference formulation and theappropriate filtering with small numerical overshoot was shown to yield the excellent ro-bustness for the most stringent problems.In this study, we assume the periodic boundary conditions in all test problems presentedin Sec. 4. The conservation property of the (magneto)hydrodynamic variables is affectedby the non-periodic boundary conditions in practical problems. For the practical imple-mentation of the non-periodic boundary conditions in the finite difference schemes with thesecondary conservation property, we refer the reader to the discussion in Morinishi [2] andDesjardins [14].
Acknowledgements
This work was supported by JSPS KAKENHI grant No. JP19K14756. A major partof numerical computations was carried out on Cray XC50 at Center for ComputationalAstrophysics, National Astronomical Observatory of Japan. This work was supported bythe computational joint research program of the Institute for Space-Earth EnvironmentalResearch (ISEE), Nagoya University. Some numerical computations were carried out on theFX100/Flow supercomputer system at the Information Technology Center, Nagoya Univer-sity. Analysis of the numerical simulations were carried out in the Center for Integrated DataScience, Institute for Space-Earth Environmental Research, Nagoya University through thejoint research program.
References [1] F. H. Harlow, J. E. Welch, Numerical Calculation of Time-Dependent Viscous Incompressible Flow ofFluid with Free Surface, Physics of Fluids 8 (12) (1965) 2182–2189. doi:10.1063/1.1761178 .[2] Y. Morinishi, T. S. Lund, O. V. Vasilyev, P. Moin, Fully Conservative Higher Order Finite Dif-ference Schemes for Incompressible Flow, Journal of Computational Physics 143 (1) (1998) 90–124. doi:10.1006/jcph.1998.5962 .[3] Y. Morinishi, O. V. Vasilyev, T. Ogi, Fully conservative finite difference scheme in cylindrical coordi-nates for incompressible flow simulations, Journal of Computational Physics 197 (2) (2004) 686–710. doi:10.1016/j.jcp.2003.12.015 .[4] A. G. Kravchenko, P. Moin, On the Effect of Numerical Errors in Large Eddy Simulations of TurbulentFlows, Journal of Computational Physics 131 (2) (1997) 310–322. doi:10.1006/jcph.1996.5597 .[5] C. A. Kennedy, A. Gruber, Reduced aliasing formulations of the convective terms within the Navier-Stokes equations for a compressible fluid, Journal of Computational Physics 227 (3) (2008) 1676–1700. doi:10.1016/j.jcp.2007.09.020 .[6] S. Pirozzoli, Numerical Methods for High-Speed Flows, Annual Review of Fluid Mechanics 43 (1) (2011)163–194. doi:10.1146/annurev-fluid-122109-160718 .
7] A. Jameson, W. Schmidt, E. Turkel, Numerical solution of the Euler equations by finite volume methodsusing Runge Kutta time stepping schemes (Jun. 1981).[8] F. Ducros, F. Laporte, T. Soul`eres, V. Guinot, P. Moinat, B. Caruelle, High-Order Fluxes for Conserva-tive Skew-Symmetric-like Schemes in Structured Meshes: Application to Compressible Flows, Journalof Computational Physics 161 (1) (2000) 114–139. doi:10.1006/jcph.2000.6492 .[9] A. Jameson, The construction of discretely conservative finite volume schemes that also globally con-serve energy or entropy, Journal of Scientific Computing 34 (2) (2008) 152–187.[10] A. Jameson, Formulation of kinetic energy preserving conservative schemes for gas dynamics and directnumerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy andkinetic energy preserving schemes, Journal of Scientific Computing 34 (2) (2008) 188–208.[11] S. Pirozzoli, Generalized conservative approximations of split convective derivative operators, Journalof Computational Physics 229 (19) (2010) 7180–7190. doi:10.1016/j.jcp.2010.06.006 .[12] J. C. Kok, A high-order low-dispersion symmetry-preserving finite-volume method for compress-ible flow on curvilinear grids, Journal of Computational Physics 228 (18) (2009) 6811–6832. doi:10.1016/j.jcp.2009.06.015 .[13] Y. Morinishi, Skew-symmetric form of convective terms and fully conservative finite difference schemesfor variable density low-Mach number flows, Journal of Computational Physics 229 (2) (2010) 276–300. doi:10.1016/j.jcp.2009.09.021 .[14] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservative finite difference scheme forvariable density low Mach number turbulent flows, Journal of Computational Physics 227 (15) (2008)7125–7159. doi:10.1016/j.jcp.2008.03.027 .[15] Y. Kuya, K. Totani, S. Kawai, Kinetic energy and entropy preserving schemes for compress-ible flows by split convective forms, Journal of Computational Physics 375 (2018) 823–853. doi:10.1016/j.jcp.2018.08.058 .[16] H. C. Yee, N. D. Sandham, M. J. Djomehri, Low-Dissipative High-Order Shock-Capturing Meth-ods Using Characteristic-Based Filters, Journal of Computational Physics 150 (1) (1999) 199–238. doi:10.1006/jcph.1998.6177 .[17] T. Shiroto, S. Kawai, N. Ohnishi, Structure-preserving operators for thermal-nonequilibriumhydrodynamics, Journal of Computational Physics 364 (2018) 1–17. arXiv:1705.03136 , doi:10.1016/j.jcp.2018.03.008 .[18] M.-J. Ni, R. Munipalli, N. B. Morley, P. Huang, M. A. Abdou, A current density conservative scheme forincompressible MHD flows at a low magnetic Reynolds number. Part I: On a rectangular collocated gridsystem, Journal of Computational Physics 227 (1) (2007) 174–204. doi:10.1016/j.jcp.2007.07.025 .[19] M.-J. Ni, R. Munipalli, P. Huang, N. B. Morley, M. A. Abdou, A current density conservative schemefor incompressible MHD flows at a low magnetic Reynolds number. Part II: On an arbitrary collocatedmesh, Journal of Computational Physics 227 (1) (2007) 205–228. doi:10.1016/j.jcp.2007.07.023 .[20] M.-J. Ni, J.-F. Li, A consistent and conservative scheme for incompressible MHD flows at a low magneticReynolds number. Part III: On a staggered mesh, Journal of Computational Physics 231 (2) (2012)281–298. doi:10.1016/j.jcp.2011.08.013 .[21] Y. Idomura, M. Ida, S. Tokuda, L. Villard, New conservative gyrokinetic full- f Vlasov code and itscomparison to gyrokinetic δ f particle-in-cell code, Journal of Computational Physics 226 (1) (2007)244–262. doi:10.1016/j.jcp.2007.04.013 .[22] T. Shiroto, N. Ohnishi, Y. Sentoku, Quadratic conservative scheme for relativistic Vlasov-Maxwell system, Journal of Computational Physics 379 (2019) 32–50. arXiv:1802.07238 , doi:10.1016/j.jcp.2018.10.041 .[23] G. T´oth, The ∇· B = 0 Constraint in Shock-Capturing Magnetohydrodynamics Codes, Journal ofComputational Physics 161 (2000) 605–652. doi:10.1006/jcph.2000.6519 .[24] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz, T. Schnitzer, M. Wesenberg, Hyperbolic Diver-gence Cleaning for the MHD Equations, Journal of Computational Physics 175 (2002) 645–673. doi:10.1006/jcph.2001.6961 .[25] A. Jameson, Origins and further development of the jameson–schmidt–turkel scheme, AIAA Journal doi:10.1080/10618569508904524 .[28] B. van Leer, Towards the Ultimate Conservation Difference Scheme. II. Monotonicity and Conserva-tion Combined in a Second-Order Scheme, Journal of Computational Physics 14 (4) (1974) 361–370. doi:10.1016/0021-9991(74)90019-9 .[29] G.-S. Jiang, C.-W. Shu, Efficient Implementation of Weighted ENO Schemes, Journal of ComputationalPhysics 126 (1996) 202–228. doi:10.1006/jcph.1996.0130 .[30] C.-W. Shu, S. Osher, Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes,Journal of Computational Physics 77 (2) (1988) 439–471. doi:10.1016/0021-9991(88)90177-5 .[31] S. Gottlieb, D. I. Ketcheson, C.-W. Shu, High order strong stability preserving time discretizations,Journal of Scientific Computing 38 (3) (2009) 251–289. doi:10.1007/s10915-008-9239-z .[32] A. Jameson, T. Baker, Solution of the euler equations for complex configurations, in: 6th ComputationalFluid Dynamics Conference Danvers, 1983, p. 1929.[33] A. Suresh, H. T. Huynh, Accurate Monotonicity-Preserving Schemes with Runge Kutta Time Stepping,Journal of Computational Physics 136 (1997) 83–99. doi:10.1006/jcph.1997.5745 .[34] D. S. Balsara, Linearized Formulation of the Riemann Problem for Adiabatic and IsothermalMagnetohydrodynamics, The Astrophysical Journal Suppliment Series 116 (1) (1998) 119–131. doi:10.1086/313092 .[35] J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley, J. B. Simon, Athena: A New Code for Astrophys-ical MHD, The Astrophysical Journal Suppliment Series 178 (1) (2008) 137–177. arXiv:0804.0402 , doi:10.1086/588755 .[36] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Com-munications on pure and applied mathematics 7 (1) (1954) 159–193.[37] A. Kurganov, E. Tadmor, New High-Resolution Central Schemes for Nonlinear Conservation Lawsand Convection-Diffusion Equations, Journal of Computational Physics 160 (1) (2000) 241–282. doi:10.1006/jcph.2000.6459 .[38] P. L. Roe, Approximate riemann solvers, parameter vectors, and difference schemes, Journal of com-putational physics 43 (2) (1981) 357–372.[39] D. Ryu, T. W. Jones, Numerical magetohydrodynamics in astrophysics: Algorithm and tests forone-dimensional flow‘, The Astrophysical Journal 442 (1995) 228–258. arXiv:astro-ph/9404074 , doi:10.1086/175437 .[40] A. Harten, J. M. Hyman, Self-Adjusting Grid Methods for One-Dimensional Hyperbolic ConservationLaws, Journal of Computational Physics 50 (2) (1983) 235–269. doi:10.1016/0021-9991(83)90066-9 .[41] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly High Order Accurate Es-sentially Non-oscillatory Schemes, III, Journal of Computational Physics 131 (1) (1997) 3–47. doi:10.1006/jcph.1996.5632 .[42] A. Mignone, P. Tzeferacos, G. Bodo, High-order conservative finite difference GLM-MHD schemes forcell-centered MHD, Journal of Computational Physics 229 (17) (2010) 5896–5920. arXiv:1001.2832 , doi:10.1016/j.jcp.2010.04.013 .[43] C.-W. Shu, S. Osher, Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes,II, Journal of Computational Physics 83 (1) (1989) 32–78. doi:10.1016/0021-9991(89)90222-2 .[44] D. S. Balsara, Second-Order-accurate Schemes for Magnetohydrodynamics with Divergence-free Reconstruction, The Astrophysical Journal Supplement Series 151 (2004) 149–184. arXiv:astro-ph/0308249 , doi:10.1086/381377 .[45] M. Dumbser, D. S. Balsara, E. F. Toro, C.-D. Munz, A unified framework for the construction ofone-step finite volume and discontinuous Galerkin schemes on unstructured meshes, Journal of Com-putational Physics 227 (2008) 8209–8253. doi:10.1016/j.jcp.2008.05.025 .
46] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatoryscheme for hyperbolic conservation laws, Journal of Computational Physics 227 (2008) 3191–3211. doi:10.1016/j.jcp.2007.11.038 .[47] G. A. Sod, Review. A Survey of Several Finite Difference Methods for Systems of Non-linear Hyperbolic Conservation Laws, Journal of Computational Physics 27 (1) (1978) 1–31. doi:10.1016/0021-9991(78)90023-2 .[48] W. Dai, P. R. Woodward, On the Divergence-free Condition and Conservation Laws in NumericalSimulations for Supersonic Magnetohydrodynamical Flows, The Astrophysical Journal 494 (1998) 317–335. doi:10.1086/305176 .[49] M. Brio, C. C. Wu, An upwind differencing scheme for the equations of ideal magnetohydrodynamics,Journal of Computational Physics 75 (1988) 400–422. doi:10.1016/0021-9991(88)90120-9 .[50] S. A. Orszag, C. M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence,Journal of Fluid Mechanics 90 (1979) 129–143. doi:10.1017/S002211207900210X .[51] D. Ryu, T. W. Jones, A. Frank, Numerical Magnetohydrodynamics in Astrophysics: Algorithm andTests for Multidimensional Flow, The Astrophysical Journal 452 (1995) 785. arXiv:astro-ph/9505073 , doi:10.1086/176347 .[52] D. Ryu, F. Miniati, T. W. Jones, A. Frank, A Divergence-free Upwind Code for Multidimensional Mag-netohydrodynamic Flows, The Astrophysical Journal 509 (1998) 244–255. arXiv:astro-ph/9807228 , doi:10.1086/306481 .[53] D. S. Balsara, D. S. Spicer, A Staggered Mesh Algorithm Using High Order Godunov Fluxes to EnsureSolenoidal Magnetic Fields in Magnetohydrodynamic Simulations, Journal of Computational Physics149 (1999) 270–292. doi:10.1006/jcph.1998.6153 .[54] D. S. Balsara, Divergence-free reconstruction of magnetic fields and WENO schemes for magne-tohydrodynamics, Journal of Computational Physics 228 (2009) 5040–5056. arXiv:0811.2192 , doi:10.1016/j.jcp.2009.03.038 .[55] D. Lee, A solution accurate, efficient and stable unsplit staggered mesh scheme for three dimensionalmagnetohydrodynamics, Journal of Computational Physics 243 (2013) 269–292. arXiv:1303.6988 , doi:10.1016/j.jcp.2013.02.049 ..