桩基础作为土木工程中最常用的基础结构形式，因其能较好地适应各种地质条件及荷载情况，且具有承载力高、沉降小等优点，被广泛应用于各类土木工程结构中^[1-2]。在桩基础研究与设计中大多只考虑其竖向静荷载承载能力，而随着近岸与海洋资源的开发，在深水环境条件下，开敞式港口码头、海上采油平台等工程结构除承受静荷载外，在极端波浪、风以及地震等环境荷载作用下，结构还承受较大的水平动力荷载^[3-7]。尤其在地震荷载作用下，饱和地基中超孔隙水压力的上升能显著加速地基承载力的下降，直至液化临界状态后完全丧失承载力，造成结构物失稳破坏^[8-10]。在过去数十年发生的主要地震中，如日本的Niigata地震，Loma Prieta地震，Kobe地震等，由地基液化所导致的港口建筑、桥梁以及各类海洋平台结构的破坏案例均有所报道^[11-12]。因此，研究地基液化产生机理，确保地基拥有足够的承载能力已成为结构抗震设计中的重要课题。

考虑到超孔隙水压的上升作为土体骨架与孔隙水耦合作用的结果^[8]，因此准确模拟土体在循环荷载作用下的非线性响应成为解决地基液化问题的必要前提。在土的弹塑性模型发展史上，基于临界状态理论^[13]发展而来的剑桥模型是第1个基于试验且能够准确描述正常固结黏性土力学特性的经典土力学本构模型^[14]，并被广泛应用于土木与岩土工程实践中^[15]，但由于剑桥模型参数均基于正常固结重塑土三轴试验得到，因此严格意义上只适用于正常固结重塑土。然而在地震荷载或其他动力循环加载作用下，土体经历弹性卸载后即进入超固结状态。作为土体应力历史的反映，其对当前土体的力学特性有着不可忽略的影响^[15]。为研究超固结土的力学特性，Nakai^[16]与Asaoka^[17]等学者通过引入下负荷面(subloading surface)概念，在剑桥模型基础上通过考虑超固结因素的影响提出了下负荷面剑桥模型，为研究循环交变荷载作用下的超固结演化规律以及在超固结状态影响下的土体力学特征提出了一个可行的理论框架^[18]。

鉴于饱和地基地震液化的巨大影响以及土体循环加载下的复杂非线性特性，本文分别通过排水以及非排水循环加载试验对超固结土体强度以及超固结状态的发展规律进行了详细分析。并在此基础上，通过引入下负荷面剑桥模型本构关系，建立了饱和地基环境条件下的桩-土耦合体系在地震作用下的三维非线性动力有限元数值模型，通过数值模拟研究了饱和地基条件下桩柱结构的地震动力响应特性，对超固结因素的影响作用进行了详细对比分析。

1 超固结下负荷面剑桥模型

在剑桥模型以及其他经典土力学弹塑性理论中，在土的卸载及再加载过程中，其应力应变关系被假定为弹性关系(图 1(a))，而实际上即使在再加载过程中也会产生塑性应变。同时考虑到正常固结土一旦卸载，即变成超固结状态(图 1(b))。在经典弹塑性理论中，只存在一个屈服面，即正常固结屈服面(图 2(a))，若某一点应力状态处在此屈服面内(超固结状态)，则假定其为弹性，这与实际试验观测现象不符。为解决这个问题很多学者提出了不同的模型。其中下负荷面模型的物理意义最为明确^[14]，可用状态变量的概念精确定义塑性势，且假定下负荷屈服面经过现有应力点并和正常屈服面几何相似，而且随应力变化而变化，因此也不需要判定应力状态是否到达屈服面，简化判定标准。

尽管基于重塑黏土三轴试验发展而来的剑桥模型的应用范围存在局限性，未能考虑超固结因素的影响，但仍然成为现代诸多弹塑性理论模型的基础。剑桥模型采用塑性体应变硬化规律，以塑性体积应变ε_V^p为硬化参数，且假定塑性变形符合相关联流动法则，其屈服函数为^[14]：

$ f = f\left( {p, q, \varepsilon _{\rm{V}}^{\rm{p}}} \right) = {C_{\rm{p}}}\ln \frac{p}{{{p_0}}} + {D^*}\eta-\varepsilon _{\rm{V}}^{\rm{p}} = 0 $

(1)

式中：ε_V^p为硬化参数；C_P=(λ-κ)/(1+e₀)，e₀为参考大气压p₀(98 kPa)下的孔隙比(图 1)，λ=0.434C_c，κ=0.434C_s，C_c和C_s分别为土体的压缩系数与膨胀系数；D^*=C_P/M^*；M^*为临界状态下的剪切应力比；p和q为平均有效应力和剪应力，η=q/p为剪切应力比。

对图 2(b)中处于超固结状态的应力点P(p, q)，由于下负荷剑桥模型假设下负荷屈服面与正常固结屈服面几何相似，则通过P点的下负荷屈服面f_s可表示为：

$ {f_{\rm{s}}} = {C_{\rm{p}}}\ln \frac{p}{{{p_{{\rm{N1}}}}}} + {D^*}\frac{q}{p} = 0 $

(2)

式中:p_N1为下负荷屈服面与p轴的交点，而p_N1e则为正常固结屈服面与p轴的交点。

根据应力关系，下负荷屈服面f_s可写为：

$ {f_{\rm{s}}} = {C_{\rm{p}}}\left[{\ln \frac{p}{{{p_0}}}-\left( {\ln \frac{{{p_{{\rm{N1e}}}}}}{{{p_0}}}-\ln \frac{{{p_{{\rm{N1e}}}}}}{{{p_{{\rm{N1}}}}}}} \right)} \right] + {D^*}\frac{q}{p} = 0 $

(3)

通过超固结状态变量$ R\left( {R = \frac{{{p_{{\rm{N1}}}}}}{{{p_{{\rm{N1e}}}}}}} \right)$的引入以及利用正常固结屈服面与下负荷面间应力关系，得到下负荷屈服面方程：

$ {f_{\rm{s}}} = \ln \frac{p}{{{p_0}}} + \frac{1}{{{M^*}}}\frac{q}{p}-\frac{1}{{{C_{\rm{p}}}}}\varepsilon _{\rm{V}}^{\rm{p}}-\ln R = 0 $

(4)

将方程变换到一般应力空间：

$ {f_{\rm{s}}} = \ln \frac{{{\sigma _m}}}{{\sigma _{m0}^*}} + \frac{{\sqrt 3 }}{{{M^*}}}\frac{{\sqrt {{J_2}} }}{{{\sigma _m}}}-\frac{1}{{{C_{\rm{p}}}}}\varepsilon _{\rm{V}}^{\rm{p}}-\ln R = 0 $

(5)

式中：σ_m=σ_ii/3，J₂=s_ijs_ij/2，s_ij=σ_ij-σ_m。对于超固结状态参数R的发展，根据Nakai和Asaoka等^[16-17]的建议，取dR=Udε_d^p，其中dε_d^p为塑性剪切应变增量，$ U =-{m_{\rm{R}}}\ln R\frac{{{M^*}}}{{{C_{\rm{p}}}}}, {m_{\rm{R}}}$为与超固结相关的材料参数。

通过协调方程联合Hooke定律可得到下负荷面剑桥模型的弹塑性本构关系(其中E_ijkl为弹性刚度张量，ε_kl，ε_kl^e和ε_kl^p分别为总应变张量、弹性应变张量和塑性应变张量)：

$ \begin{array}{l} {\rm{d}}{\sigma _{ij}} = {E_{ijkl}}{\rm{d}}\varepsilon _{kl}^{\rm{e}} = {E_{ijkl}}\left( {{\rm{d}}{\varepsilon _{kl}}-{\rm{d}}\varepsilon _{kl}^{\rm{p}}} \right)\\ \;\;\;\;\;\;\; = {E_{ijkl}}{\rm{d}}{\varepsilon _{kl}}-{E_{ijkl}}\Lambda \frac{{\partial f}}{{\partial {\sigma _{kl}}}} \end{array} $

(6)

$ \mathit{\Lambda = }\frac{{\frac{{\partial f}}{{\partial {\sigma _{ij}}}}{E_{ijkl}}{\rm{d}}{\varepsilon _{kl}}}}{{{h^{\rm{p}}} + \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{E_{ijkl}}\frac{{\partial f}}{{\partial {\sigma _{kl}}}}}}, {h^{\rm{p}}} = \frac{1}{{{C_{\rm{p}}}}}\left\{ {\frac{{\partial f}}{{\partial {\sigma _{mm}}}}-\frac{{{m_{\rm{R}}}\ln R}}{R}\frac{1}{{{\sigma _m}}}} \right\} $

通过模型本构关系的建立，得到以下加载准则：①加载时，|dε_ij^p|>0，Λ>0及${\rm{d}}{f_\sigma } = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\rm{d}}{\sigma _{ij}} > 0 $(应力应变曲线硬化过程); Λ>0及${\rm{d}}{f_\sigma } = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\rm{d}}{\sigma _{ij}} < 0 $(应力应变曲线软化过程); ②卸载时，|dε_ij^p|=0，Λ≤0。

2 模型验证及对比研究分析

为验证模型的可靠性以及研究超固结因素对土体力学和液化特性的影响，本文在下负荷面剑桥模型下，分别通过排水与非排水循环加载试验进行了数值模拟试验研究。试验土以藤森黏土为例^[17]，主要参数λ=0.09，κ=0.02，e₀=0.88，M^*=1.36，m_R=2.1，泊松比ν=0.2，γ=20.4 kN/m³。

2.1 单调加载下数值模拟验证

图 3为平均主应力一定条件下(98 kPa)的三轴压缩试验结果^[17]与下负荷面剑桥模型预测值对比结果，其中OCR为超固结比。

由图 3可见，下负荷面剑桥模型很好地模拟了超固结土的应力软化及剪胀现象。4种不同超固结状态土的屈服强度计算值分别达到1.36，1.53，1.72和1.92(图 3(a))。相较于重塑土，超固结土的屈服强度得到显著提高，屈服强度与超固结比间表现出明显的正相关性。而在加载屈服过程中，随着加载过程中超固结比的衰减(图 3c)，土体的加载响应开始软化，并在完全屈服后均趋于同一应力水平。图 3(b)给出了加载屈服过程中土体的体积压缩过程曲线。从图中压缩曲线的发展规律可以看出，在加载初始阶段，所有土体在剪切应力作用下均经历了1个明显的体积压缩过程。其中正常固结土的孔隙比压缩量达到3.5%，而对于超固结土，随着加载过程中超固结状态的不断衰减(图 3(c))，使得孔隙比朝着增长的方向发展(Δe)，并在屈服后表现出明显的剪胀特性，尤其在高超固结比条件下，其剪胀现象表现得更为显著(Δe=1.19%, 2.41%, 3.21%)。

2.2 循环加载数值模拟试验

通过单调加载试验对比研究可以看出，通过考虑超固结因素的影响，下负荷面剑桥模型对超固结土在单调加载下的屈服过程和剪胀现象能够做出准确描述。而在循环交变荷载作用下，考虑到土体的屈服过程同时伴随着超固结状态的演化发展，尤其在饱和非排水加载条件下，伴随着超孔隙水压的升高与有效应力的减小能加速土体的软化进程，直至达到液化状态后完全失去承载能力。为进一步研究土体的循环加载响应与液化特性，以及初始超固结比的影响作用，本文分别在排水和非排水条件下对3种超固结比条件下(OCR=1，2，4)土体的循环交变荷载作用下的应力应变响应关系、可压缩性以及相应的孔压、超固结比状态的发展规律进行了数值对比试验研究。试验加载幅值σ₁-σ₃为100 kPa，应变加载步长Δε_d=0.002，总试验步数为2 000步。

图 4即为三轴排水条件下的循环加载试验结果，从图中试验结果的对比可以看出，相较于正常固结土，随着超固结比的提高，超固结土在循环加载过程中表现出更小的塑性变形与非线性特性，且在循环加载过程中伴随着超固结状态逐渐达到稳定状态后(图 4(c))，其应力应变关系也随之达到稳定状态。同样从土体压缩过程曲线(图 4(b))也可以看出，在正常固结土中，孔隙比压缩量达到0.18，而在超固结土下其压缩量分别为0.12和0.06，相较于正常固结土，超固结土表现出更低的可压缩性。

图 5为非排水条件下循环加载试验结果，与排水试验条件下所有土体均被压缩至最终稳定状态不同，在非排水条件下随着超孔隙水压的上升，尤其在正常固结条件下，土体应力应变响应曲线表现为更大的回滞圈。图 5(b)给出了循环加载过程中超孔隙水压的发展过程曲线，从图中对比可见，由于正常固结土的高可压缩性，使得土体加载应力转移至孔隙水上，导致土体有效应力随之减小(图 5(c))，并最终影响土体应力应变关系响应曲线。同样从超固结比的演化过程曲线(图 5(d))的对比也可以看出，由于非排水条件加载下土体无法压缩至最终稳定状态，在循环加载过程中相应的超固结比演化在循环荷载作用下表现出更加明显的波动过程。

3 桩柱结构地震响应数值算例

通过对超固结因素的研究可以看出，土体的超固结状态能够显著影响土体的力学特性，尤其在循环非排水条件下考虑土体超固结因素后能显著影响土体的超孔隙水压上升速率、加速土体的软化过程。为进一步研究超固结因素在饱和地基下桩柱结构地震响应中的影响作用，在ADINA81有限元程序的基础上开发了剑桥模型和下负荷面剑桥模型，并以一简化四桩柱码头结构为例对结构的非线性地震响应进行了模拟研究。码头计算原型为离岸深水港岩基浅埋轻型的四桩钢管结构码头(6(a))^[19]，在本文的数值模型中，为简化计算模型，仅保留了钢管桩码头的主桩结构，并对上部结构的荷载进行了简化。其中码头结构模型上部结构距离泥面35 m，桩柱在地面以下部分长20 m，结构自重2 000 kN，流体的影响以附加水质量的形式简化考虑。钢管桩桩径为2.8 m，壁厚32 mm，弹性模量210 GPa，材料密度7 800 kg/m³，桩柱结构成18 m×24 m布置(图 6)。地基土以藤森黏土为例，计算域取100 m×40 m×30 m，在截断面上采取黏弹性边界以模拟半无限空间中地震波的运动。输入地震波采用Loma波^[9]进行模拟计算(图 7)。为更好对比研究结构的非线性地震响应，本文数值模拟中将Loma波放大至0.2g作为入射波进行对比计算。

以下负荷面剑桥模型计算结果为例给出了震后地基中超孔隙水压比的分布云图(见图 8(a))。从图中超孔隙水压分布可以看出，受桩土间动力耦合作用的影响，桩侧附近土体的超孔隙水压显著高于远场地基，且随着深度的增加，桩柱与自由场间运动差异减小，相应的超孔隙水压比差异随之降低。为进一步对比研究不同模型下超孔隙水压分布规律，图 8(b)给出了桩柱附近以及地基远场在深度方向上的超孔隙水压比分布。从图中曲线对比可以看出，超孔隙水压比在地表处达到最大值，尤其在近场处受桩土动力耦合作用的影响，地表附近已达到液化状态。而随着地基深度的加深，超孔隙水压比呈现出明显的指数衰减趋势。在5 m深处，在剑桥模型和下负荷面剑桥模型下近场超孔隙水压比分别达到0.91和0.60，随着深度的继续增加，在20 m深度位置相应的超孔隙水压比分别降低为0.36和0.30。而通过不同模型地基下超孔隙水压曲线对比可以看出，在考虑超固结影响因素后，在地震作用下地基表现出了更好的抗液化能力。

图 9给出了不同模型下结构顶点动力响应对比，为更好对比研究结构的非线性地震响应，本文同时在理想弹性地基(E=3(1-2ν)(1+e₀)p/κ)假设下对结构的地震响应进行了数值模拟，并将其计算结果作为参考背景对结构的非线性响应进行了对比分析。从图中对比结果可以看出，结构的地震动力响应峰值加速度分别达到18.36, 12.91和14.34 m/s²，相应的位移峰值为0.29, 0.41和0.34 m。相较于理想弹性地基，在考虑土体非线性因素后，随着地基土的屈服，结构模型的整体刚度降低，尤其在剑桥模型下表现得更加明显，结构加速度响应更小，而由于地基承载力的下降，结构在地震作用下表现出更高的位移响应，并在震后表现出更显著的塑性位移(0，0.07和0.04 m)。

为进一步研究地震作用下桩柱结构的受力情况，通过桩柱截面上的应力积分，图 10给出了桩柱的弯矩包络图。同样为方便对比桩柱结构中的弯矩，图中弯矩值均以弹性地基计算结果作为参照基准以相对值给出。从图中弯矩曲线可以看出，在整个桩身中弯矩包络曲线存在两个明显的峰值，且在泥面附近处达到最大值。相较于理想弹性地基模型，剑桥模型与下负荷面剑桥模型地基下最大弯矩值分别达到1.36和1.10。通过本文结构算例的地震动力响应与受力结果对比可以看出，在考虑超固结因素后，下负荷面剑桥模型模拟计算结果介于弹性模型与剑桥模型之间，相较于剑桥模型，地基表现出了更好的抗液化与承载能力。

4 结语

通过考虑超固结因素的影响，本文在下负荷面剑桥模型下对超固结土力学强度、可压缩性，以及非排水循环加载试验下超孔隙水压的发展规律进行了数值模拟研究。并在ADINA81程序基础上通过剑桥模型和下负荷面剑桥模型本构关系模块的开发，对饱和地基条件下桩柱结构的地震响应进行了有限元数值模拟研究，得到以下主要结论：

(1) 通过对超固结状态演化过程的模拟，下负荷面剑桥模型能够准确模拟出超固结土单调加载屈服过程中的应力软化以及体积剪胀等特性。相较于正常固结土，随着超固结比的增大，土体的屈服强度也相应得到显著提高，但在土体最终屈服退化为正常固结状态后，超固结土与正常固结土的屈服应力趋于一致。

(2) 在排水循环荷载作用下，随着土体被压缩至其稳定密实状态的同时达到更高超固结状态，且相较于正常固结土，由于超固结土较小的初始孔隙比，在循环加载下表现出更低的可压缩性; 在非排水循环加载下，相较于正常固结土，超固结土能够更好地维持其初始稳定状态，表现出更好的受力与抗液化特性。

(3) 在饱和地基桩柱结构地震模拟数值算例中，受桩土间动力耦合作用的影响，桩柱结构近场地基的超孔隙水压要显著高于远场地基。同时通过不同土体本构模型下模拟计算结果的对比可以看出:下负荷面剑桥模型下得到的地震响应模拟结果介于理想弹性模型与剑桥模型之间; 相较于剑桥模型，在考虑超固结因素后，由于超固结土更好的抗液化以及力学特性，随着相应模型下地基承载能力的提高，结构的地震响应与受力情况均得到明显改善。