OpenFOAM的两个燃烧求解器

cccrrryyy

根据头文件里的描述，reactingFoam是Solver for combustion with chemical reactions，而rhoReactingFoam是Solver for combustion with chemical reactions using density based thermodynamics package。很好奇OpenFOAM为什么要这样设计。两个其实都是要解压力方程的，不像是我们可能认为的一个是压力基求解器而另一个是密度基求解器，那么这两个求解器的区别到底是什么。另外，rhoReactingFoam的压力方程是如何推导出来的呢？

我还注意到一点。如果抛开化学反应和求解组分的质量分数，reactingFoam的动量和压力方程和rhoPimpleFoam完全一样，而rhoPimpleFoam是Transient solver for turbulent flow of compressible fluids for HVAC and similar applications。这里特别强调是暖通之类的应用，是表达什么意思？就是所谓的low-Mach number形式吗？

李东岳

OpenFOAM那面热动力包分为了rhoThermo和psiThermo，所以这俩个求解器应该是一个基于前者一个基于后者，应该收敛性在某些情况下不一样。我没对比过这两个求解器。但是很多算例都可以通用。rhoPimpleFoam是为了区分rhoCentralFoam，这俩算法区别很大，后者算高超更好，前者也能算但激波应该没那么尖锐，但应该跟后者区分开来，buoyantPimpleFoam算HVAC附加了浮力

cccrrryyy

@东岳感谢东岳老师回复！恩，热力学包分为rhoThermo和psiThermo两种，对于在reactingFoam和rhoReactingFoam上的应用来说，我的理解是前者调用psiThermo（和rhoPimpleFoam一样），后者调用rhoThermo。我好奇的是为什么要这样分。比如你提到的rhoPimpleFoam和rhoCentralFoam，前者是压力基求解器，后者是密度基求解器，所以分别调用不同的热力学包，可以理解。而reactingFoam和rhoReactingFoam，两者都是压力基求解器（不知道这样说是不是百分百正确），但是两者的压力方程很不一样。

可能我过于关注细节了，不过用OpenFOAM用的越多就越发现，它的代码设计很考究，任何处理都有其原因，reactingFoam和rhoReactingFoam一定有它们分别存在的理由。而且我感觉关键点就在它们的压力方程上。

李东岳

@cccrrryyy 在 OpenFOAM的两个燃烧求解器中说：

而reactingFoam和rhoReactingFoam，两者都是压力基求解器（不知道这样说是不是百分百正确），但是两者的压力方程很不一样。

方程一样吧？

是的，很多代码是为了更稳健，任何处理都有其原因的

cccrrryyy

OpenFOAM-6里面， reactingFoam的压力方程是

while (pimple.correctNonOrthogonal())
   {
       fvScalarMatrix pEqn
       (
           fvm::ddt(psi, p)
         + fvc::div(phiHbyA)
         - fvm::laplacian(rhorAUf, p)
        ==
           fvOptions(psi, p, rho.name())
       );

       pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

       if (pimple.finalNonOrthogonalIter())
       {
           phi = phiHbyA + pEqn.flux();
       }
   }

rhoReactingFoam的压力方程是直接调用了rhoPimpleFoam的压力方程

fvScalarMatrix pDDtEqn
    (
        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
      + fvc::div(phiHbyA)
     ==
        fvOptions(psi, p, rho.name())
    );

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));

        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }

如果都是常规的动量->能量->压力方程这种求解顺序，很奇怪为什么要弄出不一样的形式来呢。比较起来看，代入rho = psi* p，那么rhoReactingFoam的压力方程多出一项psi*correction(fvm::ddt(p))。这一项好像和rhoReactingFoam获取密度方式有关，在pEqn的前续代码中有体现，但我不太理解。

另外神奇的事是，从OpenFOAM-5.x开始，rhoReactingFoam的压力方程都是直接调用rhoPimpleFoam的压力方程，而reactingFoam用不同的压力方程；但是OpenFOAM-4.x之前，压力方程一样的是rhoPimpleFoam和reactingFoam，然后rhoReactingFoam用不同的压力方程。

李东岳

所以就是说，为什么一个求解器调用fvm::ddt(psi, p)，一个求解器调用fvc::ddt(rho) + psi*correction(fvm::ddt(p))，二者的区别是什么？

李东岳

没有找到任何理论解释这部分内容

如果从代码反推，fvm::ddt(psi, p)和fvc::ddt(rho) + psi*correction(fvm::ddt(p))的区别只不过是前者：$\frac{\p \psi p}{\p t}$，后者$\frac{\p \rho}{\p t}+\psi\frac{\p p}{\p t}$，第一个隐性，第二个显性

另一个我进行的一个解释是，可压缩连续性方程可以展开为
\begin{equation}
\frac{\p p}{\p t}=-c\nabla\cdot(\rho\bfU)
\end{equation}
如果可压缩性较小，$c$会很大导致$\frac{\p p}{\p t}$无解，因此对于可压缩性较小的流动，最好采用fvm::ddt(psi, p)和fvc::ddt(rho) + psi*correction(fvm::ddt(p))这种形式。我还看了一下cavitatingFoam的压力方程代码，OpenFOAM-6中是

fvScalarMatrix pEqn
        (
            fvm::ddt(psi, p)
          - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
          + fvc::div(phi, rho)
          + fvc::div(phiGradp)
          - fvm::laplacian(rhorAUf, p)
        );

OpenFOAM-7中是

fvScalarMatrix pEqn
        (
            fvc::ddt(rho)
          + psi*correction(fvm::ddt(p))
          + fvc::div(phi, rho)
          + fvc::div(phiGradp)
          - fvm::laplacian(rhorAUf, p)
        );

可以看出来这面psi-based和rho-based标准不是很统一。个人感觉应该跟具体的应用有很大关系，某些应用建议psi-based，某些应用建议rho-based。

不确定我说的对，求更多大佬分享经验

cccrrryyy

同感，希望有别的大佬来解释。从这个角度出发，也许可压缩性较小的（所谓的compressible fluid，而不是compressible gas）fvm::ddt(psi, p)可能就足够了？也许OpenFOAM-4之前的版本里面rhoPimpleFoam的描述就是这个意思，因为一般HVAC里面都是compressible fluid吧。但是OepnFOAM之后的版本里面rhoPimpleFoam的压力方程也自己换掉了，和你说的cavitatingFoam类似的情况，有点搞不懂了。
坐等大佬解答。

cccrrryyy

@东岳偶然翻看到之前的讨论，结合东岳老师提到的声速有关的内容和最近的一些新体会，贴上来和大家讨论。

对于可压缩流动的定义好像存在不同的理解，我现在倾向于认为正确的是，是否严格意义的可压缩要看密度是否随压力变化。如果密度是随压力和温度、组分等变化，那么是完全的可压缩。如果密度不随压力变化，而只随温度、组分等变化，也是可压缩（因为密度会变化），但不是完全的可压缩。后者在很多传热、燃烧应用中是广泛存在的（密度受压力影响很小，密度差主要由温度和组分变化导致），在处理的时候大多采用所谓的低马赫数近似（low Mach number approximation）。其本质是假设密度不随压力变化，不一定是字面意思上的马赫数很低才可以使用，更多的应该是指代一套数值方法（这点不是很确定）。FOAM里面有一种状态方程是incompressiblePerfectGas，和普通的perfectGas区别在于，前者的压力用一个参考压力（常数），后者用的则是计算得到的压力（非常数）。

感觉FOAM里面之所以有基于psi的热物性和基于rho的热物性就是为了能够处理不同类型的可压缩流动。传热求解器（solver/heatTransfer）比如buoyantPimpleFoam采用了基于rho的热物性，而燃烧求解器（solver/combustion）比如reactingFoam，fireFoam采用了基于psi的热物性。燃烧求解器也有用基于rho的热物性的，目前还不清楚背后的原因。可压缩求解器（solver/compressible）中，像rhoPimpleFoam本来是采用基于psi的热物性，但高版本开始也把基于rho的热物性揉进去了，很好奇如果用rhoPimpleFoam采用同一套状态方程去求解同一个问题，但一个用psiThermo一个用rhoThermo，结果是否会一样？

啰嗦了两大段，还是回到之前的问题。很明显reactingFoam和rhoReactingFoam两个求解器调用的压力方程形式不同是因为密度更新的方式不同，而密度更新方式的不同很可能来源于是否采用低马赫数假设这一套计算方法。这一点现在还不明确，所以发上来希望有大佬能给讲讲

另外有个问题，如果采用低马赫数近似，密度随压力不变，那是不是就相当于认为声速无穷大？这样假设会导致什么呢？

yuky

@cccrrryyy 我想请问一下您对这个“如果采用低马赫数近似，密度随压力不变，那是不是就相当于认为声速无穷大？这样假设会导致什么呢？”有个结论了吗，对于发动机来说，燃料进入燃烧室通常是速度较快的，但燃料开始燃烧之后速度又是正常低马赫数，请问这种用reactingFoam合适吗

myheart

@cccrrryyy 大佬，请教个问题，最近这两个求解器也困惑了我许久，拜读这个帖子后仍有疑问。
问题1：按照李老师所说，声速较大会导致rhoReactingFoam中采用的压力方程无法求解，那是不是意味着reactingFoam适用于声速可能无穷大也就是低马赫数时密度不随压力改变的情况？
问题2：同时由于rhoReactingFoam采用rhoPimpleFoam的压力方程，其为压力基求解器，捕捉激波的能力没有rhoCentralFoam强，所以其实际上也不适用于高马赫数流动？
问题3：既然这俩都不能计算高速可压流，那基于psi和基于rho的热物性带来的差异主要体现在哪？

CFD中文网

OpenFOAM的两个燃烧求解器