Numerical simulation of shock-free strong compression of 1D gas layer’s problem subject to conditions on characteristic

The article describes results of one-dimensional gas layers strong compression problem’s math modeling. The article provided numerical algorithm for solving the one-dimensional ideal gas shock-free strong compression problem in R. Mises configuration is proposed too. The method combines finite-difference method ROMB and tracking feature method. The method allows to calculate gas-dynamic characteristic (velocity, density, etc.) of ideal gas layer while time increase and provides better accuracy in comparison with other finite-difference method. The accuracy of the proposed method was demonstrated in calculations of test plane-symmetry problem. Exact solution and numerical one agree quite well. Numerical results of solving one-dimensional problems with different symmetry and gas characteristic is also shown. The main results of numerical simulations are shown in graphs and tables.


Introduction
The mathematical description of shock-free compression of ideal gas to any predefined density including infinite is of interest for controlled thermonuclear fusion [1,2]. The idea of shock-free compression is very attractive as a way of attaining arbitrary high densities at minimal energy inputs due to the absence of shocks in the flow, i.e., unchanged initial entropy [3]. Attractiveness of the idea is seen from lots of publications on the topic (a comprehensive bibliography is provided in [4]). The mathematical theory of shock-free strong gas compression is developed in [4,5]. In particular, for cylindrical (ν = 1) and spherically (ν = 2) symmetric layers of a polytropic gas with γ > 1 , it is shown that the continuous joining of two flows gives a solution to shockfree strong compression of gas with a nonzero mass to any predefined density. The authors of [6] propose an algorithm which in reverse time simulates shock-free strong compression of 1D layers of initially resting homogeneous gas with initial density ρ 0 = 1 to any predefined finite and constant density ρ * > 1 and then reconstructs, again in reverse time, the path of the gas compressing piston. Reference [7] presents calculations with the algorithm for a number of 1D problems with different symmetries, where the gas layers are compressed externally, i.e., the piston radius decreases. It should be mentioned particularly that the calculations, including the reconstruction of the compressing piston path, are done in reverse time. In [6,7], the question of whether the reconstructed path implements when time increases is not discussed at all. The algorithm relies essentially upon the local theorem of existence of a solution in the vicinity of a singular point (t = t * , r = r * ) [4,5]. For reconstructing the piston path in case of large masses of gas being compressed, it is necessary to depart rather far from (t = t * , r = r * ) and the question  [6,7] but modified for a non-decreasing piston radius. It presents direct-time calculations for the problem of interest by the ROMB finite difference scheme [9], using the piston path reconstructed in reverse-time calculation. That is, the calculations explicitly simulate the action of compression piston on resting gas that at last time gives the desired gas state -the required density in the entire mass. The direct-time ROMB calculations showed up the flow regions where the algorithm distorts solutions including the known exact solution for plane symmetry. Near the weak discontinuity the solution smears, i.e., the compression wave runs away from the exact solution. This often happens when finite difference methods are used to simulate compression and rarefaction waves. The calculation region can be reduced if take the right boundary not constant (r = r * ) but movable. In [10] another algorithm is proposed, which uses the known advantages of the ROMB scheme and the known data on the path of the weak discontinuity between gas at rest and gas in motion. The algorithm helps improve numerical accuracy near the characteristic which separates the compression wave and the background flow. However the solution of the problem under study is known to have a singularity at point (t = t * , r = r * ) which in calculation manifests itself in the following way. While the overall gradient in the solution grows and reaches about 30, the gradient of the sought function near the boundary sound characteristic becomes as large as to prohibit further calculation. The study completes with derivation of analytic formulas for the Generalized Centered Wave (GCW) characterize and their use for calculation in the region between the compression piston and the boundary sound characteristic C + [11]. With the resulted analytic formulas the calculation region at times close to t = t * can be reduced to a portion where singularities in the solution are absent, i.e., to the region between the compression piston and the boundary GCW characteristic where the desired values of density are reached. This helped us accurately calculate compression of 1D gas layers to t = t * .

Problem statement
Consider 1D symmetric isentropic (s = const) flows of a polytropic gas. Such flows are solutions of equations [4]: Here c = ρ (γ−1)/2 is sound velocity in gas; ρ is density; γ > 1 is adiabatic exponent in the u is a projection of vector V onto the radius vector in the plane x 1 Ox 2 at ν = 1 or the radius vector in the space of variables Let a background flow c(t, r) = 1, u(t, r) = 0 be defined in a vicinity of point (t = t * , r = r * ) for r * > 0. In Figure 1 it is shown in region 0, and in region 1 we are to find a flow with a predefined distribution of one of the hydrodynamic parameters which will be mated to the background flow via the weak discontinuity.
The weak discontinuity is the sound C + characteristic whose path with the above background values is defined by the equation: The region of the sought flow is limited by the straight line t = t * and the sound characteristic C + 0 of the background flow. It consists of two parts: a lower triangle (denoted 1 in Figure 2) where the Riemann generalized centered wave is defined, and an upper triangle (denoted 2 in Figure 2) where the flow with a specified constant sound velocity at t = t * is defined. We need to find the law describing the path of a piston which with no shock compresses a resting homogeneous gas of density ρ = 1 to density ρ = ρ(t * , r) = ρ * = const, which remains constant within a compressed layer of width d * , (Figure 2). The path is defined as a solution of the following Cauchy problem: Here the function u(t, r) is found from the solution of system (1). The theorems proved in [4] assert that there is a nonzero mass m 0 between the piston path r = r p (t) and the line r = r * , which can be compressed shock-free to any predefined density ρ * . But the theorems do not determine the maximal value of m 0 and hence the maximal width d * (Figure 2).

The reverse-time calculation method
The problem is numerically solved with an algorithm which exploits characteristic's method and relies upon results obtained in [4,6,12]. The main difference from the classical method of characteristics is that the solution is sought in reverse time and the singular point (r * , t * ) is treated by creating not a single but a bundle of sound characteristics C + from it. Figure 3 a schematically shows how the characteristic grid in region 1 is calculated. A detailed description of the algorithm can be found in [6]. The path is reconstructed in reverse time from the point (t = t * , r = r * − d * ) where d * is the width of the compressed layer which defines its mass. Calculation with the algorithm from [6] runs until the piston path either crosses the characteristic C + 0 , or escapes from the region where the solution exists.

The direct-time calculation method
Consider the following 1D system of hydrodynamic equations in Lagrangian mass coordinates: here t, r, m, ρ, v = 1 ρ , ε, p are, respectively, time, mass, density, specific volume, specific internal energy, pressure, and E = ε + u 2 /2 is total specific energy.
The system is closed with the equation of state for ideal gas p = (γ − 1) · ρC v T, ε = C v T, where C v is specific heat at constant volume and T is gas temperature. System (2) is completed with the boundary and initial conditions u| r=r L (t) = u p (t), where r L , r R are the left and right boundaries of the section where system (2) is being solved. At t = 0: r L = r p (0),r R = r * ,. Here u p (t) is the velocity of the compression piston from the reconstruction of its path in reverse time.
System (2), (3) is solved with the ROMB scheme which is described in rather detail in [9]: pressure and Newton specific internal energy are linearized and iterations start. New variables (pressures P and velocities U at mesh nodes) are introduced through relating equations. After appropriate manipulation we come to a system of linear algebraic equations for P and U with a four-diagonal matrix which is solved with a method similar to stream sweeping [13].

Analytic representation of GCW characteristics
The analytic formulas for the sound GCW characteristics are derived under certain assumptions. It was obtained in [11] with help of analitical method from [14]. Below are the ordinary differential equations (ODE) that are solved to obtain the functional dependencies and the final formulas.

Calculations
The reverse-time calculation algorithm and results of calculations according to the algorithm can be found in [8]. In this paper we only use the trajectories of compression piston was calculated with help of reverse-time method. Now look at direct-time calculations. Figure 4 show results of calculation by the ROMB scheme without the explicit tracking of the weak discontinuity. It is seen that the ROMB scheme adequately reproduces the compression region but significantly smears the weak discontinuity. Of special note is a flat section seen on the left in the profiles of hydrodynamic parameters. It proves that the piston path reconstructed in reverse-time calculation does compress gas to the required density.
A test calculation was done to compare calculations by ROMB without singularities and by a hybrid method which calculates only the region between the compression piston and the C 0 characteristic. It is seen from Figure 5 that the latter describes the solution in the compression wave and weak discontinuity regions very much better than ROMB.
Calculations were done for cylindrical and spherical symmetries with different adiabatic exponents and initial masses. ROMB conservatism ensures energy conservation at the end of calculation at level close to computer zero. Energy conservation at each time step tn is checked by calculating energy disbalance with the formula where m is mass, E n k and E 0 k are specific total energies at current and initial times, and W is the work done by the compression piston. Data on energy disbalance, calculated according to (9), at different times from one of the calculations are provided in Table 1. Total energy disbalance is a result of inaccuracy in the determination of the current piston velocity through linear interpolation. The data of Table 1 show it to remain insignificant and almost unchanged through the entire calculation.
Below are figures that demonstrate good accuracy of the modified method, and in the region where the compression wave and the resting gas join, too. Figure 6 shows calculated results for a cylindrical symmetric problem with γ = 1.4, ν = 1, m * = 10. The solution obtained in reverse- time calculation by the method of characteristics with recalculation is used as a reference. It should be noted that the region where the compression wave joints the resting gas is reproduced rather accurately.
The modified method and analytic formulas (7), (8) helped obtain in direct-time calculations the hydrodynamic parameters of a gas layer up to and inclusive t = t * , in the region between the compression piston and the boundary GCW characteristic. Figure 7 shows their profiles at t = t * , for γ = 1.4, ν = 2, ρ * = 10 5 , m * = 10. The solution is seen to agree well with the reference.

Conclusion
The results presented above allow the following inferences. 1. The law describing the motion of a piston with a non-vanishing radius, which compresses a 1D layer of gas of a specified mass to a required density, was obtained in the reverse-time calculation. The algorithm proposed for solving such problems was tested. Reverse-time calculations allow the piston path to be reconstructed within 1-2 minutes on standard PC. 2. With the known r = r p (t), direct-time calculations were done for a number of problems by the ROMB method. Their results show that the numerical algorithm reconstructs the solution with good accuracy up to time 0.999t * regardless of problem symmetry. Further calculations show that calculation accuracy degrades near the singular point where derivatives on the boundary sound characteristic start rapidly grow. The effect can be reduced by moving off the right boundary of the calculation region.
3. The basic error of calculation is confined within a rather narrow region about the weak discontinuity that separates the compression wave and the resting gas.
4. The implemented ROMB modification allows attaining higher accuracy of calculation up to t = t * .
5. The approach proposed not only helps reconstruct the path of compression piston for an arbitrary 1D gas layer, but also proves numerically that the piston does compress a specified mass of gas (including m * = 10) to a desired density (including ρ * = 10 5 ). Figure 8 shows laws of compressing piston motion in time and it's work. These results therefore provide specific recommendations for physical experiment.