关于在求解器中使用热物理库中的函数

huangyuhui723

各位老师好，我目前在修改一个求解器，想要在求解器中使用热物理库中的HE函数对he的边界值进行更新，我在thermoI.h中看到了这样的定义：

template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::HE(const scalar p, const scalar T) const
{
    return Type<thermo<Thermo, Type>>::HE(*this, p, T);
}

和heThermo.c里这样的使用格式：

&MixtureType::thermoMixtureType::HE

看的有点懵，不知道应该以什么样的格式来调用这个函数，希望有知道的老师能解答一下疑惑，感谢。

wangfei9088

@huangyuhui723 以of10为例。

对的，he是通过heThermo.C里这个函数构造的。

     he_
     (
         IOobject
         (
             BasicThermo::phasePropertyName
             (
                 MixtureType::thermoType::heName(),
                 phaseName
             ),
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         volScalarFieldProperty //这个函数确定cell和patch的值
         (
             "he",
             dimEnergy/dimMass,
             &MixtureType::cellThermoMixture,      //cell
             &MixtureType::patchFaceThermoMixture, //patch
             &MixtureType::thermoMixtureType::HE,
             this->p_,
             this->T_
         ),
         this->heBoundaryTypes(),
         this->heBoundaryBaseTypes()
     ),

能看到he的边界是通过这个函数确定的：

 template<class BasicThermo, class MixtureType>
 template<class PatchFaceMixture, class Method, class ... Args>
 Foam::tmp<Foam::scalarField>
 Foam::heThermo<BasicThermo, MixtureType>::patchFieldProperty
 (
     PatchFaceMixture patchFaceMixture,
     Method psiMethod,
     const label patchi,
     const Args& ... args
 ) const
 {
     tmp<scalarField> tPsi
     (
         new scalarField(this->T_.boundaryField()[patchi].size())
     );
     scalarField& psi = tPsi.ref();
 
     forAll(this->T_.boundaryField()[patchi], facei)
     {
         psi[facei] =
             ((this->*patchFaceMixture)(patchi, facei).*psiMethod)
             (
                 args[facei] ...
             );  //计算he边界上的值
     }
 
     return tPsi;
 }

关键是这个patchFaceMixture函数，调用的是coefficientMultiComponentMixture.C里的：

 template<class ThermoType>
 const typename
 Foam::coefficientMultiComponentMixture<ThermoType>::thermoMixtureType&
 Foam::coefficientMultiComponentMixture<ThermoType>::patchFaceThermoMixture
 (
     const label patchi,
     const label facei
 ) const
 {
     mixture_ =
         this->Y()[0].boundaryField()[patchi][facei]
        *this->specieThermos()[0];
 
     for (label i=1; i<this->Y().size(); i++)
     {
         mixture_ +=
             this->Y()[i].boundaryField()[patchi][facei]
            *this->specieThermos()[i];
     }
 
     return mixture_;
 }

可以看出混合物的he的边界值是通过上面这个函数更新的，即各个组分的质量分数加权求得的。如果需要更改边界值，修改这个函数。但是注意，直接修改这个函数会改变所有特性的值包括he。

speciesThermos()[i]指的是组分i的特性，OpenFOAM源代码中每个组分的焓值调用的是thermoI.H里的：

 template<class Thermo, template<class> class Type>
 inline Foam::scalar
 Foam::species::thermo<Thermo, Type>::HE(const scalar p, const scalar T) const
 {
     return Type<thermo<Thermo, Type>>::HE(*this, p, T);
 }

调用的是sensibleEnthalpy.H里的：

             scalar HE
             (
                 const Thermo& thermo,
                 const scalar p,
                 const scalar T
             ) const
             {
                 return thermo.Hs(p, T);
             }

调用的是janafThermoI.H里的：

 template<class EquationOfState>
 inline Foam::scalar Foam::janafThermo<EquationOfState>::Hs //显焓
 (
     const scalar p,
     const scalar T
 ) const
 {
     return Ha(p, T) - Hf();
 }
 
 
 template<class EquationOfState>
 inline Foam::scalar Foam::janafThermo<EquationOfState>::Ha //总焓
 (
     const scalar p,
     const scalar T
 ) const
 {
     const coeffArray& a = coeffs(T);
     return
     (
         ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
       + a[5]
     ) + EquationOfState::H(p, T);
 }
 
 
 template<class EquationOfState>
 inline Foam::scalar Foam::janafThermo<EquationOfState>::Hf() const //生成焓
 {
     const coeffArray& a = lowCpCoeffs_;
     return
     (
         (
             (((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
           + a[0]
         )*Tstd + a[5]
     );
 }

也就是说，如果要更改每个组分的特性，需要改变不同的janaf系数，这一般是不可能这么做的。当然，可以使用不同的thermo模型。

上面仅以coefficientMultiComponentMixture和janaf为例说明这些函数的调用关系。OpenFOAM0-10中很多调用关系，比如&MixtureType::thermoMixtureType::HE，理解会有点难。

如果觉得OpenFOAM0-10中这些关系难以理解，建议看看ESI版本的OpenFOAM（比如OpenFOAM-v2012）和基金会以前的版本（比如OpenFOAM-6和OpenFOAM-7），这些版本的代码直接很多，但意思是一样的，会容易理解很多。

供参考。

huangyuhui723

@wangfei9088 感谢老师，我参照您的思路捋了一下，把相关类实例化之后能够使用HE函数计算了，但是还发现一个问题，我在求解器中通过其它计算更新了T之后想对he赋值，但一直赋值不上，尝试了直接=和forAll都不行，是因为he指向的thermo.he()会在热物理库调用中自行赋值么，请问您知道对he赋值的正确方法么？

wangfei9088

@huangyuhui723 我明白你的意思了。

@huangyuhui723 在关于在求解器中使用热物理库中的函数中说：

通过其它计算更新了T之后想对he赋值

如果是求解TEqn，可以是得到温度T后更新he。OpenFOAM里的求解器一般都是求解EEqn，得到he后通过thermo.correct()更新T。

如果已知边界温度T后要对he赋值，应该在求解EEqn之前对he边界赋值，然后求解EEqn。大致代码是这样：

volScalarField& he = thermo.he();


//****************************************//加这一段
forAll(p.boundaryField(), patchi)
{
    forAll(p.boundaryField()[patchi], facei)
    {
        he.boundaryFieldRef()[patchi][facei] = function(p.boundaryField()[patchi][facei], T.boundaryField()[patchi][facei]);
    }
}
//****************************************//

     fvScalarMatrix EEqn
     (
        ......
     );

这个function函数就是T与he的关系，类似上面，混合物的he通过各组分的显焓的质量分数加权得到，显焓用janaf里的关系确定。

供参考。

huangyuhui723

@wangfei9088 感谢老师，不好意思最近没来得及回复。我后面还是尝试了直接写了TEqn，没有去折腾热物理库了。很感谢老师的指点。

CFD中文网

关于在求解器中使用热物理库中的函数