research/OF-POSTECH-1
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

SLFMFoam beta pdf related classes

Browse source
This commit is contained in:
ignis 2017-09-03 04:16:50 +09:00
parent 7b6c05c25c
commit 48e2b161d7
26 changed files with 1651 additions and 947 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

68
solvers_post/AMC/AMC.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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::scalarField> Foam::CMC::AMC(const scalarField& eta)
{
tmp<scalarField> tRes(new scalarField(eta.size(), 0.0));
scalarField& Res = tRes.ref();
for(label i=0 ; i<eta.size() ; i++)
{
Res[i] = AMC(eta[i]);
}
return tRes;
}
// ************************************************************************* //

72
solvers_post/AMC/AMC.H Normal file
View file

@ -0,0 +1,72 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
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<scalarField> 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
// ************************************************************************* //

120
solvers_post/BetaFunction/BetaFunction.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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 * * * * * * * * * * * * * * //
// ************************************************************************* //

198
solvers_post/BetaFunction/BetaFunction.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
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<scalar>& 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
// ************************************************************************* //

58
solvers_post/BetaFunction/BetaFunctionI.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

65
solvers_post/BetaFunction/BetaFunctionIO.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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;
}
// ************************************************************************* //

126
solvers_post/BetaGrid/BetaGrid.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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 * * * * * * * * * * * * * * //
// ************************************************************************* //

161
solvers_post/BetaGrid/BetaGrid.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
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
// ************************************************************************* //

58
solvers_post/BetaGrid/BetaGridI.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

65
solvers_post/BetaGrid/BetaGridIO.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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;
}
// ************************************************************************* //

161
solvers_post/BetaIntegrator/BetaIntegrator.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<scalar>& fx)
const
{
scalar evensum(0.0), oddsum(0.0), sum(0.0);
scalar h = (xh - xl)/scalar(n);
for(label i=0 ; i<fx.size() ; i++)
{
if(i%2 == 0)
{
evensum += fx[i];
}
else
{
oddsum += fx[i];
}
}
sum = -1.0*fx.first() + 2.0*evensum + 4.0*oddsum - 1.0*fx.last();
return sum*(h/3.0);
}
Foam::scalar Foam::CMC::BetaIntegrator::sumSimps(const UList<scalar>& 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 * * * * * * * * * * * * * * //
// ************************************************************************* //

166
solvers_post/BetaIntegrator/BetaIntegrator.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
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<scalar>& fx) const;
//- Calculate integral using Simpson's rule
scalar sumSimps(const UList<scalar>& 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
// ************************************************************************* //

58
solvers_post/BetaIntegrator/BetaIntegratorI.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

56
solvers_post/BetaIntegrator/BetaIntegratorIO.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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;
}
// ************************************************************************* //

312
solvers_post/SLFMFoam/BetaPDF.H
View file

@ -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<scalarField> 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<etaCut_.size()-1 ; i++)
{
scalarField tmp;
tmp.append(etaCut_[i]);
del = (etaCut_[i+1] - etaCut_[i])/N_[i];
for(label j=0 ; j<N_[i] ; j++)
{
etaSpace_.append(etaSpace_[cnt] + del);
tmp.append(etaSpace_[cnt]+del);
cnt++;
}
etaPart_[i].append(tmp);
}
AMCfine_ = AMC(etaSpace_);
pdfNum_ = scalarField(etaSpace_.size(), 0.0);
pdfDen_ = 1.0;
}
// Destructor
~BetaPDF(){}
// Member functions
// set beta-pdf parameter alpha_ and beta_
void setParameter(const scalar mf, const scalar mfVar)
{
mf_ = mf;
mfVar_ = mfVar;
fdelta_ = false;
delta_ox = false;
delta_fu = false;
scalar gamma = mf*(1.0-mf)/(mfVar+SMALL) - 1.0;
if(gamma<SMALL || (mfVar)/(mf*(1-mf)+SMALL)<0.001)
{
fdelta_ = true;
}
else
{
alpha_ = max(mf*gamma, 0.0);
beta_ = max((1.0-mf)*gamma, 0.0);
if(alpha_ < 1.0)
{
delta_ox = true;
}
if(beta_ < 1.0)
{
delta_fu = true;
}
limit_ab();
}
pdfNum_ = etaFunc(alpha_, beta_, etaSpace_);
pdfDen_ = integrate(pdfNum_)
+ Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
+ Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
}
// beta-pdf weighted integration for given mf, mfVar and f
scalar evaluate(const scalar mf, const scalar mfVar, const scalarField& etaValue, const scalarField& f)
{
scalar result(0);
setParameter(mf, mfVar);
if(fdelta_ == false)
{
scalar num(0), den(0);
for(label i=0 ; i<etaCut_.size()-1 ; i++)
{
scalarField fPart = interpolateXY(etaPart_[i], etaValue, f);
scalarField fx = fPart * etaFunc(alpha_, beta_, etaPart_[i]);
scalarField gx = etaFunc(alpha_, beta_, etaPart_[i]);
num += simps(etaCut_[i], etaCut_[i+1], N_[i], fx);
den += simps(etaCut_[i], etaCut_[i+1], N_[i], gx);
}
num += f[0]*Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
den += Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
num += f[etaValue.size()-1]*Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
den += Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
result = num/den;
}
else
{
result = interpolateXY(mf, etaValue, f);
}
return result;
}
// beta-pdf weighted integration for given mf, mfVar and f
scalar betaIntegrate(const scalarField& f) const
{
scalar result = 0.0;
if(fdelta_)
{
result = interpolateXY(mf_, etaSpace_, f);
}
else
{
scalar num = integrate (f * pdfNum_)
+ f.first()*Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
+ f.last()*Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
result = num/pdfDen_;
}
return result;
}
scalar integrate(const scalarField& f) const
{
scalar total = 0.0;
forAll(N_, i)
{
const labelList::subList prev(N_, i);
const label baseI = sum(prev);
const scalarField::subField subInterval(f, N_[i]+1, baseI);
total += simps(etaCut_[i], etaCut_[i+1], N_[i], subInterval);
}
return total;
}
scalar value(const scalar mf, const scalar mfVar, const scalar etaValue)
{
scalar result(0);
setParameter(mf, mfVar);
if(fdelta_ == false)
{
scalar den(0);
for(label i=0 ; i<etaCut_.size()-1 ; i++)
{
scalarField gx = etaFunc(alpha_, beta_, etaPart_[i]);
den += simps(etaCut_[i], etaCut_[i+1], N_[i], gx);
}
den += Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL);
den += Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
scalarField pdf = etaFunc(alpha_, beta_, etaSpace_)/den;
result = interpolateXY(etaValue, etaSpace_, pdf);
}
return result;
}
void limit_ab()
{
scalar fmax = 1.0+(beta_-1.0)/(alpha_-1.0);
fmax = 1.0/fmax;
if(alpha_ > 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<scalar>& fx) const
{
scalar evensum(0.0), oddsum(0.0), sum(0.0);
scalar h = (xh - xl)/scalar(N);
for(label i=0 ; i<fx.size() ; i++)
{
if(i%2 == 0)
{
evensum += fx[i];
}
else
{
oddsum += fx[i];
}
}
sum = -1.0*fx.first() + 2.0*evensum + 4.0*oddsum - 1.0*fx.last();
return sum*(h/3.0);
}
//Amplitude Mapping Closure (from KIVA)
//Define exp(-2*(erf^-1(2*eta - 1))^2)
scalarField AMC(const scalarField& eta) const
{
const scalar pi = 3.141592;
const scalar spi = Foam::sqrt(pi);
scalarField result(eta.size(), 0.0);
for(label i=0 ; i<eta.size() ; i++)
{
scalar a1=0.5, a0=(2.0*eta[i]-1.0), slope, 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;
}
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<varValue.size()-1 ; v++)
{
scalar var = maxVar*varValue[v];
C1table[v] = 1.0/evaluate(mf,var,x,fx);
}
}
};

193
solvers_post/SLFMFoam/FlameStructure/FlameStructure.C
View file

@ -25,13 +25,10 @@ License
#include "FlameStructure.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// const Foam::CMC::dataType Foam::CMC::FlameStructure::staticData();
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
#include "fileName.H"
#include "IFstream.H"
#include "hashedWordList.H"
#include "interpolateXY.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -41,42 +38,156 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::CMC::FlameStructure::FlameStructure(const dictionary &dict, const label nEta, const label nY)
Foam::SLFM::FlameStructure::FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList &spTable)
:
// baseClassName(),
NstDir_(dict.lookup("NstFolder")),
NstList_(dict.lookup("NstList")),
Y_(nY, scalarFieldArray1d(NstList_.size(), scalarField(nEta, 0.0))),
W_(nY, scalarFieldArray1d(NstList_.size(), scalarField(nEta, 0.0))),
T_(NstList_.size(), scalarField(nEta, 0.0)),
Q_(NstList_.size(), scalarField(nEta, 0.0)),
h_(NstList_.size(), scalarField(nEta, 0.0)),
rho_(NstList_.size(), scalarField(nEta, 0.0)),
Rgas_(NstList_.size(), scalarField(nEta, 0.0))
{}
Foam::CMC::FlameStructure::FlameStructure(const FlameStructure&)
// :
// baseClassName(),
// data_()
{}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
/*
Foam::autoPtr<Foam::CMC::FlameStructure>
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<FlameStructure>(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<etamax_SLFM ; j++)
{
scalar Yj = 0.0;
fin >> 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<<rho_temp<<endl;
rho_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
//blank line
fin.getLine(gbg);
//ENTHALPY...
fin.getLine(gbg);
//read enthalpy (from SLFM library)
for(label j = 0 ; j < etamax_SLFM ; j++)
{
fin>>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 * * * * * * * * * * * * * * //

35
solvers_post/SLFMFoam/FlameStructure/FlameStructure.H
View file

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

14
solvers_post/SLFMFoam/FlameStructure/FlameStructureIO.C
View file

@ -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;

12
solvers_post/SLFMFoam/Make/files
View file

@ -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

4
solvers_post/SLFMFoam/Make/options
View file

@ -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 \

401
solvers_post/SLFMFoam/Math.H
View file

@ -1,401 +0,0 @@
#include <cmath>
//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<NumofNst ; i++)
{
result[i] = 0.0;
}
}
}
inline void Interpolation(scalarField& eta, scalarField& eta_temp, scalar NumofNst, scalar Nspecies, scalarField& Y_SLFM_temp, scalarField& result)
{
scalar eta_size(eta.size()), eta_temp_size(eta_temp.size());
scalar lower_eta(0), higher_eta(0);
scalar lindex(0), hindex(0);
scalar frac(0), diff(0);
for(label j=0; j<eta_size ; j++)
{
for(label jj=0; jj<eta_temp_size ; jj++)
{
if(eta_temp[jj] <= eta[j])
{
lower_eta = eta_temp[jj];
lindex = jj;
}
}
for(label jj=eta_temp_size-1 ; jj>=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 ; m<NumofNst ; m++)
{
for(label i=0 ; i<Nspecies ; i++)
{
diff = Y_SLFM_temp[(m*eta_temp_size+hindex)*Nspecies+i]
- Y_SLFM_temp[(m*eta_temp_size+lindex)*Nspecies+i];
result[(m*eta_size+j)*Nspecies+i] = Y_SLFM_temp[(m*eta_temp_size+lindex)*Nspecies+i] + frac*diff;
}
}
}
}
inline void select_ij(label ij, scalar& i, scalar& j)
{
if(ij == 0) //variance for H, O2, CO, OH, CH4, T
{
i = 0;
j = 0;
}
else if(ij == 1)
{
i = 1;
j = 1;
}
else if(ij == 2)
{
i = 2;
j = 2;
}
else if(ij == 3)
{
i = 3;
j = 3;
}
else if(ij == 4)
{
i = 4;
j = 4;
}
else if(ij == 5)
{
i = 5;
j = 5;
}
else if(ij == 6) //covariance between H and O2
{
i = 0;
j = 1;
}
else if(ij == 7) //covariance between CO and OH
{
i = 2;
j = 3;
}
else if(ij == 8) //covariance between H and CH4
{
i = 0;
j = 4;
}
else if(ij == 9) //covariance between SPECIES and T
{
i = 0;
j = 5;
}
else if(ij == 10)
{
i = 1;
j = 5;
}
else if(ij == 11)
{
i = 2;
j = 5;
}
else if(ij == 12)
{
i = 3;
j = 5;
}
else if(ij == 13)
{
i = 4;
j = 5;
}
else
{
abort();
}
}
inline void select_species(label i, word& speciei)
{
if(i == 0)
{
speciei = "H";
}
else if(i == 1)
{
speciei = "O2";
}
else if(i == 2)
{
speciei = "CO";
}
else if(i == 3)
{
speciei = "OH";
}
else if(i == 4)
{
speciei = "CH4";
}
}
inline scalar select_CiCj(scalar i)
{
if(i == 1 || i == 2 || i == 4 || i == 5)
{
return 1.0;
}
else if(i == 0 || i == 3)
{
return 2.0;
}
else
{
abort();
}
}

9
solvers_post/SLFMFoam/SLFMFoam.C
View file

@ -47,6 +47,11 @@ Contact
#include "simpleControl.H"
#include "FlameStructure/FlameStructure.H"
#include "interpolateXY.H"
#include "AMC.H"
#include "BetaFunction.H"
#include "BetaGrid.H"
#include "BetaIntegrator.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@ -59,7 +64,6 @@ int main(int argc, char *argv[])
#include "createControl.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "BetaPDF.H" //make detailed integration space and perform beta-pdf integration
#include "createFvOptions.H"
#include "initContinuityErrs.H"
@ -67,9 +71,10 @@ int main(int argc, char *argv[])
turbulence->validate();
//SLFM
Info<<"Read SLFM library"<<endl;
#include "readSLFMlib.H" //read SLFM library
Info<<"Construct Beta-PDF"<<endl;
BetaPDF bpdf(SLFMdict);
#include "preIntegration.H" //AMC C1 coefficient pdf integration
Info<< "\nStarting time loop\n" << endl;

25
solvers_post/SLFMFoam/createFields.H
View file

@ -291,7 +291,6 @@ IOdictionary SLFMdict //SLFM properties dictionary
)
);
const label Ysize = Y.size();
const wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
labelList postSpecieIndices(postSpecieNames.size(),-1);
forAll(postSpecieNames, i)
@ -309,28 +308,10 @@ const label etamax(lowerN+upperN+1);
word outletName(SLFMdict.lookup("outletName"));
labelList nIntervals (SLFMdict.lookup("detailedN"));
scalarField etaGrid(sum(nIntervals)+1, 0.0);
{
scalarList etas (SLFMdict.lookup("detailedEta"));
forAll(nIntervals, i)
{
const label itv = nIntervals[i];
const labelList::subList prev(nIntervals, i);
const label baseI = sum(prev);
const scalar low = etas[i];
const scalar upp = etas[i+1];
const scalar delta = (upp - low) / scalar(itv);
for (label j = 0; j < itv; j++)
{
etaGrid[baseI+j] = low + delta * j;
}
}
etaGrid.last() = etas.last();
}
CMC::BetaGrid bg(SLFMdict);
const scalarField &etaGrid(bg.etaSpace());
const scalarField amcGridValue(CMC::AMC(etaGrid));
Info << etaGrid << endl;

24
solvers_post/SLFMFoam/preIntegration.H
View file

@ -39,28 +39,38 @@ scalarFieldArray3d Ytable
);
Info<<"Construct C1 coefficient table"<<endl;
for(label j=1 ; j<etamax-1; j++)
for(label j = 1; j < etamax - 1; j++)
{
bpdf.C1coeff(etaValue[j], varValue, C1table[j]);
for(label v = 1; v < varValue.size() - 1; v++)
{
const scalar mf = etaValue[j];
const scalar mfVar = varValue[v]*mf*(1.0-mf);
CMC::BetaFunction bf(mf, mfVar);
CMC::BetaIntegrator bi(bf, bg);
C1table[j][v] = 1.0/bi.betaIntegrate(amcGridValue);
}
}
Info<<"Construct Temp and Rgas table"<<endl;
for(label j=0 ; j<etamax ; j++)
{
const scalar mf = etaValue[j];
for(label v=0 ; v<=NVar ; v++)
{
const scalar mf = etaValue[j];
const scalar mfVar = varValue[v]*mf*(1.0-mf);
bpdf.setParameter(mf, mfVar);
CMC::BetaFunction bf(mf, mfVar);
CMC::BetaIntegrator bi(bf, bg);
for(label n=0 ; n<NstList.size() ; n++)
{
Ttable[j][v][n] = bpdf.betaIntegrate(T_SLFM[n]);
Rtable[j][v][n] = bpdf.betaIntegrate(Rgas_SLFM[n]);
Ttable[j][v][n] = bi.betaIntegrate(T_SLFM[n]);
Rtable[j][v][n] = bi.betaIntegrate(Rgas_SLFM[n]);
forAll(Y, yi)
{
Ytable[yi][j][v][n] = bpdf.betaIntegrate(Y_SLFM[yi][n]);
Ytable[yi][j][v][n] = bi.betaIntegrate(Y_SLFM[yi][n]);
}
}
}

135
solvers_post/SLFMFoam/readSLFMlib.H
View file

@ -12,7 +12,7 @@ typedef List<scalarField> scalarFieldArray1d;
typedef List<scalarFieldArray1d> scalarFieldArray2d;
typedef List<scalarFieldArray2d> 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<NstList.size() ; n++)
{
scalar N = NstList[n];
fileName fname = "sdr"+Foam::name(N)+".inp";
IFstream fin(runTime.constant()/NstFolder/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 Y_temp(etamax_SLFM, 0.0);
scalarField W_temp(etamax_SLFM, 0.0);
for(label i=0 ; i<NoSpecies ; i++)
{
//H2... (species name)
string spLine;
fin.getLine(spLine);
IStringStream spStream(spLine);
word spName;
spStream >> spName;
const label spI (composition.species()[spName]);
// Info<<spI << spName << endl;
//read species mass fraction (from SLFM library)
for(label j=0 ; j<etamax_SLFM ; j++)
{
scalar Yj = 0.0;
fin >> 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<etamax_SLFM ; j++)
{
fin>>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<etamax_SLFM ; j++)
{
scalar Yj = 0.0;
fin >> 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<etamax_SLFM ; j++)
{
fin>>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<etamax_SLFM ; j++)
{
fin>>rho_temp[j];
}
//Info<<rho_temp<<endl;
rho_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, rho_temp);
//blank line
fin.getLine(gbg);
//ENTHALPY...
fin.getLine(gbg);
//read enthalpy (from SLFM library)
scalarField h_temp(etamax_SLFM, 0.0);
for(label j=0 ; j<etamax_SLFM ; j++)
{
fin>>h_temp[j];
}
h_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, h_temp);
for(label i=0 ; i<Y.size() ; i++)
{
Rgas_SLFM[n] += Y_SLFM[i][n]/composition.W(i); // = 1/meanMW

2
solvers_post/SLFMFoam/setSDR.H
View file

@ -1,4 +1,4 @@
scalarField AMCcoeff(bpdf.AMC(etaValue));
scalarField AMCcoeff(CMC::AMC(etaValue));
//cell loop