简介一份围绕有限元弱形式主题的 doc 学习文档面向有限元初学者、物理仿真研究人员及工程分析人员。内容系统梳理 PDE 弱形式的基本概念、三种物理问题描述方式偏微分方程、能量最小化形式、弱形式之间的联系并结合弹性静力学 Navier 方程与弹性能量表达式展示从 PDE 到泛函变分再到有限元离散求解的完整思路。文档还特别说明弱形式在非线性多物理场问题中的优势以及 COMSOL、SOL Multiphysics、ANSYS 等软件中的应用场景可帮助读者理解有限元底层数学基础并为使用商业软件的自定义 PDE 建模提供参考。资料压缩包共 1 个 doc 文件大小 649KB内容集中、便于快速通读。已有 128 人学习适合作为研究生课程补充或工程师入门有限元弱形式的精炼材料。1. 有限元的弱形式从“偏微分方程解不出来”到“矩阵方程”很多做数值计算的人第一次被迫面对“弱形式”这三个字是在一份名为《有限元的弱形式.doc》的讲义里前几页还在规规矩矩推导偏微分方程翻过一页突然冒出两个积分和一堆带下标的函数空间再往后就是刚度矩阵、载荷向量和网格剖分。如果你也卡在这个位置上先把最反直觉的结论说出来弱形式不是把方程“变弱”而是把偏微分方程的逐点约束换成积分约束。换完之后解的光滑性要求从“二阶连续可导”降到“一阶导数平方可积”代价是方程从单点严格满足变成在任意测试函数的加权意义下满足。这个交换的直接收益是分片多项式现在可以进入解空间而分片多项式正是有限元刚度矩阵能组装起来的前提。做matlab有限元编程求解实例时第一步永远是从一维泊松方程开始写弱形式再把组装出来的矩阵打印出来和手算结果核对。这篇文章就把这条路完整走一遍先推弱形式再写最小可运行代码最后用误差范数验证程序没写错。适合刚接触有限元的工程师也适合边界条件总处理不干净的从业者。2. 强形式到弱形式分部积分与测试函数怎么选2.1 先写出强形式再乘一个测试函数以一维稳态热传导问题为例控制方程是-d/dx(k du/dx) fx ∈ (0, 1)边界条件取 u(0) u0k du/dx |_{x1} g。所谓“强形式”字面意思是方程在定义域内每一个点 x 上都要严格成立。这要求 u 在闭区间内二阶连续可导。对实际工程中的温度场、位移场来说这个要求相当苛刻两种材料交界面处温度连续但温度梯度可以突变集中载荷作用点附近位移的导数也会间断。强形式在数学上仍然成立但直接求解析解或者构造数值格式都不容易。有限元的思路是不要求方程逐点满足而是先选一个测试函数 v(x)把强形式两边同时乘 v再对全区间积分。测试函数可以理解成一把尺子方程不再逐点测量而是用无数把尺子去量“总体偏差”。当 v 取遍所有足够光滑且在边界上满足特定条件的函数时加权等式与原方程等价。这一步之后微分方程变成了积分方程允许被积函数出现有限个跳跃点因为跳跃点测度为零不影响积分值。2.2 分部积分导数转移给测试函数对等式 ∫(-d/dx(k u)) v dx ∫ f v dx 的左端做分部积分∫ k u v dx - [k u v]₀¹ ∫ f v dx边界项展开是 k u(1) v(1) - k u(0) v(0)。此时要区分两类边界信息右端 x1 给出的是 k u(1) g属于已知自然边界条件可以直接把 g 代入左端 x0 给出的是 u(0) u0但 k u(0) 是未知量这一项必须被消掉。消掉的办法是约束测试函数 v(0)0。本质边界条件上的测试函数取零是弱形式理论里最容易被忽略的规则。取了 v(0)0 后弱形式的最终形式是求 u使得 u(0)u0 且 u、u 平方可积对所有满足 v(0)0 且 v、v 平方可积的 v有∫ k u v dx ∫ f v dx g v(1)右端最后一项 g v(1) 就是 Neumann 边界条件进入有限元方程的入口。如果 g0这一项自然消失不需要做任何额外处理。这就是“自然”两个字的含义它自己掉进弱形式里而不是被人为塞进去的。2.3 “弱”在哪里函数空间与解的光滑性强形式要求 u ∈ C²弱形式只要求 u ∈ H¹。H¹ 是 Sobolev 空间包含所有函数本身平方可积、一阶导数也平方可积的函数。分片线性函数正好落在 H¹ 里它在每个单元内是线性函数导数在单元内为常数在单元边界处跳跃但跳跃值有限平方积分有限。对比项强形式弱形式对待定函数的要求C²二阶连续可导H¹函数和一阶导数平方可积边界条件本质与自然边界都要显式列出本质边界进解空间自然边界进积分表达式数值策略需要差分逼近二阶导数只对一阶导数积分分片多项式可直接使用物理含义每个点逐点平衡加权整体平衡“弱”不代表计算精度差。它放宽的是逐点要求换来的是数值可操作性。你选一个分片线性的 v_h代入弱形式得到的就是关于节点值的代数方程组。理解这个层次后那些带下标的函数空间定义就不再是障碍而只是说明了“允许哪些函数进来”。2.4 用具体函数验证分部积分边界项弱形式推导容易在边界项符号上出错。常见做法是用一组具体函数核对恒等式在 MATLAB 里跑通后再写正式程序% 用具体函数验证分部积分恒等式 % ∫ -u v dx ∫ u v dx - [uv]_0^1 syms x u x^2; % 测试函数 u v x * (1 - x); % 满足 v(0)0, v(1)0 的测试函数 L 1; I1 int(-diff(u, x, 2) * v, x, 0, L); % 原始弱形式左端 I2 int(diff(u, x) * diff(v, x), x, 0, L) - ... (subs(diff(u, x) * v, x, L) - subs(diff(u, x) * v, x, 0)); simplify(I1 - I2) % 输出 0 说明分部积分推导正确逻辑说明I1 是原方程左端乘测试函数后的积分I2 是分部积分后去掉边界项的结果。对满足 v(0)v(1)0 的测试函数边界项subs(diff(u,x)*v, x, L)和subs(diff(u,x)*v, x, 0)都应该为 0但代码里保留完整形式是为了验证符号运算本身没出错。把 u、v 换成分段函数也能做同样验证只是符号积分可能变慢。这个习惯值得保留任何弱形式推导先用具体函数验边界项再进离散。3. 有限元离散线性基函数下的刚度矩阵与载荷向量3.1 从弱形式到线性方程组弱形式仍然是无限维问题因为函数空间 V 中有无穷多个函数。有限元离散做的事情是把 V 换成一个有限维子空间 V_h比如所有在节点 0 到 N 上取值、在单元内部是一次多项式的连续分片线性函数。这个空间里的任意函数完全由 N1 个节点值 U₀, U₁, ..., U_N 决定。用形状函数 φ_i(x) 描述第 i 个节点上的“尖塔”φ_i(x_i)1在相邻两个单元内线性降到 0其他位置全为 0。V_h 中任意函数都能写成 u_h(x) Σ U_j φ_j(x)代入弱形式并依次取 v φ_i得到 N1 个代数方程K_ij ∫ k φ_i φ_j dxF_i ∫ f φ_i dx g φ_i(1)这就是 K U F 的来源。刚度矩阵 K 是稀疏的因为 φ_i 和 φ_j 的支撑区间只在 i、j 相同或相邻时相交每个单元只对相邻自由度有贡献。这个局部性决定了有限元组装的“单元循环”方案逐个单元计算局部矩阵再按节点编号投放进全局矩阵。3.2 单元刚度矩阵为什么要记住 1 -1 -1 1在单个单元 e 上单元长度为 h局部节点 1 和 2 对应全局节点 e 和 e1。两个线性形状函数的导数为 φ₁ -1/hφ₂ 1/h代入 K_ij 的积分k 取常数时得到局部刚度矩阵ke (k/h) * [[1, -1], [-1, 1]]这个 2×2 矩阵是整个一维有限元最常用的公式。物理意义可以从两个角度看每行元素之和为零对应刚体平移模式下单元内没有应力对角线元素为正保证系统正定。组装时局部节点 1 的贡献累加到全局第 e 行、第 e 列局部节点 2 的贡献累加到第 e1 行、第 e1 列。载荷向量的单元贡献取决于 f 的积分。最简单的近似是两点梯形积分fe h/2 * (f(x_e) f(x_{e1}))。f 是线性函数时这个公式精确f 是非线性函数时它是一阶近似。先把梯形积分跑通框架再升级成高斯积分是更稳的推进节奏。3.3 matlab有限元编程求解实例组装与求解的最小代码下面是一段完整的一维 matlab有限元编程求解实例代码对应方程 -u π² sin(πx)边界条件 u(0)u(1)0解析解 u sin(πx)。k 取 1g 取 0。% 一维泊松方程 -u f, u(0)u(1)0 % 解析解 u sin(pi*x) f (x) pi^2 .* sin(pi * x); N 20; % 单元数 x linspace(0, 1, N 1); % 节点坐标 h 1 / N; % 单元长度 K sparse(N 1, N 1); % 全局刚度矩阵 F zeros(N 1, 1); % 全局载荷向量 for e 1:N n [e, e 1]; % 当前单元对应的全局节点 ke (1 / h) * [1, -1; -1, 1]; % 单元刚度矩阵k1 fe h / 2 * [f(x(e)); f(x(e 1))]; % 单元载荷向量梯形积分 K(n, n) K(n, n) ke; % 投放进全局矩阵 F(n) F(n) fe; % 投放进全局向量 end % 本质边界条件 u(0)0, u(1)0只求解内部自由节点 free 2:N; U zeros(N 1, 1); U(free) K(free, free) \ F(free); % 与解析解对比前 4 个节点 ue sin(pi * x); disp([x(1:4), U(1:4), ue(1:4)]);逻辑说明每个单元独立计算局部刚度矩阵和局部载荷向量再依据全局节点编号把贡献“累加”进 K 和 F。sparse在 N 增大时节省大量内存避免全稠密矩阵。free 2:N排除左端节点 1 和右端节点 N1因为这两个节点值已知为零只对内部节点求解解出来后与端节点上的零合并成完整 U。参数说明N 控制网格密度N 增大时近似解更接近真解但二维问题自由度会平方增长f 是右端热源项以函数句柄传入可以让代码不绑死具体表达式h 是网格尺寸初始化后先打印确认等于 1/N能排查网格划分错误。U(free) K(free, free) \ F(free)用的是 MATLAB 内置稀疏直接求解器一维问题 N 在 10 万以内都很快。3.4 载荷向量升级两点高斯积分梯形积分对光滑但非线性的 f 会有可感知的误差。常见做法是换成两点高斯-勒让德积分积分点取 ±1/√3权重取 1。把参考区间 [-1, 1] 映射到物理单元 [x_e, x_{e1}] 后单元中点为 x_m (x_e x_{e1})/2高斯点位于 x_m ± h/(2√3)。替换前面代码里的 fe 行xm (x(e) x(e 1)) / 2; s h / (2 * sqrt(3)); fe h / 2 * [f(xm - s); f(xm s)];逻辑说明两点高斯积分对三次多项式以内的 f 精确梯形积分只对一次多项式精确。高斯点上的 f 值经过雅可比因子 h/2 映射到物理单元。对大多数光滑右端项两点高斯足够如果 f 含强烈变化先加密网格再考虑增加高斯点数量而不是只靠单元内加密。4. 弱形式中的边界条件Dirichlet与Neumann的两张面孔4.1 本质边界条件与自然边界条件的自由度差异回到弱形式边界条件分两类一类直接给定 u 的值叫本质边界条件或 Dirichlet 条件一类给定导数组合叫自然边界条件或 Neumann 条件。本质边界条件必须进入解空间u_h 在对应节点上取给定值测试函数在对应节点上取 0。自然边界条件不需要改解空间它通过弱形式中的边界积分直接进入右端项。数值上经常出错的点有两个一是在 Neumann 边界上错误地强制节点为零二是在 Dirichlet 边界上忘记做任何处理。判断方法很简单看弱形式边界积分项里 u 是否已知。u 已知是本质边界要改矩阵u 未知但它的导数组合已知是自然边界只加右端项。4.2 直接消去法把已知节点值移到右边非齐次 Dirichlet 边界 u(0)a、u(1)b 是最常见的情况。直接消去法的思路是先把已知节点值写入解向量再把已知节点对未知节点方程的影响移到右端最后只对自由节点求解。bnodes [1; N 1]; % 本质边界对应的全局节点编号 bvals [a; b]; % 边界节点值 U zeros(N 1, 1); U(bnodes) bvals; % 写入已知节点值 F2 F - K * U; % 已知节点贡献移到右端 free setdiff(1:N 1, bnodes); U(free) K(free, free) \ F2(free);逻辑说明完整方程是 K U F但 U 中部分分量已知。把 K 与已知 U 相乘得到的是边界节点对自由节点方程贡献的“力”移到右端得到 F2。如果只取 K(free, free) 和 F(free)边界节点影响被完全丢掉结果会明显偏离真解。setdiff适合教学和小规模问题大规模运算建议在组装前完成自由度编号避免每次构造 free 索引。直接消去法的优点是精确边界值以机器精度进入解。缺点是维护自由度索引稍微麻烦。自适应加密下节点编号会变化边界节点和自由节点必须分开记录。4.3 罚函数法一个参数决定收敛性罚函数法把 Dirichlet 条件作为强惩罚项放进弱形式相当于在边界上加了很大的弹簧。离散后的修改非常简单C 1e8; % 惩罚系数经验值取刚度矩阵最大对角元的 1e6 到 1e8 倍 for i 1:numel(bnodes) K(bnodes(i), bnodes(i)) K(bnodes(i), bnodes(i)) C; F(bnodes(i)) F(bnodes(i)) C * bvals(i); end U K \ F;惩罚系数 C 选小了边界约束不足表现为边界节点值偏离给定值选大了矩阵条件数迅速恶化迭代求解器可能不收敛。经验准则是让 C 与刚度矩阵对角最大元素保持至少 6 个数量级差距然后检查abs(U(bnodes) - bvals)是否在可接受范围。实际工程里直接消去法更可控罚函数法适合快速原型实现。对比项直接消去法罚函数法实现复杂度需要维护自由节点索引只加对角项实现最简单边界满足程度机器精度取决于惩罚系数条件数影响无额外恶化系数过大会恶化适用场景中小规模、精度要求高快速原型、非结构化网格4.4 自然边界条件往右端项加值的顺序不能错右端 Neumann 边界 k u(1)g弱形式载荷向量最后一个分量需要加 gF(end) F(end) g。原因是最末端节点的形状函数 φ_N 在该点取值为 1边界积分 g v(1) 在测试函数取 φ_N 时直接落到最后一个自由度上。有一个顺序问题容易忽略如果同一节点上同时存在 Dirichlet 和 Neumann 条件必须先加自然边界贡献再做本质边界消去。先消去本质边界会把 Neumann 贡献错误地带进自由节点方程结果整个右端项都会偏。推荐把“先自然、后本质”作为任何有限元程序里的固定顺序写进注释里。5. 用误差范数验证弱形式实现MATLAB算例与网格收敛5.1 构造有解析解的光滑问题验证弱形式程序是否写对的可靠方式是设计一个已知解析解的问题。取 u(x) sin(πx)则 -u π² sin(πx)。在代码里设 f (x) pi^2 .* sin(pi*x)k1边界条件 u(0)u(1)0。这个解足够光滑没有角点奇异性干扰误差只来自网格离散和数值积分方便观察收敛阶。5.2 用单元高斯积分计算误差范数误差计算不要用 interp1 对数值解重新采样那样会混入插值误差。更干净的做法是在每个单元上用两点高斯积分直接计算精确解和数值解之差err_L2 0; err_grad 0; for e 1:N xm (x(e) x(e 1)) / 2; s h / (2 * sqrt(3)); xg [xm - s, xm s]; % 高斯点 ug [U(e), U(e 1)]; % 数值解在高斯点的线性插值 ue sin(pi * xg); % 精确解 ue_x pi * cos(pi * xg); % 精确解的导数 grad_uh (U(e 1) - U(e)) / h; % 单元内数值解导数为常数 err_L2 err_L2 h / 2 * sum((ug - ue).^2); err_grad err_grad h / 2 * sum((grad_uh - ue_x).^2); end err_L2 sqrt(err_L2); err_grad sqrt(err_grad); err_H1 sqrt(err_L2^2 err_grad^2);逻辑说明高斯点上的数值解由该单元两个节点值线性插值得到精确解直接用解析公式。误差平方乘上雅可比因子 h/2 后累加最后开方得到 L2 范数。梯度误差部分利用了线性元导数为常数的特性数值导数不需要额外差分。线性元的理论收敛阶为 L2 误差约 O(h²)H1 误差约 O(h)。加密网格时可以按下面表格核对范数线性元预期收敛阶网格加密一倍后的现象L2 误差约 2误差减为原来的约 1/4H1 误差约 1误差减为原来的约 1/25.3 程序能跑通的三个快速检查第一个检查刚度矩阵对称性norm(K - K, fro)应该为 0。组装索引写错时经常会丢失非对角线的一半贡献这个语句能立刻暴露问题。第二个检查刚体模式在纯 Neumann 问题里把 U 全部置 1计算K * ones(N1,1)应该接近 0带 Dirichlet 边界时内部行的结果仍接近 0边界行会有响应。第三个检查求解残差res norm(F - K * U)如果消去边界条件时只是把行清零而不是缩减矩阵残差里会出现边界行的巨大数值。N 超过 10 万后一维三对角刚度矩阵的直接求解仍然很快但推广到二维三角形网格时带宽增大内存开销会明显上升。常见做法是切换成共轭梯度法加不完全 Cholesky 预条件pcg(K(free,free), F2(free), 1e-8, 400, ichol(K(free,free)))。此时如果还在用罚函数法处理 Dirichlet 边界巨大对角线会破坏预条件子收敛会非常慢。弱形式理论看到这里你会明白直接消去法虽然代码多一点却是最可预测的方案。把 K 和 F 的维度与网格节点数对齐打印出来能快速定位组装漏单元的常见问题下一步把一维代码推广到二维三角形单元时弱形式本身不动变的只是单元矩阵和积分规则。本文还有配套的精品资源点击获取