255 lines
6.5 KiB
C++
255 lines
6.5 KiB
C++
class BetaPDF
|
|
{
|
|
|
|
//private variable
|
|
|
|
//beta-pdf parameter alpha and beta
|
|
scalar alpha_, beta_;
|
|
|
|
//cutting point and number of eta-space
|
|
scalarField etaCut_, N_;
|
|
|
|
//detailed integration space
|
|
scalarField etaSpace_;
|
|
|
|
//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.
|
|
Info<<"Construct Beta-PDF"<<endl;
|
|
|
|
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_);
|
|
}
|
|
|
|
// Destructor
|
|
~BetaPDF(){}
|
|
|
|
// Member functions
|
|
|
|
// set beta-pdf parameter alpha_ and beta_
|
|
void setParameter(scalar& mf, scalar& 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;
|
|
}
|
|
}
|
|
}
|
|
// beta-pdf weighted integration for given mf, mfVar and f
|
|
scalar evaluate(scalar& mf, scalar& mfVar, scalarField& etaValue, scalarField& f)
|
|
{
|
|
scalar result(0);
|
|
|
|
setParameter(mf, mfVar);
|
|
|
|
if(fdelta_ == false)
|
|
{
|
|
limit_ab();
|
|
|
|
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;
|
|
}
|
|
scalar value(scalar& mf, scalar& mfVar, scalar& etaValue)
|
|
{
|
|
scalar result(0);
|
|
|
|
setParameter(mf, mfVar);
|
|
|
|
if(fdelta_ == false)
|
|
{
|
|
limit_ab();
|
|
|
|
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(scalar& a, scalar& b, scalarField& eta)
|
|
{
|
|
return pow(eta, a-1.0)*pow(1.0-eta, b-1.0);
|
|
}
|
|
scalar etaFunc(scalar& a, scalar& b, scalar& eta)
|
|
{
|
|
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(scalar& xl, scalar& xh, scalar& N, scalarField& fx)
|
|
{
|
|
scalar evensum(0.0), oddsum(0.0), sum(0.0);
|
|
|
|
scalar h = (xh - xl)/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(scalarField& eta)
|
|
{
|
|
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(scalar& eta)
|
|
{
|
|
return interpolateXY(eta, etaSpace_, AMCfine_);
|
|
}
|
|
void C1coeff(scalar& mf, 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);
|
|
}
|
|
}
|
|
};
|
|
|