多相流求解器中的曳力项实现

田畔的风

最近研究了一下OpenFOAM-v2206的twoPhaseEulerFoam求解器，注意到在文件pU/Ueqns.H中，动量方程的代码为

U1Eqn =
        (
            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
          - fvm::Sp(contErr1, U1)
          + MRF.DDt(alpha1*rho1 + Vm, U1)
          + phase1.turbulence().divDevRhoReff(U1)
         ==
          - Vm
           *(
                fvm::ddt(U1)
              + fvm::div(phi1, U1)
              - fvm::Sp(fvc::div(phi1), U1)
              - DDtU2
            )
          + fvOptions(alpha1, rho1, U1)
        );
        U1Eqn.relax();
        U1Eqn += fvm::Sp(Kd, U1);
        fvOptions.constrain(U1Eqn);
        U1.correctBoundaryConditions();
        fvOptions.correct(U1);

我注意到曳力项被写作U1Eqn += fvm::Sp(Kd, U1);

然而按照我的经验，曳力似乎应该被写作
$$
F_D = K_d (\mathbf{U}_1 - \mathbf{U}_2)
$$
也就是后面的速度应该是两相速度的差值。但是好像OpenFOAM所有的多相流求解器都是按照开头代码那样实现的，非常困惑...恳请各位指点一下...

李东岳

$$
F_D = -K_d (\mathbf{U}_1 - \mathbf{U}_2)
$$
是这个吧？

田畔的风

@李东岳李老师你好，这边主要不是纠结符号正负，主要是疑惑为什么求解器中的曳力实现是$K_d \mathbf{U}_1$而不是$K_d (\mathbf{U}_1-\mathbf{U}_2)$？

李东岳

@田畔的风另外那一项在压力方程里面，你看一下压力方程

田畔的风

@李东岳感谢李老师！现在对这个求解器理解了一些，对于$\mathbf{U}_1$的控制方程, $K_d \mathbf{U}_2$一项作为显式项被直接纳入pEqn.H的界面通量的修正。

CFD中文网

多相流求解器中的曳力项实现