不同条件下土壤压力沉陷关系有限元分析

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

(陆军工程大学石家庄校区，河北 石家庄 050003)

车辆行驶过程中，对地面产生垂直作用力，使土壤产生沉陷作用，影响车辆的行驶性能。因此研究土壤在垂直载荷的作用下，沉陷量与载荷的关系至关重要。土壤的压力沉陷特性最早应用于拖拉机在田间行驶性能的研究。经典的土壤所受压力和沉陷量之间的关系以指数函数形式为代表。贝克在贝恩斯坦提出的压力沉陷公式基础上，结合土木工程领域的地基承载力理论，提出了考虑土壤参数和载荷板宽度的土壤压力沉陷特性方程，并且得到了广泛的应用[1-3]。经过几十年的发展，土壤的压力沉陷特性模型陆续出现了英国科学家利斯[4]的幂函数模型、俄国科学家库兹可夫的双曲正切模型和日本科学家的双曲线模型[5-6]。上述经典土壤压力沉陷理论模型对研究车辆越野性能具有重要的指导意义，但是后人在此基础上的研究以静态和准静态加载条件为研究基础[7-9]，针对载荷加载速率对土壤压力沉陷特性的影响考虑较少。

近年来，有限元理论在土壤力学分析中的应用日趋广泛，土木工程领域对有限元理论的应用相对更加系统和完善[10～12]。地基极限承载力的有限元分析是土木工程中非常重要的一方面，结合太沙基的极限承载力理论并借鉴土木工程中运用有限元理论分析地基极限承载力的技巧，研究车辆地面力学领域土壤压力沉陷特性是可行的[13]。笔者基于有限元理论研究了不同载荷板尺寸和加载速率对土壤压力沉陷关系的影响。

1 有限元分析中的土壤本构模型

ABAQUS是一款大型有限元分析软件，特别是在土力学的分析中具有独特的优势[14]。ABAQUS软件内置了多种土壤本构模型，包括以广义胡克定律为基础的各向同性弹性模型、正交各向异性弹性模型和各向异性弹性模型，以及以Mohr-Coulomb模型、Drucker-Prager模型和临界状态塑性模型为代表的塑性模型[15]。

1.1 土壤的弹性模型

材料的弹性本构模型描述材料在弹性变形范围内应力应变的关系，包括两个最基本的参数杨氏模量E和泊松比μ。土壤常用的弹性本构模型包括线弹性模型、多孔介质弹性模型和线粘弹性模型。笔者采用适用最广泛的各向同性弹性模型描述土壤的弹性行为。

1.2 土壤的塑性模型

土壤的塑性模型描述了土壤的屈服、硬化、流动准则、剪胀等塑性行为[16]。经典的Mohr-Coulomb模型和Drucker-Prager模型无法解决由于等向压缩导致的土壤永远不会屈服的现象。为此，笔者将采用修正Drucker-Prager盖帽模型来描述土壤的塑性行为。修正Drucker-Prager盖帽模型在线性Drucker-Prager模型的基础上引入材料的压缩屈服，并控制剪胀的扩展。

修正Drucker-Prager盖帽模型的屈服面如图1所示。屈服面由3部分组成，包括剪切破坏面、盖帽曲面和过渡面。

盖帽屈服面方程为

R(d+patanβ)=0,

(1)

式中：Fc为屈服面;p为主应力;pa为屈服面与过渡面交点值;R为控制盖帽几何形状的参数;t为偏应力;α为控制过渡曲面形状的参数;β为p-t平面上的摩擦角;d为p-t上的粘聚力。

修正Drucker-Prager盖帽模型的剪切破坏准则使用Mohr-Coulomb准则，Mohr-Coulomb准则被认为最符合土壤的剪切破坏条件[17]，其在p-t平面上表示为一条直线，如图1所示。修正Drucker-Prager盖帽模型的剪切破坏面方程为

Fs=t-ptanβ-d=0.

(2)

此外，修正Drucker-Prager盖帽模型采用非相关联的流动准则。

1.3 土壤本构参数的确定

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

文献[18]采用上述两种土工试验对起伏路面松软土壤进行参数测定，在含水率5%的情况下，测得土壤内摩擦角φ=27.33°，粘聚力c=6.38 kPa.通过三轴压缩试验测得土壤弹性模量E=20.2 MPa，泊松比μ=0.32，总密度为1 813 kg/m3.利用三维问题中Mohr-Coulomb模型与Drucker-Prager模型参数之间的转换关系可以将上述测得的内摩擦角φ和粘聚力c转换为Drucker-Prager模型中的内摩擦角β和粘聚力d，转换关系如下：

(3)

(4)

(5)

式中：k为流应力比;σc为单轴抗压强度。

利用上述Mohr-Coulomb模型与Drucker-Prager模型参数之间的转换公式，得到修正Drucker- Prager盖帽模型中的参数分别为：弹性模量E=20.2 MPa,泊松比μ=0.32,内摩擦角β=47.31°,粘聚力d=40.15 kPa.

2 土壤平板沉陷试验有限元建模

2.1 有限元模型前处理

三维建模中，采用三维实体可变性模型，建立4个半径r分别为0.05、0.10、0.15、0.20 m的载荷板。考虑到土体模型边界对仿真结果的影响和模型计算效率，土体模型采用2 m×2 m×1.5 m的长方体。建立材料属性分别为钢和土壤的截面属性，钢的材料参数为：密度7 800 kg/m3、弹性模量E=207 GPa、泊松比μ=0.3；土壤本构模型采用各向同性线弹性模型和修正Drucker-Prager盖帽塑性模型，参数如1.3所述。

装配体及边界条件为：载荷板装配与土体上表面正中央，采用刚体约束并设置参考点；约束参考点水平方向位移自由度和3个方向转动自由度，约束土体四周表面的水平位移自由度和底部3个方向位移自由度。土体和载荷板网格划分均采用八节点线性六面体单元，设置减缩积分和沙漏控制。土体网格拓扑结构采用中间致密四周相对稀疏的方法。采用显式动态分析。每种尺寸的载荷板分别以1 592.4、796.2、318.5、159.2、106.2 kPa/s的加载速率沉陷7～15 cm。装配体及网格划分如图2所示。

2.2 仿真结果分析

取半径为0.10 m的载荷板在318.5 kPa/s的加载速率下得到的位移增量云图和压力-沉陷量关系曲线，分别如图3、4所示。

从图3可以看出，土体在压力的作用下产生下陷，并且具有明显的应力核。从图4可以看出，土壤在竖直载荷作用下产生的变形主要分为3个阶段。第1阶段，土壤随着载荷从0逐渐增大，产生弹性变形，载荷与沉陷量为线性关系，这一阶段土壤所受压力与沉陷量的关系如图中直线部分所示；第2阶段，随着载荷的进一步增大，土壤进入从弹性变形到塑性变形的过渡阶段，这一阶段在图3中表现不明显，土壤弹塑性变形区分较容易，极限承载力位置明显，说明土壤发生整体剪切破坏；第3阶段，土壤达到极限承载力，进入完全塑性变形阶段，土壤在较小载荷增量的作用下发生失稳变形是这一阶段的主要特征。曲线整体特征符合车辆实际行驶过程中压力沉陷关系[19]。

2.3 考虑加载速率的土壤压力沉陷关系

车辆在实际行驶中，车轮或履带板对地面的作用是动态的，随着车速的提高，载荷作用的时间越短，动态效应越明显。特别是高速行驶的越野车辆，其行驶机构对地面的加载过程与传统的土壤压力沉陷模型适用的加载条件存在明显差别。笔者采用4种不同尺寸的平板分别以不同加载速率进行有限元仿真分析，研究加载速率对土壤压力沉陷关系的影响。以r=0.10 m载荷板为例，加载速率对压力沉陷关系曲线的影响如图5所示。

从图5中可以看出，在r=0.1 m的载荷板作用下，土壤的极限承载力随加载速率的增大而增大，并且随着加载速率的增大，极限承载力增大效应越明显。4种尺寸载荷板压力沉陷曲线随加载速率的变化趋势基本相同。

2.4 考虑平板尺寸的土壤压力沉陷关系

载荷板尺寸是影响土壤压力沉陷关系的一个重要因素。Bekker公式在贝恩斯坦公式的基础上，考虑土壤粘聚系数、土壤摩擦系数和载荷板尺寸对土壤压力沉陷关系的影响，认为系数k与平板尺寸的倒数成线性关系。实际试验过程中，有学者发现两者之间并非线性关系[20]。土壤压力沉陷关系与载荷板几何形状、尺寸、加载条件等因素有关。

笔者以圆形载荷板为例，研究不同载荷板半径对土壤压力沉陷关系的影响。分别设置载荷板半径分别为0.05、0.10、0.15、0.20 m，以318.5 kPa/s的加载速率加载，得到压力沉陷关系曲线如图6所示。

从图6可以看出，在相同加载速率的情况下，同一沉陷量处的载荷随平板半径的增大而减小，极限承载力随着载荷板半径的减小而增大。这说明载荷板尺寸对土壤压力沉陷关系会产生一定影响，并且二者之间存在的负相关关系与Bekker公式中载荷板宽度b对压力p的影响是一致的。

3 加载条件对模型参数的影响

传统土壤压力沉陷模型参数被认为是与载荷板形状、尺寸、加载条件无关的常量，例如Bekker公式中的kc、kφ、n；Bernstein方程中的k、n等模型参数。但是，俄罗斯学者[19]利用实测数据对各种土壤压力沉陷解析模型[21-22]参数进行了解算，并与实测数据进行对比，发现模型参数不同的试验条件和载荷板形状尺寸得出的模型参数并非常数。

笔者以Bernstein方程为例，分析加载速率和载荷板尺寸对模型参数的影响。

3.1 加载速率对模型参数的影响

Bernstein方程为

p=kzn,

(6)

式中：p为土壤单位面积所受压力;z为土壤沉陷量;n为土壤变形指数;k为土壤体积压缩系数。

对Bernstein方程两边对数化得到：

lnp=nlnz+lnk.

(7)

理想条件下，式(7)在直角坐标系中为一条直线。取加载速率为318.5 kPa/s，载荷板半径为0.10 m时的压力和沉陷量数据，得到其在对数坐标系中的关系如图7所示。

从图7可以看出，仿真得出的土壤压力沉陷关系在对数坐标系中并非理想的线性关系，这一结果与文献[19]的研究结论是一致的。其余4种加载速率下，各尺寸载荷板的压力沉陷关系在对数坐标系中也存在类似特征。采用最小二乘法对lnz-lnp关系进行线性拟合，在加载速率分别为1 592.4、796.2、318.5、159.2、106.2 kPa/s的条件下，解算出不同半径载荷板采用p=kzn拟合时的参数n和k的值，如图8所示。

从图8(a)中可以看出，同一载荷板在不同加载速率下的n值并非常数，而是随着加载速率的不同成一定趋势的变化。载荷板半径为分别0.05、0.10、0.20 m时，n值随着加载速率的增大总体成减小趋势，并且减小速率随着加载速率的增大逐渐趋于平缓。载荷板半径为0.15 m时，n值随加载速率的变化趋势与其他3条曲线的变化趋势相反。对比图8中(a)、(b)两图，可以看出加载速率对参数k的影响规律与对参数n的影响规律一致。

3.2 载荷板尺寸对模型参数的影响

设置载荷板半径分别为0.05、0.10、0.15、0.20 m，并设置相同的加载速率进行仿真分析。利用最小二乘法求解出一定加载速率下不同载荷板半径对应的参数k、n的值。得到载荷板尺寸与参数k、n的关系如图9所示。

从图9(a)中可以看出，参数n的值随着载荷板半径的增大呈现整体增大的趋势，5种加载速率下得到的关系曲线均是如此。不同加载速率之间参数n与载荷板半径的关系出现差异。加载速率分别为796.2、318.5、159.2 kPa/s时，参数n与载荷板半径的关系近似为线性关系；加载速率分别为1 592.4、106.2 kPa/s时参数n随载荷板半径的增大跳动性较明显，但整体呈现越来越大的趋势。

从图9(b)中可以看出，不同加载速率下参数k与载荷板半径的关系差异性比图9(a)更明显。加载速率分别为796.2、318.5、159.2 kPa/s时，参数k随载荷板半径增大的趋势较图9(a)缓慢一些；加载速率分别为1 592.4、106.2 kPa/s时，参数k随载荷板半径的变化跳动幅度更大。比较图9中(a)、(b)两图曲线的变化趋势，发现载荷板半径对参数n、k的影响具有相似性，不同点是载荷板半径对参数k的影响更敏感。

结合3.1节对加载速率与模型p=kzn参数影响的分析可以发现，土壤压力沉陷关系模型p=kzn参数并非是与加载速率和载荷板尺寸无关的常数。对比图8和图9可以看出加载速率和载荷板尺寸对模型参数的影响规律也不相同，并且加载速率和载荷板尺寸对模型参数的影响不是相互独立的，考虑单一条件对模型参数的影响时，另外一个参数的变化也会影响模型参数的变化趋势。

4 结论

笔者利用有限元仿真，分析了加载速率和载荷板尺寸对土壤压力沉陷关系的影响，发现加载速率和载荷板尺寸的不同会影响土壤压力沉陷关系。采用p=kzn方程对仿真数据进行拟合时，不同加载速率和载荷板尺寸会得出不同的模型参数。具体结论如下：

1)不同试验条件会对土壤压力沉陷关系产生影响，载荷板尺寸相同时，加载速率越大，同一沉陷深度的压力越大，随着载荷板尺寸的增大，加载速率对压力沉陷关系的影响有减小的趋势。加载速率相同时，随着载荷板尺寸的增加，相同沉陷深度的压力减小。

2)不同试验条件会对土壤压力沉陷模型参数产生影响。p=kzn方程参数在不同加载速率和载荷板尺寸的条件下并非常数，参数n、k受单一条件影响的规律较为类似，不同条件对参数n、k的影响规律差别较大，参数与加载条件之间的关系需要进一步研究。