## Abstract

In this study, model order reduction of high-fidelity off-road mobility models is explored to address the computational intensity of nonlinear finite element deformable tire–soil interaction models. To this end, a model order reduction procedure for the tire–soil interaction model is developed with the proper orthogonal decomposition (POD), and it is integrated into the off-road mobility simulation framework, leveraging high-performance computing. The POD is, however, limited in that the modes are dependent on snapshot data collected during the running of a full order model, limiting the modes to being accurate only for the specific scenario from which they were collected. Due to this limitation, a method of mode adaptation through interpolation on a tangent space of the Grassmann manifold is investigated to allow modes to be predicted for cases in which a full order model has not been run. It is demonstrated by several numerical examples that the POD modes are effective at retaining predictive accuracy while reducing computational time. The results show that adapted POD modes are more capable of characterizing the behavior of the model than modes produced at a different value of the simulation parameter. The POD-based reduced order modeling approach is further extended to the full vehicle simulation on deformable terrain through the co-simulation coupling algorithm by leveraging the high-performance computing technique.

## 1 Introduction

The ability to predict off-road mobility is critically important not only for vehicle performance and design evaluation but also for mission planning. In order to predict vehicle behavior on deformable terrain, a number of models have been developed with varying degrees of complexity. In the simple terramechanics models, such as those developed by Bekker and Wong [1,2], the basic soil material behavior is modeled by the empirical soil pressure-sinkage curves as well as the shear stress-displacement curves, which require extensive in situ measurements. While these models provide physical insight that facilitates a basic understanding of the mechanics of tire–soil interaction, soil deformation is vastly simplified; thus complex deformable soil behavior exhibited in various vehicle maneuvering scenarios cannot be predicted. The soil contact model goes a step further to describe the spatial distribution of soil deformation using a discretized terrain surface with numerous grid points, on which the pressure-sinkage curves are defined [3]. These models are semi-empirical in nature. Therefore, physics-based high-fidelity numerical models such as the finite element (FE) method [4,5] and the discrete element (DE) method [68] are widely used to predict complex soil behavior resulting from the dynamic interaction with tires. In the DEM model, the soil is represented by a large collection of small particles that interact to describe the granular soil behavior. In the FE model, on the other hand, the soil is treated as a continuum, and the plasticity and hardening of the soil are considered using yield theories such as the Mohr–Coulomb and Drucker– Prager failure criteria.

The improvements in model accuracy with the physics-based terrain models, however, have come at the expense of computational efficiency, which limits the effectiveness of a model as the model dimensionality increases significantly to ensure accuracy. For this reason, extensive efforts have been made to effectively reduce computational costs of the high-fidelity physics-based models in varied ways, including the use of high-performance computing techniques with graphics processing units for DE models [9]; parallelized fast solvers [10]; and multiscale modeling to eliminate limitations of the existing FE and DE soil models [11,12]. Since soil deformation is produced by complex grain-scale particle motions occurring in the order of millimeters to micrometers, a scale separation technique is utilized effectively in the multiscale approach. That is, the macroscale soil deformation and the grain-scale granular soil material behaviors are modeled by separate computational models that are connected through the scale-bridging algorithm. Furthermore, the hierarchical multiscale terrain model was further enhanced by introducing data-driven neural network surrogate models in place of the lower-scale DE representative volume elements (RVEs) for computational speedup [13]. However, the computational improvement of the upper-scale FE model remains an open issue in off-road mobility simulations.

In this regard, various model order reduction (MOR) techniques have been studied in applied mathematics and computational fluid dynamics (CFD) for the sake of increased computational efficiency with the least possible loss in model accuracy [14]. In CFD, MOR using the proper orthogonal decomposition (POD) of nonlinear systems has been applied and proven to be effective in reducing the model dimensionality while ensuring accuracy [15]. POD is a method in which “snapshots” of the state of the system are taken over time during the transient simulation of a full order model (FOM) and used to develop reduced order bases consisting of completely independent modes. It has become especially popular in CFD due to its ability to produce reduced order bases and to extract underlying flow structures [16]. POD has been used to identify structures in plane Couette flow [17], reduce the order of circular cylinder vortex shedding models [18], and model flow over an airfoil for inverse design [19] among other applications. Because a POD basis is produced through a transient full order solution, its modes are dependent on the specific parameters of the original simulation. In order to make the modes more robust, bases at different operating points are generated and the sets of modes are combined through the reorthogonalization process [20]. Others have tried to make POD modes more robust through the use of artificial neural networks [21], multilevel algorithms to incorporate snapshot data incrementally [22], and interpolation between orthogonal bases through the use of a tangent space to the Grassmann manifold [23].

Despite all of the advancements in model order reduction using POD with application to CFD and other mechanical simulations, no studies have been explored in reduced order modeling for off-road vehicle mobility due to the complexity of incorporating path-dependent material behavior of deformable soil in the reduced order model, along with the effect of contact with deformable tires. It poses critical challenges in producing a reliable reduced order model for tire–soil interaction problems. This study, therefore, aims to develop a numerical procedure to integrate the POD-based MOR into a high-fidelity off-road mobility simulation framework, leveraging high-performance computing. With this approach, the potential of the POD method in producing a reduced-order tire–soil interaction model is examined by rigorously evaluating its predictive ability and computational improvements. Furthermore, a method of mode adaptation is applied through interpolation of reduced order bases on a tangent space of the Grassmann manifold as a means of basis adaptation for practical use in off-road mobility prediction.

## 2 Adaptive Model Order Reduction

In this section, the model order reduction technique using POD is briefly summarized. Furthermore, the adaptive MOR technique with a tangent space of the Grassmann manifold is overviewed to account for changes in modeling parameters of the reduced order model (ROM) with POD modes.

### 2.1 Proper Orthogonal Decomposition.

POD is a statistical method that decomposes systems with a large number of interdependent variables into a smaller number of independent variables [24]. In this approach, a set of orthogonal bases, ${φ k}k=1nϕ$, is generated by minimizing the difference between the full and reduced order solutions over time as follows:
${φ k}k=1nϕ=arg min{φ k}k=1nϕ∫t0tN‖x(t)−Px(t)‖2dt$
(1)
where $x(t)$ is a vector of the physical coordinates in the full order solution, while the matrix P is given by $P=∑k=1nϕφ kφ kT$ with the orthonormal reduced order basis vector $φ k$. The physical coordinate vector can then be approximated as
$x(t)=∑k=1nϕφ kξk=Φ ξ(t)$
(2)
where $ξ$ is a vector of the reduced order coordinates defined by $ξ=[ξ1ξ2⋯ξnϕ]T$, while $Φ$ is the reduced order matrix given by $Φ=[φ 1φ 2⋯φ nϕ]∈ℝn×nϕ$ with n being the number of physical coordinates and $nϕ$ being the number of truncated reduced order modes. The basis vectors are obtained as solutions to the optimization problem posed by Eq. (1). They are given as eigenvectors of the covariance matrix $S=XXT$ of the snapshot matrix X given by
$X=[x(t0)−x¯x(t1)−x¯⋯x(tN)−x¯]∈ℝn×N$
(3)
where $N$ is the number of snapshots. The average of the solution vectors from time $t0$ to $tN$ is given by $x¯$, and it is subtracted from the solution vector $x$ to have zero average over time. In practice, the basis vectors can be obtained directly by performing the singular value decomposition of the snapshot matrix X of the full order solution as follows [24]:
$X=UΣVT∈ℝn×N$
(4)
where the left singular matrix $U$ describes the orthogonal spatial modes, the right singular matrix $V$ defines the orthogonal temporal modes, and the singular values in the diagonal of the matrix $Σ$ represent the amount of energy stored in each mode. Accordingly, the truncated reduced order basis matrix is obtained as
$Φ={Uk}k=1nϕ∈ℝn×nϕ$
(5)
The lower energy modes are truncated such that the remaining energy comprises 99% or greater portion of the variation in the original model [24]. In this way, a tradeoff must be made between the gains in computational efficiency and the losses in accuracy. In many large-scale models, on the other hand, the number of coordinates is much larger than that of snapshots (i.e., $n≫N$), making the singular value decomposition of Eq. (4) computationally very expensive. Since the nonzero eigenvalues of the covariance matrix $S=XXT∈ℝn×n$ are the same as those of the matrix $S′=XTX∈ℝN×N$, one can obtain the reduced order bases by solving the following smaller eigenvalue problem with the matrix $S′∈ℝN×N$:
$S′uk=λkuk$
(6)
where $uk∈ℝN$ is the eigenvector for the kth mode. Accordingly, the reduced order bases can be obtained more efficiently for a large-scale model as
$Φ={φ k}k=1nϕ=Xuk/λk∈ℝn×nϕ$
(7)

### 2.2 Adaptation of Reduced Order Proper Orthogonal Decomposition Modes.

Due to the dependence of POD modes on the snapshot matrix X as dynamic data from a full order simulation, POD is weak to changes in operating parameters. A change in a simulation parameter renders the modes inappropriate. In order to address the limitation of POD, a method of mode adaptation is required to generate modes for situations for which full order models have not been run and dynamic snapshot data $X$ is not available. For a given set of POD modes corresponding to different values of some modeling parameter, a method of interpolation is desired to produce modes for additional values of the modeling parameter. A naïve approach would be to simply interpolate each element of the basis matrices directly. The problem with this method is that POD modes are orthonormal. Therefore, directly interpolating between orthonormal matrices is highly unlikely to produce a new orthonormal matrix, making the interpolated modes useless for MOR purposes. In order to interpolate between bases and ensure orthonormality, a different interpolation scheme has to be pursued.

All orthogonal bases of the same dimension belong to the same Grassmann manifold $G(nϕ,n)$, which is defined as the set of all $nϕ$-dimensional subspaces of the vector space $ℝn$. This is because, for each of these bases, $Φ∈ℝn×nϕ$, the columns span a subspace which can be represented as a point on the surface of $G(nϕ,n)$. For each of these points on the manifold, a tangent space T with the same dimensions exists. Since T is a vector space originated at the point of tangency, it is a flat space where interpolation can be conducted freely [23]. The points on T are the same dimension as the points on $G(nϕ,n)$ and will be denoted with a tilde as $Φ̃$ to differentiate them from the points on the manifold $Φ$.

In order to map a point on a tangent space to the Grassmann manifold, an explicit exponential mapping exists as follows [25]:
$Φ=[Φ0Ṽ cos(Σ̃)+Ũ sin(Σ̃)]ṼT$
(8)
where $Φ0$ defines the bases at the point of tangency as a reference point. Note that one has $Φ̃=ŨΣ̃ṼT$ with $Φ̃$ being the bases at a point on T. Conversely, an explicit logarithmic mapping from the Grassmann manifold to a tangent space for points in the neighborhood of the point of tangency exists as well and is computed as follows [23]:
$Φ̃=Ũ tan−1(Σ̃) ṼT$
(9)
where
$(I−Φ0Φ0T)Φ(Φ0TΦ)−1=ŨΣ̃ṼT$
(10)

and $Φ0$ defines the bases at the point of tangency, while $Φ$ is the original point on the manifold. Using these mapping formulae, we can map the sample set of POD modes onto a tangent space, perform interpolation to a new modeling parameter, and map the new point back onto the Grassmann manifold. Accordingly, the bases can be interpolated while retaining orthonormality, producing acceptable bases, as schematically shown in Fig. 1 [23]. The interpolation process for reduced order bases with a tangent space of the Grassmann manifold is summarized below:

Fig. 1
Fig. 1
Close modal
1. Define the sampling points $s0, s1, ⋯ snR$ associated with modeling parameters used for interpolation. Here, $s0$ is the reference point from which the tangent space will be constructed for interpolating the reduced order modes, and $nR$ is the number of additional samples.

2. Calculate the POD modes for each of the sampling points following the procedure outlined in Sec. 2.1 as:
$Φi=Φ(si)∈ℝn×nϕ for i=1,…,nR$
(11)
3. Map the POD bases onto the tangent space constructed at the reference point $s0$ as
$Φ̃i=Φ̃(si)=Ũi tan−1(Σ̃i) ṼiT for i=1,…,nR$
(12)

where $(I−Φ0Φ0T)Φi(Φ0TΦi)−1=ŨiΣ̃iṼiT$

1. Interpolate the projected POD bases at $s0, s1, ⋯ snR$ to obtain the one associated with the model parameters of interest $s$ as shown in Fig. 1
$Φ̃s=Φ̃(s)=interp(Φ̃0, Φ̃1, ⋯, Φ̃nR|s)$
(13)
2. Map the interpolated projected POD basis for parameter $s$ back onto the manifold as
$Φs=Φ(s)=[Φ0Ṽs cos(Σ̃s)+Ũs sin(Σ̃s)]ṼsT$
(14)

where $Φ̃s=ŨsΣ̃sṼsT$.

## 3 Model Order Reduction for High-Fidelity Off-Road Mobility Model

The adaptive POD MOR technique is applied to the computational algorithm for deformable tire and soil models in this section.

### 3.1 Finite Element Tire-Soil Interaction Model

#### 3.1.1 Flexible Tire Model.

The multilayer fiber-reinforced rubber structure of a pneumatic tire is modeled using shear deformable composite shell elements, which are fully integrated into the multibody vehicle dynamics simulation framework [26]. The equations of motion of the nonlinear FE tire model with the composite shell elements are written in terms of the nodal coordinates $et$ as
$Mte¨t=Qts+Qte$
(15)
where $Mt$ is the generalized constant mass matrix, $Qte$ is the generalized external force vector, and $Qts$ is the generalized elastic force vector. For the element with $m$ layers, $Qts$ is defined as
$Qts=−∑l=1m∫V0l(∂εl∂et)T∂Wl∂εldV0l$
(16)

where $W$ is the energy density function of each layer of the composite shell element incorporating both an orthotropic material model for fiber-reinforced rubber and hyperelastic rubber material. $εl$ is the compatible strain vector of the l-th layer. With the shear deformable composite shell elements, the nonlinear properties of each of the layers in the tread, shoulder, sidewall, and bead sections of the tire can be considered. In order to model tire tread blocks, many contact nodes are defined on the shell-element outer surface to describe the tread block geometry as shown in Fig. 2 [5]. These contact nodes are neither FE nodes nor FE quadrature points. These contact nodes are introduced to define contact forces on the tire tread and sidewall surfaces. By doing so, distributions of the contact pressure as well as the surface traction within the contact patch can be predicted while considering the tread block geometry without any additional finite elements for tread blocks. A penalty-based compliant contact force is defined at each contact node in the direction normal to the tread surface, while the tangential forces are calculated using a slip-dependent friction force model [26].

Fig. 2
Fig. 2
Close modal

#### 3.1.2 Deformable Soil Model.

Soil deformation is modeled with a continuum FE approach [5], which can be further extended to the multiscale soil model [11]. The equations of motion of FE soil elements can be written as
$Mse¨s=Qss+Qse$
(17)
where $Ms$ is the constant mass matrix, $Qse$ is the generalized external force vector considering contact forces with the tire model, and $Qss$ is the generalized internal soil force vector defined by
$Qss=−∫V(∂D∂e˙)Tτ dV$
(18)
where $D$ is the rate of deformation tensor and τ is the Kirchhoff stress tensor that is conjugate to the logarithmic strain tensor obtained by
$εL=12ln(B)=∑i=13ln(λi)bi⊗bi$
(19)
where $B=∑i=13(λi)2bi⊗bi$ for $λi$ ($i=1,2,3$) being the ith eigenvalue of the spatial stretch tensor obtained by the polar decomposition of the displacement gradient tensor $F$; and $bi$ is the principal direction of the stretch tensor. The stress responses of soil are modeled by phenomenological constitutive models in the FE approach. That is, nonlinear soil material behavior is modeled as an elastoplastic material experiencing large plastic deformation using the multiplicative finite strain plasticity theory [27]. To this end, the displacement gradient tensor $F$ is decomposed as $F=FeFp$, where $Fe$ and $Fp$ are the elastic and plastic portions of the displacement gradient tensor, respectively. This means that the elastic deformation and plastic deformation are coupled, which is appropriate as the large plastic strains are exhibited in off-road tire–soil interaction problems. To determine the plastic deformation, the capped Drucker–Prager is used, for which the yield functions are given by [27]
$Φa=J2(s)+ηp−ξcΦb=1b2(p−pt+a(εvp))2+(q(s)M)2−a(εvp)2}$
(20)

where $Φa$ is the function defining the conical yield surface, while $Φb$ is the function defining the elliptical cap yield surface based on a modified Cam–Clay yield criterion. In the conical yield surface function $Φa$, $J2(s)$ is the second stress invariant with the deviatoric stress tensor $s$; $p(τ)$ is the mean stress; $c$ is the cohesion; $ξ=1/3$; and $η=tan β/3$ with $β$ being the Drucker–Prager friction angle. In the elliptical cap yield surface function $Φb$, $pt$ is the tensile yield pressure; $q(s)=3J2(s)$; $M=3η$; and $a$ and $b$ are the cone and cap size parameters. The parameter $a$ is a function of the volumetric plastic strain $εvp$ to model the strain hardening characteristics. Thus, the relationship between the compaction pressure and the volumetric plastic strain needs to be experimentally identified by a one-dimensional soil consolidation test. The return mapping algorithm is used to determine the plastic strain at the FE quadrature points in every step [5]. It is worth noting that the phenomenological constitutive assumptions can be eliminated by introducing the multiscale modeling technique for complex granular materials [11].

In the tire–soil interaction simulation, the FE soil elements only in the vicinity of the tire are considered in the solution process by moving the FE solution domain according to the tire position, as shown in Fig. 2. By doing so, FE elements far from the deformation region in the current time-step can be excluded from the computation, and the degrees-of-freedom of the FE soil model can be kept constant regardless of the distance that the tire travels in various simulation scenarios. More details on the update scheme of the moving soil patches are found in the literature [12].

### 3.2 Reduction of the Coupled Equations.

By considering rigid vehicle components modeled as a multibody system subjected to joint constraints, the equations of motion of the off-road mobility model are written in the form of the coupled differential-algebraic equations as follows:
$Mrq̈r+CqrTλ=Qr(qr,et,q̇r,ėt) ⋯ Rigid vehicle components Mtët=Qt(qr,et,es,q̇r,ėt,ės,α) ⋯ Nonlinear FE tiresMsës=Qs(et,es,ėt,ės) ⋯ Deformable FE terrainC(qr,t)=0 ⋯ Joint constraints}$
(21)
where the subscript r denotes rigid bodies, t denotes FE tires, and s indicates FE soil. The set of differential-algebraic equations is integrated forward in time with variable stepsize control using an implicit time integrator. Since the tire shell elements are assembled with a rigid rim by spring and damper elements to account for compliance of the tire bead portion, no kinematic constraints are imposed on the FE tire and soil models. The deformable tire and soil models contribute the majority of the degrees-of–freedom; thus, we focus on reducing the flexible body components. The preceding equations are partitioned into the rigid-body components subjected to kinematic constraints and the flexible tire and soil components as
$Mrq̈r+CqrTλ=Qr(qr,qf,q̇r,q̇f)Mfq̈f=Qf(qr,qf,q̇r,q̇f)C(qr,t)=0}$
(22)
where $qf=[etTesT]T$, $Mf=diag(Mt,Ms)$, and $Qf=[QtTQsT]T$. Using the generalized-α method, one can discretize the equations in time to obtain the solution at time-step n +1 as
$fr=(1−αm)Mr,n+1+αq̈r,n+1+α+αmMr,n+αq̈r,n+α+(1−αf)(CqTλ−Qr)n+1+αf(CqTλ−Qr)n=0ff=(1−αm)Mf,n+1+αq̈f,n+1+α+αmMf,n+αq̈f,n+α−(1−αf)Qf,n+1−αfQf,n=0 g=1βh2C(qn+1,tn+1)=0}$
(23)
where $αm$ and $αf$ are numerical damping parameters. It is worth noting that the generalized mass matrix in Eq. (23) is evaluated at $tn+1+α$ to account for shift in time by $α$. That is, $Mn+1+α=M(qn+1+α)$ with $qn+1+α=qn+h(1+α)q̇n$ [28]. The preceding discretized equations are solved for the generalized coordinates $qn+1$ and Lagrange multipliers $λn+1$ at time-step n +1 using the Newton–Raphson method. To this end, the model order reduction is applied to flexible components (i.e., FE tire and soil models), which contribute the majority of the degrees-of-freedom of the tire–soil interaction model. That is, one has to solve the following reduced equations in every iteration step k:
$[Jrr(k)J¯rf(k)(1−αf)Cq(k)TJ¯fr(k)J¯ff(k)0(1/βh2)Cq(k)00][Δqr(k)Δξ(k)Δλ(k)]=−[fr(k)f¯f(k)g(k)]$
(24)

where $qf(k)=Φξ(k)∈ℝnf$ is defined by the POD mode matrix $Φ∈ℝnf×nϕ$ and the reduced order coordinates $ξ∈ℝnϕ$; $Jrr(k)=∂fr(k)/∂qr$; $J¯rf(k)=(∂fr(k)/∂qf)Φ$; $J¯fr(k)=ΦT(∂ff(k)/∂qr)$; $J¯ff(k)=ΦT(∂ff(k)/∂qf)Φ$; and $f¯f(k)=ΦTff(k)$. It is worth noting that the constraint Jacobian is not condensed by the POD modes as the tire and soil models are not subjected to kinematic constraints. It is known that the reduced order constraint Jacobian matrix (e.g., $CqΦ$) can become singular due to the rank deficiency [29]. Thus, special care needs to be exercised when the constraint Jacobian matrix is condensed, and the linearly dependent column vectors of the reduced order constraint Jacobian matrix and the associated Lagrange multipliers need to be eliminated in each iteration. The computational improvements by ROM are attributed to the decreased number of equations that need to be solved in every time step, which is the primary expected source of computational improvement with the POD model order reduction. In addition, the model order reduction eliminates higher frequency components from the equations of motion, thereby allowing for a larger stepsize to speed up the time integration process.

## 4 Single Tire Soil Bin Mobility Simulation

A soil bin mobility test scenario is considered to examine the predictive ability of the reduced order tire–soil interaction model generated by POD. The FOM is composed of a single flexible tire rolling over deformable dry soil at a prescribed constant velocity of 5 m/s, as shown in Fig. 2. The soil material behavior is modeled by the capped Drucker–Prager failure criterion. The tire, modeled with 1920 shear deformable composite shell elements, is inflated to a pressure of 230 kPa and subjected to a wheel load of 6 kN. It is dropped into the soil and rolls forward to travel a total of 3 m. The moving soil patch size is 1.6 m in length, 0.48 m in width, and 0.4 m in depth. More details on the test condition as well as the parameters for the FE tire and soil models are found in the literature [5]. There are a total of 34,417 equations to be solved, in which 12,000 equations come from the tire model and 22,395 equations come from the soil model. Clearly, the vast majority of the degrees-of-freedom are contributed by the FE tire and soil models in this example.

By collecting snapshots of the FOM simulation, the POD modes are generated for the ROM. Once the modes are available, the number of modes to be used in the ROM must be determined. This decision is a tradeoff between accuracy and computational efficiency. One criterion is the cumulative energy of the modes. Typically, the modes used should contribute around 99% of the total energy [24], as calculated by summing the associated singular values and dividing by the sum of all of the singular values. In Fig. 3, the cumulative energy contribution of the modes is plotted. Approximately 99% of the energy is accounted for by the first 50 modes. This suggests that at least 50 modes should be used to ensure good agreement between the FOM and ROM solutions. For each additional mode added, the results should become more accurate.

Fig. 3
Fig. 3
Close modal

The transient longitudinal and vertical tire forces predicted by the ROMs with 40, 80, and 160 modes are compared with those of the FOM in Figs. 46, respectively, while the stress distributions of the tire surface as well as the midplane of the soil patch are shown in Fig. 7. Note that the negative longitudinal tire force (Fx) means that the tire is subjected to resistance forces exerted by the soil deformation. The lateral tire force is zero due to a zero steering angle. It can be seen from Figs. 4 and 7(a) that 40 modes are not sufficient to produce accurate solutions. The vertical tire force fluctuates inaccurately, and the soil stress distribution underneath the rolling tire does not agree with that of the FOM. On the other hand, the accuracy is significantly improved as the number of modes is increased to 80 and 160. The tire forces and the soil stress distribution predicted by the ROM agree well with those of the FOM. As in the result of the cumulative energy contribution of POD modes, more than 50 modes are needed to produce the accurate solutions of the time-domain tire–soil interaction simulation. To evaluate the deviations of the tire forces of the ROM from those of the FOM, the cross plots of the tire forces for 40, 80, and 160 mode cases are shown in Figs. 8(a)8(c), respectively. In these figures, the exact match between the FOM and ROM tire forces is given if a data point lies on the 45-deg. diagonal line. It is observed from these figures that the deviations of Fx and Fz in the 80- and 160-mode ROMs, shown in Figs. 8(b) and 8(c), fall within a similar range. To quantify the accuracy of the tire forces, the root-mean-square error (RMSE) is evaluated, and the relationship between the error and the number of modes is presented in Fig. 9. It is observed from this figure that the ROM solution converges toward the full order solution as the number of modes increases. Furthermore, beyond around 80 modes, the results fall within a reasonable margin of error.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal
Fig. 9
Fig. 9
Close modal

The computational time of the FOM per 1.0-m traveling distance was 102 min (1.7 h), while those of the ROM with 80 and 160 modes were 50 min and 61 min, respectively. That is, the computational time was reduced by 51% with the 80-mode ROM, while that of the 160-mode ROM was 40%. A substantial computational time reduction was achieved with the POD MOR while retaining predictive accuracy. The rate of the computational time reduction increases as the number of modes decreases as shown in Fig. 10, where one needs 80 modes and higher to ensure the predictive accuracy as shown in Fig. 9. The computational time improvement is attributed to the decreased number of equations given by MOR, i.e., Eq. (24). It is worth noting here that the total number of iterations required in the implicit time integration was substantially decreased by MOR, as summarized in Table 1. This decrease was caused by the elimination of higher frequency components by MOR. These components contribute very little to the solution while increasing the numerical stiffness of the FE tire and soil models. Thus, MOR helped improve the numerical convergence of the iterative Newton–Raphson procedure and led to substantial computational time reduction.

Fig. 10
Fig. 10
Close modal
Table 1

Total number of iterations for the tire to travel 1.0 m distance

FOMROMROMROM
160 modes80 modes40 modes
Total number of iterations10,8426,5286,3266,158
FOMROMROMROM
160 modes80 modes40 modes
Total number of iterations10,8426,5286,3266,158

In these examples, the modeling parameters used to obtain snapshots of the FOM are assumed to be identical to those in the ROM simulation. The modes obtained with a set of parameters, however, would not be suited for a ROM with a different set of parameters. To better understand the accuracy of POD modes for changes in modeling parameters, the modes produced with the 6 kN and 7 kN wheel loads are applied to ROM simulations with the 6.5 kN wheel load while keeping the other parameters the same. The 160 modes are considered as this model is sufficiently accurate for predicting the tire forces and vertical rim displacement, which represents the soil sinkage considering the tire deflection.

Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal
Fig. 13
Fig. 13
Close modal

It is important to note here that the interpolation technique is not limited to any single variable. It can also be used for soil parameters, tire parameters, or other variables that could be considered in the design space. A multivariate interpolation can then be performed for multiple of these variables that are considered at a time. Accordingly, a FOM can be run at a strategic sampling of parameter values, and the discrete set of POD modes can be turned into a continuous range of predictable bases through the use of interpolation on the tangent space of the Grassmann manifold. Furthermore, the cumulative energy contributions of POD modes at these sampling points provide a good estimate for the number of adapted modes required for the simulation of interest.

## 5 Full Vehicle Simulation on Deformable Terrain

Since the basic performances of the POD-based ROM for the tire–soil interaction simulation were discussed in Sec. 4 with the soil bin mobility scenario, the ROM is applied to a full vehicle simulation in this section for use in vehicle mobility prediction. To this end, an off-road vehicle shown in Fig. 14 is considered. To leverage parallel computing, the high-dimensional off-road vehicle mobility model is partitioned into five solution domains, called subsystems, as shown in Fig. 14: a multibody vehicle subsystem and four deformable tire-soil subsystems, each of which uses a moving soil patch to describe soil deformation resulting from the interaction with the tire. By doing so, four tire-soil subsystems, the most computationally expensive components in the system, can be run in parallel with the co-simulation technique. Within each compute node allocated to the tire-soil subsystem, shared-memory parallelization is used to conduct the FE tire and soil simulation with openmp, while the four knuckle points of the vehicle are connected to those of the tire rims in the four separate tire-soil subsystems through the force–displacement co-simulation coupling with mpi, as shown in Fig. 14. The POD MOR algorithm can then be incorporated into each of the tire-soil subsystems in the same way as the single tire soil bin mobility model presented in Sec. 4. Such a hybrid mpi-openmp parallelization algorithm is used for computational speedup, along with the POD-based reduced-order tire–soil interaction models. To account for the multipass effect of soil deformation, FE soil coordinates/velocities and the soil plasticity parameters in the four soil patches are saved and used to keep updating terrain deformation in the global terrain domain [12]. That is, a master processor in the vehicle mpi node collects the current FE soil data from the four soil patches when the moving soil patches are updated. The updated FE soil data in the global terrain domain can then be extracted as restart information for the soil patches in the subsequent simulation. By doing so, the soil deformation produced by a front tire can be considered in the rear soil patch as the initial soil condition to account for the multipass effect.

Fig. 14
Fig. 14
Close modal

A total of 132,983 equations are used to construct the full vehicle model. However, they are separated into 315 equations for the vehicle subsystem, and four sets of 33,167 equations for each tire-soil subsystem, to which MOR is applied. In the simulation, the vehicle is dropped into the soil and travels at a constant velocity without traction. Using snapshots from the FOM, POD modes are generated for each of the four tire-soil subsystems. The right-front tire forces of the FOM are compared to those of ROM with 40 and 120 modes in Figs. 15 and 16, respectively. Whereas the 40-mode ROM can predict the vertical tire force accurately, the longitudinal tire force is oscillatory and deviates from the FOM result. Due to the effect of the vehicle suspension, the vertical tire force changes smoothly in time; thus, 40 modes appear to be sufficient to capture the vertical tire force variation in this example. Increasing the number of modes to 120 leads to a significant improvement in results, and the longitudinal tire forces of the FOM and ROM are in good agreement. This tendency is clearly observed in the cross plots of the tire forces with 40, 80, and 120 modes in Fig. 17. The tire force RMSE shown in Fig. 18 indicates that at least 120 modes are needed in this problem to predict the tire forces accurately. Since the longitudinal tire force describes the rolling resistance exerted by the soil deformation, accurate prediction of the longitudinal tire force implies a correct prediction of the soil deformation as shown Fig. 19, where the soil stress distribution of the 120-mode ROM agrees well with that of FOM. The computational time of the full-order vehicle model per 1.0-m traveling distance was 90 min, while that of the 120-mode ROM was 65 min (28% reduction). The computational time was effectively reduced with the POD MOR for the full vehicle model as well while retaining accuracy. The reduction rate is, however, not as high as the single tire soil bin mobility model presented in Sec. 4. This is because the maximum stepsize taken by the reduced order tire-soil model cannot exceed the predetermined co-simulation macrostepsize in the full vehicle model to ensure numerical stability. This resulted in a limited reduction in the number of iterations for this problem. More computational reduction can, therefore, be expected by improving the co-simulation algorithm to allow for a larger maximum stepsize taken by the ROM tire-soil subsystems.

Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal
Fig. 17
Fig. 17
Close modal
Fig. 18
Fig. 18
Close modal
Fig. 19
Fig. 19
Close modal

## 6 Summary and Conclusions

A high-fidelity tire–soil interaction model plays a critical role in predicting vehicle mobility on deformable terrain to enable reliable vehicle design as well as mission planning. However, considering complex large granular soil deformation with nonlinear material behavior leads to computationally intensive models, which has been a major stumbling block in the practical usability of these models for the sake of design processes such as optimization and stochastic analysis. In order to address this issue, a method of MOR with POD was explored for off-road mobility simulation to improve the computational efficiency while retaining as much of the accuracy as possible. The POD modes were used to condense the equations of motion of the nonlinear FE tire and soil models discretized in time by the implicit time integration scheme. Producing POD modes requires dynamic data from a FOM, which limits its use as a method of MOR because a FOM must first be run when every ROM is created. Therefore, a method of mode adaptation through interpolation on a tangent space of the Grassmann manifold was investigated to allow modes to be predicted for cases in which a FOM has not been run.

In order to demonstrate these capabilities, a single tire soil bin mobility model was used to examine the accuracy and computational time reduction of the POD-based reduced order tire–soil interaction models developed in this study. The reduced order models were shown to be effective at accurately recovering the results of the FOM with as few as 80 modes in the single tire soil bin mobility model, while the number of degrees-of-freedom of the full-order tire and soil models was 34,395. This resulted in a computational time reduction of 51% while retaining the accuracy of the dynamic contact forces exerted on the deformable tire and soil. The computational improvements brought by the ROM had two main sources: a smaller system of equations in the solution procedure and a larger permissible maximum time-step. This larger maximum time-step is attributed to the fact that much of the high frequencies in the FOM are filtered out through the truncation of the low energy modes, which is beneficial to the time integration. Furthermore, the method of interpolation on a tangent space of the Grassmann manifold was employed using modes from other cases to predict modes for a case for which POD modes had not been produced. The predicted modes were shown to be much more accurate at predicting the sinkage of the tire in the soil and the tire forces. Additionally, the numerical procedure for model order reduction of the tire–soil interaction simulation was further extended to a full-scale vehicle model, composed of interconnected multibody vehicle components and the nonlinear FE tire–soil interaction models through the co-simulation coupling algorithm with high-performance computing techniques. In the future work, more complex loading scenarios, such as cornering and braking, will be examined in order to verify that POD modes are appropriate in these scenarios. Furthermore, the POD-based reduce order modeling will be applied to the upper-scale FE terrain model in the hierarchical multiscale off-road mobility model for computational speedup.

## Acknowledgment

This research was supported by the NASA Iowa Space Grant (Award No. NNX16AL88H; Funder ID: 10.13039/100000104) and in part by the Automotive Research Center (ARC) in accordance with U.S. Army CCDC Ground Vehicle Systems Center (GVSC) (Cooperative Agreement W56HZV-19-2-0001; Funder ID: 10.13039/100008192).

## References

1.
Bekker
,
M. G.
,
1956
,
Theory of Land Locomotion
,
University of Michigan Press
, Ann Arbor, MI.
2.
Wong
,
J. Y.
,
2008
,
Theory of Ground Vehicles
,
Wiley
, Hoboken, NJ.
3.
Krenn
,
R.
, and
Gibbesch
,
A.
,
2011
, “
Soft Soil Contact Modeling Technique for Multi-Body System Simulation
,”
Trends in Computational Contact Mechanics
,
Springer
,
Berlin, Heidelberg
, pp.
135
155
.
4.
Shoop
,
S. A.
,
2001
,
Finite Element Modeling of Tire-Terrain Interaction (No. ERDC/CRREL-TR-01-16)
,
Engineering Research and Development Center, Cold Regions Research and Engineering Lab.
,
Hanover NH
.
5.
Yamashita
,
H.
,
Jayakumar
,
P.
,
Alsaleh
,
M.
, and
Sugiyama
,
H.
,
2018
, “
Physics-Based Deformable Tire–Soil Interaction Model for Off-Road Mobility Simulation and Experimental Validation
,”
ASME J. Comput. Nonlinear Dyn.
,
13
(
2
), p.
021002
.10.1115/1.4037994
6.
Smith
,
W.
, and
Peng
,
H.
,
2013
, “
Modeling of Wheel–Soil Interaction Over Rough Terrain Using the Discrete Element Method
,”
J. Terramech.
,
50
(
5–6
), pp.
277
287
.10.1016/j.jterra.2013.09.002
7.
Recuero
,
A.
,
Serban
,
R.
,
Peterson
,
B.
,
Sugiyama
,
H.
,
Jayakumar
,
P.
, and
Negrut
,
D.
,
2017
, “
A High-Fidelity Approach for Vehicle Mobility Simulation: Nonlinear Finite Element Tires Operating on Granular Material
,”
J. Terramech.
,
72
, pp.
39
54
.10.1016/j.jterra.2017.04.002
8.
Wasfy
,
T. M.
,
Jayakumar
,
P.
,
Mechergui
,
D.
, and
Sanikommu
,
S.
,
2018
, “
Prediction of Vehicle Mobility on Large-Scale Soft-Soil Terrain Maps Using Physics-Based Simulation
,”
Int. J. Veh. Perform.
,
4
(
4
), pp.
347
381
.10.1504/IJVP.2018.095753
9.
Negrut
,
D.
,
Serban
,
R.
,
Mazhar
,
H.
, and
Heyn
,
T.
,
2014
, “
Parallel Computing in Multibody System Dynamics: Why, When, and How
,”
ASME J. Comput. Nonlinear Dyn.
,
9
(
4
), p.
041007
.10.1115/1.4027313
10.
Corona
,
E.
,
Gorsich
,
D.
,
Jayakumar
,
P.
, and
Veerapaneni
,
S.
,
2019
, “
Tensor Train Accelerated Solvers for Nonsmooth Rigid Body Dynamics
,”
ASME Appl. Mech. Rev.
,
71
(
5
), p. 050804. 10.1115/1.4043324
11.
Yamashita
,
H.
,
Chen
,
G.
,
Ruan
,
Y.
,
Jayakumar
,
P.
, and
Sugiyama
,
H.
,
2019
, “
Hierarchical Multiscale Modeling of Tire–Soil Interaction for Off-Road Mobility Simulation
,”
ASME J. Comput. Nonlinear Dyn.
,
14
(
6
), p.
061007
.10.1115/1.4042510
12.
Yamashita
,
H.
,
Chen
,
G.
,
Ruan
,
Y.
,
Jayakumar
,
P.
, and
Sugiyama
,
H.
,
2020
, “
Parallelized Multiscale Off-Road Vehicle Mobility Simulation Algorithm and Full-Scale Vehicle Validation
,”
ASME J. Comput. Nonlinear Dyn.
,
15
(
9
), p.
091007
.10.1115/1.4046666
13.
Chen
,
G.
,
Yamashita
,
H.
,
Ruan
,
Y.
,
Jayakumar
,
P.
,
Knap
,
J.
,
Leiter
,
K. W.
,
Yang
,
X.
, and
Sugiyama
,
H.
,
2021
, “
Enhancing Hierarchical Multiscale Off-Road Mobility Model by Neural Network Surrogate Model
,”
ASME J. Comput. Nonlinear Dyn.
,
16
(
8
), p.
081005
.10.1115/1.4051271
14.
Schilders
,
W.
,
2008
, “
Introduction to Model Order Reduction
,”
Model Order Reduction: Theory, Research Aspects and Applications
,
Springer
,
Berlin, Heidelberg, Germany
, pp.
3
32
.
15.
Lassila
,
T.
,
Manzoni
,
A.
,
Quarteroni
,
A.
, and
Rozza
,
G.
,
2014
, “
Model Order Reduction in Fluid Dynamics: Challenges and Perspectives
,”
Reduced Order Methods for Modeling and Computational Reduction
,
Springer
,
Cham, Switzerland
, pp.
235
273
.
16.
Berkooz
,
G.
,
Holmes
,
P.
, and
Lumley
,
J. L.
,
1993
, “
The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows
,”
Annu. Review Fluid Mech.
,
25
(
1
), pp.
539
575
.10.1146/annurev.fl.25.010193.002543
17.
Moehlis
,
J.
,
Smith
,
T. R.
,
Holmes
,
P.
, and
Faisst
,
H.
,
2002
, “
Models for Turbulent Plane Couette Flow Using the Proper Orthogonal Decomposition
,”
Phys. Fluids
,
14
(
7
), pp.
2493
2507
.10.1063/1.1483300
18.
Stabile
,
G.
,
Hijazi
,
S.
,
Mola
,
A.
,
Lorenzi
,
S.
, and
Rozza
,
G.
,
2017
, “
POD-Galerkin Reduced Order Methods for CFD Using Finite Volume Discretisation: Vortex Shedding Around a Circular Cylinder
,”
Commun. Appl. Ind. Math.
,
8
(
1
), pp.
210
236
.10.1515/caim-2017-0011
19.
LeGresley
,
P.
, and
Alonso
,
J.
,
2001
, “
Investigation of Non-Linear Projection for POD Based Reduced Order Models for Aerodynamics
,”
39th Aerospace Sciences Meeting and Exhibit
, Reno, NV, Apr., p.
926
.
20.
Weickum
,
G.
,
Eldred
,
M. S.
, and
Maute
,
K.
,
2009
, “
A Multi-Point Reduced-Order Modeling Approach of Transient Structural Dynamics With Application to Robust Design Optimization
,”
Struct. Multidiscip. Optim.
,
38
(
6
), pp.
599
611
.10.1007/s00158-008-0309-5
21.
Hesthaven
,
J. S.
, and
Ubbiali
,
S.
,
2018
, “
Non-Intrusive Reduced Order Modeling of Nonlinear Problems Using Neural Networks
,”
J. Comput. Phys.
,
363
, pp.
55
78
.10.1016/j.jcp.2018.02.037
22.
,
F.
,
Helenbrook
,
B. T.
, and
,
G.
,
2018
, “
Multilevel Algorithm for Obtaining the Proper Orthogonal Decomposition
,”
AIAA J.
,
56
(
11
), pp.
4423
4436
.10.2514/1.J056807
23.
Amsallem
,
D.
, and
Farhat
,
C.
,
2008
, “
Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity
,”
AIAA J.
,
46
(
7
), pp.
1803
1813
.10.2514/1.35374
24.
Kerschen
,
G.
,
Golinval
,
J. C.
,
Vakakis
,
A. F.
, and
Bergman
,
L. A.
,
2005
, “
The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems: An Overview
,”
Nonlinear Dyn.
,
41
(
1–3
), pp.
147
169
.10.1007/s11071-005-2803-2
25.
Lu
,
Y.
,
Blal
,
N.
, and
Gravouil
,
A.
,
2018
, “
Space–Time POD Based Computational Vademecums for Parametric Studies: Application to Thermo-Mechanical Problems
,”
,
5
(
1
), p.
3
.10.1186/s40323-018-0095-6
26.
Yamashita
,
H.
,
Jayakumar
,
P.
, and
Sugiyama
,
H.
,
2016
, “
Physics-Based Flexible Tire Model Integrated With LuGre Tire Friction for Transient Braking and Cornering Analysis
,”
ASME J. Comput. Nonlinear Dyn.
,
11
(
3
), p.
031017
.10.1115/1.4032855
27.
de Souza Neto
,
E. A.
,
Peric
,
D.
, and
Owen
,
D. R.
,
2011
,
Computational Methods for Plasticity: Theory and Applications
,
Wiley
, Hoboken, NJ.
28.
Jay
,
L. O.
,
Merwine
,
B. C.
,
Sugiyama
,
H.
, and
Yamashita
,
H.
,
2020
, “
A Two-Stage Extension of the Generalized-α Method for Constrained Systems in Mechanics
,”
ASME
Paper No. DETC2020-22385.10.1115/DETC2020-22385
29.
Luo
,
K.
,
Hu
,
H.
,
Liu
,
C.
, and
Tian
,
Q.
,
2017
, “
Model Order Reduction for Dynamic Simulation of a Flexible Multibody System Via Absolute Nodal Coordinate Formulation
,”
Comput. Methods Appl. Mech. Eng.
,
324
, pp.
573
594
.10.1016/j.cma.2017.06.029