• 返回主页
  • buoyantPimpleFoam解析

    2016.11.23:统一HbyA符号;
    2016.09.04:有关correction()函数的描述;
    求解器类型:可压缩浮力驱动流求解器。

    浮力(浮升力)是流体在重力作用下由于流体中各部分密度的不均匀而引起的一种体积力。可见,浮力的产生起源于密度的不均一。buoyantPimpleFoam是OpenFOAM传热求解器之一。主要用于考虑密度差引起的浮力不可忽略的工况下的计算。相比较于pimpleFoambuoyantPimpleFoam的求解变量增加了能量$e$(或$h$)和密度$\rho$。求解框架和pimpleFoam大体相同。另外类似的求解器为boussinesqBuoyantPimpleFoam,其和buoyantPimpleFoam的区别主要是前者采用Boussinesq进行假定,并且求解的是温度方程。在实际应用中,boussinesqBuoyantPimpleFoam使用的范围存在限制。有关Boussinesq假定可以参考CFD中的能量方程

    首先我们有连续性方程: \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解析,我们对方程($\ref{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}_{\mathrm{P}}=\frac{1}{A_{\mathrm{P}}}\left( -\sum A_{\mathrm{N}}\mathbf{U}_{\mathrm{N}}+E_{\mathrm{P}} \right) \end{equation} 把方程($\ref{HbyA}$)中代入到方程($\ref{mom3}$)中,有: \begin{equation}\label{mom4} \mathbf{U}_{\mathrm{P}} = \mathbf{HbyA}_{\mathrm{P}} - \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}_{\mathrm{P}} - \frac{1}{A_{\mathrm{P}}} \left( \nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho \right) \end{equation} 把方程($\ref{mom5}$)代入到($\ref{cont}$)有: \begin{equation}\label{cont2} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \left( \mathbf{HbyA}_{\mathrm{P}} - \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 \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} 对于可压缩流体,我们有: \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} 对其去时间$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} 把方程($\ref{rhoF}$)代入到($\ref{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}
    下面我们进入代码分析。首先进入createFields.H(省略部分已经解析的内容):

    #include "readGravitationalAcceleration.H"
    #include "readhRef.H"
    #include "gh.H"//定义重力加速度矢量和位置矢量
    
    
    Info<< "Reading field p_rgh\n" << endl;
    volScalarField p_rgh
    (
        IOobject
        (
            "p_rgh",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    
    // Force p_rgh to be consistent with p
    p_rgh = p - rho*gh;//定义$p_{\mathrm{rgh}}$
    
    Info<< "Creating field kinetic energy K\n" << endl;
    volScalarField K("K", 0.5*magSqr(U));
    

    其中的volScalarField K的定义请参考:CFD中的能量方程。然后,在buoyantPimpleFoam.C中,首先求解rhoEqn.H

    {
        fvScalarMatrix rhoEqn
        (
            fvm::ddt(rho)
          + fvc::div(phi)
          ==
            fvOptions(rho)
        );
    
        fvOptions.constrain(rhoEqn);
    
        rhoEqn.solve();
    
        fvOptions.correct(rho);
    }
    

    其对应方程($\ref{cont}$)。然后求解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中的能量方程"
        }
    

    参考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} 方程($\ref{E}$)即对应EEqn.H

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

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

        fvScalarMatrix p_rghDDtEqn
        (
            fvc::ddt(rho) + //方程($\ref{p}$)第1项
          + psi*correction(fvm::ddt(p_rgh))//方程($\ref{p}$)第2项
          + fvc::div(phiHbyA)//方程($\ref{p}$)第3项
          ==
            fvOptions(psi, p_rgh, rho.name())
        );
    
    ...
    
            fvScalarMatrix p_rghEqn
            (
                p_rghDDtEqn
              - fvm::laplacian(rhorAUf, p_rgh)//方程($\ref{p}$)第4项
            );
            p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
    

    求解后,我们需要更新压力:$p=p_\mathrm{rgh}+\rho \mathbf{g}\cdot \mathbf{h}$。然后再次进去密度方程rhoEqn.H。在大部分OpenFOAM求解器中,密度方程求解了2次(原因不详)。经过自带的算例测试,求解一次密度方程和两次密度方程的结果并没有显著差别。

    另外在代码中,存在代码psi*correction(fvm::ddt(p_rgh)),其主要对应方程($\ref{p}$)第2项(离散的对角阵)。correction()实际返回的为A - (A & A.psi()),其中A为我们离散后的矩阵系数。A.psi()为$p_{\mathrm{rgh0}}$。首先,我们通过fvm::ddt(p_rgh)离散获取p_rgh的离散矩阵系数,然后通过correction函数,将已知的A & A.psi()减去(需要注意的是这一步重组后进入方程组的右方)。另外需要提及的是,correction函数在最终收敛的时候,失去作用。
    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国