企业官网建设流程全解析

简介：这份资源面向从事岩土工程数值模拟的技术人员与科研工作者，聚焦FLAC3D实体单元无法直接输出弯矩、轴力的痛点，系统讲解基于应力积分的截面内力提取方法。内容涵盖FISH函数实现思路、隧道初期支护与地铁支护桩等实战案例，并延伸至单位换算、截面选取、网格划分等易错细节，以及结果验证与精度校核的实用技巧。资源包为1个PDF文件，约104KB，内容紧凑，适合作为随查随用的技术参考。目前已有135人学习。读者可从中获得一套可复用的内力提取代码框架、典型工程场景的受力分析思路，以及判断计算结果是否可靠的验证方法，对提升数值模拟精度、支撑论文与工程报告中的数据论证具有实际帮助。

1. 为什么实体单元里“读”不出弯矩轴力，却还要硬读

用 FLAC3D 做地下结构或桩基分析的人，几乎都遇到过这个尴尬：模型跑完了，想看衬砌或桩的弯矩图，结果软件里翻遍 history、zone extra、结构单元输出，就是找不到一个叫“弯矩”的字段。原因不复杂——FLAC3D 的实体单元（zone）本质是连续介质，它算出来的是应力张量，不是内力。弯矩、轴力这些是截面内力，是应力沿截面二次积分的结果，软件不会替你自动做这一步。

所以“实体单元弯矩轴力提取”这件事，本质是自己动手把应力场换算成截面内力。它解决的是：当你用实体单元模拟梁、板、衬砌、桩这类构件时，如何从应力结果反推出工程上真正关心的弯矩和轴力。适合谁？适合已经能跑通 FLAC3D 静力分析、但卡在“结果怎么用”这一步的从业者。这篇笔记就按我实际做过的路径，把原理、命令、参数和踩过的坑一次讲清楚。

2. 从应力到内力：实体单元提取弯矩轴力的底层逻辑

2.1 实体单元没有内力，只有应力张量

FLAC3D 的 zone 每个单元存储的是应力分量，三维情况下是六个独立分量：σxx、σyy、σzz、σxy、σyz、σxz。这些是材料点的应力，不是截面上的合力。弯矩和轴力的定义是：

• 轴力 N = ∫ σ dA，沿截面法向积分
• 弯矩 M = ∫ σ · y dA，y 是到中性轴的距离

也就是说，你必须先定义一个截面，再在这个截面上对应力做面积分。FLAC3D 不会自动帮你定义截面，更不会自动积分。常见做法是：在构件上切一个截面，取该截面附近一层 zone 的应力，按面积加权求和。

这里有个容易翻车的地方：zone 的应力是存储在单元中心的，如果你直接取一个单元的值乘以面积，精度取决于网格密度。网格太粗，弯矩误差能到 30% 以上。我一般会要求截面附近至少 4 到 6 层单元，否则后面积分出来的弯矩图会像锯齿一样。

2.2 截面内力积分：离散求和代替连续积分

连续积分在 FLAC3D 里没法直接做，必须离散化。假设截面沿 y 方向，中性轴在 y0，那么：

• 轴力 N = Σ σxx_i · A_i
• 弯矩 M = Σ σxx_i · A_i · (y_i - y0)

其中 A_i 是第 i 个单元在截面上的面积，y_i 是单元中心到中性轴的距离。σxx_i 是截面法向应力。如果是梁单元，法向就是轴向；如果是板或衬砌，法向是厚度方向。

这个求和看起来简单，但有两个参数必须提前定好：中性轴位置 y0 和截面面积 A_i。中性轴一般取几何中心，但如果材料非线性或者截面不对称，中性轴会偏移。我通常先按几何中心算一遍，再根据轴力结果反推中性轴——如果 N 不为零而外荷载没有轴力，说明中性轴取错了。

2.3 用 FISH 遍历 zone 并输出应力与坐标

FLAC3D 里做这件事最直接的工具是 FISH。下面这段代码是我常用的骨架：遍历指定范围内的 zone，取出应力分量和中心坐标，写到文件里。

; 定义截面范围：x 在 [x1, x2]，z 在 [z1, z2]，沿 y 方向切截面 def extract_section x1 = 0.0 x2 = 1.0 z1 = -0.5 z2 = 0.5 y_target = 5.0 ; 截面所在的 y 坐标 y_tol = 0.6 ; 允许的 y 方向容差，取一层单元 fp = open('section_stress.txt', 1, 1) local z = zone_head loop while z # null ; 取单元中心坐标 xc = z.x(z) yc = z.y(z) zc = z.z(z) ; 判断是否在截面范围内 if xc > x1 and xc < x2 and zc > z1 and zc < z2 and abs(yc - y_target) < y_tol ; 取应力分量，FLAC3D 中 z.sxx(z) 返回单元应力 sxx = z.sxx(z) syy = z.syy(z) szz = z.szz(z) ; 取单元体积，用于估算截面面积 vol = z.vol(z) ; 输出：x, y, z, sxx, vol status = write(fp, string(xc) + ' ' + string(yc) + ' ' + string(zc) + ' ' + string(sxx) + ' ' + string(vol)) endif z = z.next(z) endloop status = close(fp) end @extract_section

逻辑说明：这段 FISH 遍历所有 zone，用中心坐标判断是否落在截面附近。y_tol控制取几层单元，一般取 0.5 到 1 倍单元尺寸。z.sxx(z)取的是单元应力，FLAC3D 默认输出的是总应力，如果用了小应变模式，注意区分有效应力和总应力。z.vol(z)是单元体积，后面用来估算截面面积——对于规则六面体，截面面积约等于体积除以单元在 y 方向的厚度。

参数说明：x1, x2, z1, z2定义截面在 x-z 平面的范围，必须覆盖整个构件截面。y_target是截面位置，y_tol太大会混入非截面单元，太小可能取不到足够单元。我一般先跑一遍看输出行数，如果少于 20 行，说明网格太粗或容差太小。

2.4 把应力积分成弯矩轴力的后处理脚本

拿到截面应力数据后，下一步是积分。我习惯用 Python 做后处理，因为 FISH 的数组操作太别扭。下面是一个最小可复现的脚本：

import numpy as np # 读取 FISH 输出的截面数据 data = np.loadtxt('section_stress.txt') x = data[:, 0] y = data[:, 1] z = data[:, 2] sxx = data[:, 3] vol = data[:, 4] # 估算单元在 y 方向的厚度，假设截面沿 y 方向切 dy = 0.1 # 需要根据实际网格尺寸修改 # 截面面积 = 体积 / dy area = vol / dy # 中性轴位置：先按几何中心 y0 = np.mean(z) # 注意这里 z 是截面高度方向 # 轴力 N = np.sum(sxx * area) # 弯矩：力臂为 z - y0 M = np.sum(sxx * area * (z - y0)) print(f'轴力 N = {N:.3f}') print(f'弯矩 M = {M:.3f}')

逻辑说明：area = vol / dy是近似，前提是单元在 y 方向厚度均匀。如果网格不均匀，最好直接用单元在截面上的投影面积，而不是用体积反推。y0取 z 的均值，对应几何中心。弯矩的力臂是z - y0，因为截面高度方向是 z。

参数说明：dy必须和实际网格尺寸一致，否则面积算错，轴力和弯矩都会成比例偏大或偏小。我一般会在 FISH 里直接输出单元在 y 方向的尺寸，而不是事后猜。另外，如果截面法向不是 x，而是 y 或 z，需要把sxx换成对应的应力分量，力臂方向也要跟着换。

3. 实战：在衬砌和桩基模型里提取内力

3.1 衬砌截面：沿厚度方向积分

衬砌是典型的薄壁构件，用实体单元模拟时，厚度方向通常只有 2 到 4 层单元。这种情况下，弯矩对网格非常敏感。我的做法是：在衬砌厚度方向至少布置 4 层单元，然后沿厚度方向积分。

假设衬砌厚度方向是 z，截面法向是 x，那么：

• 轴力 N = Σ σxx_i · A_i
• 弯矩 M = Σ σxx_i · A_i · (z_i - z0)

z0 是衬砌厚度中心。如果衬砌内外侧应力差很大，弯矩就大；如果应力均匀，弯矩接近零。这里有个血泪经验：衬砌的轴力往往很大，弯矩相对小，如果网格太粗，弯矩会被轴力的数值误差淹没。我一般会先把轴力算准，再算弯矩，两者分开验证。

3.2 桩基截面：圆形截面的离散处理

桩基是圆形截面，用实体单元模拟时，截面是六边形或八边形近似。这时候积分不能简单按矩形算，必须按实际单元面积加权。

常见做法是：在桩身切一个水平截面，取截面附近一层 zone，每个单元的面积用vol / dz估算，dz 是单元在竖向的厚度。然后：

• 轴力 N = Σ σzz_i · A_i
• 弯矩 Mx = Σ σzz_i · A_i · (y_i - y0)
• 弯矩 My = Σ σzz_i · A_i · (x_i - x0)

x0, y0 是桩截面中心。注意圆形截面的中性轴通过圆心，但离散后单元中心不一定对称，所以 x0, y0 要用面积加权平均算，而不是简单取坐标均值。

# 圆形截面桩基内力提取 x0 = np.sum(x * area) / np.sum(area) y0 = np.sum(y * area) / np.sum(area) N = np.sum(szz * area) Mx = np.sum(szz * area * (y - y0)) My = np.sum(szz * area * (x - x0))

逻辑说明：x0, y0用面积加权，是为了消除离散误差。如果直接用np.mean(x)，当单元分布不均匀时，中心会偏。Mx和My分别是绕 x 轴和 y 轴的弯矩，对应桩基的两个方向。

参数说明：szz是竖向应力，对应桩的轴向。如果桩是斜桩，需要把应力张量旋转到桩轴方向，这一步容易出错，我一般建议斜桩直接用结构单元，不要用实体单元硬做。

3.3 用 history 记录关键截面内力随加载步的变化

单次提取只能得到最终状态的内力。如果想看内力随加载过程的变化，需要在每一步都提取。FLAC3D 的 history 可以记录 zone 应力，但不能直接记录积分后的内力。我的做法是：用 FISH 写一个自定义 history，每一步调用积分函数，把结果存到 history 里。

; 定义 history 变量 def m_history m_val = 0.0 n_val = 0.0 ; 调用积分函数，返回 M 和 N m_val = calc_moment n_val = calc_axial ; 存入 history history m_val type history history n_val type history end

逻辑说明：calc_moment和calc_axial是前面积分逻辑的封装。每次 step 之后调用m_history，就能记录内力时程。注意 FISH 的 history 数量有限，不要每个单元都记，只记关键截面。

参数说明：history的采样间隔由history interval控制，默认每 10 步记一次。如果加载很快，可以改成 1。但记太密会让结果文件很大，我一般设 5 到 10。

4. 避坑：实体单元提取内力的五个常见翻车点

4.1 网格太粗，弯矩像锯齿

现象：提取出来的弯矩沿构件长度方向跳动很大，相邻截面弯矩差几倍。原因：截面附近单元太少，应力积分误差大。解决：截面附近至少 4 层单元，最好 6 层。如果已经跑完不想重算，可以用插值平滑，但那是遮丑，不是解决问题。

4.2 中性轴取错，轴力不为零

现象：外荷载没有轴力，但提取出来的 N 明显不为零。原因：中性轴位置取错，或者截面法向应力分量取错。解决：先检查应力分量是否对应截面法向，再用面积加权算中性轴。如果 N 仍然不为零，检查边界条件是否有约束反力。

4.3 应力分量搞混，弯矩符号反了

现象：弯矩图正负号和预期相反。原因：FLAC3D 的应力符号约定和材料力学不同，拉应力为正还是压应力为正，容易搞混。解决：先做一个简单悬臂梁算例，手算弯矩和 FLAC3D 提取结果对比，确认符号约定。我一般会在脚本里加一个符号开关，方便调整。

4.4 单元体积反推面积，网格不均匀时误差大

现象：轴力和弯矩整体偏大或偏小。原因：用vol / dy估算面积，当单元在 y 方向厚度不均匀时，面积算错。解决：直接在 FISH 里输出单元在截面法向的尺寸，或者用z.vol(z)除以实际厚度。更稳妥的做法是用单元在截面上的投影面积，但 FLAC3D 不直接提供，需要自己算。

4.5 动力分析中内力振荡，提取结果不可用

现象：动力时程分析中，提取的弯矩高频振荡。原因：应力波在单元之间传播，单点应力波动大。解决：在截面上做空间平均，或者用低通滤波。我一般会在截面上取多层单元平均，而不是单层。如果振荡仍然严重，说明网格太粗，需要细化。

5. 进阶：用 Python 批量后处理与内力图绘制

单截面提取只能得到一个点，实际工程需要沿构件长度画内力图。我的做法是：在 FISH 里沿构件长度方向切多个截面，每个截面输出一个数据文件，然后用 Python 批量读取、积分、绘图。

import numpy as np import matplotlib.pyplot as plt # 截面位置列表 sections = [1.0, 2.0, 3.0, 4.0, 5.0] N_list = [] M_list = [] for y_sec in sections: # 读取对应截面的数据文件 data = np.loadtxt(f'section_{y_sec:.1f}.txt') x = data[:, 0] z = data[:, 2] sxx = data[:, 3] vol = data[:, 4] dy = 0.1 area = vol / dy y0 = np.sum(z * area) / np.sum(area) N = np.sum(sxx * area) M = np.sum(sxx * area * (z - y0)) N_list.append(N) M_list.append(M) # 绘制内力图 fig, ax1 = plt.subplots() ax1.plot(sections, M_list, 'b-o', label='弯矩 M') ax1.set_xlabel('沿构件长度方向位置 (m)') ax1.set_ylabel('弯矩 (kN·m)', color='b') ax2 = ax1.twinx() ax2.plot(sections, N_list, 'r-s', label='轴力 N') ax2.set_ylabel('轴力 (kN)', color='r') plt.title('实体单元提取的内力图') plt.show()

逻辑说明：这段脚本循环读取每个截面的数据文件，分别积分得到 N 和 M，然后画双轴图。弯矩和轴力量级差很大，所以用双 y 轴。如果弯矩图有突变，检查截面位置是否落在集中荷载或支座附近。

参数说明：sections列表是截面位置，间隔根据构件长度定，一般 0.5 到 1 米一个点。dy必须和 FISH 里输出数据时的假设一致。如果截面法向不是 x，把sxx换成对应分量，力臂方向也跟着换。

一个具体技巧：如果构件是变截面的，每个截面的面积不同，积分时要用实际面积，不能用一个常数。我一般会在 FISH 里把每个截面的单元面积单独输出，Python 里直接读，避免反推误差。

最后说个我自己的教训：早期做衬砌分析时，为了省时间，网格划得很粗，结果弯矩图完全不能用，重新划网格跑了两天。后来我养成习惯，凡是需要提取内力的构件，截面附近网格至少加密一倍，宁可多算几小时，也不要返工。希望帮到你。

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