387 lines
8.8 KiB
C
387 lines
8.8 KiB
C
const int Ysize = Y.size();
|
|
|
|
//mixture fraction
|
|
Info<<"Reading field mf\n"<<endl;
|
|
volScalarField mf
|
|
(
|
|
IOobject
|
|
(
|
|
"mf",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::MUST_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh
|
|
);
|
|
|
|
//mixture fraction variance
|
|
Info<<"Reading field mfVar\n"<<endl;
|
|
volScalarField mfVar
|
|
(
|
|
IOobject
|
|
(
|
|
"mfVar",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::MUST_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh
|
|
);
|
|
|
|
//Scalar Dissipation Rate
|
|
volScalarField SDR
|
|
(
|
|
IOobject
|
|
(
|
|
"SDR",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("SDR", dimEnergy/dimTime/dimEnergy, 0.0)
|
|
);
|
|
|
|
//delta function index
|
|
volScalarField deltaftn
|
|
(
|
|
IOobject
|
|
(
|
|
"deltaftn",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("deltaftn", dimless, -1.0)
|
|
);
|
|
|
|
IOdictionary CMCdict
|
|
(
|
|
IOobject
|
|
(
|
|
"CMCdict",
|
|
runTime.constant(),
|
|
mesh,
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE
|
|
)
|
|
);
|
|
|
|
Switch uni_eta( CMCdict.lookup("uni_eta", false) );
|
|
Switch init_start_CMC( CMCdict.lookup("init_start_CMC", true) );
|
|
Switch multipleFG( CMCdict.lookup("multipleFG", false));
|
|
word latestCMCresults(CMCdict.lookup("latestCMCresultsFileName"));
|
|
word SLFMinit(CMCdict.lookup("SLFMinit_FileName"));
|
|
word fuelName(CMCdict.lookup("fuel_name"));
|
|
word CMCresults_Time(CMCdict.lookup("CMCresults_Time"));
|
|
scalar previousnf(readScalar(CMCdict.lookup("previousnf")));
|
|
scalar fuelperFG(readScalar(CMCdict.lookup("fuelperFG")));
|
|
label nf(readLabel(CMCdict.lookup("nf")));
|
|
|
|
scalar etamax( readScalar(CMCdict.lookup("etamax")) );
|
|
scalarField etaValue(etamax+1, 0.0);
|
|
scalarField deta(etamax+1 ,0.0 );
|
|
|
|
scalar section(readScalar(CMCdict.lookup("section")));
|
|
scalarField MFcut(section+1, 0.0);
|
|
scalarField Neta(section+1, 0.0);
|
|
|
|
scalarField pdf(etamax+1,0.0);
|
|
scalarField f(etamax+1, 0.0);
|
|
Info<<"Set eta space grid"<<tab<<uni_eta<<endl;
|
|
|
|
//Uniform eta space grid
|
|
if(uni_eta == true)
|
|
{
|
|
for(label etacount = 0 ; etacount <= etamax ; etacount++)
|
|
{
|
|
deta[etacount] = 1/etamax;
|
|
if(etacount == 0)
|
|
{
|
|
etaValue[0] = 0;
|
|
}
|
|
else
|
|
{
|
|
etaValue[etacount] = etaValue[etacount-1] + deta[etacount];
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
else if(uni_eta == false)
|
|
{
|
|
dictionary non_uni_eta_dict = CMCdict.subDict("non_uni_eta_dict");
|
|
|
|
MFcut = non_uni_eta_dict.lookup("MFcut");
|
|
Neta = non_uni_eta_dict.lookup("Neta");
|
|
|
|
int i = 0;
|
|
for(label etacount=1 ; etacount <= etamax ; etacount++)
|
|
{
|
|
if(etacount > Neta[i+1])
|
|
{
|
|
i = i+1;
|
|
}
|
|
if(Neta[i] < etacount && etacount <= Neta[i+1])
|
|
{
|
|
etaValue[etacount] = etaValue[etacount-1] + (MFcut[i+1]-MFcut[i])/(Neta[i+1] - Neta[i]);
|
|
}
|
|
}
|
|
|
|
for(label etacount=0 ; etacount <= etamax ; etacount++)
|
|
{
|
|
if(etacount == 0 || etacount == etamax)
|
|
{
|
|
deta[0] = etaValue[1] - etaValue[0];
|
|
deta[etamax] = etaValue[etamax] - etaValue[etamax-1];
|
|
}
|
|
else
|
|
{
|
|
deta[etacount] = ( etaValue[etacount+1] - etaValue[etacount-1] ) / 2;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
scalarField ChemDeltaT( (etamax+1)*nf, 0.0000001); //initial chemical time step;
|
|
label group(0); //flame group number to add evaporation source (default = 0), Flamegroup.H line# 29
|
|
|
|
|
|
//Peta field
|
|
PtrList<volScalarField> Peta(etamax+1); //Parameter List, number of eta space
|
|
for(label etacount = 0 ; etacount <= etamax ; etacount++)
|
|
{
|
|
word Petai = "Peta_"+Foam::name(etacount);
|
|
Peta.set (etacount,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
Petai,
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("Peta", dimless, 0.0)
|
|
)
|
|
);
|
|
}
|
|
|
|
wordList FGtypes
|
|
(
|
|
mf.boundaryField().size(),
|
|
zeroGradientFvPatchScalarField::typeName
|
|
);
|
|
|
|
forAll(mf.boundaryField(), patchi)
|
|
{
|
|
FGtypes[patchi] = mf.boundaryField()[patchi].type();
|
|
}
|
|
|
|
PtrList<volScalarField> F(nf); //flame group field
|
|
for(label k = 0 ; k < nf ; k++)
|
|
{
|
|
word Fk = "F_"+Foam::name(k);
|
|
F.set (k,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
Fk,
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::READ_IF_PRESENT,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("F", dimless, 0.0),
|
|
FGtypes
|
|
)
|
|
);
|
|
}
|
|
|
|
volScalarField F_total //summation of F[i] for each cell
|
|
(
|
|
IOobject
|
|
(
|
|
"F_total",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::READ_IF_PRESENT,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("F_total", dimless, 0.0)
|
|
);
|
|
|
|
PtrList<volScalarField> Ffrac(nf);
|
|
for(label k = 0 ; k < nf ; k++)
|
|
{
|
|
word Ffrack = "Ffrac_"+Foam::name(k);
|
|
Ffrac.set (k,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
Ffrack,
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("Ffrac", dimless, 0.0)
|
|
)
|
|
);
|
|
}
|
|
|
|
//conditional SDR
|
|
PtrList<volScalarField> CSDRCELL( etamax+1 ); //Parameter List, number of eta space
|
|
for(label j = 0 ; j <= etamax ; j++)
|
|
{
|
|
word CSDRCELLi = "CSDRCELL_"+Foam::name(j);
|
|
CSDRCELL.set (j,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
CSDRCELLi,
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("CSDRCELL", dimless, 0.0)
|
|
)
|
|
);
|
|
}
|
|
|
|
PtrList<volScalarField> SUMRR(Ysize); //Parameter List, number of eta space
|
|
for(label i = 0 ; i < Ysize ; i++)
|
|
{
|
|
word SUMRRi = "SUMRR_"+Foam::name(i);
|
|
SUMRR.set (i,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
SUMRRi,
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("SUMRR", dimless, 0.0)
|
|
)
|
|
);
|
|
}
|
|
|
|
//chemistry source term
|
|
volScalarField SUMSh
|
|
(
|
|
IOobject
|
|
(
|
|
"SUMSh",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("SUMSh", dimEnergy/dimTime/dimVolume, 0.0)
|
|
);
|
|
|
|
// Add eta-space for species Yi using vectorized array method (2D or 3D array -> 1D array)
|
|
// Qi(i , j) = Qi[i+j*(imax+1)]
|
|
// Qi(i , j , k) = Qi[ (k*(jmax+1)+j)*(imax+1) + i ]
|
|
//when i = 0~imax, j=0~jmax, k=0~kmax
|
|
// i = number of species
|
|
// j = number of eta-space
|
|
// k = number of flame group
|
|
|
|
/////////////////Qi and Qh ZONE value define//////////////////////
|
|
scalarField QiCMC( Ysize*(etamax+1)*nf , 0 );
|
|
scalarField QhCMC( (etamax+1)*nf , 0 );
|
|
scalarField QiCMCNew( Ysize*(etamax+1)*nf , 0 );
|
|
scalarField QhCMCNew( (etamax+1)*nf , 0 );
|
|
scalarField SDRCMC((etamax+1)*nf,0);
|
|
|
|
/////////////////////Conditional T field////////////////////////////
|
|
scalarField TCMC( (etamax+1)*nf, 0);
|
|
|
|
/////////////////////Conditional density field////////////////////////////
|
|
scalarField rhoCMC( (etamax+1)*nf, 0);
|
|
|
|
/////////////////////CMC equation source term/////////////////////////
|
|
scalarField CMCRR( Ysize*(etamax+1)*nf , 0 );
|
|
scalarField CMCSh( (etamax+1)*nf , 0 );
|
|
scalar Pin(p.weightedAverage(mesh.V()).value());
|
|
scalar Pin_old(Pin);
|
|
|
|
scalar o2index(0);
|
|
scalar fuelindex(0);
|
|
|
|
for(label i=0;i<Ysize;i++)
|
|
{
|
|
if(Y[i].name()=="O2")
|
|
{
|
|
o2index=i;
|
|
}
|
|
if(Y[i].name()==fuelName)
|
|
{
|
|
fuelindex=i;
|
|
}
|
|
}
|
|
|
|
scalar outi(0);
|
|
|
|
OFstream CMCresultsFile
|
|
(
|
|
("CMCresults."+ runTime.timeName() +".dat")
|
|
);
|
|
|
|
//make header of Tecplot format
|
|
CMCresultsFile<<"TITLE = "<<"\""<<"Results of CMC model calculation (Qi, Qh, T, rho, SDR)"<<"\""<<endl;
|
|
CMCresultsFile<<"VARIABLES = "<<"\"eta\"";
|
|
CMCresultsFile<<", ";
|
|
for(label i = 0 ; i<Ysize ; i++)
|
|
{
|
|
CMCresultsFile<<"\""<<Y[i].name()<<"\", ";
|
|
}
|
|
CMCresultsFile<<"\"Qh\", \"T\", \"RHO\", \"SDR\"";
|
|
CMCresultsFile<<endl;
|
|
//make header of Tecplot format
|
|
|
|
volScalarField mf_temp
|
|
(
|
|
IOobject
|
|
(
|
|
"mf_temp",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("zero",dimless,0.0)
|
|
);
|
|
|
|
|
|
|
|
if(init_start_CMC == false)
|
|
{
|
|
#include "restartCMC.H"
|
|
#include "CMCwrite.H"
|
|
}
|
|
|
|
scalar NewFGindex(0);
|
|
|