401 lines
8.2 KiB
C++
401 lines
8.2 KiB
C++
#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();
|
|
}
|
|
}
|