2 控制方程

2.1 流体与固体

日常生活中我们会遇到四种物质状态——固体、液体、气体与等离子体——它们乍看之下由对约束的响应来区分：固体既保持形状又保持体积，液体保持体积但会顺应容器的形状，气体则会膨胀充满一切可用的空间。然而对于连续介质力学而言，真正具有决定性意义的区别并非几何的，而是动力学的，它关乎物质对剪切应力的响应。

流体是一种在静止时无法承受剪切应力的物质：只要施加一个剪切力，无论它多么微小，流体都会持续地变形——它会流动——只要该力持续作用。相反，固体在足够小的剪切载荷下会以有限的静态变形作出响应，并稳定到一个新的平衡构型。因此，静止的流体只传递法向应力，即压强；剪切应力仅在运动中产生，并通过粘性得以抵抗。于是，流体一词同时涵盖了液体与气体：尽管二者在密度和可压缩性上差异极大，但它们都会在任意小的剪切作用下流动，并且在宏观层面上都服从本章将要推导的同一组平衡律。关于其背后的物理，我们参考 Batchelor [1] 以及 Landau 和 Lifshitz [23] 的经典论述。

物质状态之间的宏观差异，根源在于分子之间的作用力。两个中性分子之间的作用力 $F$ （作为其间距 $r$ 的函数，示意见图 2.1）在短程强烈排斥，在长程微弱吸引，长程处按 $F \sim - k / r^{7}$ 衰减，其中常数 $k$ 由分子本身决定；两种贡献在平衡间距 $d_{0}$ 处相互平衡。对于像水这样的极性分子，其正负电荷中心不重合，分子带有一个永久偶极矩，相应地其吸引尾部也更强。

图 2.1：两个原子之间作用力随其间距变化的示意图。

这些分子间作用力的相对强弱决定了分子活动性的程度，并进而决定了物质状态。在固体中，分子被紧密束缚，只能围绕固定、有序的位置振动；在液体中，分子仍保持紧密接触，但可以相互滑移，只保留短程有序；在气体中，分子几乎自由且无序。表 2.1 总结了这一图像。

表 2.1：固体与流体的分子性质

	分子间作用力	分子排列	统计类型
固体	强	有序	量子
液体	中	部分有序	量子 + 经典
气体	弱	无序	经典

尽管物质在分子尺度上是离散的，流体力学却把流体描述为一种连续介质。所关心的场量——密度、速度、压强与温度——都被视为位置与时间的光滑函数，它们在每一点上被定义为在某个控制体上的平均值，该控制体相对于分子间距足够大，但相对于流动的宏观特征长度 $L$ 又足够小。这种尺度上的分离即为连续介质假设。只要分子碰撞之间的平均自由程 $ℓ$ 相对于 $L$ 可以忽略，即当 Knudsen 数 $Kn = ℓ / L$ 很小时，该假设就成立；它在稀薄气体以及被约束到分子尺度的流动中失效 [1, 11]。同样的假设也构成了可变形固体力学的基础，这两个领域共同构成了连续介质力学。

最后，固体行为与流体行为之间的界限并非泾渭分明。许多材料——稠密悬浮液、聚合物溶液与熔体、膏体与凝胶——在某些条件下像液体一样流动，在另一些条件下又像固体一样抵抗，这取决于所施加载荷的大小与时间尺度。这类非牛顿行为的系统分类与建模是第 3 章的主题；而本章则发展任何此类材料（无论其具体本构响应如何）都遵守的运动学与平衡律。

2.2 运动的 Lagrange 描述与 Euler 描述

在采用了 2.1 节的连续介质假设之后，我们现在来描述流体的运动。与固体不同，流体可以经历任意大的变形，这使得描述方式的选择变得重要。最直接的方法是逐个跟踪每一个物质粒子：用其在初始时刻 $t = 0$ 的参考构型中的位置 $X$ 来标记一个粒子，则其运动由映射 $ϕ$ 给出，它把该粒子映射到当前位置 $x$ ，

\begin{matrix} (2.1) & x = ϕ (X, t) . \end{matrix}

假设 2.1

对每个 $t$ ， $ϕ (\cdot, t)$ 是一个微分同胚，即 $ϕ (\cdot, t) \in C^{1}$ 且为双射，其逆 $ϕ (\cdot, t)^{- 1} \in C^{1}$ ；
映射 $ϕ (\cdot, t)$ 对每个 $t$ 可逆，逆为 $ϕ^{- 1} (\cdot, t)$ ；
$ϕ (\cdot, 0) = Id$ 。

图 2.2 展示了所考察流体某一给定体积的一般变形。在有限变形分析中，对于位移 $u := x - X$ 的大小不作任何假设。在无穷小变形分析中，则假设位移 $u$ 比所考察流体系统的特征尺寸小若干个数量级。

图 2.2：可变形体积的一般运动。

为了描述运动的物质描述与空间描述之间的区别，我们考察一个标量，例如流体的温度 $T$ 。

物质描述： 温度 $T$ 相对于 $t = 0$ 时刻坐标为 $X$ 的物质粒子来描述，即

\begin{matrix} (2.2) & T = \hat{T} (X, t) . \end{matrix}

空间描述： 温度 $T$ 相对于空间中任一时刻 $t$ 的固定位置 $x$ 来描述，即

\begin{matrix} (2.3) & T = T (x, t) . \end{matrix}

此处及下文中，我们用符号 $\hat{\cdot}$ 标记那些相对于物质粒子表达的量。

物质描述，即 (2.2) 中给出的 Lagrange 描述，关注给定流体物质粒子温度 $T$ 的变化。而 (2.3) 表示的空间描述，也称 Euler 描述，其侧重点在于给定位置 $x$ 处温度的变化；无论是哪个物质粒子正经过位置 $x$ 。对于流体而言，Lagrange 描述显然极为繁琐，因为要理解整体运动就必须跟踪众多物质粒子的运动。这正是流体力学工作者几乎只使用 Euler 描述的原因。

物质时间导数把从 Lagrange 视角考察的物理量的时间变化率，与从 Euler 视角观测到的同一物理量的变化率联系起来。考察由此前定义的方程给出的流体运动

\begin{matrix} (2.4) & x = ϕ (X, t), \end{matrix}

由此，粒子的速度，即遵循 Lagrange 描述的速度，被定义为 $ϕ$ 的时间导数

\begin{matrix} (2.5) & \hat{v} (X, t) = \frac{\partial ϕ (X, t)}{\partial t} . \end{matrix}

尽管这里我们考察的是运动的 Lagrange 描述，速度仍是一个空间矢量。通过对方程求逆，速度可以表达为与 Euler 描述相关联的空间位置 $x$ 的函数，

\begin{matrix} (2.6) & v (x, t) = \hat{v} (ϕ^{- 1} (x, t), t) . \end{matrix}

空间中固定点 $x$ 处一般标量或张量 $T (x, t)$ 的变化率，由如下关于时间的偏导数表示：

\begin{matrix} (2.7) & \frac{\partial T (x, t)}{\partial t} = lim_{δ t \to 0} \frac{T (x, t + δ t) - T (x, t)}{δ t} . \end{matrix}

由一个以速度 $u$ 运动的观测者所测得的同一物理量在 $x$ 处的变化率，则由下式给出

\begin{matrix} (2.8) & lim_{δ t \to 0} \frac{T (x + u δ t, t + δ t) - T (x, t)}{δ t} = \frac{\partial T (x, t)}{\partial t} + u \cdot \nabla T (x, t) . \end{matrix}

若取 $u$ 为空间位置 $x$ 处、时刻 $t$ 的粒子速度，即 $v (x, t)$ ，则该变化率对应于 $\frac{d T (x, t)}{d t}$ ，即所谓的物质时间导数

\begin{matrix} (2.9) & \frac{d T (x, t)}{d t} = \frac{\partial T (x, t)}{\partial t} + \underset{对流部分}{\underset{⏟}{v (x, t) \cdot \nabla T (x, t)}} . \end{matrix}

表达式 $v (x, t) \cdot \nabla T (x, t)$ 对应于物质时间导数的对流部分。

在方程 (2.9) 中把标量 $T$ 替换为速度 $v (x, t)$ 的各个分量，我们得到加速度矢量 $a (x, t)$ ，

\begin{matrix} (2.10) & a (x, t) = \frac{\partial v (x, t)}{\partial t} + (D v (x, t)) v (x, t) . \end{matrix}

对某些应用而言，下面这个与 (2.10) 等价的加速度矢量表达式会很有用

\begin{matrix} (2.11) & a (x, t) = \frac{\partial v (x, t)}{\partial t} + \nabla \frac{| v (x, t) |^{2}}{2} + (\nabla \times v (x, t)) \times v (x, t) . \end{matrix}

在方程 (2.10) 中， $L := D v (x, t)$ 是称为速度梯度张量的二阶张量。

设采用笛卡尔基，并令 $v (x, t) = (v_{1} (x, t), v_{2} (x, t), v_{3} (x, t))^{T}$ ，我们对 $T$ 与 $v$ 的梯度、散度、旋度与 Jacobi 矩阵采用如下记号，其中 $χ \in {x, X}$ ， $d \in {2, 3}$ ：

\begin{matrix} (2.12) & \nabla_{χ} T (χ) := (\partial_{χ_{1}} T, \dots, \partial_{χ_{d}} T) \in R^{1 \times d} \end{matrix}

\begin{matrix} (2.13) & \nabla_{χ} \times v (χ) := \partial_{χ_{1}} v_{2} - \partial_{χ_{2}} v_{1} \in R 当 d = 2 \end{matrix}

\begin{matrix} (2.14) & \nabla_{χ} \cdot v (χ) := \partial_{χ_{1}} v_{1} + \dots + \partial_{χ_{d}} v_{d} \in R \end{matrix}

\begin{matrix} (2.15) & \nabla_{χ} \times v (χ) := (\begin{matrix} \partial_{χ_{2}} v_{3} - \partial_{χ_{3}} v_{2} \\ \partial_{χ_{3}} v_{1} - \partial_{χ_{1}} v_{3} \\ \partial_{χ_{1}} v_{2} - \partial_{χ_{2}} v_{1} \end{matrix}) \in R^{d \times 1} 当 d = 3 \end{matrix}

\begin{matrix} (2.16) & D_{χ} v (χ) := (\begin{matrix} \partial_{χ_{1}} v_{1} & \dots & \partial_{χ_{d}} v_{1} \\ ⋮ & \dots & ⋮ \\ \partial_{χ_{1}} v_{d} & \dots & \partial_{χ_{d}} v_{d} \end{matrix}) \in R^{d \times d} . \end{matrix}

于是，三维中的速度梯度张量为

\begin{matrix} (2.17) & L = D v = D_{x} v = (\begin{matrix} \frac{\partial v_{1}}{\partial x_{1}} & \frac{\partial v_{1}}{\partial x_{2}} & \frac{\partial v_{1}}{\partial x_{3}} \\ \frac{\partial v_{2}}{\partial x_{1}} & \frac{\partial v_{2}}{\partial x_{2}} & \frac{\partial v_{2}}{\partial x_{3}} \\ \frac{\partial v_{3}}{\partial x_{1}} & \frac{\partial v_{3}}{\partial x_{2}} & \frac{\partial v_{3}}{\partial x_{3}} \end{matrix}) . \end{matrix}

在不引起歧义时，我们在 Euler 情形下略去微分算子的下标 $x$ ，即写 $\nabla$ 而非 $\nabla_{x}$ 。速度梯度张量之于速度，正如变形梯度张量之于位移（见 2.3 节）。速度梯度张量可以分解为对称部分与反对称部分之和：

\begin{matrix} (2.18) & L = D + W, \end{matrix}

其中

\begin{matrix} (2.19) & D = \frac{1}{2} (L + L^{T}), \end{matrix}

\begin{matrix} (2.20) & W = \frac{1}{2} (L - L^{T}) . \end{matrix}

对称张量 $D$ 称为变形率张量。设采用笛卡尔基，变形率张量 $D$ 的各分量为

\begin{matrix} (2.21) & D = (\begin{matrix} \frac{\partial v_{1}}{\partial x_{1}} & \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{2}} + \frac{\partial v_{2}}{\partial x_{1}}) & \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{3}} + \frac{\partial v_{3}}{\partial x_{1}}) \\ \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{2}} + \frac{\partial v_{2}}{\partial x_{1}}) & \frac{\partial v_{2}}{\partial x_{2}} & \frac{1}{2} (\frac{\partial v_{2}}{\partial x_{3}} + \frac{\partial v_{3}}{\partial x_{2}}) \\ \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{3}} + \frac{\partial v_{3}}{\partial x_{1}}) & \frac{1}{2} (\frac{\partial v_{2}}{\partial x_{3}} + \frac{\partial v_{3}}{\partial x_{2}}) & \frac{\partial v_{3}}{\partial x_{3}} \end{matrix}) . \end{matrix}

$D$ 的对角分量之和对应于速度的散度 $div (v)$ 。 $D$ 的对角元描述沿基方向的拉伸速度，而非对角元则与相应基方向线元之间夹角的变化率有关。关于 $D$ 各分量的详细解释，参见 [1]。

反对称张量 $W$ 称为自旋张量，它描述转动速度。设采用笛卡尔基， $W$ 的各分量为

\begin{matrix} (2.22) & W = (\begin{matrix} 0 & - \frac{1}{2} (\frac{\partial v_{2}}{\partial x_{1}} - \frac{\partial v_{1}}{\partial x_{2}}) & \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{3}} - \frac{\partial v_{3}}{\partial x_{1}}) \\ \frac{1}{2} (\frac{\partial v_{2}}{\partial x_{1}} - \frac{\partial v_{1}}{\partial x_{2}}) & 0 & - \frac{1}{2} (\frac{\partial v_{3}}{\partial x_{2}} - \frac{\partial v_{2}}{\partial x_{3}}) \\ - \frac{1}{2} (\frac{\partial v_{1}}{\partial x_{3}} - \frac{\partial v_{3}}{\partial x_{1}}) & \frac{1}{2} (\frac{\partial v_{3}}{\partial x_{2}} - \frac{\partial v_{2}}{\partial x_{3}}) & 0 \end{matrix}) . \end{matrix}

$W$ 的张量分量与涡量矢量 $ω$ 直接相关

\begin{matrix} (2.23) & ω = \nabla \times v \end{matrix}

\begin{matrix} (2.24) & = {(\frac{\partial v_{3}}{\partial x_{2}} - \frac{\partial v_{2}}{\partial x_{3}}, \frac{\partial v_{1}}{\partial x_{3}} - \frac{\partial v_{3}}{\partial x_{1}}, \frac{\partial v_{2}}{\partial x_{1}} - \frac{\partial v_{1}}{\partial x_{2}})}^{T} . \end{matrix}

对任意矢量 $n$ ，乘积 $W n$ 对应于 $\frac{ω}{2}$ 与 $n$ 的叉积，即

\begin{matrix} (2.25) & W n = \frac{1}{2} ω \times n . \end{matrix}

矢量 $\frac{ω}{2}$ 描述所考察物质流体微元的转动速度。 $W$ 显然与变形率完全无关。

2.3 变形梯度与应变张量

为了从物理上理解变形梯度张量，我们考察两个物质粒子 $P$ 与 $Q$ 。 $P$ 与 $Q$ 在变形前 $t = 0$ 时刻的位置由矢量 $X_{P}$ 与 $X_{Q}$ 给出。变形后，物质粒子 $P$ 与 $Q$ 将处于如下空间位置

\begin{matrix} (2.26) & x_{P} = ϕ (X_{P}, t); 以及 x_{Q} = ϕ (X_{Q}, t) . \end{matrix}

图 2.3：变形梯度张量定义的示意图。此处假设为无穷小量。

由图 2.3 可见，距离 $d X = X_{Q} - X_{P}$ 在 $t$ 时刻变形后将变为 $d x$ ，即

\begin{matrix} (2.27) & d x = x_{Q} - x_{P}, \end{matrix}

\begin{matrix} (2.28) & = ϕ (X_{P} + d X, t) - ϕ (X_{P}, t), \end{matrix}

\begin{matrix} (2.29) & = ϕ (X_{P}, t) + {\frac{\partial ϕ (X, t)}{\partial X} |}_{X = X_{P}} d X - ϕ (X_{P}, t) + O (| d X |^{2}), \end{matrix}

\begin{matrix} (2.30) & = \hat{F} (X_{P}, t) d X, \end{matrix}

其中

\begin{matrix} (2.31) & \hat{F} (X_{P}, t) = {\frac{\partial ϕ (X, t)}{\partial X} |}_{X = X_{P}}, \end{matrix}

描述了点 $X_{P}$ 、时刻 $t$ 处的变形梯度张量。为了得到等式 (2.30)，我们假设线段 $d X$ 为无穷小量。空间元矢量 $d x$ 也称为物质矢量 $d X$ 的前推（push forward）等价物。反过来，物质矢量称为空间矢量 $d x$ 的拉回（pull back）等价物，其含义为

\begin{matrix} (2.32) & d x = \hat{F} (X_{P}, t) d X . \end{matrix}

在本书后续内容中，我们将避免显式写出变形张量的自变量，采用记号 $\hat{F}$ 。关于变形梯度张量性质的更详细描述，例如可参考 [6, 20, 22]。应变定义为无穷小线段在变形下的相对长度变化。基于变形梯度张量的定义，我们可以把长度 $∥ d x ∥$ 表达为如下形式：

\begin{matrix} (2.33) & ∥ d x ∥^{2} = (\hat{F} (X_{P}, t) d X, \hat{F} (X_{P}, t) d X) \end{matrix}

\begin{matrix} (2.34) & = d X^{T} {\hat{F}}^{T} (X_{P}, t) \hat{F} (X_{P}, t) d X \end{matrix}

\begin{matrix} (2.35) & = d X^{T} \hat{C} (X_{P}, t) d X, \end{matrix}

其中

\begin{matrix} (2.36) & \hat{C} = {\hat{F}}^{T} \hat{F}, \end{matrix}

表示右 Cauchy–Green 张量，也称为 Green 变形张量。称之为右张量，是因为在 $\hat{C}$ 的乘积定义中张量 $\hat{F}$ 位于右侧。该张量显然是对称且正定的，因为

\begin{matrix} (2.37) & (\hat{C} d X, d X) = ({\hat{F}}^{T} \hat{F} d X, d X) = (\hat{F} d X, \hat{F} d X) = ∥ \hat{F} d X ∥^{2} > 0 \forall d X \neq 0. \end{matrix}

为了定义左 Cauchy–Green 张量，我们引入 $\hat{F}$ 的 Euler 对应物

\begin{matrix} (2.38) & F (x, t) = {\hat{F}}^{- 1} (ϕ^{- 1} (x, t), t) . \end{matrix}

类似于 (2.33)，我们现在可以用张量 $F (x, t)$ 表达长度 $∥ d X ∥$

\begin{matrix} (2.39) & ∥ d X ∥^{2} = ({\hat{F}}^{- 1} (X_{P}, t) d x, {\hat{F}}^{- 1} (X_{P}, t) d x) \end{matrix}

\begin{matrix} (2.40) & = (F (x_{p}, t) d x, F (x_{p}, t) d x) \end{matrix}

\begin{matrix} (2.41) & = d x F^{T} (x_{p}, t) F (x_{p}, t) d x \end{matrix}

\begin{matrix} (2.42) & = d x b^{- 1} d x, \end{matrix}

其中 $b^{- 1} = F^{T} F = {\hat{F}}^{- T} {\hat{F}}^{- 1}$ 是左 Cauchy–Green 张量 $b$ 的逆， $b$ 定义为

\begin{matrix} (2.43) & b = \hat{F} {\hat{F}}^{T} . \end{matrix}

与 (2.36) 中 $\hat{C}$ 的定义相反，此处称之为左张量，是因为在 $b$ 的定义中 $\hat{F}$ 位于左侧。

现在我们考察线段 $d X$ 在变形下的长度平方变化：

\begin{matrix} (2.44) & \frac{1}{2} (∥ d x ∥^{2} - ∥ d X ∥^{2}) = \frac{1}{2} (d X^{T} \hat{C} d X - d X^{T} d X) \end{matrix}

\begin{matrix} (2.45) & = d X^{T} [\frac{1}{2} (\hat{C} - I)] d X \end{matrix}

\begin{matrix} (2.46) & = d X^{T} \hat{E} d X, \end{matrix}

其中 $\hat{E} = \frac{1}{2} (\hat{C} - I)$ 描述 Lagrange 应变张量（或 Green 应变张量）， $I$ 为（二阶）单位张量。基于方程 (2.46)，长度平方变化也可以用线段 $d x$ 表达：

\begin{matrix} (2.47) & \frac{1}{2} (∥ d x ∥^{2} - ∥ d X ∥^{2}) = d X^{T} \hat{E} d X, \end{matrix}

\begin{matrix} (2.48) & = d x^{T} {\hat{F}}^{- T} \hat{E} {\hat{F}}^{- 1} d x, \end{matrix}

\begin{matrix} (2.49) & = d x^{T} [\frac{1}{2} (I - {\hat{F}}^{- T} {\hat{F}}^{- 1})] d x, \end{matrix}

\begin{matrix} (2.50) & = d x^{T} [\frac{1}{2} (I - F^{T} F)] d x . \end{matrix}

与之相应的张量 $e = \frac{1}{2} (I - F^{T} F)$ 称为 Euler 应变张量（或 Almansi 应变张量）。

2.4 应力张量

设为连续介质，作用在流体任意控制体 $V (t)$ 上的力有两类。第一类是长程力，例如重力，它施加单位体积上的力。这类能够渗透到流体内部的力称为体积力或体力。第二类是与分子相互作用直接相关的短程力。这类力通常源于迁移或振荡分子造成的动量输运。这类被称为应力的力，只能作用在一层薄层上，其穿透深度相对于表面元的尺寸很小。我们用 $d s$ 记 $x$ 在表面元 $S$ 上的一个邻域。作用在 $d s$ 上的平均力用 Cauchy 牵引矢量 $t (x, t, n)$ 表示，其中 $n$ 是与表面元 $d s$ 相关联的单位法矢。所采用的约定是： $t (x, t, n)$ 是流体在表面元 $n$ 所指向一侧所施加的应力。 $t (x, t, n)$ 中方向与 $n$ 相同的法向分量表示张力（拉伸）。

Cauchy 的基本结论断言，牵引矢量对取向 $n$ 这种先验上很复杂的依赖关系，实际上是线性的，并由一个二阶张量所中介。

定理 2.2（Cauchy 应力定理）

设牵引矢量 $t (x, t, n)$ 关于 $x$ 连续，且线动量在每个物质体积上守恒。则 $t$ 对单位法矢 $n$ 的依赖是线性的：存在一个二阶张量场 $σ (x, t)$ ，即 Cauchy 应力张量，使得

\begin{matrix} (2.51) & t (x, t, n) = σ (x, t) n \end{matrix}

对每一个单位法矢 $n$ 成立。

证明（Cauchy 四面体）。 固定 $x$ 与 $t$ ，考察一个小四面体，它有三个面与坐标轴正交，还有第四个倾斜面，其外法单位矢为 $n = (n_{1}, \dots, n_{d})$ 。若倾斜面的面积为 $A$ ，则外法为 $- e_{i}$ 的面的面积为 $A_{i} = n_{i} A$ 。在该四面体上对线动量作平衡，则表面牵引、体力与惯性力满足

A t (n) + \sum_{i = 1}^{d} A_{i} t (- e_{i}) + (体积项) = 0.

若 $ε$ 记四面体的直径，则体积项（体力与惯性）按 $ε^{d}$ 标度，而表面项按 $ε^{d - 1}$ 标度。两边除以 $A \sim ε^{d - 1}$ 并令 $ε \to 0$ ，体积项消失；再利用 $t (- e_{i}) = - t (e_{i})$ 与 $A_{i} / A = n_{i}$ ，便剩下

t (n) = \sum_{i = 1}^{d} n_{i} t (e_{i}) .

定义 $σ$ 为其第 $i$ 列是法矢为 $e_{i}$ 的面上牵引 $t (e_{i})$ 的张量，即 $σ_{j i} = (t (e_{i}))_{j}$ ，则右端即为 $σ n$ ，这正是 (2.51)。 $◻$

因此，分量 $σ_{j i}$ 是作用在法向指向 $e_{i}$ 方向的表面元上的牵引的第 $j$ 个分量；对角元是法向应力，非对角元是剪切应力。此外，Cauchy 应力张量是对称的， $σ = σ^{T}$ ，这是 2.6 节所建立的角动量平衡的结果。

2.5 Reynolds 输运定理

对于给定量 $T$ ，其 Lagrange 描述与 Euler 描述之间成立如下关系：

\begin{matrix} (2.52) & T (ϕ (X, t), t) = T (x, t) = \hat{T} (X, t) = \hat{T} (ϕ^{- 1} (x, t), t) . \end{matrix}

Jacobi 矩阵 $D_{X} ϕ$ 及其行列式 $\hat{J} (X, t) := det (D_{X} ϕ (X, t))$ 对随后考察的体积分起着关键作用。

引理 2.3

成立

\begin{matrix} (2.53) & \partial_{t} \hat{J} (X, t) = \hat{J} (X, t) (\nabla_{x} \cdot v) (ϕ (X, t), t) \end{matrix}

以及

\begin{matrix} (2.54) & \hat{J} (X, t) > 0 对所有 X, t . \end{matrix}

证明。 利用链式法则、 $\hat{v}$ 的定义以及 (2.52)，我们首先得到

\begin{matrix} (2.55) & \frac{\partial}{\partial t} \frac{\partial ϕ_{i} (X, t)}{\partial X_{j}} = \frac{\partial}{\partial X_{j}} {\hat{v}}_{i} (X, t) = \frac{\partial}{\partial X_{j}} v_{i} (ϕ (X, t), t) = \sum_{k = 1}^{d} \frac{\partial}{\partial x_{k}} v_{i} (x, t) \frac{\partial}{\partial X_{j}} ϕ_{k} (X, t) . \end{matrix}

现在我们回顾一些线性代数的结论：设 $S \in R^{d \times d}$ 是一个行列式为 $det (S)$ 的矩阵。用 $S^{(i, j)} \in R^{d - 1 \times d - 1}$ 记删去 $S$ 的第 $i$ 行、第 $j$ 列后所得的矩阵，并用 $A_{i j} := (- 1)^{i + j} det (S^{(i, j)})$ 记行列式的代数余子式。用这些记号，下列公式成立：

\begin{matrix} (2.56) & (S^{- 1})_{i j} = \frac{1}{det (S)} A_{j i} \end{matrix}

\begin{matrix} (2.57) & δ_{i k} det (S) = \sum_{j = 1}^{d} A_{i j} S_{k j} . \end{matrix}

由 (2.57) 我们直接得到

\begin{matrix} (2.58) & \frac{\partial}{\partial S_{i j}} det (S) = A_{i j} . \end{matrix}

以 $S = D_{X} ϕ$ 、 $\hat{J} = det (S)$ ，利用 (2.56)–(2.58) 得到：

\begin{aligned} \frac{\partial}{\partial t} \hat{J} & = \sum_{i, j = 1}^{d} \frac{\partial \hat{J}}{\partial [\frac{\partial ϕ_{i}}{\partial X_{j}}]} \frac{\partial}{\partial t} \frac{\partial ϕ_{i} (X, t)}{\partial X_{j}} & (链式法则) \\ = \sum_{i, j = 1}^{d} A_{i j} \sum_{k = 1}^{d} \frac{\partial}{\partial x_{k}} v_{i} (x, t) \frac{\partial}{\partial X_{j}} ϕ_{k} (X, t) & (由 (2.58) 与 (2.55)) \\ = \sum_{i, j, k = 1}^{d} A_{i j} \frac{\partial}{\partial X_{j}} ϕ_{k} (X, t) \frac{\partial}{\partial x_{k}} v_{i} (x, t) \\ = \sum_{i, k = 1}^{d} δ_{i k} \hat{J} \frac{\partial}{\partial x_{k}} v_{i} (x, t) & (由 (2.57)) \\ = \hat{J} \sum_{i = 1}^{d} \frac{\partial}{\partial x_{i}} v_{i} (x, t) \\ = \hat{J} (X, t) (\nabla_{x} \cdot v) (ϕ (X, t), t) \end{aligned}

断言 (2.54) 由假设 2.1 直接得到：由于 $ϕ (\cdot, t)$ 是微分同胚， $D_{X} ϕ (X, t)$ 可逆，从而对所有 $X, t$ 有 $\hat{J} (X, t) \neq 0$ 。又因 $ϕ (X, 0) = X$ ，故 $\hat{J} (X, 0) = 1 > 0$ 。因此对所有 $X, t$ 有 $\hat{J} (X, t) > 0$ 。 $◻$

现在我们考察体积及其中所含的量。设 $\hat{V} \subset R^{d}$ 记一组物质点， $V (t) := ϕ (\hat{V}, t)$ 、 $t \geq 0$ 记这些点在 $t$ 时刻所占据的体积。由假设 2.1，有 $\hat{V} = V (0)$ 。

定理 2.4（Reynolds 输运定理）

设 $\hat{V} \subset R^{d}$ 为 Lebesgue 可测集， $Ω \subset R^{d}$ 满足对所有 $t \geq 0$ 有 $V (t) \subset Ω$ 。假设 $ϕ$ 与 $f : Ω \times (0, \infty) \to R$ 充分可微，且其导数 $D_{X} ϕ (\cdot, t)$ 与 $\nabla_{x} f (\cdot, t)$ 分别在 $\hat{V}$ 与 $Ω$ 上对每个 $t \in (0, \infty)$ 可积。则成立

\begin{matrix} (2.59) & \frac{d}{d t} \int_{V (t)} f (x, t) d x = \int_{V (t)} \partial_{t} f (x, t) + \nabla_{x} \cdot (f v) (x, t) d x . \end{matrix}

证明。 成立

\begin{aligned} \frac{d}{d t} \int_{V (t)} f (x, t) d x & = \frac{d}{d t} \int_{\hat{V}} f (ϕ (X, t), t) \hat{J} (X, t) d X & (换元法) \\ = \int_{\hat{V}} \frac{d}{d t} (f (ϕ (X, t), t) \hat{J} (X, t)) d X & (Leibniz 积分法则) \\ = \int_{\hat{V}} \partial_{t} f (ϕ (X, t), t) \hat{J} (X, t) \\ + \hat{J} (X, t) \sum_{i = 1}^{d} \frac{\partial}{\partial x_{i}} f (ϕ (X, t), t) \frac{\partial ϕ_{i}}{\partial t} (X, t) & (乘积与链式法则) \\ + f (ϕ (X, t), t) \partial_{t} \hat{J} (X, t) d X \\ = \int_{\hat{V}} \partial_{t} f (ϕ (X, t), t) \hat{J} (X, t) \\ + \hat{J} (X, t) \nabla_{x} f (ϕ (X, t), t) \cdot v (ϕ (X, t), t) \\ + f (ϕ (X, t), t) \hat{J} (X, t) (\nabla_{x} \cdot v) (ϕ (X, t), t) d X & (引理 2.3) \\ = \int_{\hat{V}} \partial_{t} f (ϕ (X, t), t) \hat{J} (X, t) \\ + (\nabla_{x} \cdot (f v)) (ϕ (X, t), t) \hat{J} (X, t) d X \\ = \int_{V (t)} \partial_{t} f (x, t) + \nabla_{x} \cdot (f v) (x, t) d x . & (换元法) \end{aligned}

$◻$

Reynolds 输运定理 2.4 使我们能就物理量的守恒得出如下结论：设 $\hat{V} \subset R^{d}$ 记一个有界域（即开且连通），它在初始时刻 $t = 0$ 含有给定量的物质粒子，并记住此前引入的映射 $ϕ (\cdot, t) : \hat{V} \to V (t)$ 。我们考察一族函数

\begin{matrix} (2.60) & f (\cdot, t) : V (t) \to R \end{matrix}

它应当刻画某物理量（例如质量、动量或能量）的（相对于体积的）密度。于是，

\begin{matrix} (2.61) & F (t) := \int_{V (t)} f (x, t) d x \end{matrix}

正是体积 $V (t)$ 中所含该量的总量。利用定理 2.4 与 Gauss 散度定理，我们得到

\begin{matrix} (2.62) & \begin{aligned} \frac{d}{d t} F (t) & = \int_{V (t)} \partial_{t} f (x, t) + \nabla_{x} \cdot (f v) (x, t) d x \\ = \int_{V (t)} \partial_{t} f (x, t) d x + \int_{\partial V (t)} f (v \cdot n) (x, t) d σ . \end{aligned} \end{matrix}

方程 (2.62) 基本上表明：量 $F$ 的时间变化率由内部密度变化，以及以速度 $v$ 穿过其表面 $\partial V (t)$ 进入体积的密度量所决定。后一项称为通量（flux）。

2.6 守恒律

现在我们可以把 2.2 节的运动学、2.4 节的应力张量以及 Reynolds 输运定理结合起来，推导流体运动的控制方程。每一条平衡律都是通过对任意物质体积 $V (t)$ 假定某一物理量的守恒，并把所得的积分恒等式局部化而得到的。

2.6.1 质量守恒

设 $ρ (x, t)$ 记质量密度。物质体积 $V (t)$ 中所含的质量为 $\int_{V (t)} ρ d x$ ，质量守恒要求它在时间上保持恒定，

\frac{d}{d t} \int_{V (t)} ρ d x = 0.

对 $f = ρ$ 应用 Reynolds 输运定理（定理 2.4）得到 $\int_{V (t)} (\partial_{t} ρ + \nabla \cdot (ρ v)) d x = 0$ 。由于物质体积是任意的且被积函数连续，该被积函数必须逐点为零，这给出连续性方程

\begin{matrix} (2.63) & \partial_{t} ρ + \nabla \cdot (ρ v) = 0. \end{matrix}

对于不可压缩流体，每个物质粒子的密度都是常数， $\frac{D ρ}{D t} = 0$ ，于是 (2.63) 化为不可压缩约束

\begin{matrix} (2.64) & \nabla \cdot v = 0. \end{matrix}

(2.63) 的一个有用推论是如下恒等式，它对任意可微（标量或矢量值）场 $ψ$ 成立：对 $f = ρ ψ$ 应用输运定理，并用 (2.63) 消去乘 $ψ$ 的那些项，

\begin{matrix} (2.65) & \frac{d}{d t} \int_{V (t)} ρ ψ d x = \int_{V (t)} ρ \frac{D ψ}{D t} d x . \end{matrix}

2.6.2 线动量平衡

对物质体积应用牛顿第二定律，把线动量的变化率等同于作用其上的总力，后者由体力 $ρ b$ 与 (2.51) 的表面牵引 $t = σ n$ 组成，

\frac{d}{d t} \int_{V (t)} ρ v d x = \int_{V (t)} ρ b d x + \int_{\partial V (t)} σ n d σ .

由散度定理，面积分被转化为体积分，

\int_{\partial V (t)} σ n d σ = \int_{V (t)} \nabla \cdot σ d x,

而左端由质量比输运恒等式 (2.65) 取 $ψ = v$ 加以变换。在任意体积 $V (t)$ 上令被积函数相等，给出 Cauchy 运动方程

\begin{matrix} (2.66) & ρ \frac{D v}{D t} = \nabla \cdot σ + ρ b . \end{matrix}

此处物质加速度 $\frac{D v}{D t} = \partial_{t} v + (\nabla v) v$ 恰好是 (2.10) 中引入的加速度场。

2.6.3 角动量平衡

在没有分布体力偶的情形下，角动量平衡对应力施加了一个逐点约束，我们在 2.4 节已作了预告。

命题 2.5（Cauchy 应力的对称性）

在没有体力偶的情形下，角动量平衡蕴含 Cauchy 应力张量是对称的， $σ = σ^{T}$ 。

证明思路。 物质体积关于原点的角动量平衡为

\frac{d}{d t} \int_{V (t)} ρ (x \times v) d x = \int_{V (t)} ρ (x \times b) d x + \int_{\partial V (t)} x \times (σ n) d σ .

用 (2.65) 变换左端、用散度定理变换面项，然后减去 $x$ 与 Cauchy 运动方程 (2.66) 的叉积，体力与加速度的所有贡献都相互抵消。剩下的是一个体积分，其被积函数是 $σ$ 反对称部分的轴矢量。由于 $V (t)$ 任意，该部分必须为零，即 $σ = σ^{T}$ 。 $◻$

2.6.4 能量平衡

为完整起见，我们记录热力学第一定律，它把内能的变化率与应力功率及热交换加以平衡。设 $e$ 为比内能， $q$ 为热流， $r$ 为单位质量的热源，则局部形式为

\begin{matrix} (2.67) & ρ \frac{D e}{D t} = σ : D - \nabla \cdot q + ρ r, \end{matrix}

其中 $σ : D$ 是应力功率， $D$ 是 2.2 节的变形率张量。能量平衡对于可压缩或非等温流动是不可或缺的；而对于本书所强调的不可压缩、等温模型，它与 (2.63) 和 (2.66) 解耦，故不再进一步讨论。

结语

方程 (2.63) 与 (2.66) 是流体运动的控制方程。它们所含的未知量——速度 $v$ 、密度 $ρ$ 与应力 $σ$ ——比标量方程的数目更多，因而并不闭合：应力还必须通过一条本构律与运动联系起来。牛顿闭合，以及由此导出的 Stokes 方程，是第 3 章的主题；而由此得到的模型的有限元离散化所依据的变分形式则将在第 4 章展开。

3 从守恒律到 Stokes 方程

第 2 章的平衡律——质量守恒与线动量平衡——只有在 Cauchy 应力张量 $σ$ 被用运动本身表达出来之后，才能确定流体的运动。本章补上这一环：为不可压缩牛顿流体给出应力的本构闭合，由此得到 Navier–Stokes 方程，再在蠕流极限下导出第 4 章所要离散的 Stokes 方程。（本章只讨论 Stokes 方程的推导；非牛顿流体的建模不在此列。）

3.1 本构闭合问题

对常密度 $ρ$ 的不可压缩流体，第 2 章的守恒律为

\begin{matrix} (3.1) & ρ \frac{D v}{D t} = div σ + ρ b ， div v = 0 ， \end{matrix}

其中 $b$ 是单位质量的给定体力。系统 (3.1) 含 $d + 1$ 个标量方程，未知量却是速度 $v$ 、压强与对称应力 $σ$ ，因而并不闭合：还需要一条把 $σ$ 与运动学联系起来的本构律。对不可压缩流体，方便的做法是把应力拆成压强部分与偏量部分（后者称为附加应力 $τ$ ）：

\begin{matrix} (3.2) & σ = - p I + τ ， tr τ = 0 ， \end{matrix}

其中压强 $p = - \frac{1}{d} tr σ$ 作为 Lagrange 乘子施行不可压缩约束， $I$ 为二阶单位张量。于是建模任务归结为给定附加应力 $τ$ 。全章中 $D$ 记变形率张量， $L = D + W$ 为速度梯度及其对称、反对称部分（见第 2 章）。

3.2 牛顿流体与 Navier–Stokes 方程

最简单的容许闭合，在变形率张量中是线性且各向同性的，它定义了牛顿流体。

定义 3.1（牛顿流体）

称不可压缩流体为牛顿流体，若其附加应力线性依赖于变形率，

\begin{matrix} (3.3) & τ = 2 μ D ， \end{matrix}

其中 $μ > 0$ 为常数动力粘度。

把 (3.2)–(3.3) 代入 (3.1)。粘性项的散度可借不可压缩性化简：

div (2 μ D (v)) = μ div (\nabla v + \nabla v^{T}) = μ (Δ v + \nabla (div v)) = μ Δ v ，

其中用到按分量的恒等式 $[div (\nabla v^{T})]_{i} = \partial_{j} \partial_{i} v_{j} = \partial_{i} (div v)$ 以及 $div v = 0$ 。于是 $div σ = - \nabla p + μ Δ v$ ，代入 (3.1) 得不可压缩 Navier–Stokes 方程

\begin{matrix} (3.4) & ρ \frac{D v}{D t} = - \nabla p + μ Δ v + ρ b ， div v = 0. \end{matrix}

它是不可压缩流动的基石模型；牛顿律 (3.3) 的特征在于粘度 $μ$ 是与流动无关的材料常数。

3.3 蠕流极限：Stokes 方程

Navier–Stokes 方程 (3.4) 左端的物质加速度 $\frac{D v}{D t} = \partial_{t} v + (v \cdot \nabla) v$ 是非线性的（源于对流项 $(v \cdot \nabla) v$ ）。以特征速度 $U$ 、特征长度 $L$ 作无量纲化，惯性项与粘性项之比由 Reynolds 数 $Re = ρ U L / μ$ 度量。当 $Re ≪ 1$ （缓慢、粘性主导的蠕流）时，对流加速度相对于粘性力可忽略；再限于定常流动 $\partial_{t} v = 0$ ，则整个物质加速度都被略去， $\frac{D v}{D t} \to 0$ 。方程 (3.4) 遂退化为线性方程

- μ Δ v + \nabla p = ρ b ， div v = 0.

把粘性项写回对称的变形率形式 $μ Δ v = div (2 μ D (v))$ （对无散度场成立），并记 $u := v$ 、 $f := ρ b$ ，再配以壁面无滑移边界条件 $u = 0$ 于 $\partial Ω$ ，便得到定常 Stokes 问题的强形式

\begin{matrix} (3.5) & - div (2 μ D (u)) + \nabla p = f ， div u = 0 在 Ω 内 ， u = 0 在 \partial Ω 上 。 \end{matrix}

这正是第 4 章将要离散的强形式（在那里重新记作 (4.2)）。之所以偏爱对称的 $div (2 μ D (u))$ 而非等价的 $μ Δ u$ ，是因为前者显式地保留了变形率张量 $D$ ——它与 Korn 不等式自然耦合，从而在第 4 章的弱形式中给出强制性。

4 变分形式

第 2 章与第 3 章所汇集的本构流动模型都是强形式的偏微分方程，它们要求逐点存在的导数，而数值逼近往往难以直接提供这类导数。它们的有限元处理转而依赖于一个等价的变分（弱）形式：把方程对光滑函数作检验，并通过分部积分放松可微性要求。本章发展这一形式，建立所得问题的适定性，并由此导出离散代数系统——从而架起通向第 6 章鞍点求解器的桥梁。

我们集中讨论定常 Stokes 问题，它是从第 3 章的 Navier–Stokes 方程 (3.4) 出发、略去对流加速度所得的线性原型；完整非线性问题所需的修改在本章末尾指出。全章中， $Ω \subset R^{d}$ （ $d \in {2, 3}$ ）记一个具有 Lipschitz 边界 $\partial Ω$ 的有界域， $(\cdot, \cdot)$ 是 $L^{2} (Ω)$ 内积；对张量， $A : B = \sum_{i, j} A_{i j} B_{i j}$ 记 Frobenius 缩并。

4.1 函数空间与弱导数

有限元方法并不处理经典的、逐点意义下的解，而是处理只在弱的、积分意义下具有导数的函数。称函数 $u \in L^{2} (Ω)$ 具有弱导数 $\partial_{i} u = g \in L^{2} (Ω)$ ，若

\int_{Ω} u \partial_{i} φ d x = - \int_{Ω} g φ d x 对所有 φ \in C_{c}^{\infty} (Ω) ，

该恒等式即分部积分公式，其边界项因 $φ$ 的紧支集而被抑制。把弱一阶导数平方可积的函数收集起来，得到 Sobolev 空间

H^{1} (Ω) = {u \in L^{2} (Ω) : \partial_{i} u \in L^{2} (Ω), i = 1, \dots, d} ，

它在范数 $∥ u ∥_{H^{1}}^{2} = ∥ u ∥_{L^{2}}^{2} + ∥ \nabla u ∥_{L^{2}}^{2}$ 下是一个 Hilbert 空间。齐次 Dirichlet 边界条件通过闭子空间 $H_{0}^{1} (Ω)$ 施加， $H_{0}^{1} (Ω)$ 是 $C_{c}^{\infty} (Ω)$ 在 $H^{1}$ 范数下的闭包。在 $H_{0}^{1} (Ω)$ 上成立 Poincaré–Friedrichs 不等式 $∥ u ∥_{L^{2}} \leq C_{Ω} ∥ \nabla u ∥_{L^{2}}$ ，因此半范数 $| u |_{H^{1}} = ∥ \nabla u ∥_{L^{2}}$ 在 $H_{0}^{1} (Ω)$ 上是一个范数，且与 $∥ \cdot ∥_{H^{1}}$ 等价。

对于不可压缩流动问题，速度是一个矢量场，其每个分量都属于 $H_{0}^{1}$ ；而压强只通过其梯度进入方程，从而只被确定到相差一个可加常数。这促使我们引入如下空间对

\begin{matrix} (4.1) & V := [H_{0}^{1} (Ω)]^{d} ， Q := L_{0}^{2} (Ω) = {q \in L^{2} (Ω) : \int_{Ω} q d x = 0} ， \end{matrix}

分别配以范数 $∥ v ∥_{V} = | v |_{H^{1}}$ 与 $∥ q ∥_{Q} = ∥ q ∥_{L^{2}}$ 。 $Q$ 中的零均值约束固定了那个未定常数，使压强唯一。

4.2 Stokes 问题的弱形式

速度 $u$ 与压强 $p$ 的定常 Stokes 问题，其强形式为

\begin{matrix} (4.2) & - div (2 μ D (u)) + \nabla p = f ， div u = 0 在 Ω 内 ， u = 0 在 \partial Ω 上 ， \end{matrix}

其中 $D (u) = \frac{1}{2} (\nabla u + \nabla u^{T})$ 是第 2 章的变形率张量， $μ > 0$ 为常数动力粘度，体力 $f \in [L^{2} (Ω)]^{d}$ 。它刻画对流加速度（见 (3.4)）可忽略的、由粘性主导的缓慢流动。

为过渡到弱形式，取检验速度 $v \in V$ ，用 $v$ 点乘动量方程并在 $Ω$ 上积分。对应力项作分部积分，利用 $v = 0$ 于 $\partial Ω$ 以及 $D$ 的对称性，可把一次求导转移到检验函数上。下面把这一步用 Gauss 散度定理逐步展开。

准备两个 Gauss 恒等式。 对二阶张量场 $S$ （其散度按分量为 $(div S)_{i} = \sum_{j} \partial_{j} S_{i j}$ ）与矢量场 $v$ ，由乘积法则 $\partial_{j} (S_{i j} v_{i}) = (\partial_{j} S_{i j}) v_{i} + S_{i j} \partial_{j} v_{i}$ ，对 $i, j$ 求和后左端为矢量场 $S^{T} v$ 的散度。对该散度应用经典 Gauss 定理，并注意 $(S^{T} v) \cdot n = \sum_{i, j} S_{i j} v_{i} n_{j} = v \cdot (S n)$ ，即得张量版分部积分公式；对压强梯度项则用标量恒等式 $div (p v) = \nabla p \cdot v + p div v$ 。

张量版与压力项的分部积分

\begin{matrix} (⋆) & \int_{Ω} (div S) \cdot v d x = \int_{\partial Ω} (S n) \cdot v d s - \int_{Ω} S : \nabla v d x ， \end{matrix}

\begin{matrix} (⋆ ⋆) & \int_{Ω} \nabla p \cdot v d x = \int_{\partial Ω} p (v \cdot n) d s - \int_{Ω} p div v d x . \end{matrix}

检验动量方程。 用 $v$ 点乘 (4.2) 的动量方程并积分，再对第一项用 $(⋆)$ （取 $S = 2 μ D (u)$ ）、对第二项用 $(⋆ ⋆)$ ：

\int_{Ω} 2 μ D (u) : \nabla v d x - \int_{\partial Ω} (2 μ D (u) n) \cdot v d s - \int_{Ω} p div v d x + \int_{\partial Ω} p (v \cdot n) d s = \int_{Ω} f \cdot v d x .

边界项消失。 因为 $v \in V = [H_{0}^{1} (Ω)]^{d}$ ，故 $v = 0$ 于 $\partial Ω$ ，两个边界积分同时为零，于是

\int_{Ω} 2 μ D (u) : \nabla v d x - \int_{Ω} p div v d x = \int_{Ω} f \cdot v d x .

利用 $D$ 的对称性。 把 $\nabla v = D (v) + W (v)$ 分解为对称与反对称部分（见 2.2 节）。由于 $D (u)$ 对称、 $W (v)$ 反对称，双点积 $D (u) : W (v) = 0$ ，故 $D (u) : \nabla v = D (u) : D (v)$ ，从而得到弱动量平衡

\int_{Ω} 2 μ D (u) : D (v) d x - \int_{Ω} p div v d x = \int_{Ω} f \cdot v d x .

检验不可压缩约束。 用检验压强 $q \in Q$ 乘以不可压缩约束并积分，得 $\int_{Ω} q div u d x = 0$ 。引入双线性形式

\begin{matrix} (4.3) & a (u, v) := \int_{Ω} 2 μ D (u) : D (v) d x ， b (v, q) := - \int_{Ω} q div v d x ， \end{matrix}

它们分别定义在 $V \times V$ 与 $V \times Q$ 上；再引入载荷泛函 $ℓ (v) := \int_{Ω} f \cdot v d x$ ，则问题取如下对称的**混合（mixed）**形式。 $b$ 定义中的负号是刻意为之：它使同一个 $b$ 既出现在动量方程里（ $b (v, p)$ ，来自 $- \int p div v$ ）又出现在约束里（ $b (u, q)$ ），从而给出对称的鞍点结构。

定义 4.1（Stokes 问题的弱形式）

求 $(u, p) \in V \times Q$ ，使得

\begin{matrix} (4.4) & \begin{aligned} a (u, v) + b (v, p) & = ℓ (v) & 对所有 v \in V ， \\ b (u, q) & = 0 & 对所有 q \in Q . \end{aligned} \end{matrix}

这两行分别表示弱动量平衡——其中压强作为约束的 Lagrange 乘子通过 $b (v, p)$ 进入——与弱不可压缩性。(4.4) 的结构，即 $V$ 上一个对称双线性形式 $a$ 与一个约束形式 $b$ 相耦合，正是第 6 章代数鞍点系统的连续原型。

4.3 适定性：Lax–Milgram 与 Babuška–Brezzi 理论

(4.4) 的适定性依赖变分学的两个经典结论。第一个针对单个强制双线性形式。

定理 4.2（Lax–Milgram）

设 $H$ 为 Hilbert 空间， $a : H \times H \to R$ 为双线性形式，满足有界性 $| a (u, v) | \leq M ∥ u ∥_{H} ∥ v ∥_{H}$ 与强制性 $a (v, v) \geq α ∥ v ∥_{H}^{2}$ （ $α > 0$ ）。则对每个有界线性泛函 $ℓ \in H^{'}$ ，问题 $a (u, v) = ℓ (v)$ （对所有 $v \in H$ ）有唯一解，且 $∥ u ∥_{H} \leq α^{- 1} ∥ ℓ ∥_{H^{'}}$ 。

对 Stokes 形式 $a$ ，有界性是显然的，而强制性由 Korn 不等式给出：在 $V$ 上，由 $2 ∥ D (v) ∥_{L^{2}}^{2} \geq | v |_{H^{1}}^{2}$ （ $v \in V$ ）可得 $a (v, v) = 2 μ ∥ D (v) ∥_{L^{2}}^{2} \geq μ ∥ v ∥_{V}^{2}$ 。若非不可压缩约束的存在，Lax–Milgram 本可立即解决问题；与压强的耦合要求把它推广到鞍点情形 [7, 16]。

于是考察抽象混合问题：给定 $V \times V$ 上有界双线性形式 $a$ 、 $V \times Q$ 上有界双线性形式 $b$ ，以及数据 $ℓ \in V^{'}$ 、 $χ \in Q^{'}$ ，求 $(u, p) \in V \times Q$ ，使得 $a (u, v) + b (v, p) = ℓ (v)$ （对所有 $v$ ）且 $b (u, q) = χ (q)$ （对所有 $q$ ）。其分析围绕核空间展开

Z := {v \in V : b (v, q) = 0 对所有 q \in Q} ，

对 Stokes 形式而言，它就是弱意义下无散度的速度所构成的空间。

定义 4.3（inf–sup 条件）

称形式 $b$ 满足 inf–sup 条件（或 Ladyzhenskaya–Babuška–Brezzi 条件），若存在常数 $β > 0$ ，使得

\begin{matrix} (4.5) & inf_{q \in Q ∖ {0}} sup_{v \in V ∖ {0}} \frac{b (v, q)}{∥ v ∥_{V} ∥ q ∥_{Q}} \geq β . \end{matrix}

定理 4.4（Babuška–Brezzi）

设 $a$ 与 $b$ 有界， $a$ 在核空间 $Z$ 上强制，且 $b$ 满足 inf–sup 条件 (4.5)。则抽象混合问题有唯一解 $(u, p)$ ，它连续依赖于数据，

∥ u ∥_{V} + ∥ p ∥_{Q} \leq C (∥ ℓ ∥_{V^{'}} + ∥ χ ∥_{Q^{'}}) ，

其中 $C$ 只依赖于强制常数 $α$ 、inf–sup 常数 $β$ ，以及 $a$ 与 $b$ 的连续性常数。

其证明及最优常数可在专著 [5, 8, 18] 中找到。应用于 Stokes 问题： $Z$ 上的强制性由上面的 Korn 估计给出（它在整个 $V$ 上成立，故在 $Z$ 上亦然），而 inf–sup 条件相当于断言散度算子把 $V$ 满射到 $Q$ 上、且有一个有界右逆：对每个 $q \in L_{0}^{2} (Ω)$ ，存在 $v \in V$ 满足 $div v = q$ 且 $∥ v ∥_{V} \leq C ∥ q ∥_{L^{2}}$ 。Lipschitz 域的这一非平凡性质 [5, 18] 补全了定理 4.4 的假设，从而得到如下结论。

定理 4.5（Stokes 问题的适定性）

对每个 $f \in [L^{2} (Ω)]^{d}$ ，弱 Stokes 问题 (4.4) 有唯一解 $(u, p) \in V \times Q$ ，它连续依赖于 $f$ 。

注 4.6（极小化与鞍点结构）。 由于 $a$ 对称，(4.4) 的解是 Lagrange 泛函 $L (v, q) = \frac{1}{2} a (v, v) - ℓ (v) + b (v, q)$ 的唯一鞍点：速度在无散度场 $Z$ 上极小化粘性耗散 $\frac{1}{2} a (v, v) - ℓ (v)$ ，而压强则是施行约束 $div v = 0$ 的 Lagrange 乘子。这是第 6 章离散 Lagrange 泛函与约束极小化的连续对应物，也是该章所专注的鞍点结构的首次出现。

4.4 Galerkin 逼近与 Céa 引理

有限元方法把 $V$ 与 $Q$ 替换为有限维子空间 $V_{h} \subset V$ 与 $Q_{h} \subset Q$ ，它们由 $Ω$ 的一个特征尺度为 $h$ 的网格 $T_{h}$ 上的分片多项式构造而成。(4.4) 的 Galerkin 逼近把试探函数与检验函数都限制到这些子空间上。

定义 4.7（离散 Stokes 问题）

求 $(u_{h}, p_{h}) \in V_{h} \times Q_{h}$ ，使得

\begin{matrix} (4.6) & \begin{aligned} a (u_{h}, v_{h}) + b (v_{h}, p_{h}) & = ℓ (v_{h}) & 对所有 v_{h} \in V_{h} ， \\ b (u_{h}, q_{h}) & = 0 & 对所有 q_{h} \in Q_{h} . \end{aligned} \end{matrix}

4.3 节的抽象理论逐字适用于 (4.6)，但有一个决定性的告诫：inf–sup 条件不会自动从连续层面继承下来，因为 (4.5) 中的上确界现在是在更小的空间 $V_{h}$ 上取的。因此稳定性要求空间对 $(V_{h}, Q_{h})$ 满足一个离散 inf–sup 条件，且其常数关于 $h$ 一致地远离零。

假设 4.8（离散 inf–sup 稳定性）

存在常数 $β_{*} > 0$ ，与 $h$ 无关，使得

\begin{matrix} (4.7) & inf_{q_{h} \in Q_{h} ∖ {0}} sup_{v_{h} \in V_{h} ∖ {0}} \frac{b (v_{h}, q_{h})}{∥ v_{h} ∥_{V} ∥ q_{h} ∥_{Q}} \geq β_{*} 对所有 h . \end{matrix}

在假设 4.8 之下，离散问题由定理 4.4 关于 $h$ 一致地适定，且逼近是拟最优的。

定理 4.9（拟最优误差估计）

设 $(u, p)$ 解 (4.4)， $(u_{h}, p_{h})$ 解 (4.6)，且假设 4.8 成立。则

\begin{matrix} (4.8) & ∥ u - u_{h} ∥_{V} + ∥ p - p_{h} ∥_{Q} \leq C (inf_{v_{h} \in V_{h}} ∥ u - v_{h} ∥_{V} + inf_{q_{h} \in Q_{h}} ∥ p - q_{h} ∥_{Q}) ， \end{matrix}

其中 $C$ 只依赖于连续性常数以及 $α, β_{*}$ ——因而与 $h$ 无关。

估计 (4.8) 是 Céa 引理的混合形式对应物：至多相差一个常数，有限元解与 $V_{h} \times Q_{h}$ 中可得的最佳逼近一样精确；再结合多项式逼近估计，它给出关于 $h$ 的收敛率 [8, 16]。决定性的假设是一致离散 inf–sup 条件 (4.7)。违反它的单元对——尤其是等阶速度–压强插值——会容许虚假压强模态，必须加以稳定化；而 inf–sup 稳定的单元对（例如 Taylor–Hood 单元）则满足它。这些问题，尤其是 (4.7) 的代数表现，将在第 6 章详细考察。

4.5 从变分问题到代数系统

选取 $V_{h}$ 的基 ${φ_{i}}_{i = 1}^{n}$ 与 $Q_{h}$ 的基 ${ψ_{k}}_{k = 1}^{m}$ ，并展开

u_{h} = \sum_{j = 1}^{n} U_{j} φ_{j} ， p_{h} = \sum_{l = 1}^{m} P_{l} ψ_{l} ，

则离散问题 (4.6) 在对基函数作检验后，成为关于系数矢量 $U \in R^{n}$ 与 $P \in R^{m}$ 的线性系统。用如下矩阵与矢量元素

\begin{matrix} (4.9) & A_{i j} = a (φ_{j}, φ_{i}) ， B_{k i} = b (φ_{i}, ψ_{k}) ， F_{i} = ℓ (φ_{i}) ， \end{matrix}

(4.6) 的两行即为 $A U + B^{T} P = F$ 与 $B U = 0$ ，也就是

\begin{matrix} (4.10) & (\begin{matrix} A & B^{T} \\ B & 0 \end{matrix}) (\begin{matrix} U \\ P \end{matrix}) = (\begin{matrix} F \\ 0 \end{matrix}) . \end{matrix}

矩阵 $A$ 是 $a$ 的 Gram 矩阵，由强制性知它对称正定；对粘性项的 Laplace 形式而言，它恰是矢量 Laplace 刚度矩阵。矩形矩阵 $B$ 表示离散散度，而离散 inf–sup 条件 (4.7) 恰好表明 $B$ 具有满行秩、且带一个关于 $h$ 一致的下界——这正是保证 (4.10) 可解的性质。

系统 (4.10) 就是第 6 章的对称鞍点系统 (6.1) 在 $C = 0$ 时的情形。对一个本不稳定的单元对作稳定化离散，会添加一个对称半正定的压强–压强项，把零块替换为 $- C$ （ $C \neq 0$ ），从而给出该章处理的一般形式 [9]。它的高效迭代求解——代数形式的 inf–sup 理论、Krylov 子空间方法与块预条件子——是下一（也是最后一）章的主题。

完整非线性不可压缩 Navier–Stokes 方程的有限元处理——对流项及其诱导的三线性形式、inf–sup 稳定单元、Picard 与 Newton 线性化，以及高 Reynolds 数下所需的稳定化——是第 5 章的主题。其每一线性化步骤同样产生形如 (4.10) 的鞍点系统，只是此时一般不再对称，其高效求解则在第 6 章展开。

本文采用署名-非商业性使用-相同方式共享 4.0 国际许可协议，转载请注明出处。

Introduction to CFD on the Basis of FEM