简介：面向从事岩土工程数值模拟的技术人员，这份PDF围绕FLAC3D实体单元无法直接输出梁单元内力的问题，系统讲解弯矩与轴力的提取思路和实现方法。内容从应力积分原理切入，结合FISH函数代码说明截面定位、高斯点遍历、单元组筛选与弯矩臂计算，并延伸到隧道初期支护、地铁支护桩等场景的受力分析；读者可获取完整的函数实现、结果记录与定时调用思路、案例受力变化解读，以及和理论值对比的验证方案。文中还提醒单位换算、截面宽度选取、网格划分等影响精度的细节，适合需要提升模拟结果可靠性的工程师与科研人员参考。包内共1个PDF文件，约104KB，篇幅集中便于按主题查阅；目前已有137人学习，可作为实体单元内力提取的实战型技术资料。

1. FLAC3D 实体单元为什么没有现成的弯矩输出

某跨平台系统里一段隧道衬砌，建模的人图省事，全部用实体单元离散，算完岩土那边要衬砌弯矩图做配筋校核。打开 FLAC3D 一看，zone 能给的只有应力张量六个分量，结构单元那套beam、pile的内力输出在这里一个都没有。实体单元本质上是连续体离散块，它只知道自己是处于受拉还是受压的应力状态，并不知道自己“承担了多少弯矩”——弯矩在这个体系里不是被储存的量，而是应力沿截面一阶矩积分出来的派生量。所以 FLAC3D 实体单元的弯矩轴力提取技术，说到底就是一件事：自己定义截面、自己还原面积、自己把应力积分成 N 和 M。这套方法适合三类人：用实体单元模拟衬砌、桩、拱、剪力墙，却拿不到内力图的；想从 zone 结果反推结构配筋依据的；以及已经用结构单元算过、需要一套独立结果做交叉验证的。

2. 从应力张量到截面内力：积分公式与符号约定

2.1 实体单元的应力存在哪、怎么取

FLAC3D 的 zone 应力在默认的线性单元里是常值，存在单元形心，zone.stress()拿到的就是这六个分量，顺序固定为 xx、yy、zz、xy、xz、yz，单位跟随模型的基本单位制。节点上的应力（gridpoint stress）是后处理阶段从相邻单元外推再平均得到的，看着更连续、画云图更漂亮，但千万别拿它做积分：一旦节点落在两种材料的分界面上，或者落在应力集中区，平均出来的值会同时被两侧单元拉扯，跨界面那一层的应力直接被抹平，积分出来的弯矩能差出一大截。取形心值虽然看上去“块状”，但每个值都对应一块确定的体积，可追溯、可复算。

这里有个容易忽略的前提：zone.stress()返回的是当前计算步的应力全量，包含初始地应力。如果你在自重平衡之后再加荷载，读到的 σ_xx 里已经掺进了岩土那部分压应力，直接积分得到的轴力会带一个很大的常数项，看起来像是“结构凭空多了一千吨轴力”。正确做法是记录平衡步的应力，做增量相减，或者用zone.stress减去zone.stress.initial一类的参考量。我一般会在初始平衡结束后立刻把所有 zone 的应力导出成一个数组存下来，后面每次提取都做差，这条规矩能省掉后面无数次返工。

还有一点：FLAC3D 沿用的是岩土符号约定，压为负、拉为正。而结构专业画弯矩图时习惯的正负跟这个不是一套。积分公式本身不挑符号，但你在写公式和画图之前必须先把约定钉死，否则同一份数据两个人能解读出相反的结论。

2.2 轴力与弯矩的截面积分式

设截面法向为 x 轴，截面形心坐标为 (y_c, z_c)，那么截面上的轴力和两个方向的弯矩写成连续形式就是：

N = ∫ σ_xx dA My = ∫ σ_xx · (z - z_c) dA Mz = -∫ σ_xx · (y - y_c) dA

离散到单元层面，每个被截面切到的 zone 贡献一块面积 A_i，形心坐标 (y_i, z_i)，轴向应力 σ_i，于是：

N = Σ σ_i · A_i My = Σ σ_i · A_i · (z_i - z_c) Mz = -Σ σ_i · A_i · (y_i - y_c)

形心由面积加权求得：z_c = Σ A_i·z_i / Σ A_i，y_c同理。这里My里的减号不是笔误，是右手坐标系下绕 y 轴转动、以 z 向力臂为正的约定；Mz前面那个负号来自 y 向力臂在右手系里方向相反。这两个符号你在自己项目里可以改成别的写法，但一旦定了，全项目所有脚本、所有图例都要统一，不然对图的时候会怀疑人生。

公式里最容易被轻视的是A_i。它不是单元体积，也不是单元某个面的面积，而是“这个单元在截面上占据的投影面积”。对于沿 x 方向规则划分的六面体网格，最省事的还原方式是A_i = V_i / dx，V_i 是单元体积，dx 是该层单元沿 x 方向的厚度。这个近似在单元接近长方体时误差很小，一旦网格扭曲严重、长宽比失控，误差就会进到弯矩里——因为弯矩本身就是应力乘以力臂的一阶矩，面积偏差会被力臂再放大一次。

2.3 符号约定与单位换算的固定写法

单位混乱是这类提取里最常见的隐性错误。FLAC3D 内部不做单位检查，你输入什么它就算什么，所以应力可能是 Pa，也可能是 kPa 或 MPa，全看你建模时怎么赋值。我的习惯是建一个固定的换算表，在脚本入口处一次性把单位锁死：

实际工程里，衬砌厚度往往是 0.3~0.8 m，模型用米制很自然；但配筋校核那边的表格是 kN、kN·m，所以脚本最后一步必须做一次除法。我吃过一次亏：模型用米、应力用 Pa，积分出来的弯矩是 N·m，直接填进配筋表当成 kN·m 用，结果差了三个数量级，配筋量算出来大得离谱，排查了大半天才发现是单位没换。从那以后我在脚本里强制写M_out = M_raw / 1e3，并且变量名一律带单位后缀，M_kNm、N_kN，看一眼就不会错。

符号方面，我统一采用“结构内力正向”约定：轴力拉为正、压为负；弯矩以使截面某侧受拉为正。由于 FLAC3D 的应力是压为负，直接积分出来的 N 在受压构件上是负值，这和结构软件的轴力符号一致，不用翻转；但弯矩的画图习惯各院不一样，我的做法是在脚本里加一个SIGN_M开关，默认取 1，出图前人工确认一次方向，确认完就固定下来，不要每次出图都临时调。

3. 用 Python 把截面内力提取跑成脚本

3.1 建一个能对得上理论解的悬臂梁基准

在把脚本往真实模型上套之前，先建一个理论解已知的基准算例，这是整套方法能不能立住的前提。选悬臂梁最合适：一根 10 m 长、0.5 m × 0.5 m 截面的弹性梁，一端固定，另一端施加竖向均布应力等效的集中力，固定端弯矩就是 P·L，轴力为零，任何位置的理论弯矩都有解析式。用这套模型跑通脚本，再去改真实模型，心里才有底。

建模命令流如下，直接贴进 FLAC3D 的命令行或者 Python 的it.command()里都能跑：

import itasca as it it.command(""" model new model title 'cantilever-section-force-benchmark' model large-strain off ; 10 m x 0.5 m x 0.5 m 悬臂梁，x 向 20 等分，y/z 各 4 等分 zone create brick size 20 4 4 ... point 0 (0,0,0) point 1 (10,0,0) ... point 2 (0,0.5,0) point 3 (0,0,0.5) zone cmodel assign elastic zone property young 3.0e10 poisson 0.2 density 2500 ; 固定端 zone gridpoint fix velocity (0,0,0) range position-x 0 ; 自由端施加竖向压应力，合力 = 1e6 Pa * 0.25 m2 = 2.5e5 N zone face apply stress-zz -1.0e6 range position-x 10 model solve ratio-local 1e-6 """)

这段命令里有几个点值得说清楚。zone create brick的point 0/1/2/3分别定义原点、x 方向端点、y 方向端点、z 方向端点，四个点定死了一个长方体，后面的size只是在这个盒子里切分。固定端用zone gridpoint fix velocity (0,0,0)把 x=0 面上所有节点的三个速度分量全部约束住，这是悬臂梁的固支条件；如果只约束法向，梁端会转动，弯矩分布就完全变了。自由端施加的是stress-zz，注意是应力不是力，合力等于应力乘以该面面积，本例是 1e6 × 0.5 × 0.5 = 2.5e5 N，算理论弯矩的时候要把这个数代进去。model solve ratio-local 1e-6把收敛判据收紧到 1e-6，弹性问题收敛很快，这个设置只是为了保证后面读应力时没有残余不平衡力干扰。

3.2 截面切片、面积还原与积分主循环

基准模型跑完之后，核心就是下面这个函数。它的思路是：给定截面位置 x0，先估出沿 x 方向的单元厚度，再挑出形心落在该层里的所有 zone，用体积除以厚度还原面积，最后套积分式。

import itasca as it import numpy as np def estimate_dx(zones): """用相邻单元形心 x 间距的中位数估计沿 x 的单元厚度""" xs = np.array([z.pos().x() for z in zones]) uniq = np.unique(np.round(xs, 6)) if len(uniq) < 2: return 1.0 return float(np.median(np.diff(uniq))) def section_forces(x0, ref_stress=None): """ x0 : 截面位置（m） ref_stress : dict{zone_id: 初始应力六分量}，用于扣除初始地应力 返回 : (N, My, Mz, A, yc, zc) """ zones = it.zone.list() dx = estimate_dx(zones) # 只取形心落在截面所在那一层的单元 layer = [z for z in zones if abs(z.pos().x() - x0) < 0.5 * dx] if not layer: raise ValueError(f"截面 x={x0} 未切到任何单元，检查网格范围") rows = [] for z in layer: s = z.stress() # (xx, yy, zz, xy, xz, yz) sxx = s[0] if ref_stress is not None: sxx -= ref_stress[z.id()][0] # 扣除初始应力 area = z.vol() / dx # 截面投影面积 rows.append((z.pos().y(), z.pos().z(), area, sxx)) arr = np.array(rows) A = arr[:, 2].sum() yc = (arr[:, 0] * arr[:, 2]).sum() / A zc = (arr[:, 1] * arr[:, 2]).sum() / A N = (arr[:, 3] * arr[:, 2]).sum() My = (arr[:, 3] * arr[:, 2] * (arr[:, 1] - zc)).sum() Mz = -(arr[:, 3] * arr[:, 2] * (arr[:, 0] - yc)).sum() return N, My, Mz, A, yc, zc

estimate_dx之所以用中位数而不是最小值，是因为网格里偶尔会有局部加密，最小值会把厚度估得偏小，导致选层时漏掉单元。section_forces里的abs(z.pos().x() - x0) < 0.5 * dx是“选一层”的判据，它的含义是：只取形心离截面半个单元厚度以内的 zone，正好构成一个单层切片。如果你的模型沿 x 方向不是等厚网格，这个判据要改成按 group 或者按层索引选，别硬套。

ref_stress这个参数是给真实模型准备的。基准模型是纯弹性、一次加载，没有初始应力，传None就行；真实岩土模型里地应力平衡完之后，每个 zone 都带着一个非零的初始应力场，这时候必须传进来做差。我一般会在平衡步结束后立刻写一个循环，把所有 zone 的 id 和应力存成字典导出，后面每次调用都读这个字典。

area = z.vol() / dx这行是整套脚本误差最大的来源。规则六面体下它几乎精确；一旦单元被压扁或者扭曲成四面体，体积和投影面积的换算关系就破了。真要高精度，得改成遍历单元的面、找到与截面共面的那个 face，取它的面积，代价是脚本复杂度上一个台阶，而且对非结构化网格的鲁棒性更差。工程上一半的精度需求，V/dx就够了。

把基准模型套进去验证：

N, My, Mz, A, yc, zc = section_forces(0.25) P = 2.5e5 # 自由端合力，N L = 10.0 M_theory = P * (L - 0.25) print(f"截面面积 A = {A:.4f} m2") # 应接近 0.25 print(f"轴力 N = {N/1e3:.2f} kN") # 应接近 0 print(f"弯矩 My = {My/1e3:.2f} kN·m") # 与 M_theory 对比 print(f"理论 M = {M_theory/1e3:.2f} kN·m")

如果 A 算出来不是 0.25 而是 0.5 或者 0.125，先查dx估得对不对，再查选层有没有选到两层。轴力应该非常接近零，如果出来一个几万牛的常数，多半是ref_stress没传或者模型有残余不平衡力。弯矩和理论值的偏差，在 20 段网格下通常能看到几个百分点，这个偏差主要不是脚本的锅，是离散本身的问题，下一章展开。

3.3 导出 CSV 与沿程弯矩包络

单个截面只能验证正确性，工程上要的是沿轴向一串截面的内力分布。把上面的函数包一层循环，逐截面扫过去，顺手导成 CSV，后面不管是给结构专业还是自己画图都方便。

import csv def scan_sections(x_start, x_end, n_sec, out_csv="section_forces.csv", ref_stress=None): xs = np.linspace(x_start, x_end, n_sec) with open(out_csv, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["x_m", "A_m2", "N_kN", "My_kNm", "Mz_kNm", "yc_m", "zc_m"]) for x in xs: try: N, My, Mz, A, yc, zc = section_forces(x, ref_stress) except ValueError: continue # 截面落在模型边界外，跳过 w.writerow([f"{x:.4f}", f"{A:.6f}", f"{N/1e3:.4f}", f"{My/1e3:.4f}", f"{Mz/1e3:.4f}", f"{yc:.6f}", f"{zc:.6f}"]) print(f"已写出 {out_csv}") scan_sections(0.25, 9.75, 39)

x_start和x_end不要取到模型的两个端点，端点那一层单元只有半边被截面穿过，V/dx会把面积算成一半，端部数据永远是错的。我的习惯是两端各留半个到一个单元厚度，0.25到9.75就是这么来的。n_sec取奇数还是偶数看你画图需要，包络图要看出极值位置，密集一点更稳，39 个截面在 10 m 长度上大约每 0.25 m 一个点，够用了。

画图用 matplotlib，两行就够：

import numpy as np import matplotlib.pyplot as plt d = np.genfromtxt("section_forces.csv", delimiter=",", names=True, encoding="utf-8") fig, ax = plt.subplots(2, 1, figsize=(7, 6), sharex=True) ax[0].plot(d["x_m"], d["N_kN"], "b-") ax[0].set_ylabel("N / kN") ax[1].plot(d["x_m"], d["My_kNm"], "r-") ax[1].set_ylabel("My / kN·m") ax[1].set_xlabel("x / m") plt.tight_layout() plt.savefig("section_force_envelope.png", dpi=150)

出图之后第一件事是拿理论解对一遍。悬臂梁的弯矩应该是一条从固定端最大值线性降到自由端零的直线，如果你画出来是锯齿状或者某一段突然跳变，别急着怀疑脚本，先看那一带的网格和材料分区。

4. 网格、截面位置、初始应力：三个决定结果可信度的参数

4.1 单元尺寸对弯矩的收敛性影响

弯矩和轴力对网格的敏感度完全不是一个量级。轴力是应力的零阶量，你把网格从 20 段加密到 40 段，N 几乎不变；弯矩是应力乘力臂的一阶矩，本质上是“截面上应力的反对称部分”，它依赖截面上的应力梯度，而梯度恰恰是有限元离散里最容易被抹平的东西。粗网格下，一个单元内部的应力被强制成常值，截面上下缘的拉压应力差被平均掉，算出来的弯矩会系统性地偏小。

我的经验是：截面尺寸方向上有 4 个单元时，弯矩误差可能在 10%~20%；加密到 8 个单元能压到 5% 上下；再往上收益递减，8 到 16 段的改进往往只有一两个百分点。所以真正要出弯矩图的时候，截面高度方向至少给 8 个单元，这是底线，不是建议。沿轴向的网格密度影响小一些，因为弯矩沿轴向变化平缓，但如果你要看的是弯矩极值位置，轴向太粗会让极值点跳到相邻截面上。

下面这张表是我在配筋校核类项目里常用的网格设置参考，具体数值要按结构尺寸缩放：

还有一点：加密网格之后弯矩变大，这是正常的，不是你算错了。很多人第一次看到加密后结果变化就去改参数，反而把对的模型调歪了。

4.2 截面位置如何避开单元界面锯齿

截面位置选在哪里，直接决定这层切片里有哪些单元。最理想的情况是截面正好落在两层单元的交界面上，此时“选一层”没有歧义；最差的情况是截面落在单元正中间，判据abs(pos.x - x0) < 0.5*dx会把前后两层都框进来一部分，切片厚度不对，面积还原出来的值也偏。

解决办法很朴素：出图之前先看一眼模型的 x 坐标列表，把要扫的截面位置取在x0 = k * dx上，k 取整数，让每个截面都落在单元界面。上面的estimate_dx已经把 dx 估出来了，scan_sections里改成按 dx 的整数倍生成 xs 就行，比等分 linspace 更稳。代价是截面位置不一定是整数米，但内力包络图对这点偏差不敏感。

另一个隐性问题是材料界面。如果截面正好穿过衬砌和围岩的交界，这一层 zone 里既有衬砌单元又有围岩单元，两者应力差一个量级，面积加权之后形心会被拉偏，弯矩直接失真。做法是在section_forces里加一个 group 过滤，只保留结构组的单元。这一点在隧道衬砌模型里是必须的，不是可选优化。

4.3 初始应力扣除与加载顺序

岩土模型和纯结构模型最大的区别就在这里。真实模型里，地应力平衡完之后，衬砌单元的应力本身就不是零，你把它当成结构内力去积分，得到的是一个巨大的假轴力。正确顺序是：地应力平衡 → 记录所有 zone 的应力快照 → 施作结构/开挖/加载 → 用当前应力减去快照 → 积分。

快照的存法有两种。简单的是在 Python 里维护一个字典，平衡步结束后遍历一次it.zone.list()存下来，模型规模在十万单元以内内存完全吃得消；规模再大就落盘成 HDF5 或者 npz，用的时候按需读。我一般写成一个独立的小函数snapshot_stress()，返回{zone_id: (sxx, syy, szz, sxy, sxz, syz)}，加载步之后传给section_forces的ref_stress参数。

有个细节要留意：如果加载过程中网格发生了大变形（比如用了model large-strain on），zone 的 id 一般不变，但形心位置会动，快照里的坐标就不能再用了。这种项目里我只对 id 做差，坐标一律取当前值，力臂也用当前形心算，逻辑上更自洽。

5. 提取结果对不上理论值时的排查顺序（避坑）

5.1 现象：轴力是个巨大的常数，弯矩曲线形状反而正常

原因：初始地应力没有扣除，或者扣了但扣的是不同计算步的快照。表现是整条 N 曲线被抬到一个很高的水平，但因为初始应力沿截面近似均匀，它对弯矩的一阶矩贡献很小，所以弯矩看着还行，轴力彻底不能用。

解决：确认快照是在哪一步存的，ref_stress传进去的是不是同一个字典。可以在section_forces里临时加一行打印，把扣除前后的 N 都打出来，一对比就知道有没有扣上。另外注意，如果模型中途改过本构或者重置过应力，快照要重存。

5.2 现象：网格加密之后弯矩变大，和“加密应该更准”的直觉相反

原因：粗网格下截面上的应力梯度被单元常值应力抹平，弯矩系统性偏小；加密之后梯度被还原出来，弯矩自然上升。这是离散误差的正常表现，不是脚本出了问题。

解决：做一次网格收敛性检查，把截面高度方向的单元数按 4、8、12 各算一遍，看弯矩是否趋于稳定值，取稳定段的值作为报告值。如果加密到 12 还在明显上涨，说明你的模型里存在应力集中或者截面附近有约束干扰，得先把这些因素排除，而不是继续加密。

5.3 现象：相邻两个截面的弯矩值跳变明显，曲线呈台阶状

原因：截面位置没有落在单元界面上，落在单元中间时，判据把前后两层各框进来一部分，或者只框进一层，导致相邻截面的“切片厚度”不一致。另一种可能是截面穿过了材料分区，一层里混了两种模量的单元。

解决：先把扫描位置改成 dx 的整数倍，让每个截面都卡在单元界面上，台阶通常会消失。如果还有，检查这一层单元的 group，把非结构组过滤掉，或者在截面筛选条件里加z.isgroup("liner")之类的判断。

5.4 现象：同一个模型重复跑两次，同一截面的结果对不上

原因：求解没有收敛到足够精度，残余不平衡力还在让应力缓慢变化；或者模型中存在未固定的刚体位移，每次求解停在不同位置。

解决：把model solve ratio-local收紧到 1e-6 或者 1e-7，重新算一次看结果是否稳定；同时用model unbalanced查一下最大不平衡力，如果量级和施加的荷载相当，说明根本没收敛。另外确认边界条件里没有漏固定方向的节点，尤其是用实体单元拼出来的组合结构，接缝处最容易漏约束。

5.5 现象：弯矩符号和结构专业给的图完全相反

原因：FLAC3D 的应力符号是压负拉正，积分式里的力臂方向也依赖坐标系约定，和结构力学的弯矩正负规则不是一套。

解决：不要试图在公式里“修”符号，那会让代码变得没法读。在脚本最后统一乘一个SIGN_M系数，取 1 或 -1，用图形和已知理论解对比一次确认，然后固定下来写进脚本头部注释。弯矩图的正负本身不影响力矩的绝对值，配筋看的是量值。

6. 进阶：从单截面提取做成沿轴扫描的自动包络

到这一步，单截面提取已经能用了，但真正在项目里跑起来，你会发现手工指定截面位置很累：一个衬砌模型几百米长，你想知道弯矩极值出现在哪，靠人眼在云图上找是不靠谱的。我后来把这套东西改成了“自动包络”，思路很简单：沿轴向按 dx 步长全扫一遍，每一截面都算 N、My、Mz，然后直接找极值位置和量值，输出一张带极值标注的表。

def auto_envelope(x_min, x_max, dx, ref_stress=None): """ 沿轴向按 dx 步长扫描，返回内力极值及其位置 """ x_list = np.arange(x_min, x_max + 0.5 * dx, dx) recs = [] for x in x_list: try: N, My, Mz, A, yc, zc = section_forces(x, ref_stress) except ValueError: continue recs.append((x, N, My, Mz)) arr = np.array(recs) if arr.size == 0: raise RuntimeError("扫描区间内没有有效截面") idx_N = np.argmax(np.abs(arr[:, 1])) idx_My = np.argmax(np.abs(arr[:, 2])) idx_Mz = np.argmax(np.abs(arr[:, 3])) print(f"N 极值 {arr[idx_N, 1]/1e3:9.2f} kN @ x = {arr[idx_N, 0]:.3f} m") print(f"My 极值 {arr[idx_My,2]/1e3:9.2f} kN·m @ x = {arr[idx_My,0]:.3f} m") print(f"Mz 极值 {arr[idx_Mz,3]/1e3:9.2f} kN·m @ x = {arr[idx_Mz,0]:.3f} m") return arr

dx直接传estimate_dx(it.zone.list())的返回值，保证扫描点全部落在单元界面上，这个问题在 4.2 节已经说过，是台阶状曲线的根源。x_min和x_max记得两头各留一个单元厚度，别扫到边界层。返回值arr是四列，可以顺手喂给前面的画图代码。

这套自动包络真正值钱的地方不在于省了手工，而在于它把“弯矩极值位置”从定性变成了定量。实体单元模型的弯矩极值位置，经常和你从结构力学直觉上判断的位置差个一两米，原因是围岩约束、材料刚度变化、施工步序都会把极值点拽走。用扫描出来的位置去回看该处的应力云图和网格，往往能发现一些建模阶段没注意到的问题，比如某个材料分区边界正好压在极值附近，或者某段网格长宽比失控。

要提醒的是，自动包络不是免费的。每扫一个截面都要遍历一次it.zone.list()，十万单元、一百个截面，跑起来是分钟级的，别放在每次调试里都执行，缓存成 CSV 增量更新更实际。另外，如果你用的本构里有软化或者损伤，内力沿程分布可能在加载后期发生重分布，包络要按施工步分别出，不能只出最终步。

我自己在这套流程上踩得最深的一脚，不是公式写错，是有一版脚本忘了把截面筛选限制在结构组里，结果围岩单元的应力混进来，弯矩图看着挺光滑、数值也“合理”，直到拿去和结构单元结果做交叉验证才发现整体偏了一个系数。从那以后我定了个习惯：任何新模型第一次出内力图，必须先和至少一种独立方法对一遍，理论解、结构单元、手算截面，随便哪一种，对上了再往下做。希望帮到你。

本文还有配套的精品资源，点击获取