A multiscale model of terrain dynamics for real-time earthmoving simulation
SServin et al.
RESEARCH
A multiscale model of terrain dynamics forreal-time earthmoving simulation
Martin Servin , Tomas Berglund and Samuel Nystedt * Correspondence:[email protected] Department of Physics, Ume˚aUniversity, SE-90187 Ume˚a,SwedenFull list of author information isavailable at the end of the article
Abstract
A multiscale model for real-time simulation of terrain dynamics is explored. Torepresent the dynamics on different scales the model combines the description ofsoil as a continuous solid, as distinct particles and as rigid multibodies. Themodels are dynamically coupled to each other and to the earthmoving equipment.Agitated soil is represented by a hybrid of contacting particles and continuumsolid, with the moving equipment and resting soil as geometric boundaries. Eachzone of active soil is aggregated into distinct bodies, with the proper mass,momentum and frictional-cohesive properties, which constrain the equipment’smultibody dynamics. The particle model parameters are pre-calibrated to the bulkmechanical parameters for a wide range of different soils. The result is acomputationally efficient model for earthmoving operations that resolve themotion of the soil, using a fast iterative solver, and provide realistic forces anddynamic for the equipment, using a direct solver for high numerical precision.Numerical simulations of excavation and bulldozing operations are performed tovalidate the model and measure the computational performance. Reference datais produced using coupled discrete element and multibody dynamics simulationsat relatively high resolution. The digging resistance and soil displacements withthe real-time multiscale model agree with the reference model up to 10-25%, andrun more than three orders of magnitude faster.
Keywords: deformable terrain; discrete element method; multibody dynamics;multiscale; real-time simulation; soil mechanics
Introduction
Physics-based simulation of earthmoving equipment and soil is an important toolfor developing smarter systems that meet the increasing demands for energy ef-ficiency, productivity and safety in the agriculture, construction, and mining in-dustries. Simulation with real-time performance is essential when developing newcontrol systems, human-machine interfaces or training operators using interactivesimulators with hardware or human in the loop. Fast simulation is vital also forapplying artificial intelligence to motion planning and control of earthmoving equip-ment. These techniques are data-hungry, requiring many repeated simulations of awide range of normal operating conditions as well as many, potentially hazardous,edge cases. Although the simulations can be run mostly in parallel, real-time per-formance, or faster, is necessary for covering a sufficiently large parameter spacewithin reasonable time. Simulating heavy equipment alone at real-time with goodaccuracy is challenging but feasible using rigid multibody dynamics and a well-optimized solver [1, 2]. It may seem out of reach to include also the environment, a r X i v : . [ c s . C E ] N ov ervin et al. Page 2 of 35 involving vast amount of soil with complex dynamics on scales that cross six ordersof magnitude [3].The first physics-based models for real-time simulation of heavy equipment anddeformable terrain appeared in mid 1990, starting with Li and Moshell assuming theMohr-Coulomb theory and volume preserving deformations [4]. Park [5] developedan elaborate model of the digging resistance in soil for the purpose of constructionexcavator simulators. The model starts from the fundamental earthmoving equation(FEE) [6] for a blade cutting a horizontal soil bed, assuming an active zone in theshape of a wedge [7]. The interface between the wedge and the passive soil representthe failure surface, stretching from the cutting edge of the blade to the free surfaceof the terrain. Park extends the model to a 3D bucket digging in sloped terrain,including the formation of a secondary separation plate by the deadload of materialin the bucket and penetration resistance from the tool’s teeth. The FEE-basedmodels take into account the strength of the soil, expressed in terms of its internalfriction and cohesion, but are limited to stationary conditions and do not describethe motion of the soil. Holz et al [8, 9, 10] combined the fundamental earthmovingequation (FEE) with contacting particle dynamics to model the cutting resistanceand motion of the soil. In this approach, static soil is adaptively converted intoparticles as the tool cut the terrain. Digging resistance, based on the FEE, is appliedas a six-degree of freedom kinematic constraint with force limits on the digging tool.A wedge-shaped active zone is assumed. The portion of terrain in the active zonethat has not yet undergone failure and particle conversion provide mass to theFEE. The particle contacts on top of the wedge contribute as surcharge mass, andthe contacts with the tool provide additional digging resistance. Later, a similarapproach is taken in [11].Unfortunately, relying on the FEE for the digging resistance has serious draw-backs. It is valid only at steady state and it suffers from unphysical singularities,as pointed out in [10]. Furthermore, the active zone possesses no inertia or momen-tum. When accelerated motion is involved, the FEE can clearly not provide correctdynamics and proper reaction forces. On the other hand, resolving the active zonewith contacting particles challenges the computational performance. If the particlesare not too many, real-time simulation is possible using nonsmooth contact dynam-ics [12] or position-based dynamics [13], with a stationary iterative solver like theGauss-Seidel algorithm. This is associated with numerical errors that manifest asartificial elasticity or insufficient friction [14, 15]. If the errors become too large thesoil will behave more like a compressive fluid rather than a stiff soil that yield onlyif the shear stress reach critical values The errors grow with the number of particlesand applied stress and decrease with the time-step and increased number of solveriterations. In the other limit, of a few coarse particles, the solver error decrease butthe spatial discretization errors grow large.To cope efficiently with the disparate length and time scales, we explore a multi-scale model for the soil dynamics. The model targets real-time simulation of earth-moving equipment and terrain with realistic reaction forces and soil deformations.A macroscale model , describing the rigid multibody dynamics, is simulated using adirect solver for high numeric precision. The active soil, resolved as particles in a mesoscale model , is simulated using an iterative solver at high-speed with large error ervin et al.
Page 3 of 35 tolerance. The key idea is that the soil dynamics is represented in both models, witha coupling that filters out the discretization and solver errors from the mesoscalemodel but capture the bulk dynamics. When an earthmoving object come in contactwith the terrain, the potential failure surfaces and zones of active soil are predicted.The soil inside the active zones is represented using a hybrid model of contact-ing particles and continuum solid, which support smooth transitioning between theresting solid, liquid and gas phases. The coupling back to the macroscale modelis mediated through aggregate bodies . They have the momentary shape, mass, andvelocity of the resolved soil in each active zone. The motion of the aggregates areconstrained relative to the equipment and the terrain failure surface, in accordancewith the Mohr-Coulomb model. For blunt objects contacting the terrain, the sub-soil stress distribution is estimated and the soil compacts if the stress reach criticalvalues. A microscale model , simulated using relatively small particles and high nu-merical precision, is used for pre-calibration of the particle parameters to match thebulk mechanical properties of a wide range of soils. The primary model parametersof the multiscale model are the bulk mechanical properties of the soil at a givenbank state: the internal friction, cohesion, dilatancy, elasticity, and mass density;as well as the equipment’s geometry, surface friction and cohesion. To test and vali-date the model, simulations of excavation and bulldozing operations in various soilsare performed. The digging resistance and soil displacements from the real-timemultiscale model are compared with those from the microscale reference model.
Modelling and simulation of soil and earthmoving equipment
Length- and timescales
Soil and granular media are strongly dissipative multiphase materials with multiplelength- and timescales [3, 16]. They consist of contacting grains with size rangingfrom clay at 10 − m to cobble and boulders at 10 m. Natural soil has a certainmoisture content, that may significantly increase the strength of the soil by cohesiveforces. At sufficient moisture levels, however, a pore pressure develops that lower theinter-particle normal forces and, consequently, the internal friction. Large deforma-tions are commonly localized in shear bands, sometimes as narrow as a few particlediameters. The soil outside the shear bands is displaced rigidly. Soil is strongly dis-sipative. Therefore, the solid phase is the natural state. If it is agitated, for instanceby an earthmoving equipment, it may transition to the liquid or gaseous phases.For sand and gravel in typical earthmoving operating conditions, the grain collisiontimescale is in the microsecond regime while the liquid timescale is around a mil-lisecond [1] . Earthmoving equipment, on the other hand, has characteristic size of10 m, operating range of around 100 m and may displace several cubic meters of soilper second. An earthmoving tool can be controlled with a spatial precision of about10 mm and its geometric features is also on this scale. A typical loading cycle, with [1] The liquid phase time scale is the characteristic time scale of particle rear-rangements, d/ (cid:112) p/ρ , with particle diameter d , mean stress p , and specific massdensity ρ . The time scale in the gaseous state is the characteristic collision time v − (5 m/ k ) ] / , assuming a strongly damped Hertzian contact model [17], for im-pact velocity v , particle mass m , and contact elasticity k = √ dE/ − ν ) , with Young’smodulus E and Poisson’s ratio ν . ervin et al. Page 4 of 35 a bucket excavator or a wheel loader, has time duration in the range between 15-30s [18, 19]. The natural timescales of the rigid body motion is about 1-10 Hz and thecontrol systems typically operate at a frequency below 100 Hz [20]. These multipleand separated length- and timescales are important to consider when modelling soiland granular media interacting with earthmoving equipment.
Distinct particles
The discrete element method (DEM) [21] describe granular media and soil as con-sisting of distinct contacting particles of finite size. It is a versatile model thatautomatically describe soil in the solid, liquid, and gaseous phases and the transi-tions between them. DEM is clearly applicable for simulating the soil dynamics inearthmoving operations [22, 23, 24]. It is, however, very computationally intensive.A common solution is to not represent the true grains with their actual distribu-tion of size, shape, and mechanical properties. Instead, the soil is represented by acollection of large, often spherical, pseudo particles with contact parameters - elas-ticity, friction, cohesion and rolling resistance - that are calibrated such that thebulk mechanical properties match the ones observed in the soil that the model ismeant to represent. The mapping between the particle parameters and macroscopicsoil parameters can be carried out using the triaxial test or cone penetration test, asdescribed in [25]. But even with pseudo-particles of size 10 to 100 mm, the numberof particles for representing a terrain the size of earthmoving equipment exceedswhat is currently possible to simulate in real-time using DEM.
Continuum
On length scales larger than the particle size, soil may be modelled using continuummechanics. The state of the soil is represented with scalar, vector and tensor fieldsfor mass density ρ ( x , t ), displacement u ( x , t ), velocity v ( x , t ), stress σ ( x , t ) andstrain (cid:15) ( x , t ). The fields obey the equations of mass continuity ∂ t ρ + ∇ · ( ρ v ) = 0 . (1)and momentum balance[ ∂ t + v · ∇ ] v = ρ − ∇ · σ + f ext , (2)where we use the short notation for partial time derivative ∂ t = ∂∂t , and f ext denoteany external force, like gravity. In the dense regime, and when the stress is belowa certain yield strength condition, the soil behaves as an elastic solid with someconstitutive law relating stress to strain, e.g., Hooke’s law for small deformations σ ( x , t ) = C : (cid:15) ( x , t ). For isotropic and homogenous materials the stiffness tensor C has only two independent parameters, often represented by the Young’s modulus E and the Poisson’s ratio ν . When the stresses reach the yield condition the solidfails and deforms plastically, rupture or flows rapidly. The simplest yield conditionfor soil is the Mohr-Coulomb criteria. It predicts that the material will fail alongany plane with normal n where the shear stress, τ n = (cid:112) ( σ · n ) − σ n , reaches thecritial value of τ n = µσ n + c, (3) ervin et al. Page 5 of 35 where the normal stress is σ n = σ αβ n β . The model has two parameters for thestrength of the material: the internal friction, µ , and the cohesion, c . The internalfriction is often represented by the angle of internal friction φ = arctan( µ ). Anal-ogously with Coulomb friction, the critical shear strength grows linearly with thenormal stress (pressure). The shear stress must also overcome any cohesive strengthof the material, which is independent of the stress and strain in the Mohr-Coulombmodel. If the material yields quasistatically it may be modelled as an elastoplasticsolid, with a plastic flow rule that accounts for strain-hardening (or softening) andthe volume expansion during shear. The latter is known as dilatancy and is definedas tr( ˙ (cid:15) ) / (cid:107) ˙¯ (cid:15) (cid:107) , where ¯ (cid:15) is the deviatoric strain. The hardening can be incorporatedin the Mohr-Coulomb law by [16, 26] by an effective internal friction µ = tan( φ + ψ ) (4)where ψ = arcsin [tr( (cid:15) ) / (cid:107) ¯ (cid:15) (cid:107) ] is the dilatancy angle. Soil may compact plasticallyunder uniform compression if the applied stress is large enough to cause particlerearrangements. This is captured by a soil’s compression index C c ≡ ∆ e ∆ ln( σ ) (5)which measure the change in void ratio, e , by a change in confining stress.The finite element method is the most widely used and versatile technique for sim-ulating deformable solids. The popular integration schemes for elastoplastic solids[27] first compute trial strain and stress fields, and then use Newton’s method forsearching the plastic strain increment that fulfill the constitutive model and plasticflow rule at each time-step. In the regime of loose soil and large shear rate, the ma-terial is better described as a viscoplastic solid or a non-Newtonian fluid with someconstitutive law between the stress and strain-rate, like the µ ( I )-rheology for gran-ular media [28]. Cemented soil, with clay or silt particles, may fracture by brittlefailure. Mesh-based numerical methods for solid dynamics have various difficultieswith large deformations and topological changes associated with soil that undergoplastic or brittle failure, and transitions between the dense solid, dilute liquid andgaseous phases. Meshfree methods, such as the material point method, are showingpromising results on 3D problems, coupled with multibody dynamics, but real-timesimulation of earthmoving operations remain out of reach [29, 30].Analytical and semi-analytical solutions, derived from the continuum theory, arevaluable for both insight and for creating fast simulation models. This include theBoussinesq type-of-equations for the stress and deformation fields underneath aload applied on the surface of a semi-infinite domain. The reaction forces on a thinblade cutting a soil bed has been thoroughly analysed [7]. A blade, or separatingplate , has two basic modes of operation, penetration and separation . Penetration isthe motion straight into the soil with relative velocity in tangential direction of theplate only. Soil cavities at the tip of the blade, often equipped with teeth, are thusforced to expand [5, 31]. The penetration resistance may be considerable althoughthe deformations are relatively small. Separation corresponds to movement normalto the plate and is the main cause of soil failure and large displacements. The edge ervin et al. Page 6 of 35 direction of motion failure surfacesoil wedgesurchargesoil floor direction of motion frictionsurcharge weightwedge weight cohesionadhesion ρρ soil resistanceseparationpenetration direction of motion failure surfacesoil wedgesurcharge soil floor direction of motion internal frictionsurcharge weightwedge weight cohesionsurfacefrictionadhesion θθ soil resistanceseparationpenetrationtool soil surface cutting edge d β Figure 1
Illustration of a blade interacting with a soil bed. There are two modes of operation,penetration, and separation (left). The failure surface form a wedge-shaped active zone (middle).When a blade cut the soil there are many forces contributing to the soil resistance (right). where the blade meet the material is referred to as the cutting edge . See Fig. 1 Theshape of the failure surface can be computed analytically in the two dimensionalcase, applicable for a wide blade, using the method of stress characteristics andassuming the Mohr-Coulomb criteria. The failure surface is often approximated witha plane. This defines an active zone with the shape of a wedge. Rankine’s theoryfor a flat soil with a blade pushing on it in the horizontal direction predicts thatthe soil fails at an angle θ = π/ − φ/ L , in the FEE is composed of four terms FL = ρgd N γ + QN Q + cdN c + c a dN a , (6)with specific soil mass density ρ , tool penetration depth d , soil cohesion c , surchargeforce Q (per tool width) and soil-tool adhesion c a . The first term is the due to theweight of the wedge, the second term is the additional (vertical) surcharge, the thirdterm the cohesive force in the failure surface, and the fourth term is the resistancedue the adhesion between the blade and the soil, see Fig. 1. The four N -factors(found in [6, 7]) depend on the geometry of the tool and the failure zone, the internalfriction and cohesion, as well as of the surface friction and adhesion. The quadraticdependency on the cutting depth d reflect that the weight depends linearly on thecross-section area of the failure zone. Note that the cohesive and adhesive forceterms are proportional to the area of the failure surface and blade contact surface,respectively. One key limitation of the FEE is that it assumes steady state and lowspeed of the blade and soil flow. Also, the N -factors suffer from singularities atcertain geoemtric configurations. Contacting rigid multibodies
The earthmoving equipment can be simulated efficiently using rigid multibody dy-namics [1, 2]. Articulated and actuated mechanisms are thus modelled as rigidbodies with kinematic constraints that represent the mechanical joints, motors, anddriveline. Dynamic contacts between the bodies are best modelled as kinematicconstraints and complementarity conditions to express unilaterality, impacts and ervin et al.
Page 7 of 35
Coulomb friction. In the current paper, we apply rigid multibody dynamics withcontact dynamics for both the equipment and the soil, and in all the levels of themultiscale model. Therefore, we describe the computational framework at greaterlevel of detail.The state of a rigid multibody system with N b bodies, N j joints and actuators and N c contacts, is represented on descriptor form in terms of the position, x ( t ) ∈ R N b ,velocity, v ( t ) ∈ R N b , and Lagrange multipliers, λ j ( t ) ∈ R N j and λ c ( t ) ∈ R N c ,that are responsible for the constraint forces in joints and contacts. The systemposition variable is a concatenation of the spatial and rotational coordinates ofthe N b bodies, x = [ x , e ], and the velocity variable holds the linear and angularvelocities, v = [ v , ω ]. The time evolution of the system state variables [ x , v , λ ] isgiven by the following set of equations M ˙ v = f ext + G Tj λ j + G Tc λ c (7) ε j λ j + η j g j + τ j G j v = u j , (8) contact law ( v , λ c , g c , G c ) , (9)where M ∈ R N b × N b is the system mass matrix, f ext is the external force, and G Tj λ j and G Tc λ c are constraint forces for joints (and motors) and for contacts,respectively. The forces have dimension R N b and is composed of linear force andtorque. Eq. (8) is a generic constraint equation, with constraint function g j ( x ),Jacobian G = ∂ g /∂ x , compliance ε j and viscous damping rate τ j . An ideal jointis represented with ε j = τ j = u j = 0, in which case Eq. (8) express a holonomicconstraint, g j ( x ) = 0. A linear or angular motor may be represented by a velocityconstraint (setting ε j = η j = 0 and τ j = 1), G j v = u j ( t ), with set speed u j ( t ). Theholonomic and nonholonomic constraints can be seen as the limit of a stiff potential, U ε = ε g T g , or a Rayleigh dissipation function, R τ = τ ( Gv ) T Gv , respectively[32]. This offer the possibility of mapping known models of viscoelasticity to thecompliant constraints. Descriptor form means that no coordinate reduction is made.The system is represented explicitly with its full degrees of freedom, although withthe presence of constraints.We consider the system to have nonsmooth dynamics [33]. That means that thevelocity and Lagrange multipliers are allowed to be time-discontinuous, reflectinginstantaneous changes from impacts, frictional stick-slip transitions or joints andactuators reaching their limits. This is unavoidable when using an implicit integra-tion scheme [2] because of the coupling between the state variables trough unilateraland frictional contacts.As contact law between particles we use a model that include cohesive-viscoelasticnormal contacts (n), tangential Coulomb friction (t) and rolling resistance (r). Theseare formulated in terms of inequality and complementarity conditions for the ve-locities, Lagrange multipliers and constraint functions. The resulting model can beseen as a time-implicit version of conventional DEM and is therefore referred to as [2] The alternative is to resolve the contact events using smooth trajectories, stiffpotentials and small time-step explicit time integration. In the limit of high stiffnessand small mass, the simulation time increase indefinitely with this approach. ervin et al.
Page 8 of 35 nonsmooth DEM [34, 14]. We use the following conditions on v and λ c = [ λ n , λ t , λ r ]as contact law :0 ≤ ε n λ n + g n + τ n G n v ⊥ ( λ n + f c ) ≥ , f c ≡ c p A p / (cid:107) G Tn (cid:107) (10) γ t λ t + G t v = 0 , (cid:107) λ t (cid:107) ≤ µ t (cid:107) G Tn λ n (cid:107) (11) γ r λ r + G r v = 0 , (cid:107) λ r (cid:107) ≤ rµ r (cid:107) G Tn λ n (cid:107) , (12)where g n is a function of the contact overlap and the Jacobians, G n , G t and G r ,are the normal, tangent and rotational directions of the contact forces [14]. Theparameters ε n , τ n , γ t in Eq. (10) control the contact stiffness and damping, and f c the cohesion. Setting these parameters to zero means that no penetration shouldoccur, g n ( x ) ≥
0, and if so the normal force should be repulsive, λ n ≥
0. Theinclusion of f c enables cohesive normal force with maximum value c P A P , where c P is the particle cohesion and A P is the particle cross section area. The cohesion isactive when the contact overlap is smaller than a certain cohesive overlap , set to afraction of the particle size, e.g., δ c = 0 . d . Eq. (11) state that contacts shouldhave zero slide velocity, G t v = 0, giving rise to a friction force that is bounded bythe Coulomb friction law with friction coefficient µ t . Similarly, Eq. (12) states that,as long as the constraint torque is no greater than the rolling resistance law, thecontacting bodies are constrained to zero relative rotational motion, G r v = 0. Here, µ r is the rolling resistance coefficient and r is the particle radius. It is a well-knownfact that the effect of particle angularity, on internal friction and angle of repose, canbe captured using spherical particles with rolling resistance. As explained in [35], a n -sided polygon can be assigned a rolling resistance coefficient µ r = (1 /
4) tan( π/ n ),which gives µ r = 0 .
05 for an eight-sided polygon, µ r = 0 . µ r = 0for a sphere ( n = ∞ ). We map the normal contact law, Eq. (10), to the Hertz-Mindlin model for contacting viscoelastic spheres, f n = k n δ / n + k n c d δ / ˙ δ n ,where δ ( x ) is the contact overlap, k n = E ∗ √ d ∗ is the contact stiffness, c d is adamping coefficient, E ∗ = [(1 − ν a ) /E a + (1 − ν b ) /E b ] − is the effective Young’smodulus, d ∗ = ( d − a + d − b ) − is the effective diameter for two contacting spheres, a and b , with Young’s modulus E a , diameter d a and Poisson’s ratio ν a etc. Themapping to Eq. (10) is accomplished by g n = δ / , ε n = 5 / k n and τ n = 5 c d / v imp . If the relative contact velocity is smaller than this value the contactsare modelled as described above. In case of impacts we apply the Newton impactlaw G n v + = − e G n v − with restitution coefficient e , while preserving all other con-straints in the system on the velocity level, Gv + = 0. This is carried out in animpact stage solve, prior to the main solve for the constrained equations of motions(7)-(9). With this division, the restitution coefficient become the key parameter formodelling the dissipative part of the normal force. For the resting contacts we cansimply enforce numerical stability using τ n = 2∆ t with little consequence of thedamping being artificially strong [36].For numerical integration we employ the SPOOK stepper [32]. It is a first order ac-curate discrete variational integrator, developed particularly for fixed time-step real-time simulation of multibody systems with non-ideal constraints and non-smoothdynamics. It has been proven to be linearly stable. The numerical time integration ervin et al. Page 9 of 35 scheme for advancing the system’s position and velocity from [ x n , v n ] at time t n to[ x n +1 , v n +1 ] at time t n +1 = t n +∆ t consist of a position update x n +1 = x n +∆ t v n +1 after having computed the new velocity x n +1 and corresponding Lagrange multiplier λ . This is done by solving the following mixed complementarity problem (MCP)[37] M − G n − ¯ G n G n Σ G n Σ v n +1 λ ¯ λ − p n q n u n = w l − w u , (13)0 ≤ z − l ⊥ w l ≥ ≤ u − z ⊥ w u ≥ p n = M v n + ∆ t M − f n , q n = − t Υ g + Υ Gv n ,and u n is the target speeds of the velocity constraints (zero for frictional con-tacts). The regularization and stabilization terms are Σ = t diag[ ε i / (1 + 4 τ / ∆ t )],¯ Σ = t diag ( γ ) and Υ = diag[1 / (1 + 4 τ / ∆ t )]. The slack variables w l and w u areonly used internally by the MCP solver. For details about the solver, see Sec. Im-plementation. Coarse-graining
The process of averaging particle kinematics and contacts forces into continuousand differentiable fields is referred to as homogenization.
Coarse-graining is oneparticular way of doing homogenization. This is useful when combining particleand continuum models of granular media. The fields are computed by sampling theparticle variables using a smoothing kernel, ζ ( x ), that is normalized (cid:82) ζ ( x )d x = 1and approach zero on a smoothing length R . The fields of mass and momentumdensity are computed ρ ( x , t ) = (cid:80) a m a ζ ( x − x a ( t )) and p ( x , t ) = (cid:80) a m a v a ζ ( x − x a ( t )), respectively, and the velocity field is simply u ( x , t ) = p ( x , t ) /ρ ( x , t ) . (16)The rate of strain tensor can thus be computed as ˙ ε αβ ( x , t ) ≡ (cid:16) ∂u α ∂x β + ∂u α ∂x β (cid:17) . Thestress tensor is the sum of the kinetic stress σ k αβ ( x , t ) = − (cid:80) a m a u aα v aβ φ ( x − x a ( t ))and the contact stress σ c αβ ( x , t ) = − (cid:80) k f abα,k x abβ (cid:82) φ ( x − x a ( t ) + s x ab ( t ))d s , wherethe summation is over the set of contacts, f abk is the contact force between particle a and b with branch vector x ab = x b − x a . Different smoothing kernels can be usedfor different purposes. The Gaussian kernel, ζ ( x ) = ( √ πR ) − exp( −(cid:107) x (cid:107) / R ),make the fields differentiable. The Heaviside function is faster to evaluate and canbe used for a discrete representation of the fields that preserve mass precisely. Fromanalytical studies and numerical experiments with contacting elastic spheres [16], itis found that the elastic bulk modulus is non-linear and change with bulk volume as ∂σ/∂ε = ZE (cid:112) ∆ V /V , where Z is the average number of contacts per grain. Thissuggest an effective Young’s modulus of the form E = E (cid:16) ± (cid:107) ρ/ρ − (cid:107) / (cid:17) , (17) ervin et al. Page 10 of 35 for a small change in density ρ/ρ = 1 + ∆ V /V relative to a reference state withmass density ρ and Young’s modulus E . The sign ± is positive for compactionand negative for expansion. The dilatancy angle also increase with the level of com-paction, and with it the internal friction by Eq. (4). Based on numerical simulationsand coarse-graining of dense granular media, Roux and Radjai [38] proposed ψ = c ϕ ( ϕ − ϕ c ) , (18)where ϕ and ϕ c are the current and critical porosity, at which the soil switch betweenpositive dilation (volume expansion) and negative dilation (volume shrinkage) uponshearing, and c ϕ is a constant that depends on the particle shape. A multiscale model of terrain dynamics
The multiscale model has three levels of abstraction, illustrated in Fig. 2, that werefer to as micro- , meso- and macroscale . In the microscopic model the terrainis fully resolved by relatively fine-grained particles with contact properties pre-calibrated to represent various type of soil, e.g., dirt, gravel, or sand. It serves asa ground-truth reference model, simulated off-line in advance with relatively smalltime-step and high numerical accuracy. The output is used for validation of thecoarser scale models and, if necessary, for calibration of parameters not known fromtheory or experiments. The meso- and macroscale models are simulated in real-timeand synchronously coupled to each other. The terrain deformations and soil flow isintegrated using the mesoscale model, with a hybrid representation of the soil thatcombine coarse particles and fields discretised with a regular grid of voxels. Inputto the mesoscale model is the motion and contact forces of the equipment at theinterfaces to the terrain. This is provided by the macroscale model, which focus onthe rigid multibody dynamics of the equipment, and any other objects interactingwith the terrain. The equipment experience the resting terrain as a quasistaticsurface and each region of agitated soil as distinct dynamic bodies, whose shape,mass velocity and mechanical strength is aggregated from the mesoscale model.It is a significant feature that each of the three models use the same parametersto characterize the soil properties. These are mapped to the model parameters forthe particle, voxel, and aggregate body dynamics. The parameter mapping relies inpart on the theory of continuum soil mechanics and in part on parameter calibrationusing the microscale model. The primary soil parameters include mass density ρ b ,internal friction angle φ b , cohesion c b , dilatancy angle ψ b , bulk elasticity modulus E b , that characterise the mechanics of the soil as found in a natural bank state .These are collected in a bulk parameter vector p b = [ ρ b , φ b , c b , ψ b , E b ]. The bankstate refers to the state at which the true soil is found in nature. It may be theresult of geological, meteorological and machined processes, that leave the soil ata particular packing density and moisture content. Two soils with identical (real)particles may thus have different bank state properties and are considered as twodistinct soil with different bulk parameters p b . Microscopic model
At highest resolution we represent the terrain as a set of relatively fine-grainfrictional-cohesive particles, N P , and the equipment as a set of rigid multibodies, ervin et al. Page 11 of 35 timelengthmicromesomacro gas solid liquid coupling granular phases σ σ σ σ calibration testshigh-resolution reference model rigid multibodiescoarse pseudo-particles resting solid aggregate body AC active zone failure surface fluidized solid compressed solid D B BABBDAC BCBC
Figure 2
Overview if the multiscale model. The microscale model is simulated off-line forvalidation and calibration of the mesoscale and macroscale models, that are simulated in real-timecoupled to each other. The mesoscale model combines a coarse particle and continuumrepresentation of the soil inside the active zone. The active soil is aggregated into a single bodywith frictional-cohesive couplings with the resting terrain and the rigid multibody equipment. N RB , with some set of joints and actuators. Let N BP denote the set of particles indirect contact with a body B ∈ N RB of the equipment. The equations of motion forthe earthmoving body and the terrain particles, labelled a , are m a ˙ v a = f ext a + G T a B λ a B + (cid:88) a (cid:48) ∈N aP G T aa (cid:48) λ aa (cid:48) , (19) m B ˙ v B = f ext B + G T B λ B + (cid:88) a ∈N BP G T a B λ a B , (20)where G T B λ B is the constraint force coupling the earthmoving body to the restof the equipment, G T B a λ a B is the contact force on the body from particle a , and G T aa (cid:48) λ aa (cid:48) is the inter-particle contact force between particle a and a (cid:48) . The body-particle interfacial force is the source for the digging resistance perceived by theequipment. The interfacial stress also alters also the internal stress in the soil,causing shear failure or brittle rupture if it reach the critical stress.The particles can roughly be divided into two domains. Either they are part ofthe active zone, N A P , which is sheared or rigidly displaced by the earthmoving body,or they remain part of the resting bed of particles, N C P . The domains are separated ervin et al. Page 12 of 35 by a failure surface AC . If the earthmoving body is a tool with penetrating teethor a sharp edge, the particles in direct contact with that constitute a third domain, N D P . The particle domains are illustrated in Fig. 2.We use a pseudo-particle representation of the soil. That means that the particlesdo not represent the true grains with their actual distribution of size, shape, and me-chanical properties. Instead, the soil is represented by a collection of relatively largespherical particles with specific mass density and contact parameters - elasticity,friction, cohesion and rolling resistance - p p = [ ρ p , µ t , µ r , c p , E p ]. These parametersare calibrated to numerical values that give the particle soil the desired bulk me-chanical properties p b . The size of the pseudo-particles is chosen small enough toresolve the important features of the equipment and precision at which it can becontrolled, i.e., around 10 −
50 mm in earthmoving applications.
Mesoscale model
The mesoscale model is a medium-resolution model of the soil dynamics using ahybrid particle-continuum approach. The soil in the active zone is primarily repre-sented by coarse particles. The continuous soil model has two phases, resting solidmass and fluidized mass . The former represents resting soil outside the active zoneand is considered an elastoplastic solid. The latter complement the use of particlesfor representing soil displacement in the active zone. The fluidized mass is convectedwith the coarse-grained velocity field of the particles, and subject to gravity. Themacroscale model supplies the equipment’s motion and contact forces at the terraininterface as input to mesoscale model. This is the basis for predicting the activezones and provide boundary conditions that drives the mesoscale soil dynamics.A regular grid is used for the discrete representation of the continuum model. Thegrid cells, or voxels , are labelled i = ( i, j, k ) by the triplet of integer positions in thegrid, aligned with the global coordinate axes. Each voxel has a centre point x i =[ x i , y i , z i ] and constant volume V = l . Each voxel has a variable mass m i , velocity v i , compaction w i and occupancy ϕ i . These are collected in a state vector s i =[ m i , v i , w i , ϕ i ]. The evolution of s i is treated by a cellular automata, with transitionrules that implements convection of fluidized mass, plastic compaction at criticalsub-soil stress, conversion between the particle and continuum representation, andsurface flow if the local slope exceeds the soil’s angle of repose δ b [39, 4]. Forcohesion-free materials it coincides well with the internal friction angle, given byEq. (4). With cohesion it may be larger [40]. Beyond this slope the terrain is not ina stable equilibrium, and will fail by avalanching into a valid state. The transitionsare constructed to preserve the total mass to machine precision.The voxel mass density, ρ i = ρ s i + ρ f i , is composed by the density of resting solid, ρ s i , and fluidized mass, ρ f i . All voxels in the terrain are fully occupied, ϕ i = 1,with solid mass and empty of fluidized mass, except for surface voxels, that mayhave occupancy less than unity, ϕ i = V s i /V ∈ [0 , ρ b but can vary locally within the range ρ i s ∈ [ ρ min , ρ max ]if subject to compaction or swelling, w i ≡ ρ i s /ρ b . The amount of solid mass in avoxel i is consequently m i s = ρ i s V i s = w i ϕ i ρ b V . The upper and lower limit on themass density imply that the compaction is bounded by the upper and lower values ervin et al. Page 13 of 35 w max = ρ max /ρ b and w min = ρ min /ρ b . We identify S ≡ w − as the swell factor ofa soil that is transformed from bank state to its maximally loose packed state.The surface voxels define a surface heightfield, z = h ( x, y ), that is used for contactswith the particles, the equipment, or any other objects in the macroscale model.It has a discrete representation h ij = h ( x i , y j ). The height value in a column withindex ( i, j ) is the centre position, z i (cid:48) of the top-most non-empty voxel, i (cid:48) = ( i, j, k (cid:48) ),plus the local mass fill ratio relative to that voxel centre, i.e., h ij = z i (cid:48) + ( ϕ i (cid:48) − / l , (21)This make the surface heightfield a continuous function of the solid occupancy, seeFig. 3. Between the grid points the surface height field is interpolated linearly.The response by the terrain is different for contacts with sharp and blunt geome-tries. The former lead to shear failure with a localized failure surface while the lattercause soil compaction. If the contacting body has a sharp cutting edge , a co-movingactive zone is predicted. The motion of the terrain inside the active zone is resolvedby particles and fluidized mass. See Fig. 3 for illustration. The basic shape of theactive zone is that of a wedge, defined by the cutting edge, the soil failure surfacethat extends from the edge to the free surface of the terrain, and the separationplane of the cutting body. The slope of the failure surface, θ , depends on the soil’sinternal friction, φ , and the orientation of the separation plane relative to the tothe terrain surface, β . From simulations with the microscale model we identify thefollowing model θ ( φ, β ) = π − (cid:18) φ + β (cid:19) , (22)This extends the classical Rankine failure angle π − φ to sloped terrain. Furthermore,to handle nonuniform distributions of material, the active zone is discretised in aset of parallel wedges in the lateral direction.As the cutting body moves through the soil, the active zone overlap with newvoxels. Overlapping soil is converted from solid mass, to particles or fluidized mass.The amount of converted mass is computed using Eq. (21), such that the surfacevoxel’s fill ratio take the exact value where the surface heightfield coincides withthe failure surface, see Fig. 3. The liberated mass is converted into new particles orgrowth of existing particles in the vicinity of the voxel, provided the local void ratioallows it. To avoid ending up with too small particles (poor performance) or toolarge particles (discretization error) the particle size is restricted to a given range[ d min , d max ]. The size range is chosen such that the number of particles remain withinthe limits of real-time simulation. Any residual mass, that cannot be converted intoparticle mass, is converted into fluidized mass in the voxel. At the next time-step,the fluidized mass is again candidate for formation and growth of particles. If itreaches a surface voxel outside the active zone, it is converted back to resting solidmass, whereby the surface heightfield is raised accordingly. The inertia of the coarseparticles is assumed to dominate over the dilute fluidized mass. Therefore, we simplymodel the flow of fluidized mass by conserved advection ∂ t ρ f + ∇ · ( ρ f u ) = r f . (23) ervin et al. Page 14 of 35 where u ( x ) is the coarse-grained particle velocity field, using Eq. (16), and r ( x , t )is a source (or sink) term for mass converted from (or to) particle mass or solidmass. Any fluidized mass that is found in a voxel with no particles is projected inthe direction of gravity towards the surface where it is converted into solid mass.When a cutting body is raised from the terrain surface a gap may arise betweenthe edge and the surface. In reality, the gap is normally closed by fine-grain soilflowing from the active zone. Mesoscale particles that are coarser than the gapcannot represent this. Instead, this is modelled by particles (and any fluid mass) inthe vicinity of the cutting edge losing mass to the surface voxels at the gap. This isillustrated in Fig. 3. The amount of mass necessary to fill the gap is computed usingEq. (21). Particles losing mass shrink correspondingly until they reach the lower sizelimit, where they are converted to fluidized mass. The process of converting soil fromsolid mass to particle and fluidised mass, and back, enable cutting and grading withhigh precision, despite the particle and voxel discretization being much coarser . Itshould also be noted that the operations for mass conversion and transport are byconstruction guaranteed to preserve the total mass to machine precision. current heightprevious heightcurrent blade tip posnew heightsolid massgrown/shrunk particlesconverted massprevious failure planecurrent failure planecurrent heightfield t n t n+1 t n+1 post fluidized mass h i m i m i sf active zone a)b) Figure 3
As a blade and its active zone moves into the terrain (a), new solid voxels are resolvedinto particles or fluidized mass that form the aggregate body. The voxel height value,corresponding to the solid occupancy, is found by projection to the failure plane of the activezone. When the blade is raised (b), the gap is filled by solid mass converted from fluidized orparticle mass in the vicinity of the edge.
The time evolution of the coarse particles, N CP , is governed by the equation ofmotion m b ˙ v b = f ext b + G T b B λ B b + G T b C λ C b + (cid:88) b (cid:48) ∈N bCP G T bb (cid:48) λ bb (cid:48) (24)and mass balance at voxel i containing the particles N i CP (cid:88) b ∈N i CP ˙ m b = − ˙ m i s − r i f , (cid:107) ˙ m b (cid:107) ≤ m b d b d (cid:48) max , d ∈ [ d min , d max ] . (25)where G T b B λ B b is the contact force from the cutting body and G T b C λ C b from theterrain surface. The change in solid mass, ˙ m i s = w i ˙ ϕ i ρ b V , is linked to the change in ervin et al. Page 15 of 35 the surface heightfield imposed by the moving cutting body. The change in particlemass depends on the change in solid mass, but is constrained by a maximum particlegrowth rate, d (cid:48) max , and size range. If the change in solid and particle mass doesnot match, the residual mass is converted to or from fluidized mass at rate r i f ,provided m i f ≥
0. Since the contact model is scale-invariant by construction, the pre-calibrated parameter p p will give the mesoscale model the desired bulk mechanicalproperties p b . Particles that come to rest outside the active zone are convertedto solid mass, whereby the surface heightfield increase correspondingly. Soil thathas undergone rapid flow tend to be more loosely packed than the original bankstate. To model this, each soil is assigned a swell factor S = ρ b /ρ min . Voxels withmass converted from particles to solid mass are thus given the compaction value w i = 1 / (1 + s ) and the local solid mass density ρ s i = ρ b / (1 + S ). The effect on thesoil’s strength by the change in compaction is addressed below.We assume that contacting blunt bodies cause soil compaction if the subsoil stressreaches critical values. Since the displacements are small it is not necessary to resolvethese deformations using particles. Instead, we operate directly on the local com-paction of solid mass w i . The contact forces f k = G T BC λ BC between the blunt body B and the terrain surface C are obtained from the macroscale model. The forces act atsome contact points x k that enclose a contact patch of area A b BC . The subsoil stress σ i in a voxel i at depth z i , and the resulting plastic deformation, can be estimatedfrom analytical solutions and numerical extensions [41, 42]. As a first order approx-imation following [43], we use the model, σ i = σ b BC (cid:20) − (cid:16) z i / (cid:112) A b BC + z i (cid:17) (cid:21) , fora circular distributed normal load, σ b BC = (cid:80) k f k /A b BC , on a semi-infinite elasticsolid. The compaction can be estimated from the compression index in Eq. (5) andby noting the relation to the void ratio w ≈ e b − e . Noting that the void ratioand compaction are related by w ≈ e b − e we obtain the following model for thecompaction in a voxel i due to subsoil stress at a depth z i beneath the contact w i = 1 + C c ln( σ i /σ b ) , (26)with clamping to the set bounds w i ∈ [ w min , w max ]. Typical values for the compres-sion index range between 0 .
001 (dense sand) and 1 (soft clay) and we use σ b = 1kPa is the consolidation stress for the bank state. A change in compaction altersthe mass density in the voxel to ρ i = ρ b w i and mass is shifted downwards columnwise to preserve the total mass. This reduce the fill ratio of the surface voxels at thecontact and the surface heightfield is reduced correspondingly. The increased levelof compaction also makes the material stiffer and stronger. From Eqs. (17) and (18)we obtain the following models for the bulk elasticity, dilatancy angle and effectiveangle of internal friction E i bulk = E (1 ± k E (cid:107) w i − (cid:107) n E ) , (27) ψ i = c ψ ( w i − w c ) , (28) φ i = φ b + ψ i , (29)where the stiffening parameters k E , n E and have default values 1 and 0 .
5, respec-tively, but can be tuned to represent more complex soil. The parameter c ψ = ∂ψ/∂w ervin et al. Page 16 of 35 control the rate of hardening, with 0 . − . w c , determines whether the soil at bank state is expanding orcompressing under shear. It is related to the bank state dilatancy angle ψ b by w c = 1 − ψ b /c ψ The compacted state is permanent until the soil is disturbed, e.g.,by earthmoving equipment interacting with it.
Macroscale model
The macroscale model focus on the rigid multibody dynamics of the equipment, B . The resting terrain is perceived as a surface heightfield, C . In each active zone,evolved in time by the mesoscale model, the soil is aggregated into a single bodywith the inertia, centre of mass and velocity of the particles and fluidized mass.These aggregate bodies , labelled A , have the degrees of freedom of rigid bodies butfinite mechanical strength. This is accomplished by compliant frictional-cohesivecontact constraints at the interfaces of the aggregate to the terrain surface, whichcurrently coincide with the failure plane, and to the equipment. The equations ofmotion are M A ˙ v A = f ext A + G T AB λ AB + G T AC λ AC , (30) M B ˙ v B = f ext B + G T B λ B + G T BA λ AB + G BC T λ BC + G T BD λ BD . (31)where the constraint forces G T AB λ AB and G T AC λ AC act on the aggregate from theequipment and the terrain failure surface, respectively. The equipment is held to-gether and actuated by G T B λ B . The separation resistance and the inertia of thesoil in the active zone is mediated by G T BA λ AB . Through G T BC λ BC , the equipmentexperience direct contact with the terrain surface, e.g., tyres, tracks or the exteriorof the bucket. Soil cutting objects may also be subject to a penetration resistanceforce G T BD λ BD . See Fig. 2 for illustration of the interaction interfaces.Each mesoscale active zone, A , is aggregated to a macroscale body with the fol-lowing mass, inertia tensor, centre of mass, linear and angular velocity M A = (cid:88) a ∈N AP m a + (cid:88) i ∈N Af m i , (32) I A αβ = (cid:88) a ∈N AP m a (cid:0) (cid:107) x a − x A (cid:107) δ αβ + x aα x aβ (cid:1) + (cid:88) i ∈N Af m i (cid:0) (cid:107) x i − x A (cid:107) δ αβ + x i α x i β (cid:1) , (33) x A = M − A (cid:88) a ∈N AP m a x a + M − A (cid:88) i ∈N Af m i x i , (34) v A = M − A (cid:88) a ∈N AP m a v a + M − A (cid:88) i ∈N Af m i v i , (35) ω A = I − A (cid:88) a ∈N AP m a ( x a − x A ) × v a + I − A (cid:88) i ∈N Af m i ( x i − x A ) × v i , (36)where N AP and N Af denote the set of particles and voxels with fluidized mass inactive zone A . ervin et al. Page 17 of 35
The contacts aggregated from the AB and AC interfaces are turned into thefollowing velocity constraint0 ≤ γ n λ n + G n v ⊥ ( λ n + f c ) ≥ , f c ≡ c b A (37) γ t λ t + G t v = 0 , (cid:107) λ t (cid:107) ≤ µ (cid:107) λ n (cid:107) . (38)Eq. (37) prevent relative motion in the normal direction, i.e., compressive or tensiledeformations. Note that, unlike model Eq. (10), this model does not include acontact overlap. This make the aggregate viscoplastic in nature, with no sense ofany reference configuration. The viscous damping is controlled by the parameter γ n = ∆ t/ ( c as E b ), where c as is the aggregate’s damping coefficient (referred to as aggregate stiffness multiplier in the implementation). If the force is tensile andreach the cohesive force limit, given by c b A , separation may occur at the interface.Eq. (38) prevent sliding motion in the failure surface and in the interface to theequipment, unless the forces reach the Coulomb condition. The friction coefficientin the failure surface is set to the soil’s internal friction µ = µ b . At the equipmentinterface the friction coefficient is set by the tool’s surface friction µ = µ tool . Thecontact points between the aggregate body A , the equipment B and the failuresurface C is the reduced set of contact points from the mesoscale model that bestapproximate the two contact areas, A AB and A AC , with four contact points each.This is illustrated in 2D in Fig. 2.If the cutting edge has a thickness or is equipped with teeth, there may be asignificant penetration force , G T BD λ BD in addition to the separation force. This ismodeled with the following velocity constraint G BD v B = 0 , G T BD λ BD ≤ f pt ≡ n t [ p t + ( c t + p t µ t ) / tan θ t ] A t , (39)where G BD = t B is a unit vector for the direction of penetration, e.g., pointing in thedirection of the teeth. For penetration to occur, the tool must overcome the forcelimit f pt [44], which depend on the tool’s friction coefficient µ t and adhesion c t , andon the number of teeth n t , with cross-section area A t , tooth angle θ t . The toothpressure, p t , is modelled using the finite cavity expansion theory [31]. In the plasticlimit the tooth pressure become p t = p +[Υ + ( α − p ] / [2 + α ] under cylindricalexpansion, where Υ = 2 c t cos φ/ (1 − sin φ ), α = (1 + sin φ ) / (1 − sin φ ). The meanpressure is computed by p = ρgz [(1 + K ) + (1 − K ) cos(2 β )], taking the lateralearth pressure and tool inclination into account [44], where K is the coefficient oflateral earth pressure, ρ is the specific soil mass density, z is the penetration depthand β is the insertion angle. A simple model for the coefficient of lateral earthpressure is given by Jaky [45] as K = 1 − sin φ . Complex digging tools
A bucket used for excavation or wheel loading is more complex than a blade. Itcan be seen as a bottom plate, with a curved back wall and sidewalls that allowmaterial to accumulate in the bucket. This deadload form a secondary separationplate , at an angle much steeper than the bottom plate. This affect the stress distri-bution, the shape of the active zone, by Eq. (22), and ultimately also the digging ervin et al.
Page 18 of 35 resistance. Furthermore, the forces from the soil that surrounds the bucket act bothtangentially (friction) and normally. This cause additional digging resistance as wellas strong resistance for motion in the lateral directions. Buckets are also used forother purposes than digging. It is a common operation to do surface leveling andcompaction using the bucket exterior. To support these features we treat a bucketas a composite model of an elementary digging tool and a set of soil deformers asdescribed below.A digging tool has a cutting edge e c , a parallel top edge e t , and a penetrationdirection t c orthogonal to these edges and to the normal of the bottom plate n c .This is the primary separation plate of the digging tool. See Fig. 4 for an illustration.A convex digging tool has an inner shape that is the void enclosed by the cuttingedge, top edge and the concave tool surface connecting them. The face that connectthe cutting edge and the top edge has normal n tc . If the bucket is full, this face actas a secondary separation plate. If the digging tool is a simple blade there is no innershape and n c = n tc . For concave tools we track the amount of material accumulatedin the inner shape (deadload) and do linear interpolation of the secondary separationplate n fill between n c and n ct . This changes the angle β in Eq. (22) and thus theslope of the failure surface. β θβ fill deadload e c e t n c n fill n ct t c β θ e c e t n c t c n d e d Figure 4
The left images illustrate a digging tool with a cutting edge e c (red line), separatingplate n c (red face), deformer edge e d (blue line), deformer face n d (blue face) and top edge e t (yellow line). The inner shape is indicated with the dashed lines. The normal n fill of the secondaryseparation plate range from n c to n ct depending on the fill ratio of the inner shape. The threeactive zones, discretized by vertical wedges, is illustrated from the rear (upper right) and fromabove (lower right). A soil deformer is simply a separation plate with no penetration resistance andno inner shape. It is defined by a deformation edge e d and a parallel top edge e t .They form a face with normal n d . The active zone from a soil deformer is notautomatically resolved using particles, unless the velocity of the deformation edgeexceed a given threshold. This avoid using particles to simulate the terrain when itremains at steady state.A bucket can thus be represented with soil deformers for the exterior faces andwith a digging tool for the cutting blade and the bucket interior. Fig. 4 shows abucket with two side walls, each assigned a soil deformer. Together with the digging ervin et al. Page 19 of 35 tool, this gives rise to a total of three active zones, which are discretised in thelateral direction. The backside of the bucket can also assigned a deformer, in whichcase the number of active zones sum to four. If not, the soil cannot be displaced bythe backside of the bucket other than by pure compression.The aggregate normal forces is modeled by Eq. (37). It is a velocity constraint,which are prone to numerical drift. The drift has no significance in soil cutting butprohibit proper resistance of soil deformers resting or being pressed onto a soil bed.A drift will cause them to sink unnaturally into the soil. This is remedied by addinga stabilizing perturbation ξg n to the constraint, i.e., γ n λ n + G n v → γ n λ n + ξg n + G n v , which turn it into a ordinary contact constraint. Transitioning between thevelocity constraint (providing smooth soil cutting) and contact constraint (avoidingartificial sinkage) is achieved by making the stabilization coefficient ξ a function ofthe orientation of the deformer relative to the terrain surface. Algorithm
The multiscale terrain model can be summarized by the Algorithm 1.
Algorithm 1
Multiscale terrain domains1: initialize terrain2: set surface heightfield h ( x, y ) and bank state soil parameters p b = [ φ b , c b , ψ b , E b ]
3: set initial voxel states s i = [ m i , v i , w i , ϕ i ] ,4: apply cellular automata to ensure the surface heightfield is soil-consistent5: initialize earthmoving equipment [ x , v , g ] B
6: set contact parameters for bodies and terrain surface [ µ, c, e, E ] BC
7: set digging tool edges and direction vectors [ e c , e t , t c ] and [ n t , A t , θ t ]
8: set deformers → [ e d , e t ] while running simulation do
10: conversion of resting and active soil11: for each body B intersecting the terrain C
12: compute active zones, discretise in wedges with failure angle θ ( φ, β )
13: convert resting soil in active zones to particle and fluidized mass14: convert resting particles outside zones to loose solid mass with w i = w min
15: apply cellular automata to re-distribute resting soil violating the angle of repose16: update the surface heightfield h ( x, y ) from the voxel state s i
17: do collision detection18: add contact constraints in the macroscale and mesoscale simulations19: for each active zone do
20: aggregate body A from voxels and particles using Eq. (32)-(36)21: add aggregate velocity constraint (37)-(38) with interfaces AB and AC end for
23: penetration resistance24: estimate pressure on edge/tooth and compute maximum penetration resistance25: add the penetration constraint Eq. (39)26: time-step Eq. (7)-(9)27: solve the macroscale MCP using direct LDLT solver → [ v , λ ] B
28: solve the mesoscale MCP using iterative PGS solver → [ v , λ ] b
29: update positions x i +1 = x i + ∆ t v i +1
30: soil compaction31: estimate the sub-soil stress from the surface contact forces BC
32: update the soil compaction w i using Eq. (26)33: displace resting solid mass vertically in compacted voxel columns34: update the soil strength parameters p ( w ) end while Implementation
The multiscale model was implemented in the C++ physics engine AGX Dynamics[46]. It supports real-time simulation of rigid multibody and particle systems with ervin et al.
Page 20 of 35 contacts, articulation joints and non-smooth dynamics. Numerical time integrationis made with the SPOOK stepper. A block-sparse LDLT solver with pivoting [47] isused as direct solver for the macroscale model, and for the equipment subsystem inthe meso- and microscale models. The solutions are exact to machine precision. Thedynamics of the contacting particles is solved to lower precision using a projectedGauss-Seidel (PGS), which is accelerated using domain decomposition for parallelprocessing and warm-starting [14, 15]. The microscale reference simulations are runwith 1 ms time-step and 500 PGS iterations. The multiscale model is run with 16 . Soil library
It is important that the bulk and particle representations are dynamically consistent.That means that the particle parameters, for each soil, need to be calibrated to thevalues that produce the desired bulk properties, e.g., internal friction, cohesion,angle of repose etc. Otherwise they do not describe the same material and thereis risk that the conversions between particle and continuum inject energy to thesystem, which can lead to numerical instability. Therefore, 15 soils were calibratedin advance to the bulk mechanical properties. These are listed in Table 1. The bulkparameters, at bank state, were chosen from tabulated values for different materials,e.g., gravel and sand. To narrow down the particle parameter search space, we wereguided by friction measurements of sand grains. The particles are given the Young’selasticity of E = 10 P a , Poisson’s ratio 0 .
15, specific mass density 2200 kg/m andcoefficient of restitution 0. Since there may exist many types of sand and gravel(with different distributions of size, shape, microscopic friction, packing densityand moisture level), each with distinct bulk properties, an integer is added to thename to distinguish between them. We also created a set of purely frictional soils (fs)and cohesive-frictional soils (cfs) with no particular real soil in mind. For testingpurposes a frictionless and a cohesive-frictionless soil (cs) were also created, butthey are not expected to be of practical use.A (virtual) triaxial test was used for parameter calibration, following the proce-dure and setup in [25], but with a lesser packing density. The model of the triaxialcell consists of frictionless rigid bodies for walls, 250 mm wide, driven by motors thatare servo-controlled to maintain a given normal stress, σ i , i = 1 , , xx, yy, zz .The particle samples have uniform size distribution with diameter ranging between27 and 33 mm, and initial porosity 0 .
42. The strength of each soil is tested by driv-ing the horizontal walls at a vertical speed of 5 mm/s while controlling the verticalwalls to maintain a given consolidation stress σ = σ . Tests are performed withdifferent consolidation stresses in the range from 1 kPa and 75 kPa. The motion ofthe walls and the stresses on them are registered during the simulation. Examplemeasurements of the stress deviator, σ − σ , as a function of the lateral strain areshown in Fig. 5 for the materials gravel-1 and wet-sand-1 . The internal friction ervin et al. Page 21 of 35 and bulk cohesion were determined from the Mohr’s circles at failure stress, seeFig. 6. The procedure was repeated, for each soil and consolidation stress, two orthree times with different particle initial states. The results are found in Table 1.Due to fluctuations and uncertainty in the Mohr method, the cohesion-free soilsshow an apparent cohesion of up to 2 kPa at most. In these cases, the bulk cohesionis explicitly set to zero. From the triaxial tests we also measure bulk elasticity asthe secant modulus, E b = ∆ σ / ∆ (cid:15) halfway before failure, and the dilatancy anglefrom ψ = arcsin [∆tr( (cid:15) ) / (cid:107) ¯ (cid:15) (cid:107) ]. The dilatancy angles are computed from the straincurves in Fig. 7. These bulk parameters are determined for a medium consolidationstress of 15 kPa. The procedure was repeated for all the 15 soils in Table 1. Table 1
A small library of soils with particle parameters pre-calibrated to desired bulk parametersusing a triaxial test. Bulk parameters Particle parametersSoil name φ b c b ψ b E b µ t µ r c p gravel-1 ◦ ◦ . . . sand-1 ◦ ◦ . . . sand-2 ◦ ◦ . . .
02 0 wet-sand-1 ◦ . ◦ . . . dirt-1 ◦ . ◦ . . . . dirt-2 ◦ . ◦ . . . . dirt-3 ◦
21 10 ◦ . . . . fs-strong ◦ ◦ . . . fs-weak ◦ ◦ . . .
05 0 cfs-strong ◦ . ◦ . . .
05 23 cfs-medium ◦ . ◦ . . .
05 12 cfs-weak ◦ . ◦ .
15 0 .
025 23 cfs-weakest ◦ . ◦ .
06 0 .
01 23 cs-weak ◦ . ◦ . . . frictionless ◦ ◦ .
25 0 . . Units deg kPa deg MPa kPa .
00 0 .
04 0 . (cid:15) [%] σ − σ [ k P a ] σ = 1 kPa σ = 5 kPa σ = 15 kPa σ = 30 kPa σ = 50 kPa σ = 75 kPa .
00 0 .
04 0 . (cid:15) [%] σ − σ [ k P a ] σ = 1 kPa σ = 5 kPa σ = 15 kPa σ = 30 kPa σ = 50 kPa σ = 75 kPa Figure 5
The principal stress as function of the vertical strain from triaxial test on gravel-1 (left) and wet-sand-1 (right). Initially the principal stress grows nearly linearly, from which thesecant modulus is determined, until shear failure occurs and the principal stress saturate.
Simulation tests
The multiscale model was tested and evaluated using two primary test systems,which are bulldozing and excavation in a flat soil bed. For reference, the same silua-tions were performed using the microscale model. The digging resistance, resultingterrain surface heightfield and computational performance of the two models wasthen compared. All materials in Table 1 was subjected to tests, see the supplemen-tary Video 1 and Video 2. The results for three frictional soils ( gravel-1 , sand-1 , ervin et al. Page 22 of 35 σ n (kPa) τ ( k P a ) c = 2.71 (kPa), φ =44.26 ◦ σ n (kPa) τ ( k P a ) c = 8.72 (kPa), φ =34.95 ◦ Figure 6
The Mohr circles for determining the internal friction and cohesion from triaxial test on gravel-1 (left) and wet-sand-1 (right). .
00 0 .
06 0 . | dev( (cid:15) ) | . . . tr( (cid:15) ) /3 triaxialTest gravel-1 σ = 1 kPa σ = 5 kPa σ = 15 kPa σ = 30 kPa σ = 50 kPa σ = 75 kPa .
00 0 .
06 0 . | dev( (cid:15) ) | . . . tr( (cid:15) ) /3 triaxialTest wet-sand σ = 1 kPa σ = 5 kPa σ = 15 kPa σ = 30 kPa σ = 50 kPa σ = 75 kPa Figure 7
The volumetric strain as function of the deviatoric strain from triaxial test on gravel-1 (left) and wet-sand-1 (right). The dilatancy angle is measured during the shear failure phase,where the slope is positive and close to linear. It was found ψ gravel-1 = 11 ◦ and ψ wet-sand = 8 ◦ for σ = 15 kPa. frictionless ) and three cohesive soils ( dirt-1 , wet-sand , cfs-weak ) are reportedmore detail. In the microscale simulations the soil bed consisted of 50e3 particles,with size uniformly distributed between 55 and 45 mm and specific mass densityof 2200 kg/m . The voxel size in the multiscale model was 0 . but the bulk mass density was set to 1474 kg/m whichcorrespond to a void ratio of 0 .
33. Default value for the aggregate stiffness multiplierwas set to 0 .
001 and 1 . sand-1 soil andpushes the material in front of it. After some distance, the material is dumped in apile, as the blade is stopped, lifted, and reversed. The cutting depth is 0 .
05 m andthe length is 5 m. The blade has mass 100 kg, sectioned vertically in two plates, 1 .
6m wide, 0 .
37 m tall, 0 .
02 m thick, and with a relative angle of 55 ◦ . It is attachedwith a lock constraint to a kinematic body that is driven with horizontal speed of0 . sand-1 and dirt-1 . Forthe penetration resistance, the blade’s edge is discretized by 80 teeth with 10 mmmaximum radius and length, and with 2 . gravel , dirt-1 and wet-sand which are given the value 20.Fig. 9 and 10 show the force from the terrain on the blade as a function ofthe horizontal position x . The force from the multiscale and reference modellargely follow each other, especially during the intermediate phase, starting afterthe blade is lowered into the ground ( x = − .
75 m) and ending where the blade ervin et al.
Page 23 of 35
Figure 8
Simulated bulldozing for comparing the multiscale terrain model (left) and themicroscale reference model (right, particles colour coded by particle height position) at time sand using material sand-1 . Second and third row show sideview of the microscale model (particlescolour coded by velocity, blue to red by 0 to 1.5 m/s) and the multiscale model, at time 0, 2, 4, 6,8, 10 and 12 s. See supplementary Video 3. gravel-1 sand-1 frictionless F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) gravel-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) sand-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) frictionless Figure 9
Force resistance from bulldozing using three frictional soils and the microscale model,with snapshot taken at 8 s. begins to rise ( x = 1 . ervin et al. Page 24 of 35 dirt-1 wet-sand-1 cfs-weak F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) dirt-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) wet-sand F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) cfs-weak Figure 10
Force resistance from bulldozing using three cohesive soils and the microscale model,with snapshot taken at 8 s. occur during the lowering and raising of the blade. The contributions to the forceon the blade from penetration and separation resistance in the multiscale model areshown in Fig. 11 for the case of sand-1 . It can be concluded that the penetrationforce is overestimated during the phase of raising the blade. It dominates the vertical − x [m] − − F x [ k N ] Terrain force F x ( x ) - bulldozing-terrain-sand-1totalpenetrationseparation − x [m] − F z [ k N ] Terrain force F z ( x ) - bulldozing-terrain-sand-1totalpenetrationseparation Figure 11
The digging resistance from bulldozing a bed of sand-1 soil simulated using themultiscale model. The force is divided in the contributions from the penetration and separationmodels. component of the force. The resting and densely packed soil is forced to expand atthe tip of the blade. When the blade is raised, the cutting edge moves through theshear zone, which should be looser and offer lesser resistance to penetration thandense soil. This is automatically captured by the resolved reference model but notin the multiscale model. In all cases, the pile of accumulated material in front ofthe blade is higher and steeper with the reference model than with the multiscalemodel. This is consistent with the larger time-step and lower PGS solver iterationcount in the multiscale model, which imply larger numerical errors that manifestthemselves as excessive slipping and rolling in the particle contacts. However, sincethe soil’s bulk properties have been calibrated in advance the errors affect the sizeof the aggregate but not on its fundamental shape or strength. In a sense, thesemesoscale errors are filtered out in the aggregation process to the macroscale modeland do not propagate into the digging resistance. As expected, the digging resistance ervin et al.
Page 25 of 35 − y [ m ] multi gravel-1 − y [ m ] ref − . . . x [m] − y [ m ] diff − y [ m ] multi sand-1 − y [ m ] ref − . . . x [m] − y [ m ] diff − y [ m ] multi frictionless − y [ m ] ref − . . . x [m] − y [ m ] diff -0.4-0.200.20.4 h ( x , y ) [ m ] Figure 12
The resulting surfaces from bulldozing in three frictional soils using the multiscalemodel (top row) and the microscale reference model (middle row). The difference is shown at thebottom row. − y [ m ] multi dirt-1 − y [ m ] ref − . . . x [m] − y [ m ] diff − y [ m ] multi wet-sand − y [ m ] ref − . . . x [m] − y [ m ] diff − y [ m ] multi cfs-weak − y [ m ] ref − . . . x [m] − y [ m ] diff -0.6-0.4-0.200.20.40.6 h ( x , y ) [ m ] Figure 13
The resulting surfaces from bulldozing in three cohesive soils using the multiscalemodel (top row) and the microscale reference model (middle row). The difference is shown at thebottom row. is much smaller for frictionless soil. It is smaller for the multiscale model thanfor the reference model. This can be understood by the fact that particles in themultiscale model slide over a frictionless plane while the motion of the particles inthe refence model are damped by the dissipative contacts (zero restitution) withthe irregular surface formed by the resting particles.The resulting surface height profiles, after a completed bulldozing cycle, are shownin Fig. 12 and 13. The multiscale and reference models are in good agreement for gravel-1 and sand-1 . The depth of the cut strips agrees within 10 millimetres andthe dimensions of the side berms are agree by roughly 10 %. The deviation in heightprofiles for the frictionless soil have the reason that is explained above. For thecohesive soils the deviation between the models is larger. The reference model with dirt-1 give wider side berms and a wider pile. The relative error is up to 50 %.For the strongly cohesive soils ( wet-sand-1 and cfs-medium ) the difference evenlarger, up to 100 %. For reference model almost all soil is accumulated in front ofthe blade. The particle cohesion appear to be much stronger in the reference modelthan with the multiscale model.The excavation test with dirt-1 can be seen in Fig. 14 and in the supplementaryVideo 4. The excavator model consists of four dynamic bodies and 4 joints: abucket on the end of an articulated arm, with an outer and inner boom, attached to ervin et al.
Page 26 of 35
Figure 14
Simulated excavation for comparing the multiscale terrain model (left) and themicroscale reference model (right, particles colour coded by particle height position) at time sand using material dirt-1 . Second and third row show sideview of the high-resolution model(particles colour coded by velocity, blue to red by 0 to 1.5 m/s) and the real-time model, at time0, 2, 4, 6, 8, 10 and 12 s. See supplementary Video 4. a revolving base on a static foundation. The test follows the setup of [24] where onlythe boom articulation joint is actuated, i.e., the inner arm is held fix. The bucketis 1 . . . . .
24 m long and weighs 300 kg. A hinge motor drives thearm to rotate at target speed 0 . . gravel-1 sand-1 frictionless F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) gravel-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) sand-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) frictionless Figure 15
Force resistance from excavation using three frictional soils and the microscale model,with snapshot taken at 5 s. soil and fil the bucket at x = 2 .
0. The horizontal force peak around x = 0 m. ervin et al. Page 27 of 35 dirt-1 wet-sand-1 cfs-weak F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) dirt-1 F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) wet-sand F x [ k N ] high-res F x ( x ) realtime F x ( x ) x [m] F z [ k N ] high-res F z ( x ) realtime F z ( x ) cfs-weak Figure 16
Force resistance from excavation using three cohesive soils and the microscale model,with snapshot taken at 5 s.
Up to this point, the force in the vertical direction has been dominated by thepenetration resistance, see Fig. 14. The weight of material inside and in front ofthe bucket increase gradually, as do the resistance to shearing the soil along thefailure plane. From x = 0 . x = 0 m these are in balance and after thisthe soil weight and the shear resistance dominate the vertical component of thedigging force. At x = − . F z at x = − .
5. There is good agreement between the models for the peak horizontalforce and the residual vertical force. Overall, the time-averaged forces agree within25 %. There is a trend for the multiscale model to overestimate the vertical diggingresistance during breakout around x = − . wet-sand-1 . With the reference model, the force fluctuates heavilybetween x = 1 . x = − . wet-sand . The force fluctuations grow large when thebeam hit the interior back wall of the bucket. The beam starts to compress, buckleand fail in an irregular manner. The aggregate in the multiscale model does notrepresent such modes of deformations and failure and does not produce large forcefluctuations. The height profiles after an excavation cycle are shown in Fig. 18 and19. For gravel-1 and sand-1 the width and height of the trench, side berms and thepile are in relatively good agreement. The depth of the trench agree within 10 mmand the other features within 10 % on average. For dirt-1 the side berms and pileis roughly 50 % larger for the reference model, and for wet-sand and cfs-weak thepile is more than twice as tall. The behaviour of the different soil materials underexcavation is shown in supplementary Video 2. It is clear that the strongly cohesivesoils are much more cohesive in the reference model compared to the multiscalemodel. ervin et al. Page 28 of 35 − x [m]0816 F x [ k N ] Terrain force F x ( x ) - excavation-terrain-dirt-1totalpenetrationseparation − x [m] − F z [ k N ] Terrain force F z ( x ) - excavation-terrain-dirt-1totalpenetrationseparation Figure 17
The digging resistance from excavation in a bed of dirt-1 soil simulated using themultiscale model. The force is divided in the contributions from the penetration and separationmodels. − y [ m ] multi gravel-1 − y [ m ] ref − − x [m] − y [ m ] diff − y [ m ] multi sand-1 − y [ m ] ref − − x [m] − y [ m ] diff − y [ m ] multi frictionless − y [ m ] ref − − x [m] − y [ m ] diff -0.4-0.200.20.4 h ( x , y ) [ m ] Figure 18
The resulting surfaces from excavation in three frictional soils using the multiscalemodel (top row) and the microscale reference model (middle row). The difference is shown at thebottom row.
The active zone model, Eq. (22), is validated using the reference model. Theinduced particle motion when digging in flat and sloping soil of the type dirt-1 with an excavator bucket can be seen in Fig. 20 and in the supplementary Video 5.The shape of the wedge formed by the mobilized soil is indicated. The angle betweenthe dynamic separation plate and the soil surface is estimated to β flat = 42 ◦ and β sloped = 45 ◦ in the two configurations. Eq. (22) imply failure planes with inclination θ modelflat = 47 ◦ and θ modelsloped = 45 ◦ , respectively. From the reference model we estimatethe failure angles at θ refflat = 53 ◦ and θ refsloped = 40 ◦ , respectively, which correspond toa 12 % error.We analysed also the computational speed of the multiscale model and the refer-ence model, aware that the results depend on implementation, optimization effortsand the hardware specification [3] The analysis was performed on the excavation testsystem shown in Fig. 14 and the result is summarized in Table 2. The number ofrigid bodies N rb and particles N p are given for the multiscale model and referencemodel. The number of equations, for the velocity update and constraint multiplier,are divided by what is treated by the direct solver, N rbeq , and by the iterative solver, N peq . The multiscale model has a real-time factor of 1 .
5, i.e., the computationaltime is 11 ms per simulation time-step 16 . [3] The tests were run on a desktop computer with Intel Core i7-7800X CPU at 3.5GHz and 32 GB RAM. ervin et al.
Page 29 of 35 − y [ m ] multi dirt-1 − y [ m ] ref − − x [m] − y [ m ] diff − y [ m ] multi wet-sand − y [ m ] ref − − x [m] − y [ m ] diff − y [ m ] multi cfs-weak − y [ m ] ref − − x [m] − y [ m ] diff -0.6-0.4-0.200.20.40.6 h ( x , y ) [ m ] Figure 19
The resulting surfaces from excavation in three cohesive soils using the multiscalemodel (top row) and the microscale reference model (middle row). The difference is shown at thebottom row.
Figure 20
Analysis of the active zone when digging in a flat (left) and sloped (right) soil bed witha complex bucket. The particles are colour coded by their velocity, with blue to red ranging from to m/s. See supplementary Video 5. for simulating a more complex vehicle at real-time speed. The reference model hasa real-time factor of 0 . N rbeq )processed by the direct solver, namely 36 additional equations (77% increase). Therelative speedup of 1500 is an effect of the 16.7 times larger time-step, 67 timesfewer equations for the iterative solver, and 50 times less iterations. Table 2
Performance analysis of the excavation test involving N rb rigid bodies and N p particles. .The direct MCP solver process N rbeq equations, and the iterative solver N peq equations. The multiscalemodel is times more efficient due to the lesser number of equations for the iterative solver, lessnumber of iterations, and larger time-step. N rb N p N rbeq N peq N pit time-step [ms] real-timemultiscale 8 950 85 1.8e4 10 16.7 1.5reference 4 50e3 48 1.2e6 500 1 0.001ratio 2 0.019 1.77 0.015 0.02 16.7 1500 Finally, the multiscale model is demonstrated in use with full vehicle models incomplex earthmoving scenarios. Fig. 21 shows a wheel loader digging in a steepwall of soil. See supplementary Video 6. The active zone in the cutting direction,discretised in five parallel wedges, is visualized in the right image. There are also ervin et al.
Page 30 of 35
Figure 21
A wheel loader digging in a step pile of soil. The resistance is too large for breaking outand instead the rear of the vehicle is raised from the ground. The right image shows the activezones, discretised in soil wedges, and the cohesive-frictional contacts between the aggregates andthe terrain. See the supplementary Video 6. active zones on each side of the bucket originating from the soil deformers. Thedifferent size and inclination of the wedges reflect the nonuniform distribution ofmass. The active zone in the digging direction is resolved with particles, but not theside deformers because of low lateral velocity. The contact points of the aggregatebodies are visualized with orange vectors. In the left image the lift cylinders areactuated to raise the bucket. However, the lift force cannot overcome the diggingresistance and, consequently, the rear of the wheel loader is lifted from the ground.This effect would not occur by the weight of the bucket and aggregate alone. It isnecessary to account also for the frictional-cohesive forces between the aggregates,the bucket, and the terrain, to capture the full resistance to breaking out from thewall.Variable soil compaction is illustrated in Fig. 22 and visualized by the intensity ofthe grey terrain. Medium grey represents nominal compaction at the bank state, atwhich the bulk strength parameters are defined. Light grey represent soil that hasdilated due to shear deformations, e.g., have been dug or pushed with the bucket.That soil has lower compaction and weaker strength according to the swell factorand dilatancy angle, respectively. Dark grey represent soil that has been compacted,e.g., due to compressive stress from the tires. It has higher compaction and strengththan the bank state value. The right images show the wheel loader driving into theleft pile of loose soil that is easily compacted.The involvement of additional rigid bodies interacting both with the earthmovingmachine and the terrain is demonstrated in Fig. 23 with a fullsized tracked excavatordigging a deep trench. The rocks contact directly with the excavator bucket, viathe direct solver, but also with the particles, through the iterative solver. The rocksforce the soil to distribute around the rock inside the bucket or pile up around therocks on the ground. The terrain is initialized with a high state of compaction suchthat the trench can sustain steep side walls. See the supplementary Video 7.The capability of representing large terrains, everywhere deformable, is demon-strated with the bulldozing example in Fig. 24 and the supplementary Video 8.Observe that the terrain is cut precisely at the dozer’s cutting blade, with a pre-cision much finer than the coarse particles, and match the vertical motion of theblade. The effect of bulldozer chassis and blade oscillations can be observed as wavepattern in the cut terrain surface. ervin et al.
Page 31 of 35
Figure 22
Demonstration of variable soil compaction and swelling. The subsoil stress from thetires cause compaction (dark grey). Particles that are converted back to resting solid becomemore loosely packed (light grey) than nominal bank state.
Figure 23
Demonstration interaction with other rigid bodies (red rocks) while digging a trenchwith a tracked excavator. See the supplementary Video 7.
Discussion
The results show good performance of the multiscale model and mostly fair agree-ment with the the reference model, needing little calibration. There are, however,a number of limitations with the method as well. The most significant deviationoccur for strongly cohesive soils, where the particle cohesion appear much strongerin the reference model than in the multiscale model. Furthermore, the shape of theactive zone is a crude approximation. It is clear from both the underlying theoryand numerical studies that the shear failure surface is not a plane and the activezone is not well approximated by a wedge, or several parallel wedges. The devi-ation in the digging resistance at breakout indicate that the slope of the failuresurface is not correct, or that the assumption of shear band localization is false.The active zone model in the present paper requires some manual work with defin-ing edges, direction vectors, teeth etc. Ideally, the dependency on such geometricfeatures should emerge from a model that merely takes the 3D geometry of a diggingtool as input. The damping coefficient of the rigid aggregate has not been derivedfrom the underlying theory or scrutinized using the microscale reference model. An ervin et al.
Page 32 of 35
Figure 24
Demonstration of high-precision bulldozing in a large terrain using a tracked vehicle.See the supplementary Video 8. alternative to improving the active zone model by analytical means is to take adata-driven approach. Following Wallin [49], it is possible to train a model, fromresolved reference simulations, to rapidly predict the digging resistance and the flowfield in a granular bed from a digging tool of particular 3D shape without explicitlydefining any cutting edges, direction vectors or teeth. Potentially this can be usedfor predicting the shape and mechanical strength of the active zone with higherprecision and generality, and possibly also the soil displacements. The drawbackof the data-driven approach is the need for running many simulations in advance,covering a wide range of soils, terrain shapes, and tool trajectories to have a usefulmodel.The computational bottleneck in the tested implementation is the mesoscale par-ticle simulation. It is run synchronously on CPU with the macroscale simulationusing fix and identical time-step. This is not necessary but made so for simplicitybecause the simulation software, AGX Dynamics, is designed for strong couplingbetween the multibody dynamics and the (nonsmooth) DEM simulation. It mightbe worth investigating alternative, possibly asynchronous, methods for simulatingthe active soil, e.g., smooth DEM on GPU.The mesoscale model support plastic compaction of the solid terrain but not sheardeformations. That is a major limitation for simulating deep ruts and vehicles orother objects sinking or being buried in the terrain.Finally, the presented method relies on a having a soil libary where particle para-maters are pre-calibrated to match the bulk mechaical paramaters. The currentlibrary, involving 15 different soils, can easily be extended to a wider range of moresoils, e.g., to include the over 100 virtual soil samples that was mapped in [25]. Fromsuch data-sets it is possible to identify mapping functions p p → p b , such that newsoils can be introduced on-the-fly. That is important for making machine learningmodels robust and transferable from one domain to another, e.g., from the simla-tion domain to the physical domain, using domain randomization [50]. However,the current model does not support inhomogeneous soil or mixing of two or moresoils. That extension is left for future development. ervin et al. Page 33 of 35
Conclusion
It has been found possible to simulate earthmoving operations in real-time with amodel that captures the rigid multibody dynamics of the equipment, the reactionforces from the terrain, and much of its deformations and flow dynamics. With amultiscale model the terrain’s active zones are represented simultaneously as a rigidbody, as particles and as a continuum. A direct solver is applied to the multibodysystem for high numerical precision and an iterative solver to the particle systemfor scalability. The models are made dynamically consistent through a soil librarythat is calibrated in advance using a reference model simulated at high-resolution.The performance, realism, and capability to represent a wide range of materials andscenarios make the solution suitable for simulation-based development control andautomation of earthmoving equipment.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author onreasonable request.
Competing interests
TB and MS are affiliated to Algoryx Simulation that develop the software used in the study.
Funding
Algoryx Simulation, eSSENCE (The Swedish e-Science Collaboration), and Ume˚a University.
Author’s contributions
The theory and algorithms was developed jointly by the authors. The implementation was made primarily by TB andSN. The simulation analysis was carried out primarily by MS. The authoring of the paper was lead by MS.
Acknowledgements
Anders Backman, Viktor Knutsson, Fredrik Nordfelth and Martin Nilsson at Algoryx Simulation has made importantcontributions in the software implementation, testing and optimization.
Author details Department of Physics, Ume˚a University, SE-90187 Ume˚a, Sweden. Algoryx Simulation AB, Kuratorv¨agen 2B,SE-90736 Ume˚a, Sweden.
References
1. Bender, J., Erleben, K., Trinkle, J.: Interactive simulation of rigid body dynamics in computer graphics (1),246–270 (2014)2. Haug, E.J.: Computer Aided Kinematics and Dynamics of Mechanical Systems, Vol 1; Basic Methods. Allynand Bacon, Boston (1989)3. Fang, Y.-g., Li, B.: Multiscale problems and analysis of soil mechanics. Mechanics of Materials , 55–67(2016). doi:10.1016/j.mechmat.2016.09.0034. Li, X., Moshell, J.M.: Modeling Soil: Realtime Dynamic Models for Soil Slippage and Manipulation. In: InComputer Graphics Proceedings, Annual Conference Series, pp. 361–368 (1993)5. Park, B.: Development of a virtual reality excavator simulator: a mathematical model of excavator digging anda calculation methodology. PhD thesis, Virginia Polytechnic Institute and State University, Blacksburg, Virginia(October 2002)6. Reece, A.R.: The fundamental equation of earth-moving mechanics. Proc. of the Institution of MechanicalEngineers (1964)7. McKyes, E.: Soil Cutting and Tillage. Developments in agricultural engineering. Elsevier, Amsterdam (1985)8. Holz, D., Beer, T., Kuhlen, T.: Soil deformation models for real-time simulation: A hybrid approach. In:Prautzsch, H., Schmitt, A., Bender, J., Teschner, M. (eds.) Workshop in Virtual Reality Interactions andPhysical Simulation ”VRIPHYS” (2009) (2009). doi:10.2312/PE/vriphys/vriphys09/021-0309. Holz, D., Azimi, A., Teichmann, M., Mercier, S.: Real-time simulation of mining and earthmoving operations: Alevel set-based model for tool-induced terrain deformations. Proceedings of the 30th ISARC, Montr´eal, Canada(2013)10. Holz, D., Azimi, A., Teichmann, M.: Advances in physically-based modeling of deformable soil for real-timeoperator training simulators. In: 2015 International Conference on Virtual Reality and Visualization (ICVRV),pp. 166–172 (2015). doi:10.1109/ICVRV.2015.6411. Jaiswal, S., Korkealaakso, P., ˚Aman, R., Sopanen, J., Mikkola, A.: Deformable terrain model for the real-timemultibody simulation of a tractor with a hydraulically driven front-loader. IEEE Access , 172694–172708(2019). doi:10.1109/ACCESS.2019.295616412. Renouf, M., Acary, V., Dumont, G.: 3d frictional contact and impact multibody dynamics. a comparison ofalgorithms suitable for real-time applications. In: Goicolea, J.M., Cuadrado, J., Orden, J.C.G. (eds.) MutlibodyDynamics 2005, ECCOMAS Thematic Conference, Madrid, Espagne (2005) ervin et al. Page 34 of 35
13. Holz, D.: Parallel Particles (P2): A Parallel Position Based Approach for Fast and Stable Simulation of GranularMaterials. In: Bender, J., Duriez, C., Jaillet, F., Zachmann, G. (eds.) Workshop on Virtual Reality Interactionand Physical Simulation. The Eurographics Association (2014). doi:10.2312/vriphys.2014123214. Servin, M., Wang, D., Lacoursi`ere, C., Bodin, K.: Examining the smooth and nonsmooth discrete elementapproach to granular matter. Int. J. Numer. Meth. Engng. , 878–902 (2014). doi:10.1002/nme.461215. Wang, D., et al. : Warm starting the projected gauss-seidel algorithm for granular matter simulation.Computational Particle Mechanics. (1), 43–52 (2016)16. Andreotti, B., Forterre, Y., Pouliquen, O.: Granular Media: Between Fluid and Solid, 1st Edition. CambridgeUniversity Press, Cambridge (2013)17. Antypov, D., Elliott, J.A.: On an analytical solution for the damped hertzian spring. EPL (Europhysics Letters) (5) (2011)18. Filla, R.: Quantifying operability of working machines. PhD thesis, Link¨oping UniversityLink¨oping University,Fluid and Mechatronic Systems, The Institute of Technology (2011)19. Du, Y., Dorneich, M.C., Steward, B.: Virtual operator modeling method for excavator trenching. Automation inConstruction , 14–25 (2016). doi:10.1016/j.autcon.2016.06.01320. Fernando, H., Marshall, J.A., Larsson, J.: Iterative learning-based admittance control for autonomousexcavation. Journal of Intelligent & Robotic Systems (3), 493–500 (2019). doi:10.1007/s10846-019-00994-321. Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. Geotechnique , 47–65(1979)22. Coetzee, C.J.: Review: Calibration of the discrete element method. Powder Technology , 104–142 (2017)23. Shmulevich, I., Asaf, Z., Rubinstein, D.: Interaction between soil and a wide cutting blade using the discreteelement method. Soil and Tillage Research (1), 37–50 (2007)24. Obermayr, M., et al. : A discrete element model and its experimental validation for the prediction of draft forcesin cohesive soil. Journal of Terramechanics , 93–104 (2014)25. Wiberg, V., Servin, M., Nordfjell, T.: Discrete element modelling of large soil deformations under heavyvehicles. Manuscript submitted for publication (2020)26. Wood, D.M., Graham, J.: Anisotropic elasticity and yielding of a natural plastic clay. International Journal ofPlasticity (4), 377–388 (1990). doi:10.1016/0749-6419(90)90009-427. Neto, E.d.S., Owen, D.R.J., Peri´c, D.: Computational Methods for Plasticity: Theory and Applications. Wiley,Chichester (2008)28. Forterre, Y., Pouliquen, O.: Flows of dense granular media. Annual Review of Fluid Mechanics (1), 1–24(2008). doi:10.1146/annurev.fluid.40.111406.10214229. Daviet, G., Bertails-Descoubes, F.: A semi-implicit material point method for the continuum simulation ofgranular materials. ACM Trans. Graph. (4) (2016). doi:10.1145/2897824.292587730. Hu, Y., Fang, Y., Ge, Z., Qu, Z., Zhu, Y., Pradhana, A., Jiang, C.: A moving least squares material pointmethod with displacement discontinuity and two-way rigid body coupling. ACM Trans. Graph. (4) (2018).doi:10.1145/3197517.320129331. Yu, H.S., Houlsby, G.T.: Finite cavity expansion in dilatant soils: loading analysis. G´eotechnique (2), 173–183(1991). doi:10.1680/geot.1991.41.2.173. https://doi.org/10.1680/geot.1991.41.2.17332. Lacoursi´ere, C.: Regularized, stabilized, variational methods for multibodies. The 48th Scandinavian conferenceon simulation and modeling (SIMS 2007), G¨oteborg, 30–31, 40–48 (2007)33. Acary, V., Brogliato, B.: Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanicsand Electronics, p. 526. Springer, Berlin (2008). http://hal.inria.fr/inria-0042353034. Radjai, F., Dubois, F.: Discrete-element Modeling of Granular Materials. Wiley-Iste (2011)35. Estrada, N., Az´ema, E., Radjai, F., Taboada, A.: Identification of rolling resistance as a shape parameter insheared granular media. Phys. Rev. E , 011306 (2011). doi:10.1103/PhysRevE.84.01130636. Wang, D., Servin, M., Berglund, T., Mickelsson, K.-O., R¨onnb¨ack, S.: Parametrization and validation of anonsmooth discrete element method for simulating flows of iron ore green pellets. Powder Technology ,475–487 (2015). doi:10.1016/j.powtec.2015.05.04037. Murty, K.G.: Linear Complementarity, Linear and Nonlinear Programming. Helderman-Verlag, Heidelberg(1988)38. Roux, S., Radja¨ı, F.: Texture-Dependent Rigid-Plastic Behavior. Physics of Dry Granular Media., pp. 229–236(1998)39. Bouchaud, J., Cates, M., Prakash, J.R., Edwards, S.: A model for the dynamics of sandpile surfaces. JournalDe Physique I , 1383–1410 (1994)40. Mitarai, N., Nori, F.: Wet granular materials. Advances in Physics , 1–45 (2006)41. Keller, T., D´efossez, P., Weisskopf, P., Arvidsson, J., Richard, G.: Soilflex: A model for prediction of soilstresses and soil compaction due to agricultural field traffic including a synthesis of analytical approaches. Soiland Tillage Research (2), 391–411 (2007). doi:10.1016/j.still.2006.05.01242. Zhu, Y., Chen, X., Owen, G.S.: Terramechanics based terrain deformation for real-time off-road vehiclesimulation. In: International Symposium on Visual Computing, pp. 431–440 (2011). Springer43. Madsen, J.: Mobility Prediction of Multi-Body Vehicle Dynamics Handling Simulations on Deformable Terrain.Doctoral thesis. University of Wisconsin-Madison. (2014).44. Bennett, N., Walawalkar, A., Heck, M., Schindler, C.: Integration of digging forces in a multi-body-systemmodel of an excavator. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-bodyDynamics (2), 159–177 (2016)45. Jaky, J.: Pressure in Silos. In: 2nd Int. Conf. Soil Mech, vol. 1. Balkema, Rotterdam, pp. 103–107 (1948)46. AGX Dynamics. . Accessed: 2020-068-12.