OpenFOAM热物理类代码详解

喜欢ヅ旅行 2022-10-27 15:55 423阅读 0赞

## 热物理类的RTS机制 ##

OpenFOAM利用RTS机制来动态创建模型实例，核心思想在于利用C++中静态成员会在类的实例化之前创建，  
以`buoyantSimpleFoam.C`对于生成热物理类的代码的为例，其中

rhoThermo::New(mesh)
    
    作者网址：https://blog.csdn.net/dsfsdffgfd

调用rhoThermo类中静态成员函数:

Foam::autoPtr<Foam::rhoThermo> Foam::rhoThermo::New
      (
          const fvMesh& mesh,
          const word& phaseName
      )
      { 
          return basicThermo::New<rhoThermo>(mesh, phaseName);
      }

返回了基类basicThermo中名为New的静态模板类，在RTS机制中也被称为Selector,其函数的定义如下:

template<class Thermo>
      Foam::autoPtr<Thermo> Foam::basicThermo::New
      (
          const fvMesh& mesh,
          const word& phaseName
      )
      { 
          IOdictionary thermoDict
          (
              IOobject
              (
                  phasePropertyName(dictName, phaseName),
                  mesh.time().constant(),
                  mesh,
                  IOobject::MUST_READ_IF_MODIFIED,
                  IOobject::NO_WRITE,
                  false
              )
          );
    
          auto cstrIter =
              lookupThermo<Thermo, typename Thermo::fvMeshConstructorTable>
              (
                  thermoDict,
                  Thermo::fvMeshConstructorTablePtr_
              );
    
          return autoPtr<Thermo>(cstrIter()(mesh, phaseName));
      }
    
      template<class Thermo, class Table>
      typename Table::iterator Foam::basicThermo::lookupThermo
      (
          const dictionary& thermoDict,
          Table* tablePtr
      )
      { 
          if (thermoDict.isDict("thermoType"))
          { 
              const dictionary& thermoTypeDict = thermoDict.subDict("thermoType");
    
              Info<< "Selecting thermodynamics package " << thermoTypeDict << endl;
    
              if (thermoTypeDict.found("properties"))
              { 
                  std::initializer_list<const char*> cmptNames
                  { 
                      "type",
                      "mixture",
                      "properties",
                      "energy"
                  };
    
                  // Construct the name of the thermo package from the components
                  const word thermoTypeName
                  (
                      thermoTypeDict.get<word>("type") + '<'
                    + thermoTypeDict.get<word>("mixture") + '<'
                    + thermoTypeDict.get<word>("properties") + ','
                    + thermoTypeDict.get<word>("energy") + ">>"
                  );
    
                  return lookupThermo<Thermo, Table>
                  (
                      thermoTypeDict,
                      tablePtr,
                      cmptNames,
                      thermoTypeName
                  );
              }
              else
              { 
                  std::initializer_list<const char*> cmptNames
                  { 
                      "type",
                      "mixture",
                      "transport",
                      "thermo",
                      "equationOfState",
                      "specie",
                      "energy"
                  };
    
                  // Construct the name of the thermo package from the components
                  const word thermoTypeName
                  (
                      thermoTypeDict.get<word>("type") + '<'
                    + thermoTypeDict.get<word>("mixture") + '<'
                    + thermoTypeDict.get<word>("transport") + '<'
                    + thermoTypeDict.get<word>("thermo") + '<'
                    + thermoTypeDict.get<word>("equationOfState") + '<'
                    + thermoTypeDict.get<word>("specie") + ">>,"
                    + thermoTypeDict.get<word>("energy") + ">>>"
                  );
    
                  return lookupThermo<Thermo, Table>
                  (
                      thermoTypeDict,
                      tablePtr,
                      cmptNames,
                      thermoTypeName
                  );
              }
          }
          else
          { 
              const word thermoTypeName(thermoDict.get<word>("thermoType"));
    
              Info<< "Selecting thermodynamics package " << thermoTypeName << endl;
    
              // Table iterator, not const_iterator
              auto cstrIter = tablePtr->find(thermoTypeName);
    
              if (!cstrIter.found())
              { 
                  FatalIOErrorInLookup
                  (
                      thermoDict,
                      Thermo::typeName,
                      thermoTypeName,
                      *tablePtr
                  ) << exit(FatalIOError);
              }
    
              return cstrIter;
          }
      }

首先定义了读写的字典`IOdictionary thermoDict`，其中`dictName`为静态变量定义为:`const Foam::word Foam::basicThermo::dictName("thermophysicalProperties");` ,下一步是调用了静态模板函数`lookupThermo`，其模板参数为`rhoThermo::fvMeshConstructorTable`，也就是返回值`rhoThermo::fvMeshConstructorTable::iterator`，参数为`Thermo::fvMeshConstructorTablePtr_`。`lookupThermo`函数先读取thermophysicalProperties文件，三种字典形式通过上述代码合为一体，例如字典经过读入写成:  
`hePsiThermo<pureMixture<const<hConst<perfectGas<specie>>,sensibleEnthalpy>>>`，或者直接读入这种形式  
`thermoType hePsiThermo<pureMixture<const<hConst<perfectGas<specie>>,sensibleEnthalpy>>>;`，把这个值赋给`thermoTypeName`后返回`rhoThermo::fvMeshConstructorTablePtr_->find(thermoTypeName)`。这段代码的HashTable是通过如下宏定义:

declareRunTimeSelectionTable
        (
            autoPtr,
            rhoThermo,
            fvMesh,
            (const fvMesh& mesh, const word& phaseName),
            (mesh, phaseName)
        )
    
      #define declareRunTimeSelectionTable(autoPtr,baseType,argNames,argList,parList)\ \ /* Construct from argList function pointer type */                         \
        typedef autoPtr<baseType> (*argNames##ConstructorPtr)argList;              \
                                                                                   \
        /* Construct from argList function table type */                           \
        typedef HashTable<argNames##ConstructorPtr, word, string::hash>            \
            argNames##ConstructorTable;                                            \
                                                                                   \
        /* Construct from argList function pointer table pointer */                \
        static argNames##ConstructorTable* argNames##ConstructorTablePtr_;         \
                                                                                   \
        /* Table constructor called from the table add function */                 \
        static void construct##argNames##ConstructorTables();                      \
                                                                                   \
        /* Table destructor called from the table add function destructor */       \
        static void destroy##argNames##ConstructorTables();                        \
                                                                                   \
        /* Class to add constructor from argList to table */                       \
        template<class baseType##Type>                                             \
        class add##argNames##ConstructorToTable                                    \
        {                                                                           \
        public:                                                                    \
                                                                                   \
            static autoPtr<baseType> New argList                                   \
            {                                                                       \
                return autoPtr<baseType>(new baseType##Type parList);              \
            }                                                                      \
                                                                                   \
            explicit add##argNames##ConstructorToTable                             \
            (                                                                      \
                const word& lookup = baseType##Type::typeName                      \
            )                                                                      \
            {                                                                       \
                construct##argNames##ConstructorTables();                          \
                if (!argNames##ConstructorTablePtr_->insert(lookup, New))          \
                {                                                                   \
                    std::cerr<< "Duplicate entry " << lookup                       \
                        << " in runtime selection table " << #baseType             \
                        << std::endl;                                              \
                    error::safePrintStack(std::cerr);                              \
                }                                                                  \
            }                                                                      \
                                                                                   \
            ~add##argNames##ConstructorToTable()                                   \
            {                                                                       \
                destroy##argNames##ConstructorTables();                            \
            }                                                                      \
                                                                                   \
            add##argNames##ConstructorToTable                                      \
                (const add##argNames##ConstructorToTable&) = delete;               \
                                                                                   \
            void operator=                                                         \
                (const add##argNames##ConstructorToTable&) = delete;               \
        };                                                                         \
                                                                                   \
        /* Class to add constructor from argList to table */                       \
        /* Remove only the entry (not the table) upon destruction */               \
        template<class baseType##Type>                                             \
        class addRemovable##argNames##ConstructorToTable                           \
        {                                                                           \
        public:                                                                    \
                                                                                   \
            const word name;  /* Lookup name for later removal */                  \
                                                                                   \
            static autoPtr<baseType> New argList                                   \
            {                                                                       \
                return autoPtr<baseType>(new baseType##Type parList);              \
            }                                                                      \
                                                                                   \
            explicit addRemovable##argNames##ConstructorToTable                    \
            (                                                                      \
                const word& lookup = baseType##Type::typeName                      \
            )                                                                      \
            :                                                                      \
                name(lookup)                                                       \
            {                                                                       \
                construct##argNames##ConstructorTables();                          \
                argNames##ConstructorTablePtr_->set(lookup, New);                  \
            }                                                                      \
                                                                                   \
            ~addRemovable##argNames##ConstructorToTable()                          \
            {                                                                       \
                if (argNames##ConstructorTablePtr_)                                \
                {                                                                   \
                    argNames##ConstructorTablePtr_->erase(name);                   \
                }                                                                  \
            }                                                                      \
                                                                                   \
            addRemovable##argNames##ConstructorToTable                             \
                (const addRemovable##argNames##ConstructorToTable&) = delete;      \
                                                                                   \
            void operator=                                                         \
                (const addRemovable##argNames##ConstructorToTable&) = delete;      \
        };

上面代码翻译一下:

/* Construct from argList function pointer type */                         
        typedef autoPtr<rhoThermo> (*fvMeshConstructorPtr)(const fvMesh& mesh, const word& phaseName);              
                                                                                   
        /* Construct from argList function table type */                           
        typedef HashTable<fvMeshConstructorPtr, word, string::hash> fvMeshConstructorTable;                                            
                                                                                   
        /* Construct from argList function pointer table pointer */                
        static fvMeshConstructorTable* fvMeshConstructorTablePtr_;         
                                                                                   
        /* Table constructor called from the table add function */                 
        static void constructfvMeshConstructorTables();                      
                                                                                   
        /* Table destructor called from the table add function destructor */       
        static void destroyfvMeshConstructorTables();                        
                                                                                   
        /* Class to add constructor from argList to table */                       
        template<class rhoThermoType>                                             
        class addfvMeshConstructorToTable                                    
        {                                                                           
        public:                                                                    
                                                                                   
            static autoPtr<rhoThermo> New (const fvMesh& mesh, const word phaseName)
            {                                                                       
                return autoPtr<rhoThermo>(new rhoThermoType (mesh, phaseName));
            }                                                                      
    
            explicit addfvMeshConstructorToTable
            (                                                                      
                const word& lookup = rhoThermoType::typeName
            )                                                                      
            {                                                                       
                constructfvMeshConstructorTables();
                if (!fvMeshConstructorTablePtr_->insert(lookup, New))
                { 
                    std::cerr<< "Duplicate entry " << lookup
                        << " in runtime selection table " << "rhoThermo"
                        << std::endl;
                    error::safePrintStack(std::cerr);
                }
            }
                                                                                   
            ~addfvMeshConstructorToTable()                                   
            {                                                                       
                destroyfvMeshConstructorTables();                            
            }                                                                      
                                                                                   
            addfvMeshConstructorToTable                                      
                (const addfvMeshConstructorToTable&) = delete;               
                                                                                   
            void operator=                                                         
                (const addfvMeshConstructorToTable&) = delete;               
        };                                                                         
                                                                                   
        /* Class to add constructor from argList to table */                       
        /* Remove only the entry (not the table) upon destruction */               
        template<class rhoThermoType>                                             
        class addRemovablefvMeshConstructorToTable                           
        {                                                                           
        public:                                                                    
                                                                                   
            const word name;  /* Lookup name for later removal */                  
                                                                                   
            static autoPtr<rhoThermo> New (const fvMesh& mesh, const word phaseName)                                   
            {                                                                       
                return autoPtr<rhoThermo>(new rhoThermoType (mesh, phaseName));              
            }                                                                      
                                                                                   
            explicit addRemovablefvMeshConstructorToTable                    
            (                                                                      
                const word& lookup = rhoThermoType::typeName                      
            )                                                                      
            :                                                                      
                name(lookup)                                                       
            {                                                                       
                constructfvMeshConstructorTables();                          
                fvMeshConstructorTablePtr_->set(lookup, New);                  
            }                                                                      
                                                                                   
            ~addRemovablefvMeshConstructorToTable()                          
            {                                                                       
                if (fvMeshConstructorTablePtr_)                                
                {                                                                   
                    fvMeshConstructorTablePtr_->erase(name);                   
                }                                                                  
            }                                                                      
                                                                                   
            addRemovablefvMeshConstructorToTable                             
                (const addRemovablefvMeshConstructorToTable&) = delete;      
                                                                                   
            void operator=                                                         
                (const addRemovablefvMeshConstructorToTable&) = delete;      
        };

上端代码主要是创建了HashTable，热物理类的独特性在于HashTable的元素不是在派生类中添加，而是通过下面宏进行定义:

makeThermos
    (
        rhoThermo,
        heRhoThermo,
        pureMixture,
        constTransport,
        sensibleEnthalpy,
        hConstThermo,
        perfectGas,
        specie
    );

这里直接对上面的宏进行翻译得到如下:

// **************************************步骤一*****************************************************
      typedef                                                                    
          constTransport                                                              
          <                                                                      
              species::thermo                                                    
              <                                                                  
                  hConstThermo                                                         
                  <                                                              
                      perfectGas                                                 
                      <                                                          
                          Specie                                                 
                      >                                                          
                  >,                                                             
                  sensibleEnthalpy                                                           
              >                                                                  
          >                                                                      
          constTransportsensibleEnthalpyhConstThermoperfectGasSpecie
                                                                                  
      typedef                                                                    
        heRhoThermo                                                                
        <                                                                      
            rhoThermo,                                                        
            pureMixture<constTransportsensibleEnthalpyhConstThermoperfectGasSpecie>                         
        >                                                                      
        heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie;
    
        // **************************************步骤二*****************************************************
        template<>
        const ::Foam::word heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie::typeName("heRhoThermo<pureMixture<constTransportsensibleEnthalpyhConstThermoperfectGasSpecie>>")
    
        // **************************************步骤三*****************************************************
        addToRunTimeSelectionTable                                                 
        (                                                                          
            basicThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,
            fvMesh
        );
        addToRunTimeSelectionTable
        (
            basicThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,
            fvMeshDictPhase
        );
        addToRunTimeSelectionTable                                                 
        (                                                                          
            fluidThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,
            fvMesh
        );
        addToRunTimeSelectionTable                                                 
        (                                                                          
            fluidThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,
            fvMeshDictPhase
        );
        addToRunTimeSelectionTable                                                 
        (                                                                          
            rhoThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,                                              
            fvMesh
        );
        addToRunTimeSelectionTable                                                 
        (                                                                          
            rhoThermo,
            heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie,                                              
            fvMeshDictPhase
        );

上面的代码步骤一先用`typedef`给一定模板参数下`constTransport`和`heRhoThermo`的别称，再在步骤二中对这一参数下的模板类特例话其静态变量，最后的步骤三利用`addToRunTimeSelectionTable`宏命令，在`rhoThermo,fluidThermo,basicThermo`的HashTable中添加上面参数的模板类，再次把上述宏命令展开:

rhoThermo::addfvMeshConstructorToTable<heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie>
        addthisTypeheRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecieConstructorTorhoThermoTable_;

终于，再次回到了`rhoThermo`中通过宏定义模板类 `addfvMeshConstructorToTable` ，实例化了一个`addfvMeshConstructorToTable<heRhoThermopureMixtureconstTransportsensibleEnthalpyhConstThermoperfectGasSpecie>`对象，在其构造函数中，首先调用了`constructfvMeshConstructorTables`函数，这一函数通过如下的宏进行定义:

#define defineRunTimeSelectionTableConstructor(baseType,argNames) \ \ /* Table constructor called from the table add function */                 \
        void baseType::construct##argNames##ConstructorTables()                    \
        {                                                                           \
            static bool constructed = false;                                       \
            if (!constructed)                                                      \
            {                                                                       \
                constructed = true;                                                \
                baseType::argNames##ConstructorTablePtr_                           \
                    = new baseType::argNames##ConstructorTable;                    \
            }                                                                      \
        }
    #define defineRunTimeSelectionTableDestructor(baseType,argNames) \ \ /* Table destructor called from the table add function destructor */       \
        void baseType::destroy##argNames##ConstructorTables()                      \
        {                                                                           \
            if (baseType::argNames##ConstructorTablePtr_)                          \
            {                                                                       \
                delete baseType::argNames##ConstructorTablePtr_;                   \
                baseType::argNames##ConstructorTablePtr_ = nullptr;                \
            }     
        }
                                                                         \
    // Create pointer to hash-table of functions
    #define defineRunTimeSelectionTablePtr(baseType,argNames) \ \ /* Define the constructor function table */                                \
        baseType::argNames##ConstructorTable*                                      \
            baseType::argNames##ConstructorTablePtr_(nullptr)

展开后写为:

/* Table constructor called from the table add function */                 
      void rhoThermo::constructfvMeshConstructorTables()                    
      {                                                                           
          static bool constructed = false;                                       
          if (!constructed)                                                      
          {                                                                       
              constructed = true;                                                
              rhoThermo::fvMeshConstructorTablePtr_                           
                  = new rhoThermo::fvMeshConstructorTable;                    
          }                                                                      
      }
    
      void rhoThermo::destroyfvMeshConstructorTables()                      
      {                                                                           
            if (rhoThermo::fvMeshConstructorTablePtr_)                          
            {                                                                       
                delete rhoThermo::fvMeshConstructorTablePtr_;                   
                rhoThermo::fvMeshConstructorTablePtr_ = nullptr;                
            }     
      }
    
      rhoThermo::fvMeshConstructorTable*                                      \
            rhoThermo::fvMeshConstructorTablePtr_(nullptr)

上面的`constructed`为局部静态变量，在计算中只会初始化一次，不会重复调用，也就是对于`rhoThermo`只会生成一个`fvMeshConstructorTablePtr_`，接下来调用了  
`fvMeshConstructorTablePtr_->insert(lookup, New)`，回到HashTable的定义:

typedef autoPtr<rhoThermo> (*fvMeshConstructorPtr)(const fvMesh& mesh, const word& phaseName);              
        typedef HashTable<fvMeshConstructorPtr, word, string::hash> fvMeshConstructorTable;

其key为’lookup’即之前特例化的模板名`heRhoThermo<pureMixture<constTransportsensibleEnthalpyhConstThermoperfectGasSpecie>>`，value为返回值`autoPtr<rhoThermo>`参数`const fvMesh& mesh, const word& phaseName`的函数指针，其值为类中的静态函数`static autoPtr<rhoThermo> New (const fvMesh& mesh, const word phaseName)`，具体调动关系可以详见表驱动模式。层层代入最后返回到`buoyantSimple.C`的代码就是:

autoPtr<rhoThermo>(autoPtr<rhoThermo>(
      new heRhoThermo<rhoThermo,
            pureMixture<constTransport<species::thermo<hConstThermo<perfectGas<Specie>>,
            sensibleEnthalpy>>>> (mesh, phaseName)))

也就是:

autoPtr<rhoThermo> pThermo(
        new heRhoThermo<rhoThermo,
            pureMixture<constTransport<species::thermo<hConstThermo<perfectGas<Specie>>,sensibleEnthalpy>>>> 
            (mesh, null));

--------------------

## 热物理类的代码实现 ##

通过热物理类的RTS机制代码，同样以字典为例最终代码得到了指向`heRhoThermo<rhoThermo,pureMixture<constTransport<species::thermo<hConstThermo<perfectGas<Specie>>,sensibleEnthalpy>>>>`实例的指针，首先看看heRhoThermo的代码实现上的区别:

--------------------

# Appendix I #

能量方程可以表示为:  
 ρ D E D t ⏞ 总能 = − ∇ ⋅ ( p U ) ⏞ 压力功 − ∇ ⋅ q ˙ ′ ′ ⏞ 热流 + ∇ ⋅ ( τ ⋅ U ) ⏞ 粘性力功 + ρ g ⋅ U ⏞ 重力功 + ρ q ˙ ′ ′ ′ ⏞ 热源 (1) \\overbrace\{\\rho\\frac\{DE\}\{Dt\}\}^\{\\textbf\{总能\}\}= \\overbrace\{-\\nabla\\cdot(p\\bold\{U\})\}^\{\\textbf\{压力功\}\} -\\overbrace\{\\nabla\\cdot\\bold\{\\dot\{q\}''\}\}^\{\\textbf\{热流\}\} +\\overbrace\{\\nabla\\cdot(\\bold\{\\tau\}\\cdot\\bold\{U\})\}^\{\\textbf\{粘性力功\}\} +\\overbrace\{\\rho\\bold\{g\}\\cdot\\bold\{U\}\}^\{\\textbf\{重力功\}\} +\\overbrace\{\\rho \\dot\{q\}'''\}^\{\\textbf\{热源\}\} \\tag\{1\} ρDtDE总能=−∇⋅(pU)压力功−∇⋅q˙′′热流\+∇⋅(τ⋅U)粘性力功\+ρg⋅U重力功\+ρq˙′′′热源(1)

(1)为能量方程的非守恒形式，通过随体导数的展开公式 D D t ≡ ∂ ∂ t + U ⋅ ∇ \\frac\{D\}\{Dt\}\\equiv\\frac\{\\partial\}\{\\partial t\}+\\bold\{U\}\\cdot\\nabla DtD≡∂t∂\+U⋅∇公式及质量守恒可以把并把任一张量 Q \\bold\{Q\} Q的非守恒形式展开为守恒形式:  
 ρ D Q D t ≡ ∂ ρ Q ∂ t + ∇ ⋅ ( ρ Q U ) \\rho\\frac\{D\\bold\{Q\}\}\{Dt\}\\equiv\\frac\{\\partial\\rho \\bold\{Q\}\}\{\\partial t\}+\\nabla\\cdot(\\rho\\bold\{Q\}\\bold\{U\}) ρDtDQ≡∂t∂ρQ\+∇⋅(ρQU)

把(1)写为守恒形式并把总能量 E E E展开为内能 e e e和机械能 K K K之和得到:

∂ ρ e ∂ t + ∇ ⋅ ( ρ e U ) + ∂ ρ K ∂ t + ∇ ⋅ ( ρ K U ) = − ∇ ⋅ ( p U ) − ∇ ⋅ q ˙ ′ ′ + ∇ ⋅ ( τ ⋅ U ) + ρ g ⋅ U + ρ q ˙ ′ ′ ′ (2) \\frac\{\\partial \\rho e\}\{\\partial t\}+\\nabla\\cdot(\\rho e\\bold\{U\}) +\\frac\{\\partial \\rho K\}\{\\partial t\}+\\nabla\\cdot(\\rho K\\bold\{U\}) = -\\nabla\\cdot(p\\bold\{U\}) -\\nabla\\cdot\\bold\{\\dot\{q\}''\} +\\nabla\\cdot(\\bold\{\\tau\}\\cdot\\bold\{U\}) +\\rho\\bold\{g\}\\cdot\\bold\{U\} +\\rho \\dot\{q\}''' \\tag\{2\} ∂t∂ρe\+∇⋅(ρeU)\+∂t∂ρK\+∇⋅(ρKU)=−∇⋅(pU)−∇⋅q˙′′\+∇⋅(τ⋅U)\+ρg⋅U\+ρq˙′′′(2)

公式(2)即为基于内能形式的能量方程。根据焓的定义:

h = e + p ρ h = e+\\frac\{p\}\{\\rho\} h=e\+ρp

把:

∂ ρ e ∂ t + ∇ ⋅ ( ρ e U ) = ∂ ρ ( h − p ρ ) ∂ t + ∇ ⋅ ( ρ ( h − p ρ ) U ) = ∂ ρ h ∂ t + ∇ ⋅ ( ρ h U ) − ∂ p ∂ t − ∇ ⋅ ( p U ) \\frac\{\\partial \\rho e\}\{\\partial t\}+\\nabla\\cdot(\\rho e\\bold\{U\}) = \\frac\{\\partial \\rho (h-\\frac\{p\}\{\\rho\})\}\{\\partial t\}+\\nabla\\cdot(\\rho (h - \\frac\{p\}\{\\rho\})\\bold\{U\}) =\\frac\{\\partial \\rho h\}\{\\partial t\}+\\nabla\\cdot(\\rho h\\bold\{U\})- \\frac\{\\partial p\}\{\\partial t\} - \\nabla\\cdot(p\\bold\{U\}) ∂t∂ρe\+∇⋅(ρeU)=∂t∂ρ(h−ρp)\+∇⋅(ρ(h−ρp)U)=∂t∂ρh\+∇⋅(ρhU)−∂t∂p−∇⋅(pU)

代入公式(2)可以得到基于焓形式的能量方程:

∂ ρ h ∂ t + ∇ ⋅ ( ρ h U ) + ∂ ρ K ∂ t + ∇ ⋅ ( ρ K U ) = ∂ p ∂ t − ∇ ⋅ q ˙ ′ ′ + ∇ ⋅ ( τ ⋅ U ) + ρ g ⋅ U + ρ q ˙ ′ ′ ′ (3) \\frac\{\\partial \\rho h\}\{\\partial t\}+\\nabla\\cdot(\\rho h\\bold\{U\}) +\\frac\{\\partial \\rho K\}\{\\partial t\}+\\nabla\\cdot(\\rho K\\bold\{U\}) = \\frac\{\\partial p\}\{\\partial t\} -\\nabla\\cdot\\bold\{\\dot\{q\}''\} +\\nabla\\cdot(\\bold\{\\tau\}\\cdot\\bold\{U\}) +\\rho\\bold\{g\}\\cdot\\bold\{U\} +\\rho \\dot\{q\}''' \\tag\{3\} ∂t∂ρh\+∇⋅(ρhU)\+∂t∂ρK\+∇⋅(ρKU)=∂t∂p−∇⋅q˙′′\+∇⋅(τ⋅U)\+ρg⋅U\+ρq˙′′′(3)

内能形式与焓形式的能量方程是完全一致的，只是形式不同，因此在理论上两者没有差别。

对动量方程:

∂ ρ U ∂ t + ∇ ⋅ ( ρ U U ) = ∇ ⋅ τ + ρ g − ∇ p \\frac\{\\partial \\rho \\bold\{U\}\}\{\\partial t\}+\\nabla\\cdot(\\rho \\bold\{U\}\\bold\{U\}) =\\nabla\\cdot\\bold\{\\tau\}+\\rho\\bold\{g\}-\\nabla p ∂t∂ρU\+∇⋅(ρUU)=∇⋅τ\+ρg−∇p  
两边同时点乘速度 U \\bold\{U\} U可以得到机械能守恒方程:

∂ ρ K ∂ t + ∇ ⋅ ( ρ K U ) = ( ∇ ⋅ τ ) ⋅ U + ρ g ⋅ U − ∇ p ⋅ U (4) \\frac\{\\partial \\rho K\}\{\\partial t\}+\\nabla\\cdot(\\rho K\\bold\{U\}) = (\\nabla\\cdot\\bold\{\\tau\})\\cdot\\bold\{U\}+\\rho\\bold\{g\}\\cdot\\bold\{U\}-\\nabla p\\cdot\\bold\{U\} \\tag\{4\} ∂t∂ρK\+∇⋅(ρKU)=(∇⋅τ)⋅U\+ρg⋅U−∇p⋅U(4)  
把公式(4)分别代入(2)(3)并且去除粘性力影响，内能形式的能量方程的源项为 − p ( ∇ ⋅ U ) -p(\\nabla\\cdot\\bold\{U\}) −p(∇⋅U):

∂ ρ e ∂ t + ∇ ⋅ ( ρ e U ) = − p ( ∇ ⋅ U ) − ∇ ⋅ q ˙ ′ ′ + ρ q ˙ ′ ′ ′ (5) \\frac\{\\partial \\rho e\}\{\\partial t\}+\\nabla\\cdot(\\rho e\\bold\{U\}) = -p(\\nabla\\cdot\\bold\{U\}) -\\nabla\\cdot\\bold\{\\dot\{q\}''\} +\\rho \\dot\{q\}''' \\tag\{5\} ∂t∂ρe\+∇⋅(ρeU)=−p(∇⋅U)−∇⋅q˙′′\+ρq˙′′′(5)  
而焓形式的能量方程则替换为 ∂ p ∂ t + U ⋅ ( ∇ p ) \\frac\{\\partial p\}\{\\partial t\} +\\bold\{U\}\\cdot(\\nabla p) ∂t∂p\+U⋅(∇p)：

∂ ρ h ∂ t + ∇ ⋅ ( ρ h U ) = U ⋅ ( ∇ p ) + ∂ p ∂ t − ∇ ⋅ q ˙ ′ ′ + ρ g ⋅ U + ρ q ˙ ′ ′ ′ (6) \\frac\{\\partial \\rho h\}\{\\partial t\}+\\nabla\\cdot(\\rho h\\bold\{U\}) = \\bold\{U\}\\cdot(\\nabla p) +\\frac\{\\partial p\}\{\\partial t\} -\\nabla\\cdot\\bold\{\\dot\{q\}''\} +\\rho\\bold\{g\}\\cdot\\bold\{U\} +\\rho \\dot\{q\}''' \\tag\{6\} ∂t∂ρh\+∇⋅(ρhU)=U⋅(∇p)\+∂t∂p−∇⋅q˙′′\+ρg⋅U\+ρq˙′′′(6)

对比公式(5)(6)，对于不可压流体或者压缩性较小的流体，(5)中的 − p ( ∇ ⋅ U ) -p(\\nabla\\cdot\\bold\{U\}) −p(∇⋅U)会很小，这在数值中具有收敛性的优势，因此对于 ψ \\psi ψ很小的流动状况，可以 ***使用内能形式的能量方程增强计算的收敛性*** ，但收敛后两者的结果没有区别。

另外不同的求解器也会选择忽略(4)(5)中的一些做功项情况，在buoyantSimpleFoam的能量方程中忽略了粘性力的做功项:

fvScalarMatrix EEqn
        (
            fvm::div(phi, he)
          + (
                he.name() == "e"
              ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
              : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
            )
          - fvm::laplacian(turbulence->alphaEff(), he) 
         ==
            rho*(U&g)
          + radiation->Sh(thermo, he)
          + fvOptions(rho, he)
        );

而在rhoSimpleFoam的能量方程中同时忽略了粘性力和重力的做功:

fvScalarMatrix EEqn
        (
            fvm::div(phi, he)
          + (
                he.name() == "e"
              ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
              : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
            )
          - fvm::laplacian(turbulence->alphaEff(), he)
         ==
            fvOptions(rho, he)
        );

在rhoCentralFoam的能量方程主要考虑公式(4)(5)中的粘性力做功:

volScalarField rhoE
        (
            IOobject
            (
                "rhoE",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            rho*(e + 0.5*magSqr(U))
        );
        surfaceScalarField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
                + fvc::dotInterpolate(mesh.Sf(), tauMC)
            )
            & (a_pos*U_pos + a_neg*U_neg)
        );
        solve
        (
            fvm::ddt(rhoE)
            + fvc::div(phiEp)
            - fvc::div(sigmaDotU)
        );

--------------------