fvmLaplacian: Added laplacianCorrection functions and updated all thermal transport implementations

Now that thermal transport is implemented as an energy implicit correction on an
explicit temperature gradient formulation it is more efficient if the implicit
correction contains only the implicit terms of the matrix and not the explicit
non-orthogonal or anisotropic correction terms which are cancelled anyway when
the evaluation of the matrix for the current state is subtracted.  The new
fvm::laplacianCorrection functions provide a convenient mechanism to efficiently
evaluate only the implicit correction to the laplacian and is now used in all
the thermophysicalTransportModels.

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C

@@ -316,11 +316,21 @@ Foam::StationaryPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
 {
     const volScalarField& alpha = *this;
 
-    return -fvm::laplacian
+    const surfaceScalarField alphaKappa
     (
-        fvc::interpolate(alpha)*fvc::interpolate(this->thermo().alphahe()),
-        he
+        alpha.name() + '*' + this->thermo().kappa().name(),
+        fvc::interpolate(alpha)*fvc::interpolate(this->thermo().kappa())
     );
+
+    // Return heat flux source as an implicit energy correction
+    // to the temperature gradient flux
+    return
+       -fvc::laplacian(alphaKappa, this->thermo().T())
+       -fvm::laplacianCorrection
+        (
+            alphaKappa/fvc::interpolate(this->thermo().Cpv()),
+            he
+        );
 }
 
 

src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C

@@ -309,7 +309,7 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField sumJ
     (

src/ThermophysicalTransportModels/fluid/laminar/Fourier/Fourier.C

@@ -130,8 +130,12 @@ Fourier<BasicThermophysicalTransportModel>::divq(volScalarField& he) const
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-        -correction(fvm::laplacian(this->alpha()*thermo.alphahe(), he))
-        -fvc::laplacian(this->alpha()*thermo.kappa(), thermo.T());
+        -fvc::laplacian(this->alpha()*thermo.kappa(), thermo.T())
+        -fvm::laplacianCorrection
+         (
+             this->alpha()*thermo.kappa()/thermo.Cpv(),
+             he
+         );
 }
 
 

src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C

@@ -310,7 +310,7 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField sumJ
     (

src/ThermophysicalTransportModels/fluid/turbulence/eddyDiffusivity/eddyDiffusivity.C

@@ -163,8 +163,8 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-        -correction(fvm::laplacian(this->alpha()*this->alphaEff(), he))
-        -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T());
+        -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T())
+        -fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 }
 
 

src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C

@@ -162,7 +162,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField hGradY
     (

src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.C

@@ -401,19 +401,17 @@ Foam::solidThermophysicalTransportModels::anisotropic::divq
 {
     const solidThermo& thermo = this->thermo();
     const volSymmTensorField Kappa(this->Kappa());
+    const surfaceVectorField& Sf = thermo.mesh().Sf();
+    const surfaceScalarField& magSf = thermo.mesh().magSf();
 
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
        -fvc::laplacian(Kappa, thermo.T())
-       -correction
+       -fvm::laplacianCorrection
         (
-            fvm::laplacian
-            (
-                Kappa/thermo.Cv(),
-                e,
-                "laplacian(alphae,e)"
-            )
+            (Sf & fvc::interpolate(Kappa/thermo.Cv()) & Sf)/sqr(magSf),
+            e
         );
 }
 

src/ThermophysicalTransportModels/solid/isotropic/isotropic.C

@@ -105,15 +105,7 @@ Foam::solidThermophysicalTransportModels::isotropic::divq
     // to the temperature gradient flux
     return
        -fvc::laplacian(kappa(), thermo.T())
-       -correction
-        (
-            fvm::laplacian
-            (
-                kappa()/thermo.Cv(),
-                e,
-                "laplacian(alphae,e)"
-            )
-        );
+       -fvm::laplacianCorrection(kappa()/thermo.Cv(), e);
 }
 
 

src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C

@@ -27,6 +27,7 @@ License
 #include "surfaceFields.H"
 #include "fvMatrix.H"
 #include "laplacianScheme.H"
+#include "gaussLaplacianScheme.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -229,9 +230,9 @@ laplacian
     const word& name
 )
 {
-    tmp<fvMatrix<Type>> Laplacian(fvm::laplacian(tgamma(), vf, name));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tgamma(), vf, name));
     tgamma.clear();
-    return Laplacian;
+    return tLaplacian;
 }
 
 
@@ -260,9 +261,9 @@ laplacian
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    tmp<fvMatrix<Type>> Laplacian(fvm::laplacian(tgamma(), vf));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tgamma(), vf));
     tgamma.clear();
-    return Laplacian;
+    return tLaplacian;
 }
 
 
@@ -325,9 +326,79 @@ laplacian
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    tmp<fvMatrix<Type>> tfvm(fvm::laplacian(tGamma(), vf));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tGamma(), vf));
     tGamma.clear();
-    return tfvm;
+    return tLaplacian;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return fvm::laplacianCorrection(fvc::interpolate(gamma), vf);
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const tmp<GeometricField<scalar, fvPatchField, volMesh>>& tgamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tLaplacianCorrection
+    (
+        fvm::laplacianCorrection(tgamma(), vf)
+    );
+    tgamma.clear();
+    return tLaplacianCorrection;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return correction
+    (
+        fv::gaussLaplacianScheme<Type, scalar>::fvmLaplacianUncorrected
+        (
+            gamma*vf.mesh().magSf(),
+            vf.mesh().nonOrthDeltaCoeffs(),
+            vf
+        )
+    );
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const tmp<GeometricField<scalar, fvsPatchField, surfaceMesh>>& tGamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tLaplacianCorrection
+    (
+        fvm::laplacianCorrection(tGamma(), vf)
+    );
+    tGamma.clear();
+    return tLaplacianCorrection;
 }
 
 

src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H

@@ -176,6 +176,36 @@ namespace fvm
         const tmp<GeometricField<GType, fvsPatchField, surfaceMesh>>&,
         const GeometricField<Type, fvPatchField, volMesh>&
     );
+
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const GeometricField<scalar, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const tmp<GeometricField<scalar, fvPatchField, volMesh>>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const GeometricField<scalar, fvsPatchField, surfaceMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const tmp<GeometricField<scalar, fvsPatchField, surfaceMesh>>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
 }