Added basic support for thermal effects in porous zones.

Currently the only supported model is fixing the temperature within the porous zone to a fixed value.

applications/solvers/compressible/rhoPorousSimpleFoam/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermalPorousZone/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
@@ -8,6 +9,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \
+    -lthermalPorousZone \
     -lspecie \
     -lcompressibleRASModels \
     -lfiniteVolume \

applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H

@@ -64,7 +64,7 @@
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
-    porousZones pZones(mesh);
+    thermalPorousZones pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     int nUCorr = 0;
@@ -84,4 +84,3 @@
             pressureImplicitPorosity = true;
         }
     }
-

applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H

@@ -9,6 +9,8 @@
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
+    pZones.addEnthalpySource(thermo, hEqn);
+
     hEqn.relax();
 
     eqnResidual = hEqn.solve().initialResidual();

applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -34,7 +34,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "RASModel.H"
-#include "porousZones.H"
+#include "thermalPorousZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "PorousZones.H"
+#include "Time.H"
+#include "volFields.H"
+#include "fvm.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ZoneType>
+template<class Type>
+void Foam::PorousZones<ZoneType>::modifyDdt(fvMatrix<Type>& m) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).modifyDdt(m);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ZoneType>
+Foam::PorousZones<ZoneType>::PorousZones
+(
+    const fvMesh& mesh
+)
+:
+    IOPtrList<ZoneType>
+    (
+        IOobject
+        (
+            "porousZones",
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        typename ZoneType::iNew(mesh)
+    ),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const oneField&,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const dimensionedScalar& rho,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const volScalarField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+template<class ZoneType>
+void Foam::PorousZones<ZoneType>::addResistance(fvVectorMatrix& UEqn) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn);
+    }
+}
+
+
+template<class ZoneType>
+void Foam::PorousZones<ZoneType>::addResistance
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    // addResistance for each zone, delaying the correction of the
+    // precessor BCs of AU
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn, AU, false);
+    }
+
+    // Correct the boundary conditions of the tensorial diagonal to ensure
+    // processor bounaries are correctly handled when AU^-1 is interpolated
+    // for the pressure equation.
+    AU.correctBoundaryConditions();
+}
+
+
+template<class ZoneType>
+bool Foam::PorousZones<ZoneType>::readData(Istream& is)
+{
+    this->clear();
+
+    IOPtrList<ZoneType> newLst
+    (
+        IOobject
+        (
+            "porousZones",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false     // Don't re-register new zones with objectRegistry
+        ),
+        typename ZoneType::iNew(mesh_)
+    );
+
+    transfer(newLst);
+
+    return is.good();
+}
+
+
+template<class ZoneType>
+bool Foam::PorousZones<ZoneType>::writeData(Ostream& os, bool subDict) const
+{
+    // Write size of list
+    os << nl << this->size();
+
+    // Write beginning of contents
+    os << nl << token::BEGIN_LIST;
+
+    // Write list contents
+    forAll(*this, i)
+    {
+        os << nl;
+        this->operator[](i).writeDict(os, subDict);
+    }
+
+    // Write end of contents
+    os << token::END_LIST << nl;
+
+    // Check state of IOstream
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::PorousZones<ZoneType>
+
+Description
+    A centralized ZoneType collection.
+
+    Container class for a set of ZoneType with the ZoneType member
+    functions implemented to loop over the functions for each ZoneType.
+
+SourceFiles
+    PorousZones.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PorousZones_H
+#define PorousZones_H
+
+#include "IOPtrList.H"
+
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+#include "dimensionedScalarFwd.H"
+#include "oneField.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+class fvMesh;
+
+
+/*---------------------------------------------------------------------------*\
+                           Class PorousZones Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ZoneType>
+class PorousZones
+:
+    public IOPtrList<ZoneType>
+{
+    // Private data
+
+        //- Reference to the finite volume mesh this zone is part of
+        const fvMesh& mesh_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        PorousZones(const PorousZones<ZoneType>&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PorousZones<ZoneType>&);
+
+
+        //- modify time derivative elements
+        template<class Type>
+        void modifyDdt(fvMatrix<Type>&) const;
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        //  with automatically constructed coordinate systems list
+        PorousZones(const fvMesh&);
+
+
+    // Member Functions
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const oneField&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const dimensionedScalar&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const volScalarField&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- Add the viscous and inertial resistance force contribution
+        //  to the momentum equation
+        void addResistance(fvVectorMatrix& UEqn) const;
+
+        //- Add the viscous and inertial resistance force contribution
+        //  to the tensorial diagonal
+        void addResistance
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU
+        ) const;
+
+        //- read modified data
+        virtual bool readData(Istream&);
+
+        //- write data
+        bool writeData(Ostream&, bool subDict = true) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PorousZones.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/porousMedia/porousZone.H

@@ -198,7 +198,6 @@ public:
             return autoPtr<porousZone>(NULL);
         }
 
-
         //- Return pointer to new porousZone created on freestore from Istream
         class iNew
         {
@@ -222,6 +221,11 @@ public:
         };
 
 
+    //- Destructor
+    virtual ~porousZone()
+    {}
+
+
     // Member Functions
 
         // Access
@@ -232,6 +236,12 @@ public:
                 return name_;
             }
 
+            //- Return mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
             //- cellZone number
             label zoneId() const
             {
@@ -275,7 +285,7 @@ public:
             }
 
 
-        //- modify time derivative elements according to porosity
+        //- Modify time derivative elements according to porosity
         template<class Type>
         void modifyDdt(fvMatrix<Type>&) const;
 
@@ -294,7 +304,7 @@ public:
         ) const;
 
         //- Write the porousZone dictionary
-        void writeDict(Ostream&, bool subDict = true) const;
+        virtual void writeDict(Ostream&, bool subDict = true) const;
 
 
     // Ostream Operator

src/finiteVolume/cfdTools/general/porousMedia/porousZones.C

@@ -25,8 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "porousZones.H"
-#include "Time.H"
-#include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -35,105 +33,4 @@ namespace Foam
     defineTemplateTypeNameAndDebug(IOPtrList<porousZone>, 0);
 }
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousZones::porousZones
-(
-    const fvMesh& mesh
-)
-:
-    IOPtrList<porousZone>
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh.time().constant(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        porousZone::iNew(mesh)
-    ),
-    mesh_(mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::porousZones::addResistance(fvVectorMatrix& UEqn) const
-{
-    forAll(*this, i)
-    {
-        operator[](i).addResistance(UEqn);
-    }
-}
-
-
-void Foam::porousZones::addResistance
-(
-    const fvVectorMatrix& UEqn,
-    volTensorField& AU
-) const
-{
-    // addResistance for each zone, delaying the correction of the
-    // precessor BCs of AU
-    forAll(*this, i)
-    {
-        operator[](i).addResistance(UEqn, AU, false);
-    }
-
-    // Correct the boundary conditions of the tensorial diagonal to ensure
-    // processor bounaries are correctly handled when AU^-1 is interpolated
-    // for the pressure equation.
-    AU.correctBoundaryConditions();
-}
-
-
-bool Foam::porousZones::readData(Istream& is)
-{
-    clear();
-
-    IOPtrList<porousZone> newLst
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false     // Don't re-register new zones with objectRegistry
-        ),
-        porousZone::iNew(mesh_)
-    );
-
-    transfer(newLst);
-
-    return is.good();
-}
-
-
-bool Foam::porousZones::writeData(Ostream& os, bool subDict) const
-{
-    // Write size of list
-    os << nl << size();
-
-    // Write beginning of contents
-    os << nl << token::BEGIN_LIST;
-
-    // Write list contents
-    forAll(*this, i)
-    {
-        os << nl;
-        operator[](i).writeDict(os, subDict);
-    }
-
-    // Write end of contents
-    os << token::END_LIST << nl;
-
-    // Check state of IOstream
-    return os.good();
-}
-
-
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/porousMedia/porousZones.H

@@ -22,152 +22,23 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Class
+Typedef
     Foam::porousZones
 
-Description
-    A centralized porousZone collection.
-
-    Container class for a set of porousZones with the porousZone member
-    functions implemented to loop over the functions for each porousZone.
-
-    The input file @c constant/porousZone is implemented as
-    IOPtrList\<porousZone\> and contains the following type of data:
-
-    @verbatim
-    1
-    (
-    cat1
-    {
-        coordinateSystem    system_10;
-        porosity    0.781;
-        Darcy
-        {
-            d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
-            f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
-        }
-    }
-    )
-    @endverbatim
-
-SourceFiles
-    porousZones.C
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef porousZones_H
 #define porousZones_H
 
+#include "PorousZones.H"
 #include "porousZone.H"
-#include "IOPtrList.H"
-
-#include "volFieldsFwd.H"
-#include "fvMatrix.H"
-#include "oneField.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
+    typedef PorousZones<porousZone> porousZones;
-/*---------------------------------------------------------------------------*\
+}
-                           Class porousZones Declaration
-\*---------------------------------------------------------------------------*/
-
-class porousZones
-:
-    public IOPtrList<porousZone>
-{
-    // Private data
-
-        //- Reference to the finite volume mesh this zone is part of
-        const fvMesh& mesh_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        porousZones(const porousZones&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const porousZones&);
-
-
-        //- modify time derivative elements
-        template<class Type>
-        void modifyDdt(fvMatrix<Type>&) const;
-
-public:
-
-    // Constructors
-
-        //- Construct from fvMesh
-        //  with automatically constructed coordinate systems list
-        porousZones(const fvMesh&);
-
-
-    // Member Functions
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const oneField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation
-        void addResistance(fvVectorMatrix& UEqn) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the tensorial diagonal
-        void addResistance
-        (
-            const fvVectorMatrix& UEqn,
-            volTensorField& AU
-        ) const;
-
-        //- read modified data
-        virtual bool readData(Istream&);
-
-        //- write data
-        bool writeData(Ostream&, bool subDict = true) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "porousZonesTemplates.C"
-#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/Allwmake

@@ -16,5 +16,6 @@ wmake libso chemistryModel
 wmake libso pdfs
 wmake libso radiation
 wmake libso barotropicCompressibilityModel
+wmake libso thermalPorousZone
 
 # ----------------------------------------------------------------- end-of-file

src/thermophysicalModels/thermalPorousZone/Make/files

@@ -0,0 +1,4 @@
+thermalPorousZone/thermalPorousZone.C
+thermalPorousZone/thermalPorousZones.C
+
+LIB = $(FOAM_LIBBIN)/libthermalPorousZone

src/thermophysicalModels/thermalPorousZone/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+
+LIB_LIBS = \
+    -lbasicThermophysicalModels \
+    -lmeshTools \
+    -lfiniteVolume

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "thermalPorousZone.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::thermalPorousZone::thermalPorousZone
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    porousZone(name, mesh, dict),
+    T_("T", dimTemperature, -GREAT)
+{
+    if (const dictionary* dictPtr = dict.subDictPtr("thermalModel"))
+    {
+        word thermalModel(dictPtr->lookup("type"));
+
+        if (thermalModel == "fixedTemperature")
+        {
+            dictPtr->lookup("T") >> T_;
+        }
+        else
+        {
+            FatalIOErrorIn
+            (
+                "thermalPorousZone::thermalPorousZone"
+                "("
+                "const word& name, "
+                "const fvMesh& mesh, "
+                "const dictionary& dict"
+                ")",
+                *dictPtr
+            )   << "thermalModel " << thermalModel << " is not supported" << nl
+                << "    Supported thermalModels are: fixedTemperature"
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::thermalPorousZone::addEnthalpySource
+(
+    const basicThermo& thermo,
+    fvScalarMatrix& hEqn
+) const
+{
+    if (zoneId() == -1 || T_.value() < 0.0)
+    {
+        return;
+    }
+
+    const labelList& cells = mesh().cellZones()[zoneId()];
+    const scalarField& V = mesh().V();
+    scalarField& hDiag = hEqn.diag();
+    scalarField& hSource = hEqn.source();
+    const scalarField& rho = thermo.rho();
+
+    scalarField hZone = thermo.h(scalarField(cells.size(), T_.value()), cells);
+    scalar rate = 1e6;
+
+    forAll (cells, i)
+    {
+        hDiag[cells[i]] += rate*V[cells[i]]*rho[cells[i]];
+        hSource[cells[i]] += rate*V[cells[i]]*rho[cells[i]]*hZone[i];
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::thermalPorousZone
+
+Description
+    Porous zone definition based on cell zones including terms for energy
+    equations.
+
+See Also
+    porousZone, thermalPorousZones and coordinateSystems
+
+SourceFiles
+    thermalPorousZone.C
+    thermalPorousZoneTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermalPorousZone_H
+#define thermalPorousZone_H
+
+#include "porousZone.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+class basicThermo;
+
+/*---------------------------------------------------------------------------*\
+                        Class thermalPorousZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermalPorousZone
+:
+    public porousZone
+{
+    // Private data
+
+        //- Temperature in the porous-zone
+        dimensionedScalar T_;
+
+
+        //- Disallow default bitwise copy construct
+        thermalPorousZone(const thermalPorousZone&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const thermalPorousZone&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        thermalPorousZone(const word& name, const fvMesh&, const dictionary&);
+
+        //- Return clone
+        autoPtr<thermalPorousZone> clone() const
+        {
+            notImplemented("autoPtr<thermalPorousZone> clone() const");
+            return autoPtr<thermalPorousZone>(NULL);
+        }
+
+        //- Return pointer to new thermalPorousZone
+        //  created on freestore from Istream
+        class iNew
+        {
+            //- Reference to the finite volume mesh this zone is part of
+            const fvMesh& mesh_;
+
+        public:
+
+            iNew(const fvMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<thermalPorousZone> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+
+                return autoPtr<thermalPorousZone>
+                (
+                    new thermalPorousZone(name, mesh_, dict)
+                );
+            }
+        };
+
+
+    //- Destructor
+    virtual ~thermalPorousZone()
+    {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the temperature in the porous-zone
+            const dimensionedScalar& T() const
+            {
+                return T_;
+            }
+
+            //- Edit access to the temperature in the porous-zone
+            dimensionedScalar& T()
+            {
+                return T_;
+            }
+
+        //- Add the thermal source to the enthalpy equation
+        void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+//#   include "thermalPorousZoneTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZoneTemplates.C

@@ -24,77 +24,57 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "porousZones.H"
+#include "porousZone.H"
-#include "volFields.H"
+#include "fvMesh.H"
-#include "fvMatrix.H"
-#include "fvm.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class Type>
+template<class RhoFieldType>
-void Foam::porousZones::modifyDdt(fvMatrix<Type>& m) const
+void Foam::porousZone::addPowerLawResistance
-{
-    forAll(*this, i)
-    {
-        operator[](i).modifyDdt(m);
-    }
-}
-
-
-// * * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    scalarField& Udiag,
-)
+    const labelList& cells,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
 {
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    const scalar C0 = C0_;
-    modifyDdt(tres());
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
-    return tres;
+
+    forAll (cells, i)
+    {
+        Udiag[cells[i]] +=
+            V[cells[i]]*rho[cells[i]]*C0*pow(magSqr(U[cells[i]]), C1m1b2);
+    }
 }
 
 
-template<class Type>
+template<class RhoFieldType>
-Foam::tmp<Foam::fvMatrix<Type> >
+void Foam::porousZone::addViscousInertialResistance
-Foam::porousZones::ddt
 (
-    const oneField&,
+    scalarField& Udiag,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    vectorField& Usource,
-)
+    const labelList& cells,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const scalarField& mu,
+    const vectorField& U
+) const
 {
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    const tensor& D = D_.value();
-    modifyDdt(tres());
+    const tensor& F = F_.value();
-    return tres;
+
+    forAll (cells, i)
+    {
+        tensor dragCoeff = mu[cells[i]]*D + (rho[cells[i]]*mag(U[cells[i]]))*F;
+        scalar isoDragCoeff = tr(dragCoeff);
+
+        Udiag[cells[i]] += V[cells[i]]*isoDragCoeff;
+        Usource[cells[i]] -=
+            V[cells[i]]*((dragCoeff - I*isoDragCoeff) & U[cells[i]]);
+    }
 }
 
 
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
-(
-    const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
-(
-    const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
-}
-
 // ************************************************************************* //

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "thermalPorousZones.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(IOPtrList<thermalPorousZone>, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::thermalPorousZones::thermalPorousZones
+(
+    const fvMesh& mesh
+)
+:
+    PorousZones<thermalPorousZone>(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::thermalPorousZones::addEnthalpySource
+(
+    const basicThermo& thermo,
+    fvScalarMatrix& hEqn
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).addEnthalpySource(thermo, hEqn);
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::thermalPorousZones
+
+Description
+    A centralized thermalPorousZone collection.
+
+    Container class for a set of thermalPorousZones with the thermalPorousZone
+    member functions implemented to loop over the functions for each
+    thermalPorousZone.
+
+    The input file @c constant/thermalPorousZone is implemented as
+    IOPtrList\<thermalPorousZone\> and contains the following type of data:
+
+    @verbatim
+    1
+    (
+    cat1
+    {
+        coordinateSystem    system_10;
+        porosity    0.781;
+        Darcy
+        {
+            d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
+            f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
+        }
+        Temperature [0 0 1 0 0] 600;
+    }
+    )
+    @endverbatim
+
+SourceFiles
+    thermalPorousZones.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermalPorousZones_H
+#define thermalPorousZones_H
+
+#include "PorousZones.H"
+#include "thermalPorousZone.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class thermalPorousZones Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermalPorousZones
+:
+    public PorousZones<thermalPorousZone>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        thermalPorousZones(const fvMesh&);
+
+
+    // Member Functions
+
+        //- Add the thermal source to the enthalpy equation
+        void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/angledDuct/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary

@@ -1,14 +1,14 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
-    format      binary;
+    format      ascii;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;

tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones

@@ -30,6 +30,13 @@ FoamFile
             d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
             f   f [0 -1 0 0 0 0 0] (0 0 0);
         }
+
+        thermalModel
+        {
+            type fixedTemperature;
+
+            T  T [0 0 0 1 0] 350;
+        }
     }
 )