移动荷载作用下正交各向异性饱和地基动力响应的有限元分析

张春丽， 康彦冰， 杨鑫鑫， 祝彦知

(中原工学院 建筑工程学院， 河南 郑州 450007)

“十四五”期间是交通运输发展的黄金时期，《交通强国建设纲要》明确了我国交通发展的总体思路和“2035年基本建成交通强国、本世纪中叶全面建成交通强国”的发展目标[1]。随着我国高速公路、铁路等交通建设的飞速发展，车辆数量也在逐年增长，长期的移动荷载引起路基的沉降和破坏。我国沿海城市存在着大量的饱和地基土，饱和地基的力学性能差，容易引发许多交通工程问题。对于移动荷载作用下饱和地基的动力响应问题，许多学者进行了研究，其中多数研究是将地基土简化为各向同性或横观各向同性，如：孙波[2]基于Biot波动方程，获得了波数-频率域以及时间-空间域的解析解，研究了移动荷载作用下饱和地基中的波动特性；梁建文等[3]采用间接边界元法研究了层状饱和地基中地铁列车荷载诱发的地基振动问题；周捡平等[4]研究了地铁列车移动荷载作用下饱和软黏土地基的动力响应；钱建固等[5]研究了承受高速移动荷载的饱和多孔地基的动力响应；Dan等[6]基于Biot理论研究了移动交通荷载作用下饱和沥青路面的水动力耦合响应；陈元杰[7]基于Biot理论与半解析法研究了饱和土体在移动荷载作用下的动力响应。实际上，土体在天然沉积过程中有明显的各向异性，而目前关于移动荷载下正交各向异性饱和地基动力响应的研究很少。本文基于ABAQUS有限元软件和FORTRAN语言，建立了正交各向异性饱和地基模型，分析了移动荷载作用下土体表面竖向位移、土体竖向应力和孔隙水压力的变化规律。

1 有限元模型建立

1.1 地基模型建立

本文研究移动荷载下的平面应变问题，通过ABAQUS软件建立有限元地基模型，模型尺寸为100 m×50 m；对模型进行网格划分，近似全局尺寸为1.0 m。模型的网格划分情况及尺寸如图1和图2所示。

图1 模型网格划分示意图Fig. 1 Schematic diagram of model meshing

在地基表面施加沿x′轴正方向、速度为c的移动谐振荷载，荷载大小为q0eiwt，如图3所示。其中，w为荷载的角频率，w=2πf，f为荷载的振动频率；t为荷载移动时间；q0为荷载幅值；2b为荷载作用长度，2b=1.5 m。

图2 地基模型图(单位：m)Fig. 2 Model chart of subgrade (unit: m)

图3 计算模型荷载分布示意图Fig. 3 Load distribution of computational model

在图2和图3中，坐标系O′x′y′是固定的，根据文献[8]，引入移动坐标系Oxy，坐标系Oxy以速度c沿x′轴的正方向移动，即坐标变换为x=x′-ct,y=y′。

1.2 子程序开发

由于ABAQUS有限元无法直接模拟移动荷载，因此在有限元模型分析中采用FORTRAN语言编写DLOAD子程序，用来模拟移动荷载，控制荷载的大小、振动频率、作用范围等。

1.3 正交各向异性参数

正交异性材料的线弹性可以通过给出9个独立的弹性刚度参数来定义，应力-应变关系如下：

(1)

对于正交各向异性材料，弹性刚度参数为：

D11=E1(1-v23v32)Υ,

(2)

D22=E2(1-v13v31)Υ,

(3)

D33=E3(1-v12v21)Υ,

(4)

D12=E1(v21-v31v23)Υ=E2(v12-v32v13)Υ,

(5)

D13=E1(v31-v21v32)Υ=E3(v13-v12v23)Υ,

(6)

D23=E2(v32-v12v31)Υ=E3(v23-v21v13)Υ,

(7)

D44=G12

(8)

D55=G13

(9)

D66=G23

(10)

其中

(11)

式中：下标1、2、3分别代表x、y、z方向；σ代表正应力，τ代表切应力，ε代表正应变；γ代表切应变；E1、E2、E3分别代表x、y、z方向的弹性模量；Gij代表ij平面方向的剪切模量；vij代表泊松比，i表示作用力的方向，j表示由作用力引起伸缩的方向。本计算模型仅考虑平面应变问题，即ε33=0，γ13=0，γ23=0。

1.4 瑞利阻尼系数

利用ABAQUS软件进行分析时，需要输入瑞利阻尼系数α和β来设置阻尼，α与质量矩阵有关，β与刚度矩阵有关，阻尼矩阵可表示为质量矩阵和刚度矩阵的线性组合，如式(12)所示。其中，[C]为阻尼矩阵，[M]为质量矩阵，[K]为刚度矩阵。

[C]=α[M]+β[K]

(12)

α、β、阻尼比的关系如式(13)所示。

(13)

式中：ξ为第i个振型的阻尼比；ωi为第i个振型的角速度。分析中通常根据第一振型确定阻尼和角速度，因此，ξi=ξ1，ωi=ω1,则：

α=ξ1ω1

(14)

(15)

式中：ξ1为第一振型的阻尼比；ω1为第一振型的角速度。取ξ1=0.025，ω1=16π rad/s，代入式(14)、式(15)中，得出：α=1.257，β=4.974×10-3。

1.5 模型参数

具体模型参数见表1和表2。

表1 土体参数

表2 荷载参数

该地基模型的表面不排水，上边界的孔隙压力取为0。

通过比例系数km(m=1，2，3，4)来控制弹性模量和泊松比的大小，令E2=k1E1，E3=k2E1，v13=k3v12，v23=k4v12。当km=1时，E1=E2=E3，此时土体为各向同性；当km≠1时，土体为正交各向异性。在以下程序分析中，如没有特别说明，默认比例系数k1=1.2，k2=0.8，k3=1.2，k4=1.6，荷载频率f=8 Hz。

2 数值分析

2.1 模型可靠度分析

地基模型采用与文献[9]相同的参数，地基土简化为各向同性，得出地基模型与文献[9]的无量纲值对比图，如图4所示。图中，|u|表示土体表面竖向位移幅值，x表示土体表面各点至荷载中心点的距离。模拟结果显示，该模型与文献[9]的无量纲值基本一致，证明了该模型的准确性。

图4 荷载移动速度为0 m/s时地基模型与文献[9]的无量纲值对比图Fig. 4 Non-dimentional comparison chart between the foundation model and literature [9] when the load moving speed is 0 m/s

2.2 弹性模量对土体竖向位移的影响

当k2=0.8时，k1分别取1.2、1.6、2.0，得出弹性模量E2为不同值时的土体表面竖向位移，如图5所示；当k1=1.2时，k2分别取0.8、1.5、2.0，得出弹性模量E3为不同值时的土体表面竖向位移，如图6所示。其中，km=1.0的曲线为各向同性土体表面的竖向位移曲线。由图5和图6可以看出，在荷载作用范围3倍以内的区域，土体表面竖向位移曲线的变化趋势基本相同；在荷载作用范围3倍以外的区域，土体表面竖向位移几乎为零。由图5可以看出，随着k1的增大，土体表面的竖向位移最大值不断减小；各向同性时土体表面的竖向位移最大。由图6可以看出，随着k2的增大，土体表面的竖向位移最大值不断增大；各向同性时土体表面的竖向位移最大值仅大于k2=0.8时土体表面的竖向位移最大值。

图5 弹性模量E2不同时土体表面竖向位移图Fig. 5 Vertical displacement diagrams of soil surface under different elastic modulus E2

图6 弹性模量E3不同时土体表面竖向位移图Fig. 6 Vertical displacement diagrams of soil surface under different elastic modulus E3

2.3 荷载对土体竖向位移的影响

2.3.1 荷载振动频率对土体竖向位移的影响

在荷载移动速度35 m/s、50 m/s、100 m/s下，荷载振动频率分别取8 Hz、16 Hz、32 Hz，得出一定荷载移动速度下不同荷载振动频率时的土体表面竖向位移，如图7、图8、图9所示。由图7-图9可以看出，当荷载移动速度不变时，在荷载作用范围3倍以内的区域，随着荷载振动频率的增大，土体表面竖向位移的最大值越来越小，位移出现负值的位置离荷载中心越来越近；在荷载作用范围3倍以外的区域，土体表面竖向位移基本无变化；当荷载振动频率为32 Hz时，在荷载中心点附近出现多个峰值。

图7 荷载移动速度35 m/s下不同振动频率时土体表面竖向位移图Fig. 7 Vertical displacement diagram of soil surface under different vibration frequencies at load moving speed of 35 m/s

图8 荷载移动速度50 m/s下不同振动频率时土体表面竖向位移图Fig. 8 Vertical displacement diagram of soil surface under different vibration frequencies at load moving speed of 50 m/s

图9 荷载移动速度100m/s下不同振动频率时土体表面竖向位移图Fig. 9 Vertical displacement diagrams of soil surface under different vibration frequencies at load moving speed of 100 m/s

2.3.2 荷载移动速度对土体竖向位移的影响

在荷载振动频率8 Hz和16 Hz下，荷载移动速度分别取35 m/s、50 m/s、100 m/s，得出一定荷载振动频率下不同荷载移动速度时的土体表面的竖向位移，如图10和图11所示。由图10和图11可以看出，在荷载作用范围3倍以内的区域，当荷载移动速度增大时，土体表面竖向位移的最大值略微增大，并且出现变形滞后现象；在荷载作用范围3倍以外的区域，土体表面竖向位移基本无变化。

图10 荷载振动频率8 Hz下不同移动速度时土体表面竖向位移图Fig. 10 Vertical displacement diagram of soil surface at different moving speeds under load vibration frequency of 8 Hz

图11 荷载振动频率16 Hz下不同移动速度时土体表面竖向位移图Fig. 11 Vertical displacement diagram of soil surface at different moving speeds under load vibration frequency of 16 Hz

2.4 荷载对土体孔隙压力的影响

在荷载振动频率8 Hz下，不同荷载移动速度对饱和地基孔隙水压力的影响情况如图12所示。由图12可以看出，随着荷载移动速度的增大，土体的最大孔隙水压力明显减小。当荷载移动速度为35 m/s时，不同荷载振动频率对饱和地基孔隙水压力的影响情况如图13所示。由图13可以看出，随着荷载振动频率的增大，土体的最大孔隙水压力略微减小。对于不同的荷载移动速度和荷载振动频率，在荷载作用中心点深度为0.5 m的位置，土体的孔隙水压力达到最大值。

图12 不同荷载移动速度下饱和地基孔隙水压力沿深度方向变化图Fig. 12 Variation diagram of pore water pressure of saturated foundation along depth under different load moving speeds

2.5 土体深度对地基应力的影响

当荷载移动速度为0 m/s和50 m/s时，土体表面、土体表面以下1 m及2 m处的地基应力如图14和图15所示，图中z表示土体深度。由图14和图15可以看出，随着土体深度的增大，土体的竖向正应力最大值有明显下降，并且竖向正应力向荷载作用点的左右两侧扩散，在荷载作用范围2倍以外区域的应力几乎为0。当荷载移动速度为0 m/s时，土体表面的竖向正应力曲线和边界条件相吻合；当荷载移动速度为50 m/s时，荷载中心点附近土体的竖向正应力出现多个峰值，且在荷载作用点的左右两侧出现拉应力，右侧的拉应力峰值比左侧大。随着荷载移动速度的增大，相同深度的土体的竖向正应力曲线震荡幅度明显增大，竖向正应力最大值明显右移，导致拉应力的点离荷载中心点越来越近。

图14 荷载移动速度为0 m/s时不同深度的地基应力图Fig. 14 Stress diagram of foundation at different depths when the load moving speed is 0 m/s

图15 荷载移动速度为50 m/s时不同深度的地基应力图Fig. 15 Stress diagram of foundation at different depths when the load moving speed is 50 m/s

3 结论

本文基于ABAQUS有限元软件和FORTRAN语言，建立了正交各向异性饱和地基模型，研究了移动荷载作用下正交各向异性饱和地基的位移、应力及孔隙水压力的变化规律，得出以下结论：

(1) 土体的弹性模量E2、E3对土体表面的竖向位移有较大影响。当弹性模量E2增加时，土体表面竖向位移的最大值减小；当弹性模量E3增加时，土体表面竖向位移的最大值增大。通过对比发现，考虑饱和土体的正交各向异性能更准确反映土体的真实情况。

(2) 荷载振动频率对土体表面的竖向位移、孔隙水压力均有影响。随着荷载振动频率的增大，土体表面竖向位移的最大值越来越小，沿深度方向的最大孔隙水压力略微减小。

(3) 荷载移动速度对土体的竖向位移、孔隙水压力、竖向正应力均有影响；随着荷载移动速度的增大，土体表面竖向位移的最大值略微增大，沿深度方向的最大孔隙水压力明显减小，相同深度的竖向正应力曲线震荡幅度明显增加。

(4) 随着土体深度的增大，土体的竖向正应力分布有很大变化，竖向正应力的最大值明显下降。