简介面向需要以 MATLAB 实现固体力学有限元分析的工程师与科研学习者这份资料包提供了一套以单元矩阵组装为核心的 FEM 基础脚本覆盖一维线性与二次单元、二维矩形域及圆柱坐标域等典型有限元示例串联几何建模、材料属性设置、边界条件施加、刚度矩阵组装与结果后处理等环节适合希望从理论过渡到编程实践的 MATLAB 用户可用于课程设计、科研入门和工程排错等场景。包体共 5 个 m 文件全部为可直接运行的 MATLAB 脚本整体仅 12KB便于快速阅读和修改目前已有 311 人学习下载。对照源码可以重点掌握一维线性/二次单元与二维平面/圆柱问题中局部刚度矩阵的构造逻辑以及全局矩阵索引组装方法对于需要扩展研究传热、弹塑性等固体力学问题的读者也提供了一套清晰可复用的代码模板。1. 为什么在 MATLAB 里从零手写 FEM 单元矩阵组装当 pde Toolbox 已经内置了比较完整的 FEA 工作流时重新用 M 脚本和循环组装一遍全局刚度矩阵看起来像是在重复造轮子。但把单元矩阵组装、边界条件索引和自由度编号一行行拆开看能解决很多黑箱报错。掌握刚度矩阵组装逻辑的工程师后续碰非线性、动态或自适应网格时定位问题会比调用包的用户快得多。压缩包里的五个文件——FEM1DL、FEM1DQ、FEM1DQnucd、FEM2DL_Box、FEM2DL_Cyl——覆盖一维线性、一维二次、二维平面四边形和轴对称圆柱的完整单元矩阵组装流程。这篇博文从局部刚度矩阵推导讲起落到自由度映射、数值积分和边界条件处理最后给出验证单元矩阵组装是否正确的方法。适合正在学有限元编程的 MATLAB 用户也是对固体力学计算流程感兴趣的工程人员的实战参考。2. 单元矩阵组装的核心局部自由度到全局编号2.1 连接矩阵决定组装位置在 MATLAB 的 FEM 代码里最容易出错的环节不是单元矩阵推导而是连接矩阵的使用。FEM1DL 中每个单元只有两个节点但全局网格往往有几十个节点局部刚度系数k(1,1)到底该加到全局矩阵的哪一行哪一列完全由连接矩阵决定。以一段长度为 0.5 m、离散成 5 个单元的杆为例节点编号按几何顺序排列连接矩阵写成% 一维杆单元5 个单元6 个节点 ne 5; nn ne 1; conn zeros(ne, 2); for e 1:ne conn(e, :) [e, e1]; % 第 e 个单元的全局节点编号 endconn(e,1)和conn(e,2)分别是第 e 个单元左、右节点在全局网格中的编号。一维问题每个节点只有一个位移自由度全局自由度编号等于全局节点编号组装时直接用节点号做矩阵索引。二维问题每个节点有 x、y 两个自由度全局自由度编号要变成2*(node-1)1和2*node索引关系完全不同。单元类型节点数每节点自由度全局自由度编号一维杆 FEM1DL21node一维杆 FEM1DQ / FEM1DQnucd31node二维平面 FEM2DL_Box422*(node-1)1、2*node轴对称 FEM2DL_Cyl422*(node-1)1、2*node2.2 局部刚度矩阵与全局组装循环杆单元的局部刚度矩阵由材料参数和几何参数唯一确定。弹性模量E200e9Pa、截面积A1e-4平方米、单元长度Le0.1m 时% 杆单元局部刚度矩阵 E 200e9; A 1e-4; Le 0.1; k_local E*A/Le * [1, -1; -1, 1];k_local是 2×2 对称矩阵物理含义是节点 1 发生单位位移而节点 2 不动时为维持单元平衡需要在两个节点上施加的力向量。对角线元素为正、副对角线为负表示刚度项抵抗相对位移的能力。这个矩阵推导自应变能积分E*A/Le是轴向刚度系数量纲为 N/m。全局组装循环把每个单元的局部矩阵累加到全局矩阵的对应位置K zeros(nn, nn); % 预分配全局刚度矩阵 F zeros(nn, 1); % 全局节点力向量 for e 1:ne n1 conn(e,1); n2 conn(e,2); K([n1 n2], [n1 n2]) K([n1 n2], [n1 n2]) k_local; endK([n1 n2],[n1 n2])是 MATLAB 的切片索引取出以n1、n2为行列的 2×2 子块。这里用累加赋值因为相邻单元在共享节点处比如节点 2 同时属于单元 1 和单元 2都会对同一个全局位置贡献刚度必须叠加否则就丢失了单元间传力路径。对单元数较大的网格建议预分配K后再组装避免循环内矩阵动态扩张。自由度规模到万级以上时把K声明为sparse(nn,nn)装配完成后用K\F求解稀疏直接法的内存占用和求解速度比稠密矩阵好一到两个数量级。2.3 组装正确性的快速自检组装写完先别急着求解做三项检查分别是对称性、连接矩阵顺序和约束后的条件数。检查项命令预期结果全局矩阵对称性norm(K - K)结果为 0连接矩阵逐行核对disp(conn(1:min(ne,10),:))每行节点号按几何相邻约束后非奇异condest(K(freeDofs,freeDofs))有限值不出现 Inf对称性检查最重要。有限元刚度矩阵理论上必须对称若norm(K-K)大于 1e-12说明组装逻辑里出现了不对称索引或重复累加。遇到这种情况我先查连接矩阵里是不是把某个单元的节点顺序写反了再查二维问题里自由度展开是否统一。3. 一维固体单元FEM1DL 与 FEM1DQ 的单元级实现3.1 FEM1DL 的形函数与 B 矩阵一维线性单元的位移场插值为 u(x) N1(x)·u1 N2(x)·u2形函数在局部坐标 ξ ∈ [-1,1] 上定义为N1(ξ) (1-ξ)/2N2(ξ) (1ξ)/2对应的几何矩阵应变-位移矩阵为B dN/dx (dN/dξ)·(dξ/dx) [-1/Le, 1/Le]其中 Le 是单元实际长度。应变 ε B·u应力 σ E·ε。单元刚度矩阵的连续形式是k_local ∫₀^Le Bᵀ·E·A·B·dx (EA/Le)·[[1, -1], [-1, 1]]这就是第 2 章里直接写出的那个 2×2 矩阵。线性单元的 B 矩阵是常量积分被解析求出不需要数值积分。这是 FEM1DL 和 FEM1DQ 在实现上最本质的区别前者一个单元就是一个常量应变场后者应变沿单元呈线性变化能捕捉更多变形细节。3.2 FEM1DQ 二次单元的形函数与数值积分FEM1DQ 用三节点二次单元局部坐标 ξ ∈ [-1,1]三个节点分别位于 ξ -1、0、1。二次形函数% 二次单元形函数xi 取 -1, 0, 1 时分别得到对应节点上的 1 N1 xi*(xi-1)/2; N2 1 - xi^2; N3 xi*(xi1)/2; % 对 xi 的导数 dN1dxi xi - 0.5; dN2dxi -2*xi; dN3dxi xi 0.5;N2在中节点处取 1在两端节点处取 0满足形函数的基本插值性质。与线性单元不同二次单元的 B 矩阵是坐标的函数刚度矩阵必须用数值积分。被积函数最高为二次多项式两点高斯积分恰好精确不需要更多积分点% 二点高斯积分计算二次单元 3x3 刚度矩阵 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; Le 0.2; E 200e9; A 1e-4; J Le/2; % 等参映射的雅可比 k zeros(3,3); for q 1:2 xi gp(q); dN [(xi-0.5), (-2*xi), (xi0.5)] / J; % 对物理坐标求导 k k gw(q) * (dN * E * A * dN) * J; enddN是 1×3 行向量dN * E * A * dN得到 3×3 局部刚度核乘上高斯权重和雅可比后累加。J Le/2是一维等参变换的映射率含义是局部坐标单位长度对应的物理长度。高斯点取在局部坐标上权重和雅可比缺一不可漏掉任一项结果都会整体偏小。3.3 非均匀网格版本的雅可比处理FEM1DQnucd 是一维二次单元的非均匀网格版本。均匀网格每个单元长度相同J可以一次算好复用非均匀网格各单元长度不同必须在循环内部对每个单元单独计算雅可比for e 1:ne Le xcoord(conn(e,2)) - xcoord(conn(e,1)); % 单元实际长度 J Le/2; % 在该单元上重新执行 3.2 的高斯积分循环 end这个文件的工程意义在于梁的固支端应力梯度大局部加密后用短单元捕捉应力变化远端用长单元控制总规模。同样的精度目标非均匀网格比均匀网格少用一半以上的单元。实际运行中注意xcoord的排序必须与连接矩阵一致否则算出的 Le 是负值刚度矩阵整体变号。3.4 一维边界条件处理与求解固定端直接划掉对应自由度即可fixedDofs 1; % 最左端固定 freeDofs setdiff(1:nn, fixedDofs); u zeros(nn, 1); u(freeDofs) K(freeDofs, freeDofs) \ F(freeDofs);setdiff返回补集把固定自由度和自由自由度分离。求解只作用于自由自由度最后把固定自由度对应的位移值赋回完整向量。多个固定节点时把fixedDofs换成行向量即可setdiff天然去重排序。用悬臂梁自由端受集中力的解析解 δ PL³/(3EI) 做验证FEM1DL 用 10 个单元能把位移误差控制在 1% 以内FEM1DQ 用 3 个单元就达到同等精度。二次单元在平滑解问题上的效率优势由此可见。4. 二维单元矩阵FEM2DL_Box 与 FEM2DL_Cyl 的组装差异4.1 双线性四边形单元的等参映射二维双线性四边形单元每个节点有 u_x、u_y 两个自由度局部坐标 (ξ,η) ∈ [-1,1]²。四个形函数% 双线性四边形单元形函数 N1 (1-xi)*(1-eta)/4; N2 (1xi)*(1-eta)/4; N3 (1xi)*(1eta)/4; N4 (1-xi)*(1eta)/4;等参变换把物理坐标表达成形函数与节点坐标的线性组合x ΣNi·xiy ΣNi·yi。雅可比矩阵 J 为J [∂x/∂ξ, ∂y/∂ξ; ∂x/∂η, ∂y/∂η]雅可比行列式 detJ 出现在面积变换关系 dx·dy detJ·dξ·dη 中。单元形状越接近正方形detJ 在单元内分布越均匀如果单元严重扭曲detJ 可能在某些高斯点变负组装出的刚度矩阵就不再正定求解直接发散。4.2 FEM2DL_Box 的数值积分与刚度核平面应力问题的应力-应变关系矩阵 D 为% 平面应力 D 矩阵 E 210e9; nu 0.3; factor E/(1-nu^2); D factor * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2];D 矩阵描述应力张量和应变张量的线性关系。平面应力和平面应变的区别只在 D 矩阵形式平面应变时把 E 替换为 E/(1-ν²)把 ν 替换为 ν/(1-ν)。用错 D 矩阵结果差一个量级表面上却看不出代码问题。在每个高斯点上先算形函数对局部坐标的导数通过雅可比逆变换得到对物理坐标的导数再组装 B 矩阵3×8并累加刚度% 2x2 高斯积分计算四边形单元局部刚度 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; k zeros(8,8); for q1 1:2 for q2 1:2 xi gp(q1); eta gp(q2); % 计算形函数导数与雅可比略 % dN_dx: 2x4 形函数对物理坐标导数 % B: 3x8 应变-位移矩阵 k k gw(q1)*gw(q2) * (B * D * B * detJ); end endB * D * B得到 8×8 局部刚度核乘上两个方向高斯权重和雅可比行列式累加。2×2 高斯积分对双线性单元已经完整捕捉 ξη 交叉项再加密积分点只会增加耗时不会提高精度这是高斯积分阶次选定的边界条件。4.3 二维组装索引与自由度展开二维组装和一维逻辑相同差别在自由度展开dofs [2*conn(e,1)-1, 2*conn(e,1), 2*conn(e,2)-1, 2*conn(e,2), ... 2*conn(e,3)-1, 2*conn(e,3), 2*conn(e,4)-1, 2*conn(e,4)]; K(dofs, dofs) K(dofs, dofs) k;2*(node-1)1是 x 方向自由度2*node是 y 方向自由度。这个编号约定必须和边界条件施加方式一致。初学时常犯的错误是用conn(e,:)直接做索引忘了每个节点有两个自由度导致全局矩阵维度不够MATLAB 直接报索引溢出。固定边界条件在二维下的施加方式fixedNodes find(xcoord 1e-12); % 找 x0 边界节点 fixedDofs []; for n fixedNodes fixedDofs [fixedDofs, 2*n-1, 2*n]; end freeDofs setdiff(1:2*nn, fixedDofs); u zeros(2*nn, 1); u(freeDofs) K(freeDofs, freeDofs) \ F(freeDofs);用坐标值判断边界节点有浮点比较隐患建议网格生成时用边界标记数组记录边界节点编号组装时直接读取避免1e-12这种魔法阈值。4.4 FEM2DL_Cyl 的轴对称假设与特殊处理轴对称问题把三维柱坐标 (r, θ, z) 缩减到 (r, z) 平面前提是几何、载荷和约束都关于 θ 对称。此时周向应变 ε_θθ u_r / r 不为零应变向量比平面问题多一项ε [ε_rr, ε_zz, ε_θθ, γ_rz]ᵀD 矩阵从 3×3 扩为 4×4B 矩阵从 3×8 扩为 4×8。刚度积分时被积函数多出径向坐标 r面积元为 2π·r·detJ·dξ·dη积分完成后乘 2π。轴对称单元的坑集中在对称轴 r0 附近。ε_θθ u_r/r 在 r 趋向 0 时数值上会爆掉需要特殊处理处理位置做法适用前提高斯点落在 r≤eps跳过该积分点并用相邻高斯点加权单元较小r0 处强制 ε_θθ 为 0直接置零轴对称问题的运动学约束对薄壁圆筒施加内压 p 时周向应力的解析解为 p·r/t。这个算例是验证 FEM2DL_Cyl 组装逻辑最快的方法计算得到的周向应力若与解析解偏差超过 2%优先检查 ε_θθ 项是否漏掉、r 因子是否乘入面积元。5. 验证单元矩阵组装从补片测试到收敛率5.1 补片测试验证组装逻辑最经典的手段是补片测试。取几个任意形状的单元组成补片施加能产生线性位移场的边界条件检查内部节点的解是否精确等于该线性场。若单元矩阵组装正确补片测试的误差收敛到机器精度。% 补片测试施加线性位移场 u_exact (x, y) 1e-3 * (2*x 3*y); % 对所有边界节点施加 u_exact % 内部节点求解后与 u_exact 比较误差应低于 1e-12 err max(abs(u_interior - u_exact_interior));补片测试对一维到二维、线性和二次单元都适用。误差在 1e-6 量级时多半是数值积分点数不够或 B 矩阵某一列写错。误差完全不收敛则说明连接矩阵或自由度编号存在系统性问题。5.2 网格加密与收敛率验证补片测试通过后用网格加密确认收敛阶。线性单元的位移 L2 范数误差随单元尺寸 h 满足二阶收敛loglog 斜率约 2二次单元是三阶收敛斜率约 3% 连续加密网格并记录 L2 误差 h_list [0.2, 0.1, 0.05, 0.025]; err_list zeros(size(h_list)); for i 1:numel(h_list) % 生成网格、组装、求解得到数值解 u_num % 误差sqrt(sum((u_num - u_exact).^2 .* element_area)) err_list(i) sqrt(sum((u_num - u_exact).^2 .* area)); end loglog(h_list, err_list, o-); hold on; loglog(h_list, h_list.^2, --); % 参考二阶斜率area是每个单元的面积乘进误差计算体现积分意义。实测斜率低于理论值优先怀疑边界条件施加不完整其次是刚度矩阵组装中漏掉某个单元。斜率对但绝对误差偏大则检查材料参数和单位制。5.3 常见报错速查表报错现象原因处理方式Matrix is singular to working precision约束缺失或自由度重复检查 fixedDofs 是否覆盖所有刚性位移结果对称但量级错误平面应力与平面应变混用核对 D 矩阵和材料参数二维结果总比解析解低每节点两个自由度被当成一个检查 dofs 索引展开轴对称 r0 附近发散ε_θθ 奇异用极限值替换或跳过轴线处积分点用厚壁圆筒的解析解对 FEM2DL_Cyl 做验证很有效薄壁圆筒在均匀内压下的周向应力解析值为 p·r/t若单元矩阵组装正确有限元结果应逼近该值。验证通过后顺手把K改成sparse存储再对比求解耗时能直观体会稀疏格式在大模型里省下的内存和时间。本文还有配套的精品资源点击获取 MATLAB有限元编程实战：从零手写单元矩阵组装与验证方法
简介面向需要以 MATLAB 实现固体力学有限元分析的工程师与科研学习者这份资料包提供了一套以单元矩阵组装为核心的 FEM 基础脚本覆盖一维线性与二次单元、二维矩形域及圆柱坐标域等典型有限元示例串联几何建模、材料属性设置、边界条件施加、刚度矩阵组装与结果后处理等环节适合希望从理论过渡到编程实践的 MATLAB 用户可用于课程设计、科研入门和工程排错等场景。包体共 5 个 m 文件全部为可直接运行的 MATLAB 脚本整体仅 12KB便于快速阅读和修改目前已有 311 人学习下载。对照源码可以重点掌握一维线性/二次单元与二维平面/圆柱问题中局部刚度矩阵的构造逻辑以及全局矩阵索引组装方法对于需要扩展研究传热、弹塑性等固体力学问题的读者也提供了一套清晰可复用的代码模板。1. 为什么在 MATLAB 里从零手写 FEM 单元矩阵组装当 pde Toolbox 已经内置了比较完整的 FEA 工作流时重新用 M 脚本和循环组装一遍全局刚度矩阵看起来像是在重复造轮子。但把单元矩阵组装、边界条件索引和自由度编号一行行拆开看能解决很多黑箱报错。掌握刚度矩阵组装逻辑的工程师后续碰非线性、动态或自适应网格时定位问题会比调用包的用户快得多。压缩包里的五个文件——FEM1DL、FEM1DQ、FEM1DQnucd、FEM2DL_Box、FEM2DL_Cyl——覆盖一维线性、一维二次、二维平面四边形和轴对称圆柱的完整单元矩阵组装流程。这篇博文从局部刚度矩阵推导讲起落到自由度映射、数值积分和边界条件处理最后给出验证单元矩阵组装是否正确的方法。适合正在学有限元编程的 MATLAB 用户也是对固体力学计算流程感兴趣的工程人员的实战参考。2. 单元矩阵组装的核心局部自由度到全局编号2.1 连接矩阵决定组装位置在 MATLAB 的 FEM 代码里最容易出错的环节不是单元矩阵推导而是连接矩阵的使用。FEM1DL 中每个单元只有两个节点但全局网格往往有几十个节点局部刚度系数k(1,1)到底该加到全局矩阵的哪一行哪一列完全由连接矩阵决定。以一段长度为 0.5 m、离散成 5 个单元的杆为例节点编号按几何顺序排列连接矩阵写成% 一维杆单元5 个单元6 个节点 ne 5; nn ne 1; conn zeros(ne, 2); for e 1:ne conn(e, :) [e, e1]; % 第 e 个单元的全局节点编号 endconn(e,1)和conn(e,2)分别是第 e 个单元左、右节点在全局网格中的编号。一维问题每个节点只有一个位移自由度全局自由度编号等于全局节点编号组装时直接用节点号做矩阵索引。二维问题每个节点有 x、y 两个自由度全局自由度编号要变成2*(node-1)1和2*node索引关系完全不同。单元类型节点数每节点自由度全局自由度编号一维杆 FEM1DL21node一维杆 FEM1DQ / FEM1DQnucd31node二维平面 FEM2DL_Box422*(node-1)1、2*node轴对称 FEM2DL_Cyl422*(node-1)1、2*node2.2 局部刚度矩阵与全局组装循环杆单元的局部刚度矩阵由材料参数和几何参数唯一确定。弹性模量E200e9Pa、截面积A1e-4平方米、单元长度Le0.1m 时% 杆单元局部刚度矩阵 E 200e9; A 1e-4; Le 0.1; k_local E*A/Le * [1, -1; -1, 1];k_local是 2×2 对称矩阵物理含义是节点 1 发生单位位移而节点 2 不动时为维持单元平衡需要在两个节点上施加的力向量。对角线元素为正、副对角线为负表示刚度项抵抗相对位移的能力。这个矩阵推导自应变能积分E*A/Le是轴向刚度系数量纲为 N/m。全局组装循环把每个单元的局部矩阵累加到全局矩阵的对应位置K zeros(nn, nn); % 预分配全局刚度矩阵 F zeros(nn, 1); % 全局节点力向量 for e 1:ne n1 conn(e,1); n2 conn(e,2); K([n1 n2], [n1 n2]) K([n1 n2], [n1 n2]) k_local; endK([n1 n2],[n1 n2])是 MATLAB 的切片索引取出以n1、n2为行列的 2×2 子块。这里用累加赋值因为相邻单元在共享节点处比如节点 2 同时属于单元 1 和单元 2都会对同一个全局位置贡献刚度必须叠加否则就丢失了单元间传力路径。对单元数较大的网格建议预分配K后再组装避免循环内矩阵动态扩张。自由度规模到万级以上时把K声明为sparse(nn,nn)装配完成后用K\F求解稀疏直接法的内存占用和求解速度比稠密矩阵好一到两个数量级。2.3 组装正确性的快速自检组装写完先别急着求解做三项检查分别是对称性、连接矩阵顺序和约束后的条件数。检查项命令预期结果全局矩阵对称性norm(K - K)结果为 0连接矩阵逐行核对disp(conn(1:min(ne,10),:))每行节点号按几何相邻约束后非奇异condest(K(freeDofs,freeDofs))有限值不出现 Inf对称性检查最重要。有限元刚度矩阵理论上必须对称若norm(K-K)大于 1e-12说明组装逻辑里出现了不对称索引或重复累加。遇到这种情况我先查连接矩阵里是不是把某个单元的节点顺序写反了再查二维问题里自由度展开是否统一。3. 一维固体单元FEM1DL 与 FEM1DQ 的单元级实现3.1 FEM1DL 的形函数与 B 矩阵一维线性单元的位移场插值为 u(x) N1(x)·u1 N2(x)·u2形函数在局部坐标 ξ ∈ [-1,1] 上定义为N1(ξ) (1-ξ)/2N2(ξ) (1ξ)/2对应的几何矩阵应变-位移矩阵为B dN/dx (dN/dξ)·(dξ/dx) [-1/Le, 1/Le]其中 Le 是单元实际长度。应变 ε B·u应力 σ E·ε。单元刚度矩阵的连续形式是k_local ∫₀^Le Bᵀ·E·A·B·dx (EA/Le)·[[1, -1], [-1, 1]]这就是第 2 章里直接写出的那个 2×2 矩阵。线性单元的 B 矩阵是常量积分被解析求出不需要数值积分。这是 FEM1DL 和 FEM1DQ 在实现上最本质的区别前者一个单元就是一个常量应变场后者应变沿单元呈线性变化能捕捉更多变形细节。3.2 FEM1DQ 二次单元的形函数与数值积分FEM1DQ 用三节点二次单元局部坐标 ξ ∈ [-1,1]三个节点分别位于 ξ -1、0、1。二次形函数% 二次单元形函数xi 取 -1, 0, 1 时分别得到对应节点上的 1 N1 xi*(xi-1)/2; N2 1 - xi^2; N3 xi*(xi1)/2; % 对 xi 的导数 dN1dxi xi - 0.5; dN2dxi -2*xi; dN3dxi xi 0.5;N2在中节点处取 1在两端节点处取 0满足形函数的基本插值性质。与线性单元不同二次单元的 B 矩阵是坐标的函数刚度矩阵必须用数值积分。被积函数最高为二次多项式两点高斯积分恰好精确不需要更多积分点% 二点高斯积分计算二次单元 3x3 刚度矩阵 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; Le 0.2; E 200e9; A 1e-4; J Le/2; % 等参映射的雅可比 k zeros(3,3); for q 1:2 xi gp(q); dN [(xi-0.5), (-2*xi), (xi0.5)] / J; % 对物理坐标求导 k k gw(q) * (dN * E * A * dN) * J; enddN是 1×3 行向量dN * E * A * dN得到 3×3 局部刚度核乘上高斯权重和雅可比后累加。J Le/2是一维等参变换的映射率含义是局部坐标单位长度对应的物理长度。高斯点取在局部坐标上权重和雅可比缺一不可漏掉任一项结果都会整体偏小。3.3 非均匀网格版本的雅可比处理FEM1DQnucd 是一维二次单元的非均匀网格版本。均匀网格每个单元长度相同J可以一次算好复用非均匀网格各单元长度不同必须在循环内部对每个单元单独计算雅可比for e 1:ne Le xcoord(conn(e,2)) - xcoord(conn(e,1)); % 单元实际长度 J Le/2; % 在该单元上重新执行 3.2 的高斯积分循环 end这个文件的工程意义在于梁的固支端应力梯度大局部加密后用短单元捕捉应力变化远端用长单元控制总规模。同样的精度目标非均匀网格比均匀网格少用一半以上的单元。实际运行中注意xcoord的排序必须与连接矩阵一致否则算出的 Le 是负值刚度矩阵整体变号。3.4 一维边界条件处理与求解固定端直接划掉对应自由度即可fixedDofs 1; % 最左端固定 freeDofs setdiff(1:nn, fixedDofs); u zeros(nn, 1); u(freeDofs) K(freeDofs, freeDofs) \ F(freeDofs);setdiff返回补集把固定自由度和自由自由度分离。求解只作用于自由自由度最后把固定自由度对应的位移值赋回完整向量。多个固定节点时把fixedDofs换成行向量即可setdiff天然去重排序。用悬臂梁自由端受集中力的解析解 δ PL³/(3EI) 做验证FEM1DL 用 10 个单元能把位移误差控制在 1% 以内FEM1DQ 用 3 个单元就达到同等精度。二次单元在平滑解问题上的效率优势由此可见。4. 二维单元矩阵FEM2DL_Box 与 FEM2DL_Cyl 的组装差异4.1 双线性四边形单元的等参映射二维双线性四边形单元每个节点有 u_x、u_y 两个自由度局部坐标 (ξ,η) ∈ [-1,1]²。四个形函数% 双线性四边形单元形函数 N1 (1-xi)*(1-eta)/4; N2 (1xi)*(1-eta)/4; N3 (1xi)*(1eta)/4; N4 (1-xi)*(1eta)/4;等参变换把物理坐标表达成形函数与节点坐标的线性组合x ΣNi·xiy ΣNi·yi。雅可比矩阵 J 为J [∂x/∂ξ, ∂y/∂ξ; ∂x/∂η, ∂y/∂η]雅可比行列式 detJ 出现在面积变换关系 dx·dy detJ·dξ·dη 中。单元形状越接近正方形detJ 在单元内分布越均匀如果单元严重扭曲detJ 可能在某些高斯点变负组装出的刚度矩阵就不再正定求解直接发散。4.2 FEM2DL_Box 的数值积分与刚度核平面应力问题的应力-应变关系矩阵 D 为% 平面应力 D 矩阵 E 210e9; nu 0.3; factor E/(1-nu^2); D factor * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2];D 矩阵描述应力张量和应变张量的线性关系。平面应力和平面应变的区别只在 D 矩阵形式平面应变时把 E 替换为 E/(1-ν²)把 ν 替换为 ν/(1-ν)。用错 D 矩阵结果差一个量级表面上却看不出代码问题。在每个高斯点上先算形函数对局部坐标的导数通过雅可比逆变换得到对物理坐标的导数再组装 B 矩阵3×8并累加刚度% 2x2 高斯积分计算四边形单元局部刚度 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; k zeros(8,8); for q1 1:2 for q2 1:2 xi gp(q1); eta gp(q2); % 计算形函数导数与雅可比略 % dN_dx: 2x4 形函数对物理坐标导数 % B: 3x8 应变-位移矩阵 k k gw(q1)*gw(q2) * (B * D * B * detJ); end endB * D * B得到 8×8 局部刚度核乘上两个方向高斯权重和雅可比行列式累加。2×2 高斯积分对双线性单元已经完整捕捉 ξη 交叉项再加密积分点只会增加耗时不会提高精度这是高斯积分阶次选定的边界条件。4.3 二维组装索引与自由度展开二维组装和一维逻辑相同差别在自由度展开dofs [2*conn(e,1)-1, 2*conn(e,1), 2*conn(e,2)-1, 2*conn(e,2), ... 2*conn(e,3)-1, 2*conn(e,3), 2*conn(e,4)-1, 2*conn(e,4)]; K(dofs, dofs) K(dofs, dofs) k;2*(node-1)1是 x 方向自由度2*node是 y 方向自由度。这个编号约定必须和边界条件施加方式一致。初学时常犯的错误是用conn(e,:)直接做索引忘了每个节点有两个自由度导致全局矩阵维度不够MATLAB 直接报索引溢出。固定边界条件在二维下的施加方式fixedNodes find(xcoord 1e-12); % 找 x0 边界节点 fixedDofs []; for n fixedNodes fixedDofs [fixedDofs, 2*n-1, 2*n]; end freeDofs setdiff(1:2*nn, fixedDofs); u zeros(2*nn, 1); u(freeDofs) K(freeDofs, freeDofs) \ F(freeDofs);用坐标值判断边界节点有浮点比较隐患建议网格生成时用边界标记数组记录边界节点编号组装时直接读取避免1e-12这种魔法阈值。4.4 FEM2DL_Cyl 的轴对称假设与特殊处理轴对称问题把三维柱坐标 (r, θ, z) 缩减到 (r, z) 平面前提是几何、载荷和约束都关于 θ 对称。此时周向应变 ε_θθ u_r / r 不为零应变向量比平面问题多一项ε [ε_rr, ε_zz, ε_θθ, γ_rz]ᵀD 矩阵从 3×3 扩为 4×4B 矩阵从 3×8 扩为 4×8。刚度积分时被积函数多出径向坐标 r面积元为 2π·r·detJ·dξ·dη积分完成后乘 2π。轴对称单元的坑集中在对称轴 r0 附近。ε_θθ u_r/r 在 r 趋向 0 时数值上会爆掉需要特殊处理处理位置做法适用前提高斯点落在 r≤eps跳过该积分点并用相邻高斯点加权单元较小r0 处强制 ε_θθ 为 0直接置零轴对称问题的运动学约束对薄壁圆筒施加内压 p 时周向应力的解析解为 p·r/t。这个算例是验证 FEM2DL_Cyl 组装逻辑最快的方法计算得到的周向应力若与解析解偏差超过 2%优先检查 ε_θθ 项是否漏掉、r 因子是否乘入面积元。5. 验证单元矩阵组装从补片测试到收敛率5.1 补片测试验证组装逻辑最经典的手段是补片测试。取几个任意形状的单元组成补片施加能产生线性位移场的边界条件检查内部节点的解是否精确等于该线性场。若单元矩阵组装正确补片测试的误差收敛到机器精度。% 补片测试施加线性位移场 u_exact (x, y) 1e-3 * (2*x 3*y); % 对所有边界节点施加 u_exact % 内部节点求解后与 u_exact 比较误差应低于 1e-12 err max(abs(u_interior - u_exact_interior));补片测试对一维到二维、线性和二次单元都适用。误差在 1e-6 量级时多半是数值积分点数不够或 B 矩阵某一列写错。误差完全不收敛则说明连接矩阵或自由度编号存在系统性问题。5.2 网格加密与收敛率验证补片测试通过后用网格加密确认收敛阶。线性单元的位移 L2 范数误差随单元尺寸 h 满足二阶收敛loglog 斜率约 2二次单元是三阶收敛斜率约 3% 连续加密网格并记录 L2 误差 h_list [0.2, 0.1, 0.05, 0.025]; err_list zeros(size(h_list)); for i 1:numel(h_list) % 生成网格、组装、求解得到数值解 u_num % 误差sqrt(sum((u_num - u_exact).^2 .* element_area)) err_list(i) sqrt(sum((u_num - u_exact).^2 .* area)); end loglog(h_list, err_list, o-); hold on; loglog(h_list, h_list.^2, --); % 参考二阶斜率area是每个单元的面积乘进误差计算体现积分意义。实测斜率低于理论值优先怀疑边界条件施加不完整其次是刚度矩阵组装中漏掉某个单元。斜率对但绝对误差偏大则检查材料参数和单位制。5.3 常见报错速查表报错现象原因处理方式Matrix is singular to working precision约束缺失或自由度重复检查 fixedDofs 是否覆盖所有刚性位移结果对称但量级错误平面应力与平面应变混用核对 D 矩阵和材料参数二维结果总比解析解低每节点两个自由度被当成一个检查 dofs 索引展开轴对称 r0 附近发散ε_θθ 奇异用极限值替换或跳过轴线处积分点用厚壁圆筒的解析解对 FEM2DL_Cyl 做验证很有效薄壁圆筒在均匀内压下的周向应力解析值为 p·r/t若单元矩阵组装正确有限元结果应逼近该值。验证通过后顺手把K改成sparse存储再对比求解耗时能直观体会稀疏格式在大模型里省下的内存和时间。本文还有配套的精品资源点击获取