A novel Arbitrary Lagrangian Eulerian Smooth Particle Hydrodynamics algorithm for nonlinear solid dynamics

This paper introduces a novel Smooth Particle Hydrodynamics (SPH) computational framework that incorporates an Arbitrary Lagrangian Eulerian (ALE) formalism, expressed through a system of first-order conservation laws. In addition to the standard material and spatial configurations, an additional (fixed) referential configuration is introduced. The ALE conservative framework is established based on the fundamental conservation principles, including mass, linear momentum and the first law of thermodynamics represented through entropy density. A key contribution of this work lies in the evaluation of the physical deformation gradient tensor, which measures deformation from material to spatial configuration through a multiplicative decomposition into two auxiliary deformation gradient tensors. Both of the deformation tensors are obtained via additional first-order conservation equations. Interestingly, the new ALE conservative formulation will be shown to degenerate into alternative mixed systems of conservation laws for solid dynamics: particle-shifting, velocity-shifting and Eulerian formulations. The framework also considers path-and/or strain rate-dependent constitutive models, such as isothermal plasticity and thermo-visco-plasticity, by integrating evolution equations for internal state variables. Another contribution of this paper is the evaluation of ALE motion (known as smoothing procedure) by solving a conservation-type momentum equation. This procedure is indeed useful for maintaining a regular particle distribution and enhancing solution accuracy in regions characterised by large plastic flows. The hyperbolicity of the underlying system is ensured and accurate wave speed bounds in the context of ALE description are presented, crucial for ensuring the stability of explicit time integrators. For spatial discretisation, a Godunov-type SPH method is employed and adapted. To guarantee stability from the semi-discretisation standpoint, a carefully designed numerical stabilisation is introduced. The Lyapunov stability analysis is carried out by assessing the time rate of the Ballistic energy of the system, aiming to ensure non-negative entropy production. In order to ensure the global conservation of angular momentum, we employ a three-stage Runge–Kutta time integrator together with a discrete angular momentum projection algorithm. Finally, a range of three dimensional benchmark problems are examined to illustrate the robustness and applicability of the framework. The developed ALE SPH scheme outperforms the Total Lagrangian SPH framework, particularly excelling in capturing plasticity regimes with optimal computational efficiency.


A B S T R A C T
This paper introduces a novel Smooth Particle Hydrodynamics (SPH) computational framework that incorporates an Arbitrary Lagrangian Eulerian (ALE) formalism, expressed through a system of firstorder conservation laws.In addition to the standard material and spatial configurations, an additional (fixed) referential configuration is introduced.The ALE conservative framework is established based on the fundamental conservation principles, including mass, linear momentum and the first law of thermodynamics represented through entropy density.A key contribution of this work lies in the evaluation of the physical deformation gradient tensor, which measures deformation from material to spatial configuration through a multiplicative decomposition into two auxiliary deformation gradient tensors.Both of the deformation tensors are obtained via additional first-order conservation equations.Interestingly, the new ALE conservative formulation will be shown to degenerate into alternative mixed systems of conservation laws for solid dynamics: particle-shifting, velocity-shifting and Eulerian formulations.The framework also considers path-and/or strain rate-dependent constitutive models, such as isothermal plasticity and thermo-visco-plasticity, by integrating evolution equations for internal state variables.Another contribution of this paper is the evaluation of ALE motion (known as smoothing procedure) by solving a conservation-type momentum equation.This procedure is indeed useful for maintaining a regular particle distribution and enhancing solution accuracy in regions characterised by large plastic flows.The hyperbolicity of the underlying system is ensured and accurate wave speed bounds in the context of ALE description are presented, crucial for ensuring the stability of explicit time integrators.For spatial discretisation, a Godunov-type SPH method is employed and adapted.To guarantee stability from the semi-discretisation standpoint, a carefully designed numerical stabilisation is introduced.The Lyapunov stability analysis is carried out by assessing the time rate of the Ballistic energy of the system, aiming to ensure non-negative entropy production.In order to ensure the global conservation of angular momentum, we employ a three-stage Runge-Kutta time integrator together with a discrete angular momentum projection algorithm.Finally, a range of three dimensional benchmark problems are examined to illustrate the robustness and applicability of the framework.The developed ALE SPH scheme outperforms the Total Lagrangian SPH framework, particularly excelling in capturing plasticity regimes with optimal computational efficiency.

Introduction
In the field of computational solid dynamics, traditional methodologies predominantly employ the Lagrangian (Total or Updated) framework for continuum representation, as extensively documented in literature [1][2][3][4][5].These methodologies ensure alignment between the computational and material domains, enabling computational nodes and material particles (carrying physics) to move concurrently.This significantly simplifies problem formulation and the development of algorithm for discrete solutions.Naturally, Lagrangian formulations are effective in directly tracking of evolving boundaries and interfaces.However, in scenarios involving very large deformations, Lagrangian-based formulations may encounter excessive distortion within the material domain, potentially compromising the accuracy of the results and, in severe instances, jeopardising the stability of the overall algorithm.In contrast, some researchers advocate for Eulerian-based computational methods, where the computational domain is anchored to a fixed spatial domain.In this case, the material continuum is permitted to convect freely within this fixed spatial domain, thus accommodating large scale distortions.The convective nature of this approach requires specific strategies, such as employing the level set method, to accurately track the evolution of boundaries and interfaces.
Given the advantages and limitations inherent to both Lagrangian and Eulerian methods, Arbitrary Lagrangian Eulerian (ALE) based techniques have been developed to harness the strengths of each approach.The fundamental principle of the ALE formulation [6][7][8][9][10][11][12][13][14][15][16][17] involves utilising a referential (fixed) domain for motion description.This referential domain is different from the material domain (used in Lagrangian descriptions) and the spatial domain (used in Eulerian descriptions).For instance, in ALE methods, the computational nodes do not strictly adhere to the material particles nor remain fixed in space.Instead, the computational nodes dynamically adjusts to follow the material to a certain extent, being continuously updated to mitigate irregular distribution of material particles.This adaptability of the ALE framework makes it particularly suited for addressing challenges in large strain solid dynamics, such as those encountered in dynamic impact/contact scenarios and crack propagation.
The ALE approach for solid dynamics was primarily developed to address challenges in modelling materials that exhibit substantial hypoelastic-plastic and hyperelastic-plastic deformations.In hypoelastic-plastic models [18,19], stress responses are quantified as rates, establishing a relationship between an objective stress rate (such as the Jaumann stress rate [20], which corrects for rigid body rotations) and the rate of deformation tensor.In contrast, hyperelastic-plastic models [15,21,22] enable stress measures to be directly obtained from the deformation gradient, utilising a strain energy potential combined with a suitable plastic projection algorithm.No stress rate equation is required.Within both frameworks, the ALE convective term within the stress update equation accommodates the differential movement between computational nodes and material particles.A widely accepted approach to handle the convective term is to employ a sequential process by dividing the ALE solution process into two successive phases: a material phase and a convection phase (known as remap stage [23][24][25]).Convection is set aside in the material phase, rendering identical to the standard Lagrangian updating process.Subsequently, stresses (along with plastic internal variables) are then adjusted to reflect the relative movement between the computational nodes and material particles.
Recently, we introduced a computational framework for solid dynamics based on a novel first-order hyperbolic ALE formulation [26], with special emphasis on hyperelastic models.Specifically, the conservation equations for mass and linear momentum were re-formulated and solved with respect to the referential domain.To more accurately capture the deformation behaviours of solids, the physical deformation gradient tensor (measured from material to spatial domains) was obtained through a multiplicative decomposition into two auxiliary deformation gradient tensors, both computed using additional geometric conservation laws.For spatial discretisation, an upwinding-based vertex-centred finite volume method was employed, as detailed in Ref. [26].
In the current work, we extend the ALE formalism described above to include thermal effects by solving an additional total energy conservation equation expressed in terms of entropy density within the reference domain.A standout feature of our ALE formulation is its adaptability to simplify into three different mixed-based sets of conservation equations pertinent to solid dynamics: the Total Lagrangian [27][28][29][30][31], Eulerian [32] and Updated Reference Lagrangian formulations [33].This adaptability can be seen as a generalisation of the velocity-shifting and particle-shifting frameworks previously established in the context of fluid dynamics.Additionally, our work also explores hyperelastic and path-dependent constitutive models [34][35][36][37] under both isothermal and thermal conditions, where the latter requires the incorporation of evolution equations for internal plastic state variables.Drawing inspiration from the pioneering work [22], we introduce a simple ALE motion strategy by solving a conservation-type momentum equation for the ALE motion.A critical aspect requiring special attention is the stability of the ALE formulation from both the continuum and numerical perspectives.To ensure continuum stability, we employ an internal energy potential that satisfies rank-one convexity [26], which guarantees the existence of plane travelling waves within the solid and the accurate propagation of strong discontinuities at physical shock speeds.More detailed analysis can be found in Ref. [26].For numerical stability, we adopt a Godunovtype Smooth Particle Hydrodynamics (SPH) method, enhanced with appropriate numerical stabilisation techniques to ensure local numerical entropy production.This aspect of stability is evaluated through Lyapunov stability analysis, by monitoring the time variation of the Ballistic energy within the coupled system discretised using the SPH method.Similar analyses have also been carried out in the context of finite volume method [38].Another contribution of the paper is to ensure the global conservation of angular momentum.This is achieved by introducing an explicit three-stage Runge-Kutta time integrator, combined with a discrete angular momentum projection algorithm.
The paper is organised as follows.Section 2 summarises the fundamental kinematic relationships in an Arbitrary Lagrangian Eulerian description.Section 3 details the conservation laws within the ALE framework, addressing both smooth and discontinuous solutions.Section 4 explores the potential to degenerate into the Total Lagrangian velocity-shifting formulation, the Updated Lagrangian particle-shifting formulation and the Eulerian formulation.Section 5 extends the discussion to incorporate models of irrecoverable plasticity.Section 6 presents the weak form of the equations and reintroduces the second law of thermodynamics in terms of the Ballistic energy within the system.The time variation of the Ballistic energy serves as a foundation for investigating non-negative entropy-producing numerical dissipation.Section 7 discusses the ALE motion procedure employed in this work to mitigate excessive particle distortions, especially in areas experiencing large plastic flows.In Section 8, we present the Godunov-type SPH discretisation procedure, with a focus on deriving entropy-stable, conservative and consistent numerical stabilisations.Furthermore, we discuss the explicit time integrator utilised for advancing the semi-discrete equations in time.Section 9 introduces a space-time discrete angular momentum projection algorithm.Section 10 provides a summary of the algorithmic flowchart for the ALE SPH framework.Section 11 presents a series of numerical examples to evaluate the capabilities of the numerical method.Section 12 summarises key conclusions and outlines prospective lines of research.Appendix A summarises the transformation rule from Total Lagrangian description to the ALE counterpart.

ALE kinematical description
Consider the motion of a continuum initially defined by a domain Ω  ⊂ ℝ 3 in its material configuration.This configuration is enclosed by a boundary, denoted as Ω  , with a corresponding unit outward normal vector   = ∑ 3 =1      .As the continuum undergoes motion, it deforms into a spatial configuration defined by a domain Ω  ⊂ ℝ 3 with a boundary Ω  and a unit outward normal vector   = ∑ 3 =1      .In an Arbitrary Lagrangian Eulerian (ALE) description, in addition to the usual material and spatial configurations described above, a third "referential" configuration is introduced.This configuration is occupied by the domain Ω  ⊂ ℝ 3 with a boundary Ω  and an outward unit normal   = ∑ 3 =1      .As a result, two new mappings are introduced, linking the referential coordinate  ∈ Ω  to: (1) the material coordinate  ∈ Ω  via  = (, ) and (2) the spatial coordinate  ∈ Ω  via  = (, ).Fig. 1 illustrates these domains and the one-to-one transformation relating the three configurations.We refer to (, ) as the spatial motion (the mapping from the referential domain to the spatial domain) and to (, ) as the material motion (the mapping from the referential domain to the material domain).With these definitions, the physical motion (, ) at a fixed time is parametrised as  = • −1 , emphasising that the three mappings {, , } are not independent.Utilising the relation (, ) = ( −1 (, ), ), the multiplicative decomposition of the physical deformation gradient tensor is expressed as with the respective deformation gradient tensors defined as It is possible to introduce the corresponding cofactors and Jacobians as follows and Here, denotes the cross product between second order tensor as discussed in Refs.[39,40].Eqs. ( 2), ( 3) and ( 4) establish line, area and volume mappings between differential elements of the different configurations as follows We now turn our attention to the time evolution of the mappings {, , }, as illustrated in Fig. 1.Specifically, we explore the velocities associated with these mappings: the (physical) velocity  of a particle , the (spatial) velocity v linked to the spatial mapping, and the (material) velocity  corresponding to the material mapping.These velocities are mathematically defined as follows where the notation The above relation illustrates that the spatial velocity is the sum of the physical velocity and the product of the deformation gradient with the material velocity.This relationship allows us to express the relative velocity , the difference between the physical velocity  and the spatial velocity v, as Notably, a relative velocity  =  implies that the referential domain coincides with the material domain, a condition that leads to the Total Lagrangian formulation.Conversely, when  = , indicating that  = − −1 , this suggests that the referential domain aligns with the spatial domain.Such alignment corresponds to the Eulerian formulation.Understanding this distinction is crucial for choosing the most suitable framework to model the behaviour of deformable solids under different conditions.

First-order conservation laws for thermal processes
In the prior work of the authors [26,38], the motion of a deformable solid during a thermal-mechanical process is described through a system of first-order conservation laws formulated in an ALE description.The governing equations are formulated as The shorthand notation [•] denotes the time differentiation holding the material coordinate  fixed, that is . This notation will be used throughout the entire manuscript.Here,   is the material density,   =    is the linear momentum per unit of material volume,   denotes the body force per unit of material volume and   is the entropy density per unit of material volume.Additionally,  represents the first Piola Kirchhoff stress tensor,  denotes the absolute temperature,   represents the heat flux defined in the reference domain and Phy represents the material time rate of physical dissipation introduced by a possible irreversible constitutive model, such as visco-elasticity or plasticity.It is important to note that this term is zero in the context of a reversible elastic model [26].The operator DIV  represents the divergence operator carried out at the reference configuration and  is the second order identity tensor.Additionally, it is crucial to ensure the satisfaction of two (initial value) differential constraints, also referred to as involutions [41]), as stated below If these constraints are satisfied initially, expressions (9b) and (9c) indicate these constraints will continue to be satisfied throughout the entire time integration process.Finally, the coordinates, including {, , }, can be updated over time using expressions (6).
In the system of equations above (9), Eq. (9d) corresponds to the classical mass continuity conservation law, whereas Eq. (9a) refers to the standard linear momentum conservation law.To address thermal effects and the presence of an irreversible constitutive model, an additional conservation Eq. (9e), denoted as the first law of thermodynamics, is included.To accurately measure the deformation of the solid from the reference configuration to both spatial and material configurations, we introduce a supplementary pair of geometric conservation laws, namely (9b) and (9c).The inclusion of these geometric conservation equations relaxes the evaluation of the fibre mappings {  ,   } from the exact gradient evaluation of both material and spatial geometries.This approach has already been shown to be advantageous in the case of Total and Updated Lagrangian descriptions, particularly in scenarios involving large deformations, as demonstrated in Refs.[33,37,[42][43][44].
In the presence of shocks or discontinuities, it is important to note that evaluating the divergence operators DIV  in system (9) becomes unfeasible.This implies that the applicability of the system (9) in its differential form is no longer valid.To address this, we must resort to the use of appropriate Rankine-Hugoniot jump conditions across a discontinuous surface.This surface is defined by a reference unit normal vector   and propagates with speed   in the reference space [26].These conditions can be expressed as Here, the symbol [[•]] represents the difference between the right and left states across a discontinuous surface.Observe that, when both the heat flux   =  and the ALE motion velocity is absent ( = ), expression (11e) simplifies to its Total Lagrangian counterpart, namely [[  ]] ≤ 0. Intuitively, this implies that the entropy behind the shock must be greater than ahead, indicating that the shock creates entropy as it propagates [45].
Remark 1.Interestingly, if we consider the material density   to be uniform1 , the mass conservation expression (9d) reduces to a geometric conservation equation for the material volume map   .The corresponding local differential equation and its jump condition are thus formulated as

Combined system for solid dynamics
Combining the equations presented above results into a system of first-order hyperbolic conservation equations written in an ALE setting (with respect to the referential domain).This can be expressed as with vector of conserved variables   , flux vector    and source term   described by where, for the sake of compactness, we define In addition to the initial and boundary (essential and natural) conditions necessary to fully define the initial boundary value problem, the closure of the system (13) requires the introduction of an appropriate constitutive law that must adhere to various physical criteria, such as thermodynamic consistency (verified through the Coleman-Noll procedure [46,47]) and compliance with the principle of objectivity [4].In the context of thermo-elasticity, this involves introducing a suitable internal energy potential, denoted as ( ,   ).This potential ensures that stress  , deformation gradient tensor  and entropy   are related through the relation  = ( ,  ) . Ensuring stability requires the internal energy potential  to demonstrate rank-one convexity (also known as ellipticity) in the sense of the Legendre-Hadamard condition [48,49].An extension to consider path-dependent constitutive models will be explored in the next section.
As for the heat flux vector   , we employ the Fourier's law for an isotropic material written in a reference description as where ℎ denotes the non-negative thermal conductivity coefficient calibrated in the spatial configuration.
Remark 2. Utilising established Calorimetry relationships [45,46]  , we obtain a simple one-to-one expression relating temperature and entropy [45] Θ( ,   ) =   exp where   is a reference temperature, η ( ) is an entropy function dependent upon the Jacobian  evaluated at the reference temperature, and   =     , where   and   represent the specific heat per unit mass and per unit material volume, respectively, which have been assumed to be constant.The symbols Θ and  are used interchangeably to denote the same temperature with different functional dependencies.
Remark 3. The evaluation of the hyperbolic nature of the system ( 13) is crucial to ensuring the existence of real-valued wave speeds across all possible deformation states.As discussed in Refs.[26,50], the condition of a (rank-one) convex internal energy functional is sufficient to guarantee hyperbolicity for a given material motion.For a more in-depth understanding, we refer the reader to Ref. [26] for a detailed explanation.
The relation between the referential wave speed   and the material wave speed   is succinctly summarised in the equation below

Alternative versions of the first-order conservation laws in solid dynamics
The application of suitable kinematic considerations can transform the ALE system of conservation laws (13) into alternative continuum conservation frameworks.In previous work [26], the authors have shown how the ALE system can degenerate into: (1) a Total Lagrangian framework, (2) an Updated Reference Lagrangian framework, and (3) an Eulerian framework.Through this process, the ALE system (13) was shown to be an elegant generalisation of its alternative continuum counterparts.This Section explores further this concept in order to explain and motivate recent velocity-shifting and particle-shifting techniques which have gained important traction within the SPH community.We will first consider in this section the simplified case of thermo-elasticity by assuming that the dissipation of the selected constitutive model is zero, that is Phy = 0, before considering a more advanced inelastic model in the following section.

Total Lagrangian formulation with velocity perturbation
The ALE system of Eqs. ( 9) can be specialised to the case when considering a mapping (, ) from the referential to the material domain such that its deformation gradient can be considered negligible (i.e.{  ,   ,   } ≈ {, , 1}) but not necessarily its velocity field  ≠ .This leads to the following particularised system of equations Notice that in this case, the geometric conservation law (9c) is unnecessary due to its redundancy and does not feature in (19).Moreover, the deformation from the referential to the spatial domain reduces to   ≈  (refer to Eq. (19b)) with the velocity  = − −1 .Finally, notice that the differential operator DIV  can be approximated as DIV  ≈ DIV  .This set of equations extends the velocity-shifting continuum formulation [51][52][53] for free surface fluid flows as originally derived in the context of SPH.Finally, by setting the perturbation velocity  ≈ , the corresponding Total Lagrangian system [37,42,43,[54][55][56][57] can be recovered, as  ≈  and the material domain is then completely aligned with the referential domain.

Updated reference Lagrangian formulation
In this specific scenario, we explore a situation where a negligible velocity  induces non-negligible material deformation   = F  .From a mathematical standpoint, this situation can be defined by the conditions {  ,   ,   } = { F  , H , J } and  ≈ , resulting into v ≈ .When these conditions are applied, the ALE system (9) degenerates to the following set of equations Notice that in this case, the geometric conservation law (9b) is unnecessary due to its redundancy, and the deformation gradient tensor from the referential to the spatial domain can be recovered as  ≈   F −1  following solution of (20b).Finally, the time derivative term can be approximated as This set of equations can be re-interpreted as the Updated Reference Lagrangian formulation as proposed by [33,37].In these works, the deformation triad { F  , H , J } was carefully chosen to continuously track the movement of the deformable solid.Moreover, this approach can also be viewed as an adaptation of the particle-shifting continuum formulation [58][59][60][61], initially developed for fluid dynamics in the context of SPH.

Eulerian formulation
The derivation of the corresponding Eulerian system is facilitated by aligning the referential and spatial domains through the enforcement of specific kinematic conditions.By setting v =  and imposing {  ,   ,   } = {, , 1}, we ensure the coincidence of the referential and spatial coordinates  = , the material velocity  = − −1  and the material mappings {  ,   ,   } = { −1 ,  −1 ,  −1 }.This set of conditions results in the simplification of the ALE formulation to its Eulerian counterpart, as discussed in Refs.[32,45,62].To facilitate this transition, we introduce the following fields: The Eulerian system is thus characterised by the following conservation equations This reformulation captures the essence of the Eulerian perspective by focusing on spatial rather than material coordinates, where div  denotes the divergence operator in the spatial configuration.

Extension to irrecoverable plasticity model
To model the strain-rate and thermally dependent behaviour, we employ a von Mises plasticity model utilising the nonlinear Johnson-Cook hardening law [37].This inelastic constitutive model is summarised in this section, with an in-depth description included in Ref. [38].The Johnson-Cook nonlinear hardening law is characterised by a yield function  , which depends on the deviatoric Kirchhoff stress  ′ and a yield stress τ .The yield stress is a function of the equivalent plastic strain ε , its plastic strain rate ̇ε  (involving the material time derivative) and the temperature .This relationship is expressed as Here, τ represents the generalised scalar von Mises equivalent stress.The nonlinear hardening law is described by where () is defined as In the above expression,  denotes the current temperature,  melting is the melting temperature of the material and  transition is the temperature at or below which there is temperature dependence of the yield stress.Furthermore, τ0  denotes the initial yield stress and ε0 is the reference strain rate calibrated based on the material.The remaining material constants include the material hardening parameter , hardening exponent , strain-rate coefficient  and temperature exposure .Notably, when  >  melting , the material is assumed to melt and behave like a fluid, offering no shear resistance since τ = 0. Interestingly, by setting the values  = 1,  = 0 and  transition ≈ ∞ (which in turn implies () = 0), the rate-independent linear hardening law can be simply recovered, that is To establish a physically meaningful flow rule, the concept of work conjugacy [4] is applied.In the context of elasto-plasticity, the material time derivative of plastic dissipation Phy (per unit material volume) is formulated in terms of the Kirchhoff stress  and the plastic rate of deformation   as [4] To ensure the principle of maximum plastic dissipation, it is crucial to establish a correlation between the plastic rate of deformation (first equality of   ) and the flow rule (second equality of   ), ensuring the material efficiently dissipates energy due to plastic deformation.Here, ̇, a proportionality factor, is referred to as the plastic multiplier (or consistency parameter).It is expressed as a rate to maintain dimensional consistency with the plastic rate of deformation   , and is obtained enforcing the yield constraint.By rearranging the above expression, and recalling the ALE transformation rule (A.4), the ALE time evolution expression of the inelastic right Cauchy-Green tensor can be obtained as The subsequent step is to establish the relationship between the consistency parameter ̇ and the evolution of the internal hardening variable ε .By recalling the von Mises equivalent stress τ from (23), the material time rate of ε is defined to be the work conjugate to τ, described as below Comparing the first term and the last term of the above equation, and again utilising the transformation rule (A.4), the ALE time evolution equation for the equivalent plastic strain is 6. Weak form and second law of thermodynamics

The conjugate variables of the underlying system
To convey a meaningful physical interpretation to the conjugate fields of interest in this work, let us consider the Ballistic energy   per unit of referential volume defined by where   =     .Both   (, ) and  (  ,   ,   ,   ,   ,  −1  ) represent alternative functional expressions of the same magnitude.In this equation, the terms on the right-hand side have different physical meanings: the first term denotes kinetic energy, the second term represents internal energy and the third term represents thermal energy.

C.H. Lee et al.
With appropriate energy conjugate fields for the two deformation measures {  ,   } defined as along with the energy conjugates to the entropy density   and to the inverse of the inelastic right Cauchy-Green strain tensor  −1 the associated work conjugates   become Note that the expression for   −1

𝑝
is obtained specifically for isotropic materials, where the dependency is expressed through   ( ,  −1  ).A detailed discussion for this expression can be found in Ref. [38].Here,  =  −   represents the temperature difference between the current temperature and the reference temperature, and the (thermal , defined as a combination of the kinetic energy () = 1  2   ( ⋅ ), the internal energy  and the heat dissipative energy   .

Variational statement
In general, a standard weak variational statement for the thermo-mechanical coupled system is established by multiplying the local differential Eqs. ( 13) (expressed in terms of the conservation variables   ) with their corresponding work conjugate virtual fields   (34).This multiplication is followed by integration over the reference domain Ω  of the body, resulting in the formulation It is interesting to notice that, under an isothermal deformation process ( =   and then  = 0), the entropy density evolution Eq. (9e) separates from the rest of the system when treated in a weak manner.This implies a degeneration of the fields (34) into the usual work conjugates utilised for isothermal plasticity, as detailed in Ref. [63].
Performing integration by parts on the last term of the right-hand side in (35), and expanding the resulting equation, yields where the normal fluxes    =       with    representing the outward normal in th referential direction.The representation (36) can be particularised to specific conservation variables, namely 2 Employing the Legendre transformation enables the Lagrangian function to be expressed in terms of Helmholtz's free energy potential , namely (,  , ) = ()−( , ).
Furthermore, it is crucial to emphasise that integrating by parts, as shown in (36), facilitates the imposition of physical boundary conditions.This is particularly useful for the linear momentum Eq. (37a) and the total energy Eq.(37e), as these expressions naturally introduce the boundary traction   , the boundary temperature   and the heat flux   .In contrast, the integration-by-parts technique is less critical for the remaining geometric conservation laws.Here, velocity boundary conditions, such as   and   , are imposed directly through essential boundary conditions.This direct approach differs from the natural enforcement achieved through surface integrals in the Eqs.(37a) and (37e).Additionally, it is essential to demonstrate that the weak form ((37a) to (37e)) ensures the satisfaction of fundamental thermodynamic principles, such as the second law.This is achieved through the addition of expressions (37a)-(37e) together with the conjugate fields specifically chosen as   =   .Subsequent algebraic manipulation leads to the following inequality Here, Πext  denotes the mechanical power associated with external forces as and  ext  accounts for the heat source and heat flux added (removed) to (from) the system as Inequality (38) represents a valid expression for the second law of thermodynamics for the system.A detailed derivation is provided in Ref. [38].
The fulfilment of inequality ( 38) is a necessary ab initio condition to ensure stability, commonly known as the classical Coleman-Noll procedure.This concept will be further explored in Section 8.2 when introducing entropy-stable numerical dissipation.

ALE material motion
In general, the primary objective of the ALE smoothing procedure is to serve as a flexible method for appropriately adjusting the position of the underlying discretisation (i.e.particles).This adjustment aims to mitigate common issues in large strain SPH simulations, such as particles clustering too closely or stretching too far apart in the spaces they occupy, especially in regions experiencing significant deformation (i.e.plastic flow).The ALE approach seeks to improve the robustness of the formulation and the accuracy of the solution fields, whilst also potentially reducing computational costs.As previously explored in Ref. [38], the material (shifting) motion of particles   can be defined by solving a conservation-type law, along with an appropriate jump condition, described as follows where   =    denotes the momentum per unit of reference volume.Additionally,   represents the first Piola-Kirchhoff stress associated with the corresponding material motion.Above conservation equation monitors how the material velocity  evolves as particles move from the referential domain to the material domain.The remaining unknown in the above expression is the definition of the stress tensor, which measures material force per unit area acting on a surface in the reference configuration.This can be achieved by introducing the scalar potential function Π through a convex combination of strain energy potential functions defined in terms of   and F  , where F  = F ( )  , with function F ( ) specifically designed to selectively incorporate certain aspects of the deformation pertinent to the problem being considered.Additionally, to prevent particles inter-penetration, it is useful to further incorporate a spring strain type of energy  spr as We refer to the numerical examples presented in Section 11 that demonstrate this concept.In the equation above, the first two terms on the right-hand side depend on the deviatoric component of a constitutive model Ψ and the last term denotes the spring model, defined as where  denotes a Young's modulus type of constant, the non-dimensional parameters  and  range from 0 to 1, and  denotes the number of dimensions of the problem.Differentiating expression (42) in time for a given deformation gradient  provides which eventually gives Remark 4. In this work, we adopt the deviatoric strain energy Ψ from the neo-Hookean model, defined as follows Here, μ =   represents a (fictitious) shear modulus defined in terms of the physical shear modulus .The dimensionless coefficient , varying within the range of [0, 1], adjusts  to reflect a fraction of the physical shear modulus.This adjustment is particularly relevant in scenarios dominated by large inertial effects, where a lower  value can prevent the fast propagation of waves and potential shock formation during ALE motion.

Smooth particle hydrodynamics spatial discretisation
Applying the standard gradient correction proposed in [64] to ensure zeroth-and first-order completeness, the SPH discretisation for the first-order system described in (13), after some algebraic manipulation, takes the following form The particle-based internal force    is defined as where the pseudo-area vectors [44] are Furthermore, the particle-based material force    , consisting of two contributions, is given by It is important to note that we derive the final term of the above equation (that is, the spring-related internal force) through a two-step process.First, we employ the skew-symmetric standard kernel gradient, which relies on the reference geometry.This means that the gradient of   at the location   is the negative of the gradient of   at   , formally expressed as     (  ) = −    (  ).Second, we assume that among the stretch values represented by   , one is significantly dominant over the others and this predominance is aligned with the direction that connects two particles, denoted by Physically, the first term on the right-hand side of equation ((50)b) depends on the deviatoric part of the constitutive model.This is used to describe how particles move and deform, focusing on changes in shape without volume change.The second term represents the restoring force preventing particles from clustering too closely together during the smoothing process.This balance is crucial for simulating realistic particle behaviour, ensuring that particles maintain an appropriate distance from each other, thus avoiding clumping.
Eqs. (47a), (47d) and (47e) correspond to the conservation of linear momentum, mass and total energy, respectively.Complementing these physical principles, Eqs.(47b) and (47c) introduce an additional set of geometric conservation equations.Lastly, Eq. (47f) denotes the ALE material motion, providing a full representation of the physical and geometric aspects within the system.
For the case of thermo-elasticity, the motion of solids is sufficiently described by the expressions from (47a) to (47d), along with the condition  Phy = 0.However, incorporating plasticity within the framework of the first-order system requires additional evolution expressions for the associated plastic state variables, expressed in (28) and (30).The corresponding semi-discrete equations at particle  are respectively.In the above expressions, the gradient at particle  for the inverse of the plastic right Cauchy-Green deformation tensor  −1  and the equivalent plastic strain ε can be approximated through the standard corrected kernel gradient evaluation, described as In the context of plasticity, the rate of plastic multiplier ̇ is determined by enforcing the yield constraint, ensuring the material response adheres to its defined yield surface.
For completeness, it is beneficial to monitor the evolution of the geometry (i.e.spatial distribution of the particles) to examine the effectiveness of the ALE material motion.In addition, geometry computation is unavoidably necessary if wishing to enforce accurate global conservation of the angular momentum (see Section 9).With that in mind, the spatial velocity v is integrated over time, and so are the material velocity  and the physical velocity  making use of As it is well known, the above mixed-based system ( 47) is prone to exhibit non-physical numerical instabilities [9,34,35,[65][66][67][68][69]] when attempting to model large strain solid dynamics problems [1,70,71].To address this issue, a Godunov-type stabilisation term is suitably added to both the pair-wise internal force    and the pair-wise material force    as The mathematical expressions for the stabilisation terms {

𝒑 𝑾
} are then obtained through the semi-discrete version of the classical Coleman-Noll procedure [72], which will be the focus of the following section.

Numerical entropy production
Based on the derivation of the second law of thermodynamics as presented in Section 6.2, we examine the semi-discrete form of inequality (38) under the assumption that   is known or has been prescribed.This stability analysis is demonstrated by multiplying expressions (47a)-(47e) with the corresponding conjugate fields (34)  ] . ( Note that the last term in the second line of (55), referred to as the rate of the plastic dissipation introduced by the plasticity model, can be expressed in various alternative ways as Next, we substitute the linear momentum Eq. (47a), the geometric conservation equation for both spatial and material deformation gradient tensors (47b) and (47c), the mass conservation Eq. (47d) into the third line of (55).Through algebraic rearrangements, we arrive at
Additionally, by considering the referential gradient of the Lagrangian function expressed as , and applying the relationship enables further manipulation of the first term within the square brackets on the fourth line of (57).This leads to Subsequently, expression (57) becomes ) . ( Notice that the first term on the right-hand side corresponds to the ALE convective term of the entropy density conservation Eq. (47e) and the second term corresponds to the ALE convective term of the evolution equation for the inverse of the inelastic right Cauchy-Green tensor (51a).When these expressions combine with the respective evolution Eqs.(47e) and (51a), the convective terms cancel out.After careful algebraic manipulation, expression (60) reduces to with the physical dissipation term Phy, featuring in the second term on the right-hand side of the above equation.Here, Πext and  ext are the semi-discrete versions of power contribution and total heat contribution, expressed as To adhere to the second law of thermodynamics (38) at the semi-discrete level, it is crucial to demonstrate all three terms on the right-hand side of (61) are non-negative.Regarding the heat conduction term (the first term on the right-hand side of (61)), the inequality is naturally satisfied due to the inherent properties of the conductive heat flux.Detailed derivation can be found in Ref. [37].As for the physical model dissipation term (the second term on the right-hand side of ( 61)), its non-positivity arises from the definition of the rate of physical plastic dissipation.Our next objective is to demonstrate the non-negativity of the numerical dissipation term, denoted as  total ≥ 0. This can be achieved by equivalently swapping indices  and  to give By simply averaging the first term and the second term of the expression above, and noting the anti-symmetric nature of the stabilisation term as  , an alternative expression for  total is shown below Sufficient conditions for  total ≥ 0 are given by where     is a positive semi-definite stabilisation matrix.Following previous work of the authors [26], the dissipation term is chosen as where . Furthermore, the dissipation term is related to the jump in velocity between interacting nodes, a characteristic feature of Godunov-type upwinding terms [63,[73][74][75].
Remark 5. Given that, in general, the ALE mesh motion is not known a priori but it will be defined via the first-order conservative-type Eq. (47f), the introduction of an appropriate stabilisation term is necessary to stabilise the convective nature inherent to the equation.In this context, and for simplicity, we propose a stabilisation term     with a mathematical structure similar to that of expression (65) as In this context, the positive-definite stabilisation tensor is in the form of

Time integration
Given the size of the semi-discrete equations, we choose a three stage Runge-Kutta explicit time integrator [33]  )) . ( Moreover, all the three geometries are also integrated using the same time integrator, resulting in a monolithic time update procedure.The unknowns   = (   ,   ,   ,   ,   )  evolve alongside the geometries (such as ,  and ) through Eq. ( 69).Additionally, the determination of the maximum time step, denoted as Δ =  +1 −   , adheres to the standard Courant-Friedrichs-Lewy (CFL) [76] criterion.This criterion is defined by the following relation Here,    represents the CFL stability number and    denotes the pressure wave speed measured in referential domain and ℎ min stands for the minimum (or characteristic) length within the computational domain.The relationship between the referential pressure wave speed    and the material pressure wave speed    is presented in expression (18).Unless specified otherwise, a value of    = 0.9 has been chosen in the subsequent examples to ensure a balance between accuracy and stability.

Plasticity
In addressing the irrecoverable elasto-plastic model presented in this paper, we utilise the implicit backward Euler time integrator as our chosen method for advancing the semi-discrete equations governing the plastic state variables, specifically  −1   and ε , as presented in ( 51)).This choice of time integration method is consistent with the traditional return mapping algorithm employed in plasticity models [4,77], thus ensuring numerical stability.The equations resulting from the time integration processes are presented below As previously discussed, the above gradient computations are approximated using the corrected kernel gradient.These approximations are expressed as For the Johnson-Cook hardening model, determining the incremental plastic multiplier Δ  requires the enforcement of the yield condition as described in (23).This enforcement typically involves the solution of nonlinear equations which would then require an iterative Newton-Raphson procedure.A more detailed discussion can be found in Ref. [37].
To solve Eqs.(71a) and (71b), we implement a predictor-corrector algorithm.In the predictor step (Lagrangian phase), the algorithm advances explicitly excluding the ALE convective effect, yielding intermediate variables { −1,int  , εint  }.In the corrector step (remapping phase), the intermediate variables are then projected to their new positions after applying the ALE smoothing procedure.The mathematical formulation for each step is detailed below, starting with the predictor step for a time increment Δ as Following this, the corrector step becomes

Discrete angular momentum preserving algorithm
The ALE algorithm presented above does not inherently ensure exactly the conservation of angular momentum, as the strain measures {  ,   } (or  ) are no longer exclusively computed from the referential gradient of both material and spatial geometries (or current geometry).Although the deviation from strict conservation is minimal, we introduce a monolithic discrete angular momentum projection algorithm for completeness.In our numerical experiments, we have observed that activation of this algorithm is not strictly necessary unless seeking "machineaccurate" strict conservation of angular momentum.Specifically, the local internal nodal force    is suitably modified in a least-square sense to preserve global angular momentum, whilst still ensuring the global conservation of linear momentum.
Building on the procedure introduced by [33], we explicitly enforce sufficient conditions for the global preservation of discrete linear and angular momentum within a time step at each of the three stages, as described below A least-square minimisation procedure is applied to obtain a modified set of internal nodal forces, denoted as T   , satisfying the specified conditions.This involves minimising the following functional (ignoring time arguments for brevity) We consider the stationary condition of the functional Π with respect to { ang ,  lin } and T   separately.The derivative of Π with respect to { ang ,  lin } leads to the constraints shown in (75).Additionally, the derivative of the functional in (76) with respect to T   yields a modified set of internal nodal forces T   , expressed as Utilising expressions (75) (with    replaced with T   ) in conjunction with (77), allows the determination of the Lagrange multipliers { ang ,  lin } as where the indicial notation

Algorithmic description
Algorithm 1 provides a summary of the ALE SPH algorithm for nonlinear solid dynamics in irreversible processes.A key advantage of the ALE approach is its capability to prevent irregular distribution of particles by dynamically adjusting their positions, especially in areas experiencing significant plastic deformation.This benefit is exploited in the example section, demonstrating its superior efficiency in comparison to traditional Lagrangian approaches.

Numerical examples
In this section, several benchmark problems are examined in order to evaluate the performance and applicability of the proposed ALE algorithm.The purpose of these examples is to demonstrate the efficiency of the ALE algorithm across various scenarios.Specifically, the ALE algorithm • ensures compliance with the Geometric Conservation Law (GCL) condition [9,19,[80][81][82], even under the scenario of fictitiously moving SPH particles, • obtains consistent second-order convergence for velocities/displacements, stresses/strains and temperature/entropy, • ensures the global conservation of angular momentum [83,84] through the incorporation of a discrete angular momentum projection algorithm, • maintains a non-negative rate of production of both numerical (i.e.Godunov-type stabilisation) and physical (e.g.plastic dissipation and irreversible heat conduction process) entropy within an insulated system, and • offers advantages in handling large plastic deformations when compared to conventional Lagrangian approaches, especially in high-speed impact and fast stretching scenarios.
In subsequent plasticity examples, we explore the ALE smoothing process by solving the conservation-type equation, as indicated in (41).However, in the context of isothermal-and thermo-elastic problems, this technique tends to results in minimal particle movement.To facilitate more significant particle movement in these instances, we instead prescribe a sinusoidal-type motion for the material particles, specifically where χ1 =  1 +   1 and χ2 =  2 +   2 denote coordinates shifted relative to the centroid of the domain   .The parameter  dictates the magnitude of particle movement, and  represents the period of the sinusoidal functions.By adopting this sinusoidal motion strategy, we achieve more pronounced and controlled particle movement.This is indeed useful for examining the robustness of the ALE algorithm in both isothermaland thermo-elastic scenarios, as it allows for an assessment of the performance of the algorithm through varying degrees of particle movement.Differentiating expression (79) with respect to time yields the material particle velocity  , expressed as Both material particle motion (79) and velocity (80) are carefully designed to ensure movement is constrained tangentially along the material boundary to prevent volume change.
To account for large thermo-elastic deformations, we utilise a Mie-Grüneisen thermo-elastic model.The corresponding internal energy density , as a function of deformation gradient  and entropy density   (per unit of material volume), is given by C.H. Lee et al.
Here, the internal energy density at reference temperature is formulated as The deviatoric ′  ( −1∕3  ) and the volumetric  ( ) contributions are where  and  denote the shear modulus and bulk modulus of the material, and Γ 0 is a positive material constant.In addition, we also incorporate a simple volumetric-based Mie-Grüneisen model [45,50,85] described in [45] by where  is a dimensionless parameter varying from zero (i.e. a perfect gas) to one (i.e.solid materials).Previous work [50] has demonstrated the polyconvexity (and thus, rank-one convexity) of this coupled model.Finally, for comparison purposes, we benchmark the ALE results against our alternative in-house, mixed-based Total Lagrangian SPH Framework (TLF) [37,44], which incorporates Godunov-type entropy-stable numerical stabilisation.This Godunov-type TLF has proven to be very effective in simulating large strain solid dynamics.For completeness, we also compare some of the results with the classical displacementbased Total Lagrangian SPH Framework (CTLF) [86].It is important to note that the primary unknown variable in CTLF is displacement, and this framework does not include any numerical dissipation.Remark 6.It is important to note that the coupled energy functional described in Eq. ( 81) can be reduced to the classical isothermal neo-Hookean model by setting the coupling constant Γ 0 = 0 and maintaining the temperature (, ) =   .This implies that   (, ) = 0, due to the direct relationship between temperature and entropy as described in (17).For further clarity, the strain energy for the neo-Hookean model can be expressed as follows

Satisfaction of geometric conservation laws: Patch test
This example involves a patch test that demonstrates the efficacy of the proposed ALE algorithm in accurately capturing rigid body translation, similar to the scenario detailed in [26].Consider a column with a square cross section of 1 m 2 and a length of  = 6 m, initially prescribed by a Fig. 3. Translation with heat.Time evolution of the  2 global errors in the velocity magnitude ‖ v − v0 ‖  2 and in temperature ‖ −   ‖  2 when using (red) the complete ALE system, (blue) a reduced system without solving variable   and (black) further reduced system without solving {  ,   }.A discretisation of 625 number of particles.A Mie-Grüneisen model is used with parameters listed in Table 1.uniform velocity field  0 = [3, 4, 0]  ms −1 throughout the entire domain.The boundaries of this structure are left unconstrained.Fig. 2(a) depicts the setup for this problem.The main aim of this example is to evaluate the ability of the algorithm to accurately simulate the behaviour of rigid body translation through arbitrary particle movement, while strictly adhering to the discrete Geometric Conservation Law (GCL) conditions.In addition, a Mie-Grüneisen model is incorporated to account for the effects of thermal variations.The material properties and the particle movement parameters for this setup are detailed in Table 1.The ALE movement over time is governed by the expression in (80).2.
Since the column is subject to rigid body translation, it is anticipated to traverse through space in the direction of the initial velocity, maintaining a constant speed and a uniform temperature distribution.For a qualitative analysis, Fig. 3 presents the temporal progression of global errors.These errors include the deviation of the current temperature  from the reference temperature   , and the differences between the current and initial velocity components.As expected, when employing the complete ALE algorithm as presented in (47), these errors fluctuate around machine precision levels.In contrast, a simplified model that omits the solution for   results in significantly larger errors, failing to meet the discrete GCL conditions associated with the volumetric strain.This discrepancy is visually evident in Fig. 4. For instance, the reduced model yields non-uniform contour profile, unlike the smooth contours achieved with the complete ALE system.Moreover, a further simplified model that neglects both   and   solutions inaccurately tracks the deformation path, leading to the eventual breakdown of the numerical algorithm.

Proof of order of convergence: manufactured problem with smooth exact solutions
As detailed in Refs.[37,50], we utilise a cubic domain with unit-length sides to assess the order of convergence of the proposed ALE framework.This example employs a specific physical mapping defined below ) sin ( with parameters  =  =  = 1 and  0 = 5 × 10 −4 m.These settings enable the exact computation of the velocity field , deformation gradient tensor  , and temperature field , as described below and ) cos ) cos Dirichlet boundary conditions are applied at the boundaries of the domain, with symmetry on the faces  1 = 0,  2 = 0 and  3 = 0 and skew-symmetric on the faces  1 = 1,  2 = 1 and  3 = 1.The simulation parameters are summarised in Table 2.
C.H. Lee et al.Fig. 5 displays the  2 norm convergence analysis at time  = 1 × 10 −3 s for the velocity components , the components of the first Piola Kirchhoff stress tensor  , and the temperature .The results demonstrate second-order convergence in space and time for all these variables, a further enhancement over traditional ALE methods [22], which typically exhibit a reduced order of convergence for the derived variables such as stresses and strains.

Satisfaction of conservation of linear and angular momenta
A column with a square cross-section and a length  of 6 m is left rotating in space with an initial velocity field described as where  0 = Ω 0 [0, 0, 1] rad/s represents the angular velocity vector with Ω 0 = 105 rad/s.This problem setup is depicted in Fig. 2(b).The objective of this example is to illustrate the conservation capabilities of the proposed ALE algorithm.The column is made of a nearly incompressible rubber material, with its material properties and parameters for particle motion summarised in Table 3.In this case, three different levels of particle refinement are analysed: {M1, M2, M3} comprising {625, 1715, 3969} particles, respectively.First, a particle refinement study simulated using the proposed ALE-SPH algorithm is investigated, as depicted in the first three rows (from top to bottom) of Fig. 6.Both the deformation pattern and pressure distribution progressively converge with increased particle refinement.For benchmarking, the performance of the proposed algorithm is also compared with that of an alternative in-house, mixed-based Godunov-type Total Lagrangian SPH framework (TLF) [44] using M3 discretisation, as shown in the last row of Fig. 6.The comparison, yielding nearly identical results, is illustrated in Fig. 6.Second, Figs.7(a) and 7(b) highlight the ability of the proposed algorithm in preserving global linear and angular momenta.The global linear momentum,  = ∫ Ω    Ω  , fluctuates close to zero value (within machine error) at all times.The global angular momentum, denoted as  = ∫ Ω   ×   Ω  , maintains its initial value throughout the simulation.Minor reduction in angular momentum is observed when the discrete angular momentum projection algorithm is omitted.Additionally, the evolution of kinetic energy , elastic energy  and the total energy , defined as the sum of  + , are illustrated.Third, and for a qualitative assessment, Fig. 7(d) monitors the time evolution of the horizontal displacement component at the position  = [0.5, 3, 0.5]  , demonstrating convergence with increased refinement.Finally, Fig. 8 showcases a series of deformed states, with the pressure distribution indicated by the colour contour plot.We observe stable solutions even after a long-term response.

Robustness of the proposed formulation in highly nonlinear scenario
This section revisits a well-documented example [1,45,86], focussing on a column subjected to an initial sinusoidal velocity field about the origin, given by where Ω 0 = 105 rad s −1 indicates the initial angular velocity magnitude (see Fig. 2(c)).A neo-Hookean material model is used.The relevant material properties and parameters related to the prescribed ALE particle motion are summarised in Table 4. Initially, the efficacy of the ALE-SPH algorithm is assessed through a particle refinement study.See the first three columns from the left in Fig. 9.The results demonstrate that structural deformations predicted using a very coarse discretisation with (M1) 625 particles are in remarkable agreement with the results from more refined discretisations, such as M2 with 1715 particles and M3 with 3969 particles.For comparison purposes, the performance of the  3.
algorithm is also evaluated against two benchmarks: a mixed-based Total Lagrangian Godunov-type SPH algorithm [44] and the recently proposed ALE vertex-centred finite volume algorithm [26].The deformed states from these algorithms are depicted in the fourth and fifth columns from the left of Fig. 9, with results that are almost identical to those obtained with the ALE-SPH algorithm.Furthermore, the robustness of the algorithm is tested across four varying magnitudes of the ALE motion, specifically  = {0, 1, 1.No discrepancies are noted when comparing the outcomes of the proposed algorithm to those of the established Total Lagrangian Godunov-type SPH algorithm [44].As depicted in Fig. 12, the ALE-SPH algorithm introduces numerical dissipation over the duration of the simulation, ensuring the discrete satisfaction of the second law of thermodynamics.

Assessment on the reliability of the algorithm under impact event
In this example, we examine a benchmark scenario involving a copper bar, initially with a length of  = 0.0324 m and a radius of  = 0.0032 m, colliding against a rigid wall at a speed of  0 = 227 m∕s.The setup for this problem is illustrated in Fig. 13(a).Utilising a Hencky-based von Mises material model along with a rate-independent linear hardening law (26), the specific parameters for this simulation are included in Table 5.To enhance computational efficiency, the simulation simplifies the scenario to consider only a quarter of the entire domain, implementing    4.  4.
as particles clumping.To prevent the excessive compression of particles near the contact zone, incorporating ALE particle motion proves advantageous.Through ALE smoothing, the particles in these critical areas can be stretched vertically upwards, thus maintaining a high quality of the spatial particles distribution.This is facilitated by solving a conservation-type law, as described in (41), with the function F specifically chosen as A series of images in Fig. 14 demonstrate the evolving state of ALE particle movement throughout the impact process.To evaluate the reliability of our ALE algorithm, a quarter section of the bar is discretised using three levels of particle refinement, namely (M1) 1560, (M2) 3744 and (M3) 7280 particles.As illustrated in Fig. 15, the deformation and plastic strain distribution patterns at time   = 80 s show remarkable consistency across all levels of particle refinement.The simulations in the left three columns are based on the ALE-SPH algorithm with increasing levels of refinement, whilst the rightmost column presents results obtained through the mixed-based Total Lagrangian Godunov-type SPH Framework (TLF) [43] Moreover, Fig. 16 presents a comparison between the outcomes achieved with the ALE algorithm (based on M3) and those obtained using the mixed-based Total Lagrangian SPH counterpart [44].The top row of this figure features snapshots that reveal the distribution of equivalent plastic strain, while the bottom row focusses on the pressure profile.Each comparison is divided into two parts: the left side represents the ALE results, and the right side displays the findings from the mixed-based Total Lagrangian approach.Despite yielding comparable results, the ALE method excels in maintaining a more regular distribution of particles around the impact zone.This advantage is evident in the detailed views provided within Figs. 17, highlighting the superiority of ALE algorithm at handling particles clumping during high-impact events.This is visible in the front view of the figure, where denser particle clustering can be observed in the right portion of the bar (representing the mixed-based Total Lagrangian approach) near the contact region, compared to the left portion of the bar (ALE-SPH algorithm).For completeness, we have also simulated the same scenario using classical displacement-based SPH approach.The comparison, illustrated in Fig. 18 at time  = 23 s, shows that the accumulation of particle clumping leads to the breakdown of the numerical scheme.Consequently, the classical displacement-based Total Lagrangian approach was unable to complete the simulation.Fig. 19(a) displays the dynamic energy transfer during the impact process, showing the conversion of a significant portion of kinetic energy into plastic dissipation and a smaller fraction into elastic strain energy.This indicates the dynamic exchange of energy that takes place during the collision.Moreover, Fig. 19(b) shows that the numerical dissipation of the system gradually decreases with particle refinement.This aligns with the second law of thermodynamics, thus affirming the physical validity of the simulation.The evolution of the radius at the position  = [0.0032,0, 0.0324]  is monitored and depicted in Fig. 19(c).Our results show excellent consistency with results from the Total Lagrangian approach described in [44].A key benefit of the ALE method is its capability to prevent particles from clumping, a critical consideration in simulations of high-speed impacts.This capability allows for the use of larger time increments in the ALE method compared to those in the Total Lagrangian method, enhancing both computational efficiency and the robustness of the algorithm under challenging scenarios.In this specific case, the ALE time steps are larger than those of the Total Lagrangian by approximately 2-2.5 times.This is shown in Fig. 19(d) where the time increment ratio is defined as the ratio between ALE and Total Lagrangian time increments.

Accuracy in the necking plastic region
Revising the initial conditions introduced in the Taylor bar problem, we now consider an altered condition where the initial velocity field of the bar is inverted, causing it to stretch on both sides.As shown in Fig. 13(b), this adjustment primarily aims to assess the effectiveness of the proposed ALE methodology in capturing large plastic flow in the vicinity of the necking region.Classical displacement-based Lagrangian approaches often prove inadequate in accurately resolving plastic deformations because particles in the necking area become excessively separated [37].A viable strategy to address this challenge involves implementing non-uniform refinement by increasing the number of particles within the necking region to better capture its deformation behaviour.
For this simulation, we employ a Hencky-based von-Mises plasticity model along with a nonlinear Johnson-Cook hardening law.The specific material properties and simulation parameters are detailed in Table 6.Considering thermal effects, we set the initial temperature profile of the bar to (,  = 0) = 573.15K.The ALE motion in this context is governed by the previously utilised function F through solving expression (41).A series of snapshots illustrating ALE particle motions, particularly emphasising the necking area, is presented in Fig. 20.Due to significant stretching, the ALE motion tends to compress, thereby ensuring that the spatial distribution of particles remain regular and of high quality -a key factor for accurately capturing the behaviour of the material under extreme conditions.
To begin with, we perform a particle refinement study utilising the proposed ALE formulation.This investigation involved three levels of particle refinement, specifically (M1) 1560, (M2) 3744 and (M3) 7280 particles.The results for both von Mises stress and temperature demonstrate convergence with increased particle refinement, as depicted in the first three columns from the left in Fig. 21.Compared to the mixed-based Total Lagrangian method, as shown in the rightmost column of Fig. 21, the ALE approach offers a notable improvement in temperature resolution.This improvement is attributed to a higher concentration of particles in the necking region.In contrast, particles in this area tend to become highly irregular under the Total Lagrangian approach, which adversely affects the accuracy of the simulation.The proposed ALE approach predicts a larger radius due to temperature softening, indicating a more precise representation of temperature effects.

Simulation of challenging coining process
The study of a coining test, depicted in Fig. 13(c), focusses on a rectangular workpiece subjected to compression by a rigid tool.Due the symmetrical nature of the problem, only a quarter of the domain, measuring 0.03 m × 0.01 m × 0.000625 m, is modelled for computational efficiency.The workpiece undergoes compression to 60% of its original height, driven by a boundary velocity of  0 =  0 [0, −1, 0]  m∕s applied  (26).The corresponding simulation parameters are tabulated in Table 5.  (92) The setup reveals a critical limitation of the Total Lagrangian finite element based algorithm in capturing the dynamic evolution of the free surface of workpiece due to severe element distortion.This study aims to highlight the enhanced performance and reliability of the proposed ALE algorithm.A hyperelastic-plastic materials combined with a linear hardening rule is utilised.The detailed material properties and simulation parameters relevant to this analysis are provided in Table 7.No thermal effects are considered.
A key challenge in the coining process is the ability of the algorithm to handle outward flow of material from the punch region [10,87].This phenomenon complicates the accurate representation of the evolving free surface, especially when a limited number of particles are available for this purpose.The proposed ALE algorithm aims to overcome this limitation by facilitating dynamic particle flow towards the punch region.This ALE approach results in a higher concentration of particles in areas of critical importance, thus improving the distribution of plastic strain, as demonstrated in Fig. 24.Compared to the Total Lagrangian approach, the ALE algorithm demonstrates superior potential in simulating the coining process with greater fidelity.This can be observed in Fig. 25.For a quantitative comparison, we monitor the time history plot of horizontal and vertical displacements of a particle at position  = [0.0125,0.01, 0.0003125]  m, as shown in Fig. 26.When the punch occurs, the particle initially bulges upward, exhibits outward flow, and then moves towards the punch region.The proposed ALE algorithm appears to successfully capture this behaviour.Unfortunately, the Godunov-type Total Lagrangian SPH algorithm cannot effectively model the actual physical behaviour of material flowing inward towards the punch region.This test case paves the way for future research into its application across complex manufacturing processes, such as cold-forming.

Conclusions
In this paper, a novel Arbitrary Lagrangian Eulerian (ALE) Smooth Particle Hydrodynamics (SPH) framework re-formulated through a system of first-order conservation laws was introduced.Physical conservation principles, including mass, linear momentum and entropy density, were supplemented with two auxiliary geometric conservation laws to facilitate the resolution of the physical deformation gradient tensor  .By exploiting the multiplicative decomposition of  into the pair {  ,   } and solving them through conservation laws, our methodology permitted accurate evaluation of strains and stresses.The hyperbolic nature of the system and the physical propagation of shocks were ensured through the rank-one convexity of the internal energy potential, which establishes the one-to-one relationship between real wave speeds in material and referential domains.Interestingly, our ALE framework can be degenerated into three simplified continuum mechanics representations for solids: Total Lagrangian, Eulerian and Updated Reference Lagrangian formulations.Both the Total and Updated Reference Lagrangian formulations can be viewed as extensions of the velocity-shifting and particle-shifting approaches, which are traditionally employed in fluid dynamics, to the context of solid dynamics.To handle irregular particle distribution, a dynamic ALE motion strategy has been employed by solving a conservationbased (material) linear momentum equation.Furthermore, for materials with history-dependent characteristic, our framework integrated additional evolution equations for internal state variables within the ALE setting.
From spatial discretisation standpoint, a Godunov-type SPH algorithm was used.Appropriate numerical stabilisation was carefully derived to fulfil the Coleman-Noll procedure at the semi-discrete level via the Ballistic energy, ensuring a consistent non-negative numerical entropy production.From temporal integration standpoint, a three-stage Runge-Kutta explicit time integrator complemented by an angular momentum projection algorithm was implemented to preserve angular momentum.The determination of wave speed limits helped in computing the time step size necessary for the stability of the time integrator.
Through extensive numerical examples, we thoroughly assessed the conservation properties, stability, and robustness of the proposed algorithm.The ALE SPH framework demonstrated significant potential in enhancing the resolution of pressure and stress fields within critical plastic regions of substantial deformation challenges, such as necking.Furthermore, an additional advantage of the ALE approach included the

Fig. 2 .
Fig. 2. Problem setup for the (a) translation test case, (b) rotating test case and (c) twisting column.

Fig. 4 .
Fig. 4. Translation with heat.A sequence of deformed states at time  = {0, 0.5, 1, 1.5} s (from left to right for each row) when using (a) the complete ALE, comparing against the reduced system without considering (b)   and (c)   and   .Colour contour indicates the absolute temperature error, denoted as abs( −   ).A discretisation of 625 number of particles.A Mie-Grüneisen model is used with parameters listed in Table 1.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4 Fig. 5 .
Fig. 5. Manufactured problem:  2 global convergence analysis at time  = 1 × 10 −3 s for (a) the components of velocity, (b) the components of the first Piola-Kirchhoff stress tensor, and (c) the temperature.A Mie-Grüneisen model is used with parameters listed in Table2.

Fig. 6 .Fig. 7 .
Fig. 6.Rotating column.Comparison of the deformed shapes at time (left)  = 0.05 s and (right)  = 0.1 s.The first three rows (top to bottom) show the model refinement of a structure simulated using the proposed ALE algorithm, whereas the last row shows a deformed structure via an alternative mixed-based Godunov-type Total Lagrangian SPH algorithm [44].Colour indicates the pressure profile.A neo-Hookean model is used.The relevant simulation parameters are summarised in Table 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 5, 2} × 10 −2 m.The simulations, presented in Fig. 10, show consistent results across all magnitudes.No instabilities are observed in any of the simulations, emphasising the importance of the tailor-made numerical stabilisation term.A qualitative assessment, shown in Fig. 11, monitors the evolution of the twist angle at two different positions, namely  = [0.5, 6, 0.5]  and  = [−0.5, 6, −0.5].

Fig. 8 .
Fig. 8. Rotating column.A series of snapshots at times  = {0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165} ms (from left to right and top to bottom).Colour indicates the pressure distribution.Results obtained using a neo-Hookean model with M3 model.The relevant simulation parameters are summarised in Table 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .Fig. 10 .Fig. 11 .
Fig. 9. Twisting column.Comparison of the deformed shapes at time  = 100 ms.The first three columns (from left to right)show the model refinement of a structure simulated using the proposed ALE algorithm.In contrast, the fourth and fifth columns illustrate the deformation achieved using an alternative in-house Godunov-type Total Lagrangian SPH algorithm[44] and an ALE vertex-centred finite volume algorithm[26].Colour indicates the pressure profile.A neo-Hookean model is used.The simulation parameters are summarised in Table4.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Twisting column.Time evolution of (a) different energy measures (e.g.kinetic energy , internal energy  and the sum of kinetic and internal energies, known as the total energy ), (b) total numerical dissipation (%) introduced by the proposed ALE method via M3 for  = {0, 0.01, 0.015, 0.02}.(c) and (d) show the time history of the numerical dissipation and its percentage for three levels of refinement, compared to the mixed-based Godunov-type Total Lagrangian SPH algorithm [44] (TLF).A neo-Hookean model is used.Parameters related to the material properties and the mesh motion are tabulated in Table4.

Fig. 13 . 14 .
Fig. 13.Problem setup for the (a) Taylor bar impact test case, (b) rapid stretching of a bar and (c) coining test case.

Fig. 15 .
Fig. 15.Taylor bar impact.Comparison of the deformed shapes at time  = 80 s.The first three columns (left to right)show the model refinement of a structure simulated using the proposed ALE algorithm, whereas the last column (at the rightmost position) shows a deformed structure via an alternative Godunov-type Total Lagrangian SPH algorithm[44].Colour indicates both von Mises (left portion of the subfigure) and equivalent plastic strain (right portion of the subfigure) profiles.A Hencky-based von Mises plasticity model using isotropic linear hardening rule is used(26).The corresponding simulation parameters are summarised in Table5.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .Fig. 17 .
Fig. 16.Taylor bar impact.A series of snapshots at times  = {20, 40, 60, 80} s (from left to right).Colour indicates (first row) the equivalent plastic strain and (second row) the pressure distribution.Each subfigure is split into two parts: the left side represents the results of ALE, and the right side illustrates those of the Total Lagrangian approach [44].Results obtained using a Hencky-based von Mises plasticity model combined with a linear hardening rule.The corresponding simulation parameters are summarised in Table 5. (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.Taylor bar impact.Comparison of the deformed shapes at time  = 23 s via M2 model.In each subfigure, left portion of the bar represents the results obtained through the proposed ALE algorithm and the right portion of bar using the Classical Total Lagrangian SPH Framework (CTLF) [86].Colour indicates pressure profile.Dense particle clustering is observed when using CTLF.A Hencky-based von Mises plasticity model using linear hardening rule is used (26).The corresponding simulation parameters are summarised in Table 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 22
compares the deformation of the bar at time  = 30 ms.As expected, the ALE method (left portion of the figure) provides superior resolution over the Total Lagrangian method (right portion of the figure) in capturing temperature near the necking areas.Furthermore, the evolution of various energy components, including Ballistic energy, aligns well between both formulations.As demonstrated in Fig.23(a), the reduction in Ballistic energy throughout the simulation process adheres to the principles of the second law of thermodynamics, thereby affirming the long-term stability of the simulation outcomes.Further exploration, as seen in Fig.23(b), monitors the development of the necking radius.

Fig. 19 .
Fig. 19.Taylor bar impact.Time evolution of (a) different energy measures (e.g.kinetic energy , internal energy  and the sum of kinetic and internal energies, known as the total energy ), (b) the accumulated numerical dissipation measured in percentage, (c) radius at the position  = [0.0032,0, 0.0324]  and (d) a ratio between the ALE and Total Lagrangian time increments.Results obtained using a Hencky-based von Mises plasticity model with isotropic linear hardening rule (26).The corresponding simulation parameters are tabulated in Table5.
to a third of the top side with  0 = 10 m∕s, denoted as  0

Fig. 20 .
Fig. 20.Rapid stretching of a bar.Two-dimensional view of a series of material mesh deformation at times  = {8, 16, 24, 30} s (from left to right).Break lines are used to indicate that a part has been shortened for representation on the figure.

Fig. 21 .
Fig. 21.Rapid stretching of a bar.Comparison of the deformed shapes at time  = 30 s.The first three columns (left to right) show the refinement of a structure simulated using the proposed ALE algorithm, whereas the last column (at the rightmost position) shows a deformed structure via an alternative Godunov-type Total Lagrangian SPH algorithm [44].Colour indicates both von Mises (left portion of the subfigure) and temperature (right portion of the subfigure) profiles.A Hencky-based von Mises plasticity model combined with a nonlinear Johnson-Cook hardening rule is used.The parameters related to material properties and ALE motion are detailed in Table 6.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 22 .Fig. 24 .
Fig. 22. Rapid stretching of a bar.Comparison of the deformed shapes at time  = 30 s via M3 model: (Left portion of the bar) the proposed ALE algorithm and (Right portion of the subfigure) Godunov-type Total Lagrangian SPH Framework (TLF) [44].Colour indicates the temperature profile.A Hencky-based von Mises plasticity model combined with a nonlinear Johnson-Cook hardening rule is used.The parameters related to material properties and ALE motion are detailed in Table 6.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 25 .
Fig. 25.Coining problem.Comparison of the deformed shapes at times  = {300, 600} s (from top to bottom) (Left portion of the block) the proposed ALE algorithm and (Right portion of the block) Godunov-type Total Lagrangian SPH Framework (TLF) [44].Colour indicates both von Mises profile.A Hencky-based von Mises plasticity model combined with a linear hardening rule is used.The corresponding simulation parameters are summarised in Table 7. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) that interconnect the internal energy density , entropy density   and temperature , and taking into account the relation  = and summing them to obtain in this work.The time integration from time step   to  +1 is governed by the following equations

Table 1
Translation with heat: material properties and parameters of the prescribed ALE particle motion.

Table 2
Manufactured problem: material properties and parameters of ALE motion.

Table 3
Rotating column: material properties and parameters of the prescribed ALE particle motion.

Table 5
Simulation parameters for isothermal Taylor bar impact.boundary conditions such as (frictionless) roller support to accomplish this.The primary objective is to demonstrate the superiority of the proposed ALE algorithm over Classical (displacement-based) Total Lagrangian SPH Framework (CTLF), particularly in situations involving substantial plastic flow in the vicinity of the contact area.CTLF approaches often suffer from significant irregular particle distributions, such symmetry

Table 6
Simulation parameters for rapid stretching of a bar.

Table 7
Simulation parameters for coining problem.