An efficient and stabilised SPH method for large strain metal plastic deformations

Due to its simplicity and robustness, smooth particle hydrodynamics (SPH) has been widely used in the modelling of solid and fluid mechanics problems. Through the years, various formulations and stabilisation techniques have been adopted to enhance it. Recently, the authors developed JST–SPH, a mixed formulation based on the SPH method. Originally devised for modelling (nearly) incompressible hyperelasticity, the JST–SPH formulation is mixed in the sense that linear momentum and a number of strain definitions, instead of the displacements, act as main unknowns of the problem. The resulting governing system of conservation laws conveniently enables the application of the Jameson–Schmidt–Turkel (JST) artificial dissipation term, commonly employed in computational fluid dynamics, to solid mechanics. Coupled with meshless SPH discretisation, this novel scheme eliminates the shortcomings encountered when implementing fast dynamics explicit codes using traditional mesh-based methods. This paper focuses on the applicability of the JST–SPH mixed formulation to the simulation of high-rate, large metal elastic–plastic deformations. Three applications—including the simulation of an industry-relevant metal forming process—are examined under different loading conditions, in order to demonstrate the reliability of the method. Results compare favourably with both data from the previous literature, and simulations performed with a commercial finite elements package. Most noticeably, these results demonstrate that the total Lagrangian framework of JST–SPH, fundamental to reduce the computational effort associated with the scheme, retains its accuracy in the presence of large distortions. Moreover, an algorithmic flow chart is included at the end of this document, to facilitate the computer implementation of the scheme.


Introduction
Displacement-based explicit dynamic codes, implemented using low-order finite element methods, are commonly used for advanced numerical simulations in aerospace, automotive, biomedical, defence and manufacturing applications. However, difficulties arise in modelling high-speed impacts or large material deformations, often leading to simulation failure. In these scenarios, most computational codes prefer to employ the 8-noded underintegrated hexahedral element to model solid components [5]. Nonetheless, for many practical applications (e. g. crashworthiness, ballistics and metal forming), the extremely large deformations still result in severe B Giorgio Greto gretog@cardiff.ac.uk 1 School of Engineering, Cardiff University, Queen's Buildings, The Parade, Cardiff CF24 3AA, UK mesh distortions. These lead to irregular and entangled elements, unless some form of adaptive remeshing is employed.
In addition, difficulties currently encountered by the widely employed, displacement-based, low-order finite element analysis include lower order of convergence for relevant field variables; excessive element distortion under large deformations, requiring periodic remeshing; locking behaviour in bending scenarios; non-physical pressure instabilities, and high-frequency noise due to Newmark-type time integrators.
Recent developments in computational methods for fast solid dynamics [1,2,[7][8][9]18,19,21,[23][24][25] recommend the representation of motion and deformation of a given body via a system of first-order, mixed formulation conservation laws. The partial differential equations (PDEs) that form this system do not present the displacement as the main unknown to be evaluated, instead yielding a set of other relevant quantities (i. e. velocity, deformation gradient, its cofactor matrix and scalar Jacobian), depending on the number of conser-vation equations being considered. This choice of the field variables will be driven on the one hand, by the need to reproduce certain features of material behaviour (e. g. near or full incompressibility) while ensuring polyconvexity of the strain energy function and, on the other hand, by the computational efficiency of the resulting scheme.
Mixed conservation laws have already been proven [1,21,25] to achieve second order of accuracy for stresses and strains. The two-field mixed formulation, composed of linear momentum p and deformation gradient F, was later augmented by including a governing equation to conserve the Jacobian J of the deformation gradient [18]. This enabled the scheme to effectively solve nearly and fully incompressible materials. Further enhancement of this { p, F, J } framework has also been recently reported in [7,19] for modelling compressible, nearly incompressible and fully incompressible materials governed by a polyconvex constitutive law. In particular, this requires the presence of the matrix of co-factors of the deformation gradient F (denoted here by H), leading to an extended system of field variables. The { p, F, H, J } system can be reformulated to an alternative description in terms of entropy conjugates, as in [7,19], due to the existence of a generalised convex entropy function. This readily facilitates the proof of existence of associated real wave speeds, which in turn establishes existence and uniqueness of real solutions.
A novel computational framework was devised to improve the numerical simulation behaviour during large material distortions, and to address the other numerical issues (e. g. locking, spurious oscillations) highlighted earlier. In a total Lagrangian perspective, the proposed methodology combines the use of smooth particle hydrodynamics (SPH), a meshless spatial discretisation technique, with an explicit, two-stage, total variation diminishing (TVD) Runge-Kutta temporal scheme.
To date, the application of mixed formulation has been largely focused on nearly incompressible hyperelasticity [1,2,21,25], where many of the aforementioned shortcomings are encountered while using linear finite element methods. The work presented in this paper aims to investigate the feasibility of using the developed mixed formulation technique in the field of metal plasticity, by exploring some applications involving large material deformations. These will include a high-velocity impact scenario, where material is subjected to a non-uniform deformation and thus develops a large geometric discontinuity; a constrained boundary problem under heavy distortion, and the simulation of severe plastic deformation in a metallurgical process of industrial relevance. This paper is organised as follows: a description of the mixed hyperbolic system of governing partial differential equations is presented in Sect. 2, while the adopted SPH discretisation methodology, that incorporates a stabilisation procedure based on a Jameson-Schmidt-Turkel (JST) artificial dissipation term, widely used in the field of com-putational fluid dynamics (CFD), will be detailed in Sect. 3. Numerical applications in metal plasticity are presented in Sect. 4, along with an assessment of the accuracy and stability of the methodology. These numerical examples are mainly chosen to demonstrate the potential and the accuracy of the developed computational strategy. Concluding remarks are then presented in Sect. 5. Finally, for the purpose of completeness, a brief outline of the chosen elasto-plastic constitutive model, including relevant computational procedures for evaluating the material deformation, is provided in "Appendix A".

Conservation laws
In this paper, a vector will be expressed with the notation a, and a tensor with A. The material and spatial coordinates are denoted by X and x, respectively. The gradient operator that refers to initial spatial coordinates will be denoted by ∇ 0 . The symbol × × × represents the tensorial cross product, a key operation utilised in obtaining a concise mathematical representation of the proposed methodology. This notation was introduced in a series of papers devoted to the development of mixed formulation techniques [7][8][9]19]. More details on the properties of the tensorial cross product can be found in the cited references.
The system of conservation laws describing the motion and deformation processes in terms of field variables p, F, H, J can be written as: In Eq. (2.1a), ρ and b, respectively, represent the material density and the body force per unit mass. System (2.1) is solved for p, F, H, J with respect to time t. It is a set of hyperbolic PDEs similar in structure to those widely used in CFD [42]. This similarity enables one to employ wellestablished artificial dissipation numerical techniques from CFD to improve the computational stability of Eqs. (2.1a) and (2.1d) [1,23,24]. It is worth noting that, in the presence of non-smooth solutions, the above system (2.1) of local conservation laws must be accompanied by suitable jump conditions, as described in [42].
The notation adopted, along with a pictorial representation of the main measures of deformation used in (2.1), is illustrated in Fig. 1.  Fig. 1 Finite deformation process on a body B in the continuum. Shown are measures of fibre, area, and volume strain-F, H, J -that define the process from the initial configuration B 0 , to B(t) at time t It is evident from the expressions (2.1b) and (2.1c), respectively, for the F and H conservation variables, that, for both physical (required compatibility of strains, [8]) and mathematical (ensure convexity, [17]) reasons, two sets of involutions need to be satisfied by the solution of (2.1): In the above equations, the components of the material curl of F in (2.2a), and of the material divergence of H in (2.2b), can be expanded with the help of the Levi-Civita permutation symbol E i jk 1 as, respectively, Using the involutions described in (2.2), Eqs. (2.1c) and (2.1d) can be further simplified into [25]: In the place of the full set of conservation variables, denoted here by U = { p, F, H, J }, reduced systems based on { p, F, J }, or only { p, F} formulations have been adopted in the past. The robustness of these reduced systems has been positively ascertained by testing them against a thorough variety of numerical regimes and external conditions [1,18,24]. However, the complete { p, F, H, J } mixed system of conservation laws (2.1) can be provided with mathematicalinstead of merely empirical-proof of existence and stability of physically relevant solutions. Such a proof may be obtained by verifying the hyperbolicity of the system, i. e. whether U can be expressed in wave-like form [18,27]: In (2.4), Z is a chosen direction, λ i is the wave velocity and an eigenvalue of the system matrix, and R i is the wave profile and the eigenvector of the system matrix corresponding to λ i . As can be seen from (2.4), in order to prove the hyperbolicity of the system, its own matrix should present real and distinct eigenvalues, and eigenvectors orthogonal to each other. Only then the hyperbolicity of a system of conservation laws can be guaranteed by specifically choosing an elastic potential function Ψ , for description of the material in use, able to satisfy certain convexity conditions. This would ensure the existence of minimum values of Ψ , for the solution of the variational problem associated with an elastic system [27].
The most stringent convexity condition is the notion of polyconvexity [3,27]. More precisely, the polyconvexity condition states that the potential strain function Ψ (F, H, J ) has to be convex in the function space formed by the components of F and H, and by J [3,8]. This also ensures a one-to-one mapping between the stress P and strain measures {F, H, J } [3,16,27,31].
The JST-SPH methodology consists of three key features: -spatial discretisation with the SPH scheme; -a JST-based numerical dissipation tool, adopted from CFD, to improve the accuracy and the stabilisation of the overall discretisation procedure; -an explicit, two-stage total variation diminishing Runge-Kutta time integrator scheme to follow the time evolution of the solutions during dynamic simulations.
The above main features of the proposed methodology will be briefly elaborated in the following sections.

The SPH method
Being a meshfree technique, SPH can be effectively employed in the simulation of high-velocity impacts and high strain-rate deformations [10,15,20,26,28,29,34,39,40]. In fact, beyond the yield point, SPH discretisation becomes attractive because large distortions make standard Lagrangian finite elements analyses difficult to pursue, without resorting to remeshing techniques that increase the computational cost [5,38]. This subsection briefly outlines the background for the SPH method and details the total Lagrangian SPH discretisation of the governing system of equations. The adoption of a total Lagrangian SPH formulation allows to avoid numerical tensile instabilities [4,11,33,40]. However, it was proven in [32] that for simulations involving material discontinuities (such as cracks and fragmentation), stabler results are achieved when switching to Eulerian kernels over the discontinuity regions while retaining Lagrangian kernels elsewhere.
In a discretised domain, the SPH interpolation of a function variable φ( x) at any arbitrary point x, denoted here by φ( x) , can be approximated by: where V b is the fractional volume assigned to a neighbouring particle in x b , h is the smoothing length, and W ( Although studies introducing adaptivity via a variable radius of support h exist in the particle methods literature (see [35]), here h and consequently D and V for in-field particles are kept constant for simplicity. For accurate interpolation of any given function, the polynomial kernel should satisfy the following reproducibility conditions: and δ( r ) represents the Dirac delta function centred at 0. The first derivative of a given function φ can be discretised in terms of SPH interpolation as Several different types of kernel functions have been employed as SPH interpolators. In the present work, to maintain the continuity for spatial derivatives of kernel function up to the fourth order, the following quintic spline polynomial is adopted as the kernel function: In (3.4), d is the number of spatial dimensions of the problem, and α is a normalising constant that depends on d as, It is widely known that the fundamental discretised form of SPH interpolation in (3.1) suffers from poor accuracies at and in the vicinity of the domain boundaries. More precisely, (3.1) does not completely satisfy the partition of unity condition at or near the boundaries, due to the truncation of kernel functions. To improve the accuracy of the SPH interpolation near boundaries, and to exactly preserve momentum, corrections must be introduced on both the kernel and the gradient of the kernel [10][11][12]. In the present work, the accuracies of the kernel and the gradient of kernel are improved by using correction methods proposed in [12]. The corrected SPH ker-nelW ( r , h) and the corrected gradient of kernel ∇ 0 W ( r , h) defined in [12] exactly fulfil first-order completeness, that is, they allow the SPH representation to exactly reproduce linear fields.
The corrected SPH interpolation of a given function f ( X ) will thus be expressed as, (3.6) In (3.6), Λ b a identifies the set of neighbours b of a particle a, at which the gradient is evaluated. The mixed { p, F, H, J } system described by Eqs. (2.1) and (2.3) can now be discretised in space using the corrected SPH interpolation as described above.
Further, to obtain the discretised form of (2.1a), the weak statement for the linear momentum evolution [10][11][12]43] can be obtained by employing work-conjugate pairs [7,19]. With the help of integration by parts, this yields In (3.7), virtual velocities are denoted by δ v. By using particle discretisation on the domain V , the above expression (3.7) is approximated to: In (3.8), A a and t a are, respectively, associated area and traction force assigned to particles a located at or near the domain boundary. The density is assumed to be homogeneous throughout the domain. Equation (3.8) can now be expressed in terms of SPH discretisation as The last term on the right-hand side of Eq. (3.9) represents the resultant internal force T a acting on particle a: Conservation laws given by Eqs. (2.3) and (2.1d), respectively, for strain measures F, H, and J , can also be discretised in a similar manner to what has been done for (3.9), with the linear momentum p evaluated through SPH: Depending on the chosen materials, the above semidiscrete formulation (3.9)-(3.11) may still suffer from accumulated numerical instabilities over a long span of time response. Therefore, it is necessary to address any such instability via suitable computational procedures.

JST artificial dissipation
In the context of SPH, various artificial dissipation techniques have been used in the past [28][29][30] to alleviate the aforementioned numerical instabilities. Alternatively, SPH can be stabilised by introducing stress points in the formulation (see [33]). In this case, however, a background mesh would be needed for the computation and assignment of fractional volumes. In the present work, an adapted nodally conservative JST stabilisation term D JST (U a ) will be incorporated into (3.9), mirroring CFD techniques. The hyperbolic nature of the first-order conservation laws in system (2.1), reflecting that of the Euler equations in fluid dynamics, makes it possible to introduce the JST term as a dissipative component [1].
The nodally conservative JST stabilisation is additively decomposed into a second-order (harmonic) operator D 2 (U a ) and a fourth-order (biharmonic) operator D 4 (U a ) as: In (3.13), c p is the pressure wave speed, x min is the particle spacing, κ (2) and κ (4) are user-defined parameters, while the∇ 2 0 symbol represents a corrected Laplacian operator, applied to kernel W as detailed in [10].

Semi-discrete governing equations
Out of the four unknowns in the vector of state U , only p and J will see JST dissipation terms appear in their respective conservation laws (3.9) and (3.11c). D JST (F) = D JST (H) = 0 are imposed in order to respect conditions (2.2), set to ensure these strain measures preserve compatibility with the domain motions. This is because F and H are now main independent variables of the problem, and are not associated with displacements x. In the context of the standard displacementbased method, it is worth recalling that F is directly computed from x, which acts as main independent variable of the problem.
In the light of the above, the spatial semi-discretisation yielded by the JST-SPH methodology reads as: (3.14) Due to the incompatibility between strains and displacements introduced by the adoption of the mixed { p, F, H, J } formulation, discussed above, the conservation of angular momentum will not in general be satisfied by the discretised system of Eq. (3.14). In addition, the conservation of linear momentum may also not be respected, due to the lack of radial symmetry introduced by kernel and gradient corrections that are applied on particles located at, or in the vicinity, of domain boundaries.
With the aim to preserve the linear and angular momenta, the internal forces T and the JST dissipation terms D JST have to be amended, in order to enforce the following conditions: To impose conditions (3.15), a procedure based on the method of Lagrangian multipliers has been implemented (see [25] for more details).

Two-stage TVD Runge-Kutta temporal integration
To obtain a fully discretised system, that is able to describe the time evolution of the solution, the set of particle equations (3.14) has to be now explicitly integrated between chosen time intervals. The solution will update from current time step t n to next step t n+1 , until the end of the simulation. An explicit, twostage TVD Runge-Kutta integrator scheme will be used for the purpose of time integration as follows: The size of the time step, t, in (3.16) is obtained using the Courant-Friedrichs-Lewy (CFL) condition [42]: In (3.17), h l is the characteristic length of the problem, here chosen to be equal to the initial distance between two particles, and α c f l is a user-defined constant assuming values

Applications
In this section, the proposed JST-SPH algorithm will be used to perform three numerical examples. The first amongst these will be the classic Taylor bar test for high-speed-impact plasticity effects: a three-dimensional cylindrical metal bar that impacts a horizontal wall at a high velocity. In the second example, the performance of the method is assessed by examining the elasto-plastic compression of a constrained rectangular body. The last example involves the simulation of the equal channel angular extrusion (ECAE) cold metal forming process, where a workpiece is forced to pass through two 90°turns inside an extrusion die. The localised distortions caused by the substantial shear stresses experienced by the workpiece at these 90°turns should constitute an ideal test for assessing the capability of the developed methodology under large strain, high-velocity conditions.

Taylor bar high-speed impact
A common benchmark test for plastic deformation at high speed is the Taylor bar impact problem [41], illustrated in Fig. 2. Elastic behaviour is governed by the material energy function presented in Eq. (A.1). Equation (A.1) is coupled with the plastic yield model described in "Appendix A", using the The bar has initial height h 0 = 32 mm and radius r 0 = 3.2 mm, and is discretised by a set of 4131 particles arranged in a regular pattern. The impact is supposed frictionless and takes place at an initial speed of the bar of v 0z = −227 m/s in the z direction, normal to the rigid wall. Due to the extensive plastic dissipation, the JST stabilising term is set to a very low value. Figure 3 presents results of the Taylor bar simulation performed by JST-SPH. As observed in experimental results [44], in Fig. 3 the plastic front appears closer to the bottom wall at the early stages of the simulation (red contour regions). It then slowly moves towards the top of the bar, as more kinetic energy dissipates into plastic strain.
To investigate the stress distribution across the bar, the von Mises equivalent stress is plotted in Fig. 4. A gradual increase in the plastic stress flow over time can be clearly observed.
In Table 1, the radius of the bottom face at time t = 80 µs is compared to the results of identical tests performed using different numerical methodologies [1]. It can be noted that the mixed formulation techniques, being locking-free, avoid excessive rigid response of the structure often demonstrated by low-order finite elements, displacement-based analyses.
Plots of the total internal energy and of energy dissipation associated with plastic deformation and JST artificial dissipation are reported in Fig. 5.
The equations used to evaluate the total internal energy U int , the total plastic dissipation W ( p) and the artificial JST dissipation W (JST) at each time step are given as follows: It is shown in Fig. 5 that artificial dissipation introduced by the JST algorithm constitutes a negligible quantity.

Highly constrained problem with moving discontinuity
An example of the use of JST-SPH in a highly constrained setting is presented in [21,25] for a fully elastic, nearly incompressible material. The effects of plasticity will be studied in this similar example, where a bidimensional rectangular block with dimensions 1 m × 0.5 m is considered under the assumption of plane strain. The material is copper, with properties chosen to be the same as those in the Taylor bar example in Sect. 4.1. An external velocity v = (0, −10) T m/s is applied throughout the simulation on the top side of the block, as These specific values of κ (2) and κ (4) were chosen in order to provide the simulation with a small amount of JST numerical dissipation, useful to suppress any spurious pressure oscillations at the initial stages of the process.
The purpose of this example is to demonstrate the capability of the total Lagrangian framework of the JST-SPH method, in the presence of high, localised distortions. The evolution of the von Mises equivalent plastic strain ( p) during the analysis is shown in the top row of Fig. 7. For comparison, the same simulation has been run using a commercial finite element method package, without any remeshing, and the results are presented in the bottom row of Fig. 7.
Results obtained by using increasing numbers of particles in the simulation are presented in Fig. 8.
It is observed from Figs. 7 and 8 that the basic finite element method is not able to handle the excessive distortions, whereas the JST-SPH method is able to simulate severe material deformation without exhibiting any numerical instability.
Plots of the total internal energy and of energy dissipation associated with plastic deformation and JST artificial dissipation are reported in Fig. 9. Energy quantities have been calculated using (4.1). It is shown in Fig. 9 that artificial dissipation introduced by the JST algorithm is a negligible quantity. Roller constraints imposed at the boundaries of the block are at the origin of irregularities appearing in the internal energy plot [calculated as in (4.1a)].

Simulation of the equal channel angular extrusion (ECAE) process
Numerous metallurgical techniques exploit the mechanics of plastic deformation to attain smaller grain size, and consequently better material properties, for metals and alloys. Amongst them, ECAE [36,37] can offer high levels of shear  strain for relatively low levels of external pressure, making it an ideal material processing technique for mass production.
In the present setting, the ECAE process will be simulated in a channel with two 90°turns, analogous to the set-up analysed in [36] using a commercial finite element software. This example will demonstrate the robustness of the JST-SPH numerical scheme under a demanding dynamical process, where the material is subjected to continuous large deformation.

Fig. 10 ECAE simulation, initial set-up
The analysis is performed in plane strain conditions, and the material is made of commercially pure aluminium (Al1100, E = 69 GPa, ν = 0.33, ρ = 2800 kg/m 3 , width l = 8 mm). The walls of the die are assumed rigid. Its rightangled corners are rounded, with 1.5 mm and 1 mm external and internal fillet radii. The initial set-up of the experiment is presented in Fig. 10.
The simulations are performed with a regular mesh of 400 particles. Contact between workpiece and die is assumed to be lubricated, and considered frictionless during the simulation. For the purpose of simplicity, a contact algorithm based on a reflective, bouncing back procedure [22] is adopted to Given that large plastic deformations are expected, it is useful to introduce numerical dissipation at the early stages of the simulation, to prevent pressure instabilities that may arise when the billet initially comes into contact with the bottom wall of the channel. To facilitate this, for each particle, the JST term k (4) was defined as function of the equivalent plastic strain ( p) , to ensure that k (4) linearly decreases, as plastic dissipation gradually sets in. This is achieved by adopting the following formula: The plastic yield stress of the billet material is assumed to follow the nonlinear hardening law given by: .3 is solved for ( p) numerically by using the Newton-Raphson method. Figures 11 and 12, respectively, depict the distributions of pressure and plastic strains during the deformation process of the billet, at various instants in time. For the purpose of comparison, Figs. 11 and 12 also present the corresponding pressure and plastic strain distributions obtained with a commercial finite element software (Abaqus/Explicit) using linear quadrilateral elements.
It can be noted that the results from the JST-SPH method correlate well with those obtained with the finite element method. Results from Figs. 11 and 12 are also consistent with observations reported in [36]. It is evident from the figures that the particle distribution in the vicinity of sharp corners of the channel is not smooth at the later stages of the simulation. This uneven distribution of the particles can be attributed to boundary effects. The plastic strain contour illustrated in Fig. 12 compares well with results reported in [36]. Further, it is clearly evident that the plastic strain distribution increases significantly as the billet material passes through each bend of the channel in the die. This feature is what is mainly exploited in the ECAE process to achieve the desired level of plastic strain in the workpiece.
Investigation of the energy patterns throughout the simulation confirms that plastic deformation is the main mechanism of energy dissipation, being on average two orders of magnitude larger than the energy dissipation introduced by activating the JST stabilisation term. Figure 13 compares the energies associated with plastic deformation, JST artificial dissipation, and the total internal energy of the bar during the simulation of the entire ECAE process. Energy quantities in Fig. 13 were calculated using (4.1).
As shown in Fig. 13, the values of total internal energy and of total plastic dissipation follow specific trends that can be linked to the motion of the billet inside the channel of the die. For example, a peak value for both energies was observed around 0.00142 s, corresponding to the instant in which the billet comes into contact with the bottom wall of the middle (horizontal) part of the channel. Following that, the energy values begin to decline steeply during the expansion of the billet material in the mid-channel. A similar, but less pronounced energy pattern can be observed once the bar reaches, and then begins to flow into the second channel, with a peak value attained around 0.035 s. Moreover, it transpires from Fig. 13 that the energy values are highly oscillatory when the billet initially passes through the angular sections of the channel. However, the frequency and the intensity of these oscillations gradually decrease, as the material moves away from the angular sections.
In addition, Fig. 13 shows the gradual decrease in energy dissipation associated with the JST term at the onset of plastic deformation. This effect follows from the definition of the JST terms in (4.2), and further demonstrates that higher artificial dissipation values are only required at the very early stages of the simulation, that is, before any plastic deformation takes place.
In order to verify the accuracy of the simulation, a convergence study was performed by varying the number of particles in the discretised domain.
The quantity under analysis is the amount of plastic strain at 0.035 s, both in local (the maximum equivalent plastic strain max ( p) amongst particles) and in global (the plastic dissipation developed in the single time step leading to instant 0.035 s) terms. Apart from the mesh of 400 particles used so far, alternative meshes with 88, 216 and 640 particles were employed. The results of the convergence analysis are presented in Fig. 14. It can be noted from the figure that the chosen parameters monotonically converge with increasing particle density. The mesh with 640 particles has also been utilised for closer comparison of the degree of similarity between the JST-SPH ECAE test conducted here with the analysis found in [36]. For this purpose, the strain contour in the direction normal to the billet movement has been investigated towards the end of the test, at simulation time t = 0.04 s. The two locations selected are "Section 1" positioned at the centre of the middle channel, and "Section 2" at 5 mm from the second corner of the die. The above locations are consistent with those of the results reported in [36].
Ten equally spaced target positions have been chosen across the two sections of interest. Local values of the plastic strain ( p) have been computed using SPH interpolated values over neighbouring particles. Results are presented in Fig. 15. Comparison of Fig. 15 with results reported in [36] reveals that plastic strain patterns of the two tests match closely along Section 2, the main difference being that the JST-SPH analysis predictably presents a smoother curve. The same can be observed for Section 1, for which there is qualitative agreement with results from [36], except for the points near the upper wall, where results obtained with JST-SPH yield larger plastic strain. This discrepancy can be attributed to the fact that these points lie in the vicinity of the upper wall, and closer to the first bend. As can be noted in Fig. 12, also the accuracy of the plastic strain at these points is affected by the same boundary effects. in [36] Parametric analyses reported in [36] provided an opportunity for further verification of the validity of the results.
In [36], geometrical properties of the die are modified in order to investigate the possibility of enhancing ECAE metallurgical performances. One of these parameters was the length of the middle channel of the die. In [36], it was noted that interesting deformed shapes develop, when the middle channel is shortened from 24 to 16 mm. The upper part of second bend has a large unfilled area, while the billet experiences more stress due to bending than due to shear. This effect is particularly severe around the inner region of the billet in the vicinity of the bottom part of the second bend. Both of these phenomena could be identified in an analogous simulation performed with JST-SPH, as presented in Fig. 16 at time t = 0.027 s. Results obtained with the finite elements commercial solver Abaqus/Explicit, using linear quadrilateral elements, are also shown in Fig. 16.
It is evident from the numerical examples that the JST-SPH scheme eliminates the excessive element distortions usually associated with mesh-based methods during the simulation of large plastic strains. In addition, the total Lagrangian framework adopted here provides an efficient modelling tool to gain better understanding of the physics underlying large material deformations.
An algorithmic diagram, illustrating the JST-SPH solution procedure for simulating elasto-plastic material deformations, is presented next.

JST-SPH algorithm for elasto-plastic simulations
Performs Total Lagrangian dynamics simulations of an elasto-plastic process using { p, F, H, J } mixed formulation: SPH space discretisation + JST dissipation term + 2 stages TVD Runge-Kutta explicit time stepping.

Initialisation
-Initiate all variables (set of particles geometry, material properties, numerical parameters, p i , F i , P i , x i for each particle i). -Compute SPH kernel smoothing length h = α(V /N ) 1/3 , with N number of particles and α a user-defined input parameter. -Locate neighbours j for each target particle i, using the alternating digital tree (ADT) search algorithm [13]. -Compute corrected kernelW , its gradient ∇ W and its Laplacian∇ 2 W . FOR RK = 1, 2 -Advance p (3.14).
-Advance the other variables F, H and J (3.14).
-Impose boundary conditions • Identify particles that have crossed the walls: IF | x i | > x W THEN · Determine normal to the wall. · Reflect p i around the normal to the wall. · Place particle position on the wall, x i = x W . -Compute first Piola-Kirchhoff stress tensor P i as in "Appendix A": • Given initial C ( p) , J and ( p) .
• Calculate pressure p.

Conclusions
This paper has focused on the application of JST-SPH mixed formulation in metal plasticity. This methodology has proved to be a valid tool for computing plastic deformations under dynamic regimes. It has been shown that there is good agreement between results reported in the previous section and data for the same tests found in the literature. The Taylor bar impact case provides an ideal application for JST-SPH, since SPH performs particularly well in high-velocity, large deformation settings. On the other hand, the ECAE test confirmed the robustness of JST-SPH: the total Lagrangian framework was able to withstand high local gradients of stresses, strains and displacements, while the presence of a constant loading in time did not lead to any kind of instability. This robustness was further evidenced in the constrained block simulation, where it can be seen that JST-SPH is able to achieve accurate results while preserving the initial connectivity under very large distortions. The same outcome cannot be obtained with standard mesh-based methods, without remeshing. It was noted, in the case of the ECAE simulation, that the results in the vicinity of the boundaries were slightly affected by sharp corners. These boundary effects can be easily eliminated by using more sophisticated boundary treatments. This will be one of the focal points in future publications. constitutive relations In the following, a stretch-based, elastic energy function will be defined as [14]: where μ and λ are the Lamé constants, and λ 1 , λ 2 , λ 3 are the principal stretches, such that the right Cauchy-Green deformation tensor, C [14], can be expressed as: In (A.2), N a , a = 1, 2, 3, are the principal directions in the reference frame.
From (A.1), it is possible to calculate the first Piola-Kirchhoff stress tensor P for a polyconvex material in the elastic realm and, in so doing, to close the system of conservation laws (2.1) [7,8,19]: In (A.3), F , H and J are the entropy conjugates of strain variables F, H and J . More detailed definitions of entropy conjugates can be found in [7,19]. Equation (A.3) provides solutions while in the elastic regime. However, once the deformation state progresses beyond the yield point and acquires a plastic component, the first Piola-Kirchhoff stress tensor P will be hereby computed as a function of the Cauchy stress σ in the current configuration, following [14]: The stress tensor σ in (A.4) can be expressed in spatial principal directions n a , a = 1, 2, 3 as σ = 3 a=1 σ a n a ⊗ n a (A.5) where σ a is given by 1  3 ln J a = 1, 2, 3 (A. 6) In (A.6), the bulk modulus κ for the material was introduced. Operating in principal directions holds the distinct advantage of making quantities computed in the current section automatically frame invariant.
The elastic principal stretches λ (e)a for a = 1, 2, 3 in (A.6), and the principal directions n a , a = 1, 2, 3 in (A.5) are obtained as the eigenvalues and eigenvectors of the elastic part b (e) of the left Cauchy-Green deformation tensor b, as b (e) = 3 a=1 λ 2 (e)a n a ⊗ n a (A.7) The plastic part of the deformation enters the algorithm just described through the preliminary evaluation of b (e) , obtained using the plastic part of the right Cauchy-Green strain tensor in (A. In (A.10), σ Y = σ Y ( ( p) ) is the yield stress, that in general will be a function of the history of accumulated plastic deformation during the process, through a material hardening rule.
In case condition (A.10) does not hold, plastic deformation develops, ( p) > 0, and the yield surface is updated through a radial return mapping procedure, whereby the normal to the yield surface m = (m 1 , m 2 , m 3 ) is expressed by m a = σ dev,a √ σ dev : σ dev a = 1, 2, 3 (A. 11) and the principal components of the deviatoric stress tensor σ dev,a are updated as [14] σ upd dev,a = σ dev,a − 2μ ( p) m a a = 1, 2, 3 ( A . 1 2 ) The resulting increment in equivalent plastic strain ( p) will be computed as the strain amount necessary to keep the updated equivalent plastic stress σ In the case of the classic hardening rule that linearly depends on ( p) , i. e.
where σ Y 0 is the initial yield stress value, and H is the hardening modulus, ( p) can be explicitly obtained as If instead a more complicated dependency of the type σ Y = σ Y ( ( p) ) is in place, then ( p) has to be extracted numerically.
Once ( p) is known, the elastic principal stretches λ (e)a can be corrected to take the plastic dissipation into account  What is left in order to complete the algorithm is then to substitute in Eq. (A.5) σ a with σ upd dev,a obtained from (A.17), and to add the dilatation pressure term κ ln J J . The plastic deformation state is also updated as The von Mises plastic algorithm described in this section is depicted in Fig. 17 in the case of isotropic hardening.