前几篇笔记我把平面三角形单元撸了一遍今天轮到平面四边形单元。很多初学者会有个错觉觉得三角形单元既然能剖分任意形状那为什么还要折腾四边形等你真正去算一个带转角的悬臂梁、一个孔边应力集中的平板或者网格稍微歪斜一点你会发现三角形单元的低阶形式真的不太够用。平面四边形单元是通向等参元、高阶单元、乃至实体单元的重要一环理解它的实现等于打通了有限元程序的任督二脉。这篇笔记我会沿着我自己写代码时的实际顺序展开先讲清楚为什么四边形值得用再把形函数、等参变换、雅可比矩阵、高斯积分、刚度矩阵组装这条线完整走一遍最后附上可直接运行的Python代码和验证算例。适合正在啃有限元源码、补课数值积分、或者想把教科书公式变成可运行程序的人。所有代码用Python加NumPy就能跑跑完你对“刚度矩阵到底怎么组出来”这件事会非常有底。1. 为什么要换平面四边形单元从三角形升级的逻辑1.1 三角单元与四边形的真实差距三角形常应变单元CST的问题在于它在每个单元内部只能描述常量应变场位移场是线性的。这意味着一个弯曲变形要用大量小三角形去逼近网格稍微粗一点结果就偏刚。你用CST算一个悬臂梁甚至会出现让人抓狂的“锁死”现象——位移解明显偏小。而平面四边形单元准确说是Q4双线性单元位移场在自然坐标下同时沿两个方向做线性插值还带一个ξ·η交叉项。交叉项的存在让单元内部能描述部分弯曲模式同样的网格密度下精度比CST高一个台阶。更重要的一点是四边形单元的映射是“等参”的它允许单元是任意直边四边形而不是必须矩形。这样在圆孔、倒角、过渡区域你可以用较少的单元贴合几何边界同时还能维持较好的单元形状。实际剖网格时四边形网格的单元数量通常能比三角形减少一半以上求解时间自然就降下来了。商业前处理软件里默认优先划分六面体网格二维问题优先四边形道理就在这里。1.2 等参单元的设计哲学等参变换听起来玄乎本质就是一句话把自然坐标系里一个规则的方形单元通过形函数映射到物理坐标系的一个任意四边形上并且位移插值函数和几何插值函数使用完全相同的形函数。你不需要为每个斜歪的四边形单独推导公式只需计算一次标准单元上的形函数及其导数其余交给雅可比矩阵完成坐标转换。这么做最大的好处是统一和通用。今天你写的是平面四节点四边形明天想升级成八节点曲边四边形只要换一套形函数改一下积分循环剩余代码框架完全不用动。我最初学有限元时就是被这种“以不变应万变”的设计给打动的真正理解了等参之后后面看实体单元、板壳单元都会顺很多。2. 形函数与等参变换推导先把数学底子打牢2.1 自然坐标下的形函数与坐标映射四节点四边形单元的母单元是自然坐标(ξ, η)下从-1到1的正方形四个节点按逆时针排列节点1对应(-1,-1)节点2对应(1,-1)节点3对应(1,1)节点4对应(-1,1)。四个形函数分别是N1 (1-ξ)(1-η)/4 N2 (1ξ)(1-η)/4 N3 (1ξ)(1η)/4 N4 (1-ξ)(1η)/4这套形函数满足两个关键性质在自身节点处取1在其他节点处取0任意一点的四个形函数之和恒等于1。第二个性质保证了刚体位移的模拟少了它单元就会有“空洞”这是自我检查时最容易踩的坑。我们再看它在物理坐标系里的作用假设单元的四个顶点坐标为(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)那么单元内任意一点坐标可以写成x N1·x1 N2·x2 N3·x3 N4·x4 y N1·y1 N2·y2 N3·y3 N4·y4这就是等参映射。因为形函数在边界上退化为线性函数比如η-1那条边只有N1和N2起作用所以直边四边形之间可以严格无缝拼接不会像高阶曲边单元那样协调性要求更高。2.2 雅可比矩阵从自然坐标到物理坐标的桥梁有了坐标映射下一步需要求位移对物理坐标的导数而形函数本身是自然坐标的函数。根据链式法则我们要通过雅可比矩阵做一次“换元”。雅可比矩阵定义为J [[∂x/∂ξ, ∂y/∂ξ], [∂x/∂η, ∂y/∂η]]每一项都可以由坐标节点值和形函数导数直接算出例如∂x/∂ξ ∂N1/∂ξ·x1 ... ∂N4/∂ξ·x4。编程时就是拿形函数导数矩阵与坐标矩阵做一次矩阵乘法。雅可比矩阵的行列式detJ有明确的几何含义它代表从母单元到物理单元的面积放大系数。高斯积分里每一项都要乘以detJ本质上是在做积分换元时补上面积变化。当单元是矩形时detJ是常数当单元扭曲成梯形或菱形时detJ在单元内部不再是常数被积函数不再是简单的多项式数值积分会引入截断误差。这也是为什么网格质量差的单元算出来精度会明显下降。我写代码时会在每个积分点强行检查detJ是否大于0。如果detJ小于等于0要么节点顺序是顺时针要么单元已经严重内凹。这种情况必须尽早暴露否则后面整个刚度矩阵是错的排查起来非常折腾。3. 单刚矩阵的高斯积分实现从公式到可运行的代码3.1 应变矩阵与本构矩阵的组装逻辑平面问题的应变分量有三个εxx、εyy、γxy。每个节点有两个位移自由度四节点单元一共有8个自由度所以几何矩阵B是3行8列。它的组装方式是每个节点对应对角位置放两个“子块”B[:, 2i] [∂Ni/∂x, 0, ∂Ni/∂y]的转置 B[:, 2i1] [0, ∂Ni/∂y, ∂Ni/∂x]的转置这个顺序必须和全局自由度编号方式保持一致我习惯每个节点先x方向自由度、再y方向自由度所有节点编号从小到大排列。B矩阵把节点位移映射成单元内应变应变 B·节点位移。本构矩阵D则把应变映射成应力应力 D·应变。平面应力问题和平面应变问题的D矩阵不同平面应力对应薄板D E/(1-ν²)乘以一个3×3矩阵平面应变对应厚截面或无限长结构D要换成带(1-2ν)分母的形式。程序里我会加一个布尔参数切换这两种模式免得算的时候记错公式。单元刚度矩阵的通用形式是k ∫ Bᵀ·D·B·t dA其中t是厚度。被积矩阵本身没有解析积分需要用数值积分逐点累加。3.2 高斯积分为什么四个点就够高斯积分的基本思想是用若干积分点上函数值的加权和近似代替连续积分。n点高斯积分最多能精确积分2n-1阶多项式所以2×2共4个积分点理论上有3阶代数精度。对于矩形单元单刚被积函数中形函数导数乘积最高是二次项所以2×2高斯积分实际上是精确的对扭曲单元由于detJ不再是常数被积函数不再是多项式但2×2仍然是很有效的近似。为什么不用1个积分点因为1点积分精度太低且会引入零能模式。为什么不用3×3被积函数的真实“困难”不在多项式阶数而在几何扭曲带来的非线性增加积分点对这种困难的改善有限反而增加计算量。实际工程里Q4单元配2×2高斯积分是绝对的主流我在代码里也默认使用2×2。3.3 element_stiffness的完整实现下面是我在实践笔记里写的单刚矩阵计算函数代码风格偏向教学但结构可以直接用在求解器里。import numpy as np def gauss_2x2(): 返回2x2高斯积分点和权重 g 1.0 / np.sqrt(3.0) return [(-g, -g, 1.0), (g, -g, 1.0), (g, g, 1.0), (-g, g, 1.0)] def shape_derivatives(xi, eta): 四节点四边形单元形函数对自然坐标的导数返回(4,2)矩阵 dN_dxi np.array([ [-(1.0 - eta) / 4.0, -(1.0 - xi) / 4.0], [(1.0 - eta) / 4.0, -(1.0 xi) / 4.0], [(1.0 eta) / 4.0, (1.0 xi) / 4.0], [-(1.0 eta) / 4.0, (1.0 - xi) / 4.0] ]) return dN_dxi def element_stiffness(coords, E2.1e11, nu0.3, thickness0.05, plane_stressTrue): 计算Q4单元单刚矩阵。coords为(4,2)数组按逆时针顺序排列 if plane_stress: D E / (1.0 - nu**2) * np.array([ [1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, (1.0 - nu) / 2.0] ]) else: D E / ((1.0 nu) * (1.0 - 2.0 * nu)) * np.array([ [1.0 - nu, nu, 0.0], [nu, 1.0 - nu, 0.0], [0.0, 0.0, (1.0 - 2.0 * nu) / 2.0] ]) k np.zeros((8, 8)) for xi, eta, w in gauss_2x2(): dN_dxi shape_derivatives(xi, eta) J dN_dxi.T coords # 2x2 Jacobian detJ np.linalg.det(J) if detJ 0: raise ValueError(detJ 0单元畸变或节点顺序错误) # 物理坐标系下的形函数导数每个积分点做一次局部变换 dN_dx np.zeros_like(dN_dxi) for i in range(4): dN_dx[i, :] np.linalg.solve(J, dN_dxi[i, :]) B np.zeros((3, 8)) for i in range(4): Nx dN_dx[i, 0] Ny dN_dx[i, 1] B[:, 2*i:2*i2] [[Nx, 0.0], [0.0, Ny], [Ny, Nx]] k B.T D B * detJ * w * thickness return k这段代码有三个地方值得说明。第一np.linalg.solve(J, dN_dxi[i, :])比直接求逆矩阵更稳定而且语义上就是解线性方程组把自然坐标导数换算成物理坐标导数。第二B矩阵组装时直接对列切片赋值利用numpy的广播避免了逐元素循环效率在单元数量不大时完全够用。第三单刚矩阵一定是8×8对称矩阵行列和近似为零调试时如果这两条不满足九成是积分点循环或形函数导数写错了。4. 全局刚度组装与算例验证让代码跑出正确结果4.1 自由度编号与总体刚度矩阵有限元程序的全局组装本质上是一个“由单元到总体”的索引累加过程。每个节点有ux和uy两个自由度所以节点编号n对应的全局自由度是2n和2n1从0开始计数在Python里最自然。我建议先把所有节点的坐标和单元表构造出来再遍历单元表对每个单元计算单刚然后把单刚里的8×8个元素“散射”进全局刚度矩阵对应的位置。为了提高可读性下面的代码没有做稀疏优化只是用稠密矩阵存储。几十万个自由度的规模不建议这样写但作为教学和中小算例完全没问题。真要追求性能把K声明成scipy.sparse.lil_matrix再转csr改造非常简单。4.2 边界条件与求解策略约束处理的常用方法有置大数法和划行划列法。置大数法实现简单但会引入数值误差大数取多少需要经验划行划列法最干净把固定自由度的行和列删掉解缩减后的方程组再把位移还原到完整向量。我用Python做这类小算例基本都选划行划列法用布尔掩码一次性取出行列索引。def assemble_global(elements, coords, E, nu, thickness, plane_stressTrue): n_nodes coords.shape[0] n_dof 2 * n_nodes K np.zeros((n_dof, n_dof)) for elem_nodes in elements: k element_stiffness(coords[elem_nodes], E, nu, thickness, plane_stress) dofs [] for n in elem_nodes: dofs.append(2 * n) dofs.append(2 * n 1) for i in range(8): for j in range(8): K[dofs[i], dofs[j]] k[i, j] return K def apply_bc(K, F, fixed_dofs): 划行划列法施加零位移约束 active np.ones(len(F), dtypebool) active[fixed_dofs] False Kaa K[np.ix_(active, active)] Fa F[active] ua np.linalg.solve(Kaa, Fa) u np.zeros(len(F)) u[active] ua return u这里有一个容易忽略的细节固定自由度上如果有外载荷划行划列后这部分载荷会被直接丢掉但求解的位移不受影响只是你没法从求解结果里直接得到支座反力。要算反力应该先解出所有自由度位移再回到原方程计算R K·u - F只看固定自由度位置的分量。4.3 均匀受拉板验证算例为了检验程序我构造一个最简单的验证算例长L2m、高H1m的矩形板厚度t0.05m左侧边完全固定右侧边施加均匀水平拉力。材料取钢材常用值弹性模量E2.1×10^11 Pa泊松比nu0.3。网格我划分成8×4个单元右端总拉力P100000N均分到右端5个节点上每个节点施加20000N。这个问题的精确解非常容易得到均匀拉伸下右端水平位移u P·L / (E·A)AH·t0.05m²算出来大约是1.9048×10^-5m。由于Q4单元对均匀应力场具有精确再现能力只要边界条件加得对数值解应该和解析解几乎完全相等。我实际跑出来的结果右端所有节点水平位移都在1.9048×10^-5m附近误差不到1e-6量级这就是单元代码正确的强信号。如果你跑出来偏差很大优先检查力向量施加是否漏了某个节点或者约束是否封住了不该封的自由度。构建网格的代码也很简单矩形区域按行列生成节点每个小格生成一个逆时针单元def build_rect_mesh(nx, ny, L, H): coords [] for j in range(ny 1): for i in range(nx 1): coords.append([i * L / nx, j * H / ny]) coords np.array(coords) elements [] for j in range(ny): for i in range(nx): n0 j * (nx 1) i n1 n0 1 n2 n1 (nx 1) n3 n0 (nx 1) elements.append([n0, n1, n2, n3]) return coords, np.array(elements)左侧固定节点是j0到ny的所有行中索引j*(nx1)的节点右端加载节点是索引j*(nx1)nx的节点。把这些过滤出来调用assemble_global和apply_bc就能完成求解。这个原型我给过不止一个同事看他们都反馈跑通这个例子后再去读那些推导细节思路会清晰很多。5. 积分方案与单元畸变隐藏的精度杀手5.1 什么是减缩积分什么是沙漏模式全积分指的是让积分点数足以精确积分刚度矩阵中的多项式项对于Q4单元就是2×2高斯积分。很多人会问“2×2只有4个点能叫全积分吗”能因为双线性单元的形函数导数是常数或线性2×2高斯积分的代数精度已经足够处理多项式部分了。真正意义上的“减缩积分”是只取1个积分点也就是单元中心点。1点积分会带来一个经典问题零能模式也叫沙漏模式。某些特定位移形态下所有积分点上的应变恰好为零但位移并不是刚体位移导致单元刚度矩阵存在额外的零特征值。表现出来就是计算结果出现波浪形的棋盘状位移场网格越密越明显。工程上为了对付剪切锁死有时候会故意用减缩积分但减缩积分单元必须配沙漏控制。自己写程序时没有沙漏控制手段建议老老实实用2×2积分。为什么全积分Q4单元在弯曲问题里会“锁死”因为纯弯曲状态的截面应变沿厚度方向应该是线性变化的而Q4单元在物理坐标下应变分量是常数或简单线性无法精确表达弯曲变形整个单元偏刚。细分网格能缓解但效率变低。我个人的经验是如果问题以弯曲为主与其纠结Q4的减缩积分不如直接换Q8八节点四边形或者升级单元阶次。5.2 单元畸变如何悄悄毁掉精度等参变换允许单元是任意四边形但不代表你可以随便画。当单元偏离矩形形状雅可比矩阵变为非常数原本精确的多项式积分不再精确误差会随畸变程度迅速上升。最典型的坏形状有两种一是单元内角接近0度或180度二是内部出现“凹”形即某个内角大于180度后者直接导致detJ在某些区域小于零程序会报错。我在实际画网格时给自己定了几条粗规则单元最长边与最短边之比尽量不超过5每个内角控制在45度到135度之间关键受力区域的单元尽量不要是严重歪斜的平行四边形。这些规则没有硬性理论门槛都是前处理经验但对计算结果的可信度影响很大。你可以在程序中算完detJ后统计最小值和最大值之比作为网格质量的一个粗估指标。5.3 工程实践中的单元质量检查拿到一个现成网格时我建议先做三个快速体检。第一遍历所有单元检查节点顺序是否逆时针一旦出现顺时针单元立刻修正为逆时针否则detJ恒为负。第二计算每个单元的雅可比行列式在四个积分点的值如果有任何一个小于等于零说明单元已经严重扭曲最好重新剖分。第三检查网格边界上有没有重复边或悬挂节点这类问题在手工网格里频繁出现结果往往是整体刚度矩阵奇异。这些检查听起来繁琐但写成脚本也就是几行循环的事。我在程序里把detJ小于某个阈值的单元列表打印出来定位到具体坐标再回到前处理软件里改网格效率非常高。6. 常见问题排查与实战心得6.1 典型问题速查表现象可能原因处理方向detJ小于等于0节点顺序为顺时针、单元内凹、节点坐标重复检查单元节点顺序插值或重构网格求解时报矩阵奇异约束不足、存在刚体位移、单元表有孤立节点检查固定边界去掉无约束刚体模式结果出现棋盘状位移场减缩积分导致的沙漏模式改回2×2全积分或添加沙漏控制位移明显偏小单元弯曲锁死、载荷漏加、约束过度检查载荷向量改用更高阶单元结果随网格加密不收敛单元畸变严重、材料参数错误、边界条件不协调检查最小detJ核对材料参数和边界单刚矩阵对角线出现负数形函数或B矩阵组装错误专门打印单刚矩阵逐列核对节点自由度顺序表里的前两条是我见过最频繁的问题。特别是初学阶段网格是自己手工生成的话节点顺序非常容易弄反。我的建议是构造单元表时一律按“左下、右下、右上、左上”的顺序写这样天然是逆时针不用每次再去判断。6.2 调试代码的几个实用技巧调试有限元程序不能像调试普通业务代码那样看输出就完事你需要分层次验证。第一步先检查形函数把节点坐标代进去看形函数是否满足“本点为1他点为0”同时检查任意点四个形函数之和是否等于1。第二步检查雅可比矩阵对矩形单元detJ应该恒等于单元面积的四分之一一旦歪斜detJ在积分点间应当有变化但始终为正。第三步检查单刚矩阵给一个简单矩形单元打印8×8单刚行和列都应该近似为零因为刚体位移对应的应变能必须为零。还有一个非常实用的“补丁测试”用一个最小网格比如2×2个单元拼成一个矩形施加一个线性位移场到边界节点上然后检查内部节点的位移解是否精确复现同一个线性场。如果程序有BUG这个测试会第一时间暴露。我做任何新的单元开发都会先把补丁测试跑过再碰实际算例这个习惯帮我省掉了大量排错时间。6.3 我个人踩过坑之后的一些体会如果你发现自己写的Q4单元程序在均匀拉伸问题上都对但一到弯曲问题就偏刚先别急着怀疑代码写错这很可能是单元本身的锁死特性。这时需要判断的问题是网格是否已经细到工程可用范围如果不想加密网格就考虑换更高阶单元或非协调模式。Q4单元本身在工程中地位这么牢固并不是因为它一定能精确模拟弯曲而是它在网格适度加密后能得到可靠结果且计算成本很低。理解这一点你才不会在调试时钻牛角尖。最后再分享一个绕不开的小技巧所有节点坐标、自由度编号、单元节点顺序在进入计算之前建议先做一次图形化检查。我用matplotlib把网格画出来把节点编号和单元号标在图上肉眼扫一遍再开始组装。很多索引错误都是在图上直接看出来的靠打印数字反推反而慢得多。写有限元程序几何和索引的代码往往比数学公式更容易出错而这些错误又最磨人。养成可视化检查的习惯后这类问题基本能在十分钟内定位。