CFD中文网

最近在做关于mixutre模型的二次开发。遇到了一个疑惑。是关于连续性方程
d(rho)/dt+div (rhoU)=0 (http://dyfluid.com/twoLiquidMixingFoam.html中的方程6)
的求解代码。
例如，在interFoam中，不需要求解此方程。而是求解alpha方程，和不可压缩形式的连续性方程
div(U)=0。事实上，我觉得这个方程也不是严格成立的。因为他对纯水和纯气的区域当然成立。但是在液面的两相混合处，我觉得理论上来说，描述方程并不可如此简化。但是毕竟液面范围很小，所以这种近似可以接受。
但是对于两相可以混合的twoLiquidMixingFoam和driftFluxFoam而言，我觉得上述方程d(rho)/dt+div (rhoU)=0则必须求解。但是我对比了下这三个求解器的代码，感觉其中的UEqn和pEqn并无太大区别。
所以我很好奇：

1. 以上求解器是否求解了d(rho)/dt+div (rhoU)=0？
2. 如果求解了，其代码段在何处？
3. 如果有Foamer愿意指点下关于不可压两相流的连续性方程的一些基础问题，我也很感兴趣。

希望大家一起讨论，希望大神指导赐教

@dzsmoglai 在 tPEF湍流模型 中说：

oPhaseEulerFoam求解器的时候，总是得

感谢@东岳 老师的答复。目前此问题已经解决。原因也很简单。因为CGKE是基于kEpsilon的。而在PhaseCompressible湍流模型类中，在编译CGKE之前，事先通过PhaseCompressible这一类中的宏makRASModel(kEpsilon)编译了PhaseCompressible版本的kEpsilon模型。因此CGKE中任何函数在调用kEpsilon时就可以调用的到。本问题的解决办法就是同样的在makePhaseCompressibles.C中make一下RNGkE就行了，即makRASModel(RNGkEpsilon)。
如果有其它Foamer对上述解释不理解，我就再补充两句。kEpsilon, RNGkEpsilon这类湍流模型实际上是通用的模型，在incompressible，compressible, phaseCompressible等不同类的湍流模型类中有不同的体现。当然，这个涉及到OpenFoam的湍流模型代码结构，集成关系最多有8层左右，尚不了解的童鞋可以去参考此博客https://marinecfd.xyz/post/openfoam-templated-turbulence-modeling-framework-part-1/。链接是第一部分，一共三个部分。

各位Foamer大家好。
最近我在修改openFoam中twoPhaseEulerFOAM中的湍流模型。主要是因为目前里边的continousGasKEpsilon（CGKE）湍流模型是基于kEpsilon的。而我想创建一个基于RNGkEpsilon的版本，叫做continousGasRNGkEpsilon(CGRNGKE)。在CGKE基础上修改代码后，在phaseCompressible湍流分类中也已经通过编译，形成了新的库文件phaseCompressibleTurbulenceModels.so。以上均没有问题。但是在重新编译twoPhaseEulerFoam求解器的时候，总是得到错误如下：
undefined reference to `Foam::RASModels::RNGkEpsilon<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::PhaseCompressibleTurbulenceModelFoam::phaseModel > > >::typeName'。
不知道哪位大神知道原因，敬请指教。谢谢

@东岳 感谢东岳老师的解答。我现在就是想修改interFoam使用的湍流模型。目前基于kEpsilon自己修改了一个，主要修改的就是把默认为0的Epsilon源项，修改为一个公式，这个公式涉及到alpha的一些运算，但是因为alpha时geometricOneField所以无法进行加减和求梯度。所以我就很困惑：
（1）CMULES算法求解的alpha肯定是个volScalarField啊（从时间文件夹中存储的alpha值就可以知道）
（2）从interFoam的Make/options里边也可以发现，它用的就是incompresibleTurbulenceModel这个库
（3）在interFoam调用turbulence->correct命令时，湍流模型的具体实现过程中，从kEpsilon.C源代码中也可以发现它调用了alpha(),那这个alpha到底时volScalarField还是geometricOneField呢？
再就是，我最想解答的
(4)如果想在epsilon source term中加入基于volScalarField类型的alpha的一些计算，该如何实现？

谢谢！

大家好！这是我在CFD中文网的第一篇帖子，希望能得到大家的帮助。
在OpenFOAM里，对于不可压流动，相分数Alpha的类型为geometricOneField。而对于可压缩流动alpha则为volScalarField.但是很奇怪的是，interFoam求解器使用了GeometricOneField这个类型来表示alpha。不难理解，interFOAM是一个采用VOF方法追踪液面的mixture模型，其alpha值在界面处肯定会取到0-1之间的值。这就带来了我的第一个问题：
（1）interFoam为什么选择geometricOneField为alpha的数据类型？
其次，我想对kEpsilon湍流模型进行些修改，然后在尝试进行一些基于alpha的加减乘除运算时，编译器就会报错。因为alpha的类型GeometricOneField类的定义里，似乎只定义了geometricOneField和自身的乘除法运算符，而没有它和常数scalar或者volScalarField类型数据的计算方法。所以根本无法通过编译。但是在构建epsilon和k方程的时候。又可以见到如下代码

   // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    volScalarField& nut = this->nut_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();

    volScalarField::Internal divU
    (
        fvc::div(fvc::absolute(this->phi(), U))().v()
    );

    tmp<volTensorField> tgradU = fvc::grad(U);
    volScalarField::Internal G
    (
        this->GName(),
        nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
    );
    tgradU.clear();

    // Update epsilon and G at the wall
    epsilon_.boundaryFieldRef().updateCoeffs();

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn
    (
        fvm::ddt(alpha, rho, epsilon_)
      + fvm::div(alphaRhoPhi, epsilon_)
      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
     ==
        ***C1_*alpha()*rho()*G*epsilon_()/k_()***
      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
      - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
      + epsilonSource()
      + fvOptions(alpha, rho, epsilon_)
    );

    epsEqn.ref().relax();
    fvOptions.constrain(epsEqn.ref());
    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
    solve(epsEqn);
    fvOptions.correct(epsilon_);
    bound(epsilon_, this->epsilonMin_);

但是我自己写的alpha()*scalar(0.5)就无法通过编译。

这就是我的第二个问题：
（2） 如何对geometricOneField类型的变量进行和常数及volScalarField进行计算？

以上两个问题，希望大神进行下解答！