控制方程与算法
首先有稳态的N-S方程:
(1)\[
\nabla\cdot\bfU=0,
\]
(2)\[
\nabla \cdot (\mathbf{U}\mathbf{U})-\nabla \cdot(\nu \nabla \mathbf{U})=-\nabla p,
\]
首先,需要对方程(2)做体积积分并隐性离散:
(3)\[
\int{\nabla \cdot \left( {\mathbf{U}\mathbf{U}} \right)\mathrm{d}V= \int {\mathbf{U}\mathbf{U}\mathrm{d}\bfS} = } \sum {{{\left( {{\mathbf{U}^n}{\mathbf{U}^*}} \right)}_f}} \bfS_f = \sum {\left( {F_f^n \mathbf{U}_f^*} \right)} ,
\]
(4)\[
{\int \nabla \cdot \left( {\nu \nabla \mathbf{U}} \right)\mathrm{d}V = \int {\nu \nabla \mathbf{U}\cdot\mathrm{d}\bfS = \sum {\left( {\nu \nabla {\mathbf{U}^*}} \right)} } _f}\cdot\bfS_f = \sum { \nu (\nabla\mathbf{U}^*)_f\cdot\bfS_f } ,
\]
(5)\[
\int \nabla p \mathrm{d}V=\int p \mathrm{d}\bfS=\sum \left(p_f^n\bfS_f\right),
\]
其中上标\(^n\)表示为当前迭代步(已知),上标\(^*\)表示预测步(待求),下标\(_f\)表示网格单元面上的值,\(V_\rP\)表示网格单元体积,\(\bfS\)表示网格单元的各个面的面矢量,\(F_f\)为通量,\(\nu\)为动力粘度。依据方程(3)-(5)有:
(6)\[
\sum F_f^n \mathbf{U}_f^* - \sum { \nu (\nabla\mathbf{U}^*)_f\cdot\bfS_f } = \sum p_f^n\bfS_f.
\]
其中的通量\(F_f^n\)采用当前已知时间步的速度\(\mathbf{U}^n\)来计算,同时保留预测速度\(\mathbf{U}^*\)为未知量。
方程(6)中\((\nabla\mathbf{U}^*)_f\cdot\bfS_f\)需要进一步处理,其中\(\nabla\mathbf{U}^*\)为定义在体心的量,即网格体心的速度梯度。\((\nabla\mathbf{U}^*)_f\)为定义在面心的量,即网格面心的速度梯度。为使用紧致基架点防止数值震荡,\((\nabla\mathbf{U}^*)_f\cdot\bfS_f\)通常表示为面法向梯度与面矢量的模的乘积:
(7)\[
(\nabla\mathbf{U}^*)_f \cdot\mathbf{S}_f=
\left((\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} \right ) \cdot|\mathbf{S}_f|
\]
其中\((\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} \)表示面法向梯度,有
(8)\[
(\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|}
= \frac{{\mathbf{U}_\mathrm{N}^* - \mathbf{U}_\mathrm{P}^*}}{{\left| \mathbf{d} \right|}},
\]
其中\(\mathbf{d}\)表示网格单元体心之间矢量差(距离矢量),下标\(_\mathrm{N}\)表示相邻网格单元的速度,下标\(_\mathrm{P}\)表示当前网格单元的速度。在这里需要注意的是,方程(8)只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性。
方程(6)中\(\mathbf{U}_f^*\)被定义为面上的预测速度,然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度,在此步可以引入各种插值格式。假设使用中心线性格式:
(9)\[
\mathbf{U}_f^* = \frac{{\mathbf{U}_\mathrm{P}^* + \mathbf{U}_\mathrm{N}^*}}{2}.
\]
(10)\[
p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}.
\]
将方程(7),(9)代入到方程(6)有
(11)\[
\sum {F_f^n \frac{{\mathbf{U}_\mathrm{P}^* + \mathbf{U}_\mathrm{N}^*}}{2}} = \sum {\nu \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^* - \mathbf{U}_\mathrm{P}^*}}{{\left| \mathbf{d} \right|}}} -\sum p_f^n\bfS_f,
\]
方程两边同时除以网格体积,并整理有:
(12)\[\begin{split}
\frac{1}{V_\rP}\left( { \sum {\frac{{F_f^n}}{2} + \sum {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} } } \right)\mathbf{U}_\mathrm{P}^*
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;
\\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
= - \frac{1}{V_\rP}\sum {\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} \right)\mathbf{U}_\mathrm{N}^*} - \frac{1}{V_\rP}\sum p_f^n\bfS_f.
\end{split}\]
将上式简写为
(13)\[
{A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}^*{\rm{ + }}\sum {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}^*} = - \frac{1}{V_\rP}\sum p_f^n\bfS_f,
\]
其中
(14)\[
A_\mathrm{P}^n= \frac{1}{V_\rP}\left(\sum {\frac{{F_f^n}}{2}} + \sum {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}}\right) ,
\]
(15)\[
A_\mathrm{N}^n= \frac{1}{V_\rP}\left({\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} \right),
\]
求解方程(13)即可获得预测速度\(\mathbf{U}^*_\mathrm{P}\)。方程(13)即为OpenFOAM中的速度方程。也即动量预测步。
注意:
OpenFOAM中的动量预测步主要是为了组建稀疏线性系统,如方程(13)。目的是获得矩阵系数\(A\)。因此并不强制进行求解,如果求解的话,可获得\(\bfU^*\),不求解的话,即为\(\bfU^n\)。
N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度,但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来,其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下,对连续性方程\(\nabla\cdot\bfU=0\)进行离散后的形式为
(16)\[
\sum \mathbf{U}_{\mathrm{P},f}^{n+1} \cdot \bfS_f=0.
\]
如果能用压力表示方程(16)中的\(\mathbf{U}_{\mathrm{P},f}^{n+1}\),是不是就可以得出压力泊松方程了呢?答案是肯定的。首先,方程(13)中的压力项不建议进行离散,这是因为方程右侧的梯度项都有可能引起数值的震荡(例如双流体模型中的湍流分散力也需要特殊处理)。同时,参考方程(13),
在收敛情况下(时间步\(t=n+1\))为:
(17)\[
{A_\mathrm{P}^{n+1}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {A_\mathrm{N}^{n+1}\mathbf{U}_\mathrm{N}^{n+1}} = -\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f,
\]
由于存在未知的\(A_\mathrm{P}^{n+1},A_\mathrm{N}^{n+1}\)。因为后续将引入迭代算法,因此可以\(A_\mathrm{P}^{n+1},A_\mathrm{N}^{n+1}\)存在一定的滞后,即\(A_\mathrm{P}^{n+1}\approx A_\mathrm{P}^{n},A_\mathrm{N}^{n+1}\approx A_\mathrm{N}^{n}\):
(18)\[
{A_\mathrm{P}^{n}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {A_\mathrm{N}^{n}\mathbf{U}_\mathrm{N}^{n+1}} = -\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f,
\]
注意:
在这里未调用任何略去临点的假定,方程(18)是严格成立的。
若将方程(18)减去(13)有:
(19)\[
{A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}'+\sum {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}'} = -\frac{1}{V_\rP}\sum p_f '\bfS_f,
\]
其中\(\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U}_\mathrm{P}^{*}\),\(\mathbf{U}_\mathrm{N}'=\mathbf{U}_\mathrm{N}^{n+1}-\mathbf{U}_\mathrm{N}^{*}\),\(p_f'=p_f^{n+1}-p_f^{n}\)。既压力的修正量和当前网格点的速度修正量以及相邻网格点的速度修正量有关。对方程(18)移项有
(20)\[
\mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{A_\mathrm{P}^n}\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f.
\]
其中我们定义
(21)\[
\mathbf{HbyA}^{n+1}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}^n}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{n}} } \right).
\]
面上的速度可以通过下面的公式获得:
(22)\[
\mathbf{U}_{\mathrm{P},f}^{n+1} = \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A^n_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f\right)_f.
\]
将方程(22)代入到方程(16)有:
(23)\[
\sum \left( \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A^n_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f\right)_f \right) \cdot \bfS_f =0,
\]
(24)\[
\sum {\mathbf{HbyA}_f^{n+1}} \cdot \bfS_f = \sum {\frac{1}{{A^n_{\mathrm{P},f}}}\left(\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f\right)_f} \cdot \bfS_f,
\]
方程(24)即下述方程的离散形式:
(25)\[
\nabla \cdot \mathbf{HbyA}^{n+1} = \nabla \cdot\left(\frac{1}{A^n_{\mathrm{P}}} \nabla p^{n+1}\right).
\]
方程(25)即为压力泊松方程。若有\(\mathbf{HbyA}^{n+1}\),则可求得收敛的压力。
注意:
这里\(\mathbf{HbyA}^{n+1}\)是未知的,因此方程(25)一个方程存在两个未知量,不可解。因此需要引入下面的迭代算法。
SIMPLE算法
可以看出,最后求解的速度\(\mathbf{U}_\mathrm{P}^{n+1}\)和压力\(p^{n+1}\)应该符合方程(13)和(25),分别对应动量方程和连续性方程。然而目前,我们只有通过动量预测步骤求出来的\(\mathbf{U}_\mathrm{P}^*\)。\(\mathbf{HbyA}_\mathrm{P}^{n+1}\)和\(p^{n+1}\)均为未知的。方程(25)不能求解。SIMPLE迭代算法则可以解决这个问题。依据SIMPLE算法,首先引入略去临点影响的假定,方程(19)可以写为:
(26)\[
{A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}' = -\frac{1}{V_\rP}\sum p_f '\bfS_f,
\]
将方程(26)和(13)相加有
(27)\[
{A^n_\mathrm{P}}\mathbf{U}_\mathrm{P}^{**}{\rm{ + }}\sum {A^n_\mathrm{N}\mathbf{U}_\mathrm{N}^*} = -\frac{1}{V_\rP} \sum p_f^*\bfS_f,
\]
其中\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^*+\mathbf{U}_\mathrm{P}'\),\(p_f^{*}=p_f^n+p_f'\)。
注意:
如果不进行动量预测步计算,则\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^n+\mathbf{U}_\mathrm{P}'\)。
相比方程(18),方程(27)由于忽略了临点的影响并不是精准的。对方程(27)进行移项有
(28)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}}\frac{1}{V_\rP} \sum p_f^*\bfS_f.
\]
(29)\[
\mathbf{HbyA}^{*}_\mathrm{P} = \frac{1}{{{A^n_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{*}} } \right).
\]
(30)\[
\mathbf{U}_{\mathrm{P},f}^{**} = \mathbf{HbyA}_f^{*} - \frac{1}{{{A^n_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP} \sum p_f^*\bfS_f\right)_f.
\]
方程(30)中可参考方程(23)-(25)的步骤组建压力泊松方程,即
(31)\[
\nabla \cdot \mathbf{HbyA}^{*} = \nabla \cdot\left(\frac{1}{A} \nabla p^{*}\right).
\]
注意:
如果不进行动量预测步计算,则\(\mathbf{HbyA}^{*}=\mathbf{HbyA}^{n}\)。
对方程(31)求解后有\(p^{*}\)。注将\(p^{*}\)与\(\mathbf{HbyA}^{*}\)回代到方程(28)的时候,有更新后的\(\mathbf{U}_{\mathrm{P}}^{**}\)。此处的\(\mathbf{U}_{\mathrm{P}}^{**}\)满足连续性方程。同时,\(\mathbf{U}_{\mathrm{P}}^{**}\)与\(p^*\)近似满足动量方程(28)。之所以近似,是因为方程(28)在推导的过程中,调用了忽略临点假定。
注意:
如果不进行动量预测步计算,则完成一次迭代后获得\(\bfU^{*}\)与\(p^{*}\)。
至此,原变量\(\bfU^n\)与\(p^n\)更新为\(\bfU^{**}\)与\(p^{*}\)。此处即完成一次迭代。我们强调,\(\bfU^{**}\)与\(p^{*}\)近似满足动量方程(28),理论上满足连续性方程。在这里,需要强调以下迭代算法的重要性。
方程(27)略掉了临点的影响,因此只能近似满足动量方程。但最终收敛的情况下,修正量会变得非常小,因此方程(27)同样会满足。
在实际计算过程中,尤其是最开始的迭代时,压力修正量过大,会导致难以收敛。因此,在\(p^n\)变化至\(p^*\)的过程中,要减少这种非常大的突变,要对压力进行亚松弛。这就导致方程(31)求解后的压力\(p\),为介于\(p^n\)变化至\(p^*\)之间的某个值。
同时,方程(31)的求解并不会将残差将至最低,也就是方程(31)求解后左右两边并不是严格相等,这也就意味着存在较大的连续性误差,\(\bfU^{**}\)并不严格的满足连续性方程。
因此,通过算例的残差、松弛必要设置,导致\(\bfU^{**}\)与\(p^{*}\)不能作为最终收敛的结果。因此需要将\(\bfU^{**}\)与\(p^{*}\)作为下一个迭代步的初始值,继续进行缓慢稳定的推进。因此,SIMPLE算法中的的迭代过程可以表示为下面几个步骤:
注意:
在这里需要强调的是,因为本身属于迭代算法,因此在求解动量方程(28)以及连续性方程(31)的时候,不需要将残差设置的非常低。
SIMPLEC算法
SIMPLEC算法认为SIMPLE的算法太过激进。这种略去邻点的行为会使得压力修正值\(p'\)过大,在某些情况下会引起发散。虽然这可以通过亚松弛来稳定计算流程,但是会减慢收敛速度。SIMPLEC算法认为任意网格点的速度修正可以近似认为是相邻网格点速度修正的平均:
(32)\[
\bfU_\rP'\approx\frac{\sum A^n_\rN\bfU_\rN'}{\sum A^n_\rN} \:\rightarrow \: \bfU_\rP'\sum A^n_\rN\approx\sum A^n_\rN\bfU_\rN'
\]
将方程(32)代入到(19)有:
(33)\[
A^n_\rP\bfU_\rP'=-\bfU_\rP'\sum A^n_\rN-\frac{1}{V_\rP} \sum p_f'\bfS_f
\]
即
(34)\[
\bfU_\rP'=-\frac{1}{A^n_\rP+\sum A^n_\rN}\frac{1}{V_\rP} \sum p_f'\bfS_f
\]
相对比于SIMPLE算法中的方程(26),对于固定的\(\bfU_\rP'\),方程(34)右侧压力修正的乘积项\(\frac{1}{A_\rP+\sum A^n_\rN}\)必然大于\(\frac{1}{A^n_\rP}\)(有限体积法离散的\(A^n_\rN<0\)),进而\(p_f'\)的值不会特别大,因此SIMPLEC算法并不需要对压力进行激进的低松弛修正。将方程(34)和方程(13)相加有:
(35)\[
{A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{**}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}} = -\frac{1}{V_\rP} \sum p_\rP^n\bfS_f- \frac{A_\mathrm{P}}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\frac{1}{V_\rP} \sum p_f'\bfS_f,
\]
既
(36)\[ \begin{align}\begin{aligned}
\begin{split}
\mathbf{U}_\mathrm{P}^{**} =
-\frac{\sum {A^n_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A^n_\mathrm{P}}} -\frac{1}{A^n_\mathrm{P}}\frac{1}{V_\rP} \sum p_\rP^n\bfS_f- \frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\frac{1}{V_\rP} \sum p_f'\bfS_f\\\begin{split}\\\\
=-\frac{\sum {A^n_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A^n_\mathrm{P}}} -\frac{1}{A^n_\mathrm{P}}\frac{1}{V_\rP} \sum p_\rP^n\bfS_f +\frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_\rP^n\bfS_f\end{split}\\\begin{split}\\\\
-\frac{1}{\left(A^n_\mathrm{P}-\sum A^n_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_\rP^n\bfS_f-\frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_f'\bfS_f\end{split}\\\begin{split}\\\\
=-\frac{\sum {A^n_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A^n_\mathrm{P}}} -\left(\frac{1}{A^n_\mathrm{P}}- \frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\right)\frac{1}{V_\rP}\sum p_\rP^n\bfS_f\end{split}\\\begin{split}\\\\
-\frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_f^*\bfS_f
\end{split}
\end{split}\end{aligned}\end{align} \]
令\(\nabla\cdot\bfU^{**}=0\)有:
(37)\[
\nabla\cdot\left(\frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\nabla p^{*}\right)
=
\nabla\cdot\mathbf{HbyA}_{simplec}^*
\]
在SIMPLEC算法中,
(38)\[
\mathbf{HbyA}_{simplec}^*=\mathbf{HbyA}^* -\left(\frac{1}{A^n_\mathrm{P}}- \frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\right)\nabla p ^n
\]
求解后有\(p^*\)。随后\(\bfU^{**}\)可以下式获得:
(39)\[
\bfU^{**}_\rP
=
\mathbf{HbyA}_{simplec}^* -\frac{1}{\left(A^n_\mathrm{P}+\sum A^n_\mathrm{N}\right)}\nabla p^{*}
\]
同样的,\(\bfU^{**}\)与\(p^{*}\)近似满足动量方程,理论上满足连续性方程。但在这里可以强调SIMPLEC与SIMPLE差别:SIMPLE略去临点的假定特别激进,导致压力修正量过高。SIMPLEC相对温和,因此压力修正量并不会太大,因此SIMPLEC不需要对压力进行场松弛。同样的,每个迭代步求解的\(\bfU^{**}\)与\(p^{*}\)更加满足动量方程。