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 - η) / 4
N₂ = (1 + ξ)(1 - η) / 4
N₃ = (1 + ξ)(1 + η) / 4
N₄ = (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 × 2N，N 为节点总数。组装时用稀疏矩阵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 型梁常用材料参数

网格尺寸对结果影响很大。我一般先用 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 kN：

P = -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 2，elements: 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 和单元位移 ue：

function 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); % 形心处 xi=0, eta=0 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 应力参考趋势

可以看到，拐角应力随网格加密持续上升，而位移基本收敛。此时不应把应力收敛作为指标，应改用拐角以外区域的应力分布或固定端反力做验证。如果要设计使用，建议在拐角处挖一个小圆角或使用子模型法，直接取角点应力会导致偏保守甚至错误的设计。

5.3 复现 topFig5.m 时的检查清单

我把这个项目复现时的排查点整理成如下清单：

1. 节点编号和单元编号是否从 1 开始；如果包含 0，Matlab 会自动当作逻辑索引，导致刚阵维度错乱。
2. 单元节点顺序是否逆时针；顺序错乱时det(J)为负，单刚矩阵奇异。
3. free_dofs是否包含了所有未约束自由度；若发现Kff条件数极大，优先检查是否有孤立节点。
4. 载荷单位是否统一；kN 和 N 混用是最常见的数值错误。
5. 绘制位移云图时，缩放系数不要直接取scale = 1；建议用max(u) / max(nodes)自动缩放，否则变形图可能小到看不见或大到完全覆盖网格。

最后还有一个实用技巧：在求解前先对K做一次condest检查。如果条件数超出 1e12，说明单位制有问题或边界约束不足，这时即使能求出位移，结果也不可信。把这一行检查写进脚本，能省下大量排查时间。

% 求解前检查整体刚阵病态程度 if condest(K(free_dofs, free_dofs)) > 1e12 warning('整体刚阵严重病态，请检查单位、材料和约束'); end

把这个condest检查放在Kff \ Ff之前，比任何调试都直接。无温度载荷 L 型梁的仿真本质上是一个经典线性静力问题，只要几何、材料、约束、载荷四项没有矛盾，求解器给出的结果就是稳定的。真正花时间的反而在网格收敛性判断和应力结果解释上。

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