BUG: heSolidThermo: initialise mu,psi. See #905.

src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,6 +36,8 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
     const scalarField& hCells = this->he_;
     const scalarField& pCells = this->p_;
     scalarField& rhoCells = this->rho_.primitiveFieldRef();
+    //scalarField& psiCells = this->psi_.primitiveFieldRef();
+    //scalarField& muCells = this->mu_.primitiveFieldRef();
     scalarField& alphaCells = this->alpha_.primitiveFieldRef();
 
     forAll(TCells, celli)
@@ -57,6 +59,8 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
         }
 
         rhoCells[celli] = volMixture_.rho(pCells[celli], TCells[celli]);
+        //psiCells[celli] = volMixture_.psi(pCells[celli], TCells[celli]);
+        //muCells[celli] = volMixture_.mu(pCells[celli], TCells[celli]);
 
         alphaCells[celli] =
             volMixture_.kappa(pCells[celli], TCells[celli])
@@ -64,26 +68,21 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
             mixture_.Cpv(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& pBf =
-        this->p_.boundaryFieldRef();
-
-    volScalarField::Boundary& TBf =
-        this->T_.boundaryFieldRef();
-
-    volScalarField::Boundary& rhoBf =
-        this->rho_.boundaryFieldRef();
-
-    volScalarField::Boundary& heBf =
-        this->he().boundaryFieldRef();
-
-    volScalarField::Boundary& alphaBf =
-        this->alpha_.boundaryFieldRef();
+    volScalarField::Boundary& pBf = this->p_.boundaryFieldRef();
+    volScalarField::Boundary& TBf = this->T_.boundaryFieldRef();
+    volScalarField::Boundary& rhoBf = this->rho_.boundaryFieldRef();
+    //volScalarField::Boundary& psiBf = this->psi_.boundaryFieldRef();
+    //volScalarField::Boundary& muBf = this->mu_.boundaryFieldRef();
+    volScalarField::Boundary& heBf = this->he().boundaryFieldRef();
+    volScalarField::Boundary& alphaBf = this->alpha_.boundaryFieldRef();
 
     forAll(this->T_.boundaryField(), patchi)
     {
         fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& prho = rhoBf[patchi];
+        //fvPatchScalarField& ppsi = psiBf[patchi];
+        //fvPatchScalarField& pmu = muBf[patchi];
         fvPatchScalarField& phe = heBf[patchi];
         fvPatchScalarField& palpha = alphaBf[patchi];
 
@@ -106,6 +105,8 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
 
                 phe[facei] = mixture_.HE(pp[facei], pT[facei]);
                 prho[facei] = volMixture_.rho(pp[facei], pT[facei]);
+                //ppsi[facei] = volMixture_.psi(pp[facei], pT[facei]);
+                //pmu[facei] = volMixture_.mu(pp[facei], pT[facei]);
 
                 palpha[facei] =
                     volMixture_.kappa(pp[facei], pT[facei])
@@ -134,6 +135,8 @@ void Foam::heSolidThermo<BasicSolidThermo, MixtureType>::calculate()
                 }
 
                 prho[facei] = volMixture_.rho(pp[facei], pT[facei]);
+                //ppsi[facei] = volMixture_.psi(pp[facei], pT[facei]);
+                //pmu[facei] = volMixture_.mu(pp[facei], pT[facei]);
 
                 palpha[facei] =
                     volMixture_.kappa(pp[facei], pT[facei])
@@ -188,8 +191,14 @@ heSolidThermo
     heThermo<BasicSolidThermo, MixtureType>(mesh, phaseName, dictName)
 {
     calculate();
+
+    // TBD. initialise psi, mu (at heThermo level) since these do not
+    // get initialised. Move to heThermo constructor?
+    this->mu_ == dimensionedScalar("zero", this->mu_.dimensions(), 0.0);
+    this->psi_ == dimensionedScalar("zero", this->psi_.dimensions(), 0.0);
 }
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class BasicSolidThermo, class MixtureType>

[Mark OLESEN]

Just as a point of interest, for zero can also use the more compact two parameter form (ref. 5632ef2d). Eg,

this->mu_ == dimensionedScalar(this->mu_.dimensions(), Zero);

also works for other dimensioned types

dimensionedVector(dimVelocity, Zero) ...

I think this is especially useful in combination with creating a volume field from components.