reactingmultiphaseEulerFoam中的表面张力实现

yfclark

各位老师好，我最近研究了一下reactingmultiphaseEulerFoam，对其中的表面张力的计算有一些困惑：
在pEqn.H中定义了：

        phigFs.set
        (
            phasei,
            (
                alpharAUfs[phasei]
               *(
                   ghSnGradRho
                 - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
                 - fluid.surfaceTension(phase)*mesh.magSf()
                )
            ).ptr()
        );

最后通过该式更新各项的速度：


            mSfGradp = pEqnIncomp.flux()/rAUf;

            forAll(phases, phasei)
            {
                phaseModel& phase = phases[phasei];

                phase.U() =
                    HbyAs[phasei]
                  + fvc::reconstruct
                    (
                        alpharAUfs[phasei]*mSfGradp
                      - phigFs[phasei]
                    );
                phase.U().correctBoundaryConditions();
                fvOptions.correct(phase.U());
            }

我个人推导的公式为：

代码里面实现的公式为：

表面张力项多了一个alpha，请大佬看看，不知道是我推导错了还是代码有bug
@东岳
reactingtwophaseEulerFoam中的式（20)和式（21）有个地方符号有问题：

最后压力项前应该为$-\alpha_{d}$

CFD中文网

reactingmultiphaseEulerFoam中的表面张力实现