运行结束后输出密度rho

[风云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];