基于地形变形的土壤力学模型在车轮运动模式中的应用：驱动、制动和侧滑

《Journal of Terramechanics》：Application of a terrain-deformation-based terramechanics model for wheel mobility modes: Driving, braking, and lateral slip

　　

# 摘要

探测车的开发和运行对于推进月球和行星探测至关重要。然而，在对应于地外重力和软质地形条件下评估车轮和车辆原型是极具挑战性的。因此，需要通过数值模拟系统性地评估车辆的移动特性。本研究提出了一种基于地形变形的土力学模型，用于分析车轮移动模式，包括驱动、制动和侧滑行为。所提出的模型考虑了土壤变形对轮-土相互作用的影响，从而能够更准确地预测各种车轮移动条件下的通行能力。为了验证该模型，进行了数值模拟并与在可变形地形上的单轮实验结果进行了比较。比较结果表明，所提出的模型有效地捕捉了地形变形对车轮性能的影响，特别是在制动和侧滑场景中，并提供了车轮移动性能的合理预测。该模型的另一显著特点在于其能够使用单一参数集处理多种移动条件。因此，我们成功地捕捉了车轮移动模式从驱动到制动以及从制动到驱动之间的过渡变化，并分析了由此产生的行驶行为。

## 1. 引言

土力学是一门研究机器与地面相互作用的交叉学科，在采矿场、灾害现场和地外地形等多种极端环境中都有应用。特别是在行星表面探测项目中，由于在局部环境中提前进行机器测试十分困难，基于土力学的模拟变得越来越重要。

正如多个航天机构发布的路线图所报告的那样（ISECG技术工作组，2019；Ishihara等，2024；Kawakatsu等，2023；NASA，2024），除着陆器的着陆和探测车的行驶外，即将开展的太空探索还将包括地外天体表面的许多操作，如钻孔、挖掘和地面采样。加压探测车的开发和运行对于由NASA主导的国际太空探索计划"阿尔忒弥斯计划"至关重要。此外，JAXA的火卫一月球样本返回任务（MMX）将搭载由CNES开发的探测车（Michel等，2022）。这些固态天体具有低重力和高真空条件，使得探测车在软质月壤上的稳定运行极具挑战性。

如前所述，数值模拟技术在探测车开发中起着重要作用。在土力学领域，基于Bekker、Wong和Reece的开创性工作（Bekker，1969；Wong和Reece，1967a，Wong和Reece，1967b；Wong，2008），半经验土力学模型已被应用于评估从小型机器人到大型矿用自卸车等各种规模的非公路车辆的行驶性能。

迄今为止，除应力评估的衍生模型外（例如：He等，2019；Horiko和Ishigami，2020；Shibly等，2005；Tsubaki和Ishigami，2021），还提出了各种扩展模型以适应制动和侧滑行为（例如：Ishigami等，2007，Ishigami等，2009；Gao等，2013；Guo等，2017；Karpman等，2024；Pavlov和Johnson，2024；Senatore和Sandu，2011）。此外，通过在多体动力学模型中实施土力学原理，已经对各种车辆模型（包括探测车）的整体行为进行了分析（例如：Arvidson等，2017；Ishigami等，2007，Ishigami等，2009；Li等，2013a；Zhou等，2014）。近年来，基于土力学的模拟已被用于利用其在短时间内生成各种条件下驾驶数据的能力来为人工智能（或机器学习）准备训练数据（Fediukov等，2024）。

传统土力学模型参考板桩贯入试验和直接剪切试验的结果，使用半经验公式根据滑转率评估轮-地接触面上的应力分布。它们基于力和力矩的平衡来预测牵引力、扭矩和沉陷等行驶性能。然而，传统土力学模型无法准确评估与车轮相互作用相关的地形表面变形。它们基于静态平衡，因此难以预测瞬态响应行为。因此，准确评估与车轮卡滞密切相关的牵引力和沉陷变化是具有挑战性的。

在这种情况下，提出了结合车轮相互作用引起的地形表面变形的模型（Guo等，2017；Jaiswal等，2019；Pavlov和Johnson，2024；Yu等，2024）。Suzuki等（2023）提出了一种扩展土力学模型（xTerramechanics模型），该模型基于元胞自动机原理纳入了地形表面的变形机制。此外，当xTerramechanics模型应用于单轮行驶分析时，它比传统土力学模型能够预测更真实的牵引力和沉陷值。还表明该模型可以表达与行驶模式相关的瞬态响应现象和卡滞现象的差异，如强制滑转行驶和自由螺旋桨行驶（Watanabe等，2023）。xTerramechanics模型也可以应用于多体动力学分析，其在分析探测车爬坡行为和刚性障碍物越障行为方面的有效性已得到证明（Kita等，2024；Watanabe等，2023）。然而，xTerramechanics模型的应用范围仅限于直线行驶。为了在各种非公路车辆（包括探测车）的开发和运行中使用，进一步推广是必要的。

在本研究中，我们将xTerramechanics模型推广以适应多种车轮移动模式，即直线驱动、制动和侧滑。为展示所提出模型的有效性，我们使用室内实验装置进行了单轮行驶测试。具体而言，我们使用单个刚性车轮进行了强制滑转测试，除控制滑转率外还控制侧偏角，以重现驱动、制动和侧滑条件。随后，我们将作用于车轮的牵引力、侧向力和沉陷与所提出模型预测的值进行比较，并证明了其有用性。所提出模型的另一特点是能够使用单一参数集处理多种行驶条件。因此，我们成功地表示了车轮移动模式从驱动到制动以及从制动到驱动的过渡变化，并讨论了行驶行为。

## 2. 模型

本节描述xTerramechanics模型。该模型由车轮表面的应力评估模型和基于元胞自动机的地形表面变形模型组成。在本研究中，采用Bekker-Wong-Reece模型（BWR模型）作为直线驱动的应力评估模型。除模型概述外，我们还描述了已作为扩展功能引入的制动模型（Gao等，2013；Guo等，2017；Wong和Reece，1967b）和侧滑模型（Ishigami等，2007，Ishigami等，2009）。选择这些模型是因为其适用性已在实验和数值研究中得到证明，从而可以独立于应力评估框架来评估所提出的地形变形逻辑的有效性。随后，我们概述了地形表面变形模型（Suzuki等，2023；Watanabe等，2023）。

### 2.1. 驱动的应力评估模型

BWR模型中刚性车轮下方的法向应力和剪应力分布分别如图1(a)和1(b)所示（Wong和Reece，1967a；Wong，2008）。在该模型中，法向应力σ由以下公式给出：

$$\sigma(\theta)=\left(\frac{c}{k_c'}+\frac{\rho g l}{k_\varphi'}\right)\left(\frac{\lambda(\theta)}{l}\right)^n \tag{1}$$

其中c为土壤凝聚力；$k_c'$和$k_\varphi'$分别为受土壤凝聚力和内摩擦角影响的沉陷模量；n表示沉陷指数；ρ为土壤堆积密度；g为重力加速度；l为车轮接触长度。由于车轮表面的应力分布取决于各位置的沉陷量λ(θ)，法向应力随接触角θ变化，λ(θ)由以下公式给出：

$$\lambda(\theta)=\begin{cases} \frac{d}{2}(\cos\theta-\cos\theta_f) & \text{for } \theta_m<\theta<\theta_f, \\ \frac{d}{2}\left[\cos\left\{\theta_f-\theta-\frac{\theta_r}{\theta_m-\theta_r}(\theta_f-\theta_m)\right\}-\cos\theta_f\right] & \text{for } \theta_r<\theta<\theta_m. \end{cases} \tag{2}$$

其中d为刚性车轮的直径，θf和θr分别为进入角和退出角。θm表示法向应力达到最大值的角度，表达式如下：

$$\theta_m=(a_0+a_1 s)\theta_f \tag{3}$$

其中a0和a1为最大应力角常数，取决于轮-土相互作用，s为滑转率，定义如下（He等，2020）：

$$s=1-\frac{v_x}{(d/2)\omega} \tag{4}$$

其中ω和vx分别为车轮的角速度和线速度（图1）。

车轮表面纵向切线方向产生的剪应力τx（图1(b)）利用土壤的剪切强度和车轮表面与土壤之间的相对滑动位移jx通过以下公式获得：

$$\tau_x(\theta)=(c+\sigma\tan\varphi)\left(1-e^{-j_x/k_x}\right) \tag{5}$$

其中φ为土壤内摩擦角，kx为剪切变形模量。jx通过积分车轮表面相对滑动速度在接触时间上获得，可以表示为θ的函数如下（Wong和Reece，1967a）：

$$j_x(\theta)=\frac{d}{2}\left\{\theta_f-\theta-(1-s)(\sin\theta_f-\sin\theta)\right\} \tag{6}$$

最后，牵引力Fx和垂直力Fz可以评估如下：

$$F_x=\frac{d}{2}b\int_{\theta_f}^{\theta_r}(\tau\cos\theta-\sigma\sin\theta)d\theta \tag{7}$$

$$F_z=\frac{d}{2}b\int_{\theta_f}^{\theta_r}(\tau\sin\theta+\sigma\cos\theta)d\theta \tag{8}$$

其中b为车轮宽度。

### 2.2. 制动的应力评估模型

制动期间滑转率的定义与驱动期间不同。制动期间车轮表面的应力分布也与驱动时不同（Wong和Reece，1967b）。本节简要描述制动期间车轮的应力分布模型。

在本研究中，制动期间的滑转率使用以下定义（Gao等，2013）：

$$s=\frac{(d/2)\omega}{v_x}-1 \tag{9}$$

注意制动期间的滑转率取负值；因此它是ISTVS标准（He等，2020）中定义的?1倍。

图2显示了制动期间刚性车轮下方切向应力分布的示意图。制动和驱动期间应力分布的一个区别是剪应力的方向发生变化。图中θt表示车轮表面产生的剪应力方向反转且法向应力达到最大值的角度（过渡点）。在本研究中，我们采用Gao等（2013）和Guo等（2017）提供的最大应力角定义，如以下公式所示：

$$\theta_t=(c_0+c_1|s|)\theta_f \tag{10}$$

其中c0和c1为定义应力过渡点角度的参数。制动期间的法向应力通过在公式(1)、(2)中设置θm = θt来评估。

虽然剪应力可以使用用于驱动的公式(5)计算，但相对滑动位移jx在驱动和制动之间存在差异。jx由以下公式给出（Gao等，2013；Guo等，2017）：

$$j_x(\theta)=\begin{cases} \frac{d}{2}(\theta_f-\theta)\left\{1+\frac{KV\sin\theta_f-\sin\theta}{(1+s)(\theta_f-\theta)}\right\} & \text{for } \theta_t<\theta\leq\theta_f, \\ \frac{d}{2}(\theta_t-\theta)\left\{\frac{\sin\theta_t-\sin\theta}{(1+s)(\theta_t-\theta)}-1\right\} & \text{for } \theta_r\leq\theta\leq\theta_t, \end{cases} \tag{11}$$

其中

$$KV=\frac{\sin\theta_f-\sin\theta_t}{\theta_f-\theta_t}-(1+s) \tag{12}$$

最后，制动期间的牵引力Fx和垂直力Fz可以评估如下：

$$F_x=\frac{d}{2}b\left\{\int_{\theta_f}^{\theta_t}\tau_x\cos\theta\,d\theta-\int_{\theta_t}^{\theta_r}\tau_x\cos\theta\,d\theta-\int_{\theta_f}^{\theta_r}\sigma\sin\theta\,d\theta\right\} \tag{13}$$

$$F_z=\frac{d}{2}b\left\{\int_{\theta_f}^{\theta_t}\tau_x\sin\theta\,d\theta-\int_{\theta_t}^{\theta_r}\tau_x\sin\theta\,d\theta-\int_{\theta_f}^{\theta_r}\sigma\cos\theta\,d\theta\right\} \tag{14}$$

其中积分范围分段以反映公式(11)。

### 2.3. 侧滑的应力评估模型

当车辆转向或穿越斜坡时，其车轮可能产生侧滑。在这种情况下，除与车轮轴向旋转相对应的纵向力外，车轮上还会产生侧向力（侧力）。Ishigami等（2007）将侧力视为两种应力的合力：车轮下方的剪应力τy和作用于车轮侧面的土压力Rb，如图3(a)所示。本研究也采用了其侧力模型，因为本研究的对象车轮具有侧板。

图3(b)显示了发生侧滑时车轮的俯视图。

车轮侧偏角β由以下方程给出，该方程使用车轮局部坐标系(x, y, z)中的速度vx和vy（Ishigami et al., 2007, Ishigami et al., 2009）：

(15) β = tan?1(vy/vx)

除侧偏角外，θz是车轮相对于世界坐标系(X, Y, Z)的X方向的角度，如图3(b)所示。该角度也在后续描述的侧向滑移过程中地形表面变形模型中使用。

剪切应力由车轮表面与地面之间的滑动产生。当发生侧偏角时，需要同时考虑车轮局部坐标系中沿x和y方向的两个相对滑动位移jx和jy。在此情况下，侧向滑移过程中车轮表面每个位置上x方向和y方向的相对滑动速度(v'x(θ), v'y(θ))之间的角度γ(θ)由以下方程给出（图3(c)）：

(16) γ(θ) = atan2(v'y(θ), v'x(θ)) = atan2(V sin β, (d/2)ω{1 - (1-s) cos θ})

其中，V是行驶方向上的车轮速度（图3(b)）。

车轮局部坐标系中的相对滑动位移jx和jy通过对车轮-土壤接触时间内的相对滑动速度进行积分获得（Ishigami et al., 2007, Ishigami et al., 2009），其表达式为：

(17) jx(θ) = (d/2){θf - θ - (1-s)(sin θf - sin θ)}

(18) jy(θ) = (d/2)(θf - θ)(1-s) tan β

其中，方程(16)、(17)、(18)仅适用于纵向滑移率s为正（驱动状态）的侧向滑移情况。需要注意的是，当将这些方程应用于制动工况时，应使用方程(9)、(11)中给出的负滑移率定义。特别是，方程(16)和(18)中出现的(1-s)项应替换为1/(1+s)。此外，由于本研究的模拟仅限于与室内实验相对应的平面地形条件，因此假设车辆速度向量V仅由水平分量组成。一般而言，当在崎岖地形上行驶、通过坡道或发生俯仰运动时，V还可能包含垂直分量。然而，本研究未考虑这些效应。

合滑动位移j由以下公式给出：

(19) j(θ) = √(jx(θ)2 + jy(θ)2)

因此，车轮下方产生的剪切应力表示为：

(20) τ(θ) = (c + σ tan φ)(1 - e^(-j/kx))

此外，由方程(16)可得：

(21) τx(θ) = τ(θ) cos γ(θ), τy(θ) = τ(θ) sin γ(θ)

接下来我们解释土压力。本研究采用阻力理论（Li et al., 2013b）来评估作用在轮侧盖板上的土压力，因为该方法简单且具有可扩展性。在本研究针对的分析中，轮侧的平面方向始终平行于z方向。因此，土压力Rb由以下方程给出：

(22) Rb = α(π/2, 0) |z?(θ)|

其中，α是当攻角为π/2且侵入角为零时单位深度的水平应力（地面土压力系数），如Li et al., 2013b所报告。z?是接触角θ处的深度。

最后，侧向力由以下公式给出：

(23) Fy = (d/2)∫(θr到θf) τy dθ + ∫(A) Rb dA

其中，A是轮侧与土壤之间的接触面积。

2.4. 地形表面变形模型

xTerramechanics模型中地形表面变形的处理基于元胞自动机（Nishimori和Ouchi, 1993; Suzuki和Ozaki, 2020; Werner, 1995）。如图4所示，整个地形表面被离散化为包含各位置高度hi,j和体密度信息的单元。这里，下标i和j分别表示与世界坐标系方向（即X和Y方向）对应的单元编号。进入角θf和离开角θr根据车轮与单元之间的几何关系自动计算。然而，在本研究中，根据车轮与单元之间的几何关系计算的离开角大于由以下方程得到的角度，因此使用下方方程对离开角进行了修正（Wong和Reece, 1967a）：

(24) θr = -cos?1(1 - κλ(0)/(d/2))

下载：下载高清图片(222KB) 下载：下载完整尺寸图片 下载：下载高清图片(163KB) 下载：下载完整尺寸图片

图4. xTerramechanics模型中车轮-土壤相互作用的示意图，其中地形表面被离散化为包含高度hi,j的单元：(a) 刚性轮下方的水平应力分布及推动作用，(b) 压缩作用，(c) 滚动作用。

其中，κ是离开角系数。

每个单元高度的变化基于三种变形作用来定义。第一种是推动作用（推土和移除土壤），即土壤运动根据车轮表面上的应力分布发生。第二种是压缩作用，考虑由于土壤压实引起的土壤性质变化。第三种是滚动作用，稳定地形表面的形状。

2.4.1. 土壤推动作用

在推动作用中，土壤运动方向使用2.1驱动应力评估模型、2.2制动应力评估模型和2.3侧向滑移应力评估模型中描述的应力评估模型获得的应力分布来定义。图4(a)展示了驱动或制动时推动作用的示意图。单元中的土壤运动定义如下（Suzuki et al., 2023; Watanabe et al., 2023）：

(25) dh_mov,i,j = -Si,j, dh_mov,i+di,j+dj = R1·Si,j

其中，dh_mov,i,j表示由推动作用引起的单元高度变化。方程(25)表明，单元中的部分高度Si,j根据车轮表面上的应力分布移动到相邻单元。需要注意的是，土壤在运动过程中预期会被压缩或膨胀；因此，高度乘以比例R1。下标di和dj根据与每个单元位置相对应的应力方向确定。如图4(a)所示，如果x（=X）方向的应力为负（σ-x），则di为1，因为土壤被推向正方向。相反，在车轮后方的区域，应力变为正（σ+x），di为-1，因为车轮通过其剪切作用将土壤向后推动。这里，x方向的应力σx由法向应力和剪切应力转换而来。Si,j根据每个车轮动力学计算步骤中车轮对地形表面的穿透（数值重叠）δi,j来确定（Suzuki et al., 2023; Watanabe et al., 2023）。在土壤变形的每个计算步骤中，单元高度移动量Si,j通过将δi,j乘以比例R2来控制，即：

(26) Si,j = R2·δi,j

接下来，我们在发生侧向滑移时推广推动作用。图5展示了侧向滑移时地形表面变形的示意图。侧向滑移时推动作用的一个特点是，轮侧也会因土压力的产生而推土。此外，根据车轮行驶方向的不同，土壤不仅在X轴方向移动，还在Y轴方向移动。因此，接下来我们将推动作用分为两部分进行说明：(a) 由车轮底部引起的推动作用和(b) 由轮侧引起的推动作用。

下载：下载高清图片(485KB) 下载：下载完整尺寸图片

图5. 侧向滑移时土壤推动作用的示意图。

(a) 车轮底部引起的推动作用

车轮底部引起的土壤运动方向与车轮上应力σx的方向相反，且该方向始终沿车轮的纵向方向，与行驶方向无关，即与侧偏角β无关。因此，X轴方向和Y轴方向土壤移动量的比值取决于角度θz。如果X轴方向的单元高度移动量为SX,i,j，Y轴方向的为SY,i,j，则它们确定如下：

(27) SX,i,j = ((π/2 - |θz|)/(π/2))·Si,j, SY,i,j = (|θz|/(π/2))·Si,j

其中，Si,j = SX,i,j + SY,i,j，单元高度移动的总量等于直线行驶时（方程(26)）的值。目标单元接收的土壤量由以下方程确定，该方程类似于方程(25)：

(28) dh_mov,i+di,j = R1·SX,i,j, dh_mov,i,j+dj = R1·SY,i,j

其中，di和dj根据水平应力σx的方向取值为1或-1。

(b) 轮侧盖板引起的推动作用

当发生侧向滑移时，轮侧产生土压力，导致推土现象，土壤移动量随侧偏角β而变化。

首先，以图5中所示的车轮局部坐标系为参考，x轴方向的单元高度移动量定义为S?x,i,j，y轴方向的为S?y,i,j。它们确定如下：

(29) S?x,i,j = ((π/2 - |β|)/(π/2))·Si,j, S?y,i,j = (|β|/(π/2))·Si,j

其中，方程(29)中的Si,j仅对应位于车轮边缘下方的单元。

由于单元编号基于世界坐标系，方程(29)必须转换到该坐标系。如果X轴方向的单元高度移动量为S?X,i,j，Y轴方向的为S?Y,i,j，它们可以表示为：

(30) S?X,i,j = ((π/2 - |θz|)/(π/2))·S?x,i,j + (|θz|/(π/2))·S?y,i,j, S?Y,i,j = (|θz|/(π/2))·S?x,i,j + ((π/2 - |θz|)/(π/2))·S?y,i,j

目标单元接收的土壤量由以下方程确定，该方程类似于方程(25)和(28)：

(31) dh_mov,i+di,j = R1·S?X,i,j, dh_mov,i,j+dj = R1·S?Y,i,j

其中，di和dj根据土压力的方向取值为1或-1。

2.4.2. 土壤压缩作用

压缩作用代表车轮下方土壤压实的效果。当地面上形成明显车辙（如雪地车辙）时，这种效应尤为显著。该作用可以表示如下（图4(b)）（Suzuki et al., 2023; Watanabe et al., 2023）：

(32) dh_comp,i,j = -R3·δi,j

与推动作用类似，地形表面高度根据δi,j降低。R2和R3之和必须小于或等于1，以避免物理矛盾。

2.4.3. 土壤滚动作用

接下来，我们解释滚动作用，该作用表达土壤从高处流向低处的现象，并促进地形表面形状的稳定性（图4(c)）。滚动作用的有效性已在风成波纹形成分析中得到证实（Nishimori和Ouchi, 1993; Suzuki和Ozaki, 2020; Werner, 1995）。该作用表示如下：

(33) dh_roll,i,j = -rh,i,j + (r/6)·h_cross + (r/12)·h_diag

其中，dh_roll,i,j是某单元中地形表面高度的变化。根据方程(33)，单元中的土壤由于重力作用滚动到相邻单元，同时该单元接收来自外部的土壤流入。流出率由系数r乘以单元高度hi,j获得。此外，rh_cross/6从i和j方向的正交单元流入，rh_diag/12从对角方向流入。每个流入项的h_cross和h_diag定义如下：

(34) h_cross = h_i+1,j + h_i-1,j + h_i,j+1 + h_i,j-1, h_diag = h_i+1,j+1 + h_i+1,j-1 + h_i-1,j-1 + h_i-1,j+1

滚动作用假设仅在相邻单元坡度超过安息角φr时发生（Suzuki和Ozaki, 2020; Suzuki et al., 2023; Watanabe et al., 2023）。此外，假设车轮内部不发生土壤流动。

2.4.4. 单元高度和体密度的演化

三种作用引起的高度变化同时并依次应用于所有单元，下一时刻的地形表面高度h'_i,j确定如下：

(35) h_i,j = h_old,i,j + dh_mov,i,j + dh_comp,i,j

(36) h'_i,j = h_i,j + dh_roll,i,j

其中，h_old,i,j是在土壤运动和压缩作用发生之前的单元高度。滚动作用计算在一个车轮动力学分析时间步内反复执行，以确保地形表面稳定变形（Watanabe et al., 2023）。因此，如果dt是车轮动力学分析的时间增量，kmax是每秒的最大滚动次数，则必须满足dt × kmax ≥ 1。此外，推动和压缩作用的间隔通过设置车轮的水平或垂直移动增量dact和旋转增量dθact来定义。这里，dact = d × dθact/2。如果满足其中任何一个增量，则推动和压缩作用发生。图6展示了使用xTerramechanics模型进行数值分析的流程图。该流程图不仅适用于单轮行驶，也适用于使用漫游车模型的多体动力学分析（Watanabe et al., 2023; Kita et al., 2024）。

下载：下载高清图片(267KB) 下载：下载完整尺寸图片

图6.

使用 xTerramechanics 模型的数值分析流程。也考虑了由于变形引起的土壤体积密度变化，因为车轮下的法向应力取决于体积密度（公式 (1)）。在 xTerramechanics 模型中，体积密度的变化是基于各自单元格的体积应变 $\epsilon_v$ 评估的。质量守恒定律得出以下方程：
(37) $\rho = \rho_0 \exp(-\epsilon_v)$
其中 $\rho_0$ 表示单元格的初始（参考）体积密度。我们使用公式 (37) 来评估土壤变形下的当前体积密度 $\rho$。

### 3. 单轮行驶实验
为了确认所提模型的有效性，作为对比，我们使用刚性车轮进行了直线行驶和制动实验，以及侧滑实验。本部分介绍了实验装置和条件。

#### 3.1. 装置
单轮行驶实验使用图 7(a) 所示的实验架进行。该装置主要由土箱、框架单元和车轮单元组成。土箱中填充了足够的硅砂（Toyoura 砂，Toyoura Keiseki Kogyo Co., Ltd., 日本）用于实验。框架单元控制车轮的纵向移动。为了进行强制滑移行驶实验，车轮纵向速度通过滚珠丝杠控制，参考车轮角速度和目标滑移率。在每次行驶实验之前，准备表面以确保均匀的地面初始条件。手动用铲子松动沙子以均匀化体积密度，并使用沿土箱导轨引导的整平板来创建平坦表面。

下载：下载高分辨率图像 (716KB) 下载：下载全尺寸图像
**图 7.** 单轮行驶测试的实验装置：(a) 概览，(b) 车轮单元，(c) 用于侧滑角的车轮角度调整。

车轮单元的结构如图 7(b) 所示。实验中使用了一个刚性车轮，在旋转轴上安装了驱动电机和六分力传感器，使其能够电子控制车轮旋转的角速度并测量作用在车轮上的力和扭矩的大小。此外，为了测量沉陷量（这是评估行驶性能的重要指标），在车轮上方安装了一个激光位移计。该计器固定在上部框架上，测量其到附加在车轮单元上的板（随车轮运动同步移动）的距离。此外，使用重力加载系统来施加车轮载荷。通过调整重量，系统可在所需重量条件下运行。

图 7(c) 显示了侧滑实验期间的车轮配置。在本研究中，行驶过程中侧滑角固定，车轮角度在车轮组件的上连接点处固定。该装置允许以 5°为间隔设置侧滑角。由于六分力传感器安装在车轮单元上，外力和扭矩是在车轮的局部坐标系 (x, y, z) 中测量的。

由六分力传感器和激光位移计测量的时间序列数据传输到用于数据记录器的 PC，并使用 Dewesoft 软件（DEWE Japan Co., Ltd., 日本）记录。

#### 3.2. 实验条件
在所有行驶实验中，为了实现强制滑移条件，控制了车轮的旋转角速度和其向前方向的平移速度。为了评估滑移率对行驶行为的影响，车轮的角速度保持恒定，而每次行驶时改变平移速度。对于每种行驶条件，至少进行三次试验，并相应地组织数据。

车轮几何形状、重量和运动条件列于表 1 中。车轮由金属制成，其表面覆盖着由金属布制成的表面材料。土箱宽度约为车轮宽度的五倍，土壤深度远大于实验期间观察到的最大沉陷量。因此，预计侧壁和底部边界对应力分布、滑移线发展、车轮力和沉陷量的影响可以忽略不计。此外，土箱尺寸与车轮尺寸的比率与先前报告的实验室车轮 - 土壤相互作用研究中采用的比率相当 (Ishigami 等，2007; Gao 等，2013; Guo 等，2017; Horiko 和 Ishigami, 2020; Suzuki 等，2023)。

**表 1.** 车轮规格和运动条件的详细信息。
| 参数 | 符号 | 单位 | 数值 |
| :--- | :--- | :--- | :--- |
| 车轮直径 | $d$ | [m] | 0.34 |
| 车轮宽度 | $b$ | [m] | 0.09 |
| 车轮重量 | $W$ | [N] | 150 |
| 角速度 | $\omega$ | [deg/s] | 7.0 (行驶)
3.5 (制动) |
| (a) 直线行驶滑移率 | $s$ | [-] | 0，0.1，0.2，0.3，0.4，0.5 |
| (b) 制动滑移率 | $s$ | [-] | -0.1，-0.2，-0.3，-0.4 |
| (c) 侧滑滑移率 | $s$ | [-] | 0，0.1，0.3，0.5 |
| 侧滑角 | $\beta$ | [deg] | 5，10，15 |

以下各节描述了本研究中进行的三种实验的条件和目标：(a) 直线行驶，(b) 制动，和 (c) 侧滑。

**(a) 直线行驶实验**
为了确定第 2 节中描述的压力评估和表面变形模型的基本参数，在六个滑移率水平下进行了强制滑移实验。每次稳态行驶期间（车轮运动稳定时）记录的数据用于评估牵引力和沉陷量的平均值。

**(b) 制动实验**
在负滑移率条件下进行了强制滑移实验。在制动实验中，如公式 (9) 所示，目标滑移率向负方向增加，导致较高的平移速度。然而，由于实验装置中滚珠丝杠的最大进给速度限制，角速度设置为行驶实验中使用值的一半，即 3.5 deg./s，滑移率的下限设置为 -0.4。与行驶实验类似，评估了每种滑移率下稳态行驶期间的平均牵引力（制动力）和沉陷量。

**(c) 侧滑实验**
如图 1 所示，通过将四种滑移率与固定的侧滑角 5°、10°和 15°组合，在 12 种条件下进行了实验。与其他实验一样，评估了稳态行驶期间的平均牵引力和沉陷量。此外，测量了每种条件下的侧向力。

### 4. 分析模型
本节描述了单轮行驶的分析条件。使用 Simulink/Simscape?模拟动态系统。使用基于 MATLAB?的 Simscape 语言实现 xTerramechanics 模型。使用 ode4 求解器求解运动方程。请注意，xTerramechanics 模型的计算流程与图 6 中所示的相同。

#### 4.1. 分析条件
图 8 显示了场模型的示意图，该模型参考第 3 节中描述的单轮行驶实验创建。地形表面使用边长 $dx = dy = 5.0$ mm 的方形单元格进行离散化。单元格尺寸应根据要分析的车轮尺寸进行选择。对于本研究中使用的车轮，5 mm 的单元格尺寸提供了适当的接触检测，并允许表示足够平滑的地形变形。进一步细化单元格分辨率并未对预测的反作用力、沉陷量或地形变形产生显著变化。因此，考虑到计算成本和模型保真度之间的权衡，采用了 5 mm 的单元格尺寸。单元格的初始高度是均匀的，且 $h_{i,j} = 300$ mm。轮廓颜色表示地形表面的高度。使用绿色表示初始高度作为基线，随着高度增加颜色变为红色，随着高度降低变为蓝色。在强制滑移条件下也进行了数值分析。建立了足够的行驶长度 (700 mm) 以获得稳态行驶行为。

下载：下载高分辨率图像 (146KB) 下载：下载全尺寸图像
**图 8.** 场模型示意图：(a) 场尺寸，(b) 初始地形表面。

分析中使用的车轮尺寸和运动条件与实验中使用的一样，如表 1 所列。与实验类似，车轮侧面由盖板组成，未考虑土壤流入车轮。表 2、表 3 分别列出了压力评估模型和地形表面变形模型的参数。此外，我们将时间增量设置为 $dt = 0.002$ s，重力加速度设置为 $g = 9.81$ m/s2。

**表 2.** 压力评估模型参数。
| 参数 | 符号 | 单位 | 数值 |
| :--- | :--- | :--- | :--- |
| 粘聚力 | $c$ | [Pa] | 0 |
| 内摩擦角 | $\varphi$ | [°] | 32 |
| 初始体积密度 | $\rho_0$ | [kg/m3] | 1490 |
| 参数 | $k'_c$ | [-] | 0 |
| 参数 | $k'_\varphi$ | [-] | 55 |
| 指数 | $n$ | [-] | 1 |
| 参数 | $k_x$ | [m] | 0.015 |
| 参数 | $k_y$ | [m] | 0.015 |
| 参数 | $a_0$ | [-] | 0.3 |
| 参数 | $a_1$ | [-] | 0.5 |
| 参数 | $c_0$ | [-] | 0.1 |
| 参数 | $c_1$ | [-] | 0.8 |
| 参数 | $\kappa$ | [-] | 0.05 |

**表 3.** 土壤变形模型参数。
| 参数 | 符号 | 单位 | 数值 |
| :--- | :--- | :--- | :--- |
| 参数 | $r$ | [-] | 0.2 |
| 参数 | $R_1$ | [-] | 0.99 |
| 参数 | $R_2$ | [-] | 0.4 |
| 参数 | $R_3$ | [-] | 0.01 |
| 参数 | $\varphi_r$ | [deg] | 32 |
| 参数 | $k_{max}$ | [-] | 2000 |
| 参数 | $d_{act}$ | [m] | 0.0005 |

xTerramechanics 模型的参数是通过反复试验确定的，使得牵引力和沉陷量都能得到合理重现。本研究中使用的土壤是硅砂（Toyoura 砂）。土壤强度参数，即粘聚力 $c$ 和内摩擦角 $\varphi$，是从直接剪切试验结果中获得的。由于硅砂的粘聚力通常非常小，为简单起见，$c$ 假设为零。压力评估模型的参数 ($k'_\varphi, k_x, a_0$, 和 $a_1$) 是使用在行驶条件下获得的实验结果确定的。此外，指数 $n$ 固定为 1。地形表面变形参数最初设置为 Suzuki 等 (2023) 在对同种 Toyoura 砂进行平刃挖掘试验中确定的值。然后对这些参数进行了微调，以考虑覆盖车轮表面的金属布的影响并重现实验观察到的剩余车辙深度。

#### 4.2. 计算流程
在数值分析中，与实验设置一样，通过控制平移速度同时保持恒定的角速度来定义滑移率。单轮模型的计算程序如下 (图 6)：

(1) **车轮坐标检测。**
Simulink 环境中的传感器模块检测世界坐标系中的车轮几何形状。基于检测到的车轮位置、方向、平移速度和角速度，计算车轮局部坐标系中的滑移率和侧滑角，以及车轮相对于 X 轴的角度 $\theta_z$。

(2) **接触计算。**
使用步骤 (1) 中获得的位置和方向评估车轮与地形单元格之间的接触条件。这使得可以计算状态，如车轮的进入和退出角度，以及进入地形的压入深度 $\delta_{i,j}$。在存在侧滑角的情况下，接触计算是在车轮的侧表面（盖板）上进行的。这使得能够识别侧表面与地形之间的接触区域，以及涉及推土的地形单元格。

(3) **应力计算。**
基于计算的滑移率和侧滑角，使用第 2 节中描述的模型计算车轮上的应力分布。同时，通过计算接触单元格高度的变化来确定每个单元格的移动作用方向。

(4) **移动和压缩作用计算。**
使用步骤 (2) 中获得的压入深度和步骤 (3) 中确定的移动方向，通过更新与车轮接触的地形单元格及其周围单元格的高度来建模地形变形。此外，根据土壤条件考虑车轮下方的压缩作用。当满足车轮水平或垂直移动增量 $d_{act}$，或旋转增量 $d\theta_{act}$ 的标准时，执行这些作用。

(5) **滚落作用计算。**
当相邻单元格之间的高度差超过休止角时，滚落作用重复 $dt \times k_{max}$ 次，导致地形表面几何形状发生变化。该过程平滑并稳定了由移动和压缩作用引起的变形地形表面。

(6) **应变计算。**
计算由于压缩作用而高度发生变化的单元格的体积应变。基于计算的应变，也估计土壤体积密度的变化，并用于下一步的应力计算。

(7) **步骤 (1) 到 (6) 的迭代**
步骤 (1)–(6) 在模拟中的每个时间步重复。

### 5. 与实验的比较与讨论
#### 5.1. 直线行驶
首先，比较了直线行驶的模拟和实验结果。如第 3 节所述，用于随后制动和侧滑分析的基本参数是通过将 BWR 模型和土壤变形参数在适当范围内拟合到获得的实验数据来确定的。图 9 比较了直线行驶 ($s \ge 0$) 期间获得的实验和分析结果。

图 9(a) 和 9(b) 分别展示了稳态条件下牵引力与滑移率、沉陷量与滑移率之间的关系。实验数据是基于稳态行驶期间获得的数值进行平均的，误差棒表示测量值的离散程度。重复测试中观察到的相对较小的变化表明，在整个测试过程中实验条件是相当一致的。
下载：下载高分辨率图像 (166KB) 下载：下载完整尺寸图像
图 9. 直线行驶和制动情况下的实验与仿真结果对比。每个图表展示了强制滑移条件下实验和仿真获得的平均稳态值：(a) 牵引力与滑移率，(b) 沉陷量与滑移率。
如图 9 所示，仿真获得的牵引力和沉陷量与实验结果相当吻合。此外，实验数据和仿真数据均表明，沉陷量随滑移率的增加而增加，而牵引力则增加并在 s ≈ 0.4–0.5 处饱和。这种行为可归因于随着滑移率增加而产生的剪切应力增大，导致车轮将更多的沙土推向后方，从而增加沉陷量。关于牵引力，虽然较大的沉陷量有助于产生更高的剪切应力以辅助牵引，但过大的沉陷会导致车轮前方堆积沙土，从而增加行驶阻力。因此，滑移极限保持在 s ≈ 0.4–0.5 左右。在仿真中，与实验相比，牵引力随滑移率增加而饱和的过程更为平缓，这被认为是受应力评估模型的影响。由于车轮后方的剪切应力分布与牵引力的大小直接相关，采用 Tsubaki 和 Ishigami (2021) 提出的高斯型法向应力分布会抑制牵引力的增加。此外，采用了 Janosi–Hanamoto 模型来表示剪切应力 - 相对滑动 - 位移关系（公式 (5)）。使用其他考虑软化行为的模型可能会产生更接近实验数据的结果（例如 He 等，2019）。
图 10 说明了仿真中滑移率分别为 s = 0.1、0.3 和 0.5 时的直线行驶后的地形表面变形。该图清楚地显示了车轮前方有明显的沙土堆积。随着滑移率的增加，残留车辙变得更加明显；实验中也有类似的观察结果。这种现象发生的原因是较高的滑移率导致较低的前进速度，导致车轮每单位行驶长度的土壤移动更频繁。
下载：下载高分辨率图像 (259KB) 下载：下载完整尺寸图像
图 10. 直线行驶后的地形表面变形快照：(a) s = 0.1, (b) s = 0.3, (c) s = 0.5。等高线颜色表示地形表面高度。
5.2. 制动
随后，对比了制动期间的实验和仿真结果。图 9(a) 和 9(b) 分别展示了各种负滑移率下的牵引力（制动力）和沉陷量的稳态值。分析中使用了表 2、表 3 中列出的、通过直线行驶期间参数拟合确定的参数。然而，针对制动条件的特定参数 c0 和 c1 被额外拟合。c0 和 c1 的拟合值列在表 2 中。
如图 9 所示，牵引力和沉陷量在实验和仿真结果中均表现出相似的趋势，其中牵引力的吻合度尤为接近。制动行为的一个特征特征是，随着滑移率变得更负，牵引力在负方向上增加。这是由于制动期间的最大应力角比行驶期间向前移动，且随着滑移率绝对值的增加，向车轮前方移动的土壤（推土）变得更加明显。此外，沉陷量随更大的负滑移率增加与法向应力分布的前移有关。与在车轮后方产生应力的行驶条件不同，制动导致后部区域应力极小。因此，车轮正下方的法向应力减小。此外，车轮与地面之间的接触面积变窄，需要更大的沉陷量来支撑车轮的重量。
为了进一步量化实验和仿真结果之间的吻合度，评估了牵引力和沉陷量的均方根误差 (RMSE) 值，总结如表 4 所示。行驶条件下牵引力的 RMSE 为 4.9 N，制动条件下为 4.0 N。行驶中略大的 RMSE 反映了在较高滑移率下观察到的差异，而仿真在制动条件下更宽的滑移率范围内与实验结果吻合良好。在整个牵引力范围内，牵引力的归一化 RMSE 小于 5%。对于沉陷量，行驶条件下的 RMSE 为 1.9 mm，制动条件下为 1.6 mm。这些误差与实验变异性相当，相对于 17.2 mm 的总体沉陷范围来说很小，表明仿真和实验之间具有可接受的吻合度。
表 4. 实验与仿真结果之间牵引力和沉陷量的均方根误差 (RMSE)。
空单元格 行驶条件 制动条件
牵引力 [N] 4.9 4.0
沉陷量 [mm] 1.9 1.6
图 11(a) 和 11(b) 展示了滑移率分别为 s = -0.1 和 -0.3 时仿真获得的制动后地形表面变形。与图 10 所示的行驶期间地形表面变形相比，制动期间可以观察到车轮前方有更大的土壤堆积。随着负滑移率的增加，这种趋势变得更加明显；实验中观察到了相似的趋势。这种现象归因于随着滑移率变得更负，最大应力角进一步向车轮前方移动，且前进速度增加。因此，更多的沙土被推向前方。图 11(c) 展示了在负滑移率 -0.3 下进行的制动测试。仿真结果成功地重现了实验中观察到的关键现象，如车轮前方的推土现象和其后方的残留车辙形成。
下载：下载高分辨率图像 (345KB) 下载：下载完整尺寸图像
图 11. 制动后地形表面变形快照：(a) s = -0.1, (b) s = -0.3, (c) 实验期间地形表面变形 (s = -0.3)。等高线颜色表示地形表面高度。
尽管沉陷量得到了定量评估，但由于实验期间未获取三维地形剖面数据，因此未对地形变形模式（包括车辙几何形状和土壤堆积）进行定量比较。因此，基于观察到的地形变形模式，对这些变形特征的吻合度进行了定性评估。
除了所提出模型的预测能力外，其计算效率也是实际应用中的重要考虑因素。值得一提的是，xTerramechanics 模型比传统地面力学模型需要更高的计算成本，但仍远低于有限元法 (FEM) 或离散元法 (DEM) 的要求。本研究中的仿真是在一台配备 Intel Core i9-13,900 CPU、16 GB RAM、运行 Windows 11 的计算机上进行的。在当前分析条件下，模拟 700 mm 的行程分别需要大约 2000 秒和 5000 秒，对应滑移率 -0.3 和 0.3。这些运行时间包括生成大量输出数据（如可视化和动画结果），且 MATLAB?/Simscape 实现未针对计算效率进行优化。相比之下，忽略地形变形的传统地面力学模型在相同条件下需要大约三分之一到一半的计算时间。尽管因此所提出的模型比传统方法产生更高的计算成本，但它仍比高保真三维 FEM 或 DEM 仿真有效得多。虽然 FEM 和 DEM 可以容纳复杂的轮子几何形状并提供详细的应力 - 应变分析，但由于相关的计算费用，它们应用于长时间、整车仿真往往是不实际的。相反，FEM 和 DEM 可以作为 xTerramechanics 框架内参数识别和模型标定的互补工具。
5.3. 侧滑
对侧滑的实验和仿真结果进行了对比。图 12(a) 和 12(b) 分别展示了各种滑移率和侧滑角下的牵引力和沉陷量的稳态值。此外，图 12(c) 展示了相应条件下的侧向力。在这些图中，图例表示侧滑角 β，空心标记代表实验结果，实心标记代表仿真结果。此对比中使用的参数与表 2、表 3 中列出的相同，这些参数是通过直线行驶期间的参数拟合确定的。
下载：下载高分辨率图像 (216KB) 下载：下载完整尺寸图像
图 12. 侧滑情况下的实验与仿真结果对比。每个图表展示了强制滑移条件下实验和仿真获得的平均稳态值：(a) 牵引力与滑移率，(b) 沉陷量与滑移率，(c) 侧向力与滑移率。
图 12(a) 和 12(b) 表明，虽然确认了侧滑的影响，但可以确认与直线行驶期间观察到的牵引力、沉陷量与滑移率关系（如图 9(a) 和 9(b) 所示）相似的趋势。图 12(c) 显示侧向力随更大的侧滑角增加，而随滑移率增加倾向于减小。侧向力随更大侧滑角的增加可归因于两个因素：如公式 (16)、(21) 所示，作用在车轮侧向的剪切应力比例随侧滑角增加，且作用在车轮侧面的土压力也变得更为显著。相反，侧向力随滑移率增加而减小可能是因为，虽然较高的滑移率导致的更大沉陷使土压力增加，但剪切应力因前进速度降低而减小。换句话说，在贡献侧向力的两个组成部分——剪切应力 ???? 和土压力 Rb 中，剪切应力被认为是当前土壤条件下的主导因素。
图 13(a) 和 13(b) 说明了滑移率 s = 0.3 和侧滑角 β = 15° 时侧滑行驶后的地形变形。除了直线行驶期间观察到的车轮前方沙土堆积外，该图显示车轮侧面有明显的沙土堆积。这归因于 xTerramechanics 模型中引入的车轮侧表面造成的推土效应。图 13(c) 说明了滑移率 0.3 和侧滑角 15° 时的行驶条件。推土沿着车轮侧面发生，符合行驶状态，仿真结果成功重现了这一现象。随着滑移率和侧滑角的增加，侧向土壤位移变得更加明显，实验和仿真中都确认了类似的行为。
下载：下载高分辨率图像 (254KB) 下载：下载完整尺寸图像
图 13. 侧滑行驶后的地形表面变形快照 (s = 0.3, β = 15°)：(a) 侧视图，(b) 顶视图，(c) 实验期间地形表面变形。等高线颜色表示地形表面高度。
基于上述结果，所提出的模型在真实条件下侧滑行驶方面与实验数据表现出良好的吻合度，特别是在 s ≈ 0.3 时，牵引力、沉陷量和侧向力与实验值非常接近。地形表面变形在仿真和实验中也显示出相似的趋势。此外，本研究中分析的所有车轮机动模式一致使用了相同的参数集，进一步验证了所提出模型的有效性和通用性。
虽然本研究未进行正式的敏感性分析，但模型公式表明土壤特性（例如粘聚力 c 和内摩擦角 ??）和轮子参数对牵引力和侧向力的影响最强，这与传统地面力学模型一致。相比之下，沉陷量不仅受土壤刚度和轮子参数影响，还受地形变形参数 R1、R2 和 R3 的影响，这些参数控制车轮下方和周围地形变形的发展。此外，推土参数 R2 对制动性能和行驶阻力有显著影响，因为它直接影响车轮前方土壤堆积产生的阻力。
5.4. 行驶 - 制动转换的瞬态响应
图 9 综合对比了行驶模式和制动模式之间的牵引力和沉陷量。

如表2和表3所列，单一组参数的确定成功复现了驱动和制动条件下的实验结果，证实了xTerramechanics模型在捕捉驱动和制动行为方面的有效性。尽管由于实验装置控制能力的限制，无法直接验证瞬态驱动-制动响应的实验结果，但所提出的模型以可接受的精度复现了稳态驱动、制动和横向滑移行为。此外，Suzuki等人（2023）使用正弦地形表面证明了xTerramechanics模型在瞬态车轮-地形相互作用条件下的预测性能。综合来看，这些发现表明所提出的模型能够合理捕捉瞬态响应，包括土壤在车轮前面积聚的情况，尽管严格的验证仍需要专门的实验。

xTerramechanics模型的一个显著特点是其能够将历史效应和瞬态响应纳入运动行为中。因此，作为最终演示，我们展示了在驱动和制动之间反复转换期间的车轮运动行为。此处使用的车轮规格和模型参数与表1、表2和表3中列出的一致。

在仿真中，如图14(a)所示，滑移率使用正弦函数从0→?0.4→0.4→0连续变化。在此条件下，使用xTerramechanics模型获得的牵引力与滑移率的关系如图14(b)所示。图中空心标记表示当车轮被迫以恒定滑移率行驶时获得的稳态值（见图9）。在初始过程（i）中，随着滑移率变为负值，制动力增大。随后，在过程（ii）中，随着滑移率趋近于零，驱动力减小。当滑移率从负值变为正值时，车轮必须越过在其前方堆积的被推土，导致牵引力暂时增加。随后，牵引力随滑移率增加而增大，达到与恒定滑移率条件下观察到的相当的值。当滑移率返回零时（iii），观察到相同的牵引力-滑移率关系。

下载：下载高分辨率图像（272KB）下载：下载全尺寸图像

图14. 驱动-制动转换瞬态响应的演示：(a) 滑移率随时间的变化，(b) 牵引力随滑移率的变化，(c) 沉陷量随滑移率的变化。

在图14(c)所示的沉陷量-滑移率关系中，可以观察到先前运动模式的历史影响更为显著。在过程（i）中，沉陷量随滑移率的变化而变化。然而，由于车轮与土壤接触至行驶的时间较短，稳态条件尚未达到。在过程（ii）的前半部分，随着负滑移率的减小，沉陷量降低。在此阶段，稳态条件得以实现，其值与恒定滑移率下获得的值（图9）紧密对应。随后，随着正滑移率的增加，沉陷量增大，但增幅小于恒定滑移条件下的情况。这归因于制动期间车轮前方推土作用的残余效应。最后，在过程（iii）中，随着滑移率减小，沉陷量也随之减小，由于对先前历史的依赖性，呈现出滞后回线。

图中还包含了使用与表2中相同的参数集所得到的传统模型的仿真结果。由于传统模型无法考虑先前地形变形的影响，牵引力和沉陷量表现出对滑移率不真实的历史无关行为。需要注意的是，虽然传统模型和xTerramechanics模型在驱动期间的牵引力预测水平相似，但传统模型低估了制动期间的牵引力，因为它没有考虑车轮前方的推土作用。此外，两个模型预测的沉陷量存在显著差异。

5.5 假设与局限

本研究中考虑的实验仅限于强制滑移条件。因此，观察到的运动行为可能与在牵引负载条件下实际行驶的漫游器的行为有所不同。尽管先前的研究（Watanabe等人，2023）表明xTerramechanics框架能够表征强制滑移和自由推进行驶条件之间的差异，但在牵引负载条件下的进一步验证仍然是有益的。此外，实验中使用的车轮为金属材质，变形可忽略不计；因此，所提出的模型假设为刚体车轮。该框架可以通过采用替代的应力评估模型来扩展到变形车轮。

虽然分层地形可以通过指定单元高度来表征，但假设每个单元内的土壤性质是均匀的。因此，该模型无法考虑随深度变化的堆积密度或取决于含水量的土壤行为。此外，所提出的模型适用于低速应用，并忽略了动态车轮-土壤相互作用引起的土壤喷射。

基于细胞自动机的地形变形模型的数值稳定性通过适当选择仿真时间步长来维持。在本研究考虑的所有运动条件下均获得了稳定的结果。然而，长时间的卡滞状态（s = 1）可能导致过度的地形变形和预测反力的波动。因此，在将模型应用于长时间 immobilized 条件时需要谨慎。

虽然所提出的模型可以通过车轮载荷和车轮-土壤相互作用计算适应不同的重力环境，但其有效性仅在地面重力条件下得到验证。因此，该模型在月球和火星等低重力环境中的适用性以及重力缩放效应的影响仍有待进一步研究。

6. 结论

本研究致力于推广xTerramechanics模型，该模型先前通过考虑地形表面变形在预测直线行驶行为方面已证明了其有效性。具体而言，通过引入现有的应力评估模型并推广土壤移动逻辑，该模型被增强以适应制动和横向滑移。为了验证所提出的模型，与使用刚体车轮的室内行驶试验进行了比较。xTerramechanics模型在强制滑移条件下以合理的精度成功预测了牵引力、沉陷量和侧向力。此外，证明单一组参数可以应用于各种车轮运动模式。最后，为了突出模型反映地形变形历史的能力，展示了驱动-制动转换期间瞬态响应的仿真示例。这些结果表明，所提出的模型为分析车轮-地形相互作用提供了一个稳健且通用的框架，并为行星漫游器和越野车辆的设计和控制提供了有价值的见解。

在 future work 中，我们计划通过将预测与在多种地形上实际越野车辆获得的实验结果进行比较，进一步验证xTerramechanics模型。将模型扩展到柔性轮胎和带有齿块的轮子仍然是一项重要挑战，特别是当齿块引起的土壤变形变得复杂时。在这种情况下，可能需要进一步开发结合先进地形变形表征技术的模型。未来的工作还将调查所识别参数集向具有不同尺寸和负载条件的车轮的可迁移性。虽然一些参数，如c、φ、ρ、R1、R2和R3，主要表征土壤条件，但其余参数取决于本研究中考虑的车轮规格。此外，未来的工作还将研究重力缩放效应，并在与行星探索相关的低重力环境中（如月球和火星）验证所提出的模型。

数据声明

支持本研究发现的数据可在合理请求下从通讯作者处获得。

作者贡献声明

Tomoya Arai：撰写–初稿、方法论、数据管理。
Mai Shimizu：可视化、数据管理。
Kei Tsuchiya：调查。
Seiji Kon：调查。
Shingo Ozaki：撰写–审阅与编辑、资金获取、概念化。

资金

本研究由国土交通省"空间建设创新项目"支持。

Tomoya Arai | Mai Shimizu | Kei Tsuchiya | Seiji Kon | Shingo Ozaki