当前位置：首页 > news >正文

【数值分析】有限元法（FEM）数值不稳定性分析：$L^2$ 与 $H_0^1$ 空间的对比实验

news 2026/6/6 10:27:31

科研入门：从数学系本科生的视角，拆解有限元法背后的力学逻辑

总起

相信每个搞数值模拟的同学，都有过类似的心理创伤：代码里模型搭得好好的，咱们吃着火锅唱着歌，正期待着漂亮的云图呢。结果突然一个网格畸形，局部应力瞬时爆表——紧接着就像多米诺骨牌一样，误差传遍全局，整个仿真瞬间炸成一团乱码。这种‘闹心’的底层原因，往往在于我们离散化的数学模型在物理和数学契合度上出了岔子。

本文就从函数空间的视角切入，聊聊为什么在H01H_0^1H01空间中建模能为有限元提供天然的稳定性，以及它是如何通过对‘能量’的刻画，让模拟不再轻易‘原地爆炸’的。”

学习路线：

理解连续介质力学中的基本对象
学习位移梯度和应变、应力
在L2L^2L2和H01H_0^1H01中建模，搞清楚为什么L2L^2L2会炸
一个直观的数值例子

1. 连续介质力学中的基本对象：位移场uuu

在开始之前，要先从高中的思维转变过来。高中物理里，一个质点就是一个点，给它一个力，它就动。但连续介质力学研究的对象是由无数个质点组成的连续体，比如一块钢板、一团橡皮泥。我们关心的不是某一个质点，而是整个连续体里每一个点在受力后跑到哪里去了。

这就引出了位移场的概念：

u(x,y,z)=[ux(x,y,z)uy(x,y,z)uz(x,y,z)]u(x,y,z) = \begin{bmatrix} u_x(x,y,z) \\ u_y(x,y,z) \\ u_z(x,y,z) \end{bmatrix}u(x,y,z)=ux(x,y,z)uy(x,y,z)uz(x,y,z)

输入是连续体里某个质点的坐标(x,y,z)(x, y, z)(x,y,z)，输出是这个质点在受力后的位移量。注意这里uuu是一个向量场，连续体里每一个点都对应一个位移向量。

扩展

知道位移就能算应变，知道应变通过本构关系能算应力
知道应力就能写平衡方程∇⋅σ+b=0\nabla \cdot \sigma + b = 0∇⋅σ+b=0
这就建立了偏微分方程，接下来要么在 Sobolev 空间里找弱解，要么用有限元数值求解

2. 从位移到应力：位移梯度、应变、本构关系

2.1 位移梯度

有了位移场uuu，下一步是对它求导，得到位移梯度：

∇u=∂u∂x=[∂ux∂x∂ux∂y∂ux∂z∂uy∂x∂uy∂y∂uy∂z∂uz∂x∂uz∂y∂uz∂z]\nabla \boldsymbol{u} = \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x}} = \begin{bmatrix} \frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} & \frac{\partial u_x}{\partial z} \\ \frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} & \frac{\partial u_y}{\partial z} \\ \frac{\partial u_z}{\partial x} & \frac{\partial u_z}{\partial y} & \frac{\partial u_z}{\partial z} \end{bmatrix}∇u=∂x∂u=∂x∂ux∂x∂uy∂x∂uz∂y∂ux∂y∂uy∂y∂uz∂z∂ux∂z∂uy∂z∂uz

这是一个3×33\times33×3的矩阵，描述的是连续体在某点附近局部的变形情况。

2.2 应变张量

取位移梯度的对称部分，得到柯西应变张量：

ε(u)=12(∇u+∇u⊤)\varepsilon(u) = \frac{1}{2}(\nabla u + \nabla u^\top)ε(u)=21(∇u+∇u⊤)

取对称部分是因为我们只关心纯变形，不关心刚体转动——反对称部分描述的是转动，不产生内力。

扩展

这里的ε(u)\varepsilon(u)ε(u)是小变形假设下的线性化应变，适用于柯西应力框架
处理大变形时，线性化失效，需要引入变形梯度张量F=I+∇u\mathbf{F} = \mathbf{I} + \nabla uF=I+∇u，并使用格林-拉格朗日应变张量

2.3 本构关系与应力

有了应变，通过本构关系推导出应力：

σ(u)=C:ε(u)\sigma(u) = \mathbb{C} : \varepsilon(u)σ(u)=C:ε(u)

C\mathbb{C}C是四阶弹性张量，就是写代码时那个包含杨氏模量EEE、泊松比ν\nuν的矩阵
:::是双点积，表示两个张量在对应维度上缩并求和

整个链条就是：u→∇u→ε(u)→σ(u)→∇⋅σ+b=0u \rightarrow \nabla u \rightarrow \varepsilon(u) \rightarrow \sigma(u) \rightarrow \nabla \cdot \sigma + b = 0u→∇u→ε(u)→σ(u)→∇⋅σ+b=0。

注意这里有个关键依赖：应力σ\sigmaσ的存在，依赖于∇u\nabla u∇u的存在。这个依赖关系，正是接下来问题的根源。

3. 在哪个函数空间建模？L2L^2L2会炸，H01H_0^1H01才稳

3.1 为什么L2L^2L2不行

L2L^2L2空间的定义只有一个要求：函数本身平方可积，即

∫∣u∣2 dx<∞\int |u|^2 \, dx < \infty∫∣u∣2dx<∞

听起来挺合理。但问题在于，L2L^2L2对导数没有任何约束。一个L2L^2L2函数完全可以在某些点上有奇异性，弱导数根本不存在。

回想一下上面的本构链：σ∼C:ε∼∇u\sigma \sim \mathbb{C} : \varepsilon \sim \nabla uσ∼C:ε∼∇u。如果∇u\nabla u∇u不存在或者不可积，那么：

应变ε\varepsilonε没有意义
应力σ\sigmaσ没有意义
系统总应变能∫σ:ε dx\int \sigma : \varepsilon \, dx∫σ:εdx直接发散

这就是为什么L2L^2L2方案在网格加密时应力会一直涨——函数空间根本没有保证导数是可控的，网格越细，奇异性就越被"看见"，算出来的应力就越大。

3.2 为什么H01H_0^1H01更稳

H1H^1H1是一个 Sobolev 空间，它的要求比L2L^2L2多了一条：一阶弱导数也必须平方可积，即

u∈L2且∇u∈L2u \in L^2 \quad \text{且} \quad \nabla u \in L^2u∈L2且∇u∈L2

下标000（即H01H_0^1H01）表示在边界上施加齐次狄利克雷条件，也就是物体边界被固定、位移为零的物理情形。

强制要求u∈H01u \in H_0^1u∈H01，就自动保证了：

∇u∈L2\nabla u \in L^2∇u∈L2，从而应变ε∈L2\varepsilon \in L^2ε∈L2，应力σ∈L2\sigma \in L^2σ∈L2
系统总应变能∫σ:ε dx\int \sigma : \varepsilon \, dx∫σ:εdx有限，有物理意义
由Lax-Milgram 定理保证弱解存在且唯一
有限元刚度矩阵条件数可控，误差分析收敛

4. 一个直观的数值例子

光说不练假把式，来看数值实验。

完整代码：https://colab.research.google.com/drive/1u8I2STixmTuEeof99jcgqoiuuT5gHxju#scrollTo=MvlRX1c-M4XH

对同一个弹性体问题，分别用L2L^2L2和H01H_0^1H01近似，随着网格加密（NNN增大，网格尺寸hhh减小），看最大应力的变化：

============================================================ 对比：L2 与 H01 近似的最大应力 ============================================================ N h L2最大应力 H01最大应力 20 0.05000 4.922136e+00 4.750000e-01 40 0.02500 6.799555e+00 4.875000e-01 80 0.01250 9.431772e+00 4.937500e-01 160 0.00625 1.314286e+01 4.968750e-01

结果很直接：