OpenFOAM中标准k-e湍流模型的一点疑问

[史浩]

对于standard k-$ϵ$ 模型，其k与$ϵ$的控制方程如下：
$$\begin{array}{}\text{(1)}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho k\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho u_{i}k\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+G_k-\rho ϵ\end{array}$$
$$\begin{array}{}\text{(2)}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho ϵ\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho u_iϵ\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[(\mu +\frac{{\mu }_t}{{\sigma }_{ϵ}})\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }{x}_j}\right]+C_{1ϵ}\frac{ϵ}{k}{G}_k-C_{2ϵ}\rho \frac{{ϵ}^2}{k}\end{array}$$
在OpenFOAM中，其实现代码在$(FOAM_SRC)/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C文件中，对应的代码如下：

     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_()      //对应epsilon方程右侧第二项
       - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)  //当U不满足连续方程时的修正项？
       - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)  //我的问题：和epsilon方程对应的右侧第三项不符
       + epsilonSource()                  //自定义源项？
       + fvOptions(alpha, rho, epsilon_)  //fvOption添加的源项
     );
 
     epsEqn.ref().relax();
     fvOptions.constrain(epsEqn.ref());
     epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
     solve(epsEqn);
     fvOptions.correct(epsilon_);
     bound(epsilon_, this->epsilonMin_);
 
     // Turbulent kinetic energy equation
     tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(alpha, rho, k_)
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
         alpha()*rho()*G  //对应k方程右侧第二项
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)  //连续性误差修正项
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)   //我的问题：与k方程右侧第三项不对应
       + kSource()                  //自定义k源项
       + fvOptions(alpha, rho, k_)  //fvOption对应的源项
     );
 
     kEqn.ref().relax();
     fvOptions.constrain(kEqn.ref());
     solve(kEqn);
     fvOptions.correct(k_);
     bound(k_, this->kMin_);

方程中，$G_k$为湍流生成项，一些书籍中，其定义如下，为什么代码中硬生生多出来一个$\rho$？
$$\begin{array}{}\text{(3)}& G_k=-\rho \overline{u_i^{{}^{\prime }}{u}_j^{{}^{\prime }}}\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}\simeq 2{\mu }_t{s}_{ij}{s}_{ij}\end{array}$$
我的问题写在了注释上了，是我的湍流方程写错了，还是OpenFOAM中对湍流方程有着特殊的处理？请各位老师指点

[yfclark]

不可压缩湍流模型：

    volScalarField G(GName(), nut_*(tgradU() && dev(twoSymm(tgradU()))));

可压缩湍流模型：

 volScalarField G(GName(), mut_*(tgradU() && dev(twoSymm(tgradU()))));

你看看是不是这样的

[史浩]

@yfclark 多谢！这样我理解了G前面的rho的含义了。麻烦老兄看一下我代码注释中的两个问题。公式补充如下：
k与e的控制方程如下：

上面的两个公式右侧第三项与OpenFOAM中代码中写的不一样，不知道是我公式的问题，还是OpenFOAM中模型的特殊处理，或者两种表达是等价的？

[yfclark]

是对应的，fvm::Sp是源项处理，fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_) 可以理解为
C2_alpha()*rho()*epsilon_()/k_()*epsilon_,只不过把C2_*alpha()*rho()*epsilon_()/k_()做为显式的系数，epsilon_是求解变量做了隐式处理，论坛里有几个帖子你搜索一下

[史浩]

@yfclark Sp是隐式源项添加方法，将源项加到对角线的元素上。
epsilon方程中的源项离散可以理解，不过我觉得写成下面的形式收敛性会不会更好一点（红色部分为隐式离散的系数）？

在k方程中，右侧第三项是-rho* epsilon，这一项对应的变量不是k，怎么隐式离散？为什么OpenFOAM中其隐式离散的系数为rho* epsilon/k？

[yfclark]

典型的教科书离散源项方法，你的离散完全正确，k里面应该时强行离散的，目的只是为了增加求解稳定性吧

[史浩]

@yfclark 多谢老师！解答了我很多疑问。
之所以对k方程右侧第三项的隐式离散比较感兴趣，是因为他是对另外一个变量的隐式离散，具体处理应该是这样的：

通过上述的方式在k方程中对epsilon进行了隐式的离散，增加求解的稳定性
如果给定任意一个传输方程组，形式如下：

两个方程的右侧第一个源项能很好的隐式离散，但是右侧第二个源项和另外一个方程的变量有关，这个怎么隐式离散？是完全凭经验，还是有一些比较好的方法进行处理？

[大喵]

大佬，我也遇到和你一样的疑问，为啥G前面多了个rho，不吝赐教万分感谢！！！

[李东岳]

@史浩 对于你${\varphi }_1$和${\varphi }_2$这种方程，已经属于耦合系统。只能显性离散。比如${\varphi }_1$的方程${\varphi }_2$显性离散

@大喵 二楼已经回答了啊

[史浩]

@东岳 多年的僵尸帖被大佬翻牌子，吓了一跳。当时求解一个耦合方程，自己想着最理想的方式就是全耦合，两个变量组装成一个矩阵，整体求解，然后自己就掉进牛角尖里出不了。后来想k-e模型中两个变量是耦合在一起的，就学习了一下OpenFOAM中湍流方程的处理方法。现在用显示耦合方法计算，然后多迭代一两次，精度并不差，就放弃了最开始整体求解的想法。
@大喵 这个对方程和代码再好好对应一下，OpenFOAM里面的和经典湍流模型方程是一致的

[大喵]

谢谢两位大佬@史浩 @东岳