企业官网建设流程全解析

做岩土数值模拟的朋友，十有八九都遇到过这个场景：FLAC3D模型辛辛苦苦算收敛，云图也出了几张，想把局部区域的主应力方向单独拎出来分析，结果发现自带的显示要么全模型密密麻麻，要么按钮找了半天。这篇文章就专门解决这个需求——用FLAC3D 7.0里的FISH脚本，把指定区域（比如某隧洞洞周）的应力张量分量导出成文件，再用MATLAB读进去，求解主应力方向并画出能直接放进报告或论文的图。

我先把话说在前面：这套流程不是银弹，但它覆盖了“局部区域筛选—数据导出—外部求主应力—定制可视化”这条很常见的后处理链路。适合正在做岩土、采矿、地下工程数值模拟的同学，也适合那些虽然用FLAC3D但不想被自带绘图限制住的人。你不需要很懂FISH，跟着流程抄一遍，改一改坐标范围和输出路径，基本就能用。

1. 为什么要把主应力方向导出来自己画

1.1 自带可视化“能看”但不“够用”

FLAC3D本身有Contour of Principal Stress，可以直接看到三个主应力的大小和方向线，日常检查模型确实方便。但真到出成果阶段，它的短板就暴露了。

第一，方向线密度没法精细控制。模型几万、几十万个zone，图里密密麻麻全是穿越符号和十字线，你想看洞室右上角那一片区域的主应力偏转，结果图上全是“毛线”。第二，FLAC3D的交互式操作支持range过滤，你可以通过界面把显示范围缩小到局部，但每一次看图都要操作一遍，没法批量产出几十个剖面的图。第三，关键是它不给你“数据”。云图只能看个大概，论文里需要统计主应力的倾角、方位角、随深度的变化规律，这些没法直接从图上量出来。

所以就有了这篇文章想解决的问题：把我们关心的那部分zone，用脚本挨个取出来，把坐标和应力张量分量写到文件里，再用MATLAB计算主应力方向并画图。这一步做完，后面无论是统计分析、方位角玫瑰图、还是跟现场实测地应力做对比，都顺畅了。

1.2 这需求到底来自哪：局部区域才是重点

有人说“FLAC3D不是能直接导出zone应力吗？何必绕一圈”。这句话对，但问题在于“部分区域”。完整模型几十万个zone，全导出来你也难处理；而且我们关注的点通常很局部。举个我实际碰到的例子，某输水隧洞开挖后，要分析洞周2倍洞径范围内的应力重分布，想看清楚最大主应力方向在洞壁两侧怎么旋转。如果用全模型导出的数据画图，洞外远处占主导的区域把颜色和向量长度直接“淹没”掉，局部细节根本看不出。所以必须先在FISH里按坐标范围、或者按到某轴线的距离筛选出那部分zone，再导出。

类似的场景还有很多：采矿工作面顶板重点关注垮落区上方一定范围；边坡稳定分析要看滑动面附近主应力方向的偏转；断层带附近要研究应力方向是否受断层影响。甚至有人需要在某个截面上取出所有zone，画一张主应力方向剖面图。核心需求是一致的：局部取数、批量处理、外部绘图。所以我建议的思路是，FISH负责“筛选和导出”，MATLAB负责“求方向和可视化”，各干各擅长的事。

2. 动手之前，先把应力数据和FISH函数摸清楚

2.1 应力张量、主应力、主方向的关系

想做主应力方向，对张量基础得稍微回顾一下。岩土单元体上的应力状态，在三维里是一个3×3对称张量，因为有对称性，独立分量只有6个：sxx、syy、szz、sxy、sxz、syz。FLAC3D里zone的应力就是按这6个分量存的。主应力是什么？就是这个3×3矩阵的特征值，对应的特征向量就是主应力方向。数学上很直接，MATLAB里一行eig就能解决。

但工程上有两个容易掉坑的点。第一，特征向量代表的“方向”只有方向没有正负，eig返回的向量可能指向相反方向，需要自己做符号规范化。第二，符号约定。FLAC3D默认以拉应力为正、压应力为负，所以它显示的σ1是代数值最大的主应力（最“拉”的那个）。很多国内工程习惯压为正，这会导致排序结果完全反过来。你不能拿着MATLAB按默认排序的结果直接跟现场地应力报告里的“最大水平主应力”对应，要先想清楚自己要用哪套约定。好消息是，应力张量整体乘以−1后，特征向量方向不变，只是特征值符号翻转，所以主应力方向本身不受影响。

2.2 FLAC3D 7.0里FISH的几个常用“取数”函数

FLAC3D 7.0的FISH已经是很现代的脚本语言了，遍历zone不用背一堆缩写，常见的几个函数要记牢：

• zone.list：得到模型中zone链表的头指针，配合z.next(zp)可以遍历全部zone。
• z.pos(zp)：返回zone质心位置的Vector，分量用->x(1)、->x(2)、->x(3)取；老一点的习惯是直接z.x(zp)、z.y(zp)、z.z(zp)。
• z.sxx(zp)、z.syy(zp)、z.szz(zp)、z.sxy(zp)、z.sxz(zp)、z.syz(zp)：读取6个应力分量。
• z.sig1(zp)、z.sig2(zp)、z.sig3(zp)：直接读取三个主应力数值，这个可以拿来和MATLAB计算结果交叉验证，非常有用。

还有file有关的一组函数：file.open、file.write、file.close，用来写文本文件。FISH写文件比较原始，你得自己拼字符串加换行，但胜在无需额外接口，模型在哪算，文件就出在哪。

注意：7.0同时保留新旧两套FISH函数写法，我下面用的是新式点函数。如果你在自己版本里遇到“unknown function”之类的报错，把z.sxx(zp)换成旧式z_sxx(zp)，坐标函数换成z_x(zp)这种形式，基本就能解决。我写这套流程时也踩过这个坑，后来统一用新式函数就顺畅了。

3. FISH脚本：筛选部分区域并导出应力张量分量

3.1 脚本总体思路

脚本要做的事很简单：遍历zone，判断质心坐标是否在指定范围内，在的话就读取6个应力分量和坐标，按行写进CSV文件。再加一个计数器，导完以后在控制台输出一共导出了多少个zone，方便你确认筛选条件没写错。

这里有个经验：筛选条件越早判断越好。我见过有人先写文件抬头、再把所有zone都写一遍，最后在MATLAB里面筛，这样也不是不行，但几十万zone的文件读写时间和文件体积都很浪费。在FISH循环里用一条if把绝大多数zone挡掉，效率完全不同。

3.2 可直接用的FISH脚本

下面这段脚本是按坐标范围筛选的，比如要导出x在20到60米、y在20到40米、z在−30到−10米范围内的所有zone：

fish define export_stress local zp = zone.list local fname = 'selected_stress.csv' local f = file.open(fname, 'write') if f = null then io.out('Cannot open file: ' + fname) return endif file.write(f, 'x,y,z,sxx,syy,szz,sxy,sxz,syz' + '\n') local count = 0 loop while zp # null local pos = z.pos(zp) local cx = pos->x(1) local cy = pos->x(2) local cz = pos->x(3) if cx >= 20.0 and cx <= 60.0 and cy >= 20.0 and cy <= 40.0 and cz >= -30.0 and cz <= -10.0 then local sxx = z.sxx(zp) local syy = z.syy(zp) local szz = z.szz(zp) local sxy = z.sxy(zp) local sxz = z.sxz(zp) local syz = z.syz(zp) file.write(f, string(cx) + ',' + string(cy) + ',' + string(cz) + ',' + string(sxx) + ',' + string(syy) + ',' + string(szz) + ',' + string(sxy) + ',' + string(sxz) + ',' + string(syz) + '\n') count = count + 1 endif zp = z.next(zp) endloop file.close(f) io.out('Exported zones: ' + string(count)) end @export_stress

写完后在FLAC3D控制台执行@export_stress，如果一切正常，会看到类似Exported zones: 15620的输出，工作目录下生成selected_stress.csv。

这是最基础的框架，但已经够用了。我一般会把坐标范围、文件名都抽出来放到脚本开头，方便一个模型用不同区域反复导出。还可以在这个框架上扩展出更复杂的筛选逻辑：比如按到某条隧洞轴线的距离小于阈值来筛、按zone是否进入塑性状态来筛，或者按某组group来筛。

提示：如果担心浮点比较误差，比如筛选边界恰好落在zone质心上，建议写成cx > 19.999 and cx < 60.001这种留一点余量的形式，避免理想化的等号判断把边界上的zone漏掉。

3.3 输出文件格式设计

输出文件我习惯用CSV，表头一行，后面每行一个zone，逗号分隔。为什么用逗号不用空格？因为MATLAB的readmatrix和readtable读取CSV最省事，空格分隔会出现字符串错位。

格式如下：

x,y,z,sxx,syy,szz,sxy,sxz,syz 23.415,30.002,-11.873,12345.6,-4567.8,...

这里有几个细节。

第一，字符串拼接时用string()把数字转成文本，FISH默认打印精度比较低，如果你需要保留足够多有效数字，建议在string()之前用round或者自己写一个格式函数，或者干脆在MATLAB里用str2double再处理，从CSV读取时精度一般不受影响。第二，文件开到几百兆以上要留意磁盘和内存，我通常会在循环里加一个if count >= 500000 then exitloop endif的保险，防止写太多文件卡死。第三，如果模型包含group，也可以在输出列里加上group名称，后续按组统计很方便。

4. MATLAB读取数据：从6个分量到主应力方向

4.1 读取CSV并组装应力张量

MATLAB读取CSV最简单的方式是readmatrix，表头自动跳过。拿到数据后，把6列应力分量拼成每个zone的3×3对称矩阵。某些初学者会在这里把矩阵拼错，尤其是sxy、sxz、syz的位置，一定要和CSV的列顺序对上。

我用下面的代码，假设CSV列顺序是x,y,z,sxx,syy,szz,sxy,sxz,syz：

data = readmatrix('selected_stress.csv'); x = data(:,1); y = data(:,2); z = data(:,3); sxx = data(:,4); syy = data(:,5); szz = data(:,6); sxy = data(:,7); sxz = data(:,8); syz = data(:,9); n = numel(x); sig = zeros(n,3); dir1 = zeros(n,3); dir2 = zeros(n,3); dir3 = zeros(n,3); for i = 1:n S = [sxx(i), sxy(i), sxz(i); sxy(i), syy(i), syz(i); sxz(i), syz(i), szz(i)]; [V, D] = eig(S); d = diag(D); [~, idx] = sort(d, 'descend'); sig(i,:) = d(idx)'; dir1(i,:) = V(:,idx(1))'; dir2(i,:) = V(:,idx(2))'; dir3(i,:) = V(:,idx(3))'; end

这段是核心。eig返回的V每一列是对应的特征向量，D的对角元素是与V列顺序一致的特征值。用sort(diag(D),'descend')把特征值从大到小排，得到的就是σ1、σ2、σ3以及它们的方向向量。

4.2 排序与符号约定：最容易错的地方

刚跑通时很多人会发现，自己用MATLAB算出来的σ1、σ2、σ3和FLAC3D里显示的数值对不上。原因十有八九是符号约定。FLAC3D默认拉应力为正，它的σ1是“代数最大”，而如果你在MATLAB里习惯先取绝对值再排序，结果自然不一样。

做法是，先明确自己项目里到底按哪套约定输出。如果你要跟岩石力学教材里“压为正”的公式对比，直接在组装矩阵时先把应力整体变号，也就是S = -[sxx ...]，特征向量不变，但特征值排序会变成压应力越大越靠前。另一种办法是不变号，排序后自己记得σ1是拉正约定下的最大主应力。两种都行，但在报告里必须写清楚，不然别人看你图里的正负号会糊涂。

符号规范化是第二个坑。eig返回的特征向量在数学上方向可以任意变号，同一个zone，可能第一次算[0.81, -0.32, 0.49]，换个机器变成[-0.81, 0.32, -0.49]，画出来的箭头朝上朝下完全相反。解决方法是给一个统一规则，比如“让每个方向向量的z分量大于0，否则反向”，或者按某一参考方向来定向。对于剖面图，我习惯让σ3的竖直分量指向z正方向，这样能直观看到巷道底板处的最小主应力方向往下指。

for i = 1:n if dir1(i,3) < 0, dir1(i,:) = -dir1(i,:); end if dir2(i,3) < 0, dir2(i,:) = -dir2(i,:); end if dir3(i,3) < 0, dir3(i,:) = -dir3(i,:); end end

这段代码就是干这个的，简单粗暴但有效。不过要注意，如果某个方向向量本来就应该指向z负方向（比如你想画矿柱内竖向主应力向下），强行规范会让所有箭头朝上，画图表达反了。我在做洞周主应力分析时一般先画出原始箭头，确认规律后再加规范逻辑，不要一开始就无脑统一符号。

4.3 用FLAC3D数值交叉验证一下

拿到MATLAB结果后，建议先抽几个zone做交叉验证。在FLAC3D里选一个zone，比如用z.sig1(zp)读出它的主应力值，再和你MATLAB计算表的对应行对比。方向也验证：在FLAC3D的plot里打开主应力方向显示，找同一个zone，看它的σ1方向线与MATLAB画出来的箭头是否大致重合。

这一步看似多余，实际能帮你提前发现坐标顺序、矩阵拼接、上下颠倒等各类问题。我第一版脚本就在这里栽过：CSV表头里把sxz和syz写反了，导致MATLAB画出的方向线在局部区域明显不对劲，交叉验证半小时才发现是列顺序错了。所以千万不要跳过验证，尤其是要往报告里放的图。

5. 画图：从能出图到能表达规律

5.1 剖面十字线图：快速看清主应力方向旋转

最常用的图是剖面图。比如隧道沿y方向开挖，我们要看xoz平面内的主应力方向，就需要筛掉所有不在剖面附近的zone，取abs(y - 30) < 0.5这样的近似平面。剖面图最适合用箭头表达方向，线长表示大小。

在MATLAB里画最大主应力σ1方向，可以用quiver：

sel = abs(y - 30) < 0.5; figure('Color', 'w'); hold on; quiver(x(sel), z(sel), dir1(sel,1), dir1(sel,3), 0.4, 'r', 'LineWidth', 1.2); quiver(x(sel), z(sel), dir3(sel,1), dir3(sel,3), 0.4, 'b', 'LineWidth', 1.0); axis equal; xlabel('x (m)'); ylabel('z (m)'); legend('\sigma_1', '\sigma_3');

这里quiver的第五个参数是箭头缩放因子，自己调整到线不重叠、长度合适为止。axis equal一定要加，否则x和z比例不一致，方向箭头画出来会被压扁，主应力方向的倾角全歪了。

更专业一点的做法是用十字线：在每个zone质心位置同时画σ1和σ3两条短线，形成类似“应力十字”的图案。这样既看到方向，又对比两个方向的大小差距。可以手动用line或者plot实现，不用quiver，因为quiver的箭头有时候太长，大小比例不好控制。

5.2 三维箭头图：全貌和局部兼得

如果只看二维剖面，会丢掉另一个主应力方向的信息。三维场景我建议画σ1方向的quiver3，点位置用zone质心，方向向量用dir1乘以应力大小比例，这样能直观看到最大主应力从远处到开挖边界附近是怎么偏转的。

figure('Color', 'w'); scale = 5e-7; % 根据实际应力量级调整 quiver3(x, y, z, dir1(:,1) .* sig(:,1) * scale, dir1(:,2) .* sig(:,1) * scale, dir1(:,3) .* sig(:,1) * scale, 0, 'r'); axis equal; view(45, 30); xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');

注意quiver3的箭头长度受坐标系范围影响很大，三维箭头经常看起来要么太短要么穿模。我后来更习惯用line把每条方向线画成一段线段：

figure; hold on; for i = 1:n L = 1.5; % 线段长度 v = dir1(i,:) * L; line([x(i) - v(1), x(i) + v(1)], ... [y(i) - v(2), y(i) + v(2)], ... [z(i) - v(3), z(i) + v(3)], 'Color', 'r'); end view(3); axis equal;

这种画法虽然慢，但方向线的长度、粗细、虚实都能完全控制，画出来像一条条应力迹线，非常适合报告用。

5.3 云图与方向线叠加，信息量翻倍

很多工程问题要同时看应力大小和方向。一个很实用的做法是用σ1的数值给点着色，再用方向线叠上去。最简单的实现是scatter，颜色表示σ1大小，向量箭头表示方向：

figure; hold on; scatter3(x, y, z, 30, sig(:,1), 'filled'); quiver3(x, y, z, dir1(:,1), dir1(:,2), dir1(:,3), 0.3, 'k'); axis equal; view(3); colormap(jet); colorbar;

如果数据点太密，可以先对σ1进行griddata插值到规则网格，再画contourf或pcolor作背景，最后叠加方向箭头。插值会光滑掉一些局部细节，所以优先用scatter直接显示离散点，等图面清楚后再考虑网格化。关键是云图和箭头用同一个坐标系，不要为了好看把箭头缩放得跟云图尺度不一致。

6. 常见问题与避坑记录

6.1 一写出来就报错，或导出文件为空

跑FISH脚本最常见的两个问题：函数不认识、筛选范围没选到zone。函数不认识多半是版本差异，按前面说的方法换旧式写法就能解决。筛选范围没选到zone，先检查单位。FLAC3D里几何尺寸按米建模，如果你脑子里想的范围是毫米、厘米，那大概率全筛空。

还有一个隐蔽问题，CSV文件是相对路径，可能写到了FLAC3D的工作目录而不是你MATLAB脚本所在目录。找文件找不到时，去模型文件所在目录找一下，几乎每次都在。最笨也最有效的排查办法是，在脚本里先不加筛选，导出前10个zone的坐标，在控制台打印出来，确认坐标提取逻辑没问题，再去加筛选条件。

6.2 单zone验证都对，但整体图画出来杂乱无章

这是方向没有规范化的典型表现。不同zone的特征向量可能在相邻位置发生正负突变，画出来就是一堆朝不同方向的乱线。重点是，这不一定是你算错了，可能是特征向量的方向符号不连续。解决办法就是前面说的符号规范化，但要注意规范化规则要连续。比如让每个方向向量的某个分量（如x分量）整体为正，如果刚好有主应力方向接近竖直，x分量接近0，符号翻转就会导致图中出现不自然的突变。稳妥的做法是先算一个参考方向，让每个点的方向向量与参考方向的夹角小于90度，小于90度就保持，大于90度就反向。

我一般在分析前先定参考：如果是隧洞周边，σ1的方向应该大致沿洞室环向；如果是地应力场，σ1通常在水平面和竖直面之间偏转。这个参考值一确定，所有箭头方向就统一了，图中也看不出剪切突变。

6.3 数据量太大，MATLAB绘图卡死怎么办

几十万个zone全导出来，MATLAB画散点、画线都会很卡。我的建议是第一层筛选放到FISH里做，只输出关注区域；第二层在MATLAB里再做抽稀，比如每隔几个点取一个，或者按剖面厚度合并。但抽稀会改变视觉密度，报告里要注明。还有一个替代思路：在FISH里直接按应力值阈值筛选，比如只输出σ1超过某个值的zone，这也是“部分区域”的另一种用法。

如果图面点还是太多，我会把方向线缩短，或者只画σ1方向不画σ3，让信息量降下来。图画出来是为了说明问题，不是为了堆数据，能表达核心规律就好。

6.4 输出文件里的科学计数法导致MATLAB读取异常

FISH默认写数字可能输出成科学计数法，用readmatrix一般能正确解析，但如果你用readtable后转数值，偶尔会碰到被识别成字符的情况。稳妥起见，导出时用string()后最好自己格式化一下。FISH里没有特别方便的printf到字符串的函数，我的土办法是先按round(x,4)控制精度，或者干脆保留原始精度，在MATLAB里统一str2double再处理。

最后再分享一个小技巧：输出文件命名时带上模型和时间戳，比如stress_zoneA_20240115.csv，并在FISH脚本开头用注释写清楚筛选范围、单位、模型名。这样后期整理数据时，不用打开文件才知道这组数据是哪来的。我自己的习惯是每个工况导一次，最终的MATLAB图全部加上工况标识，对比分析时会省很多事。

我在实际项目里用这套流程处理过某隧洞洞周应力重分布，一开始也是从FLAC3D自带云图里猜方向，后来导出数据在MATLAB里做统计，才真正看出σ1在洞壁两侧大概旋转了20多度。这个结论在报告里的说服力，比一句“方向发生偏转”强太多了。希望这套办法也能帮到你，先用小模型把这个链路跑通，再上大模型，问题会少很多。