大于治水
首发于大于治水
interPhaseChangeFoam求解器解析

interPhaseChangeFoam求解器解析

学了of半年,为了研究空化模拟问题,死磕interPhaseChangeFoam求解器代码,写了一些个人见解,与同行们分享,欢迎交流~

另,本文内容交阅于东岳老师,前辈会给出更加细致与独到的讲解,另外版本请参阅东岳流体。(已发布)

------------------分割-------------------

正文开始:~

interPhaseChangeFoam是用来求解两种有质量传递的不可压流体多相流求解器,界面处理方式是利用VOF方法进行捕捉。其控制方程主要由四部分组成,分别为:连续方程,动量方程,相方程,质量传递模型组成,其中质量传递模型取Schnerr_Sauer空化模型进行分析,具体形式如下:

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) = \frac { \dot { m } } { \rho _ { l } } \\ \nabla \cdot U = \dot { m } \left[ \frac { 1 } { \rho _ { l } } - \frac { 1 } { \rho _ { v } } \right] \\ \frac { \partial ( \rho U ) } { \partial t } + \nabla \cdot ( \rho U U ) = - \nabla p _ { r g h } + \nabla \cdot \left[ \mu _ { m } \left( \nabla U + \nabla U ^ { T } \right) \right] - g \cdot h \nabla \rho + \sigma \kappa \nabla \alpha _ { l } \\ \dot { m } = \left\{ \begin{array} { l } { \frac { 3 \rho _ { l } \rho _ { v } } { \rho } \alpha ( 1 - \alpha ) \frac { 1 } { R _ { b } } \sqrt { \frac { 2 \left( p _ { v } - p \right) } { 3 \rho _ { l } } } \left( p _ { v } > p \right) } \\ { \frac { 3 \rho _ { l } \rho _ { v } } { \rho } \alpha ( 1 - \alpha ) \frac { 1 } { R _ { b } } \sqrt { \frac { 2 \left( p - p _ { v } \right) } { 3 \rho _ { l } } } \left( p _ { v } < p \right) } \end{array} \right.

其中,R_b=\sqrt[\frac{1}{3}]{\frac{1-\alpha_l}{\alpha_l} \frac{3}{4\pi} \frac{1}{n}} 为平均泡半径,n为液相中的空泡数目, p_v为对应温度的饱和蒸汽压力 。对于其他诸如表面张力\sigma \kappa \nabla\alpha_lp_{rgh} 等请参考 interFoam求解器

对数值求解过程有个整体的认识,介绍流程如下:首先,对相方程进行求解,分为显式和隐式,更新\alpha_l及其通量和混合密度\rho;其次,由动量方程给出速度预测值;最后,求解泊松方程进行压力值更新与速度修正,并更新通量。较之interFoam求解器而言,该求解器最大的不同在于空化模型的加入和处理,所以本文主要解析如下几点:

1. 空化模型;(SchnerrSauer.C/H, phaseChangeTwoPhaseMixture.C/H)

2. 相方程的求解;(alphaControls.H, alphaEqn.H, alphaEqnSubCycle.H)

3. 泊松方程的求解。(pEqn.H)

Section 1: 空化模型

给出在SchnerrSauer.H头文件中的基本参数,主要有四个,在计算设置时,需要用户在transportProperties内指定。统计信息如表1。

在此基础上,延伸出重要的二级中间参数有三个,它们定义在SchnerrSauer.C的成员函数里,说明如表2。

最后给出质量传递源项,该项分别在求解相方程和泊松方程中加入。具体分别是:

SchnerrSauer.C中的mDotAlphal()和mDotP()的质量传递源项;

phaseChangeTwoPhaseMixture.C中的vDotAlphal()和vDotP()体积传递源项。

这四项返回值均为Pair<tmp<volScalarField>>,即为一个数组,代表了净空化传质 率和净冷凝传质率,它们具体的形式在相方程和泊松方程的求解中会继续分析,下面给出具体形式:

这里要注意,在相方程和泊松方程的求解中,对传质源项采用的都是隐式处理。

Section 2: 相方程

对相方程的求解需对原方程进行变换。在Weller所提出的MULES方法中,为了使interface面更佳清晰(sharp),处理时加入了人工压缩项\nabla \cdot (U_c \alpha_l \alpha_v),其中U_c保证界面压缩,\nabla保证守恒性, \alpha_l \alpha_v保证有界性。具体形式如Formulation 1:

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } \mathrm { U } \right) + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) - \nabla \cdot \left( \alpha _ { l } U \right) = \frac { \dot { m } } { \rho _ { l } }

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } \mathrm { U } \right) + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) - \nabla \cdot \left( \alpha _ { l } \left( \alpha _ { v } U _ { v } + \alpha _ { l } U _ { l } \right) \right) = \frac { \dot { m } } { \rho _ { l } }

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } U \right) + \nabla \cdot \left( \alpha _ { l } \alpha _ { v } U _ { c } \right) = \frac { \dot { m } } { \rho _ { l } }

这里强调一点,为什么说\nabla \cdot (U_c \alpha_l \alpha_v)是人工压缩项呢?因为在VOF假设中,各相是共享速度和压力的。也就是说,基于这个假设,U_c=0。所以U_c是人为设定的,其形式为:U_c=c_\alpha \left| U \right| \frac{\nabla \alpha}{ \left| \nabla \alpha \right|}

from alphaEqnSubCycle.H,phic代表的是c_\alpha \left| U \right|

surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));

from alphaEqn.H,phir代表的是U_c,interface.nHatf()即为\frac{\nabla \alpha}{ \left| \nabla \alpha \right|}

surfaceScalarField phir("phir", phic*interface.nHatf());

在相方程的隐式求解中,加入了速度散度项的影响,其目的和在动量方程中考虑速度散度是一样的,都是为了强调连续方程的作用。不过不同的是,由于质量传递的作用,该模型的速度散度不为0。具体演化过程如Formulation 2:

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) = \alpha _ { l } ( \nabla \cdot U ) - \alpha _ { l } ( \nabla \cdot U ) + \frac { \dot { m } } { \rho _ { l } }

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) = \alpha _ { l } ( \nabla \cdot U ) - \alpha _ { l } \left( \dot { m } \left[ \frac { 1 } { \rho _ { l } } - \frac { 1 } { \rho _ { v } } \right] \right) + \frac { \dot { m } } { \rho _ { l } }

\frac { \partial \alpha _ { l } } { \partial t } + \nabla \cdot \left( \alpha _ { l } U _ { l } \right) = \alpha _ { l } ( \nabla \cdot U ) + \dot { m } \left[ \frac { 1 } { \rho _ { l } } - \alpha _ { i } \left( \frac { 1 } { \rho _ { l } } - \frac { 1 } { \rho _ { v } } \right) \right]

在alphaControls.H中,有以下控制计算的参数:

  • nAlphaCorr,代表的是\alpha_l通量修正的次数;
  • nAlphaSubCycles,代表的是将每个时间步长均分成对应个数的子步长,循环后加和
  • MULESCorr,即Weller所提出的MULES修正,一般选为yes;
  • icAlpha,各向同性压缩系数,类似于一个松弛因子的作用,一般为0。

以上参数均在fvSolution中的alpha子字典中进行定义。

以下代码演示的是在alphaEqnSubCycle.H中,控制相方程求解的整个大循环,此举的目的在于可使用大的时间步长。

可以看出,求解相方程的次数为nAlphaSubCycles,当该参数大于1时,需要将每一个子时间步的结果进行加和。

   volScalarField divU(fvc::div(phi));

    if (nAlphaSubCycles > 1)
    {
        dimensionedScalar totalDeltaT = runTime.deltaT();
        surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);

        for
        (
            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
            !(++alphaSubCycle).end();
        )
        {
            #include "alphaEqn.H"
            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
        }

        rhoPhi = rhoPhiSum;
    }
    else
    {
        #include "alphaEqn.H"
    }

    rho == alpha1*rho1 + alpha2*rho2;

下面进入相方程的具体求解步骤。在alphaEqn.H中,对Formulation2进行离散,可以看出,在MULESCorr为yes时,第一步是预测 \alpha _l ;相方程时间项和对流项的离散方法已定,也就是说在fvScheme中不可选;另外,在源项的处理上,对\dot{m}中的\alpha_l进行隐式处理。对应的代码为:

   if (MULESCorr)
    {
        fvScalarMatrix alpha1Eqn
        (
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)
          - fvm::Sp(divU, alpha1)
         ==
            fvm::Sp(vDotvmcAlphal, alpha1)
          + vDotcAlphal
        );

        alpha1Eqn.solve();

//省略了write()函数

        talphaPhi = alpha1Eqn.flux();
    }

对于MULES通量修正方法,它是基于FCT(Flux-Corrected Transport, 文章链接) 发展而来。定义的talphaPhiCorr为相方程的Formulation1左侧的后两项;MULESCorr为yes时,通过Formulation1与2相减修正相通量(MULES算法分析在此)。如果MULESCorr为false则相方程的求解采用显式方法。

Section 3: 泊松方程

将动量方程离散成以下形式:

A _ { P } U _ { P } = \sum A _ { N } U _ { N } + S _ { P } \left( - \left( \nabla p _ { r g h } \right) - g \cdot h \nabla \rho + f _ { \sigma , P } \right) V

U _ { P } = \frac { 1 } { A _ { P } } \sum A _ { N } U _ { N } + \frac { 1 } { A _ { P } } \left( - \left( \nabla p _ { r g h } \right) - g \cdot h \nabla \rho + f _ { \sigma , P } \right)

H_P^A=\frac{1}{A_P^U}\sum{A_N U_N} ,向面上插值有:

U _ { f } = H _ { f } ^ { A } + \frac { 1 } { A _ { f } } \left( - \left( \nabla p _ { r g h } \right) _ { f } - ( g \cdot h ) _ { f } ( \nabla \rho ) _ { f } + f _ { \sigma , f } \right)

乘以面积S_f,通量为:

\phi _ { f } = U _ { f } \cdot S _ { f } = H _ { f } ^ { A } \cdot S _ { f } + \frac { 1 } { A _ { f } } \left( - ( g \cdot h ) _ { f } ( \nabla \rho ) _ { f } + f _ { \sigma , f } \right) \cdot S _ { f } - \frac { 1 } { A _ { f } } \left( \nabla p _ { r g h } \right) _ { f } \cdot S _ { f }

其中,PhiHbyA为H_f^A \cdot S_f,Phig为\frac{1}{A_f^U} (-(g \cdot h)_f (\nabla \rho)_f + f_{\sigma,f} ) \cdot S_f,代码如下:

   volScalarField rAU("rAU", 1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);

    surfaceScalarField phig
    (
        (
            interface.surfaceTensionForce()
          - ghf*fvc::snGrad(rho)
        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;//注意此处,PhiHbyA与Phig二者相加

由连续方程有:

\sum U _ { f } \cdot S _ { f } = \left [ \text { vDotP() }  [ 0 ] - \text { vDotP() } [ 1 ] \right ] \left( p - p _ { v } \right)

\sum P h i H b y A + \sum P h i g - \sum \frac { 1 } { A _ { f } } \left( \nabla p _ { r g h } \right) _ { f } \cdot S _ { f } = \left [ \text { vDotP() }  [ 0 ] - \text { vDotP() } [ 1 ] \right ] \left( p_{rgh} - \left( p _ { v } - \rho g h \right) \right)

对应代码如下:

   while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
          - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
        );

        p_rghEqn.setReference(pRefCell, pRefValue);

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + p_rghEqn.flux();

            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

可以看出在p_rghEqn方程中加入了质量传递项,与interFoam稍有不同,而其他部分基本一致。

最后,由p_{rgh}求解p代码如下,如果进出口边界都为自然边界,则需要在fvSolution中设定参考压力值。

   p == p_rgh + rho*gh;

    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rho*gh;
    }

编辑于 2018-11-01

文章被以下专栏收录