有限体积法介绍

有限差分法（FDM）与有限体积法（FVM）的根本区别

• 有限差分法 (FDM)：直接在网格节点上将偏微分方程中的导数项替换为差分近似。它关注的是节点上的值，并试图近似该点处导数的精确值。
• 有限体积法 (FVM)：首先将计算区域划分为一系列不重叠的控制体积（或单元）。然后，对每个控制体积内的偏微分方程进行积分。这使得方程从微分形式转换为了积分形式，表达了物理量在控制体积内的守恒。FVM关注的是控制体积内的平均值以及通过控制体表面的通量。

有限体积法的控制方程推导

考虑一个通用的控制体积 \(\Omega\) (在2D中，通常是一个小矩形区域 \(\Delta x \Delta y\))，其边界为 \(\partial\Omega\)。

1. x-动量方程

FDM形式:

\[\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)\]

为了方便积分，我们可以将对流项写成散度形式（利用连续性方程 \(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0\)）：

\[u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \frac{\partial (uu)}{\partial x} + \frac{\partial (vu)}{\partial y} - u(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}) = \frac{\partial (uu)}{\partial x} + \frac{\partial (vu)}{\partial y}\]

所以，x-动量方程可以写为： \(\frac{\partial u}{\partial t} + \frac{\partial (uu)}{\partial x} + \frac{\partial (vu)}{\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)\) 或者更紧凑地，使用矢量记号 \(\mathbf{u} = (u,v)\): \(\frac{\partial u}{\partial t} + \nabla \cdot (u\mathbf{u}) = -\frac{\partial p}{\partial x} + \frac{1}{Re} \nabla \cdot (\nabla u)\)

积分形式 (FVM): 对上述方程在控制体积 \(\Omega\) 上积分： \(\int_{\Omega} \frac{\partial u}{\partial t} dV + \int_{\Omega} \nabla \cdot (u\mathbf{u}) dV = \int_{\Omega} \left(-\frac{\partial p}{\partial x}\right) dV + \int_{\Omega} \frac{1}{Re} \nabla \cdot (\nabla u) dV\) 应用高斯散度定理 (\(\int_{\Omega} \nabla \cdot \mathbf{F} dV = \oint_{\partial\Omega} \mathbf{F} \cdot \mathbf{n} dS\)): \(\frac{d}{dt} \int_{\Omega} u dV + \oint_{\partial\Omega} (u\mathbf{u}) \cdot \mathbf{n} dS = \oint_{\partial\Omega} \left(-\frac{p}{\rho} \mathbf{i} \right) \cdot \mathbf{n} dS + \oint_{\partial\Omega} \frac{1}{Re} (\nabla u) \cdot \mathbf{n} dS\)

各项解释:

1.1 瞬态项:

\[\frac{d}{dt} \int_{\Omega} u dV\]

表示控制体积内 \(u\)-动量随时间的变化率。

离散化：

\[\frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} \cdot \Delta x \Delta y\]

其中\(u_{i,j}^n\)为上一步的单元中心速度，\(\Delta x, \Delta y\)为网格尺寸。二维中，\(dV=ds=\Delta x\Delta y\)。

1.2 对流项:

\[\oint_{\partial\Omega} (u\mathbf{u}) \cdot \mathbf{n} dS\]

表示通过控制体边界的 \(u\)-动量的对流通量。\(\mathbf{n}\) 是边界的外法线向量，\(dS\) 是边界的微元面积。

离散化： \(\sum_{f} (u_f \mathbf{u}_f \cdot \mathbf{n}_f) S_f\)

其中 \(f\) 代表控制体的各个面（东、西、南、北），\(u_f\) 和 \(\mathbf{u}_f\) 是面 \(f\) 上的速度值（通常通过插值得到），\(S_f\)是面 \(f\) 的面积。

在代码中，面通量采用中心插值：

• 东面通量（$F_e$）：\(F_e = 0.5 \times (u_{i,j} + u_{i,j+1}) \times \Delta y\)
• 西面通量（$F_w$）：\(F_w = 0.5 \times (u_{i,j} + u_{i,j-1}) \times \Delta y\)
• 北面通量（$F_n$）：\(F_n = 0.5 \times (v_{i,j} + v_{i+1,j}) \times \Delta x\)
• 南面通量（$F_s$）：\(F_s = 0.5 \times (v_{i,j} + v_{i-1,j}) \times \Delta x\)

对流项净通量离散为： \(\text{conv\_term} = \frac{F_e \cdot \phi_{E} - F_w \cdot \phi_{W} + F_n \cdot \phi_{N} - F_s \cdot \phi_{S}}{\Delta x \Delta y}\) 其中$\phi$为待对流的变量（\(u\)或\(v\)）。

在俺的代码中，F_e = 0.5 * (u_prev[j, i] + u_prev[j, i + 1]) * dy 是面上的法向速度分量，而 conv_term 中的 F_e * u_prev[j, i + 1] 是通量项。

1.3 压力梯度项:

\(\int_{\Omega} \left(-\frac{\partial p}{\partial x}\right)\)

表示作用在控制体积上的由于压力梯度产生的 \(x\) 方向的力。 离散化：通常近似为 \((P_W - P_E) \Delta y\)，其中 \(P_W\)和 \(P_E\) 是控制体西面和东面上的压力值（或相邻单元中心的压力）。

1.4 扩散项:

\(\oint_{\partial\Omega} \frac{1}{Re} (\nabla u) \cdot \mathbf{n} dS\)

表示通过控制体边界的 \(u\)-动量的分子扩散通量。 离散化：

\[\sum_{f} \frac{1}{Re} (\nabla u)_f \cdot \mathbf{n}_f S_f\]

扩散项同样通过面通量表达：

• 东面扩散通量：\(D_e = \frac{u_{i,j+1} - u_{i,j}}{\Delta x} \Delta y\)
• 西面扩散通量：\(D_w = \frac{u_{i,j} - u_{i,j-1}}{\Delta x} \Delta y\)
• 北面扩散通量：\(D_n = \frac{u_{i+1,j} - u_{i,j}}{\Delta y} \Delta x\)
• 南面扩散通量：\(D_s = \frac{u_{i,j} - u_{i-1,j}}{\Delta y} \Delta x\)

扩散项净通量离散为： \(\text{diff\_term} = \frac{1}{Re} \cdot \frac{D_e - D_w + D_n - D_s}{\Delta x \Delta y}\)

例如，对于东面 (e)：\(\frac{1}{Re} \left(\frac{\partial u}{\partial x}\right)_e \Delta y\)。梯度 \((\frac{\partial u}{\partial x})_e\) 通常用相邻单元中心值的差分近似，如 \(\frac{U_{E} - U_{P}}{\Delta x_{PE}}\)。 俺的代码中 D_e = (u_prev[j, i + 1] - u_prev[j, i]) / dx * dy 就是东面扩散通量乘以雷诺数倒数前的部分。

2. y-动量方程

FDM形式:

\[\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -\frac{\partial p}{\partial y} + \frac{1}{Re}\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right)\]

同样可以写为：

\[\frac{\partial v}{\partial t} + \nabla \cdot (v\mathbf{u}) = -\frac{\partial p}{\partial y} + \frac{1}{Re} \nabla \cdot (\nabla v)\]

积分形式 (FVM): \(\frac{d}{dt} \int_{\Omega} v dV + \oint_{\partial\Omega} (v\mathbf{u}) \cdot \mathbf{n} dS = \int_{\Omega} \left(-\frac{\partial p}{\partial y}\right) dV + \oint_{\partial\Omega} \frac{1}{Re} (\nabla v) \cdot \mathbf{n} dS\)

各项的离散化与x-动量方程类似，只是将 \(u\) 替换为 \(v\)，压力梯度项变为 \(y\) 方向。

3. 连续性方程

FDM形式:

\[\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0\]

或者 \(\nabla \cdot \mathbf{u} = 0\)。

积分形式 (FVM):

\[\int_{\Omega} (\nabla \cdot \mathbf{u}) dV = 0\]

应用高斯散度定理: \(\oint_{\partial\Omega} \mathbf{u} \cdot \mathbf{n} dS = 0\)

表示通过控制体边界的净体积流率（因为密度为常数，也代表净质量流率）。对于不可压缩流，此项必须为零，即流入控制体的质量等于流出的质量。

离散化：

\[\sum_{f} (\mathbf{u}_f \cdot \mathbf{n}_f) S_f = 0\]

例如，对于一个2D矩形控制体：\((u_e \Delta y) - (u_w \Delta y) + (v_n \Delta x) - (v_s \Delta x) = 0\)。 其中 \(u_e, u_w, v_n, v_s\)分别是东、西、北、南面中心处的法向速度分量。

---

上述探讨了有限体积法在x\y动量与连续性方程的形式，下面来看看程序实现：

主要思路与有限差分一致，不同在于控制方程的有限差分离散，这里将变为控制方程的有限体积通量法，即计算对流项与扩散项有所不同，后面求解压力泊松方程的算法一样(雅可比、高斯-赛尔、逐次松弛因子迭代)。

中间速度计算

预测步忽略压力项(即中间速度求解)离散方程为：

\[u_{i,j}^{*} = u_{i,j}^n + \Delta t \left[ -\text{conv\_term}（对流项） + \text{diff\_term}（扩散项） \right]\]

$v$分量同理。

可能你发现上面不是说：

\[\frac{d}{dt} \int_{\Omega} u dV = \frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} \cdot \Delta x \Delta y\]

那\(\Delta x \Delta y\)哪里去了？这是因为对流项净通量和扩散项净通量已经除以了\(\Delta x \Delta y\)

高亮显示这段话

这也就意味着后面求解压力泊松方程的程序代码可以与有限差分法的完全一样。

FVM控制方程的计算方式

1. 网格生成: 将计算域划分为有限数量的控制体积。
2. 方程积分: 对每个控制体积，将守恒方程（如动量方程、连续性方程）在体积上积分，得到积分形式的方程。
3. 通量离散:
   • 控制体表面上的通量（对流和扩散）通过插值方法，用相邻控制中心的变量值来近似。常用的插值格式有中心差分、迎风格式、混合格式、QUICK格式等。
   • 例如，东面上的值 \(u_e\) 可以是 \(0.5(U_P + U_E)\)（中心差分），或者根据流向选择 \(U_P\) 或 \(U_E\) （迎风格式）。
   • 面上的梯度也用相邻单元中心值的差分来近似。
4. 源项处理: 体积积分的源项（如压力梯度项，如果未通过散度定理转换）通常假设在控制体积内为常数（等于控制中心的值）或线性变化，然后积分。
5. 代数方程组: 对每个控制体积应用上述离散化后，会得到一组关于控制体积中心变量值的代数方程。例如，对于\(u\)-动量方程，会得到一个形如 \(a_P U_P = \sum a_{nb} U_{nb} + S_u\) 的方程，其中 \(a_P\) 是中心系数，\(a_{nb}\) 是邻点系数，\(S_u\) 是源项。
6. 求解: 联立所有控制体积的代数方程，形成一个大型稀疏矩阵方程组，然后使用迭代或直接方法求解，得到每个控制体积中心的变量值（如 \(U_P, V_P, P_P\)）。

总结FVM与FDM在控制方程处理上的不同

特性	有限差分法 (FDM)	有限体积法 (FVM)
基本思想	在节点处用差分近似导数	对控制体积积分守恒方程，计算通过表面的通量
方程形式	直接处理偏微分方程	将偏微分方程转换为积分形式
守恒性	通常不严格保证局部守恒（除非特殊构造）	严格保证每个控制体积的物理量守恒
网格	主要用于结构化网格，非结构网格较复杂	灵活适用于结构化和非结构化网格，以及复杂几何
变量位置	变量定义在网格节点上	变量通常定义在控制体积中心（单元平均值）
离散对象	导数项	通过控制体表面的通量和控制体积内的源项
关注点	节点值的精确性和导数近似的精度	控制体积内物理量的平衡和通过界面的交换