分数阶黏弹性模型在饱和软土与单桩相互作用中的应用

艾智勇，赵勇智，刘文杰

（1.同济大学土木工程学院，上海 200092；2.同济大学岩土及地下工程教育部重点实验室，上海 200092）

桩基础因其优良的承载特性，在软土地区得到了广泛的应用。目前，对于饱和地基与竖向受荷单桩的相互作用，已经取得了丰硕的成果，如文献［1-3］等。而对于软黏土，其固有的流变性质比较突出，地基的蠕变阶段较长，因此有必要考虑土体的流变特性来进行桩-土共同作用分析。王建华等［4］基于Merchant黏弹性模型，使用积分方程法得到了垂直受荷单桩与层状黏弹性地基的共同作用解答。吴文兵等［5］以广义Voigt模型来模拟土体的流变特性，研究了层状黏弹性地基中单桩变形的时间效应。Feng等［6］通过修正的Burgers模型，对多层黏弹性介质中的竖向受荷单桩进行了研究。Ai等［7］对横观各向同性黏弹性饱和土中轴向受荷桩的时变行为进行了研究。这些研究多采用Kelvin模型、Merchant模型等传统整数阶导数模型来表征地基的流变特性。然而，在某些情况下，特别是在蠕变或松弛的初始阶段，模型的应力-应变曲线与实验数据不能很好地吻合［8］。有关研究［9-10］表明，基于分数阶导数的本构模型可以更好地与实验结果相吻合，并且可以更好地模拟材料的黏弹性。为此，Ai等［11］根据分数阶导数的Merchant模型，采用精细积分法求得了多层横观各向同性黏弹性地基的固结解答，并证明了分数阶导数模型在模拟软黏土流变行为上的优势。

本文旨在引入分数阶Merchant模型，采用边界元与有限元耦合的方法，探讨分数阶黏弹性模型在饱和软土地基与单桩相互作用分析中的应用。

1 分数阶黏弹性地基的边界元解

1.1 分数阶黏弹性模型

根据Ai等［12］的研究，相比于Maxwell模型和Kelvin模型，Merchant模型更能体现软土的流变性质，因此选取分数阶的Merchant模型来描述地基的应力-应变关系。传统整数阶和改进的分数阶Merchant模型分别如图1a、1b所示，其中η和η1表示黏滞系数。二者均含有弹性模量分别为E0和E1的胡克弹性体H0和H1；不同的是，分数阶模型用Abel粘壶［13］代替了传统的牛顿粘壶。

图1 黏弹性本构模型示意图Fig.1 Sketch of viscoelastic constitutive model

对于传统Merchant模型，其的应力-应变关系为

式中：σ′s为有效应力；ε为应变；t为时间。

传统的牛顿粘壶满足牛顿摩擦定律，即σ′s(t)～d1ε(t)/dt1，它表示应力与应变对时间的一阶导数成正比。当用Abel粘壶替代传统的牛顿粘壶后，分数阶Merchant模型的应力-应变关系可表示为

式中：dαdtα为对时间t阶次为α(0＜α≤1)的导数。

对于分数阶导数，这里采取Riemann-Liouville的定义［14］。在区间(0，T)上，函数g(t)的α次微分为

式中：Γ(x)为Gamma函数，它的定义为

为简化推导，引入Laplace-Hankel变换和逆变换［15］

对式（2）进行Laplace-Hankel变换，并结合分数阶导数的定义，即式（3）-（5），可以推导得到

1.2 地基中任意点广义应力和位移的边界元解

根据弹性-黏弹性对应原理［16］，用式（8）获得的柔度系数代替相关的弹性参数，可获得横观各向同性分数阶黏弹性饱和软土地基的解［11］，即地基中任意α点受n向单位荷载时引起的s时刻β点m方向上地 基 位 移(α，β，s)及 超 孔 压(α，β，s)(m=r，z；n=z)，它们也是积分变换域内的地基柔度系数。

如图2所示将地基划分为N个层元，将桩划分为i+1个桩单元。其中，桩嵌入的土层部分划分与桩单元一致。按照边界元法的理论，把获得的地基柔度系数(α，β，s)及(α，β，s)作为边界元法的核函数，则对任意厚度为L的地基区域，在桩土边界面上沿深度方向积分，可得该区域内荷载作用下地基任意一点的位移和孔压。

图2 地基与单桩模型示意图Fig.2 Sketch of the model of the foundation and the single pile

式（9）、式（10）即为边界元法中的边界积分方程，也是地基内部沿深度方向受荷的边界元解答。

2 轴向受荷单桩的有限元解

采用2节点轴力杆单元模拟轴向受荷单桩，2节点轴力杆单元的刚度矩阵为

式中：Ep和Ap分别为桩的弹性模量和横截面积；l为杆单元长度。

根据有限元原理，可得以下有限元求解方程：

3 边界元与有限元的耦合求解

图2为多层横观各向同性黏弹性饱和地基与其中嵌入的由2节点单元所组成的单桩。地基划分为N个层元，单桩划分为i+1个单元，且桩嵌入的土层部分划分与桩单元一致。E0、E1、η和α为分数阶Merchant模型参数；νh为水平向应力引起的正交水平向应变的泊松比；νvh为竖直向应力引起的水平向应变的泊松比。

总的等效节点力矩阵p由外部载荷引起的等效节点力和桩与地基之间的相互作用力引起的等效节点力两部分组成，因此式（12）可以写成如下形式：

式中：U(t)即节点位移矩阵a，它表示其是随时间变化的；F(t)和Q(t)分别是由外部荷载和桩-土相互作用力（侧摩阻力和端阻力）引起的等效节点力。

如图2，在桩-土界面处，土内表面单元与桩单元对应，则边界单元节点与桩单元节点重合。设任意时刻边界某一单元内的任意一点的荷载与该单元2个节点上的边界力具有插值关系，即

式中：N=[N1，N2]；P ne(t)为由地基内表面单元e的2个节点上的n向边界力所构成的向量。

将积分区段L e取为边界单元长度，对式（14）进行Laplace-Hankel变换后代入式（9）和式（10）可得

式中：͂(L e，β，s)和͂(L e，β，s)分别表示L e段的所有荷载引起的β点的m向位移和超孔压。

对所有单元进行累加，可以得到由整个桩身范围的作用力引起的β点的位移及超孔压，为

式中：L为整个桩身桩长；i+1为总的单元数。

将β点依次取作单元节点，可以得到一系列关系式。在此基础上，它们统一用矩阵方程表示为

式中：͂(s)、͂(s)和͂(s)分别是由所有节点位移、超孔压和边界力所构成的总向量；͂和͂是由式（17）和式（18）构建的总柔度矩阵。

对于桩单元而言，根据有限元理论，地基反力引起的等效节点力向量为

将式（14）代入式（21），可以得到

对式（22）进行Laplace-Hankel变换，然后将2个方向的力向量合并到一起，并组装所有的桩单元，可以得到总的关系式，为

结合式（19）与式（23），可以得到

因桩-土单元划分是一致的，由桩-土界面处的位移连续条件可知

对式（13）进行Laplace-Hankel变换，并将式（24）、式（25）代入，可以得到外力与单元节点位移的关系，即

式中：͂(s)为单桩所受外力荷载向量。

将求得的͂(s)代入式（19）可得对应的边界力，再代入式（20）便得到对应时刻桩周土中的超孔压，而这些量的真实值可通过Laplace-Hankel逆变换得到。

4 数值计算与分析

因目前未见有分数阶黏弹性饱和地基与单桩共同作用的文献，为验证本文理论的正确性，将本文理论退化为基于整数阶Merchant模型的黏弹性饱和地基与单桩共同作用问题（分数阶次取1），并与王建华等［4］的结果进行对比。对比算例为嵌入3层饱和黏弹性地基中的竖向受荷单桩，长径比L p/D=25，桩顶荷载为V。第3层土为半无限空间，黏弹性土体参数选取如表1。其中，hi为第i层土的厚度，ki为渗透系数，E0i表示第i层土的Merchant模型参数E0，E1i表示第i层土的参数E1。为便于分析，定义量纲一化时间因子，量 纲 一 化 位 移w*=，量纲一化轴力

表1 饱和黏弹性地基参数Tab.1 Parameters of saturated viscoelastic soils

单桩沿深度变化的轴力和竖向位移经量纲一化处理分别绘于图3和图4。由图可见，桩身轴力和竖向位移随时间逐渐增大，但各土层对应的桩段上的轴力和竖向位移随深度的减少程度不同，这体现了土层分布的成层性对计算结果的影响。另外，本文解答与文献［4］吻合较好，这就验证了本文理论和数值计算方法的适用性。

图3 黏弹性地基中桩身轴力Fig.3 Axial force along the pile in viscoelastic soils

图4 黏弹性地基中桩身竖向位移Fig.4 Vertical displacement of the pile in viscoelas⁃tic soils

考虑到对桩身长径比、桩-土刚度比进行分析的文献已较为丰富，故本文仅针对性地探究Abel粘壶分数阶次α对地基-单桩相互作用的影响，并阐明其优越性。算例设计如下：大小为V的轴向荷载作用在双层黏弹性地基中的单桩顶部，桩的长径比，桩身模量为20GPa。双层地基厚度分别为；地基渗透系数为k1=1×10−7m⋅s−1，k2=1×10−8m⋅s−1；泊 松 比v1=0.25，v2=0.35；其余黏弹性参数如表2所示，其中，下 标i表 示 土 层 号，。为了便于分析，定义量纲一化参数如下：时间因子，位移

表2 黏弹性地基参数Tab.2 Parameters of viscoelastic soils

对于分数阶次如何取值，可照文献［9-10］中的方法，即：针对实际工程不同区域及流变特性不同的软土，获取土样进行室内试验，得到其应力-应变曲线；然后通过拟合方法选取合适的α值，以拟合其独特的流变性，并使计算结果与试验情况进一步吻合。对分数阶次的大概范围而言，根据分级加载单面排水流变固结试验与渗透试验，2组萧山软土土样分数阶次α分别为0.13和0.35［10，17］；为此，刘忠玉和杨强［10］在研究分数阶次对土体固结流变的影响时，分数阶次α分别取0.10、0.35、0.50和1.00。类似地，为研究分数阶次的影响，本文分别计算了α=0.1，0.6，1.0这3种情况下的桩身位移，并将不同τ时刻的位移曲线表示在图5中。其中，τ=0.001表示初始时刻，τ=108和τ=1015分别表示了中间时刻和固结流变完成时刻。

图5 不同分数阶次下桩身竖向位移Fig.5 Vertical displacement of the pile at different fractional orders

从图5可知，在不同的分数阶次下，无论是桩身初始位移（实线），还是最终位移（点划线）均相同，这说明分数阶次对桩身位移的初始值和最终值不产生影响；然而在中间时刻（虚线），随着流变的进行，不同分数阶次的桩身位移明显不同；其中，分数阶次α=1.0对应的桩身位移最大，而α=0.1对应的桩身位移最小。这是因为分数阶次越大，地基固结流变的速率越快，故桩身位移达到终态的时间也越短。

5 结语

为更好地描述软土地基与单桩相互作用的时变行为，引入分数阶Merchant模型，以探讨分数阶黏弹性模型在饱和横观各向同性地基与单桩相互作用分析中的应用，然后通过数值算例验证其适用性，并对分数阶次的影响进行了分析，得出了以下结论：

（1）分数阶次对桩身位移的初始值和最终值不产生影响。在相同的工况下，不同的分数阶次对应的初始沉降和最终沉降都相同。

（2）分数阶次显著影响着流变过程。分数阶次越大，地基固结流变的速率越快，故桩身位移达到终态的时间也越短。

本文所用分数阶黏弹性本构模型虽能更好地拟合土体的流变性质，但桩周土仍处于弹性阶段，并未考虑土体的塑性变形；另外，本文仅对单桩与土体的相互作用进行了分析，将来可进一步考虑桩-桩之间的相互作用，进而对群桩基础进行研究。

作者贡献声明：

艾智勇：提出理论分析方法，建立桩土相互作用模型，指导理论求解及修改，指导后续分析工作，指导文章撰写及修改工作。

赵勇智：参与理论求解及修改工作，参与算例设计及分析工作，参与完成文章初稿，进行稿件修改。

刘文杰：引入分数阶导数模型，参与理论求解工作，参与完成文章初稿。