1.3 控制方程的离散
1.3.1 高斯定理
在有限体积法中,离散的核心数学公式为高斯定理:
∫V(∇⋅a)dV=∮∂VdS⋅a(1.20.a)
∫V(∇ϕ)dV=∮∂VdSϕ(1.20.b)
∫V(∇a)dV=∮∂VdSa(1.20.c)
其中V为空间中的一个体积,∂V是控制体的表面。dS是∂V的一个微元,方向为∂V外法向。我们最常用的表达式为(1.20.a)。
对于如图1.2所示的一个控制体VP,其(1.20.a)形式的高斯定理为:
∫VP(∇⋅a)dV=∮∂VPdS⋅a(1.21)
有限体积法中的控制体表面往往是由几个平面拼接在一起的(如图1.2),此时我们可以将方程(1.21)右侧转化为几个面上计算结果的求和,即:
∮∂VPdS⋅a=f∑(∫fdS⋅a)(1.22)
在每个小平面上对∫fdS⋅a进行泰勒展开并引入有限体积法的二阶精度限制将式(1.1)式改写为面上的形式带入可得:
∫fdS⋅a=(∫fdS)⋅af+[∫fdS(x−xf)]:(∇a)f+⋯=Sf⋅af(1.23)
式中Sf=∫fdS,为一大小等于小平面面积,通过小平面中心指向外侧的向量,我们称之为面向量(如图1.2)。af为小平面中心的a的值。
联立1.22-1.23式,可得:
∫VP(∇⋅a)dV=f∑(Sf⋅af)(1.24)
同时引入式(1.3),我们可以得到有限体积法中非常重要的一组等式:
∫VP(∇⋅a)dV=(∇⋅a)PVP=f∑(Sf⋅af)(1.25)
1.3.2 标量输运方程的离散
对于形如式(1.16)的标量输运方程,暂时不考虑时间上的积分,在空间上将其写成如下的积分形式:
∫VP∂t∂(ρϕ)dV+∫VP∇⋅(ρϕu)dV=∫VP∇⋅(Γ∇ϕ)dV+∫VPSϕdV(1.27)
对于对流项和扩散项可以直接使用式(1.24)进行离散:
∫VP∇⋅(ρϕu)dV=f∑[Sf⋅(ρϕu)f]=f∑[Sf⋅(ρuf)ϕf]=f∑(Fϕf)(1.28)
上式中F=Sf⋅(ρuf)为质量通量(mass flux)。显然,我们还需要知道面上的标量值ϕf才可以计算式(1.28),但根据之前的叙述,有限体积法的值都存储在控制体中心,因此我们需要体心值去得到面上的值,这就涉及到了对流项的离散格式,关于这一点,我们将在第三章进行分析。
∫VP∇⋅(Γ∇ϕ)dV=f∑[Sf⋅(Γ∇ϕ)f]=f∑[Sf⋅Γf(∇ϕ)f](1.29)
显然,我们同样需要给定面上的扩散系数Γf和面上的梯度(∇ϕ)f才可以计算式(1.29),这就引入了扩散项的离散格式,关于这一点,我们也将在第三章进行讨论。
- 源项
源项是一个很广泛的概念,所有不能够写成扩散项,对流项和瞬态项的项都可以归结为源项,因此对于源项我们需要一个非常通用的处理方法。通常情况下,我们对源项离散前需要将源项进行线性化,即:
Sϕ=Su+Spϕ(1.30)
考虑到Sϕ是ϕ的函数,Su与Sp同样可以是ϕ的函数。结合式(1.3),对源项的离散可以写为:
∫VPSϕdV=SuVp+SpϕpVp(1.31)
关于源项离散的更多内容,我们将在第四章进行讨论。
- 瞬态项的离散
在瞬态问题中,我们还需要进行时间上的积分,根据有限体积法二阶精度的精神,对于时间上的积分我们可仿照式(1.3)推出:
∫tt+△tϕdt=ϕti△t(1.32)
其中ϕti为标量ϕ在t∼△t内的一个特征值,其可以是ϕt,也可以ϕt+△t,也可以是ϕt和ϕt+△t的函数。不同ϕti对应着不同的时间项离散格式,在这里我们默认采用t+△t时刻的值进行计算。时间离散格式还需要解决另一个问题,那就是如何表示时间项的导数,在这里,我们暂时采用最简单的表达方式:
∂t∂(ρϕ)=△t(ρϕ)t+△t−(ρϕ)t(1.33)
至此,我们可以对式1.27进行时间上的积分,并进行离散:
∫tt+△t[∫VP∂t∂(ρϕ)dV]dt+∫tt+△t[∫VP∇⋅(ρϕu)dV]dt=∫tt+△t[∫VP∇⋅(Γ∇ϕ)dV]dt+∫tt+△t[∫VPSϕdV]dt→△t(ρϕ)Pt+△t−(ρϕ)Pt△t+⎣⎡f∑(Fϕf)⎦⎤t+△t△t=⎣⎡f∑[Sf⋅Γf(∇ϕ)f]⎦⎤t+△t△t+[SuVp+SpϕpVp]t+△t△t(1.34)
整理可得:
△t(ρϕ)Pt+△t−(ρϕ)Pt+⎣⎡f∑(Fϕf)⎦⎤t+△t=⎣⎡f∑[Sf⋅Γf(∇ϕ)f]⎦⎤t+△t+[SuVp+SpϕpVp]t+△t(1.35)
1.3.3 离散后方程的求解
式1.35是我们真正求解的方程,其为一个简单的代数方程。主要到其中有两类标量值,一类是位于体心的ϕP,一类是位于控制面上的ϕf,在实际的计算中,我们通过引入离散格式将面上的值转化为当前控制体体心值ϕP与相邻控制体体心值ϕN的函数。无论采用何种离散格式,方程1.35总可以被整理为如下的形式:
aPϕpn+N∑aNϕNn=S(1.36)
我们首先解释一下上标n,上标n表示new,即新值,表示待求的未知量。同样的,我们之后用上标o表示ordinary,即旧值,表示已知量。在瞬态离散中,旧值是上一时间步的值,新值是当前待求时间步的值。如果我们按照式1.3进行离散,那唯一的旧值来源将是瞬态项中的(ρϕ)Pt。但显然我们还有其他的时间离散方案,例如我们对瞬态项采用隐式离散(即对流项的值都设为未知的),对扩散项进行显式离散(即扩散项都用目前已知的值进行计算),旧值的来源则将包括扩散项。对于稳态计算,不存在瞬态项,旧值是迭代计算中上一迭代步的值,新值是本迭代步待求的值。由于旧值是已知的,因此计算结果都归到了式(1.36)右侧的常数项S中。
已经说过,我们通过引入离散格式将面上的值转化为当前控制体体心值ϕP与相邻控制体体心值ϕ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
边界条件值是已知的,归入源项:
⎩⎨⎧B1ϕ1+C1ϕ2+0+0=S1−A1ϕ010+B2ϕ1+C2ϕ2+D2ϕ3+0=S20+0+C3ϕ2+D3ϕ3=S3−E3ϕ02(1.37)
写成矩阵形式为:
⎣⎡B1B20C1C2C30D2D3000⎦⎤⎣⎡ϕ1ϕ2ϕ3⎦⎤=⎣⎡S1−A1ϕ01S2S3−E3ϕ02⎦⎤(1.38)
即:
[A][ϕ]=[S](1.39)
在实际的计算中,我们通过迭代的方法求解矩阵方程(1.39)以得到我们需要的解,关于矩阵求解的部分我们放到第四章进行分析。对比分析式1.36,式1.37和式1.38可得很容易发现,所有式1.36中的aP最终转化为了系数矩阵中的对角项,所有aN最终转化为了系数矩阵中的非对角项。因此,式1.36除了可以表示一个单独控制体的离散结果,还可以表示整个系统的离散结果。具体来说,式(1.36)相当于将式1.39改写为:
[Adiag][ϕ]+[Aoffdiag][ϕ]=[S](1.40)
在今后的章节中,当我们再次列写形如式(1.36)的离散方程时,我们一般指其矩阵形式(1.40)。
系列说明:
接触有限体积法有一段时间了,也看了一些资料,但是有时候总觉得看过一遍之后什么也记不住。老话说得好,眼过千遍不如手过一遍,久而久之我就有了写一些比较像样子的笔记的想法。初学的时候曾经写过一本叫“OpenFOAM编程笔记:单相不可压缩流动”的册子,但当时基础太差(现在基础也不好),错误太多,倒不如推倒重来。
本系列将持续更新,欢迎各位挑错交流,挑错交流可以直接留言也可以联系邮箱cloudbird7@foxmail.com。等到预定内容全部写完后我将集结所有内容为一个独立的文件,开放下载。