企业官网建设流程全解析

简介：面向ANSYS有限元初学者的刚度矩阵专题资料包，围绕结构静力、动力分析中的核心概念，覆盖了从刚阵定义、单元刚阵形成到全局组装、矩阵提取与后处理应用的完整链路。包内共39个文件，以MATLAB的.m程序为主，附带Fortran源程序、C文件、文本说明、动态库及可执行程序等，合计仅380KB，既可阅读理论文档也可直接运行脚本验证计算流程。已有358人学习下载，适合需要结合代码理解刚度矩阵原理的工程师与研究生使用。内容包含ANSYS Workbench提取刚度矩阵的操作指引、框架结构刚阵组装脚本、坐标变换、静力凝聚与动力响应等模块，能帮助读者加深对Ku=F方程、边界条件及网格质量影响的理解，是一份轻量实用的入门参考资料。

1. Stiffness Matrix 在 ANSYS 里到底是个什么角色

做结构分析的工程师大多有过这种体验：在 Workbench 里点了 Solve，几秒钟后看到云图，但中间发生了什么却像黑箱。真正把有限元当工具而不是当神的人，都会去碰那个核心的矩阵——刚度矩阵。它本质上就是线性系统 (K u = F) 里的 (K)，所有静力学、模态、谐响应分析的解都是从这一个矩阵出发的。ANSYS 里你设置的单元类型、材料弹性模量、截面尺寸、边界条件，最后都会被折算成这个矩阵里的系数。它可以很大，比如一个十万节点的模型，(K) 就是十万阶的稀疏矩阵，但它的组装逻辑和你在课本上学到的三自由度弹簧系统没有任何区别。这套资源里的 km_form.m、SpaceFrameAssemble 脚本、matrixout.f90，刚好把这套逻辑从理论到代码串了一遍。适合正在做二次开发、想导出 ANSYS 矩阵做减缩或者想自己写单元的人。

2. 单元刚度矩阵的形成：km_form.m 与坐标变换矩阵

2.1 有限元里“单元贡献”是怎么落地的

整体刚度矩阵不是凭空出现的，它由每个单元的局部刚度矩阵 (k^e) 经过坐标变换后，按照自由度编号“投递”到全局矩阵中。局部坐标下的单元刚度矩阵只和单元几何、材料属性有关，比如一个平面梁单元，在局部坐标系下是四阶或六阶矩阵，里面只有 (E)、(I)、(A)、(L) 这些参数。但实际结构中单元朝向千变万化，必须在全局坐标系里组装。这时用到的就是坐标变换矩阵 (T)，局部矩阵到全局矩阵的映射关系是 (K^e = T^T k^e T)。资源里的 example_coordinate_transformation.m 就是专门做这个事的。

2.2 km_form.m 里做了什么

km_form.m 这个名字一看就是“k matrix formation”，它负责根据节点坐标和材料参数生成单元刚度矩阵。以空间框架单元为例，每个节点有 6 个自由度（三个平动、三个转动），单元局部刚度矩阵是 12×12。常见做法是先算出轴向、扭转、两个平面内弯曲的刚度分量，再按自由度顺序填入矩阵。伪代码逻辑是：

function k_local = km_form(E, A, Iy, Iz, G, J, L) % 空间梁单元局部刚度矩阵，自由度顺序: [ux v w rx ry rz] x2 % 轴向刚度 k_axial = E * A / L * [1 -1; -1 1]; % 扭转刚度 k_torsion = G * J / L * [1 -1; -1 1]; % 绕z轴的弯曲刚度（x-y平面内） k_bend_z = E * Iz / L^3 * [12 6*L -12 6*L; ... 6*L 4*L^2 -6*L 2*L^2; ... -12 -6*L 12 -6*L; ... 6*L 2*L^2 -6*L 4*L^2]; % 然后按自由度顺序散开 k_local = zeros(12, 12); % ... 按节点1的6个自由度和节点2的6个自由度填块 end

这段代码里最关键的是自由度顺序约定。有的程序按[u v w rx ry rz]排，有的按[ux uy uz rotx roty rotz]排，顺序错了整个矩阵就是错的。你把 k_local 打印出来，对角线元素应该都是正数，而且每行之和为零（刚体位移模式下内力为零，这是检验单元矩阵的最低标准）。

参数说明：E是弹性模量，A是截面积，Iy和Iz是两个主轴惯性矩，G是剪切模量，J是扭转常数，L是单元长度。空间梁单元如果不考虑剪切变形和翘曲，这六个参数就能定死局部矩阵。如果你用的是欧拉-伯努利梁而不是铁木辛柯梁，程序里通常就没有剪切修正系数，这一点看代码里的12*E*Iz/L^3前面的系数就能判断——若是12*E*Iz/(L^3*(1+phi))，那就考虑了剪切变形。

2.3 坐标变换的常见坑位

实际工程里，单元局部坐标系的 x 轴通常沿单元轴向，但 y 轴和 z 轴的方向需要人为指定主方向，否则变换矩阵不唯一。ANSYS 里你设置梁的截面方向时，其实就是在定这个主方向。example_coordinate_transformation.m 里一般用一个三点定义法：给定单元两个端点坐标和第三个参考点，来构造局部坐标系的方向余弦矩阵。用方向余弦拼出 (T)（12×12 块对角矩阵）后，K_global_e = T' * k_local * T这一步在 MATLAB 里写起来很简单，但很多新手会漏掉转置。(T) 是正交矩阵，理论上 (T^{-1}=T^T)，但如果你在构造 (T) 时因为方向余弦没有归一化导致不正交，那么T'和inv(T)就不一样，结果自然出错。

我一般做完坐标变换后会做两个自检：一是把变换前后的矩阵特征值对比，坐标变换不应该改变特征值（因为 (T^T k T) 是相似变换）；二是计算整个结构的刚体位移模态，应该得到六个接近零的频率或零特征值。如果特征值差得远，十有八九是方向余弦矩阵算错了。

3. 全局刚度矩阵组装：SpaceFrameAssemble 与自由度编号

3.1 组装不是矩阵加法，是按编号投递

单元矩阵算完之后，要把它们组装成全局矩阵 (K)。这个过程的本质是“直接刚度法”：每个单元的两个节点在全局自由度列表里都有对应的全局编号，单元矩阵里的每个元素 (k_{ij}) 要加到全局矩阵的 (K_{row, col}) 上，其中row和col是由单元节点自由度映射出来的全局行、列号。资源里的 SpaceFrameAssemble 这个 MATLAB 程序，注释里写得很清楚，它就是做这种投递的。它的典型输入是单元节点连接矩阵elems、节点坐标nodes和单元刚度矩阵的 cell 数组。

function K = SpaceFrameAssemble(node_coord, elem_node, elem_mat) % node_coord: 节点坐标矩阵，每行一个节点 [x y z] % elem_node: 单元连接定义，每行两个节点编号 % elem_mat: 单元局部或全局刚度矩阵 cell 数组 nn = size(node_coord, 1); % 节点数 ndof = 6; % 每个节点自由度 K = zeros(nn*ndof); % 预分配全局矩阵 for e = 1:size(elem_node, 1) n1 = elem_node(e, 1); n2 = elem_node(e, 2); dof1 = (n1-1)*ndof + (1:ndof); dof2 = (n2-1)*ndof + (1:ndof); dof_index = [dof1, dof2]; K(dof_index, dof_index) = K(dof_index, dof_index) + elem_mat{e}; end end

这里的核心参数是ndof。对于平面刚架，ndof 取 3；对于空间刚架，ndof 取 6；对于桁架，ndof 取 2 或 3。如果你把 ndof 取错，矩阵维度直接对不上，MATLAB 会立刻报错。但维度对上了不代表编号正确——常见错误是单元节点顺序反了，导致单元矩阵的“第一个节点”和“第二个节点”互换，而单元矩阵本身是分块对称的，如果两个节点自由度块完全一样（比如对称截面），结果恰好不报错，但实际上自由度映射错位，全局矩阵会有物理上说不通的耦合。

3.2 半带宽优化和稀疏存储

直接生成一个 (N \times N) 的满矩阵，在节点数超过几千时就会耗光内存。实际 ANSYS 内部用的是稀疏矩阵存储，只保存非零元素。MATLAB 里把K = zeros(...)改成K = sparse(...)就能省掉大量内存。而带宽优化则是一个经典问题：节点编号顺序直接影响矩阵的轮廓大小，编号不好的矩阵半带宽大，求解慢。ANSYS 的 Wavefront 求解器和现在默认的稀疏直接求解器都会自动重排自由度来减小带宽。你自己写程序时，可以用 MATLAB 里的symrcm或amd函数对全局自由度编号做重排。

% 组装完成后，用 Cuthill-McKee 算法减小带宽 p = symrcm(K); K_banded = K(p, p);

不过要注意，重排之后自由度编号变了，后续施加边界条件和读取结果时，位移向量也要按照p做对应的映射。否则你解出来的u(p)和真实的节点位移对不上。我一般会把映射向量存下来，并在注释里写明原始自由度编号和重排后编号的对应关系。

3.3 边界条件：去掉奇异，而不是置零

组装好的 (K) 在没有施加边界条件时是奇异的，因为结构存在刚体位移自由度。常见做法是把固定约束自由度对应的行和列删掉，形成缩小后的 (K_{red}) 和载荷向量 (F_{red})。另一种做法是惩罚法，在约束自由度上叠加一个很大的数，比如 (K_{ii} + 10^{12})，但惩罚法的精度取决于你这个数取得多大，取小了约束不严，取大了引入数值病态。资源里的 exam3_2.m 里多半就是用的缩减自由度法，因为它更接近有限元教材里的标准流程。关键坑在于自由度编号是 1 到 (N) 的连续整数，如果用 MATLAB 的setdiff筛选，别搞混节点编号和自由度编号。比如固定节点 5 的所有平动自由度，就需要把(5-1)*6+1: (5-1)*6+3这些自由度编号从集合里删掉。

4. 从 ANSYS 里把刚度矩阵提取出来：Static Condensation 与文件交换

4.1 ANSYS 导出矩阵的两种路线

很多人以为 ANSYS 里看不到刚度矩阵，其实经典 ANSYS（Mechanical APDL）可以直接用HBMAT命令导出刚度矩阵、质量矩阵和阻尼矩阵到外部文件。Workbench 则要借助 Solution 里的 Output Controls 或者用命令对象插入一段 APDL 代码。常见做法是：在 Workbench 的 Model 下插入 Commands（APDL），写一段：

/SOLU ! 计算刚度矩阵并写入文件 WRFULL,1 HBMAT,'K_matrix','txt',' ', '', 'K', 1, YES, YES

HBMAT后面第一个参数是文件名，第二个是扩展名，第三个是路径，第四个留空，第五个是矩阵类型标识（K表示刚度矩阵），第六个是格式（1 表示 text，2 表示 binary），第七个YES表示是否用 Harwell-Boeing 格式写出，第八个YES表示是否写出排序信息。导出的文件是文本格式，但它是 Harwell-Boeing 稀疏格式，不是普通矩阵的样子。要把它读进 MATLAB，需要自己写解析函数，或者用资源里的 matrixout.f90 去转换。

4.2 Fortran 程序 matrixout.f90 的转换逻辑

资源里的 matrixout.f90 和 BINLIB.LIB 应该是从 ANSYS 子结构分析二次开发里来的老代码。它的功能大概率是把 ANSYS 输出的二进制或 Harwell-Boeing 文件读出来，再转成 Fortran 顺序存储的稠密矩阵。Harwell-Boeing 格式的行列索引从 1 开始，且非零元素指针数组比实际非零数多一个末尾标记，读的时候很容易偏移一位。如果你要自己写读取器，核心逻辑是：

! 读取 HBMAT 导出的 .txt 文件（text 格式） read(unit, '(A)') tail_line ! 跳过前 4 行头部信息 ! 第 5 行开始是列指针, 然后行索引, 然后数值

空间有限，实际你更推荐在 MATLAB 里写。网上也有现成的hb_read函数，但如果你不想依赖外部包，可以用下面的思路：前四行分别是文件头、格式描述、矩阵维度、非零数，之后的行数是列指针，再后面是行号，最后是实部数值。需要注意的是 ANSYS 导出的矩阵可能是对称的，只保存上三角或下三角，读取后要K = K + K' - diag(diag(K))补全。

4.3 Static Condensation 有什么用

Static Condensation（静力凝聚）在资源里出现了两次，一个文档一个 m 文件，说明这包东西不只是提取矩阵那么简单。静力凝聚的核心思想是消去不需要的内部自由度，只保留边界自由度或主自由度。对于线性静力问题，把节点分为主自由度 (m) 和从自由度 (s)，刚度矩阵分块为：

[ \begin{bmatrix} K_{mm} & K_{ms} \ K_{sm} & K_{ss} \end{bmatrix} \begin{bmatrix} u_m \ u_s \end{bmatrix}

\begin{bmatrix} F_m \ F_s \end{bmatrix} ]

当从自由度上没有外力时，(u_s = -K_{ss}^{-1} K_{sm} u_m)，代入后可得到减缩后的刚度矩阵：

[ K_{cond} = K_{mm} - K_{ms} K_{ss}^{-1} K_{sm} ]

这个公式在 MATLAB 里的实现非常直接，但数值上要小心：如果 (K_{ss}) 的条件数很差，求逆会放大误差。常见做法是对 (K_{ss}) 做 Cholesky 分解再回代，而不是直接inv(K_ss) * K_sm。资源里的 static condensation.m 我猜测就是用 chol 或反斜杠运算符做的。从 ANSYS 提取出完整矩阵后，再做一次静力凝聚，就能把模型自由度减到几百个，用于后续子结构分析或试验相关性分析。

5. 用提取出来的矩阵做模态缩减与参数验证

5.1 组装修正后的质量矩阵和阻尼矩阵

刚度矩阵提出来之后，光有它无法做动力学分析，还需要质量矩阵。ANSYS 可以同时导出质量矩阵，HBMAT命令将矩阵类型改为M即可。在 MATLAB 里，你可以用提取的 (K) 和 (M) 求解广义特征值问题来验证模型：求 (\det(K - \omega^2 M)=0) 的根。如果导出的矩阵是减缩的，自由度数量对不上，特征值结果可能偏移。一个快速验证方法：把导出的 (K) 用于一个只有三个自由度的悬臂梁模型，手算前两阶频率，再和 ANSYS 的模态分析结果对比，误差应该在 1% 以内。超过这个范围，先检查边界条件是否在矩阵里体现——ANSYS 导出的矩阵有时是未施加边界条件的自由矩阵，需要你自己删除约束自由度。

5.2 灵敏度分析和矩阵扰动检查

在做结构优化时，刚度矩阵的偏导是不可或缺的量。你不需要去解析地推导 (\partial K/\partial t_i)（构件厚度），可以在 MATLAB 里用有限差分法逼近：

% 对厚度 t 做 0.1% 扰动，观察某个特征值的灵敏度 delta = 1e-3 * t; K_pert = update_stiffness(t + delta); omega_pert = sqrt(eig(K_pert, M)); sens = (omega_pert - omega) / delta;

这个做法的前提是你的刚度矩阵函数必须是由你自己组装的，比如前面 mx_form.m 那一套。如果你只是从 ANSYS 导出一次矩阵，那就只能做一次性分析，没法做参数化扰动。所以我一般会把 ANSYS 当成“高精度数值参考”，把可参数化的 MATLAB 组装程序当“设计迭代引擎”，两边对同一模型做交叉验证。

5.3 实际输出后的自检清单

分析做完之后，有一件事值得做：把你组装的全局矩阵的稀疏模式画出来，和 ANSYS 导出的矩阵对比。用spy(K)看非零元素分布，如果两个矩阵的自由度排序方式一致，图案应该几乎一样。由于 ANSYS 内部可能对自由度重新排序，图案会不同，但非零元总数应该接近。若非零数差异巨大，多半是单元连接表出错，某个单元的两个节点编号写反，或者遗漏了某个单元。资源里的EXTRACT.plg是经典 ANSYS 的宏文件，里面定义了提取矩阵的命令流，你可以直接把它拖进 ANSYS 里跑，然后对照生成的 .full 文件和你的 MATLAB 组装结果。只要这一步对上了，后续用减缩矩阵做子结构、固定界面模态综合，就都有了可信的底子。跑通一次之后，把整套流程封装成函数，以后换个模型只需改输入文件，不用再重写一遍提取逻辑。

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