diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C new file mode 100644 index 0000000000..73ca510c59 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C @@ -0,0 +1,229 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "normBasicXiSubXiEq.H" +#include "addToRunTimeSelectionTable.H" +#include "ignition.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + defineTypeNameAndDebug(normBasicSubGrid, 0); + addToRunTimeSelectionTable(XiEqModel, normBasicSubGrid, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiEqModels::normBasicSubGrid::normBasicSubGrid +( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su +) +: + XiEqModel(XiEqProperties, modelType,thermo, turbulence, Su), + Cxpe1_(XiEqModelCoeffs_.get("Cxpe1")), + Cxpe2_(XiEqModelCoeffs_.get("Cxpe2")), + Cxpe3_(XiEqModelCoeffs_.get("Cxpe3")), + Cxpe4_(XiEqModelCoeffs_.get("Cxpe4")) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiEqModels::normBasicSubGrid::~normBasicSubGrid() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::XiEqModels::normBasicSubGrid::XiEq() const +{ + const fvMesh& mesh = Su_.mesh(); + const volVectorField& U = mesh.lookupObject("U"); + const volScalarField& b = mesh.lookupObject("b"); + const volScalarField& Nv = mesh.lookupObject("Nv"); + const volSymmTensorField& nsv = + mesh.lookupObject("nsv"); + const volSymmTensorField& Bv = + mesh.lookupObject("Bv"); + + volScalarField magU(mag(U)); + + const scalarField Cw = pow(mesh.V(), 2.0/3.0); + + tmp tN + ( + new volScalarField + ( + IOobject + ( + "tN", + mesh.time().constant(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar(Nv.dimensions(), Zero) + ) + ); + + volScalarField& N = tN.ref(); + + N.primitiveFieldRef() = Nv.primitiveField()*Cw; + + tmp tns + ( + new volSymmTensorField + ( + IOobject + ( + "tns", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedSymmTensor(nsv.dimensions(), Zero) + ) + ); + volSymmTensorField& ns = tns.ref(); + + tmp tB + ( + new volSymmTensorField + ( + IOobject + ( + "tB", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedSymmTensor(Bv.dimensions(), Zero) + ) + ); + volSymmTensorField& B = tB.ref(); + + //calculating flame normal + + volVectorField flNormal + ( + "flNormal", + fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()) + ); + + volScalarField mgb("mgb", mag(flNormal)); + + dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL); + + const volScalarField bc(b*(1.0-b)); + + dMgb += 1.0e-8* + (bc*mgb)().weightedAverage(mesh.V()) + /(bc.weightedAverage(mesh.V()) + SMALL); + + mgb += dMgb; + flNormal /= mgb; + + B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw); + volScalarField Ntemp("Ntemp", N); + volScalarField Np("Np", max(N - (flNormal & ns & flNormal), scalar(1))); + + // B_ is Bv*sqrt(Cw) + volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np)); + bl.min(1.0); + + volScalarField up(sqrt((2.0/3.0)*turbulence_.k())); + + IOdictionary combustionProperties + ( + IOobject + ( + "combustionProperties", + mesh.time().constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + + ignition ign(combustionProperties, mesh.time(), mesh); + + dimensionedVector ignLoc("ignLoc", dimLength, ign.sites()[0].location()); + + dimensionedScalar filtRad2 + ( + "filtRad2", + dimLength, + 6.0*ign.sites()[0].diameter() + ); + + const volScalarField filDist(mag(mesh.C() - ignLoc)); + + const volScalarField filterMult + ( + pos(filDist - filtRad2)*neg(bl - 0.99)*pos(N - 1e-3) + ); + + tmp XiSubEq + ( + scalar(1) + + min( min(Cxpe1_, Cxpe2_*magU/up)*sqrt(bl), Cxpe3_) + * filterMult + ); + + return XiSubEq; +} + + +bool Foam::XiEqModels::normBasicSubGrid::read(const dictionary& XiEqProperties) +{ + XiEqModel::read(XiEqProperties); + + XiEqModelCoeffs_.readEntry("Cxpe1", Cxpe1_); + XiEqModelCoeffs_.readEntry("Cxpe2", Cxpe2_); + XiEqModelCoeffs_.readEntry("Cxpe3", Cxpe3_); + XiEqModelCoeffs_.readEntry("Cxpe4", Cxpe4_); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.H b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.H new file mode 100644 index 0000000000..b17d2401d5 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.H @@ -0,0 +1,117 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiEqModels::normBasicSubGrid + +Description + + +SourceFiles + normBasicSubGrid.C + +\*---------------------------------------------------------------------------*/ + +#ifndef normBasicSubGrid_H +#define normBasicSubGrid_H + +#include "XiEqModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + +/*---------------------------------------------------------------------------*\ + Class normBasicSubGrid Declaration +\*---------------------------------------------------------------------------*/ + +class normBasicSubGrid +: + public XiEqModel +{ + // Private data + + // Constants in the equilibrium Xp equation + scalar Cxpe1_; + scalar Cxpe2_; + scalar Cxpe3_; + scalar Cxpe4_; + + + // Private Member Functions + + //- Disallow copy construct + normBasicSubGrid(const normBasicSubGrid&); + + //- Disallow default bitwise assignment + void operator=(const normBasicSubGrid&); + + +public: + + //- Runtime type information + TypeName("normBasicSubGrid"); + + + // Constructors + + //- Construct from components + normBasicSubGrid + ( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su + ); + + + //- Destructor + virtual ~normBasicSubGrid(); + + + // Member Functions + + //- Return the flame-wrinking XiEq + virtual tmp XiEq() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiEqProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiEqModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C new file mode 100644 index 0000000000..a5eddfa34b --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C @@ -0,0 +1,300 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "normBasicXiSubG.H" +#include "addToRunTimeSelectionTable.H" +#include "zeroGradientFvPatchField.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiGModels +{ + defineTypeNameAndDebug(normBasicSubGrid, 0); + addToRunTimeSelectionTable(XiGModel, normBasicSubGrid, dictionary); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiGModels::normBasicSubGrid::normBasicSubGrid +( + const dictionary& XiGProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su +) +: + XiGModel(XiGProperties, modelType, thermo, turbulence, Su), +// Bv_ +// ( +// IOobject +// ( +// "Bv", +// Su.mesh().facesInstance(), +// Su.mesh(), +// IOobject::MUST_READ, +// IOobject::NO_WRITE +// ), +// Su.mesh() +// ), + k1_(XiGModelCoeffs_.get("k1")), + kb1_(XiGModelCoeffs_.get("kb1")), + kbe_(XiGModelCoeffs_.get("kbe")), + kbx_(XiGModelCoeffs_.get("kbx")), + k2_(XiGModelCoeffs_.get("k2")), + LOverCw_(XiGModelCoeffs_.get("LOverCw")) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiGModels::normBasicSubGrid::~normBasicSubGrid() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::XiGModels::normBasicSubGrid::G() const +{ + const objectRegistry& db = Su_.db(); + const fvMesh& mesh = Su_.mesh(); + const volVectorField& U = db.lookupObject("U"); + const volScalarField& b = db.lookupObject("b"); + const volScalarField& Nv = db.lookupObject("Nv"); + const volScalarField& St = db.lookupObject("St"); + const volSymmTensorField& nsv = db.lookupObject("nsv"); + const volScalarField& Lobs = db.lookupObject("Lobs"); + const volSymmTensorField& Bv = db.lookupObject("Bv"); + + const scalarField Cw(pow(Su_.mesh().V(), 2.0/3.0)); + volScalarField CwVol + ( + IOobject + ( + "CwVol", + mesh.time().timeName(), + mesh + ), + mesh, + dimensionSet(dimLength), + Cw, + zeroGradientFvPatchField::typeName + ); + CwVol.correctBoundaryConditions(); + + if (!db.foundObject("Ep")) + { + FatalErrorIn + ( + "Foam::tmp Foam::XiGModels::" + "normBasicSubGrid::G() const" + ) + << "Looking for Ep in db which does not exist " + << Foam::abort(FatalError); + } + + const volScalarField& Ep = db.lookupObject("Ep"); + const volScalarField& Xp = db.lookupObject("Xp"); + const volScalarField& Xi = db.lookupObject("Xi"); + + //tmp tGtot = XiGModel_->G(); + tmp tGtot + ( + new volScalarField + ( + IOobject + ( + "tGtot", + Su_.mesh().time().timeName(), + Su_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + Su_.mesh(), + dimensionedScalar(inv(dimTime), Zero) + + ) + ); + + volScalarField& Gtot = tGtot.ref(); + + + //calculating flame normal + + volVectorField flNormal(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf())); + volScalarField mgb("mgb", mag(flNormal)); + dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL); + { + volScalarField bc(b*(1.0-b)); + + dMgb += 1.0e-8* + (bc*mgb)().weightedAverage(mesh.V()) + /(bc.weightedAverage(mesh.V()) + SMALL); + } + mgb += dMgb; + flNormal /= mgb; + + + tmp tN + ( + new volScalarField + ( + IOobject + ( + "tN", + Su_.mesh().time().timeName(), + Su_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + Su_.mesh(), + dimensionedScalar(Nv.dimensions(), Zero) + + ) + ); + + volScalarField& N = tN.ref(); + + N.primitiveFieldRef() = Nv.primitiveField()*Cw; + + tmp tns + ( + new volSymmTensorField + ( + IOobject + ( + "tns", + Su_.mesh().time().timeName(), + Su_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + Su_.mesh(), + dimensionedSymmTensor(nsv.dimensions(), Zero) + ) + ); + + volSymmTensorField& ns = tns.ref(); + + ns.primitiveFieldRef() = nsv.primitiveField()*Cw; + + tmp tB + ( + new volSymmTensorField + ( + IOobject + ( + "tB", + Su_.mesh().time().timeName(), + Su_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + Su_.mesh(), + dimensionedSymmTensor(Bv.dimensions(), Zero) + ) + ); + + volSymmTensorField& B = tB.ref(); + + B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw); + + volScalarField Np(max(N - (flNormal & ns & flNormal), scalar(1))); + + // B_ is Bv*sqrt(Cw) + volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np)); + bl.min(1.0); + + volScalarField flSpeed("flSpeed", ((U & flNormal) + St)*b/(b+SMALL)) ; + + volScalarField up("up", sqrt((2.0/3.0)*turbulence_.k())); + + const volScalarField Gtot1 + ( + "Gtot1", + ( + k1_ + kb1_*min(pow(bl, kbe_), kbx_) + )*mag(flSpeed)/(max(Lobs, LOverCw_*CwVol)) + ); + + const volScalarField Gtot2("Gtot2", k2_*Ep*Su_*Xi/(Xp - 0.999)); + + const volScalarField value(pos(N - 1.e-3)*neg(bl - 0.99)); + + Gtot = value*Gtot1+(1.0 - value)*Gtot2; + + + //if (Xi.mesh().time().outputTime()) + { + //Gtot.write(); + //bl.write(); + //Lobs.write(); + //flSpeed.write(); + //N.write(); + } + + return tGtot; +} + + +Foam::tmp Foam::XiGModels::normBasicSubGrid::Db() const +{ + // Not used // + const objectRegistry& db = Su_.db(); + const volScalarField& Xi = db.lookupObject("Xi"); + const volScalarField& rho = db.lookupObject("rho"); + const volScalarField& mgb = db.lookupObject("mgb"); + const volScalarField& Lobs = db.lookupObject("Lobs"); + const volScalarField& Db = db.lookupObject("Db"); + + //return turbulence_.muEff() + return Db + + rho*Su_*(Xi - 1.0)*mgb*(0.5*Lobs)*Lobs/(mgb*Lobs + 1.0); +} + + +bool Foam::XiGModels::normBasicSubGrid::read(const dictionary& XiGProperties) +{ + XiGModel::read(XiGProperties); + + XiGModelCoeffs_.readEntry("k1", k1_); + XiGModelCoeffs_.readEntry("kb1", kb1_); + XiGModelCoeffs_.readEntry("kbe", kbe_); + XiGModelCoeffs_.readEntry("kbx", kbx_); + XiGModelCoeffs_.readEntry("k2", k2_); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.H b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.H new file mode 100644 index 0000000000..080203f55b --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.H @@ -0,0 +1,132 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiGModel::normBasicSubGrid + +Description + + + +SourceFiles + normBasicSubGrid.C + +\*---------------------------------------------------------------------------*/ + +#ifndef normBasicSubGrid_H +#define normBasicSubGrid_H + +#include "XiGModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiGModels +{ + +/*---------------------------------------------------------------------------*\ + Class normBasicSubGrid Declaration +\*---------------------------------------------------------------------------*/ + +class normBasicSubGrid +: + public XiGModel +{ + // Private data + + //- Sub-grid generation rate coefficient + scalar k1_; + + //- Sub-grid generation rate coefficient - * sqrt(b) + scalar kb1_; + + //- Sub-grid generation rate coefficient - * b + scalar kbe_; + + //- Sub-grid generation rate upper limit coefficient - * b + scalar kbx_; + + //- Sub-grid generation rate coefficient + scalar k2_; + + //- Maximum Lobs/CellWidth + scalar LOverCw_; + + // Private Member Functions + + //- Disallow copy construct + normBasicSubGrid(const normBasicSubGrid&); + + //- Disallow default bitwise assignment + void operator=(const normBasicSubGrid&); + + +public: + + //- Runtime type information + TypeName("normBasicSubGridG"); + + + // Constructors + + //- Construct from components + normBasicSubGrid + ( + const dictionary& XiGProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su + ); + + + //- Destructor + virtual ~normBasicSubGrid(); + + + // Member Functions + + //- Return the flame-wrinking generation rate + virtual tmp G() const; + + //- Return the flame diffusivity + virtual tmp Db() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiGProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiGModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.C b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.C new file mode 100644 index 0000000000..39c86f748c --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.C @@ -0,0 +1,319 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "basicSch.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace PDRDragModels +{ + defineTypeNameAndDebug(basicSch, 0); + addToRunTimeSelectionTable(PDRDragModel, basicSch, dictionary); +} +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::PDRDragModels::basicSch::basicSch +( + const dictionary& PDRProperties, + const compressible::RASModel& turbulence, + const volScalarField& rho, + const volVectorField& U, + const surfaceScalarField& phi +) +: + PDRDragModel(PDRProperties,turbulence, rho, U, phi), + Csu("Csu", dimless, PDRDragModelCoeffs_), + Csk("Csk", dimless, PDRDragModelCoeffs_), + Aw_ + ( + IOobject + ( + "Aw", + U_.mesh().facesInstance(), + U_.mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + U_.mesh() + ), + + CR_ + ( + IOobject + ( + "CR", + U_.mesh().facesInstance(), + U_.mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + U_.mesh() + ), + nrCoef_(PDRDragModelCoeffs_.get("nrCoef")), + nrExp2_(PDRDragModelCoeffs_.get("nrExp2")), + lCoef_(PDRDragModelCoeffs_.get("lCoef")), + maxSchFac_(PDRDragModelCoeffs_.get("maxSchFac")), + subGridSchelkin_(PDRDragModelCoeffs_.get("subGridSchelkin")) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::PDRDragModels::basicSch::~basicSch() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::PDRDragModels::basicSch::Dcu() const +{ + tmp tDragDcu + ( + new volSymmTensorField + ( + IOobject + ( + "tDragDcu", + U_.mesh().time().constant(), + U_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + U_.mesh(), + dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero) + + ) + ); + + volSymmTensorField& DragDcu = tDragDcu.ref(); + + if (on_) + { + const volScalarField& betav = + U_.db().lookupObject("betav"); + + DragDcu = + (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*sqr(Aw_); + } + + return tDragDcu; +} + + +Foam::tmp Foam::PDRDragModels::basicSch::Gk() const +{ + tmp tGk + ( + new volScalarField + ( + IOobject + ( + "tGk", + U_.mesh().time().constant(), + U_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + U_.mesh(), + dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero) + + ) + ); + + volScalarField& Gk = tGk.ref(); + + if (on_) + { + const volScalarField& betav = + U_.db().lookupObject("betav"); + + const volSymmTensorField& CT = + U_.db().lookupObject("CT"); + + Gk = + (0.5*rho_)*mag(U_)*(U_ & CT & U_) + + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_); + + if (subGridSchelkin_) + { + Gk *= schFac(); + } + } + + + return tGk; +} + +Foam::tmp Foam::PDRDragModels::basicSch::schFac() const +{ + const volScalarField& Su_ = U_.db().lookupObject("Su"); + const volScalarField& rhou_ = U_.db().lookupObject("rhou"); + const volScalarField& muu_ = U_.db().lookupObject("muu"); + + tmp tfac + ( + new volScalarField + ( + IOobject + ( + "tfac", + U_.mesh().time().constant(), + U_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + U_.mesh(), + dimensionedScalar(dimless, Zero) + + ) + ); + + volScalarField& schFac = tfac.ref(); + + const volScalarField& k = turbulence_.k(); + const volScalarField& epsilon = turbulence_.epsilon(); + + const volScalarField up(sqrt((2.0/3.0)*k)); + const volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon); + + volScalarField Rs(Su_*l*rhou_/muu_); + + if (subGridSchelkin_) + { + schFac = max + ( + 1.0, + min + ( + maxSchFac_, + pow(Rs, 2.0 * SchelkinExponent(nrCoef_, nrExp2_, Su_)) + ) + ); + } + + return tfac; +} + + +Foam::tmp Foam::PDRDragModels::basicSch::SchelkinExponent +( + const scalar nrCoef, + const scalar nrExp, + const volScalarField& Su +) const +{ + const fvMesh& mesh = Su.mesh(); + + const volVectorField& U = mesh.lookupObject("U"); + + const volScalarField& Nv = mesh.lookupObject("Nv"); + const volSymmTensorField& nsv = + mesh.lookupObject("nsv"); + + tmp tN + ( + new volScalarField + ( + IOobject + ( + "tN", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedScalar(Nv.dimensions(), Zero) + + ) + ); + + volScalarField& N = tN.ref(); + + N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0); + + tmp tns + ( + new volSymmTensorField + ( + IOobject + ( + "tns", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedSymmTensor(nsv.dimensions(), Zero) + ) + ); + + volSymmTensorField& ns = tns.ref(); + + ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0); + + const volVectorField Uhat + ( + U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4)) + ); + + const volScalarField nr(sqrt(max(N - (Uhat & ns & Uhat), scalar(1.0)))); + + //Re use tN + N.primitiveFieldRef() = + nrCoef*((scalar(1.0) - pow(nrExp, nr))/(1.0 - nrExp) - scalar(1.0)); + + return tN; + +} + + +bool Foam::PDRDragModels::basicSch::read(const dictionary& PDRProperties) +{ + PDRDragModel::read(PDRProperties); + + PDRDragModelCoeffs_.readEntry("Csu", Csu.value()); + PDRDragModelCoeffs_.readEntry("Csk", Csk.value()); + + return true; +} + + +void Foam::PDRDragModels::basicSch::writeFields() const +{ + Aw_.write(); + CR_.write(); +} + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.H b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.H new file mode 100644 index 0000000000..0a129b6df5 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.H @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::PDRDragModels::basicSch + +Description + + +SourceFiles + basicSch.C + +\*---------------------------------------------------------------------------*/ + +#ifndef basicSch_H +#define basicSch_H + +#include "PDRDragModel.H" +#include "XiEqModel.H" + + +namespace Foam +{ +namespace PDRDragModels +{ + +/*---------------------------------------------------------------------------*\ + Class basicSch Declaration +\*---------------------------------------------------------------------------*/ + +class basicSch +: + public PDRDragModel +{ + // Private data + + dimensionedScalar Csu; + dimensionedScalar Csk; + volScalarField Aw_; + volSymmTensorField CR_; + + //- Schelkin effect Model constants + const scalar nrCoef_; + const scalar nrExp2_; + const scalar lCoef_; + const scalar maxSchFac_; + + //- Use sub-grid Schelkin effect + bool subGridSchelkin_; + + + // Private Member Functions + + //- Disallow copy construct + basicSch(const basicSch&); + + //- Disallow default bitwise assignment + void operator=(const basicSch&); + + +public: + + //- Runtime type information + TypeName("basicSch"); + + + // Constructors + + //- Construct from components + basicSch + ( + const dictionary& PDRProperties, + const compressible::RASModel& turbulence, + const volScalarField& rho, + const volVectorField& U, + const surfaceScalarField& phi + ); + + + //- Destructor + virtual ~basicSch(); + + + // Member Functions + + //- Return the momentum drag coefficient + virtual tmp Dcu() const; + + //- Return the momentum drag turbulence generation rate + virtual tmp Gk() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& PDRProperties); + + //- Write fields + void writeFields() const; + + //- Return the schelkin factor for drag turbulence generation rate + tmp schFac() const; + + //- Return the sub-grid Schelkin effect exponent + tmp SchelkinExponent + ( + const scalar, + const scalar, + const volScalarField& + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace PDRDragModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C new file mode 100644 index 0000000000..2d68189c4e --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C @@ -0,0 +1,261 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "BLMgMaXiEq.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + defineTypeNameAndDebug(BLMgMaXiEq, 0); + addToRunTimeSelectionTable(XiEqModel, BLMgMaXiEq, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiEqModels::BLMgMaXiEq::BLMgMaXiEq +( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su +) +: + XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su), + kaCoef_(XiEqModelCoeffs_.get("kaCoef")), + lowK0_(XiEqModelCoeffs_.get("lowK0")), + lowKg_(XiEqModelCoeffs_.get("lowKg")), + XiEqCoef_(XiEqModelCoeffs_.get("XiEqCoef")), + alphaCoefP_(XiEqModelCoeffs_.get("alphaCoefP")), + betaCoefP_(XiEqModelCoeffs_.get("betaCoefP")), + alphaCoefN_(XiEqModelCoeffs_.get("alphaCoefN")), + betaCoefN_(XiEqModelCoeffs_.get("betaCoefN")), + maLim_(XiEqModelCoeffs_.get("maLim")), + maLim1_(XiEqModelCoeffs_.get("maLim1")), + quenchCoef_(XiEqModelCoeffs_.get("quenchCoef")), + quenchExp_(XiEqModelCoeffs_.get("quenchExp")), + quenchM_(XiEqModelCoeffs_.get("quenchM")), + quenchRate1_(XiEqModelCoeffs_.get("quenchRate1")), + quenchRate2_(XiEqModelCoeffs_.get("quenchRate2")), + lCoef_(XiEqModelCoeffs_.get("lCoef")), + SuMin_(0.01*Su.average()), + uPrimeCoef_(XiEqModelCoeffs_.get("uPrimeCoef")), + nrExp_(XiEqModelCoeffs_.get("nrExp")), + subGridSchelkin_(XiEqModelCoeffs_.get("subGridSchelkin")), + MaModel + ( + IOdictionary + ( + IOobject + ( + "combustionProperties", + Su.mesh().time().constant(), + Su.mesh(), + IOobject::MUST_READ + ) + ), + thermo + ) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiEqModels::BLMgMaXiEq::~BLMgMaXiEq() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::XiEqModels::BLMgMaXiEq::XiEq() const +{ + const volScalarField& k = turbulence_.k(); + const volScalarField& epsilon = turbulence_.epsilon(); + volScalarField up("up", sqrt((2.0/3.0)*k)); + if (subGridSchelkin_) + { + up.primitiveFieldRef() += + calculateSchelkinEffect(uPrimeCoef_, nrExp_); + } + + volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon); + volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu()); + + volScalarField upBySu("upBySu", up/(Su_ + SuMin_)); + volScalarField K("K", kaCoef_*upBySu*upBySu/sqrt(Rl)); + volScalarField Ma("Ma", MaModel.Ma()); + + volScalarField regime("regime", MaModel.Ma()*scalar(0.0)); + + tmp tXiEq + ( + new volScalarField + ( + IOobject + ( + "XiEq", + epsilon.time().timeName(), + epsilon.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + epsilon.mesh(), + dimensionedScalar(dimless, Zero) + ) + ); + + const objectRegistry& db = Su_.db(); + + const volScalarField& b = db.lookupObject("b"); + + const volScalarField multiP1(0.0*pos(b - 0.99) + 1.0*neg(b - 0.99)); + const volScalarField multiP2(1.0*pos(b - 0.01) + 0.0*neg(b - 0.01)); + + volScalarField& xieq = tXiEq.ref(); + + forAll(xieq, celli) + { + scalar alpha; + scalar beta; + scalar gulderMa; + + if (Ma[celli]>= 0) + { + gulderMa = + 1.0 + + (0.1402 - 0.007*Ma[celli]) + * upBySu[celli]*sqrt(upBySu[celli]/K[celli]); + + regime[celli] = multiP1[celli]*multiP2[celli]; + + } + else + { + gulderMa = + 1.0 + + (0.005*Ma[celli]*Ma[celli]+0.01*Ma[celli] + 0.125) + * upBySu[celli]*sqrt(upBySu[celli]/K[celli]); + + regime[celli] = 2*multiP1[celli]*multiP2[celli]; + } + + + if (K[celli] < (lowK0_ + lowKg_*Ma[celli]) ) + { + xieq[celli] = gulderMa; + } + else + { + if (Ma[celli] >= 0.0) + { + alpha = alphaCoefP_*(maLim_ - Ma[celli]); + beta = betaCoefP_*(maLim_ - Ma[celli]); + regime[celli] = 3*multiP1[celli]*multiP2[celli]; + } + else + { + alpha = alphaCoefN_*(maLim1_ - Ma[celli]) ; + beta = betaCoefN_*(maLim_ + Ma[celli]); + regime[celli] = 4*multiP1[celli]*multiP2[celli]; + } + xieq[celli] = XiEqCoef_*alpha*pow(K[celli], beta)*upBySu[celli]; + } + + if (Ma[celli] > -3.0 && Ma[celli] < 11.0) + { + scalar K0p8 = quenchCoef_*pow( Ma[celli] - quenchM_, quenchExp_); + scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[celli]; + if (K[celli] > (K0p8 - 0.223/quenchRate)) + { + xieq[celli] *= 0.8*exp(-quenchRate*(K[celli] - K0p8)); + regime[celli] = 5*multiP1[celli]*multiP2[celli]; + } + } + } + + forAll(xieq.boundaryField(), patchi) + { + scalarField& xieqp = xieq.boundaryFieldRef()[patchi]; + const scalarField& Kp = K.boundaryField()[patchi]; + const scalarField& Map = Ma.boundaryField()[patchi]; + const scalarField& upBySup = upBySu.boundaryField()[patchi]; + + forAll(xieqp, facei) + { + scalar alpha; + scalar beta; + + if (Map[facei] > 0.0) + { + alpha = alphaCoefP_*(maLim_ - Map[facei]); + beta = betaCoefP_*(maLim_ - Map[facei]); + } + else + { + alpha = alphaCoefN_*(maLim_ - Map[facei]); + beta = betaCoefN_*(maLim_ + Map[facei]); + } + xieqp[facei] = + XiEqCoef_*alpha*pow(Kp[facei], beta)*upBySup[facei]; + + if (Map[facei] > -3.0 && Map[facei] < 11.0) + { + scalar K0p8 = quenchCoef_*pow(Map[facei] - quenchM_, quenchExp_); + scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[facei]; + + if (Kp[facei] > (K0p8 - 0.223/quenchRate)) + { + xieqp[facei] *= 0.8*exp(-quenchRate*(Kp[facei] - K0p8)); + } + } + else + { + Info<< + "Markstein Number out of range for Quench Formulation" << endl; + } + } + } + + return tXiEq; +} + + +bool Foam::XiEqModels::BLMgMaXiEq::read(const dictionary& XiEqProperties) +{ + XiEqModel::read(XiEqProperties); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.H b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.H new file mode 100644 index 0000000000..04adbeb649 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.H @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiEqModel::BLMgMaXiEq + +Description + Model for XiEq based on Bradley, Lawes and Mansour (2011) + Cobustion and Falme, 158, 123 correlation + with a linear correction function to give a plausible profile for XiEq. + See \link SCOPELaminarFlameSpeed.H \endlink for details on the SCOPE laminar + flame speed model. + +SourceFiles + BLMgMaXiEq.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BLMgMaXiEq_H +#define BLMgMaXiEq_H + +#include "XiEqModel.H" +#include "SCOPELaminarFlameSpeed.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + +/*---------------------------------------------------------------------------*\ + Class BLMgMaXiEq Declaration +\*---------------------------------------------------------------------------*/ + +class BLMgMaXiEq +: + public XiEqModel +{ + // Private data + + // Model constants + + scalar kaCoef_; + scalar lowK0_; + scalar lowKg_; + scalar XiEqCoef_; + scalar alphaCoefP_; + scalar betaCoefP_; + scalar alphaCoefN_; + scalar betaCoefN_; + scalar maLim_; + scalar maLim1_; + scalar quenchCoef_, quenchExp_, quenchM_; + scalar quenchRate1_, quenchRate2_; + scalar lCoef_; + + //- Minimum Su + dimensionedScalar SuMin_; + + //- Schelkin effect Model constants + scalar uPrimeCoef_; + scalar nrExp_; + + //- 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. + laminarFlameSpeedModels::SCOPE MaModel; + + + // Private Member Functions + + //- Disallow copy construct + BLMgMaXiEq(const BLMgMaXiEq&); + + //- Disallow default bitwise assignment + void operator=(const BLMgMaXiEq&); + + +public: + + //- Runtime type information + TypeName("BLMgMaXiEq"); + + + // Constructors + + //- Construct from components + BLMgMaXiEq + ( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su + ); + + + //- Destructor + virtual ~BLMgMaXiEq(); + + + // Member Functions + + //- Return the flame-wrinking XiEq + virtual tmp XiEq() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiEqProperties); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiEqModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.C new file mode 100644 index 0000000000..a94e085ee9 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.C @@ -0,0 +1,141 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "instability2XiEq.H" +#include "addToRunTimeSelectionTable.H" +#include "IFstream.H" +#include "fvCFD.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + defineTypeNameAndDebug(instability2XiEq, 0); + addToRunTimeSelectionTable(XiEqModel, instability2XiEq, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiEqModels::instability2XiEq::instability2XiEq +( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su +) +: + XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su), + saModel_ + ( + IOdictionary + ( + IOobject + ( + "combustionProperties", + Su.mesh().time().constant(), + Su.mesh(), + IOobject::MUST_READ + ) + ), + thermo + ), + CIn_(saModel_.CIn()), + defaultCIn_(XiEqModelCoeffs_.get("defaultCIn")), + XiEqInFade_(XiEqModelCoeffs_.get("XiEqInFade")), + XiEqModel_ + ( + XiEqModel::New(XiEqModelCoeffs_, modelType, thermo, turbulence, Su) + ) +{ + if (CIn_ <= 0.0) + { + CIn_ = defaultCIn_; + } +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiEqModels::instability2XiEq::~instability2XiEq() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::XiEqModels::instability2XiEq::XiEq() const +{ + IOdictionary combustionProperties + ( + IOobject + ( + "combustionProperties", + Su_.mesh().time().constant(), + Su_.mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh()); + + //const scalar ignTim = ign.sites()[0].tmIgn(); + scalar curTime = Su_.mesh().time().value(); + scalar deltaT = Su_.mesh().time().deltaTValue(); + const scalar ignTim = curTime - deltaT - ign.sites()[0].time(); + + volScalarField turbXiEq(XiEqModel_->XiEq()); + + volScalarField XiEqIn1("XiEqIn1", 0.0*turbXiEq); + + dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_); + dimensionedScalar ignTm("ignTm", dimTime, ignTim); + XiEqIn1 = exp(CIn*Su_*Su_*ignTm) - 1.0; + + return + ( + 1.0 + sqrt(XiEqInFade_*sqr(XiEqIn1) + sqr(turbXiEq - 1.0)) + ); +} + + +bool Foam::XiEqModels::instability2XiEq::read(const dictionary& XiEqProperties) +{ + XiEqModel::read(XiEqProperties); + + XiEqModelCoeffs_.readEntry("defaultCIn", defaultCIn_); + XiEqModelCoeffs_.readEntry("XiEqInFade", XiEqInFade_); + + return XiEqModel_->read(XiEqModelCoeffs_); +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.H b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.H new file mode 100644 index 0000000000..faba197089 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.H @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiEqModels::instability2XiEq + +Description + + +SourceFiles + instability2XiEq.C + +\*---------------------------------------------------------------------------*/ + +#ifndef instability2XiEq_H +#define instability2XiEq_H + +#include "laminarFlameSpeed.H" +#include "SCOPELaminarFlameSpeed.H" +#include "ignitionSite.H" +#include "ignition.H" +#include "Time.H" +#include "fvMesh.H" +#include "XiEqModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiEqModels +{ + +/*---------------------------------------------------------------------------*\ + Class instability2 Declaration +\*---------------------------------------------------------------------------*/ + +class instability2XiEq +: + public XiEqModel +{ + // Private data + + // -Laminar burning speed + laminarFlameSpeedModels::SCOPE saModel_; + + //- GIn (initial instability G)divided by Su^2. Read from fuel file + scalar CIn_; + + //- Default CIn if not in fuel file + scalar defaultCIn_; + + //- Determines how fast XiEqIn fades out as turbulence comes in + scalar XiEqInFade_; + + //- Equilibrium Xi model due to all other effects + autoPtr XiEqModel_; + + + // Private Member Functions + + //- Disallow copy construct + instability2XiEq(const instability2XiEq&); + + //- Disallow default bitwise assignment + void operator=(const instability2XiEq&); + + +public: + + //- Runtime type information + TypeName("instability2XiEq"); + + + // Constructors + + //- Construct from components + instability2XiEq + ( + const dictionary& XiEqProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su + ); + + + //- Destructor + virtual ~instability2XiEq(); + + + // Member Functions + + //- Return the flame-wrinking XiEq + virtual tmp XiEq() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiEqProperties); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiEqModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.C b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.C new file mode 100644 index 0000000000..63c5ea7d3d --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.C @@ -0,0 +1,171 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "IFstream.H" +#include "instability2G.H" +#include "addToRunTimeSelectionTable.H" +#include "fvCFD.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiGModels +{ + defineTypeNameAndDebug(instability2G, 0); + addToRunTimeSelectionTable(XiGModel, instability2G, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiGModels::instability2G::instability2G +( + const dictionary& XiGProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su +) +: + XiGModel(XiGProperties, modelType, thermo, turbulence, Su), + saModel_ + ( + IOdictionary + ( + IOobject + ( + "combustionProperties", + Su.mesh().time().constant(), + Su.mesh(), + IOobject::MUST_READ + ) + ), + thermo + ), + CIn_(saModel_.CIn()), + defaultCIn_(XiGModelCoeffs_.get("defaultCIn")), + GInFade_(XiGModelCoeffs_.get("GInFade")), + GInMult_(XiGModelCoeffs_.get("GInMult")), + lambdaIn_("lambdaIn", XiGModelCoeffs_), + XiGModel_ + ( + XiGModel::New(XiGModelCoeffs_,modelType,thermo, turbulence, Su) + ) +{ + if (CIn_ <= 0.0) + { + CIn_ = defaultCIn_; + } +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiGModels::instability2G::~instability2G() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::XiGModels::instability2G::G() const +{ + IOdictionary combustionProperties + ( + IOobject + ( + "combustionProperties", + Su_.mesh().time().constant(), + Su_.mesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh()); + + scalar curTime = Su_.mesh().time().value(); + scalar deltaT = Su_.mesh().time().deltaTValue(); + const scalar ignTim = curTime - deltaT - ign.sites()[0].time(); + + volScalarField turbXiG(XiGModel_->G()); + + volScalarField GIn("GIn", 0.0*turbXiG); + + + forAll (GIn, i) + { + GIn[i] = CIn_*Su_[i]*Su_[i]*exp(CIn_*Su_[i]*Su_[i]*ignTim)*GInMult_; + } + + dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_); + dimensionedScalar ignTm("ignTm", dimTime, ignTim); + + GIn = CIn*Su_*Su_*exp(CIn*Su_*Su_*ignTm)*GInMult_; + + GIn *= + GIn + / + ( + GIn + + GInFade_*turbXiG + + dimensionedScalar("GSmall", inv(dimTime), SMALL) + ); + + return (GIn + turbXiG); +} + + +Foam::tmp Foam::XiGModels::instability2G::Db() const +{ + const objectRegistry& db = Su_.db(); + const volScalarField& Xi = db.lookupObject("Xi"); + const volScalarField& rho = db.lookupObject("rho"); + const volScalarField& mgb = db.lookupObject("mgb"); + const volScalarField& Db1 = db.lookupObject("Db"); + + //return turbulence_.muEff() + return Db1 + + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_); +} + + +bool Foam::XiGModels::instability2G::read(const dictionary& XiGProperties) +{ + XiGModel::read(XiGProperties); + + XiGModelCoeffs_.readEntry("defaultCIn", defaultCIn_); + XiGModelCoeffs_.readEntry("GInFade", GInFade_); + XiGModelCoeffs_.readEntry("GInMult", GInMult_); + XiGModelCoeffs_.readEntry("lambdaIn", lambdaIn_); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.H b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.H new file mode 100644 index 0000000000..3a5b0600a6 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.H @@ -0,0 +1,141 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiGModels::instability2G + +Description + + +SourceFiles + instability2G.C + +\*---------------------------------------------------------------------------*/ + +#ifndef instability2G_H +#define instability2G_H + +#include "laminarFlameSpeed.H" +#include "SCOPELaminarFlameSpeed.H" +#include "XiGModel.H" +#include "ignitionSite.H" +#include "ignition.H" +#include "Time.H" +#include "fvMesh.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiGModels +{ + +/*---------------------------------------------------------------------------*\ + Class instability2G Declaration +\*---------------------------------------------------------------------------*/ + +class instability2G +: + public XiGModel +{ + // Private data + + laminarFlameSpeedModels::SCOPE saModel_; + + // GIn (inituial instability G)divided by Su^2. Read from fuel file + scalar CIn_; + + //- Default CIn if not in fuel file + scalar defaultCIn_; + + // Determine how fast GIn fades out as turbulence starts + scalar GInFade_; + + // Set GIn large so that XiEq determines Xi value. Son increase byfactor: + scalar GInMult_; + + //- instability2G length-scale + dimensionedScalar lambdaIn_; + + //- Xi generation rate model due to all other processes + autoPtr XiGModel_; + + + // Private Member Functions + + //- Disallow copy construct + instability2G(const instability2G&); + + //- Disallow default bitwise assignment + void operator=(const instability2G&); + + +public: + + //- Runtime type information + TypeName("instability2G"); + + + // Constructors + + //- Construct from components + instability2G + ( + const dictionary& XiGProperties, + const word& modelType, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su + ); + + + //- Destructor + virtual ~instability2G(); + + + // Member Functions + + //- Return the flame-wrinking generation rate + virtual tmp G() const; + + //- Return the flame diffusivity + virtual tmp Db() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiGProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiGModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.C b/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.C new file mode 100644 index 0000000000..4dcfada06c --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.C @@ -0,0 +1,253 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2015 OpenFOAM Foundation + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "transportTwoEqs.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiModels +{ + defineTypeNameAndDebug(transportTwoEqs, 0); + addToRunTimeSelectionTable(XiModel, transportTwoEqs, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::XiModels::transportTwoEqs::transportTwoEqs +( + const dictionary& XiProperties, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su, + const volScalarField& rho, + const volScalarField& b, + const surfaceScalarField& phi +) +: + XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi), + XiShapeCoef_(XiModelCoeffs_.get("XiShapeCoef")), + CpfiDot_(XiModelCoeffs_.get("CpfiDot")), + CpfiCross_(XiModelCoeffs_.get("CpfiCross")), + GEtaExp_(XiModelCoeffs_.get("GEtaExp")), + LOverCw_(XiModelCoeffs_.get("LOverCw")), + XiEqModel_ + ( + XiEqModel::New(XiProperties, "XiEqModel", thermo, turbulence, Su) + ), + XiGModel_(XiGModel::New(XiProperties, "XiGModel", thermo, turbulence, Su)), + XpEqModel_ + ( + XiEqModel::New(XiProperties, "XpEqModel", thermo, turbulence, Su) + ), + XpGModel_ + ( + XiGModel::New(XiProperties, "XpGModel", thermo, turbulence, Su) + ), + Ep_ + ( + IOobject + ( + "Ep", + b.time().timeName(), + b.db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + b.mesh() + ) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::XiModels::transportTwoEqs::~transportTwoEqs() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + + +Foam::tmp Foam::XiModels::transportTwoEqs::Db() const +{ + return XiGModel_->Db(); +} + + +void Foam::XiModels::transportTwoEqs::correct +( + const fv::convectionScheme& mvConvection +) +{ + const volScalarField XiEqEta(XiEqModel_->XiEq()); + volScalarField GEta(XiGModel_->G()); + + GEta *= max( 1.0, exp( GEtaExp_*(1.0 - (Xi_ - 1.0)/(XiEqEta - 0.999)))) ; + + const volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999)); + + const volScalarField XiEqStar(R/(R - GEta)); + + const volScalarField XiEq + ( + 1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XiEqStar - 1.0) + ); + + const volScalarField G(R*(XiEq - 1.0)/XiEq); + + + const objectRegistry& db = b_.db(); + const volScalarField& betav = db.lookupObject("betav"); + const volScalarField& p = db.lookupObject("p"); + const volScalarField& mgb = db.lookupObject("mgb"); + const surfaceScalarField& phiSt = + db.lookupObject("phiSt"); + const volScalarField& Db = db.lookupObject("Db"); + const surfaceScalarField& nf = db.lookupObject("nf"); + + surfaceScalarField phiXi + ( + "phiXi", + phiSt + + ( + - fvc::interpolate(fvc::laplacian(Db, b_)/mgb)*nf + + fvc::interpolate(rho_) + * fvc::interpolate(Su_*(1.0/(Xi_*Xp_) - (Xi_*Xp_)))*nf + ) + ); + + + dimensionedScalar zero + ( + "zero", + dimensionSet(2, -6, -2, 0, 0, 0, 0), + scalar(0.0) + ); + + const volScalarField Gpfi + ( + CpfiDot_*sqrt(max(fvc::grad(rho_)&fvc::grad(p), zero ))/rho_*b_*(1.0-b_) + + CpfiCross_*sqrt(mag(fvc::grad(rho_)^fvc::grad(p) ))/rho_*b_*(1.0-b_) + ); + + fvScalarMatrix XiEqn_ + ( + betav*fvm::ddt(rho_, Xi_) + + mvConvection.fvmDiv(phi_, Xi_) + + fvm::div(phiXi, Xi_) + - fvm::Sp(fvc::div(phiXi), Xi_) + == + betav*rho_*(R + Gpfi ) + - fvm::Sp(betav*rho_*(R - G), Xi_) + ); + + XiEqn_.relax(); + XiEqn_.solve(); + + // Correct boundedness of Xi + // ~~~~~~~~~~~~~~~~~~~~~~~~~ + Xi_.max(1.0); + Xi_ = min(Xi_, 2.0*XiEq); + + // Calculation of Xp generated by obstacles + volScalarField XpEqEta("XpEqEta",XpEqModel_->XiEq()); + + const volScalarField GpEta("GpEta", XpGModel_->G()); + + const volScalarField Rp("Rp", GpEta*XpEqEta/(XpEqEta - 0.999)); + + const volScalarField XpEq + ( + "XpEq", + 1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XpEqEta - 1.0) + ); + + const volScalarField Gpp("Gpp", Rp*(XpEq - 1.0)/XpEq); + + + fvScalarMatrix XpEqn_ + ( + betav*fvm::ddt(rho_, Xp_) + + mvConvection.fvmDiv(phi_, Xp_) + + fvm::div(phiXi, Xp_) + - fvm::Sp(fvc::div(phiXi), Xp_) + == + betav*rho_*Rp + - fvm::Sp(betav*rho_*(Rp - Gpp), Xp_) + ); + + XpEqn_.relax(); + XpEqn_.solve(); + + Xp_.max(1.0); + Xp_ = min(Xp_, 20.0*XpEq); + + // Calculate Ep + const volScalarField& Lobs = db.lookupObject("Lobs"); + const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0); + volScalarField LI(Lobs); + + LI.primitiveFieldRef() = max(LI.primitiveField(),LOverCw_*sqrt(Cw)); + + fvScalarMatrix EpEqn_ + ( + betav*fvm::ddt(rho_, Ep_) + + mvConvection.fvmDiv(phi_, Ep_) + + fvm::div(phiXi, Ep_) + - fvm::Sp(fvc::div(phiXi), Ep_) + == + betav*rho_*Gpp*Xp_/LI + - fvm::Sp(betav*rho_*Rp, Ep_) + ); + + EpEqn_.relax(); + EpEqn_.solve(); + + Ep_.max(0.0); + Ep_.min(100000.0); +} + + +bool Foam::XiModels::transportTwoEqs::read(const dictionary& XiProperties) +{ + XiModel::read(XiProperties); + + XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef_); + XiModelCoeffs_.readEntry("CpfiDot", CpfiDot_); + XiModelCoeffs_.readEntry("CpfiCross", CpfiCross_); + XiModelCoeffs_.readEntry("GEtaExp", GEtaExp_); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.H b/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.H new file mode 100644 index 0000000000..24f1b3dcd8 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2015 OpenFOAM Foundation + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::XiModels::transportTwoEqs + +Description + +SourceFiles + transportTwoEqs.C + +\*---------------------------------------------------------------------------*/ + +#ifndef transportTwoEqs_H +#define transportTwoEqs_H + +#include "XiModel.H" +#include "XiEqModel.H" +#include "XiGModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace XiModels +{ + +/*---------------------------------------------------------------------------*\ + Class transportTwoEqs Declaration +\*---------------------------------------------------------------------------*/ + +class transportTwoEqs +: + public XiModel +{ + // Private data + + scalar XiShapeCoef_; + scalar CpfiDot_; + scalar CpfiCross_; + scalar GEtaExp_; + + //- Maximum Lobs/CellWidth + scalar LOverCw_; + + //- Equilibrium for Xi (turbulence) + autoPtr XiEqModel_; + + //- Generation for Xi (turbulence) + autoPtr XiGModel_; + + //- Equilibrium for Xp (obstacles) + autoPtr XpEqModel_; + + //- Generation for Xp (obstacles) + autoPtr XpGModel_; + + //- Dissipation lenght scale for subgrid obstacles + volScalarField Ep_; + + + // Private Member Functions + + //- Disallow copy construct + transportTwoEqs(const transportTwoEqs&); + + //- Disallow default bitwise assignment + void operator=(const transportTwoEqs&); + + +public: + + //- Runtime type information + TypeName("transportTwoEqs"); + + + // Constructors + + //- Construct from components + transportTwoEqs + ( + const dictionary& XiProperties, + const psiuReactionThermo& thermo, + const compressible::RASModel& turbulence, + const volScalarField& Su, + const volScalarField& rho, + const volScalarField& b, + const surfaceScalarField& phi + ); + + + //- Destructor + virtual ~transportTwoEqs(); + + + // Member Functions + + //- Return the flame diffusivity + virtual tmp Db() const; + + //- Correct the flame-wrinking Xi + virtual void correct() + { + NotImplemented; + } + + //- Correct the flame-wrinking Xi using the given convection scheme + virtual void correct(const fv::convectionScheme& mvConvection); + + //- Update properties from given dictionary + virtual bool read(const dictionary& XiProperties); + + //- Write fields of the XiEq model + virtual void writeFields() + { + XiEqModel_().writeFields(); + XpEqModel_().writeFields(); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace XiModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/combustion/PDRFoam/readTimeControls.H b/applications/solvers/combustion/PDRFoam/readTimeControls.H new file mode 100644 index 0000000000..b35ebd100a --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/readTimeControls.H @@ -0,0 +1,47 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2015 OpenFOAM Foundation + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Global + readTimeControls + +Description + Read the control parameters used by setDeltaT + +\*---------------------------------------------------------------------------*/ + +const bool adjustTimeStep = + runTime.controlDict().getOrDefault("adjustTimeStep", false); + +scalar maxCo = runTime.controlDict().getOrDefault("maxCo", 1.0); + +scalar maxDeltaT = + runTime.controlDict().getOrDefault("maxDeltaT", GREAT); + +scalar maxDeltaTRatio = + runTime.controlDict().getOrDefault("maxDeltaTRatio", 1.2); + + +// ************************************************************************* //