CFD 多相VOF算法勘误

学流体的小明

最近在写论文，因为用的是interFoam求解器，所以在写控制方程的时候一致参考李东岳老师的这个笔记 http://dyfluid.com/interFoam.html 。最近仔细看了下密度、运动粘度、动力粘度之间的关系，发现一直以来都没注意到的错误：
错误在于公式24

$$
\nabla\cdot\boldsymbol{\tau}=\nabla\cdot\left(\nu\left(\nabla\bfU+\nabla\bfU^{\mathbf{T}}\right)\right)
$$

在多相流的控制方程里面，密度$\rho$是和速度放一起了，所以上面公式当中的 $\nu$ 应该为 $\mu$ 。

interFOAM类的求解器用的是不可压缩不相容的混合物模型，在代码中是这样定义的

immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);

它要求输入每一个相的密度$\rho$和运动粘度$\nu$，而$\nu$是在$\mu$的基础上计算的。
在immiscibleIncompressibleTwoPhaseMixture 这个类的父类incompressibleTwoPhaseMixture中有计算类成员nu_的函数：

void Foam::incompressibleTwoPhaseMixture::calcNu()
{
    nuModel1_->correct();
    nuModel2_->correct();

    const volScalarField limitedAlpha1
    (
        "limitedAlpha1",
        min(max(alpha1_, scalar(0)), scalar(1))
    );

    // Average kinematic viscosity calculated from dynamic viscosity
    nu_ = mu()/(limitedAlpha1*rho1_ + (scalar(1) - limitedAlpha1)*rho2_);
}

而$\mu$的计算公式为 $\mu = \alpha_1 \rho_1 \nu_1 + (1-\alpha_1) \rho_2 \nu_2$：

Foam::tmp<Foam::volScalarField>
Foam::incompressibleTwoPhaseMixture::mu() const
{
    const volScalarField limitedAlpha1
    (
        min(max(alpha1_, scalar(0)), scalar(1))
    );

    return tmp<volScalarField>::New
    (
        "mu",
        limitedAlpha1*rho1_*nuModel1_->nu()
      + (scalar(1) - limitedAlpha1)*rho2_*nuModel2_->nu()
    );
}

也就是说$\nu$公式和东岳老师笔记 http://dyfluid.com/interFoam.html 的公式25不一致，正确的公式是
$$\nu = \frac{{{\alpha _1}{\rho _1}{\nu _1} + \left( {1 - {\alpha _1}} \right){\rho _2}{\nu _2}}}{{{\alpha _1}{\rho _1} + \left( {1 - {\alpha _1}} \right){\rho _2}}}$$

那么粘性应力在程序中是怎样计算的？
interFOAM的动量方程在UEqn.H当中，是

fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
    );

可以看到粘性应力项和湍流粘度项混在一起，直接由湍流模型的对象turbulence给出。但是只有当给出具体湍流模型的时候，才能知道究竟是怎样的表达式。
我再往下去就有点晕了，不知道要找哪个源代码文件当中的具体函数了，求大佬指点一下。就是没找到

不去探求源代码的话，我的理解是，湍流模型给出的都是nut，在各个算例文件当中都给的是动力湍流粘度。
来自湍流的动力粘度为$\rho \nu_t$是要用当地的密度乘以湍流模型给出的运动粘度。
也就是说：这里的turbulence->divDevRhoReff(rho, U)在WALE模型下，或者其它的一些模型下，计算粘性应力 + 湍流应力在方程中的那一项，对应的表达式为
$$\nabla\cdot\boldsymbol{\tau}=\nabla\cdot\left( \mu_{eff} \left(\nabla\bfU+\nabla\bfU^{\mathbf{T}}\right)\right)$$
其中
$${\mu _{eff}} = {\mu} + {\mu _t} = {\mu} + \rho {\nu _t} = {\alpha _1}{\rho _1}{\nu _1} + \left( {1 - {\alpha _1}} \right){\rho _2}{\nu _2} + \left[ {{\alpha _1}{\rho _1} + \left( {1 - {\alpha _1}} \right){\rho _2}} \right]{\nu _t}$$

李东岳

感谢勘误！我回法国更新一下。谢谢

Amadeus

太有帮助了! 感谢

CFD中文网

CFD 多相VOF算法 勘误

CFD 多相VOF算法勘误