FDM-3:压力泊松方程
这是顶盖驱动方腔流的有限差分方法-压力泊松方程概述.
1. 无量纲控制方程（二维情形）
我们从无量纲化后的Navier-Stokes方程和连续性方程开始：
- 无量纲动量方程（以 \(u,v\) 分量表示）：
- \[\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{\partial p}{\partial x} + \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)\]
这里，\(u, v, p, x, y, t\) 均为无量纲变量，\(Re\) 为雷诺数。
- 无量纲连续性方程（不可压缩条件）：
- \[\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0\]
2. 速度修正公式的由来（无量纲形式）
解释为何速度修正采用如下形式：
\[\frac{u^{n+1} - u^*}{\Delta t} = -\frac{\partial p^{n+1}}{\partial x} \quad\]1. 回顾完整的无量纲动量方程
\(x\) 方向动量方程的时间离散形式为：
\[\frac{u^{n+1} - u^n}{\Delta t} = \underbrace{- \left( u^n\frac{\partial u^n}{\partial x} + v^n\frac{\partial u^n}{\partial y} \right) + \frac{1}{Re}\left(\frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2}\right)}_{H_u^n} - \frac{\partial p^{n+1}}{\partial x}\]我们的目标是求解满足连续性方程的 \(u^{n+1}\) 和 \(v^{n+1}\)，而 \(p^{n+1}\) 是保证这一点的关键。
2. “临时速度场” \(u^*\) 的计算
在投影法（分步法）的第一步中，计算临时速度 \(u^*\)（和 \(v^*\)）时，只考虑对流和扩散项（以及旧压力梯度，若有），即式子(1)：
\[\frac{u^* - u^n}{\Delta t} = H_u^n\]计算 \(u^*\) 时未使用第 \(n+1\) 步的压力梯度 \(\frac{\partial p^{n+1}}{\partial x}\)，因此 \(u^*\) 是一个“预测”速度，尚未满足不可压缩条件。
3. 速度修正的推导
将式子(1)带入完整的动量方程可写为：
\[\frac{u^{n+1} - u^n}{\Delta t} = \frac{u^* - u^n}{\Delta t} - \frac{\partial p^{n+1}}{\partial x}\]整理得速度修正公式：
\[\frac{u^{n+1} - u^*}{\Delta t} = -\frac{\partial p^{n+1}}{\partial x}\]\(v\) 分量同理。
4. 物理和数学意义
-
物理意义：该修正表示从临时速度 \(u^*\) 到最终无散度速度 \(u^{n+1}\) 的变化，是由新的压力梯度在时间步长 \(\Delta t\) 内产生的加速度。压力梯度作为力，调整速度场以消除速度场中的“源”或“汇”，即保证不可压缩性。
-
数学意义（投影）：该过程等价于将临时速度场 \(\mathbf{u}^*\) 投影到无散度速度场的子空间。压力梯度场 \(\nabla p^{n+1}\) 与无散度场正交（在 \(L^2\) 意义下），修正步骤去除速度场中梯度场的分量，得到满足连续性的速度场 \(\mathbf{u}^{n+1}\)。
5. 为什么是这种形式？
此形式源于：
- 无量纲动量方程的结构：压力梯度项是其固有部分。
- 时间步进策略（分步法）：先计算不含新压力梯度的临时速度，再用压力梯度修正。
- 满足连续性方程的需求：压力泊松方程的推导基于此修正形式，求解出的压力场正是使速度无散度的关键。
补充：
无量纲压力泊松方程的具体推导
第一步：计算临时速度场
首先对动量方程进行时间离散。设 \(u^n, v^n\) 为第 \(n\) 步的速度，\(u^*, v^*\) 为不考虑第 \(n+1\) 步压力梯度时计算得到的临时速度场。
定义 \(H_u^n\) 为 \(x\) 方向动量方程中除时间导数和压力梯度外的所有项：
\[H_u^n = - \left( u^n\frac{\partial u^n}{\partial x} + v^n\frac{\partial u^n}{\partial y} \right) + \frac{1}{Re}\left(\frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2}\right)\]同理，\(y\) 方向动量方程为：
\[H_v^n = - \left( u^n\frac{\partial v^n}{\partial x} + v^n\frac{\partial v^n}{\partial y} \right) + \frac{1}{Re}\left(\frac{\partial^2 v^n}{\partial x^2} + \frac{\partial^2 v^n}{\partial y^2}\right)\]然后，利用显式时间步进（如前向欧拉）计算临时速度：
\[\frac{u^* - u^n}{\Delta t} = H_u^n\] \[\frac{v^* - v^n}{\Delta t} = H_v^n\]即：
\[u^* = u^n + \Delta t \cdot H_u^n\] \[v^* = v^n + \Delta t \cdot H_v^n\]此时的临时速度场 \((u^*, v^*)\) 一般不满足连续性方程。
第二步：用压力修正速度
时间步 \(n+1\) 的速度场 \((u^{n+1}, v^{n+1})\) 通过压力梯度对临时速度进行修正得到：
\[\frac{u^{n+1} - u^*}{\Delta t} = -\frac{\partial p^{n+1}}{\partial x}\] \[\frac{v^{n+1} - v^*}{\Delta t} = -\frac{\partial p^{n+1}}{\partial y}\]即：
\[u^{n+1} = u^* - \Delta t \frac{\partial p^{n+1}}{\partial x}\] \[v^{n+1} = v^* - \Delta t \frac{\partial p^{n+1}}{\partial y}\]第三步：强制速度场无散度
修正后的速度场 \((u^{n+1}, v^{n+1})\) 必须满足无量纲连续性方程：
\[\frac{\partial u^{n+1}}{\partial x} + \frac{\partial v^{n+1}}{\partial y} = 0\]将第二步的表达式代入：
\[\frac{\partial}{\partial x} \left( u^* - \Delta t \frac{\partial p^{n+1}}{\partial x} \right) + \frac{\partial}{\partial y} \left( v^* - \Delta t \frac{\partial p^{n+1}}{\partial y} \right) = 0\]展开得：
\[\frac{\partial u^*}{\partial x} - \Delta t \frac{\partial^2 p^{n+1}}{\partial x^2} + \frac{\partial v^*}{\partial y} - \Delta t \frac{\partial^2 p^{n+1}}{\partial y^2} = 0\]整理为压力项：
\[\Delta t \left( \frac{\partial^2 p^{n+1}}{\partial x^2} + \frac{\partial^2 p^{n+1}}{\partial y^2} \right) = \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y}\]最终得到无量纲压力泊松方程：
\[\frac{\partial^2 p^{n+1}}{\partial x^2} + \frac{\partial^2 p^{n+1}}{\partial y^2} = \frac{1}{\Delta t} \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right)\]求解压力泊松方程的方法
二维无量纲压力泊松方程：
\[\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2} = \frac{1}{\Delta t} \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right)\]- 右侧的 \(\frac{1}{\Delta t} \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right)\) 是源项，通常在代码中表示为
b[j, i]，代表临时速度场的散度乘以 \(1/\Delta t\)。 - 求解该方程得到的压力场 \(p^{n+1}\)，用于速度修正步骤，确保最终速度场 \((u^{n+1}, v^{n+1})\) 满足连续性方程（无散度）。
采用五点中心差分格式离散拉普拉斯算子（网格间距为 \(dx, dy\)）：
\[\frac{p_{j,i+1} - 2p_{j,i} + p_{j,i-1}}{dx^2} + \frac{p_{j+1,i} - 2p_{j,i} + p_{j-1,i}}{dy^2} = b_{j,i}\]整理求解 \(p_{j,i}\)：
\[p_{j,i} \left( \frac{2}{dx^2} + \frac{2}{dy^2} \right) = \frac{p_{j,i+1} + p_{j,i-1}}{dx^2} + \frac{p_{j+1,i} + p_{j-1,i}}{dy^2} - b_{j,i}\]两边乘以 \(dx^2 dy^2\)：
\[p_{j,i} (2dy^2 + 2dx^2) = dy^2(p_{j,i+1} + p_{j,i-1}) + dx^2(p_{j+1,i} + p_{j-1,i}) - dx^2 dy^2 b_{j,i}\]因此，迭代求解压力的更新公式为：
\[p_{j,i} = \frac{dy^2(p_{j,i+1}^{old} + p_{j,i-1}^{old}) + dx^2(p_{j+1,i}^{old} + p_{j-1,i}^{old}) - dx^2 dy^2 b_{j,i}}{2(dx^2 + dy^2)}\]这与我代码中 p_old_iter 的表达式一致。