diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files b/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files deleted file mode 100644 index 412c5852c..000000000 --- a/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files +++ /dev/null @@ -1,4 +0,0 @@ -NOxFoam.C - -EXE = $(FOAM_USER_APPBIN)/NOxFoam_thermalNOx - diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/Make/options b/applications/solvers/combustion/NOxFoam_thermalNOx/Make/options deleted file mode 100644 index 21381519d..000000000 --- a/applications/solvers/combustion/NOxFoam_thermalNOx/Make/options +++ /dev/null @@ -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 diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/NOxFoam.C b/applications/solvers/combustion/NOxFoam_thermalNOx/NOxFoam.C deleted file mode 100644 index a66da6afd..000000000 --- a/applications/solvers/combustion/NOxFoam_thermalNOx/NOxFoam.C +++ /dev/null @@ -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 . - -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; -} - - -// ************************************************************************* // diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/SnoCalc.H b/applications/solvers/combustion/NOxFoam_thermalNOx/SnoCalc.H deleted file mode 100644 index ca0d43097..000000000 --- a/applications/solvers/combustion/NOxFoam_thermalNOx/SnoCalc.H +++ /dev/null @@ -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(); - - } - - diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/createFields.H b/applications/solvers/combustion/NOxFoam_thermalNOx/createFields.H deleted file mode 100644 index f37042cc4..000000000 --- a/applications/solvers/combustion/NOxFoam_thermalNOx/createFields.H +++ /dev/null @@ -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)); diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files index 7772265f4..af4e7a7c3 100644 --- a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files +++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files @@ -1,5 +1,5 @@ NOxFoam.C -EXE = $(FOAM_APPBIN)/NOxFoam_thermFuel_2 +EXE = $(FOAM_APPBIN)/NOxFoam diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C index 4fe2bb280..5000ece86 100644 --- a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C +++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C @@ -22,11 +22,9 @@ License along with OpenFOAM. If not, see . 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" diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H index abb6f3a27..8092672b6 100644 --- a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H +++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H @@ -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; + } diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H index 28097cd1d..116a4683f 100644 --- a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H +++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H @@ -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 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)); diff --git a/tutorials/combustion/NOxFoam_thermal/constant/.chem_FOAM.swp b/tutorials/combustion/NOxFoam_thermal/constant/.chem_FOAM.swp deleted file mode 100644 index 59d02ce57..000000000 Binary files a/tutorials/combustion/NOxFoam_thermal/constant/.chem_FOAM.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermal/constant/.g.swp b/tutorials/combustion/NOxFoam_thermal/constant/.g.swp deleted file mode 100644 index cc596dead..000000000 Binary files a/tutorials/combustion/NOxFoam_thermal/constant/.g.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermal/constant/.thermo.incompressiblePoly.swp b/tutorials/combustion/NOxFoam_thermal/constant/.thermo.incompressiblePoly.swp deleted file mode 100644 index b6f71dae1..000000000 Binary files a/tutorials/combustion/NOxFoam_thermal/constant/.thermo.incompressiblePoly.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/fuelNitrogenMassFraction b/tutorials/combustion/NOxFoam_thermal/constant/NOxProperties similarity index 60% rename from tutorials/combustion/NOxFoam_thermalAndFuel/constant/fuelNitrogenMassFraction rename to tutorials/combustion/NOxFoam_thermal/constant/NOxProperties index 68813842e..493108704 100644 --- a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/fuelNitrogenMassFraction +++ b/tutorials/combustion/NOxFoam_thermal/constant/NOxProperties @@ -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; // ************************************************************************* // diff --git a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.chem_FOAM.swp b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.chem_FOAM.swp deleted file mode 100644 index 59d02ce57..000000000 Binary files a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.chem_FOAM.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.g.swp b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.g.swp deleted file mode 100644 index cc596dead..000000000 Binary files a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.g.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.thermo.incompressiblePoly.swp b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.thermo.incompressiblePoly.swp deleted file mode 100644 index b6f71dae1..000000000 Binary files a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/.thermo.incompressiblePoly.swp and /dev/null differ diff --git a/tutorials/combustion/NOxFoam_thermal/constant/modelParameter b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/NOxProperties similarity index 58% rename from tutorials/combustion/NOxFoam_thermal/constant/modelParameter rename to tutorials/combustion/NOxFoam_thermalAndFuel/constant/NOxProperties index bb08e1ba4..493108704 100644 --- a/tutorials/combustion/NOxFoam_thermal/constant/modelParameter +++ b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/NOxProperties @@ -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; + +// ************************************************************************* // diff --git a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/modelParameter b/tutorials/combustion/NOxFoam_thermalAndFuel/constant/modelParameter deleted file mode 100644 index bb08e1ba4..000000000 --- a/tutorials/combustion/NOxFoam_thermalAndFuel/constant/modelParameter +++ /dev/null @@ -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;