A first-order hyperbolic arbitrary Lagrangian Eulerian conservation formulation for non-linear solid dynamics

The paper introduces a computational framework using a novel Arbitrary Lagrangian Eulerian (ALE) formalism in the form of a system of first‐order conservation laws. In addition to the usual material and spatial configurations, an additional referential (intrinsic) configuration is introduced in order to disassociate material particles from mesh positions. Using isothermal hyperelasticity as a starting point, mass, linear momentum and total energy conservation equations are written and solved with respect to the reference configuration. In addition, with the purpose of guaranteeing equal order of convergence of strains/stresses and velocities/displacements, the computation of the standard deformation gradient tensor (measured from material to spatial configuration) is obtained via its multiplicative decomposition into two auxiliary deformation gradient tensors, both computed via additional first‐order conservation laws. Crucially, the new ALE conservative formulation will be shown to degenerate elegantly into alternative mixed systems of conservation laws such as Total Lagrangian, Eulerian and Updated Reference Lagrangian. Hyperbolicity of the system of conservation laws will be shown and the accurate wave speed bounds will be presented, the latter critical to ensure stability of explicit time integrators. For spatial discretisation, a vertex‐based Finite Volume method is employed and suitably adapted. To guarantee stability from both the continuum and the semi‐discretisation standpoints, an appropriate numerical interface flux (by means of the Rankine–Hugoniot jump conditions) is carefully designed and presented. Stability is demonstrated via the use of the time variation of the Hamiltonian of the system, seeking to ensure the positive production of numerical entropy. A range of three dimensional benchmark problems will be presented in order to demonstrate the robustness and reliability of the framework. Examples will be restricted to the case of isothermal reversible elasticity to demonstrate the potential of the new formulation.

In addition, with the purpose of guaranteeing equal order of convergence of strains/stresses and velocities/displacements, the computation of the standard deformation gradient tensor (measured from material to spatial configuration) is obtained via its multiplicative decomposition into two auxiliary deformation gradient tensors, both computed via additional first-order conservation laws.Crucially, the new ALE conservative formulation will be shown to degenerate elegantly into alternative mixed systems of conservation laws such as Total Lagrangian, Eulerian and Updated Reference Lagrangian.Hyperbolicity of the system of conservation laws will be shown and the accurate wave speed bounds will be presented, the latter critical to ensure stability of explicit time integrators.For spatial discretisation, a vertex-based Finite Volume method is employed and suitably adapted.To guarantee stability from both the continuum and the semi-discretisation standpoints, an appropriate numerical interface flux (by means of the Rankine-Hugoniot jump conditions) is carefully designed and presented.Stability is demonstrated via the use of the time variation of the Hamiltonian of the system, seeking to ensure the positive production of numerical entropy.A range of three dimensional benchmark problems will be presented in order to demonstrate the robustness and reliability of the framework.Examples will be restricted to the case of isothermal reversible elasticity to demonstrate the potential of the new formulation.

INTRODUCTION
Classical formulations used in computational solid dynamics are typically based on purely Lagrangian (Total or Updated) descriptions of the continuum, [1][2][3][4][5] in which the computational mesh (i.e.nodes) coincides with the underlying material discretisation (i.e., particles).In this description, particles and nodes move concurrently, which greatly simplifies both the formulation of the problem and the resulting discrete solution algorithm.Naturally, Lagrangian formulations enable direct tracking of moving boundaries and interfaces.However, for some problems involving very large deformations, Lagrangian type formulations can suffer from excessive mesh distortion which can inevitably lead to inaccurate results and, in extreme cases, may also lead to the breakdown of the overall algorithm due to mesh inversion.As a result, the use of remeshing strategies 6 becomes unavoidable if the simulation is to successfully progress to completion.The design of remeshing approaches can sometimes be non-trivial and render smoothing projection errors in the unknown fields.On the other hand, some authors have favoured the use of Eulerian based computational methods, in which the computational mesh remains fixed in space and the background material continuum is allowed to freely convect through the fixed mesh, facilitating thus the treatment of very large distortions.As a result of this convective effect, specific treatment is necessary in order to follow moving boundaries and interfaces.Motivated by the above pros and cons, Arbitrary Lagrangian Eulerian (ALE) based techniques emerge in order to combine the advantages of purely Lagrangian and Eulerian approaches.The basic idea of the ALE formulation [7][8][9][10][11][12][13][14][15] is the use of a referential (fixed) domain for the description of the motion, different from the material domain (Lagrangian description) and the spatial domain (Eulerian description).Specifically, the computational mesh is neither attached to the material nor kept fixed in space.Indeed, the computational mesh partially follows the material points and is continuously evolved with respect to the material in order to reduce elements distortion.Thus, the use of an ALE description of the continuum has the potential to cope with large strain solid dynamics problems concerning, for instance, dynamic impact/contact and crack propagation.
The ALE formulation in the context of fast solid dynamics was originally conceived for highly deformable hypoelastic-plastic and hyperelastic-plastic processes.In hypoelastic-plastic models, 16,17 stresses are described in rate form relating an objective stress rate with the rate of deformation tensor.On the contrary, in hyperelastic-plastic models, [18][19][20] stresses are computed from the deformation gradient via a strain energy potential combined with a suitable plastic projection algorithm (i.e.no rate equation for stresses is then needed).In any of the above two approaches, a key issue is the ALE convective term featuring in the stress update equation, which reflects the relative motion between the computational mesh and the material particles.One popular approach to handle the convective term is the use of a sequential strategy whereby splitting the ALE solution process into two consecutive stages, namely, a material stage and a convection stage (also known as remap stage [21][22][23] ).Convection is neglected in the material stage, rendering identical to the standard Lagrangian updating process.Subsequently, stresses (and plastic internal variables) are then convected to account for the relative mesh motion in the remap stage.
In the context of hyperelasticity 20,24,25 (for instance, no plastic strains involved), Reference 20 developed an ALE finite element method without the need to consider convective terms.This can be achieved by suitably linking the referential, material and spatial domains with appropriate mappings, and so the deformation gradient tensor describing the fibre map between spatial and material domains (and thus, stresses) can be computed from the referential domain by applying the chain rule.Even in this case, if large plastic strains are considered, the convection of the internal plastic variables describing the plastic response is still needed. 19n the current work, a new first-order hyperbolic ALE computational framework for solid dynamics is presented, with emphasis on a hyperelastic model.Specifically, mass, linear momentum and the total energy conservation equations are re-formulated and solved with respect to the referential domain.In order to guarantee equal order of convergence of strains/stresses and velocities/displacements, the update of the physical deformation gradient tensor (measured from material to spatial domains) is obtained via its multiplicative decomposition into two auxiliary deformation gradient tensors, both computed via additional geometric conservation laws.One attractive feature of the presented ALE formulation is its possibility to degenerate into three different mixed-based sets of conservation equations, namely Total Lagrangian formulation, [26][27][28][29] Eulerian formulation 30 and Updated Reference Lagrangian formulation. 31This paper sets the foundations of a robust ALE computational framework capable of handling path-dependent constitutive models [32][33][34] (such as large strain elasto-visco-plasticity 35 ), along with the consideration of thermal effects 36 induced through strong thermo-mechanical coupling due to shocks.This will be the focus of our next publication on the topic.
An important aspect that requires special attention is that of the stability of the proposed ALE formulation, from both the continuum and numerical standpoints.With respect to the former, the use of a strain energy potential satisfying rank-one convexity guarantees as a consequence both the existence of plane travelling waves within the solid and the propagation of strong discontinuities at physical shock speeds.As for the spatial discretisation, an acoustic Riemann based (upwinding) approach is pursued where appropriate numerical stabilisation is introduced, ensuring local production of numerical entropy.The latter is shown by monitoring the time variation of the Hamiltonian of the system.
The paper is broken-down as follows.Section 2 introduces the general description of conservation laws in a Total Lagrangian description, both at integral and at local level, allowing for both smooth and discontinuous solutions.Section 3 introduces the fundamental kinematic relationships in an Arbitrary Lagrangian Eulerian description.Section 4 presents the transformation of conservation laws from the Total Lagrangian to the ALE description, detailing integral, smooth local and jump conditions.Section 5 presents the in-depth derivation of all the conservation laws used in this work for the simulation of solid dynamics, including both physical and geometric conservation laws.Section 6 demonstrates the hyperbolicity of the resulting system of conservation laws.Section 7 studies the time variation of the Hamiltonian of the system, used as a preamble to formulate non-negative entropy producing numerical dissipation.Section 8 details the spatial discretisation technique pursued in this paper, that is, a vertex centred finite volume method, paying special attention to the derivation of the conservative and consistent numerical fluxes.Section 9 presents the explicit time integrator used to advance the semi-discrete equations in time.Section 10 summarises the algorithmic flowchart of the resulting computational framework.Section 11 presents a series of numerical examples to assess the capabilities of the numerical method.Section 12 summarises some key conclusions and lines of prospective research.

TOTAL LAGRANGIAN CONSERVATION LAWS DESCRIPTION
Consider the motion of a continuum which in its material or initial configuration is defined by a domain Ω X ⊂ R 3 of boundary Ω X with unit outward normal vector N X = ∑ 3 I=1 N I X E I .After the motion, the continuum occupies a spatial (or current) configuration defined by a domain Ω x ⊂ R 3 of boundary Ω x with unit outward normal vector n x = ∑ 3 i=1 n i x e i .The motion is described by a time-dependent mapping (X, t) which links a material particle X ∈ Ω X to the spatial configuration x ∈ Ω x according to x = (X, t).As shown in previous works by the authors, 26,28,29,37 the motion of a deformable solid can be described by a system of first-order conservation laws formulated in a Total Lagrangian description as follows where  X denotes a vector of unknown conservation variables,  X represents the corresponding flux vector and  X is the vector of source terms, the symbol d dt represents the material time derivative and dΩ X and dA X = dA X N X represent a suitable differential volume element and outward area vector, respectively.For the case of sufficiently smooth solutions, where spatial derivatives of fluxes are well-defined, the above integral expression (1) is equivalent to a set of first-order differential equations described as where  t | | |X denotes the total or material time derivative (understood by holding X fixed) and DIV X is the divergence operator defined with respect to the material configuration, namely, DIV X ( X ) = ∑ 3 . In addition, for discontinuous solutions, the integral conservation laws (1) also lead to an appropriate set of jump conditions (known as Rankine-Hugoniot conditions) across a discontinuity surface of outward unit normal N X moving with wave speed U X [38][39][40][41][42][43][44] as In an Arbitrary Lagrangian Eulerian (ALE) description, in addition to the material and spatial configurations as presented in the previous section, a third "referential" configuration is introduced, occupied by the domain Ω  ⊂ R 3 of boundary Ω  with outward unit normal N  = ∑ 3 I=1 N I   I .As a result, two new mappings are introduced linking the referential coordinate  ∈ Ω  to: (1) the material configuration X ∈ Ω X via X = (, t) and (2) the spatial configuration x ∈ Ω x via x = (, t). Figure 1 shows these domains and the one-to-one transformation relating the three configurations.We refer to (, t) as the spatial motion (i.e. the mapping from the referential domain to the spatial domain), and refer to (, t) as the material motion (i.e., the mapping from the referential domain to the material domain).With these definitions in place, the physical motion (X, t) at a fixed time is parametrised as  = • −1 , emphasising that the three mappings {, , } are not independent.Utilising the relation (X, t) = ( −1 (, t), t), the multiplicative decomposition of the physical deformation gradient tensor is expressed as with the respective deformation gradient tensors defined as Here, F  denotes the material deformation gradient tensor measured from reference to material configuration, whereas F  is the spatial deformation gradient tensor measured from reference to spatial configuration.It is possible to introduce the corresponding Jacobians and cofactors as follows and where ⨯ denotes a cross product between second order tensor as discussed in References 45,46.Formulae ( 5), ( 6) and ( 7) allow to establish line, area and volume mappings between differential elements of the different configurations as follows Next, let us consider the time rates of the above mappings {, , } depicted in Figure 1.For instance, the (physical) velocity v of a particle X, the (spatial) velocity v associated with the spatial map, and the (material) velocity W associated to the material map are given by

THE ALE TRANSFORMATION
The starting point for deriving the ALE conservation form of the fundamental physical principle is to establish the relation between the usual Total Lagrangian conservation laws (expressed in the material domain) and the ALE conservation laws counterpart (expressed in the referential domain).In order to achieve this, recall first the global Total Lagrangian conservation laws as presented in (1).Re-parametrisation of X = (, t), together with the use of the standard Jacobian transformation dΩ X = J  dΩ  and Nanson's rule 4 mathematically defined as dA X = H  dA  , enables both the flux term and the source term on the right hand side of the above expression to be re-expressed as respectively.Insofar as the material domain is time-dependent (due to its re-parametrisation with respect to the referential domain), the application of Reynolds transport theorem (B1) to the material time derivative on the left hand side of (1) is necessary, which results in With the substitution of ( 11) and ( 12) into (1), the equivalent ALE conservation equations over a fixed referential domain become which leads to the corresponding set of local conservation equations where DIV  denotes the divergence operator defined with respect to the referential configuration.In addition, a suitable set of jump conditions across any moving shock in direction N  and speed U  can also be introduced as The above transformation will be applied to Lagrangian conservation variables such as mass, linear momentum and energy in the following section.
Remark 1. Equations ( 13), ( 14) and ( 15) can be re-expressed in terms of a set of reference conservation variables   and fluxes   as 47 t and where the reference conservation variables, fluxes and source terms are

ALE FIRST ORDER CONSERVATION EQUATIONS FOR SOLIDS
Utilising the transformation expressions described in (13), ( 14) and (15), both the global and local ALE conservation equations for mass, linear momentum and total energy, along with their appropriate Rankine-Hugoniot conditions, can be readily obtained.

Conservation of mass
In a Total Lagrangian setting, the global version of the classical mass conservation equation states that the mass of the domain Ω X is constant in time, that is where  X is the material density of the continuum (i.e.mass per unit material volume).With the use of the transformation Equation ( 13), the global ALE mass conservation equation therefore emerges as Using the Gauss divergence theorem on the boundary integral of (21), it is possible to obtain the differential (strong) form counterpart together with its jump condition described as respectively.Remark 2. Interestingly, if we consider the material density  X to be uniform*, the above mass conservation expression (21) reduces to a geometric conservation equation for the Jacobian map J  (also known as material volume map) presented as Again, application of Gauss divergence theorem on the boundary integral renders the corresponding local differential equation and its jump condition as It is also useful to present an alternative non-conservative form of the differential equation for J  .This is achieved by enforcing the divergence-free conditions DIV  H  = 0 (i.e Piola identity) in (24), and after some simple algebraic manipulations, this gives

Conservation of linear momentum
In order to describe the motion of a solid, we consider the usual global conservation of linear momentum p X =  X v (per unit material volume) within an arbitrary material domain Ω X , with boundary Ω X , formulated as where f X are the body forces per unit material volume and P is the first Piola Kirchhoff two-point stress tensor.Making use of the global transformation relation derived in (13), the ALE linear momentum conservation equation written in integral form becomes Using the Gauss theorem on the boundary integral of (27) gives together with the jump conditions shown as

Conservation of total energy or first law of thermodynamics
It is possible to include the total energy conservation law into our system of conservation laws as 40 where E X represents the total energy per unit of material volume.Even though in the current work we only consider the case of adiabatic deformation where both the heat flux and the heat source terms are neglected, the above energy is still useful in order to evaluate accurately numerical dissipative effects † .Utilising the transformation rule previously derived in (13), the equivalent ALE global conservation becomes Consequently, the corresponding local form of the total energy conservation law and its jump conditions are and

Evaluation of deformation gradient tensor
In the standard ALE solid formulation, 19,36 the physical deformation gradient is generally computed based upon the referential gradient of the respective geometries, mathematically defined as F = )(

𝜕𝚿 𝜕𝝌
) −1 .This implies that the accuracy of the physical deformation gradient tensor relies heavily on the resolution of the respective referential gradient computations.9][50] This approach has been shown to be very effective when simulating solids undergoing large deformations, 51,41,52 with a possibility of experiencing very large deformations and shocks. 44n order to derive the conservation law for the spatial deformation F  , we can first integrate Equation ( 5) over an arbitrary reference volume, followed by the use of the Gauss theorem to give where I is the second order identity tensor.Differentiating in time gives a global conservation equation for the spatial deformation gradient as The local differential equation and jump conditions are Following the exact same procedure discussed above, the local conservation law for the material deformation gradient tensor F  as well its associated jump conditions are

A unified system of first-order conservation laws for solid dynamics
Combining the physical conservation Equations ( 28), ( 24) and ( 32) and the supplementary geometric conservation Equations ( 36) and ( 37), a complete set of first-order conservation equations written in a compact manner in an ALE description (with respect to the referential domain) can be established as with vector of conserved variables   , flux vector  I  and source term   described by where, for the sake of compactness, we define Above equations (40) summarise the relationship between material and referential fields.As a result, the corresponding flux vector associated with the referential unit normal N  can be expressed as In addition to initial and boundary (essential and natural) conditions required for the complete definition of the initial boundary value problem, closure of system (38) requires the introduction of a suitable constitutive law relating stress P and deformation gradient tensor F. This constitutive law must fulfill a series of physical requirements, including thermodynamic consistency (via the Coleman-Noll procedure 53,54 ) and the principle of objectivity. 4In the case of reversible elasticity, this can be achieved via the introduction of a suitable strain energy potential Ψ such that stress P and deformation gradient tensor F are related via P = Ψ(F) F .Moreover, stability requires that the strain energy potential Ψ is rank-one convex (or elliptic) in the sense of Legendre-Hadamard. 55,56ery importantly, enforcement of suitable kinematic restrictions permits the degeneration of above ALE system (38) into three alternative systems of first-order conservation equations suitable for solid dynamics, namely, Total Lagrangian formulation, Eulerian formulation and the recently proposed Updated Reference Lagrangian formulation, 31 which involves the concept of incremental kinematics.Details of their derivations are given in Remarks 3,4, and 5.The ALE system (38) can thus be seen as an elegant generalisation of other existing continuum conservation laws descriptions.
Remark 3. In order to recover the Total Lagrangian description, it is important to ensure both referential and material domains overlap by enforcing W = 0 and F  = I, which in turn yields  = X, v = v, F  = F.These kinematic restrictions enables the ALE system (38) to become Remark 4. The respective Eulerian system can also be obtained by enforcing the referential domain to coincide with the spatial domain by setting v = 0 and F  = I, leading to  = x, W = −F −1 v and F  = F −1 .These conditions then imply the ALE system (38) to be degenerated into References 57,58,30 where the variables per unit of spatial volume are and the Cauchy stress tensor  = PH −1 , with div x denoting the divergence operator defined with respect to the spatial configuration.
Remark 5. Another recently proposed Updated Reference Lagrangian formulation, 31,35 which requires the introduction of an intermediate configuration for the multiplicative decomposition of the conservation variables, can be simply recovered.This is achieved by setting W = 0 but now with the given (or prescribed) material mapping functions such as . With these at hand, the ALE system (38) reduces to the so-called Updated Reference Lagrangian system summarised as

HYPERBOLICITY: TRAVELLING PLANE WAVES AND JUMP CONDITIONS
The stability of above system (38) of conservation laws (and of any of its degenerate systems counterparts) is fundamentally dependent upon the specific choice of constitutive model (i.e.hyperelastic strain energy function Ψ).As stated above, the requirement of rank-one convexity is sufficient to ensure stability, guaranteeing as a consequence both the existence of plane travelling waves within the solid and the propagation of strong discontinuities at physical shock speeds ‡ .These two considerations are presented in depth below.

Rank-one convexity: Legendre-Hadamard condition or hyperbolicity
The study of the hyperbolicity 60,61 of the system described in (38) is crucial in order to guarantee the existence of real wave speeds for the entire range of deformation states.In this work, it is instructive to examine the hyperbolicity of system (38) particularised for the case of a given material mesh motion, that is, F  and J  will be assumed to be known ab initio or prescribed.This implies that the unknown variables {F  , J  } and their corresponding geometric conservation laws become redundant in (38), thus can be removed from the hyperbolicity analysis presented below.The eigenvalues (or wave speeds) of the reduced conservation system (i.e.comprised of unknown variables p  and F  ) can be determined by identifying possible plane wave solutions (in the absence of source terms) of the type 26 given by where  represents a scalar real valued function, U   represents the wave speed (measured with respect to the referential configuration) corresponding to the eigenmode and N  is the direction of propagation of the wave in the referential configuration.Substitution of above expression (45) into system (38) leads to the following characteristic equations where  N  can be identified with the so-called flux Jacobian matrix.Alternatively, above Equation ( 46) can be re-written by utilising the concept of directional derivative 4 as follows where D(•) [] represent the directional derivative of a field (•) with respect to a given variation . 4 Considering each of the conservation laws for variables {p  , F  } in system (38), it yields Post-multiplication of (48b) by F −1

𝚿 and introducing
where the kinematic relationship between the measures of deformation has been used in the third line of above Equation ( 49) and the velocity relationship (10) has been used in the final line of (49).It is possible to show that where the velocity term v  has been introduced for convenience.Furthermore, introducing the Nanson's rule expression Λ H  N X = H  N  , and substituting it along with (50) into (49), results in In order to further simply above expression (51), recall 40 that a similar plane wave type substitution into Equation (42c) renders where U X denotes the wave speed measured with respect to the material configuration.This above expression enables to realise that which when substituted into the last term on the right hand side of Equation ( 51) results in Comparison of Equations ( 52) and ( 54) permits to infer that which summarises the relationship between the wave speed measured in referential U   and material U  X configurations.Crucially, above formula (55) allows to straightforwardly make use of the closed-form expressions for material wave speeds U  X presented by the authors in previous publications. 40e turn now our attention to Equation (48a) by expanding the first term on its right hand side as where the last of the equations in (40) has been used in the first line of (58) and the Nanson's rule in the second line.
Similarly, expansion of the second term on the right hand side of (48a) yields where  =  2 Ψ FF is the constitutive tensor.Substitution of above two Equations ( 56), ( 57) into (48a) gives Re-arrangement of above Equation ( 58) renders which represents the well-known wave propagation material eigenvalue problem and  NN is the so-called acoustic tensor.Naturally, existence of real wave speeds requires the acoustic tensor to be positive definite, which is ensured if the strain energy potential is rank-one convex or elliptic, which is known as the Legendre-Hadamard condition. 59,55,56

Rank-one convexity: Rankine-Hugoniot or jump conditions
As in Section 6.1, the material motion is assumed to be known ab initio.The Rankine-Hugoniot jump conditions for the total energy (33) can be recalled and re-written in a material description as To proceed, it is necessary to explore the algebraic expression for the jump of the product of any two arbitrary variables a and b such that ⟦ab⟧ = ⟦a⟧b Ave + a Ave ⟦b⟧ with 57 where [•] L and [•] R represent values at either side of the discontinuity.Applying the above jump property onto the terms containing velocity jump (i.e.⟦v ⋅ v⟧) and the jump between the first Piola-Kirchoff stress tensor and the velocity (mathematically expressed as ⟦P T v⟧), gives Observing also that ⟦•⟧ = ⟦•⟧ L + ⟦•⟧∕2 and removing the subscript index (•) L for simplicity, above expression can be re-written as Following the notation used in Reference 56, introducing −⟦v⟧∕U X = u, that is ⟦F⟧ = u ⊗ N X , and making use of expression (3), it is easy to show that ⟦P⟧ ∶ ⟦F⟧ =  R U 2 X (u ⋅ u).With this at hand, above expression becomes which is equivalently written as The above inequality (64) holds for a stored energy function satisfying rank-one convexity.Equations ( 63) and ( 64) imply that strong discontinuities in the medium travel with real speeds, in the same way as in Section 6.1, acoustic disturbances (plane waves) also travel at physical wave speeds.

HAMILTONIAN AND SECOND LAW OF THERMODYNAMICS
In order to provide a proper physical meaning to the conjugate fields of the underlying system, consider the Hamiltonian   per unit of referential volume defined by References 62,35,63 with   = J   R and   (, t) and  (p  , F  , F  , J  ) represent alternative functional representations of the same magnitude.In above equation, the first term on the right hand side represents the kinetic energy and the second term represents the strain energy potential per unit of volume measured in the referential configuration.For the case of isothermal reversible elasticity (which is the case of the current work), the notion of Hamiltonian is simply the total energy per unit reference volume.It is now possible to introduce energy conjugate fields to the two deformation measures {F  , F  } defined as and as a result, Here,  X (v, F) = (v) − Ψ(F) denotes the usual material Lagrangian function, defined as the difference between the kinetic energy ) and the strain energy potential measured in the material configuration.It is instructive to revisit the global version of the second law of thermodynamics when written in terms of the Hamiltonian energy density  .The time derivative of the Hamiltonian is obtained via the chain rule as where, Equations ( 66), ( 67) and ( 25) have been substituted in the second and third lines of (68), respectively, and  = − ( F T P +  X I ) denotes the classical Eshelby stress tensor. 55Pairs said to be dual or energy conjugates with respect to the referential volume in the sense that their inner product yields energy rate per unit of referential volume.For instance, the energy conjugate field to the time rate of the material-based deformation gradient is the Eshelby stress tensor.Consequently, we can substitute the linear momentum Equation ( 28) and the geometric conservation Equations ( 36) and ( 37) into the fifth line of ( 68), and after some algebraic manipulations, it gives Moreover, noticing that ) , above equation after some arrangement reduces to Applying the Gauss theorem on the divergence term above, expression (70) becomes where Πext  denotes the external work given by Note that in Equation ( 71), the material surface integral vanishes as the mesh velocity W has to be tangential to the material surface to avoid changes in physical volume, that is W ⋅ dA X = 0. Above equality (71) represents the global statement of the first law of thermodynamics.Furthermore, when considering irreversible processes (e.g., irreversible thermo-mechanics or inelasticity), above equality transforms into the following inequality This represents a valid expression for the second law of thermodynamics of a system.Satisfaction of inequality ( 73) is a necessary ab initio condition to ensure stability, otherwise referred to as the classical Coleman-Noll procedure.This concept will be further exploited in Section 8.2 when introducing appropriate numerical dissipation to the finite volume spatial discretisation employed in this work.

Semi-discrete formulation for ALE solid dynamics
The vertex centred finite volume spatial discretisation presented in this work requires the introduction of a median dual mesh 41,51,64,65 for the definition of control volumes (see Figure 2).With this in mind, and making use of Here,  a  and  a  are the average values of both the conservation variables and source term vector within the fixed referential control volume, respectively.Moreover, the surface integral of ( 74) is approximated via a second-order central difference scheme, resulting in Here, the representation of the internal nodal force is defined as C ab  and the boundary traction t B a = P a  N a  .The above mixed-based system (76a-g) is prone to exhibit non-physical numerical instabilities 10,66,67,68,69,32,33,61 when attempting to model large strain solid dynamics problems. 1,70,71To remedy this, an upwind based stabilisation term must be suitably derived (via the semi-discrete version of the classical Coleman-Noll procedure 72 ) and introduced to the internal nodal force T a  , presented as Such stabilisation term will be discussed in the following section.

Numerical entropy production
The semi-discrete version of inequality ( 73) is where suitable conjugate fields (66) and ( 67) have been substituted in the first line of (78).We can now substitute the linear momentum Equation (76a) and the geometric conservation equation for both spatial and material deformation gradient tensors (76b) and (76c) into the second line of ( 78), and after some algebraic manipulations, to give Here, Πext represents the semi-discrete version of power contribution described as Finally, the semi-discrete version of the Hamiltonian (79) reduces to To ensure the satisfaction of the second law of thermodynamics at semi-discrete level, it is the objective to prove that the numerical dissipation term on the right-hand side of ( 81) is non-negative, that is  total ≥ 0. This can be achieved by equivalently swapping indices a and b to give Simply averaging the first term and the second term of the expression above, and noting the anti-symmetric nature of the stabilisation term as  ba p  = − ab p  , an alternative expression for  total is Sufficient conditions guaranteeing  total ≥ 0 are where S ab is a positive semi-definite stabilisation matrix.Following the previous work of the authors, 40,44 the dissipation term used is chosen as where It is very interesting to observe how the dissipation term is related to the difference in velocity between interacting nodes, typical of Riemann solver based upwinding terms. 44,73,74,75

TIME INTEGRATION
Insofar as the resulting set of semi-discrete equations is rather large, it will only be appropriate to employ an explicit type of time integrator.In this work, a three stage Runge-Kutta explicit time integrator 31 is used.This is described by the following time update equations from time step t n to t n+1 ) , ) ) , ) ) . ( Additionally, both material and spatial geometries are also updated via the same time integrator presented above, resulting in a monolithic time update procedure in which the unknowns   = ( p  , F  , F  , J  ) T are updated along with the material geometry  and the spatial geometry  through Equation ( 86).Moreover, the maximum time step Δt = t n+1 − t n is governed by a standard Courant-Friedrichs-Lewy (CFL) 76 condition given by where  CFL is the CFL stability number, c p, is the pressure wave speed measured in referential domain and h min represents the minimum (or characteristic) length within the computational domain.Notice that the relationship between the referential pressure wave speed c p, and the material pressure wave speed c p,X is already presented in expression (55).The closed-form expression for c p,X can be obtained, and which has been reported in Reference 28.Unless otherwise stated, a value of  CFL = 0.9 has been chosen in the following examples to ensure both accuracy and stability.

ALGORITHMIC DESCRIPTION
For ease of implementation, Algorithm 1 summarises the algorithmic description of the proposed Arbitrary Lagrangian Eulerian vertex-based finite volume methodology for hyperelasticity at finite strains.One attractive feature of the ALE algorithm is the ability to address excessive element distortion by suitably re-adapting the mesh.This is especially the case when extending to plasticity problems, and which will be explored in forthcoming publications.In the current work, we provide analytically the material mesh motion at every time step of the time integration process.The objective is to assess the effectiveness and robustness of the overall algorithm.

NUMERICAL EXAMPLES
In this section, a wide spectrum of three dimensional benchmarked examples are presented in order to assess the performance, effectiveness and applicability of the proposed ALE algorithm described above.With these examples the aim is to show that the overall ALE formulation • achieves second order convergence for both velocities (or displacements) and stresses (or strains), ) ENFORCE strong boundary conditions on velocities: v a and va (9) EVALUATE first Piola P a end (10) EXPORT results for this time step (11) ADVANCE in time end • ensures the satisfaction of the Geometric Conservation Law (GCL) condition, 10,47,65,17,77 • guarantees a non-negative rate of production of entropy in an isolated system, and • circumvents locking difficulties and hour-glassing. 78,79 the following numerical computations, and to examine the robustness of the proposed algorithm, we prescribe the sinusoidal-type material mesh motion on a cuboid described by ) sin where χ1 =  1 +  C 1 and χ2 =  2 +  C 2 represent the coordinates shifted with respect to the centroid of the domain  C ,  is a parameter determining the magnitude of the mesh motion and T denotes the period of the sinusoidal functions.Differentiating expression (88) in time gives the material mesh velocity W expressed as ) sin ) cos ) sin Both material mesh motion (88) and material mesh velocity (89) are carefully designed to only allow movement tangentially on the material boundary surface in order to avoid changes in physical volume.This is illustrated in Figure 3.

F I G U R E 3
Visualisation of the analytical material motion (88) on a unit cube with a magnitude  = 10 −2 m.A line of eight nodes in one of the faces is highlighted in red so that their displacements can be easily followed.
To account for large and reversible deformations, a neo-Hookean hyperelastic model is used.The associated strain energy functional Ψ(F) can be conveniently decomposed into the summation of a deviatoric component Ψ(J −1∕3 F) and a volumetric component U(J) given by (90) where  and  are the shear modulus and bulk modulus.The Lamé first parameter  can be subsequently computed as  =  − 2 3 .In previous publications, 45 the authors have shown that this model is polyconvex and, thus, rank-one convex.

Satisfaction of geometric conservation laws: Patch test
The first example corresponds to a standard rigid body translation test.A column of unit 1 m squared cross section and of length H = 6 m is initially prescribed by a uniform velocity field v 0 = [3, 1, 0] T ms −1 across the entire domain.All the boundaries of the structure are left free.See Figure 4 for the problem setup.The primary aim of this example is to check the capability of the proposed ALE algorithm in accurately reproducing rigid body translation behaviour by arbitrarily moving the mesh, thus satisfying the so-called discrete Geometric Conservation Law (GCL) conditions.For instance, we expect the components of the first Piola Kirchhoff stress tensor remain null and the velocity components remain exactly the same during the entire simulation.In this example, a neo-Hookean material is used.Its corresponding material properties and the mesh motion parameters are presented in Table 1.
To begin with, it is interesting to show graphically the mesh motion according to the expression described in (89).This is displayed in Figures 5 and 6.For qualitative comparison purposes, we monitor the time history of: (i) the global error related to the difference between the Jacobian and unity (i.e.J − 1), and (ii) the global errors related F I G U R E 4 Problem setup for the translation test case: Geometry and its dimension.to the difference between the three components of the velocity and the initial velocity field v 0 .These can be seen in Figure 6.As expected, these two errors described above fluctuate around machine accuracy when using the complete ALE algorithm presented in (76a-g).Relatively large errors are observed when considering a reduced system by neglecting some of the conservation variables related to the material mesh motion.For instance, not solving J  in system (76a-g) would violate the GCL condition associated with the volumetric behaviour.Clearly, a non-smooth profile for J − 1 is observed in Figure 7. Numerical errors can be further magnified (such as incorrect deformation path) when not solving the variables {F  , J  } of the system (76a-g), which indeed can potentially lead to the breakdown of the numerical scheme.

Mesh convergence: Low dispersion swinging cube
The objective of this example is to assess the order of convergence of the proposed ALE framework.As already explored in References 48,50,69,38,64,26, a cubic shaped domain of unit length is considered.The exact field of the physical mapping function (X, t) is specifically chosen as ) cos ) cos ) sin ) cos ) cos ) sin when using (red) the complete ALE system, (green) a reduced system without solving variable J  and (blue) further reduced system without solving {F  , J  }.A discretisation of 4 × 24 × 4 linear tetrahedral elements.A neo-Hookean model is used with parameters listed in Table 1.

F I G U R E 7
Translation: A sequence of deformed states at time t = {0, 0.8, ..., 4.8} s (from left to right for each row) when using (A) the complete ALE, comparing against the reduced system without considering (B) J  and (C) F  and J  .Colour contour indicates the error of ||J − 1|| L 2 .A discretisation of 4 × 24 × 4 linear tetrahedral elements.A neo-Hookean model is used with parameters listed in Table 1.
with the parameters A = B = C = 1 to ensure the existence of a non-zero pressure distribution.For values of U 0 below 0.0001 m, the solution can be considered to be linear.Consequently, the exact velocity field v and deformation gradient tensor F can be computed as Dirichlet boundary conditions compatible with the exact field  ex (92) are applied on the boundary of the domain.In particular, the cube has symmetric boundary conditions (i.e.restricted to tangential movement) at the faces X 1 = 0, X 2 = 0 and X 3 = 0 and skew symmetric boundary conditions (i.e.restricted to normal movement) at the faces X 1 = 1, X 2 = 1 and X 3 = 1.The problem is initialised with appropriate velocity and deformation gradient distributions such that Use of the mesh velocity described in (89) enables the spatial velocity field to be expressed as A list of parameters used for this simulation is summarised in Table 2.As compared to the closed form solutions described in (93), Figures 8 and 9 show the L 1 and L 2 norm convergence analysis at time t = 1 × 10 −3 s for the components of the velocity v and the components of the first Piola-Kirchhoff stress tensor P. As expected, the proposed computational framework achieves equal second order convergence for all the variables solved, namely velocity and the stress tensor.This equal order convergence for all variables (especially, stresses/strains) is one of the advantages of the proposed framework.For completeness, we have also included their corresponding detailed results for both L 1 and L 2 norm errors in Tables 3  and 4.

Conservation properties and the discrete satisfaction of second law: L-shaped block
4][85] Another interesting objective of this problem is to examine if the proposed algorithm ensures the discrete satisfaction of the second law of thermodynamics by measuring the total numerical dissipation introduced by the algorithm.The geometry of the problem is depicted in Figure 10.The L-shaped block is subjected to a pair of time-varying boundary forces applied on two boundary faces, which are described as

F I G U R E 10
Problem setup for L-shaped block: Geometry and its dimension.
A neo-Hookean model is considered.The associated material parameters used in the simulation can be found in Table 5.Three different levels of progressive mesh refinement are considered: {M1, M2, M3} comprising {5616, 18954, 44928} number of linear tetrahedral elements, respectively.Regarding the mesh motion, we artificially move the material mesh using the mapping described in (89).For clarity, a sequence of material mesh motion is displayed in Figure 11.
First, a mesh refinement study is carried out.This is seen in the first three columns (from left to right) of Figure 12.The deformation pattern simulated using a relatively coarse mesh (M1), is in good agreement with the results obtained using finer discretisations (M2 and M3 models).Improved pressure resolution is naturally observed.For benchmarking purposes, results obtained with an alternative in-house Total Lagrangian vertex-based finite volume algorithm 41,44 with M3 discretisation are also included and compared.Comparing the results of the proposed ALE algorithm and those of the Total Lagrangian finite volume algorithm, near-identical results are seen (see Figure 12).  5.
times.The components of the total angular momentum slightly decrease after the loading phase t > 5 s, although this is very difficult to notice § .In addition, Figure 13C depicts the evolution of various energy measures.This include kinetic energy, elastic strain potential energy and the total energy (which is the summation of kinetic energy and elastic strain energy).A slight decrease in the total energy is unavoidable after the loading phase due to the addition of upwinding-based numerical dissipation (77) to the system.The total numerical dissipation continuously increases throughout the entire simulation, and hence ensuring the discrete satisfaction of the second law of thermodynamics.This is seen in Figure 13D.Third, and for qualitative comparison purposes, Figure 14 monitors the time evolution of the horizontal velocity component v x at position X = [0, 10, 0] T and position X = [6, 0, 0] T .The solution converges with a successive level of  ALE results are obtained using three different meshes (namely M1, M2 and M3), comparing against alternative Total Lagrangian vertex-based finite volume algorithm. 41A neo-Hookean model is used.Parameters related to material properties and the mesh motion are summarised in Table 5.
refinement.Finally, a sequence of deformed states are depicted Figure 15, where the colour contour indicates the pressure profile.Stable solutions are observed even after a relatively long-term response.

11.4
Robustness of the proposed formulation

Bending column
Another well documented example 1,86,57 is considered in this section.The deformation of a column with a unit squared cross section and of length H = 6 m is initiated with a linear velocity profile across the entire length of the column given by where V = 10 m∕s.The column is clamped at the bottom surface and is left free on the rest of the boundaries (see Figure 16).The primary aim of this example is to illustrate the performance of the proposed ALE algorithm in bending dominated scenario, regardless of the magnitude of the mesh motion.With respect to the material model, a neo-Hookean constitutive law is used, with material parameters being tabulated in Table 6.In this example, three different model refinements are explored, namely (M1) 2304, (M2) 18,432 and (M3) 147,456 linear tetrahedral elements.
Figure 17 shows how the material mesh deforms at time t = 0.2 s (according to expression (89)), with four different magnitudes  = {0, 1, 1.5, 2} of the mesh motion.The Total Lagrangian formulation can be simply recovered by enforcing the value of  = 0.As depicted in Figure 18, even with a different mesh motion magnitude, almost identical results (e.g., deformed shapes and pressure contour) are observed.In Figure 19, we monitor the time history of the horizontal velocity component v x and also one of the component of the first Piola P xX .The ALE results show extremely good agreement with the results obtained using the alternative Total Lagrangian counterpart.Figure 20A illustrates the time history of F I G U R E 15 L-shaped block: A sequence of deformed structures with pressure distribution at times t = 0, 1, 2, ..., 19 s (from left to right and top to bottom), respectively.Results obtained using a neo-Hookean model with M3 mesh.Parameters related to material properties and the mesh motion are summarised in Table 5.

F I G U R E 16
Bending column: Problem setup for bending column: geometry and its dimension.

TA B L E 6
Bending column: Material properties and parameters of the mesh motion.5.  5.

F I G U R E 18
the kinetic energy, strain potential energy and total energy (which is the summation of kinetic and strain energies).The total numerical dissipation introduced by the overall algorithm is non-negative and accumulates over time.This is shown in Figure 20B, indicating the discrete satisfaction of second law of thermodynamics.The dissipation of the solution is analysed in Figure 20C,D.Notice that the difference between the numerical total energy and the conserved total energy is the actual dissipation of the algorithm.It can be seen that the numerical dissipation is clearly reduced as the mesh is refined.A neo-Hookean model is used with material properties and the mesh motion parameters summarised in Table 6.

(A) (B) (C) (D)
F I G U R E 20 Bending column: Time history of (A) various energy measures and (B) the numerical dissipation introduced by the algorithm using three level of mesh refinements using three values of .(C) and (D) shows the time history of numerical dissipation and its percentage, respectively.A neo-Hookean model is used with material properties and its mesh motion related parameters are tabulated in Table 6.

F I G U R E 21
Twisting column: Problem setup, geometry and its dimension.

F I G U R E 22
Twisting column: Comparison of the deformed shapes at time t = 100 ms.The first three columns (left to right) show the mesh 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 in-house Total Lagrangian vertex-based finite volume algorithm. 41Colour indicates the pressure profile.A neo-Hookean model is used.The parameters related to material properties and mesh motion are summarised in Table 6.6.

F I G U R E 24
Twisting column: Comparison of the deformed structures at t = 140 ms for  = 0, 1, 1.5, 2 (from left to right) via (M3).
Colour indicates the pressure distribution.A neo-Hookean model is used.The parameters related to material properties and mesh motion are summarised in Table 6.

F I G U R E 25
Twisting column: A series of snapshots at times t = 100,300, 500 ms from the bottom view of the deformed column for (First row) Total Lagrangian vertex-centred finite volume algorithm 41 and (Second row)  = 2. Time evolution of the angles at positions A, B, C, D are also monitored.

F I G U R E 26
Twisting column: Time evolution of (A) different energy measures, (B) total numerical dissipation (%) introduced by the proposed ALE method via M3 for  = 1, 1.5, 2, compared to a in-house Total Lagrangian vertex centred algorithm.(C) and (D) show the time history of the numerical dissipation and its percentage for three levels of refinement.A neo-Hookean model is used.Parameters related to the material properties and the mesh motion are tabulated in Table 6.

Twisting column
Similar to the previous bending example, the column is now twisted with an initial sinusoidal velocity field relative to the origin given by where Ω 0 = 105rads −1 represents the magnitude of initial angular velocity (see Figure 21).Identical material properties and mesh motion parameters presented in Table 6 for bending example are used.First, a mesh refinement study for the proposed ALE algorithm is carried out.In Figure 22, the deformation pattern of the structure predicted using a small number of elements (M1) agrees extremely well with the results obtained using finer discretisations (M2 and M3).For benchmarking purpose, a Total Lagrangian vertex-centred finite volume algorithm 41 discretised with M3 is also included and compared.Practically identical results are observed.This is seen in Figure 23.Additionally, in order to examine the robustness of the algorithm, we run the problem using four different magnitudes of the mesh motion, namely  = {0, 1, 1.5, 2}.Almost identical results are obtained (see Figure 24).No instabilities are observed in any of the simulations, emphasising the importance of the tailor-made numerical stabilisation term.For qualitative comparison, Figure 25 monitors the evolution of the accumulated angle at four different positions.Comparing the proposed algorithm with the published Total Lagrangian finite volume algorithm, 41 no noticeable difference is observed.As shown in Figure 26, the overall numerical dissipation introduced by the ALE algorithm accumulates over time throughout the entire simulation, ensuring the discrete satisfaction of the second law of thermodynamics.

CONCLUSIONS
This paper has introduced a novel Arbitrary Lagrangian Eulerian (ALE) description for the modelling of solid dynamics using conservation laws.Physical conservation laws (i.e., mass, linear momentum and total energy density) are supplemented with two geometric conservation laws for the solution of the deformation gradient tensors {F  , F  } between the referential and material and between referential and spatial configurations, respectively.Exploiting the multiplicative decomposition of the physical deformation gradient tensor F into the pair {F  , F  } permits the evaluation of the strains and stresses with equal order convergence that those of velocities and displacements.The hyperbolicity of the system, as well as the physical propagation of shocks, is demonstrated by requiring rank-one convexity of the underpinning hyperelastic strain energy potential, establishing the relationship of the wave speeds between material and referential configurations.The framework elegantly degenerates into alternative continuum descriptions, namely, Total Lagrangian, Eulerian and Incremental Updated Lagrangian formulations.From the discretisation standpoint, a vertex centred Finite Volume algorithm is used, with interface numerical fluxes carefully defined to ensure non-negative numerical entropy production.This is demonstrated using the discrete variation of the Hamiltonian of the system and monitored numerically by solving the total energy density as an additional decoupled conservation law.A three-stage Runge-Kutta explicit time integrator is used to evolve the system of discrete equations in time.The accurate computation of wave speed bounds is exploited when establishing the time step size for the stability of the time integrator.
A range of numerical examples is presented in order to demonstrate the conservation, stability and order of convergence of the resulting algorithm.Examples are restricted to the case of reversible hyperelasticity, in order to showcase the potential of the formulation before exploiting more complex irreversible inelasticity constitutive models in forthcoming publications.A user-designed material mesh motion is pre-defined (i.e., the material mesh motion is not solved) and ALE results are benchmarked against an alternative Total Lagrangian formulation, demonstrating the robustness of results.This paper paves the way for the simulation of large deformation inelastic processes where the ALE formulation can be exploited to its maximum potential.The extension of our framework into thermo-visco-plastic deformation processes has already shown excellent results and will be reported in our next publication.

1 2 [
where b ∈ Λ a represents the set of neighbouring control volumes b associated with the control volume a and C   = A  3 N  represents the (tributary) boundary area vector.For a given edge connecting nodes a and b, the mean undeformed area vector C ab  satisfies the reciprocal relation, that is C ab  = −C ba  .The terms within the parenthesis in (75) correspond to the evaluation of the control volume internal interface flux  I N ab  and boundary fluxes  B ,a .This evaluation is comprised of a summation over edges (first term in the parenthesis) and a summation over boundary faces (second term in the parenthesis).The internal flux  I N ab  =  N  ( a ) +  N  ( b ) ] denotes the average states of the left and the right control volumes of a given edge ab.The boundary flux  B a is enforced through either Neumann or Dirichlet boundaries.Expression (75) is now particularised for each individual component of   under consideration, yielding

F I G U R E 5
Translation: Time evolution of the prescribed mesh motion at time t = {0, 40, 80,120, 160,200} ms (from left to right).(A) (B) F I G U R E 6 Translation: Time evolution of the L 2 global errors in ||J − 1|| L 2 and velocity components ||v i − v 0,i || L 2 (with i = {x, y, z})

− 2 F
Second, Figure 13A,B show the capability of the proposed ALE algorithm in preserving the conservation of global linear and angular momenta.The global linear momentum is expected to oscillate around zero value (machine error) at all 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License TA B L E 5 L-shaped block: Material properties and parameters of material mesh motion.I G U R E 11 L-shaped block: For visualisation purposes, we display the two-dimensional view of a series of material mesh deformation described in (89) at time t = 0.1, … 0.4 s (from left to right).F I G U R E 12 L-shaped block: Comparison of deformed shapes at time t = 24 s.The first three columns (left to right) show the mesh refinement of a structure simulated using the proposed ALE algorithm, whereas the last column (on the right) shows results via an alternative in-house Total Lagrangian vertex-based finite volume algorithm.41Colour indicates the pressure profile.A neo-Hookean model is used.The parameters related to material properties and mesh motion are summarised in Table5.13 L-shaped block: Time evolution of (A) global linear momentum, (B) global angular momentum, (C) different energy measures and (D) total numerical dissipation (%) introduced in the algorithm.A neo-Hookean model is used.Parameter related to the material properties and the mesh motion are tabulated in Table

Young− 2 F
I G U R E 17 Bending column: Snapshots of the mesh movement at time t = 0.2 s using four different values of  = {0, 1, 1.5, 2} (from left to right).A neo-Hookean model is used with material properties and the mesh motion related parameters summarised in Table Bending column comparison of the deformed structures at (Top) t = 0.2 s and (Bottom) t = 1.1 s.Different values of  are studied (via M3).The first column corresponds to  = 0, the second column denotes  = 1, the third column  = 1.5 and the fourth column  = 2. Colour indicates the pressure distribution.A neo-Hookean model is used with material properties and the mesh motion parameters presented in Table

19
Bending column: Time history of the (A) horizontal velocity component v x and (B) the component of the stress tensor P xX at point A (i.e.X = (0.5, 6, 0.5) T m).
Twisting column: A sequence of deformed structures at times t = 20, 40, 60, 80,100 ms (from left to right) for (Top row) Total Lagrangian vertex centred finite volume algorithm 41 and (Bottom row)  = 2, respectively.Colour indicated pressure distribution.A neo-Hookean model is used.Parameters related to the material properties and the mesh motion are tabulated in Table 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

)
10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License , 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Control volumes for (A) an interior node and (B) a boundary node for a median dual tessellation visualised in 2D.10970207 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Algorithm 1. ALE vertex-based finite volume algorithm Input : initial geometry  a and initial states of p a  , F a  , F a  , J a  , E a  and W a for all vertices Output: current material geometry  a , current spatial geometry  a , physical velocity v a , material velocity W a , and current states of F a , P a and E a Translation: Material properties and parameters of material mesh motion.
Swinging cube: Material properties and parameters of material mesh motion.Swinging cube: L 1 and L 2 global convergence analysis at time t = 1 × 10 −3 s for the components of velocity.Results obtained using a neo-Hookean model and the material properties used are summarised in Table2.10970207,0,Downloadedfromhttps://onlinelibrary.wiley.com/doi/10.1002/nme.7467byUniversityOfGlasgow,Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)onWileyOnlineLibraryfor rules of use; OA articles are governed by the applicable Creative Commons License Swinging cube: L 1 and L 2 global convergence analysis at time t = 1 × 10 −3 s for the components of the first Piola Kirchhoff stress tensor.Results obtained using a neo-Hookean model and the material properties used are summarised in Table2.Swinging cube: Numerical values for the L 1 norm error of the velocity components v x , v y , v z and the diagonal components of the first Piola Kirchhoff stress P xX , P yY , P zZ .Convergence rates are calculated using the results of the two finest meshes.A neo-Hookean model is used with parameters summarised in Table2.
TA B L E 3Note: Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Swinging cube: Numerical values for the L 2 norm error of the velocity components v x , v y , v z and the diagonal components of the first Piola Kirchhoff stress P xX , P yY , P zZ .Convergence rates are calculated using the results of the two finest meshes.A neo-Hookean model is used with parameters summarised in Table 2. Note: 10970207, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7467by University Of Glasgow, Wiley Online Library on [24/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License