自然界中的土体普遍具有流变性，流变使土的应力-应变关系产生时间依赖性(松弛、蠕变)。土体流变分为流体流变和土骨架流变，流体流变是指土体中孔隙水和胶结物在应力作用下发生的黏性流动变形，而土骨架流变是指固体颗粒在应力作用下的错动和重新排列^[1]。实际中的土体流变，是流体流变和土骨架流变的耦合作用结果。

饱和土的一维固结理论自从TERZAGHI于1925年首次提出，如今已有近百年的发展^[2]。期间，BIOT^[3]从弹性理论出发，考虑土体中应力、应变的相容条件和渗流连续条件，提出比较完善的三维固结理论。Terzaghi和Biot固结理论为研究土体固结行为奠定了基础，学者们在此基础上围绕土体非均匀性、各向异性、边界条件、排水条件及加载方式等方面做了拓展性研究，取得显著成果^{[4 − 6]}。

TAYLOR和MERCHANT^[7]最早提出用Kelvin流变模型考虑土骨架变形的时间效应。陈宗基^[8]将流变理论应用于土体固结分析，推导出了固结方程及其解答，指出流变效应在固结初期就已发生。XIE等^[9]采用三元件Merchant模型考虑土的流变特性，对循环荷载条件下的一维固结问题做了研究，指出循环荷载对固结速率具有显著影响。上述研究均为基于Terzaghi固结理论的一维固结问题，而且仅考虑了土骨架流变，未考虑孔隙水的流体流变。近些年，学者们通过在Biot固结理论中引入孔隙水的动力黏度，以考虑流体流变，结果表明：流体流变对土体固结响应的影响非常明显^{[10 − 11]}。

传统流变本构模型是关于时间的整数阶微分方程，研究中发现：整数阶流变模型在拟合试验结果时需要引入很多拟合参数才能实现土体固结曲线的有效拟合，导致力学模型非常复杂^{[12 − 13]}。近几年的研究表明：SCOTT BLAIR^[14]和HEYMANS等^[15]提出的分数阶导数流变本构对土体流变具有更好的适用性^{[16 − 18]}。由于分数阶微分算子为卷积积分，其具有全局相关性，能够准确地描述力学与物理过程中的历史记忆性。最新报道的研究中，学者们利用Caputo分数导数修正Kelvin模型，重新推导了饱和土的流变固结方程，并通过Hankel变换和Laplace变换给出了相应的解答，为分数阶导数模型在土体固结理论中的应用作出了重要贡献^{[19 − 21]}。

目前关于地基土三维固结问题的研究，大多针对圆形分布荷载作用下的轴对称情形，在柱坐标系下轴对称问题可得到简化处理。但实际工程中，荷载大多数情况下都并非圆形分布，因此轴对称固结模型的解答在工程中适用性有限。此外，现有研究中采用的分数阶流变模型多为Kelvin模型，该模型无法反映应力-应变关系的瞬时流变。基于上述文献综述，本文采用分数阶导数Merchant流变本构描述土骨架应力-应变关系，并将外部荷载考虑为矩形分布荷载；在三维直角坐标系中构建饱和土地基的非轴对称流变固结模型，采用双重Fourier变换和Laplace变换推导出变换域解析解；最后，通过数值反演得到地基中任意位置处的时域响应。

1   分数阶导数流变本构

为了反映土骨架的流变性，采用流变本构模型描述土骨架的应力-应变关系。流变理论中，描述流变效应的本构模型主要有Maxwell模型、Kelvin模型和Merchant流变模型，其中Maxwell模型为Hooke弹簧和Newton黏壶的两元件串联模型，在时间趋于无穷时，Maxwell模型的松弛模量趋于0，造成蠕变变形无限增大，因此Maxwell模型能很好地反映材料的松弛而不能真实反映蠕变；Kelvin模型为Hooke弹簧和Newton黏壶的两元件并联模型，在时间趋于无穷时，Kelvin模型的松弛模量趋于材料的静模量，尽管蠕变变形能够收敛于稳定值，但应力松弛过程中应力无法恢复至0，所以Kelvin模型能很好地反映材料的蠕变而不能真实反映松弛；Merchant模型是由Kelvin单元和一个Hooke弹簧串联而成的三元件模型，既能真实地反映松弛也能真实地反映蠕变。传统Newton黏壶用应变的一阶整数阶导数描述应变速率，KOELLER^[22]将黏壶应变速率替换为分数阶导数，进而提出Abel黏壶。本文将传统Merchant模型中的Newton黏壶替换为Abel黏壶，则得到分数阶导数Merchant流变固体模型，如图1所示。图1中：σ和ε为总应力和总应变；ε₁为单个弹簧元件的应变；ε₂为Kelvin单元的应变；k₁为单个弹簧元件的弹性模量，k₂为Kelvin单元的弹性模量；η_s和a为Abel黏壶的黏度和分数阶导数的阶数(0≤a≥1)。分数阶导数Merchant模型的应力-应变关系描述为：

式中：τ_σ=η_s/(k₁+k₂)为应力松弛时间；τ_ε=η_s/k₂为蠕变弛豫时间；E₀=k₁k₂/(k₁+k₂)为静弛豫模量；t为时间变量。从式(1)可以看出，当a=0时模型退化为线弹性体，a=1时模型退化为传统整数阶流变模型。

对任意时域函数f(t)引入如下Lapace变换：

式中：$\hat f\left( s \right)$为函数$f\left( t \right)$的Laplace变换，s为时间变量t的变换参数，上标“^”为对应物理量已进入Laplace变换域；ζ为任意实数；$\mathit{\mathit{\text{i}}}=\sqrt{-1}$为虚数单位。利用式(2)对式(1)执行Laplace变换，得到：

式(3)反映了分数阶导数流变模型在Laplace变换域的应力-应变关系，由式(3)得出复模量$\hat{E}_{\mathrm{r}}$的表达式为：

为讨论分数阶导数Merchant流变模型的模量随时间的变化规律，对式(4)进行Laplace逆变换。根据积分变换理论，有以下关系：

式中：符号L⁻¹为Laplace逆变换；Γ( )为伽玛函数。利用式(2)和式(5)对式(4)进行Laplace逆变换，得到E_r的时域表达式为^[23]：

式中，M_α,β( )为Mittag-Leffler函数，其定义为：

图  1  分数阶导数Merchant流变模型

Figure  1.  Fractional derivative Merchant rheological model

下载: 全尺寸图片 幻灯片

利用式(6)绘制出弛豫模量E_r随时间的变化曲线，如图2所示。图2反映出弛豫模量的瞬时响应以及随时间的递减性和收敛性。图2中曲线反映出，黏性阶数a=0时，流变模型退化为线弹性模型，此时弛豫模量不随时间变化，恒等于未弛豫模量E_u=E₀(τ_ε/τ_σ)；当a=1时分数阶导数流变模型退化为整数阶流变模型，该情况下弛豫模量E_r由E_u突变至E₀，无法反映弛豫模量随时间的递减过程，进一步表明分数阶导数流变模型具有反映瞬时流变的优点(历史记忆性)，而整数阶流变模型无此特点。对应于实际土体，导数阶数a越大，表征土体越软，流变性越强。

图  2  弛豫模量随时间的变化

Figure  2.  The variation of relaxation modulus over time

下载: 全尺寸图片 幻灯片

利用广义胡克定律，将式(3)中Laplace变换域的流变本构模型拓展到三维应力状态，写为：

式中：${\hat \sigma _x}$、${\hat \sigma _y}$、${\hat \sigma _{\textit{z}}}$、${\hat \tau _{xy}}$、${\hat \tau _{y{\textit{z}}}}$、${\hat \tau _{x{\textit{z}}}}$为Laplace变换域中的土骨架应力分量(有效应力)；${\hat \varepsilon _x}$、${\hat \varepsilon _y}$、${\hat \varepsilon _{\textit{z}}}$、${\hat \gamma _x}$、${\hat \gamma _y}$、${\hat \gamma _{\textit{z}}}$为应变分量；$\hat \theta = {\hat \varepsilon _x} + {\hat \varepsilon _y} + {\hat \varepsilon _{\textit{z}}}$为土骨架体应变；$\hat G$为剪切模量；$\hat \lambda$为Lamé常数，根据定义有：

式中，μ为泊松比。

2   半空间饱和黏土地基固结模型

如图3所示，在三维直角坐标系x-y-z中构建饱和黏土地基的无限半空间模型，其表面作用矩形分布压力q(x,y,t)，荷载分布长度为l_a，宽度为l_b。局部压力q(x,y,t)为t=0时刻突加的恒定荷载，q(x,y,t)的表达式为：

式中，H( )为Heaviside阶跃函数，其定义为：

图  3  饱和土地基非轴对称固结模型

Figure  3.  Non axisymmetric consolidation for saturated ground

下载: 全尺寸图片 幻灯片

2.1   饱和黏土地基固结控制方程

1) 饱和黏土固结平衡方程

根据Biot固结理论，写出固结平衡方程为：

式中：p为超静孔隙水压力，后文中均简称为孔压；$\gamma ' = \left( {{\rho _{{\mathrm{sat}}}} - {\rho _{\rm{f}}}} \right)g$为土的有效重度，其中${\rho _{{\text{sat}}}}$为土体饱和密度，${\rho _{\rm{f}}}$为水的密度，g为重力加速度。

2) 土骨架变形的几何方程：

式中，${u_x}$、${u_y}$、${u_{\textit{z}}}$为土骨架位移分量。

3) 达西渗流连续方程

对于饱和土，孔隙水的流出量即等于土骨架的体应变，根据达西定律得到渗流连续方程为^[24]：

式中：${\nabla ^2} = \dfrac{{{\partial ^2}}}{{\partial {x^2}}} + \dfrac{{{\partial ^2}}}{{\partial {y^2}}} + \dfrac{{{\partial ^2}}}{{\partial {{\textit{z}}^2}}}$为Laplace算子；k_f为渗透系数。

2.2   边值条件

1) 初始条件

2) 边界条件

考虑地表为透水边界，则矩形分布荷载作用下的地表边界条件为：

无限半空间无穷深度处的边界条件为：

3   模型求解

饱和黏土地基流变固结模型的控制方程，由一系列偏微分方程组成。为实现偏微分方程组的求解，对空间坐标x、y引入如下双重Fourier变换：

式中：$\tilde {\tilde f}\left( {\xi ,\eta } \right)$为函数$f\left( {x,y} \right)$的双重Fourier变换；$\xi$、$\eta$分别为x、y坐标的变换参数；上标“≈”为对应物理量已进入双重Fourier变换域。

利用式(2)和式(16)对式(11)、式(12)同时进行双重Fourier变换和Laplace变换，得到：

经过双重Fourier变换和Laplace变换后，平衡方程和几何方程都变为关于z坐标的常微分方程。将式(8)和式(19)代入式(18)后得到：

根据体应变的散度性质和Fourier变换的微分性质，有：

对式(20)~式(22)取散度，并考虑式(23)，即$\dfrac{\partial }{{\partial x}}\left( {20} \right) + \dfrac{\partial }{{\partial y}}\left( {21} \right) + \dfrac{\partial }{{\partial {\textit{z}}}}\left( {22} \right)$，整理后得到：

同理，对式(13)进行双重Fourier变换和Laplace变换，得到：

将式(25)代入式(24)消去$\hat {\tilde{ \tilde p}}$后得到：

式(26)为$\hat {\tilde{ \tilde \theta}}$关于z坐标的一元二阶常微分方程，直接利用常微分方程的基本理论，得出式(26)的通解为：

式中：${r_1} = \sqrt {{\xi ^2} + {\eta ^2} - \varphi }$，其中$\varphi = {{{\rho _{\rm{f}}}gs} /{[ {{k_{\rm{f}}}( {\hat \lambda + 2\hat G} )} ]}}$；C₁、C₂为待定系数。经过分析可知，${{\mathrm{Re}}} ( {{r_1}} ) {\geqslant} 0$。根据边界条件式(16)，${\textit{z}} \to \infty$时的变形为0，因此得出式(27)中的C₂=0。

将式(27)代入式(25)，整理后得到$\hat {\tilde{ \tilde p}}$的常微分方程：

根据非齐次常微分方程理论得出式(28)的通解为：

式中：${r_3} = \sqrt {{\xi ^2} + {\eta ^2}}$；C₃为待定系数。

将式(27)、式(29)代入式(20)~式(22)，依次解得位移分量${\hat {\tilde{ \tilde u}}_x}$、${\hat {\tilde{ \tilde u}}_y}$、${\hat {\tilde{ \tilde u}}_{\textit{z}}}$的通解：

式中，C₄、C₅、C₆为待定系数。同理，将式(27)和式(30)代入式(8)和式(19)，进一步解出应力分量的通解为：

式(27)中已得出C₂=0，接下来只需利用地表边界条件确定出待定系数C₁、C₃、C₄、C₅、C₆问题即可得解。将式(29)、式(30)、式(31)代入边界条件式(15)及式(23)，同时取z=0，得到关于C₁、C₃、C₄、C₅、C₆的5个线性方程组：

式中，$\hat {\tilde {\tilde q}} = {q_0}\dfrac{{\sin \left( {{l_{\rm{a}}}\xi } \right)\sin \left( {{l_{\rm{b}}}\eta } \right)}}{{{l_{\rm{a}}}{l_{\rm{b}}}\xi \eta s}}$为荷载函数q(x,y,t)的Fourier-Laplace变换。

求解式(32)，得出C₁、C₃、C₄、C₅、C₆的结果为：

式中：$\varLambda = 2{r_1}{r_3} + {{\hat \lambda \varphi } / {2\hat G}}$；${F_1} = 2r_3^2[ {( {\hat G - \hat \lambda } )\varphi - \hat G{{( {{r_1} - {r_3}} )}^2}} ]$；${F_2} = 2\hat G{r_3}{( {{r_1} - {r_3}} )^2}$；${F_3} = \varphi [ {\hat \lambda ( {r_3^2 + {r_3}} ) - 2\hat G} ]$。

经过上述推导，得出了饱和土流变固结的Fourier-Laplace变换域解析解，利用IFFT和Crump法^{[25 − 26]}进行双重Fourier和Laplace数值反演，则进一步得到地基中任意点处的时域固结响应。

4   结果分析

4.1   方法验证与对比分析

由式(1)分析可知，在本文模型基础上取分数阶导数a=0时，模型退化为Biot固结模型，MCNAMEE和GIBSON^[27]已给出了Biot三维固结的经典解。基于MATLAB平台进行编程计算，图4给出了本文退化结果与经典解的对比，同时图中也与Terzaghi一维固结的结果做了对比，模型计算参数如表1所示。为了实现标准化对比，计算结果中对物理量做了如下无量纲处理：无量纲孔压${p^*} = {p /{{q_0}}}$，无量纲时间${t^*} = {{{k_{\rm{f}}}{G_0}t} / {\left[ {{\rho _{\rm{f}}}g\left( {l_{\rm{a}}^2 + l_{\rm{b}}^2} \right)} \right]}}$，其中${G_0} = {{{E_0}} / {\left[ {2\left( {1 + \mu } \right)} \right]}}$。图4的对比结果显示，本文退化解与经典Biot固结的结果完全吻合。从图4中能够清晰看出，Biot三维固结能够反映出土体固结初期的Mandol-Cryer效应，而Terzaghi一维固结无法反映该物理现象。此外，Terzaghi固结和Biot固结在加载瞬间，孔压随即达到高水平位置，无递增过程。

图  4  与经典固结理论的对比结果

Figure  4.  Comparison results with classical consolidation theory

下载: 全尺寸图片 幻灯片

表  1  模型计算参数

Table  1.  Calculation parameters

参数	单位	数值
静弛豫模量E0	N·m−2	1.16×108
泊松比μ	−	0.29
土体饱和密度ρsat	kg·m−3	2600
水密度ρf	kg·m−3	1000
渗透系数kf	m·s−1	1×10−7
黏壶黏度ηs	Pa·s	1×108
恒定荷载q0	N·m−2	100
荷载分部长度la	m	1
荷载分部宽度lb	m	1

下载: 导出CSV 

| 显示表格

图5进一步对比了Biot固结与本文流变固结的差异性。图5中显示，本文流变模型取a=0的结果与Biot固结的结果完全相同，随着a的取值增大，加载时的瞬时孔压降低，加载初期出现一段孔压递增过程，a的取值越大，加载时的瞬时孔压越低，且递增时间越长。而经典Biot固结理论，无法反映加载时的瞬时孔压，也不能反映加载初期的孔压递增过程。此外，图5的对比结果表明，土骨架流变性只影响孔压递增和消散的时间路径，但不改变固结过程中的孔压峰值。图6为饱和黏土固结的室内试验结果，通过对比图5的计算结果能够发现，本文模型a=0.6和a=0.8的理论计算曲线与试验曲线的形态非常相似，文献[28]中的试验土样为南京地区的饱和软黏土，具有较强的流变性，因而a值较大。

图  5  流变固结与Biot固结的差异性对比

Figure  5.  Comparison with Biot consolidation

下载: 全尺寸图片 幻灯片

图  6  饱和土固结试验曲线^[28]

Figure  6.  Test curves of saturated soil consolidation

下载: 全尺寸图片 幻灯片

上述对比验证表明：本文分数阶导数流变固结模型能够完备的退化到Biot固结理论，而且能够反映加载瞬间的瞬时孔压以及加载初期的孔压递增过程，所反映的孔压响应规律与试验结果一致^{[29 − 30]}。

4.2   算例分析

图7给出了地基不同深度处的孔压消散曲线，反映了地基不同深度处土体固结响应的差异性。图7中显示，随着深度增加，土体中孔压开始消散的时间越为滞后。图8的对比结果进一步表明，土骨架的流变阶数a增大，加载时地基不同深度处的瞬时孔压差减小(${\left. {\Delta {p^*}} \right|_{a = 0.2}} \gt {\left. {\Delta {p^*}} \right|_{a = 0.6}}$)。

图9为深度z=l_a平面上不同水平距离处的孔压响应曲线，由图9可知，地基中水平方向不同位置处的土体孔压响应同样具有滞后性，距离荷载作用位置越远，孔压开始消散的时间越为延后。图10显示，当土骨架的黏性阶数a增大时，加载时地基水平方向的瞬时孔压差也将减小。图11给出了地表沉降的时间发展曲线，图11中$u_{\textit{z}}^*$为无量纲竖向位移，$u_{\textit{z}}^* = {{{G_0}{u_{\textit{z}}}} /{( {{q_0}\sqrt {l_{\rm a}^2 + l_{\rm b}^2} } )}}$。图11中显示：土骨架流变性越强加载时的瞬时沉降越小，沉降发展时间越长。

图  7  不同深度处土体的固结响应

Figure  7.  Consolidation response of soil at different depths

下载: 全尺寸图片 幻灯片

图  8  流变对不同深度处土体固结的影响

Figure  8.  Influence of rheology on soil at different depths

下载: 全尺寸图片 幻灯片

图  9  不同水平位置处土体的固结响应

Figure  9.  Consolidation response at different horizontal positions

下载: 全尺寸图片 幻灯片

图  10  流变对不同水平位置处土体固结的影响

Figure  10.  Influence of rheology on different horizontal position

下载: 全尺寸图片 幻灯片

图  11  流变性对地表沉降的影响

Figure  11.  Influence of rheology on surface subsidence

下载: 全尺寸图片 幻灯片

5   结论

本文利用分数阶导数Merchant流变本构描述土骨架应力-应变关系，结合Biot固结理论构建饱和土地基在矩形分布荷载作用下的三维非轴对称流变固结模型，采用双重Fourier变换和Laplace变换推导模型控制方程的变换域解析解，并利用数值反演得出时空域解。通过与经典固结理论和试验结果进行对比，验证了模型和算法的适用性，分析了分数阶导数模型的历史记忆性描述土体流变的优越性。本文基于Biot固结理论构建分析模型，仍属于弹性框架下的理论，所以只适用于小变形情况下的土体固结分析，对于包含塑性变形的大变形固结则不再适用。通过本文分析，得出以下几条结论：

(1) 传统固结理论无法反映土体在加载时的瞬时固结响应，施加荷载后，孔压和变形即刻达到高水平位置；而流变固结能够反映加载后的孔压递增过程，流变性越强，初始瞬时孔压和变形越小且递增时间越长。

(2) 流变固结考虑了土骨架变形的时间效应，因此地基不同位置处的土体固结响应出现明显的时间滞后性，离荷载作用点越远，滞后性越明显。

(3) 土骨架流变性使加载时不同位置处的瞬时孔压差减小，在孔压递增阶段孔压差逐渐增大。流变性只改变了孔压响应的时间路径，但不改变固结过程中的孔压峰值。

(4) 土骨架流变性对地表瞬时沉降和沉降的发展时间都具有显著影响，流变性越强，瞬时沉降量越小，沉降发展时间越长。