reactingmultiphaseEulerFoam中的表面张力实现


  • OpenFOAM讲师

    各位老师好,我最近研究了一下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}$


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2