简介：一份面向ANSYS APDL与Matlab用户的有限元矩阵后处理脚本。脚本针对结构工程与科研分析场景，可直接接收APDL导出的文本格式刚度矩阵和质量矩阵，在Matlab中完成数据读取、矩阵维度重构与数据整理，并提供特征值分解示例，便于读者进一步开展模态分析、频率响应和结构动力学研究，适合需要将有限元模型导出至Matlab做二次开发的工程师与科研人员。压缩包仅包含一个脚本文件，大小仅有六百一十二字节，内容精简、步骤完整，可直接放入工作目录运行，也可按需修改自由度参数适应不同模型规模。已有三千九百一十二人学习，使用热度较高，说明它能够有效解决APDL与Matlab之间的数据衔接问题。借助这段代码，读者可以跳过繁琐的手工赋值环节，快速获得可计算的刚度、质量矩阵，并在此基础之上扩展绘图展示与振动特性分析，显著提升后处理效率与建模分析灵活性。
1. 从黑匣子到 Matlab：为什么要把 ANSYS APDL 的 K、M 矩阵导出来
做结构动力学分析的时候，我经常遇到一个需求：ANSYS 的模态结果只是一个黑匣子，模态频率和振型都有了，但做轴承稳定性判断、模型降阶、传感器优化布置或者动刚度修正时，我真正需要的是刚度矩阵 K 和质量矩阵 M。先把 K、M 从有限元模型里完整取出来，再交给 Matlab 做二次计算，这条路比在 APDL 里绕来绕去要快得多。ANSYS APDL 输出有限元模型刚度矩阵和质量矩阵，配合 Matlab 后处理代码，解决的正是这个需求：把求解器手里的全局矩阵变成你自己能操作的数据。适合做转子动力学、振动控制、模型修正和联合仿真的工程师，也适合想把 ANSYS 的结果从黑匣子变成白盒子的所有人。
2. APDL 侧导出：HBMAT 命令与最小可跑通脚本
2.1 为什么选 HBMAT：三条导矩阵路径怎么选
要拿到全局刚度矩阵和质量矩阵，APDL 里常见的做法有三条：HBMAT 命令、WRFULL 写满文件、子结构法输出。三条路径我都试过，最终长期使用的是 HBMAT，原因很简单：它直接输出文本格式的稀疏矩阵，Matlab 读起来最不折腾。
WRFULL 输出的是二进制 .full 文件，里面包含矩阵和自由度信息，但解析二进制格式要自己对字节流，工作量大且容易翻车。子结构法也能生成矩阵文件，但它的设计初衷是超单元减缩，会顺带做自由度凝聚，如果你不想被凝聚逻辑干扰，配置反而更啰嗦。HBMAT 则直接在求解完成后把 K 或 M 写成一个 Harwell-Boeing 格式的文本文件，也就是行业里常说的 HB 格式。这种格式是稀疏矩阵的标准交换格式，Matlab 侧只需要解析三个数组段就能恢复完整矩阵，逻辑非常直观。
HBMAT 命令在我常用的写法里是这样构成的：
| 参数 | 作用 | 我常用的设置 |
|---|---|---|
| 文件名 | 输出文件的主名 | KMAT、MMAT 这类见名知意的名字 |
| 扩展名 | 文件扩展名，可留空 | 留空，后面自己补 .HB |
| 导出对象 | STIFF 或 MASS | 刚度矩阵写 STIFF，质量矩阵写 MASS |
| 矩阵格式 | HB 或 CCM | 选 HB，Matlab 解析最简单 |
| 存储模式 | 0 只存下三角，2 存完整对称矩阵 | 默认 0，省空间 |
| 打印选项 | 是否回显矩阵 | 小模型用 YES，大模型用 NO |
关键点在于，HBMAT 必须在完成一次 SOLVE 之后执行。矩阵是求解过程中组装好的，没求解就调用，拿不到全局矩阵。
2.2 最小可跑通脚本：先求解、再导出 K 和 M
我用一个最简单的两单元梁模型来演示完整流程，这套脚本可以直接抄去改模型。单位制要提前统一，这里用的是毫米、牛、吨这套常用组合，密度对应 7.85e-9 t/mm³。
FINISH /CLEAR, NOSTART /PREP7 ET,1,BEAM188 MP,EX,1,2.1e5 MP,PRXY,1,0.3 MP,DENS,1,7.85e-9 SECTYPE,1,BEAM,RECT SECDATA,50,50 N,1,0,0,0 N,2,1000,0,0 N,3,2000,0,0 E,1,2 E,2,3 FINISH /SOLU D,1,ALL,0 SOLVE FINISH /SOLU HBMAT,KMAT,,,STIFF,HB,0,,YES HBMAT,MMAT,,,MASS,HB,0,,YES FINISH先做一次静力求解，再把刚度矩阵写入 KMAT.HB，把质量矩阵写入 MMAT.HB。模型本身没有任何载荷，SOLVE 只是在形式上完成一次组装，矩阵已经在这个阶段生成。D,1,ALL,0 把节点 1 完全约束住，这很重要，因为 HBMAT 导出的是参与求解的活动自由度对应的矩阵，被约束掉的自由度不会出现在输出文件里。
不是所有单元都一定能导出质量矩阵。如果你用的单元没有密度，或者单元类型本身不带质量，MMAT.HB 里面就会是空矩阵或直接报错。这也是为什么我在脚本里提前定义 MP,DENS。做完这一步，工作目录下应该出现两个 HB 文件，文件大小通常远小于全量稠密矩阵，因为存的是稀疏数据。
2.3 导出文件头部长什么样：HB 格式的直观印象
拿到 KMAT.HB 之后，直接用文本编辑器打开，或者用命令行 head 看前几行，能直观看到 HB 格式的布局。我用 bash 看一眼：
head -n 8 KMAT.HB看到的典型内容大致是这样：
Harwell-Boeing 文件: ANSYS 导出的稀疏矩阵 STIFFNESS MATRIX FROM ANSYS (10I6) 8 8 14 0 0 1 3 6 10 15 ...前几行是注释和格式说明，第 4 行的三个关键数字依次是：矩阵行数、矩阵列数、非零元素个数。后面的整数数组是列指针，再后面是行索引，最后是数值数组。理解这个布局是 Matlab 解析的前提，代码逻辑也完全围绕这三段数据展开。
值得注意，ANSYS 默认输出的是对称矩阵的下三角部分，非零元个数约等于满矩阵非零数的一半多一点，而不是你直觉里那个 nnz。如果你把输出的 nnz 和模型总自由度数直接对比，会发现数量对不上，这是正常的，因为只有下三角被写进去了。
3. Matlab 侧读取：不靠第三方，用 30 行代码解析 HB 文本
3.1 HB 文件的布局：三个数组段
HB 格式虽然带了一个看起来挺正式的名字，但落到文件层面就三件事：列指针数组、行索引数组、数值数组。这三个数组配合矩阵阶数，就能还原出完整的稀疏矩阵。
列指针数组的长度是 列数 + 1，它的作用是把非零元按列分组。假设第 j 列的非零元位于列指针数组的第 j 和第 j+1 个元素之间，那么这一列就对应行索引数组和数值数组中同一区间的内容。行索引数组给出每个非零元的行号，数值数组给出对应位置的值。
我整理成一张表方便对照：
| 数组段 | 长度 | 含义 |
|---|---|---|
| 列指针 | ncol+1 | 每列非零元在后续两个数组中的起始位置 |
| 行索引 | nnz | 每个非零元所在行号 |
| 数值 | nnz | 每个非零元的实际数值 |
在 ANSYS 导出的 HB 文件里，列指针通常是 Fortran 风格的 1 基索引，也就是第一个指针通常是 1。如果用 0 基索引去解析，矩阵第一列会直接丢数据。这个细节我踩过，后面避坑章节会专门说。
3.2 read_hb 解析函数：返回稀疏 K、M 的代码
下面这个函数是我现在一直在用的简化版，只处理 ANSYS 最常见的实对称矩阵导出场景。文件头注释里写清了入参和出参，直接复制就能用。
function [K, n] = read_hb(fname, sym_flag) % read_hb 读取 ANSYS APDL HBMAT 导出的 HB 文本文件 % 输入: % fname : 文件名，如 'KMAT.HB' % sym_flag : 1 表示只存了下三角(默认)，0 表示完整矩阵 % 输出: % K : Matlab 稀疏矩阵 % n : 矩阵阶数 if nargin < 2 sym_flag = 1; end fid = fopen(fname, 'rt'); if fid < 0 error('无法打开: %s', fname); end % 跳过前 4 行说明行 fgetl(fid); % 注释 fgetl(fid); % 标题 fgetl(fid); % 格式字符串 line4 = fgetl(fid); % 维度行 dims = sscanf(line4, '%d'); n = dims(1); % 行数 ncol = dims(2); % 列数 nnz = dims(3); % 非零元个数 % 读取列指针，长度为 ncol+1 col_ptr = fscanf(fid, '%d', ncol+1); % 读取行索引 row_ind = fscanf(fid, '%d', nnz); % 读取数值 val = fscanf(fid, '%f', nnz); fclose(fid); % 由列指针还原每个非零元的列号 col_ind = zeros(nnz, 1); for j = 1:ncol seg = col_ptr(j) : col_ptr(j+1)-1; col_ind(seg) = j; end K = sparse(row_ind, col_ind, val, n, n); % 如果文件只存下三角，补全为完整对称矩阵 if sym_flag K = K + K.' - diag(diag(K)); end end这段代码的核心逻辑就是先读维度行，再依次读三个数组，最后用 sparse 组装。最容易被忽略的是 seg 区间计算，col_ptr(j) 到 col_ptr(j+1)-1 这一段才属于第 j 列。如果这里把端点开闭搞反，矩阵的非零元会整体错位。
补全对称矩阵那行用了 K + K.' - diag(diag(K))，因为对角元加了两次，要减回去一次。如果原始文件已经存了完整矩阵，sym_flag 传 0 就不要补全，否则对角元会放大一倍。
3.3 参数说明：sym_flag 与存储模式如何对应
调用这个函数的时候，sym_flag 的设置要和 APDL 里 HBMAT 命令的存储模式对齐。我在 2.1 表格里提过，HBMAT 的存储模式参数选 0 时只存下三角，这时 Matlab 侧必须传 sym_flag=1。如果 APDL 侧选了完整矩阵存储，Matlab 侧就传 0。
实际使用中我建议 APDL 侧保持默认的只存下三角模式，文件小，读写快，Matlab 补全只需要一行代码，没有理由在 APDL 侧让文件体积翻倍。
调用方式很直接：
K = read_hb('KMAT.HB', 1); M = read_hb('MMAT.HB', 1); % 检验维度是否一致 assert(size(K,1) == size(M,1), 'K 和 M 阶数不一致');拿到 K 和 M 之后，你的自由度顺序问题就浮出水面了。矩阵本身是稀疏的，但顺序是 ANSYS 内部的自由度顺序，不是你在前处理里建模时的节点顺序。这正是下一章要处理的事情。
4. 做自由度排序映射：让矩阵行号对上物理节点
4.1 为什么矩阵行号是“乱序”的
把 KMAT.HB 读进 Matlab 后，第一件事不是算特征值，而是搞清楚矩阵的第 1 行到底对应哪个物理自由度。ANSYS 在求解器内部会对自由度重新编号，这个过程对用户不可见，所以你看到的矩阵行号只是内部方程号，不等于节点号也不需要等于节点号。
这个排序黑匣子，是做矩阵二次开发时最容易踩坑的地方。你如果直接把矩阵行号当成节点号去画振型，画出来的结果大概率是杂乱的。在做模态综合或 Guyan 减缩时，我一般会定义一组主自由度，用 M 命令把关心的节点自由度挑出来。这种情况下，HBMAT 导出的矩阵会以这组主自由度的顺序作为行列顺序，这个顺序通常就是你定义主自由度时的顺序。
所以我的习惯是：在 APDL 里做矩阵导出的同时，把主自由度顺序同步记到一张文本表里。矩阵是谁的顺序，表里就记录谁。后面 Matlab 侧做映射时，以这个表为准，而不是猜测。
4.2 在 APDL 里把自由度顺序写进一张映射表
定义主自由度时，顺手把节点号和方向写进文件。最可靠的方式不是事后用命令从模型里反查，而是在定义的时候就记录。下面是配合 HBMAT 使用的主自由度映射表输出脚本：
! 选择主自由度 M,2,UX M,2,UY M,2,ROTZ M,3,UX M,3,UY M,3,ROTZ ! 把主自由度顺序写入 DOF_MAP.TXT *CFOPEN,'DOF_MAP','TXT' ! 每一行: 节点号, 方向编号(1=UX,2=UY,6=ROTZ) *VWRITE,2,1 (F8.0,1X,I2) *VWRITE,2,2 (F8.0,1X,I2) *VWRITE,2,6 (F8.0,1X,I2) *VWRITE,3,1 (F8.0,1X,I2) *VWRITE,3,2 (F8.0,1X,I2) *VWRITE,3,6 (F8.0,1X,I2) *CFCLOS这个表必须和你定义主自由度的顺序完全一致，顺序错一个，整个映射就全部错位。方向编号里 1 代表 UX，2 代表 UY，3 代表 UZ，4 代表 ROTX，5 代表 ROTY，6 代表 ROTZ，这是 ANSYS 的方向标识法，不要在编号上玩创意。
如果你直接用全部活动自由度导出而不用主自由度，映射表的生成会更复杂，因为自由度数量是节点数乘方向上自由度数的结果。我的建议是：导出矩阵的场景，尽量配合主自由度集使用，既缩小矩阵规模，也让顺序可控。这个做法在模型缩聚里天然成立，一般的 K、M 提取也够用。
4.3 Matlab 里完成行号到节点的定位
拿到 DOF_MAP.TXT 之后，Matlab 侧就能建立矩阵行号和物理自由度之间的对应关系。用 readtable 读这个文本，然后给每个自由度算一个全局编号。
% 读取自由度映射表 tp = readtable('DOF_MAP.TXT', 'ReadVariableNames', false); node = tp.Var1; dof = tp.Var2; % 方向编号转换: 1->UX, 2->UY, 6->ROTZ dof_name = {'UX', 'UY', 'UZ', 'ROTX', 'ROTY', 'ROTZ'}; % 构造全局自由度编号, 用于和 K 矩阵行号对齐 n_dof_total = 6; % 每个节点按 6 个自由度考虑 global_id = (node - 1) * n_dof_total + dof; fprintf('K 矩阵第 1 行对应: 节点 %d, 方向 %s\n', ... node(1), dof_name{dof(1)});这套映射有一个前提：你的模型里所有节点自由度编号都是按 6 个自由度连续排列的。对于大多数实体单元和梁单元模型，这个前提成立。壳单元或某些缩减单元可能只有部分自由度，此时要把 n_dof_total 改成实际参与的自由度类型数。
5. 矩阵导出与后处理避坑：5 条实际踩过的记录
5.1 质量矩阵导出为 0 或提示不可用
现象：HBMAT 导出质量矩阵时，得到的 MMAT.HB 里全是零，或者命令直接提示质量矩阵不可用。
原因：单元没有密度。MP,DENS 没定义，或者定义了但单元类型本身不产生质量矩阵。另一个可能的原因是还没有执行 SOLVE，全局质量矩阵还没组装。
解决：检查 APDL 里是否定义了 MP,DENS，单元类型是否支持质量矩阵。先执行一次 SOLVE 再调用 HBMAT，顺序不能反。
5.2 读出来矩阵行列数跟预期差很多
现象：Matlab 读出的矩阵阶数，和节点数乘自由度算出来的结果完全对不上。
原因：HBMAT 导出的是活动自由度矩阵，被约束掉的自由度不会进入输出。你如果按全部节点自由度去算规模，自然就偏大。另外，如果模型里还有耦合自由度或主从约束，参与求解的自由度数会进一步变化。
解决：用约束施加后的活动自由度数量作为预期值。可以通过 APDL 里的 *GET 命令拿到实际自由度数量，再和读出的 n 比对。
5.3 矩阵行号顺序前后对不上
现象：同一次模型，先后导出了两次 K、M，Matlab 里按相同映射去取某一行的值，两次结果不对应同一个自由度。
原因：两次导出之间，网格顺序、主自由度定义顺序或者约束施加顺序被改动过，内部编号就变了。如果两个矩阵本身来自两次不同的求解，顺序不一致是必然的。
解决：K 和 M 必须在同一次求解流程里连续导出，导出的同时把映射表也写出来。不要复用上次分析产生的 DOF_MAP.TXT，每次都以当前导出的映射表为准。
5.4 全稠密存储导致内存直接打爆
现象：Matlab 读取 K 矩阵时，内存占用瞬间飙升，甚至直接卡死退出。
原因：用了 full(K) 把稀疏矩阵转成稠密矩阵，或者在某些操作中无意间触发了稠密化。HB 文件本身是稀疏存储，几千自由度规模的模型转成稠密矩阵后，内存需求会按自由度的平方增长。
解决：全程保持稀疏矩阵形态。后续特征值计算用 eigs，不要用 eig。如果必须做矩阵加减，确保参与运算的另一个矩阵也是 sparse。我对这个问题的处理原则是：任何一步操作前，先确认 whos 里变量类型是 sparse。
5.5 eig 算出的频率数量级离谱
现象：用 K 和 M 算出的一阶频率，和 ANSYS 模态分析结果差了上千倍。
原因：单位制不一致。如果 APDL 用的是毫米、牛、吨单位制，模型里 E 的单位是 MPa，密度是 t/mm³，长度是 mm。直接把这些矩阵拿到国际单位制下计算，K 和 M 的比例关系就会失真，频率自然不对。
解决：在 APDL 侧统一单位，并在 Matlab 侧做换算。我常用的换算方式是：长度 mm 转 m 时，K 除以 1e6，M 除以 1e3，再求特征值。换算完成后和 ANSYS Modal 模块的结果对比，通常能对上。
6. 用 eigs 做一次自检：拿固有频率验证 K、M 没读错
6.1 从矩阵到频率的完整代码
矩阵读进来了，映射表也建立好了，接下来必须验证一件事：读出来的 K 和 M 没有错位。最简单可靠的方法是用 eigs 计算固有频率，和 ANSYS 模态分析结果做对比。这是一个自检动作，却能省掉后面大量无用功。
% 读取 K 和 M K = read_hb('KMAT.HB', 1); M = read_hb('MMAT.HB', 1); % 单位换算: 假设 APDL 使用 mm/N/t 单位制 K_si = K / 1e6; % N/mm -> N/m M_si = M / 1e3; % t -> kg % 计算前 10 阶特征值 opts.tol = 1e-10; [D, V] = eigs(K_si, M_si, 10, 'smallestabs', opts); % 固有频率，单位 Hz f_hz = sqrt(real(diag(D))) / (2 * pi); % 排序输出 [f_hz_sorted, idx] = sort(f_hz); disp(f_hz_sorted);eigs 的 'smallestabs' 选项求的是幅值最小的特征值，正好对应最低几阶模态。tol 设置到 1e-10 是为了让特征值精度足够，否则后面对比时误差可能来自求解器而不是矩阵本身。
6.2 与 ANSYS 模态结果对比的排查顺序
把算出来的频率和 ANSYS Modal 分析结果放在一起对比，允许的偏差应该在小数点后第四位以内。如果偏差明显，按这个顺序排查：先查单位换算，再查自由度映射表，最后查求解器设置。
我自己的习惯是：拿到 K、M 后，第一件事不是做降阶或优化，而是先把一阶频率算出来和 ANSYS 对一次。对上了，矩阵就是可信的；对不上，后面所有后处理都是在错误的地基上盖楼。
这条流程我重复过很多次，后来养成了固定习惯：APDL 脚本里永远把 HBMAT 和映射表输出放在同一个求解流程里，Matlab 侧永远先用频率验证再往下走。这套做法让我少走了太多弯路，希望帮到你。
本文还有配套的精品资源，点击获取
