The piecewise parabolic method for Riemann problems in nonlinear elasticity

We present the application of Harten-Lax-van Leer (HLL)-type solvers on Riemann problems in nonlinear elasticity which undergoes high-load conditions. In particular, the HLLD (“D” denotes Discontinuities) Riemann solver is proved to have better robustness and efficiency for resolving complex nonlinear wave structures compared with the HLL and HLLC (“C” denotes Contact) solvers, especially in the shock-tube problem including more than five waves. Also, Godunov finite volume scheme is extended to higher order of accuracy by means of piecewise parabolic method (PPM), which could be used with HLL-type solvers and employed to construct the fluxes. Moreover, in the case of multi material components, level set algorithm is applied to track the interface between different materials, while the interaction of interfaces is realized through HLLD Riemann solver combined with modified ghost method. As seen from the results of both the solid/solid “stick” problem with the same material at the two sides of contact interface and the solid/solid “slip” problem with different materials at the two sides, this scheme composed of HLLD solver, PPM and level set algorithm can capture the material interface effectively and suppress spurious oscillations therein significantly.

Nonlinear elastic deformation of solid material undergoing high-load conditions commonly occurs in industrial application areas, such as the design of automobile anti-collision device, the evaluation of the capability of spacecraft structural materials against hypervelocity impact, and etc. For such a phenomena of technological interest, the research strategy can be classified into three categories: Lagrangian view 1,2 , Eulerian view [3][4][5] and the combination of the above two views 6 . The Lagrangian view is more popular in engineering area as the mesh point for this view is always coincident with material point in the deformation process of material and no material point crosses the mesh boundary, leading to no necessity for processing material interface. While the material has a large deformation, the computational mesh will warp or distort. Although this problem can be overcome through mesh remeshing and deleting the abnormal meshes, the effects on the computational accuracy cannot be ignored yet. While for Eulerian view, fixed meshes are adopted and materials are allowed to pass through mesh interface, implying that mesh warp or distortion induced in the case of large deformation could be avoided in a large extent. Thus, it is relatively suitable for studying the problems featured with large material deformation. Despite of this, it is still a great challenge to apply the Eulerian reference frame as the capturing of the interfaces between different materials is quite complicated therein.
Elastic deformation of solid material is a common phenomenon in engineering, but to describe large deformation process 7 accurately is still a tough task. In recent years, a variety of model equations and numerical methods have been developed under Eulerian view to simulate nonlinear elastic deformation behaviors of solid material, as reviewed by Benson 8 ,Wilkins 9 . With regard to model equation, the ones of elastic strain law of materials should be added in Eulerian reference frame on the basis of the governing equations in Lagrangian reference frame, which are composed of mass, momentum and energy equations. The model equations describing elastic strain law can be divided into two major groups. The first one, whose variable is deviator stress tensor [10][11][12] , is constructed in non-conservative form and thermodynamically inconsistent. It is extensively applied in engineering as it could also account for plastic effect of solid material through Maxwell relaxation model, despite that its disadvantages of not ensuring the strict numerical resolution due to the non-conservative form and being unable to

Results
In this section, we first present a comparative study of the HLL-type Riemann solvers which are applied in single material cases. (Particularly, we select test cases that are known to cause severe difficulties for numerical computations to assess the accuracy and robustness of the schemes.) With the state equation given in Equation (10) employed, Equation (8) is solved in each test case under different initial and boundary condition. The common material parameters, which appear in state equation, are shown in Table 1. Moreover, first-order Godunov's scheme is extended to PPM, which is featured with higher order of accuracy and used with different types of Riemann solvers for single material problem. Further, the validation of PPM combined with different types of Riemann solvers (especially the HLLD solver) are examined in one or two dimensional cases with multi materials by numerical experiments. In the test cases below, all one-dimensional initial value problems are solved in a computational domain 0 m ≤ x ≤ 1 m, and the initial discontinuity is placed at the location of x 0 = 0.5 m. Transmissive boundary conditions are applied at the computational domain boundaries and the CFL number is set to be 0.8 unless otherwise specified.
Test case 1: Contact discontinuity problem. In this case, we solve the Riemann problem with initial density and velocity distributions corresponding to a contact discontinuity. At both sides of this discontinuity, the materials are both copper and the initial distributions are as follows: Here we solve the case within the dimensional time t = 0.7 ms on a mesh of 500 cells and make a quantitative comparison with the exact solution W(x, t) = W(x − 0.01t, 0). Figure 1 presents the profiles of density ρ and velocity component u computed by first-order Godunov's scheme with three kinds of Riemann solvers, i.e. HLL, HLLC and HLLD (In the following we denote the methods of first-order Godunov's scheme combined with HLL, HLLC and HLLD Riemann solvers by '1st + HLL' , '1st + HLLC' and '1st + HLLD' , respectively.) It is clearly shown that the HLLC and HLLD solvers produce nearly accurate density profiles which are almost coincident, while the density profile derived by HLL is rather different from the exact solution. For the velocity distribution, the results from three methods, among which only that from HLLD is nearly coincident with the exact solution, are considerably different. The L 1 norm errors and convergence orders of density ρ and deformation gradient F 11 at time t = 0.7 ms for test case 1 are presented in Table 2 for 1st + HLL, 1st + HLLC and 1st + HLLD. It is observed that all first-order methods are converging at the order of Test case 2: Five-wave shock-tube problem. In the present case, the materials at both sides of the discontinuity interface still are copper. The initial distributions satisfy compatibility condition and are given as  For this case, the solution structure comprises five waves (from left to right): two left-travelling rarefaction waves, a right-travelling contact discontinuity wave, a right-travelling rarefaction wave and a right-travelling shock wave, respectively. The solution derived by high-order scheme with a fine mesh is used as the reference solution for comparison.  We solve this problem within the dimensional time t = 0.06 ms on a mesh of 500 cells by using first order scheme with HLL, HLLC and HLLD solvers. The density profiles are shown in Fig. 2, where the exact solution is denoted by solid line. It is seen that HLL has the most diffusive result in the region 0.4 m ≤ x ≤ 0.6 m, where HLLC and HLLD provide more accurate results. Also, it is found that HLLD performs better than HLLC as shown in the local enlarged graph Fig. 2(b). To check grid convergence, the cell number is increased from 500 to 1000 and 2000 with the same time interval and the results from different gird resolutions are compared in Fig. 3. As shown in Fig. 3(a), the convergence rate is relatively low for HLL, and the density profile with 2000 grid cells is even very diffusive around x = 0.5 m. With increasing cell number, the density profiles derived from HLLC and HLLD are improved in the vicinity of x = 0.5 m, although the convergence rates are still very low due to their low-order nature. Relative CPU-times, which are the corresponding multiples relative to the CPU time taken by first order scheme coupled with HLL with the coarsest grid (N = 500), are given in Table 3 for different combinations of computational scheme and Riemann solver. While first order Godunov's scheme is used, the CPU time spent by HLLC is approximately 1.7 times as long as that by HLL, and the time spent by HLLD is approximately 2.2 times as long as that by HLL. In summary, the first-order scheme can only give relatively reasonable solution even with very fine meshes and the use of HLLC and HLLD, having low convergence rate.
In the present framework, a simple way to increase the resolution of shock and discontinuity is to employ high-order reconstruction. Herein, we accomplish such reconstruction by applying Piecewise Parabolic Method (PPM), which is coupled with Riemann solvers, including HLL, HLLC and HLLD. Figure 4 plots the profiles of density derived from different combinations of PPM and Riemann solver: the one of PPM and HLL (denoted as 'PPM + HLL'), the one of PPM and HLLC (denoted as 'PPM + HLLC') and the one of PPM and HLLD (denoted as 'PPM + HLLD'). As shown in Fig. 4(a), the profiles can be resolved reasonably using a mesh containing only 100 cells. When the number of cells is increased to 500, the contact discontinuities, shock waves and the rarefaction wave in the region of 0.1 m ≤ x ≤ 0.2 m are all captured accurately, as shown in Fig. 4(b-d). Compared with exact solution, the scheme PPM + HLLD is seen to have the most accurate solution, particularly in the region of 0.4 m ≤ x ≤ 0.5 m where the solution exhibits very slow convergence rate. The CPU time spent by PPM + HLLC is approximately 1.3 times as long as that by PPM + HLL, while the time spent by PPM + HLLD is approximately 1.8 times as long as that by PPM + HLL. The L 1 -norm errors at time t = 0.06 ms and convergence rates for all combinations of PPM and Riemann solver are presented in Table 4 for test case 2.  in the five-wave shock-tube example, we further test the accuracy and robustness of the high-order scheme PPM + HLLD in solving a more complex problem, i.e. the seven-wave shock-tube problem, by adding an additional degree of shear deformation. The initial conditions of this problem are The solution structure for this case is composed of three left-travelling rarefaction waves, a right-travelling contact discontinuity, two right-travelling rarefaction waves and a right-travelling shock wave. The results shown in Fig. 5 are obtained by the scheme PPM + HLLD on a mesh of 500 cells with the dimensional computational time t = 0.06 ms. As shown in Fig. 5, all the seven waves are captured precisely, implying that the PPM + HLLD scheme indeed has superiority and good robustness in processing the solid material problem with complex multi-wave structure. CPU times for test case 3 are comparable to those for test case 2 shown in Table 3. The L 1 -norm errors at time t = 0.06 ms and convergence rates of PPM + HLLD are presented in Table 5 for test case 3.

Test case 4: Solid/solid 'stick' problem.
With the purpose of examining the performance of the method PPM + HLLD in the multi material problem, we will solve the solid/solid 'stick' problem and solid/solid 'slip' problem in the following two cases (test case 4 and 5) by using it. For the former, the materials at two sides   Table 3. Computation times for test case 2. CPU times are shown relative to the first order scheme coupled with HLL Riemann solver on the coarsest grid.
of contact interface where the 'stick' condition is satisfied are the same and thus the Riemann solvers for both single-material and multi-material situations can be utilized for solution. When multi-material HLLD solver is used, level-set method for interface tracking must be applied in conjunction. This technique of PPM combined with multi-material HLLD solver and level-set method is referred to as 'PPM + HLLD(M)' . And, the one of PPM coupled with single-material HLLD solver is denoted as 'PPM + HLLD(S)' . The initial conditions for test case 4 are    Table 4. L 1 norm errors and orders of convergence for test case 2.
In the solving process, the CFL number is set to be 0.6 and a mesh of 1000 cells is used in the computational domain. The profiles of density and entropy are shown in Fig. 6 at time t = 0.06 ms. Compared with the exact solution, the ones from PPM + HLLD(S) and PPM + HLLD(M) techniques are both featured by sharp peaks across shocks. Further, it is also found that spurious overshoots occur in both density and entropy profiles in the vicinity of x = 0.55, due to the fact that variables are not conserved across linearly degenerate field. Compared with PPM + HLLD(S), PPM + HLLD(M) can effectively suppress the spurious overshoot phenomenon due to the use of entropy fix technique, as illustrated in the Fig. 6(b,d). Nevertheless, it is noteworthy that this kind of error could not be eliminated completely in the present method framework, which still needs to be improved further in future. L 1 norm error and its convergence order for selected variables are shown in Table 6 with different cell numbers and method combinations such as 1st + HLLD, PPM + HLLD(S) and PPM + HLLD(M). It is found that L 1 error norm for primitive variables exhibits approximately first-order convergence trend for both PPM + HLLD(S) and PPM + HLLD(M) techniques on account of the discontinuities (longitudinal shock, contact discontinuity and transverse shock) present in the solution. Further, the errors for all schemes decrease with increasing cell number, implying that the solutions are converging. For tangential velocity u 2 and deformation gradient F 31 , the errors obtained by applying PPM on the coarsest mesh are lower than those by 1st-order Godunov's scheme on the finest mesh, which indicates PPM indeed has higher order of reconstruction accuracy once again. In summary, the results for this case demonstrate that the material interface could be tracked accurately by level-set method, and highly precise solutions could be obtained by applying PPM + HLLD(M) with the spurious overshoots being suppressed simultaneously.
Test case 5: solid/solid 'slip' multi-material problem. The initial conditions for this case are different from that for test case 5. In detail, the material at the left side of interface is now aluminum, while the one at the right side is still copper, with the slip condition satisfied at the material interface in the middle. The solution structure of this problem comprises six nonlinear waves, which are three left-travelling shocks inside the aluminum   medium, three right-travelling shocks inside the copper medium, and a right-travelling contact discontinuity. We solve the problemm by means of the technique PPM + HLLD(M) for multi material situation within the dimensional time t = 0.05 ms. CFL number is fixed to be 0.6 and two meshes, including 100 cells and 1000 cells respectively, are used. Various variable profiles are shown in Fig. 7. With only 100 cells, the wave structure is captured reasonably by PPM + HLLD(M), although the numerical solution seems to be diffusive near the discontinuities.
As the cell number is increased to 1000, the profiles of density, stress components and velocities near shocks become sharp without significant spurious oscillations, implying that PPM + HLLD(M) could be applied in multi material problems successfully. The L 1 norm errors at time t = 0.05 ms and convergence rates of PPM + HLLD(M) are presented in Table 7 for test case 5.
Test case 6: impact of a projectile on a solid plate. In this two-dimensional case, the technique PPM + HLLD(M) for multi material situation is utilized to simulate the impact problem of a projectile on a solid plate surrounded by vacuum. The initial configuration is shown in Fig. 8.   initial time all materials are assumed to be in a stress free state: F = I and S = 0. Further, the cooper target is set to be static, while the copper projectile is initialized with a non-zero velocity component: u 1 = 800 m/s. The solution of this test case was performed using a similar Eulerian model as that developed by Favrie 34 and Gorsse 24 , where the projectile and plate are surrounded by air.
The schlieren images of the density are shown at time t = 1.0 × 10 −5 s, t = 2.0 × 10 −5 s, t = 2.3 × 10 −4 s and t = 6.5 × 10 −4 s in Fig. 9. It is found that elastic shock propagates away from the impact surface, reaches the free surface and is subsequently reflected to form a rarefaction wave. Further, the elastic material is deformed and its surface oscillates with time, having strong similarity to the phenomena depicted by Favrie 34 Fig. 10 for this test case, for which the differences between the results from different grid resolutions are nearly invisible. The L 1 norm errors for density ρ, normal velocity u 1 and deformation gradient F 11 at the output time t = 1 ms and the convergence rate of PPM + HLLD scheme are presented in Table 8, from which we could conclude that the results converge with the order of larger than 3, implying that PPM + HLLD is at least third-order accurate for smooth problems.

Discussion
In this study, the governing equations developed by Godunov & Romenski 14,31 are utilized to describe the elastic deformation of solid materials in Eulerian reference frame. The existing Godunov-type shock-capturing schemes have been applied in conjunction with HLL family of Riemann solver to solve Riemann problems in nonlinear elasticity. The two-state HLL Riemann solver 7 has been widely used as the standard shock-capturing scheme due to its simplicity and high effectiveness. However, numerical experiments show that HLL solver is of strong dissipation particularly in the cases with contact surfaces and strong shock waves, leading to inaccurate resolution of physical features and unacceptable numerical smearing. HLLC 35,36 method assumes a three-wave model and thus it could lead to better resolution of intermediate waves. Nevertheless, for systems with the eigenstructures containing more than three distinct characteristic fields, HLLC seems to be inadequate and tends to behave like HLL with inaccurate resolution of intermediate waves, particularly when these waves move slowly relative to the mesh used, as illustrated in test case 1. Therefrom, HLLC solver has been improved and developed to be HLLD solver by admitting the correct number of characteristic field 37,38 . HLLD solver involves five wave structure and has been successfully applied in MHD 39,40 . In this paper we extend the usage of HLLD solver to the impact problem of nonlinear elasticity, considering the similarity of the wave structures in nonlinear elasticity with those appearing in MHD. And, a comparative study shows that 1st-order Godunov's scheme coupled with HLLD solver has advantages in capturing multi-wave structures in solid material with large deformation, in comparison with that coupled with HLL and HLLC solvers.
With the purpose of increasing the convergence rate of solution, we apply the well-known piecewise parabolic method, an extension of 1st-order Godunov's scheme to higher order accuracy, and couple it with HLL-type Riemann solvers to solve the problems featured with complex nonlinear wave structures. The results of numerical experiments show that the technique PPM + HLLD provides solutions with highest resolution for single material problems compared with other techniques. Moreover, as HLLD Riemann solver can deal with boundary conditions at the interface between different solid materials, we develop the technique PPM + HLLD(M) by coupling PPM with multi-material HLLD Riemann solver for the treatment of multi material problems, where level-set algorithm is used for tracking material interfaces. As shown in the solid/solid 'stick' problem, the PPM + HLLD(M) technique in conjunction with level-set algorithm can suppress greatly 'overshoot' or 'heating errors' , which may occur in the vicinity of contact discontinuity as observed in the work of Titarev 20 . For the solid/solid 'slip' multi-material problem, it is demonstrated that the PPM + HLLD(M) technique is able to capture nonlinear wave structure accurately even in the presence of strong shocks. While for the two-dimensional impact problem (test case 6), the PPM + HLLD(M) technique can reproduce similar dynamic behaviors with those given by Favrie 34 and Gorsse 24 , proving the robustness of this technique in two-dimensional cases. Further,   compared with the method adopted by Barton et al. 30 on the basis of characteristic relation of invariants, the PPM + HLLD technique is relatively easier to be implemented in the existing code. The accuracy analysis on our PPM scheme confirms that it exhibits approximately third-order convergence and even achieves fourth-order accuracy with time step Δt tending to 0 for smooth problem. The relevant results of accuracy and error analysis imply that for the Riemann problem for elasticity with discontinuity, linearly degenerate discontinuous waves produce sub-linear convergence rates which eventually dominate the global convergence rates. In our examples, the convergence rates with first-order Godunov's scheme are close to 0.5, while those with PPM are close to 0.75, being consistent with the theoretical rates for contact discontinuities presented by J.W. Banks 41 , where the authors reported that for a general scheme of order p, the order of convergence for unlimited stable schemes is established as p/(p + 1). In summary, the PPM + HLLD technique has been proved to be a robust tool in solving Riemann problems in nonlinear elasticity for both single material problem and multi material problem, and it could be natural to be extended to plastic problems.

Methods
Governing equations of nonlinear elasticity. The model developed by Godunov & Romenski 14,31 is used to describe the deformation process of solid material in the Eulerian reference frame. The physical variables used contain velocities u i , deformation gradient tensor F ij = ∂x i /∂X j (where x i and X j denote the fixed Eulerian coordinates and Lagrangian coordinates respectively) and specific entropy S. Following the notations used by Titarev 20 and Barton 21 , a hyperbolic equation system depicting momentum, strain, and energy conservation laws in Cartesian coordinates can be written as Here ρ denotes the material density, σ is the Cauchy stress, E = ε + u i u i /2 is the total energy, ε is the specific internal energy and the Einstein summation convention over repeated indices is implied (i, j, k = 1, 2, 3). The material density, stress tensor, specific internal energy, strain tensor G ij and temperature T can be represented as functions of the variables mentioned above: where ρ 0 is the density of the initially unstressed medium. Further, the continuity equation is expressed as Following the treatment by Barton 42 , we also use the continuity equation (5) to replace the strain conversation equation for the deformation gradient component ρF 11 for convenience. Thereafter, the governing equations can be expressed in matrix form: where U is the vector composed of the conservative variables, F i is the corresponding flux vector and S C is the vector of source terms associated with the constraints for the deformation tensor: Introducing the vector of primitive variables W = (u, F T e 1 , F T e 2 , F T e 3 , S), we could rewrite Equation (6) as a hyperbolic quasi-linear system t x where the Jacobian matrix A k is given by Barton 21 for details on the eigenstructure of this system.
The equation of state (EOS) for elastic solid employed is the formula about the specific internal energy ε, which is expressed in terms of three independent invariants of the Finger tensor (I 1 , I 2 , I 3 ) 20,43 : Finite Volume Methods. Considering now the one-dimensional system of elastic solid and taking k = 1 in Equation (6), the following conservative equation is obtained: where the source term S c becomes equal to zero and is neglected in numerical computation for one-dimensional cases. Note that for two-dimensional and three-dimensional cases, the source item S c is not neglected and treated as source terms in the numerical computation for the purpose of capturing the correct wave speed in the quasi-linear system (Equation (8)) and ensuring the correctness of numerical results. Equation (6) is not conservative in nature for multi-dimensional case. Nevertheless, we follow standard practice 21,22 and solve the two-dimensional impact problem (test case 6) by using Equation (6) sequentially. With grid spacing Δx = x i+1/2 − x i−1/2 and time step Δt = t n+1 − t n , the finite volume scheme for solving this hyperbolic Equation (12) can be written as where U i n is an approximation to the average of spatial integral in the cell [x i+1/2 , x i−1/2 ] and F i+1/2 is the numerical flux yet to be defined. As proposed by Godunov 16 , the numerical flux F i+1/2 is derived by solving the Riemann problem with the initial data In this study, we use HLL (Harten, Lax and van Leer) family of solvers to obtain approximate solutions of the Riemann problem in nonlinear elasticity and carry out a comparative study on their robustness and accuracy. It is worthwhile that for HLL-type solvers, the wave propagation speed (λ k ) and the wave decomposition of Δq into Δ k q are not rigorously derived from the Jacobian matrix, leading to their usability for simulating the problems with complicated Jacobian matrix. In the following, we will briefly introduce HLLD Riemann solver. While to see the details of HLL and HLLC Riemann solver, please refer to Harten 7 , Gavrilyuk 4 and Toro 15 .
HLLD Riemann Solver. HLLD Riemann solver gives a nonlinear approximate solution. Its central idea is to assume a wave configuration for the solution that consists of five waves (two slow waves, two fast waves and a contact discontinuity), which separates six constant states. As shown in Fig. 11, there are four intermediate states: ∼ − U , U *− , U *+ , and ∼ + U . The fastest (longitudinal) waves between U ± and ∼ ± U are denoted S L ± and the slow shear waves S S ± separate the states U ∼ ± and U *± . Each wave is considered to be a discontinuity and the Rankine-Hugoniot relation is satisfied across each wave ( ± S S and ± S L ): From Eqs (15a-15c) one can find that there are more unknowns than equations, thus some other conditions must be imposed. In order to obtain the unknown intermediate state vectors U ∼ ± , ⁎± U , F ±  and ± F ⁎ , the following conditions need to be satisfied: • Tangential velocities u 2 , u 3 and tangential stresses σ 12 , σ 13 are continuous across fast (longitudinal) waves and may jump across slow (shear) waves. • Density ρ, normal velocity u 1 and normal stress σ 11 are continuous across slow waves and may jump across fast waves. • Normal stress σ 11 and normal velocity u 1 are continuous across contact discontinuity; for the 'stick' multi-material problem, shear stress and tangential velocity are equal for different materials at the interface; while for the 'slip' multi-material problem, tangential component of the stress vector are zero.
Data availability statement. The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.