From 414f2d1d489eec6b6d4bd0f4ef392587ec426fa6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 13 Nov 2014 14:32:44 +0000
Subject: [PATCH] BUG: turbulentTemperatureRadCoupledMixed: incorrect
 multi-layer averaging (#1440)

---
 ...mperatureRadCoupledMixedFvPatchScalarField.C | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
index ac6d9c66eb..36c6106c73 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C
@@ -123,14 +123,16 @@ turbulentTemperatureRadCoupledMixedFvPatchScalarField
 
         if (thicknessLayers_.size() > 0)
         {
+            // total thermal transmittance by harmonic averaging
             forAll (thicknessLayers_, iLayer)
             {
                 const scalar l = thicknessLayers_[iLayer];
                 if (l > 0.0)
                 {
-                    contactRes_ += kappaLayers_[iLayer]/l;
+                    contactRes_ += l/kappaLayers_[iLayer]; // inverse sum
                 }
             }
+            contactRes_ = 1.0/contactRes_; // new total inverse
         }
     }
 
@@ -210,17 +212,16 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
 
 
     // Swap to obtain full local values of neighbour K*delta
-    tmp<scalarField> KDeltaNbr(new scalarField(TcNbr.size(), 0.0));
+    scalarField KDeltaNbr;
     if (contactRes_ == 0.0)
     {
-        // Swap to obtain full local values of neighbour K*delta
-        KDeltaNbr() = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();
+        KDeltaNbr = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();
     }
     else
     {
-        KDeltaNbr() = contactRes_;
+        KDeltaNbr.setSize(nbrField.size(), contactRes_);
     }
-    mpp.distribute(KDeltaNbr());
+    mpp.distribute(KDeltaNbr);
 
     scalarField KDelta(kappa(*this)*patch().deltaCoeffs());
 
@@ -237,11 +238,11 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
         mpp.distribute(QrNbr);
     }
 
-    scalarField alpha(KDeltaNbr() - (Qr + QrNbr)/Tp);
+    scalarField alpha(KDeltaNbr - (Qr + QrNbr)/Tp);
 
     valueFraction() = alpha/(alpha + KDelta);
 
-    refValue() = (KDeltaNbr()*TcNbr)/alpha;
+    refValue() = (KDeltaNbr*TcNbr)/alpha;
 
     mixedFvPatchScalarField::updateCoeffs();