写在前面
在岩土工程的数值模拟分析中，初始地应力的平衡是必须予以重视的问题。
对复杂尤其是存在复杂约束条件及相互作用的模型中，仅使用传统的地应力平衡方法很难满足平衡要求。
因此，本文提出了一种分步平衡法，将ABAQUS生死单元技术与自动平衡法结合起来，来解决较复杂岩土性质和接触条件下的地应力平衡问题。
基于生死单元技术分步平衡法的应用
以吸力式桶型基础为例来验证此方法的地应力平衡效果。考虑有限元模型的对称性及数值模拟的计算效率，取地基土体和桶体的一半进行建模分析。吸力桶基础外径D0为3m，壁厚ts为0.1m，长度L为6m。为了避免边界效应对数值模拟结果产生影响，地基土半径D取为7D0，深度H取3.33L。
图1为划分的有限元网格。
在地基底部边界上，约束三个方向的自由度，在侧面边界上对水平方向两个自由度进行约束，而在对称面上则对法向自由度进行约束。
土体采用符合Mohr-Coulomb屈服准则的理想弹塑性本构模型。目前基于ABAQUS有限元软件采用Mohr-Coulomb本构模型对砂土地基进行数值模拟计算分析时，地基土的变形模量多采用恒定值。而实际上，砂土的变形模量随着围压水平而改变，并不是一个恒定值，这种压硬性可以采用Janbu提出的幂函数公式来描述，如公式(1)所示。
式中：Es是砂土地基的侧限压缩模量；大气压强σat取100 kN/m2；κ决定了基准应力状态下的土体刚度；λ反映了土体刚度的应力依赖性程度；σm为平均主应力。
由于ABAQUS中的Mohr-Coulmb弹塑性模型默认采用变形模量，所以根据以下土力学公式将侧限压缩模量Es转换为变形模量E，如式(2)所示。
式中：E为变形模量；ν为土体泊松比。参考相关规范以及Achmus等的建议，本次数值模拟采用κ=600，λ=0.55，具体模型参数见表1所示。
通过Fortran语言编写子程序进行二次开发，实现土体模量的应力相关性。
桶土界面接触行为采用“摩擦接触对”算法模拟，法向上采用硬接触，切向上采用Coulomb定律描述其力学性质，摩擦系数f=tan(0.6φ)，其中φ为砂土内摩擦角。
对此，分别采用自动平衡法和本文提出的分步平衡法来进行地应力平衡。使用地应力自动平衡方法，是在一个地应力分析步完成地应力平衡，吸力桶与地基土体直接在地应力分析步中设置接触条件及相互作用;而分步平衡法则需两个地应力分析步来分步完成初始地应力平衡，吸力桶基础与土体的接触及相互作用只有在吸力桶结构激活时才起作用，能够更加真实地反映实际工况。
平衡效果对比
分别对砂土地基定模量和变模量两种情况下的地应力平衡效果进行评价。评价时，主要通过竖向应力大小和竖向位移归零程度两方面进行衡量。
定模量模拟效果对比
对于定模量情况，砂土地基的变形模量E取值为60MPa。分别采用自动平衡法和本文提出的分步平衡法来进行地应力平衡，在分析步中施加体力加载。
图2为竖向应力S33分布的比较，图中应力单位为Pa。
通过图2可以看出，竖向应力随着土体深度的增加而增大。竖直方向的自重应力S33要满足S33=γ·H。在地基底面，竖向自重应力的大小理论上约等于220kN/m2，分步平衡法得到的结果比其约低18%，而自动平衡法的误差则达到了28%。从自重应力分布来看，分步平衡法得到的结果更接近于理论值，随后进行的外力加载分析结果也会更接近于真实情况。
在自动平衡法中，为了能够收敛，吸力桶的重度一般假设为土体的重度，这对于固定式海上风机桶基工程的受力分析可能影响不大，但是对于桶形基础用于承受拉拔荷载时会造成一定误差。而对于分步平衡法，在第二步激活桶体单元时，可以施加真实的钢材料的重度，从而更加真实地反映实际工况。
另外，从图2还可以看出，分步平衡法所得到的竖向应力分布更加均匀，而自动平衡法得到的竖向应力在吸力桶周围产生了明显的局部错位情况。
初始地应力平衡效果的优劣程度还取决于其竖向位移U3归零是否能够满足要求。一般岩土工程领域对地应力分析步的竖向位移U3归零的要求是小于10-4m量级。
图3为两种方法在定模量模拟情况下的竖向位移的对比，单位为m。
从图3中可以看出，无论是采用哪种方法对定模量进行数值模拟，其竖向位移归零均能满足量级要求。但是就位移分布规律来看，自动平衡法模拟结果中，桶体内部土体产生了类似条纹状的彩带，出现了不符合规律的位移变化，效果不是特别理想。而分步平衡法模拟出来的位移变化情况，桶基础中土体因为桶的贯入而发生沉降，桶盖周围会因为贯入而产生部分隆起，更加符合模型试验和原位测试观察到的实际现象
变模量模拟效果对比
采用相同的模型尺寸及子程序对两种方法的地应力平衡效果进行对比分析。在ABAQUS应用子程序时，土体变形模量是在给定的场变量FV1区间内进行线性插值赋予具体数值的。图4为子程序在两种方法中调用后的模拟情况，图中模量单位为Pa。
根据式(1)和式(2)可以估计地基土体底面的变形模量理论值。由于砂土地基一般满足K0固结条件，其侧压力系数约等于：
因此地基底面的平均应力可以计算为:
将式(4)代入式(1)可得地基底面的侧限压缩模量Es=68607kPa，进而根据式(2)得变形模量E=50965kPa≈5.1×107Pa。两种方法得到的地基底面变形模量值相差不大，比理论值约低6%。
但对于地基顶面，此处的平均应力σm=0，根据式(1)和式(2)其变形模量理论值应该约为0。从图4看出，分步平衡法准确反映了地基浅层土体的低模量水平，而自动平衡法得到结果约为3961kPa，错误高估了地基浅层土体的变形模量，导致地基土体的变形模量分布规律与式(1)的幂函数关系相差较大，对后续桶体在拉拔荷载作用下的位移预测带来很大影响。
图5为变模量情况下竖向应力S33分布的比较，图中应力单位为Pa。与定模量模拟情况类似，从自重应力分布来看，分步平衡法得到的结果更接近于理论值，随后进行的外力加载分析结果也会更接近于实际情况。
图6为两种方法在变模量模拟情况下的竖向位移的对比，单位为m。
通过比较可以看到，采用自动平衡法得到的竖向位移量级已经达到10-4m，这对于固定式海上风机桶基工程的位移分析影响不大，但对于桶形基础承受拉拔荷载时的位移会产生一定误差。而且，自动平衡法所得到的浅层土体竖向位移场产生了较大的波浪形趋势，位移分布极不均匀，这与实际情况有较大的差异。
采用分步平衡法的地应力平衡效果明显提高，其平衡后的竖向位移在10-5次方量级，比自动平衡法结果低一个量级，因此更加满足地应力平衡的要求。与对应的定模量情况模拟得到的竖向位移分布规律类似，位移场仅在吸力桶基础临近区域受到影响，桶内土体因贯入的影响产生沉降，桶盖周围土体受到挤压而发生少量隆起，与实际情况吻合度较高。
通过上述分析可以看到，分析过程中存在结构与土体相互作用或者岩土体材料性质比较复杂时，仅采用自动平衡法，其平衡效果往往不理想或不能满足平衡要求，甚至有时会出现收敛困难而导致数值计算任务中断。本文通过吸力桶工程实例的地应力平衡效果的对比，可以看出分步平衡法可有效改善地应力平衡效果，得到符合实际情况的数值模拟结果，对类似数值模拟具有一定的参考价值。
