条形荷载下梯度非均匀非饱和土的动力响应分析1)

施力维 马 强 ,†, 舒进辉

* (青海大学土木工程学院,西宁 810016)

† (青海省建筑节能材料与工程安全重点实验室,西宁 810016)

引言

动荷载作用下地基动力响应问题的研究一直是土木工程界一个重要的研究课题.诸多学者对动荷载作用下均质单相弹性地基[1-2]、两相饱和土地基[3-4]和三相非饱和土地基[5-10]的动力响应问题进行了一系列研究,有关动荷载作用下均质土地基动力响应的研究已有很多重要的成果.

然而需要注意的是,上述对土体动力响应问题的研究绝大部分工作是基于土体具有统计各向同性且均匀等假设为前提条件的.而在自然环境下的土体绝大多数都为固-液-气三相组成的非饱和多孔介质,且由于沉积年代的不同表现出不同的物理力学性质,土体沿深度方向具有非均匀性.因此采用层状非饱和土土地基模型更加符合实际工程情况.目前关于单相层状非均匀介质的计算模型较多,主要有薄层单元法[11-12]、传递矩阵法[13-15]、刚度矩阵法[16-18]、传递和反射矩阵方法[19-21]等.除了上述方法,由于回传射线矩阵法(reverberation-ray matrix method,RRMM)推导的矩阵方程中不包含正指数项,因此较好地防止了在高速、高频动载荷作用下因土层厚度或参数相差较大而出现的数值计算问题,被广泛运用于层状模型计算中[22].文献[23]扩展和改进了回传射线矩阵的方法,用于分析各向异性层状弹性介质中自由波的传播,证明了回传射线矩阵法的高精度.柳伟等[24]将回传射线矩阵法推广至桩土系统的振动分析中,分析了外露长度、埋置深度、桩端约束情况对埋置结构自振特性的影响.文献[25]将回传射线矩阵法推广到地基梁自振特性的研究中,对不同边界条件下的自振频率、衰减系数及模态函数进行了分析.随后文献[26]建立了在动载荷下梯度非均匀饱和土体的回传射线矩阵法计算列式,提出回传射线矩阵法对非均匀两相饱和土动力响应分析也具有很好的适用性.Ma 等[27]基于回传射线矩阵法研究了饱和多孔功能梯度材料对饱和土地基中移动荷载引起的振动的隔振效果.

综上可知目前关于层状弹性地基和饱和土地基波动响应问题的研究较多,有关层状非饱和土地基动力响应问题的研究还鲜有报道.因此,鉴于非饱和土是土体在自然界更为普遍的一种存在状态,考虑土体的非均匀性,本文基于多孔介质混合物理论,采用回传射线矩阵法研究了条形荷载作用下梯度非均匀非饱和土地基的动力响应问题.通过傅里叶积分变换和Helmholtz 矢量分解原理,结合边界条件,建立了在条形荷载作用下非均匀非饱和土体的计算列式.分析了非饱和土物理力学性质在深度方向呈梯度变化时应力、位移以及孔隙压力等物理量的变化规律.

1 非饱和土控制方程及通解

1.1 非饱和土控制方程

基于非饱和土多孔介质混合物理论,三相非饱和土体动力控制方程如下[8,10]

式中,pw和pa分别为孔隙水压力和孔隙气压力;ρ=(1-n)ρs+nSrρw+n(1-Sr)ρa,其中 ρ 代表非饱和土介质的总密度,ρm(m=s,w,a)分别代表固、液、气三相的密度,Sr代表饱和度,n代表孔隙率;顶标·表示对时间t求导,wi和vi分别代表孔隙水和孔隙气沿i方向相对于土骨架ui的位移.χ代表有效应力系数,受含水量等因素影响,在本文中根据文献[28]将χ的值等同于Sr的值;λ和 μ代表Lame 常数;a=1-Kb/Ks,Kb和Ks分别代表土骨架和土颗粒的压缩模量,其中Kb=λ+2μ/3 ;Kw和Ka分别代表孔隙水和孔隙气的体积压缩模量;e=∇·u代表土骨架的体积应变,u代表土骨架位移张量,w和v代表孔隙水和孔隙气的位移张量;系数A11～A24表达式详见文献[10].g代表重力加速度,kw和ka分别代表孔隙水和孔隙气的渗透系数.根据Fredlund[29]理论,水和空气的渗透系数可表示为

式中,ηw,ηa分别代表孔隙水和孔隙气的黏性系数;krw,kra分别代表孔隙水和孔隙气的相对渗透系数;κ代表土的固有渗透率.相对渗透系数可由土-水特征曲线导出,本文采用V-G 模型[30]可得

式中,Se代表有效饱和度;α,m,d代表拟合参数,m=1-1/d;Se=(Sr-Sw0)/(1-Sw0),Sw0代表水的束缚饱和度;pc=pa-pw为基质吸力.

1.2 方程的求解

在简谐振动中,所有变量可表示为f=的形式,其中表示变量f的幅值(简化起见,下文公式推导略去～号),ω 为圆频率.对空间变量x进行傅里叶变换及其逆变换

式中,ξ 为x方向上的波数;f是空间域变量,是相应频率域变量;i为虚数,t代表时间.

式(1a)和(1b)可改写为

对式(5)两边取散度后代入式(1c)和(1d)中可得

式中b11=baρaA11,b12=baρaA12,b13=baρaA13+ω2ρa,b21=bwρwA21,b22=bwρwA22,b23=bwρwA23+ω2ρw.

将式(5)代入到式(1e)中可得

式中,b1=-aχ-ω2/bw,b2=-a(1-χ)-ω2/ba,b3=ρω2+ω4ρw/bw+ω4ρa/ba.

同样,对式(5)两边取散度后代入到式(7)并联立式(6)可得

对方程式(6)和式(8)中的空间变量z进行傅里叶变换

求解常微分方程组(9)可得通解如下

将式(7)进行傅里叶变换可重写为

结合式(10)和式(11)可得两个二阶线性非齐次微分方程,求解可得其通解为

对土骨架体积应变e=∇·u进行傅里叶变换并联立式(10c)和式(12)可得

将式 εij=代入到非饱和土应力-应变关系式 σij=λeδij+2μεij-δijap中,式中P为等效孔隙流体压力,表达为p=χpw+(1-χ)pa.通过对应力-应变关系式和式(5)进行傅里叶变换后可得

结合方程式(10)和式(12),最后可得应力、孔隙水和孔隙气的相对位移在傅里叶变换域的通解为

方程式(10a)、式(10b)、式(12)和式(15)为非饱和土体动力响应问题的应力、位移和孔压等物理量在频率域的解答.

2 回传射线矩阵法求解

本文考虑如图1 所示的条形简谐荷载作用下梯度非均匀非饱和土地基的动力响应模型.非饱和土厚度为H,表面受到条形简谐荷载,其中荷载幅值为q0,分布长度为 2l.

图1 条形荷载作用下梯度非均匀非饱和土地基示意图Fig.1 Physical model of graded non-homogeneous unsaturated soil under strip load

对于非均匀非饱和土,由于土体的非均匀性导致土体动力问题的控制方程为变系数偏微分方程,一般情况下很难获得解析解.由于回传射线矩阵法(RRMM)对非均匀单相弹性介质和非均匀饱和多孔介质都具有很好的适用性,因此本文尝试采用RRMM 法对非均匀非饱和多孔介质动力控制方程进行求解.为此,通常将非均匀材料在材料物理力学性质变化方向上被简化为多个均匀材料薄层,层数N的数量取决问题所需精度.回传射线矩阵法的关键在于建立一对对偶坐标系,图2 给出了第i层的对偶坐标系,(xi(i+1),zi(i+1))和 (xi(i-1),zi(i-1))中上标i表示局部坐标系位于第i层,上标 (i+1)和 (i-1)分别表示指向第i+1和i-1 层,通过建立局部坐标系,可以有效避免计算过程中的数值问题.

图2 局部坐标示意图Fig.2 Dual local coordinates at the interface

2.1 散射关系

在局部坐标系 (xi(i-1),0)下,式(10a)、式(10b)、式(12)以及式(15)可改写为

式中X,Mn,Nn,Ln,Hn的表达式见附录A.

节点i处的应力、位移以及孔压在对偶坐标系下的连续性条件为

利用式(16)和式(17)整理后得

式中入射波和出射波的波幅向量分别为节点i处的局部散射矩阵为Si,散射矩阵中各元素详见附录B.

考虑非饱和土地基表面受条形均布荷载作用,底面固定的边界条件为

在z=0 处

在z=H处

利用式(16)、式(19)和式(20),整理后可得

式中

矩阵S0和SN中各元素详见附录C.

合并方程式(18)、式(21)和式(22)可得

方程(23)中提供的 8N个方程不满足向量A和B中 16N个未知量解的要求,因此须再寻找一组方程.

2.2 相位关系和回传射线矩阵

将式(16d)代入式(24)中可得

引入新矢量

结合式(23)、式(25)和式(26)可得

式中P为整体相位矩阵,定义为

U为整体置换矩阵,定义为

将式(27)代入式(23),最终可得

式中R=SPU称为回传射线矩阵,I为单位矩阵.

通过式(30)求出C和D后,代入式(16),就可以获得非饱和土地基中任一点处的应力、位移和孔压等物理量,再通过傅里叶逆变换即可求得相应的各物理量在空间域上的表示.

3 数值算例

3.1 验证

为了验证本文算法的正确性,将本文非均匀非饱和土地基退化为均匀非饱和土地基,并与Shi[10]求得的条形荷载下非饱和土地基的解答进行比较.验证时取计算参数为[10]:H=20 m ,μs=1.91×107Pa,ρs=2700 kg/m3,Sr=0.8,υ=0.2,n=0.6,Ks=3.6×1010Pa,k=1×10-12m2,q0=1 kPa,ω=1 rad/s,γ=0 .图3 给出了本文数值解与文献解的计算结果,从图中可以看出,本文计算结果与文献结果[10]比较接近,从而验证了本文算法的有效性.

图3 本文数值解与文献解的比较Fig.3 Comparison between present work and analytical

3.2 数值结果和讨论

考虑沿深度方向非均匀非饱和土层的物理力学参数按幂函数变化,则有

式中,G(x3)为在x3处的 λ,μ,n,ρ,k,υ,Sr等物理力学参数;γ 为梯度因子,γ 值的不同代表非饱和土体的非均匀程度不同.可以看出,当 γ=0 时,非饱和土物理力学参数为土层表面处的材料参数;当γ→∞时,非饱和土物理力学参数趋向于土层底面处的材料参数.

作为本文解答的具体应用,取土层厚度H=40 m,层数N=40 ,均布荷载幅值q0=1 kPa,频率ω=1 rad/s,l=1 m .非饱和土地基的物理力学参数见表1所示[31].简化起见,假设沿深度方向各物理力学参数按同样的梯度因子变化.

表1 非均匀地基的物理力学参数Table 1 Physico-mechanical properties of non-homogeneous foundation

为了研究梯度因子对非均匀非饱和土地基动力响应的影响规律,图4 给出了梯度因子 γ 对竖向位移沿深度z的变化曲线.从图中可以看出,随着梯度因子 γ 的增大,竖向位移幅值刚开始随之显著减小,但随着 γ 的继续增大竖向位移降低幅度开始变小.这是由于梯度因子 γ 的增大,使得抵抗变形能力强的基底材料增多,随着梯度因子 γ 的继续增大导致土的整体刚度较大,因此竖向位移幅值增量较小.

图4 梯度因子 γ 对竖向位移沿深度z 变化的影响曲线Fig.4 Influence of the gradient factor γ on vertical displacement variations against with depth z

图5 和图6 给出了梯度因子 γ 对孔隙水压和孔隙气压沿深度z的变化曲线.从图5 中可以看出,随着梯度因子 γ 的增大,孔隙水压的幅值先增大后减小,并且土体非均匀程度越高,孔隙水压的幅值越大,当非饱和土体趋于均质土时(γ=0,50),孔隙水压幅值大小相接近.从图6 中可以发现,随着梯度因子γ的增大,孔隙气压幅值不断增大,并且其沿深度方向的振动频率增大,波峰值不断靠近地表处附近.

图5 梯度因子 γ 对孔隙水压沿深度z 变化的影响曲线Fig.5 Influence of the gradient factor γ on pore water pressure variations against with depth z

图6 梯度因子 γ 对孔隙气压沿深度z 变化的影响曲线Fig.6 Influence of the gradient factor γ on pore air pressure variations against with depth z

图7 给出了梯度因子 γ 对正应力沿深度z的变化曲线.图中可以看出,正应力幅值随着梯度因子γ的增大先增大后减小,并且土体非均匀程度越高,正应力的幅值越大,当非饱和土体趋于均质土时(γ=0,50),正应力的幅值大小相接近,这与孔隙水压幅值大小受土体非均匀程度的影响相同.综合图4～图7还可以发现梯度因子 γ 显著地改变了竖向位移、正应力、孔隙水压和孔隙气压沿深度z方向上的振动模态.

图7 梯度因子 γ 对正应力沿深度z 变化的影响曲线Fig.7 Influence of the gradient factor γ on normal stress variations against with horizontal depth z

4 结论

本文基于多孔介质混合物理论,研究了条形荷载作用下梯度非均匀非饱和土地基的动力响应问题,建立了非均匀非饱和土地基在条形荷载作用下回传射线矩阵法的计算列式,分析了非饱和土体物理力学性质沿深度方向按梯度变化时应力、位移以及孔隙压力等物理量的变化规律.结果表明:(1)土体非均匀性对非饱和土的动力响应具有显著影响,梯度因子改变了竖向位移、正应力和孔隙压力在深度方向上的振动模态,其中孔隙气压在深度方向的振动频率随着梯度因子的增加而不断增大,峰值不断靠近地表处附近；(2)竖向位移随梯度因子的增大不断减小,正应力和孔隙压力随着梯度因子的增大先增大后减小,并且土体非均匀程度越高,其相应幅值越大;而孔隙气压随梯度因子的增大不断增大.

附录A

附录B

附录C

附录D