CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

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

    OpenFOAM
    4
    11
    2696
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 史
      史浩 讲师 最后由 李东岳 编辑

      对于standard k-$\epsilon$ 模型,其k与$\epsilon$的控制方程如下:
      \begin{equation}
      \frac{\partial }{\partial t}\left ( \rho k \right )+\frac{\partial }{\partial x_{i}}\left ( \rho u_{i}k\right )=\frac{\partial }{\partial t}\left [ \left ( \mu +\frac{\mu_t}{\sigma_{k}} \right ) \frac{\partial k}{\partial x_{j}} \right ]+G_{k}-\rho \epsilon
      \end{equation}
      \begin{equation}
      \frac{\partial }{\partial t}\left ( \rho \epsilon \right )+\frac{\partial }{\partial x_{i}}\left ( \rho u_{i} \epsilon \right )=\frac{\partial }{\partial t}\left [ \left ( \mu +\frac{\mu_{ t }}{\sigma_{ \epsilon }} \right ) \frac{\partial \epsilon }{\partial x_{j}} \right ]+C_{1\epsilon}\frac{\epsilon}{k}G_{k}- C_{2\epsilon}\rho \frac{\epsilon^2 }{k}
      \end{equation}
      在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{equation}
      G_{k}=-\rho \overline{u_{i}^{'}u_{j}^{'}}\frac{\partial u_{j}}{\partial x_{i}}\simeq 2\mu_{t}s_{ij}s_{ij}
      \end{equation}
      我的问题写在了注释上了,是我的湍流方程写错了,还是OpenFOAM中对湍流方程有着特殊的处理?请各位老师指点

      让我们随波逐流

      1 条回复 最后回复 回复 引用
      • Y
        yfclark 讲师 最后由 编辑

        不可压缩湍流模型:

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

        可压缩湍流模型:

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

        你看看是不是这样的

        史 1 条回复 最后回复 回复 引用
        • 史
          史浩 讲师 @yfclark 最后由 编辑

          @yfclark 多谢!这样我理解了G前面的rho的含义了。麻烦老兄看一下我代码注释中的两个问题。公式补充如下:
          k与e的控制方程如下:
          943938f6-644c-4fa6-9920-c2e908106852-image.png
          上面的两个公式右侧第三项与OpenFOAM中代码中写的不一样,不知道是我公式的问题,还是OpenFOAM中模型的特殊处理,或者两种表达是等价的?

          让我们随波逐流

          1 条回复 最后回复 回复 引用
          • Y
            yfclark 讲师 最后由 编辑

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

            史 1 条回复 最后回复 回复 引用
            • 史
              史浩 讲师 @yfclark 最后由 李东岳 编辑

              @yfclark Sp是隐式源项添加方法,将源项加到对角线的元素上。
              epsilon方程中的源项离散可以理解,不过我觉得写成下面的形式收敛性会不会更好一点(红色部分为隐式离散的系数)?
              3c093154-7e72-471f-9ee4-38c4c87e942b-image.png
              在k方程中,右侧第三项是-rho* epsilon,这一项对应的变量不是k,怎么隐式离散?为什么OpenFOAM中其隐式离散的系数为rho* epsilon/k?

              让我们随波逐流

              1 条回复 最后回复 回复 引用
              • Y
                yfclark 讲师 最后由 编辑

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

                史 1 条回复 最后回复 回复 引用
                • 史
                  史浩 讲师 @yfclark 最后由 编辑

                  @yfclark 多谢老师!解答了我很多疑问。
                  之所以对k方程右侧第三项的隐式离散比较感兴趣,是因为他是对另外一个变量的隐式离散,具体处理应该是这样的:
                  63bd1e8e-6987-46d9-89a7-1486548a160e-image.png
                  通过上述的方式在k方程中对epsilon进行了隐式的离散,增加求解的稳定性
                  如果给定任意一个传输方程组,形式如下:
                  fe33f579-2424-4d05-8cb3-ad21dfca2477-image.png
                  两个方程的右侧第一个源项能很好的隐式离散,但是右侧第二个源项和另外一个方程的变量有关,这个怎么隐式离散?是完全凭经验,还是有一些比较好的方法进行处理?

                  让我们随波逐流

                  1 条回复 最后回复 回复 引用
                  • 大
                    大喵 最后由 编辑

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

                    1 条回复 最后回复 回复 引用
                    • 李东岳
                      李东岳 管理员 最后由 编辑

                      @史浩 对于你$\phi_1$和$\phi_2$这种方程,已经属于耦合系统。只能显性离散。比如$\phi_1$的方程$\phi_2$显性离散

                      @大喵 二楼已经回答了啊

                      线上CFD课程开始报名:http://www.dyfluid.com/class.html

                      CFD高性能服务器 http://dyfluid.com/servers.html

                      史 1 条回复 最后回复 回复 引用
                      • 史
                        史浩 讲师 @李东岳 最后由 编辑

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

                        让我们随波逐流

                        1 条回复 最后回复 回复 引用
                        • 大
                          大喵 最后由 编辑

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

                          1 条回复 最后回复 回复 引用
                          • First post
                            Last post