前阵子帮人排查一份二维平面应力的自研代码现象很奇怪单元刚度矩阵对称、量纲也对可一算悬臂梁挠度就比解析解小了三成多。最后定位到的原因跟高斯积分、材料参数都没关系——三个节点的局部编号本来是逆时针装配到某个区域时被写成了顺时针形函数里的坐标差项整体变号单元的 B 矩阵虽然还是常数矩阵但方向全反了。改一个符号结果立刻对上。这件事之后我一直在想为什么形函数这么容易写错答案很朴素——绝大多数人是把形函数当公式背下来的背的时候记住了N1 1 - ξ - η这种样子却说不清楚N1在三角形里到底长什么样、为什么在节点 1 处等于 1、在对边上等于 0。一旦坐标顺序、单元类型、节点数量发生变化背书式的记忆立刻失效。这篇要聊的“划线法”就是把这个东西从“背公式”拉回到“看得见”的一种思路。它的核心动作只有一个在三角形内部划线。对线性单元划的是平行于对边的直线对二次单元划的是双曲线。形函数的值就是这些线在单元里走的“位置刻度”。理解了这个几何图像三节点单元、六节点单元、甚至更高阶的单元你都能自己在纸上把它推出来而不是去翻手册。内容会从几何原理讲到代码落地包含完整的形函数推导、B 矩阵计算、刚度矩阵装配以及我在实际使用中总结的几套自查手段。适合正在手写有限元程序、或者被形函数编号问题反复折磨的读者。1. 形函数为什么不能靠背从一次符号错位说起1.1 公式背得越熟出错时越难定位先复盘开头那个案例。原始的坐标处理代码里形函数的导数项是这样写的dNdx np.array([y1-y2, y2-y0, y0-y1]) / (2*A)这行代码本身没错但它隐含了一个前提节点局部编号0→1→2必须构成逆时针回路。当某片区域的网格在导出时被工具反向排列2*A变成了负数dNdx整体反号B矩阵也跟着反号。单元刚度Ke t·A·BᵀDB里A是负的B是负的BᵀDB里的负号被平方抵消所以Ke看上去依然“对称、正定、量纲正常”——这正是最难查的地方。如果当时脑子里有清晰的几何图像问题会简单得多∇N₁的方向应该垂直于对边边 2-3并指向节点 1模长等于1/h₁h₁是节点 1 到对边的垂距。写完之后拿这个判据一测方向反了立刻就能发现。这就是几何视角相对代数视角的价值它给了你一个独立于公式的物理判据。1.2 形函数必须同时满足的三条硬约束不管用什么方法构造一套合法的形函数必须同时满足下面三条。这三条不是建议是定义本身约束名称数学表达物理含义违反后的后果节点性Kronecker 性Nᵢ(xⱼ, yⱼ) δᵢⱼ节点 i 的单位位移只影响节点 i 本身边界条件加不上支座处位移不收敛单位分解性Σ Nᵢ ≡ 1刚体平移不产生任何应变单元无法表示刚体位移结果整体偏硬线性完备性Σ Nᵢxᵢ xΣ Nᵢyᵢ y常应变状态可被精确表示分片检验失败网格加密不收敛前两条容易理解第三条经常被忽略。它的本质要求是形函数空间必须包含所有一次多项式。这条性质直接决定了单元能不能通过分片检验patch test而分片检验是收敛性的必要条件——过不了这个检验网格再密也是白搭。六节点二次单元之所以“更高级”是因为它在满足这三条之外还能精确表示二次位移场所以对弯曲主导的问题表现更好。但注意高次完备性不是必须的一次完备性才是。1.3 “划线”划的到底是什么这里说清楚我的用法所谓划线法指的是把形函数Nᵢ理解成一个标量场然后去考察它的等值线。因为在三角形单元里Nᵢ的等值线形状极其规整而且有明确的几何意义。对三节点线性单元N₁的等值线是一族平行于对边2-3的直线。N₁ 0这条线就是整条对边N₁ 1退化成节点 1 一个点中间每一条平行线N₁的取值就是“当前线到对边的距离”除以“节点 1 到对边的距离”。这就是“划线”的字面含义——你在三角形里画一条平行于底边的线这条线上所有点的N₁值都一样。对六节点二次单元的顶点形函数N₁ 0的线不止一条而是两条平行于对边的直线一条是对边本身另一条穿过两条邻边的中点。这两条线把三角形切成三段中间一段N₁为负。这个图像一出来为什么顶点形函数要写成L₁(2L₁-1)就完全说得通了。对六节点二次单元的边中点形函数N₄ 0的线退化成两条边等值线则是以这两条边为渐近线的双曲线族整个场像一道沿边 1-2 隆起的山脊。不同教材对这套方法叫法不一有的叫几何法有的叫图形法有的直接归到面积坐标法里。名字不重要重要的是它把“形函数的数值”和“单元内的几何位置”绑定了。下面几章会把这套图像完整落地。2. 一条平行线切出的比例划线法的几何内核2.1 三角形里的相似关系给了一切考虑三角形123在内部任取一点P从P作一条平行于边 2-3 的直线交边 1-2 于Q、交边 1-3 于R。设P到边 2-3 的垂距为d₁节点 1 到边 2-3 的垂距为h₁。由相似三角形关系有|PQ| / |23| (h₁ - d₁) / h₁更重要的是P沿这条平行线移动时d₁完全不变。换句话说“到对边 2-3 的垂距”这个量在平行于 2-3 的整条直线上是常数。这正是一个标量场等值线的定义。于是我们可以直接定义N₁(P) d₁(P) / h₁三个性质随手验证P落在节点 1 时d₁ h₁N₁ 1P落在边 2-3 上时d₁ 0N₁ 0P沿任意平行线移动时N₁不变所以等值线确实是直线族。而且由于d₁(P)是P点坐标的一次函数点到直线的有向距离公式是关于坐标的线性式N₁天然就是坐标的线性函数——线性完备性自动满足。N₂、N₃同理分别用对边 3-1 和 1-2 的垂距来定义。三个形函数就这么“划”出来了一行代数推导都没用。2.2 梯度方向与模长最好用的验算判据上面的定义立刻给出一个特别实用的结论。因为Nᵢ的等值线平行于对边所以梯度∇Nᵢ必然垂直于对边又因为Nᵢ从对边上的 0 线性增长到节点上的 1走完全程的距离恰好是hᵢ所以|∇Nᵢ| 1 / hᵢ这两句话是免费的单元测试。写代码时算完dNdx和dNdy之后随手验证import numpy as np def check_gradient(coords): coords: (3,2) 逆时针节点坐标 dNdx, dNdy, A tri3_derivatives(coords) for i in range(3): j, k (i1) % 3, (i2) % 3 # 对边向量与高 e coords[k] - coords[j] L np.linalg.norm(e) h 2.0 * A / L # 梯度方向应垂直对边 grad np.array([dNdx[i], dNdy[i]]) cos_angle abs(np.dot(grad, e)) / (np.linalg.norm(grad) * L) assert cos_angle 1e-12, fN{i1} 梯度不垂直于对边 assert abs(np.linalg.norm(grad) - 1.0/h) 1e-12, fN{i1} 梯度模长错误这段检查只花几微秒但能挡住坐标顺序错、节点索引错、面积公式写错这三类最常见的问题。我在好几个项目里都把它放进了单元级单元测试里跑一次网格就全部覆盖。还有一个更隐蔽的用处|∇Nᵢ| 1/hᵢ说明单元越扁hᵢ越小梯度越大。当三角形某条高接近零时形函数梯度会急剧放大B 矩阵的条件数爆炸刚度矩阵出现病态。这就是畸形三角形必须控制纵横比的根本原因跟“网格质量”这种含糊说法比起来这个解释有明确的数学来源。2.3 面积坐标与划线法其实是同一件事有限元教材里最常见的写法是面积坐标Lᵢ Aᵢ / AAᵢ是P点与对边围成的子三角形面积A是单元总面积。这个定义和上面的垂距比是等价的因为A₁ (1/2)·|23|·d₁ A (1/2)·|23|·h₁两式相除A₁/A d₁/h₁ N₁。所以面积坐标的几何含义就是“垂距比”只不过借了面积这个中介。之所以教材偏爱面积坐标是因为它天然满足L₁ L₂ L₃ 1三个子三角形拼成整个三角形单位分解性一眼可见。我个人在做推导时用垂距比因为直观写代码时用面积坐标因为它对称、好扩展到高阶。两套语言随时切换不冲突。顺便把几种常见的构造方法摆在一起对比方便你按场景选择方法操作方式优势局限适用场景划线法几何法考察等值线形状与位置图像直观高次单元也能画出来手算表达式偏慢理解原理、验算、教学中推导待定系数法设N a bx cy代入节点条件解方程机械可靠不易漏项高次单元方程数暴涨一次单元、小规模单元的快速推导面积坐标法直接用Lᵢ构造多项式高度对称易于程序化需要额外记忆坐标变换公式二次及以上单元的代码实现张量积法一维形函数相乘程序通用性强三角形单元不适用四边形、六面体单元3. 三节点常应变单元把 N₁、N₂、N₃ 逐个划出来3.1 三次划线三个形函数按 2.1 的定义N₁由“到边 2-3 的垂距比”给出。写出坐标形式设节点坐标为(x₁,y₁)、(x₂,y₂)、(x₃,y₃)单元面积的两倍记为2A (x₂-x₁)(y₃-y₁) - (x₃-x₁)(y₂-y₁)逆时针编号时2A 0。边 2-3 的方程可以写成(y₂-y₃)x (x₃-x₂)y (x₂y₃ - x₃y₂) 0把节点 1 的坐标代进去得到常数项再整理就能得到标准形式N₁ [(y₂-y₃)x (x₃-x₂)y (x₂y₃ - x₃y₂)] / (2A)N₂、N₃按1→2→3→1循环替换下标即可。三个式子形式完全一样只有下标在轮转——这种对称性是划线法给出的直接结果因为三个形函数的构造过程本来就是一回事只是换了一条对边。把三式相加分子里的x项系数是(y₂-y₃)(y₃-y₁)(y₁-y₂) 0y项系数是(x₃-x₂)(x₁-x₃)(x₂-x₁) 0常数项是循环和也恰好等于2A。所以N₁N₂N₃ ≡ 1单位分解性自动成立。3.2 用待定系数法做交叉验证几何法推出来的式子我习惯用待定系数法再走一遍两路结果对上才敢往代码里塞。步骤很死板设N₁ a bx cy代入三个节点条件a b·x₁ c·y₁ 1 a b·x₂ c·y₂ 0 a b·x₃ c·y₃ 0这是一个 3×3 线性方程组。系数矩阵的行列式恰好等于2A也就是两倍面积几何上对应三个节点围成的平行四边形的面积。只要三点不共线行列式非零解存在且唯一。解得b (y₂-y₃)/(2A)c (x₃-x₂)/(2A)a (x₂y₃-x₃y₂)/(2A)和几何法的结果完全一致。这个交叉验证的价值不在于“确认公式对不对”而在于确认你的几何直觉没有跑偏。我在早期推导六节点单元时就是靠这个方式发现自己把某个边中点形函数的系数写成了 2 而不是 4。3.3 从形函数到 B 矩阵常数应变是怎么来的有了形函数的显式表达式求导是机械操作。因为Nᵢ关于x、y都是一次的导数就是常数∂N₁/∂x (y₂-y₃)/(2A) ∂N₁/∂y (x₃-x₂)/(2A) ∂N₂/∂x (y₃-y₁)/(2A) ∂N₂/∂y (x₁-x₃)/(2A) ∂N₃/∂x (y₁-y₂)/(2A) ∂N₃/∂y (x₂-x₁)/(2A)**注意这些导数里没有x和y。**这就是“常应变三角形单元CST”这个名字的来历B 矩阵在整个单元内是常数矩阵所以应变在单元内均匀分布。这个结论是形函数一次性的直接推论不需要额外证明。把导数按平面问题的标准排布组装成 3×6 的 B 矩阵B [ ∂N₁/∂x 0 ∂N₂/∂x 0 ∂N₃/∂x 0 ] [ 0 ∂N₁/∂y 0 ∂N₂/∂y 0 ∂N₃/∂y ] [ ∂N₁/∂y ∂N₁/∂x ∂N₂/∂y ∂N₂/∂x ∂N₃/∂y ∂N₃/∂x ]单元刚度矩阵就是Ke t·A·BᵀDB——因为被积函数是常数积分直接退化成“乘以面积”连高斯积分都不需要。这是 CST 单元最快的计算路径我在做大规模网格快速预估时经常用它先跑一遍量级。3.4 用 Python 把这条链路完整跑通下面这段代码是我常用的最小实现包含形函数导数、B 矩阵、刚度矩阵以及三条自检import numpy as np def tri3_derivatives(coords): coords: (3,2) 逆时针排列的节点坐标 返回 dNdx, dNdy, 面积A x coords[:, 0] y coords[:, 1] twoA (x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]) if twoA 0: raise ValueError(f单元面积非正 ({twoA/2:.4e})检查编号方向) dNdx np.array([y[1]-y[2], y[2]-y[0], y[0]-y[1]]) / twoA dNdy np.array([x[2]-x[1], x[0]-x[2], x[1]-x[0]]) / twoA return dNdx, dNdy, twoA / 2.0 def tri3_B(coords): dNdx, dNdy, A tri3_derivatives(coords) B np.zeros((3, 6)) B[0, 0::2] dNdx B[1, 1::2] dNdy B[2, 0::2] dNdy B[2, 1::2] dNdx return B, A def tri3_Ke(coords, D, thickness1.0): B, A tri3_B(coords) return thickness * A * (B.T D B) def D_plane_stress(E, nu): return E / (1 - nu**2) * np.array([ [1.0, nu, 0.0], [nu, 1.0, 0.0], [0.0, 0.0, (1 - nu) / 2] ]) def self_check(coords): dNdx, dNdy, A tri3_derivatives(coords) # 1. 单位分解梯度之和为零 assert abs(dNdx.sum()) 1e-12 and abs(dNdy.sum()) 1e-12 # 2. 线性完备Σ Nᵢxᵢ x assert abs((dNdx * coords[:, 0]).sum() - 1.0) 1e-12 assert abs((dNdy * coords[:, 1]).sum() - 1.0) 1e-12 # 3. 梯度模长等于高度的倒数 for i in range(3): j, k (i1) % 3, (i2) % 3 h 2.0 * A / np.linalg.norm(coords[k] - coords[j]) assert abs(np.hypot(dNdx[i], dNdy[i]) - 1.0/h) 1e-12 return True if __name__ __main__: coords np.array([[0.0, 0.0], [2.0, 0.0], [0.5, 1.5]]) print(自检通过:, self_check(coords)) D D_plane_stress(2.1e11, 0.3) Ke tri3_Ke(coords, D, thickness0.02) print(单元刚度矩阵对称性:, np.allclose(Ke, Ke.T)) print(对角元全正:, np.all(np.diag(Ke) 0))这三条自检的覆盖范围值得说一下单位分解检查能抓住系数写错和求和下标写错线性完备检查能抓住面积符号错误因为符号错了两边会同时变号但归一化后可能仍通过所以还需要配合第三条梯度模长检查覆盖了所有几何相关错误。实话说我后来所有单元的测试都是从这个模板扩展出去的。4. 六节点二次单元划线从直线变成双曲线4.1 边中点形函数的图像沿边隆起的一道山脊六节点三角形在三条边上各加一个中点节点编号习惯是 1、2、3 为顶点4、5、6 分别是 1-2、2-3、3-1 的中点。面积坐标下六个形函数是N₁ L₁(2L₁-1) N₂ L₂(2L₂-1) N₃ L₃(2L₃-1) N₄ 4L₁L₂ N₅ 4L₂L₃ N₆ 4L₃L₁先看N₄ 4L₁L₂。把垂距比代进去N₄ 4 · (d₁/h₁) · (d₂/h₂)也就是说N₄的值正比于“到边 2-3 的垂距”与“到边 3-1 的垂距”的乘积。N₄ 0的条件是d₁ 0或d₂ 0即点落在边 2-3 或边 3-1 上——这两条边正是节点 4 所在的边 1-2 的“邻边”。而N₄的等值线d₁·d₂ 常数在几何上是一族以边 2-3 和边 3-1 为渐近线的双曲线。当你从边 2-3 出发向边 1-2 移动时d₁增大N₄随之增大从边 3-1 出发也是同理。两条路径的乘积场在边 1-2 的中点节点 4达到峰值 1然后向内陆迅速衰减。整体形状像一道沿着边 1-2 隆起的山脊越远离这条边山脊高度越低。这个图像解释了一件常被误解的事N₄的影响范围不局限于边 1-2 附近它在整个单元内部都不为零只是衰减得很快。所以边中点节点的位移变化会通过N₄影响到远方顶点附近的位移场这是二次单元比线性单元“传得更远”的原因也是它刚度表现得“更柔”的原因。4.2 4L₁L₂ 里的系数 4 是怎么划出来的如果是第一次看这个表达式很容易疑惑为什么是 4 而不是别的数答案藏在“划线”的峰值条件里。我们要构造一个在边 1-2 上为抛物线、在节点 4 处等于 1 的形状函数。边 1-2 上L₃ 0所以L₁ L₂ 1。取最简单的乘积形式c·L₁L₂代入节点 4 的位置L₁ L₂ 1/2c · (1/2) · (1/2) 1 → c 4系数 4 就是这么来的。如果你换成L₁L₂L₃之类的形式在边上会恒为零因为L₃ ≡ 0整个函数在边上消失了无法满足节点性。所以乘积因子里只能包含在边上不为零的那两个坐标。同样的逻辑推N₅ 4L₂L₃节点 5 在边 2-3 上L₁ 0所以用L₂L₃N₆ 4L₃L₁。三个边中点形函数的形式完全对称只是取哪两个坐标在轮换。顺便验一下节点性N₄在节点 4 处等于 1在节点 1L₂ 0处等于 0在节点 2L₁ 0处等于 0在节点 5L₁ 0处等于 0在节点 6L₂ 0处等于 0在节点 3L₁ L₂ 0处等于 0。六个节点全对。这种“一眼验六个”的顺畅感只有在图像清楚的时候才有。4.3 顶点形函数为什么要“减一项”再看顶点形函数N₁ L₁(2L₁-1)。如果只写L₁它在顶点 1 处等于 1在边 2-3 上等于 0看起来已经满足了“顶点 1 的形函数”。但问题在于边中点节点 4、6 落在L₁ 1/2的位置上此时L₁ 1/2 ≠ 0节点性被破坏。所以必须引入一个修正让N₁在L₁ 1/2处也归零。二次多项式L₁(2L₁-1)的两个零点正好是L₁ 0和L₁ 1/2完美符合要求。回到划线图像N₁ 0的等值线有两条一条是对边 2-3L₁ 0另一条是平行于对边、穿过边 1-2 和边 3-1 中点的中位线L₁ 1/2。两条线把三角形切成三带靠近顶点 1 的窄带N₁ 0中间带N₁ 0靠近对边的带N₁ 0等等这个说法要修正一下——实际是L₁ 1/2时为正0 L₁ 1/2时为负L₁ 0处为 0。严格说f(L₁) 2L₁² - L₁在L₁ 1/4处取到最小值-1/8然后在L₁ 1/2处回到零之后单调上升到顶点处的 1。所以形函数在单元内部存在负值区最小值是 -0.125。这个负值区不是数学缺陷而是二次插值的必然结果要让形函数在边中点处精确归零它就必须在中间某段“凹下去”。负值区的代价是单元刚度矩阵可能出现负的非对角元虽然整体仍正定某些迭代求解器在这个区域收敛会慢一些。实际使用中不用刻意处理但心里要有数。4.4 二次单元的收益与代价把线性单元和二次单元放在一起对比能看得更清楚项目三节点线性单元六节点二次单元形函数等值线平行直线族直线族 双曲线族位移场次数一次二次应变场常数单元内均匀线性变化边界描述只能直线边可拟合曲线边边中点可偏移单元刚度计算常数被积函数一次乘法需要至少 3 点数值积分弯曲问题表现严重偏刚需极密网格明显改善但仍有剪切锁定倾向典型陷阱纵横比敏感单点积分导致沙漏、边中点误放二次单元真正值钱的地方在于能表示线性应变这直接改善了对弯曲主导问题的模拟精度。在线性单元里一个单元只能给一个常应变值弯曲梁必须靠很多单元“阶梯式”逼近曲率二次单元则可以线性变化用少得多的单元就能逼近。但二次单元不是免费的。最典型的坑是数值积分阶次Ke的被积函数是BᵀDB而B是线性的因为二次形函数求导后是一次所以被积函数是二次多项式至少需要二次精度的积分公式。如果用形心单点积分会直接导致刚度矩阵秩亏出现零能模式沙漏结构看起来软塌塌地变形而没有任何应变能增长。这个错误在自研代码里非常常见因为单点积分算得快跑小算例时还不容易暴露。5. 落到代码从形函数公式到可用的单元刚度5.1 节点编号约定要写死在文档里写单元代码的第一件事不是写公式是定编号约定。我的习惯是把下面三条写在文件头注释里任何改动都要同步更新局部节点 0、1、2 逆时针排列面积恒为正六节点单元的局部节点 3、4、5 分别位于边 (0,1)、(1,2)、(2,0) 的中点全局自由度排布为[u₀, v₀, u₁, v₁, ...]节点内u在前v在后。第二条尤其重要。六节点单元的中点在边上的位置组合方式有若干种不同教材和不同软件的习惯不一。一旦和网格生成器输出不一致形函数会整体错位——最典型的表现是刚度矩阵依然对称正定但计算结果整体偏移且误差随网格加密不收敛。我见过的最隐蔽案例是两个模块对“节点 4 在哪条边”的理解差了半个单元前面几百行代码全部正确只是接错了三个索引。5.2 用面积坐标统一生成形函数与导数六节点单元的代码我建议全部在面积坐标下写然后统一转换到物理坐标。这样形函数表达式最简洁且便于扩展到更高阶。def tri6_shape(L): L: (n,3) 面积坐标返回 (n,6) 形函数 L1, L2, L3 L[:, 0], L[:, 1], L[:, 2] return np.column_stack([ L1*(2*L1-1), L2*(2*L2-1), L3*(2*L3-1), 4*L1*L2, 4*L2*L3, 4*L3*L1 ]) def tri6_dshape_dL(L): 返回 (n,6,3)形函数对三个面积坐标的偏导 L1, L2, L3 L[:, 0], L[:, 1], L[:, 2] n L.shape[0] d np.zeros((n, 6, 3)) d[:, 0, 0] 4*L1 - 1 d[:, 1, 1] 4*L2 - 1 d[:, 2, 2] 4*L3 - 1 d[:, 3, 0] 4*L2; d[:, 3, 1] 4*L1 d[:, 4, 1] 4*L3; d[:, 4, 2] 4*L2 d[:, 5, 2] 4*L1; d[:, 5, 0] 4*L3 return d从面积坐标到物理坐标的转换靠的是∂Lᵢ/∂x和∂Lᵢ/∂y。这三个导数恰好是 3.3 节里线性单元形函数的导数因为面积坐标本身就是线性单元的形函数而且对六节点单元来说同样是常数∂L₁/∂x (y₂-y₃)/(2A)等等——这里要特别注意对于六节点单元如果边中点节点允许偏移三条边就是曲线2A不再是常数。标准的等参六节点三角形通过面积坐标的有理变换处理这种情况但如果不做等参映射把中点固定在几何边中点上那三条边就退化成直线此时∂Lᵢ/∂x才是常数。这两种实现方式的结果差别很大选择前要清楚自己在做什么。5.3 数值积分至少三点宁可七点三角形上的三点积分公式二次精度点1: (1/6, 1/6, 2/3) 权重 1/3 点2: (1/6, 2/3, 1/6) 权重 1/3 点3: (2/3, 1/6, 1/6) 权重 1/3这里的权重已经归一化积分时乘上单元面积即可GAUSS3 np.array([[1/6, 1/6, 2/3], [1/6, 2/3, 1/6], [2/3, 1/6, 1/6]]) W3 np.array([1/3, 1/3, 1/3]) def tri6_Ke(coords, D, thickness1.0): dLdx, dLdy, A tri3_derivatives(coords[:3]) # 六节点单元的几何边仍按顶点三角形处理 Ke np.zeros((12, 12)) for gp, w in zip(GAUSS3, W3): N tri6_shape(gp.reshape(1, 3))[0] dNdL tri6_dshape_dL(gp.reshape(1, 3))[0] # (6,3) dNdx dNdL dLdx dNdy dNdL dLdy B np.zeros((3, 12)) B[0, 0::2] dNdx B[1, 1::2] dNdy B[2, 0::2] dNdy B[2, 1::2] dNdx Ke w * A * thickness * (B.T D B) return Ke需要七点积分四次精度的场景不多主要是当你把材料非线性、几何非线性引入后被积函数次数上升三点就不够了。七点公式的坐标和权重在常用手册里都能查到我的建议是开发阶段先上七点跑通并把结果和三点对比确认差异在可接受范围后再降阶省时间。这个顺序比反过来安全得多因为降阶之后的误差非常隐蔽不会报错只会让结果悄悄偏软。5.4 分片检验验证形函数的最快手段比起一个个查系数更高效的办法是直接做分片检验。取一小片网格比如四个三角形拼成的方形区域施加与线性位移场一致的边界条件u(x,y) a₀ a₁x a₂y v(x,y) b₀ b₁x b₂y如果形函数满足线性完备性那么解出来的节点位移必须精确等于这个线性场在节点处的取值误差在机器精度量级。任何一个形函数写错、编号接错、积分不对这个检验都会立刻炸掉。def patch_test_shape(coords, n_rand200, seed0): 在单元内随机取点验证形函数能精确重现线性场 rng np.random.default_rng(seed) # 生成面积坐标样本 r rng.random((n_rand, 2)) mask r.sum(axis1) 1.0 r[mask] 1.0 - r[mask] # 反射回单位三角形 L np.column_stack([1 - r.sum(axis1), r]) # 取节点上的线性场值 a np.array([0.3, 1.7, -0.9]) u_node a[0] a[1]*coords[:, 0] a[2]*coords[:, 1] # 插值 N tri6_shape(L) u_interp N u_node u_exact a[0] a[1]*(N coords[:, 0]) a[2]*(N coords[:, 1]) return np.max(np.abs(u_interp - u_exact))如果这个返回值在 1e-14 量级形函数就是对的。这个检验不需要边界条件、不需要求解线性方程组几十行代码就能写完是我眼里性价比最高的验证手段。6. 划线法在实战中容易翻车的几个位置6.1 顶点顺序与雅可比符号第 1 章讲过的那次返工根子在这里。补充一点具体做法我的单元构造函数第一步永远是算行列式符号不对直接抛异常绝不让它悄悄往下走。有些人担心网格生成器偶尔会输出顺时针单元抛异常会导致大批量计算中断于是选择“自动纠正”。我的建议是不纠正直接报错——因为顺时针单元往往意味着网格生成环节有更系统的问题掩盖它只会让后面更难查。如果确实需要在后处理阶段容忍方向不一致比如读取外部网格可以在读入时统一做一次方向检查并把所有单元规整为逆时针但这一步必须是显式的、有日志记录的不能藏在Ke的计算里。6.2 畸形三角形让线性单元照样失准线性单元的应变是常数很多人因此认为它对形状不敏感。这个想法是错的。形函数梯度模长是1/hᵢ当某个顶点到对边的高很小时对应的梯度会剧烈放大B 矩阵元素之间量级差好几个数量级BᵀDB的条件数随之恶化求解出来的应力在这类单元里会出现明显振荡。判断标准不用太复杂看两个量就够最小内角低于 15° 就要警惕和纵横比最长边比最短高超过 10 就要处理。这两个量在网格生成后批量扫一遍成本很低。真正麻烦的是畸形单元往往出现在几何过渡区恰好是应力集中的地方结果错得最厉害的位置正好是你最关心的位置。6.3 高阶单元的边中点节点不要随手放在几何中点这是二次单元最容易被忽略的一点。如果边中点节点严格放在两个顶点的几何中点那么三条边就是直线单元退化成“直边二次单元”几何描述精度和线性单元一样——只有位移场是二次的。这在曲边界问题比如含孔的板、圆管截面上会导致边界离散误差成为主导误差无论你把位移场阶次提到多高都没用。正确的做法是用等参映射把边中点节点放在真实曲线的中点位置上通过面积坐标的有理变换把参考单元映射到曲边单元。这一步的代价是雅可比矩阵不再是常数每个积分点都要单独算代码复杂度上一个台阶。但如果你的问题里几何曲率是主要矛盾这个投入必须付。顺带一提等参映射会把∂Lᵢ/∂x变成非常数5.2 节里那句“对六节点单元来说同样是常数”就不再成立了。写代码时如果混用了两种假设会得到一个看起来能跑但精度不升反降的结果非常难查。6.4 我自己的一套固定自查流程把这几年踩过的坑浓缩成一套固定动作每次新增单元类型或者改动形函数代码我都会从头走一遍几何层算单元面积、最小高、最小角确认都是正数且在合理范围形函数层在单元内随机采 200 个点验证Σ Nᵢ 1、Nᵢ(xⱼ) δᵢⱼ、Σ ∇Nᵢ 0梯度层验证|∇Nᵢ| 1/hᵢ且方向垂直于对边线性单元专用高阶单元可以对每个顶点局部近似验证完备性层跑 5.4 节的分片检验误差应在机器精度刚度层验证Ke对称、对角元为正、特征值中恰好有 3 个零二维刚体模态六节点单元还要确认没有多余的零特征值有就说明积分阶次不够整体层做一个有解析解的小算例比如单向拉伸、纯弯曲比较收敛阶。第 5 条里的“恰好 3 个零特征值”是个特别锋利的判据。多一个零特征值就说明出现了沙漏模式少一个就说明约束加错了。用numpy.linalg.eigvalsh扫一遍几毫秒的事。最后分享一个我自己觉得挺有用的小习惯每推导出一种新单元的形函数我会先在纸上把等值线画出来——线性单元画三条平行线二次单元画那条中位线和高高低低的双曲线。画完再回头看公式如果图像和公式对不上一定是公式写错了。这个动作花不了两分钟但比我盯着矩阵看半天管用得多。形函数这东西说到底是一张画在三角形上的地形图公式只是这张图的坐标化描述先看见图再写公式顺序反过来就容易出事。