运行结束后输出密度rho

[风云5091]

@wangfei9088 前辈您好，最近在学习reactingFoam源码，请问reactingFoam求解器的不可压缩性体现在什么方面呀？它的连续性方程采取的也是可压缩的那一套，我理解的不可压缩性可能指的是某种物质的密度不随温度变化，比如甲烷等等？而每个网格的密度由于网格中不同组分的质量分数不同所以必然会发生改变，不知我的猜想是否对吗，没有在源码中找到相关的说明

[风云5091]

@wangfei9088 前辈还有就是openfoam中如果要实现输出单个反应的反应速率应该怎么实现呀，我看很多文献中的敏感性分析都用到了这个，从算完的结果中能直接导出来吗

[wangfei9088]

@风云5091 可压缩指的是处理间断吧，超音速，激波等等，可以看东岳老师写的解析。
reactingFoam里混合物密度的处理是每个组分的理想气体状态方程得到的密度之和。
$$
\begin{equation}
\rho = \sum_{k=1}^N \rho_k
\end{equation}
$$
以OpenFOAM-v2012为例，perfectGasI.H 计算每个组分的密度：

 template<class Specie>
 inline Foam::scalar Foam::perfectGas<Specie>::rho(scalar p, scalar T) const
 {
     return p/(this->R()*T);
 }

混合物的特性（包括密度）是在multiComponentMixture.H 里计算的，四个函数里分别乘以质量分数：

         const ThermoType& cellMixture(const label celli) const;
  
         const ThermoType& patchFaceMixture
         (
             const label patchi,
             const label facei
         ) const;
  
         const ThermoType& cellVolMixture
         (
             const scalar p,
             const scalar T,
             const label celli
         ) const;
  
         const ThermoType& patchFaceVolMixture
         (
             const scalar p,
             const scalar T,
             const label patchi,
             const label facei
         ) const;

混合物中每个组分的密度也可以在SpecieMixture.C 中获取。

 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::rho
 (
     const label speciei,
     const scalar p,
     const scalar T
 ) const
 {
     return this->getLocalThermo(speciei).rho(p, T);
 }

最后，网格和边界里的密度是在哪儿计算并更新的呢？hePsiThermo.C或heRhoThermo.C。

这些都是源码里提供的，供参考。

[wangfei9088]

@风云5091 OpenFOAM-v2012中有个reactionsSensitivityAnalysis ，具体是在controlDict里加一段：

functions
{
    sensitivityAnalysis
    {
        libs                (fieldFunctionObjects);
        type                psiReactionsSensitivityAnalysis;
        writeControl        writeTime;
    }
}

可能对你有用，供参考。

[风云5091]

@wangfei9088 好的好的，谢谢前辈！！

[风云5091]

@wangfei9088 前辈这个我之前说的有点问题，我之前以为敏感性分析就是分析每个基元反应的速率，结果好像不是！我目前就是想得到计算收敛后每个基元反应的速率，然后和流场结果对比一下，确定在哪块流场区域哪个基元反应相对更重要

[wangfei9088]

@风云5091 我的知识也有限，没见过后处理有这个计算基元反应速率的功能。提供个思路。
StandardChemistryModel.H 里有计算某个反应中某个组分反应速率的函数：

             //- Return reaction rate of the speciei in reactionI
             virtual tmp<volScalarField::Internal> calculateRR
             (
                 const label reactionI,
                 const label speciei
             ) const;

specieReactionRates.C 用到了上面那个函数，得到的是总的反应速率：

             volScalarField::Internal RR
             (
                 chemistryModel_.calculateRR(ri, si)
             );

我想应该可以仿照它写一个满足要求的functionObjects吧。供参考。

[风云5091]

@wangfei9088 好的，感谢前辈提供的思路，基元反应速率在我研究的这个领域还挺重要的，至少我看的大部分文献里都提到了这个，暑假可能要多学习学习openfoam代码了

[风云5091]

@wangfei9088 如果可以的话，前辈您可以提供一些您当时学习燃烧这块openfoam源代码的思路/资料吗，感觉这块李老师研究的也比较少，网上能搜索到的资源也比较少。

[wangfei9088]

@风云5091 啊，我也没有特别的感悟，燃烧模拟的理论部分当时看的比较多的书是Theoretical and numerical combustion，方程和公式推导很详细。有大佬在网站上写了OpenFOAM里关于燃烧的一些代码解析，结合多看OpenFOAM源代码学习。如果学习中有疑问，可以随时在CFD中文网里提呀。论坛里卧虎藏龙，肯定有大佬帮忙解答的。

[风云5091]

@wangfei9088 好的好的，谢谢前辈！！

[李东岳]

@wangfei9088

大佬，我跑个题。我看你那面做化学反应比较多。有没有兴趣用PINN搞。化学反应用PINN，其他部分还是用传统CFD。
这个路子你们组里面搞了没。没搞的话感兴趣搞么。

[wangfei9088]

@李东岳 哈哈，李老师，我不是大佬，我们组里没人用PINN研究化学反应。看论坛里讨论得这么火热，我也想学习学习PINN的，就是还不熟悉。搞这个要向李老师多学习呀！

[尚善若水]

@wangfei9088 请教大佬，参照specieReactionRates的写法，先尝试输出每个反应步的体积平均Qdot(命名sumVQdotRi/V)，计算Qdot时需要物种的比焓h[speciei]，这个如何访问呢？

const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？

bool Foam::functionObjects::reactionSpecieQdots::write()
{
    logFiles::write();

    const label nSpecie = chemistryModel_.nSpecie();
    const label nReaction = chemistryModel_.nReaction();

    // Region volume
    const scalar V = this->V();

    for (label reactioni=0; reactioni<nReaction; reactioni++)
    {
        if (Pstream::master())
        {
            writeTime(file());
            file() << token::TAB << reactioni;
        }

        const PtrList<volScalarField::Internal> RR
        (
            chemistryModel_.reactionRR(reactioni)
        );

        scalar sumVQdotsR = 0;

        for (label speciei=0; speciei<nSpecie; speciei++)
        {
            scalar sumVQdotRi = 0;

            const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？

            if (all())
            {
                sumVQdotRi = fvc::domainIntegrate(RR[speciei]*(-h[speciei])).value();
            }
            else
            {
                sumVQdotRi =
                    gSum
                    (
                        scalarField
                        (
                            fvMeshFunctionObject::mesh_.V()*RR[speciei]*(-h[speciei]),
                            cells()
                        )
                    );
            }

            if (Pstream::master())
            {
                file() << token::TAB << sumVQdotRi/V;
            }

            sumVQdotsR += sumVQdotRi;
        }

        if (Pstream::master())
        {
            file() << token::TAB << sumVQdotsR/V;
        }

        if (Pstream::master())
        {
            file() << nl;
        }
    }

    if (Pstream::master())
    {
        file() << nl << endl;
    }

    return true;
}

[wangfei9088]

@尚善若水 没太看懂const scalar h[speciei] = specieThermos_h[speciei].hf(); //如何访问？这是啥意思。
如果需要得到每个组分的化学焓（生成焓），可以这么写：

const scalar hi = specieThermo_[speciei].Hc(); // 化学焓 J/kg

OpenFOAM版本不同，可能有些许差异，我用的是OpenFOAM-v2012。

[尚善若水]

@wangfei9088 这是我看到基金会版本里面用来算Qdot的，我直接贴上去了。我的问题是specieThermo_这个可以直接访问吗？const label nSpecie = chemistryModel_.nSpecie();可以这么写是前面定义了

我看到specieThermo_是在别的地方定义的。我弄不明白这个。

[wangfei9088]

@尚善若水 看了一下，开发版和of2012差别很大，开发版里不能用specieThermo_了。依据要不要单位，组分的生成焓在开发版里可以这么写：

dimensionedScalar hi = chemistryModel_.thermo().hfi(speciei); //- Enthalpy of formation [J/kg]
或者
scalar hi = chemistryModel_.thermo().hfiValue(speciei); //- Enthalpy of formation [J/kg]

[尚善若水]

@wangfei9088 好的，我仔细看下。谢谢。

[尚善若水]

@wangfei9088 感谢大佬，程序正常输出，正确性还没测试，也许可以和cantera的结果比一下，代码就是上面发的加上根据大佬说的改的。另外，我在想可不可以更进一步，写出每一个反应步每个组分的Qdot，命名为Qdot.Ri.Speciei，类似输出RR.Speciei?

[wangfei9088]

@尚善若水 应该可以吧。RR[speciei][celli]是celli网格里speciei的反应速率，那该网格里该组分的Qdot应该可以这么写：

scalar Qdot-speciei-celli = hi * RR_[speciei][celli];