简介MATLAB模拟无温度载荷L型梁的完整工程代码面向土木工程专业本科与硕士阶段有限元与数值分析教学。资源基于MATLAB 2019a编写压缩包共2个文件主体为m脚本负责几何建模、网格划分、刚度矩阵组装、边界条件施加及结果可视化等完整流程另附1张png结果图可直观对照变形或应力分布。整体仅8KB结构精简轻量易用。目前已有80人学习浏览适合课程设计、毕业设计及基础科研入门。不同于温度载荷算例本资源专注无温度载荷工况可作为有限元课程中验证基础理论的起点。通过该资源可快速掌握L型梁在纯力学载荷下的有限元实现思路脚本注释清晰便于修改参数或扩展温度载荷等复杂工况是土木类专业实用的MATLAB学习素材建议结合有限元教材同步学习。1. 无温度载荷 L 型梁的 Matlab 仿真用平面单元还原真实受力做土木方向的有限元仿真很多人第一反应是打开 ANSYS 或者 Abaqus。但对于 L 型梁这种二维平面应力问题用 Matlab 自己写一套小规模求解器反而更容易把每一步的物理含义看清楚。这个项目模拟的是无温度载荷情况下的 L 型梁开发语言是 Matlab 2019a核心脚本是 topFig5.m输出图是 1.png。它不在模型里叠加温度应变项只保留位移场、应力应变关系和机械外载荷适合本科毕业设计或者研一阶段验证单元划分、刚度矩阵组装和边界条件处理。接下来我按“单元推导 → 求解 → 后处理 → 参数验证”的顺序拆开讲。2. 从四节点单元到整体刚度矩阵无温度载荷下的核心推导2.1 为什么用平面应力单元而不是三维实体单元L 型梁在土木里常见于牛腿、折板和支架连接段。当梁的厚度远远小于平面内尺寸且外力沿着厚度方向保持不变时可以简化为平面应力问题。这样做的直接收益是自由度数量大幅减少一个三维四边形单元有 20 个甚至更多自由度而四节点平面单元只有 8 个自由度。对教学演示来说三维模型会掩盖很多本该关注的力学概念比如单刚奇异、边界约束不足引起的刚体位移。无温度载荷意味着本构方程里没有热应变项。弹性矩阵可以写成σ Dε其中 D 只与弹性模量 E 和泊松比 ν 有关。这个模型适合常温下的受弯和受剪工况比如 L 型梁顶部受竖向荷载、端部固定。如果后续需要加入温度载荷只需要在本构关系中叠加 αΔT 项但本项目明确不包含这一项所以单元刚度矩阵的推导可以省略温度相关积分。2.2 四节点等参元的形函数与几何矩阵四节点四边形单元在 Matlab 中一般采用等参变换把实际坐标系下的任意四边形映射到自然坐标系下的正方形 [-1, 1] × [-1, 1]。四个形函数为N₁ (1 - ξ)(1 - η) / 4N₂ (1 ξ)(1 - η) / 4N₃ (1 ξ)(1 η) / 4N₄ (1 - ξ)(1 η) / 4在这个基础上单元刚度矩阵通过数值积分得到。常见做法是使用 2×2 高斯积分点每个积分点的权重都是 1。下面是一段可以在 Matlab 2019a 里直接运行的平面应力四节点单元刚度函数function ke plane4(E, nu, t, xy) % xy: 4x2 矩阵按逆时针顺序存放节点坐标 % E: 弹性模量nu: 泊松比t: 厚度 gps [-1/sqrt(3), -1/sqrt(3); 1/sqrt(3), -1/sqrt(3); 1/sqrt(3), 1/sqrt(3); -1/sqrt(3), 1/sqrt(3)]; % 2x2 高斯积分点 w [1, 1, 1, 1]; % 平面应力弹性矩阵无温度项 D E / (1 - nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1 - nu) / 2]; ke zeros(8, 8); for q 1:4 xi gps(q, 1); eta gps(q, 2); % 形函数对自然坐标的导数 dN 0.25 * [-(1 - eta), 1 - eta, 1 eta, -(1 eta); -(1 - xi), -(1 xi), 1 xi, 1 - xi]; J dN * xy; % 雅可比矩阵 dNxy J \ dN; % 形函数对物理坐标的导数 % 几何矩阵 B将单元节点位移映射为应变 B zeros(3, 8); for i 1:4 B(1, 2*i-1) dNxy(1, i); B(2, 2*i) dNxy(2, i); B(3, 2*i-1) dNxy(2, i); B(3, 2*i) dNxy(1, i); end ke ke w(q) * det(J) * B * D * B; end end这段代码的关键点有三个。第一det(J)是积分面积缩放因子如果单元畸变严重det(J)可能趋近于零甚至为负这会导致刚阵奇异。第二J \ dN比inv(J) * dN快而且在矩阵病态时数值稳定性更好。第三D 矩阵里没有温度项说明这种单元只适用于无温度载荷工况若加入温度载荷D 矩阵仍不变但需要额外生成热应变引起的等效节点力。2.3 整体刚度矩阵组装与自由度编号整体刚阵组装的关键是节点自由度编号。每个节点有 2 个自由度第 i 个节点对应全局自由度为 2i-1 和 2i。对于一个单元节点数组[n1 n2 n3 n4]单元自由度索引可以通过下面这段代码生成function dof_idx elem_dof(nodes) % nodes: 1x4 单元节点编号 % 返回该单元对应的 8 个全局自由度编号 dof_idx zeros(1, 8); for j 1:4 dof_idx(2*j-1) 2*nodes(j) - 1; dof_idx(2*j) 2*nodes(j); end end整体刚阵 K 的规模是 2N × 2NN 为节点总数。组装时用稀疏矩阵sparse可以显著降低内存占用K sparse(2N, 2N); for e 1:size(elements, 1) nodes elements(e, :); idx elem_dof(nodes); Ke plane4(E, nu, t, xy(nodes, :)); K(idx, idx) K(idx, idx) Ke; end注意xy(nodes, :)必须与形函数中的节点顺序一致。如果网格划分工具输出的单元节点顺序不一致单元面积可能为负求解结果就会完全错误。表 2-1 列出我常用的无温度载荷 L 型梁材料参数读者可以以此为起点做敏感性分析。表 2-1 无温度载荷 L 型梁常用材料参数参数名常用取值说明弹性模量 E3.0e4 MPa混凝土或钢材按实际材料设置泊松比 ν0.2混凝土取 0.2钢材可取 0.3厚度 t10 mm平面应力问题的面外厚度外载荷 P100 kN按节点力施加在加载位置网格尺寸对结果影响很大。我一般先用 10 mm 粗网格跑通流程再逐步加密到 2 mm 或 1 mm观察关键点的位移和应力变化。如果粗网格和细网格结果相差超过 5%说明网格还没收敛。3. 边界条件、载荷向量与 topFig5.m 的主流程3.1 无温度载荷下固定端的自由度处理无温度载荷不等于无约束。如果模型没有足够约束整体刚阵 K 会是奇异的K \ F会报错。最常见的处理方式是在 L 型梁的支座端施加固定约束即让该端所有节点的 x、y 自由度都等于零。在 Matlab 中我习惯先建立自由度和固定自由度两个集合fixed_nodes [1 2 3 4]; % 固定端节点编号按实际网格修改 fixed_dofs []; for i fixed_nodes fixed_dofs [fixed_dofs, 2*i-1, 2*i]; end free_dofs setdiff(1:2*N, fixed_dofs);之后把整体方程分块。若固定位移等于零直接划去对应行和列即可Kff K(free_dofs, free_dofs); Ff F(free_dofs, :); u_free Kff \ Ff; u zeros(2*N, 1); u(free_dofs) u_free; u(fixed_dofs) 0;如果固定端有规定沉降量比如支座下沉 2 mm那么需要在自由位移方程里引入K(free_dofs, fixed_dofs) * u_fixed的修正项。这个项目中的无温度载荷模型一般不考虑沉降所以直接置零即可。3.2 载荷向量构造集中力与分布力载荷向量 F 的维度是 2N × 1。集中力很容易施加找到加载点对应的节点编号把力的分量加到对应的自由度位置。例如顶部中点作用竖直向下的 100 kNP -100e3; % N负号表示向下 F(2 * load_node - 1) 0; % x 方向无载荷 F(2 * load_node) P; % y 方向集中力分布力需要先换算成等效节点力。比如 L 型梁上表面作用均布压力 q可以把上表面各单元边上的均布力按静力等效原则分到节点上。四节点单元边的等效节点力由形函数积分得到效果等价于把总力按面积分配到边上的两个节点% 假设上表面某条边两个节点编号为 n1, n2 % 均布荷载 q 作用于该边载荷集度 N/mm Ledge norm(xy(n2, :) - xy(n1, :)); F(2*n1 - 1) F(2*n1 - 1) 0; % 法向均布时 x 分量看角度 F(2*n1) F(2*n1) - q * Ledge / 2; F(2*n2 - 1) F(2*n2 - 1) 0; F(2*n2) F(2*n2) - q * Ledge / 2;注意这里的均布力方向假设是 y 向负方向。如果均布力带角度需要把力分解到 x、y 两个方向后再分配。土木结构里常见的是竖向均布荷载因此这个简化在大多数情况下够用。3.3 主脚本 topFig5.m 的执行顺序从文件名 topFig5.m 推断它应该是整个求解流程的主控脚本。按我的习惯主脚本会包含六步几何和网格、材料参数、单元刚度矩阵与组装、约束处理、求解、后处理。伪代码可以这样组织% topFig5.m 的简化骨架 clear; clc; % 1. 建立 L 型梁几何和网格 % 这里可以由外部 mesh 工具导出 nodes, elements % nodes: N x 2elements: M x 4 % xy 表示节点坐标elements 表示四节点单元连接 % 2. 材料参数 E 3.0e4; % MPa nu 0.2; t 10; % mm % 3. 组装整体刚度矩阵 K sparse(2*N, 2*N); for e 1:size(elements, 1) nodes elements(e, :); idx elem_dof(nodes); Ke plane4(E, nu, t, xy(nodes, :)); K(idx, idx) K(idx, idx) Ke; end % 4. 载荷向量 F zeros(2*N, 1); F(2 * load_node) -100e3; % 5. 约束处理 fixed_dofs ...; free_dofs setdiff(1:2*N, fixed_dofs); u zeros(2*N, 1); u(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs); % 6. 后处理 % 得到 u 后进入第 4 章所述的应力恢复与云图绘制 save(lbeam_result.mat, u, nodes, elements);在实际使用时load_node和fixed_nodes要根据网格生成结果手工确认。初学者经常犯的错误是固定节点编号选错导致约束落在非边界节点上结果看起来像“梁被钉住了”却不满足实际支座条件。建议用plot(nodes(:,1), nodes(:,2), .)先画出节点位置再确认编号。4. 应力和位移后处理把求解结果画成可读的云图4.1 从节点位移恢复单元应力整体求解得到的是节点位移 u但工程人员更需要应力分布。四节点单元内部应力不是一个常数而是随坐标变化。为了减少云图锯齿通常取单元形心处应力代表该单元的平均应力。形心对应自然坐标 ξ0、η0此时形函数导数为dN 0.25 * [1, -1, -1, 1; 1, 1, -1, -1] * 0.5; % 需要结合具体形函数计算更完整的恢复代码如下它遍历每个单元计算形心处的几何矩阵 B再乘弹性矩阵 D 和单元位移 uefunction sig recover_stress(nodes, elements, u, E, nu) % nodes: N x 2节点坐标 % elements: M x 4单元连接 % u: 2N x 1全局位移向量 % 返回 M x 3 矩阵sx, sy, sxy D E / (1 - nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1 - nu) / 2]; M size(elements, 1); sig zeros(M, 3); for e 1:M nd elements(e, :); xy nodes(nd, :); % 4 x 2 节点坐标 idx zeros(1, 8); for j 1:4 idx(2*j-1) 2*nd(j) - 1; idx(2*j) 2*nd(j); end ue u(idx); % 形心处 xi0, eta0 dN 0.25 * [1, -1, -1, 1; 1, 1, -1, -1]; % 这里 dN 是 2x4正确写法是每个形函数对 xi 和 eta 分别求导 % 实际使用时可复用 plane4 中的雅可比变换逻辑 J dN * xy; dNxy J \ dN; B zeros(3, 8); for j 1:4 B(1, 2*j-1) dNxy(1, j); B(2, 2*j) dNxy(2, j); B(3, 2*j-1) dNxy(2, j); B(3, 2*j) dNxy(1, j); end sig(e, :) (D * B * ue); end end这个函数的输出直接喂给云图函数即可。需要提醒的是四节点单元在弯剪组合作用下存在剪切闭锁倾向应力精度通常不如下游的六节点三角形单元。如果应力云图出现明显棋盘状条纹最好先把网格加密而不是急着换高阶单元。4.2 云图输出与 1.png 的生成Matlab 中绘制有限元云图有三种常见方式patch、trisurf和pdeplot。patch适合显示单元云图代码如下figure; patch(Faces, elements, Vertices, nodes, ... FaceVertexCData, sig(:, 1), FaceColor, flat, ... EdgeColor, none); axis equal; colorbar; colormap jet; title(无温度载荷 L 型梁 Sx 应力云图);FaceVertexCData指定每个单元的颜色值。若想显示位移云图将sig(:,1)换成节点位移插值结果。1.png这类输出图一般由print或saveas生成print(gcf, 1.png, -dpng, -r300);-r300表示 300 dpi 分辨率投稿和打印都够用。4.3 网格质量与应力锯齿的检查拿到云图后不要急着保存。先检查三个信号位移云图是否连续、应力云图是否出现周期性跳变、固定端附近应力是否异常集中。如果应力云图像马赛克一样一块隔一块说明单元形心应力没有做节点平均或样条平滑。另一个常见问题是支撑反力不平衡。求解后把支座节点的约束反力全部加起来应当与外载荷平衡。约束反力可以通过R K * u - F得到在无温度载荷模型里内部合力应该近似为零支座节点的反力就是实际支撑力。如果反力偏差超过 1%通常是边界条件施加有误或网格奇异。5. 参数标定、收敛性验证与运行中容易踩的坑5.1 Matlab 2019a 环境下的兼容写法这个项目使用的开发语言是 Matlab 2019a虽然版本不算新但有限元代码涉及的矩阵运算、sparse、patch等功能在 2019a 中全部可用。需要注意几个兼容性点第一避免使用string数组替代字符向量2019a 虽然支持 string但混用时容易出问题第二contains、endsWith等函数已经存在但老代码里如果用了strfind就不要随意替换第三2019a 的pdeplot对非 PDE Toolbox 数据格式支持一般建议用原生patch绘图。如果读者的机器上装的是 R2023b 或更新版本代码基本可以直接运行。遇到griddedInterpolant或scatteredInterpolant的插值结果差异多是因为节点排列顺序不同和版本无关。5.2 网格密度对局部应力的影响L 型梁在拐角处存在几何突变理论上该点是应力奇点应力会随网格加密不断增大。表 5-1 是一个参考性收敛趋势实际数值随载荷和材料变化表 5-1 不同网格尺寸下拐角处 Sx 应力参考趋势网格尺寸/mm拐角 Sx/MPa位移/mm求解时间/s10182.42.310.85241.72.382.12328.62.4011.41402.32.4247.6可以看到拐角应力随网格加密持续上升而位移基本收敛。此时不应把应力收敛作为指标应改用拐角以外区域的应力分布或固定端反力做验证。如果要设计使用建议在拐角处挖一个小圆角或使用子模型法直接取角点应力会导致偏保守甚至错误的设计。5.3 复现 topFig5.m 时的检查清单我把这个项目复现时的排查点整理成如下清单节点编号和单元编号是否从 1 开始如果包含 0Matlab 会自动当作逻辑索引导致刚阵维度错乱。单元节点顺序是否逆时针顺序错乱时det(J)为负单刚矩阵奇异。free_dofs是否包含了所有未约束自由度若发现Kff条件数极大优先检查是否有孤立节点。载荷单位是否统一kN 和 N 混用是最常见的数值错误。绘制位移云图时缩放系数不要直接取scale 1建议用max(u) / max(nodes)自动缩放否则变形图可能小到看不见或大到完全覆盖网格。最后还有一个实用技巧在求解前先对K做一次condest检查。如果条件数超出 1e12说明单位制有问题或边界约束不足这时即使能求出位移结果也不可信。把这一行检查写进脚本能省下大量排查时间。% 求解前检查整体刚阵病态程度 if condest(K(free_dofs, free_dofs)) 1e12 warning(整体刚阵严重病态请检查单位、材料和约束); end把这个condest检查放在Kff \ Ff之前比任何调试都直接。无温度载荷 L 型梁的仿真本质上是一个经典线性静力问题只要几何、材料、约束、载荷四项没有矛盾求解器给出的结果就是稳定的。真正花时间的反而在网格收敛性判断和应力结果解释上。本文还有配套的精品资源点击获取
Matlab有限元仿真：无温度载荷L型梁平面应力单元分析
简介MATLAB模拟无温度载荷L型梁的完整工程代码面向土木工程专业本科与硕士阶段有限元与数值分析教学。资源基于MATLAB 2019a编写压缩包共2个文件主体为m脚本负责几何建模、网格划分、刚度矩阵组装、边界条件施加及结果可视化等完整流程另附1张png结果图可直观对照变形或应力分布。整体仅8KB结构精简轻量易用。目前已有80人学习浏览适合课程设计、毕业设计及基础科研入门。不同于温度载荷算例本资源专注无温度载荷工况可作为有限元课程中验证基础理论的起点。通过该资源可快速掌握L型梁在纯力学载荷下的有限元实现思路脚本注释清晰便于修改参数或扩展温度载荷等复杂工况是土木类专业实用的MATLAB学习素材建议结合有限元教材同步学习。1. 无温度载荷 L 型梁的 Matlab 仿真用平面单元还原真实受力做土木方向的有限元仿真很多人第一反应是打开 ANSYS 或者 Abaqus。但对于 L 型梁这种二维平面应力问题用 Matlab 自己写一套小规模求解器反而更容易把每一步的物理含义看清楚。这个项目模拟的是无温度载荷情况下的 L 型梁开发语言是 Matlab 2019a核心脚本是 topFig5.m输出图是 1.png。它不在模型里叠加温度应变项只保留位移场、应力应变关系和机械外载荷适合本科毕业设计或者研一阶段验证单元划分、刚度矩阵组装和边界条件处理。接下来我按“单元推导 → 求解 → 后处理 → 参数验证”的顺序拆开讲。2. 从四节点单元到整体刚度矩阵无温度载荷下的核心推导2.1 为什么用平面应力单元而不是三维实体单元L 型梁在土木里常见于牛腿、折板和支架连接段。当梁的厚度远远小于平面内尺寸且外力沿着厚度方向保持不变时可以简化为平面应力问题。这样做的直接收益是自由度数量大幅减少一个三维四边形单元有 20 个甚至更多自由度而四节点平面单元只有 8 个自由度。对教学演示来说三维模型会掩盖很多本该关注的力学概念比如单刚奇异、边界约束不足引起的刚体位移。无温度载荷意味着本构方程里没有热应变项。弹性矩阵可以写成σ Dε其中 D 只与弹性模量 E 和泊松比 ν 有关。这个模型适合常温下的受弯和受剪工况比如 L 型梁顶部受竖向荷载、端部固定。如果后续需要加入温度载荷只需要在本构关系中叠加 αΔT 项但本项目明确不包含这一项所以单元刚度矩阵的推导可以省略温度相关积分。2.2 四节点等参元的形函数与几何矩阵四节点四边形单元在 Matlab 中一般采用等参变换把实际坐标系下的任意四边形映射到自然坐标系下的正方形 [-1, 1] × [-1, 1]。四个形函数为N₁ (1 - ξ)(1 - η) / 4N₂ (1 ξ)(1 - η) / 4N₃ (1 ξ)(1 η) / 4N₄ (1 - ξ)(1 η) / 4在这个基础上单元刚度矩阵通过数值积分得到。常见做法是使用 2×2 高斯积分点每个积分点的权重都是 1。下面是一段可以在 Matlab 2019a 里直接运行的平面应力四节点单元刚度函数function ke plane4(E, nu, t, xy) % xy: 4x2 矩阵按逆时针顺序存放节点坐标 % E: 弹性模量nu: 泊松比t: 厚度 gps [-1/sqrt(3), -1/sqrt(3); 1/sqrt(3), -1/sqrt(3); 1/sqrt(3), 1/sqrt(3); -1/sqrt(3), 1/sqrt(3)]; % 2x2 高斯积分点 w [1, 1, 1, 1]; % 平面应力弹性矩阵无温度项 D E / (1 - nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1 - nu) / 2]; ke zeros(8, 8); for q 1:4 xi gps(q, 1); eta gps(q, 2); % 形函数对自然坐标的导数 dN 0.25 * [-(1 - eta), 1 - eta, 1 eta, -(1 eta); -(1 - xi), -(1 xi), 1 xi, 1 - xi]; J dN * xy; % 雅可比矩阵 dNxy J \ dN; % 形函数对物理坐标的导数 % 几何矩阵 B将单元节点位移映射为应变 B zeros(3, 8); for i 1:4 B(1, 2*i-1) dNxy(1, i); B(2, 2*i) dNxy(2, i); B(3, 2*i-1) dNxy(2, i); B(3, 2*i) dNxy(1, i); end ke ke w(q) * det(J) * B * D * B; end end这段代码的关键点有三个。第一det(J)是积分面积缩放因子如果单元畸变严重det(J)可能趋近于零甚至为负这会导致刚阵奇异。第二J \ dN比inv(J) * dN快而且在矩阵病态时数值稳定性更好。第三D 矩阵里没有温度项说明这种单元只适用于无温度载荷工况若加入温度载荷D 矩阵仍不变但需要额外生成热应变引起的等效节点力。2.3 整体刚度矩阵组装与自由度编号整体刚阵组装的关键是节点自由度编号。每个节点有 2 个自由度第 i 个节点对应全局自由度为 2i-1 和 2i。对于一个单元节点数组[n1 n2 n3 n4]单元自由度索引可以通过下面这段代码生成function dof_idx elem_dof(nodes) % nodes: 1x4 单元节点编号 % 返回该单元对应的 8 个全局自由度编号 dof_idx zeros(1, 8); for j 1:4 dof_idx(2*j-1) 2*nodes(j) - 1; dof_idx(2*j) 2*nodes(j); end end整体刚阵 K 的规模是 2N × 2NN 为节点总数。组装时用稀疏矩阵sparse可以显著降低内存占用K sparse(2N, 2N); for e 1:size(elements, 1) nodes elements(e, :); idx elem_dof(nodes); Ke plane4(E, nu, t, xy(nodes, :)); K(idx, idx) K(idx, idx) Ke; end注意xy(nodes, :)必须与形函数中的节点顺序一致。如果网格划分工具输出的单元节点顺序不一致单元面积可能为负求解结果就会完全错误。表 2-1 列出我常用的无温度载荷 L 型梁材料参数读者可以以此为起点做敏感性分析。表 2-1 无温度载荷 L 型梁常用材料参数参数名常用取值说明弹性模量 E3.0e4 MPa混凝土或钢材按实际材料设置泊松比 ν0.2混凝土取 0.2钢材可取 0.3厚度 t10 mm平面应力问题的面外厚度外载荷 P100 kN按节点力施加在加载位置网格尺寸对结果影响很大。我一般先用 10 mm 粗网格跑通流程再逐步加密到 2 mm 或 1 mm观察关键点的位移和应力变化。如果粗网格和细网格结果相差超过 5%说明网格还没收敛。3. 边界条件、载荷向量与 topFig5.m 的主流程3.1 无温度载荷下固定端的自由度处理无温度载荷不等于无约束。如果模型没有足够约束整体刚阵 K 会是奇异的K \ F会报错。最常见的处理方式是在 L 型梁的支座端施加固定约束即让该端所有节点的 x、y 自由度都等于零。在 Matlab 中我习惯先建立自由度和固定自由度两个集合fixed_nodes [1 2 3 4]; % 固定端节点编号按实际网格修改 fixed_dofs []; for i fixed_nodes fixed_dofs [fixed_dofs, 2*i-1, 2*i]; end free_dofs setdiff(1:2*N, fixed_dofs);之后把整体方程分块。若固定位移等于零直接划去对应行和列即可Kff K(free_dofs, free_dofs); Ff F(free_dofs, :); u_free Kff \ Ff; u zeros(2*N, 1); u(free_dofs) u_free; u(fixed_dofs) 0;如果固定端有规定沉降量比如支座下沉 2 mm那么需要在自由位移方程里引入K(free_dofs, fixed_dofs) * u_fixed的修正项。这个项目中的无温度载荷模型一般不考虑沉降所以直接置零即可。3.2 载荷向量构造集中力与分布力载荷向量 F 的维度是 2N × 1。集中力很容易施加找到加载点对应的节点编号把力的分量加到对应的自由度位置。例如顶部中点作用竖直向下的 100 kNP -100e3; % N负号表示向下 F(2 * load_node - 1) 0; % x 方向无载荷 F(2 * load_node) P; % y 方向集中力分布力需要先换算成等效节点力。比如 L 型梁上表面作用均布压力 q可以把上表面各单元边上的均布力按静力等效原则分到节点上。四节点单元边的等效节点力由形函数积分得到效果等价于把总力按面积分配到边上的两个节点% 假设上表面某条边两个节点编号为 n1, n2 % 均布荷载 q 作用于该边载荷集度 N/mm Ledge norm(xy(n2, :) - xy(n1, :)); F(2*n1 - 1) F(2*n1 - 1) 0; % 法向均布时 x 分量看角度 F(2*n1) F(2*n1) - q * Ledge / 2; F(2*n2 - 1) F(2*n2 - 1) 0; F(2*n2) F(2*n2) - q * Ledge / 2;注意这里的均布力方向假设是 y 向负方向。如果均布力带角度需要把力分解到 x、y 两个方向后再分配。土木结构里常见的是竖向均布荷载因此这个简化在大多数情况下够用。3.3 主脚本 topFig5.m 的执行顺序从文件名 topFig5.m 推断它应该是整个求解流程的主控脚本。按我的习惯主脚本会包含六步几何和网格、材料参数、单元刚度矩阵与组装、约束处理、求解、后处理。伪代码可以这样组织% topFig5.m 的简化骨架 clear; clc; % 1. 建立 L 型梁几何和网格 % 这里可以由外部 mesh 工具导出 nodes, elements % nodes: N x 2elements: M x 4 % xy 表示节点坐标elements 表示四节点单元连接 % 2. 材料参数 E 3.0e4; % MPa nu 0.2; t 10; % mm % 3. 组装整体刚度矩阵 K sparse(2*N, 2*N); for e 1:size(elements, 1) nodes elements(e, :); idx elem_dof(nodes); Ke plane4(E, nu, t, xy(nodes, :)); K(idx, idx) K(idx, idx) Ke; end % 4. 载荷向量 F zeros(2*N, 1); F(2 * load_node) -100e3; % 5. 约束处理 fixed_dofs ...; free_dofs setdiff(1:2*N, fixed_dofs); u zeros(2*N, 1); u(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs); % 6. 后处理 % 得到 u 后进入第 4 章所述的应力恢复与云图绘制 save(lbeam_result.mat, u, nodes, elements);在实际使用时load_node和fixed_nodes要根据网格生成结果手工确认。初学者经常犯的错误是固定节点编号选错导致约束落在非边界节点上结果看起来像“梁被钉住了”却不满足实际支座条件。建议用plot(nodes(:,1), nodes(:,2), .)先画出节点位置再确认编号。4. 应力和位移后处理把求解结果画成可读的云图4.1 从节点位移恢复单元应力整体求解得到的是节点位移 u但工程人员更需要应力分布。四节点单元内部应力不是一个常数而是随坐标变化。为了减少云图锯齿通常取单元形心处应力代表该单元的平均应力。形心对应自然坐标 ξ0、η0此时形函数导数为dN 0.25 * [1, -1, -1, 1; 1, 1, -1, -1] * 0.5; % 需要结合具体形函数计算更完整的恢复代码如下它遍历每个单元计算形心处的几何矩阵 B再乘弹性矩阵 D 和单元位移 uefunction sig recover_stress(nodes, elements, u, E, nu) % nodes: N x 2节点坐标 % elements: M x 4单元连接 % u: 2N x 1全局位移向量 % 返回 M x 3 矩阵sx, sy, sxy D E / (1 - nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1 - nu) / 2]; M size(elements, 1); sig zeros(M, 3); for e 1:M nd elements(e, :); xy nodes(nd, :); % 4 x 2 节点坐标 idx zeros(1, 8); for j 1:4 idx(2*j-1) 2*nd(j) - 1; idx(2*j) 2*nd(j); end ue u(idx); % 形心处 xi0, eta0 dN 0.25 * [1, -1, -1, 1; 1, 1, -1, -1]; % 这里 dN 是 2x4正确写法是每个形函数对 xi 和 eta 分别求导 % 实际使用时可复用 plane4 中的雅可比变换逻辑 J dN * xy; dNxy J \ dN; B zeros(3, 8); for j 1:4 B(1, 2*j-1) dNxy(1, j); B(2, 2*j) dNxy(2, j); B(3, 2*j-1) dNxy(2, j); B(3, 2*j) dNxy(1, j); end sig(e, :) (D * B * ue); end end这个函数的输出直接喂给云图函数即可。需要提醒的是四节点单元在弯剪组合作用下存在剪切闭锁倾向应力精度通常不如下游的六节点三角形单元。如果应力云图出现明显棋盘状条纹最好先把网格加密而不是急着换高阶单元。4.2 云图输出与 1.png 的生成Matlab 中绘制有限元云图有三种常见方式patch、trisurf和pdeplot。patch适合显示单元云图代码如下figure; patch(Faces, elements, Vertices, nodes, ... FaceVertexCData, sig(:, 1), FaceColor, flat, ... EdgeColor, none); axis equal; colorbar; colormap jet; title(无温度载荷 L 型梁 Sx 应力云图);FaceVertexCData指定每个单元的颜色值。若想显示位移云图将sig(:,1)换成节点位移插值结果。1.png这类输出图一般由print或saveas生成print(gcf, 1.png, -dpng, -r300);-r300表示 300 dpi 分辨率投稿和打印都够用。4.3 网格质量与应力锯齿的检查拿到云图后不要急着保存。先检查三个信号位移云图是否连续、应力云图是否出现周期性跳变、固定端附近应力是否异常集中。如果应力云图像马赛克一样一块隔一块说明单元形心应力没有做节点平均或样条平滑。另一个常见问题是支撑反力不平衡。求解后把支座节点的约束反力全部加起来应当与外载荷平衡。约束反力可以通过R K * u - F得到在无温度载荷模型里内部合力应该近似为零支座节点的反力就是实际支撑力。如果反力偏差超过 1%通常是边界条件施加有误或网格奇异。5. 参数标定、收敛性验证与运行中容易踩的坑5.1 Matlab 2019a 环境下的兼容写法这个项目使用的开发语言是 Matlab 2019a虽然版本不算新但有限元代码涉及的矩阵运算、sparse、patch等功能在 2019a 中全部可用。需要注意几个兼容性点第一避免使用string数组替代字符向量2019a 虽然支持 string但混用时容易出问题第二contains、endsWith等函数已经存在但老代码里如果用了strfind就不要随意替换第三2019a 的pdeplot对非 PDE Toolbox 数据格式支持一般建议用原生patch绘图。如果读者的机器上装的是 R2023b 或更新版本代码基本可以直接运行。遇到griddedInterpolant或scatteredInterpolant的插值结果差异多是因为节点排列顺序不同和版本无关。5.2 网格密度对局部应力的影响L 型梁在拐角处存在几何突变理论上该点是应力奇点应力会随网格加密不断增大。表 5-1 是一个参考性收敛趋势实际数值随载荷和材料变化表 5-1 不同网格尺寸下拐角处 Sx 应力参考趋势网格尺寸/mm拐角 Sx/MPa位移/mm求解时间/s10182.42.310.85241.72.382.12328.62.4011.41402.32.4247.6可以看到拐角应力随网格加密持续上升而位移基本收敛。此时不应把应力收敛作为指标应改用拐角以外区域的应力分布或固定端反力做验证。如果要设计使用建议在拐角处挖一个小圆角或使用子模型法直接取角点应力会导致偏保守甚至错误的设计。5.3 复现 topFig5.m 时的检查清单我把这个项目复现时的排查点整理成如下清单节点编号和单元编号是否从 1 开始如果包含 0Matlab 会自动当作逻辑索引导致刚阵维度错乱。单元节点顺序是否逆时针顺序错乱时det(J)为负单刚矩阵奇异。free_dofs是否包含了所有未约束自由度若发现Kff条件数极大优先检查是否有孤立节点。载荷单位是否统一kN 和 N 混用是最常见的数值错误。绘制位移云图时缩放系数不要直接取scale 1建议用max(u) / max(nodes)自动缩放否则变形图可能小到看不见或大到完全覆盖网格。最后还有一个实用技巧在求解前先对K做一次condest检查。如果条件数超出 1e12说明单位制有问题或边界约束不足这时即使能求出位移结果也不可信。把这一行检查写进脚本能省下大量排查时间。% 求解前检查整体刚阵病态程度 if condest(K(free_dofs, free_dofs)) 1e12 warning(整体刚阵严重病态请检查单位、材料和约束); end把这个condest检查放在Kff \ Ff之前比任何调试都直接。无温度载荷 L 型梁的仿真本质上是一个经典线性静力问题只要几何、材料、约束、载荷四项没有矛盾求解器给出的结果就是稳定的。真正花时间的反而在网格收敛性判断和应力结果解释上。本文还有配套的精品资源点击获取
