基于ABAQUS生死单元技术的平衡地应力方法研究

肖国丰，马思琪，刘 洋，张兴瑞，范庆来

(鲁东大学 海洋岩土工程研究所， 山东 烟台 264025)

随着科学技术的不断进步以及电子信息产业的发展，计算机技术也取得了前所未有的革新与突破，数值模拟也随之得以快速发展、广泛应用。ABAQUS作为一套功能强大的工程模拟有限元软件，被广泛应用于岩土工程的实践与研究，为岩土工程领域线性以及非线性问题的研究做出了突出贡献。同时ABAQUS对岩土材料的模拟具有很强的适用性，经过长时间的版本升级其功能也愈加完善。在实际岩土工程中，土体经过长年累月的沉淀，已形成相对稳定的地质状态。而使用有限元软件对岩土工程问题进行求解的时候，为了使有限元软件尽可能真实可靠的模拟实际工况，需在加载之前使其地应力达到一种平衡状态，即土体内部应力与施加的重力荷载相平衡。只有初始地应力满足平衡要求，后续工况模拟才能得到符合实际工况的解。

在岩土工程的数值模拟分析中，初始地应力的平衡是必须予以重视的问题。针对ABAQUS对初始地应力平衡提出的多种方法，目前已有一些成果对其进行分析探讨。于新等[1]应用自动平衡法与ODB导入法分别进行计算，认为对简易模型来说两种方法差别不大；杨金尤等[2]则对带帽刚性桩复合地基进行了计算，其中定义关键字法受条件约束较大，而ODB导入法对复杂模型存在需反复导入的弊端；而对简易模型来说，导入ODB方法同样可以明显改善地应力平衡效果[3]；代汝林等[4]同时对不同复杂程度的模型分别采用多种方法进行计算，对于地质条件复杂的岩土体采用自动平衡法和应力提取法效果明显优于其他方法；郭亚然等[5]则提出一种粘弹性边界模型来解决动力响应分析中地应力平衡问题。

以上研究多对简易模型进行了分析探讨，而对复杂尤其是存在复杂约束条件及相互作用的模型中，仅使用传统的地应力平衡方法很难满足平衡要求。因此，本文提出了一种分步平衡法，将ABAQUS生死单元技术与自动平衡法结合起来，来解决较复杂岩土性质和接触条件下的地应力平衡问题。

1 基于生死单元技术的分步平衡法的提出

进行初始地应力平衡的目的是忽略工程结构物的施工贯入对地基土体带来的影响，直接对已经安装就位的工程结构物的受荷情况进行分析，而此时就要求其工程结构物在地基土中已经处于一个只存在应力而不存在位移的相对稳定状态[6]。在简单的岩土工程问题中，其初始地应力平衡有着多种方法，其中自动地应力平衡方法作为“一步平衡法”的代表，由于其操作简单得到了广泛应用。但如果实际模拟工况较为复杂，工程结构物与土体接触较多且存在用户定义场时，仅使用自动地应力平衡方法有时很难满足地应力平衡要求，甚至会遇到收敛困难的问题[7-8]。

本文基于ABAQUS中的生死单元功能结合自动平衡法对地基初始地应力场进行平衡。该方法主要分为两步，因此本文称之为“分步平衡法”。在分别建立地基土体和结构物部件后，第一步“杀死”结构物单元，仅对地基土体进行地应力平衡，模拟工程结构物未贯入状态；第二步“杀死”结构物占据区域内的土体单元，同时激活装配好的工程结构物部件，并设置满足实际工况的接触条件，进行结构与地基相互作用后的地应力平衡，模拟工程结构物安装就位后的稳定状态。分步平衡法相较于自动平衡法来说，能够有效缓解收敛性困难，同时地应力平衡与位移归零精度更高。

2 分步平衡法的应用

以吸力式桶型基础为例来验证此方法的地应力平衡效果。考虑有限元模型的对称性及数值模拟的计算效率，取地基土体和桶体的一半进行建模分析。吸力桶基础外径D0为3 m，壁厚ts为0.1 m，长度L为6 m。为了避免边界效应对数值模拟结果产生影响，地基土半径D取为7D0，深度H取3.33L。图1为划分的有限元网格。在地基底部边界上，约束三个方向的自由度，在侧面边界上对水平方向两个自由度进行约束，而在对称面上则对法向自由度进行约束。

图1 有限元网格

土体采用符合Mohr-Coulomb屈服准则的理想弹塑性本构模型。目前基于ABAQUS有限元软件采用Mohr-Coulomb本构模型对砂土地基进行数值模拟计算分析时，地基土的变形模量多采用恒定值[9-10]。而实际上，砂土的变形模量随着围压水平而改变，并不是一个恒定值，这种压硬性可以采用Janbu提出的幂函数公式来描述[11]，如公式(1)所示。

(1)

式中：Es是砂土地基的侧限压缩模量;大气压强σat取100 kN/m2;κ决定了基准应力状态下的土体刚度;λ反映了土体刚度的应力依赖性程度;σm为平均主应力。

由于ABAQUS中的Mohr-Coulmb弹塑性模型默认采用变形模量，所以根据以下土力学公式将侧限压缩模量Es转换为变形模量E，如式(2)所示。

(2)

式中：E为变形模量；ν为土体泊松比。参考相关规范以及Achmus等的建议[11]，本次数值模拟采用κ=600，λ=0.55，具体模型参数见表1所示。通过Fortran语言编写子程序进行二次开发，实现土体模量的应力相关性。

表1 材料参数

桶土界面接触行为采用“摩擦接触对”算法模拟，法向上采用硬接触，切向上采用Coulomb定律描述其力学性质，摩擦系数f=tan(0.6φ)，其中φ为砂土内摩擦角。

对此，分别采用自动平衡法和本文提出的分步平衡法来进行地应力平衡。使用地应力自动平衡方法，是在一个地应力分析步完成地应力平衡，吸力桶与地基土体直接在地应力分析步中设置接触条件及相互作用；而分步平衡法则需两个地应力分析步来分步完成初始地应力平衡，吸力桶基础与土体的接触及相互作用只有在吸力桶结构激活时才起作用，能够更加真实地反映实际工况。

3 平衡效果对比

分别对砂土地基定模量和变模量两种情况下的地应力平衡效果进行评价。评价时，主要通过竖向应力大小和竖向位移归零程度两方面进行衡量。

3.1 定模量模拟效果对比

对于定模量情况，砂土地基的变形模量E取值为60 MPa。分别采用自动平衡法和本文提出的分步平衡法来进行地应力平衡，在分析步中施加体力加载。

图2为竖向应力S33分布的比较，图中应力单位为Pa。通过图2可以看出，竖向应力随着土体深度的增加而增大。竖直方向的自重应力S33要满足S33=γ·H。在地基底面，竖向自重应力的大小理论上约等于220 kN/m2，分步平衡法得到的结果比其约低18%，而自动平衡法的误差则达到了28%。从自重应力分布来看，分步平衡法得到的结果更接近于理论值，随后进行的外力加载分析结果也会更接近于真实情况。

图2 竖向自重应力分布的对比

另外，在自动平衡法中，为了能够收敛，吸力桶的重度一般假设为土体的重度，这对于固定式海上风机桶基工程的受力分析可能影响不大，但是对于桶形基础用于承受拉拔荷载时会造成一定误差。而对于分步平衡法，在第二步激活桶体单元时，可以施加真实的钢材料的重度，从而更加真实地反映实际工况。

另外，从图2还可以看出，分步平衡法所得到的竖向应力分布更加均匀，而自动平衡法得到的竖向应力在吸力桶周围产生了明显的局部错位情况。

初始地应力平衡效果的优劣程度还取决于其竖向位移U3归零是否能够满足要求。一般岩土工程领域对地应力分析步的竖向位移U3归零的要求是小于10-4m量级。图3为两种方法在定模量模拟情况下的竖向位移的对比，单位为m。

图3 竖向位移的对比

从图3中可以看出，无论是采用哪种方法对定模量进行数值模拟，其竖向位移归零均能满足量级要求。但是就位移分布规律来看，自动平衡法模拟结果中，桶体内部土体产生了类似条纹状的彩带，出现了不符合规律的位移变化，效果不是特别理想。而分步平衡法模拟出来的位移变化情况，桶基础中土体因为桶的贯入而发生沉降，桶盖周围会因为贯入而产生部分隆起，更加符合模型试验和原位测试观察到的实际现象[12-13]。

3.2 变模量模拟效果对比

采用相同的模型尺寸及子程序对两种方法的地应力平衡效果进行对比分析。在ABAQUS应用子程序时，土体变形模量是在给定的场变量FV1区间内进行线性插值赋予具体数值的。图4为子程序在两种方法中调用后的模拟情况，图中模量单位为Pa。

图4 土体模量分布的对比K0=1-sinφ=0.37

根据式(1)和式(2)可以估计地基土体底面的变形模量理论值。由于砂土地基一般满足K0固结条件，其侧压力系数约等于：

(3)

因此地基底面的平均应力可以计算为：

(4)

将式(4)代入式(1)可得地基底面的侧限压缩模量Es=68 607 kPa，进而根据式(2)得变形模量E=50 965 kPa≈5.1×107Pa。两种方法得到的地基底面变形模量值相差不大，比理论值约低6%。但对于地基顶面，此处的平均应力σm=0，根据式(1)和式(2)其变形模量理论值应该约为0。从图4看出，分步平衡法准确反映了地基浅层土体的低模量水平，而自动平衡法得到结果约为3 961 kPa，错误高估了地基浅层土体的变形模量，导致地基土体的变形模量分布规律与式(1)的幂函数关系相差较大，对后续桶体在拉拔荷载作用下的位移预测带来很大影响。

图5为变模量情况下竖向应力S33分布的比较，图中应力单位为Pa。与定模量模拟情况类似，从自重应力分布来看，分步平衡法得到的结果更接近于理论值，随后进行的外力加载分析结果也会更接近于实际情况。

图5 竖向应力分布的对比

图6为两种方法在变模量模拟情况下的竖向位移的对比，单位为m。通过比较可以看到，采用自动平衡法得到的竖向位移量级已经达到10-4m，这对于固定式海上风机桶基工程的位移分析影响不大，但对于桶形基础承受拉拔荷载时[14]的位移会产生一定误差。而且，自动平衡法所得到的浅层土体竖向位移场产生了较大的波浪形趋势，位移分布极不均匀，这与实际情况有较大的差异。

图6 竖向位移的对比

采用分步平衡法的地应力平衡效果明显提高，其平衡后的竖向位移在10-5次方量级，比自动平衡法结果低一个量级，因此更加满足地应力平衡的要求。与对应的定模量情况模拟得到的竖向位移分布规律类似，位移场仅在吸力桶基础临近区域受到影响，桶内土体因贯入的影响产生沉降，桶盖周围土体受到挤压而发生少量隆起，与实际情况吻合度较高[12-13]。

通过上述分析可以看到，分析过程中存在结构与土体相互作用或者岩土体材料性质比较复杂时，仅采用自动平衡法，其平衡效果往往不理想或不能满足平衡要求，甚至有时会出现收敛困难而导致数值计算任务中断。本文通过吸力桶工程实例的地应力平衡效果的对比，可以看出分步平衡法可有效改善地应力平衡效果，得到符合实际情况的数值模拟结果，对类似数值模拟具有一定的参考价值。

4 结 论

在岩土工程数值计算中，对于复杂工况，要获得一个只存在地应力而不存在位移的平衡状态，仅采用自动平衡法等一步平衡法往往不理想，甚至有时会出现收敛困难。本文基于ABAQUS软件中的生死单元功能，提出了一种分步平衡法，通过分步创建桶土之间的接触以更好地模拟地基在桶土相互作用过程中的实际响应。通过分析，得到如下结论：

(1) 对于定模量情况下的两种方法对比，分步平衡法得到的竖向应力场更接近于理论解，竖向位移量级尽管与自动平衡法相差不大，但是得到的竖向位移分布规律更接近于已有的模型试验结果。

(2) 对于变模量情况下的两种方法对比，分步平衡法不仅得到的竖向应力场更接近于理论解，而且竖向位移量级也比自动平衡法低1个量级。

(3) 对于压硬性砂土地基中桶土相互作用问题，分步平衡法能够合理反映地基浅层土体的低模量特征，而自动平衡法得到的结果则错误高估了地基浅层土体的变形模量，导致地基土体的变形模量分布规律与Janbu公式相差较大。