前阵子我写一个简单的平面刚架求解器自编了一个梁单元算悬臂梁自由端挠度。手算答案是 PL³/(3EI)程序却差了一截。我以为是积分点不够换成更高阶的高斯积分还是不对。查到最后问题出在单元刚度矩阵的推导上——我随手套了一个网上流传的“标准公式”却没搞明白它是从虚功原理推出来的公式里的符号方向和我代码里的节点自由度定义根本不一致。从那次之后我意识到做结构力学程序算法第一课不是编程技巧而是把虚功原理吃透。虚功原理这个词搞结构的人都不陌生但很多人对它的理解停留在“考试会做求位移的题”这个层面。一旦要自己写算法、自己推导单元刚度矩阵、自己处理边界条件就会发现虚功原理其实贯穿了所有东西平衡方程、几何方程、物理方程、位移法、力法、有限元。它不是“一个公式”而是一套逻辑框架。这篇文章是结构力学程序算法理论基础系列的第一篇我打算用做程序的角度把虚功原理彻底拆开讲清楚重点讲它在算法落地的三个场景单元刚度矩阵怎么推、位移怎么用单位荷载法算、程序里容易踩哪些概念坑。这篇内容适合三类人一类是准备写结构分析程序的学生或工程师一类是天天用有限元软件但想知道黑箱里发生了什么的人还有一类是考注册结构、想把虚功原理真正弄懂而不是死记公式的人。读完你至少能自己推导一个杆系单元的刚度矩阵能说清楚单位荷载法为什么可行也知道为什么程序报错“刚度矩阵奇异”时往往是虚位移空间出了问题。1. 为什么结构程序绕不开虚功原理1.1 我自己的一个翻车现场程序输出与手算不符先说前面那个悬臂梁的事。我用的是经典欧拉-伯努利梁单元理论上刚度矩阵是k EI/L³ × [[12, 6L, -12, 6L], [6L, 4L², -6L, 2L²], [-12, -6L, 12, -6L], [6L, 2L², -6L, 4L²]]这个矩阵我背得滚瓜烂熟但那次写程序时用了不同的节点自由度顺序——先所有平动自由度再所有转动自由度。矩阵里那一堆 6L、12 之类的项对应关系完全错位了。程序没报错结果就是不对。这让我意识到一件事如果你只背公式不掌握公式是怎么从虚功原理推出来的你连“检查公式是否适用于自己的自由度排列”这种最基本的事都做不到。后来我老老实实从虚功方程把刚度矩阵推了一遍所有疑惑都消失了。你想一个连单元矩阵都需要验证的程序员虚功原理是绕不过去的。1.2 虚功原理在程序算法中的地位不只是推导工具在写结构程序的时候我们面对的核心问题是已知结构几何、材料、荷载求节点位移和杆端力。这个问题用矩阵位移法来解流程大致是划分单元给节点编号。对每个单元建立刚度矩阵。把单元刚度矩阵组装成整体刚度矩阵。引入边界条件。求解线性方程组得到节点位移。回代得到杆端力。这个流程里第 2 步是关键中的关键。那单元刚度矩阵从哪里来虚功原理。更准确地说是从虚位移原理推导出来的。不止是杆系单元平面应力单元、板壳单元、实体单元的刚度矩阵清一色都是这么来的。此外还有一类场景求结构某一点的位移或转角。手算时我们用单位荷载法而单位荷载法的理论基础就是虚功原理中的虚力原理。写程序时如果你想输出“某节点在某种荷载组合下的位移”完全可以通过虚功原理做后处理校验而不必重新求解整体方程。所以虚功原理在程序算法里承担的是“母体”角色一切公式从它而出一切结果可以用它来验证。你把它当成一个必须背的公式那就太亏了。1.3 一句话讲清楚虚功是什么很多人被“虚”这个字绊住。虚功的“虚”不是说这个功是假的、不存在的而是说做功的位移是“虚设的、几何可能的、不一定真实发生的位移”。我常用一个类比你想测一个弹簧的刚度正常做法是压它一下量出力和位移的比值。但“虚”的做法是假设这个弹簧被压缩了一个微小量 δ这个 δ 可以不是你真正施加的位移然后计算力在这个假设位移上做的功。如果这个弹簧本身处于平衡状态那么这个“虚功”就一定满足某种关系。这里的关键点有三个虚位移必须满足结构的约束条件比如固定端处虚位移为零。虚位移必须是微小的这样才能用一阶变分处理。虚位移是人为假想的与真实加载过程无关。程序算法处理的是离散化后的节点自由度所以虚位移在程序里就是一个虚拟的节点位移向量 δU。后面你会发现所有平衡方程都可以写成“把 δU 乘到方程两边利用 δU 的任意性”这样的形式。虚功原理就是把微分形式的平衡方程转化成积分或求和形式的变分方程。这个转化正是数值方法能够落地的基础。2. 虚位移原理与虚力原理两个方向两种用途2.1 虚位移原理从位移场出发找平衡条件虚位移原理说的是如果一个结构处于平衡状态那么所有外力包括支座反力在任意满足约束的虚位移上所做的虚功等于结构内力在相应虚应变上所做的虚功。用数学式子表达就是δW_ext δW_int其中δW_ext Σ Pᵢ × δuᵢ ∫ q × δv dxδW_int ∫ σ × δε dV这个式子的妙处在于它把“平衡”从“一个点的力为零”扩展成了“整个结构的能量关系”。在程序里我们不用去管每个微元体的受力平衡只需要保证节点力与内力在这个积分关系下平衡就够了。这就是有限元方法里“弱形式”的由来——不强求每一点都严格满足平衡微分方程而只要求在积分意义下满足。如果你在悬臂梁自由端加一个向下的集中力 P真实位移场是那个三次曲线。现在你假想一个虚位移 δv(x)它满足固定端位移为零。虚位移原理就会告诉你P × δv(自由端) ∫ M × δκ dx其中 δκ 是虚曲率。正是这个方程最终导出了梁单元刚度矩阵。2.2 虚力原理从力场出发找变形协调虚力原理是另一个方向。它说的是如果一个结构的变形是协调的就是几何上连续、没有裂缝或重叠那么虚设的平衡力系在外力作用点上的虚余功等于虚设内力在真实变形上做的虚余功。它的用途和虚位移原理正好相反。虚位移原理用来求力或者建立刚度方程虚力原理用来求位移。手算结构力学里面的单位荷载法用的就是虚力原理。你可以这样理解这两个原理的分工虚位移原理已知可能的位移虚位移求力是否平衡 → 对应程序里的刚度法位移法。虚力原理已知可能的力虚力求位移是否协调 → 对应程序里的柔度法力法。写程序的人用到虚力原理最多的地方就是求某点位移在需求位移的地方加一个单位虚力用这个虚力产生的内力场去和真实荷载产生的内力场做积分就得到真实位移。这个方法在结构力学里叫单位荷载法后面我会专门讲它在程序里怎么实现。2.3 两种原理在程序中的分工我整理了一个对照表写程序前看这个表思路会很清晰比较项虚位移原理虚力原理虚设的量虚位移满足约束虚力满足平衡已知量真实外力、真实内力真实位移/应变求解目标平衡条件、刚度方程位移、变形协调程序应用推导单元刚度矩阵单位荷载法求位移对应方法位移法、刚度法、有限元力法、柔度法这两种原理并不矛盾它们其实是同一个能量互等定理的两个侧面。程序里先用力法求出位移影响系数或者用单位荷载法校验位移结果都是常见操作。理解了这一点你就不容易把虚功原理和“只能用来求位移”划等号了。3. 从连续体到离散化虚功方程如何变成矩阵运算3.1 插值形函数连续位移场的“降维打击”写程序时我们不能把梁上每一个点的位移都当成未知数那样未知量是无穷多个没法算。所以我们引入“形函数”这个概念把单元内任意一点的位移用节点位移的插值来表示。比如一个两节点梁单元每个节点有挠度和转角两个自由度单元内任意一点的挠度可以写成v(x) N₁(x)v₁ N₂(x)θ₁ N₃(x)v₂ N₄(x)θ₂其中 Nᵢ(x) 就是形函数v₁、θ₁、v₂、θ₂ 是节点位移。对欧拉-伯努利梁这组形函数是三次 Hermite 插值它保证了挠度连续、转角连续。这一步看起来只是数学处理但它和虚功原理结合之后就非常关键了虚位移同样可以用同一组形函数来插值。也就是说δv(x) N₁(x)δv₁ N₂(x)δθ₁ N₃(x)δv₂ N₄(x)δθ₂为什么必须是同一组形函数因为只有这样才能让虚功方程两边统一到节点自由度上。这就是伽辽金法Galerkin法的核心思想真实位移和虚位移在同一个函数空间里。3.2 单元刚度矩阵的虚功推导过程现在我们把虚位移原理应用到离散单元上。外力虚功可以写成节点力向量 F 与虚节点位移 δU 的点积δW_ext δUᵀ F内力虚功呢需要先把虚应变和虚位移联系起来。对梁单元曲率 κ 与挠度的二阶导数有关δκ d²(δv)/dx²用形函数表示就是 δκ B × δU这个 B 矩阵是形函数二阶导数的组合。于是内力虚功为δW_int ∫ δκᵀ × EI × κ dx ∫ (B δU)ᵀ × EI × (B U) dx δUᵀ × (∫ Bᵀ EI B dx) × U根据虚位移原理δW_ext δW_int所以δUᵀ F δUᵀ × (∫ Bᵀ EI B dx) × U因为 δU 是任意的所以括号里的积分式就等于刚度矩阵k ∫ Bᵀ EI B dx没错这个积分就是单元刚度矩阵的来源。我之前翻车就是因为背了 k 的最终形式却没意识到它是形函数二阶导、弹性模量和截面惯性矩在单元长度上的积分。一旦自由度顺序变了B 矩阵里各个列的顺序也要跟着变光背公式是背不住的。类似的推导对其他单元也成立。平面应力单元的刚度矩阵是 k ∫ Bᵀ D B t dA其中 D 是弹性矩阵t 是厚度。实体单元则是 k ∫ Bᵀ D B dV。所有位移类单元刚度矩阵都是同一个套路Bᵀ D B 的积分。3.3 为什么高斯积分不是可有可无的上文出现了积分 ∫ Bᵀ EI B dx。程序里这个积分怎么算对等截面直杆Bᵀ EI B 被积函数是多项式理论上可以精确积分但程序里我们一般统一用数值积分——最常见的是高斯积分。高斯积分的思路是在单元内取若干积分点计算每个积分点上的被积函数值加权求和。对两节点梁单元因为形函数是三次的B 矩阵是形函数的二阶导所以 Bᵀ EI B 最多是二次多项式理论上取 2 个高斯积分点即可精确积分。但不要因为这个例子简单就轻视积分点的选择。在平面单元里如果积分点不足会出现“零能量模式”——单元看起来在变形但应变能为零整个矩阵就奇异了。低阶单元配合减缩积分还可能产生沙漏模式。这些都是程序写出来结果不对的常见原因。我的经验是在新单元代码里做一次能量校验用数值积分的结果和精确积分的结果对比。如果两者差别较大要么积分点数不够要么形函数写错了。用这种方式排错比对着输出结果瞎猜快得多。3.4 组装与边界条件虚功方程在整体层面的含义每个单元的虚功方程都成立对所有单元求和就得到整体虚功方程δU_totalᵀ F_total Σ ∫ δκᵀ EI κ dx注意内部节点上的内力在单元之间是成对出现、相互抵消的这是组装过程为什么可行的重要原因。在程序里组装就是把单元刚度矩阵的元素按照节点自由度编号“对号入座”放进整体刚度矩阵。边界条件的处理也很有深意。固定端节点的虚位移必须为零所以你需要在整体方程中把这些自由度划掉或者用“乘大数法”把对应的对角元素变得极大。从虚功原理的角度看这不过是在说那个自由度上的虚位移不允许出现所以它对应的平衡方程就从方程组里退出了。如果你忘记施加边界条件整体刚度矩阵会有刚体位移模式——相当于结构可以在空间里自由平移、转动而不产生应变刚度矩阵奇异线性方程组的解不唯一。程序会报“矩阵奇异”或者“除数接近零”。这时候第一反应应该是检查约束条件而不是查积分点。4. 单位荷载法求位移虚功原理最经典的算法落地4.1 单位荷载法的物理逻辑现在说回虚力原理最常见的应用——单位荷载法。这个方法的陈述很简单要求结构上某一点沿某一方向的位移就在该点沿该方向施加一个单位虚力然后计算这个虚力引起的内力在真实位移上做的虚功这个虚功在数值上就等于所求位移。为什么因为虚力是单位力它在真实位移上做的外力虚功就是 1 × Δ这个 Δ 就是所求位移。根据虚力原理1 × Δ ∫ (内力虚功的共轭形式) dx对梁结构来说真实荷载产生弯矩 M(x)单位虚力产生虚弯矩 M̄(x)那么Δ ∫ M̄(x) × M(x) / (EI) dx如果是桁架就是Δ Σ N̄ × N × L / (EA)其中 N̄ 是单位虚力引起的轴力N 是真实荷载引起的轴力L 是杆长EA 是轴向刚度。4.2 莫尔积分手算公式与程序积分的对应上面的式子就是结构力学教材里的莫尔积分也叫单位荷载法积分式。手算时我们会分段积分因为弯矩表达式在荷载突变处要分段。程序里做这件事就更直接了沿杆件取积分点在每个积分点分别计算真实弯矩 M(x) 和虚弯矩 M̄(x)然后做数值积分。我在程序里实现单位荷载法时通常不重新求解整体方程而是复用一个已经求好的单元内力结果先做一次完整的结构分析得到所有单元在真实荷载下的内力M、N、V。在需求位移处施加单位虚力再做一次完整的结构分析得到所有单元的虚内力M̄、N̄。对每个单元沿长度做积分累加得到位移。对线弹性结构这两次分析都可以快速完成代价很小。但它的价值很大它是独立于刚度矩阵求解的一条验证路径可以有效检验第一遍求解是否出错。4.3 算例简支梁跨中挠度的程序化实现举个例子。一个跨度 L 的简支梁承受均布荷载 q求跨中挠度。手算结果是 5qL⁴/(384EI)这是结构力学教材里的经典结论。用单位荷载法程序化实现步骤是这样第一次分析真实荷载下跨中弯矩表达式是M(x) qLx/2 - qx²/2第二次分析在跨中施加向下的单位虚力跨中虚弯矩表达式是M̄(x) x/2 当 0 ≤ x ≤ L/2 M̄(x) (L-x)/2 当 L/2 ≤ x ≤ L然后积分Δ ∫₀ᴸ M̄(x) × M(x) / (EI) dx把两段分别积分再相加结果就是 5qL⁴/(384EI)。我在程序里写过这段逻辑核心伪代码如下def midspan_deflection(L, q, EI): # 真实弯矩函数 M lambda x: q * L * x / 2 - q * x * x / 2 # 单位虚力作用下的弯矩函数 M_bar lambda x: x / 2 if x L / 2 else (L - x) / 2 # 分段高斯积分每段取 4 个积分点 def gauss_integral(a, b, n4): # 标准高斯点和权重 points, weights gauss_legendre(n) mapping lambda xi: 0.5 * (a b) 0.5 * (b - a) * xi return 0.5 * (b - a) * sum( w * M_bar(mapping(xi)) * M(mapping(xi)) / EI for xi, w in zip(points, weights) ) return gauss_integral(0, L / 2) gauss_integral(L / 2, L) print(midspan_deflection(6.0, 20.0, 2.0e4)) # 输出 1.875e-3 之类的结果这里分段是因为 M̄(x) 在跨中有一阶不连续。若不分段数值积分也能给出近似值但分段更精确、更稳妥。这个细节就是程序实现和手算思路相互印证的地方。4.4 虚力选取的注意事项用单位荷载法时虚力的方向、位置必须与所求位移的方向、位置严格对应。求绝对竖向位移就加竖向单位力求转角就加单位力偶求两点的相对位移就在两个点各加一个方向相反的单位力。这些规则说起来简单但写程序时很容易因为坐标方向约定不同而出错。我推荐一个排错方法用单位荷载法计算一个已知解析解的简单问题比如悬臂梁自由端挠度、简支梁跨中挠度。如果连这个都对不上先别急着查积分代码检查虚力的施加方向和真实内力的符号约定是否一致。符号问题在虚功计算里是排第一位的 bug 来源。5. 程序实现中的关键细节与常见坑5.1 虚功不是“假功”符号约定要统一写虚功相关的代码最容易踩的坑是符号约定。虚功原理本身对符号不挑剔——只要你保持一致正负号最后会自己摆平。但程序里不同模块的符号约定很容易互相冲突。比如我前面那个反面教材我用了一个经典的梁单元刚度矩阵但我的自由度排列是 [v₁, v₂, θ₁, θ₂]而不是 [v₁, θ₁, v₂, θ₂]。矩阵里的 6L、-6L、2L²、4L² 这些项跟着形函数 B 矩阵的列顺序走自由度排列一变矩阵必须重新排。如果不重新推一遍只是机械地按记忆填数必错无疑。处理办法很简单单元刚度矩阵永远不要手写死而是从形函数出发实时计算。即使为了效率要缓存也要在建矩阵前先打印一个 2×2 或 4×4 的小矩阵验算它是否对称、对角项是否为正再投入到组装流程里。5.2 忽略剪切变形时哪些构件会出问题欧拉-伯努利梁单元只考虑了弯曲应变能没有考虑横向剪切应变能。在虚功方程里这就等于把内部虚功那一项只保留了弯曲贡献忽略了剪切贡献。对细长梁这个假设足够好。但如果梁的高跨比大于 1/10 左右比如深梁、短柱、剪力墙连梁剪切变形的影响就不可忽略了。这时候再用欧拉-伯努利梁单元算出来的位移偏小、刚度偏大程序结果和实测对不上。解决方法是改用铁木辛柯梁单元Timoshenko梁单元。它在虚功方程里额外加入一项剪切应变能δW_int ∫ δκ × EI × κ dx ∫ δγ × k_s × GA × γ dx其中 γ 是剪切应变k_s 是截面剪切修正系数矩形截面取 5/6圆形截面取 6/7 左右。加了这一项之后单元刚度矩阵就不再是原来那个 4×4 的简单形式积分点选择也要重新验证。我的建议是如果你的程序要用于实际工程至少要同时实现欧拉-伯努利梁和铁木辛柯梁两种单元根据构件的长细比自动选择。5.3 刚度矩阵奇异虚位移自由度与约束的关系整体刚度矩阵奇异几乎每个写结构程序的人都遇到过。现象是求解器报“矩阵奇异”或“零主元”。从虚功原理角度理解这是因为整体虚位移中还存在一个非零的刚体位移模式它在任何真实外力下都不产生内部虚功导致虚功方程无法给出唯一解。这是很好记的判断方法如果你的结构有 6 个刚体自由度平面问题 3 个两个平动、一个转动你至少要约束 6 个自由度才能让整体刚度矩阵非奇异。在程序里最简单有效的检查是计算整体刚度矩阵的秩或者看特征值如果有接近零的特征值就说明有未约束的刚体模式。有些程序还会用“最小约束工况”来调试一个被三个不同方向的支座约束住的简支梁如果求解成功再把约束逐步放松就能定位到哪根单元、哪个自由度上的约束漏了。5.4 验证程序的三个实用手段写完成一个分析程序后怎么确认它对我自己的习惯是“三件套”验证第一件单元素精解测试。取一个单元一端固支一端加力位移手算和程序解逐项对比。这一步通过说明单元刚度矩阵和荷载向量没问题。第二件整体能量校验。求解完成后计算外力做的总功 W_ext 0.5 × FᵀU再计算所有单元的应变能总和 U_int 0.5 × Σ(U_eᵀ k_e U_e)。理论上二者必须完全相等。如果不等八成是组装或者边界处理有 bug。能量校验是抓组装错误的好手段。第三件与商业软件对拍。建一个简单的算例比如三跨连续梁、平面刚架用通用有限元软件跑一遍对比节点位移和杆端弯矩。这个测试不需要太复杂几个典型工况就够了。对拍通过之后程序才算真正可以交付使用。除了这三件套还有一个从虚功原理直接引出的检查方法改变单元的划分数量看结果是否收敛。理论上加密网格后位移应该单调逼近精确解。如果结果不收敛或者乱跳通常说明单元本身有问题——比如不满足分片试验或者存在零能量模式。6. 从虚功原理出发还能走到哪里6.1 虚功方程其实是很多高级算法的出发点写到这里有人可能会问虚功原理只是用来推刚度矩阵、求个位移吗当然不止。虚功原理的变分形式是很多更深层理论的出发点。把虚位移原理中“虚位移任意”这一条和“总势能取驻值”联系起来就得到了最小势能原理。在程序里这相当于把求解过程转化为一个优化问题在所有满足约束的位移场中真实位移使总势能最小。有限元方法里的很多误差估计、自适应加密算法都是从最小势能原理出发做的。把虚功方程拓展到动力问题加入惯性力项和阻尼力项就得到了 Hamilton 原理的雏形。这就是模态分析、瞬态动力学程序的根基。我后来写动力分析模块的时候发现只要把静力虚功方程里加上 ρ × ü 项整个程序框架几乎不用大改只是刚度矩阵旁边多了质量矩阵和阻尼矩阵。再说远一点几何非线性问题里平衡方程需要在变形后的构型上建立虚功方程也要写成当前构型下的形式这就导出了更新的拉格朗日格式UL 格式和完全的拉格朗日格式TL 格式。这些内容看起来很高深但底层的那句话还是不变的外力的虚功等于内力的虚功。6.2 我的个人体会把虚功原理当成“第一性原理”做了这些年结构程序我的体会是虚功原理在结构分析中的地位相当于能量守恒在物理中的地位。它不告诉你“某一点的应力是多少”但它给出了“整个结构必须满足的全局约束”。恰恰是这种全局约束在计算机上最容易表达、最容易离散化、最容易做数值求解。如果你正在计划写自己的分析程序我给的建议是先花一周把虚功原理和它的两个分支虚位移原理、虚力原理彻底吃透把梁单元的刚度矩阵自己推一遍把单位荷载法的积分代码写一遍。这两件事一旦完成后面无论遇到什么单元类型、什么复杂工况你都有了从源头推理的能力而不是永远在别人的公式库里找拼图。这也就是这个“结构力学程序算法理论基础”系列的第一篇的目的——把地基打牢。后面再聊单元库的组织、平面问题的等参变换、非线性求解策略都会反复用到这里的概念。到时候你会感谢现在认真推导过的每一行虚功方程。