各位老师好，我最近研究了一下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());
}
我个人推导的公式为：
78b21ff7-7858-4542-a777-3f2e993f3957-image.png
代码里面实现的公式为：
efa31825-beb7-4ee1-8483-980d46400369-image.png
表面张力项多了一个alpha，请大佬看看，不知道是我推导错了还是代码有bug
@东岳
reactingtwophaseEulerFoam中的式（20)和式（21）有个地方符号有问题：
d07eb439-d134-41ab-862d-944dc3a3efd9-image.png
最后压力项前应该为$-\alpha_{d}$
