有限元如何把连续场变成矩阵？从强形式到弱形式的完整主线

有限元如何把无穷维的连续场问题，转换为有限维矩阵方程，同时保留物理平衡，并判断计算结果是否可信？

全文只回答这一个问题。它可以拆成三问：

 

1. 为什么要从强形式变成弱形式？ ——强形式要求每一点都严格平衡，而我们想用的分段多项式在接缝处根本谈不上"平不平衡"：一阶导已经不连续，强形式却还要求二阶导。
2. 为什么弱形式能够变成矩阵方程？ ——把"每点受力平衡"换成"做的功对得上"。功能累加，接缝就不再是障碍；再规定解只由几个节点值定住，方程也就变成有限条。
3. 为什么求解器收敛不等于结果可信？ ——求解器只保证那几条方程被解出来，不保证方程描述的是真实的结构。判据要另外建立：模型、表述、离散、代数、物理验证五层证据，见第 8 章。

有限元最后交给计算机的是一组矩阵方程，但它的起点并不是矩阵，而是连续介质中每一点的物理守恒。两者之间的路线是：

全文分四个部分、九章：

 

• 第一部分（1–3 章）　从物理守恒到弱形式。 把"每一点都要平衡"的微分方程，改写成积分意义下的等价陈述，并交代什么样的函数才有资格参与。这一部分不含任何近似。
• 第二部分（4–7 章）　从弱形式到矩阵方程。 做唯一的一次近似——把解限制在有限维空间里，得到       ；随后用一个两单元算例从强形式一路算到边界反作用量，再把同一条流水线换到结构力学上验证一遍。
• 第三部分（8 章）　结果的可信度。 求解器收敛之后，结果凭什么可信：模型、表述、离散、代数、物理验证五层证据。
• 第四部分（9 章）　主线之外的扩展。 动力学、非线性、不连续 Galerkin、等几何分析、稳定化格式等，都留到主线走通之后再给一张地图。

符号约定 （全文统一）：

第一部分　从物理守恒到弱形式

微分方程要求每一点都严格平衡，这个要求对分段多项式太苛刻。本部分就从这一点出发，说明为什么必须改用弱形式：第 1 章交代要解的是什么问题，第 2 章把方程换成积分意义下的等价陈述并说明这样做换来了什么，第 3 章交代换过之后什么样的函数才有资格参与。到本部分结束，方程已经改头换面，但还没有引入任何近似。

1. 有限元从什么问题出发

连续问题的未知量是一个场 ——数学上就是一个定义在整个区域上的函数。结构力学的位移场     、传热的温度场     、流体的速度与压力     、电磁的电势、声学的声压——这些量在区域内有无穷多个取值，而计算机只能存有限个数。这就是全部困难的来源。

 

有限元的应对。 把物体切成许多小块（单元），规定未知量在每块内部只能是低次多项式——最常用的就是线性，于是整个解由有限个自由度 ——那些多项式前面的系数——完全确定。最常用的一类单元里，这些系数恰好等于节点上的场变量值（原因见第 4 章），因此习惯上直接叫"节点自由度"。至于方程如何跟着改写，是第 2 章的事。

本文使用的模型问题。 全文围绕同一个数学原型展开—— 稳态扩散（势场）问题 ，也就是泊松型方程（一般      可随位置变化，比经典 Poisson 方程更广）。

 

最好理解的例子是稳态传热 ：温度分布不再随时间变化时，流进任意一小块材料的热量，必须等于流出的热量加上它自身产生的热量。 渗流 （未知量是水头）、 静电场 （电势）、 稳态物质扩散 （浓度）写出来是同一个式子——凡是"通量正比于某个量的梯度、且稳态下守恒"的现象，都长这样。

式中各符号：

关于通量的符号 ：传热、扩散、静电等问题中， 物理上向外的通量是    （热从高温流向低温，故带负号）。本文推导一律采用      这个组合，并把      称为广义通量 ——它与物理通量差一个负号，但全程自洽。若要代入实测的物理通量，按第 2 章末尾的约定改号即可。一维杆问题没有这个困扰：     就是轴力本身。

三个式子分管三处：

• 第一式在区域内部      ，要求每一点守恒——净流出等于该点产生的量；
• 第二式在边界的一段      ，直接给定未知量的值       ；
• 第三式在边界的另一段      ，给定穿过边界的通量       （       就是通量在外法向上的分量）。

其中     、    。一维时它退化为

这一个式子同时是两个物理问题，第 6 章的算例会用到它的双重身份：

这个提法称为强形式（strong form） ——"强"指它要求方程在每一个点上严格成立，也因此要求      具有足够高阶的导数。这个要求正是第 2 章要处理的麻烦。

 

两类边界条件提供的信息不同。

这两类在弱形式中的地位也完全不同：     上的条件必须显式施加，     上的条件会在分部积分后自动进入。这一点第 2、3 章会说清。

2. 为什么要把强形式改写成弱形式

强形式对解的要求太高。 式中含二阶导数，意味着      必须处处二阶可导。但现实中材料会分界、截面会突变、载荷会集中，解本来就不那么光滑；更要紧的是，我们打算用的分段线性函数根本达不到这个要求。

 

残差：还有多少没平衡。 把一个候选解代回强形式， 本应处处得零 ——那是真解才有的性质。候选解够不着真解时，代回去总要剩下一点：

物理上，残差就是"净流出减去源"，即某点上还剩多少没抵消掉的量；在结构力学里它就是不平衡力 。

为什么不能逐点检查残差。 取一个分段线性的候选解：它在单元内部是直线（二阶导为零），在节点处一阶导跳变。于是强形式残差不能作为普通函数在单元界面上逐点定义 ——界面两侧通量不等，那里集中着一份没抵消的力，而它挤在零厚度的面上，不是"单位体积的量"，强形式的提问方式装不下它。也就是说，"逐点检查残差是否为零"这个问题对这类函数提不出来 。

提高单元阶次并不能解决这一点。二次单元的二阶导在单元内部存在了，但相邻单元之间仍只有      连续，界面上一阶导照样跳变。 症结在界面连续性，不在多项式次数 （详见附录 2）。

那就换个问法。 逐点行不通，是因为"某一点上平不平衡"这个问题本身依赖二阶导。既然如此，就不要盯着单点，改成把残差摊到整个区域上去称量 ——要求它对测试空间中的每一个允许测试函数，称量结果都为零。注意是"每一个"，不是某一次称出零就算过关。

 

这一步用两个动作完成：先把残差乘上一个函数再积分（ 换问法 ），再对得到的式子做分部积分（ 降低对      的要求 ）。前者解决"怎么问"，后者解决"谁有资格回答"。

动作一：乘      并积分，把"力"变成"功"。

这里的      叫测试函数 （test function），是我们自己挑的一个已知函数 ，不是待求量。它充当探针：拿它去和残差作内积，就是从某一个角度问一次"平不平衡"。     可以任意取（唯一的限制在 Dirichlet 边界上，见第 3 章），取得越多，检验得越彻底。

在力学里      有个现成的名字： 虚位移 ——一个假想的、微小的、不违反约束的位移。这样一来，残差是单位体积上的不平衡力，于是

于是"每一点合力为零"就换成了" 对任意虚位移，总虚功为零 "。两者等价——若某处合力不为零，总能挑一个集中在那里的      把功做出来。这一步之所以有意义，是因为功可以累加、可以积分，而"逐点为零"不能 。

动作二：分部积分，把虚功拆成内、外两部分。 从

 

出发，用散度的乘积法则

移项、对      积分，再对最后一项用散度定理：

边界项按      拆开。

• 在        上，测试函数按定义取零（理由见第 3 章），该项整体消失—— 这正是        上未知的边界反作用量不必事先知道的原因 （结构中即支反力，传热中即维持恒温壁面所需的热流）；
• 在        上，       已由边界条件给定，直接代入。

于是得到弱形式 ：

左边是内部项——结构力学里就是应变乘应力，即内力虚功 ；右边是体力与边界力做的外力虚功 。弱形式读作内力虚功 = 外力虚功 ，这正是虚功原理。

"降阶"降的是什么。 不是方程的阶——它还是二阶。变的是两个导数的归属：

于是对      只要求一阶弱导数，     分段多项式合法。这是主要目的，另外白拿两样：边界项自动出现，成为自然边界条件；     与      地位对称。代价是      也必须可导，且方程只在积分意义下成立。

通量符号的约定。 若把物理向外通量定义为     ，则边界项应写成     。正负号取决于通量定义，推导时必须保持一致。

其他方程（线弹性、不可压、对流扩散、四阶的梁板问题）走的是同一条路线，但恒等式、边界项含义与分部积分次数各不相同，见附录 1。

3. 什么函数有资格进入弱形式

弱形式里反复出现"对所有      成立"。要把这句话说清楚，得先讲明白      和      各自能从哪些函数里挑。

能量有限：    。

即      且其一阶弱导数    。泊松问题的能量含     ，在结构力学中      对应应变，这一项就是应变能。所以      的物理含义是： 位移不能到处发散，应变不能到处发散，储存的总能量必须有限。

它管的是积分而不是每一点，所以不要求经典导数处处存在——这正是分段多项式能够进入的原因。但允许到什么程度需要说准 ：

•       连续、       在单元界面上有限跳变： 可以属于       ；
•       本身跨界面发生跳跃： 通常不属于       ，因为这样的间断会产生分布型导数。

孤立点上修改取值不影响函数的等价类，但跨界面的间断不是"孤立点"。

试函数空间：谁有资格当候选解。 以一维问题、两端给定     、     为例：

 

上面式子的意思是：     是一大把曲线的集 合——左端全都钉在 0、右端全都钉在 1（这是两个边界条件），中间形状各异，但能量都必须有限（这是     ）。真实解就是这堆曲线中的一条。一般区域上写成

测试函数空间：允许怎样去"拨动"它。 设     ，给它一个微小扰动     。在      上扰动后仍须满足     ，而原来的      已满足     ，所以

物理上：     是真实位移，     是虚位移。固定端的真实位移被约束住，固定端的虚位移也必须为零——你不能"虚拟地"去挪动一个焊死的支座。

为什么连"微小变化"也不允许。     是给定的数据，不是待求的量，这是硬约束，没有"近似满足"这一档。只要     ，对任何     都有     ——它不是"稍差一点的候选"，而是直接掉出集 合     。     只能让偏离变小，不能让它变成零。

用那把曲线来想：     里所有曲线的左端都钉死在同一点，允许的移动只能是"保持钉死、改变中间形状"，方向必然取自     。（数学上说，     不是线性空间而是仿射集 ：取一个满足边界条件的     ，则     。）

限制只在     。 在      上和区域内部，     完全自由——那里给定的是通量而非位移，边界值本来就是待求的。这个不对称直接带来第 2 章的结果：     上的边界项因      而消失，边界反作用量不必事先知道就能推进求解；反过来， 想把它反求出来，恰恰要用一个在      上不为零的      做后处理 （第 6 章算例会用到）。

连续弱问题。 至此可以完整写出：寻找     ，使得

 

读作： 对每一种允许的虚拟变化     ，内部响应都等于外部作用。 不同的      是不同的探针——有的只看局部一小块，有的看整体走势；对每一根探针都平衡，     就是连续弱解。

注意这一步还没有做任何离散近似。 从给定强问题建立与之相容的连续弱问题，不是有限维近似；在适当的正则性与边界条件下，强解与弱解相容。有限元的离散误差要到第 4 章、把无限维空间替换成有限维空间时才引入。

为什么普通有限元只要求     。 两个一维线性单元拼成的函数

在公共节点处函数值连续，但左右斜率可以不同（    ），于是     、    。

物理上： 位移不能断开，但应变可以跳。 材料不能撕裂，所以      必须连续；而相邻单元的应变本就允许不同（想想两段刚度不同的材料交界处），所以斜率跳变可以接受。只要各单元内梯度有限，就有     ，故     —— 能量有限，够格参赛 。但其经典二阶导数在节点处不存在，一般不属于     。

这正是降阶的意义： 强形式要      的二阶导数，折线给不出；弱形式只要一阶弱导数，折线就合法了。

第二部分　从弱形式到矩阵方程

全文唯一的一次近似在这里发生：把解限制在有限维空间里。第 4 章说明这个空间怎么造、凭什么判断解得好不好；第 5 章把它落到实处——单元矩阵、组装、边界条件、求解；第 6 章用一个两单元算例从强形式一路算到边界反作用量，闭合整条主线；第 7 章把同一条流水线原样搬到结构力学上，验证它不依赖于具体的物理问题。

4. 有限维近似与 Galerkin 方法

网格、单元与节点。 计算区域被划分成有限个单元：    。常见单元包括一维杆与梁单元，二维三角形、四边形单元，三维四面体、六面体单元，以及板、壳、界面和接触单元。

形函数。 以长度      的一维两节点线性单元为例：

它们满足节点插值条件     ——正是这个条件保证了自由度      就是节点      处的场变量值本身。形函数做三件事：在单元内部插值场变量、建立自由度与连续场的联系、决定有限元近似空间及其收敛能力。

 

有限元空间。 第 3 章的试函数空间     （候选解住的地方）与测试函数空间     （探针取自的地方）都是无限维的——候选解有无穷多条，探针有无穷多根。有限元各取其中有限维的一部分：

这里的      要求      里每个函数都真正够格：能量有限、在      上取到规定值。做到这一点的叫协调有限元 ，本文主线就是这一类。非协调元与不连续 Galerkin 故意放弃这一条以换取别的好处，见第 9 章与附录 4。

    由分段多项式张成，两者的差别仍在 Dirichlet 边界上：     中的函数满足     ，     中的函数在那里取零。离散近似写作

Galerkin 方法：一个需要说明的选择。 试函数空间是"解的候选形状"，测试函数空间是"检验用的探针"——本是两件不同的事，原则上可以各选各的。探针取 Dirac 函数（只在若干个点上检查）就成了配点法；取     （残差对各自由度的导数）则给出最小二乘法。

 

Galerkin 法的选择是：探针就用形函数本身，两者取同一组     。 三个理由：

• 数量正好 ：       个形函数当探针，得到        个方程，配        个未知数；
• 力学上本就如此 ：       是虚位移，而虚位移是位移的变分，只能住在位移住的地方（第 3 章）；
• 结果好 ：       与        地位对称；当双线性型        本身对称时离散矩阵也对称，且解具有最优投影性质（第 8 章）。

需要限定： 仅仅采用 Galerkin 方法并不能保证任意问题的矩阵都对称 。矩阵对称来自      本身对称、且试函数与测试函数取自同一空间这两个条件的叠加。对流占优问题会故意打破这种对称（Petrov–Galerkin／SUPG，见附录 4）。

离散问题。 现在可以写出有限元真正在解的问题了。回忆第 3 章的连续弱问题——寻找     ，使      对一切      成立，其中      是内部响应、     是外部作用。

把其中的两个空间换成有限维的，就得到离散问题 ：寻找     ，使得

两者逐项对照，只改了三个地方：

第三行是关键： 双线性型与载荷项一个字都没动 ，物理内容原封不动。变的只是"在哪里找解"和"用多少根探针检验"。

而"对一切      成立"这句话，实际上只是      个条件。因为任何      都是形函数的线性组合，只要在每个形函数上分别成立，对整个空间就自动成立：

 

把      代进去，未知数正是      个自由度     。 这已经是一个代数方程组了 ，第 5 章只是把它按单元具体算出来而已。

换成"残差"的说法。 定义弱残差泛函

它衡量"外部作用减去内部响应"还差多少，是一个作用在测试函数上的泛函，而不是一个逐点取值的函数——这样定义正是为了避开第 2 章那个麻烦：强形式残差      在单元界面上并不是普通函数。于是离散问题可以写成

它与      一般并不相等——后者要用到      的二阶导，而      全程只用一阶导（两者的确切差别见附录 3）。习惯上仍把上式读作"残差与测试空间正交"。

正交不等于处处为零。 这是全文最重要的理论结论之一：上式只保证残差在有限维测试空间上的投影为零。举个最简单的：在      上取

虽然     ，但      显然不处处为零——只是正负部分在积分中抵消了。     这根探针只能称出"总量"，称不出分布。

探针有限，就一定有看不见的东西，这正是离散误差的来源 （误差究竟有多大，第 8 章用 Céa 引理给出定量的界）。反过来，只有探针足够丰富才能推出残差为零：若      对一切      成立，则可取     （需     ），得     ，从而      几乎处处。

三个层次必须区分。

后两层差在两处：

空间离散误差就在最后这一步引入 ，它同时截断了两样东西：解的搜索范围从      缩到     ，检验的角度从无穷多根缩到      根。网格加密或阶次提高，等于两条腿一起扩张，离散解才有可能收敛到连续弱解。（工程计算中除此之外还有模型误差、积分误差、迭代误差和数据误差，见第 8 章。）

5. 从单元方程到整体求解

把      代入弱形式，依次取     ：

但这个全局积分不便直接计算。实际流程是：

第一步：单元矩阵。 积分按单元拆开，在每个单元上计算

 

其中      是形函数的梯度矩阵，作用是把节点自由度换算成梯度（第 6 章 ③ 有一维线性单元的具体推导）。这个      的形式贯穿所有问题，只是中间的材料量随物理问题改变（第 7 章）。

第二步：等参映射与数值积分。 实际单元形状不规则，把参考单元映射到物理单元：

 

几何与场变量使用同一组形函数，称为等参思想 。导数经雅可比矩阵变换：

单元积分用高斯积分：

（取绝对值以不依赖单元定向；实际网格通常约定各单元定向一致、    ，此时绝对值可省。     变号意味着单元翻转，属于网格错误。）

积分方案与单元稳定性密切相关：积分点不足会产生零能模态（沙漏），某些问题中全积分会导致锁死，减缩积分与选择性积分需与具体单元理论配套（机理见附录 4）。

第三步：组装。 按局部自由度与全局自由度的对应关系累加：

 

    是自由度映射矩阵 ，记录"这个单元的第几号自由度，是全局的第几号"。以三节点两单元为例，单元 2 连接全局节点 2、3，于是

也就是把      的单元矩阵搬到      中正确的位置，其余补零。求和时，公共节点（这里是节点 2）上两个单元的贡献自动叠加。

实际程序里不会真的构造    ——那太浪费。代码只存一张编号表，按索引把      的元素累加到      的对应位置。上面的写法只是把这件事表达成矩阵运算，便于推导。

物理意义是：相邻单元共享节点自由度、公共节点处满足位移兼容、各单元在公共节点上的内力累加、最终形成整体平衡。由于一个节点只与邻近单元的节点直接耦合，     是稀疏矩阵。

第四步：边界条件。 Neumann 条件（自然边界条件）不需要像 Dirichlet 那样去约束自由度，但边界积分    仍要照算并组装进    ——"自然"指的是它自动出现在弱形式里，不是指它不用算。Dirichlet 条件（本质边界条件）则根本没出现在弱形式里 。对照弱形式就能看清两者的差别：

 

    是方程右端的一项，算进      就完事；     却写在候选解的范围里（第 3 章），式子本身不含它。既然如此，解出来的      就没有任何机制保证边界取值正确，必须在代数层面显式施加：直接消元、矩阵行列修改、罚函数法、拉格朗日乘子法、Nitsche 方法。

对本文这类纯扩散问题，未施加 Dirichlet 条件前      是奇异的——零空间是常数场（整体温度可以任意平移）；未约束的线弹性问题同理，零空间是刚体模态。这不是错误，而是物理事实的正确反映：没有参考点，解只能确定到相差一个常数或一次刚体运动。（并非所有问题都如此，例如含 Robin 边界或反应项时刚度阵本身就可能非奇异。）

第五步：求解    。矩阵性质全部继承自前面的选择：对称来自      对称加 Galerkin，正定来自应变能加足够约束，稀疏来自形函数的局部支撑。中小规模用直接法（Cholesky／    ，一次分解可多次回代），大规模用迭代法（预条件共轭梯度、多重网格）。

 

这一步的展开不在本文范围内。直接法与迭代法的原理、预条件的作用、条件数为何随网格加密而恶化，《矩阵求解》系列有系统讨论，可以参看。

6. 贯穿算例：两个单元，从强形式算到边界反作用量

前五章把整条路线讲完了，但还没有算出过一个数。这一章用第 1 章那个一维稳态扩散问题走完全程——从强形式出发，经弱形式、单元矩阵、组装、边界条件，一直算到节点值、单元通量和边界反作用量，最后与解析解逐项对照。

 

取常数     、常数分布载荷     ，左端固定、右端受给定通量 ：

按前面的对照表，这既是"左端恒温、右端给定热流的导热杆"，也是"左端固定、右端受轴向力     、承受均布载荷的拉杆"。

① 弱形式。 按第 2 章，寻找     ，使得对一切     ：

这一步还没有离散 ：     仍是无穷维空间      中的真解，     也遍历整个     。式子与原方程等价，只是换了个说法。注意右端的     ——Neumann 条件已经自动进来了。

② 离散。 取两个长度      的线性单元，节点位于     ，自由度     。

③ 单元矩阵。 先算     。单元内     ，两个形函数是     、    ，求导得     、    ，于是

    就是把节点自由度换算成梯度的那个矩阵 （结构力学中即应变—位移矩阵）。这里它是常数，因为线性单元的斜率在单元内不变。代入     ：

均布载荷平分到两端节点。

④ 组装。 两个单元的      完全相同，区别只在它们占据全局矩阵的哪几行几列。单元 1 连节点 1–2，单元 2 连节点 2–3，各自摊到      里：

节点 2 被两个单元共用，所以    位置上    ；节点 1 与 3 各只属于一个单元，保持为 1。载荷向量同理，节点 2 处两个单元各分到     ，合为     ：

注意      只出现在节点 3——它来自弱形式右端的      一项，而     、其余形函数在该点为零。 Neumann 条件就这样自动变成了一个节点载荷。

⑤ 施加 Dirichlet 条件。    ，划去第一行第一列，剩下

⑥ 求解。

⑦ 与解析解比较。 精确解为

代入      与     ，得到的正是上面的     、    —— 节点值一分不差 。这是一维常系数问题的特殊性质（此时有限元解在节点处与真解重合），二三维不会这么巧，但足以说明流程没跑偏。

 

⑧ 回算通量。 由      求导得单元内通量（结构中即轴力），每个单元内为常数：

精确通量为     ，是线性变化的。用单元中点去比：     处精确值为     ，     处为     ——与上面两个单元常数值完全吻合。

这不是巧合，而是超收敛 ：在合适的条件下，导数类结果在单元内某些特定点上的精度高于其他位置。但要注意限定——本例是常系数、均匀一维网格、真解为二次多项式且积分精确，条件相当理想； 一般问题中超收敛点的存在与位置都依赖于单元类型、网格规则性与解的光滑性，不能想当然地推广 。

顺带澄清一个常见说法：商业软件在积分点上算应力，首要原因是本构积分本来就在那里执行 （塑性、损伤等都以积分点为状态载体），而不单是精度考虑。节点应力通常由外推、恢复或多单元平均得到，其可靠性取决于单元、网格与所用的恢复方法（如 ZZ 恢复，见附录 3）。

 

⑨ 回算边界反作用量。     上那个未知量（本例中即左端固定处的量：结构读法下是支反力 ，传热读法下是维持该处温度所需的热流 ）对应被划掉的那一行。回到组装后的完整方程，第一行为

代入      与     ，得

校验整体守恒 ：外部输入总量为分布源      加端部通量     ，     恰好等大反号——热学上读作"流进多少就得从恒温端带走多少"，力学上读作"外载荷与支反力平衡"。这正是第 3 章那句话的兑现—— 求解时不需要知道它，求解后用一个在      上不为零的检验函数（这里就是被划掉的那一行）把它算出来 。这个校验也是第 8 章"守恒检查"最基本的一条。

7. 结构力学有限元是同一框架的向量版本

第 1 到 6 章全程用的是一个标量问题。而工程中最常见的有限元分析是结构力学，未知量是向量位移场。那么前面那套流水线要不要推倒重来？

不要。 换的只是符号，流程一步不改。 下面这张表就是全部的过渡工作：