Invariant finite-difference schemes with conservation laws preservation for one-dimensional MHD equations

Invariant finite-difference schemes are considered for one-dimensional magnetohydrodynamics (MHD) equations in mass Lagrangian coordinates for the cases of finite and infinite conductivity. For construction these schemes previously obtained results of the group classification of MHD equations are used. On the basis of the classical Samarskiy-Popov scheme new schemes are constructed for the case of finite conductivity. These schemes admit all symmetries of the original differential model and have difference analogues of all of its local differential conservation laws. Among the conservation laws there are previously unknown ones. In the case of infinite conductivity, conservative invariant schemes constructed as well. For isentropic flows of a polytropic gas proposed schemes possess the conservation law of energy and preserve entropy on two time layers. This is achieved by means of specially selected approximations for the equation of state of a polytropic gas. Also, invariant difference schemes with additional conservation laws are proposed. A new scheme for the case of finite conductivity is tested numerically for various boundary conditions which shows accurate preservation of difference conservation laws.


Introduction
Magnetic hydrodynamics equations describe the flows of electrically conductive fluids such as plasma, liquid metals, and electrolytes and are widely used in modeling processes in various fields from engineering to geophysics and astrophysics.
In the present publication, we restrict ourselves to considering plane one-dimensional MHD flows under the assumption that the medium is inviscid and thermally non-conducting. A group classification of the MHD equations under the above conditions was carried out recently in [1] (for some particular results see also [2][3][4][5]). The group classification splits into four essentially different cases according to whether the conductivity of the medium is finite or infinite, and the longitudinal component of the magnetic field vector is zero or a non-zero constant.
The MHD equations are nonlinear, so that even in the one-dimensional case only their particular solutions are known [6][7][8][9][10]. Therefore, numerical modeling in magnetohydrodynamics is of great practical interest. There are many approaches to numerically modeling MHD equations, including finite-difference, finite element and finite volume methods (see, e.g. [11][12][13][14][15][16][17][18][19]). Further we consider finite-difference schemes taking as a starting point the classical Samarsky-Popov schemes [12,13] for the MHD equations for the case of finite conductivity. The main properties of the considered schemes are invariance, i.e. preservation of the symmetries of the original differential equations, and the presence of difference analogues of local differential conservation laws. It is known that there is a connection between the invariance of equations and the presence of conservation laws [20][21][22][23].

arXiv:2112.03118v2 [math.NA] 26 Feb 2022
Invariant schemes have been studied for a long time [24][25][26][27][28], and over the past decades, significant progress has been made in the development of methods for their construction and integration. For schemes for ordinary differential equations with Lagranagian or Hamiltonian functions, a number of methods [29][30][31][32] have been developed that make it possible to decrease the order or even integrate the schemes. A method based on the Lagrangian identity has also been developed for the case when the equations do not admit a variational formulation [33,34].
For partial difference schemes, the main methods used are the method of differential invariants [28,35,36] and the difference analogue of the direct method [28,37]. Using these methods, the authors have constructed invariant schemes for various shallow water models [38][39][40][41]. Also, some previously known schemes have been investigated from a group analysis point of view. In particular, symmetries and conservation laws of the Samarskiy-Popov schemes for the one-dimensional gas dynamics equations of a polytropic gas have been investigated in [42][43][44]. Based on the results of the group classification [1] and Samarskiy-Popov schemes for the MHD equations, we further construct invariant finite-difference schemes possessing conservation laws. The set and number of conservation laws depend on the conductivity, the form of the magnetic field vector and the equation of state of the medium. This paper is organized as follows. In Section 2 the simplest version of the finiteconductivity MHD equations in mass Lagrangian coordinates in case of one-dimensional plane flows is considered. Electric and magnetic fields are represented by one-component vectors, which greatly simplifies the form of the equations. This was the main case considered in Samarsky and Popov's publications [12,13]. The section also provides basic notation and definitions. Then, symmetries and conservation laws of the Samarsky-Popov scheme for the MHD equations are investigated. In addition to the previously known conservation laws, the center-of-mass conservation law is given, as well as new conservation laws for the specific conductivity function, obtained on the basis of the group classification [1]. Section 3 is devoted to various generalizations of the scheme of Section 2. In Section 3.1 the scheme for arbitrary electric and magnetic fields is considered. Its symmetries are investigated and conservation laws are given. The case of infinite conductivity is considered in Section 3.2. It is shown that in this case the Samarsky-Popov scheme requires some additional modifications in order to possess the conservation law of angular momentum. In the case of a polytropic gas, it turns out to be possible to preserve not only energy, but also entropy along pathlines. This can be done using a specially selected equation of state for a polytropic gas. At the end of the section, an example of an invariant scheme is given that does not possess a conservation law of energy, but preserves entropy and has additional conservation laws in the case of isentropic flows. In Section 4, one of the invariant schemes for the case of finite conductivity is numerically implemented for the example of plasma bunch deceleration by crossed electromagnetic fields. The results are discussed in the Conclusion.

Conservative schemes for MHD equations with finite conductivity
Problems of continuum mechanics and plasma physics are often considered in mass Lagrangian coordinates [13,45] since for them the formulation of boundary conditions is greatly simplified. In particular, the conservative Samarsky-Popov schemes for the equations of gas dynamics and magnetohydrodynamics have been constructed in mass Lagrangian coordinates.
In mass Lagrangian coordinates the MHD equations, describing the plane one-dimensional MHD flows, are [1,13] 1 where t is time, s is mass Lagrangian coordinate, x is Eulerian coordinate, ρ is density, p is pressure, ε is internal energy, u = (u, v, w) is the velocity of a particle, E = (E x , E y , E z ) is the electric field vector, H = (H x , H y , H z ) is the magnetic field vector, and i = (i x , i y , i z ) is the electric current. The conductivity σ is some function of p and ρ, i.e., σ = σ(p, ρ).
Following [13] we firstly consider the simplest case of one-component electric and magnetic fields. Here we also introduce the notation and some basic concepts. In the next section some generalizations are considered, including the case of infinite conductivity. For simplicity, the longitudinal component of the magnetic field H is set to zero, and the coordinate system is chosen in such a way that H = (0, H, 0). Consequently, the electric current i and the electric field E are also one-component vectors, i.e., i = (0, 0, i), E = (0, 0, E). Electromagnetic force f = ( f , 0, 0) acts in the x-direction, and the velocity is u = (u, 0, 0) (see Figure 1).
Given the above, the system of the one-dimensional MHD equations with a finite conductivity |σ| < ∞ in mass Lagrangian coordinates can be written as [13] 1 x t = u, where κ = 1/(4π) and q is Joule heating per unit mass.
In particular, we consider a polytropic gas for that the following relation holds Equation (2e) for the energy evolution can be rewritten in the semi-divergent form or in the divergent form Notice that the electromagnetic force f = −iH/ρ can be represented in the divergent form f = −κ(H 2 /2) s , and equation (2b) can be rewritten as Further we assume κ = 1 since it can be discarded by means of the scaling transformatioñ s = κs,p = κ p,ρ = κρ,σ = κσ.
Here and further φ t , φˇt and φ s , φs denote finite-difference derivatives of some quantity φ = φ(t n , s m , u n m , ...) which are defined with the help of the finite-difference right and left shifts along the time and space axes correspondingly  The indices n and m are respectively changed along time and space axes t and s. The time and space steps are defined as follows Following the Samarskiy-Popov notation throughout the text we denote and Notice that on a uniform lattice h i = h = const in its integral nodes (13) becomes Remark 1. The energy equation (8e) can be reduced to one of the three following forms [13] using equivalent algebraic transformations: These different forms of equation reflect the balance of certain types of energy, i.e. they express the different physical aspects of energy conservation. To emphasize this property, such schemes are also called completely conservative.

Invariance of Samarskiy-Popov's schemes
System (2) can be rewritten in the following form that is more convenient for symmetry analysis 1

Remark 2.
Notice that for the polytropic gas with the state equation (3) one can rewrite the energy evolution equation (8e) of the Samarskiy-Popov scheme as In this form the energy evolution equation corresponds to equation (18e).
Calculations show [1] that the Lie algebra admitted by the system for an arbitrary σ = σ(p, ρ) is 1 The group generator is prolonged to the finite-difference space as follows [24,28] The scheme of the form defined on a uniform orthogonal mesh is invariant if the following criterion of invariance holds [28]X 1 Here and further the notation To preserve uniformness and orthogonality of the mesh it is also needed [24,28] where D ±τ and D ±s are finite-difference differentiation operators One can verify that scheme (8) is indeed invariant with respect to the generators X 1 , ..., X 4 and all the generators (20) satisfy the mesh orthogonality and uniformness conditions (25), (26). Hence, one can use an orthogonal uniform mesh (23b) which is invariant one.

Conservation laws possessed by Samarskiy-Popov's scheme
All the conservation laws of system (18) have their finite-difference counterparts for the Samarskiy-Popov scheme. They are given in Table 1.

An additional conservation law
only occurs in case σ = ρ [1]. In this case, system (18) admits two more symmetries, namely Conservation law (27) has its finite-difference counterpart which can be found by direct calculations.
Notice that neither conservation law (27) nor the center-of-mass law was mentioned in [13]. Perhaps, the authors of [13] have known the finite-difference analogue of (29).

The case of finite conductivity
We consider a more general case H = (H 0 , H y , H z ), E = (0, E y , E z ), i = (0, i y , i z ), u = (u, v, w), x = (x, y, z), and H 0 = const. Here we used the fact that the coordinate system can always be chosen in such a way that the first component of the vector field E is equal to zero.
Further we consider equations (1), where, by analogy with (18), the energy evolution equation (1h) is written as A generalization of scheme (8) for E = (0, E y , 0) and H = (H 0 , 0, H z ) is given in [13]. Since the MHD equations are almost symmetric in terms of the components E y , E z and H y , H z , one can extend the scheme proposed in [13] as follows where 0 (α, β 1 , β 2 ) 1 and Notice that this generalization of the scheme was discussed in [13] but it was not given explicitly.

1.
If H 0 = 0 and σ is arbitrary, then the admitted Lie algebra is In case σ = ρ, there are two additional generators are admitted, namely There are also two more conservation laws in the latter case (see Table 2). Additional conservation laws do not occur for any other forms of the function σ.

2.
If H 0 = 0 and σ is arbitrary then the admitted Lie algebra is where h 1 and h 2 are arbitrary functions of s. Additional conservation laws do not occur for any specific σ.
In both the cases above, scheme (31) is invariant. The rotation generator X 5 is only admitted for β 1 = β 2 . The remaining generators are admitted by the scheme for any set of parameters α, β 1 , β 2 and γ > 1.
The conservation laws possessed by system (31) and their finite-difference counterparts are given in Table 2. Here and further, conservation laws whose fluxes vanish for H 0 = 0 are marked with †. In case H 0 = 0, their densities preserve along the pathlines.

The case of infinite conductivity σ → ∞
In this case, system (1) reduces to where the internal energy is given by (3).
In addition to the analogues of conservation laws presented in the previous section, system (36) possesses the conservation law of angular momentum, namely As σ → ∞, scheme (31) becomes Center of mass law Energy conservation One can verify that scheme (38) is invariant one. As the symmetries of (36) and the corresponding difference schemes are reviewed in Section 3.2.3, we defer our discussion until then.

Conservation of angular momentum and energy
Apparently, the latter scheme does not preserve angular momentum, i.e. it does not possess a difference analogue of the conservation law (37). One can verify it by algebraic manipulations with the scheme or with the help of the finite-difference analogue of the direct method [37,38]. We overcome this issue by modifying the latter scheme as follows This allows one to obtain the whole set of finite-difference analogues of the conservation laws of equation (36) excluding the conservation of the entropy along the pathlines. The conservation laws are presented in Table 3. Notice that the three-layer conservation law of energy given in the table can be rewritten in the following two-layer form by means of (78c) Also, in order to verify the conservation law (37), one has to consider the following equations which can be obtained by integration of (78c) We also notice that the modified scheme (39) is still invariant and a completely conservative one.

Conservation of the entropy along the pathlines
This represents the conservation of the entropy S along pathlines which is a crucial difference between the finite and infinite conductivity cases.
It is known [42] that the Samarskiy-Popov scheme for polytropic gas does not preserve the entropy S for arbitrary γ. However, the following relation holds on solutions of the system which approximates the differential relation The latter relation holds along trajectories of the particles up to O(τ) for α = 0.5 or up to O(τ 2 ) for α = 0.5. In [42], an entropy preserving invariant scheme for gas dynamics equations in case of polytropic gas with γ = 3 was proposed. This scheme conserves the entropy along the pathlines but has only one conservation law, namely the conservation law of mass. It seems that the conservation of entropy by the difference scheme usually leads to the "loss" of some other conservation laws.
Here we propose a way of preserving the entropy along the pathlines for polytropic gas with integer values of adiabatic exponent γ 2 for scheme (39). We show that this can be done by choosing appropriate approximations of the state equation (3).
We notice that by means of (3) and (36a), equation (36g) can be represented as the identity In the finite-difference case, the rules of differentiation are different. As a result, not every approximation of the latter identity is a finite-difference identity. For a proper discrete analogue of (45) the right hand side of the identity should also be expressed in the divergent form as well as the left hand side. Choosing the difference approximation for the scheme in the case of a polytropic gas, one has an additional "degree of freedom": the choice of approximation for the state equation (3). This should be done so that both the left and right hand sides of the resulting approximation for (45) are divergent expressions. Notice that this does not affect the conservativeness of the total energy conservation law equation since it does not depend on any specific form of the equation of state.
Further we consider the shifted version of the equation (38d) First, we choose the following approximation of the state equation (3) for γ = 2, Substituting (47) into (46), one derives Solving with resect top (α) , one getsp The latter equation can be rewritten as This means conservation of entropy S along pathlines for γ = 2 on two time layers. We have achieved the integrability of the difference analogue of equation (36g) by choosing a suitable approximation for the state equation.
In a similar way one can arrive at the conservation of entropy for γ = 3, namely Similarly, for γ = 4 etc. Thus, by induction, one establishes the following general formula for an arbitrary natural γ 2 Entropy preservation formula (54) are presented in Table 3 among the other conservation laws.

Remark 4.
From the preservation of entropy S t = 0 in the differential case it follows (for simplicity, we consider the specific case γ = 2) Since the constant can be omitted, this means p(0, s) In the finite-difference case, by means of (51), one derives the following analogue of (55) where N = T/τ and we recall that p (α) = αp + (1 − α)p. Similar to the differential case, the latter gives which means entropy preservation for a given liquid particle.
Remark 5. The approach described above can also lead to entropy conservation for rational values of γ. Without proof of the existence of a general formula we present the result for γ = 5/3 which occurs for one-atomic ideal gas. One can verify that for γ = 5/3 the approximation for the internal energy ε leads to the following preservation of entropy 2 3p Center of mass law In case H 0 = 0, the admitted Lie algebra is In case γ = 2, there are two additional generators, namely Here q 3 , q 4 and q 5 are arbitrary functions of s. Scheme (39), (54) admits all the generators (64). However, the scheme does not admit the generators X 9 and X 10 .
There are the following additional conservation laws for system (36).

a)
In case H 0 = 0, there is an additional conservation law which corresponds to the generator ∂ s u ρ • The conservation law corresponding to the generator ∂ s is The latter follows from system (36). When conductivity of the medium tends to infinity, the phenomenon of frozen-in magnetic field is observed (see, e.g. [46]). In this case, in the absence of the longitudinal component H 0 of the magnetic field, the quantity B 0 which is proportional to the magnetic pressure turns out to preserve along the pathlines. • In case γ = 2, the admitted generator corresponds to the conservation law which follows from system (36). (70) is a basis one. Its partial derivative with respect to s is equivalent to (67), and its partial derivative with respect to t is
By virtue of the content of Remark 6, one can verify that the finite-difference analogues of (68) and (71) hold along the pathlines for scheme (39), namely and (74) Scheme (39) also admits the generators ∂ s and (69) under the same conditions as for the differential case.
Analyzing scheme (39), one can conclude that for the additional conservation laws (66), (67) and (70) there is no approximations in terms of rational expressions. This means that construction of finite difference analogues of the mentioned conservation laws is extremely hard.
Further we restrict ourselves to the case γ = 2 and S = S 1 = const, and consider another invariant scheme on an extended finite difference stencil.
We introduce the pressure for the polytropic gas as Then, the conservation law of entropy is defined by the following invariant expression The scheme under consideration is based on scheme (39) and it has the following form One can verify that the latter scheme is invariant. It admits the same symmetries as scheme (39), (54).
The following quantities hold for (78) Scheme (78) possesses the difference analogues of (66) and (67), namely To construct the latter conservation law one should use the following relation Remark 9. The angular momentum and center of mass conservation laws are The remaining conservation laws of mass, momentum, magnetic flux and entropy follow directly from the scheme since it is written in a divergent form.

Numerical experiments
In this section, we consider the problem of deceleration of a plasma bunch in a crossed electromagnetic field under the presence and absence of a longitudinal component of magnetic field. We use scheme (31), and consider how the conservation laws hold on the solutions of this scheme. In addition to the transverse component H y of the magnetic field, we also consider the case of the presence of a longitudinal magnetic field H 0 = 0.
A plasma bunch is considered, which moves from left to right in a railgun channel. The channel is filled with a relatively cold weakly conducting gas. With the help of an external electric circuit, a strong transverse magnetic field is generated in the channel, which causes the bunch to decelerate. During its motion, the plasma bunch closes the electric circuit and moves along the background gas; therefore, the magnetic field and pressure at the left boundary of the computational domain are considered equal to zero. The differential boundary conditions are as follows where 0 s S and 0 t t max , S is the total mass of the gas, J and V are current and voltage, and C 0 , L 0 , R 0 are the external circuit parameters. The boundary conditions (88c) and (88d) are approximated in the same way as in [13], namely where M = S/h . All calculations were carried out using the dimensionless version of scheme (31) with the value of the coefficient κ = 4π. For the dimensionless form of the scheme, the initial conditions are: ρ 0 = 1.0, p 0 = 0.0056, R 0 = 1.17, C 0 = 1.64, L 0 = 0.0035, S = 4.0, the temperature of the plasma T 0 = 3.0, the initial speed of the plasma bunch u 0 = 0.75. The gas is considered polytropic with γ = 5/3. The uniform mesh steps are h = 0.067 and τ = 0.003, and t max = 0.7. The initial voltage V 0 is varied between 1.67 and 2.6 which approximately correspond to the voltage 650 and 1000 V. In experiments where the longitudinal magnetic field H 0 is present, a value close to 1 is taken for H 0 . In the calculations, a linear artificial viscosity is used, with a viscosity coefficient ν = 2h.
The problem under consideration is close to the problem described in [47] in which, however, tabulated real plasma parameters, including electrical conductivity, were used. In our problem, we used the ideal gas equation and an exponential conductivity function.
Scheme (31) is implemented using the iterative methods described in [13]. In this case, the scheme equations are divided into two parts, dynamic and magnetic. The dynamic part is preliminarily linearized using the Newton method, and for the magnetic part a flow version of the sweep method is used [48], which is well suited for the case of finite conductivity, especially when its values are small. The bunch motion is modeled by a shock wave. Conductivity σ of the plasma bunch is proportional to T 3/2 , and the conductivity function σ = σ(ρ, T) is very sensitive to the density ρ in such a way that in the rarefied background gas region it has values close to zero.
Three essentially different cases are considered: 1.
The bunch is decelerated using a transverse magnetic field H y at a relatively low voltage in the circuit.

2.
The bunch is decelerated using a transverse magnetic field at a high voltage in the circuit.

3.
A rather strong longitudinal magnetic field H 0 is added to the previous case. (Calculations show that a weak longitudinal magnetic field has little effect on the experimental results.) In all cases, at the initial moment of time, the gas particles are given a small constant transverse velocity v > 0. This is necessary in order to track the influence of the longitudinal magnetic field on the transverse component of the particle velocity, which should be observed only in the third numerical experiment. Figure 2 shows the evolution of the magnetic field and plasma temperature in the first experiment. The magnetic field is not strong enough to stop the bunch. If the bunch reaches the right boundary of the computational domain, the reflection of the wave can be observed due to the boundary condition u(S, t) = 0. Figure 3 shows the second case where the transverse magnetic field is strong enough. The plasma bunch is decelerated by the magnetic field and after a short period of time begins to move backward. Adding a sufficiently strong longitudinal magnetic field H 0 to the previous experiment leads to an intermediate picture: the magnetic field is "smeared" over the computational domain, the plasma deceleration process is not as intense as in the previous case, and is inhomogeneous along the mass coordinate, which leads to a kind of fragmentation of the temperature profile (see Figure 4).

Conclusion
Finite-difference schemes for MHD equations in the case of plane one-dimensional flows are considered. The Samarsky-Popov classical scheme for the case of finite conductivity is taken as a starting point. Symmetries and conservation laws of this scheme are investigated. It is shown that the scheme admits the same symmetries as the original differential model. It also has difference analogues of the conservation laws of the original model. In addition to the conservation laws previously known for the scheme, new conservation laws are given, which are obtained on the basis of the group classification recently carried out in [1].
The classical Samarskiy-Popov scheme is generalized to the case of arbitrary vectors of electric and magnetic fields, as well as to the case of infinite conductivity. In the case of finite conductivity the scheme possesses difference analogues of all differential local conservation laws obtained in [1], some of which were not previously known. In the case of infinite conductivity, straightforward generalization of the scheme leads to a scheme that does not preserve angular momentum. The proposed modification makes it possible to obtain an invariant scheme that also possesses the conservation law of angular momentum. In addition, it is shown how to approximate the equation of state for a polytropic gas to preserve the entropy along the pathlines on the extended stencil for two time layers.
A numerical implementation of the generalized Samarsky-Popov scheme for the case of finite conductivity is performed for the problem of deceleration of a plasma bunch by crossed electromagnetic fields. Various cases of the action of fields on a plasma are considered. Calculations show that the finite-difference conservation laws are preserved on the solutions of the scheme quite accurately. Funding: This research was supported by Russian Science Foundation Grant No. 18-11-00238 "Hydrodynamicstype equations: symmetries, conservation laws, invariant difference schemes".

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.