企业官网建设流程全解析

1. 为什么结构单元只能看主应力，却拿不到主应变

做FLAC3D隧道或边坡支护设计的朋友，十有八九都撞上过同一堵墙：软件里结构单元（壳单元、衬砌单元、土工格栅这些）的后处理结果，应力路径给得特别全，最大主应力、最小主应力、中间主应力，再加上个direction cosine，想看哪个看哪个。可一旦你想看看主应变，对不起，没有。后处理面板里翻遍了，fish里把结构单元的所有属性名都过了一遍，就是找不到跟应变相关的分量。

这个问题看着不起眼，真做项目的时候挺膈应人。尤其是做监测反分析、做支护结构安全评估、或者写论文需要画应力应变曲线的时候，你手里只有应力没有应变，整个分析链条就断了一截。更烦的是，FLAC3D自带的文档对这件事基本上是一笔带过，既没有明确说“结构单元不计算应变”，也没有告诉你怎么绕过去，就让你自己在里面瞎折腾。

我最早碰到这个问题是在做一段浅埋暗挖隧道的二次衬砌安全性评估，当时需要在报告中给出“初支主应力/主应变随开挖步的变化曲线”。结果埋着头在fish里翻了一下午，把struct get property能列出的属性都快背下来了，愣是没找到一个跟strain相关的词。后来跟同门一聊，发现大家都在这儿卡过，只是最后各自偷偷摸摸用不同的方式凑了个结果出来。

所以这篇文章想做的事就一个：把FLAC3D结构单元主应变缺失这件事彻底说透。先说清楚软件为什么不肯给，再说通过应力反推应变的几条可行路线。整个过程我会按我自己实际验证过的思路来讲，避免给一堆理论上成立、实际上根本跑不通的花架子。

2. FLAC3D结构单元的应力输出机制：它到底算了什么，又漏了什么

2.1 结构单元的本构关系不是“空壳”，应力是真算出来的

先说清楚一个基本前提：FLAC3D里的结构单元，无论是壳单元（shell）、衬砌单元（liner）还是土工格栅单元（geogrid），在计算过程中是实打实走了一遍本构模型的。也就是说，每一步的应力增量不是后处理时硬凑出来的，而是由材料模型根据当前的变形状态计算出来的。既然有应力增量，就必然存在对应的应变增量，否则本构方程根本推不动。

这里可以拿线弹性本构来说：Δσ = D * Δε，D是弹性矩阵。每一步迭代，程序先通过节点的位移增量算出应变增量，再代入本构算出应力增量，然后更新应力。整个过程和zone（连续体单元）的计算逻辑是完全一致的。

既然计算过程里明明有应变，为什么后处理不给？这里就得说FLAC3D的设计思路了。

2.2 软件为什么不暴露应变结果：一个典型的“够用就好”设计

做数值模拟的人都知道，FLAC3D有个很明显的性格特点：它给你开放的内容，都是它认为你在工程判断中必须用的；而那些它觉得“不常用”的中间量，往往就直接不更新到内存里。

结构单元的应变就是这个命运。程序在计算过程中确实产生了应变增量，但它完成应力更新之后，并不会把这部分应变数据持久化保存到结构单元的属性列表里。这跟zone单元是截然不同的——zone单元有专门的一整套应变输出机制（包括zone.strain.xx等分量），而结构单元从一开始就被设计成一个“面向力与应力”的对象，位移和应力是它的对外接口，应变被当成中间计算垃圾直接扔掉了。

还有一个细节可以佐证这个设计思路：结构单元的破坏判断，全部是基于应力比和受力状态做的。比如衬砌单元里头的混凝土屈服准则，用的是拉应力/压应力判断，没有哪个准则需要直接把应变作为触发条件。程序的设计者自然认为没有输出应变的必要。

所以结论很简单：不是“FLAC3D算不出”，而是“FLAC3D觉得没必要给你看”。但这不代表我们拿不到，因为应力已经摆在那了，应力和应变的关系就写在材料的本构方程里，反过来推就是了。

2.3 结构单元已有的应变替代项：你真的全找到了吗

在往下走之前，我建议读者先去把结构单元的所有属性再翻一遍。翻的方法很简单，在FLAC3D里随便建一个结构单元，然后用fish遍历它的属性名：

; 建一个最简单的壳单元示例 struct shell create by-size ... ; 遍历属性 local s = struct.near(0,0,0) local attrs = struct.attribute.list(s) loop foreach a (attrs) io.out(a) endloop

我当年就是这么干的，最后在长串列表里发现了一个叫bending_moment或force之类的力类属性，以及一大组应力分量。但真的没有strain这类属性。

有人可能会说：shell单元不是有strain结果吗？zone里的zone.strain.invariant是一堆的呀。你注意区分——结构单元的属性和zone单元的属性是完全不同的两套命名空间，很多人在fish脚本里写zone.strain去访问结构单元，结果返回null，就是这个原因。

但这篇文章真正想讲的不是抱怨设计，而是给你几条真正能落地的主应变求解方案。

3. 方案一：基于广义胡克定律，从主应力直接反推主应变

3.1 原理篇：线弹性本构下主应力和主应变的关系

在结构单元处于线弹性阶段时，应力和应变的关系完全由弹性模量E和泊松比ν决定。假设你已经从FLAC3D里拿到了三个主应力σ1、σ2、σ3（这个很简单，后处理直接能读），那对应的主应变ε1、ε2、ε3可以通过广义胡克定律求得：

• ε1 = (σ1 - ν*(σ2 + σ3)) / E
• ε2 = (σ2 - ν*(σ1 + σ3)) / E
• ε3 = (σ3 - ν*(σ1 + σ2)) / E

这个公式对大家都熟悉，关键问题不在公式，在于三个使用前提：

第一个前提：材料处于线弹性阶段。如果你的衬砌已经进入塑性（FLAC3D默认混凝土模型虽然不是摩尔库伦，但同样存在拉裂、压碎破坏），那么线弹性的应力-应变对应关系就已经失效了。因为塑性阶段的应变包含不可恢复的塑性部分，而应力无法唯一决定总应变。

第二个前提：你取到的三个主应力值是总应力，不是增量。只要FLAC3D的应力输出是基于总应力更新的，那用它反推出来的就是总应变，这点不用担心。

第三个前提：材料是各向同性的。如果是加了钢筋的混凝土衬砌，用复合材料等效参数的时候，E和ν取的是等效值，这时候反推出来的应变是宏观等效应变，不是混凝土基质或钢筋各自的应变了。对工程分析来说，这个等效结果往往是够用的。

3.2 实操篇：fish里怎么把主应变算出来并输出曲线

有了公式，剩下就是编程实现的问题。我的做法是在FLAC3D里写一段fish函数，遍历目标结构单元，读取三个主应力，然后带入上式。

下面这段代码是我实际用过的简化版本，大家可以根据自己的模型微调：

; 从结构单元读取主应力并计算主应变 ; 假设结构单元类型为shell，且材料为线弹性 def calc_principal_strain local s = struct.near(0,0,0) ; 找目标结构单元 local sig1 = struct.stress.1(s) local sig2 = struct.stress.2(s) local sig3 = struct.stress.3(s) local E = 3.0e10 ; 弹性模量 local nu = 0.2 ; 泊松比 eps1 = (sig1 - nu*(sig2 + sig3)) / E eps2 = (sig2 - nu*(sig1 + sig3)) / E eps3 = (sig3 - nu*(sig1 + sig2)) / E end calc_principal_strain io.out('eps1 = ' + string(eps1))

说明几个地方：

• struct.stress.1返回的是最大主应力，struct.stress.2是中间主应力，struct.stress.3是最小主应力。这里的顺序FLAC3D内部已经排好了，不需要你自己排序。
• 这个方法对所有结构单元类型（shell、liner、geogrid）理论上都适用，因为广义胡克定律跟单元类型无关，只跟材料有关。
• 如果你要输出随时间/开挖步变化的曲线，可以在solve循环的每一个阶段调用这段函数，把算出的eps1存到一个历史变量里，用history记录即可。

3.3 方案一的局限：一旦材料屈服，误差有多大

我在做某个围岩级别比较差的隧道案例时，初期支护的混凝土衬砌其实已经局部进入了塑性状态。这时候用线弹性公式算出来的应变，和真实的总应变（包含塑性应变）会存在不小的偏差，尤其是最大主应变方向如果有拉裂缝产生，弹性反推的结果会明显偏小。

那怎么办？有两条路：一是只在弹性阶段进行主应变分析，塑性区明确标记“本方法不适用”；二是使用下一节讲的增量叠加法，从计算过程的应变增量入手。

方案一的适用场景很明确：材料没有进入塑性，或者塑性区范围很小不影响整体趋势判断。像衬砌设计中的正常使用状态验算（裂缝宽度控制之前），一般还在这个框架内。

4. 方案二：用增量法重构应变历史（适用于塑性阶段）

4.1 核心思路：既然每一步都算了应变增量，那就把它攒起来

前面提过，FLAC3D在计算过程中每一步都产生了应变增量，只是没保存。但有一个东西它是保存了的——每一步的应力增量（或者说是更新后的应力状态）。如果我们能从相邻两步的应力状态差反推出应变增量，再把这些增量累加起来，就能重构出总应变，这个思路在处理塑性问题时也基本成立。

这里的分寸在于：塑性阶段的应力-应变关系不再是简单的线性映射，但如果我们把步长取得足够小（FLAC3D的力学时步本来就很小），那么每两个邻近输出点之间，可以近似认为材料处于线性加卸载的路径上。用该时刻的切线刚度（或当前屈服状态对应的弹性刚度）来建立应力增量和弹性应变增量之间的桥梁：

Δε_e = D^(-1) * Δσ

然后：

ε_total = Σ Δε_e + ε_plastic

问题是ε_plastic怎么来？通常情况下，如果你只关心弹性应变部分（工程上很多时候关心的就是弹性的那一部分，比如判断裂缝宽度），那么上面这个增量反推法已经够用了。如果你想得到包含塑性应变的完整总应变，那你需要做一个很麻烦的事情：在fish里挂钩struct的每一步计算，将应变增量在本地累加。

4.2 关键接口：fish回调与结构单元应力读取

FLAC3D支持fish回调（callback），我们可以利用它在每一步力学计算结束之后，读取当前应力状态，做差分，然后累加。

这里给出一个简化框架：

global eps1_acc = 0.0 global eps2_acc = 0.0 global eps3_acc = 0.0 global sig1_prev = 0.0 global sig2_prev = 0.0 global sig3_prev = 0.0 global first_flag = 1 fish define update_strain local s = struct.near(0,0,0) local sig1 = struct.stress.1(s) local sig2 = struct.stress.2(s) local sig3 = struct.stress.3(s) if first_flag = 1 then sig1_prev = sig1 sig2_prev = sig2 sig3_prev = sig3 first_flag = 0 else local dsig1 = sig1 - sig1_prev local dsig2 = sig2 - sig2_prev local dsig3 = sig3 - sig3_prev local E = 3.0e10 local nu = 0.2 ; 主应力方向在增量过程中可能发生旋转，这里先做简化近似 eps1_acc += (dsig1 - nu*(dsig2 + dsig3)) / E eps2_acc += (dsig2 - nu*(dsig1 + dsig3)) / E eps3_acc += (dsig3 - nu*(dsig1 + dsig2)) / E sig1_prev = sig1 sig2_prev = sig2 sig3_prev = sig3 endif end

然后把update_strain挂到callback里面：

fish callback add update_strain

或者你在solve循环里自己控制，每走n步手动调一次。

4.3 这个方案的坑：主应力方向旋转的问题

这个方案最大的坑在于：主应力的方向不是固定的。每一步三个主应力对应的方向可能都在变（尤其在开挖卸荷、支护施作后的应力重分布阶段），而你用struct.stress.1这种接口只能拿到数值，拿不到方向。

如果你的模型应力路径比较温和、主应力方向基本不变，那这个增量反推法还是比较准的。但要是主应力方向出现了明显旋转（比如大变形隧道的分步开挖），那你在dsig层面做的“三个分量分别差分”就是在一个混合坐标系里操作，算出来的增量方向意义就会打折扣。

更严格的做法是：每一步都同时读取主应力方向和主应力大小，然后把前后两步的应力张量在同一个固定坐标系下做差分，再求主应变。这需要用到struct.stress.direction之类的接口，工作量会大不少，但对精度要求高的项目还是值得做的。

4.4 方案二 vs 方案一：怎么选

我的经验是分情况：

• 如果你只是做常规的支护结构受力分析，且构件基本处于弹性状态，直接用方案一，简单、快、不折腾。
• 如果你的结构局部进入了塑性，但你想看的只是某个特定位置的总体应变趋势，方案二也够用，只要别太纠结塑性应变的绝对精度。
• 如果你要做严谨的塑性应变解耦（弹性应变和塑性应变的分离），那方案二还不够，建议走fish应变积分方案（下一章详述）。

5. 方案三：自己动手，用fish对节点位移做几何方程求解

5.1 原理：应变的本质是位移的导数，和应力无关

前面两个方案本质上都是“借力打力”，通过本构关系从应力反推应变。但不要忘了应变的原始定义：应变是位移场的空间梯度。FLAC3D的结构单元是离散的有限单元，它的每个节点都有位移（这个是可以直接读取的）。只要我们能拿到一个结构单元所有节点的位移，再用有限元形函数求导，就可以得到单元的应变，完全绕开应力-应变本构关系。

这个方案的好处非常明显：它不依赖材料是否线弹性，不依赖主应力方向是否旋转，因为它是直接从运动学（几何方程）出发的，物理意义最干净。

坏处是：编程量比前两个方案大不少，而且需要你对有限元形函数有一定基础。

5.2 实现思路：壳单元局部坐标系下的应变计算

以三节点的壳单元（triangular shell element）为例。FLAC3D的shell单元的几何在局部坐标系下是一个平面三角形（或者四边形），我们可以在fish里取出每个节点的坐标和位移。

具体做法：

1. 获取结构单元的节点ID列表。
2. 对每一个节点，读取其全局坐标和全局位移。
3. 建立一个局部坐标系（以单元平面为xy平面，法向为z）。
4. 将位移转换到局部坐标系。
5. 利用三角形单元的形函数偏导数，计算面内应变分量εx、εy、γxy。
6. 结合壳理论，将面内应变和弯曲应变叠加，得到上下表面的应变。
7. 最后做特征值分解求主应变。

这里给出一个fish片段思路（以三角形单元为例，只展示面内应变的核心计算）：

; 伪代码框架 ; 节点坐标: p1, p2, p3 ; 节点位移: d1, d2, d3 ; 计算局部坐标下的几何矩阵B，然后 B * 节点位移 = 应变 ; 这里需要先建立局部坐标系，做坐标变换

完整的代码比较长，我不全贴在这里。但思路务必记住：节点位移是FLAC3D一定会输出的量，结构单元的位移接口包括struct.node.pos和struct.node.disp。有了位移，一切皆有可能。

5.3 一种更省事的替代：把结构单元转成zone来做应变分析

如果你只是偶尔需要某个位置的应变，实在不想写形函数求导的代码，我给你一个土办法：在结构单元附近建立一个非常薄的zone层，紧贴在壳单元或衬砌单元上，然后观察这个zone层的应变。

原理很简单：如果zone层足够薄、刚度适当，它的变形基本上和结构单元是一致的，那么zone层的主应变就非常接近结构单元的主应变。这个方法的优点在于FLAC3D对zone的应变输出非常完整——zone.strain.xx、zone.strain.yy、zone.strain.xy等全都有，直接后处理就能看主应变方向和大小的云图。

缺点也很明显：zone层的存在会略微改变模型的受力状态（因为增加了额外刚度），而且zone层太薄时会引发网格畸变或计算效率下降。所以这个方法仅推荐作为快速验证手段，不推荐作为最终结果。我自己一般是在做初步的参数敏感性分析时用这个招数快速建立“应力-应变”对应关系，等确定了关键位置再用方案一或方案二做精细化提取。

6. 从主应变到工程决策：别算完就完事

6.1 拉应变临界值：什么时候该担心

很多朋友费了老劲把主应变算出来，结果不知道怎么用。这里说几个常用的工程判断思路，供大家参考。

混凝土类支护结构，最关心的通常是最大主拉应变。它跟裂缝直接相关。按照混凝土结构设计规范里关于裂缝宽度验算的思路，混凝土极限拉应变一般在100~150微应变（1e-4到1.5e-4）之间。超过这个值，意味着混凝土很可能开裂。

我做支护设计时一般按这个经验来判断：

• ε1 < 0.5e-4：安全，不需要额外关注。
• 0.5e-4 < ε1 < 1.2e-4：处于临界区，结合应力比一起看。
• ε1 > 1.2e-4：大概率开裂，需要检查是否存在超过抗拉强度的区域。

但这只是静态判断，真正要小心的是主应变的方向。如果你算出来的最大主应变方向和结构单元的环向（如果是隧道衬砌）夹角超过30度，那裂缝形态会从横向缝变成斜裂缝甚至纵向裂缝，治理思路完全不同。

6.2 主应变塑性区识别：一个比应力更直观的指标

另外一个非常实用的用途是判断塑性区的分布范围与贯通趋势。

用主应力判断塑性区有个问题：压应力区也可能出现高应力但并未破坏（混凝土三向受压时强度会大幅提高），而主应变对破坏形态的刻画更直接——拉应变集中区往往就是潜在开裂区，剪应变集中区往往就是剪切破坏带。

我在边坡锚杆（土工格栅）的模拟中经常用最大剪应变γmax = ε1 - ε3来识别潜在滑动面。FLAC3D自带的zone塑性区指示器用的是应力状态，但如果你把结构单元（比如土工格栅加筋层）的ε1 - ε3画出来，往往能比应力云图更早看到塑性应变局部化的苗头。

6.3 双向受力状态下的主应变方向应用

主应变方向另外一个重要的应用是判断结构单元的主要受力方向。比如土工格栅加筋土挡墙中，格栅的主应变方向通常可以指示加筋体的主受力方向。如果你的格栅布置方向和主拉应变方向夹角太大，说明你的布筋方向没对准，加筋效率在打折扣。

这个分析如果是用方案一算出来的，还需要额外做一步：从FLAC3D读取主应力方向，用主应力方向作为主应变方向（在线弹性条件下，各向同性材料的主应力和主应变方向是完全一致的）。如果是用方案三直接求的应变张量，那特征向量直接就是主应变方向，更直接准确。

7. 三种方案对比总结与我的最终建议

表格看得比较直观：

如果是普通项目，赶工期、只需要曲线趋势，直接上方案一。如果材料局部进入塑性且你希望结果更可靠些，上方案二，注意控制主应力旋转的影响。如果项目本身就需要应变张量、主应变方向的精确值，或者你在做学术研究，直接上方案三，一劳永逸。

说一个我踩过的坑：方案一里的弹性模量和泊松比设置，FLAC3D结构单元本构里你可能用了struct shell property young=... poisson=...，但fish里读取的时候要确认读的是不是同一个属性。有时候你建模时用的prop命令和fish接口访问的默认属性名不一致，导致算出离谱的应变。建议读者在计算前先手算一个简单例子验证——建立一个单轴拉伸的壳单元，施加一个已知力，然后手动算一下理论应变，跟fish输出比对，误差在1%以内说明参数通路是通的。

最后再分享一个小经验：结构单元的应变分析，不要孤立地只看单元本身的应力和应变，还要结合周围zone的变形云图一起看。结构单元和zone之间的协调关系是否合理，往往是判断模型是否“假收敛”的一把尺子。我遇到过几次结构单元应力看起来很光滑、但fish算出来的应变在局部区域出现跳变的情况，最后排查下来都是网格过渡区刚度不匹配导致的，属于模型问题，不是算法问题。所以无论用哪个方案，先画一遍应变云图，看看空间分布是否连续，这个习惯能帮你省掉大量后期排查的麻烦。