INT: reacting*EulerFoam: correct handling of K for energy equation (#2033)

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2015-2021 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -70,13 +70,6 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class BasePhaseModel>
-Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasePhaseModel>
@@ -100,23 +93,32 @@ Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
+    const volScalarField& rho = this->rho();
 
-    const volVectorField U(this->U());
+    const tmp<volVectorField> tU(this->U());
-    const surfaceScalarField alphaPhi(this->alphaPhi());
+    const volVectorField& U(tU());
-    const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
 
-    const volScalarField contErr(this->continuityError());
+    const tmp<surfaceScalarField> talphaPhi(this->alphaPhi());
-    const volScalarField K(this->K());
+    const surfaceScalarField& alphaPhi(talphaPhi());
+
+    const tmp<surfaceScalarField> talphaRhoPhi(this->alphaRhoPhi());
+    const surfaceScalarField& alphaRhoPhi(talphaRhoPhi());
+
+    const tmp<volScalarField> tcontErr(this->continuityError());
+    const volScalarField& contErr(tcontErr());
+
+    tmp<volScalarField> tK(this->K());
+    const volScalarField& K(tK());
 
     volScalarField& he = this->thermo_->he();
 
     tmp<fvScalarMatrix> tEEqn
     (
-        fvm::ddt(alpha, this->rho(), he)
+        fvm::ddt(alpha, rho, he)
       + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
 
-      + fvc::ddt(alpha, this->rho(), K) + fvc::div(alphaRhoPhi, K)
+      + fvc::ddt(alpha, rho, K) + fvc::div(alphaRhoPhi, K)
       - contErr*K
 
       - fvm::laplacian
@@ -135,7 +137,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
         tEEqn.ref() += filterPressureWork
         (
             fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
-          + this->thermo().p()*fvc::ddt(alpha)
+          + (fvc::ddt(alpha) - contErr/rho)*this->thermo().p()
         );
     }
     else if (this->thermo_->dpdt())

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,7 +55,7 @@ class AnisothermalPhaseModel
 :
     public BasePhaseModel
 {
-    // Private member functions
+    // Private Member Functions
 
         //- Optionally filter the pressure work term as the phase-fraction -> 0
         tmp<volScalarField> filterPressureWork
@@ -67,6 +68,7 @@ public:
 
     // Constructors
 
+        //- Construct from phase system and phase name
         AnisothermalPhaseModel
         (
             const phaseSystem& fluid,
@@ -76,7 +78,7 @@ public:
 
 
     //- Destructor
-    virtual ~AnisothermalPhaseModel();
+    virtual ~AnisothermalPhaseModel() = default;
 
 
     // Member Functions

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2015-2021 OpenFOAM Foundation
     Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -212,13 +212,6 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class BasePhaseModel>
-Foam::MovingPhaseModel<BasePhaseModel>::~MovingPhaseModel()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasePhaseModel>
@@ -255,8 +248,7 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
 
     if (K_.valid())
     {
-        K_.clear();
+        K_.ref() = 0.5*magSqr(this->U());
-        K();
     }
 }
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,7 +66,7 @@ class MovingPhaseModel
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Velocity field
         volVectorField U_;
@@ -103,7 +104,7 @@ protected:
 
 private:
 
-    // Private static member functions
+    // Private Member Functions
 
         //- Calculate and return the flux field
         tmp<surfaceScalarField> phi(const volVectorField& U) const;
@@ -113,6 +114,7 @@ public:
 
     // Constructors
 
+        //- Construct from phase system and phase name
         MovingPhaseModel
         (
             const phaseSystem& fluid,
@@ -122,7 +124,7 @@ public:
 
 
     //- Destructor
-    virtual ~MovingPhaseModel();
+    virtual ~MovingPhaseModel() = default;
 
 
     // Member Functions