research/OpenFOAM-4.x-lab
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

removed thermal NOx solver and add fuel NOx switch to thermal-fuel NOx solver

Browse source
This commit is contained in:
ignis 2018-02-01 14:37:28 +09:00
parent bfbb887274
commit 5ac3bc47ea
18 changed files with 241 additions and 732 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
applications/solvers/combustion/NOxFoam_thermalNOx/Make/files
View file

@ -1,4 +0,0 @@
NOxFoam.C
EXE = $(FOAM_USER_APPBIN)/NOxFoam_thermalNOx

10
applications/solvers/combustion/NOxFoam_thermalNOx/Make/options
View file

@ -1,10 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-lsampling

97
applications/solvers/combustion/NOxFoam_thermalNOx/NOxFoam.C
View file

@ -1,97 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-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/>.
Application
NOxFoam_thermalNOx
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFvOptions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
if(instantaneousRadicals==false)
{
Info<< "Partial Equilibrium Approach is selected for O, OH radicals" << endl;
}
else
{
Info<< "Instantaneous massfraction field will be used for O, OH radicals" << endl;
}
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "SnoCalc.H"
mut=rho*nut;
while (simple.correctNonOrthogonal())
{
fvScalarMatrix NOEqn
(
fvm::ddt(rho, NO)
+ fvm::div(phi, NO)
- fvm::laplacian(mut, NO)
==
Sno
);
NOEqn.relax();
fvOptions.constrain(NOEqn);
NOEqn.solve();
fvOptions.correct(NO);
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

44
applications/solvers/combustion/NOxFoam_thermalNOx/SnoCalc.H
View file

@ -1,44 +0,0 @@
const scalar Wno = 0.030006; //Molecular weight of NO [kg/mol]
const scalar Wn2 = 0.028013; // N2 [kg/mol]
const scalar Wo2 = 0.031998; // O2 [kg/mol]
const scalar Wh2o = 0.018015; // H2O[kg/mol]
const scalar Wo = 0.015999; // O [kg/mol]
const scalar Woh = 0.017008; // OH [kg/mol]
forAll (mesh.cells(),celli)
{
Xno[celli]=rho[celli]*NO[celli]/Wno; //Molar concentration of NO [mol/m^3]
Xn2[celli]=rho[celli]*N2[celli]/Wn2; // N2 [mol/m^3]
Xo2[celli]=rho[celli]*O2[celli]/Wo2; // O2 [mol/m^3]
Xh2o[celli]=rho[celli]*H2O[celli]/Wh2o; // H2O[mol/m^3]
// O,OH molar concentrations
if(instantaneousRadicals==false)
{
Xo[celli]=36.64*Foam::pow(T[celli],0.5)*Foam::pow(Xo2[celli],0.5)*Foam::exp(-27123/T[celli]);
//Molar concentration of O2 [mol/m^3], Partial Equilibrium Approach
Xoh[celli]=2.129E+02*Foam::pow(T[celli],-0.57)*Foam::exp(-4595/T[celli])*Foam::pow(Xo[celli],0.5)*Foam::pow(Xh2o[celli],0.5);
//Molar concentration of OH [mol/m^3], Partial Equilibrium Approach
}
else
{
Xo[celli]=rho[celli]*O[celli]/Wo; //Molar concentration of O [mol/m^3]
Xoh[celli]=rho[celli]*OH[celli]/Woh; // OH[mol/m^3]
}
//Reaction constants
kf1[celli]=1.8E+08*Foam::exp(-38370/T[celli]);
kf2[celli]=1.8E+04*T[celli]*Foam::exp(-4680/T[celli]);
kf3[celli]=7.1E+07*Foam::exp(-450/T[celli]);
kr1[celli]=3.8E+07*Foam::exp(-425/T[celli]);
kr2[celli]=3.81E+03*T[celli]*Foam::exp(-20820/T[celli]);
Sno[celli] = Wno*2*kf1[celli]*Xo[celli]*Xn2[celli]*
((1-(kr1[celli]*kr2[celli]*Xno[celli]*Xno[celli])/(kf1[celli]*Xn2[celli]*kf2[celli]*Xo2[celli]))/
(1+(kr1[celli]*Xno[celli])/(kf2[celli]*Xo2[celli]+kf3[celli]*Xoh[celli])))/runTime.time().deltaTValue();
}

340
applications/solvers/combustion/NOxFoam_thermalNOx/createFields.H
View file

@ -1,340 +0,0 @@
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rho*U) & mesh.Sf()
);
volScalarField nut
(
IOobject
(
"nut",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField mut
(
IOobject
(
"mut",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("mut",dimensionSet(1,-1,-1,0,0,0,0),0.0)
);
volScalarField Xno
(
IOobject
(
"Xno",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xno",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField Xn2
(
IOobject
(
"Xn2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xn2",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField Xo2
(
IOobject
(
"Xo2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xo2",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField Xh2o
(
IOobject
(
"Xh2o",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xh2o",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField Xo
(
IOobject
(
"Xo",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xo",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField Xoh
(
IOobject
(
"Xoh",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xoh",dimensionSet(0,-3,0,0,1,0,0),0.0)
);
volScalarField NO
(
IOobject
(
"NO",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField N2
(
IOobject
(
"N2",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField O2
(
IOobject
(
"O2",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField H2O
(
IOobject
(
"H2O",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField O
(
IOobject
(
"O",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField OH
(
IOobject
(
"OH",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField Sno
(
IOobject
(
"Sno",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Sno",dimensionSet(1,-3,-1,0,0,0,0),0.0)
);
volScalarField kf1
(
IOobject
(
"kf1",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf1",dimensionSet(0,3,-1,0,-1,0,0),0.0)
);
volScalarField kf2
(
IOobject
(
"kf2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf2",dimensionSet(0,3,-1,0,-1,0,0),0.0)
);
volScalarField kf3
(
IOobject
(
"kf3",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf3",dimensionSet(0,3,-1,0,-1,0,0),0.0)
);
volScalarField kr1
(
IOobject
(
"kr1",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kr1",dimensionSet(0,3,-1,0,-1,0,0),0.0)
);
volScalarField kr2
(
IOobject
(
"kr2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kr2",dimensionSet(0,3,-1,0,-1,0,0),0.0)
);
IOdictionary modelParameter
(
IOobject
(
"modelParameter",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
bool instantaneousRadicals(modelParameter.lookupOrDefault("instantaneousRadicals", true));

2
applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files
View file

@ -1,5 +1,5 @@
NOxFoam.C
EXE = $(FOAM_APPBIN)/NOxFoam_thermFuel_2
EXE = $(FOAM_APPBIN)/NOxFoam

16
applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C
View file

@ -22,11 +22,9 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
simpleReactingParcelFoam
NOxFoam
Description
Steady state solver for compressible, turbulent flow with reacting,
multiphase particle clouds and optional sources/constraints.
\*---------------------------------------------------------------------------*/
@ -49,19 +47,19 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
if(instantaneousRadicals==false)
if(instantaneousRadicals)
{
Info<< "Instantaneous mass fraction field will be used for O, OH radicals" << endl;
}
else
{
Info<< "Partial Equilibrium Approach is selected for O, OH radicals" << endl;
}
else if(instantaneousRadicals==true)
{
Info<< "Instantaneous massfraction field will be used for O, OH radicals" << endl;
}
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
mut=rho*nut;
#include "SourceCalc.H"
#include "NH3Eqn.H"
#include "NOEqn.H"

114
applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H
View file

@ -1,94 +1,42 @@
const scalar Wno = 0.030006; //Molecular weight of NO [kg/mol]
const scalar Wn2 = 0.028013; // N2 [kg/mol]
const scalar Wo2 = 0.031998; // O2 [kg/mol]
const scalar Wh2o = 0.018015; // H2O[kg/mol]
const scalar Wo = 0.015999; // O [kg/mol]
const scalar Woh = 0.017008; // OH [kg/mol]
const scalar Wnh3 = 0.017031; // NH3[kg/mol]
const scalar Wn = 0.014007; // N [kg/mol]
const scalar RR = 8.314; //Universal gas constant [J/(mol*K)]
const scalar A1 = 4.0E+06; // [1/s]
const scalar A2 = 1.8E+08; // [1/s]
const scalar E1 = 133947.2; // [J/mol]
const scalar E2 = 113017.95; // [J/mol]
Xno = rho*NO/Wno; // Molar concentration of NO [mol/m^3]
// Calculation of source term. (Thermal NO)
forAll (mesh.cells(),celli)
{
// molar concentrations
Xno[celli]=rho[celli]*NO[celli]/Wno; //Molar concentration of NO [mol/m^3]
Xn2[celli]=rho[celli]*N2[celli]/Wn2; // N2 [mol/m^3]
Xo2[celli]=rho[celli]*O2[celli]/Wo2; // O2 [mol/m^3]
Xh2o[celli]=rho[celli]*H2O[celli]/Wh2o; // H2O[mol/m^3]
Xnh3[celli]=rho[celli]*NH3[celli]/Wnh3; // nh3[mol/m^3]
// O,OH molar concentrations
if(instantaneousRadicals==false)
{
Xo[celli]=36.64*Foam::pow(T[celli],0.5)*Foam::pow(Xo2[celli],0.5)*Foam::exp(-27123/T[celli]);
//Molar concentration of O2 [mol/m^3], Partial Equilibrium Approach
Xoh[celli]=2.129E+02*Foam::pow(T[celli],-0.57)*Foam::exp(-4595/T[celli])*Foam::pow(Xo[celli],0.5)*Foam::pow(Xh2o[celli],0.5);
//Molar concentration of OH [mol/m^3], Partial Equilibrium Approach
}
else
{
Xo[celli]=rho[celli]*O[celli]/Wo; //Molar concentration of O [mol/m^3]
Xoh[celli]=rho[celli]*OH[celli]/Woh; // OH[mol/m^3]
}
// NH3
//Oxygen Reaction Order, a
if (Xo2[celli] <= 4.1E-03) {
a[celli] = 1.0;
}
else if ( Xo2[celli] > 4.1E-03 and Xo2[celli] <= 1.11E-02 ) {
a[celli] = -3.95-0.9*Foam::log(Xo2[celli]);
}
else if ( Xo2[celli] > 1.11-02 and Xo2[celli] <= 0.03 ) {
a[celli] = -0.35-0.1*Foam::log(Xo2[celli]);
}
else if ( Xo2[celli] > 0.03 ) {
a[celli] = 0;
}
//Conversion rate of NH3
R1[celli] = A1*Xnh3[celli]*Foam::pow(Xo2[celli],a[celli])*Foam::exp(-E1/(RR*T[celli])); // [1/s]
R2[celli] = A2*Xnh3[celli]*Xno[celli]*Foam::exp(-E2/(RR*T[celli])); // [1/s]
//Source and sink of NH3
Snh3_p[celli] = (Sfuel1[celli]+Sfuel2[celli])*Yn[celli]*Wnh3/Wn/mesh.V()[celli]; //NH3 production
Snh3_1[celli] = -R1[celli]*Wnh3*p[celli]/(RR*T[celli]); //NH3 consumption -> NO (oxidation)
Snh3_2[celli] = -R2[celli]*Wnh3*p[celli]/(RR*T[celli]); //NH3 consumption -> N2 (reduction)
//Sum of Sources (NH3 production, consumption 1, and consumption 2 )
Snh3[celli] = Snh3_p[celli] +Snh3_1[celli] +Snh3_2[celli];
//Thermal NO
//Reaction rate constants [m^3/(mol*s)]
kf1[celli]=1.8E+08*Foam::exp(-38370/T[celli]);
kf2[celli]=1.8E+04*T[celli]*Foam::exp(-4680/T[celli]);
kf3[celli]=7.1E+07*Foam::exp(-450/T[celli]);
kr1[celli]=3.8E+07*Foam::exp(-425/T[celli]);
kr2[celli]=3.81E+03*T[celli]*Foam::exp(-20820/T[celli]);
//Calculation of source term. (Thermal NO)
SthermNO[celli] = Wno*2*kf1[celli]*Xo[celli]*Xn2[celli]*
SthermNO[celli] = Wno.value()*2*kf1[celli]*Xo[celli]*Xn2[celli]*
((1-(kr1[celli]*kr2[celli]*Xno[celli]*Xno[celli])/
(kf1[celli]*Xn2[celli]*kf2[celli]*Xo2[celli]))/
(1+(kr1[celli]*Xno[celli])/
(kf2[celli]*Xo2[celli]+kf3[celli]*Xoh[celli])))
/runTime.time().deltaTValue();
// Fuel NO
//Calculation of source term. (Fuel NO)
SfuelNO_1[celli] = R1[celli]*Wno*p[celli]/(RR*T[celli]); // NH3 + O2 -> NO (source)
SfuelNO_2[celli] =-R2[celli]*Wno*p[celli]/(RR*T[celli]); // NH3 + NO -> N2 (sink)
SfuelNO[celli] = SfuelNO_1[celli] + SfuelNO_2[celli]; // Sum of fuel NO source terms [kg/(m^3*s)]
// Sum of NO Sources (source of thermal NO and fuel NO)
Sno[celli] = SthermNO[celli] + SfuelNO[celli]; // [kg/(m^3*s)]
}
Sno = SthermNO;
if (calculateFuelNOx)
{
Xnh3 = rho*NH3/Wnh3; // Molar concentration of nh3[mol/m^3]
forAll (mesh.cells(),celli)
{
// Conversion rate of NH3
R1[celli] = A1*Xnh3[celli]*Foam::pow(Xo2[celli],a[celli])*Foam::exp(-E1/(RR*T[celli])); // [1/s]
R2[celli] = A2*Xnh3[celli]*Xno[celli]*Foam::exp(-E2/(RR*T[celli])); // [1/s]
}
// NH3 sinks
Snh3_1 = -R1*Wnh3*p/(Rgas*T); // NH3 consumption -> NO (oxidation)
Snh3_2 = -R2*Wnh3*p/(Rgas*T); // NH3 consumption -> N2 (reduction)
// net NH3 source (NH3 production, consumption 1, and consumption 2)
Snh3 = Snh3_p + Snh3_1 + Snh3_2;
// Fuel NO
// Calculation of source term. (Fuel NO)
SfuelNO_1 = R1*Wno*p/(Rgas*T); // NH3 + O2 -> NO (source)
SfuelNO_2 = -R2*Wno*p/(Rgas*T); // NH3 + NO -> N2 (sink)
SfuelNO = SfuelNO_1 + SfuelNO_2; // Sum of fuel NO source terms
// Sum of NO Sources (source of thermal NO and fuel NO)
Sno += SfuelNO;
}

289
applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H
View file

@ -1,4 +1,52 @@
IOdictionary NOxProperties
(
IOobject
(
"NOxProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
wordList fuelSourceNames(NOxProperties.lookup("fuelSourceNames"));
PtrList<volScalarField::Internal> fuelSources(fuelSourceNames.size());
forAll (fuelSources, si)
{
fuelSources.set
(
si,
new volScalarField::Internal
(
IOobject
(
fuelSourceNames[si],
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
);
}
scalarField nitrogenMassInFuels(NOxProperties.lookup("nitrogenMassInFuels"));
bool instantaneousRadicals
(
NOxProperties.lookupOrDefault("instantaneousRadicals", true)
);
bool calculateFuelNOx
(
NOxProperties.lookupOrDefault("calculateFuelNOx", false)
);
volScalarField T
(
IOobject
@ -12,6 +60,19 @@ volScalarField T
mesh
);
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField p
(
IOobject
@ -26,6 +87,7 @@ volScalarField p
);
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
@ -39,19 +101,7 @@ volVectorField U
mesh
);
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
@ -65,7 +115,6 @@ surfaceScalarField phi
linearInterpolate(rho*U) & mesh.Sf()
);
volScalarField nut
(
IOobject
@ -89,11 +138,9 @@ volScalarField mut
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("mut",dimensionSet(1,-1,-1,0,0,0,0),0.0)
rho * nut
);
volScalarField Xno
(
IOobject
@ -104,8 +151,8 @@ volScalarField Xno
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xno",dimensionSet(0,-3,0,0,1,0,0),0.0)
mesh,
dimensionedScalar("Xno",dimMoles/dimVolume,0.0)
);
volScalarField Xn2
@ -119,7 +166,7 @@ volScalarField Xn2
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xn2",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xn2",dimMoles/dimVolume,0.0)
);
volScalarField Xo2
@ -133,7 +180,7 @@ volScalarField Xo2
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xo2",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xo2",dimMoles/dimVolume,0.0)
);
volScalarField Xh2o
@ -147,7 +194,7 @@ volScalarField Xh2o
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xh2o",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xh2o",dimMoles/dimVolume,0.0)
);
volScalarField Xnh3
@ -161,7 +208,7 @@ volScalarField Xnh3
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xnh3",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xnh3",dimMoles/dimVolume,0.0)
);
volScalarField Xo
@ -175,7 +222,7 @@ volScalarField Xo
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xo",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xo",dimMoles/dimVolume,0.0)
);
volScalarField Xoh
@ -189,7 +236,7 @@ volScalarField Xoh
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Xoh",dimensionSet(0,-3,0,0,1,0,0),0.0)
dimensionedScalar("Xoh",dimMoles/dimVolume,0.0)
);
volScalarField NO
@ -277,12 +324,18 @@ volScalarField NH3
"NH3",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh
mesh,
dimensionedScalar("NH3",dimless,0.0)
);
if (calculateFuelNOx)
{
NH3.writeOpt() = IOobject::AUTO_WRITE;
}
volScalarField Sno
(
IOobject
@ -290,11 +343,11 @@ volScalarField Sno
"Sno",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Sno",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("Sno",dimDensity/dimTime,0.0)
);
volScalarField SthermNO
@ -304,11 +357,11 @@ volScalarField SthermNO
"SthermNO",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("SthermNO",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("SthermNO",dimDensity/dimTime,0.0)
);
volScalarField SfuelNO
@ -322,7 +375,7 @@ volScalarField SfuelNO
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("SfuelNO",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("SfuelNO",dimDensity/dimTime,0.0)
);
volScalarField SfuelNO_1
@ -332,11 +385,11 @@ volScalarField SfuelNO_1
"SfuelNO_1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("SfuelNO_1",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("SfuelNO_1",dimDensity/dimTime,0.0)
);
volScalarField SfuelNO_2
@ -350,7 +403,7 @@ volScalarField SfuelNO_2
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("SfuelNO_2",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("SfuelNO_2",dimDensity/dimTime,0.0)
);
volScalarField a
@ -364,7 +417,7 @@ volScalarField a
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("a",dimensionSet(0,0,0,0,0,0,0),0.0)
dimensionedScalar("a",dimless,0.0)
);
volScalarField R1
@ -378,7 +431,7 @@ volScalarField R1
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("R1",dimensionSet(0,0,-1,0,0,0,0),0.0)
dimensionedScalar("R1",dimless/dimTime,0.0)
);
volScalarField R2
@ -392,7 +445,7 @@ volScalarField R2
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("R1",dimensionSet(0,0,-1,0,0,0,0),0.0)
dimensionedScalar("R1",dimless/dimTime,0.0)
);
volScalarField Snh3_p
@ -406,7 +459,7 @@ volScalarField Snh3_p
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Snh3_p",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("Snh3_p",dimDensity/dimTime,0.0)
);
volScalarField Snh3_1
@ -420,7 +473,7 @@ volScalarField Snh3_1
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Snh3_1",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("Snh3_1",dimDensity/dimTime,0.0)
);
volScalarField Snh3_2
@ -434,7 +487,7 @@ volScalarField Snh3_2
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Snh3_2",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("Snh3_2",dimDensity/dimTime,0.0)
);
volScalarField Snh3
@ -448,11 +501,9 @@ volScalarField Snh3
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Snh3",dimensionSet(1,-3,-1,0,0,0,0),0.0)
dimensionedScalar("Snh3",dimDensity/dimTime,0.0)
);
volScalarField kf1
(
IOobject
@ -464,7 +515,7 @@ volScalarField kf1
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf1",dimensionSet(0,3,-1,0,-1,0,0),0.0)
dimensionedScalar("kf1",dimVolume/dimTime/dimMoles,0.0)
);
volScalarField kf2
@ -478,7 +529,7 @@ volScalarField kf2
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf2",dimensionSet(0,3,-1,0,-1,0,0),0.0)
dimensionedScalar("kf2",dimVolume/dimTime/dimMoles,0.0)
);
volScalarField kf3
@ -492,7 +543,7 @@ volScalarField kf3
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kf3",dimensionSet(0,3,-1,0,-1,0,0),0.0)
dimensionedScalar("kf3",dimVolume/dimTime/dimMoles,0.0)
);
volScalarField kr1
@ -506,7 +557,7 @@ volScalarField kr1
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kr1",dimensionSet(0,3,-1,0,-1,0,0),0.0)
dimensionedScalar("kr1",dimVolume/dimTime/dimMoles,0.0)
);
volScalarField kr2
@ -520,72 +571,90 @@ volScalarField kr2
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("kr2",dimensionSet(0,3,-1,0,-1,0,0),0.0)
dimensionedScalar("kr2",dimVolume/dimTime/dimMoles,0.0)
);
const dimensionedScalar Wno ("Wno", dimMass/dimMoles, 0.030006); // Molecular weight of NO [kg/mol]
const dimensionedScalar Wn2 ("Wn2", dimMass/dimMoles, 0.028013); // Molecular weight of N2 [kg/mol]
const dimensionedScalar Wo2 ("Wo2", dimMass/dimMoles, 0.031998); // Molecular weight of O2 [kg/mol]
const dimensionedScalar Wh2o ("Wh2o", dimMass/dimMoles, 0.018015); // Molecular weight of H2O[kg/mol]
const dimensionedScalar Wo ("Wo", dimMass/dimMoles, 0.015999); // Molecular weight of O [kg/mol]
const dimensionedScalar Woh ("Woh", dimMass/dimMoles, 0.017008); // Molecular weight of OH [kg/mol]
const dimensionedScalar Wnh3 ("Wnh3", dimMass/dimMoles, 0.017031); // Molecular weight of NH3[kg/mol]
const dimensionedScalar Wn ("Wn", dimMass/dimMoles, 0.014007); // Molecular weight of N [kg/mol]
volScalarField Sfuel1
(
IOobject
(
"Sfuel1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
const dimensionedScalar Rgas (Foam::constant::physicoChemical::R);
const scalar RR = Rgas.value(); // Universal gas constant [J/(mol*K)]
const scalar A1 = 4.0E+06; // [1/s]
const scalar A2 = 1.8E+08; // [1/s]
const scalar E1 = 133947.2; // [J/mol]
const scalar E2 = 113017.95; // [J/mol]
volScalarField Sfuel2
(
IOobject
(
"Sfuel2",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// molar concentrations
Xn2 = rho*N2/Wn2; // Molar concentration of N2 [mol/m^3]
Xo2 = rho*O2/Wo2; // Molar concentration of O2 [mol/m^3]
Xh2o = rho*H2O/Wh2o; // Molar concentration of H2O[mol/m^3]
IOdictionary fuelNitrogenMassFraction
(
IOobject
(
"fuelNitrogenMassFraction",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// O,OH molar concentrations
if(instantaneousRadicals)
{
Xo = rho*O/Wo; // Molar concentration of O [mol/m^3]
Xoh = rho*OH/Woh; // Molar concentration of OH[mol/m^3]
}
else
{
forAll (mesh.cells(),celli)
{
// Molar concentration of O [mol/m^3], Partial Equilibrium Approach
Xo[celli] = 36.64 * Foam::pow(T[celli],0.5)
* Foam::pow(Xo2[celli],0.5) * Foam::exp(-27123.0/T[celli]);
// Molar concentration of OH [mol/m^3], Partial Equilibrium Approach
Xoh[celli] = 2.129E+02 * Foam::pow(T[celli],-0.57) * Foam::exp(-4595.0/T[celli])
* Foam::pow(Xo[celli],0.5) * Foam::pow(Xh2o[celli],0.5);
}
}
volScalarField Yn
(
IOobject
(
"Yn",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(fuelNitrogenMassFraction.lookup("Yn"))
);
forAll (mesh.cells(),celli)
{
// Thermal NO
// Reaction rate constants [m^3/(mol*s)]
kf1[celli] = 1.8E+08*Foam::exp(-38370/T[celli]);
kf2[celli] = 1.8E+04*T[celli]*Foam::exp(-4680/T[celli]);
kf3[celli] = 7.1E+07*Foam::exp(-450/T[celli]);
kr1[celli] = 3.8E+07*Foam::exp(-425/T[celli]);
kr2[celli] = 3.81E+03*T[celli]*Foam::exp(-20820/T[celli]);
}
IOdictionary modelParameter
(
IOobject
(
"modelParameter",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// a, Order of Reaction "NH3 + O2 -> NO"
forAll (mesh.cells(),celli)
{
if (Xo2[celli] <= 4.1E-03)
{
a[celli] = 1.0;
}
else if ( Xo2[celli] > 4.1E-03 and Xo2[celli] <= 1.11E-02 )
{
a[celli] = -3.95 - 0.9*Foam::log(Xo2[celli]);
}
else if ( Xo2[celli] > 1.11E-02 and Xo2[celli] <= 0.03 )
{
a[celli] = -0.35 - 0.1*Foam::log(Xo2[celli]);
}
else // if ( Xo2[celli] > 0.03 )
{
a[celli] = 0.0;
}
}
// NH3 source
forAll (mesh.cells(),celli)
{
scalar nh3Source = 0.0;
forAll (fuelSources, si)
{
nh3Source += fuelSources[si][celli] * nitrogenMassInFuels[si];
}
Snh3_p[celli] = nh3Source*Wnh3.value()/Wn.value()/mesh.V()[celli]; // NH3 production
}
bool instantaneousRadicals(modelParameter.lookupOrDefault("instantaneousRadicals", true));

BIN
tutorials/combustion/NOxFoam_thermal/constant/.chem_FOAM.swp
View file

Binary file not shown.

BIN
tutorials/combustion/NOxFoam_thermal/constant/.g.swp
View file

Binary file not shown.

BIN
tutorials/combustion/NOxFoam_thermal/constant/.thermo.incompressiblePoly.swp
View file

Binary file not shown.

18
tutorials/combustion/NOxFoam_thermalAndFuel/constant/fuelNitrogenMassFraction → tutorials/combustion/NOxFoam_thermal/constant/NOxProperties
View file

@ -1,23 +1,25 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object fuelNitrogenMassFraction;
version 2;
format ascii;
class dictionary;
location "constant";
object SOxProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
fuelSourceNames ( reactingCloud1:rhoTrans_C7H16 );
Yn Yn [ 0 0 0 0 0 0 0 ] 0.01;
nitrogenMassInFuels (0.01);
instantaneousRadicals false; //true;
calculateFuelNOx true;
// ************************************************************************* //

BIN
tutorials/combustion/NOxFoam_thermalAndFuel/constant/.chem_FOAM.swp
View file

Binary file not shown.

BIN
tutorials/combustion/NOxFoam_thermalAndFuel/constant/.g.swp
View file

Binary file not shown.

BIN
tutorials/combustion/NOxFoam_thermalAndFuel/constant/.thermo.incompressiblePoly.swp
View file

Binary file not shown.

20
tutorials/combustion/NOxFoam_thermal/constant/modelParameter → tutorials/combustion/NOxFoam_thermalAndFuel/constant/NOxProperties
View file

@ -7,13 +7,19 @@
\*---------------------------------------------------------------------------*/
FoamFile
{
version 4.x;
format ascii;
class dictionary;
location "constant";
object fuelNitrogenMassFraction;
version 2;
format ascii;
class dictionary;
location "constant";
object SOxProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
fuelSourceNames ( reactingCloud1:rhoTrans_C7H16 );
instantaneousRadicals false; //true;
nitrogenMassInFuels (0.01);
instantaneousRadicals false; //true;
calculateFuelNOx true;
// ************************************************************************* //

19
tutorials/combustion/NOxFoam_thermalAndFuel/constant/modelParameter
View file

@ -1,19 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 4.x;
format ascii;
class dictionary;
location "constant";
object fuelNitrogenMassFraction;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
instantaneousRadicals false; //true;