simpleFoam解析


李东岳


1. 引言

simpleFoam为OpenFOAM中一个基于有限体积法的,用于求解稳态不可压缩牛顿流体N-S方程的求解器。在simpleFoam中,压力和速度的耦合通过SIMPLE/SIMPLEC算法进行计算,且没有考虑其他体积力(如重力等)。类似于icoFoam,simpleFoam可用于理解 OpenFOAM中压力速度耦合的框架策略。对于流场内温度变化较大,但可忽略浮力以及重力的可压缩稳态求解器,可选用rhoSimpleFoam。附加重力的可压缩浮力驱动流稳态求解器可选用buoyantSimpleFoam。本文从最基本的N-S方程在笛卡尔网格下的有限体积离散开始推导,介绍simpleFoam求解器中植入的SIMPLE/SIMPLEC算法。需要注意的是,由于为不可压缩流体求解器,下文中$p$表示压力除以密度。

2. 控制方程与算法
2.1 有限体积法离散
首先有稳态的N-S方程: \begin{equation} \nabla\cdot\bfU=0, \label{C} \end{equation} \begin{equation} \nabla \cdot (\mathbf{U}\mathbf{U})-\nabla \cdot(\nu \nabla \mathbf{U})=-\nabla p, \label{mom} \end{equation} 对方程\eqref{mom}进行对速度$\bfU$隐性离散有:
\begin{equation} \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)} , \label{mom2} \end{equation}
\begin{equation} {\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 } , \label{mom3} \end{equation}
\begin{equation} \int \nabla p \mathrm{d}V=\int p \mathrm{d}\bfS=\sum \left(p_f^n\bfS_f\right), \label{gradp} \end{equation} 其中上标$^n$表示为当前迭代步(已知),上标$^*$表示预测步(待求),下标$_f$表示网格单元面上的值,$V_\rP$表示网格单元体积,$\bfS$表示网格单元的各个面的面矢量,$F_f$为通量,$\nu$为动力粘度。方程\eqref{mom3}中$(\nabla\mathbf{U}^*)_f\cdot\bfS_f $需要进一步处理,其中$\nabla\mathbf{U}^*$为定义在体心的量,即网格体心的速度梯度。$(\nabla\mathbf{U}^*)_f$为定义在面心的量,即网格面心的速度梯度。为使用紧致基架点防止数值震荡,$(\nabla\mathbf{U}^*)_f\cdot\bfS_f$通常表示为面法向梯度与面矢量的模的乘积: \begin{equation} (\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| \label{sngrad} \end{equation} 其中$(\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} $表示面法向梯度,有 \begin{equation} (\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|}}, \label{sngrad2} \end{equation} 其中$\mathbf{d}$表示网格单元体心之间矢量差(距离矢量),下标$_\mathrm{N}$表示相邻网格单元的速度,下标$_\mathrm{P}$表示当前网格单元的速度。在这里需要注意的是,方程\eqref{sngrad2}只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性。将方程\eqref{mom2}-\eqref{sngrad}代入到方程\eqref{mom}有: \begin{equation} \sum F_f^n \mathbf{U}_f^* = \sum { \nu (\nabla\mathbf{U}^*)_f\cdot\bfS_f } -\sum p_f^n\bfS_f. \label{momF} \end{equation}

其中的通量$F_f^n$采用当前已知时间步的速度$\mathbf{U}^n$来计算,同时保留预测速度$\mathbf{U}^*$为未知量。

2.2. 动量预测
方程\eqref{momF}中$\mathbf{U}_f^*$被定义为面上的预测速度,然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度,在此步可以引入各种插值格式。假设使用中心线性格式:
\begin{equation} \mathbf{U}_f^* = \frac{{\mathbf{U}_\mathrm{P}^* + \mathbf{U}_\mathrm{N}^*}}{2}. \label{linear} \end{equation}
\begin{equation} p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}. \label{linear2} \end{equation}
将方程\eqref{sngrad},\eqref{linear}代入到方程\eqref{momF}有
\begin{equation} \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, \end{equation}
\begin{equation} \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{equation}
将上式简写为 \begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^*{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^*} = - \frac{1}{V_\rP}\sum p_f^n\bfS_f, \label{apanmom} \end{equation} 其中 \begin{equation} A_\mathrm{P}= \frac{1}{V_\rP}\left(\sum {\frac{{F_f^n}}{2}} + \sum {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}}\right) , \end{equation} \begin{equation} A_\mathrm{N}= \frac{1}{V_\rP}\left({\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} \right), \end{equation}

求解方程\eqref{apanmom}即可获得预测速度$\mathbf{U}^*_\mathrm{P}$。方程\eqref{apanmom}即为OpenFOAM中的动量预测方程。需要注意的是,动量预测步骤并不是必须的,这主要和速度压力耦合求解策略有关。如果不调用动量预测步骤,则$\mathbf{U}^*_\mathrm{P}=\mathbf{U}^n_\mathrm{P}$。

2.3. 压力泊松方程
N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度,但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来,其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下,对连续性方程$\nabla\cdot\bfU=0$进行离散后的形式为 \begin{equation} \sum \mathbf{U}_{\mathrm{P},f}^{n+1} \cdot \bfS_f=0. \label{cont} \end{equation} 如果能用压力表示方程\eqref{cont}中的$\mathbf{U}_{\mathrm{P},f}^{n+1}$,是不是就可以得出压力泊松方程了呢?答案是肯定的。首先,方程\eqref{apanmom}中的压力项不建议进行离散,这是因为方程右侧的梯度项都有可能引起数值的震荡(例如双流体模型中的湍流分散力也需要特殊处理)。同时,参考方程\eqref{apanmom}, 在收敛情况下(时间步$t=n+1$)为: \begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{n+1}} = -\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f, \label{apanmomNew} \end{equation} 若将方程\eqref{apanmomNew}减去\eqref{apanmom}有:
\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}'+\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}'} = -\frac{1}{V_\rP}\sum p_f '\bfS_f, \label{simple} \end{equation}
其中$\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U}_\mathrm{P}^{*}$(其他同理)。在这里未调用任何略去临点的假定,方程\eqref{apanmomNew}是严格成立的(精准的)。既压力的修正量和当前网格点的速度修正量以及相邻网格点的速度修正量有关。对方程\eqref{apanmomNew}移项有 \begin{equation} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{A_\mathrm{P}}\frac{1}{V_\rP}\sum p_f^{n+1}\bfS_f. \end{equation} \begin{equation} \mathbf{HbyA}^{n+1}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{n+1}} } \right). \label{hbya} \end{equation} 面上的速度可以通过下面的公式获得: \begin{equation} \mathbf{U}_{\mathrm{P},f}^{n+1} = \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP}\sum p_f^*\bfS_f\right)_f. \label{uf} \end{equation} 将方程\eqref{uf}代入到方程\eqref{cont}有: \begin{equation}\label{eg} \sum \left( \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP}\sum p_f^*\bfS_f\right)_f \right) \cdot \bfS_f =0, \end{equation} \begin{equation} \sum {\mathbf{HbyA}_f^{n+1}} \cdot \bfS_f = \sum {\frac{1}{{A_{\mathrm{P},f}}}\left(\frac{1}{V_\rP}\sum p_f^*\bfS_f\right)_f} \cdot \bfS_f, \label{jim} \end{equation} 方程\eqref{jim}即下述方程的离散形式: \begin{equation} \nabla \cdot \mathbf{HbyA}^{n+1} = \nabla \cdot\left(\frac{1}{A_{\mathrm{P}}} \nabla p^{n+1}\right). \label{poss} \end{equation}

方程\eqref{poss}即为压力泊松方程。若有$\mathbf{HbyA}^{n+1}$,则可求得收敛的压力。

2.4. SIMPLE算法
可以看出,最后求解的速度$\mathbf{U}_\mathrm{P}^{n+1}$和压力$p^{n+1}$应该符合方程\eqref{apanmom}和\eqref{poss},分别对应动量方程和连续性方程。然而目前,我们只有通过2.1节动量预测步骤求出来的$\mathbf{U}_\mathrm{P}^*$。$\mathbf{HbyA}_\mathrm{P}^{n+1}$和$p^{n+1}$均为未知的。方程\eqref{poss}不能求解。SIMPLE迭代算法则可以解决这个问题。依据SIMPLE算法,首先引入略去临点影响的假定,方程\eqref{simple}可以写为: \begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}' = -\frac{1}{V_\rP}\sum p_f '\bfS_f, \label{uprime} \end{equation} 将方程\eqref{uprime}和\eqref{apanmom}相加有 \begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{**}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^*} = -\frac{1}{V_\rP} \sum p_f^*\bfS_f, \label{pisoapan} \end{equation} 其中$\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^*+\mathbf{U}_\mathrm{P}'$,$p_f^{*}=p_f^n+p_f'$。相比方程\eqref{apanmomNew},方程\eqref{pisoapan}由于忽略了临点的影响并不是精准的。对方程\eqref{pisoapan}进行移项有 \begin{equation} \mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\frac{1}{V_\rP} \sum p_f^*\bfS_f. \end{equation} \begin{equation} \mathbf{HbyA}^{*}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{*}} } \right). \label{hbya2} \end{equation} \begin{equation} \mathbf{U}_{\mathrm{P},f}^{**} = \mathbf{HbyA}_f^{*} - \frac{1}{{{A_{\mathrm{P},f}}}}\left(\frac{1}{V_\rP} \sum p_f^*\bfS_f\right)_f. \label{pisouf} \end{equation} 方程\eqref{pisouf}中可参考方程\eqref{eg}-\eqref{poss}的步骤组建压力泊松方程,即 \begin{equation} \nabla \cdot \mathbf{HbyA}^{*} = \nabla \cdot\left(\frac{1}{A} \nabla p^{*}\right). \label{pisoposs} \end{equation}

对方程\eqref{pisoposs}求解后有$p^{*}$。将$p^*$回代到方程\eqref{pisouf}的时候,有$\mathbf{U}_{\mathrm{P}}^{**}$,继续求解有$p^{**}$以及$\mathbf{U}_{\mathrm{P}}^{***}$。总而言之,SIMPLE算法中的的迭代过程可以表示为下面几个步骤:

  1. 依据当前时间步的条件求解速度$\mathbf{U}^*$,组建$\mathbf{HbyA}^*$;
  2. 求解方程\eqref{pisoposs}求解压力$p^*$;
  3. 通过方程\eqref{pisouf}依据压力$p^*$以及$\mathbf{HbyA}^*$更新速度有$\mathbf{U}^{**}$;
  4. 回到第一步进行下一次迭代;
2.5. SIMPLEC算法

SIMPLEC算法认为SIMPLE的算法太过激进。这种略去邻点的行为会使得压力修正值$p'$过大,在某些情况下会引起发散。于是在SIMPLE算法中,需要调用松弛技术来减小每一次的压力修正。然而松弛技术却减慢了收敛速率。SIMPLEC算法认为任意网格点的速度修正可以近似认为是相邻网格点速度修正的平均: \begin{equation} \bfU_\rP'\approx\frac{\sum A_\rN\bfU_\rN'}{\sum A_\rN} \:\rightarrow \: \bfU_\rP'\sum A_\rN\approx\sum A_\rN\bfU_\rN' \label{simplecnew} \end{equation} 将方程\eqref{simplecnew}代入到\eqref{simple}有: \begin{equation} A_\rP\bfU_\rP'=-\bfU_\rP'\sum A_\rN-\frac{1}{V_\rP} \sum p_f'\bfS_f \label{simplecnew2} \end{equation} 即 \begin{equation} \bfU_\rP'=-\frac{1}{A_\rP+\sum A_\rN}\frac{1}{V_\rP} \sum p_f'\bfS_f \label{simplec2} \end{equation} 相对比于SIMPLE算法中的方程\eqref{uprime},对于固定的$\bfU_\rP',方程$\eqref{simplec2}右侧压力修正的乘积项$\frac{1}{A_\rP+\sum A_\rN}$必然大于$\frac{1}{A_\rP}$(有限体积法离散的$A_\rN<0$),进而$p_f'$的值不会特别大,因此SIMPLEC算法并不需要对压力进行激进的低松弛修正。将方程\eqref{simplec2}和方程\eqref{apanmom}相加有:

\begin{equation} {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, \label{simplec3} \end{equation}
\begin{equation} \begin{split} \mathbf{U}_\mathrm{P}^{**} = -\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A_\mathrm{P}}} -\frac{1}{A_\mathrm{P}}\frac{1}{V_\rP} \sum p_\rP^n\bfS_f- \frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\frac{1}{V_\rP} \sum p_f'\bfS_f \\\\ =-\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A_\mathrm{P}}} -\frac{1}{A_\mathrm{P}}\frac{1}{V_\rP} \sum p_\rP^n\bfS_f +\frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_\rP^n\bfS_f \\\\ -\frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_\rP^n\bfS_f-\frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_f'\bfS_f \\\\ =-\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{*}}}{{A_\mathrm{P}}} -\left(\frac{1}{A_\mathrm{P}}- \frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\right)\frac{1}{V_\rP}\sum p_\rP^n\bfS_f \\\\ -\frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\frac{1}{V_\rP}\sum p_f^*\bfS_f \end{split} \label{simplec4} \end{equation}
令$\nabla\cdot\bfU^{**}=0$有:
\begin{equation} \nabla\cdot\left(\frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\nabla p^{*}\right) = \nabla\cdot\left(\bfHbyA^* -\left(\frac{1}{A_\mathrm{P}}- \frac{1}{\left(A_\mathrm{P}+\sum A_\mathrm{N}\right)}\right)\nabla p ^n \right) \label{simplec6} \end{equation}
求解后有$p^*$。SIMPLEC算法中的的迭代过程总结如下:

  1. 依据当前时间步求解速度$\mathbf{U}^*$;
  2. 求解方程\eqref{simplec6}求解压力$p^*$;
  3. 通过方程\eqref{simplec4}更新速度有$\mathbf{U}^{**}$;
  4. 回到第一步进行下一次迭代;

3. 验证算例
3.1. 后向台阶流

simpleFoam自带若干算例。在OpenFOAM-7中,自带的算例主要为以下几个:

上述算例在此不做过多介绍。本文提供的算例对标相关文献的研究工作 (Armaly et al., 1983)


算例网格为采用blockMesh生成的2D纯六面体结构网格。左侧给定非均一的速度进口,右侧为速度出口,上下为用于附加剪切的壁面。算例的雷诺数为$\mathrm{Re}_S\approx 389$。

constant文件夹的物性主要设置如下:

相关算例运行首先需要执行网格转换命令:

blockMesh

网格生成完毕,后直接运行直接模拟求解器运行即可:

simpleFoam

大约迭代1175次收敛。收敛后,可运行下行命令提取相关数据:

postProcess -func sampleDict -latestTime

相关数据后,可运行下行命令对数据进行处理(实验结果在算例文件中已经包含)

./data/makeGraph

数据处理采用的是gnuplot ,结果如下所示:


点击下载算例

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网