做基坑、隧道、边坡支护数值模拟的同行应该都有过这种经历模型算完了云图也导出来了领导或审图单位问你要桩身弯矩图、轴力图你盯着屏幕上的应力云图一时语塞——实体单元明明只有应力哪来的弯矩和轴力这个问题在 FLAC3D 里特别典型。结构单元如 beam、pile 可以直接用sel recover或者 FISH 拿到截面内力但很多人建模时为了跟地层贴合、考虑接触或者模拟实际几何形态用的是实体单元。实体单元输出的只有应力张量没有现成的内力和弯矩。偏偏设计环节最关心的就是内力于是“实体单元提取弯矩轴力”就成了每个用 FLAC3D 做结构分析的岩土工程师绕不开的活。这篇东西我就只讲一件事在 FLAC3D 6.0 里怎么用 FISH 把实体单元的应力积成弯矩和轴力。梁、隧道衬砌、桩三类场景我都会拆开说代码是 6.0 版本的5.0 的朋友可以看思路但别直接抄新版很多东西变了。文章会涉及截面法的基本原理、建模时的分组埋点、完整的 FISH 代码框架以及我踩过的一些坑。适合已经会用 FLAC3D 建模计算、但还没亲手写过应力积分代码的读者也适合那些想从 5.0 迁到 6.0、正被新 API 搞懵的人。1. 思路先行实体单元的内力怎么从应力里变出来1.1 实体单元只有应力截面内力全靠“积”先想明白一个基本问题实体单元计算结束后每个 zone 上存的是什么是应力张量六个分量——三个正应力、三个剪应力。FLAC3D 的云图能给你看应力场、位移场但它不会替你回答“这根桩在埋深 5 米处承受多大轴力”这种问题。原因是实体单元本身没有“截面”的概念。单元只是空间里的一块材料受力状态用应力表示而轴力、弯矩属于“截面内力”二者差了本质上的一步截面内力需要对截面上的应力做积分。公式其实很朴素以截面法向为 x 方向为例轴力 N ∫σxx dA弯矩绕 z 轴 Mz ∫σxx · y dA这里 y 是应力作用点到中性轴的距离。也就是说把截面“切开”把面上每个微元的法向正应力乘上微元面积加起来就是轴力再乘上离中性轴的力臂加起来就是弯矩。这就是材料力学里最基础的截面法放到 FLAC3D 里不过是把这个积分改成对所有 zone 做累加。关键点在于“怎么把离散的 zone 攒成一个连续截面”。三维模型里没有天然的切割面你需要在代码里指定一个几何平面然后找出这个平面附近参与积分的单元把它们的应力贡献累加起来。1.2 公式就三行但符号约定最容易被忽视轴力拉为正还是压为正弯矩的绕向怎么定义中性轴取在哪里这三个问题看起来小但实际做数据的时候能把人绕晕。我的习惯是统一按国内结构设计的常用约定走轴力以受拉为正受压为负弯矩的正负不强行约定算完直接给曲线让看图的人按工程习惯自己判读只要出图时代码里的符号保持一致就行。更需要注意的是截面的“法向方向”。同样一个截面法向取正 x 和负 x积出来的轴力符号是完全相反的。代码里用的法向量其实就是你心里定义的“截面朝向”。我在下面给的示例代码里法向是用向量(nx, ny, nz)明确写出来的改这个向量就可以控制符号和提取方向这个小参数别忽略。中性轴位置也是个容易含糊的点。截面有几何形心和应力中性轴严格意义上弯矩积分的力臂应该取到截面几何形心而不是应力为零的那个轴。因为实体单元应力里包含了轴力分量应力零点会偏离形心如果你用应力零点去当转轴算出来的弯矩会混入轴力效应。工程上取几何形心就够了除非你做的是精细组合截面分析。1.3 梁、隧道、桩的提取差异一张表说清这三种结构虽然都叫“实体单元提取内力”但切面方向和关注的内力分量差异很大。我做了个表方便对照结构类型切面方向核心内力典型提取方式桩 / 立柱垂直于桩轴的横截面沿桩长分布轴力 N、水平双向弯矩 My / Mz沿桩长间隔切 N 个截面插值绘制内力沿深度图隧道衬砌横断面内的径向截面沿环向分布环向轴力 N、环向弯矩 M每隔 10°~15° 取一个径向切面输出环向内力包络图实体梁 / 支撑垂直于梁轴的截面沿梁长分布轴力 N、弯矩 M剪力 V沿梁长切截面与 beam 单元结果互相验证从代码角度三者共用同一套“应力积分”内核差异只在于截面位置怎么定义、法向量怎么旋转、哪些 zone 参与积分。理解了这一层后面改代码就很顺了。2. 建模时就埋好伏笔分组、切面与初始应力2.1 提前分组比事后用坐标去筛靠谱得多提取内力的第一步是找到参与积分的 zone。最简单粗暴的办法是用坐标范围去圈比如x 10 且 x 10.5。但真实模型里桩、衬砌和土体犬牙交错坐标范围经常圈出非目标单元尤其是在弯折、倾斜的结构附近。我在建模阶段就会做一件事把需要提取内力的结构单独分成一个组。桩归pile衬砌归lining支撑归strut。FLAC3D 6.0 里分组的命令很成熟支持 group 槽位和命名组zone group pile zone group lining分组之后FISH 里判断单元归属就干净了if zone.inGroup(zp, pile) then这比坐标筛选可靠性高一个量级尤其是当结构物本身是倾斜的坐标范围根本没办法写。我做过的几个倾斜桩项目如果不用分组提取代码基本没法写。2.2 切面位置网格规则处才是可靠的“截面”截面应该切在哪里很多人会直接在关心的位置切比如桩顶下 2 米处、衬砌的拱顶位置这没问题。但有个容易被忽略的点别把截面切在网格尺寸突变的地方或者应力集中区。比如桩与承台交界处网格往往加密、应力紊乱在那附近积分出来的轴力是震荡的。我一般在应力集中区外扩 2~3 个单元再切截面宁可损失一点精度也别把局部应力尖峰算进内力里。截面厚度这个概念也要提前想清楚。严格数学意义上的截面是零厚度的面但 FLAC3D 的单元是离散的你不可能让单元刚好躺在一个面上。实际做法是设定一个薄层把落在该薄层内的单元都纳入积分。这个层厚直接决定结果的光滑程度层太薄参与积分的单元少曲线毛糙层太厚结果被“平均化”峰值会被抹平。我在后面代码里用变量d_thick控制这个厚度。2.3 初始应力场不处理提取的轴力直接失真这是实战中最大的坑没有之一。很多人模型是这样做的先建立包含衬砌或桩在内的全模型跑初始地应力平衡然后开挖、施作支护继续计算。结算后直接提取衬砌应力积分。这时候你会发现轴力大得离谱——因为衬砌单元的应力里包含了初始地应力场也就是埋深带来的围压也一并被积分进去了。设计要的是“结构内力”是结构承担的那部分力不是地应力背景。处理办法主要有两种第一种在结构单元激活时把初始应力清零。常见做法是给衬砌、桩单元单独分组激活后统一初始化应力zone initialize stress ... zone group pile但对于算完之后的模型这一招不好使因为你不能随便重置应力场。第二种用 FISH 记录增量应力。在计算开始前遍历目标 zone把初始应力分量存入 extra 变量计算结束提取时用当前应力减去初始应力。FLAC3D 6.0 里 extra 变量是可以给 zone 挂自定义数值的; 计算前 zone.extra(zp, 1) st.xx zone.extra(zp, 2) st.yy zone.extra(zp, 3) st.zz结算后提取时local sx st.xx - zone.extra(zp, 1) local sy st.yy - zone.extra(zp, 2) local sz st.zz - zone.extra(zp, 3)我自己更常用第二种因为不需要改初始化的顺序而且能回溯应力路径。3. FLAC3D 6.0 里的 FISH 实现从单截面到批量输出3.1 6.0 读取应力的正确姿势zone.stress 是结构体如果你以前用的是 5.0应该很熟悉z_stress(zp, 1)这种索引式的应力访问方式。6.0 改成了结构体返回一行代码直接拿六个分量local st zone.stress(zp)访问分量的语法是st.xx st.yy st.zz st.xy st.xz st.yz这个改动让代码可读性好了很多但也坑了一批从 5.0 迁移的人。我一开始也照着旧习惯写z_stress直接报错。6.0 里遍历 zone 的写法也变了不再是loop zp (1, znum)配合z_id而是loop foreach zp zone.list这两种改动是整个提取脚本迁移到 6.0 时最核心的两处。3.2 单截面提取代码拿去改改就能用下面这个代码是一个最小可用的单截面提取框架目标是从实体桩或实体梁的某个截面上提取轴力和两个方向的弯矩。我加了不少注释方便你改成自己的模型。; 实体单元截面内力提取FLAC3D 6.0 / FISH ; 适用桩、梁、柱等实体结构 ; 思路指定一个截面位置与法向找出该截面附近一个薄层内的 ; 目标 zone对法向正应力做面积积分得到轴力 ; 对中性轴的力矩积分得到弯矩。 function get_section_force ; --- 用户输入参数 --- local sec_pos 10.0 ; 截面位置本例为 x 坐标 local nx 1.0 ; 截面法向 x 分量 local ny 0.0 ; 截面法向 y 分量 local nz 0.0 ; 截面法向 z 分量 local d_thick 0.1 ; 截面薄层厚度参与积分的单元范围 local y_neutral 0.0 ; 中性轴 y 坐标几何形心 local z_neutral 0.0 ; 中性轴 z 坐标 ; --- 累加变量 --- local force_n 0.0 ; 轴力 local moment_my 0.0 ; 绕 y 方向轴的弯矩 local moment_mz 0.0 ; 绕 z 方向轴的弯矩 local area_sum 0.0 ; 截面累计面积校验用 ; --- 遍历所有 zone --- loop foreach zp zone.list ; 只处理目标分组后的单元避免误伤土体 if zone.inGroup(zp, pile) then local p zone.centroid(zp) ; 判断单元是否位于截面薄层范围内 local dot (p.x - sec_pos) * nx p.y * ny p.z * nz if math.abs(dot) d_thick then ; 应力张量 local st zone.stress(zp) ; 法向正应力按一般公式只取 n^T σ n 的完全式 ; 本例法向为 x所以退化为 st.xx local sig_n st.xx * nx * nx st.yy * ny * ny st.zz * nz * nz sig_n 2.0 * st.xy * nx * ny 2.0 * st.xz * nx * nz 2.0 * st.yz * ny * nz ; 单元体积折算为面积贡献做法详情见正文说明 local vol zone.vol(zp) local dA vol / d_thick ; 单元中心到中性轴的距离 local ry p.y - y_neutral local rz p.z - z_neutral ; 累加内力 force_n sig_n * dA moment_mz sig_n * ry * dA ; 绕整体 z 轴的弯矩 moment_my sig_n * rz * dA ; 绕整体 y 轴的弯矩 area_sum dA end_if end_if end_loop ; --- 输出结果 --- io.out(截面位置 x string(sec_pos)) io.out(轴力 N string(force_n)) io.out(绕 z 轴弯矩 Mz string(moment_mz)) io.out(绕 y 轴弯矩 My string(moment_my)) io.out(累计面积 string(area_sum)) end ; 调用 get_section_force这里面有两个地方要重点说明。第一法向正应力的完整公式是二次型形式涉及剪应力分量。很多简化代码只取正应力里对应的那个分量比如法向为 x 就直接拿st.xx这在截面法向与坐标轴完全平行时没问题但遇到斜截面会出错。我上面写的是完整公式直接支持任意倾斜截面代价是代码里多了几行乘法实测开销可以忽略。第二体积折算成面积的思路离散模型中截面是“薄片”薄层厚度是d_thick。一个 zone 的体积除以薄层厚度近似等于这个 zone 对截面的面积贡献。严格说这个近似依赖于单元在法向方向上的尺寸与d_thick的匹配关系。我的经验是d_thick取该截面位置附近单元法向尺寸的 1~1.5 倍结果最稳。3.3 批量切面与内力曲线输出单个截面提取只能给一个数实际工程要的是完整的内力分布曲线。比如桩身弯矩沿深度的变化、衬砌环向弯矩包络图都需要批量切面。批量做法就是把上面那段函数包进一个循环截面位置逐步移动每次调用后把结果存进 table; 批量提取桩身内力从 x5 到 x15每隔 0.5 米切一个截面 table.create(N_table) table.create(Mz_table) table.create(My_table) loop for i 1 to 20 local sec_x 5.0 i * 0.5 ; 这里把前面函数里的 sec_pos 改为 sec_x ; 为节省篇幅假设已经封装成函数 get_internal_at(sec_x) ; 返回值为 n_force, mz, my 三个值 local nf 0 local mz 0 local myv 0 ; get_internal_at(sec_x, nf, mz, myv) ; 伪代码示意按实际修改 table.push(N_table, sec_x, nf) table.push(Mz_table, sec_x, mz) table.push(My_table, sec_x, myv) end_loop ; 导出 txt方便画图 table.export(N_table, 桩身轴力.txt) table.export(Mz_table, 桩身弯矩_z.txt) table.export(My_table, 桩身弯矩_y.txt)实际封装时要注意 FISH 的局部变量作用域函数之间传参最好通过全局变量或 table 传递避免嵌套函数变量遮蔽的问题。我踩过这个坑用一个函数里定义的局部变量直接传给另一个函数结果怎么都取不到值检查半天才发现是变量名被外层覆盖了。3.4 弯矩方向分解绕两个局部轴的弯矩怎么算桩和梁在三维空间里受力弯矩不是一个数而是两个方向的弯矩分量。上面代码直接算了绕整体 y 轴和整体 z 轴的两个分量这在结构主轴恰好与整体坐标轴平行时是够用的。但斜桩、斜梁就麻烦了。你需要把整体坐标系下的My、Mz转换到截面局部坐标系下变成绕局部强轴和弱轴的弯矩。做法是先定义截面的局部坐标系法向为n再取一个截面内的方向作为局部 y 轴用叉乘得到局部 z 轴。然后把整体坐标下的应力积分结果投影到局部轴上。这本质上是坐标变换我建议直接用 FISH 里的向量运算函数别手写旋转变换矩阵容易算错。FLAC3D 6.0 的 math 库里有向量的叉乘、点乘函数可以省不少事。代码框架上你只需要把 3.2 节里累加的moment_mz、moment_my从整体坐标分量改成局部坐标分量即可local dot_local_y sig_n * dot(ry_vector, local_y_vec) local dot_local_z sig_n * dot(rz_vector, local_z_vec)这里的ry_vector、rz_vector是单元中心到中性轴的向量在截面内的投影。原理不复杂但要注意向量方向的一致性投影之前先做归一化。4. 实战对比桩、隧道衬砌、实体梁的内力提取4.1 桩身内力沿深度每 0.5 米切一刀桩用实体单元建模时最关心的是桩身轴力和水平弯矩沿竖向的分布。以竖向桩为例截面法向就是竖向切面位置就是不同的深度z。我一般沿桩长每隔 0.5 米切一个截面太密了没有必要曲线也跳。如果是端承桩靠近桩端位置的轴力会明显衰减这时候可以在桩端附近加密切面到 0.25 米能更精细地捕捉轴力传递规律。这里有个容易犯的错如果桩身截面不是圆形而是异形桩比如方桩、T 型桩几何形心不好拍脑袋定。我建议先在建模阶段用fish自动计算截面形心把计算结果存到一个变量里避免在提取代码里手工填。几何形心错了弯矩整体都会带一个系统误差且这个误差不容易被发现因为曲线形态看起来依然合理。4.2 隧道衬砌沿环向切面画内力包络图衬砌和桩不一样的地方在于内力方向。桩提取的是竖向截面上的竖向内力衬砌要提取的是横断面内、沿衬砌环向切线方向的内力。我是这么处理的以隧道中心为原点建一个局部极坐标系每隔 15° 取一个径向截面。这里的“径向截面”就是过隧道中心、沿该角度的半径方向延伸的面。在这个面上积分应力的环向切向分量得到该角度位置的环向轴力 N再对径向方向取矩得到环向弯矩 M。代码上法向量不再固定为(1,0,0)而是随角度变化local theta 15.0 * math.pi / 180.0 local nx math.cos(theta) local ny math.sin(theta) local nz 0.0nx、ny就是该角度处衬砌环向切线方向的分量。然后把 3.2 节的函数里的法向替换掉就是衬砌内力提取脚本。衬砌厚度很小的情况下参与积分的就是沿径向剖开的一小层单元d_thick要调小一点我一般取衬砌单元沿环向尺寸的 1~1.5 倍。输出之后你会得到一圈 24 组(N, M)把 N 和 M 沿角度展开画图就是经典的衬砌内力分布图。拱顶、仰拱、拱肩位置的弯矩极值一目了然拿去对比监测数据或者规范验算都够用。4.3 实体梁跟 beam 单元结果对一对精度才放心实体梁的提取思路和桩几乎一致但我特别建议做一件事在模型里同时建一根 beam 单元作为对照或者单独做一个简支梁小算例用 beam 单元结果和实体积分结果互相验证。我试过一个悬臂梁算例一端固支端部加集中力分别用 beam 单元和实体单元建模然后用 3.2 的代码提取实体梁根部弯矩结果跟 beam 单元查出来的弯矩对比误差在 3% 以内。这个误差来源主要是截面薄层的离散化和中性轴位置误差。对照验证的真正价值是帮你发现代码里隐含的符号错误。实体单元积分出来的轴力如果是负号而 beam 单元查出来是正号不要怀疑是模型问题先查你的法向量方向和中性轴定义。符号问题在单一模型里看不出来一旦对照就原形毕露。5. 常见问题与排查技巧实录5.1 结果大得离谱先查初始应力扣了没有如果提取出来的轴力数量级超出正常范围第一个怀疑对象不是代码而是初始应力场。我前面已经说过这个坑这里再展开一句判断方法很简单在计算结束后的zone.stress(zp)里打印一两个目标单元的应力和单纯的初始地应力对比一下。如果接近说明你积进去的大部分是地应力背景不是结构内力。处理方式是 2.3 节提到的增量应力提取。还有一个更省事的思路如果你只是要内力分布曲线的形态不太在意绝对数值可以用“模型位移清零”前后的应力差来近似。但这个方法仅在地层和结构都处于弹性阶段时有效塑性变形大的模型不适用。5.2 弯矩符号总是“反”的法向和中性轴位置再捋一遍弯矩符号反了十有八九是截面法向方向定义反了。截面法向翻 180°轴力和弯矩全翻号。还有一种是中性轴位置填错比如忘了考虑截面偏心导致力臂的方向反了。我的排查习惯是先在目标截面上挑一个单元手动算一遍该单元对轴力和弯矩的贡献看正负号是否符合预期。比如竖向桩顶部受压那顶部单元对轴力的贡献应该是负的压为负如果代码里出来是正的说明法向取反了。5.3 应力云图连续积分结果却跳变切面别擦着网格边界提取结果出现锯齿状跳变网格边界是最常见的原因。如果切面位置刚好落在网格边界上部分单元被计算进截面部分被排除积分面积会忽大忽小。解决办法是适当增大d_thick让参与积分的单元数量稳定下来。我一般做法是先试一遍调整d_thick让累计面积area_sum保持稳定再用这个参数跑批量切面。area_sum就是提取消错时最好的诊断变量每次提取都应该顺手打印出来检查。5.4 5.0 迁移到 6.0 踩过的坑旧命令别直接搬这个问题太典型了单独列出来说。5.0 时代的 FISH 用z_stress(zp, i)拿应力6.0 改成了zone.stress(zp).xx这种结构体访问方式。遍历方式从loop i (1, znum)改成了loop foreach zp zone.list。分组判断从z_group(zp, 1)变成了zone.inGroup(zp, name)。如果你在 6.0 里直接跑 5.0 的脚本前几行就会报错。我重新封装了一套 6.0 的提取函数之后才发现新旧版本不只是语法差异连遍历顺序都变了后者会影响你输出数据的顺序。如果你要把旧数据跟新结果对比注意单位格点的排列顺序别只看数值不看对应关系。5.5 量纲与单位自查清单最后列一个自查清单每次提取完数据都过一遍检查项典型错误自查方法应力单位模型是 Pa输出按 kPa 处理打印一个单元应力对照手册长度单位米 vs 厘米面积dA的量纲是否对体积折算厚度d_thick忘记按模型尺寸调整累计面积是否接近真实截面面积轴力符号拉压约定不统一与 beam 单元结果对照弯矩单位力·长度长度是米还是毫米数量级是否合理单位问题在数值模型里最容易阴沟翻船。我吃过一次亏模型用米和 Pa 建的dA是平方米应力是 PaN/m²乘出来轴力是 N 没错但我想当然除以 1000 当成了 kN结果一个量级直接错掉。每个中间量都过一遍量纲是省不掉的功夫。做内力提取这一整套流程做下来我个人最大的体会是FLAC3D 里没有“内力云图”这个一键功能但只要你理解截面法的本质用 FISH 把应力积分出来并不难。难的是在你模拟的全过程里始终清楚自己到底想要哪个平面的哪个内力量以及这个量在工程约定里怎么定义正负。我每次写新的提取脚本都会先用一个解析解已知的简单模型验证代码符号和量级确认无误后再套到实际工程模型上。最后再分享一个小技巧提取结果导出成 txt 之后我一般都直接在 Python 里读出来画图不用 FLAC3D 自带的绘图功能。因为批量截面提取之后进一步画包络图、跟规范上限值比线用 Python 处理要方便得多。如果你只是要快速出图table.export之后在表格软件里插入散点图也够用了。