diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
index 24316a1685..e916b7a56a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H
@@ -1,30 +1,30 @@
+if (mesh.changing())
 {
-    if (mesh.changing())
+    forAll(U.boundaryField(), patchi)
     {
-        forAll(U.boundaryField(), patchi)
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].initEvaluate();
-            }
-        }
-
-        forAll(U.boundaryField(), patchi)
-        {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].evaluate();
-
-                phi.boundaryField()[patchi] =
-                    rho.boundaryField()[patchi]
-                   *(
-                       U.boundaryField()[patchi]
-                     & mesh.Sf().boundaryField()[patchi]
-                    );
-            }
+            U.boundaryField()[patchi].initEvaluate();
         }
     }
 
+    forAll(U.boundaryField(), patchi)
+    {
+        if (U.boundaryField()[patchi].fixesValue())
+        {
+            U.boundaryField()[patchi].evaluate();
+
+            phi.boundaryField()[patchi] =
+                rho.boundaryField()[patchi]
+               *(
+                   U.boundaryField()[patchi]
+                 & mesh.Sf().boundaryField()[patchi]
+               );
+        }
+    }
+}
+
+{
     volScalarField pcorr
     (
         IOobject
@@ -40,13 +40,13 @@
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divrhoU
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divrhoU
         );
 
         pcorrEqn.solve();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 1d70a0f354..40cf109457 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -22,7 +22,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )/fvc::interpolate(rho)
     );
 
@@ -55,7 +55,7 @@ else
     (
         "phiHbyA",
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
+      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index dac6980917..aae3285186 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
@@ -56,13 +56,10 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "createPcorrTypes.H"
+    #include "createRhoUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    // Create old-time absolute flux for ddtCorr
-    surfaceScalarField phiAbs("phiAbs", phi);
-
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -72,12 +69,6 @@ int main(int argc, char *argv[])
         #include "readControls.H"
         #include "compressibleCourantNo.H"
 
-        // Make the fluxes absolute before mesh-motion
-        fvc::makeAbsolute(phi, rho, U);
-
-        // Update absolute flux for ddtCorr
-        phiAbs = phi;
-
         #include "setDeltaT.H"
 
         runTime++;
@@ -86,20 +77,23 @@ int main(int argc, char *argv[])
 
         {
             // Store divrhoU from the previous time-step/mesh for the correctPhi
-            volScalarField divrhoU(fvc::div(phi));
+            volScalarField divrhoU(fvc::div(fvc::absolute(phi, rho, U)));
 
             // Do any mesh changes
             mesh.update();
 
             if (mesh.changing() && correctPhi)
             {
+                // Calculate absolute flux from the mapped surface velocity
+                phi = mesh.Sf() & rhoUf;
+
                 #include "correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, rho, U);
             }
         }
 
-        // Make the fluxes relative to the mesh-motion
-        fvc::makeRelative(phi, rho, U);
-
         if (mesh.changing() && checkMeshCourantNo)
         {
             #include "meshCourantNo.H"
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
index ebc62e7e59..b693bb7860 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H
@@ -37,13 +37,13 @@ if (mesh.changing())
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi)
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -54,6 +54,6 @@ if (mesh.changing())
             phi -= pcorrEqn.flux();
         }
     }
-}
 
-#include "continuityErrs.H"
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
index 2dca530c14..c3fde6b3a1 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/cavitatingDyMFoam.C
@@ -47,16 +47,15 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readThermodynamicProperties.H"
-    #include "readControls.H"
-    #include "createFields.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
 
-    surfaceScalarField phivAbs("phivAbs", phiv);
-    fvc::makeAbsolute(phivAbs, U);
-
+    #include "readThermodynamicProperties.H"
+    #include "readControls.H"
+    #include "createFields.H"
+    #include "createUf.H"
+    #include "createPcorrTypes.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -75,18 +74,8 @@ int main(int argc, char *argv[])
 
         scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        {
-            // Calculate the relative velocity used to map relative flux phiv
-            volVectorField Urel("Urel", U);
-
-            if (mesh.moving())
-            {
-                Urel -= fvc::reconstruct(fvc::meshPhi(U));
-            }
-
-            // Do any mesh changes
-            mesh.update();
-        }
+        // Do any mesh changes
+        mesh.update();
 
         if (mesh.changing())
         {
@@ -94,7 +83,16 @@ int main(int argc, char *argv[])
                 << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                 << " s" << endl;
 
-            #include "correctPhi.H"
+            if (correctPhi)
+            {
+                // Calculate absolute flux from the mapped surface velocity
+                phiv = mesh.Sf() & Uf;
+
+                #include "correctPhi.H"
+
+                // Make the flux relative to the mesh motion
+                fvc::makeRelative(phiv, U);
+            }
         }
 
         // --- Pressure-velocity PIMPLE corrector loop
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
index 8115b791e9..c05d664838 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/correctPhi.H
@@ -1,18 +1,27 @@
+if (mesh.changing())
 {
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll (p.boundaryField(), i)
+    forAll(U.boundaryField(), patchI)
     {
-        if (p.boundaryField()[i].fixesValue())
+        if (U.boundaryField()[patchI].fixesValue())
         {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+            U.boundaryField()[patchI].initEvaluate();
         }
     }
 
+    forAll(U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].fixesValue())
+        {
+            U.boundaryField()[patchI].evaluate();
+
+            phiv.boundaryField()[patchI] =
+                U.boundaryField()[patchI]
+              & mesh.Sf().boundaryField()[patchI];
+        }
+    }
+}
+
+{
     volScalarField pcorr
     (
         IOobject
@@ -29,7 +38,7 @@
     );
 
     surfaceScalarField rhof(fvc::interpolate(rho, "div(phiv,rho)"));
-    dimensionedScalar rAUf("(1|A(U))", dimTime, 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index e7f078d1a6..0f0d576841 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -18,7 +18,7 @@
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + rhorAUf*fvc::ddtCorr(U, phivAbs);
+         + rhorAUf*fvc::ddtCorr(U, Uf);
     fvc::makeRelative(phiv, U);
 
     surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
@@ -43,7 +43,6 @@
         if (pimple.finalNonOrthogonalIter())
         {
             phiv += (phiGradp + pEqn.flux())/rhof;
-            phivAbs = fvc::absolute(phiv, U);
         }
     }
 
@@ -82,4 +81,10 @@
     U.correctBoundaryConditions();
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
+
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += mesh.Sf()*(phiv - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    }
 }
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H
new file mode 100644
index 0000000000..216e9b4247
--- /dev/null
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/readControls.H
@@ -0,0 +1,9 @@
+#include "readTimeControls.H"
+
+scalar maxAcousticCo
+(
+    readScalar(runTime.controlDict().lookup("maxAcousticCo"))
+);
+
+bool correctPhi =
+    pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
index 76e9217b65..97c63e6a3e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H
@@ -1,28 +1,26 @@
+if (mesh.changing())
 {
-    if (mesh.changing())
+    forAll(U.boundaryField(), patchi)
     {
-        forAll(U.boundaryField(), patchi)
+        if (U.boundaryField()[patchi].fixesValue())
         {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].initEvaluate();
-            }
-        }
-
-        forAll(U.boundaryField(), patchi)
-        {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].evaluate();
-
-                phi.boundaryField()[patchi] =
-                U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
-            }
+            U.boundaryField()[patchi].initEvaluate();
         }
     }
 
-    #include "continuityErrs.H"
+    forAll(U.boundaryField(), patchi)
+    {
+        if (U.boundaryField()[patchi].fixesValue())
+        {
+            U.boundaryField()[patchi].evaluate();
 
+            phi.boundaryField()[patchi] =
+            U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
+        }
+    }
+}
+
+{
     volScalarField pcorr
     (
         IOobject
@@ -38,13 +36,13 @@
         pcorrTypes
     );
 
-    dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divU
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
         );
 
         pcorrEqn.solve();
@@ -54,8 +52,8 @@
             phi -= pcorrEqn.flux();
         }
     }
+
+    phi.oldTime() = phi;
+
+    #include "continuityErrs.H"
 }
-
-phi.oldTime() = phi;
-
-#include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
index 57a79b6091..e291b7dbcd 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
@@ -1,20 +1,6 @@
 {
     #include "continuityErrs.H"
 
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll (p.boundaryField(), i)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
     volScalarField pcorr
     (
         IOobject
@@ -30,9 +16,7 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    //adjustPhi(phi, U, pcorr);
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index 0c5c8a2a86..d0285dc357 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
     #include "createMRFZones.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
+    #include "createPcorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
index 7120a09919..2df4bf7ffb 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
@@ -55,15 +55,15 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 {
     patchType() = ptf.patchType();
 
-    // Map value. Set unmapped values and overwrite with mapped ptf
-    if (&iF && iF.size())
-    {
-        fvPatchField<scalar>::operator=(patchInternalField());
-    }
-    map(ptf, mapper);
     // Map gradient. Set unmapped values and overwrite with mapped ptf
     gradient() = 0.0;
     gradient().map(ptf.gradient(), mapper);
+
+    // Evaluate the value field from the gradient if the internal field is valid
+    if (&iF && iF.size())
+    {
+        evaluate();
+    }
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
index f3bf760877..5c53f64fdd 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -174,7 +174,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvc::absolute
 {
     if (tphi().mesh().moving())
     {
-        return tphi + fvc::meshPhi(rho, U);
+        return tphi + fvc::interpolate(rho)*fvc::meshPhi(rho, U);
     }
     else
     {
diff --git a/src/regionModels/regionCoupling/Make/files b/src/regionModels/regionCoupling/Make/files
index 846df33e9e..c1ebaf7990 100644
--- a/src/regionModels/regionCoupling/Make/files
+++ b/src/regionModels/regionCoupling/Make/files
@@ -1,4 +1,5 @@
 derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
 derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C
+derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libregionCoupling
diff --git a/src/regionModels/regionCoupling/Make/options b/src/regionModels/regionCoupling/Make/options
index 28569ccda6..741293c78e 100644
--- a/src/regionModels/regionCoupling/Make/options
+++ b/src/regionModels/regionCoupling/Make/options
@@ -6,8 +6,15 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude\
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
@@ -23,6 +30,8 @@ LIB_LIBS = \
     -lpyrolysisModels \
     -lsurfaceFilmModels \
     -lsolidChemistryModel \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
     -lfiniteVolume \
     -lmeshTools \
     -lcompressibleRASModels \
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
new file mode 100644
index 0000000000..ed3d509d60
--- /dev/null
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
@@ -0,0 +1,372 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mappedPatchBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::filmModelType&
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+filmModel() const
+{
+    const regionModels::regionModel& model =
+        db().time().lookupObject<regionModels::regionModel>
+        (
+            "surfaceFilmProperties"
+        );
+
+    return dynamic_cast<const filmModelType&>(model);
+}
+
+
+const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+pyrolysisModelType&
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+pyrModel() const
+{
+    const regionModels::regionModel& model =
+        db().time().lookupObject<regionModels::regionModel>
+        (
+            "pyrolysisProperties"
+        );
+
+    return dynamic_cast<const pyrolysisModelType&>(model);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    temperatureCoupledBase(patch(), "undefined", "undefined-K"),
+    TnbrName_("undefined-Tnbr"),
+    QrNbrName_("undefined-QrNbr"),
+    QrName_("undefined-Qr"),
+    convectiveScaling_(1.0),
+    filmDeltaDry_(0.0),
+    filmDeltaWet_(0.0)
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 1.0;
+}
+
+
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+(
+    const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField& psf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(psf, p, iF, mapper),
+    temperatureCoupledBase(patch(), psf.KMethod(), psf.kappaName()),
+    TnbrName_(psf.TnbrName_),
+    QrNbrName_(psf.QrNbrName_),
+    QrName_(psf.QrName_),
+    convectiveScaling_(psf.convectiveScaling_),
+    filmDeltaDry_(psf.filmDeltaDry_),
+    filmDeltaWet_(psf.filmDeltaWet_)
+{}
+
+
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    temperatureCoupledBase(patch(), dict),
+    TnbrName_(dict.lookup("Tnbr")),
+    QrNbrName_(dict.lookup("QrNbr")),
+    QrName_(dict.lookup("Qr")),
+    convectiveScaling_(dict.lookupOrDefault<scalar>("convectiveScaling", 1.0)),
+    filmDeltaDry_(readScalar(dict.lookup("filmDeltaDry"))),
+    filmDeltaWet_(readScalar(dict.lookup("filmDeltaWet")))
+{
+    if (!isA<mappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::"
+            "filmPyrolysisRadiativeCoupledMixedFvPatchScalarField\n"
+            "(\n"
+            "    const fvPatch& p,\n"
+            "    const DimensionedField<scalar, volMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << mappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << dimensionedInternalField().name()
+            << " in file " << dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+
+    if (dict.found("refValue"))
+    {
+        // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+    }
+}
+
+
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+(
+    const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField& psf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(psf, iF),
+    temperatureCoupledBase(patch(), psf.KMethod(), psf.kappaName()),
+    TnbrName_(psf.TnbrName_),
+    QrNbrName_(psf.QrNbrName_),
+    QrName_(psf.QrName_),
+    convectiveScaling_(psf.convectiveScaling_),
+    filmDeltaDry_(psf.filmDeltaDry_),
+    filmDeltaWet_(psf.filmDeltaWet_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::
+updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Get the coupling information from the mappedPatchBase
+    const mappedPatchBase& mpp =
+        refCast<const mappedPatchBase>(patch().patch());
+
+    const label patchI = patch().index();
+    const label nbrPatchI = mpp.samplePolyPatch().index();
+
+    const polyMesh& mesh = patch().boundaryMesh().mesh();
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+    const fvPatch& nbrPatch =
+        refCast<const fvMesh>(nbrMesh).boundary()[nbrPatchI];
+
+    scalarField intFld(patchInternalField());
+
+    const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField&
+        nbrField =
+        refCast
+        <
+            const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+        >
+        (
+            nbrPatch.lookupPatchField<volScalarField, scalar>(TnbrName_)
+        );
+
+    // Swap to obtain full local values of neighbour internal field
+    scalarField nbrIntFld(nbrField.patchInternalField());
+    mpp.distribute(nbrIntFld);
+
+    scalarField& Tp = *this;
+
+    const scalarField K(this->kappa(*this));
+    const scalarField nbrK(nbrField.kappa(*this));
+
+    // Swap to obtain full local values of neighbour K*delta
+    scalarField KDeltaNbr(nbrK*nbrPatch.deltaCoeffs());
+    mpp.distribute(KDeltaNbr);
+
+    scalarField myKDelta(K*patch().deltaCoeffs());
+
+    scalarList Tfilm(patch().size(), 0.0);
+    scalarList htcwfilm(patch().size(), 0.0);
+    scalarList filmDelta(patch().size(), 0.0);
+
+    const pyrolysisModelType& pyrolysis = pyrModel();
+    const filmModelType& film = filmModel();
+
+    label myPatchINrbPatchI = -1;
+
+    // Obtain Rad heat (Qr)
+    scalarField Qr(patch().size(), 0.0);
+    if (QrName_ != "none") //region0
+    {
+        Qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
+        myPatchINrbPatchI = nbrPatch.index();
+    }
+
+    if (QrNbrName_ != "none") //pyrolysis
+    {
+        Qr = nbrPatch.lookupPatchField<volScalarField, scalar>(QrNbrName_);
+        mpp.distribute(Qr);
+        myPatchINrbPatchI = patchI;
+    }
+
+    const label filmPatchI =
+        pyrolysis.nbrCoupledPatchID(film, myPatchINrbPatchI);
+
+    const scalarField htcw(film.htcw().h()().boundaryField()[filmPatchI]);
+
+    // Obtain htcw
+    htcwfilm =
+        const_cast<pyrolysisModelType&>(pyrolysis).mapRegionPatchField
+        (
+            film,
+            myPatchINrbPatchI,
+            filmPatchI,
+            htcw,
+            true
+        );
+
+
+    // Obtain Tfilm at the boundary through Ts.
+    // NOTE: Tf is not good as at the boundary it will retrieve Tp
+    Tfilm = film.Ts().boundaryField()[filmPatchI];
+    film.toPrimary(filmPatchI, Tfilm);
+
+    // Obtain delta
+    filmDelta =
+        const_cast<pyrolysisModelType&>(pyrolysis).mapRegionPatchField<scalar>
+        (
+            film,
+            "deltaf",
+            myPatchINrbPatchI,
+            true
+        );
+
+     // Estimate wetness of the film (1: wet , 0: dry)
+     scalarField ratio
+     (
+        min
+        (
+            max
+            (
+                (filmDelta - filmDeltaDry_)/(filmDeltaWet_ - filmDeltaDry_),
+                scalar(0.0)
+            ),
+            scalar(1.0)
+        )
+     );
+
+    scalarField qConv(ratio*htcwfilm*(Tfilm - Tp)*convectiveScaling_);
+
+    scalarField qRad((1.0 - ratio)*Qr);
+
+    scalarField alpha(KDeltaNbr - (qRad + qConv)/Tp);
+
+    valueFraction() = alpha/(alpha + (1.0 - ratio)*myKDelta);
+
+    refValue() = ratio*Tfilm + (1.0 - ratio)*(KDeltaNbr*nbrIntFld)/alpha;
+
+    mixedFvPatchScalarField::updateCoeffs();
+
+    if (debug)
+    {
+        scalar Qc = gSum(qConv*patch().magSf());
+        scalar Qr = gSum(qRad*patch().magSf());
+        scalar Qt = gSum((qConv + qRad)*patch().magSf());
+
+        Info<< mesh.name() << ':'
+            << patch().name() << ':'
+            << this->dimensionedInternalField().name() << " <- "
+            << nbrMesh.name() << ':'
+            << nbrPatch.name() << ':'
+            << this->dimensionedInternalField().name() << " :" << nl
+            << "     convective heat[W] : " << Qc << nl
+            << "     radiative heat [W] : " << Qr << nl
+            << "     total heat     [W] : " << Qt << nl
+            << "     walltemperature "
+            << " min:" << gMin(*this)
+            << " max:" << gMax(*this)
+            << " avg:" << gAverage(*this)
+            << endl;
+    }
+
+}
+
+
+void filmPyrolysisRadiativeCoupledMixedFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    mixedFvPatchScalarField::write(os);
+    os.writeKeyword("Tnbr")<< TnbrName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("QrNbr")<< QrNbrName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("convectiveScaling") << convectiveScaling_
+    << token::END_STATEMENT << nl;
+    os.writeKeyword("filmDeltaDry") << filmDeltaDry_ <<
+    token::END_STATEMENT << nl;
+    os.writeKeyword("filmDeltaWet") << filmDeltaWet_ <<
+    token::END_STATEMENT << endl;
+    temperatureCoupledBase::write(os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.H b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.H
new file mode 100644
index 0000000000..734b19b364
--- /dev/null
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisRadiativeCoupledMixed/filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.H
@@ -0,0 +1,232 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+
+    Foam::
+    compressible::
+    filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+
+Description
+    Mixed boundary condition for temperature, to be used in the flow and
+    pyrolysis regions when a film region model is used.
+
+
+    Example usage:
+        myInterfacePatchName
+        {
+            type            filmPyrolysisRadiativeCoupledMixed;
+            Tnbr            T;
+            kappa           fluidThermo;
+            QrNbr           none;
+            Qr              Qr;
+            kappaName       none;
+            filmDeltaDry    0.0;
+            filmDeltaWet    3e-4;
+            value           $internalField;
+        }
+
+    Needs to be on underlying mapped(Wall)FvPatch.
+    It calculates local field as
+         ratio = (filmDelta - filmDeltaDry)/(filmDeltaWet - filmDeltaDry),
+
+    when ratio = 1 is considered wet and the film temperarture is fixed at
+    the wall. If ratio = 0 (dry) it emulates the normal radiative solid BC.
+
+    In between ratio 0 and 1 the gradient and value contributions are
+    weighted using the ratio field in the followig way:
+
+    qConv = ratio*htcwfilm*(Tfilm - *this)*convectiveScaling_;
+    qRad = (1.0 - ratio)*Qr;
+
+    Then the solid can gain or loose energy through radiation or conduction
+    towards the film.
+
+    Note: kappa : heat conduction at patch.
+        Gets supplied how to lookup/calculate kappa:
+    - 'lookup' : lookup volScalarField (or volSymmTensorField) with name
+    - 'basicThermo' : use basicThermo and compressible::RASmodel to calculate K
+    - 'solidThermo' : use basicSolidThermo K()
+
+    Qr is the radiative flux defined in the radiation model.
+
+
+SourceFiles
+    filmPyrolysisRadiativeCoupledMixedFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmPyrolysisRadiativeCoupledMixedFvPatchScalarField_H
+#define filmPyrolysisRadiativeCoupledMixedFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+#include "temperatureCoupledBase.H"
+#include "thermoSingleLayer.H"
+#include "pyrolysisModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+   Class filmPyrolysisRadiativeCoupledMixedFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+:
+    public mixedFvPatchScalarField,
+    public temperatureCoupledBase
+{
+public:
+
+    typedef Foam::regionModels::surfaceFilmModels::thermoSingleLayer
+        filmModelType;
+
+    typedef Foam::regionModels::pyrolysisModels::pyrolysisModel
+        pyrolysisModelType;
+
+
+private:
+
+    // Private data
+
+        //- Name of field on the neighbour region
+        const word TnbrName_;
+
+         //- Name of the radiative heat flux in the neighbout region
+        const word QrNbrName_;
+
+        //- Name of the radiative heat flux in local region
+        const word QrName_;
+
+        //- Convective Scaling Factor (as determined by Prateep's tests)
+        const scalar convectiveScaling_;
+
+        //- Minimum delta film to be consired dry
+        const scalar filmDeltaDry_;
+
+        //- Maximum delta film to be consired wet
+        const scalar filmDeltaWet_;
+
+        //- Retrieve film model from the database
+        const filmModelType& filmModel() const;
+
+        //- Retrieve pyrolysis model from the database
+        const pyrolysisModelType& pyrModel() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("filmPyrolysisRadiativeCoupledMixed");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  turbulentTemperatureCoupledBaffleMixedFvPatchScalarField onto a
+        //  new patch
+        filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+        (
+            const
+            filmPyrolysisRadiativeCoupledMixedFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+        (
+            const filmPyrolysisRadiativeCoupledMixedFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new filmPyrolysisRadiativeCoupledMixedFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Get corresponding K field
+        tmp<scalarField> K() const;
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
index a9aef8f49f..0327dd1da7 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
@@ -69,8 +69,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     q_(p.size(), 0.0),
     h_(p.size(), 0.0),
     Ta_(p.size(), 0.0),
-    thicknessLayer_(0),
-    kappaLayer_(0)
+    thicknessLayers_(),
+    kappaLayers_()
 {
     this->refValue() = 0.0;
     this->refGrad() = 0.0;
@@ -93,8 +93,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     q_(ptf.q_, mapper),
     h_(ptf.h_, mapper),
     Ta_(ptf.Ta_, mapper),
-    thicknessLayer_(ptf.thicknessLayer_),
-    kappaLayer_(ptf.kappaLayer_)
+    thicknessLayers_(ptf.thicknessLayers_),
+    kappaLayers_(ptf.kappaLayers_)
 {}
 
 
@@ -112,8 +112,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     q_(p.size(), 0.0),
     h_(p.size(), 0.0),
     Ta_(p.size(), 0.0),
-    thicknessLayer_(dict.lookupOrDefault<scalar>("thicknessLayer", 0.0)),
-    kappaLayer_(dict.lookupOrDefault<scalar>("kappaLayer", 0.0))
+    thicknessLayers_(),
+    kappaLayers_()
 {
     if (dict.found("q") && !dict.found("h") && !dict.found("Ta"))
     {
@@ -125,6 +125,11 @@ externalWallHeatFluxTemperatureFvPatchScalarField
         mode_ = fixedHeatTransferCoeff;
         h_ = scalarField("h", dict, p.size());
         Ta_ = scalarField("Ta", dict, p.size());
+        if (dict.found("thicknessLayers"))
+        {
+            dict.lookup("thicknessLayers") >> thicknessLayers_;
+            dict.lookup("kappaLayers") >> kappaLayers_;
+        }
     }
     else
     {
@@ -176,8 +181,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     q_(tppsf.q_),
     h_(tppsf.h_),
     Ta_(tppsf.Ta_),
-    thicknessLayer_(tppsf.thicknessLayer_),
-    kappaLayer_(tppsf.kappaLayer_)
+    thicknessLayers_(tppsf.thicknessLayers_),
+    kappaLayers_(tppsf.kappaLayers_)
 {}
 
 
@@ -194,8 +199,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     q_(tppsf.q_),
     h_(tppsf.h_),
     Ta_(tppsf.Ta_),
-    thicknessLayer_(tppsf.thicknessLayer_),
-    kappaLayer_(tppsf.kappaLayer_)
+    thicknessLayers_(tppsf.thicknessLayers_),
+    kappaLayers_(tppsf.kappaLayers_)
 {}
 
 
@@ -252,7 +257,19 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         }
         case fixedHeatTransferCoeff:
         {
-            q = (Ta_ - Tp)*(1.0/h_ + thicknessLayer_/(kappaLayer_ + VSMALL));
+            scalar totalSolidRes = 0.0;
+            if (thicknessLayers_.size() > 0)
+            {
+                forAll (thicknessLayers_, iLayer)
+                {
+                    const scalar l = thicknessLayers_[iLayer];
+                    if (l > 0.0)
+                    {
+                        totalSolidRes += kappaLayers_[iLayer]/l;
+                    }
+                }
+            }
+            q = (Ta_ - Tp)*(1.0/h_ + totalSolidRes);
             break;
         }
         default:
@@ -308,9 +325,6 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::write
 {
     mixedFvPatchScalarField::write(os);
     temperatureCoupledBase::write(os);
-    os.writeKeyword("thicknessLayer")<< thicknessLayer_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("kappaLayer")<< kappaLayer_ << token::END_STATEMENT << nl;
 
     switch (mode_)
     {
@@ -323,6 +337,10 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::write
         {
             h_.writeEntry("h", os);
             Ta_.writeEntry("Ta", os);
+            os.writeKeyword("thicknessLayers")<< thicknessLayers_
+                << token::END_STATEMENT << nl;
+            os.writeKeyword("kappaLayers")<< kappaLayers_
+                << token::END_STATEMENT << nl;
             break;
         }
         default:
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H
index d8c4bc9a12..000bc8473c 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H
@@ -26,8 +26,9 @@ Class
 
 Description
     This boundary condition supplies a heat flux condition for temperature
-    on an external wall.Optional thin thermal layer resistance can be
-    specified through thicknessLayer and kappaLayer entries.
+    on an external wall. Optional thin thermal layer resistances can be
+    specified through thicknessLayers and kappaLayers entries for the
+    fixed heat transfer coefficient mode.
 
     The condition can operate in two modes:
         \li fixed heat transfer coefficient: supply h and Ta
@@ -52,8 +53,8 @@ Description
             q               uniform 1000;      // heat flux / [W/m2]
             Ta              uniform 300.0;     // ambient temperature /[K]
             h               uniform 10.0;      // heat transfer coeff /[W/Km2]
-            thicknessLayer  0.001              // thickness of layer [m]
-            kappaLayer      0.0                // thermal conductivity of
+            thicknessLayers (0.1 0.2 0.3 0.4); // thickness of layer [m]
+            kappaLayers     (1 2 3 4)          // thermal conductivity of
                                                // layer [W/m/K]
             value           uniform 300.0;     // initial temperature / [K]
             kappaName       none;
@@ -118,11 +119,11 @@ private:
         //- Ambient temperature / [K]
         scalarField Ta_;
 
-        //- Thickness of the thin wall
-        scalar thicknessLayer_;
+        //- Thickness of layers
+        scalarList thicknessLayers_;
 
-        //- Thermal conductivity of the thin wall
-        scalar kappaLayer_;
+        //- Conductivity of layers
+        scalarList kappaLayers_;
 
 
 public: