CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    两类pEqn.H求解可压缩流动计算时稳定性问题

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

      有对pimple算法求解pEqn.H有研究的老哥吗,关于这个算法有一些困惑,(1)https://github.com/OpenFOAM/OpenFOAM-2.3.x/blob/master/applications/solvers/combustion/reactingFoam/reactingFoam.C
      中密度方程为什么要求解两次,一次出现在runtime++后面,一次包含在pEqn.H里面。
      (2)reactingFoam和rhoreactingFoam中的pEqn.H求解的压力方程并不相同,而且reactingFoam的pEqn.H和rhoPimpleFoam里面的pEqn.H基本上是一致的。有使用过rhoreactingFoam的大佬吗?rhoreactingFoam在求解可压缩流动是是否比reactingFoam有优势?我现在使用的是reactingFoam这一类的pEqn.H,计算过程中一个时间步内压力变化很大,导致温度出现了震荡现象,导致物性也发生了震荡,无法收敛。
      下面贴出这两中pEqn.H
      rhoPimpleFoam/reactingFoam

      if (pimple.transonic())
      {
          surfaceScalarField phid
          (
              "phid",
              fvc::interpolate(psi)
             *(
                  (fvc::interpolate(HbyA) & mesh.Sf())
                + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
              )
          );
      
          fvOptions.makeRelative(fvc::interpolate(psi), phid);
      
          while (pimple.correctNonOrthogonal())
          {
              fvScalarMatrix pEqn
              (
                  fvm::ddt(psi, p)
                + fvm::div(phid, p)
                - fvm::laplacian(rhorAUf, p)
                ==
                  fvOptions(psi, p, rho.name())
              );
      
              fvOptions.constrain(pEqn);
      
              pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
      
              if (pimple.finalNonOrthogonalIter())
              {
                  phi == pEqn.flux();
              }
          }
      }
      

      rhoreactingFoam

       thermo.rho() -= psi*p;
      
          volScalarField rAU(1.0/UEqn.A());
          surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
      
          volVectorField HbyA("HbyA", U);
          HbyA = rAU*UEqn.H();
      
          if (pimple.transonic())
          {
              surfaceScalarField phiHbyA
              (
                  "phiHbyA",
                  (
                      (fvc::interpolate(rho*HbyA) & mesh.Sf())
                    + rhorAUf*fvc::ddtCorr(rho, U, phi)
                  )/fvc::interpolate(rho)
              );
      
              fvOptions.makeRelative(phiHbyA);
      
              surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
      
              phiHbyA *= fvc::interpolate(rho);
      
              fvScalarMatrix pDDtEqn
              (
                  fvc::ddt(rho) + fvc::div(phiHbyA)
                + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
              );
      
              while (pimple.correctNonOrthogonal())
              {
                  fvScalarMatrix pEqn
                  (
                      pDDtEqn
                    - fvm::laplacian(rhorAUf, p)
                   ==
                      fvOptions(psi, p, rho.name())
                  );
      
                  fvOptions.constrain(pEqn);
      
                  pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
      
                  if (pimple.finalNonOrthogonalIter())
                  {
                      phi = phiHbyA + pEqn.flux();
                  }
              }
          }
          thermo.rho() += psi*p;
      

      只贴了transonic部分的

      1 条回复 最后回复 回复 引用
      • C
        cccrrryyy 教授 最后由 编辑

        我最近也在思考这个问题,不知道你这边有没有什么进展。我看了OpenFOAM-4和OpenFOAM-6还有一些其他版本的,分享一下我目前的一些看法。

        就你发的这个2.3.x而言,主程序里面runtime++后面的#include "rhoEqn.H"是一个根据速度分布求解密度的方程,应该只是为了有一个初始的密度场,不然后续的计算无法进行。在真正开始pimple循环后,EEqn.H末尾有一个thermo.correct()。这一步根据计算后的能量场,更新了温度、粘度、可压缩比(psi)等等物性,所以在pEqn.H的开头,出现了一步rho = thermo.rho(),实际上进行了更新密度的操作(rho = psi * p,注意这个时候psi已经变了),之后再进行压力的求解。注意,在pEqn.H的79行和80行,再次调用"rhoEqn.H"(计算根据速度场应该得到的密度场,或者说是没有更新之前的密度场),以及 "compressibleContinuityErrs.H"(计算前面所说的密度场和由热物性所得到的密度场之间的区别,即密度场更新前和更新后的差异)。不知道我以上的理解对不对。

        我还没弄明白的问题和你一样,reactingFoam和rhoReactingFoam的压力方程为什么代码上有差异。看看有没有别人回答吧。

        I don't want to survive, I want to thrive.

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

          乍一看,三年过去了

          http://dyfluid.com/buoyantSimpleFoam.html

          在真正开始pimple循环后,EEqn.H末尾有一个thermo.correct()。这一步根据计算后的能量场,更新了温度、粘度、可压缩比(psi)等等物性,所以在pEqn.H的开头,出现了一步rho = thermo.rho(),实际上进行了更新密度的操作(rho = psi * p,注意这个时候psi已经变了),之后再进行压力的求解。

          最近更新了这个,尤其是关于两个密度更新的问题的讨论。不可压缩流速度修正与压力变量直接相关。但是附加浮力之后,速度修正与压力变量以及密度变量(浮力)直接相关。我赞同上面的理解。同时,能量方程之后的密度更新可有可无,只不过是收敛性的问题。

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

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

          1 条回复 最后回复 回复 引用
          • Moved from Algorithm by  李东岳 李东岳 
          • First post
            Last post