很多刚接触有限元的同学都有同一种感觉看一维杆单元的刚度矩阵推导还能跟上一到二维平面问题就懵了。尤其是遇上四边形单元——等参变换、Jacobian矩阵、高斯积分这些名词一下子砸过来教材上公式又多又长对着书看半天代码里还是不知道从哪一行开始写。我当年跟着教材手推第一个四节点等参单元的单刚矩阵时也花了整整两天踩了不少坑。这篇内容要做的就是把手把手带你把平面等参四边形单元的刚度矩阵用 MATLAB 完整实现出来。我会先讲清楚等参变换到底在做什么再给出可以直接运行的单刚计算函数附上完整的组装和求解脚本最后用简单算例验证理论解。适合两类人正在为有限元课程作业发愁的学生以及刚接手结构分析、想搞懂商业软件背后到底在算什么的新工程师。跟着走一遍你会发现“单元刚度矩阵”这件事并没有想象中那么神秘。1. 等参四边形的底层逻辑从“规则”到“任意”1.1 为什么四边形单元要做等参变换有限元里最早接触的三角形单元形函数可以直接在物理坐标里显式写出来比如面积坐标那一套推导还算直观。可一旦换成四边形麻烦就来了任意四边形在物理坐标里并没有一个简单统一的“双线性多项式”可以当作形函数而且不同单元的倾斜角度、边长比例都不一样每换一个单元都要重写一套函数这显然不现实。解决办法是曲线救国先把真实四边形“映射”到一个标准的正方形上去。这个正方形放在自然坐标里横坐标叫 ξ纵坐标叫 η范围都是 [-1, 1]。在这个标准正方形上四个形函数是固定的、特别简洁的双线性函数。然后利用坐标映射把标准正方形压扁、拉伸、旋转变成真实网格里的那个四边形。关键在于真实坐标和自然坐标之间的映射关系也采用同一套形函数来插值比如x N₁x₁ N₂x₂ N₃x₃ N₄x₄y N₁y₁ N₂y₂ N₃y₃ N₄y₄这里的 N₁ 到 N₄ 就是标准正方形上的形函数。位移场同样用这四个形函数插值几何插值和位移插值用的是同一套函数所以叫“等参”。这个词在英文里是 isoparametric“iso”就是“相同”的意思。这个思路的价值在于只要写好一次形函数和高斯积分就能处理任意形状的四边形单元。无论单元是被压扁的梯形、倾斜的平行四边形还是四条边都不规则的凸四边形核心计算流程完全不用改。商业有限元软件里 CPS4、Q4 这类单元能处理复杂网格靠的就是这套等参逻辑。1.2 四个双线性形函数到底有什么性质四节点等参单元的形函数写出来是四个双线性多项式N₁ (1-ξ)(1-η)/4N₂ (1ξ)(1-η)/4N₃ (1ξ)(1η)/4N₄ (1-ξ)(1η)/4注意它们有几个重要性质。第一每个形函数在自己的节点上等于 1在其他三个节点上等于 0。这一点保证了“插值性”也就是说单元节点的位移值就是真实结点位移不会出现奇奇怪怪的偏移。第二四个形函数在任何位置加起来恒等于 1这保证刚体位移模式下单元不会产生虚假应变——也就是说如果四个节点位移全是 1单元内任意点的插值位移也是 1应变自然为零。这个性质很多人会忽略但当你在调试代码时发现单元刚度矩阵有奇异模式多半就是形函数不满足常数完备性造成的。这里是“双线性”而不是“完全二次”每一张形函数包含常数项、一次项和 ξη 乘积项。这是四边形单元区别于三角形单元的另一个特色。ση 乘积项使得单元内部位移场可以出现轻微“扭曲”但沿单元边界的位移分布是一次的因此相邻等参四边形单元共享边上的位移是协调的不会开裂、重叠。2. 刚度矩阵的数学推导三步走2.1 从位移插值到应变矩阵每一步都很机械先写位移场再求应变最后组装成刚度矩阵。单元内任一点的位移被四个节点的位移插值出来u(ξ,η) Σ Nᵢ(ξ,η)·uᵢv(ξ,η) Σ Nᵢ(ξ,η)·vᵢ这里的 uᵢ、vᵢ 是节点 i 的两个位移分量。平面问题的应变有三个分量εx、εy、γxy。把位移代进去ε [∂u/∂x ; ∂v/∂y ; ∂u/∂y ∂v/∂x] B·d其中 d 是单元全部节点位移组成的 8×1 向量d [u₁, v₁, u₂, v₂, u₃, v₃, u₄, v₄]ᵀ。B 矩阵是 3×8 的应变-位移矩阵它的每一列对应一个节点自由度B [∂N₁/∂x 0 ∂N₂/∂x 0 ∂N₃/∂x 0 ∂N₄/∂x 00 ∂N₁/∂y 0 ∂N₂/∂y 0 ∂N₃/∂y 0 ∂N₄/∂y∂N₁/∂y ∂N₁/∂x ∂N₂/∂y ∂N₂/∂x ∂N₃/∂y ∂N₃/∂x ∂N₄/∂y ∂N₄/∂x]看到这里你可能会问形函数明明写在自然坐标里怎么对 x、y 求导这就轮到 Jacobian 矩阵出场了。2.2 Jacobian 矩阵连接两套坐标的桥梁自然坐标 (ξ,η) 和物理坐标 (x,y) 之间靠链式法则沟通。对任意一个形函数 Nᵢ有∂Nᵢ/∂ξ ∂Nᵢ/∂x · ∂x/∂ξ ∂Nᵢ/∂y · ∂y/∂ξ∂Nᵢ/∂η ∂Nᵢ/∂x · ∂x/∂η ∂Nᵢ/∂y · ∂y/∂η写成一个矩阵关系[dNᵢ/dξ] J · [dNᵢ/dx]转置形式具体看排列这里的 J 就是 2×2 的 Jacobian 矩阵。它的元素可以这样算既然 x Σ Nᵢ(ξ,η)·xᵢ那么 ∂x/∂ξ Σ (∂Nᵢ/∂ξ)·xᵢ其它元素同理。写成矩阵形式就是J [∂x/∂ξ ∂y/∂ξ∂x/∂η ∂y/∂η] [dN_dξ ; dN_dη] · nodes其中 nodes 是 4×2 的节点坐标矩阵。这个式子看起来抽象但在代码里异常直观就是两行导数向量乘坐标矩阵而已。J 有两个作用。第一个作用是把自然坐标下的导数换算成物理坐标下的导数方法是求逆dN/dx 和 dN/dy 由 J 的逆矩阵乘出来。第二个作用是提供面积映射因子。自然坐标系里一个微元面积 dξ·dη映射到物理坐标后面积变成 |J|·dξ·dη。所以说白了|J| 就是那个拉伸系数积分的时候自然而然地要乘上去。2.3 单刚积分与高斯数值积分平面问题单元刚度矩阵的标准形式是Kₑ ∫∫ Bᵀ·D·B·t·|J| dξ dη积分区域是自然坐标下 [-1,1]×[-1,1]。D 是材料弹性矩阵t 是厚度。这个积分的被积函数通常不是简单多项式尤其单元形状不规则时更复杂所以实际代码里几乎都用高斯数值积分代替解析积分。高斯积分的思路是把积分近似成被积函数在若干特定点上的加权求和。对二维问题就是对两个方向分别做高斯积分Kₑ ≈ Σᵢ Σⱼ wᵢ·wⱼ·B(ξᵢ,ηⱼ)ᵀ·D·B(ξᵢ,ηⱼ)·t·|J(ξᵢ,ηⱼ)|对四节点双线性单元最常用的是 2×2 高斯积分。四个积分点在 ξ±1/√3、η±1/√3权重全部是 1。为什么两个点就够了因为一维 n 点高斯积分对 2n-1 次多项式精确而双线性单元的被积函数最高次数不超过二次2 点高斯刚好能精确积分。对于矩形、平行四边形这类形状规则的单元2×2 就是精确解不是近似解。单元一旦严重畸变被积函数变复杂可以考虑升到 3×3但结果差异通常只在网格质量差的时候明显。3. MATLAB代码实现一个函数搞定单元刚度3.1 函数接口设计写代码不是把公式翻译成语法那么简单接口设计能省掉你后面一大堆调试时间。我建议把单刚计算独立成一个函数输入输出尽量少调用逻辑清晰。函数输入包括四部分节点坐标、材料参数、问题类型、积分阶数。节点坐标是一个 4×2 矩阵每行是一个节点的 x、y 坐标顺序必须是单元局部顺序通常逆时针排列。材料参数是弹性模量 E 和泊松比 nu。问题类型用字符串区分平面应力还是平面应变。积分阶数默认用 2表示每个方向两个高斯点。这样一个函数就能覆盖两种常见工况而不需要每换一种材料或单元形状就重新写一遍。3.2 完整代码展示下面这段是完整代码保存为Quad4Stiffness.m就能用。我把注释写得很详细方便你对照公式一行行看。function Ke Quad4Stiffness(nodes, E, nu, t, problemType, ngp) % Quad4Stiffness 计算平面四节点等参四边形单元刚度矩阵 % nodes: 4x2 矩阵每行为一个节点坐标 [x, y]按逆时针排列 % E: 弹性模量 % nu: 泊松比 % t: 厚度 % problemType: planeStress 或 planeStrain % ngp: 每个方向高斯积分点数取2或3 % Ke: 8x8 单元刚度矩阵 % % 示例: % nodes [0 0; 1 0; 1 1; 0 1]; % Ke Quad4Stiffness(nodes, 200e9, 0.3, 0.01, planeStress, 2); % ---------- 1. 材料弹性矩阵 D ---------- if strcmpi(problemType, planeStress) D E/(1 - nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; elseif strcmpi(problemType, planeStrain) fac E*(1-nu)/((1nu)*(1-2*nu)); D fac * [1, nu/(1-nu), 0; nu/(1-nu), 1, 0; 0, 0, (1-2*nu)/(2*(1-nu))]; else error(problemType 必须是 planeStress 或 planeStrain); end % ---------- 2. 高斯积分点和权重 ---------- if ngp 2 gpt [-0.577350269189626, 0.577350269189626]; gw [1.0, 1.0]; elseif ngp 3 gpt [-0.774596669241483, 0.0, 0.774596669241483]; gw [0.555555555555556, 0.888888888888889, 0.555555555555556]; else error(ngp 目前只支持 2 或 3); end % ---------- 3. 初始化单元刚度矩阵 ---------- Ke zeros(8, 8); % ---------- 4. 双重循环完成高斯积分 ---------- for i 1:ngp for j 1:ngp xi gpt(i); eta gpt(j); % 形函数对自然坐标的偏导数结果是 1x4 行向量 dN_dxi 0.25 * [-(1-eta), (1-eta), (1eta), -(1eta)]; dN_deta 0.25 * [-(1-xi), -(1xi), (1xi), (1-xi)]; % Jacobian 矩阵: 2x4 乘 4x2 得到 2x2 J [dN_dxi; dN_deta] * nodes; detJ det(J); % 检查行列式防止节点顺序错误 if detJ 0 error([det(J)%f 0请检查节点顺序是否为逆时针且单元形状合法], detJ); end invJ inv(J); % 形函数对物理坐标的偏导数 dN_dx invJ(1,1)*dN_dxi invJ(1,2)*dN_deta; dN_dy invJ(2,1)*dN_dxi invJ(2,2)*dN_deta; % 构造应变-位移矩阵 B (3x8) B zeros(3, 8); for n 1:4 B(1, 2*n-1) dN_dx(n); B(2, 2*n) dN_dy(n); B(3, 2*n-1) dN_dy(n); B(3, 2*n) dN_dx(n); end % 累加高斯积分 Ke Ke B * D * B * detJ * gw(i) * gw(j); end end % ---------- 5. 乘厚度 ---------- Ke Ke * t; end3.3 关键代码逐行解读这段代码里最容易被忽略的是 dN_dx 那两行。你可能觉得形函数对自然坐标的导数已经算出来了乘以一个逆矩阵自然是物理坐标导数但写代码时很容易把 invJ 的行列弄反。dN_dx 和 dN_dy 应该这样对应∂N/∂x J⁻¹(1,1)·∂N/∂ξ J⁻¹(1,2)·∂N/∂η∂N/∂y J⁻¹(2,1)·∂N/∂ξ J⁻¹(2,2)·∂N/∂η我第一次写的时候习惯性地把 invJ 的行和列写反结果出来的刚度矩阵怎么验都不对最后盯了半天才发现这里错了。建议你在调试阶段加一句disp(J); disp(invJ);拿规则正方形单元手动算一下 J再对照代码输出能很快确认这部分有没有写对。B 矩阵的构造也容易出错。我们的自由度排列是每个节点先 x 后 y即 d [u1, v1, u2, v2, u3, v3, u4, v4]。因此 B 矩阵第一行只在奇数编号的列上有值对应 εx 与水平位移的关系第二行在偶数编号的列上有值对应 εy第三行是剪应变两种自由度都参与。代码里的循环把每对自由度都同步填进去这样总自由度顺序就统一了。还有一个细节我特意让 detJ 小于等于零时直接报错。这在调试初期特别有用。你在命令行用Ke Quad4Stiffness(...)如果节点顺序写成了顺时针程序会立刻告诉你是单元定义问题而不是给你一个看起来“能算但结果全错”的矩阵。等熟悉以后可以考虑把这个检查去掉以节省一点运行时间但个人建议保留因为这个错误太容易犯了。4. 会组装才算会应用把单刚拼成整体刚度矩阵4.1 自由度编号与组装逻辑单个四边形单元有两个自由度所以单刚是 8×8。但实际模型有成百上千个节点总自由度是 2 乘以节点总数。把每个单元的 8×8 矩阵放进整个 K 矩阵的合适位置这个“放”的过程就是集成或者叫组装。关键是自由度编号。如果你的全局节点编号是已知的比如第 n 个节点的两个自由度分别是 2n-1 和 2n那么某单元第 k 个节点在全局里的自由度就是 2·elem(k)-1 和 2·elem(k)。组装时把单刚矩阵的每一行每一列对应到这些全局自由度上相加进去就行。我提供一段通用的组装函数它遍历所有单元逐个调用单刚函数再累加到全局矩阵。保存为assembleGlobal.mfunction K assembleGlobal(nodes, elems, E, nu, t, problemType, ngp) % 组装总体刚度矩阵 % nodes: Nn x 2全局节点坐标 % elems: Ne x 4每个单元四个节点编号按逆时针排列 % 其余参数同 Quad4Stiffness % K: 2*Nn x 2*Nn 总体刚度矩阵 nNode size(nodes, 1); K zeros(2*nNode, 2*nNode); nElem size(elems, 1); for e 1:nElem idx elems(e, :); Ke Quad4Stiffness(nodes(idx, :), E, nu, t, problemType, ngp); % 当前单元8个自由度在全局编号中的位置 dofMap zeros(1, 8); for k 1:4 dofMap(2*k-1) 2*idx(k) - 1; dofMap(2*k) 2*idx(k); end K(dofMap, dofMap) K(dofMap, dofMap) Ke; end end4.2 组装代码与边界条件处理组装好了 K接下来就是施加边界条件和求解。这一步和标题的“单元刚度矩阵”关系不大但如果不做很多人会卡在这里到底怎么验证单刚算得对不对。给一个小例子。一个 2×2 的四边形网格共 9 个节点右端受拉左侧固定。节点布局是每边 0.5 个单位% 单轴拉伸验证2x2 网格 nodes [0, 0; 0.5, 0; 1, 0; 0, 0.5; 0.5, 0.5; 1, 0.5; 0, 1; 0.5, 1; 1, 1]; elems [1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8]; E 200e9; nu 0.3; t 0.01; K assembleGlobal(nodes, elems, E, nu, t, planeStress, 2); % 边界条件左侧节点 1,4,7 固定 x、y fixedDofs [1, 2, 7, 8, 13, 14]; freeDofs setdiff(1:18, fixedDofs); % 右端节点 3,6,9 各施加 F/3 的 x 向集中力 Fvec zeros(18, 1); Fvec([5, 11, 17]) 1000/3; % 求解 u zeros(18, 1); u(freeDofs) K(freeDofs, freeDofs) \ Fvec(freeDofs); % 理论解应力 F/(t*L)应变 应力/E位移 应变*L theoryUx 1000 / (t * 200e9); % L 1 fprintf(右端节点x位移 %.6e m\n, u(5)); fprintf(理论x位移 %.6e m\n, theoryUx);如果你跑一下会发现数值结果和理论解一致差别在浮点误差范围内。这里我特意用了 2×2 网格而不是单单元网格就是为了演示组装过程和多个单刚叠加的效果。每一步都符合直觉左边固定右边拉单元只在 x 方向被拉长y 方向因为泊松效应略微收缩但左端约束了 y 方向位移结果中的 y 位移表现正常。边界条件处理的关键是只对 freeDofs 求解避免 K 奇异。任何有限元模型如果没有施加足够的约束K 都会是奇异的MATLAB 解出来要么是 NaN要么是巨大无比的数字。这个问题很常见不要怀疑你的刚度矩阵算错了先检查是不是边界条件漏了。5. 算例验证贴数据说话5.1 单轴拉伸与材料力学解对照前面代码里的单轴拉伸例子已经足够说明问题了。我再把数字整理一下方便你理解。模型尺寸是 1m×1m厚度 0.01m弹性模量 200GPa总拉力 1000N。材料力学解应力 σ F/(t·L) 1000/(0.01×1) 100000 Pa 0.1MPa。应变 ε σ/E 100000/200e9 5e-7。因为杆长 1m位移就是 5e-7m。用程序跑出来的右端节点 x 位移就是这个值浮点误差不超过 1e-18 量级。这个验证虽然简单但非常有价值。它证明了几件事形函数写对了Jacobian 变换没有算错B 矩阵构造正确高斯积分权重和积分点没有选错组装逻辑无误。如果你在这上面卡住了比如结果差 100 倍先检查是不是忘了乘厚度 t如果结果差一个等于 1e-6 左右的常数十有八九是单位制混乱造成的。有限元代码里单位制不统一是一切错误的根源我在写代码前会先把单位写在文件头注释里避免三天后自己都忘了这套模型用的是米制还是毫米制。5.2 结果可视化与后处理单元刚度矩阵本身看不见摸不着但求解之后的位移场可以画出来这能帮你直观确认单元没有翻转、没有畸变、边界条件作用方式对不对。MATLAB 里用patch函数画变形前后的网格非常方便scale 1e8; % 位移放大系数让变形可见 figure; hold on; for e 1:size(elems, 1) idx elems(e, :); % 未变形 patch(nodes(idx, 1), nodes(idx, 2), w, EdgeColor, k); % 变形后 un zeros(4,1); vn zeros(4,1); for k 1:4 un(k) nodes(idx(k), 1) scale * u(2*idx(k)-1); vn(k) nodes(idx(k), 2) scale * u(2*idx(k)); end patch(un, vn, none, EdgeColor, r, LineWidth, 1.5); end axis equal; grid on; title(变形前后对比红色为放大后的变形);跑完上述例子你会看到四个白色方形组成原来的网格红色变形后网格在 x 方向整体伸长y 方向因为泊松效应略微变窄。这个画面能给你很强的信心代码逻辑没问题。后处理里还有一个常用动作由节点位移计算单元应力。公式是 σ D·B·d_e其中 d_e 是单元位移向量。提取出来算一下就能输出每个单元的应力分量。有兴趣的话自己加上去不建议在单刚函数里做后处理保持单刚函数的纯净性这样它只负责返回 Ke其它功能模块化更清晰。6. 常见问题与排查技巧实录6.1 节点编号错乱导致行列式为负这是我见过最多的错误。等参四边形单元要求四个节点按逆时针顺序传入。如果你按顺时针传或中间两个节点交换位置Jacobian 行列式会变成负数甚至为零单刚矩阵全部错乱。最典型的错误是用“先下到上、再左到右”的直觉去写单元定义。网格里读出来的编号未必是单元局部顺序必须显式检查每个单元四个节点是不是在单元内部按逆时针排列。我加在代码里的if detJ 0 error(...)就是为此设计的。如果报了这个错先别改代码去排查该单元的节点顺序。排查方法可以这样画一下该单元的四个节点和连线顺序。代码里加一段plot(nodes([1:4,1],1), nodes([1:4,1],2), o-);如果是顺时针你会看到路径首尾相连但方向不对交换第 2 和第 4 个节点就能修正。6.2 平面应力与平面应变选错平面应力和平面应变的 D 矩阵不一样这是初学者容易忽略的坑。平面应力假设厚度方向应力为零适用于薄板类结构比如钢板、膜结构。平面应变假设厚度方向应变为零适用于厚截面或长条形结构比如大坝、隧道截面、长管的横截面分析。两种 D 矩阵形式在代码里已经分别写清楚了。如果你把一个平面应变模型当作平面应力算刚度会明显偏小位移就会偏大。反过来把薄板当平面应变算结果又偏保守。工程里必须一开始就确定清楚不能只改材料参数而忘了改问题类型。判断方法很朴素结构 z 方向尺寸远小于另外两个方向用 planeStress远大于另外两个方向用 planeStrain。拿不准的时候就查你参考的商业软件里用的是哪种单元比如平面应力单元和平面应变单元通常是不同的单元类型。6.3 高斯积分点数该选多少四节点等参四边形单元用 2×2 高斯积分是标准做法对规则形状单元是精确积分对轻微畸变单元精度也够。3×3 高斯积分能进一步降低畸变单元带来的误差但会增加约 2.25 倍的计算量。如果你在调试一个特别扭曲的网格可以先用 2×2 算再用 3×3 算一遍比较两个结果的差异。如果差异很小说明积分点足够差异很大说明网格质量太差问题不在于积分点数量而在于该修复网格了。我踩过的坑一度认为把网格加密就能提高精度但忽略积分点数不够这个前提。用 2×2 高斯积分加密网格时每个单元内部应变是线性变化的加密网格会让结果趋向收敛。但如果单元本身畸变严重比如一个内角接近 180 度的四边形再加密也只是让坏单元更小不会让映射关系变好。此时更明智的做法是重新划分网格或者在畸变区域改用三角形单元。6.4 整体刚度矩阵奇异怎么办K 奇异的表现是求解命令返回 NaN 或者矩阵接近奇异警告。最常见原因是约束不足。平面问题至少要约束三个自由度而且这三个自由度不能都在一条直线上。比如你只约束了左边两个节点的 x 和 y 自由度但其实约束总数是不够的——如果约束的三个自由度出现共线模型仍可能发生“转绕”刚体模式K 还是奇异。排查方法是检查rank(K)。如果等于 2 倍节点数减 3说明没有刚体模式可以求解。如果小于这个数会有刚体位移。还有一个经验把稀疏矩阵可视化一下spy(K)可以看到刚度矩阵的带宽结构如果出现某些行完全没有非零元素多半是某个节点没有参与任何单元也就是孤立节点。这也会导致奇异。另外如果你的模型确实约束足够却仍然奇异检查单元之间是否真正连在了一起。比如组装时节点编号写错两个相邻单元根本没有共享正确节点导致连接处自由度不连通。这个问题在大型网格里极难靠肉眼发现我的习惯是每次组装完抽查几个单元的 dofMap确认相邻单元的公共边对应同一个全局节点编号。7. 最后分享一点我的经验这套代码本身并不长但吃透它等于打通了有限元编程最关键的环节。我个人的经验是不要满足于“能跑出结果”而是拿它不断做对照实验把网格从 1×1 加密到 4×4观察位移收敛趋势把单元某个节点手动挪动一点看刚度矩阵和结果如何响应把平面应力换成平面应变估计结果会向哪个方向变化。这些实验每个都花不了几分钟但能帮你建立非常扎实的“有限元直觉”后面学更高阶单元、非线性分析都会顺得多。一个额外的小技巧把 Quad4Stiffness 和 assembleGlobal 分别保存成独立的 .m 文件再用一个测试脚本统一调用。不要把所有代码塞进一个脚本里否则改单元参数要反复滚动页面调试效率会明显下降。函数越独立越容易复用。等参四边形的这条路走通之后下一步可以尝试八节点曲边四边形单元或者六面体等参单元逻辑完全相同只是形函数和自由度数量变多了。到时候再翻有限元教材你会发现自己已经能主动理解它而不是被动接受它了。