twoPhaseEulerFoam中的dilatation rates



  • 最近在推导体积分数alpha的方程。理论上来说compressibleInterFoam和twoPhaseEulerFoam中的体积分数方程应该是相同的(不考虑相间的质量传递时),即“东岳流体”中的公式:
    0_1509951529273_Capture.PNG
    在compressibleInterFoam中,dilatation项定义为:

    dgdt =
                (
                    pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
                  - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
                );
    

    这个和东岳流体中的定义是一致的。
    而在twoPhaseEulerFoam中,dilatation项的定义为:

                fluid.dgdt() =
                (
                    alpha1*(pEqnComp2 & p_rgh)
                  - alpha2*(pEqnComp1 & p_rgh)
                );
    

    这个显然和公式推导出的不一致,但是在后续的体积分数方程计算时,我并没有看到程序对这项有什么特殊处理,这是否是错误的?
    另外还想问问各位大神,pEqnComp2 & p_rgh是代表什么意思?pEqnComp2 是fvScalarMatrix类型,而p_rgh是volScalarField类型,中间的操作符&表示什么,是OpenFoam进行算符重构了吗?



  • dgdt要和MULES求解的相方程匹配。单看twoPhaseEulerFoam,在OpenFOAM-2.2.x中定义为:

    dgdt =
                (
                    pos(alpha2)*(pEqnComp2 & p)/rho2
                  - pos(alpha1)*(pEqnComp1 & p)/rho1
                );
    

    在OpenFOAM-2.3.x中定义为:

    fluid.dgdt() =
                (
                    alpha1*(pEqnComp2 & p_rgh)
                  - alpha2*(pEqnComp1 & p_rgh)
                );
    

    我并没有看到程序对这项有什么特殊处理

    你应该看的是OpenFOAM-2.3.x以后的版本。你可以看一下alphaEqn的处理。

    &是进行重构了。矩阵稀疏×场,在OpenFOAM中很常见,类似一种反离散 ,还原回去。



  • @李东岳 多谢东岳。还有几处不太明白。
    (1)OpenFoam对 &的重构是在哪定义的呢?具体到当前的这个&,类似于反离散?具体返回什么值呢?
    (2)我看的是OpenFOAM-4.0.x。无论dgdt的定义如何,最后应该都能回归到体积分数的求解公式,也就是:0_1509957228618_Capture.PNG
    OpenFOAM-2.3中twoPhaseEulerFoam里dgdt的定义,与OpenFOAM4.0.x里compressibleInterFoam的定义是相同的,最后在显式和隐式源项相加后还能恢复到公式。
    但是OpenFOAM4.0.x里twoPhaseEulerFoam的定义中,显式和隐式源项相加后并不能得到体积分数公式。在twoPhaseSystem中,其源项定义为:

     forAll(dgdt_, celli)
            {
                if (dgdt_[celli] > 0.0)
                {
                    Sp[celli] -= dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
                    Su[celli] += dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
                }
                else if (dgdt_[celli] < 0.0)
                {
                    Sp[celli] += dgdt_[celli]/max(alpha1[celli], 1e-4);
                }
            }
    

    也就是得到,alpha*Sp+Su=dgdt + Su这个结果,这个公式显然与体积分数公式不同。体积分数公式中右端项为0_1509957915146_1.PNG
    而现在得到的结果是:0_1509957939111_2.PNG



  • 是的,你说的对。我又看了下,在23x以后,dgdtpEqnComp的定义为:

    pEqnComp1 =
                (
                    fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
                  - (fvOptions(alpha1, rho1)&rho1)
                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
                )/rho1
              + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
    
            pEqnComp2 =
                (
                    fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
                  - (fvOptions(alpha2, rho2)&rho2)
                  - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
                )/rho2
              + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
    

    在这俩个定义中,和分别包含了alpha,因此再乘以时候,就是alpha1*alpha2



  • @李东岳 多谢多谢,是我自己看程序不仔细,pEqnComp1确实已经包含了alpha1。


登录后回复
 

与 CFD中文网 的连接断开,我们正在尝试重连,请耐心等待