Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.C

@@ -81,7 +81,7 @@ bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 {
     PDRDragModelCoeffs_ = PDRProperties.subDict(type() + "Coeffs");
 
-    PDRDragModelCoeffs_.lookup("PDRDragModel") >> on_;
+    PDRDragModelCoeffs_.lookup("drag") >> on_;
 
     return true;
 }

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C

@@ -91,24 +91,80 @@ Foam::PDRDragModels::basic::~basic()
 
 Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
 {
-    const volScalarField& betav =
+    tmp<volSymmTensorField> tDragDcu
-        U_.db().lookupObject<volScalarField>("betav");
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "tDragDcu",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U_.mesh(),
+            dimensionedSymmTensor
+            (
+                "zero",
+                dimMass/dimTime/pow(dimLength, 3),
+                pTraits<symmTensor>::zero
+            ),
+            zeroGradientFvPatchSymmTensorField::typeName
+        )
+    );
 
-    return (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*sqr(Aw_);
+    volSymmTensorField& DragDcu = tDragDcu();
+
+    if (on_)
+    {
+        const volScalarField& betav =
+            U_.db().lookupObject<volScalarField>("betav");
+
+        DragDcu =
+            (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*sqr(Aw_);
+    }
+
+    return tDragDcu;
 }
 
 
 Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
 {
-    const volScalarField& betav =
+    tmp<volScalarField> tGk
-        U_.db().lookupObject<volScalarField>("betav");
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "tGk",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U_.mesh(),
+            dimensionedScalar("zero", dimMass/dimLength/pow(dimTime, 3), 0.0),
+            zeroGradientFvPatchVectorField::typeName
+        )
+    );
 
-    const volSymmTensorField& CT =
+    volScalarField& Gk = tGk();
-        U_.db().lookupObject<volSymmTensorField>("CT");
 
-    return
+    if (on_)
-        (0.5*rho_)*mag(U_)*(U_ & CT & U_)
+    {
-      + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
+        const volScalarField& betav =
+            U_.db().lookupObject<volScalarField>("betav");
+
+        const volSymmTensorField& CT =
+            U_.db().lookupObject<volSymmTensorField>("CT");
+
+        Gk =
+            (0.5*rho_)*mag(U_)*(U_ & CT & U_)
+          + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
+    }
+
+    return tGk;
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C

@@ -50,7 +50,12 @@ Foam::XiEqModels::Gulder::Gulder
 :
     XiEqModel(XiEqProperties, thermo, turbulence, Su),
     XiEqCoef_(readScalar(XiEqModelCoeffs_.lookup("XiEqCoef"))),
-    SuMin_(0.01*Su.average())
+    SuMin_(0.01*Su.average()),
+    uPrimeCoef_(readScalar(XiEqModelCoeffs_.lookup("uPrimeCoef"))),
+    subGridSchelkin_
+    (
+        readBool(XiEqModelCoeffs_.lookup("subGridSchelkin"))
+    )
 {}
 
 
@@ -67,9 +72,9 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
     const volScalarField& epsilon = turbulence_.epsilon();
 
-    if (subGridSchelkin())
+    if (subGridSchelkin_)
     {
-        up.internalField() += calculateSchelkinEffect();
+        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.H

@@ -61,6 +61,12 @@ class Gulder
         //- Minimum laminar burning velocity
         const dimensionedScalar SuMin_;
 
+        //- Schelkin effect Model constant
+        scalar uPrimeCoef_;
+
+        //- Use sub-grid Schelkin effect
+        bool subGridSchelkin_;
+
 
     // Private Member Functions
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C

@@ -53,6 +53,11 @@ Foam::XiEqModels::SCOPEXiEq::SCOPEXiEq
     XiEqExp_(readScalar(XiEqModelCoeffs_.lookup("XiEqExp"))),
     lCoef_(readScalar(XiEqModelCoeffs_.lookup("lCoef"))),
     SuMin_(0.01*Su.average()),
+    uPrimeCoef_(readScalar(XiEqModelCoeffs_.lookup("uPrimeCoef"))),
+    subGridSchelkin_
+    (
+        readBool(XiEqModelCoeffs_.lookup("subGridSchelkin"))
+    ),
     MaModel
     (
         IOdictionary
@@ -84,9 +89,9 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
     const volScalarField& epsilon = turbulence_.epsilon();
 
     volScalarField up(sqrt((2.0/3.0)*k));
-    if (subGridSchelkin())
+    if (subGridSchelkin_)
     {
-        up.internalField() += calculateSchelkinEffect();
+        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
@@ -152,6 +157,8 @@ bool Foam::XiEqModels::SCOPEXiEq::read(const dictionary& XiEqProperties)
     XiEqModelCoeffs_.lookup("XiEqCoef") >> XiEqCoef_;
     XiEqModelCoeffs_.lookup("XiEqExp") >> XiEqExp_;
     XiEqModelCoeffs_.lookup("lCoef") >> lCoef_;
+    XiEqModelCoeffs_.lookup("uPrimeCoef") >> uPrimeCoef_;
+    XiEqModelCoeffs_.lookup("subGridSchelkin") >> subGridSchelkin_;
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.H

@@ -70,6 +70,12 @@ class SCOPEXiEq
         //- Minimum Su
         dimensionedScalar SuMin_;
 
+        //- Schelkin effect Model constant
+        scalar uPrimeCoef_;
+
+        //- Use sub-grid Schelkin effect
+        bool subGridSchelkin_;
+
         //- The SCOPE laminar flame speed model used to obtain the
         //  Marstein number.  Note: the laminar flame speed need not be
         //  obtained form the same model.

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C

@@ -77,11 +77,6 @@ Foam::XiEqModel::XiEqModel
             IOobject::NO_WRITE
         ),
         Su.mesh()
-    ),
-    uPrimeCoef_(XiEqModelCoeffs_.lookupOrDefault<scalar>("uPrimeCoef", 0.0)),
-    subGridSchelkin_
-    (
-        XiEqModelCoeffs_.lookupOrDefault<bool>("subGridSchelkin", false)
     )
 {}
 
@@ -98,10 +93,6 @@ bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 {
     XiEqModelCoeffs_ = XiEqProperties.subDict(type() + "Coeffs");
 
-    uPrimeCoef_ = XiEqModelCoeffs_.lookupOrDefault<scalar>("uPrimeCoef", 0.0);
-    subGridSchelkin_ =
-        XiEqModelCoeffs_.lookupOrDefault<bool>("subGridSchelkin", false);
-
     return true;
 }
 
@@ -119,7 +110,7 @@ void Foam::XiEqModel::writeFields() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::XiEqModel::calculateSchelkinEffect() const
+Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
 {
     const fvMesh& mesh = Su_.mesh();
 
@@ -188,7 +179,7 @@ Foam::XiEqModel::calculateSchelkinEffect() const
 
     const scalarField cellWidth(pow(mesh.V(), 1.0/3.0));
 
-    const scalarField upLocal(uPrimeCoef_*sqrt((U & CT & U)*cellWidth));
+    const scalarField upLocal(uPrimeCoef*sqrt((U & CT & U)*cellWidth));
 
     const scalarField deltaUp(upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0));
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.H

@@ -80,12 +80,6 @@ protected:
         //
         volSymmTensorField nsv_;
 
-        //- Schelkin effect Model constant
-        scalar uPrimeCoef_;
-
-        //- Use sub-grid Schelkin effect
-        bool subGridSchelkin_;
-
 
 private:
 
@@ -162,14 +156,9 @@ public:
             return turbulence_.muEff();
         }
 
-        //- Return state of the sub-grid Schelkin effect
-        bool subGridSchelkin() const
-        {
-            return subGridSchelkin_;
-        }
 
         //- Return the sub-grid Schelkin effect
-        tmp<volScalarField> calculateSchelkinEffect() const;
+        tmp<volScalarField> calculateSchelkinEffect(const scalar) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties) = 0;

src/regionModels/regionModel/regionModel/regionModel.C

@@ -44,19 +44,22 @@ namespace regionModels
 void Foam::regionModels::regionModel::constructMeshObjects()
 {
     // construct region mesh
-    regionMeshPtr_.reset
+    if (!time_.foundObject<fvMesh>(regionName_))
-    (
+    {
-        new fvMesh
+        regionMeshPtr_.reset
         (
-            IOobject
+            new fvMesh
             (
-                lookup("regionName"),
+                IOobject
-                time_.timeName(),
+                (
-                time_,
+                    regionName_,
-                IOobject::MUST_READ
+                    time_.timeName(),
+                    time_,
+                    IOobject::MUST_READ
+                )
             )
-        )
+        );
-    );
+    }
 }
 
 
@@ -66,19 +69,22 @@ void Foam::regionModels::regionModel::constructMeshObjects
 )
 {
     // construct region mesh
-    regionMeshPtr_.reset
+    if (!time_.foundObject<fvMesh>(regionName_))
-    (
+    {
-        new fvMesh
+        regionMeshPtr_.reset
         (
-            IOobject
+            new fvMesh
             (
-                dict.lookup("regionName"),
+                IOobject
-                time_.timeName(),
+                (
-                time_,
+                    regionName_,
-                IOobject::MUST_READ
+                    time_.timeName(),
+                    time_,
+                    IOobject::MUST_READ
+                )
             )
-        )
+        );
-    );
+    }
 }
 
 
@@ -109,11 +115,21 @@ void Foam::regionModels::regionModel::initialise()
 
             nBoundaryFaces += regionPatch.faceCells().size();
 
-            const mappedWallPolyPatch& mapPatch =
+            const mappedPatchBase& mapPatch =
-                refCast<const mappedWallPolyPatch>(regionPatch);
+                refCast<const mappedPatchBase>(regionPatch);
 
-            const label primaryPatchI = mapPatch.samplePolyPatch().index();
+            if
-            primaryPatchIDs.append(primaryPatchI);
+            (
+                primaryMesh_.time().foundObject<polyMesh>
+                (
+                    mapPatch.sampleRegion()
+                )
+            )
+            {
+
+                const label primaryPatchI = mapPatch.samplePolyPatch().index();
+                primaryPatchIDs.append(primaryPatchI);
+            }
         }
     }
 
@@ -197,7 +213,8 @@ Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
     regionMeshPtr_(NULL),
     coeffs_(dictionary::null),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_()
+    intCoupledPatchIDs_(),
+    regionName_("none")
 {}
 
 
@@ -228,7 +245,8 @@ Foam::regionModels::regionModel::regionModel
     regionMeshPtr_(NULL),
     coeffs_(subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_()
+    intCoupledPatchIDs_(),
+    regionName_(lookup("regionName"))
 {
     if (active_)
     {
@@ -273,7 +291,8 @@ Foam::regionModels::regionModel::regionModel
     regionMeshPtr_(NULL),
     coeffs_(dict.subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
-    intCoupledPatchIDs_()
+    intCoupledPatchIDs_(),
+    regionName_(dict.lookup("regionName"))
 {
     if (active_)
     {
@@ -312,13 +331,6 @@ void Foam::regionModels::regionModel::evolve()
 {
     if (active_)
     {
-        if (primaryMesh_.changing())
-        {
-            FatalErrorIn("regionModel::evolve()")
-                << "Currently not possible to apply " << modelName_
-                << " model to moving mesh cases" << nl << abort(FatalError);
-        }
-
         Info<< "\nEvolving " << modelName_ << " for region "
             << regionMesh().name() << endl;
 

src/regionModels/regionModel/regionModel/regionModel.H

@@ -117,6 +117,9 @@ protected:
             //- List of patch IDs internally coupled with the primary region
             labelList intCoupledPatchIDs_;
 
+        //- Region name
+        word regionName_;
+
 
     // Protected member functions
 

src/regionModels/regionModel/regionModel/regionModelI.H

@@ -60,7 +60,11 @@ inline const Foam::word& Foam::regionModels::regionModel::modelName() const
 
 inline const Foam::fvMesh& Foam::regionModels::regionModel::regionMesh() const
 {
-    if (!regionMeshPtr_.valid())
+    if (time_.foundObject<fvMesh>(regionName_))
+    {
+        return time_.lookupObject<fvMesh>(regionName_);
+    }
+    else if (!regionMeshPtr_.valid())
     {
         FatalErrorIn
         (
@@ -75,7 +79,14 @@ inline const Foam::fvMesh& Foam::regionModels::regionModel::regionMesh() const
 
 inline Foam::fvMesh& Foam::regionModels::regionModel::regionMesh()
 {
-    if (!regionMeshPtr_.valid())
+    if (time_.foundObject<fvMesh>(regionName_))
+    {
+        return const_cast<fvMesh&>
+        (
+            time_.lookupObject<fvMesh>(regionName_)
+        );
+    }
+    else if (!regionMeshPtr_.valid())
     {
         FatalErrorIn
         (

src/regionModels/regionModel/regionModel1D/regionModel1D.C

@@ -39,7 +39,6 @@ namespace regionModels
 
 void Foam::regionModels::regionModel1D::constructMeshObjects()
 {
-    const fvMesh& regionMesh = regionMeshPtr_();
 
     nMagSfPtr_.reset
     (
@@ -49,11 +48,11 @@ void Foam::regionModels::regionModel1D::constructMeshObjects()
             (
                 "nMagSf",
                 time().timeName(),
-                regionMesh,
+                regionMesh(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            regionMesh,
+            regionMesh(),
             dimensionedScalar("zero", dimArea, 0.0)
         )
     );
@@ -163,7 +162,7 @@ bool Foam::regionModels::regionModel1D::read(const dictionary& dict)
 {
     if (regionModel::read(dict))
     {
-        moveMesh_.readIfPresent("moveMesh", dict);
+        moveMesh_.readIfPresent("moveMesh", coeffs_);
 
         return true;
     }

src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.C

@@ -27,6 +27,7 @@ License
 #include "thermoBaffle2D.H"
 
 #include "fvm.H"
+#include "fvcDiv.H"
 #include "addToRunTimeSelectionTable.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fvMatrices.H"
@@ -58,6 +59,13 @@ bool thermoBaffle2D::read()
 }
 
 
+bool thermoBaffle2D::read(const dictionary& dict)
+{
+    this->solution().lookup("nNonOrthCorr") >> nNonOrthCorr_;
+    return regionModel1D::read(dict);
+}
+
+
 void thermoBaffle2D::solveEnergy()
 {
     if (debug)
@@ -128,6 +136,19 @@ void thermoBaffle2D::solveEnergy()
         Q
     );
 
+    if (moveMesh_)
+    {
+        surfaceScalarField phiMesh
+        (
+            fvc::interpolate(rhoCp*T_)*regionMesh().phi()
+        );
+
+        TEqn -= fvc::div(phiMesh);
+    }
+
+    Info<< "T gas min/max   = " << min(T_).value() << ", "
+        << max(T_).value() << endl;
+
     TEqn.relax();
     TEqn.solve();
 
@@ -157,7 +178,7 @@ thermoBaffle2D::thermoBaffle2D
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         regionMesh(),
         dimensionedScalar
@@ -175,7 +196,7 @@ thermoBaffle2D::thermoBaffle2D
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         regionMesh(),
         dimensionedScalar
@@ -209,7 +230,7 @@ thermoBaffle2D::thermoBaffle2D
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         regionMesh(),
         dimensionedScalar
@@ -227,7 +248,7 @@ thermoBaffle2D::thermoBaffle2D
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         regionMesh(),
         dimensionedScalar

src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.H

@@ -102,9 +102,13 @@ protected:
 
     // Protected member functions
 
-        //- Read control parameters from dictionary
+        //- Read control parameters IO dictionary
         virtual bool read();
 
+        //- Read control parameters from dictionary
+        virtual bool read(const dictionary& dict);
+
+
         // Equations
 
             //- Solve energy equation

src/regionModels/thermoBaffleModels/thermoBaffleModel/thermoBaffleModel.C

@@ -53,6 +53,13 @@ bool thermoBaffleModel::read()
 }
 
 
+bool thermoBaffleModel::read(const dictionary& dict)
+{
+    regionModel1D::read(dict);
+    return true;
+}
+
+
 void thermoBaffleModel::init()
 {
     if (active_)

src/regionModels/thermoBaffleModels/thermoBaffleModel/thermoBaffleModel.H

@@ -91,9 +91,12 @@ protected:
 
     // Protected Member Functions
 
-        //- Read control parameters from dictionary
+        //- Read control parameters from IO dictionary
         virtual bool read();
 
+        //- Read control parameters from dictionary
+        virtual bool read(const dictionary&);
+
 
 public:
 

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C

@@ -230,7 +230,7 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         else //out
         {
             this->refGrad()[i] = 0.0;
-            this->refValue()[i] = KDelta[i]*q[i] + patchInternalField()()[i];
+            this->refValue()[i] = q[i]/KDelta[i] + patchInternalField()()[i];
             this->valueFraction()[i] = 1.0;
         }
     }

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/PDRProperties

@@ -50,26 +50,26 @@ instabilityCoeffs
             XiEqModelL
             {
                 XiEqModel       Gulder;
-                uPrimeCoef      1.0;
-                subGridSchelkin true;
 
                 GulderCoeffs
                 {
                     XiEqCoef   0.62;
+                    uPrimeCoef      1.0;
+                    subGridSchelkin true;
                 }
             }
 
             XiEqModelH
             {
                 XiEqModel       SCOPEXiEq;
-                uPrimeCoef      1.0;
-                subGridSchelkin true;
 
                 SCOPEXiEqCoeffs
                 {
                     XiEqCoef   1.6;
                     XiEqExp    0.33333;
                     lCoef      0.336;
+                    uPrimeCoef      1.0;
+                    subGridSchelkin true;
                 }
             }
         }

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes

@@ -60,7 +60,7 @@ laplacianSchemes
     laplacian(muEff,U) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
-    laplacian((rho*inv((((1)*A(U))+((((0.5*rho)*CR)*mag(U))+((((Csu*(1))*betav)*muEff)*sqr(Aw)))))),p) Gauss linear limited 0.333;
+    laplacian((rho*inv((((1)*A(U))+tDragDcu))),p) Gauss linear limited 0.333;
     laplacian(Db,b) Gauss linear limited 0.333;
     laplacian(Db,ft) Gauss linear limited 0.333;
     laplacian(Db,h) Gauss linear limited 0.333;