企业官网建设流程全解析

很多刚接触有限元的同学都有同一种感觉：看一维杆单元的刚度矩阵推导还能跟上，一到二维平面问题就懵了。尤其是遇上四边形单元——等参变换、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-η)/4
N₂ = (1+ξ)(1-η)/4
N₃ = (1+ξ)(1+η)/4
N₄ = (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 0
0 ∂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)/((1+nu)*(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), (1+eta), -(1+eta)]; dN_deta = 0.25 * [-(1-xi), -(1+xi), (1+xi), (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; end

3.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.m：

function 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 end

4.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 文件，再用一个测试脚本统一调用。不要把所有代码塞进一个脚本里，否则改单元参数要反复滚动页面，调试效率会明显下降。函数越独立，越容易复用。

等参四边形的这条路走通之后，下一步可以尝试八节点曲边四边形单元，或者六面体等参单元，逻辑完全相同，只是形函数和自由度数量变多了。到时候再翻有限元教材，你会发现自己已经能主动理解它，而不是被动接受它了。