1.3 控制方程的离散（OpenFOAM理论笔记系列）

最新推荐文章于 2024-07-11 09:45:31 发布

CloudBird07

最新推荐文章于 2024-07-11 09:45:31 发布

阅读量4.4k

点赞数 7

分类专栏： OpenFOAM理论笔记文章标签： openfoam cfd

本文链接：https://blog.csdn.net/CloudBird07/article/details/107450545

版权

OpenFOAM理论笔记专栏收录该内容

10 篇文章 163 订阅

订阅专栏

1.3 控制方程的离散

1.3.1 高斯定理

在有限体积法中，离散的核心数学公式为高斯定理：
$\int_V(\nabla\cdot\vec a )dV=\oint_{\partial V}d\vec S\cdot\vec a \tag{1.20.a}$
$\int_V(\nabla\phi )dV=\oint_{\partial V}d\vec S\phi \tag{1.20.b}$
$\int_V(\nabla\vec a )dV=\oint_{\partial V}d\vec S\vec a \tag{1.20.c}$
其中 $V$ 为空间中的一个体积， $\partial V$ 是控制体的表面。 $d\vec S$ 是 $\partial V$ 的一个微元，方向为 $\partial V$ 外法向。我们最常用的表达式为(1.20.a)。

对于如图1.2所示的一个控制体 $V_P$ ，其(1.20.a)形式的高斯定理为：
$\int_{V_P}(\nabla\cdot\vec a )dV=\oint_{\partial V_P}d\vec S\cdot\vec a \tag{1.21}$
有限体积法中的控制体表面往往是由几个平面拼接在一起的(如图1.2)，此时我们可以将方程(1.21)右侧转化为几个面上计算结果的求和，即：
$\oint_{\partial V_P}d\vec S\cdot\vec a=\sum_f\left(\int_fd\vec S\cdot\vec a\right) \tag{1.22}$
在每个小平面上对 $\int_fd\vec S\cdot\vec a$ 进行泰勒展开并引入有限体积法的二阶精度限制将式(1.1)式改写为面上的形式带入可得：
$\int_fd\vec S\cdot\vec a=\left(\int_{f} d \vec{S}\right) \cdot \vec{a}_{f}+\left[\int_{f} d \vec{S}\left(\vec{x}-\vec{x}_{f}\right)\right]:(\nabla \vec{a})_{f}+\cdots \\=\vec{S_f}\cdot \vec{a_f} \tag{1.23}$
式中 $\vec{S_f}=\int_{f} d \vec{S}$ ,为一大小等于小平面面积，通过小平面中心指向外侧的向量，我们称之为面向量(如图1.2)。 $\vec{a_f}$ 为小平面中心的 $\vec a$ 的值。

联立1.22-1.23式，可得：
$\int_{V_P}(\nabla\cdot\vec a )dV=\sum_f\left(\vec{S_f}\cdot \vec{a_f}\right) \tag{1.24}$
同时引入式(1.3)，我们可以得到有限体积法中非常重要的一组等式：
$\int_{V_P}(\nabla\cdot\vec a )dV={\left(\nabla\cdot\vec a\right)}_PV_P=\sum_f\left(\vec{S_f}\cdot \vec{a_f}\right) \tag{1.25}$

1.3.2 标量输运方程的离散

对于形如式(1.16)的标量输运方程，暂时不考虑时间上的积分，在空间上将其写成如下的积分形式：
$\int_{V_P}\frac{\partial(\rho\phi)}{\partial t} dV+\int_{V_P}\nabla\cdot(\rho\phi\vec u) dV=\int_{V_P}\nabla\cdot(\Gamma\nabla\phi) dV+\int_{V_P}S_\phi dV \tag{1.27}$
对于对流项和扩散项可以直接使用式(1.24)进行离散：

对流项：
$\int_{V_P}\nabla\cdot(\rho\phi\vec u) dV=\sum_f\left[ {S_f}\cdot\left(\rho\phi\vec u\right)_f\right]=\sum_f\left[ {S_f}\cdot(\rho\vec {u_f})\phi_f\right]\\=\sum_f(\mathbf{F}\phi_f) \tag{1.28}$
上式中 $\mathbf{F}={S_f}\cdot(\rho\vec {u_f})$ 为质量通量(mass flux)。显然，我们还需要知道面上的标量值 $\phi_f$ 才可以计算式(1.28),但根据之前的叙述，有限体积法的值都存储在控制体中心，因此我们需要体心值去得到面上的值，这就涉及到了对流项的离散格式，关于这一点，我们将在第三章进行分析。
扩散项：
$\int_{V_P}\nabla\cdot(\Gamma\nabla\phi) dV=\sum_f\left[ {S_f}\cdot\left(\Gamma\nabla\phi\right)_f\right]=\sum_f\left[ {S_f}\cdot\Gamma_f\left(\nabla\phi\right)_f\right] \tag{1.29}$
显然，我们同样需要给定面上的扩散系数 $\Gamma_f$ 和面上的梯度 $\left(\nabla\phi\right)_f$ 才可以计算式(1.29)，这就引入了扩散项的离散格式，关于这一点，我们也将在第三章进行讨论。
源项
源项是一个很广泛的概念，所有不能够写成扩散项，对流项和瞬态项的项都可以归结为源项，因此对于源项我们需要一个非常通用的处理方法。通常情况下，我们对源项离散前需要将源项进行线性化，即：
$S_\phi=S_u+S_p\phi \tag{1.30}$
考虑到 $S_\phi$ 是 $\phi$ 的函数， $S_u$ 与 $S_p$ 同样可以是 $\phi$ 的函数。结合式(1.3)，对源项的离散可以写为：
$\int_{V_P}S_\phi dV=S_uV_p+S_p\phi_pV_p \tag{1.31}$
关于源项离散的更多内容，我们将在第四章进行讨论。
瞬态项的离散
在瞬态问题中，我们还需要进行时间上的积分，根据有限体积法二阶精度的精神，对于时间上的积分我们可仿照式(1.3)推出：
$\int_t^{t+\triangle t}\phi dt=\phi^{t_i}\triangle t \tag{1.32}$
其中 $\phi^{t_i}$ 为标量 $\phi$ 在 $t\sim\triangle t$ 内的一个特征值，其可以是 $\phi^{t}$ ，也可以 $\phi^{t+\triangle t}$ ，也可以是 $\phi^{t}$ 和 $\phi^{t+\triangle t}$ 的函数。不同 $\phi^{t_i}$ 对应着不同的时间项离散格式，在这里我们默认采用 $t+\triangle t$ 时刻的值进行计算。时间离散格式还需要解决另一个问题，那就是如何表示时间项的导数，在这里，我们暂时采用最简单的表达方式：
$\frac{\partial(\rho\phi)}{\partial t}=\frac{(\rho\phi)^{t+\triangle t}-(\rho\phi)^{t}}{\triangle t} \tag{1.33}$
至此，我们可以对式1.27进行时间上的积分，并进行离散：
$\int_t^{t+\triangle t}\left[\int_{V_P}\frac{\partial(\rho\phi)}{\partial t} dV\right]dt+\int_t^{t+\triangle t}\left[\int_{V_P}\nabla\cdot(\rho\phi\vec u) dV\right]dt\\=\int_t^{t+\triangle t}\left[\int_{V_P}\nabla\cdot(\Gamma\nabla\phi) dV\right]dt+\int_t^{t+\triangle t}\left[\int_{V_P}S_\phi dV\right]dt \\ \\\rightarrow \\ \frac{(\rho\phi)_P^{t+\triangle t}-(\rho\phi)_P^{t}}{\triangle t}\triangle t+\left[\sum_f(\mathbf{F}\phi_f)\right]^{t+\triangle t}\triangle t\\=\\\left[ \sum_f\left[ {S_f}\cdot\Gamma_f\left(\nabla\phi\right)_f\right]\right]^{t+\triangle t}\triangle t+\left[S_uV_p+S_p\phi_pV_p \right]^{t+\triangle t}\triangle t \tag{1.34}$
整理可得：
$\frac{(\rho\phi)_P^{t+\triangle t}-(\rho\phi)_P^{t}}{\triangle t}+\left[\sum_f(\mathbf{F}\phi_f)\right]^{t+\triangle t}\\=\\\left[ \sum_f\left[ {S_f}\cdot\Gamma_f\left(\nabla\phi\right)_f\right]\right]^{t+\triangle t}+\left[S_uV_p+S_p\phi_pV_p \right]^{t+\triangle t} \tag{1.35}$

1.3.3 离散后方程的求解

式1.35是我们真正求解的方程，其为一个简单的代数方程。主要到其中有两类标量值，一类是位于体心的 $\phi_P$ ,一类是位于控制面上的 $\phi_f$ ，在实际的计算中，我们通过引入离散格式将面上的值转化为当前控制体体心值 $\phi_P$ 与相邻控制体体心值 $\phi_N$ 的函数。无论采用何种离散格式，方程1.35总可以被整理为如下的形式：
$a_P\phi_p^n+\sum_Na_N\phi_N^n=S \tag{1.36}$
我们首先解释一下上标n，上标n表示new，即新值，表示待求的未知量。同样的，我们之后用上标o表示ordinary，即旧值，表示已知量。在瞬态离散中，旧值是上一时间步的值，新值是当前待求时间步的值。如果我们按照式1.3进行离散，那唯一的旧值来源将是瞬态项中的 $(\rho\phi)_P^{t}$ 。但显然我们还有其他的时间离散方案，例如我们对瞬态项采用隐式离散(即对流项的值都设为未知的)，对扩散项进行显式离散(即扩散项都用目前已知的值进行计算)，旧值的来源则将包括扩散项。对于稳态计算，不存在瞬态项，旧值是迭代计算中上一迭代步的值，新值是本迭代步待求的值。由于旧值是已知的，因此计算结果都归到了式(1.36)右侧的常数项 $S$ 中。

已经说过，我们通过引入离散格式将面上的值转化为当前控制体体心值 $\phi_P$ 与相邻控制体体心值 $\phi_N$ 的函数。对于方程左侧，第一项表示当前控制体体心的部分，第二项表示相邻控制体体心值的部分，由于相邻的控制体不止一个，因此我们使用加和符号表示所有相邻控制体的值。

式(1.36)对每个控制体都可以列一个，对于含有多个控制体的待求解系统来说，多个形如式(1.36)的方程将构成一个方程组。这里我们举一个简单的例子，假设一个一维系统存在三个控制体，从左至右依次记为1,2,3号控制体，1号控制体左侧的边界条件记为01，3号控制体右侧的边界条件记为02(目前可以把边界条件看做一个隐藏的无需求解的控制体)，对每个控制体列写形如(1.36)的方程：
$\left\{\begin{array}{l} \text{控制体1：}A_{1}\phi_{01}+\underbrace{B_{1}\phi_1}_{\text{当前控制体的值}}+\underbrace{C_{1}\phi_2}_{\text{相邻控制体的值}}+0+0=S_1 \\ \\ \text{控制体2：}0+\underbrace{B_{2}\phi_1}_{\text{相邻控制体的值}}+\underbrace{C_{2}\phi_2}_{\text{当前控制体的值}}+\underbrace{D_{2}\phi_3}_{\text{相邻控制体的值}}+0=S_2 \\ \\ \text{控制体3：}0+0+\underbrace{C_{3}\phi_2}_{\text{相邻控制体的值}}+\underbrace{D_{3}\phi_3}_{\text{相邻控制体的值}}+E_{3}\phi_{02}=S_3 \end{array} \right.$
边界条件值是已知的，归入源项：
$\left\{\begin{array}{l} B_{1}\phi_1+C_{1}\phi_2+0+0=S_1-A_{1}\phi_{01} \\ 0+B_{2}\phi_1+C_{2}\phi_2+D_{2}\phi_3+0=S_2 \\ 0+0+C_{3}\phi_2+D_{3}\phi_3=S_3-E_{3}\phi_{02} \end{array} \right. \tag{1.37}$
写成矩阵形式为：
$\begin{bmatrix} B_1 & C_1 & 0 & 0 \\ B_2 & C_2 & D_2 & 0 \\ 0 & C_3 & D_3 & 0 \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\\phi_3 \end{bmatrix}= \begin{bmatrix} S_1-A_{1}\phi_{01} \\ S_2 \\ S_3-E_{3}\phi_{02} \end{bmatrix} \tag{1.38}$
即：
$[A][\phi]=[S] \tag{1.39}$
在实际的计算中，我们通过迭代的方法求解矩阵方程(1.39)以得到我们需要的解，关于矩阵求解的部分我们放到第四章进行分析。对比分析式1.36，式1.37和式1.38可得很容易发现，所有式1.36中的 $a_P$ 最终转化为了系数矩阵中的对角项，所有 $a_N$ 最终转化为了系数矩阵中的非对角项。因此，式1.36除了可以表示一个单独控制体的离散结果，还可以表示整个系统的离散结果。具体来说，式(1.36)相当于将式1.39改写为：
$[A_{diag}][\phi]+[A_{offdiag}][\phi]=[S] \tag{1.40}$
在今后的章节中，当我们再次列写形如式(1.36)的离散方程时，我们一般指其矩阵形式(1.40)。

系列说明：
接触有限体积法有一段时间了，也看了一些资料，但是有时候总觉得看过一遍之后什么也记不住。老话说得好，眼过千遍不如手过一遍，久而久之我就有了写一些比较像样子的笔记的想法。初学的时候曾经写过一本叫“OpenFOAM编程笔记：单相不可压缩流动”的册子，但当时基础太差（现在基础也不好），错误太多，倒不如推倒重来。
本系列将持续更新，欢迎各位挑错交流，挑错交流可以直接留言也可以联系邮箱cloudbird7@foxmail.com。等到预定内容全部写完后我将集结所有内容为一个独立的文件，开放下载。

CloudBird07

关注

7
点赞
踩
15

收藏

觉得还不错? 一键收藏
7
评论
1.3 控制方程的离散（OpenFOAM理论笔记系列）

1.3 控制方程的离散1.3.1 高斯定理在有限体积法中，离散的核心数学公式为高斯定理：∫V(∇⋅a⃗)dV=∮∂VdS⃗⋅a⃗(1.20.a)\int_V(\nabla\cdot\vec a )dV=\oint_{\partial V}d\vec S\cdot\vec a\tag{1.20.a}∫V(∇⋅a)dV=∮∂VdS⋅a(1.20.a)∫V(∇ϕ)dV=∮∂VdS⃗ϕ(1.20.b)\int_V(\nabla\phi )dV=\oint_{\partial V}d\vec S\
复制链接

扫一扫

专栏目录