动态加载条件下土壤承压沉陷特性仿真研究

冯文选，马吉胜，何 健，吴大林

(1.陆军工程大学石家庄校区，石家庄 050003； 2.西京学院， 西安 710123)

车辆越野行驶过程中行驶机构与地面之间的力学相互作用分为垂直应力应变关系和水平应力应变关系。其中，垂直应力应变关系可以用车辆作用于地面垂直方向的载荷与相应载荷下地面沉陷量之间的关系来描述[1]。研究车辆—地面压力沉陷特性关系对研究车辆在特定地域的通过性和车辆的行驶阻力具有重要作用。

从德国学者Bernstein建立车轮陷入地面深度与接地压力之间的关系开始，各国学者均对车辆—地面压力沉陷关系进行了不同的研究。苏联学者Goriatchkin将Bernstein公式推广为指数形式。美国学者bekker在上述二人研究的基础上，结合土木工程中地基下陷规律提出了bekker公式，并设计出测试地面承压沉陷特性的试验方法和仪器。英国科学家Reece[2]考虑土体密度和黏聚力的影响，对bekker方程进行了改进。俄国科学家库兹可夫提出了双曲正切模型，日本科学家的双曲线模型[3—4]。俄罗斯学者Modest Lyasko在总结其他学者实验结果的基础上，提出了参数与实验条件无关的LSA模型[5]。

中国学者在土壤的承压沉陷特性方面也做了许多研究工作。吉林工业大学庄继德等人，以新疆沙漠沙为研究对象，提出了描述沙土压力沉陷特性的模型[6]。杨启梁等[7]在研究水田载荷—下陷特性中，采用形积当量作为反映测板形状、尺寸特征的参数，推导出了水田土壤压力沉陷公式。南京农业大学姚艳等人采用二项式对室内重塑土壤压力沉陷关系进行拟合，得到了较高的拟合精度，并进行了农田实测数据的验证[8]。原军械工程学院赵家丰等人结合土壤承载极限理论，提出一种改进的土壤承载模型，并利用不同土质实验参数进行验证，取得较好的拟合精度[9]。

总结国内外研究可以发现，各国学者对车辆—地面压力沉陷特性的研究主要建立在基于准静态加载的平板载荷试验之上。为更加贴近车辆行驶过程于地面相互作用的实际情况，本文采用有限元仿真的方法分析加载速率对车辆—地面压力沉陷特性的影响，并试图利用数值逼近的方法建立动态加载条件下车辆—地面压力沉陷关系模型。

1 土壤力学特性

1.1 土壤本构关系

土壤的应力—应变关系可以用固体力学中的弹塑性理论来描述[10]。土壤的力学行为可以分为弹性行为和塑性行为。弹性行为是土壤在弹性变形范围内的应力—应变关系，可以用弹性本构模型来描述。弹性模型主要包含两个参数，弹性模量E和泊松比μ。常见的土壤弹性本构模型包括线弹性模型、多孔介质弹性模型和线粘弹性模型。本文采用适用最广泛的各向同性弹性模型描述土壤的弹性行为。

土壤的塑性行为主要包括屈服、硬化、剪胀和流动性，可以用塑性本构模型来描述。经典的土壤塑性本构模型主要有Mohr-Coulomb模型、Drucker-Prager模型和临界状态塑性模型[11]。本文采用被广泛应用的Mohr-Coulomb模型描述土壤的塑性行为。Mohr-Coulomb模型中的剪切屈服函数为：

F=Rmcq-ptanφ-c=0

(1)

其中：φ为q-p平面上的摩擦角，c为黏聚力，Rmc为屈服函数在π平面上的形状参数。

在ABAQUS软件中，为了避免π平面上屈服面存在尖角导致流动方向不唯一的现象，Mohr-Coulomb模型采用非相关联的流动准则。并通过指定黏聚力c与等效塑性应变之间的关系，控制土壤的硬化规律。

1.2 土壤参数的确定

土壤本构参数可以通过土工试验来确定。土体抗剪强度试验是测定土体剪切特性的试验的一种，分为室内试验和现场试验。室内试验包括直剪试验和三轴剪切试验，直剪试验仪器设备简单、操作方便、试验用土少，可以测量土壤的抗剪强度参数粘聚力c和内摩擦角φ。三轴压缩试验作为另一种常见的室内土工试验，与直剪试验相比，除了可以测定抗剪强度指标粘聚力c和内摩擦角φ，还可以测定土壤的弹性模量E、泊松比μ和压缩系数等参数[12]。

采用上述两种土工试验对松软粘土进行参数测定[13]，在含水率5%的情况下，测得土壤内摩擦角φ=27.33°，粘聚力c=6.38 kPa。通过三轴压缩试验测得土壤弹性模量E=20.2 MPa，泊松比μ=0.32。总密度为1 932 kg/m3。

2 土壤压力沉陷试验有限元分析

2.1 有限元模型建立

有限元仿真采用大型非线性有限元分析软件ABAQUS为仿真平台，采取三维实体建模。载荷板为半径R=0.1 m的圆板。为简化计算，载荷板采用刚体约束，并设置参考点。综合考虑土体模型边界对仿真结果的影响和仿真计算效率，土体采用1 m×1 m×1 m的正方体模型。建立材料分别为钢和土壤的截面属性，钢的材料参数为：密度7 800 kg/m3、弹性模量E=207 GPa、泊松比μ=0.3；土壤本构模型采用各向同性线弹性模型和Mohr-Coulomb模型，参数如1.2所述。

约束土体模型四个侧面水平方向位移自由度，底面设置为完全固定约束；约束载荷板参考点水平方向位移自由度和三个方向转动自由度。整个模型采用八节点线性六面体单元，设置减缩积分、沙漏控制。土体网格拓扑结构采用中间致密四周相对稀疏的方法。设置载荷板参考点竖直方向载荷。装配体模型如图1所示。

2.2 动、静态加载仿真结果对比

ABAQUS软件提供的Standard求解器可以忽略部件惯性效应的影响，对稳态过程和准静态过程等静力学问题具有良好的模拟效果；Explict求解器对时间、质量惯性等动态因素敏感，可以模拟动力学问题[14]。基于上述两种求解器分别对准静态加载和动态加载条件下的土壤压力沉陷特性进行分析。

在2.1节模型的基础上，设置载荷板参考点竖直方向载荷大小50 kN，设置静力通用分析步模拟土壤压力沉陷准静态加载。设置载荷大小不变，加载速率为100 kPa/s，采用显示动力分析步模拟土壤压力沉陷动态加载条件。取载荷板参考点竖直方向的位移、载荷数据，并采用经典土壤压力沉陷公式p=kzn对仿真数据进行拟合，拟合曲线如图2所示。

从图2中可得：

1) 对比准静态加载与动态加载条件下土壤压力沉陷曲线，可以看出土壤在两种加载条件下的压力沉陷特性存在明显差异，特别是土壤进入塑性变形之后，动态加载时土壤的极限承载力明显大于准静态加载时。

2) 从Bernstein方程对土壤压力沉陷关系拟合效果来看，Bernstein方程对准静态加载条件下土壤压力沉陷关系拟合良好，拟合决定系数R2=0.923；但是对动态加载条件下土壤压力沉陷关系的拟合存在较大误差，拟合决定系数R2=0.852，拟合效果明显不如前者。

总结上述，可见加载速率会对土壤压力沉陷特性产生明显影响；以准静态加载为试验条件的传统土壤压力沉陷模型对符合准静态加载条件的土壤压力沉陷仿真结果具有良好的拟合效果，但是对动态加载时的土壤压力沉陷关系不能较好的拟合。

2.3 不同加载速率下土壤承压沉陷特性分析

车辆越野行驶过程中，随着车速的不断提高，车辆对地面施加载荷的速率可以达到几千千帕每秒。以美军M4坦克为例，其最大行驶速度42 km/h，总重32 t。文献[3]中对各国主战坦克平均最大接地压力Pmm进行了总结，其中M4坦克的平均最大接地压力Pmm=282.2 kPa。按照M4坦克履带接地长度为5 m，宽0.25 m，采用下式计算不同行驶速度对应的对地加载速率：

(2)

式中:P为加载速率(kPa/s)；v为车辆行驶速度(m/s)；l为履带接地长度(m)。

计算结果如表1所示：

表1 M4坦克行驶速度、加载速率对比

以表1中四种行驶速度所对应的平板加载速率为例，在2.1节建模的基础上，分别设置不同的加载速率，得到p1～p4四种加载速率下的土壤压力沉陷关系曲线，如图3所示。图3上还有准静态加载时的土壤压力沉陷曲线。

从图3可以看出，首先四种加载速率下得到的土壤压力沉陷曲线均与准静态加载时的土壤压力沉陷曲线存在明显差异，与1.4节得到的结果一致。其次，动态加载条件下，四种加载速率之间的土壤压力沉陷关系也存在不同。沉陷量相同时，土壤的承载力随着加载速率的提高而增大，在土壤进入塑性变形区开始阶段尤为明显。

3 土壤压力沉陷关系Prony指数函数拟合

上节通过对动态加载条件下土壤压力沉陷关系与准静态加载条件下土壤压力沉陷关系的对比，得出加载速率是影响土壤压力沉陷关系的重要因素，并且经典土壤压力沉陷拟合模型对动态加载条件下的土壤压力沉陷关系具有较大误差。因此，本节将从数值逼近的角度采用Prony指数型函数逼近方法对土壤压力沉陷关系进行拟合。

3.1 Prony指数型函数逼近方法

Prony法是以构造形如式2的指数型函数为目的的一种非线性数值逼近算法[15]。

(2)

其中：A，s为待定参数，e为自然底数，t为自变量，f为函数。

Prony法的主要原理为：令t=iT,使得式3成立。

f(iT)=fi，i=0,1,…,2n-1

(3)

其中：T为步长，{fi}为给定型值。

定义变量zj、αi，且分别满足如下两式：

zj=esjT,j=1,…，n

(4)

(5)

联立式(2)、(3)、(4)可得：

(6)

由式(5)、(6)可得：

k=0,1,…,n-1

(7)

由于αn=1，可将式(7)写为：

(8)

从而求出αi(i=0,1,…,n-1)。然后依据式(5)求解袋代数方程：

zn+αn-1zn-1+…+α1z+α0=0

(9)

得到n个根z1,z2,…,zn。最后按照式(10):

(10)

求解指数sj，再由式(6)求解系数Aj。

3.2 土壤压力沉陷关系拟合

当n=2时，Prony指数型函数等价于下式：

P=aebz+cedz

(11)

其中：P为土壤承受压力，z为沉陷量，a、b、c、d为常数。

采用式(11)对图3中不同加载速率下得到的土壤压力沉陷曲线进行拟合，拟合结果如表2所示。

表2 土壤压力沉陷曲线指数函数拟合参数

从表2可以看出，Prony指数型函数对不同加载速率下的土壤压力沉陷特性曲线均具有很好的拟合效果。其拟合精度R2最低为0.988。

以加载速率p1为例，分别采用Prony指数型函数和Bernstein方程对土壤压力沉陷关系进行拟合。拟合结果如图4所示。

由图4可以看出Prony指数函数对动态加载条件下的土壤压力沉陷关系拟合精度明显高于Bernstein方程。对其他加载速率下的土壤压力沉陷关系进行拟合，也能取得类似的结果。

如图5所示为准静态加载条件下上述两种拟合结果。

从图5中可得，Prony指数函数与Bernstein方程对准静态加载时的土壤压力沉陷关系均有良好的拟合效果，Prony指数函数拟合精度有进一步的提高。

综上所述，Prony指数函数对准静态和动态加载条件下的土壤压力沉陷关系均能取得较高的拟合精度。

3.3 Prony指数函数试验验证

Wills BDM[16]在湿黏土中采用尺寸分别为0.050 8×0.304 8、0.076 2×0.407 5、0.101 6×0.609 6 m的矩形平板进行土壤压力沉陷试验，并得到土壤压力沉陷试验数据。俄罗斯学者Modest Lyasko[5]引用其试验结果分别采用bekker模型、LSA模型、Kacigin-Guskeo[17]提出来的Kacigin-Guskeo函数进行拟合分析。在此基础上，采用Prony法对上述湿黏土压力沉陷试验数据拟合，得到压力沉陷曲线如图6所示。

如图6(a)所示，采用Prony指数函数对三种尺寸矩形载荷板的压力沉陷数据进行拟合，按照载荷板尺寸由小到大，其决定系数R2分别为0.996 6、0.999 7、0.999 4，试验验证表明土壤压力沉陷实测结果与采用Prony法模拟的结果一致。与图6(b)中各模型拟合结果相比，双指数模型具有更好的拟合效果。

注：图(a)中●为0.101 6×0.609 6 m试验值、◆为0.076 2×0.407 5 m试验值、■为0.050 8×0.304 8 m寸试验值

4 结论

1) 基于有限元仿真方法，分析了加载速率对土壤压力沉陷特性关系的影响。发现动态加载时的土壤压力沉陷特性与准静态加载时存在明显差异。不同加载速率对土壤的压力沉陷特性也存在影响，同一沉陷深度时土壤的承载力随加载速率的提高而增大。

2) 以Bernstein方程为例，分析了加载速率对经典模型对土壤压力沉陷特性拟合效果的影响。Bernstein方程对准静态加载条件下的土壤压力沉陷关系拟合效果较好，对动态加载条件下的土壤压力沉陷关系拟合误差较大。

3) 从数值逼近的角度对土壤压力沉陷关系进行拟合，发现Prony指数型函数对不同加载速率下的土壤压力沉陷关系均具有较好的拟合效果。