multiphaseEulerFoam: Corrected handling of K for energy equation

Resolves bug-report https://bugs.openfoam.org/view.php?id=3634
OpenFOAM · Feb 24, 2021 · a060c63 · a060c63
1 parent 9c82a94
commit a060c63
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 8 deletions.
diff --git a/...ltiphaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/...ltiphaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,12 +102,20 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
 
-    const volVectorField U(this->U());
-    const surfaceScalarField alphaPhi(this->alphaPhi());
-    const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
+    const tmp<volVectorField> tU(this->U());
+    const volVectorField& U(tU());
 
-    const volScalarField contErr(this->continuityError());
-    const volScalarField K(this->K());
+    const tmp<surfaceScalarField> talphaPhi(this->alphaPhi());
+    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();
 
@@ -119,6 +127,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
       + fvc::ddt(alpha, this->rho(), K) + fvc::div(alphaRhoPhi, K)
       - contErr*K
+
       + this->divq(he)
      ==
         alpha*this->Qdot()

diff --git a/...ultiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/...ultiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -264,8 +264,7 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
 
     if (K_.valid())
     {
-        K_.clear();
-        K();
+        K_.ref() = 0.5*magSqr(this->U());
     }
 }