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} \otimes\mathbf{U})-\nabla \cdot(\nu \nabla \mathbf{U})=-\nabla p, \label{mom} \end{equation}

对方程\eqref{mom}进行对速度$\bfU$隐性离散有:

\begin{equation} \int \int {\nabla \cdot \left( {\mathbf{U}\otimes\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\mathbf{U}\otimes\mathbf{U}\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\mathbf{U}^n}\otimes{\mathbf{U}^r}} \right)}_f}} \bfS_f\Delta t = \sum {\left( {F_f^n \mathbf{U}_f^r} \right)} \Delta t, \label{mom2} \end{equation}
\begin{equation} {\int \int \nabla \cdot \left( {\nu \nabla \mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\nu \nabla \mathbf{U}\mathrm{d}\bfS\mathrm{d}t = \sum {\left( {\nu \nabla {\mathbf{U}^r}} \right)} } _f}\bfS_f\Delta t = \sum {\left( {\nu \left( {\nabla _f^ \bot {\mathbf{U}^r}} \right)} \right)} \Delta t, \label{mom3} \end{equation}
\begin{equation} \int \int \nabla p \mathrm{d}V\mathrm{d}t=\int\int p \mathrm{d}\bfS\mathrm{d}t=\sum \left(p_f^n\bfS_f\right)\Delta t, \label{gradp} \end{equation}

其中上标$^n$表示为当前的时间步(已知),上标$^r$表示预测时间步(待求),下标$_f$表示网格单元面上的值,$V_\rP$表示网格单元体积,$\bfS$表示网格单元的各个面的面矢量,$\Delta t$表示时间步长,$F_f$为通量,$\nu$为动力粘度。方程\eqref{mom3}中$\nabla _f^ \bot {\mathbf{U}^r}$(面法向梯度)在正交网格中进一步可以表示为:

\begin{equation} \nabla _f^ \bot {\mathbf{U}^r}{\rm{ = }}{\left( {\nabla {\mathbf{U}^r}} \right)_f}\bfS_f = \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^r - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}, \label{sngrad} \end{equation}

其中$d$表示网格单元体心之间的距离,下标$_\mathrm{N}$表示相邻网格单元的速度,下标$_\mathrm{P}$表示当前网格单元的速度。$\nabla _f^ \bot {\mathbf{U}^n}$在三维情况下为一个矢量。在这里需要注意的是,方程\eqref{sngrad}只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性。将方程\eqref{mom2}-\eqref{sngrad}代入到方程\eqref{mom}有:

\begin{equation} \sum {\left( {F_f^n \mathbf{U}_f^r} \right)} = \sum {\left( {\nu \nabla _f^ \bot {\mathbf{U}^r}} \right)} -\sum \left(p_f^n\bfS_f\right). \label{momF} \end{equation}

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

2.2. 动量预测

方程\eqref{momF}中$\mathbf{U}_f^r$被定义为面上的预测速度,然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度,在此步可以引入各种插值格式。假设使用中心线性格式:

\begin{equation} \mathbf{U}_f^r = \frac{{\mathbf{U}_\mathrm{P}^r + \mathbf{U}_\mathrm{N}^r}}{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 {\left( {F_f^n \frac{{\mathbf{U}_\mathrm{P}^r + \mathbf{U}_\mathrm{N}^r}}{2}} \right)} = \sum {\left( {\nu \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^r - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}} \right)}-\sum \left(p_f^n\bfS_f\right), \end{equation}
\begin{equation} \left( { \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right)\mathbf{U}_\mathrm{P}^r = - \sum {\left( {\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)\mathbf{U}_\mathrm{N}^r} \right)} -\sum \left(p_f^n\bfS_f\right). \end{equation}

将上式简写为

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^r{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^r} = -\sum \left(p_f^n\bfS_f\right), \label{apanmom} \end{equation}

其中

\begin{equation} A_\mathrm{P}=\left( { \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right), \end{equation}
\begin{equation} A_\mathrm{N}=\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right), \end{equation}

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

2.3. 压力泊松方程

N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度,但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来,其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下,对连续性方程$\nabla\cdot\bfU=0$进行离散后的形式为

\begin{equation} \sum \left(\mathbf{U}_{\mathrm{P},f}^{n+1} \cdot \bfS_f\right)=0. \label{cont} \end{equation}

如果能用压力表示方程\eqref{cont}中的$\mathbf{U}_{\mathrm{P},f}^{n+1}$,是不是就可以得出压力泊松方程了呢?答案是肯定的。首先,方程\eqref{apanmom}中的压力项不建议进行离散,这是因为方程右侧的梯度项都有可能引起数值的震荡(例如双流体模型中的湍流分散力也需要特殊处理)。这样方程\eqref{apanmom}可以写为:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{r}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}} = -\nabla p_\mathrm{P}^{n}, \label{apanmomr} \end{equation}

在收敛情况下(时间步$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}} = -\nabla p_\mathrm{P}^{n+1}, \label{apanmomNew} \end{equation}

若将方程\eqref{apanmomNew}减去\eqref{apanmomr}有:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}'+\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}'} = -\nabla p_\mathrm{P}', \label{simple} \end{equation}

其中$\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U}_\mathrm{P}^{r}$(其他同理)。在这里未调用任何略去临点的假定,方程\eqref{apanmomNew}是严格成立的(精准的)。既压力的修正量和当前网格点的速度修正量以及相邻网格点的速度修正量有关。对方程\eqref{apanmomNew}移项有

\begin{equation} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\nabla {p^{n+1}}. \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}

通过Rhie-Chow插值,面上的速度可以通过下面的公式获得:

\begin{equation} \mathbf{U}_{\mathrm{P},f}^{n+1} = \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}}. \label{uf} \end{equation}

将方程\eqref{uf}代入到方程\eqref{cont}有:

\begin{equation}\label{eg} \sum \left( \left( \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}} \right) \cdot \bfS_f \right)=0, \end{equation}
\begin{equation} \sum {\left( {\mathbf{HbyA}_f^{n+1}} \cdot \bfS_f\right) } = \sum {\left( {\frac{1}{{A_{\mathrm{P},f}}}{\nabla _f}{p^{n+1}}} \cdot \bfS_f\right) }, \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}^r$。$\mathbf{HbyA}_\mathrm{P}^{n+1}$和$p^{n+1}$均为未知的。方程\eqref{poss}不能求解。SIMPLE迭代算法则可以解决这个问题。依据SIMPLE算法,首先引入略去临点影响的假定,方程\eqref{simple}可以写为:

\begin{equation} A_\mathrm{P}\bfU_\mathrm{P}'=- \nabla p_\mathrm{P}', \label{uprime} \end{equation}

将方程\eqref{uprime}和\eqref{apanmomr}相加有

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}} = -\nabla p_\mathrm{P}^{*}, \label{pisoapan} \end{equation}

其中$\mathbf{U}_\mathrm{P}^{*}=\mathbf{U}_\mathrm{P}^{r}+\mathbf{U}_\mathrm{P}'$,$p_\mathrm{P}^{*}=p_\mathrm{P}^{r}+p_\mathrm{P}'$。相比方程\eqref{apanmomNew},方程\eqref{pisoapan}由于忽略了临点的影响并不是精准的。对方程\eqref{pisoapan}进行移项有

\begin{equation} \mathbf{U}_\mathrm{P}^* = \mathbf{HbyA}^{r}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\nabla {p^{*}}. \end{equation}
\begin{equation} \mathbf{HbyA}^{r}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{r}} } \right). \label{hbya2} \end{equation}
\begin{equation} \mathbf{U}_{\mathrm{P},f}^{*} = \mathbf{HbyA}_f^{r} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{*}}. \label{pisouf} \end{equation}

方程\eqref{pisouf}中可参考方程\eqref{eg}-\eqref{poss}的步骤组建压力泊松方程,即

\begin{equation} \nabla \cdot (\mathbf{HbyA}^{r}) = \nabla \cdot(\frac{1}{A_{\mathrm{P},f}} \nabla p^{*}). \label{pisoposs} \end{equation}

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

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

SIMPLEC算法认为SIMPLE的算法太过激进。这种略去邻点的行为会使得压力修正值$p'$过大,在某些情况下会引起发散。于是在SIMPLE算法中,需要调用松弛技术来减小每一次的压力修正。然而松弛技术却减慢了收敛速率。SIMPLEC算法并不想略去邻点的修正$\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}'}$,在方程\eqref{simple}左右两侧减去$\sum {A_\mathrm{N}\mathbf{U}_\mathrm{P}'}$:

\begin{equation} \left(A_\mathrm{P}-\sum A_\mathrm{N}\right)\mathbf{U}_\mathrm{P}'=-\sum A_\mathrm{N}\left(\mathbf{U}_\mathrm{N}'+\bfU_\rP '\right)- \nabla p_\rP', \label{simplec} \end{equation}
略去右侧第一项:
\begin{equation} A_\mathrm{P}\mathbf{U}_\mathrm{P}'=- \frac{A_\mathrm{P}}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP', \label{simplec2} \end{equation}
将方程\eqref{simplec2}和方程\eqref{apanmomr}相加有:
\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}} = -\nabla p_\mathrm{P}^{n}- \frac{A_\mathrm{P}}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP', \label{simplec3} \end{equation}
\begin{equation} \mathbf{U}_\mathrm{P}^{*} =-\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}}}{{A_\mathrm{P}}} -\frac{1}{A_\mathrm{P}}\nabla p_\mathrm{P}^{n}- \frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP' \\\\ =-\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}}}{{A_\mathrm{P}}} -\frac{1}{A_\mathrm{P}}\nabla p_\mathrm{P}^{n} +\frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP ^n \\\\ -\frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP ^n-\frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP' \\\\ =-\frac{\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}}}{{A_\mathrm{P}}} -\left(\frac{1}{A_\mathrm{P}}- \frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\right)\nabla p_\rP ^n \\\\ -\frac{1}{\left(A_\mathrm{P}-\sum A_\mathrm{N}\right)}\nabla p_\rP ^* \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(-\frac{\sum {A_\mathrm{N}\mathbf{U}^{r}}}{{A_\mathrm{P}}} -\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}^r$;
  2. 求解方程\eqref{simplec6}求解压力$p^*$;
  3. 通过方程\eqref{simplec4}更新速度有$\mathbf{U}^{*}$;
  4. 回到第一步进行下一次迭代;

3. 代码分析
...
//方程\eqref{mom}
tmp tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
	
	...
	
	volScalarField rAU(1.0/UEqn.A());
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));//方程\eqref{hbya}
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

    tmp rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
		//方程\eqref{simplec6}的右半部分
		
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }
	
	fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );
	//方程\eqref{simplec6}

东岳流体 版权所有
勘误、讨论、补充内容请前往CFD中文网