diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 58c0b576737cd27884fa4063967ebeab8ab871d8..b8b7e2f1d4068011822bab887aef7c3777a062ad 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -365,8 +365,12 @@ $(solutionControl)/simpleControl/simpleControl.C
$(solutionControl)/pimpleControl/pimpleControl.C
porousMedia = $(general)/porousMedia
-$(porousMedia)/porousZone.C
-$(porousMedia)/porousZones.C
+$(porousMedia)/porosityModel/porosityModel/porosityModel.C
+$(porousMedia)/porosityModel/porosityModel/porosityModelNew.C
+$(porousMedia)/porosityModel/porosityModel/porosityModelList.C
+$(porousMedia)/porosityModel/porosityModel/IOporosityModelList.C
+$(porousMedia)/porosityModel/DarcyForchheimer/DarcyForchheimer.C
+$(porousMedia)/porosityModel/powerLaw/powerLaw.C
MRF = $(general)/MRF
$(MRF)/MRFZone.C
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C
new file mode 100644
index 0000000000000000000000000000000000000000..f292cb63c85741f2f12496502b9baff933ffa945
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C
@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "addToRunTimeSelectionTable.H"
+#include "DarcyForchheimer.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ namespace porosityModels
+ {
+ defineTypeNameAndDebug(DarcyForchheimer, 0);
+ addToRunTimeSelectionTable(porosityModel, DarcyForchheimer, mesh);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
+(
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ porosityModel(name, modelType, mesh, dict),
+ coordSys_(coeffs_, mesh),
+ D_("D", dimless/sqr(dimLength), tensor::zero),
+ F_("F", dimless/dimLength, tensor::zero),
+ rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
+ muName_(coeffs_.lookupOrDefault<word>("mu", "mu")),
+ nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
+{
+ // local-to-global transformation tensor
+ const tensor& E = coordSys_.R();
+
+ dimensionedVector d(coeffs_.lookup("d"));
+ if (D_.dimensions() != d.dimensions())
+ {
+ FatalIOErrorIn
+ (
+ "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
+ "("
+ "const dictionary&, "
+ "const coordinateSystem&, "
+ "const keyType&"
+ ")",
+ coeffs_
+ ) << "incorrect dimensions for d: " << d.dimensions()
+ << " should be " << D_.dimensions()
+ << exit(FatalIOError);
+ }
+
+ adjustNegativeResistance(d);
+
+ D_.value().xx() = d.value().x();
+ D_.value().yy() = d.value().y();
+ D_.value().zz() = d.value().z();
+ D_.value() = (E & D_ & E.T()).value();
+
+ dimensionedVector f(coeffs_.lookup("f"));
+ if (F_.dimensions() != f.dimensions())
+ {
+ FatalIOErrorIn
+ (
+ "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer"
+ "("
+ "const dictionary&, "
+ "const coordinateSystem&, "
+ "const keyType&"
+ ")",
+ coeffs_
+ ) << "incorrect dimensions for f: " << f.dimensions()
+ << " should be " << F_.dimensions()
+ << exit(FatalIOError);
+ }
+
+ adjustNegativeResistance(f);
+
+ // leading 0.5 is from 1/2*rho
+ F_.value().xx() = 0.5*f.value().x();
+ F_.value().yy() = 0.5*f.value().y();
+ F_.value().zz() = 0.5*f.value().z();
+ F_.value() = (E & F_ & E.T()).value();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::DarcyForchheimer::~DarcyForchheimer()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+ fvVectorMatrix& UEqn
+) const
+{
+ const vectorField& U = UEqn.psi();
+ const scalarField& V = mesh_.V();
+ scalarField& Udiag = UEqn.diag();
+ vectorField& Usource = UEqn.source();
+
+ if (UEqn.dimensions() == dimForce)
+ {
+ const volScalarField& rho =
+ mesh_.lookupObject<volScalarField>(rhoName_);
+ const volScalarField& mu =
+ mesh_.lookupObject<volScalarField>(muName_);
+
+ apply(Udiag, Usource, V, rho, mu, U);
+ }
+ else
+ {
+ const volScalarField& nu =
+ mesh_.lookupObject<volScalarField>(nuName_);
+
+ apply(Udiag, Usource, V, geometricOneField(), nu, U);
+ }
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+) const
+{
+ const vectorField& U = UEqn.psi();
+ const scalarField& V = mesh_.V();
+ scalarField& Udiag = UEqn.diag();
+ vectorField& Usource = UEqn.source();
+
+ apply(Udiag, Usource, V, rho, mu, U);
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::correct
+(
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU
+) const
+{
+ const vectorField& U = UEqn.psi();
+
+ if (UEqn.dimensions() == dimForce)
+ {
+ const volScalarField& rho =
+ mesh_.lookupObject<volScalarField>(rhoName_);
+ const volScalarField& mu =
+ mesh_.lookupObject<volScalarField>(muName_);
+
+ apply(AU, rho, mu, U);
+ }
+ else
+ {
+ const volScalarField& nu =
+ mesh_.lookupObject<volScalarField>(nuName_);
+
+ apply(AU, geometricOneField(), nu, U);
+ }
+}
+
+
+void Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
+{
+ os << indent << name_ << endl;
+ dict_.write(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H
new file mode 100644
index 0000000000000000000000000000000000000000..928422938a1da0357a1563335e93f20b618409f4
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H
@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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::DarcyForchheimer
+
+Description
+ Darcy-Forchheimer law porosity model, given by:
+
+ \f[
+ S = - (\mu d + \frac{\rho |U|}{2} f) U
+ \f]
+
+ where
+ \vartable
+ d | Darcy coefficient [1/m2]
+ f | Forchheimer coefficient [1/m]
+ \endvartable
+
+ Since negative Darcy/Forchheimer parameters are invalid, they can be used
+ to specify a multiplier (of the max component).
+
+ The orientation of the porous region is defined with the same notation as
+ a co-ordinate system, but only a Cartesian co-ordinate system is valid.
+
+SourceFiles
+ DarcyForchheimer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DarcyForchheimer_H
+#define DarcyForchheimer_H
+
+#include "porosityModel.H"
+#include "coordinateSystem.H"
+#include "dimensionedTensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+ Class DarcyForchheimer Declaration
+\*---------------------------------------------------------------------------*/
+
+class DarcyForchheimer
+:
+ public porosityModel
+{
+private:
+
+ // Private data
+
+ //- Local co-ordinate system
+ coordinateSystem coordSys_;
+
+ //- Darcy coefficient [1/m2]
+ dimensionedTensor D_;
+
+ //- Forchheimer coefficient [1/m]
+ dimensionedTensor F_;
+
+ //- Name of density field
+ word rhoName_;
+
+ //- Name of dynamic viscosity field
+ word muName_;
+
+ //- Name of kinematic viscosity field
+ word nuName_;
+
+
+ // Private Member Functions
+
+ //- Apply
+ template<class RhoFieldType>
+ void apply
+ (
+ scalarField& Udiag,
+ vectorField& Usource,
+ const scalarField& V,
+ const RhoFieldType& rho,
+ const scalarField& mu,
+ const vectorField& U
+ ) const;
+
+ //- Apply
+ template<class RhoFieldType>
+ void apply
+ (
+ tensorField& AU,
+ const RhoFieldType& rho,
+ const scalarField& mu,
+ const vectorField& U
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ DarcyForchheimer(const DarcyForchheimer&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const DarcyForchheimer&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("DarcyForchheimer");
+
+ //- Constructor
+ DarcyForchheimer
+ (
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Destructor
+ virtual ~DarcyForchheimer();
+
+
+ // Member Functions
+
+ //- Add resistance
+ virtual void correct(fvVectorMatrix& UEqn) const;
+
+ //- Add resistance
+ virtual void correct
+ (
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+ ) const;
+
+ //- Add resistance
+ virtual void correct
+ (
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU
+ ) const;
+
+
+ // I-O
+
+ //- Write
+ void writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "DarcyForchheimerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..184206a010dccd166878507f0b002e1e8d35f9b9
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::DarcyForchheimer::apply
+(
+ scalarField& Udiag,
+ vectorField& Usource,
+ const scalarField& V,
+ const RhoFieldType& rho,
+ const scalarField& mu,
+ const vectorField& U
+) const
+{
+ const tensor& D = D_.value();
+ const tensor& F = F_.value();
+
+ forAll(cellZoneIds_, zoneI)
+ {
+ const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+ forAll(cells, i)
+ {
+ const label cellI = cells[i];
+
+ const tensor Cd = mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+
+ const scalar isoCd = tr(Cd);
+
+ Udiag[cellI] += V[cellI]*isoCd;
+ Usource[cellI] -= V[cellI]*((Cd - I*isoCd) & U[cellI]);
+ }
+ }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::DarcyForchheimer::apply
+(
+ tensorField& AU,
+ const RhoFieldType& rho,
+ const scalarField& mu,
+ const vectorField& U
+) const
+{
+ const tensor& D = D_.value();
+ const tensor& F = F_.value();
+
+ forAll(cellZoneIds_, zoneI)
+ {
+ const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+ forAll(cells, i)
+ {
+ const label cellI = cells[i];
+ AU[cellI] += mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C
new file mode 100644
index 0000000000000000000000000000000000000000..2a655eefddec9d7fedb8efc01ac4ba1f16c986d9
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.C
@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "IOporosityModelList.H"
+#include "fvMesh.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::IOobject Foam::IOporosityModelList::createIOobject
+(
+ const fvMesh& mesh
+) const
+{
+ IOobject io
+ (
+ "porosityProperties",
+ mesh.time().constant(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
+
+ if (io.headerOk())
+ {
+ Info<< "Creating porosity model list from " << io.name() << nl << endl;
+
+ io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+ return io;
+ }
+ else
+ {
+ Info<< "No porosity models present" << nl << endl;
+
+ io.readOpt() = IOobject::NO_READ;
+ return io;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::IOporosityModelList::IOporosityModelList
+(
+ const fvMesh& mesh
+)
+:
+ IOdictionary(createIOobject(mesh)),
+ porosityModelList(mesh, *this)
+{}
+
+
+bool Foam::IOporosityModelList::read()
+{
+ if (regIOobject::read())
+ {
+ porosityModelList::read(*this);
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+// ************************************************************************* //
+
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H
new file mode 100644
index 0000000000000000000000000000000000000000..e9e534926ab338adf3107a0a13bfab93f3ec368c
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/IOporosityModelList.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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::IOporosityModelList
+
+Description
+ List of porosity models with IO functionality
+
+SourceFiles
+ IOporosityModelList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOporosityModelList_H
+#define IOporosityModelList_H
+
+#include "IOdictionary.H"
+#include "porosityModelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class IOporosityModelList Declaration
+\*---------------------------------------------------------------------------*/
+
+class IOporosityModelList
+:
+ public IOdictionary,
+ public porosityModelList
+{
+private:
+
+ // Private Member Functions
+
+ //- Create IO object if dictionary is present
+ IOobject createIOobject(const fvMesh& mesh) const;
+
+ //- Disallow default bitwise copy construct
+ IOporosityModelList(const IOporosityModelList&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const IOporosityModelList&);
+
+
+public:
+
+ // Constructors
+
+ //- Construct from mesh
+ IOporosityModelList(const fvMesh& mesh);
+
+
+ //- Destructor
+ virtual ~IOporosityModelList()
+ {}
+
+
+ // Member Functions
+
+ //- Read dictionary
+ virtual bool read();
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..f0997c822211a116787e1990364ad82955b952f0
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "porosityModel.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(porosityModel, 0);
+ defineRunTimeSelectionTable(porosityModel, mesh);
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
+
+void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
+{
+ scalar maxCmpt = max(0, cmptMax(resist.value()));
+
+ if (maxCmpt < 0)
+ {
+ FatalErrorIn
+ (
+ "void Foam::porosityModel::adjustNegativeResistance"
+ "("
+ "dimensionedVector&"
+ ")"
+ ) << "Negative resistances are invalid, resistance = " << resist
+ << exit(FatalError);
+ }
+ else
+ {
+ vector& val = resist.value();
+ for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
+ {
+ if (val[cmpt] < 0)
+ {
+ val[cmpt] *= -maxCmpt;
+ }
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::porosityModel::porosityModel
+(
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ name_(name),
+ mesh_(mesh),
+ dict_(dict),
+ coeffs_(dict.subDict(modelType + "Coeffs")),
+ active_(readBool(dict_.lookup("active"))),
+ zoneName_(dict_.lookup("cellZone")),
+ cellZoneIds_(mesh_.cellZones().findIndices(zoneName_))
+{
+ Info<< " creating porous zone: " << zoneName_ << endl;
+
+ bool foundZone = !cellZoneIds_.empty();
+ reduce(foundZone, orOp<bool>());
+
+ if (!foundZone && Pstream::master())
+ {
+ FatalErrorIn
+ (
+ "Foam::porosityModel::porosityModel"
+ "("
+ "const word&, "
+ "const word&, "
+ "const fvMesh&, "
+ "const dictionary&"
+ ")"
+ ) << "cannot find porous cellZone " << zoneName_
+ << exit(FatalError);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::porosityModel::~porosityModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+ fvVectorMatrix& UEqn
+) const
+{
+ if (cellZoneIds_.empty())
+ {
+ return;
+ }
+
+ this->correct(UEqn);
+}
+
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+) const
+{
+ if (cellZoneIds_.empty())
+ {
+ return;
+ }
+
+ this->correct(UEqn, rho, mu);
+}
+
+
+void Foam::porosityModel::porosityModel::addResistance
+(
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU,
+ bool correctAUprocBC
+) const
+{
+ if (cellZoneIds_.empty())
+ {
+ return;
+ }
+
+ this->correct(UEqn, AU);
+
+ if (correctAUprocBC)
+ {
+ // Correct the boundary conditions of the tensorial diagonal to ensure
+ // processor boundaries are correctly handled when AU^-1 is interpolated
+ // for the pressure equation.
+ AU.correctBoundaryConditions();
+ }
+}
+
+
+bool Foam::porosityModel::read(const dictionary& dict)
+{
+ active_ = readBool(dict.lookup("active"));
+ coeffs_ = dict.subDict(type() + "Coeffs");
+ dict.lookup("cellZone") >> zoneName_;
+ cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..56131ba3dd97cd5c287f699acc9eb3f9fcce153c
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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::porosityModel
+
+Description
+ Top level model for porosity models
+
+SourceFiles
+ porosityModel.C
+ porosityModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef porosityModel_H
+#define porosityModel_H
+
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "fvMatricesFwd.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedVector.H"
+#include "keyType.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class porosityModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class porosityModel
+{
+private:
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ porosityModel(const porosityModel&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const porosityModel&);
+
+
+protected:
+
+ // Protected data
+
+ //- Porosity name
+ word name_;
+
+ //- Reference to the mesh database
+ const fvMesh& mesh_;
+
+ //- Dictionary used for model construction
+ const dictionary dict_;
+
+ //- Model coefficients dictionary
+ dictionary coeffs_;
+
+ //- Porosity active flag
+ bool active_;
+
+ //- Name(s) of cell-zone
+ keyType zoneName_;
+
+ //- Cell zone Ids
+ labelList cellZoneIds_;
+
+
+ // Protected Member Functions
+
+ //- Adjust negative resistance values to be multiplier of max value
+ void adjustNegativeResistance(dimensionedVector& resist);
+
+ virtual void correct(fvVectorMatrix& UEqn) const = 0;
+
+ virtual void correct
+ (
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+ ) const = 0;
+
+ virtual void correct
+ (
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU
+ ) const = 0;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("porosityModel");
+
+ //- Selection table
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ porosityModel,
+ mesh,
+ (
+ const word& modelName,
+ const word& name,
+ const fvMesh& mesh,
+ const dictionary& dict
+ ),
+ (modelName, name, mesh, dict)
+ );
+
+ //- Constructor
+ porosityModel
+ (
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Return pointer to new porosityModel object created on the freestore
+ // from an Istream
+ class iNew
+ {
+ //- Reference to the mesh database
+ const fvMesh& mesh_;
+ const word& name_;
+
+ public:
+
+ iNew
+ (
+ const fvMesh& mesh,
+ const word& name
+ )
+ :
+ mesh_(mesh),
+ name_(name)
+ {}
+
+ autoPtr<porosityModel> operator()(Istream& is) const
+ {
+ const dictionary dict(is);
+
+ return autoPtr<porosityModel>
+ (
+ porosityModel::New
+ (
+ name_,
+ mesh_,
+ dict
+ )
+ );
+ }
+ };
+
+ //- Selector
+ static autoPtr<porosityModel> New
+ (
+ const word& name,
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Destructor
+ virtual ~porosityModel();
+
+
+ // Member Functions
+
+ //- Return const access to the porosity model name
+ inline const word& name() const;
+
+ //- Return const access to the porosity active flag
+ inline bool active() const;
+
+ //- Add resistance
+ virtual void addResistance(fvVectorMatrix& UEqn) const;
+
+ //- Add resistance
+ virtual void addResistance
+ (
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+ ) const;
+
+ //- Add resistance
+ virtual void addResistance
+ (
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU,
+ bool correctAUprocBC
+ ) const;
+
+
+ // I-O
+
+ //- Write
+ virtual void writeData(Ostream& os) const = 0;
+
+ //- Read porosity dictionary
+ virtual bool read(const dictionary& dict);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "porosityModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H
new file mode 100644
index 0000000000000000000000000000000000000000..1df2d1197298852ceb8523be29a096f87085b6c5
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H
@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::word& Foam::porosityModel::name() const
+{
+ return name_;
+}
+
+
+inline bool Foam::porosityModel::active() const
+{
+ return active_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C
new file mode 100644
index 0000000000000000000000000000000000000000..14ab39d96436521ad20c9db8cdeef36972d4fc86
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C
@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "porosityModelList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
+/*
+void Foam::porosityModelList::XXX()
+{}
+*/
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::porosityModelList::porosityModelList
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ PtrList<porosityModel>(),
+ mesh_(mesh)
+{
+ reset(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::porosityModelList::~porosityModelList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::porosityModelList::active() const
+{
+ bool a = false;
+ forAll(*this, i)
+ {
+ a = a || this->operator[](i).active();
+ }
+
+ if (!a)
+ {
+ Info<< "No porosity models active" << endl;
+ }
+
+ return a;
+}
+
+
+void Foam::porosityModelList::reset(const dictionary& dict)
+{
+ label count = 0;
+ forAllConstIter(dictionary, dict, iter)
+ {
+ if (iter().isDict())
+ {
+ count++;
+ }
+ }
+
+ this->setSize(count);
+ label i = 0;
+ forAllConstIter(dictionary, dict, iter)
+ {
+ if (iter().isDict())
+ {
+ const word& name = iter().keyword();
+ const dictionary& modelDict = iter().dict();
+
+ this->set
+ (
+ i++,
+ porosityModel::New(name, mesh_, modelDict)
+ );
+ }
+ }
+}
+
+
+bool Foam::porosityModelList::read(const dictionary& dict)
+{
+ bool allOk = true;
+ forAll(*this, i)
+ {
+ porosityModel& pm = this->operator[](i);
+ bool ok = pm.read(dict.subDict(pm.name()));
+ allOk = (allOk && ok);
+ }
+ return allOk;
+}
+
+
+bool Foam::porosityModelList::writeData(Ostream& os) const
+{
+ forAll(*this, i)
+ {
+ os << nl;
+ this->operator[](i).writeData(os);
+ }
+
+ return os.good();
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+ fvVectorMatrix& UEqn
+) const
+{
+ forAll(*this, i)
+ {
+ this->operator[](i).addResistance(UEqn);
+ }
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+) const
+{
+ forAll(*this, i)
+ {
+ this->operator[](i).addResistance(UEqn, rho, mu);
+ }
+}
+
+
+void Foam::porosityModelList::addResistance
+(
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU,
+ bool correctAUprocBC
+) const
+{
+ forAll(*this, i)
+ {
+ this->operator[](i).addResistance(UEqn, AU, correctAUprocBC);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const porosityModelList& models
+)
+{
+ models.writeData(os);
+ return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H
new file mode 100644
index 0000000000000000000000000000000000000000..e580a2d0fadd7c46a0df99d1acc83f4c00e203a0
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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::porosityModelList
+
+Description
+ List container for porosity models
+
+SourceFiles
+ porosityModelList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef porosityModelList_H
+#define porosityModelList_H
+
+#include "fvMesh.H"
+#include "dictionary.H"
+#include "fvMatricesFwd.H"
+#include "porosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+class porosityModelList;
+Ostream& operator<<(Ostream& os, const porosityModelList& models);
+
+/*---------------------------------------------------------------------------*\
+ Class porosityModelList Declaration
+\*---------------------------------------------------------------------------*/
+
+class porosityModelList
+:
+ PtrList<porosityModel>
+{
+private:
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ porosityModelList(const porosityModelList&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const porosityModelList&);
+
+
+protected:
+
+ // Protected data
+
+ //- Reference to the mesh database
+ const fvMesh& mesh_;
+
+
+public:
+
+ //- Constructor
+ porosityModelList(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ ~porosityModelList();
+
+
+ // Member Functions
+
+ //- Return active status
+ bool active() const;
+
+ //- Reset the source list
+ void reset(const dictionary& dict);
+
+ //- Add resistance
+ void addResistance(fvVectorMatrix& UEqn) const;
+
+ //- Add resistance
+ void addResistance
+ (
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+ ) const;
+
+ //- Add resistance
+ void addResistance
+ (
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU,
+ bool correctAUprocBC = true
+ ) const;
+
+
+ // I-O
+
+ //- Read dictionary
+ bool read(const dictionary& dict);
+
+ //- Write data to Ostream
+ bool writeData(Ostream& os) const;
+
+ //- Ostream operator
+ friend Ostream& operator<<
+ (
+ Ostream& os,
+ const porosityModelList& models
+ );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..1da8c91c4dfb583e022ce27059a0abaee880e63e
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelNew.C
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "porosityModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::porosityModel> Foam::porosityModel::New
+(
+ const word& name,
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+{
+ const word modelType(dict.lookup("type"));
+
+ Info<< "Porosity region " << name << ":" << nl
+ << " selecting model: " << modelType << endl;
+
+ meshConstructorTable::iterator cstrIter =
+ meshConstructorTablePtr_->find(modelType);
+
+ if (cstrIter == meshConstructorTablePtr_->end())
+ {
+ FatalErrorIn
+ (
+ "porosityModel::New"
+ "("
+ "const word& name,"
+ "const fvMesh&, "
+ "const dictionary&"
+ ")"
+ )
+ << "Unknown " << typeName << " type " << modelType << nl << nl
+ << "Valid " << typeName << " types are:" << nl
+ << meshConstructorTablePtr_->sortedToc()
+ << exit(FatalError);
+ }
+
+ return autoPtr<porosityModel>(cstrIter()(name, modelType, mesh, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C
new file mode 100644
index 0000000000000000000000000000000000000000..15fbae39dfb78486fe133002bb94ef1990956ec1
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 "addToRunTimeSelectionTable.H"
+#include "powerLaw.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ namespace porosityModels
+ {
+ defineTypeNameAndDebug(powerLaw, 0);
+ addToRunTimeSelectionTable(porosityModel, powerLaw, mesh);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLaw::powerLaw
+(
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ porosityModel(name, modelType, mesh, dict),
+ C0_(readScalar(coeffs_.lookup("C0"))),
+ C1_(readScalar(coeffs_.lookup("C1"))),
+ rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLaw::~powerLaw()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::porosityModels::powerLaw::correct
+(
+ fvVectorMatrix& UEqn
+) const
+{
+ const vectorField& U = UEqn.psi();
+ const scalarField& V = mesh_.V();
+ scalarField& Udiag = UEqn.diag();
+
+ if (UEqn.dimensions() == dimForce)
+ {
+ const volScalarField& rho =
+ mesh_.lookupObject<volScalarField>(rhoName_);
+
+ apply(Udiag, V, rho, U);
+ }
+ else
+ {
+ apply(Udiag, V, geometricOneField(), U);
+ }
+}
+
+
+void Foam::porosityModels::powerLaw::correct
+(
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+) const
+{
+ const vectorField& U = UEqn.psi();
+ const scalarField& V = mesh_.V();
+ scalarField& Udiag = UEqn.diag();
+
+ apply(Udiag, V, rho, U);
+}
+
+
+void Foam::porosityModels::powerLaw::correct
+(
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU
+) const
+{
+ const vectorField& U = UEqn.psi();
+
+ if (UEqn.dimensions() == dimForce)
+ {
+ const volScalarField& rho =
+ mesh_.lookupObject<volScalarField>(rhoName_);
+
+ apply(AU, rho, U);
+ }
+ else
+ {
+ apply(AU, geometricOneField(), U);
+ }
+}
+
+
+void Foam::porosityModels::powerLaw::writeData(Ostream& os) const
+{
+ os << indent << name_ << endl;
+ dict_.write(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H
new file mode 100644
index 0000000000000000000000000000000000000000..070a00e9a021d4dc49d743f240f9e08b5edf03e5
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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::powerLaw
+
+Description
+ Power law porosity model, given by:
+
+ \f[
+ S = - \rho C_0 |U|^{(C_1 - 1)} U
+ \f]
+
+ where
+ \vartable
+ C_0 | model coefficient
+ C_1 | model coefficient
+ \endvartable
+
+
+SourceFiles
+ powerLaw.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef powerLaw_H
+#define powerLaw_H
+
+#include "porosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+ Class powerLaw Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLaw
+:
+ public porosityModel
+{
+private:
+
+ // Private data
+
+ //- C0 coefficient
+ scalar C0_;
+
+ //- C1 coefficient
+ scalar C1_;
+
+ //- Name of density field
+ word rhoName_;
+
+
+ // Private Member Functions
+
+ //- Apply resistance
+ template<class RhoFieldType>
+ void apply
+ (
+ scalarField& Udiag,
+ const scalarField& V,
+ const RhoFieldType& rho,
+ const vectorField& U
+ ) const;
+
+ //- Apply resistance
+ template<class RhoFieldType>
+ void apply
+ (
+ tensorField& AU,
+ const RhoFieldType& rho,
+ const vectorField& U
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ powerLaw(const powerLaw&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const powerLaw&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("powerLaw");
+
+ //- Constructor
+ powerLaw
+ (
+ const word& name,
+ const word& modelType,
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Destructor
+ virtual ~powerLaw();
+
+
+ // Member Functions
+
+ //- Add resistance
+ virtual void correct(fvVectorMatrix& UEqn) const;
+
+ //- Add resistance
+ virtual void correct
+ (
+ fvVectorMatrix& UEqn,
+ const volScalarField& rho,
+ const volScalarField& mu
+ ) const;
+
+ //- Add resistance
+ virtual void correct
+ (
+ const fvVectorMatrix& UEqn,
+ volTensorField& AU
+ ) const;
+
+
+ // I-O
+
+ //- Write
+ void writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "powerLawTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..0e1d7d2343f9e558013e4fab619203e862ec0876
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLaw::apply
+(
+ scalarField& Udiag,
+ const scalarField& V,
+ const RhoFieldType& rho,
+ const vectorField& U
+) const
+{
+ const scalar C0 = C0_;
+ const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+ forAll(cellZoneIds_, zoneI)
+ {
+ const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+ forAll(cells, i)
+ {
+ const label cellI = cells[i];
+
+ Udiag[cellI] +=
+ V[cellI]*rho[cellI]*C0*pow(magSqr(U[cellI]), C1m1b2);
+ }
+ }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLaw::apply
+(
+ tensorField& AU,
+ const RhoFieldType& rho,
+ const vectorField& U
+) const
+{
+ const scalar C0 = C0_;
+ const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+ forAll(cellZoneIds_, zoneI)
+ {
+ const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
+
+ forAll(cells, i)
+ {
+ const label cellI = cells[i];
+
+ AU[cellI] =
+ AU[cellI] + I*(rho[cellI]*C0*pow(magSqr(U[cellI]), C1m1b2));
+ }
+ }
+}
+
+
+// ************************************************************************* //