基于非局部Biot 理论的循环荷载下饱和土地基动力特性研究

丁海滨，管凌霄，童立红，徐长节,3，颜建伟

(1.华东交通大学轨道交通基础设施性能监测与保障国家重点实验室，江西，南昌 330013；2.江西省地下空间技术开发工程研究中心，江西，南昌 330013；3.华东交通大学江西省岩土工程基础设施安全与控制重点实验室，江西，南昌 330013)

探究循环动荷载下饱和土地基的动力响应及固结沉降对实际工程具有重要意义，如厂房机械振动对地基的影响，列车荷载下饱和地基动力特性及软土地区地基加固[1]等。饱和孔隙介质的动力控制方程最早是由BIOT[2- 3]提出，该方程能够较准确地反映饱和土的动力特性。之后，许多学者[4-5]通过室内动力实验验证了Biot 理论的正确性。Biot理论也因其形式简单且其模型参数易通过试验获得，而在实际工程得到了广泛的应用。Biot 动力控制方程的原始形式中的变量为u(土骨架位移)和w(流体相对土骨架位移)，即u-w格式控制方程。后续学者们根据所研究问题的不同，将原始的控制方程变换为u-p和u-w-p等格式的控制方程[6-7]。

由于基于Biot 理论的二维或者三维固结问题的解析解难以获得，目前学者们所求得的解析解都是对所研究的问题进行简化后获得的，如：简化为轴对称问题[8-9]，或者将该问题简化为一维问题，求得其解析解[10-12]。简化所得到的解析解与实际情况存在一定的差异，为此，学者大多采用数值方法求解Biot 理论控制方程。PREVOST[13]采用有限元法和隐式积分算法，研究了波在饱和土两相介质中的传播问题。CARTER 等[14]采用有限元数值方法求解了弹塑性地基的固结问题。SIMON等[15]采用有限元法研究了u-w，u-w-p及u-p三种形式的控制方程在饱和孔隙介质的复杂边值中的应用。FERRONATO等[16]提出了三维混合有限元方程，求解了Biot 固结方程，以减轻数值计算引起的孔隙水压不稳定问题。MENENDEZ等[17]采用有限元法求解了不可压缩流体和渗透系数变化的固结模型，分析了饱和土介质的非线性固结问题。基于Arbitrary Lagrangian-Eulerian 法和摩尔库伦(MC)及修正剑桥本构模型(MCC)，SABETAMAL等[18]采用Generalized-α 积分法对u-p耦合控制方程时间进行离散，求解了荷载作用下饱和孔隙介质的大变形问题。基于Biot 固结理论，WU 等[19-20]采用数值流变模型建立了饱和土动力固结的有限元方程，编制了有限元计算程序，分析了外荷载作用下低渗透性饱和黏土边坡的动力固结问题。NAVAS等[21-23]基于u-w格式的控制方程，采用局部最大熵原理(LME)构建无网格法形函数，研究了饱和土地基的固结问题。NAZEM等[24]提出最大熵无网格法，用于研究饱和土介质的固结问题，结果显示所提出的模型可给出固结沉降的稳定的解。TORABI等[25]提出了一种求解固结方程的稳定时间积分算法。

然而，以上针对饱和孔隙介质动力固结的研究都是基于经典Biot 理论，虽然Biot 理论在工程中得到了广泛的应用，其假设饱和土中的波长远大于孔隙尺寸，但在高频情况下该假设已不再成立[26]。为此，TONG 等[27]通过引入了非局部参数考虑了孔隙尺寸效应的影响，提出了非局部Biot理论。随后，利用非局部Biot 理论，TONG等[28]研究了饱和土介质中Rayleigh 波传播特性；XU等[29]，DING 等[30]及徐长节等[31]研究饱和土介质中圆形衬砌的动力响应问题；TONG 等[32- 33]研究移动荷载下饱和土地基的动力响应问题。

由以上分析可知，高频下，孔隙尺寸效应对结构动力特性的影响显著，为此，本文基于TONG等[27]所构建的u-w形式的非局部Biot 理论，采用有限单元法，结合虚位移原理及Newmark 积分法，得到了非局部Biot 理论的有限元离散方程，编制对应的有限元计算程序，分析孔隙尺寸效应对循环动荷载作用下饱和土地基动力特性的影响。本文得到了非局部Biot 理论的控制方程数值解，使其能用于解决复杂的工程问题，以期为实际工程中循环动荷载下饱和土地基动力响应及固结沉降问题提供参考。

1 非局部Biot 控制方程的有限元空间离散

1.1 非局部Biot 理论控制方程

采用有限元求解一个系统，首先需建立并求解该系统的线性方程组，称之为有限元方程。有多种方法可用于求解有限元方程，如加权参数法、虚位移法及变分原理等。本章采用虚位移法获得系统的有限元方程。非局部Biot 理论控制方程式为[27]：

式中：u、w分别为土骨架位移和流体相对于土骨架位移； α =1-Kb/Ks；M为Biot 参数； ∇2为Laplace算子； σ′为总Cauchy 应力； ρf为孔隙流体的体积密度， ρ =(1-n0)ρs+n0ρf，n0为初始孔隙比， ρs为土颗粒密度； ℓ0为非局部参数，该参数可由室内波速试验测定，其包含孔隙尺寸效应及由于波动所引起的孔隙动力效应两部分[27]，主要反映土颗粒之间相互作用的影响范围；为曲度因子；η为流体粘滞系数；κ为渗透系数；F(ζ)为粘性修正系数，ω为 圆频率，；a为孔隙半径，圆孔状孔隙时：裂缝状孔隙时：ξ 为 弯曲因子；ber 和 bei分别为第一类零阶开尔文函数的实部和虚部。

1.2 边界条件

设有限元求解域为Ω，求解域的边界为，则式(1)应满足边界条件：

强制边界条件：

自然边界条件：

且各边界条件应满足如下关系：

1.3 空间离散

本章采用虚位移原理得到式(1)的空间离散方程，假设 δu为 土骨架的虚位移矢量， δw为流体相对固体虚位移矢量，则式(1)可表示为：

对变量场u和w引入同阶等参数单元差值函数：

式中：ue(x,t)和we(x,t)分别为单元土骨架和孔隙流体的位移场；x为空间坐标矢量；t为时间；分别为土骨架和孔隙流体节点位移矢量；Ni(x)为单元形函数；ne为单元节点数目；I为单位矩阵。和仅表示矩阵的乘积运算。

将式(7)中第一式和第二式分别代入式(5)和式(6)，采用Galerkin 法，可将式(5)和式(6)写成等效积分弱形式[34]：

其中：

2 非局部Biot 理论控制方程的时间离散

式(8)为非局部Biot 理论离散后的控制方程，本文考虑线弹性材料，因此则式(8)简化为：

为方便后续推导，将(Kss+K′ss)合并为一项，并令为其中弹性矩阵变为，，则式(9)可简化为：

为编程方便将式(10)写为：

为求解式(11)，将时域划分为时间间隔为Δt的小区间， Δt应足够小，以保证计算的收敛性和精度。由Newmark 积分可知，假设当前时间为n+1，并假设前一时刻n的位移、速度和加速度已知。则n+1时刻的位移、速度和加速度可表示为：

式(17)即为所求的离散方程，β1、β2为Newmark-积分常数，为保证计算稳定，应该满足β2≥β1≥0.5。本文计算所选取的计算参数值为 β1=0.6和β2=0.605，以引入一定的数值阻尼提高计算的稳定性及收敛性[35]。

3 循环荷载下饱和土动力固结分析

3.1 模型验证

为验证所编写u-w形式的非局部Biot 理论有限元程序的正确性，利用该程序分析饱和土一维固结问题，并将其与Terzaghi 一维固结解析解结果进行对比。为此，假设一个深度为H0=10.0m，宽度为1.0m 的一维固结模型(模型深度远大于宽度)，模型上边界设置为排水边界，顶部施加P0=100 kPa 的均布荷载，底部及两侧为不透水边界，计算模型及荷载形式如图1所示。为采用本计算程序模拟上述一维情况，将上述土柱采用4×40个4节点矩形单元离散，并固定所有节点的土骨架及孔隙流体相对土骨架的水平自由度(x方向)。为模拟土柱底部及两侧的不排水情况，将模型底部和两侧节点分别设置为wz=0及wx=0；顶部为排水边界，即wz所对应的自然边界条件；顶部外荷载P0全部施加到土骨架竖向(uz方向)自由度上。同时本文理论中的非局部参数τ=0.0 mm，以将非局部Biot 理论退化至经典Biot 理论控制方程。饱和土计算参数按JEREM IC等[36]论文选取，Terzaghi一维固结中假设了土颗粒及孔隙水压力不可压缩，因此将土颗粒及孔隙水的体积模量取一个大值，如表1所示。

图1 一维固结土柱及荷载形式Fig.1 One-dimension consolidation column and loading condition

表1 饱和土体参数Table1 parametersof saturated soil

图2为荷载作用下饱和土土柱一维固结的数值解与解析解的对比曲线，图中横坐标为无量纲超静孔压p/p0，纵坐标为无量纲深度h/h0。如前所述，Newmark 积分参数取为 β1=0.6 和 β2=0.605，即引入一定的数值阻尼以消除由于外荷载的突然施加而导致的数值波动。由图2可知，本文所编制程序的计算结果与Terzaghi 解析解结果吻合的很好，由此说明了本文所编制的u-w格式的饱和土体动力固结分析程序的正确性。

图2 静载下一维饱和土柱超静孔压的数值解与解析解对比Fig.2 Comparing the numerical solution and analytical resultsof the excess pore-pressure for one-dimension saturated column under static load

由于Biot 控制方程的复杂性，目前仍未有二维及三维问题饱和土动力响应问题的解析解，因此为进一步验证本文所编制u-w格式饱和土控制方程的正确性，将本文的饱和土双相介质退化为单相介质(通过固定孔隙流体相对固体位移自由度wi，同时取Biot 参数α =M=0)，并与已有数值结果对比[37]。为此，建立高10m，宽100m 的平面应变地基模型，该模型左上角表面1m 范围内施加均布荷载，均布荷载在时间T内成三角形变化，如图3所示。退化后土体参数E=100 MPa，ν=0.3 及 ρ=2000 kg/m3，荷载持续时间为T=0.16 s。模型左边为自由边界，底面为固定边界条件，计算结果分析了地基右边固定及右边自由两种工况下地基顶部中点A处的位移，并将计算结果与JU 和WANG[37]的论文结果对比。

图3 平面应变地基模型示意图/m Fig.3 Foundation model of plane strain

图4为模型右边界固定及自由情况下本文计算结果与JU 和WANG[37]的计算结果对比曲线，由图可知，右端固定和自由时，本文结果与JU 和WANG[37]的计算结果均吻合的很好，由此验证本文计算程序的正确性。

图4 模型顶面中点A水平位移ux 时程曲线Fig.4 The time-history curve of horizontal displacement ux of point A at the top surface of themodel

3.2 算例分析

为研究孔隙尺寸效应对循环动荷载下饱和土动力响应的影响，本节将选取几个典型的实例进行参数分析，以研究孔隙尺寸效应对饱和土地基的动力固结及动力响应的影响。具体如下：

3.2.1循环荷载下一维动力固结

循环荷载下饱和土的一维固结计算模型尺寸如图1所示，但模型顶部作用正弦荷载P=P0sin(ωt)，其中：P0=100 kPa为 正弦荷载幅值；ω=2πf为圆频率；f=1/T为正弦荷载频率；T为正弦荷载周期；位移观测点位于土柱中心，孔隙水压力观测点位于模型底面，数值积分采用Newmark 积分，时间间隔 Δt应随所施加荷载频率的变化而变化，如：频率较小时，时间间隔 Δt可选取较大值，而频率较大时选取的时间间隔 Δt应较小(本文计算时满足Δt≤T/8，即1/4荷载区域内应有包含两个时间积分点，以保证计算的精度)；模型左右边界约束固体颗粒及流体的法向位移，底部为固定边界和不透水边界条件；饱和土的物理力学参数见表2中的材料1。

表2 饱和土物理力学参数Table2 physical and mechanical parametersof saturated soil

图5为荷载频率不同时，观测点O的位移时域响应，表征孔隙尺寸效应的非局部参数取值分别为0.00 m、0.04m 和0.08 m。由图可知，荷载频率为500Hz 时，观测点O的位移响应随非局部参数的增加，变化不大。而由图5(b)～图5(d)可以看出，随非局部参数的增加，位移响应的峰值向后移动，即饱和土中波的相位差发生了变化。这是由于频率较小时，循环荷载在饱和土内产生的波的波长较长，此时可以认为饱和土中的波长远大于饱和土介质的孔隙尺寸，由此导致孔隙尺寸效应的影响并不明显。而随着荷载振动频率的增加，循环荷载在饱和土中产生的波的频率增加，此时，饱和土中的波长接近甚至小于饱和土介质的孔隙尺寸，由TONG 等[27]研究结果可知，考虑孔隙尺寸效应时，随波频率的增加，饱和土介质的波速减小，由此导致随非局部参数的增加，观测点O的位移响应向后移动。

图6为不同荷载频率时，观测点B的孔隙水压力随时间变化曲线。由图可知，荷载频率较低时，非局部Biot 理论的计算结果与经典Biot 理论计算结果几乎一致，而频率较高时，此时，随非局部参数的增加，改变了饱和土介质中波速的相位差，由此导致孔隙水压力响应稍有滞后。其原因与图5相同，在此不在赘述。此外，对比图6中的四幅图可看出，随入射频率的增加，整体孔隙水压力变小。

图5 不同荷载频率下一维饱和土地基位移响应时程曲线Fig.5 The displacement time-history curve of the onedimensional saturated foundation for different load frequency

图6 不同荷载频率下一维饱和土地基孔压响应时程曲线Fig.6 The pore-pressure time-history curveof the onedimensional saturated foundation for different load frequency

3.2.2循环荷载下二维饱和土地基动力响应

为充分研究孔隙尺寸效应对饱和土地基动力响应的影响，本节采用如图3的循环荷载下饱和土地基的计算模型尺寸，但在模型左端1m 范围内作用有均匀分布的正弦荷载，荷载周期为T，荷载频率取为1000Hz。模型底面为固定且不透水边界，左右边界为自由且排水边界(不对土骨架位移及流体相对土骨架位移施加任何约束)。观测点位于距离模型右侧10m 处的A点(地表)和模型表面中心B点。本实例饱和土所采用物理力学参数见表2中的材料2。

图7和图8分别为循环动荷载下观测点B和观测点A的固体位移及流体相对固体位移的时程曲线，由图可知，时间小于1.6×10-2s时，观测点B的位移响应为0，这是由于左端均布荷载产生的波还未到达观测点B，随时间的增加，观测点B开始波动，但考虑非局部效应，B点出现波动的时间稍有延迟，且随着时间的增加位移响应的延迟愈加明显。该现象也可从观测点A看出，所不同的是A点产生波动的时间大约在3×10-2s时刻。从图中还可以看出，考虑孔隙尺寸效应对饱和土介质中位移场的幅值影响不大，只是对饱和土介质中波速影响较大，由此导致位移场的相位差发生变化。由此说明，孔隙尺寸效应对循环荷载所引起地基中的能量影响较小。此外，由图中的插图可知，随时间的增加，位移场波动幅值略有增加的趋势，这是由于所施加的动载的持续作用，在介质中的能量积累所致。

图7 观测点A的位移时域响应Fig.7 Displacement response in time domain at point A

图8 观测点B 位移时域响应Fig.8 Displacement response in time domain at point B

3.2.3循环荷载下饱和土地基动力响应

本节例子所模拟的是一个饱和土地基，在其中心处作用一个均匀分布荷载P，荷载随时间呈正弦变化。考虑到该问题的对称性，仅考虑一半模型，如图9所示。约束模型两侧及底部的法向位移，且模型两侧及底面为不透水边界条件。模型的几何参数及网格划分如图9所示，物理力学参数见表3.2中的材料3。正弦荷载幅值为100 kPa，振动频率为1000 Hz，Newmark 积分步长选取为Δt=1.0×10-5s，即所分析的时长为3×10-2s。本节主要分析，非局部参数对循环荷载下饱和土地基位移及孔压响应的影响，观测点位于模型的四个顶点，顺时针方向依次为点1、点2、点3及点4(如图9所示)。

图9 条形荷载作用下地基响应计算模型Fig.9 Calculating modelof the responseof foundation subjected to a strip load

图10为循环荷载下观测点1和2处位移时程曲线，表征孔隙尺寸效应的非局部参数分别取为ℓ0=0.00 m 和ℓ0=0.08m。由图可知，观测点1的位移响应随时间的增加呈现出波动上升的趋势，而考虑非局部参数后，位移响应相对于经典Biot 理论(非局部参数取0)的计算结果有向后延迟现象，其原因是由于观测点1距离荷载作用区域有一定的距离，而考虑非局部效应后饱和土中波速有所降低所致。而观测点2位于荷载的下方，因此随时间增加位移场的相位差相差不大，但其振动幅值有一定程度的减小。

图10 位移时程曲线Fig.10 Displacement time-history curve

图11为循环荷载下，上述四个观测点的孔压随时间变化曲线，由图可看出，观测点1、观测点2和观测点4孔压响应开始时间有不同程度的延迟，其原因为各观测点与荷载的距离不同，导致其产生响应的时间不同。由于观测点2正好位移荷载下方，因此，孔压响应没有延迟现象，且考虑及不考虑孔隙尺寸效应所引起孔压响应的相位基本一致。此外，考虑孔隙尺寸效应后，观测点1的孔压振动幅值有一定的增大，而观测点2、观测点3和观测点4处孔压幅值则呈现出不同程度减小，且孔压响应均有滞后的现象。

图11 孔压时程曲线Fig.11 Pore pressure time-history curve

4 结论

本文基于非局部Biot 理论，采用虚位移法构建了有限元空间离散方程，采用Newmark 积分法对时间进行了离散，并采用MATLAB软件编制了饱和土地基二维动力固结的有限元计算程序，通过与解析解及已有的单相介质数值解的对比，验证了本文有限元计算程序的正确性。进而，采用该有限元程序分析了三种典型的工程实例，得出如下结论：

(1)荷载频率振动频率较低时，非局部Biot 理论与经典Biot 理论所预测的结果差别很小，即非局部参数的影响很小；

(2)当输入高频荷载时，孔隙尺寸效应的影响逐渐显现，且非局部参数对位移和应力响应的相位产生显著影响，使得位移场及孔压的响应时间有所延迟。这是由于考虑孔隙尺寸效应后，饱和土中的波速有所降低，由此导致荷载在饱和土中引起的波传播至观测点的时间有所延长，从而引起位移及孔压延迟响应的现象。

本文的数值模型可以用于研究如厂房机械振动、高速铁路、高速公路地基的动力响应及固结沉降分析等饱和土动力响应及固结沉降问题。