An entropy-stable updated reference Lagrangian smoothed particle hydrodynamics algorithm for thermo-elasticity and thermo-visco-plasticity

This paper introduces a novel upwind Updated Reference Lagrangian Smoothed Particle Hydrodynamics (SPH) algorithm for the numerical simulation of large strain thermo-elasticity and thermo-visco-plasticity. The deformation process is described via a system of first-order hyperbolic conservation laws expressed in referential description, chosen to be an intermediate configuration of the deformation. The linear momentum, the three incremental geometric strains measures (between referential and spatial domains), and the entropy density of the system are treated as conservation variables of this mixed coupled approach, thus extending the previous work of the authors in the context of isothermal elasticity and elasto-plasticity. To guarantee stability from the SPH discretisation standpoint, appropriate entropy-stable upwinding stabilisation is suitably designed and presented. This is demonstrated via the use of the Ballistic free energy of the coupled system (also known as Lyapunov function), to ensure the satisfaction of numerical entropy production. An extensive set of numerical examples is examined in order to assess the applicability and performance of the algorithm. It is shown that the overall algorithm eliminates the appearance of spurious modes (such as hour-glassing and non-physical pressure fluctuations) in the solution, typical limitations observed in the classical Updated Lagrangian SPH framework.


Introduction
Total Lagrangian Smoothed Particle Hydrodynamics (SPH) [2][3][4] is a well-established numerical method for the simulation of explicit solid dynamics.One of the most attractive features of SPH is its mesh-free nature, thus not relying on the use of an underlying mesh.The absence of mesh, and the calculation of pair-wise interactions among particles based exclusively on their separation, allow ease of computation for problems involving large deformation.Due to its ability to handle large distortions within reasonable accuracy and stability [5][6][7], the SPH method has been shown to be competitive in comparison to alternative mesh-based methods, where latter would typically require expensive remeshing strategies.
However, for problems experiencing severe distortions, a Total Lagrangian SPH formulation [8,9] will unavoidably require updates of the material (or initial) configuration.Nonphysical zero-energy modes [10][11][12][13][14] are highly likely to be activated when performing such updates (i.e.re-calculation of kernel gradient and its correction).For instance, when very few updates are performed during the entire simulation, the accumulated errors may potentially remain small and unnoticed.However, when updates are frequently performed (for example, at every time step of the time integration process), the solution can be negatively affected resulting in spurious oscillations.
Significant efforts have been undertaken over the years to rectify undesirable zero-energy modes [15] in the SPH solution.Specifically, two types of techniques are widely used.First, the instability is precluded by introducing artificial viscous stress [16][17][18][19], or other related (Laplacian-type) techniques [20][21][22][23], conceptually similar to the treatment of hour-glass modes used in finite element method [10].Another common practice, in the framework of SPH, is the introduction of higher order derivatives [8,[24][25][26][27][28][29] to stabilise numerical computations.In [24], the approach is shown to be relatively expensive as it requires the evaluation of a third-order tensor for stabilisation based on Hessian difference.Second, the introduction of a staggered SPH approach [12,13,30,31] by incorporating a secondary set of particles (known as stress points) for stress computation, allowing the variables and their derivatives to be computed at different positions.Despite significant development in the field, the ab initio stability of SPH schemes still remains an open problem.
Moreover, even though the developments described above has greatly improved the current state of SPH methods, there is still a need for a robust Updated Lagrangian SPH framework, especially when attempting to model problems undergoing large geometry distortions, such as high-speed impact or high-speed stretching.Under this circumstance, consideration of thermal effects becomes necessary so as to attain a realistic representation of stresses.With this in mind, the aim of this paper is to further extend the recent SPH work [1] presented in the context of isothermal elasticity and plasticity models to account for possible strongly thermally-coupled scenarios, through the consideration of thermo-elasticity and thermo-visco-plasticity.Specifically, and by adopting referential configuration as an intermediate configuration during the deformation process, an extra conservation equation corresponding to the first law of thermodynamics (written in terms of the entropy density of the system) is solved in addition to the conservation of linear momentum and the three incremental geometric conservation laws (measured from referential domain to spatial domain).Interestingly, the methodology can indeed be degenerated into either a mixed-based set of Total [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] or Updated Lagrangian formulations [47] provided certain conditions are met.One key aspect that requires careful consideration is the overall stability of the algorithm.In the current work, an upwinding (Riemann-based) approach is exploited where a consistently derived numerical stabilisation is introduced guaranteeing the production of total numerical entropy over the entire simulation.The latter is shown by monitoring of the Ballistic energy of the system via the semi-discrete entropy analysis known as Coleman-Noll procedure.Another objective of the paper is to show that the overall SPH algorithm is capable of frequently carrying out updates of the reference configuration without introducing non-physical modes.In the numerical examples presented, unless otherwise stated, updates of the reference configuration are performed at every time step of the time integration process.Obviously this is not necessary but it has been done in order to assess whether the algorithm triggers possible instabilities.
The outline of the paper is organised as follows.Section 2 summarises the first-order system of Updated Reference Lagrangian conservation laws for isothermal hyperelasticity.Section 3 begins by introducing fundamental concepts of thermodynamics, necessary for the remainder of the paper.The section then presents two commonly used thermomechanical models, namely Mie-Grüneisen-based thermoelastic model and thermo-visco-plastic Johnson Cook model.Section 4 presents the variational formulation of the problem and the second law of thermodynamics written in terms of the Ballistic free energy.Section 5 presents the Smoothed Particle Hydrodynamics spatial discretisation where special focus is paid to the upwinding numerical dissipation employed.A proof of total entropy production (which is a summation of the physical dissipation introduced by the inelastic model and the numerical dissipation introduced by the SPH scheme) is included as a necessary condition for stability at the semi-discrete level.For completeness, Sect.6 illustrates the algorithmic flowchart of the resulting numerical scheme.Section 7 presents a number of challenging numerical examples with the objective to assess the robustness of the algorithm, where comparisons will be performed against an alternative mixed-based Total Lagrangian SPH implementation already benchmarked.Finally, Sect.8 presents some concluding remarks.

Updated reference Lagrangian conservation laws for reversible processes
Consider the deformation of a solid from an initial undeformed configuration V , of boundary ∂ V and outward unit normal N, to a current deformed configuration v , of boundary ∂ v and outward unit normal n, at time instant t.Consider an additional configuration χ , of boundary ∂ χ and outward unit normal N χ , corresponding to an intermediate configuration of the solid.This additional intermediate configuration can be adopted as a reference configuration, leading to what we refer to as an Updated Reference Lagrangian description [1].As a result, an Updated Reference Lagrangian system of first-order hyperbolic conservation laws can be used to describe the motion of a solid x = φ χ (χ, t) as follows (refer to [1] for an extended presentation): In the system above, p χ = ρ χ v is the linear momentum per unit reference volume, ρ χ represents the reference density, v represents the velocity field, f χ is the body force per unit reference volume, the triplet { f , h, j} represents the incremental deformation gradient tensor, its co-factor and its Jacobian (accounting for deformations from reference to spatial configurations).The referential stress tensor σ χ , obtained from the push forward (or push back) of the first Piola-Kirchhoff P (or Cauchy σ ) stress tensor to the referential configuration, is described as σ χ = J −1 χ P F T χ = jσ f −T .Symbol ⊗ represents the standard dyadic outer product, whilst denotes the tensor cross product [40] between vectors and/or second-order tensors.In addition, div χ and curl χ represent the divergence and curl operators carried out with respect to the referential configuration, and the respective identity tensor i is defined as i = 3 i=1 e i χ ⊗ e i χ with The incremental deformation tensor and its co-factor { f , h} (1b) and (1c) must satisfy appropriate compatibility conditions [1], namely Once equations (1b)-(1d) are solved and the triplet of incremental deformations { f , h, j} is obtained, the triplet of deformation measures {F, H, J } (mapped from material to spatial configurations) can then be obtained via multiplicative decomposition as where {F χ , H χ , J χ } denote the triplet of (known) deformation measures between the material configuration and the reference configuration.Notice that, if we update {F χ , H χ , J χ } continuously throughout the time integration process, a purely Updated Lagrangian first-order system [47] of conservation laws is retrieved.On the other hand, the Total Lagrangian formulation [32][33][34][35][36][37][38][39][40][41][42][43][44][45] is recovered if {F χ , H χ , J χ } are strongly enforced at the origin (that is, the reference configuration coincides with the material configuration).Detailed explanation of the transformations between the various formulations can be found in Reference [1].
3 Extension to irreversible processes

First law of thermodynamics in terms of total energy, internal energy and entropy
The system described in (1) can be generalised to take into account thermal effects [35,40], as is the case in thermoelasticity or thermo-visco-plasticity scenarios.The resulting process is irreversible and requires an additional conservation law (with corresponding unknown) describing the total energy balance.This is known as the first law of thermodynamics and, in the current work, is expressed in referential description as where q χ represents the heat flux per unit reference area and r χ is the heat source term per unit reference volume.Notice that the total energy density E χ and heat flux vector q χ can be related to the standard Lagrangian measures as , where E R and Q represents the total energy density per unit undeformed volume and material flux vector (per unit undeformed area), respectively.The total energy density E χ (χ , t) in the above equation includes both kinetic energy and internal energy contributions.Multiplying the linear momentum evolution equation (1a) by v and subtracting it from the above energy expression (5), the internal energy evolution equation becomes where e χ (χ , t) represents the internal energy per unit reference volume and the symbol ∇ χ represents the gradient operator evaluated at the referential domain, which is defined in indicial notation as [∇ χ ] I = ∂ ∂χ I .The internal energy density e χ (χ , t) is postulated to be a function of the incremental deformation variables X = { f , h, j}, the entropy density (per unit of reference volume) η χ and a set of state variables [48][49][50][51] (i.e.plastic deformation in this case) collected in the form of a tensor α, namely e χ (χ, t) = E χ (X, η χ , α). (7) Notice that E χ denotes the same internal energy density as e χ but with a different functional dependency.The entropy density field η χ (χ, t) is defined as the (energy) dual conjugate variable to the temperature θ(χ , t) described by Again, the pull back equivalents of both the referential entropy and the referential internal energy density are defined as Here, η and E represent the Lagrangian entropy and the internal energy per unit undeformed volume, whereas α R represents a set of state variables measured with respect to the material configuration.
Similarly, stress conjugate fields with respect to the incremental deformation measures { f , h, j} are defined as [1] Comparing the time derivative of the internal energy density e χ (χ, t) (6) to that of its equivalent re-expression E χ (X, η χ , α), and using the tensor cross product properties already presented in [41] together with expressions (8), ( 9) and (1b-1d), it is possible to relate the incremental stress tensor σ χ to the conjugate stresses defined in (9) as which leads to the following relationship It is also possible to re-express the first law of thermodynamics in terms of the entropy density η χ (χ, t) by combining (10) and ( 6) to give with ḊPhy representing the rate of physical dissipation introduced by the constitutive model, such as due to plasticity.This term is indeed zero when considering a reversible elastic model.Alternatively, and noting that 1 θ div χ q χ = div χ a conservation-type of law for the entropy density emerges as [40] Regarding the heat flux vector q χ , we consider a simple Fourier's law of heat conduction for an isotropic material as (14) with h representing the thermal conductivity coefficient calibrated in the spatial configuration.

General thermal relationship
In general, the Calorimetry relationships between internal energy density E χ , entropy density η χ and temperature θ can be derived [40] from the definition of the specific heat at constant reference volume c χ v .Specifically, where the specific heat can be alternatively expressed as c Here, ρ R represents the material density and C v and c v , respectively, represent the specific heat per unit mass and the specific heat per unit undeformed volume.Expression (15) can be re-written using the chain rule to yield Given the fact that ∂E χ /∂η χ = θ (8), a relationship between the temperature θ and the entropy density η χ at constant elastic deformation can be established after rearranging renders With expression above, and for simplicity assuming constant specific heat coefficient c χ v (such that it does not depend upon the elastic deformation and temperature), the relationship between entropy and temperature can be integrated exactly as [35,40] ηχ Reversing the above equation yields the expression for temperature field but now written in terms of {X, η χ , α} to give As shown above, the notation and θ is used to describe the same temperature with different functional dependency.
In addition, it is also possible to write a relationship for the internal energy density as functions of deformation X, entropy density η χ and a set of state variables α.This can be achieved by integrating expression (15) with respect to the temperature field between the limits θ R and a given temperature θ , and noting that θ(χ, t) = (X, η χ , α) (19) and Ẽχ (X, θ, α) = E χ (X, η χ , α) (15), to give

Stress evaluation
To complete the definition of the incremental stress tensor σ χ (11), and given that the constitutive relation in general depends on the standard deformation maps {F, H, J } (from material domain to spatial domain), it is convenient to utilise the previously described pull back equivalent of the internal energy density ẼR χ ( f , h, j, α) = J −1 χ ẼR (F, H, J , α R ) and of the entropy density function ηR χ ( f , h, j, α) = J −1 χ ηR (F, H, J , α R ).For instance, consider a simple volumetric-based Mie-Grüneisen model described by where q is a dimensionless parameter varying from zero (i.e. a perfect gas) to one (i.e.solid materials) and 0 is a (positive) material constant.It is also convenient to relate the stress conjugate fields { f , h , j } with those of a Total Lagrangian description { F , H , J } [35] defined as Indeed, for the conjugate stresses f , it yields Similarly, for the conjugate stresses h and j Remark 1 It is particularly useful to obtain stress expressions in terms of the symmetric Kirchhoff stress tensor {τ F , τ H , τ J } since it is usually needed when considering plasticity models.To achieve this, substitution of ( 4) into ( 23) and (24) gives alternative expressions for { f , h , j } to be described by with the stress relations being defined as In the current work, two well-established thermomechanical models, namely thermo-elasticity and thermaland rate-dependent Johnson-cook plasticity, will be presented and summarised in the following section.

Thermo-elastic model
For the case of a Mooney-Rivlin model, a standard distortional-volumetric internal energy density formulated at reference temperature θ R is described as [40] ẼMR where {ζ R , ξ R , χ R } are material parameters.These parameters can then be calibrated against those of linear thermoelasticity, namely, shear modulus μ, bulk modulus κ and thermal expansion coefficient α as [35] The material conjugate stresses (22) can now follow by taking the derivative of expression (27) with respect to {F, H, J }, and Notice that when the value of ξ R = 0, the above material degenerates to a neo-Hookean type of thermo-elastic model [39].It is now straightforward to obtain the components of the referential stress { f , h , j } by a direct substitution of the material stresses { F , H , J } (( 29) and ( 30)) into (23) and (24).

Thermo-visco-plastic model
Many engineering applications often exhibit some irrecoverable (or permanent) strain and thermal-dependent plastic deformation.To describe this behaviour, a von Mises plasticity model incorporating Johnson Cook hardening law [2] is considered and summarised here for completeness.In the context of large strains, it is customary to decompose the deformation gradient tensor F multiplicatively into an elastic component F e and a permanent deformation component F p as [52,53] This would subsequently lead to the evaluation of the elastic left Cauchy Green strain tensor b e , which is written in terms of the incremental deformation gradient tensor f and the inverse of the right Cauchy Green strain c −1 p measured at reference domain, described by As shown in Reference [52], the formulations developed to describe von Mises plasticity models are greatly simplified by operating in principal directions.For this reason, the left Cauchy Green strain tensor described in (32) can now be alternatively obtained by evaluating the principal directions of b e , that is n α , to give where λ e,α represents the elastic principal stretches.Recalling the distortional stretches being λe,α = J −1/3 λ e,α [52], it is now instructive to introduce the Hencky-based internal energy functional in terms of the elastic logarithmic stretches The volumetric and distortional components are and respectively.With these at hand, the principal components of symmetric Kirchhoff stress tensor are obtained in a standard manner [52] as with the deviatoric components ταα and pressure p emerging as which, after some algebra, gives 3 ; Attention is now focussed on the evaluation of the radial return mapping algorithm in order to ensure that τ stays on the yield surface function.In the current work, we consider the von Mises-based Johnson Cook hardening rule defined by a yield function of the deviatoric Kirchhoff stress τ and a yield stress τy (depending on the equivalent plastic strain εp , its plastic strain rate εp and temperature θ ) as The nonlinear strain-rate dependent hardening law is in the form of where Here, θ is the current temperature, θ melt is the melting temperature of the material and θ transition is the temperature at or below which there is no temperature dependence of the yield stress.Moreover, τ 0 y is the initial yield stress and ε0 represent the reference strain-rate.The remaining material constants are material hardening parameter H , hardening exponent N , strain-rate coefficient C and temperature exponent M. Notice that when θ > θ melt , the material is assumed to melt and behave like a fluid, offering no shear resistance since τy = 0.
An algorithmic representation of the one-step discrete time integration process (i.e. from n to n + 1) of the Henckybased model with von Mises rate-dependent plasticity model described above is summarised in Algorithm 1.In the case of Johnson Cook hardening rule, the plastic multiplier γ has to be obtained via the enforcement of the yield condition (40).This generally leads to the solution of nonlinear equations which would require an iterative Newton-Raphson method.In order to prevent singularities potentially arising from the derivative of a function within each Newton-Raphson iterative process, it is instructive to apply a change of variable by defining a new variable β in terms of plastic strain εp , plastic multiplier γ and Johnson Cook hardening exponent N , that is β = εp + γ N .For clarity, the iterative solution procedure is summarised in Algorithm 2. (5) Obtain yield criterion:

Combined equations
Combining equations ( 1) and ( 13) into a first-order hyperbolic system expressed in the reference configuration gives Algorithm 2: Evaluation of plastic multiplier through Newton-Raphson algorithm Here, U χ is the vector of conservation variables (per unit of reference configuration), F χ i is the flux vector in i-th direction at reference domain and S χ is the source term (per unit of reference configuration).Their corresponding components are

Variational formulation
In order to provide a proper physical implication to the conjugate fields of the first-order system (43) at hand, we introduce the Ballistic energy density B χ [40] per unit of reference volume (also known as the Lyapunov function of the thermomechanical process) defined by with B χ (χ , t) and Bχ ( p χ , f , h, j, η χ , α) being alternative functional representations of the same magnitude.In the above equation, the first term of the right-hand side represents the kinetic energy, the second term represents the internal energy density and third term represents the thermal heat component.Recalling the definition of conjugate stresses (9), it is now possible to obtain the associated work conjugates V χ as [35] where ϑ = −θ R denotes the temperature change.With this, the standard weak statement [54,55] of the underlying system is established by multiplying the differential equations (43) with their suitable 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 product of work conjugate pairs.In order to introduce physical boundary contributions acting on the body, the flux term on the right hand side of ( 47) is now integrated by parts and the resulting equation yields Here, the normal fluxes are defined as with N χ i being the outward unit normal to the reference domain in the i-th direction.Above general representation (48) can be particularised to the conservation equations under consideration, namely the linear momentum p χ , the triplet of incremental geometric deformation measures { f , h, j} and the entropy density η χ as Notice here that δϑ = δ is the virtual conjugate field of the entropy η χ .The objective of integrating by parts as shown above is to enable the enforcement of boundary conditions via physical boundary fluxes.This is especially useful for the momentum update (49a) and the entropy density equation (50) as both expressions naturally introduce the boundary tractions t B , boundary heat flux q B and boundary temperature θ B .

Second law of thermodynamics
It is instructive to revisit the global version of the second law of thermodynamics when written in terms of the Ballistic energy density B. Taking the derivatives with respect to its arguments, the time derivative of the Ballistic energy density is obtained via the chain rule as where, equations ( 46), (1b-1d) and (11) have been substituted in the second, third and fourth lines of (51), respectively.Consequently, we can substitute the linear momentum equation (1a) into (51) to give Moreover, let us now focus on the term associated with the time rate of the entropy density.By replacing δϑ with ϑ in expression (50), and noticing that the term involving entropy density rate on the right-hand side of (53) becomes Combining ( 55) and ( 53), and carrying out integration by parts of the div χ term in equation ( 53), it after some rearrangement renders where ˙ ext χ denotes the mechanical power associated with external forces, defined as and Q ext χ represents both heat source and heat flux added (removed) to (from) the system, defined as Recalling the Fourier's law of heat conduction ( 14), the first term on the right-hand side of ( 56) is non-positive, which is demonstrated as below Additionally, consider the case of elasto-plasticity [49] where the elastic potential energy (34) is expressed in terms of elastic left Cauchy-Green tensor b e = f c −1  p f T (with . Under this circumstance, the state variable is in fact the inverse of the plastic right Cauchy Green tensor (with respect to reference configuration), that is With this, the rate of plastic (physical) dissipation ḊPhy described in (56) becomes Insofar as the time rate of plastic strain ε p χ has been defined as the work conjugate to the von Mises equivalent stress τ [48], equation above can then be recast as [49] where τ is the deviatoric component of the Kirchhoff stress and the transformation ε p χ = J −1 χ εp .Observing that in the above expression the rate of dissipation is always nonnegative, that is ḊPhy ≥ 0, equation ( 56) can be transformed into the following inequality This represents a valid expression for the second law of thermodynamics [48] of a system.Satisfaction of inequality ( 62) is a necessary ab initio condition to ensure stability, otherwise referred to as the classical Coleman-Noll procedure [56].This fundamental concept will be further exploited in Sect.5.2 when introducing consistently derived numerical dissipation to the SPH discretisation.

SPH semi-discrete equations
Combining the use of nodal (or particle) integration for approximating the weak form integrals (49) and the standard corrected gradient evaluation for ∇ χ δV χ [35,57] to ensure zeroth-and first-order completeness, that is the SPH discretisation for the system { p χ , f , h, j, η χ } described in (49a-49d) and ( 50), after some algebraic manipulation becomes where the pair-wise internal force T χ ab is defined by with the pseudo area operators, see Reference [34], being defined, respectively, as In the above expressions, b a represents the set of neighbouring particles b belonging to the domain of influence of a given radius 2h of particle a, A χ a and V χ a represent the referential tributary area and the volume.The boundary traction t a B is directly computed from the given (Neumann) boundary conditions, whereas the heat flux q B and the boundary temperature θ B are the prescribed thermal boundary conditions.Note that A χ a = 0 for those particles not located on the boundary.It is also worthwhile pointing out that even with the use of kernel gradient correction (69), expression (63) still ensures the global conservation of linear momentum due to the pair-wise nature of internal force representation (68).
Finally, in order to address non-physical zero-energy modes due to the rank-deficiency inherent to the use of nodal integration (e.g.collocation), appropriate numerical dissipation terms {D p χ ab , D j ab } (refer to expressions (63) and (66)) must be introduced.These terms, being locally conservative by construction, can be suitably derived utilising the semi-discrete version of the second law of thermodynamics written in terms of the Ballistic energy (62), guaranteeing non-negative entropy production.This will be demonstrated in the following section.It is interesting to note that the stabilisation term incorporated to the linear momentum evolution (63) addresses the appearance of hourglass modes due to rank deficiency, whereas the stabilisation in the Jacobian evolution (66) would be used to remove pressure fluctuations.
Since the resulting set of equations is rather large, it will be suitable to employ an explicit type of time marching scheme.
In this work, a three-stage Runge-Kutta explicit time integrator presented in [1] is used.

Numerical entropy production
In this section, inequality ( 62) is assessed for the set of SPH semi-discrete equations described in (63)-(67).The semidiscrete form of ( 51) is where, equations (65), (66) and (11) have been substituted in the first and second lines of (70), respectively.Consequently, we can substitute the evolution of linear momentum equation (63) and of the first law (67) into (70) and, after some simple algebra, gives Here, ˙ ext and Q ext denote the semi-discrete power contribution and total heat contribution, respectively, expressed as Under the framework of variational consistency [33,58], the first term on the first line of (71) must be zero.This can be easily proved if the discretisation of velocity gradient for internal work is consistent with the discretisation for the incremental geometric conservation equations, which indeed is our case as shown below It is now the objective to prove that the remaining terms on the right-hand side of (71) are non-positive (to be in agreement with inequality (62)).With respect to the heat conduction term (second term on the first line of (71)), where the last inequality is fulfilled due to the nature of the conductive heat flux.With respect to the physical model dissipation term (the first term on the second line of (71)), it is again non-positive due to the definition of the rate of plastic dissipation, that is As for the numerical dissipation term (second term on the second line of (71)), this term can be equivalently written by swapping indices a and b to give Simple averaging the first and second terms of the expression above, and noting the local conservation nature of the stablisation terms such as D Dissipation terms remained to be defined in order to guarantee non-negative total entropy production (and thus, the fulfilment of the second law of thermodynamics).Sufficient conditions to ensure this, namely D total > 0, are given by Interestingly, the dissipation terms are directly related to the jump (or difference) in velocity and stresses between pairwise particles, typical upwind method [49] of first-order Godunov-type scheme.In order to ensure second-order accuracy in space, and following our previous work [34], a linear reconstruction procedure based on the use of corrected kernel gradient operator is used for the reconstruction of the left and right states at the mid-edge connecting particle a and particle b.In addition, we have also implemented the classical monotonicity-preserving Venkatakrishnan slope limiter [59] to better handle spurious oscillations in the region of shocks.

Algorithmic description
For ease of implementation, Algorithm (3) recaps the algorithmic description of the Updated Reference Lagrangian SPH methodology for thermo-mechanical coupled problems at finite strains.One interesting feature of the proposed SPH algorithm is the ability to suitably update the reference configuration when certain criteria are met.This indeed will be explored in forthcoming publications, especially in the area of dynamic fracture in brittle materials where the principal stress criterion of Rankine will be employed.However, in this work, we decide to update the reference domain at every time step of the time integration process.This is to check whether the proposed algorithm is capable of removing undesirable spurious modes, a typical shortcoming of the standard Updated Lagrangian SPH formulation [24].

Numerical examples
In this section, a series of three-dimensional numerical examples is presented in order to assess the performance, effectiveness and applicability of the proposed Updated Reference Lagrangian Smoothed Particle Hydrodynamics (URL-SPH) algorithm described above.It is crucial to show that the overall URL-SPH formulation • achieves equal second-order convergence for velocities, stresses and temperature, • alleviates spurious oscillations in the region of shocks (or discontinuities), • circumvents zero-energy modes (under dynamic stretching) and pressure instabilities, • preserves the total linear and angular momenta over a long term response, and • guarantees a non-negative rate of production of total entropy within a coupled system.
In the following numerical computations, the global a posteriori angular momentum projection algorithm as shown in [1] is activated.Moreover, the kernel function as well as its gradient evaluation must be expressed in terms of the reference configuration.To achieve this, we first map the domain of interest from reference domain to material domain, perform the necessary calculations (e.g. both kernel and its gradient evaluation), and then push forward to Mie-Grüneisen coefficients q 1 0 8.5889 the reference domain followed by the application of appropriate corrections for completeness.Details can be found in Sect.4.1 of Reference [1].In terms of the temporal stability of the algorithm, the Courant-Friedrichs-Lewy number of 0.9 has been chosen [1].In addition, comparisons are also carried out against the results simulated using an alternative in-house Total Lagrangian SPH algorithm [35], which can indeed be recovered by ensuring no updates of the configuration takes place over the entire simulation, that is It is not the premise of the paper to claim that URL-SPH algorithm outperforms our previously developed Total Lagrangian SPH algorithm [35], but to demonstrate that the current URL-SPH algorithm can indeed be equally compelling and competitive in the applications of solid mechanics.This in general is not necessarily the case for standard Updated Lagrangian SPH algorithms [24], where non-physical zero-energy modes can accumulate in the solution over time and eventually lead to the breakdown of the numerical scheme.

Swinging cube
To check the convergence pattern of the proposed URL-SPH algorithm, we consider a unit cube subjected to both mapping and temperature profiles described, respectively, as and with the parameter β being defined as With these equations at hand, and to guarantee the structure is in equilibrium state, the exact profile for the heat source term r R becomes whilst the body force term remains zero, that is f R = 0.The complete derivation procedure was detailed in [40].In this example, a linear thermo-elastic model1 is considered with the parameters summarised in Table 1.
Regarding boundary conditions of the cube, we enforce symmetry boundary conditions (i.e.restricted to tangential movement) at the faces X = 0, Y = 0 and Z = 0 and anti-symmetry boundary conditions (i.e.restricted to normal movement) at the faces X = 1, Y = 1 and Z = 1.Additionally, reference temperature θ R must be enforced at every time step of the time integration process at those three boundary faces, namely X = 1, Y = 1 and Z = 1.
Comparing with the exact expressions provided in (80) and (81), Fig. 1 illustrates the L 2 global convergence analysis of the overall SPH algorithm at time t = 8 × 10 −4 s.Indeed, one crucial advantage of the proposed SPH formulation over the standard (displacement-based) SPH algorithm is the ability to achieve equal second-order convergence for all the variables solved, namely linear momentum (or velocity), the stress tensor (or strain) and the temperature (or entropy).

Cable with step function loading
A wave propagation of a cable under the influence of shocks is considered.The purpose of this test case is to show the shock capturing capability of the proposed SPH algorithm.Similar type of problems were also explored in References [60][61][62].A cable of length L = 10 m, with a unit cross section A = 1 m 2 , is fixed at the left end (X = 0), whilst a step traction loading is enforced at the right end (X = 10 m) given as with T 0 = 0.001 Pa.Roller support (also known as symmetry boundary conditions) are applied on the remaining boundaries.For simplicity, no thermal effects are included.First, we consider a linear elastic model where the material properties are Young's modulus E = 1 Pa, material density ρ R = 1 kg/m 3 and Poisson's ratio ν = 0.The closed-form expression for the horizontal displacement of the cable is given as a function of position X and time t as where the natural frequencies ω n are given by Upon the sudden application of external force at the right end of the cable, a shock stress wave propagates towards the left fixed end and then gets reflected back.Figure 2a and b monitor the time history of both the horizontal velocity and the axial stress wave measured at the middle of a cable, that is when X = 5 m.As expected, first-order URL-SPH algorithm (without kernel gradient evaluation for reconstruction) introduces slightly more numerical dissipation to the solution.To further improve the solution accuracy, a linear reconstruction procedure is used.Minor under-and over-shoot oscillations are seen in the presence of shocks.The oscillatory behaviour however can be addressed with the introduction of appropriate slope limiter.Moreover, a particle refinement study is also carried out.Three model refinements are used, namely (M1) 606, (M2) 1206 and (M3) 2406 number of SPH particles.
As displayed in Fig. 2c,d, improved representation of shock profile is clearly seen by increasing the number of particles.Second, we examine the exact same problem but now with Johnson-Cook plasticity model (41).Again, let us consider isothermal condition which requires g(θ ) = 0 to be enforced in the yield function (41).The remaining parameters used in the Johnson-Cook model are N = 1, C = 0, τ 0 y = 0.0015 Pa and the hardening parameter H = 0.05 Pa.For benchmarking purposes, the results obtained using the in-house vertex-based finite volume algorithm [49] is also plotted and compared.The proposed URL-SPH algorithm combined with slope limiter removes unwanted oscillatory behaviour and, more crucially, is in excellent agreement with the finite volume results [49].This is seen in Fig. 3.

L-shaped block
As documented in References [1,35,40], the primary aim of this test case is to assess the capability of the proposed URL-SPH algorithm in the preservation of both the linear and angular momenta of a system.The geometry of the problem is illustrated in Fig. 4. The block is subjected to a pair of (time-varying) boundary forces F 1 (t) and F 2 (t) which can be mathematically described as When considering thermo-mechanical coupled physics, we also need to prescribe an initial temperature distribution Mie-Grüneisen coefficients q 1 0 0.0255 across the structure defined as elsewhere.
This can be equivalently done by enforcing the associated entropy profile via the relation between the entropy and temperature (refer to ( 18)), that is A Mie-Grüneisen-based thermo-elastic model as detailed in Sect.3.3.1 is considered.The associated material parameters used in the simulation can be found in Table 2.For completeness, three different levels of particle refinement are considered: {M1, M2, M3} comprising {828, 5445, 13950} number of particles, respectively.
First, a particle refinement study is carried out.This can be seen in the first three columns (from left to right) of Fig. 5.The deformation pattern, together with pressure and temperature profiles, simulated using a relatively small number of particles (M1), is in good agreement with those results obtained using finer discretisations (M2 and M3 models).For benchmarking purposes, an alternative in-house Total Lagrangian SPH algorithm [35] with M3 discretisation is also included and compared.Comparing the results of the proposed algorithm and those of the Total Lagrangian SPH algorithm, near-identical results are observed (see Fig. 5).
Second, Figs.6a and b shows the ability of the proposed algorithm in ensuring the conservation of global linear and angular momenta.The global linear momentum, L total = χ p χ d χ , is expected to oscillate around zero value (machine error) at all times.The global angular momentum, A total = χ x × p χ d χ , is indeed conserved after the loading phase t > 5 s.Another interesting variable to be monitored is the global entropy η total χ = χ η χ d χ , which increases over time throughout the entire simulation.This is seen in Fig. 6c, indicating the discrete satisfaction of second law of thermodynamics.In addition, Fig. 6d depicts the evolution of various energy measures.These include kinetic energy In this case, the external power only arises as a result of external boundary traction, that is ψext = ∂ χ t B •v B d A χ .With these at hand, the total energy E total = K total +E total Mech +E total Ther −ψ ext can now be computed.This consequently leads to an alter-Fig.5 L-shaped block: comparison of deformed shapes at time t = 15 s.The first three columns (left to right) show the particle refinement of a structure simulated using the URL-SPH algorithm, whereas the last column (on the right) shows a deformed structure via alternative in-house Total Lagrangian SPH algorithm [35].The first row depicts the pressure contour and the second row illustrates temperature contour.A Mie-Grüneisen thermo-elastic constitutive model as described in Sect.3.3.1 is used.The corresponding material parameters are summarised in Table 2 native energy measure known as Ballistic energy, that is B total = E total −θ R η total χ .A slight decrease in the total energy is unavoidable after the loading phase due to the incorporation of upwinding-based numerical dissipation (76) to the system.
Third, and for qualitative comparison purposes, Figs.6e and f monitor the time evolution of the velocity component v x and the temperature at position X = [6, 0, 0] T .The solution converges with a successive level of refinement.Finally, a sequence of deformed states are depicted in Figs.7 and  8, where the colour contour plot indicates the pressure and temperature distributions, respectively.Stable solutions are observed even after a relatively long-term response.

Punch test
We consider a block with 3×3 of vertical holes with diameter D (see Fig. 9a).The block is left free on its top surface and is constrained with roller supports (i.e.symmetric boundary conditions) applied to the rest of the surfaces.The objective of this test case [40] is to check if the proposed URL-SPH algorithm is capable of alleviating the appearance of spurious pressure in a highly constrained scenario.The deformation of the block is initiated with a compressive velocity profile applied in a quarter of the domain (X ≥ 0 and Y ≥ 0), together with a linear temperature profile, described as    In Fig. 14a, the time history of the kinetic energy, internal energies (e.g.heat and mechanical contributions), total energy and Ballistic energy is monitored.Specifically, the difference between the total energy (cyan dashed line) and the Ballistic energy (green dashed line) is regarded as the global entropy associated with irreversible heat conduction, which must be positive in this case.This is shown in Fig. 14b as the value of global entropy is non-negative and increases over time.Moreover, the standard displacement-based SPH algorithm [24] was reported to trigger possible instabilities by carrying our frequent updates of the reference configuration.This would then lead to the breakdown of the overall scheme.
To highlight this issue, we update the reference configuration at every 1, 10, 20 and 30 time steps.We then monitor the temperature evolution at position X = [0, 0, 0.5] T m and also the evolution of global dissipation.Almost identical results (i.e.stable and smooth) are observed (refer to Fig 14c and d).
For visualisation purposes, Figs. 15 and 16 display a series of deformed states without exhibiting locking.

Taylor bar impact
This is another benchmark problem where a copper bar of initial length L = 0.0324 m and of initial radius r = 0.0032 m impacts against a rigid wall with a velocity of 227 m/s.The initial temperature profile of the bar is set to the reference temperature, that is θ(χ , t = 0) = 298.15K.The geometry of the problem is illustrated in Fig. 17.Its objective is to assess the performance of the proposed SPH algorithm under high speed impact.Specifically, a von Mises material with strain rate-and thermal-dependent Johnson-Cook hardening law is chosen.The material parameters used in the simulation are    3 tabulated in Table 4.For ease of computation, we perform the simulation of the bar impact by considering only a quarter of the domain via appropriate symmetry boundary condition such as roller support.
Aiming to demonstrate the consistency of the algorithm, we discretise the quarter of a bar using three different levels of particle refinement, namely (M1) 1560, (M2) 3744 and (M3) 7280 number of SPH particles.In Fig. 18a, the evolution in time of various energy measures is monitored.When impact occurs, most of the kinetic energy is converted into irrecoverable heat dissipation and plastic dissipation whilst part of the kinetic energy is transferred into elastic strain energy.Fig-ure 18b shows the reduction of total numerical dissipation when increasing the particle density.The global numerical entropy increases over time for the entire simulation, hence guaranteeing long term stability.We also monitor both radius and length evolution at position X = [0, 0.0032, 0] T and X = [0.0032,0, 0.0324] T , respectively.This is seen in Fig. 18c and d.Our results agree extremely well with the solutions obtained via the Total Lagrangian SPH counterpart [35].This is usually not the case when employing the standard displacement-based Updated Lagrangian SPH algorithm [24] as it was reported to introduce particle clumping at the contact boundary.3 In addition, and as shown in Fig. 19, the deformation pattern of the structure together with its temperature and von Mises contour matches extremely well across all the three particle refinements.For visualisation purposes, Figs.20 and 21 show a sequence of deformed states of the bar for a rela- Finally, we can further assess the robustness of the algorithm by substantially increasing the value of the initial temperature profile θ(χ, t = 0) from 298.15 K to 573.15 K .Due to the softening behaviour caused by the high temperature accumulated at the contact plane, the time history of the length and radius are relatively larger in comparison to that of the previous case.Their plots are shown in Figs.22c and  d.Finally, Fig. 23 displays a series of snapshots for the bar impact in terms of temperature distribution, simulated via the proposed method and the Total Lagrangian SPH algorithm.Practically identical results are observed.

Necking bar
Similar to the Taylor bar previously explored in Sect.7.5, we now stretch the bar on both sides by reversing its initial velocity field.In this example, no fracture is considered.To account for thermal effects, the initial temperature profile of the bar is prescribed as θ(χ, t = 0) = 573.15K.The primary interest of this necking problem is to demonstrate that the proposed URL-SPH methodology is capable of alleviating spurious modes even for problems involving massive stretching, a persistent shortcoming typically encountered in the classical Updated Lagrangian SPH algorithm [24].The material properties used in the simulation are exactly the same as those reported in Table 4.A strain rate-and thermal-dependent Johnson-Cook hardening rule is chosen.Given the presence of symmetry planes, only one-eighth of the structure is discretised with appropriate boundary condi- Figure 24 illustrates a comparison of the proposed SPH algorithm against our in-house Total Lagrangian SPH algorithm [35] at time t = {20, 40} ms.Both formulations yield similar results in terms of deformed shape and pressure field.For qualitative comparison, we also monitor the radius reduction of the bar in the necking region as a function of the elongation.Comparing with the Total Lagrangian counterpart [35], it is interesting to notice that the proposed SPH method is able to capture post-necking deformation behaviour with a smaller number of particles.This is seen in Fig. 25b. Figure 26 compares the pressure resolution with and without the dissipation term appearing in the incremental Jacobian conservation equation (66).Observe that, with the incorporation of sufficient dissipation D j ab , the overall algorithm removes the non-physical mechanisms similar to hour-glassing in the solution.Figure 27 shows a series of deformed states where unwanted low-energy modes are not present.No temperature and/or plastic strain instabilities are observed.4 plasticity has been presented.From the continuum viewpoint, the methodology is built upon a suitable multiplicative decomposition of the conservation variables by introducing an intermediate (or incremental) configuration during the thermally coupled deformation process.This requires the re-formulation of a system of first-order hyperbolic conservation laws, usually expressed in material (or initial) configuration, to this new intermediate configuration.In addition to conservation laws for the linear momentum and the three incremental geometric conservation laws (for the deformation gradient, its co-factor and its determinant) previously used in isothermal process [1], a further conservation law representing the first law of thermodynamics written in terms of the entropy density is incorporated to extend the range of applications into thermally coupled hyperelasticity and strain-rate dependent plasticity.
From the spatial discretisation standpoint, a Smoothed Particle Hydrodynamics method utilising the standard kernel gradient corrections is presented.In order to address spurious energy modes inherent to the collocation nature of SPH, appropriate upwinding numerical dissipation is introduced.Such numerical dissipation is specifically designed via the use of the Coleman-Noll procedure at the semi-discrete level, demonstrated by monitoring the so-called Ballistic energy of the system.From the time integration standpoint, a standard explicit three-stage Runge-Kutta time marching scheme is employed.With the aim of demonstrating the reliability of the methodology, a wide spectrum of numerical examples is presented and compared.It has been shown that the resulting SPH algorithm addresses several numerical artefacts posed by standard Updated Lagrangian SPH methods, namely spurious pressure fluctuations, hour-glassing and numerical errors related to global conservation, completeness and long-term instability.4 Fig. 27 Necking bar: a sequence of deformed structures at times t = {0, 8, 16, . . ., 72} μs (left to right and top to bottom).In each subfigure, and in terms of colour representation, top half of the bar represents plastic strain and bottom half of the bar represents temperature profile.Results obtained using a Hencky-based Johnson-Cook hardening rule as described in Sect.3.3.2.The corresponding material parameters are summarised in Table 4

end ( 7 )
Update deviatoric Kirchoff stress tensor: alternative expression for D total is

1 2
defined as positive semi-definite stabilisation matrices[1] Ave p,ab n ab ⊗ n ab +c Ave s,ab (I −n ab ⊗n ab ) ; •] a + [•] b ) and the direction vector is given by n ab = x b −x a x b −x a .c s and c p correspond to the shear and volumetric wave speeds obtained via the classical wave propagation theory[49].In addition, 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 are, respectively, defined as C χ ,Skew ab = 1 2 C χ ab − C χ ba and c Skew ab = h a C χ ab − h b C χ ba .

Fig. 1
Fig. 1 Swinging cube: L 2 global convergence analysis at time t = 8 × 10 −4 s for (a) the components of linear momentum (or velocity), (b) the components of the stress tensor and (c) the temperature.Results obtained using a linear thermo-elastic model and the material properties used are summarised in Table 1

Fig. 2 2 Fig. 3
Fig. 2 Shock dominated problem: the first row illustrates the time history of (a) horizontal velocity component v 1 and (b) axial stress σ 11 χ measured at the middle of a structure.The second row shows the particle refinement analysis for (c) velocity component v 1 and (d) axial

Fig. 4
Fig. 4 Problem setup for L-shaped block: a geometry and b initial temperature profile

Fig. 6 L
Fig. 6 L-shaped block: time evolution of (a) global linear momentum, (b) global angular momentum, (c) global entropy, (d) different energy measures plotted with two different scales, (e)velocity component v x at position X = [6, 0, 0] T and (f) temperature at position X = [6, 0, 0] T .In terms of the energy plot, the magnitudes of internal heat energy and

Fig. 10
Fig.10Punch block: the first three columns (M1, M2 and M3, from left to right) show the a sequence of particle refinement of a structure simulated using URL-SPH algorithm, whereas the last column shows a deformed structure via in-house Total Lagrangian SPH algorithm[35].The first and third rows depict the pressure and temperature contour

Fig. 11
Fig. 11 Punch block: comparison of deformed shapes at time t = 0.25 s with M3 model: (a) URL-SPH algorithm, (b) URL-SPH algorithm without dissipation (i.e. by setting the values of D pχ ab = 0 and D j ab = 0)

Fig. 12
Fig. 12 Punch block: comparison of deformed shapes at time t = 0.065 s with M3 model: (a) URL-SPH algorithm, (b) URL-SPH algorithm without dissipation (i.e. by setting the values of D pχ ab = 0 and D j ab = 0)

Fig. 14
Fig. 14 Punch block: time evolution of (a) different energy measures, (b) global entropy, (c) temperature at position X = [0, 0, 0.5] T and (d) global numerical dissipation.A Mie-Grüneisen thermo-elastic model as described in Sect.3.3.1 is used.The corresponding material parameters are summarised in Table 3

Fig. 19
Fig.19 Taylor bar impact: comparison of deformed shapes at time t = 80 ms.The first three columns (M1, M2 and M3, from left to right) show the particle refinement of a structure simulated using the URL-SPH algorithm, whereas the last column shows a deformed structure via in-house Total Lagrangian SPH algorithm[1].The first row depicts

Fig. 23
Fig. 23 Taylor bar impact: a sequence of deformed structures with temperature profile at times t = {15, 30, 45, 60, 80} ms (from left to right) using (top row) the proposed URL-SPH algorithm and (bottom

Fig. 24
Fig. 24 Necking bar: comparison of deformed shapes at time (a) t = 20 ms and (b) t = 40 ms.The first three columns (left to right representing M1, M2 and M3) show the particle refinement of a structure simulated using the proposed URL-SPH algorithm, whereas the last column shows a deformed structure via the mixed-based Total Lagrangian SPH algorithm [1] (via M3).Colour contour indicates pressure field.A Hencky-based Johnson-Cook hardening rule as described in Sect.3.3.2 is used.The corresponding material parameters are summarised in Table 4

Fig. 25
Fig. 25 Necking bar: (a) time evolution of different energy measures and (b) radius reduction as function of elongation of the bar.A Hencky-based Johnson-Cook hardening rule as described in Sect.3.3.2 is used.The corresponding material parameters are summarised in Table 4

Fig. 26
Fig. 26 Necking bar: comparison of deformed shapes at time t = 60 ms via M3 model: (Left) URL-SPH algorithm without adding numerical dissipation term D j ab in the conservation equation for j (66), and (Right) URL-SPH algorithm incorporating sufficient numerical dissipation term D j ab .Zoom-in view of region near necking is included.A Hencky-based Johnson-Cook hardening rule as described in Sect.3.3.2 is used.The corresponding material parameters are summarised in Table 4

Algorithm 3 :
Updated Reference Lagrangian SPH AlgorithmInput : initial geometry X a and initial states of p a χ , f a , h a , j a , η a

χ
Output: current geometry x a , particle velocity v a and current states of F a , H a , J 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 χ , ḟa , ḣa , ja and ηa χ (14) ENSURE conservation of angular momentum (15) COMPUTE smoothed velocities using the corrected kernel (16) EVOLVE p a χ , f a , h a , j a , η a χ and x a (17) COMPUTE σ a

Table 1
Linear thermo-elasticity: material parameters used in the simulation

Table 2 L
-shaped block: material parameters used in the simulation

Table 3
Punch block: geometry and material parameters used in the simulation

Table 4
Taylor bar: material parameters used in the simulation tions.For completeness, three levels of particle refinement for the model are studied.Model M1 contains a number of 1740 SPH particles, model M2 comprises 4108 particles and model M3 contains 8160 particles.To accurately capture the onset of necking, it is important to place more particles within the necking region.