return 0;
wmake

buoyantPimpleFoam解析

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

buoyantPimpleFoam是OpenFOAM的传热求解器之一,其用于求解瞬态的、由于温度变化导致的密度变化、浮力驱动流动。相对于pimpleFoam,buoyantPimpleFoam的求解变量增加了能量$e$(或$h$)和密度$\rho$。求解思想和pimpleFoam大体相同。另外一个容易混淆的求解器为boussinesqBuoyantPimpleFoam,二者的区别在于后者采用Boussinesq进行假定,且仅仅求解温度方程。在实际应用中,boussinesqBuoyantPimpleFoam使用的范围存在限制。有关Boussinesq假定可以参考CFD中的能量方程

OpenFOAM中并没有附加重力的恒密度单相流求解器。例如在恒密度的icoFoam和simpleFoam中,是压力差引起的流体流动。但在某些情况下,用户想要在不可压缩流体内添加重力场。考虑一个管道内的静止的液体,在进出口没有压差的情况下,附加重力或不附加重力预测的速度是相同的:液体均不会产生流动。但附加重力会引起压力自上而下的一个均匀增强的渐进变化。因此在附加重力的不可压缩求解器中,可以预测这种压力自顶而下越来越大的情况。buoyantPimpleFoam可以通过将beta(热膨胀系数)设置为常量,达到类似植入重力的pimpleFoam求解器。在这个情况下,重力引起均一的压力梯度变化,且无对流速度。

方程求解和离散

首先有连续性方程:

\begin{equation}\label{cont} \frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \mathbf{U})=0 \end{equation}

动量方程:

\begin{equation}\label{mom} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U})-\nabla \cdot (\mu \nabla \mathbf{U})=-\nabla p + \rho \mathbf{g} \end{equation}

其中$\rho$为流体的密度,$\mathbf{U}$为速度矢量,$t$表示时间,$p$表示压力,$\mathbf{g}$表示重力加速度。参考icoFoam解析,对方程\eqref{mom}半离散化之后有:

\begin{equation}\label{mom2} A_{\mathrm{P}}\mathbf{U}_{\mathrm{P}}+\sum A_{\mathrm{N}}\mathbf{U}_{\mathrm{N}}-E_{\mathrm{P}}=-\nabla p + \rho \mathbf{g} \end{equation}
\begin{equation}\label{mom3} \mathbf{U}_{\mathrm{P}} =\frac{1}{A_{\mathrm{P}}}\left( -\sum A_{\mathrm{N}}\mathbf{U}_{\mathrm{N}}+E_{\mathrm{P}} \right)- \frac{1}{A_{\mathrm{P}}} \left( \nabla p - \rho \mathbf{g} \right) \end{equation}

其中$A_{\mathrm{P}}$,$A_{\mathrm{N}}$为离散后的矩阵系数。令

\begin{equation}\label{HbyA} \mathbf{HbyA}=\frac{1}{A_{\mathrm{P}}}\left( -\sum A_{\mathrm{N}}\mathbf{U}_{\mathrm{N}}+E_{\mathrm{P}} \right) \end{equation}

把方程\eqref{HbyA}中代入到方程\eqref{mom3}中,有

\begin{equation}\label{mom4} \mathbf{U}_{\mathrm{P}} = \mathbf{HbyA}-\frac{1}{A_{\mathrm{P}}} \left( \nabla p - \rho \mathbf{g} \right) \end{equation}

参考interFoam解析中的推导:$\nabla p - \rho \mathbf{g}=\nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho$,其中$\mathbf{h}$为位置矢量,有

\begin{equation}\label{mom5} \mathbf{U}_{\mathrm{P}} = \mathbf{HbyA} - \frac{1}{A_{\mathrm{P}}} \left( \nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho \right) \end{equation}

把方程\eqref{mom5}代入到\eqref{cont}有

\begin{equation}\label{cont2} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \left( \mathbf{HbyA}-\frac{1}{A_{\mathrm{P}}} \left( \nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho \right) \right) \right)=0 \end{equation}
\begin{equation}\label{cont3} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \mathbf{HbyA}_{\mathrm{P}} - \frac{\rho}{A_{\mathrm{P}}} \mathbf{g} \mathbf{h} \nabla \rho \right) -\nabla \cdot \left( \frac{\rho}{A_{\mathrm{P}}} \nabla p_{\mathrm{rgh}} \right)=0 \end{equation}

对于可压缩流体有

\begin{equation}\label{rho} \rho_1 = \rho_0+\mathrm{psi}\left(p_1-p_0 \right) = \rho_0+\mathrm{psi}\left(p_{\mathrm{rgh1}}-p_{\mathrm{rgh0}} \right) \end{equation}

其中$\mathrm{psi}$表示可压缩性。其对时间$t$的微分为

\begin{equation}\label{rhoF} \frac{\partial \rho_1}{\partial t} = \frac{\partial \rho_0}{\partial t}+\mathrm{psi}\frac{\partial \left(p_{\mathrm{rgh1}}-p_{\mathrm{rgh0}} \right)}{\partial t} \end{equation}

把方程\eqref{rhoF}代入到\eqref{cont3}有

\begin{equation}\label{p} \frac{\partial \rho_0}{\partial t}+\mathrm{psi}\frac{\partial \left(p_{\mathrm{rgh1}}-p_{\mathrm{rgh0}} \right)}{\partial t}+\nabla \cdot \left(\rho \left( \mathbf{HbyA}_{\mathrm{P}} - \frac{1}{A_{\mathrm{P}}} \mathbf{g} \mathbf{h} \nabla \rho \right) \right) -\nabla \cdot \left( \frac{1}{A_{\mathrm{P}}} \nabla p_{\mathrm{rgh}} \right)=0 \end{equation}

下面分析能量方程,主要的推导过程请参考CFD中的能量方程

\begin{equation} \frac{\partial \rho h}{\partial t}+\nabla \cdot (\rho \mathbf{U} h) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K) - \frac{\partial p}{\partial t}= \rho r - \nabla \cdot q + \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{hEqn} \end{equation}

在buoyantPimpleFoam中,忽略了$\nabla \cdot(\tau \cdot \mathbf{U})$和$\rho r$。并且认为热通量$q$可以指定为$\alpha_\mathrm{eff} h$。因此,能量方程变为

\begin{equation}\label{E} \frac{\partial \rho h}{\partial t}+\nabla \cdot (\rho \mathbf{U} h) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K) - \frac{\partial p}{\partial t}= - \nabla \cdot \left(\alpha_\mathrm{eff} \nabla h \right) + \rho \mathbf{g} \cdot \mathbf{U} \end{equation}

方程\eqref{E}即最终的能量方程。

代码分析

下面分析主要的代码。首先进入createFields.H

// Force p_rgh to be consistent with p
p_rgh = p - rho*gh; // 声明p_rgh场

然后求解密度方程rhoEqn.H,其对应方程\eqref{cont}:

    fvScalarMatrix rhoEqn
    (
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
        fvOptions(rho)
    );

然后然后求解UEqn.H:

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        fvOptions(rho, U)//方程(7)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()//方程(7)右端的两项
            )
        );

        fvOptions.correct(U);
        K = 0.5*magSqr(U);//K的定义参考"CFD中的能量方程"
    }

能量方程EEqn.H

{
    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)//方程(14)第1,2项
      + fvc::ddt(rho, K) + fvc::div(phi, K)//方程(14)第3,4项
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )//当求解内能的时候,参考"CFD中的能量方程"中的方程14中的第8项
          : -dpdt//当求解内能的时候,方程(14)第5项
        )
      - fvm::laplacian(turbulence->alphaEff(), he)//方程(14)第6项
     ==
        rho*(U&g)//方程(14)第7项
      + radiation->Sh(thermo)//其他辐射源项
      + fvOptions(rho, he)//其他源项
    ); 
}

接下来求解压力方程。压力方程的组建思想和其他求解器相同。即通过连续性方程组建压力泊松方程,在此不再赘述。仅摘取压力方程如下:

    fvScalarMatrix p_rghDDtEqn
    (
        fvc::ddt(rho)  //方程(12)第1项
      + psi*correction(fvm::ddt(p_rgh))//方程(12)第2项
      + fvc::div(phiHbyA)//方程(12)第3项
      ==
        fvOptions(psi, p_rgh, rho.name())
    );

...

        fvScalarMatrix p_rghEqn
        (
            p_rghDDtEqn
          - fvm::laplacian(rhorAUf, p_rgh)//方程(12)第4项
        );

求解后需要更新压力:$p=p_\mathrm{rgh}+\rho \mathbf{g}\cdot \mathbf{h}$。然后再次进去密度方程rhoEqn.H。 另外在代码中,psi*correction(fvm::ddt(p_rgh))主要对应方程\eqref{p}的第2项(离散的对角阵)。correction()实际返回的为A - (A & A.psi()),其中A为我们离散后的矩阵系数。A.psi()为$p_{\mathrm{rgh0}}$。首先,通过fvm::ddt(p_rgh)离散获取p_rgh的离散矩阵系数,然后通过correction函数,将已知的A & A.psi()减去(需要注意的是这一步重组后进入方程组的右方)。另外需要提及的是,correction函数在最终收敛的时候,失去作用。

2018.04.13重修页面排版

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