企业官网建设流程全解析

做岩土、隧道、边坡数值模拟的朋友，早晚都会撞上“初始地应力场”这个词。Abaqus里如果不把这个场设对，后面算隧道开挖、基坑卸荷、边坡稳定性，第一步就可能出现几十厘米甚至几米的虚假位移，整个模型约等于作废。我见过不少人在论坛问“为什么我开挖出来的位移这么大”“地应力平衡永远不收敛”，十有八九就是初始应力没给对。这篇文章我就把自己用Abaqus设置初始地应力场的几种常用方法、完整操作步骤、以及调试经验一次性讲清楚，新手看完能照着做，老手也可以对照检查自己的设置流程。

1. 为什么要设置初始地应力场：不设置会出什么乱子

1.1 初始地应力场是什么：岩体本来就是“带应力上班”的

很多没接触过岩土方向的人第一次听到“初始地应力场”会觉得抽象，其实道理很简单：地面以下几十米甚至几百米的岩土体，并不是一块刚从Abaqus里新建的、无应力无变形的“白纸材料”。它在上覆岩土自重、地下水、历史构造运动等因素下，早就处于一个受力平衡的应力状态。这个“与生俱来”的应力状态，就是初始地应力场。它由竖向自重应力和水平构造应力组成，通常用竖向应力σv = ρgh、水平应力σh = K0σv来近似。K0是侧压力系数，对于正常固结土可以用K0 = ν/(1-ν)估算，岩石工程里更多是根据实测或经验取0.5到1.5之间的一个值。

Abaqus默认模型是从零应力状态开始的。如果不人为设置初始地应力，直接施加重力，相当于把原本已经稳定平衡了几千上万年的岩体，从“无应力”状态强行压到“有应力”状态。这个加载过程会产生一个明显的压缩变形，位移量级往往远大于你关心的工程位移。也就是说，你第一步算出来的“沉降”“变形”根本不是工程要看的增量变形，而是一个由“人为把重力加上去”造成的虚假压缩量，后续所有开挖、支护、卸载的分析结果都会被这个虚假位移污染。

1.2 不设置初始应力，隧道开挖算出来全是“假位移”

举个具体例子。一个埋深50米的隧道，上覆岩层平均密度按2000 kg/m³算，重力加速度取9.8，那么隧道所在位置的竖向应力大约是1MPa。如果岩体的弹性模量只有100MPa量级，应变大约就是1%。从地表到隧道深度这50米厚的岩体，在自重作用下理论上要被压缩掉数十厘米。这个位移是地壳在漫长地质年代里早已完成的变形，不是工程开挖引起的位移。你要是把它留在计算结果里，后面看隧道拱顶沉降云图时，会发现整个模型都在往下“沉”，数值大得离谱，但这不是隧道开挖造成的，而是初始压缩变形没有扣除。

地应力平衡要做的事，就是“先让这个初始应力场在重力和其他外载下自洽平衡，并把由此产生的位移尽可能清零”。平衡做完之后，模型处于“有应力、无位移”的状态。后续不管是开挖、加支护还是施加超载，新得到的位移才是真正的增量位移，才能拿去和现场监测数据对比。这个逻辑是整个岩土数值模拟的地基，地基不正，楼盖得再漂亮也没用。

1.3 哪些工程必须做地应力平衡，哪些可以偷懒

只要你模拟的对象是岩土体，而且关心的是“开挖/加载之后相对初始状态的变形和应力变化”，那就必须做初始地应力平衡。典型场景包括：深埋隧道与地下洞室开挖、边坡稳定性分析、基坑开挖、桩基与地基沉降、矿山采动、油气井井壁稳定等。这些场景里，初始应力不仅影响位移，还会直接影响破坏判据，比如岩体是否进入塑性、节理是否张开，都和围压水平密切相关。初始应力场给错了，破坏模式都会变。

也有可以简化的情况。比如你只做地表浅层的一个小型填方工程，材料强度很低，自重应力影响本来就小；又比如你根本不关注岩土体本身的初始位移，只关心结构构件在外部荷载下的内力响应，那这类问题可以只把重力当作普通荷载加载，不必严格做地应力平衡。但我个人的建议是，凡是模型里出现了“岩土体+重力+开挖或加载”这三个要素，就老老实实把地应力平衡写上。很多期刊审稿人和工程评审专家对这个步骤有明确要求，你模型里第一分析步不是Geostatic，或者初始位移没清零，很容易被人一句话打回来。

2. 四种主流设置方法怎么选：自动平衡、SIGINI、ODB导入、分步法

2.1 方法一：GEOSTATIC自动平衡法，适合自重应力场

Abaqus/Standard里专门提供了Geostatic分析步，配合关键字*Initial Conditions, Type=Stress, Geostatic，可以在一个分析步内自动完成初始应力场的平衡。它的基本逻辑是：你告诉Abaqus“这个区域的初始应力随深度按某个梯度分布”，Abaqus在Geostatic分析步里把重力加上去，反复迭代修正应力，使模型在重力下达到平衡，并把由初始应力引起的位移收敛到极小。

这个方法的优点是非常省事，尤其适合水平成层、地表水平、边界规则的模型。你只需要在CAE里给一个Geostatic分析步，然后在Predefined Field里定义初始应力，或者直接改关键字输入几行数据就行。缺点是它假设初始应力是以水平分层为基本规律的，对于起伏地表、强烈构造应力、复杂地形条件下，纯靠这种方式给出来的应力场不一定符合实际情况。而且Geostatic自动平衡对网格质量和边界约束要求比较高，有时候会不收敛。

2.2 方法二：SIGINI用户子程序，适合复杂构造应力场

如果你要模拟的初始应力场不是简单随深度线性变化，而是随坐标有更复杂的关系，比如考虑了褶皱、断层、水平构造应力非均匀分布，或者你只是想写一个自定义的侧压力系数表达式，那最好用SIGINI子程序。SIGINI是Abaqus专门用于定义初始应力场的用户子程序，Abaqus在计算开始前会调用它，在每个积分点上给你当前单元的坐标和积分点信息，你再把应力分量填进SIGMA数组就行。

用SIGINI的好处是灵活，你可以在程序里写任意函数，读取外部数据文件，甚至按照不同材料区域做不同处理。缺点是，First你得会一点Fortran或者Python风格的程序思维，Second调试要比CAE里点几下鼠标麻烦一些。实际上SIGINI的Fortran模板很固定，把公式填对，编译通过，再用一个小模型验证结果，后面就能稳定复用了。我在3.2节会直接给一个可以套用的模板。

2.3 方法三：ODB/文件导入法，适合已有模型或复杂初始场

有时候你不想手动写公式，只想“用Abaqus算出来的应力场作为下一次分析的初始应力场”。最典型的做法是：先建一个模型，只施加重力和边界条件，跑一遍静力分析得到应力场，确定这个应力场是自己满意的、符合实测规律的，然后把这份应力场作为初始条件导入到正式计算模型里。

在Abaqus里可以通过两种方式实现。一种是用Initial Conditions, Type=Stress, File=job.odb，直接把之前分析得到的ODB文件作为初始应力来源；另一种是先把应力分量提取出来写成数据文件，再用Initial Conditions, Type=Stress, Input=xxx.dat读入。这种方式尤其适合“先算一个小子模型得应力，再映射到大模型”“地质体很复杂，用数值方法先算地应力场再做工程分析”这类工作流。操作上稍微繁琐一点，但精度和可控性都不错。

2.4 方法四：分步重力加载法，做不了子程序时的备用方案

如果你既不想写SIGINI，又担心GEOSTATIC自动平衡不收敛，还有一招“土办法”：先用一个Static, General分析步把重力加载上去，得到一个包含自重应力且包含自重位移的结果；然后在下一步分析前，把位移场清零（不是把应力清零），这样相当于人为抹掉了自重引起的位移，保留自重应力。实现起来可以通过重启动或者*Restart，也可以在后处理里把位移场导出再减去初始位移。

这个方法思路直白，很多老工程师会用它做初步试探。它的缺点是，位移清零不是一个严格的力学操作，如果模型里有塑性、接触等非线性因素，直接清零位移可能会破坏应力-应变关系的一致性，导致后续结果出现不协调。因此我一般把它当作备选方案，或者只用来快速验证整体量级，正式的科研和工程分析还是优先用前面三种。

2.5 选型对比表

3. 新手必看：CAE里设置初始地应力场的完整实操步骤

3.1 用GEOSTATIC自动平衡：从建模到关键字修改的完整流程

先说最常用的GEOSTATIC自动平衡法。假设你要做一个水平地层的隧道开挖模型，地表水平，模型范围200m×100m，隧道埋深30m。整体流程是：建几何、赋材料、装配、设分析步、加荷载和边界、定义初始应力、提交计算。

材料参数里必须包含密度和弹性模量、泊松比。密度是地应力计算的第一要素，没有密度，重力就无从谈起。分析步方面，第一步必须设置为Geostatic，而不是默认的Static, General。在Abaqus/CAE里，Step模块下选择Create Step，Procedure type选General，然后找到Geostatic，点开之后一般保持默认设置即可。求解过程中允许迭代修正初始应力，所以建议把增量步数设为一个较大的值，比如100，防止第一次迭代就报错。

边界条件方面，推荐的做法是：模型底部约束竖向位移，左右两侧约束水平位移，前后两面（如果是二维模型就是平面应变约束）约束对应自由度。这样模型在重力作用下不会整体刚体移动，又能自由产生侧向变形。如果模型很大，你也可以用“底部竖向约束+两侧法向约束”的常规岩土约束组合。

荷载方面，在Load模块里创建重力荷载，施加重力加速度，方向沿Y轴负向，大小9.8。注意Abaqus里体积力的单位取决于你用的单位系统，如果用国际单位m·kg·s，重力加速度就是9.8；如果用mm·t·s单位制，重力加速度要写成9800。单位不一致是地应力平衡不收敛的第一大原因，务必先确认。

关键一步是设置初始应力。在CAE中可以通过Load模块的Predefined Field创建，也可以直接修改inp文件。我更推荐在inp文件里增加关键字，因为看得清楚，也方便后期批量调整。在*Step, name=Geostatic之前插入：

*Initial Conditions, type=stress, geostatic Eall, 0., 1000000., 0., -100., 0.65, 0.65

这行的含义是：单元集Eall，在深度坐标y=0处的竖向应力为0，在y=-100处的竖向应力为1MPa，模型顶面坐标是0，底面坐标是-100，水平侧压系数K0在面内为0.65，面外也为0.65。Abaqus会根据这两个深度点的竖向应力，按线性关系插值出整个模型每一点的竖向应力，再乘以K0得到水平应力。这里有个细节容易搞错：应力值必须带正负号，Abaqus默认压应力为负，但Geostatic这种输入格式它内部会自动按土压力习惯处理。更规范的做法是参考手册里的符号规定，最好先在简单模型上试一次，确保应力的正负方向符合预期。

提交计算后，打开ODB看第一步的位移云图。如果初始地应力设置正确，位移量级应该非常小，理想情况下达到10⁻⁴m以下，很多模型甚至能到10⁻⁶m。如果位移云图整体是红彤彤向下沉的，说明初始应力与重力不匹配，需要检查K0、密度、边界条件和单位。

3.2 用SIGINI子程序：一个可直接套用的Fortran模板

SIGINI用起来其实不难，它的核心逻辑是：Abaqus在每个积分点开始计算前调用你写的子程序，你根据传入的坐标COORDS，把该点的初始应力分量赋值给SIGMA数组。先记住Fortran模板：

SUBROUTINE SIGINI(SIGMA,COORDS,NTENS,NCRDS,NOEL,NPT, * LAYER,KSPT,LREBAR,NAMES) INCLUDE 'ABA_PARAM.INC' DIMENSION SIGMA(NTENS), COORDS(NCRDS) CHARACTER*80 NAMES(2) REAL rho, g, depth, K0 rho = 2000.0 g = 9.8 depth = -COORDS(2) K0 = 0.65 SIGMA(1) = -rho*g*depth SIGMA(2) = -K0*rho*g*depth SIGMA(3) = -K0*rho*g*depth SIGMA(4) = 0.0 RETURN END

这段代码默认你的重力方向是Y负向，所以COORDS(2)是Y坐标，depth取负号后变成正值深度。SIGMA(1)是Y方向的正应力，也就是竖向应力；SIGMA(2)和SIGMA(3)是两个水平方向的正应力；SIGMA(4)是剪切分量。对于平面应变模型，NTENS=3或4，需要根据实际的应力分量顺序调整。

写完子程序后，需要在模型关键字里声明使用SAF。在*Initial Conditions里面把type改成user：

*Initial Conditions, type=stress, user

然后在Job模块提交任务时，在Edit Job的General选项卡里，把Fortran子程序文件添加进去，或者用命令行提交：

abaqus job=jobname user=sigini.for

调试SIGINI时有个很实用的技巧：在子程序里临时加一段文件输出代码，把COORDS和SIGMA的值打印到一个txt文件里。这样提交一个小模型后，直接打开txt看各点的应力是否正确。不要一上来就跑大模型，先搞一个10×10的简单模型验证，等应力分布符合预期了再上线。这样调试速度快，也避免被Abaqus的各种报错信息绕晕。

3.3 用ODB文件导入法：从已有模型无缝传递应力场

ODB导入法比较适合“地应力场很复杂，已经算好了一个稳定应力场，要在它的基础上接着做工程分析”的情况。我常用的流程是这样：先用一个不带开挖的完整地质模型，在Static, General分析步里只施加重力和边界条件，算出稳定状态下的应力场。这个模型可以包含起伏地形、多层地层、断层影响，只要你觉得它足够真实就行。

算完之后，正式工程模型的开挖部分通常要在这个地质模型上“切”出来。为了省去重新设置初始应力的麻烦，我会在正式模型的关键字里加入：

*Initial Conditions, type=stress, file=geostatic.odb

这个写法的意思是：从geostatic.odb这个输出数据库里读取应力场，作为正式模型的初始应力条件。注意两个模型的几何位置和坐标系必须完全一致，否则应力场映射会出错。如果你的正式模型网格和地质模型网格不完全一致，Abaqus会根据网格节点坐标做插值，通常问题不大，但网格差异过大会导致应力场不光滑。

还有一种常见做法是把应力场导出成数据文件，再用*Initial Conditions, type=stress, input=xxx.dat读入。数据文件格式一般是：单元号或单元集名，然后跟着S11、S22、S33、S12等应力分量。这个方法的好处是你可以在导入前对数据进行后处理，比如人为调整K0、滤掉某些奇异点的应力值。缺点是数据文件可能很大，手动编辑不现实，最好通过Python脚本自动生成。更详细的格式建议参考Abaqus Keywords Reference Manual，不同版本之间稍微有点差异。

3.4 判断地应力平衡成功的3个硬指标

很多朋友做完地应力平衡后不确定自己到底算没算对，就盯着云图颜色瞎猜。我总结了三个可量化的硬指标，满足这三条基本就算平衡成功。

第一条，第一步分析能收敛。Geostatic分析步如果一直不收敛，或者每步都疯狂迭代，说明初始应力与荷载或者边界条件不匹配。此时先不要急着往下算，赶紧检查材料参数、单位、约束和应力输入。

第二条，位移量级足够小。在ODB里查看第一个分析步结束时的U magnitude，好的平衡结果是10⁻⁴m以下，稍差一些也要在10⁻³m量级。如果你的模型尺寸是几百米，位移超过0.01m基本就是不合格的，需要在后面分析里人为减去初始位移或者重新修正初始应力。

第三条，应力场分布合理。查看S22（竖向应力）云图，应该基本随深度线性增加，且最大值接近ρgh理论值；水平应力S11大约等于K0倍竖向应力。每条深度的应力曲线拉出来应该是一条平滑直线。如果应力云图里面出现斑块状、锯齿状，多半是网格质量或者初始应力插值出了问题。

4. 初始地应力场设置中的常见报错与排查技巧

4.1 自动平衡不收敛、负特征值，先看这6个原因

做初始地应力场时最常见的报错是Geostatic分析步不收敛，或者出现负特征值警告。我自己排查这类问题基本按下面这个顺序来，命中率很高。

第一，单位不一致。密度、尺寸、弹性模量、重力加速度，任何一个单位没统一，应力场就会错得离谱。检查方式很简单，算一下模型最深处的理论自重应力ρgh，再对比初始条件里输入的应力值，量级不应该差太多。

第二，缺边界条件或约束不足。模型如果缺少必要的约束，在重力和初始应力平衡过程中会出现刚体移动，Abaqus会报零主元或负特征值。记住岩土模型的标准配置：底边固定竖向，左右两侧约束法向，必要时还应在前后方向加约束。

第三，材料参数有问题。弹性模量太小、泊松比取值异常、密度没赋上，都会导致收敛困难。特别是有些模型用了线弹性材料，弹性模量低到几十MPa，又刚好处于高应力区，变形量太大，平衡就很困难。

第四，初始应力输入方向或数值符号不对。Geostatic数据行的应力值正负搞反、侧压力系数填得太大，都会让初始应力与重力不匹配。

第五，网格质量太差。长细比夸张的单元、严重扭曲的单元，在应力平衡时会产生局部奇异，Abaqus计算不收敛的概率会明显上升。

第六，模型里有不该参与初始平衡的接触或边界条件。比如你设置了接触对、弹簧、阻尼器，它们会干扰Geostatic分析步的平衡过程。对于这类组件，建议在初始应力平衡阶段通过Model Change将它们暂时移除，或者不在此阶段激活，等平衡完成后再激活。

4.2 初始应力与塑性屈服同时出现怎么办

深埋高应力区做地应力平衡时，还有一个让人头疼的报错：initial stress exceeds yield stress，或者说初始应力已经超过材料屈服强度。这在高埋深软岩、高地应力区很常见。Abaqus在力平衡前要检查初始应力是否在屈服面内，如果不在，会产生大量塑性应变，平衡就乱了。

处理办法有三条路。第一条，把第一步平衡分析改为弹性模型。也就是说，先用线弹性材料跑地应力平衡，让应力场稳定下来；然后在中途切换到弹塑性材料，通过Field或材料状态变量把泊松比、屈服强度等参数更新成真实值。第二条，如果材料本来就是弹塑性，可以在初始应力设置时把应力水平整体调低一点，确保初始状态处于弹性范围内，再在后续分析中通过荷载逐步增加到真实应力。第三条，使用自动平衡并配合Abaqus的初始应力修正，让Abaqus在迭代过程中自动调整应力使其回归屈服面。这种方法需要特别小心，因为Abaqus可能会将超出屈服面的应力投影回屈服面，导致初始应力场与目标应力场产生偏差。

我个人的建议是：对于深埋高应力岩体，优先采用“弹性试算+塑性切换”的方式。先算出一个满足平衡条件的弹性初始应力场，确认位移清零后再通过重启动或者Field切换材料参数。这样做既保证了初始应力场稳定，又允许后续分析充分反映塑性行为。

4.3 环境类问题速查：libpng error、GPU加速、中断卡死

除了模型本身的问题，Abaqus运行环境也会在初始地应力调试阶段捣乱。按你搜到的热词，我整理几个常见的环境坑。

一是libpng error。这个错误通常在Abaqus启动或者打开CAE、ODB时弹出，表现为一个带“libpng error”字样的警告框，有些版本会直接导致图形界面异常。多数情况是显卡驱动与Abaqus自带的图形库不兼容。解决办法：更新显卡驱动；在环境文件abaqus_v6.env里设置相关图形选项；如果还不行，可以在命令行提交计算，完全绕开图形界面。这个问题不影响inp模型的求解，所以遇到时不用太慌。

二是GPU加速。Abaqus/Explicit支持GPU加速，Abaqus/Standard从部分版本开始也能用GPU加速某些求解器。启用GPU之前，先确认你的显卡型号、驱动版本、CUDA版本和Abaqus版本匹配。如果不匹配，最直接的表现就是计算中途报错或者速度反而更慢。地应力平衡这类小模型通常用不到GPU，建议关闭GPU加速，用CPU多核跑反而更稳。

三是运行中中断不了。Job运行时点Stop没反应，或者卡在“Writing ODB”这一步。常见原因是系统资源占用过高，或者ODB文件被其他程序锁定。可以先尝试等一会儿，如果还不行就打开任务管理器结束Abaqus相关进程。写ODB时被杀进程容易留下损坏的ODB文件，下次计算前建议把原ODB删掉或者另存一个新名称。

四是“节点没有连接到任何单元”的警告。某些网格操作或删除单元后，模型里会残留孤立节点。这类节点不参与计算，但会在输出诊断信息里反复出现，干扰你判断真正的报错。找孤立节点可以用Mesh模块的Verify功能检查，也可以用Python脚本遍历网格，把没有归属单元的节点ID列出来，然后通过Edit Mesh或者重新建模清理掉。

我把这些环境坑放进速查表，便于对照：

5. 进阶实战：焊接仿真、cohesive单元和Voronoi模型中的应力场处理

5.1 焊接仿真为什么不能直接照搬地应力平衡的思路

焊接仿真在Abaqus里越来越多见，但要注意，焊接中的“应力”和岩土中的“初始地应力”并不完全是一回事。焊接模拟的核心是热-力耦合，材料经历快速升温、局部熔化、冷却收缩，最终形成残余应力场。岩土里所谓的初始地应力是为了在计算开挖前让模型处于自平衡的天然应力状态，而焊接模拟里，你通常不是先给整个工件一个“初始应力”，而是通过移动热源逐步把热应力算出来。

真正和“初始应力场”沾边的是多道焊模拟。焊接完第一道之后，工件里已经存在残余应力，第二道焊要在这个残余应力基础上继续计算。这时候就可以把第一道焊接算出来的应力场通过ODB导入或者重启动的方式作为第二道焊的初始状态。方法上可以参考第3.3节ODB导入法，但要注意焊接模型里还有温度场、材料状态、单元生死等额外变量，导入时必须把温度和相关状态变量一起传递，不能只传应力。否则第二道焊的温度场和应力场对不上，计算结果虚得没法看。

5.2 cohesive单元搭配Voronoi模型做岩石破裂时，初始地应力怎么给

cohesive单元和Voronoi模型的组合，现在很多做岩石破裂、混凝土断裂、多晶材料损伤的朋友都在用。Voronoi模型把材料划分成很多不规则的多边形“块体”，cohesive单元则铺在块体边界上，用来模拟裂缝的萌生和扩展。这种模型在引入初始地应力时，会踩一个很典型的坑：初始应力平衡阶段，cohesive单元在还不需要开裂的时候就已经提前损伤甚至破坏了。

原因是，地应力平衡阶段单元之间会有很大的压应力或者剪应力，如果cohesive单元的损伤初始阈值设置得比较低，或者初始刚度比较小，它可能在平衡过程中就被“压坏”了。等后续正式加载时，模型里全是已经损伤的cohesive单元，裂纹还没加载就出现了，完全失真。

我有两个比较实用的处理思路。一个是在初始地应力平衡阶段，暂时不让cohesive单元参与计算。可以通过Model Change功能把cohesive单元所在的set在Geostatic分析步开始时移除，等平衡完成后的下一个分析步再重新激活。重新激活时，cohesive单元虽然没有继承初始应力，但对于裂缝模拟来说，只要块体单元已经处于正确的应力状态，cohesive的初始应力可以通过界面本构的初始间隙间接体现，很多研究都是这样简化的。另一个思路是把cohesive单元的损伤起始位移在初始平衡阶段设得非常大，同时保持弹性刚度足够大，让它在这个阶段“坚不可摧”，等平衡完成后再通过材料参数切换把真实损伤参数换回来。

这两种方法我都试过，Model Change方式更干净，但对单元重激活时的数值稳定要求更高；参数切换方式操作起来直观，但要注意切换瞬间可能带来应力突变。具体选哪种，要看你研究问题的重点。如果是做岩石破裂过程，我推荐用Model Change，把cohesive单元的影响留到真正加载阶段。

5.3 初始地应力场与后续动力分析、开挖卸载的配合

最后再说一个经常被忽略的衔接问题。初始地应力平衡完之后，后续分析可能是静力开挖，也可能是地震动力响应，这两者对初始应力场的要求不完全一样。静力开挖相对简单，平衡完直接进入开挖步即可，位移云图会从接近于零的初始状态重新变化。但动力分析时要特别注意，初始应力场必须能平稳地转入动力分析步，否则在第一个动力增量步会产生巨大的不平衡力，相当于给模型来了一记瞬间冲击。

处理方法是，在转入动力分析之前，先加一个Static, General稳态分析步，让地应力平衡后的应力场平稳过渡到动力分析的初始状态；或者在动力分析中使用*Initial Conditions续传应力场，并结合阻尼设置吸收可能出现的数值振荡。另外，开挖卸载模拟中如果要用到单元生死，被移除的单元里的初始应力也要按顺序释放，不能一下子全去掉，否则会在开挖边界上产生剧烈的应力重分布，导致周围单元瞬间进入塑性。更合理的做法是通过多个分析步分级降低被挖单元的模量，模拟应力逐步释放的过程，再移除单元。

我自己做这类项目时有个习惯：无论用什么方法设置初始应力场，都会在正式计算前单独跑一个“地质模型+初始应力平衡”的小版本，把平衡结果和理论值核对一遍。这一步工作看起来多花了几分钟，却能避免后面整个工程模型因为一个初始应力错误而白跑几天。尤其是模型里同时有Voronoi、cohesive、热力耦合这些复杂要素时，前期的地应力平衡越扎实，后面的问题越少。