1.3 控制方程的离散

1.3 控制方程的离散

1.3.1 高斯定理

在有限体积法中,离散的核心数学公式为高斯定理:

(1.20.a)V(a)dV=VdSa\int_V(\nabla\cdot\vec a )dV=\oint_{\partial V}d\vec S\cdot\vec a \tag{1.20.a}

(1.20.b)V(ϕ)dV=VdSϕ\int_V(\nabla\phi )dV=\oint_{\partial V}d\vec S\phi \tag{1.20.b}

(1.20.c)V(a)dV=VdSa\int_V(\nabla\vec a )dV=\oint_{\partial V}d\vec S\vec a \tag{1.20.c}

其中VV为空间中的一个体积,V\partial V是控制体的表面。dSd\vec SV\partial V的一个微元,方向为V\partial V外法向。我们最常用的表达式为(1.20.a)。

对于如图1.2所示的一个控制体VPV_P,其(1.20.a)形式的高斯定理为:

(1.21)VP(a)dV=VPdSa\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)右侧转化为几个面上计算结果的求和,即:

(1.22)VPdSa=f(fdSa)\oint_{\partial V_P}d\vec S\cdot\vec a=\sum_f\left(\int_fd\vec S\cdot\vec a\right) \tag{1.22}

在每个小平面上对fdSa\int_fd\vec S\cdot\vec a进行泰勒展开并引入有限体积法的二阶精度限制将式(1.1)式改写为面上的形式带入可得:

(1.23)fdSa=(fdS)af+[fdS(xxf)]:(a)f+=Sfaf\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}

式中Sf=fdS\vec{S_f}=\int_{f} d \vec{S},为一大小等于小平面面积,通过小平面中心指向外侧的向量,我们称之为面向量(如图1.2)。af\vec{a_f}为小平面中心的a\vec a的值。

联立1.22-1.23式,可得:

(1.24)VP(a)dV=f(Sfaf)\int_{V_P}(\nabla\cdot\vec a )dV=\sum_f\left(\vec{S_f}\cdot \vec{a_f}\right) \tag{1.24}

同时引入式(1.3),我们可以得到有限体积法中非常重要的一组等式:

(1.25)VP(a)dV=(a)PVP=f(Sfaf)\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)的标量输运方程,暂时不考虑时间上的积分,在空间上将其写成如下的积分形式:

(1.27)VP(ρϕ)tdV+VP(ρϕu)dV=VP(Γϕ)dV+VPSϕdV\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)进行离散:

  • 对流项:

(1.28)VP(ρϕu)dV=f[Sf(ρϕu)f]=f[Sf(ρuf)ϕf]=f(Fϕf)\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}

上式中F=Sf(ρuf)\mathbf{F}={S_f}\cdot(\rho\vec {u_f})为质量通量(mass flux)。显然,我们还需要知道面上的标量值ϕf\phi_f才可以计算式(1.28),但根据之前的叙述,有限体积法的值都存储在控制体中心,因此我们需要体心值去得到面上的值,这就涉及到了对流项的离散格式,关于这一点,我们将在第三章进行分析。

  • 扩散项:

(1.29)VP(Γϕ)dV=f[Sf(Γϕ)f]=f[SfΓf(ϕ)f]\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}

显然,我们同样需要给定面上的扩散系数Γf\Gamma_f和面上的梯度(ϕ)f\left(\nabla\phi\right)_f才可以计算式(1.29),这就引入了扩散项的离散格式,关于这一点,我们也将在第三章进行讨论。

  • 源项
    源项是一个很广泛的概念,所有不能够写成扩散项,对流项和瞬态项的项都可以归结为源项,因此对于源项我们需要一个非常通用的处理方法。通常情况下,我们对源项离散前需要将源项进行线性化,即:

(1.30)Sϕ=Su+SpϕS_\phi=S_u+S_p\phi \tag{1.30}

考虑到SϕS_\phiϕ\phi的函数,SuS_uSpS_p同样可以是ϕ\phi的函数。结合式(1.3),对源项的离散可以写为:

(1.31)VPSϕdV=SuVp+SpϕpVp\int_{V_P}S_\phi dV=S_uV_p+S_p\phi_pV_p \tag{1.31}

关于源项离散的更多内容,我们将在第四章进行讨论。

  • 瞬态项的离散
    在瞬态问题中,我们还需要进行时间上的积分,根据有限体积法二阶精度的精神,对于时间上的积分我们可仿照式(1.3)推出:

(1.32)tt+tϕdt=ϕtit\int_t^{t+\triangle t}\phi dt=\phi^{t_i}\triangle t \tag{1.32}

其中ϕti\phi^{t_i}为标量ϕ\phittt\sim\triangle t内的一个特征值,其可以是ϕt\phi^{t},也可以ϕt+t\phi^{t+\triangle t},也可以是ϕt\phi^{t}ϕt+t\phi^{t+\triangle t}的函数。不同ϕti\phi^{t_i}对应着不同的时间项离散格式,在这里我们默认采用t+tt+\triangle t时刻的值进行计算。时间离散格式还需要解决另一个问题,那就是如何表示时间项的导数,在这里,我们暂时采用最简单的表达方式:

(1.33)(ρϕ)t=(ρϕ)t+t(ρϕ)tt\frac{\partial(\rho\phi)}{\partial t}=\frac{(\rho\phi)^{t+\triangle t}-(\rho\phi)^{t}}{\triangle t} \tag{1.33}

至此,我们可以对式1.27进行时间上的积分,并进行离散:

(1.34)tt+t[VP(ρϕ)tdV]dt+tt+t[VP(ρϕu)dV]dt=tt+t[VP(Γϕ)dV]dt+tt+t[VPSϕdV]dt(ρϕ)Pt+t(ρϕ)Pttt+[f(Fϕf)]t+tt=[f[SfΓf(ϕ)f]]t+tt+[SuVp+SpϕpVp]t+tt\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}

整理可得:

(1.35)(ρϕ)Pt+t(ρϕ)Ptt+[f(Fϕf)]t+t=[f[SfΓf(ϕ)f]]t+t+[SuVp+SpϕpVp]t+t\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是我们真正求解的方程,其为一个简单的代数方程。主要到其中有两类标量值,一类是位于体心的ϕP\phi_P,一类是位于控制面上的ϕf\phi_f,在实际的计算中,我们通过引入离散格式将面上的值转化为当前控制体体心值ϕP\phi_P与相邻控制体体心值ϕN\phi_N的函数。无论采用何种离散格式,方程1.35总可以被整理为如下的形式:

(1.36)aPϕpn+NaNϕNn=Sa_P\phi_p^n+\sum_Na_N\phi_N^n=S \tag{1.36}

我们首先解释一下上标n,上标n表示new,即新值,表示待求的未知量。同样的,我们之后用上标o表示ordinary,即旧值,表示已知量。在瞬态离散中,旧值是上一时间步的值,新值是当前待求时间步的值。如果我们按照式1.3进行离散,那唯一的旧值来源将是瞬态项中的(ρϕ)Pt(\rho\phi)_P^{t}。但显然我们还有其他的时间离散方案,例如我们对瞬态项采用隐式离散(即对流项的值都设为未知的),对扩散项进行显式离散(即扩散项都用目前已知的值进行计算),旧值的来源则将包括扩散项。对于稳态计算,不存在瞬态项,旧值是迭代计算中上一迭代步的值,新值是本迭代步待求的值。由于旧值是已知的,因此计算结果都归到了式(1.36)右侧的常数项SS中。

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

式(1.36)对每个控制体都可以列一个,对于含有多个控制体的待求解系统来说,多个形如式(1.36)的方程将构成一个方程组。这里我们举一个简单的例子,假设一个一维系统存在三个控制体,从左至右依次记为1,2,3号控制体,1号控制体左侧的边界条件记为01,3号控制体右侧的边界条件记为02(目前可以把边界条件看做一个隐藏的无需求解的控制体),对每个控制体列写形如(1.36)的方程:

{控制体1:A1ϕ01+B1ϕ1当前控制体的值+C1ϕ2相邻控制体的值+0+0=S1控制体2:0+B2ϕ1相邻控制体的值+C2ϕ2当前控制体的值+D2ϕ3相邻控制体的值+0=S2控制体3:0+0+C3ϕ2相邻控制体的值+D3ϕ3相邻控制体的值+E3ϕ02=S3\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.

边界条件值是已知的,归入源项:

(1.37){B1ϕ1+C1ϕ2+0+0=S1A1ϕ010+B2ϕ1+C2ϕ2+D2ϕ3+0=S20+0+C3ϕ2+D3ϕ3=S3E3ϕ02\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}

写成矩阵形式为:

(1.38)[B1C100B2C2D200C3D30][ϕ1ϕ2ϕ3]=[S1A1ϕ01S2S3E3ϕ02]\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}

即:

(1.39)[A][ϕ]=[S][A][\phi]=[S] \tag{1.39}

在实际的计算中,我们通过迭代的方法求解矩阵方程(1.39)以得到我们需要的解,关于矩阵求解的部分我们放到第四章进行分析。对比分析式1.36,式1.37和式1.38可得很容易发现,所有式1.36中的aPa_P最终转化为了系数矩阵中的对角项,所有aNa_N最终转化为了系数矩阵中的非对角项。因此,式1.36除了可以表示一个单独控制体的离散结果,还可以表示整个系统的离散结果。具体来说,式(1.36)相当于将式1.39改写为:

(1.40)[Adiag][ϕ]+[Aoffdiag][ϕ]=[S][A_{diag}][\phi]+[A_{offdiag}][\phi]=[S] \tag{1.40}

在今后的章节中,当我们再次列写形如式(1.36)的离散方程时,我们一般指其矩阵形式(1.40)。

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