The piecewise parabolic method for elastic-plastic flow in solids

A numerical technique of high-order piecewise parabolic method in combination of HLLD (”D” denotes Discontinuities) Riemann solver is developed for the numerical simulation of elastic-plastic flow. The introduction of the plastic effect is realized by decomposing the total deformation gradient tensor as the product of elastic and plastic deformation gradient tensors and adding plastic source term to the conservation law model equation with the variable of the elastic deformation gradient tensor. For the solution of the resulting inhomogeneous equation system, a temporal splitting strategy is adopted and a semi-implicit scheme is performed to solve the ODES in the plastic step, which is conducted to account for the contributions from plastic source terms. As seen from the results of test cases involving large deformation and high strain rate, the computational model used can reflect the characteristics of constitutive relation of material under strong impact action and our numerical method can realize the exact simulation of the elastic-plastic behavior of solid material, especially the accurate capture of the elastic-plastic waves. Further, it could also deal with high-speed impact problems with multi-material components, catching material interfaces correctly and keeping the interfaces sharp, when combined with interface tracking technique such as the level-set algorithm.

of plastic behavior of the material under impact action. However, this description way is only applicable to perfect plasticity and has some difficulty in the achievement of the reflection of elastic-plastic shock wave on free interfaces. Another way to depict plastic behavior based on hyperelastic model is to decompose the total deformation tensor (deformation gradient tensor or its inverse) into the two parts of elastic and plastic deformation by using multiplicative decomposition method and to realize the characterization of plastic effect by means of the additional source term. Plohr 12 had constructed both rate-independent and rate-dependent plasticity model based on the inverse of deformation gradient tensor and obtained a larger equation system with the total equation number of 21. And, Miller 13 and Barton 14 had taken the second way, i.e. added the plastic source term to the control equation for elasticity based on deformation gradient tensor, to realize the depiction of plastic effect.
Considering that elastic-plastic shock wave and discontinuousness exist in the elastic-plastic flow problem widely and the numerical method of Godunov type has been proven to have apparent superiority for capturing both shock wave and discontinuousness in fluid mechanics, it is natural to extend the use of high-order scheme of Godunov type based on the solution of local Riemann problem to solid mechanics in Eulerian framework, in case that the control equations are constructed to be of properly conservative form. Accordingly, several schemes of Godunov type, such as piecewise linear scheme, first-order Godunov scheme, MUSCL-Hancock scheme and MUSCL scheme, have been applied in the elastic-plastic deformation problems (see Miller et al. 13 ; Kluth et al. 11 ; Hank et al. 15 ; Favrie et al. 16 ; Howell et al. 1 ; Brauer et al. 4 ). Note that the above works are mostly focused on modifying the piecewise-constant distribution of original variables as the piecewise-linear distribution in space, achieving second-order accuracy of numerical computation. To our knowledge, the numerical scheme of Godunov type with higher order accuracy, such as Piecewise Parabolic Method (PPM) 17 which has third-order accuracy and stronger ability of processing discontinuousness, has not been utilized in the problem of elastic-plastic deformation.
In addition to scheme precision, the solution of local Riemann problem is another important and noticeable respect for the schemes of Godunov type. In the nonlinear hyperbolic system of fluid mechanics, it is commonly realized by linearizing the Jacobian matrix, e.g. Roe scheme 18 . While for solid mechanics, it is considerably difficult to derive the linearized Jacobian matrix analytically due to the complexity of equation of state. Thus, two kinds of numerical method for the computational treatment of this problem are raised. The first choice is to do the arithmetical average of Jacobian matrices at the left and right cell centre states simply 19 . However, this way failed to solve the Riemann problem when applied in the impact case of nonlinear elasticity 19 . The second choice is to make use of Harten-Lax-van Leer (HLL) solver based on two-wave assumption 20 and HLLC solver ('C' denotes Contact) based on three-wave assumption 21,22 . Howell et al. 1 had used HLL solver to simulate two-dimensional elastic-perfectly plastic deformation problem of solid material and the application of HLLC solver in the elastic-plastic deformation problem had been discussed in Ndanou et al. 23 , Cheng et al. 24 and Brauer et al. 4 . However, HLL and HLLC solvers based on two-wave and three-wave assumptions may be inadequate when applied to nonlinear elasticity, in view that the exact wave pattern of Riemann problem may involve more than five distinct waves therein. Thus, we have developed HLLD ('D' denotes Discontinuities) solver based five-wave assumption and found that it has obvious advantages in capturing wave patterns containing more than five waves when applied to nonlinear elasticity (see Zhang et al. 25 . In this paper, the numerical technique of PPM reconstruction combined with HLLD solver is chosen to simulate the elastic-plastic deformation behavior of solid materials, with the purpose of both the derivation of high order numerical accuracy and the accurate capture of multi-wave patterns.

Results
Test case 1: Accuracy test. To test the accuracy of our method in a plasticity-dominated problem, we perform the computation for a one-dimensional nontrivial and shockless problem similar to the purely elastic problem presented in Zhang 25 . Two initially Gaussian distributions with width 5 are used to initialize F 22 and F 33 as functions of coordinate x i : The other components of velocity u and F are all set to 0. The parameters for equation of state 15 and constitutive model are as follows: ρ 0 = 2.7 g/cm 3 , p 01 = 73 GPa, p 02 = 172 GPa, p 03 = 40 GPa, μ = 24.8 GPa, σ y = 0.2976 GPa, The computational region is 0 m ≤ x ≤ 40 m and a transmissive boundary condition is used. We consider the solution derived by PPM + HLLD scheme with a fine mesh as the reference solution for accuracy and error analysis. Here, the grid spacing is 0.1 m and the CFL number is set to be 0.1 so that the dominant error is the spatial error. Figure 1 plots the initial conditions and computed results of this test case at t = 1 ms. In this scale, the differences between the results at 50, 100, and 200 grid points are nearly invisible. The convergence rate of PPM + HLLD scheme are presented in Table 1. From Table 1, we could conclude that the numerical results converge with the order of 3 and thus the PPM + HLLD scheme is third-order accurate for smooth problems.
Test case 2: Elastic-plastic piston-like problem. This case, from Maire et al. 26 , is similar to the one-dimensional hydrodynamic piston test. At the initial time, the copper medium with the length interval of 0 m ≤ x ≤ 1 m is at rest. The velocity at its left boundary is u 1 = 20 m/s and its right boundary satisfies transmissive boundary condition. The material parameters for copper medium are: ρ 0 = 8.93 g/cm 3 , p 01 = 140.7 GPa, p 02 = 287.1 GPa, p 03 = 233.5 GPa, μ = 45 GPa, σ y = 0.09 GPa. The analytical solution of this problem consists of a plastic shock wave and an elastic precursor shock wave. This test case is run until the termination time t = 0.15 ms.
The convergence of the numerical results derived by the PPM + HLLD method with meshes of 100, 200 and 400 cells is presented in Fig. 2. As seen from Fig. 2, the results with different grid resolutions are all converged to the exact solution from Maire et al. 26 Table 2. The numerical accuracy of the whole scheme is first order due to the existence of shock waves in this case.
Test case 3: Elastic-plastic impact problem. This test case, similar to the impact test problem in Titarev et al. 19 , describes the impact situation of two semi-infinite slabs of aluminum traveling with the same velocity and in opposite directions. The computational region is 0 m ≤ x ≤ 10 m and the initial discontinuousness lies at the location of x 0 = 5 m. Further, the parameters for the state equation and constitutive model of aluminum medium are as follows: ρ 0 = 2.7 g/cm 3 , p 01 = 73 GPa, p 02 = 172 GPa, p 03 = 40 GPa, μ = 24.8 GPa, σ y = 0.2976 GPa, τ 0 = 10 −9 s. Figure 3 gives the distributions of density ρ, the velocity in the x direction u 1 as well as the stress components σ 11 and σ 22 at t = 0.6 ms with different impact speeds. The results imply that with lower impact speed, the aluminum material is of elastic-plastic nature and the symmetric impact of aluminum slabs generates symmetric elastic-plastic shock waves. While as stated in Titarev et al. 19 , the aluminum material is elasticity and only the symmetric elastic waves arise in the resulting wave structures. With further increase of impact speed, both the intensity of shock wave and plastic effect are enhanced, and the stress and density behind the shock wave become   larger. As PPM + HLLD uses flattening algorithm to add a little dissipation at the locations of strong shock waves, which is introduced by an additional slope limiter via a 'flattening' parameter and suppresses spurious post-shock oscillations for strong shock waves, and there is no non-physical undershoots and overshoots as in the computations of Titarev et al. 19 , our numerical method is shown by this case to have good effectiveness for solving high speed impact problem.

Test case 4: Wilkins's problem.
This test was originally proposed by Wilkins 10 and deals with the one-dimensional impact of a thin aluminum flyer into a half space of aluminum plate. The initial one-dimensional computational region 0 mm ≤ x ≤ 50 mm consists of two states which are separated by an interface at x 0 = 5 mm.
The left state is moving with a non-zero velocity while the right state is initially at rest. The materials of the flyer and plate considered are both aluminium, and other material parameters are the same as in the test case 3. At the right boundary of aluminum plate, the transmissive boundary condition is satisfied. The left side of flyer is surrounded by vacuum. Further, the interface location is tracked via the level set field and a solid-vacuum Riemann problem is solved to determine the values at the ghost cells at the left boundary.
The simulation starts at the moment that the flyer plate is just in touch with the target. We use 500 cells to simulate the problem up to the termination time t = 5 μs and CFL is equal to 0.6. Since there is no exact solution for comparison, we generate a reference solution by running a high-order method on fine-mesh. Two impact speeds u 1 = 0.8 km/s and 2 km/s are considered. The plots of density ρ and normal stress σ 11 at different times are given with the impact speed u 1 = 0.8 km/s in Figs 4 and 5, respectively. Squares indicate numerical solution while the solid line represents the fine-mesh reference solution. For the lower impact speed, the impact of the aluminum  Table 2. L 1 -norm errors and convergence rates as well as relative CPU times for elastic-plastic piston-like problem. flyer into the plate results in two elastic-plastic shock waves travelling to both the left and right sides. With increasing time, the right-travelling elastic shock wave is distinguishable from the right-travelling plastic shock wave in the right-travelling shock waves due to the larger speed of elastic wave, forming speed and stress steps. The left-travelling shock wave reaches the left free surface and is subsequently reflected to form a right-travelling rarefaction wave. Afterwards, this rarefaction wave is split into elastic rarefaction wave and plastic rarefaction wave obviously due to different wave speeds, which is the so-called elasticity precursor phenomena, and the right-travelling rarefaction wave begins to pursue the initial right-travelling waves. When larger impact speed u 1 = 2 km/s is imposed, the plastic phenomena is more obvious and the peak values of density and stress become even larger, as shown in Figs 6 and 7. As the elastic limit has been exceeded at this time, only the left-travelling and right-travelling plastic shock waves are generated after the impact event. Subsequently the left-travelling plastic shock wave reaches the left free surface and is reflected to form right-travelling rarefaction wave, which will be split into elastic rarefaction wave and plastic rarefaction wave some time later and catch the initial right-travelling plastic shock wave eventually. As no analytic solution exists for this problem leading to incapability of quantitative comparison analysis, but our results of the peak values of both wavespeed and stress are seen to be in good qualitative agreement with wilkins et al. 10 . Further, this case approves both the ability of our numerical technique on processing multi material interfaces and its strong robustness for the plastic-dominated case.
Test case 5: Two-dimensional impact. In this two-dimensional case, we utilize the PPM + HLLD scheme for the multi-material situation to investigate the impact problem of a projectile on a solid plate surrounded by vacuum. The initial configuration of the impact problem is depicted in Fig. 8, in which the projectile is a square with sides of length 0.1 m and the plate is 0.5 m long and 0.1 m wide. The materials of the projectile and plate are assumed to be aluminium. At t = 0, the projectile gets in touch with the plate and all materials are assumed to be in a stress free state: F = I and S = 0. Further, the aluminium target is set to be static, while the aluminium projectile is initialized with The mesh sizes are set as Δx = Δy = 1/800 m, and the CFL number is fixed to be 0.4. Figure 9 gives the density schlieren (top) and the von Mises criteria σ σ ′ − ( ) y 2 2 3 2 (bottom) for times t = 1.0 × 10 −5 s, t = 1.5 × 10 −5 s, t = 2.5 × 10 −5 s and t = 3.0 × 10 −5 s. These schlieren images are all obtained by plotting |▽ρ| field using a logarithmic scale. As seen from Fig. 9(a), the longitudinal shock wave and transverse shear wave are generated in the impactor and plate after the impact. As the speed of the latter, which induces plastic deformation, is lower than that of the former, the wave front of von Mises criteria is slower than that of the longitudinal shock wave in the density plot. With increasing time, both waves spread to the right boundary of plate and are reflected back to the plate ( Fig. 9(b,c)). Accompanied with the propagation of waves in the plate and their multiple reflections, the deformation of the plate is observed clearly and a bulge is formed at the back part. Further, the plastic effects are quite obvious at both the location of contact and the back part of plate, where large deformation exists, whereas other parts bear a relatively small plastic strain. The robustness of the technique PPM + HLLD in two-dimensional cases is well verified based on these qualitative results.

Discussion
This paper has constructed the PPM method suitable for the elastic-plastic deformation problem of solid material under the impact action. This method is based on the elastic-plastic model equation under Eulerian framework, which is formed by the addition of plastic source term into the elastic model equation with deformation gradient tensor as the original variables. In our previous work, the HLL-type solvers are applied to the Riemann problem of nonlinear elasticity and HLLD Riemann solver is shown to be more effective compared with HLL and HLLC solvers for the problems with solution structure of more than five waves. Furthermore, the numerical technique of PPM combined with HLLD solver is seen to have strong accuracy and reliability for the solution of Riemann problem. For the elastic-plastic problem in this study, the solution way is to apply the splitting method to decompose the elastic-plastic problem into the solutions of elasticity control equation and plastic source term. With the purpose of extending the numerical technique of PPM + HLLD to the elastic-plastic problem, the difficulty lies in the fact that the plastic source term may range from zero to positive infinity (i.e. the nature of solid material changes from elasticity to fluid), leading to very small time scale with potentially quite high strain rate. To prevent the possible numerical rigidity, we use the semi implicit method, which is of less computational cost and easier to be realized compared with full implicit method, to solve the ODES for the plastic source term. The high order and effectiveness of our method are confirmed by one-dimensional nontrivial and shockless problem (test case 1) and elastic-plastic piston-like problem (test case 2), respectively. For solving the impact problem with multi solid materials involving different material interfaces and even the interface between solid material and vacuum, we have applied the modified ghost fluid method in fluid mechanics, where the level-set method is used to track material interface and the interaction of the different materials at the two sides of the interface is simulated by the HLLD Riemann solver for multi material cases. The effectiveness of the above method is fully reflected by the computational results of Wilkins's problem and two-dimensional elastic-plastic impact problem.

Methods
Governing equations. We follow the formulation presented in Zhang 25 , where the model equations proposed by Godunov & Romenski 27,28 are used in Eulerian framework for nonlinear elasticity, and extend it to plasticity modelling, so as to accomplish the study of elastic-plastic deformation behavior of solid materials. The hyperbolic partial differential equations that depict mass, momentum and energy conservation in Cartesian coordinates read where ρ denotes the material density, u i denotes velocity, σ 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 strain equation of nonlinear elasticity is described in terms of deformation gradient tensor F ij = ∂x i /∂X j (x i and X j represent the fixed Eulerian coordinates and Lagrangian coordinates respectively): is an artificial vector variable that provides the conservative form of equations for F ij . In order to capture the correct wave speed in the quasi-linear system (12) and to ensure the correctness of numerical results, the term −u i β j is treated as a source term in the numerical computation. To extend elastic strain equation to account for plasticity, we utilize a multiplicative decomposition of the total deformation gradient tensor 13 T e p  where the superscript 'e' denotes elastic and 'p' means plasticity. The evolution of F p , the plastic part of deformation gradient tensor, satisfies the general expression  where L p is called the plastic deformation rate tensor. After substituting (6,7) into (5), the following formula can be obtained: is the constitutive law defining the plastic deformation rate, μ is the shear modulus, σ σ σ δ ′ = − ij ij kk ij is the tensor of deviatoric stress, δ ij is Kronecker delta, and τ is the relaxation time, the introduction of which makes the material deform in a way phenomenologically consistent with Maxwell solid model. Note that when F is mentioned in the system accounting for inelastic deformation (e.g. (8)), it refers to the elastic deformation gradient tensor in fact, rather than total deformation gradient tensor. Further, the relaxation time τ used can be determined from Favrie et al. 16 as 2 is the yield function of von Mises, σ y is the yield limit and τ 0 is a characteristic relaxation time. We make use of the Von Mises criterion for the actual processing of (8). That is, we decide whether the plasticity effect should be introduced in (8) or not according to the comparison results of the equivalent stress derived its prescribed value which can be found in Drumheller et al. 29 and Steinberg 30 . In detail, if the equivalent stress is less than the prescribed value, the plasticity effect is ignored and the elastic model is used directly. On the contrary, the plasticity is introduced via source terms in the equations for elastic deformation gradient tensor F.
It is noted that in conjunction with the equations for the deformation gradient tensor F one has the mass continuity equation (2). Similar to the treatment by Barton 31 , we use the continuity equation (2) rather than the strain conversation equation for the deformation gradient component ρF 11 for maintaining a fully determined system. Then the governing equations can be expressed in matrix form: where U is the vector formed by the conservative variables and F i is the corresponding flux vector. The two terms, S C and S p , on right hand side of (10), are associated with the constraints for the deformation gradient tensor and the plastic deformation, respectively. Here, be the vector of primitive variables. Equation (10) is rewritten as a hyperbolic quasi-linear system where the Jacobian matrix A k is given by and the vector of source terms is The equation of state used in this paper for the specifc internal energy ε is constructed by Miller et al. 13 where ρ 0 is the density of the initially unstressed medium, V = 1/ρ is the specific volume, p(V) is the hydrostatic pressure (Wilkins et al. 10 ) with the expression as and l 1 , l 2 , l 3 are the three independent invariants of the right Cauchy-Green strain tensor C = F T F: Further, the material density, stress tensor, specific internal energy ε and temperature T can be represented as functions of the variables mentioned above: Numerical scheme. To solve the Equation (10), a splitting strategy is used at each time step (Strang et al. 32 ).
Let H Δt be the solution operator corresponding to the homogeneous part and the compatibility part of the problem and S Δt be the operator describing the contribution of the plastic source step over the time step Δt. Starting from the initial data U n , which is the vector of volume averaged conserved variables stored at the cell centres, the solution for the next time level is computed using operator splitting = . The operator H Δt is derived based on the spatially unsplit method is the spatial coordinates, Δx i is the grid spacing in the direction x i and F i±1/2,j,k , F i,j±1/2,k , F i,j,k±1/2 are the numerical fluxes evaluated at the cell boundaries. Each single-component region is given a boundary of "ghost cell" data, determined from a boundary condition, or obtained by a Modified Ghost Method (MGM). The single-component regions are then advanced in time with high-order Godunov methods i.e. piecewise parabolic method (PPM) [33][34][35] . The approach we use for elastic-plastic solid is basically identical the approach described in Zhang 25 . Multi-dimensional non-split method is implemented for PPM. The fluxes in (20) are constructed by approximate one-dimensional Riemann solvers in the direction orthogonal to the cell sides of the Cartesian mesh. For the sake of simplicity we report here the main ideas of the scheme. While to see the details of PPM reconstruction and PPM characteristic tracing, please refer to Miller et al. 33 .
HLLD Riemann Solver. The main idea of HLLD Riemann solver that gives a nonlinear approximate solution is to assume a wave configuration for the solution that consists of five waves separating six constant states. The five waves include two slow waves, two fast waves and a contact discontinuity. Figure 10 shows four intermediate states: ∼ − U , U *− , U *+ , and ∼ + U . We denote the fastest (longitudinal) waves between U ± and ± U as ± S L and the slow shear waves as ± S S that separate the states ± U and U *± . Each wave is considered to be a discontinuity and across each wave the Rankine-Hugoniot relation is satisfied: Since there are more unknowns than equations based on (21a-c), we need to impose some other conditions.
In order to obtain the unknown intermediate state vectors ±  F and ± ⁎ F , we assume the following conditions are 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, while density ρ, normal velocity u 1 and normal stress σ 11 are just the opposite. • Normal stress σ 11 and normal velocity u 1 are continuous across contact discontinuity; shear stress and tangential velocity are equal at the material interface for the 'stick' multi-material problem, while tangential component of the stress vector are zero for the 'slip' multi-material problem.
The left and right fastest wave speeds are computed as L L R L L R 1 1 13 13 and the slow wave speeds are estimated by Within HLL approximation 20 , the intermediate wave speed S C in the present solver is evaluated as where express the i-th components of the vectors ± Q L and ± Q R , respectively. From the wave speeds above, ∼ U state can be described as:  For the 'stick' multi-material problem, the following formulas are satisfied: (3) (4) (4) ; (2) ( 4); (3) ( 4) For the 'slip' multi-material problem, we obtain 2  2  12  3  3  13  12  13 This expression could be used in the particular case of a solid-inviscid fluid or solid-vacuum interaction. Expressions of other variables are described as follows ; ( (7) ( ) )/ (4); ( (8) ( ) )/ (4); ( (9) ( ) )/ (4); ( (10) ( ) )/ (4); ( (11) ( ) )/ (4); ( (12) ( ) )/ (4); ( (13) )/ (4) If we use an explicit Euler scheme, the global time step will be dictated by a CFL condition. However, if this limit way for time step is considered and the strain rates are high enough, the timescales associated with the relaxation operator could be relatively small and hence the system of ODEs Equation (33) will become stiff. Thus, we use the following semi-implicit scheme to prevent the potential of stiff ODEs instead: where (ρF) n+1 is the solution at time t n+1 , the term σ′ is taken explicitly at the previous step t * and the term ρF is chosen implicitly.

Material interface evolution. The Level-Set Equation.
The level-set algorithm is useful for interface tracking in multi-material problem. The central idea of the level-set technique is to represent the moving interface Γ(t) as the level set {ϕ = 0} of a function ϕ. At any moment t, ϕ(x, t) equals to zero on the front Γ(t) and we obtain the evolution equation In general, ϕ(x, t) will not satisfy the signed distance over time. Therefore, we need to adopt the reinitialization algorithm that transforms ϕ(x, t) to make it be the signed distance from location x to interface Γ(t). This transformation is achieved by solving the following equation 0 0 to steady state, which is the desired signed distance function. The 5th-order weighted essentially non-oscillatory (WENO) scheme coupled with the 3rd-order Runge-Kutta time integration scheme is applied to solve Equation (35).
HLLD Multi-Material Riemann Solver. The solution procedure for multi-material problems is similar to HLLD Riemann solver for the single material problem except solid-solid or solid-vacuum boundary conditions at the material interface. The known states W L = W(x i , t n ) and W R = W(x i+1 , t n ) passing an interface of two materials at time t n are used to pose a multi-material Riemann problem. The fluxes on the left and right sides of the interface are given by We define the states in the respective materials' ghost cells as W i+1 = W i+2 = W i+3 = W *− and W i−2 = W i−1 = W i = W *+ . We also adopt an entropy fix technique 36 to suppress 'heating errors' . In practical applications, the initial values of W L and W R are set to be W *− and W *+ , respectively. With these settings, about 10 iterations reaches the convergence precision.
Multi-dimensional implementation of MGFM. For multi-dimensional problems, the implementation of MGFM lies in the fact that the material interfaces are not perpendicular to the coordinate axis, but at a certain angle with them. Thus, it is required to extrapolate the relevant physical quantity g in the direction normal to the material interface by solving the following transport equation: n is the unit vector in the direction normal to material interface. When we perform extrapolation from the region of φ < 0 to that of φ > 0, the '+' sign is used in (38) and the values of physical quantities in the ScIEntIfIc RepoRTs | (2018) 8:9989 | DOI:10.1038/s41598-018-28182-7 region of φ < −ε are kept; while we perform extrapolation from the region of φ > 0 to that of φ < 0, the '−' sign is adopted and the values of physical quantities in the region of φ > ε are maintained. ε is generally set to 3Δx. Then the extrapolation value for material 1 and the real value for material 2 coexist on the ghost meshes at the boundaries of material 1. In order to solve the problem of the angle between the material interface and coordinate axis, one coordinate axis could be rotated to the normal direction of material interface, leading to the transformation of multi dimensional problem to one dimensional problem with the following transformation relationship formula: ROT ROT