diff --git a/solvers_post/AMC/AMC.C b/solvers_post/AMC/AMC.C new file mode 100644 index 0000000..af2d84b --- /dev/null +++ b/solvers_post/AMC/AMC.C @@ -0,0 +1,68 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "AMC.H" + +#include "error.H" + +// * * * * * * * * * * * * * Functions * * * * * * * * * * * * // + +//Amplitude Mapping Closure (from KIVA) +//Define exp(-2*(erf^-1(2*eta - 1))^2) +inline Foam::scalar Foam::CMC::AMC(const scalar eta) +{ + const scalar a0 = (2.0*eta - 1.0); + scalar a1 = 0.5; + scalar slope; + scalar da = GREAT; + + while(mag(da) > SMALL) + { + slope = (2.0/spi)*Foam::exp(-1.0*Foam::pow(a1,2.0)); + da = (a0 - Foam::erf(a1))/slope; + a1 = a1 + da; + } + + return Foam::exp(-2.0*Foam::pow(a1,2.0)); +} + +//Amplitude Mapping Closure (from KIVA) +Foam::tmp Foam::CMC::AMC(const scalarField& eta) +{ + tmp tRes(new scalarField(eta.size(), 0.0)); + + scalarField& Res = tRes.ref(); + + for(label i=0 ; i. + +Function + Foam::CMC::AMC + +Description + +SourceFiles + AMC.C + +\*---------------------------------------------------------------------------*/ + +#ifndef AMC_H +#define AMC_H + +#include "tmp.H" +#include "scalarField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace CMC +{ + +/*---------------------------------------------------------------------------*\ + Function AMC Declaration +\*---------------------------------------------------------------------------*/ + + const scalar pi = 3.141592; + + const scalar spi = Foam::sqrt(pi); + + //Amplitude Mapping Closure (from KIVA) + tmp AMC(const scalarField& eta); + + + //Amplitude Mapping Closure (from KIVA) + //Define exp(-2*(erf^-1(2*eta - 1))^2) + inline scalar AMC(const scalar eta); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CMC +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/solvers_post/BetaFunction/BetaFunction.C b/solvers_post/BetaFunction/BetaFunction.C new file mode 100644 index 0000000..b6a13ae --- /dev/null +++ b/solvers_post/BetaFunction/BetaFunction.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaFunction.H" + +#include "error.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// const scalar Foam::CMC::BetaFunction::staticData(); + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CMC::BetaFunction::_setParameters() +{ + scalar gamma = mf_*(1.0 - mf_)/(mfVar_ + SMALL) - 1.0; + + if (gamma < SMALL || (mfVar_)/(mf_*(1 - mf_) + SMALL) < 0.001) + { + fdelta_ = true; + } + else + { + a_ = max(mf_*gamma, 0.0); + b_ = max((1.0 - mf_)*gamma, 0.0); + + if(a_ < 1.0) + { + delta0_ = true; + } + if(b_ < 1.0) + { + delta1_ = true; + } + + _limitAB(); + } +} + + +void Foam::CMC::BetaFunction::_limitAB() +{ + scalar fmax = 1.0/(1.0 + (b_ - 1.0)/(a_ - 1.0)); + + if(a_ > 500.0) + { + a_ = 500.0; + b_ = (a_ - 1.0 - fmax*(a_ - 2.0))/fmax; + } + else if(b_ > 500.0) + { + b_ = 500.0; + a_ = (1.0 + fmax*(b_ - 2.0))/(1.0 - fmax); + } +} + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CMC::BetaFunction::BetaFunction(const scalar mf, const scalar mfVar) +: + mf_(mf), + mfVar_(mfVar), + a_(0.0), + b_(0.0), + fdelta_(false), + delta0_(false), + delta1_(false) +{ + _setParameters(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::CMC::BetaFunction::~BetaFunction() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaFunction/BetaFunction.H b/solvers_post/BetaFunction/BetaFunction.H new file mode 100644 index 0000000..ca353c6 --- /dev/null +++ b/solvers_post/BetaFunction/BetaFunction.H @@ -0,0 +1,198 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +Class + Foam::BetaFunction + +Description + +SourceFiles + BetaFunctionI.H + BetaFunction.C + BetaFunctionIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BetaFunction_H +#define BetaFunction_H + +#include "scalar.H" +#include "Switch.H" +#include "scalarField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class Istream; +class Ostream; + +// Forward declaration of friend functions and operators +namespace CMC +{ + class BetaFunction; +} +Istream& operator>>(Istream&, CMC::BetaFunction&); +Ostream& operator<<(Ostream&, const CMC::BetaFunction&); + + +/*---------------------------------------------------------------------------*\ + Class BetaFunction Declaration +\*---------------------------------------------------------------------------*/ + +namespace CMC +{ + +class BetaFunction +{ + // Private data + + //- Mean mixture fraction + scalar mf_; + + //- Mixture fraction variance + scalar mfVar_; + + //- Beta pdf parameter alpha + scalar a_; + + //- Beta pdf parameter beta + scalar b_; + + //- Flag beta pdf special case: forced-delta function + Switch fdelta_; + + //- delta function at pure oxidizer stream + Switch delta0_; + + //- delta function at pure fuel stream + Switch delta1_; + + + // Private Member Functions + + //- Calculate alpha and beta and set special case flags + void _setParameters(); + + //- Limit too large value of beta pdf parameters + void _limitAB(); + + //- Disallow default bitwise copy construct + // BetaFunction(const BetaFunction&); + + //- Disallow default bitwise assignment + // void operator=(const BetaFunction&); + + +public: + + // Static data members + + //- Static data staticData + // static const scalar staticData; + + + // Constructors + + //- Construct from components + BetaFunction(const scalar mf, const scalar mfVar); + + //- Construct from Istream + BetaFunction(Istream&); + + //- Construct as copy + // BetaFunction(const BetaFunction&); + + + //- Destructor + ~BetaFunction(); + + + // Member Functions + + scalarField etaFunc(const UList& eta) const + { + return pow(eta, a_-1.0)*pow(1.0-eta, b_-1.0); + } + + scalar etaFunc(const scalar eta) const + { + return Foam::pow(eta, a_-1.0)*Foam::pow(1.0-eta, b_-1.0); + } + + + // Access + + Switch fdelta() const {return fdelta_;} + + Switch delta0() const {return delta0_;} + + Switch delta1() const {return delta1_;} + + scalar mf() const {return mf_;} + + scalar mfVar() const {return mfVar_;} + + scalar alpha() const {return a_;} + + scalar beta() const {return b_;} + + // Check + + // Edit + + // Write + + + // Member Operators + + // void operator=(const BetaFunction&); + + + // Friend Functions + + // Friend Operators + + // IOstream Operators + + friend Istream& ::Foam::operator>>(Istream&, BetaFunction&); + friend Ostream& ::Foam::operator<<(Ostream&, const BetaFunction&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CMC +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "BetaFunctionI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/solvers_post/BetaFunction/BetaFunctionI.H b/solvers_post/BetaFunction/BetaFunctionI.H new file mode 100644 index 0000000..eca8d9e --- /dev/null +++ b/solvers_post/BetaFunction/BetaFunctionI.H @@ -0,0 +1,58 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaFunction/BetaFunctionIO.C b/solvers_post/BetaFunction/BetaFunctionIO.C new file mode 100644 index 0000000..9d568fa --- /dev/null +++ b/solvers_post/BetaFunction/BetaFunctionIO.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaFunction.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CMC::BetaFunction::BetaFunction(Istream& is) +{ + // Check state of Istream + is.check("Foam::CMC::BetaFunction::BetaFunction(Foam::Istream&)"); +} + + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaFunction&) +{ + // Check state of Istream + is.check + ( + "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaFunction&)" + ); + + return is; +} + + +Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaFunction&) +{ + // Check state of Ostream + os.check + ( + "Foam::Ostream& Foam::operator<<(Foam::Ostream&, " + "const Foam::CMC::BetaFunction&)" + ); + + return os; +} + + +// ************************************************************************* // diff --git a/solvers_post/BetaGrid/BetaGrid.C b/solvers_post/BetaGrid/BetaGrid.C new file mode 100644 index 0000000..71b4211 --- /dev/null +++ b/solvers_post/BetaGrid/BetaGrid.C @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaGrid.H" + +#include "error.H" +#include "dictionary.H" +#include "SubList.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// const scalar Foam::CMC::BetaGrid::staticData(); + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CMC::BetaGrid::BetaGrid(const dictionary& dict) +: + etaCut_(dict.lookup("detailedEta")), + nSpacings_(dict.lookup("detailedN")), + etaSpace_(sum(nSpacings_)+3) +{ + // Check for input list size + if (etaCut_.size() != nSpacings_.size() + 1) + { + FatalErrorInFunction + << "Number of grid interval points does not match number of grid sizes" + << abort(FatalError); + } + + // Check for evenness of number of spacings in each interval + forAll (nSpacings_, i) + { + if ((nSpacings_[i] % 2) == 1) + { + FatalErrorInFunction + << "Number of spacings in intervals should be even" + << abort(FatalError); + } + } + + forAll(nSpacings_, i) + { + const label ns = nSpacings_[i]; + const labelList::subList prev(nSpacings_, i); + const label baseI = sum(prev) + 1; + + const scalar low = etaCut_[i]; + const scalar upp = etaCut_[i+1]; + const scalar delta = (upp - low) / scalar(ns); + + for (label j = 0; j < ns; j++) + { + etaSpace_[baseI+j] = low + delta * j; + } + } + etaSpace_[0] = 0.0; + etaSpace_[etaSpace_.size()-2] = etaCut_.last(); + etaSpace_[etaSpace_.size()-1] = 1.0; +} + + +Foam::CMC::BetaGrid::BetaGrid(const BetaGrid&) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::CMC::BetaGrid::~BetaGrid() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // + +void Foam::CMC::BetaGrid::operator=(const BetaGrid& rhs) +{ + // Check for assignment to self + if (this == &rhs) + { + FatalErrorInFunction + << "Attempted assignment to self" + << abort(FatalError); + } +} + +// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaGrid/BetaGrid.H b/solvers_post/BetaGrid/BetaGrid.H new file mode 100644 index 0000000..fc17e0a --- /dev/null +++ b/solvers_post/BetaGrid/BetaGrid.H @@ -0,0 +1,161 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +Class + Foam::BetaGrid + +Description + +SourceFiles + BetaGridI.H + BetaGrid.C + BetaGridIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BetaGrid_H +#define BetaGrid_H + +#include "scalarField.H" +#include "labelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class Istream; +class Ostream; + +// Forward declaration of friend functions and operators +namespace CMC +{ + class BetaGrid; +} +class dictionary; + +Istream& operator>>(Istream&, CMC::BetaGrid&); +Ostream& operator<<(Ostream&, const CMC::BetaGrid&); + + +/*---------------------------------------------------------------------------*\ + Class BetaGrid Declaration +\*---------------------------------------------------------------------------*/ + +namespace CMC +{ + +class BetaGrid +{ + // Private data + + //- Piecewise uniform interval points + scalarField etaCut_; + + //- Sizes of piecewise uniform intervals + labelList nSpacings_; + + //- Mixture fraction eta grid points + scalarField etaSpace_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + BetaGrid(const BetaGrid&); + + //- Disallow default bitwise assignment + void operator=(const BetaGrid&); + + +public: + + // Static data members + + //- Static data staticData + // static const scalar staticData; + + + // Constructors + + //- Construct from components + BetaGrid(const dictionary& dict); + + //- Construct from Istream + BetaGrid(Istream&); + + //- Construct as copy + // BetaGrid(const BetaGrid&); + + + //- Destructor + ~BetaGrid(); + + + // Member Functions + + // Access + const labelList& N() const {return nSpacings_;} + + const scalarField& etaCut() const {return etaCut_;} + + const scalarField& etaSpace() const {return etaSpace_;} + + // Check + + // Edit + + // Write + + + // Member Operators + + // void operator=(const BetaGrid&); + + + // Friend Functions + + // Friend Operators + + // IOstream Operators + + friend Istream& ::Foam::operator>>(Istream&, BetaGrid&); + friend Ostream& ::Foam::operator<<(Ostream&, const BetaGrid&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CMC +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "BetaGridI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/solvers_post/BetaGrid/BetaGridI.H b/solvers_post/BetaGrid/BetaGridI.H new file mode 100644 index 0000000..eca8d9e --- /dev/null +++ b/solvers_post/BetaGrid/BetaGridI.H @@ -0,0 +1,58 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaGrid/BetaGridIO.C b/solvers_post/BetaGrid/BetaGridIO.C new file mode 100644 index 0000000..d567d45 --- /dev/null +++ b/solvers_post/BetaGrid/BetaGridIO.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaGrid.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CMC::BetaGrid::BetaGrid(Istream& is) +{ + // Check state of Istream + is.check("Foam::CMC::BetaGrid::BetaGrid(Foam::Istream&)"); +} + + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaGrid&) +{ + // Check state of Istream + is.check + ( + "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaGrid&)" + ); + + return is; +} + + +Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaGrid&) +{ + // Check state of Ostream + os.check + ( + "Foam::Ostream& Foam::operator<<(Foam::Ostream&, " + "const Foam::CMC::BetaGrid&)" + ); + + return os; +} + + +// ************************************************************************* // diff --git a/solvers_post/BetaIntegrator/BetaIntegrator.C b/solvers_post/BetaIntegrator/BetaIntegrator.C new file mode 100644 index 0000000..07ca064 --- /dev/null +++ b/solvers_post/BetaIntegrator/BetaIntegrator.C @@ -0,0 +1,161 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaIntegrator.H" + + +#include "error.H" +#include "SubList.H" +#include "SubField.H" +#include "interpolateXY.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// const scalar Foam::CMC::BetaIntegrator::staticData(); + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +//extended Simpson's rule (Numerical recipes, 2nd Ed. p.128) +//for equally spaced and even N intervals (or odd N+1 points) +Foam::scalar Foam::CMC::BetaIntegrator::simps +(const scalar xl, const scalar xh, const label n, const UList& fx) + const +{ + scalar evensum(0.0), oddsum(0.0), sum(0.0); + + scalar h = (xh - xl)/scalar(n); + + for(label i=0 ; i& f) const +{ + scalar total = 0.0; + + const labelList &N(bg_.N()); + const scalarField &etaCut(bg_.etaCut()); + + forAll(N, i) + { + const labelList::subList prev(N, i); + const label baseI = sum(prev) + 1; + + const scalarField::subField subInterval(f, N[i]+1, baseI); + + total += simps(etaCut[i], etaCut[i+1], N[i], subInterval); + } + + return total; +} + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CMC::BetaIntegrator::BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg) +: + etaSpace_(bg.etaSpace()), + pdfNum_(bg.etaSpace().size(), 0.0), + pdfDen_(0.0), + bf_(bf), + bg_(bg) +{ + const scalarField::subField domain(etaSpace_, etaSpace_.size()-2, 1); + scalarField::subField range(pdfNum_, etaSpace_.size()-2, 1); + range = bf_.etaFunc(domain); + pdfDen_ = sumSimps(pdfNum_) + + Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL) + + Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::CMC::BetaIntegrator::~BetaIntegrator() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +// beta-pdf weighted integration for given mf, mfVar and f +Foam::scalar Foam::CMC::BetaIntegrator::betaIntegrate(const scalarField& f) const +{ + scalar result = 0.0; + + if(bf_.fdelta()) + { + result = interpolateXY(bf_.mf(), etaSpace_, f); + } + else + { + scalar num = sumSimps (f * pdfNum_) + + f.first()*Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL) + + f.last()*Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL); + result = num/pdfDen_; + } + + return result; +} + +// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // + +void Foam::CMC::BetaIntegrator::operator=(const BetaIntegrator& rhs) +{ + // Check for assignment to self + if (this == &rhs) + { + FatalErrorInFunction + << "Attempted assignment to self" + << abort(FatalError); + } +} + +// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaIntegrator/BetaIntegrator.H b/solvers_post/BetaIntegrator/BetaIntegrator.H new file mode 100644 index 0000000..a0ee5e2 --- /dev/null +++ b/solvers_post/BetaIntegrator/BetaIntegrator.H @@ -0,0 +1,166 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +Class + Foam::CMC::BetaIntegrator + +Description + +SourceFiles + BetaIntegratorI.H + BetaIntegrator.C + BetaIntegratorIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BetaIntegrator_H +#define BetaIntegrator_H + +#include "scalarField.H" +#include "BetaFunction.H" +#include "BetaGrid.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class Istream; +class Ostream; + +// Forward declaration of friend functions and operators +namespace CMC +{ + class BetaIntegrator; +} +Istream& operator>>(Istream&, CMC::BetaIntegrator&); +Ostream& operator<<(Ostream&, const CMC::BetaIntegrator&); + + +/*---------------------------------------------------------------------------*\ + Class BetaIntegrator Declaration +\*---------------------------------------------------------------------------*/ + +namespace CMC +{ + +class BetaIntegrator +{ + // Private data + + //- Beta pdf numerators + scalarField pdfNum_; + + //- Beta pdf denominator + scalar pdfDen_; + + //- Quadrature points + const scalarField &etaSpace_; + + //- Description of data_ + BetaFunction bf_; + + //- Description of data_ + const BetaGrid &bg_; + + + // Private Member Functions + + //- Calculate integral using Simpson's rule + scalar simps(const scalar xl, const scalar xh, const label n, const UList& fx) const; + + //- Calculate integral using Simpson's rule + scalar sumSimps(const UList& f) const; + + //- Disallow default bitwise copy construct + BetaIntegrator(const BetaIntegrator&); + + //- Disallow default bitwise assignment + void operator=(const BetaIntegrator&); + + +public: + + // Static data members + + //- Static data staticData + // static const scalar staticData; + + + // Constructors + + //- Construct from components + BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg); + + //- Construct as copy + // BetaIntegrator(const BetaIntegrator&); + + + //- Destructor + ~BetaIntegrator(); + + + // Member Functions + + scalar betaIntegrate(const scalarField& f) const; + + // Access + + // Check + + // Edit + + // Write + + + // Member Operators + + // void operator=(const BetaIntegrator&); + + + // Friend Functions + + // Friend Operators + + // IOstream Operators + + friend Istream& ::Foam::operator>>(Istream&, BetaIntegrator&); + friend Ostream& ::Foam::operator<<(Ostream&, const BetaIntegrator&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CMC +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "BetaIntegratorI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/solvers_post/BetaIntegrator/BetaIntegratorI.H b/solvers_post/BetaIntegrator/BetaIntegratorI.H new file mode 100644 index 0000000..eca8d9e --- /dev/null +++ b/solvers_post/BetaIntegrator/BetaIntegratorI.H @@ -0,0 +1,58 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/solvers_post/BetaIntegrator/BetaIntegratorIO.C b/solvers_post/BetaIntegrator/BetaIntegratorIO.C new file mode 100644 index 0000000..8818197 --- /dev/null +++ b/solvers_post/BetaIntegrator/BetaIntegratorIO.C @@ -0,0 +1,56 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "BetaIntegrator.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaIntegrator&) +{ + // Check state of Istream + is.check + ( + "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaIntegrator&)" + ); + + return is; +} + + +Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaIntegrator&) +{ + // Check state of Ostream + os.check + ( + "Foam::Ostream& Foam::operator<<(Foam::Ostream&, " + "const Foam::CMC::BetaIntegrator&)" + ); + + return os; +} + + +// ************************************************************************* // diff --git a/solvers_post/SLFMFoam/BetaPDF.H b/solvers_post/SLFMFoam/BetaPDF.H deleted file mode 100644 index 127dc25..0000000 --- a/solvers_post/SLFMFoam/BetaPDF.H +++ /dev/null @@ -1,312 +0,0 @@ -class BetaPDF -{ - - //private variable - - //raw parameter mixture fraction and mf variance - scalar mf_, mfVar_; - - //beta-pdf parameter alpha and beta - scalar alpha_, beta_; - - //cutting points - scalarField etaCut_; - - //number of eta-space - labelList N_; - - //detailed integration space - scalarField etaSpace_; - - //pdf numerators - scalarField pdfNum_; - - //pdf denominator - scalar pdfDen_; - - //detailed integration space, part - typedef List scalarFieldArray1d; - scalarFieldArray1d etaPart_; - - //AMC, exp(-2*(erf^-1(2*eta - 1))^2) field - //for detailed integration space - scalarField AMCfine_; - - //flag to check forced-delta ftn - //or delta ftn at oxidizer or fuel - bool fdelta_ = false, delta_ox = false, delta_fu = false; - -public: - - // Constructor - BetaPDF(IOdictionary& SLFMdict) - : - alpha_(0), - beta_(0), - etaCut_(SLFMdict.lookup("detailedEta")), - N_(SLFMdict.lookup("detailedN")), - etaPart_(N_.size()) - { - //Ref. F.Liu et al., INT. J. THERM. SCI. 41 (2002) 763-772. - scalar del(0); - label cnt(0); - - etaSpace_.append(etaCut_[0]); - for(label i=0 ; i 500.0) - { - alpha_ = 500.0; - beta_ = (alpha_-1.0-fmax*(alpha_-2.0))/fmax; - } - else if(beta_ > 500.0) - { - beta_ = 500.0; - alpha_ = (1.0+fmax*(beta_-2.0))/(1.0-fmax); - } - } - scalarField etaFunc(const scalar a, const scalar b, const scalarField& eta) const - { - return pow(eta, a-1.0)*pow(1.0-eta, b-1.0); - } - scalar etaFunc(const scalar a, const scalar b, const scalar eta) const - { - return Foam::pow(eta, a-1.0)*Foam::pow(1.0-eta, b-1.0); - } - //extended Simpson's rule (Numerical recipes, 2nd Ed. p.128) - //for equally spaced and even N intervals (or odd N+1 points) - scalar simps(const scalar xl, const scalar xh, const label N, const UList& fx) const - { - scalar evensum(0.0), oddsum(0.0), sum(0.0); - - scalar h = (xh - xl)/scalar(N); - - for(label i=0 ; i SMALL) - { - slope = 2.0/spi*Foam::exp(-1.0*Foam::pow(a1,2.0)); - da = (a0 - Foam::erf(a1))/slope; - a1 = a1+da; - } - - result[i] = Foam::exp(-2.0*Foam::pow(a1,2.0)); - } - - return result; - } - scalar AMC(const scalar eta) const - { - return interpolateXY(eta, etaSpace_, AMCfine_); - } - void C1coeff(const scalar mf, const scalarField& varValue, scalarField& C1table) - { - scalar maxVar = mf*(1.0-mf); - scalarField x, fx; - - x.append(0.0); - x.append(etaSpace_); - x.append(1.0); - - fx.append(0.0); - fx.append(AMCfine_); - fx.append(0.0); - - for(label v=1 ; v -Foam::CMC::FlameStructure::New() + NstDir_(NstDir), + NstList_(NstList), + eta_(etaGrid), + species_(spTable), + Y_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))), + W_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))), + T_(NstList_.size(), scalarField(eta_.size(), 0.0)), + Q_(NstList_.size(), scalarField(eta_.size(), 0.0)), + h_(NstList_.size(), scalarField(eta_.size(), 0.0)), + rho_(NstList_.size(), scalarField(eta_.size(), 0.0)), + Rgas_(NstList_.size(), scalarField(eta_.size(), 0.0)) { - return autoPtr(new FlameStructure); + for (label n = 0; n < NstList_.size(); n++) + { + scalar N = NstList_[n]; + + fileName fname = "sdr"+Foam::name(N)+".inp"; + IFstream fin(NstDir_/fname); + + string gbg; + label etamax_SLFM, NoSpecies; + + //INPUT FILE... + fin.getLine(gbg); + + //No.GRID... + fin.getLine(gbg); + + //Number of eta-point and species + fin >> etamax_SLFM >> NoSpecies; + + //blank line + fin.getLine(gbg); + + //PRESSURE... + fin.getLine(gbg); + + //1.00... + fin.getLine(gbg); + + //MIXTURE... + fin.getLine(gbg); + + //read eta space (from SLFM library) + scalarField etaValue_SLFM(etamax_SLFM, 0.0); + for(label j = 0 ; j < etamax_SLFM ; j++) + { + fin>>etaValue_SLFM[j]; + } + + //blank line + fin.getLine(gbg); + + //INITIAL... + fin.getLine(gbg); + + //species loop + scalarField temp(etamax_SLFM, 0.0); + for(label i = 0 ; i < NoSpecies ; i++) + { + word spName; + //H2... (species name) + fin >> spName; + const label spI (species_[spName]); + + // Info << spI << spName << endl; + + //read species mass fraction (from SLFM library) + for(label j=0 ; j> Yj; + temp[j] = max(0.0, Yj); + } + + Y_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp); + + //read reaction rate (from SLFM library) + for(label j = 0 ; j < etamax_SLFM ; j++) + { + fin>>temp[j]; + } + + W_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp); + + //blank line + fin.getLine(gbg); + } + + //read temperature + { + //H2... (species name) + fin.getLine(gbg); + + //read temperature + for(label j = 0 ; j < etamax_SLFM ; j++) + { + scalar Yj = 0.0; + fin >> Yj; + temp[j] = max(0.0, Yj); + } + + T_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp); + + //read heat source + for(label j = 0 ; j < etamax_SLFM ; j++) + { + fin>>temp[j]; + } + + Q_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp); + + //blank line + fin.getLine(gbg); + } + + //DENSITY... + fin.getLine(gbg); + + //read density (from SLFM library) + for(label j = 0 ; j < etamax_SLFM ; j++) + { + fin>>temp[j]; + } + //Info<>temp[j]; + } + + h_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp); + } } -*/ // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -Foam::CMC::FlameStructure::~FlameStructure() +Foam::SLFM::FlameStructure::~FlameStructure() {} @@ -85,16 +196,6 @@ Foam::CMC::FlameStructure::~FlameStructure() // * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // -void Foam::CMC::FlameStructure::operator=(const FlameStructure& rhs) -{ - // Check for assignment to self - if (this == &rhs) - { - FatalErrorInFunction - << "Attempted assignment to self" - << abort(FatalError); - } -} // * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // diff --git a/solvers_post/SLFMFoam/FlameStructure/FlameStructure.H b/solvers_post/SLFMFoam/FlameStructure/FlameStructure.H index 864f257..96a5ba8 100644 --- a/solvers_post/SLFMFoam/FlameStructure/FlameStructure.H +++ b/solvers_post/SLFMFoam/FlameStructure/FlameStructure.H @@ -22,7 +22,7 @@ License along with OpenFOAM. If not, see . Class - Foam::CMC::FlameStructure + Foam::SLFM::FlameStructure Description @@ -48,16 +48,18 @@ namespace Foam class Istream; class Ostream; -namespace CMC +class hashedWordList; + +namespace SLFM { -class FlameStructure; + class FlameStructure; } // Forward declaration of friend functions and operators -Istream& operator>>(Istream&, CMC::FlameStructure&); -Ostream& operator<<(Ostream&, const CMC::FlameStructure&); +Istream& operator>>(Istream&, SLFM::FlameStructure&); +Ostream& operator<<(Ostream&, const SLFM::FlameStructure&); -namespace CMC +namespace SLFM { /*---------------------------------------------------------------------------*\ @@ -77,7 +79,13 @@ class FlameStructure word NstDir_; //- Species mass fraction - scalarList NstList_; + const scalarList NstList_; + + //- Species mass fraction + const scalarField eta_; + + //- Species mass fraction + const hashedWordList &species_; //- Species mass fraction scalarFieldArray2d Y_; @@ -107,13 +115,10 @@ class FlameStructure FlameStructure(); //- Disallow default bitwise copy construct - FlameStructure(const FlameStructure&); + // FlameStructure(const FlameStructure&); //- Disallow default bitwise assignment - void operator=(const FlameStructure&); - - //- limit value of alpha and beta - void limitAB(); + // void operator=(const FlameStructure&); public: @@ -127,10 +132,10 @@ public: // Constructors //- Construct from components - FlameStructure(const dictionary &dict, const label nEta, const label nY); + FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList& spTable); //- Construct from Istream - FlameStructure(Istream&); + // FlameStructure(Istream&); //- Construct as copy // FlameStructure(const FlameStructure&); @@ -196,7 +201,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace CMC +} // End namespace SLFM } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/solvers_post/SLFMFoam/FlameStructure/FlameStructureIO.C b/solvers_post/SLFMFoam/FlameStructure/FlameStructureIO.C index 12dbe5c..75b6748 100644 --- a/solvers_post/SLFMFoam/FlameStructure/FlameStructureIO.C +++ b/solvers_post/SLFMFoam/FlameStructure/FlameStructureIO.C @@ -28,41 +28,41 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::CMC::FlameStructure::FlameStructure(Istream& is) /* +Foam::SLFM::FlameStructure::FlameStructure(Istream& is) : base1(is), base2(is), member1(is), member2(is) - */ { // Check state of Istream - is.check("Foam::CMC::FlameStructure::FlameStructure(Foam::Istream&)"); + is.check("Foam::SLFM::FlameStructure::FlameStructure(Foam::Istream&)"); } + */ // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -Foam::Istream& Foam::operator>>(Istream& is, CMC::FlameStructure&) +Foam::Istream& Foam::operator>>(Istream& is, SLFM::FlameStructure&) { // Check state of Istream is.check ( - "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::FlameStructure&)" + "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::SLFM::FlameStructure&)" ); return is; } -Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::FlameStructure&) +Foam::Ostream& Foam::operator<<(Ostream& os, const SLFM::FlameStructure&) { // Check state of Ostream os.check ( "Foam::Ostream& Foam::operator<<(Foam::Ostream&, " - "const Foam::CMC::FlameStructure&)" + "const Foam::SLFM::FlameStructure&)" ); return os; diff --git a/solvers_post/SLFMFoam/Make/files b/solvers_post/SLFMFoam/Make/files index 93044ed..2df083b 100644 --- a/solvers_post/SLFMFoam/Make/files +++ b/solvers_post/SLFMFoam/Make/files @@ -1,5 +1,17 @@ FlameStructure/FlameStructure.C FlameStructure/FlameStructureIO.C + +../AMC/AMC.C + +../BetaFunction/BetaFunction.C +../BetaFunction/BetaFunctionIO.C + +../BetaGrid/BetaGrid.C +../BetaGrid/BetaGridIO.C + +../BetaIntegrator/BetaIntegrator.C +../BetaIntegrator/BetaIntegratorIO.C + SLFMFoam.C EXE = $(FOAM_USER_APPBIN)/SLFMFoam diff --git a/solvers_post/SLFMFoam/Make/options b/solvers_post/SLFMFoam/Make/options index 5239a88..365bb2f 100644 --- a/solvers_post/SLFMFoam/Make/options +++ b/solvers_post/SLFMFoam/Make/options @@ -1,6 +1,10 @@ DEV_PATH=../../libs EXE_INC = \ + -I../AMC \ + -I../BetaFunction \ + -I../BetaGrid \ + -I../BetaIntegrator \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \ diff --git a/solvers_post/SLFMFoam/Math.H b/solvers_post/SLFMFoam/Math.H deleted file mode 100644 index 87a24e5..0000000 --- a/solvers_post/SLFMFoam/Math.H +++ /dev/null @@ -1,401 +0,0 @@ -#include - -//Define ln(gamma_ftn(xx)) -inline double gammaln(double XX) -{ - double STP = 2.5066282746310005; - const double HALF = 0.5, ONE = 1, FPF = 5.5; - double X, Y, TMP, SER, ga; - - static double COF[] = { - 76.18009172947146, - -86.50532032941677, - 24.01409824083091, - -1.231739572450155, - 0.001208650973866179, - -0.000005395239384953 - }; - - X=XX; - Y=X; - TMP = X + FPF; - TMP = (X+HALF)*std::log(TMP)-TMP; - SER = 1.000000000190015; - - for(int j = 0; j <= 5 ; j++) - { - Y = Y + ONE; - SER = SER+COF[j] / Y ; - } - - ga = TMP + std::log(STP*SER/X) ; - - return ga; -} - -//Define exp(-2*(erf^-1(2*eta - 1))^2) -inline double AMC(double eta) //from kiva -{ - - const double pi = 3.141592; - const double spi = std::sqrt(pi); - double a1=0.5, a0=(2*eta - 1), slope, da=0, result=0; - if(eta < 0.000001 || eta > 0.999999) - { - result = 1e-30; - } - else - { -AMC_while:; - slope = 2/spi*std::exp(-1* std::pow(a1,2) ); - da = (a0 - Foam::erf(a1))/slope; - a1 = a1+da; - if( std::fabs(da) > 1e-12 ) //abs function for double variable - { - goto AMC_while; - } - result = std::exp(-2* std::pow(a1,2) ); - } - return result; -} - -//Define TDMA function -inline void TDMA(scalarField& a, scalarField& b, scalarField& c, scalarField& Qi, scalarField& QiNew, scalar Numofetaspace) -{ - scalar bet; - scalarField gam(Numofetaspace,0); - - bet = b[0]; - QiNew[0] = Qi[0]/bet; - - for(label j=1 ; j< Numofetaspace ; ) - { - gam[j] = c[j-1]/bet; - bet = b[j] - a[j]*gam[j]; - QiNew[j] = ( Qi[j] - a[j]*QiNew[j-1] ) / bet; - j = j+1; - } - - for(label j = Numofetaspace - 2 ; j >= 0 ;) - { - QiNew[j] = QiNew[j] - gam[j+1] * QiNew[j+1]; - j = j-1; - } -} - -inline double simpson(scalar min, scalar max, scalar interval, scalarField& ftn) -{ - double result_odd=0, result_even=0, result=0; - for(label i=min+1 ; i < max ; i=i+2) - { - //odd points - result_odd += ftn[i]; - } - result_odd = result_odd*(4.0/3.0)*interval; - - for(label i=min+2 ; i < max ; i=i+2) - { - //even points - result_even += ftn[i]; - } - result_even = result_even*(2.0/3.0)*interval; - - //start and end points - result = result_odd + result_even + (ftn[min]+ftn[max])*interval*(1.0/3.0); - - return result; -} - -inline double integration(scalar deltaftn, scalarField& MFcut, scalarField& Neta, scalarField& pdf, scalarField& f) -{ - double result = 0.0; - double min, max, interval; - if(deltaftn < 0.5) //one delta function (deltaftn == 0) - { - for(label i=0; i<= (Neta[MFcut.size()-1]);i++) - { - if(pdf[i] > 1e-30) - { - result = f[i]; - } - } - } - else if(deltaftn > 0.5 && deltaftn < 1.5) //two delta function (deltaftn == 1) - { - for(label i=0; i<= (Neta[MFcut.size()-1]);i++) - { - if(pdf[i] > 1e-30) - { - result += 0.5*f[i]; - } - } - } - else if(deltaftn > 1.5) //not a delta function (deltaftn == 2) - { - scalarField ftn = pdf*f; - for(label i=0; i<(MFcut.size()-1); i++ ) - { - min = Neta[i]; - max = Neta[i+1]; - interval = (MFcut[i+1]-MFcut[i])/(max-min); - result += simpson(min, max, interval, ftn); - } - } - else - { - //Integration error - } - return result; -} - -//inline double trapz() -//{ -//} - -//inline double rectang() -//{ -//} - -//Returns lognormal PDF of the Nst -//Ref) S.H.Kim and Kang Y. Huh, Combust Flame, vol.138, pp.336-352, 2004 -// Equation 17 -inline double lognormalPDF(scalar Nst, scalar Nst_local) -{ - double result, numerator; - const double pi = 3.141592; - const double sigmaN = 1.0; - - numerator = sqr(std::log(Nst/Nst_local)+0.5*sqr(sigmaN)); - result = std::exp(-0.5*numerator/sqr(sigmaN)); - result = result/std::sqrt(2.0*pi)/Nst/sigmaN; - - return result; -} - -//Returns the integral of the lognormalPDF function between a and b, by -//ten-point Gauss-Legendre integration, Numerical Recipe 2nd Ed. p.140 -inline double qgaus(scalar a, scalar b, scalar Nst_local) -{ - double dx, xm, xr, result; - static double w[] = { - 0.2955242247, - 0.2692667193, - 0.2190863625, - 0.1494513491, - 0.0666713443 - }; - static double x[] = { - 0.1488743389, - 0.4333953941, - 0.6794095682, - 0.8650633666, - 0.9739065285 - }; - - xm = 0.5*(b+a); - xr = 0.5*(b-a); - result = 0.0; - - for(label i=0 ; i<5 ; i++) - { - dx = xr*x[i]; - result = result + w[i]*(lognormalPDF(xm+dx, Nst_local) - + lognormalPDF(xm-dx, Nst_local)); - } - result = xr*result; - - return result; -} - -//Returns the P_N(Nst)*dNst -//Ref) S.H.Kim and Kang Y. Huh, Combust Flame, vol.138, pp.336-352, 2004 -// Equation 16 -inline void calculate_PdNst(scalar NumofNst, scalarField& Nst, scalar Nst_local, scalarField& result) -{ - const double alarge = 1000.0; - scalar a,b; - - if(Nst_local > 0.5*Nst[0]) - { - a = 0.0; - b = 0.5*(Nst[0]+Nst[1]); - result[0] = qgaus(a, b, Nst_local); - - for(label i=1 ; i<(NumofNst-1) ; i++) - { - a = 0.5*(Nst[i]+Nst[i-1]); - b = 0.5*(Nst[i]+Nst[i+1]); - result[i] = qgaus(a, b, Nst_local); - } - - a = 0.5*(Nst[NumofNst-1]+Nst[NumofNst-2]); - b = alarge; - result[NumofNst-1] = qgaus(a, b, Nst_local); - } - else - { - result[0] = 1.0; - for(label i=1 ; i=0 ; jj--) - { - if(eta_temp[jj] > eta[j]) - { - higher_eta = eta_temp[jj]; - hindex = jj; - } - } - - frac = (eta[j]-lower_eta)/(higher_eta-lower_eta); - - for(label m=0 ; mvalidate(); //SLFM + Info<<"Read SLFM library"< scalarFieldArray1d; typedef List scalarFieldArray2d; typedef List scalarFieldArray3d; -CMC::FlameStructure slfmLibrary(SLFMdict, etaGrid.size(), Y.size()); +SLFM::FlameStructure slfmLibrary(runTime.constant()/NstFolder, NstList, etaGrid, composition.species()); scalarFieldArray2d& Y_SLFM (slfmLibrary.Y()); scalarFieldArray2d& W_SLFM (slfmLibrary.W()); @@ -24,139 +24,6 @@ scalarFieldArray1d& Rgas_SLFM (slfmLibrary.Rgas()); for(label n=0 ; n>etamax_SLFM>>NoSpecies; - - //blank line - fin.getLine(gbg); - - //PRESSURE... - fin.getLine(gbg); - - //1.00... - fin.getLine(gbg); - - //MIXTURE... - fin.getLine(gbg); - - //read eta space (from SLFM library) - scalarField etaValue_SLFM(etamax_SLFM, 0.0); - for(label j=0 ; j>etaValue_SLFM[j]; - } - - //blank line - fin.getLine(gbg); - - //INITIAL... - fin.getLine(gbg); - - //species loop - scalarField Y_temp(etamax_SLFM, 0.0); - scalarField W_temp(etamax_SLFM, 0.0); - for(label i=0 ; i> spName; - const label spI (composition.species()[spName]); - // Info<> Yj; - Y_temp[j] = max(0.0, Yj); - } - - Y_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp); - - //read reaction rate (from SLFM library) - for(label j=0 ; j>W_temp[j]; - } - - W_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp); - - //blank line - fin.getLine(gbg); - } - - //read temperature - { - //H2... (species name) - fin.getLine(gbg); - - //read temperature - for(label j=0 ; j> Yj; - Y_temp[j] = max(0.0, Yj); - } - - T_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp); - - //read heat source - for(label j=0 ; j>W_temp[j]; - } - - Q_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp); - - //blank line - fin.getLine(gbg); - } - - //DENSITY... - fin.getLine(gbg); - - //read density (from SLFM library) - scalarField rho_temp(etamax_SLFM, 0.0); - for(label j=0 ; j>rho_temp[j]; - } - //Info<>h_temp[j]; - } - - h_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, h_temp); - for(label i=0 ; i