diff --git a/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files b/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files
new file mode 100644
index 000000000..412c5852c
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermalNOx/Make/files
@@ -0,0 +1,4 @@
+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
new file mode 100644
index 000000000..21381519d
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermalNOx/Make/options
@@ -0,0 +1,10 @@
+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
new file mode 100644
index 000000000..a66da6afd
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermalNOx/NOxFoam.C
@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
new file mode 100644
index 000000000..ca0d43097
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermalNOx/SnoCalc.H
@@ -0,0 +1,44 @@
+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
new file mode 100644
index 000000000..f37042cc4
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermalNOx/createFields.H
@@ -0,0 +1,340 @@
+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
new file mode 100644
index 000000000..7772265f4
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files
@@ -0,0 +1,5 @@
+NOxFoam.C
+
+EXE = $(FOAM_APPBIN)/NOxFoam_thermFuel_2
+
+
diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/options b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/options
new file mode 100644
index 000000000..81352f91b
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I${LIB_SRC}/meshTools/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lfvOptions \
+ -lsampling
diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NH3Eqn.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NH3Eqn.H
new file mode 100644
index 000000000..355accc2d
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NH3Eqn.H
@@ -0,0 +1,17 @@
+// NH3 transport equation
+while (simple.correctNonOrthogonal())
+{
+ fvScalarMatrix NH3Eqn
+ (
+ fvm::ddt(rho, NH3)
+ + fvm::div(phi, NH3)
+ - fvm::laplacian(mut, NH3)
+ ==
+ Snh3
+ );
+
+ NH3Eqn.relax();
+ fvOptions.constrain(NH3Eqn);
+ NH3Eqn.solve();
+ fvOptions.correct(NH3);
+}
diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOEqn.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOEqn.H
new file mode 100644
index 000000000..94fc44e6f
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOEqn.H
@@ -0,0 +1,17 @@
+// NO transport equation
+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);
+}
diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C
new file mode 100644
index 000000000..4fe2bb280
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C
@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ simpleReactingParcelFoam
+
+Description
+ Steady state solver for compressible, turbulent flow with reacting,
+ multiphase particle clouds and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#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 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"
+
+ 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_thermal_fuelNOx/SourceCalc.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H
new file mode 100644
index 000000000..abb6f3a27
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H
@@ -0,0 +1,94 @@
+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]
+
+ 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]*
+ ((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)]
+
+ }
+
+
diff --git a/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H
new file mode 100644
index 000000000..28097cd1d
--- /dev/null
+++ b/applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H
@@ -0,0 +1,591 @@
+
+volScalarField T
+(
+ IOobject
+ (
+ "T",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField p
+(
+ IOobject
+ (
+ "p",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+
+volVectorField U
+(
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField rho
+(
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+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 Xnh3
+(
+ IOobject
+ (
+ "Xnh3",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Xnh3",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 NH3
+(
+ IOobject
+ (
+ "NH3",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField Sno
+(
+ IOobject
+ (
+ "Sno",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Sno",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField SthermNO
+(
+ IOobject
+ (
+ "SthermNO",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("SthermNO",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField SfuelNO
+(
+ IOobject
+ (
+ "SfuelNO",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("SfuelNO",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField SfuelNO_1
+(
+ IOobject
+ (
+ "SfuelNO_1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("SfuelNO_1",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField SfuelNO_2
+(
+ IOobject
+ (
+ "SfuelNO_2",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("SfuelNO_2",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField a
+(
+ IOobject
+ (
+ "a",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("a",dimensionSet(0,0,0,0,0,0,0),0.0)
+);
+
+volScalarField R1
+(
+ IOobject
+ (
+ "R1",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("R1",dimensionSet(0,0,-1,0,0,0,0),0.0)
+);
+
+volScalarField R2
+(
+ IOobject
+ (
+ "R2",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("R1",dimensionSet(0,0,-1,0,0,0,0),0.0)
+);
+
+volScalarField Snh3_p
+(
+ IOobject
+ (
+ "Snh3_p",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Snh3_p",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField Snh3_1
+(
+ IOobject
+ (
+ "Snh3_1",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Snh3_1",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField Snh3_2
+(
+ IOobject
+ (
+ "Snh3_2",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Snh3_2",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField Snh3
+(
+ IOobject
+ (
+ "Snh3",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Snh3",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)
+);
+
+
+volScalarField Sfuel1
+(
+ IOobject
+ (
+ "Sfuel1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField Sfuel2
+(
+ IOobject
+ (
+ "Sfuel2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+IOdictionary fuelNitrogenMassFraction
+(
+ IOobject
+ (
+ "fuelNitrogenMassFraction",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+);
+
+volScalarField Yn
+(
+ IOobject
+ (
+ "Yn",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar(fuelNitrogenMassFraction.lookup("Yn"))
+);
+
+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/SOxFoam/Make/files b/applications/solvers/combustion/SOxFoam/Make/files
new file mode 100644
index 000000000..81fc43a64
--- /dev/null
+++ b/applications/solvers/combustion/SOxFoam/Make/files
@@ -0,0 +1,4 @@
+SOxFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/SOxFoam
+
diff --git a/applications/solvers/combustion/SOxFoam/Make/options b/applications/solvers/combustion/SOxFoam/Make/options
new file mode 100644
index 000000000..75d5f6e93
--- /dev/null
+++ b/applications/solvers/combustion/SOxFoam/Make/options
@@ -0,0 +1,11 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I${LIB_SRC}/meshTools/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lfvOptions \
+ -lsampling
diff --git a/applications/solvers/combustion/SOxFoam/SOxFoam.C b/applications/solvers/combustion/SOxFoam/SOxFoam.C
new file mode 100644
index 000000000..1069f167f
--- /dev/null
+++ b/applications/solvers/combustion/SOxFoam/SOxFoam.C
@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ simpleReactingParcelFoam
+
+Description
+ Steady state solver for compressible, turbulent flow with reacting,
+ multiphase particle clouds and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#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;
+
+ while (simple.loop())
+ {
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ #include "Sso2Calc.H"
+ mut=rho*nut;
+
+ while (simple.correctNonOrthogonal())
+ {
+ fvScalarMatrix SO2Eqn
+ (
+ fvm::ddt(rho, SO2)
+ + fvm::div(phi, SO2)
+ - fvm::laplacian(mut, SO2)
+ ==
+ Sso2
+ );
+
+ SO2Eqn.relax();
+ fvOptions.constrain(SO2Eqn);
+ SO2Eqn.solve();
+ fvOptions.correct(SO2);
+ }
+
+ 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/SOxFoam/Sso2Calc.H b/applications/solvers/combustion/SOxFoam/Sso2Calc.H
new file mode 100644
index 000000000..286a89652
--- /dev/null
+++ b/applications/solvers/combustion/SOxFoam/Sso2Calc.H
@@ -0,0 +1,11 @@
+const scalar Wso2 = 0.064066; //Molecular weight of SO2[kg/mol]
+const scalar Ws = 0.032065; // S [kg/mol]
+
+ forAll (mesh.cells(),celli)
+ {
+
+ Sso2[celli] = (Sfuel1[celli]+Sfuel2[celli])*Ys[celli]*Wso2/Ws/mesh.V()[celli];
+
+ }
+
+
diff --git a/applications/solvers/combustion/SOxFoam/createFields.H b/applications/solvers/combustion/SOxFoam/createFields.H
new file mode 100644
index 000000000..51680499c
--- /dev/null
+++ b/applications/solvers/combustion/SOxFoam/createFields.H
@@ -0,0 +1,176 @@
+
+volScalarField rho
+(
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ //thermo.rho()
+ mesh
+);
+
+volVectorField U
+(
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+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)
+);
+
+
+
+//#include "compressibleCreatePhi.H"
+//
+//mesh.setFluxRequired(p.name());
+//
+//dimensionedScalar rhoMax
+//(
+// dimensionedScalar::lookupOrDefault
+// (
+// "rhoMax",
+// simple.dict(),
+// dimDensity,
+// GREAT
+// )
+//);
+//
+//dimensionedScalar rhoMin
+//(
+// dimensionedScalar::lookupOrDefault
+// (
+// "rhoMin",
+// simple.dict(),
+// dimDensity,
+// 0
+// )
+//);
+
+
+volScalarField SO2
+(
+ IOobject
+ (
+ "SO2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField Sso2
+(
+ IOobject
+ (
+ "Sso2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("Sso2",dimensionSet(1,-3,-1,0,0,0,0),0.0)
+);
+
+volScalarField Sfuel1
+(
+ IOobject
+ (
+ "Sfuel1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+volScalarField Sfuel2
+(
+ IOobject
+ (
+ "Sfuel2",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+IOdictionary fuelSulfurMassFraction
+(
+ IOobject
+ (
+ "fuelSulfurMassFraction",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+);
+
+volScalarField Ys
+(
+ IOobject
+ (
+ "Ys",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar(fuelSulfurMassFraction.lookup("Ys"))
+);
+