A physics-based deformable tire–soil interaction simulation capability that can be fully integrated into the monolithic multibody dynamics computer algorithm is developed by extending a deformable tire model based on the flexible multibody dynamics approach to off-road mobility simulations with a moving soil patch technique and it is validated against test data. A locking-free nine-node brick element is developed for modeling large plastic soil deformation using the multiplicative finite strain plasticity theory along with the capped Drucker–Prager failure criterion. To identify soil parameters including cohesion and friction angle, the triaxial compression test is carried out, and the soil model developed is validated against the test data. In addition to the component level validation for the tire and soil models, the tire–soil interaction simulation capability developed in this study is validated against the soil bin mobility test results. The tire forces and rolling resistance coefficients predicted by the simulation model agree well with the test results. It is shown that effect of the wheel loads and tire inflation pressures is well captured in the simulation model. Furthermore, it is demonstrated that the moving soil patch technique, with which soil behavior only in the vicinity of the rolling tire is solved to reduce the soil model dimensionality, leads to a significant reduction in computational time, thereby enabling use of the high-fidelity physics-based tire–soil interaction model in the large-scale off-road mobility simulation.

## Introduction

An accurate prediction of mobility capability on deformable terrains is crucial to off-road vehicle performance evaluation including the drawbar pull and rolling resistance; thus, high-fidelity physics-based simulation has been pursued to account for highly nonlinear structural and material behaviors of tires and soils to enable quantitative evaluations of off-road mobility capabilities. While the analytical and semi-empirical terramechanics models have gained acceptance in characterizing fundamental tire–soil interaction behaviors, complex soil deformation and highly nonlinear geomaterial behavior are, in general, difficult to capture [1–3]. The tire deformation effect is also either neglected or vastly simplified in many of the terramechanics models. Therefore, to enable quantitative evaluations of off-road mobility capabilities on deformable terrains under various speeds, loading, and soil conditions, physics-based computational models have been proposed using either the finite element (FE) [4–8] or discrete-element (DE) approaches [9–11]. In the finite element approach, soil behavior is modeled as a continuum. The large deformation soil behavior is described by phenomenological constitutive models such as Coulomb–Mohr and Drucker–Prager failure criteria. Continuum soil behavior including the strain hardening effect, which is important to characterize the soil sinkage and the normal contact pressure on deformable terrain can be effectively modeled using phenomenological constitutive models in the FE approach. On the contrary, the granular material behavior can be modeled by the interparticle dynamics using the DE method. Since soil behavior is modeled by a collection of moving particles, large deformation of noncohesive granular materials can be modeled using large-scale DE models by exploiting high performance computing techniques [10].

Although various FE deformable tire–soil interaction simulation models have been utilized with commercial FE software [6–8], use of the FE tire and soil models in multibody vehicle simulation requires the integration of the nonlinear finite element approach to the general multibody dynamics computer algorithm. Due to the essential difference in formulations and solution procedures adopted in the classical finite element and multibody-dynamics computer algorithms [12–14], finite element tire–soil models and multibody-dynamics vehicle models have to be either run separately using different computational models in different computational frameworks or run in parallel using co-simulation techniques. However, use of co-simulation techniques, in general, requires small step size to assure the numerical stability especially when discontinuous and sudden changes in the state of motion occur in the two different systems. In off-road mobility simulations, the deformable tire–soil interface is subjected to transient changes in forces due to severe vehicle traction and cornering maneuvers, thus special care needs to be taken in handling co-simulation algorithms. Furthermore, there are various commercial flexible tire models such as FTire and CDTire utilized in multibody vehicle dynamics simulations through co-simulation interface, and their tire–soil simulation capabilities are being developed [15,16]. Recently, a nonlinear FE tire model based on the absolute nodal coordinate formulation [17] was further integrated with DE deformable terrain model in Chrono software using a three-way co-simulation technique to enable high-fidelity off-road vehicle simulations [18].

While FE approaches have been pursued in the area of computational geomechanics, use of FE soil models in multibody dynamics vehicle simulations poses challenges concerning the finite element and multibody-dynamics coupling problem. Thus, a monolithic flexible multibody dynamics approach is proposed in this study. To this end, the physics-based FE tire simulation capability based on the absolute nodal coordinate formulation, which accounts for the coupling of the structural tire dynamics and transient tire friction under hard braking and cornering maneuvers, is further extended to off-road mobility simulations using monolithic multibody dynamics solver. This allows for integration of a high-fidelity FE tire–soil interaction model into an off-road vehicle model in the multibody dynamics analysis. The main contribution of this paper lies in the following points:

- (1)
A physics-based deformable tire–soil interaction simulation capability that can be fully integrated into the monolithic multibody dynamics computer algorithm is developed by extending the deformable tire model based on the flexible multibody dynamics approach to off-road mobility simulations with a moving soil patch technique. With this approach, the finite element soil model only in the vicinity of the rolling tire enters into the equations of motion to keep the soil model dimensionality relatively small regardless of the size of terrain and traveled distance considered in the off-road mobility simulation. Furthermore, impacts of the moving soil patch technique on the computational efficiency and accuracy are quantitatively evaluated.

- (2)
For modeling the large deformable continuum soil behavior, a nine-node brick element is proposed by introducing an additional center node to the standard eight-node brick element, which defines the second derivative of the global position vector, to alleviate element lockings without special techniques such as the enhanced assumed strain (EAS) approach. This allows for a straightforward implementation of the multiplicative finite strain plasticity theory along with the capped Drucker–Prager failure criterion to model the large plastic soil deformation exhibited as sinkage on deformable terrains.

- (3)
The deformable soil model is validated against the test data by the triaxial compression test, while an off-road flexible tire model considering its multilayer fiber-reinforced rubber (FRR) structure is experimentally validated for the contact pressure, the load-deflection curve, and the contact patch lengths under different wheel loads and inflation pressures. In addition to the component level validation, the tire–soil interaction simulation capability developed in this study is validated against the soil bin mobility test results.

## Physics-Based Flexible Tire Model and Validation

### Modeling of Tires Using Composite Shell Element.

In the flexible tire model utilized in this study for off-road mobility simulations, the complex FRR structure of a tire [19] is modeled using the shear deformable composite shell element based on the absolute nodal coordinate formulation, which allows for a straightforward implementation of nonlinear finite element models in general multibody dynamics computer algorithms without resorting to co-simulation techniques [17]. As shown in Fig. 1(a), the global position vector of a material point in the shell element is defined by [17,20]

*k*th node of element

*i*, one has

*N*layers, one has [17]

*l*th layer of the composite shell element, with which an orthotropic material model for the fiber-reinforced rubber as well as the hyperelastic material model for rubber are incorporated in the composite shell element [17,21]. $\epsilon il$ in Eq. (4) is the compatible strain vector of the

*l*th layer obtained from the assumed displacement field. To alleviate element lockings for the shell element, the EAS [22,23] and assumed natural strain [24,25] approaches are used to define the modified strain field as follows:

where $\epsilon \u0303$ is the covariant strain vector, for which the assumed natural strain is applied to the transverse shear and normal strain components, while the EAS is utilized for the transverse normal strain and the in-plane strain components. The constant transformation matrix $T$ is as given in the literature [20,23]. For more details on the shell element formulation and implementation, one can refer to the literature [17,20].

Using the shear deformable composite shell element, a tire model for a commercial off-road tire of 235/75R15 is developed. The tire cross section geometry and the layer property are shown in Fig. 2(a). The tire cross section is divided into the tread, shoulder, sidewall, and bead sections. An orthotropic material in each layer is stacked layer by layer to calculate the generalized elastic forces for a multilayered composite material, thereby allowing for incorporating various layer properties of the fiber reinforced rubber tire structure across the cross section. More details on the modeling of fiber reinforce rubber structure of tires using the present composite shell element are given in the literature [17]. The tire geometric data are extracted from the tire cut section. Inflation pressure is modeled by the distributed loads applied to the inner surface of the tire model. To model the contact with deformable soil, the outer surface of the finite element tire model is discretized into numerous contact nodes defined on the shell element surface as shown in Figs. 1(b) and 1(c). These contact nodes are different from the finite element nodes, and they are used to describe the tire tread geometry and to define contact forces with deformable terrain. To detect contact points between the tire and terrain surfaces, interpenetration between the contact nodes on the tire tread and terrain surfaces is evaluated. For contact node *i* on the tire tread surface shown in Fig. 1(b), the nearest contact node *j* on the terrain is detected, and then the distance $\delta ij$ normal to the terrain surface at contact node *j* is evaluated by $\delta ij=(ri\u2212rj)\u22c5nj$, where $nj$ is a unit normal to the surface at contact node *j*. If sign of the distance function is negative, two contact nodes are in contact and the normal contact force $Fn$ is defined by a compliant force function using $Fn=(k\delta ij+c\delta \u02d9ij\Vert \delta ij\Vert )nj$, where $\delta \u02d9ij$ is the deflection rate along the normal, $k$ is a spring contact, and $c$ is a damping coefficient. The contact spring constant *k* is determined such that the distributed contact stiffness within the contact patch remains constant regardless of the number of nodal points in contact. A proportional damping to the spring constant is assumed at each contact point for determining the damping coefficient *c*. The tangential forces $Ft$ at the contact point are calculated using either discretized distributed-parameter LuGre tire friction model as proposed in the previous study [17,26] or simplified slip-dependent friction force model [27]. Use of LuGre tire friction model allows for considering history-dependent friction behavior [28,29] and is suited for predicting detailed contact stress distribution over the contact patch under transient vehicle maneuvers [17,26].

### Comparison With Test Results.

To validate the flexible off-road tire model developed in this study, the tire is vertically loaded on a flat rigid ground to measure the normal contact pressure, the tire deflection, and the contact patch lengths. The tire is discretized into 24 shell elements in width and 80 elements in circumference, considering the multilayer fiber-reinforced rubber properties described in Fig. 2(a). The deformed cross section shape under a 4 kN wheel load is compared with the undeformed shape in Fig. 2(b). The normal contact pressure distribution along the longitudinal and lateral axes of the contact patch coordinate system at 4 kN wheel load are, respectively, compared with the test results in Figs. 3(a) and 3(b) for the inflation pressure of 230 kPa. The contact pressures predicted by the tire model agree well with the test results in magnitude and shape. Furthermore, the contact patch lengths in the longitudinal and lateral axes of the tire under various vertical wheel loads in the range of 1.0–8.0 kN are, respectively, compared with the test results in Figs. 4(a) and 4(b) for two different inflation pressures, 180 kPa and 230 kPa. It is observed from these figures that the simulation results agree well with the test results. Figures 5(a) and 5(b), respectively, present the lateral and vertical tire deflections under various vertical wheel loads in the range of 1.0–8.0 kN for two different inflation pressures (180 kPa and 230 kPa). The effect of the inflation pressure on the tire’s lateral and vertical stiffness is well captured in the tire model, and the load-deflection curves are in good agreement with the test data.

## Continuum Finite Element Soil Model

### Brick Element for Modeling Soils.

In this section, a continuum-based finite element soil model is introduced for application to multibody off-road mobility simulations using the physics-based flexible tire model discussed in Sec. 2. For this purpose, a nine-node brick element is developed to alleviate element lockings exhibited in the standard eight-node brick element. In this element, an additional node defining the second derivative of the global position vector is introduced to the center of the brick element, as shown in Fig. 6(a), to account for the quadratic polynomials in the shape function. This leads to the following global position vector of brick element *i*:

*k*(

*k*= 1,...,9) are defined by

where $\xi i$, $\eta i$, $\zeta i$ are the element natural coordinates associated with $xi$, $yi$, and $zi$, respectively. While the use of locking remedies such as selective reduced integration and the EAS method can alleviate the element lockings of the eight-node brick element, use of the nine-node element leads to a straightforward implementation of finite plasticity models for large deformable soils. This is because use of EAS requires additional iterative solution procedures to determine the internal EAS parameters for each element at every time step, making the solution procedure for large deformation soil plasticity model computationally intensive. Furthermore, use of the selective reduced order integration can lead to numerical instability when the element experiences large distortion. It is also important to notice here that the present nine-node brick element can be equivalent to the EAS-9 eight-node brick element [23] in accuracy, in which the linear polynomials are introduced to the assumed strain field using the additional nine EAS parameters. On the contrary, this is achieved in the present nine-node brick element by introducing the additional nine nodal coordinates defining the second derivative of the global position vector at the center node associated with the quadratic polynomials at the displacement level. While extra degrees-of-freedom (DOF) are needed, there is no element-level iterative solution procedure required to determine the enhanced strain field; thus, the finite strain elastoplastic forces and consistent elastoplastic tangent moduli for a continuum soil model can be implemented in a straightforward manner. Furthermore, this feature facilitates seamless integration of the continuum soil model in general multibody dynamics computer algorithms.

The preceding equation indicates that the plastic state, which is defined as the intermediate stress free configuration, can be conceptually recovered by an elastic unloading from a current configuration as illustrated in Fig. 6(b). In the large plastic strain regime, elastic and plastic deformations are coupled, resulting in the classical additive decomposition of elastic and plastic strains no longer valid for finite plasticity problems as in off-road tire–soil interaction problems. Furthermore, it is known that this formulation overcomes drawbacks of the hypo-plasticity formulation including erroneous dissipated elastic responses exhibited in incremental constitutive equations [31–33].

*i*th eigenvalue of spatial elastic stretch tensor $Ve$ defined by the polar decomposition $Fe=VeR$ with the rotation tensor $R$. The vector $bi$ defines the principal stretch vector. For an isotropic material, the Kirchhoff stress tensor that is conjugate to the logarithmic strain tensor is defined using the elastic coefficient matrix

**E**as

where $M$ is the constant generalized mass matrix and $Qe$ is the generalized external force vector that accounts for the contact forces with the tire model.

### Capped Drucker–Prager Failure Criterion and Return Mapping.

where $\Phi a$ is the Drucker–Prager yield function that defines the yield surface, while $\Phi b$ defines the elliptical cap yield function based on a modified Cam–Clay yield criterion, as illustrated in Fig. 7 [36]. For the equation of $\Phi a$, $J2(s)=s:s/2$ with $s(\tau )=\tau \u2212p(\tau )I$ and $p(\tau )=tr(\tau )/3$; *d* is the cohesion; and the parameters $\eta $ and $\xi $ are defined using the Drucker–Prager friction angle $\beta $ as $\eta =tan\u2009\beta /3$ and $\xi =1/3$, respectively [36]. In the equation of $\Phi b$, *a* and $b$ are parameters for the cone and cap sizes; $pt=(\xi /\eta )d$ defines the tensile yield pressure; $q(s)=3J2(s)$; and $M=3\eta $. It is important to note here that the paramour *a* in Eq. (15) is given as a function of the volumetric plastic strain $\epsilon vp$ to describe the strain hardening characteristics of soil and it is defined as [36]

where *C _{c}* and

*C*are identified experimentally using the one-dimensional consolidation test data, while

_{s}*e*

_{0}and

*p*

_{c}_{0}are, respectively, the void ratio and the mean effective stress at the initial state.

*n*+ 1 as follows:

*n*. Using the trial Kirchhoff stress tensor evaluated with Eq. (18), the yield functions are calculated to determine whether or not plastic loading occurs. If not, the stress and the plastic right Cauchy–Green tensor are updated as $\tau n+1e=\tau n+1trial$ and $(Cn+1p)\u22121=(Cnp)\u22121$ for the next time step. If plastic loading occurs, plastic strains need to be determined by bringing the trial stress state back to either (1) the capped Drucker–Prager yield surface, (2) the cone apex point, (3) the cap surface, or (4) the intersection as illustrated in Fig. 7(b). The detailed return mapping algorithm implemented in this study is found in the literature [39]. Accordingly, the incremental plastic multipliers can be determined and then the deviatoric Kirchhoff stress tensor $sn+1$ and the hydrostatic pressure $pn+1$ are updated. The Kirchhoff stress tensor is updated by $\tau n+1=sn+1+pn+1I$ to define the generalized elastoplastic force vector using Eq. (13). To update the plastic strain information, the total elastic strain tensor is updated as

*G*and

*K*are, respectively, the shear modulus of rigidity and bulk modulus. The elastic left Cauchy–Green tensor can then be updated as

*i*th principal logarithmic strain. Using Eq. (20), the inverse of the plastic right Cauchy–Green tensor at $t=tn+1$ is updated for the next step as

In other words, the cumulative plastic strain information is saved through the inverse of the plastic right Cauchy–Green tensor $(Cn+1p)\u22121$. It is important to notice here that the consistent elastoplastic tangent moduli have to be evaluated based on the current yield surface or point that the stress lies on to ensure the convergence of finite element solutions in the time integration process [39].

## Numerical Examples of Soil Model and Validation

In order to demonstrate the capabilities of the nine-node large deformable continuum soil element using the capped Drucker–Prager failure criterion, the triaxial compression test is carried out for validation. Furthermore, a pressure sinkage test model is developed to discuss the accuracy and convergence of the finite element solutions.

### Triaxial Compression Test and Validation.

The triaxial compression test is used to determine the strength and the stress–strain behavior of soil under uniform confining pressures [38]. As shown in Fig. 8, a cylindrical soil specimen subjected to a uniform confining pressure is vertically compressed until failure. The soil specimen is fully saturated first by supplying water from the bottom tube as shown in Fig. 8. The water circulates through the soil and is drained from the top tube. After that, a uniform confining pressure is applied to consolidate the cylindrical soil specimen and excessive water is drained from both tubes. The deviator stress is applied vertically such that the stress element in the soil is subjected to the major principal stress in the vertical direction. Three confining pressures (25, 50, and 200 kPa) are considered in this test to identify parameters of the continuum soil model with the capped Drucker–Prager failure criterion, including the Young’s modulus, shear modulus of rigidity, friction angle, and cohesion. The specification of the soil specimen is summarized in Table 1.

Specimen diameter | 71.10 mm |

Specimen height | 142.27 mm |

Water content | 11.20% |

Initial area | 39.70 cm^{2} |

Initial volume | 564.86 cm^{3} |

Density | 2.149 g/cm^{3} |

Dry density | 1.933 g/cm^{3} |

Confining pressure | 25 kPa, 50 kPa, and 200 kPa |

Specimen diameter | 71.10 mm |

Specimen height | 142.27 mm |

Water content | 11.20% |

Initial area | 39.70 cm^{2} |

Initial volume | 564.86 cm^{3} |

Density | 2.149 g/cm^{3} |

Dry density | 1.933 g/cm^{3} |

Confining pressure | 25 kPa, 50 kPa, and 200 kPa |

#### Test Results and Parameter Identification.

The deviator stress obtained in the test is presented as a function of the axial strain for three different confining pressures in Fig. 9. It is observed from this figure that the soil strength increases as the confining pressure increases due to the increase in soil consolidation. Furthermore, since the soil is partially saturated in the low confining pressure case, difference in the strength for 25 and 50 kPa confining pressures is small. To identify the cohesion and friction angle at failure, the deviator stress, principal stress, and mean effective stress at failure are summarized in Table 2 for the three confining pressures. Using this test result, the deviator stress and the mean effective stress at failure are plotted in Fig. 10, and the linear regression is carried out to identify the cohesion from the intercept and the friction angle from the slope as shown in this figure. The cohesion *d* and friction angle *β* used in the capped Drucker–Prager yield criterion are, respectively, identified as $d0=211$ kPa and $\beta =51.8\u2009deg$. Notice here that the corresponding Coulomb–Mohr friction angle is $\varphi =24.4\u2009deg$ using the equation: $tan\u2009\beta =6\u2009sin\u2009\varphi /(3\u2212sin\u2009\varphi )$ [38]. Young’s modulus and shear modulus of rigidity are identified using the slope of the stress–strain curve. Parameters for the strain hardening behavior defined by Eq. (17) are identified by the one-dimensional consolidation test of the same soil specimen, for which the change in the void ratio, defined as a ratio of the volume of void space to the solid volume, is shown as a function of the mean effective stress in Fig. 11. In this figure, the loading result is shown by black rectangular symbols, while the unloading is shown by red circular symbols. Using this test result, the linear regression is performed to determine the slopes at the loading and unloading stages that correspond to *C _{c}* and

*C*in Eq. (17), respectively, as shown in Fig. 11.

_{s}*C*and

_{c}*C*are then identified as 0.2035 and 0.00268.

_{s}Confining pressure (kPa) | Deviator stress at failure (kPa) | Major principal stress at failure (kPa) | Mean effective stress at failure (kPa) | |
---|---|---|---|---|

Test 1 | 25 | 444.4 | 469.4 | 173.1 |

Test 2 | 50 | 450.0 | 500.0 | 200.0 |

Test 3 | 200 | 808.6 | 1008.6 | 469.5 |

Confining pressure (kPa) | Deviator stress at failure (kPa) | Major principal stress at failure (kPa) | Mean effective stress at failure (kPa) | |
---|---|---|---|---|

Test 1 | 25 | 444.4 | 469.4 | 173.1 |

Test 2 | 50 | 450.0 | 500.0 | 200.0 |

Test 3 | 200 | 808.6 | 1008.6 | 469.5 |

#### Triaxial Soil Test Validation.

Using the soil model parameters identified by the triaxial and the one-dimensional consolidation tests, a numerical triaxial compression test model is developed using the nine-node brick soil element with the capped Drucker–Prager failure criterion for validation. The model parameters identified by the test for this soil specimen are summarized in Table 3. The deviator stress and axial strain relationship obtained using the nine-node brick soil element is compared with the test data at 200 kPa confining pressure in Fig. 12. It is observed from this figure that the nonlinear plastic failure and hardening behavior of the soil under triaxial compression test condition is well captured using the soil model developed and implemented in this study for off-road mobility simulation, and the result agrees reasonably well with the test result.

Density, $\rho $ | 2149 kg/m^{3} |

Young’s modulus, $E$ | 54.1 MPa |

Poisson’s ratio, $\nu $ | 0.293 |

Cohesion, $d$ | 210.9 kPa |

Drucker–Prager friction angle, $\beta $ | 51.78 deg |

Mohr–Coulomb friction angle, ϕ | 24.37 deg |

Compression index, $Cc$ | 0.20351 |

Swelling index, $Cs$ | 0.00268 |

Initial void ratio, $e0$ | 0.424 |

Initial mean effective stress, $p0$ | 0.358 MPa |

Density, $\rho $ | 2149 kg/m^{3} |

Young’s modulus, $E$ | 54.1 MPa |

Poisson’s ratio, $\nu $ | 0.293 |

Cohesion, $d$ | 210.9 kPa |

Drucker–Prager friction angle, $\beta $ | 51.78 deg |

Mohr–Coulomb friction angle, ϕ | 24.37 deg |

Compression index, $Cc$ | 0.20351 |

Swelling index, $Cs$ | 0.00268 |

Initial void ratio, $e0$ | 0.424 |

Initial mean effective stress, $p0$ | 0.358 MPa |

### Pressure Sinkage Test Model and Comparison.

To further discuss the accuracy of the continuum soil model using the nine-node brick element with the capped Drucker–Prager failure criterion and the multiplicative finite strain plasticity theory, a numerical pressure sinkage test of the same soil is carried out and the result is compared with an ABAQUS model. A square plate ($0.25\xd70.25\xd70.01$ m) on the top surface of the rectangular soil specimen ($0.8\xd70.8\xd70.6$ m) is pressed vertically as shown in Fig. 13(a) to obtain the pressure–sinkage curve. The contour plot in this figure describes the distribution of the norm of soil deflection in the specimen for visualization. The top surface of soil is subjected to contact force exerted by a rigid plate. A quarter model is used by imposing the symmetric boundary condition as shown in the same figure. In the ABAQUS model for comparison, the trilinear brick element using the reduced integration with hourglass control (C3D8R) is used together with the capped Drucker–Prager yield criteria. The material properties of the soil are same as those identified by the triaxial compression test (Table 3). The pressure–sinkage curves are shown in Fig. 13(b), exhibiting good agreement in the yielding behavior with the ABAQUS reference solution as shown in Fig. 13(b).

## Tire–Soil Interaction Simulation

### Tire–Soil Interaction Simulation With Moving Soil Patch Technique.

where subscript *r*, *t*, and *s* refer, respectively, to rigid bodies, tire elements, and soil elements. $C(qr,et,t)=0$ defines kinematic constraints to impose mechanical joint constraints and prescribed motion trajectories for vehicle and tire components. $\alpha $ is the vector of internal EAS parameters for the composite shell element [20]. It is important to notice here that kinematic constraints of the finite element tire model are imposed through intermediate coordinates that map the rigid body reference coordinates to the finite element nodal coordinates at the constraint definition point such that existing joint constraint library developed for the rigid body reference coordinates can be utilized for flexible bodies without modification. For more details on imposing joint constraints on finite elements based on the absolute nodal coordinate formulation, one can refer to the literature [40]. The preceding equations of motions are integrated forward in time simultaneously using the implicit HHT time integrator [41] with variable step size control.

It is important to note here that the high-fidelity tire–soil interaction simulation, in general, leads to high computational cost due to the large number of soil elements required to model a wide range of terrains. To address this issue, the moving soil patch technique [11] is introduced to the continuum finite element tire–soil interaction simulation algorithm for multibody off-road mobility simulation. Due to energy dissipation resulting from plastic deformation of soil, soil far from the contact region has no contribution to the tire–soil interaction behavior. For this reason, finite elements of soil far behind the tire are removed from the equations of motion after the tire passes, and then new soil elements are added to the front of the tire such that the number of degrees-of-freedom associated with the soil model remains small regardless of the traveled distance of a vehicle in the off-road mobility simulation. Having removed the soil elements behind the tire from the equations of motion, the nodal coordinates and plastic strains of the removed soil elements are saved such that the state of soil right before the element removal can be recovered any time when the other tire comes close and rolls over this soil portion. In other words, the soil element freezes temporarily and then defreezes to allow for multiple pass scenarios on deformable terrains without loss of generality. Furthermore, in order for the tire to be positioned around the center of the moving soil patch, distances between the center of the tire rim and edges of the moving soil patch are monitored such that the soil patch can be updated when the minimum distance becomes less than a threshold to eliminate the boundary effect.

### Soil Bin Off-Road Mobility Test.

In order to validate the physics-based tire–soil interaction simulation capability for multibody off-road mobility simulation, indoor soil bin mobility test shown in Fig. 14 was carried out. Soil used in the test is same as the one presented in the triaxial compression test discussed in Sec. 4, while the tire used in the test is same as the one discussed in Sec. 2 (235/75R15 commercial off-road tire). The speed of the moving carriage, to which the single tire is attached, is controlled to 1 m/s and the tire is free to rotate about its spinning axis without external tractive and braking torques. The test was operated at two different wheel loads (6 and 8 kN) produced by an actuator and three different inflation pressures (180, 230, and 280 kPa). The steering angle to the direction of travel (i.e., slip angle) is set to 3 deg, which is constant in time. The three-axis tire forces (Fx, Fy, and Fz) are measured by force transducers embedded in the tire rim and compared with the simulation results for validation. The soil sinkage was measured at the middle of the groove by a depth sensor by hand as shown in Fig. 14. The test is run twice for each test scenario: one on the right track in the soil bin and another on the left track as shown in Fig. 14. To produce noticeable soil sinkages, the soil is not compacted by a roller after loosening the soil by a tiller and flattening the soil surface for each test. The soil density, void ratio, and water content are calculated using soil samples randomly collected from the soil bin. Their mean values are used as parameters for the soil model. The soil density is 1556 kg/m^{3}; the water content is 8.21%; and the void ratio defined as a ratio of the volume of void space to the volume of solids is 0.893, which is higher than that of the triaxial soil test due to the lack of soil compaction by a roller. The variability of the measured void ratio is shown in Fig. 15. The mean value is used as the initial void ratio to identify the strain hardening curve. The cohesion is then identified as 12.80 kPa. The other elastic parameters are assumed to be same as those in the triaxial compression test.

### Numerical Results.

Using the numerical procedure developed in this study, the off-road mobility soil bin test model is developed with tire and soil model parameters identified individually as presented in previous and this sections. The tire model considering a simplified tread pattern is inflated by applying distributed pressure load to the inner surface first and then dropped on the deformable terrain by applying the vertical load to the center of the rim. The longitudinal velocity is prescribed as a driving constraint imposed on the rigid rim to move the tire forward at a constant speed. The soil bin width is selected to be 0.4 m such that the boundary effect becomes small enough. The deformed shapes of the tire and soil using the moving soil patch approach are presented in Fig. 16, where the contour plot describes distribution of the vertical soil deflection on the terrain surface. It is observed that the soil patch is updated as the tire travels. The symmetric boundary condition is imposed on the interface planes of the moving soil patch. The tire forces (Fx, Fy and Fz) obtained by the test and the simulation are compared for the wheel load of 6 and 8 kN in Figs. 17(a) and 17(b), respectively. The tire inflation pressure is 230 kPa for both cases, while the coefficient of friction is assumed to be 0.25. The tire force test data exhibit small fluctuation due to the tire test rig carriage movement and actuation, yet the time-domain nominal tire force test and simulation results are in good agreement under different vertical wheel load (Fz) conditions. While Fx and Fy are small due to free rolling condition, the figure clearly demonstrates that effect of the vertical force Fz on Fx and Fy is well captured in the simulation model and they are in good agreement with the test data.

In order to discuss the effect of the length of the moving soil patch on accuracy, two soil patch models with different lengths (0.8 m and 1.25 m) are considered and compared with the 2.75-m-long complete soil model. Each finite element size in all soil models remains the same to ensure a fair comparison. The side view of the tire and soil deformed shapes is shown in Fig. 18 for these two cases, where the contour plot describes distribution of the von Mises stress in the central plane of the soil bin model for visualization. As shown in this figure, a shorter soil patch model leads to more frequent update of the soil patch. The soil patch is updated when the distance between the rim center and the forward edge of the soil patch becomes less than 0.4 m to eliminate the boundary effect in this problem. This results in the 0.8-m-long moving soil patch being updated at every traveled distance of 0.05 m, while the 1.25-m-long moving soil patch being updated at every traveled distance of 0.5 m as illustrated in Fig. 18. The tire rim vertical position for the traveled distance of 2.0 m is presented in Fig. 19, and they are compared with the complete soil model as the reference solution. In this simulation, the tire is dropped on the soil from the initial state and then the steady-state rolling condition is achieved, resulting in an approximately 0.06 m sinkage of the soil. It is observed from this figure that the 0.8 m and 1.25-m-long moving soil patch models lead to the same result as the reference complete soil model.

To discuss the computational cost of the high-fidelity physics-based tire–soil interaction simulation model using the moving soil patch approach, CPU times of the complete and moving soil patch models are summarized in Fig. 20 for different traveled distances and terrain lengths. Intel Xeon, CPU E5-2687W v3, 3.10 GHz is used with OpenMP parallel computing capability for ten threads. The tire traveling speed is assumed to be 5 m/s. In the complete soil model, the total soil bin length is determined by the total traveled distance required in the simulation. Thus, length of the complete soil model is selected to be 0.75 m longer than the traveled distance of the tire in the simulation to eliminate the boundary effect. For instance, in the case of 1.0-m-long tire rolling simulation, the soil bin length is selected to be 1.75 m. The numbers of degrees-of-freedom of the soil elements for the soil bin length of 1.75 m, 2.75 m, and 3.75 m are, respectively, 91,203, 142,923 and 194,643, while the numbers of degrees-of-freedom of the 0.8-m and 1.25-m-long moving soil patch model are 42,069 and 65,343 regardless of the traveled distance. The number of degrees-of-freedom of the tire model remains 12,000. It is observed from Fig. 20 that the CPU time increases exponentially as the traveled distance increases in the complete soil bin model due to the significant increase in the soil model dimensionality, while CPU time increases at a rate proportional to the traveled distance in the moving soil patch models. This is attributed to the fact that the soil model dimensionality remains constant regardless of the traveled distance; thus, the computational time depends solely on the traveled distance (i.e., simulation time) in this case. Furthermore, use of the smaller-size moving soil patch can reduce the total CPU time effectively. For instance, CPU times of two moving soil patch models (0.8 and 1.25 m long) are, respectively, 13.9 h and 22.3 h for the traveled distance of 2.0 m as shown in Fig. 20, and this difference increases as the traveled distance (i.e., simulation time) increases. Accordingly, significant reduction in computational time can be achieved in long terrain simulation using the moving soil patch model, which is important to enable the use of the high-fidelity physics-based tire–soil interaction model in large-scale off-road mobility simulations. An optimum size of the moving soil patch and threshold for an update of the soil patch, on the other hand, depend on simulation scenarios (i.e., soil properties, wheel loads, and slips) and desired accuracy. Furthermore, to further improve the computational efficiency, distributed memory parallel computing using message passing interface (MPI) could be used as demonstrated for tire–soil interaction simulations using the present tire model along with DE soil model [18].

### Comparison With Test Results.

To ensure the accuracy of the tire–soil interaction simulation capability under different wheel loads and tire inflation pressures, the soil sinkages obtained by the test and simulation are compared in Fig. 21, where simulation results are taken from those at the end of simulation since the steady-state motion is achieved. It is observed from the simulation results that increase in the wheel load and the tire inflation pressure leads to an increase in the soil sinkage. Since sinkages were measured by hand, there exists relatively high variability in the test data. However, the effect of the wheel load on the sinkage is clearly observed. The increase in the sinkage by higher inflation pressure is caused by higher contact pressure on deformable soil. Furthermore, an increase in the wheel load and inflation pressure leads to an increase in the magnitude of the longitudinal tire force as observed in Fig. 22. Since the tire has larger sinkage as the inflation pressure increases, the tire experiences larger resistance forces from soil to maintain a constant speed and this results in the larger longitudinal tire forces as observed in both simulation and test results. This phenomenon is clearly observed as the rolling resistance coefficient (Fx/Fz) in Fig. 23, in which the rolling resistance coefficient increases as the inflation pressure increases, whereas it remains nearly constant regardless of the wheel loads. This trend is well captured in the simulation results. The similar trend is observed for the lateral forces in Fig. 24. The lateral tangential force coefficient (Fy/Fz) remains nearly constant for different wheel loads, while the higher inflation pressure slightly increases the lateral tangential force coefficient as shown in Fig. 25. These simulation results agree well with the test results in magnitude and trend, considering the variability of soil data. It is evident from this comparison and individual validation of the tire and soil models discussed in Secs. 2 and 4 that the fundamental tire, soil, and their interaction behavior are well predicted using the simulation model developed in this study. Further validation study would be pursued in the future for traction and braking scenarios as well as transient cornering maneuvers.

## Summary and Conclusions

In the analysis of an off-road vehicle on deformable terrain, the tire and soil models undergo large deformation with highly nonlinear material behavior. The tire structure consists of layers of plies and belts embedded in rubber and the modeling of multilayer composite tire structure is critically important for predicting the normal contact pressure and slip distribution in the contact patch. In other words, the tire structural dynamics, the normal contact pressure, and tangential forces (i.e., shear contact stresses) are fully coupled in the nonlinear tire dynamics simulation. In the off-road vehicle simulation, modeling of the dynamic coupling between the large deformable nonlinear tire and terrain makes the physics-based off-road simulation technically challenging. Soil exhibits highly nonlinear yielding behavior with large permanent deformation and it is modeled by either FE or DE approaches. While FE approaches have been pursued in the area of computational geomechanics, use of FE soil models in multibody dynamics vehicle simulations poses challenges concerning the finite element and multibody-dynamics coupling problem. Thus, a monolithic flexible multibody dynamics approach was proposed in this study. For this purpose, the physics-based finite element tire simulation capability based on the absolute nodal coordinate formulation was further extended to off-road mobility simulations. For modeling the large plastic soil deformation exhibited as sinkage on deformable terrains, a nine-node brick element that has an additional center node for the second derivative of the global position vector was proposed using the multiplicative finite strain plasticity theory along with the capped Drucker–Prager failure criterion. The deformable soil model was validated against the test data by the triaxial compression test, while an off-road flexible tire model considering its multilayer fiber-reinforced rubber structure was validated for the contact pressure, the load-deflection curve, and the contact patch lengths under different wheel loads and inflation pressures. With these models, the finite element tire–soil interaction simulation model was developed and implemented in the monolithic multibody dynamics algorithm with a moving soil patch technique, with which the soil elements only in the vicinity of the rolling tire enter into the equations of motion to keep the soil model dimensionality relatively small regardless of the terrain coverage and the traveled distance considered in the off-road mobility simulation. The tire–soil simulation capability developed was validated against the indoor soil bin mobility test. It was demonstrated that tire forces and rolling resistance coefficients predicted by the simulation model agree well with the test results and effect of the wheel loads and tire inflation pressures was well captured in the simulation model developed in this study.

## Acknowledgment

This research was supported by the Automotive Research Center (ARC) in accordance with Cooperative Agreement W56HZV-04-2-0001 U.S. Army Tank Automotive Research, Development and Engineering Center (TARDEC), Caterpillar Inc., and Yokohama Rubber Co., Ltd.

## Disclaimer

Reference herein to any specific commercial company, product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or the Department of the Army (DoA). The opinions of the authors expressed herein do not necessarily state or reflect those of the U.S. Government or the DoA, and shall not be used for advertising or product endorsement purposes.

## Funding Data

Tank Automotive Research, Development and Engineering Center (Grant No. W56HZV-04-2-0001).