各位老师好,我最近研究了一下reactingmultiphaseEulerFoam,对其中的表面张力的计算有一些困惑:
在pEqn.H中定义了:
最后通过该式更新各项的速度:
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}$