Particle-based solid for nonsmooth multidomain dynamics

A method for simulation of elastoplastic solids in multibody systems with nonsmooth and multidomain dynamics is developed. The solid is discretised into pseudo-particles using the meshfree moving least squares method for computing the strain tensor. The particle’s strain and stress tensor variables are mapped to a compliant deformation constraint. The discretised solid model thus fit a unified framework for nonsmooth multidomain dynamics simulations including rigid multibodies with complex kinematic constraints such as articulation joints, unilateral contacts with dry friction, drivelines, and hydraulics. The nonsmooth formulation allows for impact impulses to propagate instantly between the rigid multibody and the solid. Plasticity is introduced through an associative perfectly plastic modified Drucker–Prager model. The elastic and plastic dynamics are verified for simple test systems, and the capability of simulating tracked terrain vehicles driving on a deformable terrain is demonstrated.


Introduction
We address the modelling and simulation of multidomain systems with nonsmooth dynamics and deformable solids of particulate nature. This includes mechatronic multibody systems with nonsmooth dynamics, such as vehicles, robots, and processing machinery that handle or operate on soil or B Martin Servin martin.servin@umu.se 1 Umeå University, Umeå, Sweden bulk materials. Fast multidomain simulation is useful for concept design exploration, development of control strategies, and for interactive real-time simulators, e.g. for operator training, human-machine interaction studies, and hardwarein-the-loop testing.

Nonsmooth multidomain dynamics
Multidomain system simulation requires integration of multiple heterogeneous subsystems into a single full-system model. Mechatronic systems are typically composed by rigid and flexible multibodies coupled by kinematic constraints for modelling of joints and sets of differential algebraic equation (dae) models for electronics, hydraulics, and powertrain dynamics [1]. If the subsystems are not strongly coupled, the full-system simulation can be realised by means of co-simulation [2]. At coarse timescales, however, the dynamics need to be treated as nonsmooth [3][4][5]. This means that velocities are allowed to be arbitrarily discontinuous and impulses-from impacts, joint limits, and frictional stick-slip transitions-can propagate instantly throughout the system. This enables implicit time integration using large timesteps. To ensure numerical stability and computational efficiency, the full system must be approached with a consistent mathematical framework and adequate numerical solvers. Current techniques for co-simulation do not support this for strongly coupled systems with nonsmooth dynamics.

Particulate-based solids
The versatile and widely used discrete element method (DEM) [6] for simulating particulate solids have been extended to the framework of nonsmooth dynamics [4,5,7,8]. When the lengthscale of deformations is much larger than the particle size, it is more computationally favourable to Fig. 1 Illustration of the idea of using particle-based discretisation and compliant deformation constraints to obtain a unified multibody system model for both deformable solids and mechatronic systems model the solid as a continuum with elastoplastic constitutive law [9]. There are several arguments for using a meshfree discretisation [10,11] over traditional mesh-based methods such as the finite element method (FEM). Meshfree methods allow large deformations and topological changes without need for remeshing [12]. This makes them well suited for fracturing materials and transitions from solid into viscous and free-surface flow, or from continuum to particulate level of detail when using multi-scale methods [13,14] or model reduction techniques [15].

Our contribution
The main contribution in this paper, illustrated in Fig. 1, is a particle-based method for modelling and simulation of elastoplastic solids as a multibody system on descriptor form [2] and with nonsmooth dynamics. This enable numerical integration with large timestep size of multidomain systems containing deformable solids, mechatronics, and frictional contacts. The objective is to achieve real-time performance or better.
The solid is modelled as a particle system with strain tensor constraints introduced as the stiff limit of energy and dissipation potentials for elastic deformations, with constraint regularisation and stabilisation based on a discrete variational approach [16,17]. An associative perfectly plastic Drucker-Prager model is employed using an elastic predictor-plastic corrector strategy to detect yielding and compute the plastic flow. In its basic form, this model does not yield in hydrostatic compression in contrast to many real materials. Many soils yield under hydrostatic compression by failure in the microscopic structures whereby air and fluid are released. Therefore, the Drucker-Prager model is extended to a capped version [18] that model also plastic compaction. A meshfree method [10,11] is chosen in order to handle large deformations and, for future development, support fracturing and transitions to viscous flow or to granular media represented by contacting discrete elements. The displacement gradient and strain tensor is approximated by the method of moving least squares (mls) [19]. A stabilisation constraint is introduced to suppress spurious short wavelength modes associated with the mls approximation [20].
The dynamics of other subsystems and their potentially strong coupling is treated within the same multibody dynamics framework using a variational time integrator [16,17]. Each timestep involves solving a block-sparse mixed linear complementarity problem (mlcp) that also support modelling of frictional contacts and secondary constraints for joints and motors. The present paper is motivated by the exploration of ground vehicles on deformable terrain [21]. Current solutions enable simulation of complex vehicle models with real-time performance but not with dynamic terrain models firmly based on solid mechanics in three dimensions. Usually, empirical terrain models of Bekker-Wong type are used [22,23]. On the other hand, there exist many solutions for more accurate offline simulations with fine spatial and temporal resolution for elastoplastic solids coupled with tire models [24,25]. But few support inclusion of multibody systems and not in real time or faster.

Notations
Matrices and vectors are represented in bold face in capital and lower case, A and x, respectively. Latin superscripts indicate the index of a specific body i, j, k = 1, . . . , N , where N is the total number of bodies in the system. Greek subscripts indicate a specific scalar component of a vector or matrix. The Einstein summation convention is used where repeated indices imply summation over them, e.g. representing matrix-vector multiplication as Ax = A αβ x β . Subscripts x, y or z indicate a specific component of a vector or matrix assuming a Cartesian coordinate system. As an example of these notations, the position vector of a body i is denoted q i . The relative position of particle j to i is q i j = q j − q i . The y component of this vector is q i j y = q j y − q i y . For rigid bodies, the position vector includes also the rotation variables of the body. The system position vector is q = q 1 , q 2 , . . . , q i , . . . , q N . The gradient is denoted ∇ x = ∂ ∂ x and ∇ α = ∂ ∂ x α . The dot notation is used for the time derivative˙≡ d dt . The Cauchy stress tensor is represented by σ C and the second Piola Kirchoff stress tensor by σ .

Discrete multidomain dynamics
This section provides a summary of discrete dynamics and of the numerical integration method that is used.

Lagrangian formulation
The Lagrangian of a constrained mechanical system with position q and velocityq is where T (q,q) = 1 2q T Mq is the kinetic energy, U (q) is the potential energy, and R(q,q) is the Rayleigh dissipation function. The system state vector (q,q) is constrained by holonomic constraints 0 = g(q), with Jacobian G = ∂ g/∂q, and nonholonomic (or nonintegrable kinematic) constraints 0 =ḡ(q,q, t). The corresponding Lagrange multipliers are λ andλ, respectively. We will restrict to Pfaffian nonholonomic constraints, 0 =ḡ(q, t) ≡Ḡq−w(t). Nonholonomic constraints can, by definition, not be deduced by differentiation of a holonomic constraint [26].
To ensure numerical stability, it is common to regularise the constraints. This may be done by treating them as the limit of strong potentials and dissipation functions U ε (q) = 1 2ε g T g and R γ (q,q) = 1 2γḡ Tḡ with ε, γ → 0 [27]. In fact, any stiff force can be transformed into a regularised constraint through a Legendre transform [17] The corresponding Euler-Lagrange equations of motion are where f ≡ −∇ q U (q) − ∇qR(q,q) are the explicit forces from nonstiff potentials. This constitutes a system of dae of index 1, which are easier to integrate numerically than the corresponding higher index dae in the absence of regularisation and stabilisation. ε = γ = 0. A term for dissipation of motion orthogonal to the holonomic constraint surface can be added to improve the convergence. Also dissipation can be physically based by introducing it as a Legendre transform of a Rayleigh dissipation function R τ (q,q) = τ 2εġ Tġ → τλ Tġ − τ ε 2λ Tλ , with damping parameter τ . This modifies Eq. (5) to ελ + ετλ + g(q) + τ Gq = 0.
The compliance and damping factors ε, τ and γ are not restricted to being scalar or diagonal. In what follows these are assumed to be matrices.

Time discretisation
Variational integration [28] provides a systematic approach to derive time integration schemes with good properties, e.g. conservation of discrete symplectic forms and discrete momentum maps. Rather than discretising the Euler-Lagrange equations of motion directly, the Lagrangian and principle of least action are defined in time-discrete form.
Employing semi-implicit Euler discretisation and linearising the constraint as g(q + q) = g(q) + G q leads to the following stepping scheme [16,17] q n+1 = q n + hq n+1 , where p n = Mq n + h f n , v n = − 4 h Υ g + Υ Gq n and regularisation and stabilisation matrices . Equation (9) is a linear system of N = dim(q) + dim(g) + dim(ḡ) equations. The matrix on the left-hand side is blocksparse, positive definite, and nonsymmetric. The regularisation appearing as the diagonal perturbation matrices Σ andΣ is needed for handling otherwise ill-posed or ill-conditioned problems, e.g. systems with constraint degeneracy and large mass ratios. The stabilisation terms − 4 h Υ g + Υ Gq n on the right-hand side counteract constraint violations, e.g. sudden and large contact penetrations at impacts or small numerical constraint drift. The presented stepping scheme, referred to as spook, has been proved to be linearly stable [16], and numerical simulations suggest a large domain of nonlinear stability.
The combined effect of the regularisation and stabilisation terms is to bring elastic and viscous properties to motion orthogonal to the constraint surface g(q) = 0, e.g. for modelling elasticity in mechanical joints. The parameters ε and γ need not be chosen arbitrarily, as for the conventional approaches to constraint regularisation and stabilisation, but can be based on physics models using parameters that can be derived from first principles, found in the literature or be identified by experiments. This becomes straightforward when the regularisation is introduced by potential energy as quadratic functions in g, i.e. U ε (q) = 1 2 g T ε −1 g. This has been exploited in previous work to constraint-based mod-elling of lumped element beams [29], wires [30], meshfree fluids [31], and granular material [8].

Nonsmooth dynamics
In discrete time, some dynamics is best treated as nonsmooth [3], meaning that the velocity may change discontinuously in accordance with some impact law, expressed in terms of inequality constraints or complementarity conditions, in addition to the equations of motion. This is an efficient way of modelling dynamical contacts, dry friction, joint limits, electric and hydraulic circuit switching and actuator dynamics in discrete time with large step size. The alternative is to use fine enough temporal resolution where the dynamics appear smooth and may be modelled by daes or ordinary differential equations alone. In general, this is an intractable approach for full-system simulations and interactive applications. The nonsmooth dynamics can be formally treated as differential variational inequalities [32]. The discrete equations of motion (9) can be linearised into a mixed linear complementarity problem (mlcp), when nonsmooth dynamics is included where H z corresponds to the left-hand side of Eq. (9), r is the negated right-hand side, and w ± are slack variables. The terms u and l correspond to the upper and lower limits on the solution vector z, respectively. The original linear system of equations is recovered by assigning ±∞ to the upper and lower limits, i.e. no limit. Rigid body contacts are modelled by the Signorini-Coulomb law for unilateral dry frictional contacts [5], denoting the penetration depth by δ, contact normal and tangent by n and t. The Signorini-Coulomb law states that if the nonpenetration constraint, g n ≡ δ ≤ 0, is violated, then the normal contact velocity and constraint force are complementary to ensure separation, 0 ≤ġ n ⊥ λ n ≥ 0, and the tangential friction force that acts to maintain zero slip, g t ≡ t Tq = 0, is bound by the Coulomb friction cone, |λ t | ≤ μλ n . The latter may be linearised by approximating the cone with a polyhedral or a box. Impacts are treated post facto in a separate impact stage succeeding the update of velocities and positions. At the impact stage, an impulse transfer is applied enabling discontinuous velocity changes, i.e. fromq − toq + . Newtons' impact law, n Tq+ = −en Tq− , is used with restitution coefficient e in conjunction with preserving all other constraints, Gq = 0. This amounts to solving the same mlcp (10) but with p n = Mq n and v n = 0 on the right-hand side of Eq. (9).
A fixed timestep approach is preferred when aiming for fast full-system simulation with many nonsmooth events per unit time, which is the case with many dynamic contacts. The alternative, using variable timestep and exact event location, makes the numerical integration computationally intractable and may fail by occurrence of Zeno points.

Heterogeneous multidomain dynamics
In multidomain simulation with strongly coupled subsystems, e.g. an elastoplastic solid and mechatronic system, the dynamics propagate instantly throughout the full system. The nonsmooth multidomain dynamics approach provides a general framework for implicit integration of such systems that enable stable simulation at large timesteps. This avoid the loss in performance and numerical stability associated with explicit co-simulation algorithms. A system composed of two strongly coupled subsystems, A and B, takes the following form where , Σ AB are the multiplier, Jacobian, and regularisation of the subsystem coupling constraints. The individual subsystem matrices, H A and H B , take the generic form of a saddle point matrix as in Eq. (9).

Constraint-based meshfree elastoplastic solid
This section describes the elastoplastic model that is assumed, the spatial discretisation method and the compliant deformation constraint.

Elastoplastic solid
The solid is assumed to sustain geometrically large deformations. The St. Venant-Kirchhoff elasticity model is used in combination with the Drucker-Prager plasticity model [9]. The material strain is expressed by the Green-Lagrange strain tensor where u(x) is the displacement field mapping a reference coordinate x to displaced positionx(x) = x + u(x). Observe that the Green-Lagrange strain tensor transforms under large rotations and no co-rotation procedure is needed. It is convenient to use the Voigt notation for representing the strain tensor on vector format, = x x , yy , zz , 2 yz , 2 xz , 2 xy T , and the second Piola Kirchoff stress tensor σ = σ x x , σ yy , σ zz , σ yz , σ xz , σ xy T , related to the Cauchy stress by σ C = J −1 F T σ F with displacement gradient F = ∇u and J = detF. The linear constitutive law reads σ = C , with stiffness matrix where λ and μ are the first and second Lamé parameters which are related to the Young's modulus and Poisson's ratio, E and ν by λ The corresponding strain energy density is Plastic deformation occurs when the stress of a material reaches the critical yield strength, Φ(σ ) = 0. The choice of yield function, Φ(σ ), is based on the type of material being modelled. For simplicity, we restrict ourself to small plastic deformation and use the second Piola Kirchoff stress to evaluate the yield function and plastic increment. The strain tensor is additively decomposed into an elastic and a plastic component, = e + p , where the latter store the permanent plastic deformation. In ideal plasticity, the plastic deformation occurs instantly according to a flow rule d p = dλ p ∂Ψ ∂σ , where Ψ (σ ) is the plastic potential and λ p is the plastic multiplier. If Ψ (σ ) = Φ(σ ), the model is said to be associative and otherwise nonassociative. The plastic flow lasts as long as the plastic multiplier is positive dλ p > 0, incrementally reducing the stress dσ p , until it reaches the elastic regime, Φ(σ ) < 0. This constitutes a nonlinear complementarity problem known as Karush-Kuhn-Tucker conditions The plastic multiplier is computed from the constitutive law, which in the plastic flow phase is dσ = C p d , where the elastoplastic tangent stiffness matrix is Predictor-corrector algorithms are conventionally used to integrate the plastic flow. In the absence of plastic hardening or softening, the plastic deformation, plastic multiplier, and total stress can be computed easily using the radial return algorithm [9] summarised in Algorithm 2. The assumed plasticity model is a capped Drucker-Prager model, following Dolarevic [18], with a compaction variable κ. The yield function Φ (σ , κ) is C 1 continuous consisting of the Drucker-Prager yield surface with a tension and a compression cap function according to where I 1 = tr(σ ) is the first invariant of the stress tensor and J 2 = 0.5 tr(σ 2 ) is the second invariant of the deviatoric stress tensorσ = σ − 1 3 I 1 1. The expressions for the tension cap Φ T (I 1 , J 2 ) and compression cap Φ C (I 1 , J 2 , κ) are given in "Appendix", and the Drucker-Prager surface is defined as where φ is the internal friction angle and c the cohesion parameter such that The capped yield surface is illustrated in Fig. 2. The tension cap is fixed and is used to regularise the plastic flow behaviour in the corner region. The compression surface cap, on the other hand, is dynamic, and the maximum hydrostatic pressure, κ, is used as the main variable.
The conventional Drucker-Prager model is a common model for the plastic deformation dynamics of soils, e.g. wet or dry sand. These materials are weak under tensile stress (I 1 > 0) and become stronger under compression (I 1 < 0) where it may support large shear stresses ( √ J 2 = 0). The capped Drucker-Prager model is a generalisation that includes plastic compaction that occurs in many soils. The compaction mechanism is the failure of individual grains whereby air or water is displaced from the soil. The compaction saturates at a maximum level, where all voids vanish. The compaction dynamics is modelled by a variable cap on compressive side of the Drucker-Prager yield surface, intersecting the hydrostatic axis at −κ. The chosen compaction hardening law is [18] where W is the maximum volume compaction, D is the hardening rate, and κ 0 is the initial cap position where compressive failure first occurs. Observe that when the plastic volume compaction − tr( p ) approaches W , the cap variable κ goes to infinity. In this regime, the material no longer compact plastically and behaves like the standard Drucker-Prager model. The elastoplastic material parameters are summarised in Table 1. Expressions for the detailed shape of the compression cap and the derivatives of the yield functions are found in "Appendix".

Meshfree method
The continuous solid of mass m and reference volume V is discretised into N p particles. The particles have mass m i = m/N p , volume V i = V /N p , position q i , and velocityq i . The particle displacement is u i = q i − x i with reference position x i as illustrated in Fig. 1. A continuous and differentiable displacement field is approximated using the mls method [19,33]. The displacement gradient field defines the strain tensor field in any point by Eq. (12). The particles can thus be understood as pseudo-particles having both particle and field properties. The mls approximation of the displacement field is where the shape function and moment matrix use a quadratic basis and weight function W (x) = 315 64πl 9 (l 2 − x T x) 3 if |x| is smaller than the influence radius l, and otherwise zero. The base function notation p j = p(x j ) is used to simplify the expressions. The gradient of the interpolated displacement field is with The strain tensor field, (x), can thus be approximated by applying Eqs. (12)- (22). Observe that u j α depends on current particle positions q while ∇ β Ψ j (x) depends only on the reference positions.

Compliant deformation constraint
The strain energy density, U (x) = 1 2 T C , is transformed into a compliant constraint via a Legendre transform as presented in Sect.2. For each particle i, a deformation constraint is imposed with particle strain components i , q), computed by evaluating the field at x = x i . Note that the particle strain depends on the displacement of all particles within the domain of influence with radius l. The compliance parameter becomes ε i = (V i C i ) −1 , which depends on the material parameters, through Eq. (13), and on the spatial resolution through the particle volume factor V i that appears when integrating from energy density to particle energy. As expected, softer materials will have larger compliance parameter and thus deform more easily. The Jacobian of the deformation constraint, G = ∂ g/∂q, can be expanded through the chain rule into The Jacobian for constraint i has block structure where Λ i j Observe that only F i αβ (q) depends on the current particle positions while Λ i j β depends only on the reference positions and may thus be pre-computed. In the derivation of the Jacobians, the following partial derivatives are used with Kronecker notation δ ηκ γ τ ≡ δ γ τ δ ηκ .

Stabilisation constraint
The mls approximation has the deficiency of underestimating the strain for deformation modes on short length scales. This makes the method unstable, manifested by spurious short wavelength deformation modes. Following [20,34], we eliminate this instability and associated discretisation error by adding the following deformation energỹ where α is a dimensionless stabilisation parameter and f i are the forces on the pseudo-particle other than elastic ones. The corresponding stabilisation constraint becomesg i = (l/E) (∇ · σ ) x=x i with regularisationε i = (αV i E) −1 and after a suitable normalisation. The stress divergence is computed using a mls approximation of the stress tensor field The stabilisation constraint becomes where the stress and strain are represented using Voigt notation such that the divergence operator isD The Jacobian of the stabilisation constraint becomes where G jk τ σ = ∂ j τ /∂q k σ is already derived for the compliant constraint.

Contacts
Contacts between elastoplastic solids and rigid or static bodies are modelled by assigning a rigid spherical shape of radius l to each pseudo-particle and imposing the Signorini-Coulomb law for unilateral dry frictional contacts [5], as described in Sect. 2.3 and in greater detail in [8]. A contact situation between a solid and three rigid bodies is illustrated in Fig. 3. The contact overlap and normal are indicated as well as the contact point, which is defined at half distance between the overlapping geometries. Contacts between neighbouring pseudo-particles are ignored. Using spherical shapes to approximate the contact boundary of the solid has obvious drawbacks, e.g. perfectly smooth surfaces cannot be represented. This should be considered as an intermediate solution only, practical to implement in prototype code. More sophisticated representations of boundaries and contact detection algorithms for meshfree methods exist [35,36] and should replace the use of spherical shapes.

Simulation
The implementation details and results from nonsmooth multidomain dynamic simulation including elastoplastic solids is presented in this section. The main simulation algorithm is given in Sect. 4.1. The elastoplastic solid model is implemented in C++ making use of the simulation software library

Main algorithm
Algorithm 1 lists the main steps in running a multidomain dynamics simulation including both elastoplastic solids and contacting rigid multibodies. The elastoplastic solid is initialised by defining a solid volume and assigning material parameters. The solid is discretised into N p pseudo-particles according to a given spatial resolution. For simplicity, the particles are positioned in a regular cubic pattern. This defines the reference state with vanishing strain and stress. The elastoplastic constraints are initialised and connectivity data listing the particles involved in each constraint is created. Similarly, rigid bodies and kinematic constraints are defined, for instance, to form an articulated vehicle and powertrain. Each body is assigned a geometric 3D shape. Contact material properties are assigned to each body, including coefficient of restitution, elasticity, and friction. These are used to generate contact constraint data included in the mlcp when triggered by contact detection during simulation. Certain quantities are computed once before time integration and then reused for the sake of optimisation, e.g. the inverse moment matrix A −1 in the mls approximation, ∇ β Φ and Λ. The simulation is run with a fixed timestep, following the stepper in Sect. 2.2. Each timestep begins with reading input signals, e.g. from operator and control system, and computation of explicit forces. Next, contact detection algorithms produce a set of contacts for intersecting geometries. Each contact position and velocity, penetration depth, normal and tangent is stored. Contacts are classified as either impacting, continuous or separating, depending on the sign of the relative contact velocity. The strain and stress fields are computed as described in Sect. 3.1 using the mls approximation of the displacement field for the spatial discretisation described in Sect. 3.2. The trial stress is tested against one of the three yield surfaces decided by Algorithm 3. If it is outside the elastic domain, the radial return Algorithm 2 computes the plastic deformation that returns the stress to the yield surface with an error threshold ε tol . The compaction variable κ is also updated by this. Constraint violation and Jaco-

Overview of the MLCP solver
The main computational task for each timestep is solving the mlcp in Eq. (10) with saddle point matrix structure given in Eq. (9). Each sub-matrix is block-sparse and the routines for building, storing, and solving are tailored for this and exploit blas3 operations as much as possible to gain computational speed. In order to achieve high performance, e.g. real-time simulation of complex multidomain systems, splitting is applied [38]. Each subproblem can then be treated with different mlcp solvers that best fit the requirements of accuracy, stability, and scalability. For the elastoplastic solid and for jointed rigid multibodies, with relatively few complementarity conditions, a direct mlcp solver is used that is described in Sect. 4.2.2.

Split solve
Assume three sets of constraints, labelled A, B, and AB for two subsystem and their coupling. The linear system (11) is split into two subproblems by duplicating the variablesq and λ AB and reordering the augmented system. A stationary iterative update procedure can then be formed based on the four matrix blocks with and the same for H B A , z B A , and r B A with A and B permuted. The iterative procedure converges if the spectral radius of the augmented system H fulfils ρ( The key point is that each subproblem can be approached using different solvers. In the case of a machine interacting with an elastoplastic terrain, the machine and solid terrain constraints, A, are split from the tangential contact forces related to dry friction, B, both systems sharing the normal contact constraints (nonpenetration), AB. The subsystem (32), with the solid, machine, and contact normals is solved first using a direct solver, while the subsystem (33), containing friction constraints, is solved using an iterative projected Gauss-Seidel solver. The split solve may be terminated after the first iteration, accepting potential errors from the friction forces, or continued with a final stage 1 solve. Experiments show that further iterations do not necessarily reduce the errors. For large contact systems (>1000 contacts), with many complementarity conditions, both normal and friction constraints may be moved to an iterative projected Gauss-Seidel solver [8].

Direct solver
The direct solver is a block pivot method [38] for mlcps based on Newton-Rhapson iterations applied to nonsmooth formulation. The detailed algorithm is an adaptation of Murty's principal pivoting method, and is found in the reference [39] as Algorithm 21.2. Before the solve step, the matrix H is factorised by a permuted ldl-factorisation, H = P L DL T P T , where the permutation P is used to minimise fill-ins and may include leaf-swapping to avoid pivoting on Σ since it is often close to zero. The factorisation requires symmetric indefinite matrices, but this is not the case with H initially as shown in Eq. (11). This is fixed by absorbing a negative sign in the vector z.

Verification tests
The rotational invariance of the elastic solid model is verified in simulations confirming that the strain invariants are unaffected by rigid rotations of the body.

Elasticity
The validity of the model is confirmed by a number of tests of macroscopic relations between an applied load and the resulting displacement. Uniaxial stretching (u) and hydrostatic compression (c) are investigated, as described in Fig. 4. In the uniaxial test, the boundary condition in the tensile direction is a plane constraint with free slip in plane. The other boundaries are free. In the hydrostatic compression test, the boundary conditions are dynamically generated unilateral contact constraints. The boundary force is increased linearly with time, by regulating the position of the boundary geometry, while measuring the stress and strain components of Eqs. (36) and (37). All tests are simulated on a uniformly discretised homogeneous isotropic three-dimensional solid with parameters given in Table 2.
The simulation result is compared with exact solutions for from elasticity theory for uniform deformations of St. Venant-Kirchhoff materials. A test cube with initial rest length l 0 and cross-sectional area A 0 is considered. The deformed side length is denoted l, and the stretch ratio in that direction is λ = l/l 0 . The following exact relations can be derived where p = F/A 0 is the applied boundary pressure, and c u = E and c c = E (1−2ν) are the Young's and bulk modulus. The simulation results for E = 10 6 Pa and ν = 0.25 are presented in Fig. 5. The influence radius was chosen experimentally to be as low as possible, reducing the number of neighbours and the simulation time, while maintaining physically accurate results. The chosen influence radius results in 27 neighbours for any evaluation point that is at least half an influence radius away from the boundary.
The simulation result with stabilisation constraint (α = 100) agrees very well with the exact solution for large deformations. Without the stabilisation constraint, the simulated stress-strain relation deviates by 5-10% from the theory at large deformations |λ − 1| 0.15. The error decrease with finer resolution. Observe that for infinitesimal deformations 1 2 λ λ 2 − 1 ≈ l/l 0 , Eqs. (36) and (37) approximate to the well-known expressions from Hooke's law. Increasing the stiffness or resolution further requires a smaller timestep for numerical stability.
The strain field under hydrostatic compression is displayed in Fig. 6. With stabilisation constraint, the strain field is nearly uniform through the material. Without the stabilisation constraint, the strain deviates by as much as 40% with the largest deviations at the boundary.

Plasticity
The Drucker-Prager cap model is tested under hydrostatic and triaxial loading and unloading, following the tests made  Fig. 7. The result is a permanent deformation of the cube by tr = 0.1 compared to tr = 0.12 in [18]. The evolution of the stress invariants, I 1 and J 2 , in the yield space shows that the deviatoric stress is negligible and stress evolves purely along the hydrostatic axis during loading and unloading.
In the triaxial test, a hydrostatic pressure of σ h = 100 kPa is first established. The load along the z-axis is then increased up to a level where the material yields while the side wall pressure is kept constant. The material is finally unloaded. The simulated evolution of the deviatoric stress as function of strain is displayed in Fig. 8 for a Drucker-Prager material with and without cap model. The simulated Drucker-Prager material is found to yield at a critical stress of 270 kPa. This is in good agreement with the analytical prediction from Eq. (18). As expected, the yield plateau is lower and less pronounced in the case of the capped Drucker-Prager material, because the stress reaches the compression cap surface before the Drucker-Prager surface. The triaxial test can be used to determine the value of the plastic hardening parameter D, whereas the maximum compaction parameter, W , can be determined from the hydrostatic compression test.

Dynamic contacts
Dynamic contact handling is demonstrated by dropping a beam to rest on two thin cylinders and then deforming it by pressing a larger cylinder down on the beam. The result for an elastic and elastoplastic beam is displayed in Fig. 9.

Multidomain demonstration
Two example systems are simulated to demonstrate multidomain capability. The first system is an articulated terrain vehicle with tracked bogies driving over a deformable terrain. Images from simulation of a bogie and of a full vehicle are shown in Fig. 10. b Stress evolution for the capped Drucker-Prager model. c Stress evolution for the conventional Drucker-Prager model bogie frames, wheels and track elements. These are interconnected with kinematic constraints to model chassis articulation, bogie and wheel axes, and tracks covering the wheeled bogies. The vehicle drivetrain consists of an engine, torque converter, drive shafts, and differentials that distribute the engine power to the front and rear parts and further between left and right bogies to the wheels. The terrain is modelled which represents a weak and soft forest terrain. The tracked bogie and vehicle create a rutting with permanent deformation and the rut depth can be measured. In the full vehicle simulation, the solid is discretised into N p = 2100 pseudoparticles. The coupled system of vehicle and terrain thus have about 20.000 variables. With 10 ms timestep and using a single CPU on a conventional desktop computer, the simulation is of the order 1/100 from real time. It should be emphasised that this is measured on prototype code with no attempt for optimisation and parallelisation. Also, the implementation uses high-fidelity direct solver that delivers higher accuracy than what is typically required in terrain-vehicle simulations. Fig. 9 Deformation of an elastic and elastoplastic beam (middle and right from initial state to the left). The colour codes the Frobenius norm of the strain tensor, ranging from 0 (green) to 0.20 (red) Fig. 11 Image sequence from a simulated cone penetration test on an elastoplastic terrain with an embedded rock The second demonstration example is a dynamic cone penetration test, where a cylindrical weight is dropped repeatedly on a cone measuring the penetration depth. This is a common way of measuring the mechanical properties of ter-rain. The penetration depth of the cone in simulated for two different cases. The first case is a homogeneous material, and the second is material with an embedded rock, represented by a rigid body, as illustrated by Fig. 11. The colours indicate the magnitude of displacement. The simulated cone penetration is presented in Fig. 12. The resistance is higher in the ground with an embedded rock, before the cone actually comes in contact with the rock.

Discussion
A meshfree elastoplastic solid model is made compatible with nonsmooth multidomain dynamics. The solid appears as a particle system in a constrained multibody system, on the same footing as articulated rigid multibodies and power transmission systems. The pseudo-particles carry field variables, e.g. stress and strain tensor fields, approximated using the moving least squares method. The dynamic interaction between the deformable solid, rigid multibodies, and other geometric boundaries is modelled by unilateral contact constraints with dry friction. The full system can be processed using the same numerical integrator and solver framework without introducing additional coupling equations with unknown parameters. Impulsive behaviour can automatically be transmitted instantly through the full and strongly coupled system. This enables fast and stable simulation of complex mechatronic systems interacting with elastoplastic materials. The MLS approximation underestimates the strain for short wavelength deformation modes. A stabilisation constraint is developed to compensate for this and it improves the accuracy and stability significantly. The explicit form of the deformation constraint and its Jacobian is given in Eqs. (24)- (26), and the stabilisation constraint in Eqs. (29)- (31). The Jacobians can be factored to a constant matrix that can be pre-computed, multiplied with current particle displacements. The verification tests show good agreement with analytical reference solutions for large deformations of elastic and elastoplastic test cubes. The stabilisation constraint is essential for the accuracy and stability. Demonstration is made with a tracked terrain vehicle driving over deformable terrain using the capped Drucker-Prager model and a cone penetrometer test on terrain with and without an embedded rock.
The computational bottleneck of the simulation lies in solving the block-sparse mixed linear complementarity problem (10). With dedicated hardware, parallel factorisation algorithms, and iterative solver for a more approximate solution of the terrain dynamics, the performance can be expected to increase by several orders in magnitude. Exploring this is left for future work. Another area of improvement is to integrate the plastic yield and flow computations with the mixed complementarity problem for the constraint forces and velocity update. This will enable integration with larger timestep for strongly coupled problems than the current predictorcorrector method allow.
The presented method rest on the St. Venant-Kirchhoff elasticity model in combination with the Drucker-Prager plasticity model and the additive decomposition of the elastic and plastic strain. These assumptions are in general only valid under small deformations and alternatives should be considered in future work.