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里取出每个节点的坐标和位移。具体做法获取结构单元的节点ID列表。对每一个节点读取其全局坐标和全局位移。建立一个局部坐标系以单元平面为xy平面法向为z。将位移转换到局部坐标系。利用三角形单元的形函数偏导数计算面内应变分量εx、εy、γxy。结合壳理论将面内应变和弯曲应变叠加得到上下表面的应变。最后做特征值分解求主应变。这里给出一个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. 三种方案对比总结与我的最终建议表格看得比较直观方案原理适用阶段编程难度精度广义胡克定律反推应力-应变本构线弹性低弹性范围高塑性失效增量叠加法应力增量-应变增量-累加弹塑性中较好主应力旋转时误差大位移几何法节点位移-形函数导数-应变通用高最准确无本构依赖贴层zone法间接近似通用低近似可能扰动模型如果是普通项目赶工期、只需要曲线趋势直接上方案一。如果材料局部进入塑性且你希望结果更可靠些上方案二注意控制主应力旋转的影响。如果项目本身就需要应变张量、主应变方向的精确值或者你在做学术研究直接上方案三一劳永逸。说一个我踩过的坑方案一里的弹性模量和泊松比设置FLAC3D结构单元本构里你可能用了struct shell property young... poisson...但fish里读取的时候要确认读的是不是同一个属性。有时候你建模时用的prop命令和fish接口访问的默认属性名不一致导致算出离谱的应变。建议读者在计算前先手算一个简单例子验证——建立一个单轴拉伸的壳单元施加一个已知力然后手动算一下理论应变跟fish输出比对误差在1%以内说明参数通路是通的。最后再分享一个小经验结构单元的应变分析不要孤立地只看单元本身的应力和应变还要结合周围zone的变形云图一起看。结构单元和zone之间的协调关系是否合理往往是判断模型是否“假收敛”的一把尺子。我遇到过几次结构单元应力看起来很光滑、但fish算出来的应变在局部区域出现跳变的情况最后排查下来都是网格过渡区刚度不匹配导致的属于模型问题不是算法问题。所以无论用哪个方案先画一遍应变云图看看空间分布是否连续这个习惯能帮你省掉大量后期排查的麻烦。