基于耗散能的海洋K0固结土统一修正剑桥模型

2021-02-02吴有平伏亮明周轩漾黄礼胜

长江科学院院报 2021年1期

吴有平，伏亮明，周轩漾，胡达，黄礼胜

（1.中国电建集团中南勘测设计研究院有限公司研发中心，长沙 410014；2.武汉大学土木建筑工程学院，武汉430072；3.湖南城市学院土木工程学院，湖南益阳 413000；4.广东省有色金属地质局，广州 510080）

1 研究背景

在海上风电桩基计算中采用较广的计算方法是p-y曲线法，该方法主要分为2类：一类是传统的半经验公式法（如API、DNV等）［1-2］；另一类是由土体的应力-应变曲线通过比例系数转换成p-y曲线，如Bouzid等［3］、Zhang和Andersen［4］推导的公式，研究表明，后者计算获得的结果较前者更优越。正如Bouzid等［3］认为p-y曲线由依赖于半经验法趋向有合理的理论基础的方法，从本构关系研究入手，获得的应力-应变曲线将更具严谨的理论基础。但现有的本构关系主要为数值计算提供服务，相应的计算所需要的参数较多，为方便获得应力-应变关系曲线，需要对现有模型进行改进。

自Roscoe等［5］于1963年提出被公认的弹塑性模型——剑桥黏土模型以来，不少模型在此基础上对其进行改进和修正。热力学理论因其物理意义清晰而被引入土体本构模型研究中，Houlsby［6］最早在土力学本构研究中采用热力学理论，后来Collins和Houlsby［7］、Collins和Hilder［8］、王秋生等［9］将该理论用于剑桥模型的改进中，获得了新的土力学模型，其理论更加严谨。为使建立的模型既适用于黏土又适用于砂土，国内外学者对此开展了不少研究。Lade［10］、Lade和Duncan［11］根据砂土的试验结果，提出了砂土的弹塑性模型。Lade和Musante［12］、Lade和Kim［13］通过修改屈服面函数将模型推广至黏性土，在1994年建立了适用于砂土和黏土的统一单屈服面模型。Crouch和Wolf［14］在1994年建立了统一三维临界状态塑性模型。Yao等［15-18］结合SMP强度准则、Lade强度准则及广义非线性强度理论提出了适用于砂土与黏土的统一硬化模型，并在下加载面的基础上推导了适用于超固结土体的统一硬化模型。此外，孔亮等［19］、李学丰等［20］引入统一硬化参量对Collins的热力学土体模型进行改进，从而将该模型扩展为既能反映砂土剪胀等变形特性又适用于黏土的统一热力学模型。Yu［21］引入应力状态系数，在剑桥模型基础上构建了一个既可应用于砂土，又可应用于黏土的本构模型，即CASM模型。Yao等［22］假设砂土CSL（Critical State Line）随洛德角而改变，扩展以前的统一硬化模型，并通过应力转换将模型扩展为三维模型。上述部分统一硬化模型采用塑性体应变作为硬化参量值得商榷，如Yu［21］、王秋生等［9］构建的模型。在这些统一硬化模型中皆需要选用较多计算参数来保证计算结果的精度，但过多的参数选取提高了模型的应用难度。为方便模型的选用，针对特定的地质条件开展了模型的简化研究。

海底地层经过长期地质作用后，在无海底滑坡、火山喷发等情况下，其地层处于K0固结状态，且为饱水状态。对于海洋工程中的桩周土体而言，在受波浪、洋流等荷载下，单次荷载作用历时短，土体孔隙水排出较慢，可近似为不排水条件，只有在持续的风荷载作用下，土体才处于排水固结过程中。因此，海洋工程中基础所处的土体通常为饱和的原状土。根据文献［23］对我国海洋底质研究成果及中南勘测设计研究院有限公司的海洋工程实践，我国沿海地层以黏土、粉土、细砂、中粗砂为主，且各地层交错出现，有别于欧洲海上风电基础的单一地层［24］。因而传统针对单一地层的p-y曲线法在我国风电桩基设计中有待改进。根据海洋工程中桩周土体特性，在修正剑桥模型的基础上，引入相变参数，推导出适用于剪胀性土体的屈服面函数。引入耗散功作为硬化参量，结合传统弹塑性理论，对Collins的热力学土体模型进行改进，提出了一种同时适用于砂土与黏土的新的统一修正剑桥模型。

2 K0固结条件下修正剑桥模型研究

2.1 修正剑桥模型的塑性势面

修正剑桥模型的塑性功表达式为［5］

式中：dWp为塑性功增量；p′、q分别为平均主应力和广义剪应力分别为塑性体应变与塑性剪应变；M为临界状态线斜率。

整理式（1）可得

求解式（3）可得

式中C为积分常数。式（4）为椭圆曲线方程，在修正剑桥模型特征点q/p′＝M处，塑性体应变在q＝0时

当q＝0时，p′＝p′c，p′c为土体前期固结压力；代入塑性势方程后可以解得

由于采用关联流动法则，所以屈服面函数与塑性势函数相等。在修正剑桥模型参数中引入剪胀指标，修正后的应力比为Mp。采用关联流动法则，则修正后的屈服面函数及塑性势函数为

式中p′c为前期土体固结压力，在p-q图中为椭圆与p轴的交点。

具有剪胀性土体的屈服面方程可以表示为

2.2 修正剑桥模型的不足与改进

当加载超过了原来的屈服面，土体产生新的塑性变形，如果屈服应力得到提高，这种特性称为硬化。用于计算从一个屈服面到另外一个屈服面塑性应变增量大小的内参量称为硬化参量。土体硬化的塑性应变可能与应力路径有关或者无关，但作为弹塑性模型的硬化参量必须是应力路径无关的。对于正常固结的黏性土，大量的试验已证实可以采用塑性体应变作为硬化参量。Nakai［25］通过对中等密度的Toyoura砂进行不同应力路径下的三轴试验，发现塑性体应变与应力路径有关，应力路径越长，其塑性应变则越大（如图1（a）所示）；显然作为硬化参量是不合适的。同时，Nakai［25］、Yao等［17］发现在不同应力路径下塑性功近似相等（如图1（b）所示）。ADF、ABEF、ABCF、ABF为三轴压缩路径，详见文献［25］。

图1 塑性体应变、单位体积塑性功W p与p的关系Fig.1 Relations of and W p against p

文献［25］中塑性功采用式（8）形式计算，即

式中：tij为剪应力张量为塑性应变张量；aij为主应力方向与空间滑动面余弦值的对称张量详见文献［25］；tN、tS分别为剪应力在空间滑动面法向分量与切线分量。

Collins和Hilder［8］认为一般假设塑性功是耗散的，而弹性功是可以恢复（存储）的，但实际很难满足。其推导的塑性功表达式为

图2为p′-q应力空间中实际应力屈服面图，若在耗散应力空间中，则屈服面的中心点位于原点处。图中CSL为临界状态线。真实应力与耗散应力主要区别是体应力p′相差一个迁移应力ρ′，在修正剑桥模型中，此迁移力为p′c/2。

图2 p′-q应力空间中实际应力屈服面［8］Fig.2 True stress yield surface in p′-q stress space［8］

在耗散应力空间中，应力耗散可以采用式（10）表达。

式中：π′为耗散压力；τ为耗散剪应力；πc为耗散应力空间中的正常固结压力。

将式（10）和式（11）消去应变增量，可以获得在耗散应力空间的屈服轨迹方程为

在各向同性压缩条件下引入相同的迁移应力，πc＝p′c/2，π′＝p′-p′c/2，屈服轨迹在真实应力空间的轨迹如式（13）或式（14）所示。

采用关联流动法则，则耗散函数的应力增量如式（15）所示。

式中δμ为流动法则中比例系数。

由式（14）和式（15）可以得塑性剪胀Δ的表达式为

代入式（10），结合式（11）、式（14）可得耗散功为

2.3 K0固结条件下修正剑桥模型屈服面函数

对于K0固结的土体，在耗散空间的旋转角度与真实应力空间的旋转角度大小相同，只是旋转点不同，如图3所示。图3中，DFL为破坏状态线，NCL为正常固结线。

图3 土体各向异性时的耗散应力空间Fig.3 Dissipative stress space for soil anisotropy

式中：A、B分别为等向压缩与剪切时描述耗散的无量纲函数；θn为倾斜椭圆最右端与原点连线即NCL的倾角为NCL的斜率，K0为静止土压力系数。

对式（18）求偏导，可以获得耗散应力分量，即

塑性应变增量采用式（20）表示，即

真实应力空间中应力轨迹如图4所示。

迁移应力分量可以表示为

图4 真实应力空间屈服轨迹Fig. 4 Yield locus in true stress space

式中ρ′、ζ分布为p′、q方向迁移应力分量。

因而屈服面函数可以表示为

对于砂土及应力硬化模型，体应变增量消失时的状态线称为临界状态线，即图3和图4中的CSL线。在耗散应力空间中，令由式（20）可得

式中θCSL为屈服面与CSL的交点与中心点连线的倾角。

假设屈服轨迹与CSL相交于（p′CSL，qCSL），则p′CSL，qCSL满足式（22）且同时满足式（24），即

令M0＝B/A＝tanφ0，联立式（22）和式（24），可以获得式（25），即

又M＝tanφCSL，A＝p′c/2，可得

代入式（22），可得屈服面函数为

剑桥模型沿p轴等向塑性体变硬化，对于初始各向异性土体，在日本广泛采用关口-太田模型沿初始固结线（K0线）不等向塑性体变硬化。试验数据表明，2个模型结果与试验值皆有出入，对此，孙德安等［26］提出了介于两者之间的修正。本研究基于该思想对屈服面的p轴引入系数β进行修正，此时旋转角tanθβ＝βtanθn，即相当于βηk0。为了便于理解，在公式中不区分tanθn与tanθβ，后续部分计算K0固结土体时直接引入系数β。M0中的B和A为倾斜后椭圆的竖向与水平轴长度，对于给出的是无倾斜的椭圆M值时，需要按M0＝M/cosθβ2进行换算。对于砂土，当η＞MP时，进入相变阶段。η为倾斜椭圆的斜率，MP为三轴试验中砂土体积由压缩突变为膨胀时的应力比。

2.4 K0固结条件下修正剑桥模型流动法则

流动法则是确定塑性应变增量方向的一条规定。如将主应力空间中同量塑性能的点连接起来形成塑性势面，用函数g表示。

塑性应变增量与应力之间存在如下关系，即

式中：dΛ为塑性应变增量系数；σij为应力张量。

采用p′、q应力表示时，式（28）可表示为

如果塑性势函数g与屈服面函数f相同，则称为关联流动法则，否则为非关联流动法则。

2.5 K0固结条件下修正剑桥模型硬化规律

根据Collins和Hilder［8］对迁移应力的定义，其表达式为

式中：Ψ2为自由能函数第2部分（第1部分为弹性能）为采用Hencky应变测量法获得的塑性应变张量。结合式（21）可得

式中ρ′v、ρ′s分别为正应力与剪应力方向对应的迁移应力。

整理式（31）得

考虑塑性部分的自由能函数Ψ2的微分可以表示为

塑性部分的自由能函数Ψ2为常量，因而dΨ2＝0，设结合式（32）可得

对式（34）进行积分，不考虑θn变化时，即不考虑旋转硬化，可得

式中C为常数。

根据应变空间的特征曲线，可得

式中：p′0为初始有效压力；λ、κ分别为正常固结线与弹性加载线的斜率；e0为初始孔隙比。

正常固结压力可由式（37）表示，即

屈服面函数（27）结合式（37）经过整理后可表示为

其中：

式中H为硬化参数。

参考文献［8］，以耗散功Wp作为硬化参量，结合式（18）、式（27），耗散功表达式为

对于K0固结的屈服面函数可以采用式（40）表达，即

p0′F(Wp)是pc′的函数，对式（40）两边微分，可得

将式（40）取微分，结合式（39）可得

加工硬化规律是指在给定应力增量下引起塑性应变增量变化的规律。假定塑性应变增量系数dΛ采用式（43）表示，即

式中：f为屈服条件函数，f(σij，H)＝0；A为硬化参数的函数。采用耗散功Wp来计算A的数值，即

式中α为由应力与应变计算功时的调整系数。

由式（43）可知

由流动法则可得

两边各乘ασij得

结合式（45）可得

f值的确定采用等向压缩路径，通过对应的应力应变计算获得f-Wp关系图。黄文熙［27］从常规三轴压缩试验的（σ1-σ3）-ε1曲线与ε1-εv曲线，可以计算2根曲线上某一应力状态的f值与Wp值，绘制1根f-Wp曲线，这根曲线相当于σ3等于某一定值。用不同的σ3做试验可得1组曲线。将这组曲线规格化，可得f＝F（Wp）函数，求得相当于某一f值的F′＝∂F/∂Wp。本文中F（Wp）选用a（Wp）b的形式表示，a、b为系数。

对于弹性部分的计算，可以由e-lnp关系中压缩和回弹曲线推算出，即

图5 计算流程Fig.5 Flow chart of computation

3 计算流程

计算流程如图5所示。计算主要过程如下：

（1）通过基本试验与三轴试验，获得剑桥模型参数λ、κ、M、θ值。

（2）根据三轴试验获得σ、ε值，代入式（39）计算Wp值；令式（40）中F（Wp）＝0，计算f值。

（3）在同一σ3下获得的不同f和Wp值，绘制f-Wp关系曲线。

（4）将不同σ3下f-Wp进行归一化，获得f-Wp关系式，即f＝F（Wp），同时获得系数a、b的值。结合式（44），获得α值。

（5）由f＝F（Wp）计算∂f/∂Wp。

（6）采用关联流动法则，令g＝f，由式（48）与式（43）计算dΛ。

（7）将荷载进行分级。

（8）由式（27）计算屈服状态。

（9）当计算加载面进入屈服时，由式（29）计算塑性应变，由式（49）计算弹性应变。

（10）当计算加载面未进入屈服时，由式（49）计算弹性应变。

（11）计算获得σ、ε、W值。

4 模型验证

为验证本研究提出的计算方法可行性，分别选取砂土和黏土进行计算，具体结果如下。

4.1 承德中密砂

承德中密砂是一种天然均匀细砂，参考有关该砂的材料属性及计算参数，结合试验得到的应力-应变关系曲线，获得的计算参数见表1。计算结果如图6所示。

表1 承德中密砂计算参数Table 1 Computational parameters of moderately-dense sand in Chengde

从图6可以发现：①主应变随剪应力的增加而近似指数函数增加，计算值与试验值吻合较好，且随围压的增大，吻合程度增高；②体应变在3种试验条件下都存在剪胀，低围压时，计算值与试验值吻合较好。同时，从图6还可看出，主应变在材料软化阶段的数值不能反映出来，这是本模型的一个缺陷。对于不需要计算破坏阶段的工程（如风电桩基）而言，该模型可以满足工程计算需要。

图6 计算与试验剪切应力-应变曲线Fig. 6 Calculated and test shear stress-strain curves

4. 2 Bothkennar 黏土

表2 Bothkennar 黏土计算参数Table 2 Parameters of Bothkennar clay for calculation

从图7可以发现：①体积应力-应变本文计算值在体应变＜5.2%时，与试验值吻合较好；②剪切应力-应变在试验条件下计算值与试验值除个别点外，都能较好吻合。而参考文献［29］的体积应力-应变在应变＜1.6%时，能与试验值较好吻合；剪切应力-应变关系与试验值存在偏差。本文计算所需参数为10，少于文献［29］所需数目。因而，本研究提出的本构模型在应用于K0固结的黏土计算中具有较高精度，且计算所需参数少。

综合前述砂土与黏土的案例计算分析可以认为，本研究提出的统一修正剑桥模型具有计算精度高、选用计算参数少且物理意义明确的优点，且能同时适用于黏土与砂土本构关系计算。新本构模型的参数减少有利于工程应用，计算精度的提高能减少安全冗余，同时可作为有限元的本构模型，为提高数值模拟结果的精度提供基础。近年来海上风电的崛起，新提出的土体本构模型为改进现有基于单一地层小直径桩基计算的p-y曲线法提供理论基础，以满足桩径日益增大且地质条件复杂条件下的桩基计算需要。