A New Updated Reference Lagrangian Smooth Particle Hydrodynamics algorithm for isothermal elasticity and elasto-plasticity

This paper presents a new Updated Reference Lagrangian Smooth Particle Hydrodynamics (SPH) algorithm for the analysis of large deformation isothermal elasticity and elasto-plasticity. Taking as point of departure a Total Lagrangian setting and considering as referential configuration an intermediate configuration of the deformation process, the equation of conservation of linear momentum and three geometric conservation laws (for the deformation gradient, its cofactor and its determinant) are re-written leading to a very generic system of first order conservation laws. The key feature of the formulation is a suitable multiplicative decomposition of the conservation variables, leading to a very simple final set of equations with striking similarities to the conventional Total Lagrangian system albeit re-written in terms of alternative Referential Updated conservation variables which are evolved in time. Taking advantage of this new Updated Reference Lagrangian formalism, a second order entropy-stable SPH upwinding stabilisation method will be introduced. With respect to previous publications by the group, a new three-stage Runge–Kutta time integration method is implemented in order to increase the CFL stability restriction. Finally, and to demonstrate the robustness and applicability of the methodology, a wide spectrum of challenging problems will be presented and compared, including some benchmark three-dimensional large deformation elasto-plasticity problems. To facilitate its ease of implementation, the paper explores the use of a series of novel expressions for the evaluation of kernels and the gradients of kernels to provide the SPH user the freedom to choose amongst various options, all leading to equally convincingly robust results.


Introduction
Numerical methods that do not require connection between nodes or the creation of sub-regions (elements) in the process of discretisation comprise a large family of methods called meshless (or meshfree) methods, see e.g.Belytschko et al. [1], Huerta et al. [2] and Chen et al. [3] for in-depth studies.The main objective of such techniques is to eliminate some of the limitations related to mesh dependence of conventional computational methods.Since the approximations are constructed in terms of nodes (particles) and not based on elements (their distortion may degrade accuracy), meshfree methods can address a range of problems much wider than mesh-based methods.Moreover, meshfree methods [4] have the potential of circumventing some difficulties usually faced in the pre-processing stage, i.e., meshfree methods eliminate the standard mesh generation process, making the integration between Computer Aided Engineering (CAE) and Computer Aided Design (CAD) easier [5].
One of the earliest meshfree methods is the Smooth Particle Hydrodynamics (SPH), first presented in 1977 by Gingold et al. [6] and Lucy [7].Libersky et al. were amongst the first to use the SPH method in the context of strength of materials [8] and dynamic solid mechanics [9].The SPH method gives numerical solutions to initialboundary value problems defined by conservation laws of continuum mechanics combined with constitutive relations for the materials involved.It is important to emphasise that SPH is not based on discrete particles colliding with each other or presenting cohesive-like behaviour.Instead, SPH is a fully Lagrangian modelling scheme that permits the discretisation of continuum partial differential equations through the approximation of properties at a set of points distributed over the domain, without the need of defining a spatial mesh.The Lagrangian nature of the method and its ability to keep accurate history of events, associated with its meshless characteristics, represent the main advantages of SPH.A massive number of publications regarding the SPH method is available in the literature.The works by Monaghan [10], Swegle et al. [11] and Vignjevic et al. [12] and references therein provide comprehensive reviews of the topic.
The classical Eulerian and Updated Lagrangian displacement-based SPH methods are known to suffer from a number of well-known drawbacks, such as numerical instabilities in the form of spurious zero-energy modes and tensile instability (the latter only in the case of Eulerian SPH) [13][14][15][16][17][18][19][20][21][22][23][24][25][26], lack of consistency [18,20,27,28], loss of conservation [24,29,30] and reduced order of convergence for derived variables such as stresses and strains [24].Distinct research groups [31][32][33][34][35][36][37][38][39][40][41][42][43][44] have come up with various solutions for the different issues faced in SPH.Some interesting work has also been reported in [45,46] where a hourglass control based stabilisation algorithm is employed for the description of elasto-and visco-plastic continuum.It is however still not yet clear how to introduce appropriate numerical viscosity through the classical Coleman-Noll procedure in order to ensure the production of non-negative numerical entropy.Additionally, the kernel function as well as the gradient of its kernel used in [45,46] are evaluated based on standard updated isotropic (spherical) kernels, which may not be optimal in the presence of strong anisotropic changes in volume.In the recent work by Lee et al. [26], a Riemann solver based stabilisation strategy for a mixed-based SPH framework for large strain explicit solid dynamics was presented.The framework distinguishes from standard SPH techniques due to the fact that the conservation of the linear momentum p is solved along with the conservation equations for the deformation gradient F and its minors.Conservation laws of geometric quantities such as the deformation gradient F, the volume map (Jacobian) J and the area map (cofactor) H are introduced aiming at the establishment of a mixed set of conservation laws.In previous works, but in different contexts, Lee et al. [47,48] proposed other methodologies which had as unknowns the linear momentum, the deformation gradient and the total energy.With both the linear momentum and the deformation gradient tensor being primary variables of the problem, stresses converge at the same rate as the velocities and displacements.Moreover, the new formulation was shown to be efficient in nearly incompressible and bending dominated scenarios.However, in case of extreme deformations in the incompressible limit, the { p, F} formulation lacks robustness.In [49], in the context of Finite Element Method (FEM), Gil et al. enhanced the formulation for nearly and truly incompressible deformations with the novelty of introducing a conservation law for the Jacobian J of the deformation gradient, providing extra flexibility to the scheme.In [50], a new geometric conservation law for the co-factor H of the deformation gradient was also added to the framework, leading to an enhanced mixed-formulation, very amenable for constitutive laws strongly dependent on the cofactor of the deformation.The development of this new conservation equation was possible due to the simplification introduced by the use of a tensor cross product operation, presented for the first time in the context of solid mechanics in [51] and further explained in [52].
The evolution of the above mentioned works resulted in a mixed methodology in the form of a system of Total Lagrangian first order conservation laws { p, F, H, J }.Following these developments, a new stabilised total Lagrangian SPH methodology was introduced in [26], aiming at the removal of spurious zero-energy modes and improved order of convergence of variables such as stresses and strains.Moreover, in this approach, linear consistency is fulfilled by performing kernel and gradient corrections, as proposed in [28].Global conservation of angular momentum is strictly ensured with the introduction of a monolithic angular momentum projection algorithm, first presented in [24].The stabilisation methodology presented in [26] does not require user-defined artificial stabilisation parameters and is based on physical pressure and shear wave speeds.
Even though the combination of all these developments has greatly improved the current state of SPH methods in fast solid dynamics, there is still the lack of a robust updated Lagrangian SPH framework, which is necessary when topological changes (i.e.fracture, separation) in the geometry need to be taken into account along the solution.With this idea in mind, the present work puts forward a new stabilised Updated Reference Lagrangian SPH framework, aimed to be a robust alternative tool for computer simulations in fast solid dynamics problems, prepared for the future endeavour of handling problems with topological changes.Specifically, and by adopting as referential configuration an intermediate configuration of the deformation process, the equation of conservation of linear momentum and three geometric conservation laws (for the deformation gradient, its cofactor and its determinant) are re-written leading to a very generic (incremental) system of first order conservation laws.One attractive feature of the methodology is its possibility to degenerate into either Total or purely Updated Lagrangian formulations.The actual use of the methodology in the context of fracture mechanics [53] is the next step of this work, beyond the current scope, along with the consideration of thermal effects induced by strong thermo-mechanical coupling due to shocks.
The paper is broken-down into the following sections.Section 2 starts by presenting the Total Lagrangian Formulation of the conservation laws to be solved, encompassing the linear momentum and three geometric conservation laws.The section continues deriving the new Updated Reference Lagrangian set of conservation laws in terms of a new set of alternative conservation variables and introducing their appropriate (entropy) work conjugates.Section 3 presents the variational statements of the problem.Section 4 presents the Smooth Particle Hydrodynamics discretisation scheme.Special attention is paid to two aspects: first, the Riemann based (upwinding) numerical dissipation employed; second, the introduction of four different options for the evaluation of kernels and gradients of kernels needed for SPH implementation.Section 5 briefly presents the three-stage explicit Runge-Kutta scheme used to evolve the semi-discrete equations in time, improving the numerical CFL restriction previously published in the literature.Section 6 presents a simple projection technique to ensure strict total conservation of angular momentum.Section 7 includes the algorithmic flowchart of the resulting numerical scheme and how to degenerate into their Total and purely Updated Lagrangian algorithmic counterparts.Section 8 presents a set of numerical examples to assess the convergence, conservation and stability of the computational framework, with detailed comparison with existing implementations available in the literature.Section 9 presents some concluding remarks.

A summary of the total Lagrangian formalism
Consider the three dimensional deformation of an isothermal body of material density ρ R moving from its initial undeformed configuration Ω V , with boundary ∂Ω V defined by an outward unit normal N, to a current deformed configuration Ω v at time t, with boundary ∂Ω v defined by an outward unit normal n (see Fig. 1).The time dependent motion φ(X, t) of the body can be described by a system of first order Total Lagrangian conservation laws [24][25][26]54,55], summarised as Here, U is the vector of conservation variables, F I is the flux column vector in I th material direction and S is the source term, described as (2) Fig. 1.Deformation of a solid body.
In above system of equations, p = ρ R v is the linear momentum per unit material volume, v represents the velocity field, f R is the body force per unit material volume, {F, H, J }1 represent the deformation gradient tensor, its co-factor and its Jacobian, P represents the first Piola-Kirchhoff stress tensor and the material identity tensor I is expressed as I = ∑ 3 I =1 E I ⊗ E I with the Cartesian material coordinate basis being defined as DIV and CURL represent the material divergence and curl operators carried out with respect to the material configuration, and the symbol represents the tensor cross product between vectors and/or second order tensors, see [52].
As the system of conservation laws presented above has more equations than needed, suitable compatibility relationships (also known as involutions [47,48,56,57]) are necessary, namely For closure of the system above, a suitable constitutive relationship for the stresses P compliant with the principle of objectivity (frame invariance) and the second law of Thermodynamics (Coleman-Noll procedure) is presented in Remark 1. Finally, for the complete definition of the initial boundary value problem, both initial and boundary (essential and natural) conditions must also be specified as appropriate.
Remark 1.In this work, and without loss of generality, a Mooney-Rivlin hyperelastic constitutive model has been considered, where the strain energy density is defined as a convex multi-variable function W of the triplet of deformation measures {F, H, J } [50,58,59] as where ζ , ξ and κ (bulk modulus) are positive material parameters.By comparison of the tangent elasticity operator at the initial configuration with that of classical linear elasticity, appropriate values for the material parameters ζ and ξ can be defined in terms of the shear modulus µ, that is, 2ζ + 3 √ 3ξ = µ [58].Following Refs.[24][25][26], the first Piola Kirchhoff stress tensor can be expressed as where the conjugate stresses {Σ F , Σ H , Σ J } are defined by and

Updated reference Lagrangian formalism
Taking inspiration from the work presented in Ref. [23,45,46], an equivalent system of Updated Reference Lagrangian conservation equations can be written in terms of a new configuration, additional to the well established material (initial) and spatial (current) configurations.This new configuration is defined as an intermediate or incremental configuration through which the continuum has passed (convected) during the deformation process at a previous time instant or load increment (refer to the so-called "Fixed" in time configuration in Fig. 1).This extra configuration can be adopted as a new reference configuration, leading to what we will refer to as a Updated Reference Lagrangian formulation.As we will see, this additional referential (intrinsic) configuration leads to an enhanced formulation, capable of degenerating into the classical Total Lagrangian or Updated Lagrangian formalisms.
By doing this, the triplet of deformation measures {F, H, J } (from material to spatial configurations) can be obtained via a multiplicative decomposition as where { f , h, j} represents the set of (incremental) deformations measured from the reference to the spatial configuration, whereas {F χ , H χ , J χ } denotes the set of (known) deformations measured from the material to the reference configuration.
In order to derive the Updated Reference Lagrangian formalism, we can utilise appropriate push forward transformations for vector and flux quantities [60], that is, multiplying material vector quantities {U , S} by J −1 χ and post-multiplying the Total Lagrangian flux F by H −1 χ , making use of the Piola identity Eq. (1) thus becomes which can then be particularised for the individual components of U , yielding Additionally, we introduce the following intermediate variables defined as With these, and with use of (9), above system (( 12)-( 15)) after some algebraic manipulation reduces to Here, div χ and curl χ represent the divergence and curl operators carried out with respect to the referential configuration, and the referential identity tensor i is being defined as i = ∑ 3 i=1 e i χ ⊗ e i χ with It is worth emphasising that the primary unknown variables of the Updated Reference Lagrangian system introduced herein are { p χ , f , h, j}.Crucially, the updates of incremental geometric deformations { f , h} (18) and ( 19) must ensure the satisfaction of appropriate involutions, namely Combining above equations into a system of conservation laws formulated in the reference configuration gives Here, U χ is the vector of conservation variables (per unit of reference configuration), F i χ is the flux vector in ith direction at reference domain and S χ is the source term (per unit of reference configuration).Their respective components are It is important to remark the simplicity of the above system of Updated Reference Lagrangian conservation equations (( 24)) with striking similarities to the previously presented Total Lagrangian system ((2)), which will naturally facilitate its computational implementation.

Hamiltonian (or total energy density) and constitutive relationship for stresses
In order to provide a proper physical meaning to the conjugate fields (with respect to the conservation variables U χ ) of the system at hand, consider the total energy density (or Hamiltonian H χ [61][62][63] per unit of reference volume in case of reversible processes) defined by and H χ (χ , t) and Ĥχ ( p χ , f , h, j) represent alternative functional representations of the same magnitude.For the isothermal case under consideration, the notion of Hamiltonian is simply the total energy per unit reference volume.Making use of Eq. ( 25), and introducing the conjugate stresses with respect to { f , h, j} to be the associated work conjugates V χ can then be obtained as Upon the use of the tensor cross product property shown in [50], the above conjugate stresses {Σ f , Σ h , Σ j } can indeed be related to the standard Cauchy stress tensor σ χ via the rate of internal work per unit reference volume, described as where the superimposed dot represents the usual material time derivative describing the change in f associated with a specific particle initially located at position X.Comparing the left hand side of Eq. (28a) with Eq. (28d) leads to the following relationship Remark 2. To complete the definition of the Cauchy stress tensor (29), one plausible option is to rewrite the referential conjugate stresses {Σ f , Σ h , Σ j } of the system in terms of their material counterparts {Σ F , Σ H , Σ J } (previously presented in Remark 1).This is achieved by first utilising the pull back equivalent of the convex multivariable function W χ ( f , h, j) = J −1 χ W (F, H, J ), followed by the use of the multiplicative decomposition for the triplet of deformation measures (9) and the chain rule.For instance, the conjugate stress of f gives Similarly, the conjugate stresses of h and j are now given by Remark 3.Many practical engineering applications often exhibit some permanent inelastic deformation.In order to describe this irrecoverable behaviour, the simple case of rate-independent von Mises plasticity in conjunction with either linear or nonlinear hardening rule will be considered in this paper.Within the context of large strains, it is customary to decompose the deformation gradient tensor F into an elastic component F e and a permanent deformation component F p , namely F = F e F p .In addition, and utilising the decomposition of F = f F χ , the elastic left Cauchy Green tensor b e can now be rewritten in terms of both incremental deformation gradient tensor f and referential plastic strain c p as The standard return mapping algorithm to update the (referential) plastic strain c p in order to ensure that the Kirchhoff stress satisfies, for instance, a von Mises type of plastic constraint can be seen in Ref. [26].

Weak form statements
In general, a standard weak variational statement for the isothermal mechanical system is established by multiplying the local differential Eqs.(23) (written in terms of the conservation variables U χ ) with their appropriate work conjugate virtual fields δV χ , and integrating over the reference domain Ω χ of the body, to give where the symbol • is used to denote the inner (dual) product of work conjugate pairs.Integrating by parts the last term on the right hand side of ( 33), and re-arranging the resulting equation yields ∫ where the normal fluxes are defined as being the outward unit normal to the reference domain in the ith direction.Above representation (34) can be particularised to the case of the linear momentum p χ and the triplet of incremental geometric deformation measures { f , h, j} as The main purpose of integrating by parts as shown above is to enable the imposition of the boundary conditions via boundary fluxes.This is indeed useful for the momentum update (35a) as it introduces naturally the boundary tractions t B , but less so in the case of geometric conservation Eqs.(35b)-(35d).

A point of departure: corrected SPH approximation at material domain
Consider the elastic body discretised by a cloud of particles as shown in Fig. 2. In the context of corrected SPH methods [20,28], any arbitrary vector function g is in general interpolated at any given position (quadrature point) via corrected (scatter) kernel functions W with a given compact support of radius 2h around every particle.For a given position X a , the vector g can be approximated as where ∥r ab ∥ = √ r ab • r ab being the distance to its centre and r ab = X a − X b the position vector.Typically, the material smoothing length h b is in the form of where f h is the (scalar) coefficient used to scale the size of the compact support.Here, Λ b a (or Λ j b ) represents the set of neighbouring particles b (or j) belonging to the domain of influence (or compact support) of a given radius 2h of particle a (or b), V b and g b represent the material volume and time-varying vector function g stored at particle b.In this paper, we employ the commonly used quadratic kernel function described in Ref. [64], with the value of scaling coefficient f h = 0.6.
To ensure the exact satisfaction of both constant and linear completeness, above corrected kernel approximation (scalar) operator C * [•] is defined according to [28] as with the parameters and In addition, for the evaluation of the material gradient of any arbitrary vector function g, we employ the following approximation introduced in [28 Specifically, the gradient correction (vector) operator C * [•] is defined as In expression (41), the term −g a is included in order to ensure that the gradient vanishes for a uniform field [28].
In addition, the use of the kernel gradient correction ∇0 ensures that the gradient of any linear field is exactly evaluated [10].
Given the fact that the kernel is usually a function of the distance ∥r∥ to its centre, the gradient evaluation can now follow via the chain rule, namely

Corrected SPH approximation at reference domain
For the case of the Updated Reference Lagrangian formalism, the kernel function as well as its gradient evaluation must be expressed in terms of the reference configuration.Since we aim at applications involving large deformations, standard updated isotropic (spherical) kernels would not be optimal (see e.g.[65]).Due to the possible presence of strong anisotropic changes in volume, and consequently in the particle spacing, the use of anisotropic kernels is better suited.Moreover, in such approach, provided that no severe topological changes take place, the initial list of neighbouring particles can be kept constant throughout the simulation, providing exceptional accuracy in a cost-efficient manner.This can be achieved by relating quantities in the material and reference configurations via the concepts of push forward φ * [•] and pull back φ −1 * [•] operations.For example, the reference kernel function W b (χ a ) can be considered as the push-forward equivalent of the material kernel function.This can be expressed in terms of the operation In the same spirit of Eq. ( 44), the push-forward equivalent of the material kernel gradient is2 Notice that the evaluations for the reference kernel function (44) and the associated gradient (45) still require geometrical information attached to the material domain, that is, r ab .In order to ensure that the kernel function depends solely on the information at the reference domain, one viable option is to approximate the material vector used in (44) as the pull back equivalent of its reference counterpart, that is where r χ ,ab = χ a − χ b denotes the reference position vector whereas the average (second order) approximation is . With this at hand, the reference kernel function can now be rewritten as Making use of ( 43) and the pull back quantities (46), we can further approximate the material gradient described in (45) to be In the current work, we exploit four different options to obtain reference kernel functions and their corresponding gradient evaluations, which are summarised as follows: • Option#1.Both the kernel function and its gradient are first corrected (ensuring zeroth-and first-order completeness) in the material domain and then are pushed forward to the reference domain.Mathematically, these are expressed as • Option#2.On the contrary, this particular option first pushes forward the material kernel function and its gradient to the reference domain, followed by the application of appropriate corrections at reference domain, that is • Option#3.As an alternative to option#2, we can rewrite the material kernel function (47) (in terms of the reference geometrical information) via the pull back equivalent of the reference vector (46) to yield • Option#4.Alternatively, we can re-define the push-forward operation of the gradient operator (45) using the average of the deformation mapping between pairwise interacting particles, that is φ For the numerical examples examined in this paper, the robustness and accuracy of the overall algorithm is not adversely affected regardless of the option used for SPH kernel approximation.However, we believe it is extremely informative to explore these various possibilities and provide the user with several options, in order to choose the most convenient one from the point of view of their computational implementation.

SPH semi-discrete equations
Upon the use of corrected SPH kernel approximations for U χ and δV χ , along with the corrected gradient evaluation for ∇ χ δV χ , in expression (34), the SPH discretisation for the system { p χ , f , h, j} described in ( 17)-( 20) becomes Here, the representation for the (algorithmic) internal nodal force T χ a is defined as with the (pseudo) area operators, see [26], defined as In the above expressions, V χ a and A χ a represent the volume and the referential tributary area and t a B its traction vector computed directly from the given traction boundary conditions.Notice that A χ a = 0 for those particles not placed on the boundary.
Finally, the remaining terms {D p χ ab , D j ab } in Eqs. ( 53) and ( 56) are the so-called numerical dissipation terms.These terms can be introduced and derived using the semi-discrete version of the classical Coleman-Noll procedure in order to ensure the production of non-negative numerical entropy.Following a similar procedure presented in Ref. [66], and after some algebraic manipulations, the Godunov-type numerical dissipations are summarised here for completeness where ) .Notice here that the stabilisation term applied to the linear momentum evolution (53) alleviates the appearance of spurious zero energy modes due to rank deficiency inherent to the use of nodal particle integration, whereas the stabilisation in the volume map update (56) may potentially be used to address pressure instabilities (especially in near incompressibility) if necessary.For all the numerical examples presented below, it is worth pointing out that no dissipation has been included in (56) by strongly enforcing the value of D j ab equal to zero.

Time discretisation
Insofar as the resulting set of semi-discrete equations is rather large, it will only be appropriate to employ an explicit type of time integrator.In this work, a three stage Runge-Kutta explicit time integrator [67] is used.This is described by the following time update equations from time step t n to t n+1 Additionally, the geometry is also updated via the same time integrator presented above, resulting in a monolithic time update procedure in which the unknowns U χ = ( p χ , f , h, j ) T are updated along with the (spatial) geometry x through Eqs.(61).Moreover, the maximum time step ∆t = t n+1 − t n is governed by a standard Courant-Friedrichs-Lewy (CFL) [47] condition given by where α C F L is the CFL stability number and c Ave p,ab the average pressure wave speed between particles a and b.Unless otherwise stated, a value of α C F L = 0.9 has been chosen in the following examples to ensure both accuracy and stability.
It is worth noticing that in expression (62), the characteristic particle spacing ∥x a − x b ∥ is computed with respect to the current (or spatial) configuration.The evaluation for pressure wave speed is however approximated by dividing the standard linear wave speed c Lin p by the average (minimum) stretch between the pairwise interacting particles a and b, described as where µ is the shear modulus, λ is the Lamé first parameter and λ a and λ b are the minimum stretch of particles a and b, respectively.

Conservation of angular momentum
The resulting SPH algorithm does not exactly fulfil conservation of angular momentum, since the strain measures { f , h, j} are no longer exclusively obtained from the current geometry.Although deviation from strict conservation is minimal, a monolithic discrete angular momentum projection algorithm is presented, following the work of [24,68].In our numerical experiments, we have observed that activation of this algorithm is not strictly necessary unless "machine accurate" strict angular momentum conservation is sought.Specifically, the local internal nodal force T χ a is suitably modified (in a least-square sense) in order to preserve the global angular momentum, whilst still ensuring the global conservation of linear momentum.
Adapting the procedure presented in [68] to the Runge Kutta time integrator considered herein, sufficient conditions for the global preservation of the discrete linear and angular momentum within a time step are explicitly enforced at each of the three stages described as A least-square minimisation procedure is used to obtain a modified set of internal nodal forces T χ a that satisfy the above conditions.This can be achieved by computing the minimum of the following functional (ignoring time arguments for brevity) After some simple algebra, a modified set of internal nodal forces T a is obtained in the form of where the Lagrange multipliers {λ ang , λ lin } are the solution of the following system of equations with the indicial notation

Algorithmic description
For ease of understanding, this section summarises the algorithmic implementation of the SPH methodology presented in previous sections.We start by including a graphical representation of the method, as illustrated in Fig. 3.
Algorithm (1) summarises the complete algorithmic description of the proposed Updated Reference Lagrangian SPH methodology.Notice that, under certain conditions, the proposed methodology degenerates to the case of either Total Lagrangian SPH [26] or pure Updated Lagrangian SPH.For example, if the conditional statement in Algorithm 1 is set to "FALSE" for every time integration process, the reference domain would remain unchanged throughout the entire simulation.In this case, the reference domain coincides with the material domain, implying that Fig. 3.The proposed Updated Reference Lagrangian algorithm comprises a series of multiplicative incremental configurations as illustrated above.For instance, from time t n to t 2n , the complete deformation gradient at time t 2n , namely F 2n , is evaluated via the previously computed deformation gradient F n (now being stored and assigned as a deformation gradient at intermediate configuration Ω v (t n )) and a series of incremental deformation gradients { f ∆ 1 , f ∆ 2 , . . ., f ∆n } between the reference configuration Ω v (t n ) and the new configuration Ω v (t 2n ).
the corrected kernel and gradient approximations are only carried out once.On the other hand, if the conditional statement returns "TRUE" at every time step of the time integration process, pure Updated Lagrangian SPH is retrieved.Specifically, corrections for kernel function and its gradient are performed at every time step, which may significantly increase the overall computational cost.One key feature of the proposed approach is the ability to suitably update the reference domain when certain criteria (e.g.internal state variable, maximum value of stretch) are met.This will be explored in forthcoming publications, especially in the application of fast dynamic fracture.However, in this work, we choose to update the reference domain at a fixed number n of time steps (pre-defined by user) with the aim to check the performance of this approach.Its graphical representation is illustrated in Fig. 3 for clarity.

Column
As a first example, a unit squared cross section column of length L = 6 m is studied.To assess the performance of the proposed method on various deformation modes (e.g.stretching, bending and twisting), the column is initiated with three different dynamic loading scenarios (see Fig. 4).These include: • Scenario I: tensile mode.The column is pulled with an initial (linear) velocity profile described as ] T with the value of V 0 = 50 ms −1 .The bottom surface of the column is clamped, the top surface is left free and the remaining boundaries are enforced with symmetric roller support.
• Scenario II: bending mode.The column is subjected to bending by the application of an initial linearly varying velocity profile given by v 0 = [ ] T where V 0 = 10 ms −1 is the maximum velocity applied.The column is clamped at the bottom and free on all other sides.

Algorithm 1: Updated Reference Lagrangian SPH Algorithm
Input : initial geometry X a and initial states of p χ ,a , f a , h a j a Output: current geometry x a , particle velocity v a and current states of F a , H a and J a (1) INITIALISE F χ,a = H χ ,a = I, J χ ,a = 1 and x a = χ a = X a (2) FIND neighbouring particles within a given support size (Λ b a ) (3) COMPUTE corrected kernel and gradient approximations for Time t 0 to Time t do if update at this step = TRUE then (4) COMPUTE the velocity as Aiming to display convergence of the algorithm, four successively refined (uniform) models are used.These include (M1) 5 × 5 × 25, (M2) 7 × 7 × 35, (M3) 9 × 9 × 49 and (M4) 17 × 17 × 97 number of particles.In terms of constitutive model, a neo-Hookean material is used.A summary of the data simulation used is also presented in Table 1.This problem is specifically designed to check if the proposed method can circumvent spurious modes in stretching dominated problems.In this case, four different options for gradient corrections (summarised in Section 4.2) ensuring zeroth and first order completeness are employed and compared.As shown in Fig. 5, all four options described above yield extremely similar results in comparison to the published Total Lagrangian result [26], almost undistinguishable.Moreover, the global numerical dissipation introduced in the algorithm can be assessed by computing the difference between the total conserved energy and the summation of kinetic and elastic strain energies, as reported for model M1 in Fig. 6.It is interesting to notice that all four possible options ensure the generation of global numerical dissipation over the entire duration of the simulation.Only slight variations are observed.Regardless of which option is used, the total numerical dissipation of the algorithm is reduced when increasing the number of particles.This is shown in Fig. 7.It is worth pointing out that, for problems accompanied with large topological changes (such as material separation in dynamic fracture), it is instructive to use Option#3 and Option#4 where the computations of kernel function and its gradient are carried out at (post-fractured) referential configuration.This is in contrast to Option#1 and Option#2 where the computations are at (pre-fractured) material configuration without considering separation process, hence not suitable for scenario involving severe topological changes.

Scenario II : Bending mode
The main objective of this problem [24,69,70] is to assess the applicability of the proposed algorithm by carrying out frequent updates of the reference configuration.Doing this with the classical SPH algorithm [23] would activate spurious mechanisms, and eventually lead to the breakdown of the scheme.For this reason, updates of the reference configuration are carried out at every time step.Obviously this is unnecessary but it has been done to check whether the algorithm triggers possible instabilities.Fig. 8 illustrates the deformed shape of a bending column where colour indicates pressure contour.Remarkably, even with a small number of particles (M1), the obtained result agrees very well with the Total Lagrangian result.No instabilities are observed.As shown in Fig. 9a, the amount of numerical dissipation reduces with an increasing number of particles.For completeness, we also simulate the same problem by updating the reference configuration at every 3, 30 and 300 time steps and the same results are obtained (see Fig. 9b).

Scenario III : Twisting mode
In order to examine the robustness of the algorithm, a well documented twisting example introduced in [24-26,71] is considered.Fig. 10 compares the deformation process of the structure at time t = 0.1 s using the four model refinements previously described.Similar results in terms of deformed shape and pressure field are observed.Finally, and for qualitative comparison, Fig. 11 monitors the evolution of the accumulated angle at four different positions, 3 namely N 11 , N 12 , N 13 and N 14 .Comparing the proposed algorithm with the published Total Lagrangian SPH algorithm [26], no noticeable difference is observed.

L-shaped structure
As previously explored in Refs.[25,59,66,72,73], the main objective of this classical benchmark problem is to examine the capability of the proposed algorithm of preserving both the linear and angular momenta of a system.An L-shaped structure (see Fig. 12) is subjected to an external torque induced by a pair of time-varying tractions acting on two of its boundary faces F 1 (t) and F 2 (t), described as In this example, a neo-Hookean model is considered.The values of all the parameters and the associated SPH ingredients used in the simulation can be found in Table 2.For completeness, three different levels of particle refinements are considered.More precisely, {M1, M2, M3} comprise {828, 5445, 13950} number of SPH particles, respectively.To assign the volume of each SPH particle, a tributary volume distribution is used.First, a particle refinement study for the proposed SPH algorithm is carried out.In Fig. 13, the deformation pattern of the structure predicted using a small number of particles (M1) agrees extremely well with the results obtained using finer discretisations (M2 and M3).For comparison purpose, a Total Lagrangian Upwind SPH [26] discretised with M3 is also included and compared.Practically identical pressure profiles are observed (see Fig. 13).
Second, Figs.14a and 14b examine the ability of the proposed algorithm (discretised with M3) of preserving both the linear momentum and angular momentum of the system.Specifically, the global linear momentum is close to (and oscillates around) zero machine accuracy at all times as no movement of the centre of mass is appreciated.The global angular momentum is however expected to be conserved after the loading phase, that is when time t > 5 s.Additionally, Fig. 14c illustrates the time histories of different forms of energy obtained using the proposed algorithm and the Total Lagrangian SPH.These include kinetic energy, strain energy and total energy, where the latter is defined as the sum of kinetic and strain energies.Identical results are observed.Given the fact that some consistent upwinding stabilisation terms are introduced in our SPH algorithm, a slight decrease in the total energy is unavoidable after the loading phase.As expected, numerical dissipation decreases with a progressive level of refinement, as reported in Fig. 14d.

Impact bar
In this example, a copper bar of length L = 0.0324 m and of radius R = 0.0032 m impacts against a rigid wall with a dropping velocity given by v 0 = [0, 0, −V 0 ] T , where V 0 = 227 m s −1 .The main purpose of the example is to assess the performance of the proposed algorithm in capturing large plastic flows under high speed impact [26,74].The geometry of the problem is illustrated in Fig. 15.A (Hencky-based) von Mises hyperelastic-plastic material, in conjunction with isotropic linear hardening rule, is chosen for the simulation of this problem.A summary of data simulation and material properties used are reported in Table 3.To assign the volume for every particle, a tributary volume is used.
Due to the symmetry of the problem, only a quarter of the geometry is represented and discretised using three levels of refinement, namely (M1) 1560, (M2) 3744 and (M3) 7280 particles.Fig. 16 compares the deformation process of the structure at time t = {40, 60, 80} µs using these three models.Similar deformed shapes are observed, with colour contours indicating equivalent von Mises stress and plastic strain profiles.In terms of the energy plot (see Fig. 17a), the kinetic energy of the system decreases upon impact.Such energy first transforms into elastic strain energy (prior to plasticity), followed by the addition of plastic dissipation (within the plastic regime).Moreover, given the nature of the proposed algorithm, a very small amount of the kinetic energy also converts to monotonic decreasing numerical dissipation during the deformation process, as displayed in Fig. 17b.For verification purposes, it is worth noticing that the results obtained using the proposed SPH algorithm match extremely well with those of the Total Lagrangian SPH counterpart.The final radius of the copper bar at time t = 80 µs predicted by the proposed algorithm is shown in Fig. 18b, benchmarked against other published numerical results [26,75]. 4As shown in Refs.[75], the solutions obtained using the standard linear 4-noded tetrahedra (being widely used in commercial software) typically suffers from volumetric locking and pressure instabilities.The proposed meshfree method clearly circumvents these issues.

Necking bar
The following benchmark example [76][77][78][79][80] illustrates necking of a cylindrical bar, of radius R = 0.006413 m and length L = 0.053334 m, under tension.The complete geometry of the bar can be viewed in Fig. 19.Notice that a geometric imperfection (one percent reduction in the radius) is introduced to induce necking in the central region   of the bar.The top surface of the cylindrical bar is stretched with smooth time-varying velocity profile described as with Here, the value of the parameters are t 0 = V 0 = V 2 = 0, t 1 = 0.0007 s, V 1 = 10 ms −1 and t 2 = 0.0014 s.Graphical representation of this bell-shaped velocity profile (69) is depicted in Fig. 19c.The material properties of the (Hencky-based) von Mises material and the simulation parameters are summarised in Table 4. Nonlinear hardening rule as proposed in [76] is chosen.
Due to the presence of symmetry planes, only a slice (of π 16 rad) of the upper half of the cylinder is simulated with appropriate boundary conditions.For completeness, two levels of particle refinement for the model are studied.Model M1 contains a number of 1428 SPH particles and model M2 comprises 5535 particles.The particle distributions are non-uniform and each particle carries a volume tributarily distributed from the associated    neighbouring particles.In order to accurately capture the onset of necking, more particles are placed within the crucial necking zone, see Fig. 20.Fig. 21 shows the time evolution of the deformation pattern with both plastic strain and equivalent von Mises stress plotted.The proposed method can clearly capture the formation of necking at the central region of the bar.Additionally, Fig. 22a depicts a comparison of the proposed SPH algorithm against the recently published Total Lagrangian SPH formulation [26] at the end of the simulation (that is, when the total elongation of the bar is of 14 mm).Remarkably, both formulations produce very similar results in terms of deformed shape and pressure field.For qualitative comparison, we also monitor the ratio between current and initial radius of the bar in necking regime as a function of the elongation.Indeed, our results agree extremely well with those of published experimental and numerical results [79,80].This can be seen in Fig. 22b.
For visualisation purposes, Fig. 23 displays a three dimensional view of the cylindrical bar at different levels of deformation.No spurious hourglassing is observed in the necking zone.

Conclusions
In this paper, a new Updated Reference Lagrangian Smooth Particle Hydrodynamics algorithm for the analysis of large deformation isothermal elasticity and elasto-plasticity has been presented.From the continuum point of view, by exploiting as referential configuration an intermediate configuration of the deformation process, the equation of conservation of linear momentum and three geometric conservation laws (for the deformation gradient, its cofactor and its determinant) have been rewritten leading to a very generic (incremental) system of first order conservation laws.Moreover, this can be degenerated into a Total Lagrangian system (previously pursued by the authors in [24][25][26]) or into a purely Updated Lagrangian system.Technically speaking, only the geometric conservation law for the deformation gradient (out of the three geometric conservation laws) is strictly necessary.An appropriate multiplicative decomposition of the conservation variables has resulted into an amenable mixed set of conservation equations with striking similarities to the conventional Total Lagrangian system, apart from the use of alternative referential (incremental) conservation variables and corresponding entropy (work) conjugate fields.
From the computational implementation standpoint, a second order (in space and time) entropy-stable SPH upwinding stabilisation method is employed by means of the use of the Rankine Hugoniot jump conditions, without resorting to any ad-hoc algorithmic regularisation.A new three-stage Runge-Kutta time integration method is implemented resulting in values of the CFL stability restriction close to one.With the aim of demonstrating the robustness and potential of the methodology, an ample suite of challenging problems (in elasticity and elastoplasticity) have been presented and thoroughly analysed.The amount of numerical dissipation introduced by the algorithm has been strictly monitored through the extra implementation of a conservation equation for the total energy density, though not strictly required from the modelling point of view.The use of alternative kernels and gradients of kernels as well the use of different strategies to select the referential (incremental) configuration has been explored.The next step of our work is the adaptation of the current Updated Reference Lagrangian SPH framework to handle dynamic fracture problems typically accompanied by complex topological changes caused by material separation.

Fig. 2 .
Fig.2.Illustration of the domain of influence (also known as compact support) for particle a.
[•] L ,R ab are the left and the right states of the reconstructed values obtained via linear reconstruction procedure [26].The positive definite stabilisation tensors are expressed in terms of the average pressure wave speed c Ave p,ab and the average shear wave speed c Ave s,ab as Ave p,ab n ab ⊗ n ab + c Ave s,ab (I − n ab ⊗ n ab ) ab = 1 2 ([•] a + [•] b ) and the unit vector is given by n ab = x b −x a ∥x b −x a ∥ .The pseudo-area vector C χ ,Skew ab (along with its norm magnitude ∥C χ ,Skew ab ∥) and its push forward equivalent (spatial) vector c Skew ab

( 8 )
COMPUTE corrected kernel and gradient approximations (9) COMPUTE σ χ ,a end (10) EVALUATE p and s-wave speeds: c p , c s (11) COMPUTE time increment: ∆t for RK time integrator = 1 to 3 do (12) COMPUTE slope of linear reconstruction procedure (13) COMPUTE right-hand-side of the mixed-based system: ṗχ,a , ḟ a , ḣa and ja (14) ENSURE conservation of angular momentum (15) COMPUTE smoothed velocities using the corrected kernel (16) EVOLVE p χ ,a , f a , h a j a and x a (17) COMPUTE σ χ,a end (18) COMPUTE smoothed variables using the corrected kernel (19) EXPORT results for this time step (20) ADVANCE in time end • Scenario III: twisting mode.The column is twisted with an initial sinusoidal velocity field relative to the origin given by v 0 = [0, Ω 0 sin ( π Y /2L), 0] T where Ω 0 = 105 rad s −1 represents the magnitude of initial angular velocity.Enforcement of boundary conditions is identical to that of Scenario II (bending mode).

Fig. 5 .
Fig. 5. Tensile mode.Pressure distribution for model M4 at two different time instants, namely t = 0.215 s (first row) and t = 0.3 s (second row).Results are obtained using the four options for the computation of the kernel function and its gradient presented in Section 4.2.Update of the reference domain is carried out at every time step of the time integration process.

Fig. 6 .
Fig. 6.Tensile mode.Evolution of total numerical dissipation for model M1 considering different options of kernel evaluations.Results are obtained by updating the reference domain at every time step of the time integration process.

Fig. 7 .
Fig. 7. Tensile mode.Evolution of total numerical dissipation for different model refinements using (a) Option#1, (b) Option#2, (c) Option#3 and (d) Option#4.Results are obtained by updating the reference domain at every time step of the time integration process.

Fig. 8 .Fig. 9 .
Fig. 8. Bending mode.Pressure distribution at time t = 0.45 s using four different model refinements, namely (M1) 5 × 5 × 25, (M2) 7 × 7 × 35, (M3) 9 × 9 × 49 and (M4) 17 × 17 × 97 number of particles.First four columns show the results obtained using the proposed algorithm with updates performed at every time step of the time integration process.The last column shows Total Lagrangian SPH results for comparison purposes.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 13 .
Fig. 13.L-shaped block.Pressure distribution at time t = 7.8 s using three different model refinements, namely (M1) 828, (M2) 5445 and (M3) 13 950 number of particles.First three columns show the results obtained using the proposed SPH algorithm with updates performed at every time step of the time integration process.The last column shows Total Lagrangian SPH results for comparison purposes.

Fig. 14 .
Fig. 14.L-shaped block.Evolution of (a) the components of global linear momentum, (b) the components of global angular momentum and (c) different forms of energy (e.g.kinetic energy K , strain energy Ψ and total energy E = K + Ψ ) using M3 model.In addition, the time evolution of (d) global numerical dissipation for three different model refinements is also monitored.Notice that [•] U L F refers to the results obtained via the proposed Updated Reference Lagrangian SPH algorithm and [•] T L F refers to the results from the published Total Lagrangian SPH algorithm [26].

Fig. 16 .Fig. 17 .
Fig. 16.Impact bar.A sequence of deformed structures at time t = {40, 60, 80} µs (from top to bottom) using {M1, M2, M3} (left to right).In terms of contour plot, von Mises stresses (left side) and equivalent plastic strain (right side).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .
Fig. 18.Impact bar.(a) Pressure comparison at 80 µs for M3.Results obtained via Total Lagrangian SPH (left side) and the proposed SPH algorithm (right side), and (b) evolution of the radius of the bar measured at the impact interface for three different model refinements.Notice that [•] U L F refers to the results obtained via the proposed Updated Reference Lagrangian SPH algorithm and [•] T L F refers to the results from the published Total Lagrangian SPH algorithm [26].

Fig. 21 .
Fig. 21.Necking bar.A sequence of deformed structures when the total elongation of the bar is of {10, 12, 14} mm (from top to bottom) using {M1, M2} (left to right).In terms of contour plot, von Mises stress profile (left side) and equivalent plastic strain (right side).

Fig. 22 .
Fig. 22. Necking bar.(a) Pressure comparison when the total elongation of the bar is 14 mm for model M2.Results obtained via Total Lagrangian SPH (left side) and the proposed SPH algorithm (right side), and (b) radius reduction in the necking region compared against published results [76].Notice that [•] U L F refers to the results obtained via the proposed Updated Reference Lagrangian SPH algorithm and [•] T L F refers to the results from the published Total Lagrangian SPH algorithm [26].

Table 1
Column: material and data simulation used.

Table 2 L
-shaped structure: material and algorithmic parameters used in the simulation.

Table 3
Impact bar: material and data simulation used.

Table 4
Necking bar: material and data simulation used.