读研那会儿第一次做有限元课设导师丢给我一本教材让我用MATLAB把平面四节点等参单元跑通。教材上的公式每一个字我都认识可一落到程序里就卡壳Jacobian矩阵到底按行组织还是按列组织、形函数对物理坐标求导时链式法则该正向用还是反向用、Gauss积分点循环应该先写ξ还是先写η。最气人的是程序好不容易跑通了拿结果跟教材例题对却对不上一度怀疑是自己数学没学好。后来把每个中间矩阵都打印出来逐项核对才彻底弄明白问题出在哪。这篇文章就按我当时踩坑的顺序把平面等参四边形单元刚度矩阵的计算原理和MATLAB完整代码讲清楚所有代码直接复制就能跑适合正在啃有限元教材、准备做课程设计或者单纯想把经典单元亲手实现一遍的朋友。1. 搞清楚等参四节点单元刚度矩阵到底在算什么先明确一个基本概念单元刚度矩阵是连接节点力与节点位移的线性关系表达式是 f Ke × d。对平面四节点四边形单元来说每个节点有 u、v 两个自由度四个节点加起来就是 8 个自由度所以 Ke 是一个 8×8 矩阵。矩阵第 i 行第 j 列元素的物理意义是当第 j 个自由度产生单位位移时在第 i 个自由度上需要施加的力。这个概念先摆在这里后面写代码时所有矩阵维度都是围着它转的。那等参两个字是什么意思简单说就是描述单元几何形状的插值函数和描述单元内位移场的插值函数用的是同一组形函数。单元的坐标由节点坐标插值得到位移也由节点位移插值得到两者共用形函数 Ni所以叫等参isoparametric。这个设计的直接好处是不管单元形状怎么变只要形函数能映射坐标就能同样映射位移程序不用为不同形状的单元单独写逻辑。计算一个四节点等参单元的刚度矩阵需要先准备以下输入数据输入参数含义矩阵/大小coord4 个节点的物理坐标 (x, y)4×2E弹性模量标量nu泊松比标量t单元厚度标量probType平面应力或平面应变字符串整个计算流程可以压缩成五步第一步在自然坐标 (ξ, η) 上构造形函数并求其对自然坐标的偏导数第二步通过坐标映射关系求出 Jacobian 矩阵第三步用链式法则把形函数对自然坐标的导数转换到物理坐标第四步组装应变矩阵 B 和弹性矩阵 D第五步在若干 Gauss 积分点上对 BDB×det(J)×t 做加权求和。下面每一章对应其中一步代码也是按这个结构写的。2. 从自然坐标到物理坐标形函数与Jacobian的映射关系2.1 四节点等参元的形函数表达式四节点四边形单元定义在自然坐标系 (ξ, η) 上取值范围是 [-1, 1]×[-1, 1]。四个节点按逆时针排列它们的自然坐标是节点ξη1-1-121-13114-11对应的四个形函数是N1 0.25×(1-ξ)×(1-η) N2 0.25×(1ξ)×(1-η) N3 0.25×(1ξ)×(1η) N4 0.25×(1-ξ)×(1η)这组形函数是二维拉格朗日插值的直接结果它有两个关键性质。第一形函数 Ni 在自身节点 i 处取值 1在其余三个节点处取值 0。第二四个形函数之和在单元内任意位置恒等于 1也就是 ΣNi ≡ 1。第二个性质非常重要它保证了单元的刚体位移模式可以被精确描述——如果所有节点都有相同的位移分量单元内任意一点的位移也会等于这个常数值不会出现不该有应变时算出力的怪事。2.2 Jacobian矩阵怎么算有了形函数就可以建立物理坐标和自然坐标的映射关系x Σ Ni(ξ,η) × xi y Σ Ni(ξ,η) × yi其中 xi、yi 是第 i 个节点的物理坐标。这个公式的意思是单元内任意一点 (x, y) 都可以由四个节点坐标加权得到权重就是形函数。接下来关键问题来了形函数 Ni 是自然坐标的函数但应变计算需要的是形函数对物理坐标 x、y 的偏导数怎么办答案是 Jacobian 矩阵。Jacobian 矩阵 J 定义为物理坐标对自然坐标的偏导矩阵J [ ∂x/∂ξ ∂y/∂ξ ] [ ∂x/∂η ∂y/∂η ]它的每个元素都可以由形函数求导再和节点坐标加权得到。比如 ∂x/∂ξ Σ (∂Ni/∂ξ)×xi∂y/∂ξ Σ (∂Ni/∂ξ)×yi其余类似。对矩形或平行四边形单元J 是常数矩阵对任意形状的四边形J 内部随 ξ、η 变化这就意味着必须在每个积分点位置重新计算 J。2.3 链式法则把自然坐标导数转换到物理坐标有了 J 之后用多元函数链式法则可以把两组导数联系起来∂Ni/∂ξ (∂Ni/∂x)(∂x/∂ξ) (∂Ni/∂y)(∂y/∂ξ) ∂Ni/∂η (∂Ni/∂x)(∂x/∂η) (∂Ni/∂y)(∂y/∂η)写成矩阵形式就是[∂Ni/∂ξ] [∂Ni/∂x] [∂Ni/∂η] J × [∂Ni/∂y]所以只要对 J 求逆就能得到[∂Ni/∂x] [∂Ni/∂ξ] [∂Ni/∂y] J⁻¹ × [∂Ni/∂η]这一段是初学者最容易搞反的地方。我当初写代码时惯性思维以为求导矩阵直接倒过来就行实际上必须通过 J 的逆来转换。MATLAB 里可以直接用 inv(J)也可以手写 2×2 逆矩阵公式对于 2×2 矩阵两者差别不大。这里还要强调节点编号顺序。为了保证 J 的行列式 det(J) 严格大于 0四个节点必须按逆时针顺序输入。反过来的话坐标映射会发生折叠det(J) 变负算出来的 Ke 会出现负能量模式结果完全错误。程序里加一句 det(J)0 的检查是成本最低的防错手段。3. 应变矩阵B和弹性矩阵D连接位移场与应力场3.1 B矩阵为什么是3×8平面问题上每个点的应变状态有三个分量εxx 方向正应变、εyy 方向正应变、γxy剪应变。这三个分量可以由位移场 u(x,y)、v(x,y) 求偏导得到εx ∂u/∂x εy ∂v/∂y γxy ∂u/∂y ∂v/∂x而单元内的位移场由节点位移插值得到u Σ Ni×uiv Σ Ni×vi。把插值关系代入应变表达式就能把应变写成节点位移的线性组合组合系数构成的矩阵就是 B 矩阵。每个节点贡献两个列分别对应 u 自由度和 v 自由度所以总列数是 8行数是 3B 就是 3×8 矩阵。自由度排列顺序是 [u1, v1, u2, v2, u3, v3, u4, v4]。在这个排列下B 矩阵组装逻辑是第 1 行εx 方向把 ∂Ni/∂x 放在第 2i-1 列第 2 行εy 方向把 ∂Ni/∂y 放在第 2i 列第 3 行剪应变 γxy∂Ni/∂y 放在第 2i-1 列∂Ni/∂x 放在第 2i 列这个排列规则是写代码时最容易出错的细节。剪应变那一行特别容易漏掉某一项或写成 ∂Ni/∂x 和 ∂Ni/∂y 交换位置一旦写错刚度矩阵行列元素全盘错位而且单独看矩阵很难发现。验证方法后面会单独讲。3.2 弹性矩阵D平面应力与平面应变的区别弹性矩阵 D 描述材料本构关系也就是应力和应变之间的线性关系。对平面问题D 是 3×3 矩阵但平面应力和平面应变两种情况表达式不一样。平面应力假设面外应力 σz 0适用于薄板类结构比如薄板受面内载荷D E/(1-ν²) × [1, ν, 0; ν, 1, 0; 0, 0, (1-ν)/2]平面应变假设面外应变 εz 0适用于长柱体或厚截面结构的横截面分析D E/((1ν)(1-2ν)) × [1-ν, ν, 0; ν, 1-ν, 0; 0, 0, (1-2ν)/2]这里有一个容易忽略的点ν 很接近 0.5 时平面应变公式里的分母 (1-2ν) 趋近于零D 矩阵数值会变得非常大这是材料近似不可压缩带来的数值特征不是程序 bug但实际使用时要留意。在代码里我建议用 switch 语句根据 probType 参数选择 D 矩阵而不是写两个独立函数。这样单元函数接口更干净调用方不容易出错。4. 2×2 Gauss积分把刚度积分变成离散求和单元刚度矩阵的完整表达式是Ke t × ∫∫ BDB × det(J) dξdη积分区域就是自然坐标系下的正方形 [-1, 1]×[-1, 1]。理论上可以用符号积分求出解析表达式但对任意四边形J 是变量B 又含 J 的逆解析展开非常冗长。工程上标准做法是数值积分其中用得最多的就是 Gauss 积分。两点 Gauss 积分在 [-1, 1] 区间上取积分点 ±1/√3 ≈ ±0.57735权重都是 1可以精确积分不超过三次的多项式。二维情况下在两个方向分别取点就得到 2×2 共四个积分点积分点ξη权重1-0.57735-0.57735120.57735-0.57735130.577350.5773514-0.577350.577351于是连续积分变成二重求和Ke t × Σᵢ Σⱼ wᵢwⱼ × [B(ξᵢ,ηⱼ)DB(ξᵢ,ηⱼ)] × det(J(ξᵢ,ηⱼ))讲一点上面的数学背景对矩形单元J 是常数矩阵B 的分量是 ξ、η 的一次式因此 BDB 是二次多项式2×2 Gauss 积分结果是精确的。对任意形状的四边形B 里含 1/det(J)严格说被积函数不是多项式2×2 积分只是工程近似。但对正常形状的单元误差很小这也是有限元教材里四节点单元默认用 2×2 积分的原因。为什么不考虑 1×1 积分因为 1 点积分会引入伪零能模式也就是常说的沙漏模式。8×8 的单元刚度矩阵理论上有 3 个零特征值对应 2 个平动和 1 个转动三个刚体自由度。如果改用 1×1 积分会出现第 4 个接近 0 的特征值单元虽然看起来是满了但实质上存在不需要消耗能量的变形路径网格一复杂就出大问题。这个经验在检查程序结果时非常有用后面验证部分还会再提。5. MATLAB完整代码与逐段验证5.1 单元刚度矩阵主程序现在把前面四章的内容全部浓缩成一个 MATLAB 函数。输入是材料参数、厚度、节点坐标和问题类型输出是 8×8 的单元刚度矩阵。为了照顾新手代码刻意写成最朴素的循环不搞矢量化炫技因为这里计算量很小清晰比高效重要。function Ke Quad4Stiffness(E, nu, t, coord, probType) % 计算平面4节点等参四边形单元的8×8刚度矩阵 % 输入 % E - 弹性模量 % nu - 泊松比 % t - 厚度 % coord - 4×2 节点坐标按逆时针排列 % probType - stress 平面应力 或 strain 平面应变 % 输出 % Ke - 8×8 单元刚度矩阵 if nargin 5 probType stress; end % 1. 弹性矩阵 D switch probType case stress D E/(1-nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; case strain D E/((1nu)*(1-2*nu)) * [1-nu, nu, 0; nu, 1-nu, 0; 0, 0, (1-2*nu)/2]; otherwise error(probType 只能是 stress 或 strain); end % 2. Gauss积分点与权重 gp [-1/sqrt(3), 1/sqrt(3)]; w [1, 1]; Ke zeros(8, 8); % 3. 双重循环遍历四个积分点 for i 1:2 for j 1:2 xi gp(i); eta gp(j); % 形函数对自然坐标的导数 dN_dxi 0.25 * [-(1-eta), (1-eta), (1eta), -(1eta)]; dN_deta 0.25 * [-(1-xi), -(1xi), (1xi), (1-xi)]; % Jacobian矩阵 J [dN_dxi; dN_deta] * coord; detJ det(J); if detJ 0 error(Jacobian行列式非正请检查节点编号是否按逆时针排列); end invJ inv(J); % 组装B矩阵 B zeros(3, 8); for node 1:4 dN_dx invJ(1,1) * dN_dxi(node) invJ(1,2) * dN_deta(node); dN_dy invJ(2,1) * dN_dxi(node) invJ(2,2) * dN_deta(node); B(1, 2*node-1) dN_dx; B(2, 2*node ) dN_dy; B(3, 2*node-1) dN_dy; B(3, 2*node ) dN_dx; end % 累加积分贡献 Ke Ke w(i) * w(j) * (B * D * B) * detJ * t; end end end代码里最关键的就是 B 矩阵组装那段内层循环。节点循环从 1 到 4每个节点分配两个自由度所以第 node 个节点对应的 u 自由度列号是 2×node-1v 自由度列号是 2×node。这个索引规律一旦写死就不会再出排列错位的问题。5.2 三个不变量验证程序是否正确程序写完不能直接信一定要验证。我常用的方法有三个每个都能在五分钟内定位大部分错误。第一个验证是检查刚体平移模式。把 8×8 刚度矩阵的所有行相加结果应该接近零向量。因为如果所有节点都有相同的平动位移单元不会产生应变也就不应该有节点力。MATLAB 里执行 max(abs(sum(Ke, 2)))理想结果是 1e-12 量级。如果不是说明 B 矩阵组装或积分有问题。第二个验证是特征值分析。对无约束单元8 个自由度对应 3 个刚体模式所以 Ke 应该有 3 个接近 0 的特征值剩下 5 个为正。如果出现 4 个接近 0 的特征值多半是积分阶次降到了 1×1出现了沙漏模式。用 eig(Ke) 排序后检查前 3 个数量级在 1e-13 左右属于正常。第三个验证是能量校验这个最直观也最能说服我。取一个边长为 2 的正方形单元坐标从 (-1,-1) 到 (1,1)令 E1、ν0、厚度 t1施加位移场 ux、v0。理论上这是单轴拉伸应变 εx1应变能密度是 0.5×E×ε²0.5单元体积是 4总应变能应该是 2。用程序算 Ke再算 0.5×d×Ke×d其中 d 是节点位移列向量结果应该非常接近 2。这个验证可以一次抓出 D 矩阵公式、B 矩阵组装、积分权重等几乎所有错误。% 验证脚本 E 200e9; nu 0.3; t 0.01; coord [0, 0; 1, 0; 1, 1; 0, 1]; % 逆时针编号 Ke Quad4Stiffness(E, nu, t, coord, stress); % 验证1刚体平移行和应为零向量 disp(行和最大值); disp(max(abs(sum(Ke, 2)))); % 验证2特征值前3个应接近0 e sort(eig(Ke)); disp(前4个特征值); disp(e(1:4)); % 验证3单轴拉伸应变能 E1 1; nu1 0; t1 1; coord1 [-1, -1; 1, -1; 1, 1; -1, 1]; K1 Quad4Stiffness(E1, nu1, t1, coord1, stress); d [-1, 0; 1, 0; 1, 0; -1, 0]; % ux, v0 dVec reshape(d, 8, 1); U 0.5 * dVec * K1 * dVec; disp(单轴拉伸应变能理论值2); disp(U);我在本地跑验证 3 的结果是 2.0000浮点误差在 1e-14 量级。每次改动代码之后我都重新跑一遍这三个验证一旦哪一步出现数量级异常马上能定位到是哪部分出了问题。5.3 从单元刚度矩阵到总刚拿到单元刚度矩阵后下一步通常是组装全局刚度矩阵。这个函数输出的自由度排列是 [u1, v1, u2, v2, u3, v3, u4, v4]总装时只需要构造一个单元自由度映射表把每个节点编号映射到全局自由度数然后把 Ke 的元素散到总刚对应位置。这一步在四节点单元里非常机械但对后面写求解器的人来说自由度顺序一致是减少 bug 的隐形保障。6. 实战中容易埋进去的坑和排查经验6.1 节点编号顺序逆时针还是顺时针这是新手最常见的报错来源。我见过很多次程序跑出来的 Ke 全是负的行列式、能量为负数的情况一查就是节点坐标按顺时针输了。解决方案其实很简单一是程序里写 detJ0 的检查二是输入坐标时养成画图确认的习惯。在 MATLAB 里用 patch 函数把四个节点按输入顺序连起来看一眼一秒就能发现顺序问题。6.2 单位制混用单位制这个坑不是程序问题但能让结果莫名其妙差 10 的若干次方。比如几何用毫米弹性模量却用 Pa算出来的刚度矩阵数值和位移结果完全对不上。我的习惯是全程采用一致单位制长度用 m弹性模量用 Pa厚度用 m力自动是 N或者长度用 mm弹性模量用 MPa厚度用 mm力还是 N。写计算书时把单位写清楚比事后检查省心得多。6.3 单元畸变对精度的影响四节点单元对形状畸变比较敏感。如果某个单元内角接近 180 度或者某两条边严重不等长Jacobian 行列式在积分点附近可能趋近于零B 矩阵数值暴涨刚度矩阵条件数变得极差。这种单元算出来的结果即使不报错精度也很差。网格划分时尽量让四边形接近矩形避免长条或尖角单元。如果问题几何实在复杂建议改用三角形单元过渡或者加密网格。6.4 平面应力还是平面应变选择要趁早什么时候用平面应力什么时候用平面应变很多人记不住。薄板受面内载荷面外可以自由变形用平面应力。长柱或厚板分析横截面面外被约束住用平面应变。如果实际结构介于两者之间两种结果可以看作上下界。程序实现上只有 D 矩阵不同但选错会让结果偏离真实情况很远这一点在选择 probType 参数时要想清楚。6.5 后续扩展八节点单元和三维单元这个程序结构稍微改改就能扩展。换成八节点四边形单元时形函数从 4 个变成 8 个每个节点仍然是 2 个自由度B 矩阵变成 3×16Gauss 积分阶次建议升到 3×3。三维六面体单元则是把自然坐标从 (ξ,η) 变成 (ξ,η,ζ)形函数变 8 个B 矩阵变 6×24积分阶次也是 3×3 起步。核心流程完全一致都是形函数导数→Jacobian→B矩阵→Gauss 加权求和这条主线所以把这套思路吃透之后迁移成本很低。最后再分享一个实际调试心得我后来每次写完类似的单元程序都会刻意构造一个简单载荷工况跟解析解或教材例题做对比。能量校验是性价比最高的手段如果应变能都对不上那基本不用看后面的应力结果。单元刚度矩阵是有限元程序的最小细胞这个细胞是对的往上一层总装、边界条件、求解器才有意义。
