update NOxFoam tutorial cases according to solver changes

NOxFoam, SOxFoam tutorial cases were uploaded.

See merge request combustion-postech/OpenFOAM-4.x!15

.gitlab-ci.yml

@@ -1,23 +1,9 @@
+before_script:
+  - "[[ -d ../ThirdParty-4.x ]] || git clone https://github.com/OpenFOAM/ThirdParty-4.x.git ../ThirdParty-4.x"
+  - source etc/bashrc
+  - rm -rf platforms
+  - wclean all
+
 allmake:
-  image: park0d/of4builder
-  tags:
-    - openfoam4
-  before_script:
-    - "[[ -d ../ThirdParty-4.x ]] || git clone https://github.com/OpenFOAM/ThirdParty-4.x.git ../ThirdParty-4.x"
-    - source etc/bashrc
-    - rm -rf platforms
-    - wclean all
   script:
     - ./Allwmake -j
-  artifacts:
-    paths:
-      - ./
-
-cantera:
-    image: park0d/of4builder
-    tags:
-        - openfoam4
-    script:
-        - whereis libcantera
-    
-  

applications/solvers/combustion/NOxFoam_thermal_fuelNOx/Make/files

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

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

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);
+}

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);
+}

applications/solvers/combustion/NOxFoam_thermal_fuelNOx/NOxFoam.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#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)
+    {
+        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;
+    }
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #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;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/NOxFoam_thermal_fuelNOx/SourceCalc.H

@@ -0,0 +1,42 @@
+    Xno = rho*NO/Wno;    // Molar concentration of NO [mol/m^3]
+
+    // Calculation of source term. (Thermal NO)
+    forAll (mesh.cells(),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();
+    }
+
+    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;
+    }

applications/solvers/combustion/NOxFoam_thermal_fuelNOx/createFields.H

@@ -0,0 +1,660 @@
+
+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
+    (
+        "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
+);
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        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
+    ),
+    rho * nut
+);
+
+volScalarField Xno
+(
+    IOobject
+    (
+        "Xno",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xno",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xn2
+(
+    IOobject
+    (
+        "Xn2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xn2",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xo2
+(
+    IOobject
+    (
+        "Xo2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xo2",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xh2o
+(
+    IOobject
+    (
+        "Xh2o",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xh2o",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xnh3
+(
+    IOobject
+    (
+        "Xnh3",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xnh3",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xo
+(
+    IOobject
+    (
+        "Xo",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xo",dimMoles/dimVolume,0.0)
+);
+
+volScalarField Xoh
+(
+    IOobject
+    (
+        "Xoh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Xoh",dimMoles/dimVolume,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::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("NH3",dimless,0.0)
+);
+
+if (calculateFuelNOx)
+{
+    NH3.writeOpt() = IOobject::AUTO_WRITE;
+}
+
+volScalarField Sno
+(
+    IOobject
+    (
+        "Sno",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Sno",dimDensity/dimTime,0.0)
+);
+
+volScalarField SthermNO
+(
+    IOobject
+    (
+        "SthermNO",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("SthermNO",dimDensity/dimTime,0.0)
+);
+
+volScalarField SfuelNO
+(
+    IOobject
+    (
+        "SfuelNO",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("SfuelNO",dimDensity/dimTime,0.0)
+);
+
+volScalarField SfuelNO_1
+(
+    IOobject
+    (
+        "SfuelNO_1",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("SfuelNO_1",dimDensity/dimTime,0.0)
+);
+
+volScalarField SfuelNO_2
+(
+    IOobject
+    (
+        "SfuelNO_2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("SfuelNO_2",dimDensity/dimTime,0.0)
+);
+
+volScalarField a
+(
+    IOobject
+    (
+        "a",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("a",dimless,0.0)
+);
+
+volScalarField R1
+(
+    IOobject
+    (
+        "R1",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("R1",dimless/dimTime,0.0)
+);
+
+volScalarField R2
+(
+    IOobject
+    (
+        "R2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("R1",dimless/dimTime,0.0)
+);
+
+volScalarField Snh3_p
+(
+    IOobject
+    (
+        "Snh3_p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Snh3_p",dimDensity/dimTime,0.0)
+);
+
+volScalarField Snh3_1
+(
+    IOobject
+    (
+        "Snh3_1",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Snh3_1",dimDensity/dimTime,0.0)
+);
+
+volScalarField Snh3_2
+(
+    IOobject
+    (
+        "Snh3_2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Snh3_2",dimDensity/dimTime,0.0)
+);
+
+volScalarField Snh3
+(
+    IOobject
+    (
+        "Snh3",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Snh3",dimDensity/dimTime,0.0)
+);
+
+volScalarField kf1
+(
+    IOobject
+    (
+        "kf1",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("kf1",dimVolume/dimTime/dimMoles,0.0)
+);
+
+volScalarField kf2
+(
+    IOobject
+    (
+        "kf2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("kf2",dimVolume/dimTime/dimMoles,0.0)
+);
+
+volScalarField kf3
+(
+    IOobject
+    (
+        "kf3",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("kf3",dimVolume/dimTime/dimMoles,0.0)
+);
+
+volScalarField kr1
+(
+    IOobject
+    (
+        "kr1",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("kr1",dimVolume/dimTime/dimMoles,0.0)
+);
+
+volScalarField kr2
+(
+    IOobject
+    (
+        "kr2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    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]
+
+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]
+
+// 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]
+
+// 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);
+    }
+}
+
+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]);
+}
+
+// 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
+}
+

applications/solvers/combustion/SOxFoam/Make/files

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

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

applications/solvers/combustion/SOxFoam/SOxFoam.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    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 "calculateSso2.H"
+
+        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;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/SOxFoam/calculateSso2.H

@@ -0,0 +1,14 @@
+const scalar Wso2 = 0.064066;   // Molecular weight of SO2[kg/mol]
+const scalar Ws   = 0.032065;   // Molecular weight of S  [kg/mol]
+
+forAll (mesh.cells(),celli)
+{
+    scalar sulfurSource = 0.0;
+
+    forAll (fuelSources, si)
+    {
+        sulfurSource += fuelSources[si][celli] * sulfurMassInFuels[si];
+    }
+
+    Sso2[celli] = sulfurSource*Wso2/Ws/mesh.V()[celli];
+}

applications/solvers/combustion/SOxFoam/createFields.H

@@ -0,0 +1,130 @@
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    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
+    ),
+    rho * nut
+);
+
+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",dimDensity/dimTime,0.0)
+);
+
+IOdictionary SOxProperties
+(
+    IOobject
+    (
+        "SOxProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+wordList fuelSourceNames(SOxProperties.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 sulfurMassInFuels(SOxProperties.lookup("sulfurMassInFuels"));

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

@@ -148,7 +148,7 @@ irreversibleReactionDelimiter {space}"=>"{space}
 startPDependentSpecie {space}"("{space}"+"{space}
 pDependentSpecie      {specieName}")"{space}
 reactionCoeffs        {space}{floatNum}{some_space}{floatNum}{some_space}{floatNum}{space}
-reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*((\+({some_space}))?){space}
+reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*{space}
 reactionKeywordSlash  {reactionKeyword}"/"{space}
 thirdBodyEfficiency   {space}{floatNum}{space}"/"{space}
 startReactionCoeffs   {space}"/"{space}
@@ -637,8 +637,7 @@ bool finishReaction = false;
                     (
                         currentSpecieName,
                         1.0,
-                        molecularWeight(currentSpecieComposition),
+                        molecularWeight(currentSpecieComposition)
-                        chargeNumber(currentSpecieComposition)
                     ),
                     currentLowT,
                     currentHighT,

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

@@ -144,28 +144,6 @@ Foam::scalar Foam::chemkinReader::molecularWeight
 }
 
 
-Foam::scalar Foam::chemkinReader::chargeNumber
-(
-    const List<specieElement>& specieComposition
-) const
-{
-    scalar nElemCharges = 0.0;
-
-    forAll(specieComposition, i)
-    {
-        label nAtoms = specieComposition[i].nAtoms;
-        const word& elementName = specieComposition[i].elementName;
-
-        if (elementName[0] == 'e')
-        {
-            nElemCharges -= nAtoms;
-        }
-    }
-
-    return nElemCharges;
-}
-
-
 void Foam::chemkinReader::checkCoeffs
 (
     const scalarList& reactionCoeffs,

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H

@@ -254,11 +254,6 @@ private:
             const List<specieElement>& specieComposition
         ) const;
 
-        scalar chargeNumber
-        (
-            const List<specieElement>& specieComposition
-        ) const;
-
         void finishElements(labelList& currentAtoms);
 
         void checkCoeffs

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C

@@ -68,70 +68,6 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
 }
 
 
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::z
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).z();
-}
-
-
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Qc
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).Qc();
-}
-
-
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Qc2
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).Qc2();
-}
-
-
-template<class MixtureType>
-void Foam::SpecieMixture<MixtureType>::Cp
-(
-    scalarField& Cps,
-    const scalar p,
-    const scalar T
-) const
-{
-    forAll(this->species(), speciei)
-    {
-        Cps[speciei] = this->getLocalThermo(speciei).Cp(p, T);
-    }
-
-    return;
-}
-
-
-template<class MixtureType>
-void Foam::SpecieMixture<MixtureType>::Cv
-(
-    scalarField& Cvs,
-    const scalar p,
-    const scalar T
-) const
-{
-    forAll(this->species(), speciei)
-    {
-        Cvs[speciei] = this->getLocalThermo(speciei).Cv(p, T);
-    }
-
-    return;
-}
-
-
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
 (

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H

@@ -81,34 +81,6 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const;
 
-            //- Number of charges of the given specie []
-            virtual scalar z(const label speciei) const;
-
-            //- Specific charge of the given specie []
-            virtual scalar Qc(const label specieI) const;
-
-            //- Absolute specific charge of the given specie []
-            virtual scalar Qc2(const label specieI) const;
-
-
-        // All species thermo properties
-
-            //- Heat capacity at constant pressure [J/(kg K)]
-            virtual void Cp
-            (
-                scalarField& Cps,
-                const scalar p,
-                const scalar T
-            ) const;
-
-            //- Heat capacity at constant volume [J/(kg K)]
-            virtual void Cv
-            (
-                scalarField& Cvs,
-                const scalar p,
-                const scalar T
-            ) const;
-
 
         // Per specie thermo properties
 

src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C

@@ -81,72 +81,4 @@ Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::W() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc() const
-{
-    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
-
-    tmp<volScalarField> tQc
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Qc", Y[0].group()),
-                Y[0].time().timeName(),
-                Y[0].mesh()
-            ),
-            Y[0].mesh(),
-            dimensionedScalar("zero", dimCurrent*dimTime/dimMass, 0)
-        )
-    );
-
-    volScalarField& Qc = tQc.ref();
-
-    dimensionedScalar Qci("Qci", dimCurrent*dimTime/dimMass, 0);
-
-    forAll(Y, i)
-    {
-        Qci.value() = this->Qc(i);
-        Qc += Y[i]*Qci;
-        // Qc += Y[i]*this->Qc(i);
-    }
-
-    return tQc;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc2() const
-{
-    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
-
-    tmp<volScalarField> tQc
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Qc2", Y[0].group()),
-                Y[0].time().timeName(),
-                Y[0].mesh()
-            ),
-            Y[0].mesh(),
-            dimensionedScalar("zero", sqr(dimCurrent*dimTime)/dimMass, 0)
-        )
-    );
-
-    volScalarField& Qc = tQc.ref();
-
-    dimensionedScalar Qc2i("Qc2i", sqr(dimCurrent*dimTime)/dimMass, 0);
-
-    forAll(Y, i)
-    {
-        Qc2i.value() = this->Qc2(i);
-        Qc += Y[i]*Qc2i;
-        //Qc += Y[i]*this->Qc2(i);
-    }
-
-    return tQc;
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H

@@ -88,43 +88,9 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const = 0;
 
-            //- Number of charges of the given specie []
-            virtual scalar z(const label speciei) const = 0;
-
-            //- Specific charge of the given specie [C/kg]
-            virtual scalar Qc(const label speciei) const = 0;
-
-            //- Absolute specific charge of the given specie [C/kg]
-            virtual scalar Qc2(const label speciei) const = 0;
-
             //- Molecular weight of the mixture [kg/kmol]
             tmp<volScalarField> W() const;
 
-            //- Specific charge of the mixture [C/kg]
-            tmp<volScalarField> Qc() const;
-
-            //- Absolute specific charge of the mixture [C^2/kg]
-            tmp<volScalarField> Qc2() const;
-
-
-        // All species thermo properties
-
-            //- Heat capacity at constant pressure [J/(kg K)]
-            virtual void Cp
-            (
-                scalarField& Cps,
-                const scalar p,
-                const scalar T
-            ) const = 0;
-
-            //- Heat capacity at constant volume [J/(kg K)]
-            virtual void Cv
-            (
-                scalarField& Cvs,
-                const scalar p,
-                const scalar T
-            ) const = 0;
-
 
         // Per specie thermo properties
 

src/thermophysicalModels/specie/specie/specie.C

@@ -29,24 +29,6 @@ License
 
 /* * * * * * * * * * * * * * * public constants  * * * * * * * * * * * * * * */
 
-//- Universal gas constant (default in [J/(kmol K)])
-const Foam::scalar Foam::specie::RR = constant::physicoChemical::R.value()*1000;
-
-//- Standard pressure (default in [Pa])
-const Foam::scalar Foam::specie::Pstd = constant::standard::Pstd.value();
-
-//- Standard temperature (default in [K])
-const Foam::scalar Foam::specie::Tstd = constant::standard::Tstd.value();
-
-//- Elementary charge (default in [C])
-const Foam::scalar Foam::specie::e = constant::electromagnetic::e.value();
-
-//- Avogadro number (default in [1/mol])
-const Foam::scalar Foam::specie::NA = constant::physicoChemical::NA.value()*1000;
-
-//- Boltzmann constant (default in [J/K])
-const Foam::scalar Foam::specie::k = constant::physicoChemical::k.value();
-
 namespace Foam
 {
     defineTypeNameAndDebug(specie, 0);
@@ -59,8 +41,7 @@ Foam::specie::specie(Istream& is)
 :
     name_(is),
     nMoles_(readScalar(is)),
-    molWeight_(readScalar(is)),
+    molWeight_(readScalar(is))
-    nCharges_(readScalar(is))
 {
     is.check("specie::specie(Istream& is)");
 }
@@ -70,8 +51,7 @@ Foam::specie::specie(const dictionary& dict)
 :
     name_(dict.dictName()),
     nMoles_(readScalar(dict.subDict("specie").lookup("nMoles"))),
-    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight"))),
+    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight")))
-    nCharges_(dict.subDict("specie").lookupOrDefault("nCharges", 0.0))
 {}
 
 
@@ -82,7 +62,6 @@ void Foam::specie::write(Ostream& os) const
     dictionary dict("specie");
     dict.add("nMoles", nMoles_);
     dict.add("molWeight", molWeight_);
-    dict.add("nCharges", nCharges_);
     os  << indent << dict.dictName() << dict;
 }
 
@@ -93,8 +72,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const specie& st)
 {
     os  << st.name_ << tab
         << st.nMoles_ << tab
-        << st.molWeight_ << tab
+        << st.molWeight_;
-        << st.nCharges_;
 
     os.check("Ostream& operator<<(Ostream& os, const specie& st)");
     return os;

src/thermophysicalModels/specie/specie/specie.H

@@ -77,9 +77,6 @@ class specie
         //- Molecular weight of specie [kg/kmol]
         scalar molWeight_;
 
-        //- Number of elementary charges of specie
-        scalar nCharges_;
-
 
 public:
 
@@ -87,29 +84,6 @@ public:
     ClassName("specie");
 
 
-    // Public constants
-
-        // Thermodynamic constants
-
-            //- Universal gas constant [J/(kmol K)]
-            static const scalar RR;
-
-            //- Standard pressure [Pa]
-            static const scalar Pstd;
-
-            //- Standard temperature [K]
-            static const scalar Tstd;
-
-            //- Elementary charge [C]
-            static const scalar e;
-
-            //- Avogadro number [1/kmol]
-            static const scalar NA;
-
-            //- Boltzmann constant [J/K]
-            static const scalar k;
-
-
     // Constructors
 
 
@@ -124,23 +98,6 @@ public:
             const scalar molWeight
         );
 
-        //- Construct from components without name
-        inline specie
-        (
-            const scalar nMoles,
-            const scalar molWeight,
-            const scalar nCharges
-        );
-
-        //- Construct from components with name
-        inline specie
-        (
-            const word& name,
-            const scalar nMoles,
-            const scalar molWeight,
-            const scalar nCharges
-        );
-
         //- Construct as copy
         inline specie(const specie&);
 
@@ -170,27 +127,6 @@ public:
             //- Gas constant [J/(kg K)]
             inline scalar R() const;
 
-            //- Charge number
-            inline scalar z() const;
-
-
-        // Electric charge function
-
-
-            // Mole specific properties
-
-                //- Electric charge [C/kmol]
-                inline scalar qc() const;
-
-
-            // Mass specific properties
-
-                //- Electric charge [C/kg]
-                inline scalar Qc() const;
-
-                //- Electric charge squared [C^2/kg]
-                inline scalar Qc2() const;
-
 
         // I-O
 

src/thermophysicalModels/specie/specie/specieI.H

@@ -41,8 +41,7 @@ inline specie::specie
 :
     name_(name),
     nMoles_(nMoles),
-    molWeight_(molWeight),
+    molWeight_(molWeight)
-    nCharges_(0.0)
 {}
 
 
@@ -53,36 +52,7 @@ inline specie::specie
 )
 :
     nMoles_(nMoles),
-    molWeight_(molWeight),
+    molWeight_(molWeight)
-    nCharges_(0.0)
-{}
-
-
-inline specie::specie
-(
-    const word& name,
-    const scalar nMoles,
-    const scalar molWeight,
-    const scalar nCharges
-)
-:
-    name_(name),
-    nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(nCharges)
-{}
-
-
-inline specie::specie
-(
-    const scalar nMoles,
-    const scalar molWeight,
-    const scalar nCharges
-)
-:
-    nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(nCharges)
 {}
 
 
@@ -92,8 +62,7 @@ inline specie::specie(const specie& st)
 :
     name_(st.name_),
     nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_),
+    molWeight_(st.molWeight_)
-    nCharges_(st.nCharges_)
 {}
 
 
@@ -101,8 +70,7 @@ inline specie::specie(const word& name, const specie& st)
 :
     name_(name),
     nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_),
+    molWeight_(st.molWeight_)
-    nCharges_(st.nCharges_)
 {}
 
 
@@ -132,30 +100,6 @@ inline scalar specie::R() const
 }
 
 
-inline scalar specie::z() const
-{
-    return nCharges_;
-}
-
-
-inline scalar specie::qc() const
-{
-    return nCharges_ * NA * e;
-}
-
-
-inline scalar specie::Qc() const
-{
-    return nCharges_ * NA * e / molWeight_;
-}
-
-
-inline scalar specie::Qc2() const
-{
-    return nCharges_ * nCharges_ * NA * e * e / molWeight_;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline void specie::operator=(const specie& st)
@@ -163,7 +107,6 @@ inline void specie::operator=(const specie& st)
     //name_ = st.name_;
     nMoles_ = st.nMoles_;
     molWeight_ = st.molWeight_;
-    nCharges_ = st.nCharges_;
 }
 
 
@@ -171,10 +114,6 @@ inline void specie::operator+=(const specie& st)
 {
     scalar sumNmoles = max(nMoles_ + st.nMoles_, SMALL);
 
-    nCharges_ =
-        nMoles_/sumNmoles*nCharges_
-      + st.nMoles_/sumNmoles*st.nCharges_;
-
     molWeight_ =
         nMoles_/sumNmoles*molWeight_
       + st.nMoles_/sumNmoles*st.molWeight_;
@@ -191,10 +130,6 @@ inline void specie::operator-=(const specie& st)
         diffnMoles = SMALL;
     }
 
-    nCharges_ =
-        nMoles_/diffnMoles*nCharges_
-      - st.nMoles_/diffnMoles*st.nCharges_;
-
     molWeight_ =
         nMoles_/diffnMoles*molWeight_
       - st.nMoles_/diffnMoles*st.molWeight_;
@@ -219,9 +154,7 @@ inline specie operator+(const specie& st1, const specie& st2)
     (
         sumNmoles,
         st1.nMoles_/sumNmoles*st1.molWeight_
-      + st2.nMoles_/sumNmoles*st2.molWeight_,
+      + st2.nMoles_/sumNmoles*st2.molWeight_
-        st1.nMoles_/sumNmoles*st1.nCharges_
-      + st2.nMoles_/sumNmoles*st2.nCharges_
     );
 }
 
@@ -238,9 +171,7 @@ inline specie operator-(const specie& st1, const specie& st2)
     (
         diffNmoles,
         st1.nMoles_/diffNmoles*st1.molWeight_
-      - st2.nMoles_/diffNmoles*st2.molWeight_,
+      - st2.nMoles_/diffNmoles*st2.molWeight_
-        st1.nMoles_/diffNmoles*st1.nCharges_
-      - st2.nMoles_/diffNmoles*st2.nCharges_
     );
 }
 
@@ -250,8 +181,7 @@ inline specie operator*(const scalar s, const specie& st)
     return specie
     (
         s*st.nMoles_,
-        st.molWeight_,
+        st.molWeight_
-        st.nCharges_
     );
 }
 

src/thermophysicalModels/specie/thermo/thermo/thermo.H

@@ -336,7 +336,6 @@ public:
             ) const;
 
 
-
         // I-O
 
             //- Write to Ostream

tutorials/combustion/NOxFoam_thermal/constant/NOxProperties

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

tutorials/combustion/NOxFoam_thermal/constant/polyMesh/boundary

@@ -0,0 +1,70 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "1/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    INLET
+    {
+        type            patch;	
+        nFaces          279;
+        startFace       88864;
+    }
+    WALL
+    {
+        type            wall;
+        nFaces          360;
+        startFace       89143;
+    }
+    SLIPWALL
+    {
+        type            wall;
+        nFaces          1062;
+        startFace       89503;
+    }
+    OUTLET
+    {
+        type            patch;	
+        nFaces          477;
+        startFace       90565;
+    }
+    CYCLIC_half0
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       91042;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half1;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+    CYCLIC_half1
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       93051;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half0;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+)
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermal/constant/polyMesh/pointZones

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       regIOobject;
+    location    "1/polyMesh";
+    object      pointZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0
+()
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermal/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     NOxFoam;       
+
+startFoam       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2000; 	
+
+deltaT          1;
+
+writeControl    timeStep;	
+
+writeInterval   100;	
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermal/system/decomposeParDict

@@ -0,0 +1,78 @@
+// The FOAM Project // File: decomposeParDict
+/*
+-------------------------------------------------------------------------------
+ =========         | dictionary
+ \\      /         | 
+  \\    /          | Name:   decomposeParDict
+   \\  /           | Family: FoamX configuration file
+    \\/            | 
+    F ield         | FOAM version: 2.1
+    O peration     | Product of Nabla Ltd.
+    A and          | 
+    M anipulation  | Email: Enquiries@Nabla.co.uk
+-------------------------------------------------------------------------------
+*/
+// FoamX Case Dictionary.
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+  //  root            "/home/warhol/chris/foam/chris2.1/run/tutorials/interFoam";
+  //  case            "damBreak";
+   // instance        "system";
+  //  local           "";
+
+    class           dictionary;
+    object          decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 24;  //12;
+
+method          scotch; //metis;
+
+simpleCoeffs
+{
+    n               (1 1 2);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (1 1 1);
+    delta           0.001;
+    order           xyz;
+}
+
+metisCoeffs
+{
+    processorWeights 
+    (
+        1
+        1
+1
+1
+1
+1
+1
+
+    );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           
+(
+);
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermal/system/fvSchemes

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,NO)      bounded Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermal/system/fvSolution

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+
+    NO
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       0;
+        relTol          0.1;
+        maxIter         20;
+    }
+}
+
+SIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    rhoMin          0.1;
+    rhoMax          1.5;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"            0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/0/NH3

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "2100";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;           
+
+boundaryField
+{
+    INLET
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+    WALL
+    {
+        type            zeroGradient;
+    }
+    SLIPWALL
+    {
+        type            zeroGradient;
+    }
+    OUTLET
+    {
+        type            zeroGradient;
+    }
+    CYCLIC_half0
+    {
+        type            cyclic;
+    }
+    CYCLIC_half1
+    {
+        type            cyclic;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/constant/NOxProperties

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

tutorials/combustion/NOxFoam_thermalAndFuel/constant/polyMesh/boundary

@@ -0,0 +1,70 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "1/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    INLET
+    {
+        type            patch;	
+        nFaces          279;
+        startFace       88864;
+    }
+    WALL
+    {
+        type            wall;
+        nFaces          360;
+        startFace       89143;
+    }
+    SLIPWALL
+    {
+        type            wall;
+        nFaces          1062;
+        startFace       89503;
+    }
+    OUTLET
+    {
+        type            patch;	
+        nFaces          477;
+        startFace       90565;
+    }
+    CYCLIC_half0
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       91042;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half1;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+    CYCLIC_half1
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       93051;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half0;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+)
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/constant/polyMesh/pointZones

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       regIOobject;
+    location    "1/polyMesh";
+    object      pointZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0
+()
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     NOxFoam;      
+
+startFoam       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         50000; 	
+
+deltaT          1;
+
+writeControl    timeStep;	
+
+writeInterval   1000;	
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/system/decomposeParDict

@@ -0,0 +1,57 @@
+// The FOAM Project // File: decomposeParDict
+/*
+-------------------------------------------------------------------------------
+ =========         | dictionary
+ \\      /         | 
+  \\    /          | Name:   decomposeParDict
+   \\  /           | Family: FoamX configuration file
+    \\/            | 
+    F ield         | FOAM version: 2.1
+    O peration     | Product of Nabla Ltd.
+    A and          | 
+    M anipulation  | Email: Enquiries@Nabla.co.uk
+-------------------------------------------------------------------------------
+*/
+// FoamX Case Dictionary.
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    object          decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 24;  
+
+method          scotch; 
+
+simpleCoeffs
+{
+    n               (1 1 2);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (1 1 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           
+(
+);
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/system/fvSchemes

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,NO)      bounded Gauss linearUpwind grad(NO);
+    div(phi,NH3)      bounded Gauss linearUpwind grad(NH3);
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/combustion/NOxFoam_thermalAndFuel/system/fvSolution

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "(NO|NH3)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       0;
+        relTol          0.1;
+        maxIter         20;
+    }
+}
+
+SIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    rhoMin          0.1;
+    rhoMax          1.5;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        NH3             0.001;
+        NO              0.01;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/constant/SOxProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "constant";
+    object          SOxProperties;
+}
+
+fuelSourceNames         ( reactingCloud1:rhoTrans_C7H16 );
+
+sulfurMassInFuels       ( 0.01 );
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/constant/polyMesh/boundary

@@ -0,0 +1,70 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "1/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    INLET
+    {
+        type            patch;	
+        nFaces          279;
+        startFace       88864;
+    }
+    WALL
+    {
+        type            wall;
+        nFaces          360;
+        startFace       89143;
+    }
+    SLIPWALL
+    {
+        type            wall;
+        nFaces          1062;
+        startFace       89503;
+    }
+    OUTLET
+    {
+        type            patch;	
+        nFaces          477;
+        startFace       90565;
+    }
+    CYCLIC_half0
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       91042;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half1;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+    CYCLIC_half1
+    {
+        type            cyclic;
+        inGroups        1(cyclic);
+        nFaces          2009;
+        startFace       93051;
+        matchTolerance  0.0001;
+        transform       rotational;
+        neighbourPatch  CYCLIC_half0;
+        rotationAxis    (1 0 0);
+        rotationCentre  (0 0 0);
+    }
+)
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/constant/polyMesh/pointZones

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       regIOobject;
+    location    "1/polyMesh";
+    object      pointZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+0
+()
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     SOxFoam;      
+
+startFoam       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         3000; 	
+
+deltaT          1;
+
+writeControl    timeStep;	
+
+writeInterval   300;	
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/system/decomposeParDict

@@ -0,0 +1,78 @@
+// The FOAM Project // File: decomposeParDict
+/*
+-------------------------------------------------------------------------------
+ =========         | dictionary
+ \\      /         | 
+  \\    /          | Name:   decomposeParDict
+   \\  /           | Family: FoamX configuration file
+    \\/            | 
+    F ield         | FOAM version: 2.1
+    O peration     | Product of Nabla Ltd.
+    A and          | 
+    M anipulation  | Email: Enquiries@Nabla.co.uk
+-------------------------------------------------------------------------------
+*/
+// FoamX Case Dictionary.
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+  //  root            "/home/warhol/chris/foam/chris2.1/run/tutorials/interFoam";
+  //  case            "damBreak";
+   // instance        "system";
+  //  local           "";
+
+    class           dictionary;
+    object          decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 24;  //12;
+
+method          scotch; //metis;
+
+simpleCoeffs
+{
+    n               (1 1 2);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (1 1 1);
+    delta           0.001;
+    order           xyz;
+}
+
+metisCoeffs
+{
+    processorWeights 
+    (
+        1
+        1
+1
+1
+1
+1
+1
+
+    );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           
+(
+);
+
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/system/fvSchemes

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,SO2)      bounded Gauss linearUpwind grad(SO2);
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/combustion/SOxFoam/system/fvSolution

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.x                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    SO2
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       0;
+        relTol          0.1;
+        maxIter         20;
+    }
+
+}
+
+SIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    rhoMin          0.1;
+    rhoMax          1.5;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"            0.7;
+    }
+}
+
+
+// ************************************************************************* //