简介面向ANSYS有限元初学者的刚度矩阵专题资料包围绕结构静力、动力分析中的核心概念覆盖了从刚阵定义、单元刚阵形成到全局组装、矩阵提取与后处理应用的完整链路。包内共39个文件以MATLAB的.m程序为主附带Fortran源程序、C文件、文本说明、动态库及可执行程序等合计仅380KB既可阅读理论文档也可直接运行脚本验证计算流程。已有358人学习下载适合需要结合代码理解刚度矩阵原理的工程师与研究生使用。内容包含ANSYS Workbench提取刚度矩阵的操作指引、框架结构刚阵组装脚本、坐标变换、静力凝聚与动力响应等模块能帮助读者加深对KuF方程、边界条件及网格质量影响的理解是一份轻量实用的入门参考资料。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*(1phi))那就考虑了剪切变形。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)*61: (5-1)*63这些自由度编号从集合里删掉。4. 从 ANSYS 里把刚度矩阵提取出来Static Condensation 与文件交换4.1 ANSYS 导出矩阵的两种路线很多人以为 ANSYS 里看不到刚度矩阵其实经典 ANSYSMechanical APDL可以直接用HBMAT命令导出刚度矩阵、质量矩阵和阻尼矩阵到外部文件。Workbench 则要借助 Solution 里的 Output Controls 或者用命令对象插入一段 APDL 代码。常见做法是在 Workbench 的 Model 下插入 CommandsAPDL写一段/SOLU ! 计算刚度矩阵并写入文件 WRFULL,1 HBMAT,K_matrix,txt, , , K, 1, YES, YESHBMAT后面第一个参数是文件名第二个是扩展名第三个是路径第四个留空第五个是矩阵类型标识K表示刚度矩阵第六个是格式1 表示 text2 表示 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 组装结果。只要这一步对上了后续用减缩矩阵做子结构、固定界面模态综合就都有了可信的底子。跑通一次之后把整套流程封装成函数以后换个模型只需改输入文件不用再重写一遍提取逻辑。本文还有配套的精品资源点击获取
ANSYS 与 MATLAB 中的刚度矩阵：组装、提取与静力凝聚
简介面向ANSYS有限元初学者的刚度矩阵专题资料包围绕结构静力、动力分析中的核心概念覆盖了从刚阵定义、单元刚阵形成到全局组装、矩阵提取与后处理应用的完整链路。包内共39个文件以MATLAB的.m程序为主附带Fortran源程序、C文件、文本说明、动态库及可执行程序等合计仅380KB既可阅读理论文档也可直接运行脚本验证计算流程。已有358人学习下载适合需要结合代码理解刚度矩阵原理的工程师与研究生使用。内容包含ANSYS Workbench提取刚度矩阵的操作指引、框架结构刚阵组装脚本、坐标变换、静力凝聚与动力响应等模块能帮助读者加深对KuF方程、边界条件及网格质量影响的理解是一份轻量实用的入门参考资料。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*(1phi))那就考虑了剪切变形。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)*61: (5-1)*63这些自由度编号从集合里删掉。4. 从 ANSYS 里把刚度矩阵提取出来Static Condensation 与文件交换4.1 ANSYS 导出矩阵的两种路线很多人以为 ANSYS 里看不到刚度矩阵其实经典 ANSYSMechanical APDL可以直接用HBMAT命令导出刚度矩阵、质量矩阵和阻尼矩阵到外部文件。Workbench 则要借助 Solution 里的 Output Controls 或者用命令对象插入一段 APDL 代码。常见做法是在 Workbench 的 Model 下插入 CommandsAPDL写一段/SOLU ! 计算刚度矩阵并写入文件 WRFULL,1 HBMAT,K_matrix,txt, , , K, 1, YES, YESHBMAT后面第一个参数是文件名第二个是扩展名第三个是路径第四个留空第五个是矩阵类型标识K表示刚度矩阵第六个是格式1 表示 text2 表示 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 组装结果。只要这一步对上了后续用减缩矩阵做子结构、固定界面模态综合就都有了可信的底子。跑通一次之后把整套流程封装成函数以后换个模型只需改输入文件不用再重写一遍提取逻辑。本文还有配套的精品资源点击获取
