1. 为什么程序里要用虚功原理结构力学这门课很多人学的时候觉得虚功原理是个很玄的东西虚位移、虚应变、内力虚功、外力虚功……概念一堆考试的时候会套公式就行。但等你真的动手去写一个结构分析程序或者去读开源代码的时候会发现虚功原理不是玄学它是整个程序算法从“物理世界”跨到“数字世界”的那座桥。说白了程序没法直接解“梁在荷载作用下的变形”这种物理问题。程序能干的只有一件事——解方程组。那关键问题就变成了怎么把力学问题翻译成方程虚功原理就是干这个翻译工作的。我记得自己第一次看有限元程序源码的时候看到别人直接算单刚矩阵、总刚矩阵觉得这些矩阵就像从天上掉下来的一样。后来追到源头才发现所有那些看起来很“矩阵化”的东西追根溯源都站在虚功原理这个起点上。那虚功原理到底说了个啥一句话版假如一个结构处于平衡状态那么对任意满足约束条件的虚位移外力做的虚功等于结构内力做的虚功应变能增量。这玩意儿表面上是为了求平衡条件但它比直接把平衡方程写出来要高明得多。为什么因为直接写平衡方程你得把每个节点、每个杆件的力都拆开分析那叫“矢量力学”的思路。而虚功原理是“能量法”的思路它不管中间那些复杂的力传递路径只看“输入的能量”和“存进结构里的能量”是不是相等。这种思路特别适合计算机——你不需要考虑结构几何形状有多复杂只需要建立一个统一的能量表达式然后让它对某个东西的变分为零方程自然而然就出来了。从程序算法的角度看虚功原理还有一层更实际的价值它天然就是“离散化”的起点。有限元法的思想是把连续体切成有限个小单元每个单元假设一个位移场然后组装起来。虚功原理在这里的用法就是先在每个单元上写出虚功方程再把它组装成整个结构的方程组。这一步一步的操作就是程序的骨架。所以说这篇“理论基础”系列的第一篇我把虚功原理放在最前面不是因为它最简单恰恰是因为它最难替代。后面不管是要讲直接刚度法、能量原理还是有限元、动力分析、非线性分析虚功原理都是地基。地基没打牢后面盖的楼越高越危险。2. 虚功原理的核心概念拆解2.1 虚位移的理解很多人刚接触虚位移的时候总是纠结那个“虚”字。它不是说这个位移是假的、不存在的而是说这个位移是“虚拟给的”不代表真实的变形过程。我打个比方。你手里拿着一根竹竿竹竿上放着一个球。现在竹竿是弯的球在某个位置保持平衡。这时候你在脑子里“想象”一下——注意只是在脑子里——把球往旁边挪一点点然后松手。这个“想象的微小挪动”就是虚位移。关键点在哪里这个虚位移是在满足约束条件的前提下的任意微小位移。也就是说你不能往墙里面挪也不能挪到天上只能在结构允许的自由度方向上挪。而且因为是“虚”的它不需要满足力学的平衡方程也不需要满足应力应变关系它只需要满足几何协调条件——也就是不能把结构扯断。从数学上讲虚位移就是位移场的变分记作δu。这里的δ是变分符号不是delta。它表示“一个无限小的、任意的位移增量”。那我们为什么需要这样一个虚拟的东西因为现实世界里结构在真实荷载作用下的位移是固定不变的——是多少就是多少你没法让它变也就没法拿它来做“对比实验”。但虚位移是任意的我们可以任意指定然后看“如果我给这个结构这样一个小位移结构的能量会怎么变化”。这种“思想实验”给了我们一个探测结构状态的手段。2.2 外力虚功与内力虚功有了虚位移接下来就是给它配对。虚功 力 × 虚位移就这么简单。外力虚功结构上所有外力包括集中力、分布荷载、支座反力在虚位移上做的功的总和。内力虚功结构内部应力在对应虚应变上做的功的总和。这里有个容易搞混的地方。很多人觉得“内力虚功”这个词很奇怪——内力不是成对出现的吗它们自己互相抵消还能做功没错内力是成对出现的它们确实平衡但它们分别在各自那一侧的虚位移上做功这功是对“变形”而言的不是对“刚体位移”而言的。再打个比方。你双手拉一根橡皮筋拉到一定长度保持平衡。这时候你在脑子里想象把皮筋再拉长一点点你的手做了功同时皮筋内部的分子链也“储存”了更多的弹性势能。外力做的功 储存进去的势能增量。这就是虚功原理最朴素的理解。在公式上它长这样δW_ext δW_int其中δW_ext是外力虚功δW_int是内力虚功通常也叫虚应变能。这个等式成立的条件就是结构处于平衡状态。2.3 虚位移原理与虚力原理的区别其实虚功原理有个“双胞胎兄弟”——虚力原理。虚位移原理给结构一个虚位移外力虚功等于内力虚功等价于满足平衡条件。虚力原理给结构一个虚力虚拟的平衡力系虚功等式等价于满足变形协调条件。这在程序算法里是个要紧的选择题。有限元法主要用到的是虚位移原理所以后续推导基本都围绕它进行。但有的时候比如做杂交单元、平衡单元、或者某些柔性体分析的时候虚力原理也有用武之地。虚位移原理的数学根源是“最小势能原理”——真实的位移让总势能取最小值。虚力原理的根源是“最小余能原理”——真实的应力让总余能取最小值。这俩在拓扑优化、结构优化领域也非常重要因为你要么用位移变量要么用力变量选哪个直接决定优化算法长什么样。对初学者记住这句话就够了虚位移看平衡虚力看协调程序里绝大多数情况用虚位移原理。3. 从虚功方程到单元刚度矩阵的完整推导3.1 平面梁单元的虚功方程写到这里我们开始动真格的了。你最终是要写代码的所以我们来走一遍从虚功原理到刚度矩阵的完整推导这是所有有限元程序的“第一性原理”。以最简单的平面梁单元为例。假设一个梁单元有两个节点每个节点有竖向位移w和转角θ两个自由度一共四个自由度。其位移场用插值函数表示通常用Hermite插值多项式因为梁的位移场需要保证位移和转角都连续w(x) N1*w1 N2*θ1 N3*w2 N4*θ2这四个形函数分别是N1 1 - 3*x^2/L^2 2*x^3/L^3 N2 x - 2*x^2/L x^3/L^2 N3 3*x^2/L^2 - 2*x^3/L^3 N4 -x^2/L x^3/L^2这些形函数可能看起来有点复杂它们的本质就是“保证单元端点位移和转角精确匹配节点自由度”的插值基函数。现在给结构一个虚位移场也就是给节点的四个自由度一个虚拟增量δw(x) N1*δw1 N2*δθ1 N3*δw2 N4*δθ2梁的虚应变这里是虚曲率就是位移场的二阶导δκ d²(δw)/dx² B * δd其中B矩阵几何矩阵是形函数二阶导组成的行向量δd是节点虚位移向量。梁的曲率与挠度的关系其实就是材料力学里面的小变形假设。外力虚功包括节点力做的虚功 分布荷载做的虚功。写成公式δW_ext δd^T * F ∫ δw(x) * q(x) dx其中F是节点力向量q(x)是分布荷载。内力虚功虚应变能等于应力弯矩在虚曲率上做的功δW_int ∫ δκ * M dx ∫ (B*δd)^T * (EI * B*d) dx这里用到了弯矩和曲率的关系M EI * κ EI * B * d。EI是抗弯刚度。3.2 刚度矩阵的形成根据虚位移原理δW_ext δW_int即δd^T * F ∫ δw * q dx ∫ (B*δd)^T * (EI*B*d) dx因为虚位移δd是任意的所以两边可以约掉δd^T就得到了单元的平衡方程F [∫ B^T * EI * B dx] * d - ∫ N^T * q dx这个方括号里的积分就是单元刚度矩阵k ∫ B^T * EI * B dx对等截面直梁单元EI是常数B是x的简单函数这个积分可以直接算出来得到经典的4×4梁单元刚度矩阵k EI/L^3 * [ 12 6L -12 6L ] [ 6L 4L² -6L 2L² ] [ -12 -6L 12 -6L ] [ 6L 2L² -6L 4L² ]是不是很眼熟结构力学教材里直接给的梁单元刚度矩阵就是这么来的。只不过教材里直接给结果而程序里你需要理解它是“从虚功方程积分得来的”。第二步等效节点荷载就是后面那个积分项f_eq ∫ N^T * q dx它可以算出来比如均布荷载q作用下等效节点力就是ql/2的剪力和ql²/12的弯矩分别加在两端节点上。这就是为什么程序里处理分布荷载时不需要把荷载细分到很多节点上直接等效到节点就行了——虚功原理保证了这种等效在能量意义上是精确的。这一整套流程从虚位移、应变、应力、虚功方程到最终得到刚度矩阵和等效荷载向量——这就是有限元程序的“最小核心循环”。你写任何一个线性静力有限元程序核心步骤都是这一套只是单元的几何形状变复杂了形函数变多了而已。3.3 程序里为什么不用解析积分你可能注意到刚才算了积分好像不太难。但实际问题中单元的形函数可能非常复杂比如高阶单元、等参单元、板壳单元这时候被积函数已经复杂到没法手算解析解了更别说还要考虑材料非线性、几何非线性。所以程序里几乎无一例外地采用数值积分最常见的是高斯积分。核心思想很简单不试图算出积分的解析表达式而是在单元内选取几个采样点高斯点算这些点的函数值加权求和逼近积分值。对上面的梁单元刚度矩阵理论上一维高斯积分取两个积分点就足够精确了。对二维四边形单元通常取2×2个高斯点对三维六面体单元取2×2×2个高斯点。程序里其实就是几层循环初始化 k 0 对于每个高斯点 (ξi, ωi): 计算该点的雅可比行列式 |J| 计算该点的 B 矩阵 计算该点的材料矩阵 D (各向同性弹性就是常量) 累加 k B^T * D * B * |J| * ωi这种做法的好处是通用性极强。你只要换一个单元类型改一下形函数和积分规则刚度矩阵的计算流程一个字都不用改。这就是商业有限元软件能够在一个框架里塞进几十种单元的秘密。4. 虚功原理在程序算法中的实际地位4.1 直接刚度法里怎么用上虚功原理说到结构力学程序很多人第一个想到的是直接刚度法——也就是先建立单元刚度矩阵再对号入座组装成总刚矩阵然后引入边界条件求解。这套流程也是用虚功原理打底的。组装的过程其实不需要再回到虚功方程本身只需要按节点自由度编号“对号入座”就行。但你要明白为什么可以直接组装因为虚功的叠加性——整个结构的总虚功等于各个单元的虚功之和。这就好理解了。整个结构的虚功方程可以拆成每个单元分别写虚功方程然后加起来。单元之间相邻节点的力是内力对它们的虚功互相抵消所以最后剩下的就是外力虚功和所有单元的虚应变能相等。也就是说结构总的刚度矩阵就等于把所有单元刚度矩阵按自由度映射叠加起来——程序里那一堆“组装循环”的背后还是虚功原理。如果你写程序的时候组装逻辑写错了比如自由度编号对应错了、单元定位向量弄反了那总刚矩阵就是错的。这时候调试起来很痛苦但如果你脑子里有“虚功原理要求所有单元的贡献按自由度对齐叠加”这个意识就很容易定位问题——组装的本质是能量的叠加不是力的叠加。4.2 位移边界条件的处理虚功方程有一个隐含的前提虚位移必须满足位移边界条件也就是支座处的虚位移为零。不然的话外力虚功里就会多出支座反力做的功方程形式就变了。程序里怎么处理最粗暴也最常见的做法是划行划列法——把总刚矩阵中约束自由度的行列直接删掉对应的荷载项也删掉。这种做法实现简单、求解效率高适合处理“零位移”约束。但有时候约束条件是“支座沉降”——给定位移不为零。这时候就不能简单删行删列了而要用罚函数法或者拉格朗日乘子法。罚函数法在总刚矩阵上直接加一个很大的数罚数来逼近约束条件实现起来很简单但会影响矩阵的条件数拉格朗日乘子法引入额外的未知量乘子物理意义更清晰实际效果也更好但要增加系统自由度。处理在程序中是一个基础但容易出错的环节。我见过不少人总刚矩阵算得对但一处理边界条件就报错或者约束写错了导致结构出现机构位移——这时候你回头想想虚功原理就会意识到约束错了等于虚位移空间选错了整个能量等式自然就不成立了刚度矩阵奇异也就是必然的。4.3 组装与求解的整体流程把前面的步骤串起来一个线性静力结构分析程序的主流程就清楚了1. 定义节点和单元几何数据 2. 计算单元刚度矩阵 k_e虚功积分 3. 计算单元等效节点荷载 f_e 4. 组装总刚矩阵 K对号入座 5. 引入边界条件位移约束、荷载 6. 求解线性方程组 K*U F 7. 后处理计算单元应力、弯矩、轴力等其中第2、3步最依赖虚功原理第4步的物理意义也是虚功的叠加性第5、6步则是在解虚功方程离散化之后的线性方程组。整个流程环环相扣一环错了后面全错。所以我一直觉得学结构力学程序算法第一件事不是背流程而是把虚功原理这一环的物理意义彻底想清楚。5. 常见问题与研习建议5.1 新手最容易踩的几个坑第一个坑把虚功原理和能量守恒搞混。能量守恒说的是“总能量不变”虚功原理说的是“平衡状态下外虚功等于内虚功”。两者有关系但不是一个东西。虚功原理不是一个普适的能量守恒律而是平衡条件的另一种数学表达。第二个坑数值积分阶次选低导致刚阵奇异。这是程序里特别容易出的问题。对梁单元有些新手用1个高斯点算刚度矩阵结果发现刚阵秩不够结构出现了零能模式。原因就是积分点太少没法精确积分形函数的二次项。经验法则是积分阶次必须能精确积分刚度矩阵被积函数中形函数导数的最高阶次。对梁单元至少2个高斯点对平面4节点等参单元至少2×2个高斯点否则你调试到天亮都找不到问题。第三个坑分布荷载漏算等效节点力。新手在算框架的时候经常只算了节点集中荷载忘记把梁上分布荷载等效到节点上。结果就是结构位移偏小很多而且位移模式完全不对。虚功原理告诉你分布荷载必须通过∫N^T*q积分等效到节点上才能保证能量等价。第四个坑忘了检查虚位移空间。如果一个结构没有足够的约束总刚矩阵是奇异的程序直接崩掉。很多时候不是你程序写错了而是你的“虚位移”里包含了刚体位移——这在物理上是不允许的在数学上就是K奇异。排查方法是用特征值分析看看刚度矩阵有没有零特征值。5.2 调试与验证的建议我写结构程序有一个习惯每写完一个功能必须先跑标准算例验证。比如拿一个两端固支梁受均布荷载的标准解析解去对比程序输出结果。这种“基准测试”benchmark的思想我强烈建议每个人都要养成。具体来说你可以建立这样几个测试用例悬臂梁端部受集中力校核挠度和转角简支梁受均布荷载校核跨中挠度和弯矩两端固支梁受均布荷载校核固端弯矩和跨中挠度这些算例在结构力学教材里都有解析解算完拿来和程序结果对比误差在1%以内基本就可以认为刚度矩阵和荷载向量没写错。如果误差大优先检查单元刚度矩阵的符号和数值——这是最常见的错误源。另外验证虚功方程本身是否满足也是一个强大的调试手段。程序算完之后可以输出每个单元的内力虚功和外力虚功之和检查它们是否相等。如果不相等说明要么是单元刚度不对要么是荷载向量不对要么是组装有错。这种“能量检验”比单纯看位移值靠谱得多因为能量是全局量有一点小错都能放大看出来。5.3 学习路径的展望虚功原理学的不仅仅是虚功原理本身更重要的是通过它你建立了一种“能量思维”。这种思维会让你在后面的学习中越来越顺。后续会在这个系列里继续展开的包括最小势能原理与变分法基础——虚功原理的“升级版”让你明白为什么结构总是“选择”让势能最小的那条变形路径。这是理解有限元、也理解各种优化算法的必要基础。有限元方法与形函数理论——从虚功方程到任意单元的构造方法深入理解为什么有限元解是“近似解”以及如何提升精度。动力分析与模态分析——虚功原理加上达朗贝尔惯性力项就得到了动力学方程再往后做模态分析、谱分析都是从这里出发。几何非线性和材料非线性——虚功方程在非线性情形下的增量形式你学完就能看懂那些强大的非线性有限元程序在做什么。我个人在实际使用中最大的体会就是很多看起来“编程味儿”很重的问题追到根子上都还是力学概念的问题。虚功原理是你手里第一把钥匙把它用透了后面的路会顺很多。