Reduced Order Modeling of Deformable Tire-Soil Interaction With Proper Orthogonal Decomposition

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 [6–8] 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.

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.

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 $Φ$⁠.

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:

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)=Ũi tan−1(Σ̃i) ṼiT for i=1,…,nR$

   (12)

where $(I−Φ0Φ0T)Φi(Φ0TΦi)−1=ŨiΣ̃iṼ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)=[Φ0Ṽs cos(Σ̃s)+Ũs sin(Σ̃s)]ṼsT$

   (14)

where $Φ̃s=ŨsΣ̃sṼsT$⁠.

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

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].

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].

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.

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.

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. 4–6, 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.

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.

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.

The vertical tire rim displacements under the 6.5 kN wheel load predicted with the 6 and 7 kN POD modes are compared with the FOM result in Fig. 11. It is observed from this figure that the vertical rim displacement is underpredicted with the 6 kN modes due to the modes generated with a lower wheel load and overpredicted with the 7 kN modes due to the modes generated with a higher wheel load. A similar tendency is observed for the tire forces shown by the cross plots in Figs. 12(a) and 12(b), where Fz results of the nonadapted 6-kN modes are distributed below the 45-deg line, while those of the 7-kN modes are distributed above the 45-deg line. To address this issue, the POD modes for the 6.5 kN wheel load are produced by the adaptive POD method as explained in Sec. 2. That is, the modes for a 6.5 kN wheel load can be generated using the 6 kN modes as the point of tangency of the Grassmann manifold and the 7 kN modes as a sample point. Using the adapted POD modes for a 6.5 kN ROM, the dynamic simulation is conducted. The vertical tire displacement is compared with that of FOM in Fig. 13, while the cross plot of the tire forces is presented in Fig. 12(c). It is clear from these figures that the vertical rim displacement and the tire forces are predicted much more accurately with the adapted modes than either of the nonadapted modes. The tire force data of the adapted ROM lie on or around the 45-deg line, corresponding to the FOM solutions as shown in Fig. 12(c). Furthermore, it is observed from Fig. 13 that the adapted modes obtained by interpolating over a 3 kN range using the 5-kN and 8-kN modes still produce an accurate result. However, an interpolation range of 5 kN using the 4-kN and 9-kN modes appears to be too large to obtain accurate adaptive modes for this example, and sampling points closer to the one being predicted are required.

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.

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.

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.

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.

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).