chtMultiRegionFoam, heSolidThermo: Moved the solid heat flux model in…

…to heSolidThermo

and changed to be an energy implicit correction to a temperature gradient
based heat-flux.  This formulation is both energy conservative and temperature
consistent.

The wallHeatFlux functionObject has been updated to use a consistent heat-flux
from the heSolidThermo.

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -1,10 +1,7 @@
 // Initialise solid field pointer lists
-PtrList<coordinateSystem> coordinates(solidRegions.size());
 PtrList<solidThermo> thermos(solidRegions.size());
 PtrList<radiationModel> radiations(solidRegions.size());
 PtrList<fv::options> solidHeatSources(solidRegions.size());
-PtrList<volScalarField> betavSolid(solidRegions.size());
-PtrList<volSymmTensorField> aniAlphas(solidRegions.size());
 
 // Populate solid field pointer lists
 forAll(solidRegions, i)
@@ -24,84 +21,4 @@ forAll(solidRegions, i)
         i,
         new fv::options(solidRegions[i])
     );
-
-    if (!thermos[i].isotropic())
-    {
-        Info<< "    Adding coordinateSystems\n" << endl;
-        coordinates.set
-        (
-            i,
-            coordinateSystem::New(solidRegions[i], thermos[i].properties())
-        );
-
-        tmp<volVectorField> tkappaByCv =
-            thermos[i].Kappa()/thermos[i].Cv();
-
-        aniAlphas.set
-        (
-            i,
-            new volSymmTensorField
-            (
-                IOobject
-                (
-                    "Anialpha",
-                    runTime.timeName(),
-                    solidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                solidRegions[i],
-                dimensionedSymmTensor
-                (
-                    "zero",
-                    tkappaByCv().dimensions(),
-                    Zero
-                ),
-                zeroGradientFvPatchSymmTensorField::typeName
-            )
-        );
-
-        aniAlphas[i].primitiveFieldRef() =
-            coordinates[i].R().transformVector(tkappaByCv());
-        aniAlphas[i].correctBoundaryConditions();
-
-    }
-
-    IOobject betavSolidIO
-    (
-        "betavSolid",
-        runTime.timeName(),
-        solidRegions[i],
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    );
-
-    if (betavSolidIO.typeHeaderOk<volScalarField>(true))
-    {
-        betavSolid.set
-        (
-            i,
-            new volScalarField(betavSolidIO, solidRegions[i])
-        );
-    }
-    else
-    {
-        betavSolid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "betavSolid",
-                    runTime.timeName(),
-                    solidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                solidRegions[i],
-                dimensionedScalar(dimless, scalar(1))
-            )
-        );
-    }
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -3,31 +3,8 @@ solidThermo& thermo = thermos[i];
 tmp<volScalarField> trho = thermo.rho();
 const volScalarField& rho = trho();
 
-tmp<volScalarField> tCv = thermo.Cv();
-const volScalarField& Cv = tCv();
-
-tmp<volSymmTensorField> taniAlpha;
-if (!thermo.isotropic())
-{
-    volSymmTensorField& aniAlpha = aniAlphas[i];
-    tmp<volVectorField> tkappaByCv = thermo.Kappa()/Cv;
-    const coordinateSystem& coodSys = coordinates[i];
-
-    aniAlpha.primitiveFieldRef() =
-        coodSys.R().transformVector(tkappaByCv());
-    aniAlpha.correctBoundaryConditions();
-
-    taniAlpha = tmp<volSymmTensorField>
-    (
-        new volSymmTensorField(aniAlpha)
-    );
-}
-
-
 volScalarField& e = thermo.he();
 
-const volScalarField& betav = betavSolid[i];
-
 const fv::options& fvOptions = solidHeatSources[i];
 
 solidNoLoopControl& pimple = pimples.solid(i);

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -3,12 +3,8 @@
     {
         fvScalarMatrix eEqn
         (
-            fvm::ddt(betav*rho, e)
-          - (
-                thermo.isotropic()
-              ? fvm::laplacian(betav*thermo.alpha(), e, "laplacian(alpha,e)")
-              : fvm::laplacian(betav*taniAlpha(), e, "laplacian(alpha,e)")
-            )
+            fvm::ddt(rho, e)
+          + thermo.divq(e)
           ==
             fvOptions(rho, e)
         );

src/functionObjects/field/wallHeatFlux/wallHeatFlux.C

@@ -104,6 +104,52 @@ Foam::functionObjects::wallHeatFlux::calcWallHeatFlux(const volVectorField& q)
 }
 
 
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::wallHeatFlux::calcWallHeatFlux
+(
+    const surfaceScalarField& q
+)
+{
+    tmp<volScalarField> twallHeatFlux
+    (
+        volScalarField::New
+        (
+            type(),
+            mesh_,
+            dimensionedScalar(dimMass/pow3(dimTime), 0)
+        )
+    );
+
+    volScalarField::Boundary& wallHeatFluxBf =
+        twallHeatFlux.ref().boundaryFieldRef();
+
+    const surfaceScalarField::Boundary& qBf = q.boundaryField();
+
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        const label patchi = iter.key();
+
+        wallHeatFluxBf[patchi] = -qBf[patchi];
+    }
+
+    if (foundObject<volScalarField>("qr"))
+    {
+        const volScalarField& qr = lookupObject<volScalarField>("qr");
+
+        const volScalarField::Boundary& radHeatFluxBf = qr.boundaryField();
+
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            const label patchi = iter.key();
+
+            wallHeatFluxBf[patchi] -= radHeatFluxBf[patchi];
+        }
+    }
+
+    return twallHeatFlux;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjects::wallHeatFlux::wallHeatFlux
@@ -206,22 +252,14 @@ bool Foam::functionObjects::wallHeatFlux::execute()
                 thermophysicalTransportModel::typeName
             );
 
-        return store
-        (
-            name,
-            calcWallHeatFlux(ttm.q())
-        );
+        return store(name, calcWallHeatFlux(ttm.q()));
     }
     else if (foundObject<solidThermo>(solidThermo::dictName))
     {
         const solidThermo& thermo =
             lookupObject<solidThermo>(solidThermo::dictName);
 
-        return store
-        (
-            name,
-            calcWallHeatFlux(-thermo.alpha()*fvc::grad(thermo.he()))
-        );
+        return store(name, calcWallHeatFlux(thermo.q()));
     }
     else
     {

src/functionObjects/field/wallHeatFlux/wallHeatFlux.H

@@ -75,6 +75,7 @@ SourceFiles
 #include "writeLocalObjects.H"
 #include "HashSet.H"
 #include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -110,6 +111,9 @@ protected:
         //- Calculate the heat-flux
         tmp<volScalarField> calcWallHeatFlux(const volVectorField& q);
 
+        //- Calculate the heat-flux
+        tmp<volScalarField> calcWallHeatFlux(const surfaceScalarField& q);
+
 
 public:
 

src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C

@@ -24,7 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "heSolidThermo.H"
-#include "volFields.H"
+#include "fvmLaplacian.H"
+#include "fvcLaplacian.H"
+#include "coordinateSystem.H"
+#include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -267,4 +270,81 @@ Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa
 }
 
 
+template<class BasicSolidThermo, class MixtureType>
+Foam::tmp<Foam::volSymmTensorField>
+Foam::heSolidThermo<BasicSolidThermo, MixtureType>::KappaLocal() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volSymmTensorField> tKappaLocal;
+
+    if (!isotropic())
+    {
+        coordinateSystem coordinates
+        (
+            coordinateSystem::New(mesh, this->properties())
+        );
+
+        const tmp<volVectorField> tKappa(Kappa());
+
+        tKappaLocal =
+        (
+            volSymmTensorField::New
+            (
+                "KappaLocal",
+                mesh,
+                dimensionedSymmTensor(tKappa().dimensions(), Zero),
+                zeroGradientFvPatchSymmTensorField::typeName
+            )
+        );
+        volSymmTensorField& KappaLocal = tKappaLocal.ref();
+
+        KappaLocal.primitiveFieldRef() =
+            coordinates.R().transformVector(tKappa());
+        KappaLocal.correctBoundaryConditions();
+    }
+
+    return tKappaLocal;
+}
+
+
+template<class BasicSolidThermo, class MixtureType>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::heSolidThermo<BasicSolidThermo, MixtureType>::q() const
+{
+    return
+      - (
+            isotropic()
+          ? fvm::laplacian(this->kappa(), this->T_)().flux()
+          : fvm::laplacian(KappaLocal(), this->T_)().flux()
+        );
+}
+
+
+template<class BasicSolidThermo, class MixtureType>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::heSolidThermo<BasicSolidThermo, MixtureType>::divq
+(
+    volScalarField& e
+) const
+{
+    return
+      - (
+            isotropic()
+          ?   fvc::laplacian(this->kappa(), this->T_)
+            + correction(fvm::laplacian(this->alpha(), e))
+          :   fvc::laplacian(KappaLocal(), this->T_)
+            + correction
+              (
+                  fvm::laplacian
+                  (
+                      KappaLocal()/this->Cv(),
+                      e,
+                      "laplacian(" + this->alpha().name() + ",e)"
+                  )
+              )
+        );
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.H

@@ -56,6 +56,9 @@ class heSolidThermo
         //- Calculate the thermo variables
         void calculate();
 
+        //- Return the anisotropic Kappa in the local coordinate system
+        tmp<volSymmTensorField> KappaLocal() const;
+
 
 public:
 
@@ -96,23 +99,23 @@ public:
         //- Update properties
         virtual void correct();
 
+        //- Return true if thermal conductivity is isotropic
+        virtual bool isotropic() const
+        {
+            return MixtureType::thermoType::isotropic;
+        }
 
-        // Derived thermal properties
-
-            //- Anisotropic thermal conductivity [W/m/K]
-            virtual tmp<volVectorField> Kappa() const;
-
-            //- Return true if thermal conductivity is isotropic
-            virtual bool isotropic() const
-            {
-                return MixtureType::thermoType::isotropic;
-            }
+        //- Anisotropic thermal conductivity [W/m/K]
+        virtual tmp<volVectorField> Kappa() const;
 
+        //- Anisotropic thermal conductivity [W/m/K]
+        virtual tmp<vectorField> Kappa(const label patchi) const;
 
-        // Per patch calculation
+        //- Return the heat flux
+        virtual tmp<surfaceScalarField> q() const;
 
-            //- Anisotropic thermal conductivity [W/m/K]
-            virtual tmp<vectorField> Kappa(const label patchi) const;
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 
 
     // Member Operators