简介这份《计算力学——有限元编程实现》资源面向工科学生与编程初学者以C实现二维有限元分析涵盖几何建模、三角形三节点单元、四边形四节点等参元、八节点四边形等参元、刚度矩阵组装、节点与线性荷载处理、方程组求解及后处理等关键流程。资源共五十二个文件主要包含C源程序、工程配置文件、编译调试信息以及多套三角形和四边形网格数据文本压缩包大小为二十二点五一兆字节。目前已有1005人学习下载。借助完整的主程序与几何形状类实现可对三节点三角形、四节点及八节点四边形单元进行求解直接编译运行查看应力位移结果便于理解有限元编程的完整脉络适合课程设计、毕业设计或自学入门。1. 有限元编程实现把连续体方程变成可执行代码最值得先跑通的路线做结构分析的人大多有个共同经历软件用得越顺手越想亲手把刚度矩阵拼一遍。计算力学这条路上最硬核也最有长期回报的部分就是有限元编程实现。它的目标不是在键盘上复刻某个商业求解器而是用几百行代码把连续体问题离散成可求解的代数方程组走完从网格划分、单元刚度计算、整体组装到边界条件处理的全流程。本文适合刚完成弹性力学与数值分析学习、想亲手实现第一个平面应力求解器的开发者也适合天天点鼠标但始终觉得求解器像黑匣子的工程师。整条路线只用 Python 加开源稀疏线性代数库不依赖任何商业闭源组件跑通之后你对“单元”“节点”“刚度”的理解会和只看公式完全不同。2. 有限元求解的底层逻辑动手写代码前必须先定的 4 个关键决策在写第一行代码前先回答四个问题方程形式上用什么、位移怎么插值、积分怎么算、每个自由度代表什么。这四个答案直接决定后面所有数组的形状和循环方式也决定了你在程序报错时该怎么排查。2.1 为什么有限元不直接离散微分方程弱形式与虚功原理弹性静力学的起点是平衡方程 ∂σij/∂xj bi 0经典写法要求域内每一点都满足。编程时把这个强形式逐点做差分并不现实应力场不连续、边界形状不规则、载荷集中都会让差分格式难以收敛。有限元的做法是给方程乘一个任意的虚位移场 δu再对整个求解域积分得到虚功方程∫Ω σ:δϵ dΩ ∫Ω b·δu dΩ ∫Γ t·δu dΓ。这就是弱形式它把“每点平衡”弱化为“在一组加权积分意义下平衡”。实际编程中你根本不用关心这个方程的连续解长什么样只需要把它离散成一组线性代数方程。弱形式每次都对一个试函数空间积分这就是为什么刚度矩阵总是对称的——只要材料矩阵对称双线性形式自然对称。后面写组装的时候对称性可以帮助你只填上三角减少一半循环工作量。另一个容易被忽略的事实是弱形式天然包含边界条件自然边界条件已经在右端项里只有本质边界条件需要手动处理。很多第一次写代码的人把全部边界都当成“固定”去处理结果把自由边也约束住了算出来的位移场完全变形。2.2 形函数与等参映射局部坐标下的插值如何回到物理坐标Q4 四边形单元四节点双线性单元是最容易上手的单元。每个节点有 2 个位移自由度单元内部位移由四个节点位移插值得到。等参思想的关键是几何坐标也用同一组形函数插值。也就是说物理坐标 x、y 和位移 u、v 都由同一组 Ni(ξ,η) 表达。形函数在局部坐标里写死了四个节点分别是 (-1,-1), (1,-1), (1,1), (-1,1)对应插值 N1(1-ξ)(1-η)/4、N2(1ξ)(1-η)/4、N3(1ξ)(1η)/4、N4(1-ξ)(1η)/4。这带来的编程影响是你只需要在局部坐标 (ξ,η) 上写形函数形函数对局部坐标的导数固定不变再通过坐标变换矩阵 J 把导数和积分区域换算回物理坐标。换元后单元刚度矩阵写成 ke ∫_{-1}^1∫_{-1}^1 B^T D B t |J| dξdη其中 B 是应变-位移矩阵D 是材料弹性矩阵t 是厚度。J 在每个积分点上都要重新算因为等参单元的映射是随单元形状变化的。一个常见误区是以为 J 是常数只在规则矩形网格里 J 才恒定一旦网格畸变J 随位置变化偷懒把 J 提出积分号就会得到错误刚度。2.3 数值积分方案怎么选平面 Q4 的 2×2 规则单元刚度矩阵里的 B^T D B 通常是多项式解析积分存在但很繁琐。工程上一律用数值积分。Q4 单元形函数导数是 (ξ,η) 的一次函数B^T D B 的最高次为二次二维情形下用 2×2 积分点就可以精确积分。积分点取局部坐标 ±1/√3权重均为 1。这是有限元教材里的标准结论但真正写程序时会发现另一个细节积分点的取法直接影响单元刚度是否包含零能模式即所谓的沙漏模式。很多人习惯性加密到 3×3这对 Q4 单元属于浪费增大计算量但不提升收敛精度。反过来如果单元几何扭曲严重2×2 也够用真正要小心的是单元内 J 不能出现负值。对于八节点 Q8 单元B 矩阵含二次项通常需要 3×3 积分才能保证刚度矩阵不出现零能量模式。这个选型在后面讲扩展路线时绕不开。我在日常调试中还有一个习惯先写一个函数打印每个积分点的 detJ单元有问题时最先看到的就是这个数变成负数或零。2.4 单元刚度矩阵的物理意义每一行是给节点的力单刚中第 i 行第 j 列 Kij 的含义是第 j 个自由度发生单位位移时在第 i 个自由度上需要施加的力。对角线为正矩阵半正定。因为单元本身没有约束允许刚体位移单刚一定是奇异的。你看到单刚行列式为零不要紧张这是正常现象只有整体矩阵在施加足够约束后才应变为非奇异。调试时这一点非常关键。如果整体刚度矩阵依然奇异说明刚体模态没有被约束住最常见的错误是漏了某个方向的位移约束。把这个物理含义记在心里你就能理解为什么有限元代码的主干始终是同一个套路组装刚度、施加约束、求解、回代。不管换什么单元、什么维数自由度编号方式变化了但这个骨架不会变。3. 用 Python 实现 Q4 平面应力求解器从网格生成到约束求解的完整代码这一章直接给可复现代码。场景选为矩形平板单轴拉伸左端固定右端受均匀拉力材料线弹性平面应力假设。程序只依赖 NumPy 和 SciPy环境方面只需要保证能 import 这两个库。3.1 先把网格和单元表定义出来节点编号决定组装循环网格生成的第一要务是节点编号规则。常见做法是“先 x 后 y”x 方向从 0 到 nxy 方向从 0 到 ny节点编号为 n jy * (nx1) ix。这个编号顺序直接关系到单元表怎么写也会影响后面组装时带宽的大小。下面是生成矩形网格的最小代码import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla def make_mesh(nx, ny, Lx, Ly): 生成矩形域 Q4 网格返回节点坐标和单元表 nn (nx 1) * (ny 1) # 节点总数 nodes np.zeros((nn, 2)) for jy in range(ny 1): for ix in range(nx 1): n jy * (nx 1) ix nodes[n, 0] ix * Lx / nx nodes[n, 1] jy * Ly / ny elements [] for jy in range(ny): for ix in range(nx): n0 jy * (nx 1) ix # 逆时针左下、右下、右上、左上 ele [n0, n0 1, n0 nx 2, n0 nx 1] elements.append(ele) return nodes, np.array(elements)网格参数nx、ny决定单元数量Lx、Ly是物理尺寸。单元表按逆时针排列节点这是为了确保坐标变换矩阵的行列式为正。如果你按顺时针给节点单元的“面积”为负后面算出来的 detJ 会变成负数应力结果会彻底乱掉。节点编号的规律建议写在注释里因为组装自由度索引时要反复用到这个映射关系。3.2 单元刚度矩阵B、D、坐标变换与 2×2 积分的实现单元刚度是整段代码的核心。它做的事情是取出该单元四个节点的坐标构造坐标变换矩阵在每个积分点上计算 B 矩阵和 D 矩阵然后累加成 8×8 的单刚。平面应力问题的 D 矩阵里包含杨氏模量E和泊松比nudef element_stiffness(xy, E, nu, t): 平面应力 Q4 单元刚度矩阵2x2 积分 D E / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2] ]) gp [-1.0 / np.sqrt(3), 1.0 / np.sqrt(3)] wt [1.0, 1.0] ke np.zeros((8, 8)) for i, xi in enumerate(gp): for j, eta in enumerate(gp): # 形函数对局部坐标的导数4x2 dN np.array([ [-(1 - eta) / 4, (1 - eta) / 4, (1 eta) / 4, -(1 eta) / 4], [-(1 - xi) / 4, -(1 xi) / 4, (1 xi) / 4, (1 - xi) / 4] ]) J dN xy # 坐标变换矩阵 2x2 detJ np.linalg.det(J) if detJ 0: raise ValueError(detJ 0单元节点顺序或几何有误) Jinv np.linalg.inv(J) dNdx Jinv dN # 形函数对物理坐标的导数 2x4 B np.zeros((3, 8)) for a in range(4): B[0, 2*a] dNdx[0, a] B[1, 2*a1] dNdx[1, a] B[2, 2*a] dNdx[1, a] B[2, 2*a1] dNdx[0, a] ke wt[i] * wt[j] * B.T D B * detJ * t return ke参数t是板厚平面应力问题默认取 1 即可。np.linalg.det(J)检查负 Jacobian 这一步不能省它是最便宜的几何自检手段。dNdx Jinv dN是关键换元把形函数对局部坐标的导数变换到物理坐标。B 矩阵的组装规律是每个节点对应两列第一列放 ∂N/∂x 到 εxx 行和 ∂N/∂y 到 γxy 行第二列放 ∂N/∂y 到 εyy 行和 ∂N/∂x 到 γxy 行。这和三自由度剪应变定义直接相关写错一列就会得到不对称矩阵。3.3 稀疏组装与自由度映射COO 格式为什么适合有限元整体刚度矩阵的规模和自由度总数有关。比如 8×8 网格有 81 个节点、162 个自由度满矩阵是 162×162勉强能用但换成 100×100 网格满矩阵就是 20402×20402接近 3.3 GB 内存直接爆掉。有限元刚度矩阵是稀疏的每个单元的贡献只会落在一个 8×8 的小窗口里所以用 COO 格式收集三元组再转 CSR 是最自然的做法def assemble_K(nodes, elements, E, nu, t): nn len(nodes) ndof nn * 2 rows, cols, vals [], [], [] for ele in elements: xy nodes[ele] ke element_stiffness(xy, E, nu, t) # 该单元 4 个节点对应的 8 个自由度索引 dof np.empty(8, dtypeint) for a in range(4): dof[2*a] 2 * ele[a] dof[2*a1] 2 * ele[a] 1 for a in range(8): for b in range(8): rows.append(dof[a]) cols.append(dof[b]) vals.append(ke[a, b]) K sp.coo_matrix((vals, (rows, cols)), shape(ndof, ndof)).tocsr() return Kdof数组把局部自由度映射到全局自由度公式是“节点号 × 2 方向”。COO 格式允许同一个位置出现重复值tocsr()会把相同下标的值自动累加这正是整体组装需要的语义。初学者最容易在这里犯的错误是把单元自由度偏移写成了ele[a]*2和ele[a]*21但局部编号a与单元表顺序不一致。建议先打印一个单元的dof和单刚手工核对一下第一个非零元素落的位置。3.4 固定边界与等效节点荷载左端约束、右端受拉的完整求解约束施加直接采用“自由度消去法”而不是罚函数法。消去法没有罚系数需要调也不引入额外条件数适合教学代码。等效节点荷载方面右端均匀拉力在 Q4 单元边界上可以按“平均分配”处理因为常应力状态下分布力对每个节点的贡献是相等的。完整求解代码如下def solve_plane_stress(nx, ny, Lx, Ly, E, nu, t, sigma): nodes, elements make_mesh(nx, ny, Lx, Ly) ndof len(nodes) * 2 K assemble_K(nodes, elements, E, nu, t) # 左端固定x0 的所有节点两个方向都约束 fixed [] for n in range(len(nodes)): if abs(nodes[n, 0]) 1e-12: fixed [2*n, 2*n1] fixed np.array(sorted(set(fixed))) # 右端受拉xLx 的节点按节点数等分总拉力 right [n for n in range(len(nodes)) if abs(nodes[n, 0] - Lx) 1e-12] total_force sigma * Ly * t f np.zeros(ndof) for n in right: f[2*n] total_force / len(right) # 消去固定自由度求解自由部分 free np.setdiff1d(np.arange(ndof), fixed) Kff K[free][:, free].tocsr() u_free spla.spsolve(Kff, f[free]) # 还原完整位移向量 u np.zeros(ndof) u[free] u_free return nodes, elements, u, K, fsigma是施加在右端的拉应力值total_force等于应力乘以截面高度再乘以厚度单位必须和 E 一致。固定自由度只包含左端节点这样刚体位移被完全消除。要注意K[free][:, free]这种写法会对矩阵做两次索引复制中等规模没问题如果自由度超过 50 万更推荐在组装时通过自由度映射只组装自由部分否则内存会翻倍。3.5 从位移解里恢复应力一段 15 行的后处理函数位移解本身是近似值工程上更关心应力。应力的标准恢复方式是回到单元内部在积分点用 B 矩阵乘材料矩阵再乘位移。下面的函数在单元中心 (ξ0, η0) 处计算应力这个位置往往是单元内精度最好的点def recover_stress(nodes, elements, u, E, nu): D E / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2] ]) stresses [] for ele in elements: xy nodes[ele] # 单元中心处 dN 是常数 dN np.array([[-0.25, 0.25, 0.25, -0.25], [-0.25, -0.25, 0.25, 0.25]]) J dN xy Jinv np.linalg.inv(J) dNdx Jinv dN B np.zeros((3, 8)) for a in range(4): B[0, 2*a] dNdx[0, a] B[1, 2*a1] dNdx[1, a] B[2, 2*a] dNdx[1, a] B[2, 2*a1] dNdx[0, a] ue np.empty(8) for a in range(4): ue[2*a] u[2*ele[a]] ue[2*a1] u[2*ele[a]1] stress D B ue stresses.append(stress) return np.array(stresses)这段代码的输出是每个单元一个应力向量顺序为 [σxx, σyy, τxy]。ue的提取顺序必须和单元表局部编号一致否则 B 矩阵和位移对不上。用这段程序跑左侧固定、右侧受拉的例子会发现σxx 稳定在 sigma 附近σyy 和 τxy 接近零说明平面应力状态恢复正确。4. 结果怎么看应力后处理、网格收敛与力平衡验证代码能跑出位移只是第一步结果对不对、能不能写进报告是另一层问题。这一章讲三件事应力应该在哪个位置取、网格加密后结果怎么变、以及一个不用看云图就能判断对错的力平衡检查。4.1 节点应力和积分点应力的区别不要直接拿节点位移去画应力云图Q4 单元的位移场是双线性的应变由位移导数得到所以应变和应力在单元内是线性变化的。但节点处的位移本身是整体求解的近似值直接拿节点位移代入 B 矩阵求节点应力相当于在单元边界处用精度最低的导数做估计会造成应力在单元之间不连续云图出现明显锯齿。这种现象在应力梯度大的区域尤其刺眼。标准做法有两个。一是取积分点应力就是上一章的单元中心应力二是把相邻单元在同一节点处的应力做平均。后者做起来很简单维护一个“节点应力累加器”每算出一个单元应力就累加到该单元四个节点上最后除以节点被共享的次数。这个平均操作能让云图变得连续但严格讲它只是后处理平滑不会提高精度。对方案件里积分点应力才是“真值”节点平均应力只是“可视化值”两者在报告里最好分开说明。4.2 网格收敛性验证加密一倍位移和应力各怎么变有限元近似解的收敛性判断方法是逐步加密网格观察目标量是否稳定。对均匀拉伸这种应力场恒定的问题Q4 单元在 1×1 网格下就能给出精确位移加密不会带来变化这种情况说明问题太简单验证不出来。要测收敛性应该换成带应力集中的问题比如悬臂梁自由端受集中力或含圆孔平板受拉。悬臂梁是很好的验证基准左端固定右端自由端受竖向力 F理论挠度 δ FL³/(3EI)其中 I t·H³/12L 是梁长H 是截面高度。把网格从 4×1 加密到 8×2、16×4看自由端挠度与理论值之差的绝对值和相对误差。比较时需要取梁理论挠度公式的适用范围细长梁误差小短粗梁会包含剪切变形需要修正。这个对比能直接暴露单元太粗、边界约束过强或载荷等效出错的问题。我一般会记录三列网格规模、自由端位移、相对理论误差误差随加密稳步缩小才说明程序框架是健康的。4.3 力平衡检查支反力与外力之和为什么必须接近零位移云图可能看起来正常但数值可能整体偏了。最直接的自检是求解完成后计算 r f - K u这个向量在内部自由度上应该是数值零在固定自由度上就是支反力。整体支反力总和应该与外力总和平衡即 sum(r) 应当接近零。如果这个数大得离谱说明有节点力漏加或约束施加不对称。r f - K u print(整体力残差:, np.sum(r)) reaction np.zeros(len(nodes)) for i, n in enumerate(sorted(set(fixed // 2))): reaction[n] np.sqrt(r[2*n]**2 r[2*n1]**2)fixed // 2把自由度索引还原成节点号。力平衡检查不依赖任何理论解是程序自身的一致性条件只要刚度矩阵、载荷向量和约束处理中有一处不对称错误这个残差就会明显偏离零。建议把这段代码放在求解函数里每次跑完自动打印省掉大量肉眼排查。5. 有限元编程排雷手册新手最容易翻车的 5 个经典问题这一章是血泪经验的集合。每个问题按“现象 → 原因 → 解决”的顺序写对应一个典型的调试场景。5.1 求解器报错矩阵奇异问题不一定出在求解器现象spsolve抛错说矩阵 singular或者位移解巨大到 10^15。原因整体刚度矩阵缺少足够的约束刚体平移或转动模态没有被消除。解决先检查固定自由度列表。把固定节点打印出来确认 x 和 y 两个方向都约束了。一个常见的低级错误是只固定了节点号忘记了每个节点有两个自由度。至少要约束三个自由度且不能全在一个方向上才能消除平面问题的刚体位移。用下面的代码在求解前检查约束数量ndof len(nodes) * 2 print(自由度总数:, ndof, 固定自由度:, len(fixed)) rank_check sp.linalg.spsolve(K[free][:, free], np.zeros(len(free)))如果矩阵尺寸远小于自由度总数先检查free索引是否越界再看固定自由度的判断条件里浮点比较是否用了abs(x - 0) 1e-12而不是 0。5.2 单元体积为负、雅可比行列式为零节点编号顺序是第一嫌疑现象单刚计算抛detJ 0异常或者程序没报错但应力云图在某个单元附近出现不合理的正负交替。原因单元表节点顺序不是逆时针或者网格生成时节点编号映射错了导致坐标变换矩阵行列式为负。解决最直接的办法是用叉积检查每个单元的有向面积。对单元的前三个节点计算 (x1-x0)(y2-y0) - (x2-x0)(y1-y0)符号为正才是逆时针。检查脚本如下for k, ele in enumerate(elements): x0, y0 nodes[ele[0]] x1, y1 nodes[ele[1]] x2, y2 nodes[ele[2]] area2 (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0) if area2 0: print(单元, k, 定向错误)如果单元本身没问题那要看网格是否产生过度畸变比如某个单元被压成了接近零面积的细条。Q4 单元对畸变容忍度一般网格质量太差时即使 detJ 为正精度也会明显下降。5.3 计算结果的量级莫名其妙十有八九是单位制不统一现象算出来的应力是 10^11和理论值差了好几个量级挠度是 10^-8看起来像没加载。原因单位制混用。常见组合是把几何尺寸用毫米、材料弹性模量用 PaN/m²直接乘在一起或者把力的单位 kN 与长度的单位 m 混搭。解决程序内部只用一套单位制最常用的是“米-牛-帕”或“毫米-牛-兆帕”。输入参数前先列一个单位对照量方案 A米制方案 B毫米制长度mmm力NN弹性模量Pa (N/m²)MPa (N/mm²)应力PaMPa位移mmm检查方法很简单算一个简单拉伸或悬臂梁的小例子看位移量级是否合理。如果程序能复现理论解再批量跑其他算例。不要在代码里做单位换算函数那会让单位问题更难排查不如在入口处强制统一。5.4 应力云图出现锯齿把高斯点应力直接当成了节点应力现象云图颜色在单元边界处呈锯齿状同一节点在相邻单元里的应力值差一大截。原因后处理时直接用节点位移和单元角点处的形函数导数求应力等于在每个单元内独立外推应力不连续。解决改取单元中心或积分点应力再做节点平均。实现方式在第四章已经给过。如果坚持要节点应力至少先做相邻单元平均否则只看单元中心应力绘图。绘图时用matplotlib.tri的三角化可以把四边形拆成两个三角形云图会更干净但这属于可视化细节不影响物理量精度。5.5 等效节点力算错了分布载荷不能一律按节点均分现象右端受均匀拉力的例子里支反力总和与外力的差不为零换成线性变化的分布载荷后自由端位移明显偏小。原因把分布力简单除以边界节点数再分配只对常值拉力恰好成立载荷沿边界变化时每个节点应该承担的值并不相等。解决对分布载荷按照形函数做边界积分。Q4 边界上的形函数是一维线性函数比如边上两个节点 p、q边界线载荷为 q(s)贡献到 p 和 q 的节点力分别是 ∫ q·Np ds 和 ∫ q·Nq ds。对线性变化的载荷两个节点其实应该按 1:2 和 2:1 分配而不是各一半。最稳妥的做法是写一个通用的边界载荷积分循环沿每条载荷边做 2 点数值积分把结果累加到对应的自由度上。6. 把玩具求解器改造成可扩展版本先验证后提速的一次升级技巧代码跑通之后第一件事不是加功能而是确认“新版本结果和旧版本完全一致”。我一般保留一个 8×8 网格的旧输出作为基准任何性能优化都必须让最大位移差小于 10^-10才认为没有改坏逻辑。一个能显著提升组装速度的技巧是利用刚度矩阵的对称性。组装循环里单元刚度ke[a,b]和ke[b,a]相等因此只需要遍历上三角部分for a in range(8): for b in range(a, 8): v ke[a, b] rows.append(dof[a]); cols.append(dof[b]); vals.append(v) if a ! b: rows.append(dof[b]); cols.append(dof[a]); vals.append(v)这能把组装时间缩短近一半。配合上k_elements的数量级在十万以下时差别不大但单元数过万后体感明显。再做一层优化是把每个单元的dof索引预先算好存成数组避免每次组装都重新计算这一层改动不大但能减少大量 Python 整数运算。从 Q4 扩展到 Q8 或三维六面体单元时组装框架可以原样保留需要改的只有三处形函数和形函数导数、单元的节点数Q8 是 8 个节点 16 个自由度、积分点方案Q8 常用 3×3。材料非线性则需要把 D 矩阵换成随应力状态变化的切线矩阵并加上 Newton 迭代外循环那是下一步的工程话题。我最早自己写的时候也栽在约束漏了没发现后来养成了“先跑通、再提速、提速后立刻验证”的习惯希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/17 10:17:41 · 阅读 1.2k 从零实现有限元求解器：Q4平面应力分析与Python代码详解
简介这份《计算力学——有限元编程实现》资源面向工科学生与编程初学者以C实现二维有限元分析涵盖几何建模、三角形三节点单元、四边形四节点等参元、八节点四边形等参元、刚度矩阵组装、节点与线性荷载处理、方程组求解及后处理等关键流程。资源共五十二个文件主要包含C源程序、工程配置文件、编译调试信息以及多套三角形和四边形网格数据文本压缩包大小为二十二点五一兆字节。目前已有1005人学习下载。借助完整的主程序与几何形状类实现可对三节点三角形、四节点及八节点四边形单元进行求解直接编译运行查看应力位移结果便于理解有限元编程的完整脉络适合课程设计、毕业设计或自学入门。1. 有限元编程实现把连续体方程变成可执行代码最值得先跑通的路线做结构分析的人大多有个共同经历软件用得越顺手越想亲手把刚度矩阵拼一遍。计算力学这条路上最硬核也最有长期回报的部分就是有限元编程实现。它的目标不是在键盘上复刻某个商业求解器而是用几百行代码把连续体问题离散成可求解的代数方程组走完从网格划分、单元刚度计算、整体组装到边界条件处理的全流程。本文适合刚完成弹性力学与数值分析学习、想亲手实现第一个平面应力求解器的开发者也适合天天点鼠标但始终觉得求解器像黑匣子的工程师。整条路线只用 Python 加开源稀疏线性代数库不依赖任何商业闭源组件跑通之后你对“单元”“节点”“刚度”的理解会和只看公式完全不同。2. 有限元求解的底层逻辑动手写代码前必须先定的 4 个关键决策在写第一行代码前先回答四个问题方程形式上用什么、位移怎么插值、积分怎么算、每个自由度代表什么。这四个答案直接决定后面所有数组的形状和循环方式也决定了你在程序报错时该怎么排查。2.1 为什么有限元不直接离散微分方程弱形式与虚功原理弹性静力学的起点是平衡方程 ∂σij/∂xj bi 0经典写法要求域内每一点都满足。编程时把这个强形式逐点做差分并不现实应力场不连续、边界形状不规则、载荷集中都会让差分格式难以收敛。有限元的做法是给方程乘一个任意的虚位移场 δu再对整个求解域积分得到虚功方程∫Ω σ:δϵ dΩ ∫Ω b·δu dΩ ∫Γ t·δu dΓ。这就是弱形式它把“每点平衡”弱化为“在一组加权积分意义下平衡”。实际编程中你根本不用关心这个方程的连续解长什么样只需要把它离散成一组线性代数方程。弱形式每次都对一个试函数空间积分这就是为什么刚度矩阵总是对称的——只要材料矩阵对称双线性形式自然对称。后面写组装的时候对称性可以帮助你只填上三角减少一半循环工作量。另一个容易被忽略的事实是弱形式天然包含边界条件自然边界条件已经在右端项里只有本质边界条件需要手动处理。很多第一次写代码的人把全部边界都当成“固定”去处理结果把自由边也约束住了算出来的位移场完全变形。2.2 形函数与等参映射局部坐标下的插值如何回到物理坐标Q4 四边形单元四节点双线性单元是最容易上手的单元。每个节点有 2 个位移自由度单元内部位移由四个节点位移插值得到。等参思想的关键是几何坐标也用同一组形函数插值。也就是说物理坐标 x、y 和位移 u、v 都由同一组 Ni(ξ,η) 表达。形函数在局部坐标里写死了四个节点分别是 (-1,-1), (1,-1), (1,1), (-1,1)对应插值 N1(1-ξ)(1-η)/4、N2(1ξ)(1-η)/4、N3(1ξ)(1η)/4、N4(1-ξ)(1η)/4。这带来的编程影响是你只需要在局部坐标 (ξ,η) 上写形函数形函数对局部坐标的导数固定不变再通过坐标变换矩阵 J 把导数和积分区域换算回物理坐标。换元后单元刚度矩阵写成 ke ∫_{-1}^1∫_{-1}^1 B^T D B t |J| dξdη其中 B 是应变-位移矩阵D 是材料弹性矩阵t 是厚度。J 在每个积分点上都要重新算因为等参单元的映射是随单元形状变化的。一个常见误区是以为 J 是常数只在规则矩形网格里 J 才恒定一旦网格畸变J 随位置变化偷懒把 J 提出积分号就会得到错误刚度。2.3 数值积分方案怎么选平面 Q4 的 2×2 规则单元刚度矩阵里的 B^T D B 通常是多项式解析积分存在但很繁琐。工程上一律用数值积分。Q4 单元形函数导数是 (ξ,η) 的一次函数B^T D B 的最高次为二次二维情形下用 2×2 积分点就可以精确积分。积分点取局部坐标 ±1/√3权重均为 1。这是有限元教材里的标准结论但真正写程序时会发现另一个细节积分点的取法直接影响单元刚度是否包含零能模式即所谓的沙漏模式。很多人习惯性加密到 3×3这对 Q4 单元属于浪费增大计算量但不提升收敛精度。反过来如果单元几何扭曲严重2×2 也够用真正要小心的是单元内 J 不能出现负值。对于八节点 Q8 单元B 矩阵含二次项通常需要 3×3 积分才能保证刚度矩阵不出现零能量模式。这个选型在后面讲扩展路线时绕不开。我在日常调试中还有一个习惯先写一个函数打印每个积分点的 detJ单元有问题时最先看到的就是这个数变成负数或零。2.4 单元刚度矩阵的物理意义每一行是给节点的力单刚中第 i 行第 j 列 Kij 的含义是第 j 个自由度发生单位位移时在第 i 个自由度上需要施加的力。对角线为正矩阵半正定。因为单元本身没有约束允许刚体位移单刚一定是奇异的。你看到单刚行列式为零不要紧张这是正常现象只有整体矩阵在施加足够约束后才应变为非奇异。调试时这一点非常关键。如果整体刚度矩阵依然奇异说明刚体模态没有被约束住最常见的错误是漏了某个方向的位移约束。把这个物理含义记在心里你就能理解为什么有限元代码的主干始终是同一个套路组装刚度、施加约束、求解、回代。不管换什么单元、什么维数自由度编号方式变化了但这个骨架不会变。3. 用 Python 实现 Q4 平面应力求解器从网格生成到约束求解的完整代码这一章直接给可复现代码。场景选为矩形平板单轴拉伸左端固定右端受均匀拉力材料线弹性平面应力假设。程序只依赖 NumPy 和 SciPy环境方面只需要保证能 import 这两个库。3.1 先把网格和单元表定义出来节点编号决定组装循环网格生成的第一要务是节点编号规则。常见做法是“先 x 后 y”x 方向从 0 到 nxy 方向从 0 到 ny节点编号为 n jy * (nx1) ix。这个编号顺序直接关系到单元表怎么写也会影响后面组装时带宽的大小。下面是生成矩形网格的最小代码import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla def make_mesh(nx, ny, Lx, Ly): 生成矩形域 Q4 网格返回节点坐标和单元表 nn (nx 1) * (ny 1) # 节点总数 nodes np.zeros((nn, 2)) for jy in range(ny 1): for ix in range(nx 1): n jy * (nx 1) ix nodes[n, 0] ix * Lx / nx nodes[n, 1] jy * Ly / ny elements [] for jy in range(ny): for ix in range(nx): n0 jy * (nx 1) ix # 逆时针左下、右下、右上、左上 ele [n0, n0 1, n0 nx 2, n0 nx 1] elements.append(ele) return nodes, np.array(elements)网格参数nx、ny决定单元数量Lx、Ly是物理尺寸。单元表按逆时针排列节点这是为了确保坐标变换矩阵的行列式为正。如果你按顺时针给节点单元的“面积”为负后面算出来的 detJ 会变成负数应力结果会彻底乱掉。节点编号的规律建议写在注释里因为组装自由度索引时要反复用到这个映射关系。3.2 单元刚度矩阵B、D、坐标变换与 2×2 积分的实现单元刚度是整段代码的核心。它做的事情是取出该单元四个节点的坐标构造坐标变换矩阵在每个积分点上计算 B 矩阵和 D 矩阵然后累加成 8×8 的单刚。平面应力问题的 D 矩阵里包含杨氏模量E和泊松比nudef element_stiffness(xy, E, nu, t): 平面应力 Q4 单元刚度矩阵2x2 积分 D E / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2] ]) gp [-1.0 / np.sqrt(3), 1.0 / np.sqrt(3)] wt [1.0, 1.0] ke np.zeros((8, 8)) for i, xi in enumerate(gp): for j, eta in enumerate(gp): # 形函数对局部坐标的导数4x2 dN np.array([ [-(1 - eta) / 4, (1 - eta) / 4, (1 eta) / 4, -(1 eta) / 4], [-(1 - xi) / 4, -(1 xi) / 4, (1 xi) / 4, (1 - xi) / 4] ]) J dN xy # 坐标变换矩阵 2x2 detJ np.linalg.det(J) if detJ 0: raise ValueError(detJ 0单元节点顺序或几何有误) Jinv np.linalg.inv(J) dNdx Jinv dN # 形函数对物理坐标的导数 2x4 B np.zeros((3, 8)) for a in range(4): B[0, 2*a] dNdx[0, a] B[1, 2*a1] dNdx[1, a] B[2, 2*a] dNdx[1, a] B[2, 2*a1] dNdx[0, a] ke wt[i] * wt[j] * B.T D B * detJ * t return ke参数t是板厚平面应力问题默认取 1 即可。np.linalg.det(J)检查负 Jacobian 这一步不能省它是最便宜的几何自检手段。dNdx Jinv dN是关键换元把形函数对局部坐标的导数变换到物理坐标。B 矩阵的组装规律是每个节点对应两列第一列放 ∂N/∂x 到 εxx 行和 ∂N/∂y 到 γxy 行第二列放 ∂N/∂y 到 εyy 行和 ∂N/∂x 到 γxy 行。这和三自由度剪应变定义直接相关写错一列就会得到不对称矩阵。3.3 稀疏组装与自由度映射COO 格式为什么适合有限元整体刚度矩阵的规模和自由度总数有关。比如 8×8 网格有 81 个节点、162 个自由度满矩阵是 162×162勉强能用但换成 100×100 网格满矩阵就是 20402×20402接近 3.3 GB 内存直接爆掉。有限元刚度矩阵是稀疏的每个单元的贡献只会落在一个 8×8 的小窗口里所以用 COO 格式收集三元组再转 CSR 是最自然的做法def assemble_K(nodes, elements, E, nu, t): nn len(nodes) ndof nn * 2 rows, cols, vals [], [], [] for ele in elements: xy nodes[ele] ke element_stiffness(xy, E, nu, t) # 该单元 4 个节点对应的 8 个自由度索引 dof np.empty(8, dtypeint) for a in range(4): dof[2*a] 2 * ele[a] dof[2*a1] 2 * ele[a] 1 for a in range(8): for b in range(8): rows.append(dof[a]) cols.append(dof[b]) vals.append(ke[a, b]) K sp.coo_matrix((vals, (rows, cols)), shape(ndof, ndof)).tocsr() return Kdof数组把局部自由度映射到全局自由度公式是“节点号 × 2 方向”。COO 格式允许同一个位置出现重复值tocsr()会把相同下标的值自动累加这正是整体组装需要的语义。初学者最容易在这里犯的错误是把单元自由度偏移写成了ele[a]*2和ele[a]*21但局部编号a与单元表顺序不一致。建议先打印一个单元的dof和单刚手工核对一下第一个非零元素落的位置。3.4 固定边界与等效节点荷载左端约束、右端受拉的完整求解约束施加直接采用“自由度消去法”而不是罚函数法。消去法没有罚系数需要调也不引入额外条件数适合教学代码。等效节点荷载方面右端均匀拉力在 Q4 单元边界上可以按“平均分配”处理因为常应力状态下分布力对每个节点的贡献是相等的。完整求解代码如下def solve_plane_stress(nx, ny, Lx, Ly, E, nu, t, sigma): nodes, elements make_mesh(nx, ny, Lx, Ly) ndof len(nodes) * 2 K assemble_K(nodes, elements, E, nu, t) # 左端固定x0 的所有节点两个方向都约束 fixed [] for n in range(len(nodes)): if abs(nodes[n, 0]) 1e-12: fixed [2*n, 2*n1] fixed np.array(sorted(set(fixed))) # 右端受拉xLx 的节点按节点数等分总拉力 right [n for n in range(len(nodes)) if abs(nodes[n, 0] - Lx) 1e-12] total_force sigma * Ly * t f np.zeros(ndof) for n in right: f[2*n] total_force / len(right) # 消去固定自由度求解自由部分 free np.setdiff1d(np.arange(ndof), fixed) Kff K[free][:, free].tocsr() u_free spla.spsolve(Kff, f[free]) # 还原完整位移向量 u np.zeros(ndof) u[free] u_free return nodes, elements, u, K, fsigma是施加在右端的拉应力值total_force等于应力乘以截面高度再乘以厚度单位必须和 E 一致。固定自由度只包含左端节点这样刚体位移被完全消除。要注意K[free][:, free]这种写法会对矩阵做两次索引复制中等规模没问题如果自由度超过 50 万更推荐在组装时通过自由度映射只组装自由部分否则内存会翻倍。3.5 从位移解里恢复应力一段 15 行的后处理函数位移解本身是近似值工程上更关心应力。应力的标准恢复方式是回到单元内部在积分点用 B 矩阵乘材料矩阵再乘位移。下面的函数在单元中心 (ξ0, η0) 处计算应力这个位置往往是单元内精度最好的点def recover_stress(nodes, elements, u, E, nu): D E / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2] ]) stresses [] for ele in elements: xy nodes[ele] # 单元中心处 dN 是常数 dN np.array([[-0.25, 0.25, 0.25, -0.25], [-0.25, -0.25, 0.25, 0.25]]) J dN xy Jinv np.linalg.inv(J) dNdx Jinv dN B np.zeros((3, 8)) for a in range(4): B[0, 2*a] dNdx[0, a] B[1, 2*a1] dNdx[1, a] B[2, 2*a] dNdx[1, a] B[2, 2*a1] dNdx[0, a] ue np.empty(8) for a in range(4): ue[2*a] u[2*ele[a]] ue[2*a1] u[2*ele[a]1] stress D B ue stresses.append(stress) return np.array(stresses)这段代码的输出是每个单元一个应力向量顺序为 [σxx, σyy, τxy]。ue的提取顺序必须和单元表局部编号一致否则 B 矩阵和位移对不上。用这段程序跑左侧固定、右侧受拉的例子会发现σxx 稳定在 sigma 附近σyy 和 τxy 接近零说明平面应力状态恢复正确。4. 结果怎么看应力后处理、网格收敛与力平衡验证代码能跑出位移只是第一步结果对不对、能不能写进报告是另一层问题。这一章讲三件事应力应该在哪个位置取、网格加密后结果怎么变、以及一个不用看云图就能判断对错的力平衡检查。4.1 节点应力和积分点应力的区别不要直接拿节点位移去画应力云图Q4 单元的位移场是双线性的应变由位移导数得到所以应变和应力在单元内是线性变化的。但节点处的位移本身是整体求解的近似值直接拿节点位移代入 B 矩阵求节点应力相当于在单元边界处用精度最低的导数做估计会造成应力在单元之间不连续云图出现明显锯齿。这种现象在应力梯度大的区域尤其刺眼。标准做法有两个。一是取积分点应力就是上一章的单元中心应力二是把相邻单元在同一节点处的应力做平均。后者做起来很简单维护一个“节点应力累加器”每算出一个单元应力就累加到该单元四个节点上最后除以节点被共享的次数。这个平均操作能让云图变得连续但严格讲它只是后处理平滑不会提高精度。对方案件里积分点应力才是“真值”节点平均应力只是“可视化值”两者在报告里最好分开说明。4.2 网格收敛性验证加密一倍位移和应力各怎么变有限元近似解的收敛性判断方法是逐步加密网格观察目标量是否稳定。对均匀拉伸这种应力场恒定的问题Q4 单元在 1×1 网格下就能给出精确位移加密不会带来变化这种情况说明问题太简单验证不出来。要测收敛性应该换成带应力集中的问题比如悬臂梁自由端受集中力或含圆孔平板受拉。悬臂梁是很好的验证基准左端固定右端自由端受竖向力 F理论挠度 δ FL³/(3EI)其中 I t·H³/12L 是梁长H 是截面高度。把网格从 4×1 加密到 8×2、16×4看自由端挠度与理论值之差的绝对值和相对误差。比较时需要取梁理论挠度公式的适用范围细长梁误差小短粗梁会包含剪切变形需要修正。这个对比能直接暴露单元太粗、边界约束过强或载荷等效出错的问题。我一般会记录三列网格规模、自由端位移、相对理论误差误差随加密稳步缩小才说明程序框架是健康的。4.3 力平衡检查支反力与外力之和为什么必须接近零位移云图可能看起来正常但数值可能整体偏了。最直接的自检是求解完成后计算 r f - K u这个向量在内部自由度上应该是数值零在固定自由度上就是支反力。整体支反力总和应该与外力总和平衡即 sum(r) 应当接近零。如果这个数大得离谱说明有节点力漏加或约束施加不对称。r f - K u print(整体力残差:, np.sum(r)) reaction np.zeros(len(nodes)) for i, n in enumerate(sorted(set(fixed // 2))): reaction[n] np.sqrt(r[2*n]**2 r[2*n1]**2)fixed // 2把自由度索引还原成节点号。力平衡检查不依赖任何理论解是程序自身的一致性条件只要刚度矩阵、载荷向量和约束处理中有一处不对称错误这个残差就会明显偏离零。建议把这段代码放在求解函数里每次跑完自动打印省掉大量肉眼排查。5. 有限元编程排雷手册新手最容易翻车的 5 个经典问题这一章是血泪经验的集合。每个问题按“现象 → 原因 → 解决”的顺序写对应一个典型的调试场景。5.1 求解器报错矩阵奇异问题不一定出在求解器现象spsolve抛错说矩阵 singular或者位移解巨大到 10^15。原因整体刚度矩阵缺少足够的约束刚体平移或转动模态没有被消除。解决先检查固定自由度列表。把固定节点打印出来确认 x 和 y 两个方向都约束了。一个常见的低级错误是只固定了节点号忘记了每个节点有两个自由度。至少要约束三个自由度且不能全在一个方向上才能消除平面问题的刚体位移。用下面的代码在求解前检查约束数量ndof len(nodes) * 2 print(自由度总数:, ndof, 固定自由度:, len(fixed)) rank_check sp.linalg.spsolve(K[free][:, free], np.zeros(len(free)))如果矩阵尺寸远小于自由度总数先检查free索引是否越界再看固定自由度的判断条件里浮点比较是否用了abs(x - 0) 1e-12而不是 0。5.2 单元体积为负、雅可比行列式为零节点编号顺序是第一嫌疑现象单刚计算抛detJ 0异常或者程序没报错但应力云图在某个单元附近出现不合理的正负交替。原因单元表节点顺序不是逆时针或者网格生成时节点编号映射错了导致坐标变换矩阵行列式为负。解决最直接的办法是用叉积检查每个单元的有向面积。对单元的前三个节点计算 (x1-x0)(y2-y0) - (x2-x0)(y1-y0)符号为正才是逆时针。检查脚本如下for k, ele in enumerate(elements): x0, y0 nodes[ele[0]] x1, y1 nodes[ele[1]] x2, y2 nodes[ele[2]] area2 (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0) if area2 0: print(单元, k, 定向错误)如果单元本身没问题那要看网格是否产生过度畸变比如某个单元被压成了接近零面积的细条。Q4 单元对畸变容忍度一般网格质量太差时即使 detJ 为正精度也会明显下降。5.3 计算结果的量级莫名其妙十有八九是单位制不统一现象算出来的应力是 10^11和理论值差了好几个量级挠度是 10^-8看起来像没加载。原因单位制混用。常见组合是把几何尺寸用毫米、材料弹性模量用 PaN/m²直接乘在一起或者把力的单位 kN 与长度的单位 m 混搭。解决程序内部只用一套单位制最常用的是“米-牛-帕”或“毫米-牛-兆帕”。输入参数前先列一个单位对照量方案 A米制方案 B毫米制长度mmm力NN弹性模量Pa (N/m²)MPa (N/mm²)应力PaMPa位移mmm检查方法很简单算一个简单拉伸或悬臂梁的小例子看位移量级是否合理。如果程序能复现理论解再批量跑其他算例。不要在代码里做单位换算函数那会让单位问题更难排查不如在入口处强制统一。5.4 应力云图出现锯齿把高斯点应力直接当成了节点应力现象云图颜色在单元边界处呈锯齿状同一节点在相邻单元里的应力值差一大截。原因后处理时直接用节点位移和单元角点处的形函数导数求应力等于在每个单元内独立外推应力不连续。解决改取单元中心或积分点应力再做节点平均。实现方式在第四章已经给过。如果坚持要节点应力至少先做相邻单元平均否则只看单元中心应力绘图。绘图时用matplotlib.tri的三角化可以把四边形拆成两个三角形云图会更干净但这属于可视化细节不影响物理量精度。5.5 等效节点力算错了分布载荷不能一律按节点均分现象右端受均匀拉力的例子里支反力总和与外力的差不为零换成线性变化的分布载荷后自由端位移明显偏小。原因把分布力简单除以边界节点数再分配只对常值拉力恰好成立载荷沿边界变化时每个节点应该承担的值并不相等。解决对分布载荷按照形函数做边界积分。Q4 边界上的形函数是一维线性函数比如边上两个节点 p、q边界线载荷为 q(s)贡献到 p 和 q 的节点力分别是 ∫ q·Np ds 和 ∫ q·Nq ds。对线性变化的载荷两个节点其实应该按 1:2 和 2:1 分配而不是各一半。最稳妥的做法是写一个通用的边界载荷积分循环沿每条载荷边做 2 点数值积分把结果累加到对应的自由度上。6. 把玩具求解器改造成可扩展版本先验证后提速的一次升级技巧代码跑通之后第一件事不是加功能而是确认“新版本结果和旧版本完全一致”。我一般保留一个 8×8 网格的旧输出作为基准任何性能优化都必须让最大位移差小于 10^-10才认为没有改坏逻辑。一个能显著提升组装速度的技巧是利用刚度矩阵的对称性。组装循环里单元刚度ke[a,b]和ke[b,a]相等因此只需要遍历上三角部分for a in range(8): for b in range(a, 8): v ke[a, b] rows.append(dof[a]); cols.append(dof[b]); vals.append(v) if a ! b: rows.append(dof[b]); cols.append(dof[a]); vals.append(v)这能把组装时间缩短近一半。配合上k_elements的数量级在十万以下时差别不大但单元数过万后体感明显。再做一层优化是把每个单元的dof索引预先算好存成数组避免每次组装都重新计算这一层改动不大但能减少大量 Python 整数运算。从 Q4 扩展到 Q8 或三维六面体单元时组装框架可以原样保留需要改的只有三处形函数和形函数导数、单元的节点数Q8 是 8 个节点 16 个自由度、积分点方案Q8 常用 3×3。材料非线性则需要把 D 矩阵换成随应力状态变化的切线矩阵并加上 Newton 迭代外循环那是下一步的工程话题。我最早自己写的时候也栽在约束漏了没发现后来养成了“先跑通、再提速、提速后立刻验证”的习惯希望帮到你。本文还有配套的精品资源点击获取