diff --git a/.gitignore b/.gitignore
index a897b8121..8297918ba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,8 @@
*.bak
*.bak[0-9][0-9]
\#*\#
+.*.swp
+.*.swo
# File-browser settings - anywhere
.directory
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H b/applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H
new file mode 100644
index 000000000..3ccd83fec
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H
@@ -0,0 +1,40 @@
+{
+ volScalarField& he = thermo.he();
+
+ fvScalarMatrix EEqn
+ (
+ fvm::ddt(alphac,rho, he) + mvConvection->fvmDiv(alphacPhi, he)
+ + fvc::ddt(alphac,rho, K) + fvc::div(alphacPhi, K)
+ + (
+ he.name() == "e"
+ ? fvc::div
+ (
+ fvc::absolute(phi/fvc::interpolate(rho), U),
+ p,
+ "div(phiv,p)"
+ )
+ : -dpdt
+ )
+ - fvm::laplacian(alphac*turbulence->alphaEff(), he)
+ ==
+ alphac*rho*(U&g)
+ + parcels.Sh(he)
+ + radiation->Sh(thermo)
+ + combustion->Sh()
+ + fvOptions(alphac,rho, he)
+ );
+
+ EEqn.relax();
+
+ fvOptions.constrain(EEqn);
+
+ EEqn.solve();
+
+ fvOptions.correct(he);
+
+ thermo.correct();
+ radiation->correct();
+
+ Info<< "T gas min/max = " << min(T).value() << ", "
+ << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/Make/files b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/files
new file mode 100644
index 000000000..dcdf1f14b
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/files
@@ -0,0 +1,3 @@
+ThermoMPPICFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/ThermoMPPICFoam
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/Make/options b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/options
new file mode 100644
index 000000000..8596ce3bb
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/options
@@ -0,0 +1,59 @@
+DEV_PATH=../../libs
+
+EXE_INC = \
+ -I. \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I${LIB_SRC}/meshTools/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+ -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+ -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+ -I$(LIB_SRC)/ODE/lnInclude \
+ -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+ -I$(DEV_PATH)/combustionModels_POSTECH/lnInclude \
+ -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude \
+ -I$(DEV_PATH)/radiationModels_POSTECH/lnInclude \
+ -I$(DEV_PATH)/intermediate_POSTECH/lnInclude
+
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lturbulenceModels \
+ -lcompressibleTurbulenceModels \
+ -llagrangian \
+ -llagrangianTurbulence \
+ -lspecie \
+ -lcompressibleTransportModels \
+ -lfluidThermophysicalModels \
+ -lliquidProperties \
+ -lliquidMixtureProperties \
+ -lsolidProperties \
+ -lsolidMixtureProperties \
+ -lthermophysicalFunctions \
+ -lreactionThermophysicalModels \
+ -lSLGThermo \
+ -lODE \
+ -lregionModels \
+ -lsurfaceFilmModels \
+ -lfvOptions \
+ -lsampling \
+ -L$(FOAM_USER_LIBBIN)\
+ -lcombustionModels_POSTECH \
+ -llagrangianIntermediate_POSTECH \
+ -lchemistryModel_POSTECH \
+ -lradiationModels_POSTECH
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C b/applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C
new file mode 100644
index 000000000..9ef3a3a76
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-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
+ ThermoMPPICFoam
+
+Description
+ Transient solver for compressible, turbulent flow with a reacting,
+ multiphase particle cloud, and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "basicThermoKinematicMPPICCloud.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "fvOptions.H"
+#include "SLGThermo.H"
+#include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "postProcess.H"
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createMesh.H"
+ #include "createControl.H"
+ #include "createTimeControls.H"
+ #include "createRDeltaT.H"
+ #include "createFields.H"
+ #include "createFieldRefs.H"
+ #include "createFvOptions.H"
+ #include "initContinuityErrs.H"
+
+ turbulence->validate();
+
+ #include "compressibleCourantNo.H"
+ #include "setInitialDeltaT.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+ #include "readTimeControls.H"
+ #include "compressibleCourantNo.H"
+ #include "setDeltaT.H"
+
+ runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ parcels.evolve();
+
+ // Update continuous phase volume fraction field
+ alphac = max(1.0 - parcels.theta(), alphacMin);
+ alphac.correctBoundaryConditions();
+ alphacf = fvc::interpolate(alphac);
+ alphacPhi = alphacf*phi;
+
+ fvVectorMatrix cloudSU(parcels.SU(U));
+ volVectorField cloudVolSUSu
+ (
+ IOobject
+ (
+ "cloudVolSUSu",
+ runTime.timeName(),
+ mesh
+ ),
+ mesh,
+ dimensionedVector
+ (
+ "0",
+ cloudSU.dimensions()/dimVolume,
+ vector::zero
+ ),
+ zeroGradientFvPatchVectorField::typeName
+ );
+
+ cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
+ cloudVolSUSu.correctBoundaryConditions();
+ cloudSU.source() = vector::zero;
+
+ #include "rhoEqn.H"
+
+ // --- Pressure-velocity PIMPLE corrector loop
+ while (pimple.loop())
+ {
+ #include "UEqn.H"
+ #include "YEqn.H"
+ #include "EEqn.H"
+
+ // --- Pressure corrector loop
+ while (pimple.correct())
+ {
+ #include "pEqn.H"
+ }
+
+ if (pimple.turbCorr())
+ {
+ turbulence->correct();
+ }
+ }
+
+ rho = thermo.rho();
+
+ runTime.write();
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H b/applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H
new file mode 100644
index 000000000..4d4bdcb65
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H
@@ -0,0 +1,43 @@
+ fvVectorMatrix UEqn
+ (
+ fvm::ddt(alphac, rho, U)
+ + fvm::div(alphacPhi, U)
+ - fvm::laplacian(alphac*turbulence->muEff(),U)
+ - fvc::div(alphac*turbulence->muEff()*dev2(fvc::grad(U)().T()))
+ - fvm::Sp(fvc::ddt(alphac,rho) + fvc::div(alphacPhi), U)
+ ==
+ cloudSU
+ + fvOptions(alphac,rho, U)
+ );
+
+ UEqn.relax();
+
+ volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+
+ surfaceScalarField phicForces
+ (
+ (fvc::interpolate(rAU*cloudVolSUSu) & mesh.Sf())
+ + alphacf*rhorAUf*(g & mesh.Sf())
+ );
+
+ fvOptions.constrain(UEqn);
+
+ if (pimple.momentumPredictor())
+ {
+ solve
+ (
+ UEqn
+ ==
+ fvc::reconstruct
+ (
+ phicForces/rAUf
+ )
+ - fvc::grad(p)
+ );
+
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
+ }
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H b/applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H
new file mode 100644
index 000000000..e83392103
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H
@@ -0,0 +1,53 @@
+tmp> mvConvection
+(
+ fv::convectionScheme::New
+ (
+ mesh,
+ fields,
+ phi,
+ mesh.divScheme("div(phi,Yi_h)")
+ )
+);
+
+{
+ combustion->correct();
+ dQ = combustion->dQ();
+ label inertIndex = -1;
+ volScalarField Yt(0.0*Y[0]);
+
+ forAll(Y, i)
+ {
+ if (Y[i].name() != inertSpecie)
+ {
+ volScalarField& Yi = Y[i];
+
+ fvScalarMatrix YEqn
+ (
+ fvm::ddt(alphac,rho, Yi)
+ + mvConvection->fvmDiv(alphacPhi, Yi)
+ - fvm::laplacian(alphac*turbulence->muEff(), Yi)
+ ==
+ combustion->R(Yi)
+ + fvOptions(alphac,rho, Yi)
+ );
+
+ YEqn.relax();
+
+ fvOptions.constrain(YEqn);
+
+ YEqn.solve(mesh.solver("Yi"));
+
+ fvOptions.correct(Yi);
+
+ Yi.max(0.0);
+ Yt += Yi;
+ }
+ else
+ {
+ inertIndex = i;
+ }
+ }
+
+ Y[inertIndex] = scalar(1) - Yt;
+ Y[inertIndex].max(0.0);
+}
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H b/applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H
new file mode 100644
index 000000000..451ff4793
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H
@@ -0,0 +1,26 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicThermoKinematicMPPICCloud parcels
+(
+ "reactingCloud",
+ rho,
+ U,
+ g,
+ slgThermo
+);
+
+
+// Continuouse phase volume fraction lower limit
+scalar alphacMin
+(
+ 1.0
+ - readScalar
+ (
+ parcels.particleProperties().subDict("constantProperties").lookup("alphaMax")
+ )
+);
+
+// Update alphac from the particle locations
+alphac = max(1.0 - parcels.theta(), alphacMin);
+alphac.correctBoundaryConditions();
+alphacf = fvc::interpolate(alphac);
+alphacPhi = alphacf*phi;
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H b/applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H
new file mode 100644
index 000000000..5842906a6
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H
@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H b/applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H
new file mode 100644
index 000000000..b6f3dfdef
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H
@@ -0,0 +1,158 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Creating combustion model\n" << endl;
+
+autoPtr combustion
+(
+ combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = combustion->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList& Y = composition.Y();
+
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+ FatalErrorInFunction
+ << "Specified inert specie '" << inertSpecie << "' not found in "
+ << "species list. Available species:" << composition.species()
+ << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ thermo.rho()
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+dimensionedScalar rhoMax
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMax",
+ pimple.dict(),
+ dimDensity,
+ GREAT
+ )
+);
+
+dimensionedScalar rhoMin
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMin",
+ pimple.dict(),
+ dimDensity,
+ 0
+ )
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr turbulence
+(
+ compressible::turbulenceModel::New
+ (
+ rho,
+ U,
+ phi,
+ thermo
+ )
+);
+
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+ IOobject
+ (
+ "dpdt",
+ runTime.timeName(),
+ mesh
+ ),
+ mesh,
+ dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme::fieldTable fields;
+
+forAll(Y, i)
+{
+ fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField dQ
+(
+ IOobject
+ (
+ "dQ",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);
+
+volScalarField alphac
+(
+ IOobject
+ (
+ "alphac",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("0", dimless, 1)
+);
+
+surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
+surfaceScalarField alphacPhi("alphacPhi", alphacf*phi);
+
+volScalarField alphacrho("alphacrho", alphac*rho);
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H b/applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H
new file mode 100644
index 000000000..053fdcf16
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H
@@ -0,0 +1,81 @@
+{
+ rho = thermo.rho();
+
+ // Thermodynamic density needs to be updated by psi*d(p) after the
+ // pressure solution - done in 2 parts. Part 1:
+ thermo.rho() -= psi*p;
+
+ rhorAUf= fvc::interpolate(rho*rAU);
+ surfaceScalarField rhof("rhof", fvc::interpolate(rho));
+
+ volVectorField HbyA("HbyA", U);
+ HbyA = rAU*UEqn.H();
+
+ surfaceScalarField phiHbyA
+ (
+ "phiHbyA",
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
+ + phicForces*rhof
+
+ )
+ );
+
+ MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+ // Update the pressure BCs to ensure flux consistency
+ constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+ fvScalarMatrix pDDtEqn
+ (
+ fvc::ddt(alphac,rho) + alphac*psi*correction(fvm::ddt(p))
+ + fvc::div(alphacf*phiHbyA)
+ ==
+ fvOptions(alphac*psi, p, rho.name())
+ );
+
+ while (pimple.correctNonOrthogonal())
+ {
+ fvScalarMatrix pEqn
+ (
+ pDDtEqn
+ - fvm::laplacian(alphacf*rhorAUf, p)
+ );
+
+ pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+ if (pimple.finalNonOrthogonalIter())
+ {
+
+ phi = phiHbyA + pEqn.flux()/alphacf;
+
+ alphacPhi = alphacf*phi;
+
+ p.relax();
+
+ thermo.rho() += psi*p;
+
+ #include "rhoEqn.H"
+ #include "compressibleContinuityErrs.H"
+
+ U = HbyA + rAU*fvc::reconstruct(phicForces/rAUf) - rAU*fvc::grad(p);
+
+ }
+ }
+
+ U.correctBoundaryConditions();
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
+
+ if (thermo.dpdt())
+ {
+ dpdt = fvc::ddt(p);
+ }
+
+ rho = thermo.rho();
+ rho = max(rho, rhoMin);
+ rho = min(rho, rhoMax);
+
+ Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H b/applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H
new file mode 100644
index 000000000..36203c48b
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2015 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 .
+
+Global
+ rhoEqn
+
+Description
+ Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+ fvScalarMatrix rhoEqn
+ (
+ fvm::ddt(alphac,rho)
+ + fvc::div(alphacPhi)
+ ==
+ fvOptions(alphac,rho)
+ );
+
+ rhoEqn.solve();
+
+ fvOptions.correct(rho);
+
+ Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+ << endl;
+
+ rho = max(rho, rhoMin);
+ rho = min(rho, rhoMax);
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H b/applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H
new file mode 100644
index 000000000..63398fb00
--- /dev/null
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-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 .
+
+\*---------------------------------------------------------------------------*/
+
+{
+ volScalarField& rDeltaT = trDeltaT.ref();
+
+ const dictionary& pimpleDict = pimple.dict();
+
+ // Maximum flow Courant number
+ scalar maxCo(readScalar(pimpleDict.lookup("maxCo")));
+
+ // Maximum time scale
+ scalar maxDeltaT(pimpleDict.lookupOrDefault("maxDeltaT", GREAT));
+
+ // Smoothing parameter (0-1) when smoothing iterations > 0
+ scalar rDeltaTSmoothingCoeff
+ (
+ pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.1)
+ );
+
+ // Damping coefficient (1-0)
+ scalar rDeltaTDampingCoeff
+ (
+ pimpleDict.lookupOrDefault("rDeltaTDampingCoeff", 0.2)
+ );
+
+ // Maximum change in cell temperature per iteration
+ // (relative to previous value)
+ scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
+
+
+ Info<< "Time scales min/max:" << endl;
+
+ // Cache old reciprocal time scale field
+ volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+ // Flow time scale
+ {
+ rDeltaT.ref() =
+ (
+ fvc::surfaceSum(mag(phi))()()
+ /((2*maxCo)*mesh.V()*rho())
+ );
+
+ // Limit the largest time scale
+ rDeltaT.max(1/maxDeltaT);
+
+ Info<< " Flow = "
+ << gMin(1/rDeltaT.primitiveField()) << ", "
+ << gMax(1/rDeltaT.primitiveField()) << endl;
+ }
+
+ // Reaction source time scale
+ {
+ volScalarField::Internal rDeltaTT
+ (
+ mag
+ (
+ parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+ + combustion->Sh()()
+ )
+ /(
+ alphaTemp
+ *rho()
+ *thermo.Cp()()()
+ *T()
+ )
+ );
+
+ Info<< " Temperature = "
+ << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+ << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+
+ rDeltaT.ref() = max
+ (
+ rDeltaT(),
+ rDeltaTT
+ );
+ }
+
+ // Update tho boundary values of the reciprocal time-step
+ rDeltaT.correctBoundaryConditions();
+
+ // Spatially smooth the time scale field
+ if (rDeltaTSmoothingCoeff < 1.0)
+ {
+ fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+ }
+
+ // Limit rate of change of time scale
+ // - reduce as much as required
+ // - only increase at a fraction of old time scale
+ if
+ (
+ rDeltaTDampingCoeff < 1.0
+ && runTime.timeIndex() > runTime.startTimeIndex() + 1
+ )
+ {
+ rDeltaT = max
+ (
+ rDeltaT,
+ (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+ );
+ }
+
+ Info<< " Overall = "
+ << gMin(1/rDeltaT.primitiveField())
+ << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H
new file mode 100644
index 000000000..d4d686e35
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H
@@ -0,0 +1,32 @@
+{
+ volScalarField& he = thermo.he();
+
+ fvScalarMatrix EEqn
+ (
+ mvConvection->fvmDiv(phi, he)
+ + (
+ he.name() == "e"
+ ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+ : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+ )
+ - fvm::laplacian(turbulence->alphaEff(), he)
+ ==
+ parcels.Sh(he)
+ + radiation->Sh(thermo)
+ + combustion->Sh()
+ + fvOptions(rho, he)
+ );
+
+ EEqn.relax();
+
+ fvOptions.constrain(EEqn);
+
+ EEqn.solve();
+
+ fvOptions.correct(he);
+ thermo.correct();
+ radiation->correct();
+
+ Info<< "T gas min/max = " << min(T).value() << ", "
+ << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files
new file mode 100644
index 000000000..a48b2dd78
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files
@@ -0,0 +1,3 @@
+simpleCoalcombustionFoam.C
+
+EXE = $(FOAM_APPBIN)/simpleCoalcombustionFoam
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options
new file mode 100644
index 000000000..f39888db2
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options
@@ -0,0 +1,56 @@
+EXE_INC = \
+ -I. \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I${LIB_SRC}/meshTools/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+ -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+ -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+ -I$(LIB_SRC)/ODE/lnInclude \
+ -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+ -I$(LIB_SRC)/combustionModels/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lturbulenceModels \
+ -lcompressibleTurbulenceModels \
+ -llagrangian \
+ -llagrangianIntermediate \
+ -llagrangianTurbulence \
+ -lspecie \
+ -lcompressibleTransportModels \
+ -lfluidThermophysicalModels \
+ -lliquidProperties \
+ -lliquidMixtureProperties \
+ -lsolidProperties \
+ -lsolidMixtureProperties \
+ -lthermophysicalFunctions \
+ -lreactionThermophysicalModels \
+ -lSLGThermo \
+ -lchemistryModel \
+ -lradiationModels \
+ -lODE \
+ -lregionModels \
+ -lsurfaceFilmModels \
+ -lcombustionModels \
+ -lfvOptions \
+ -lsampling
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H
new file mode 100644
index 000000000..a29d4b8c5
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H
@@ -0,0 +1,21 @@
+ MRF.correctBoundaryVelocity(U);
+
+ tmp tUEqn
+ (
+ fvm::div(phi, U)
+ + MRF.DDt(rho, U)
+ + turbulence->divDevRhoReff(U)
+ ==
+ rho()*g
+ + parcels.SU(U)
+ + fvOptions(rho, U)
+ );
+ fvVectorMatrix& UEqn = tUEqn.ref();
+
+ UEqn.relax();
+
+ fvOptions.constrain(UEqn);
+
+ solve(UEqn == -fvc::grad(p));
+
+ fvOptions.correct(U);
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H
new file mode 100644
index 000000000..a16f67ee3
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H
@@ -0,0 +1,53 @@
+tmp> mvConvection
+(
+ fv::convectionScheme::New
+ (
+ mesh,
+ fields,
+ phi,
+ mesh.divScheme("div(phi,Yi_h)")
+ )
+);
+
+{
+ combustion->correct();
+ dQ = combustion->dQ();
+ label inertIndex = -1;
+ volScalarField Yt(0.0*Y[0]);
+
+ forAll(Y, i)
+ {
+ if (Y[i].name() != inertSpecie)
+ {
+ volScalarField& Yi = Y[i];
+
+ fvScalarMatrix YEqn
+ (
+ mvConvection->fvmDiv(phi, Yi)
+ - fvm::laplacian(turbulence->muEff(), Yi)
+ ==
+ parcels.Srho(i)
+ + combustion->R(Yi)
+ + fvOptions(rho, Yi)
+ );
+
+ YEqn.relax();
+
+ fvOptions.constrain(YEqn);
+
+ YEqn.solve(mesh.solver("Yi"));
+
+ fvOptions.correct(Yi);
+
+ Yi.max(0.0);
+ Yt += Yi;
+ }
+ else
+ {
+ inertIndex = i;
+ }
+ }
+
+ Y[inertIndex] = scalar(1) - Yt;
+ Y[inertIndex].max(0.0);
+}
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H
new file mode 100644
index 000000000..f80b43e82
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H
@@ -0,0 +1,8 @@
+Info<< "\nConstructing reacting cloud" << endl;
+coalCloudList parcels
+(
+ rho,
+ U,
+ g,
+ slgThermo
+);
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H
new file mode 100644
index 000000000..5842906a6
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H
@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H
new file mode 100644
index 000000000..03413049b
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H
@@ -0,0 +1,123 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Creating combustion model\n" << endl;
+
+autoPtr combustion
+(
+ combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = combustion->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList& Y = composition.Y();
+
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+ FatalErrorInFunction
+ << "Specified inert specie '" << inertSpecie << "' not found in "
+ << "species list. Available species:" << composition.species()
+ << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ thermo.rho()
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+#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
+ )
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr turbulence
+(
+ compressible::turbulenceModel::New
+ (
+ rho,
+ U,
+ phi,
+ thermo
+ )
+);
+
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme::fieldTable fields;
+
+forAll(Y, i)
+{
+ fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField dQ
+(
+ IOobject
+ (
+ "dQ",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H b/applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H
new file mode 100644
index 000000000..055eff6f0
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H
@@ -0,0 +1,57 @@
+{
+ // Thermodynamic density needs to be updated by psi*d(p) after the
+ // pressure solution - done in 2 parts. Part 1:
+ thermo.rho() -= psi*p;
+
+ volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+ tUEqn.clear();
+ surfaceScalarField phiHbyA
+ (
+ "phiHbyA",
+ fvc::interpolate(rho)*fvc::flux(HbyA)
+ );
+
+ MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+ // Update the pressure BCs to ensure flux consistency
+ constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+ while (simple.correctNonOrthogonal())
+ {
+ fvScalarMatrix pEqn
+ (
+ fvc::div(phiHbyA)
+ - fvm::laplacian(rhorAUf, p)
+ ==
+ parcels.Srho()
+ + fvOptions(psi, p, rho.name())
+ );
+
+ pEqn.solve();
+
+ if (simple.finalNonOrthogonalIter())
+ {
+ phi = phiHbyA + pEqn.flux();
+ }
+ }
+
+ p.relax();
+
+ // Second part of thermodynamic density update
+ thermo.rho() += psi*p;
+
+ #include "compressibleContinuityErrs.H"
+
+ U = HbyA - rAU*fvc::grad(p);
+ U.correctBoundaryConditions();
+ fvOptions.correct(U);
+
+ rho = thermo.rho();
+ rho = max(rho, rhoMin);
+ rho = min(rho, rhoMax);
+ rho.relax();
+
+ Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C b/applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C
new file mode 100644
index 000000000..49ec91e64
--- /dev/null
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ simpleCoalcombustionFoam
+
+Description
+ Steady state solver for compressible, turbulent flow with reacting,
+ multiphase particle clouds and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "coalCloudList.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "IOporosityModelList.H"
+#include "fvOptions.H"
+#include "SLGThermo.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "postProcess.H"
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createMesh.H"
+ #include "createControl.H"
+ #include "createFields.H"
+ #include "createFieldRefs.H"
+ #include "createFvOptions.H"
+ #include "initContinuityErrs.H"
+
+ turbulence->validate();
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (simple.loop())
+ {
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ parcels.evolve();
+
+ // --- Pressure-velocity SIMPLE corrector loop
+ {
+ #include "UEqn.H"
+ #include "YEqn.H"
+ #include "EEqn.H"
+ #include "pEqn.H"
+ }
+
+ turbulence->correct();
+
+ runTime.write();
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
index 6aa801f69..4e454f228 100644
--- a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
+++ b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
@@ -34,6 +34,7 @@ License
#include "COxidationKineticDiffusionLimitedRate.H"
#include "COxidationHurtMitchell.H"
#include "COxidationMurphyShaddix.H"
+#include "COxidationKD_CO.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -47,6 +48,7 @@ License
); \
makeSurfaceReactionModelType(COxidationIntrinsicRate, CloudType); \
makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType); \
+ makeSurfaceReactionModelType(COxidationKD_CO, CloudType); \
makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C
new file mode 100644
index 000000000..3982f1ed9
--- /dev/null
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "COxidationKD_CO.H"
+#include "mathematicalConstants.H"
+#include "thermodynamicConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template
+Foam::COxidationKD_CO::
+COxidationKD_CO
+(
+ const dictionary& dict,
+ CloudType& owner
+)
+:
+ SurfaceReactionModel(dict, owner, typeName),
+ Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+ C1_(readScalar(this->coeffDict().lookup("C1"))),
+ C2_(readScalar(this->coeffDict().lookup("C2"))),
+ E_(readScalar(this->coeffDict().lookup("E"))),
+ CsLocalId_(-1),
+ O2GlobalId_(owner.composition().carrierId("O2")),
+ COGlobalId_(owner.composition().carrierId("CO")),
+ WC_(0.0),
+ WO2_(0.0),
+ HcCO_(0.0)
+{
+ // Determine Cs ids
+ label idSolid = owner.composition().idSolid();
+ CsLocalId_ = owner.composition().localId(idSolid, "C");
+
+ // Set local copies of thermo properties
+ WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+ const scalar WCO = owner.thermo().carrier().W(COGlobalId_);
+ WC_ = WCO - 0.5*WO2_;
+
+ HcCO_ = owner.thermo().carrier().Hc(COGlobalId_);
+
+ const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+ const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+ Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+}
+
+
+template
+Foam::COxidationKD_CO::
+COxidationKD_CO
+(
+ const COxidationKD_CO& srm
+)
+:
+ SurfaceReactionModel(srm),
+ Sb_(srm.Sb_),
+ C1_(srm.C1_),
+ C2_(srm.C2_),
+ E_(srm.E_),
+ CsLocalId_(srm.CsLocalId_),
+ O2GlobalId_(srm.O2GlobalId_),
+ COGlobalId_(srm.COGlobalId_),
+ WC_(srm.WC_),
+ WO2_(srm.WO2_),
+ HcCO_(srm.HcCO_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template
+Foam::COxidationKD_CO::
+~COxidationKD_CO()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template
+Foam::scalar Foam::COxidationKD_CO::calculate
+(
+ const scalar dt,
+ const label cellI,
+ const scalar d,
+ const scalar T,
+ const scalar Tc,
+ const scalar pc,
+ const scalar rhoc,
+ const scalar mass,
+ const scalarField& YGas,
+ const scalarField& YLiquid,
+ const scalarField& YSolid,
+ const scalarField& YMixture,
+ const scalar N,
+ scalarField& dMassGas,
+ scalarField& dMassLiquid,
+ scalarField& dMassSolid,
+ scalarField& dMassSRCarrier
+) const
+{
+ // Fraction of remaining combustible material
+ const label idSolid = CloudType::parcelType::SLD;
+ const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
+
+ // Surface combustion active combustible fraction is consumed
+ if (fComb < SMALL)
+ {
+ return 0.0;
+ }
+
+ const SLGThermo& thermo = this->owner().thermo();
+
+ // Local mass fraction of O2 in the carrier phase
+ const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+ // Diffusion rate coefficient
+ const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
+
+ // Kinetic rate
+ const scalar Rk = C2_*exp(-E_/(constant::thermodynamic::RR*T)); //
+
+ // Particle surface area
+ const scalar Ap = constant::mathematical::pi*sqr(d);
+
+ // Change in C mass [kg]
+ scalar dmC = Ap*rhoc*constant::thermodynamic::RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
+
+ // Limit mass transfer by availability of C
+ dmC = min(mass*fComb, dmC);
+
+ // Molar consumption
+ const scalar dOmega = dmC/WC_;
+
+ // Change in O2 mass [kg]
+ const scalar dmO2 = dOmega*Sb_*WO2_;
+
+ // Mass of newly created CO [kg]
+ const scalar dmCO = dOmega*(WC_ + Sb_*WO2_);
+
+ // Update local particle C mass
+ dMassSolid[CsLocalId_] += dOmega*WC_;
+
+ // Update carrier O2 and CO mass
+ dMassSRCarrier[O2GlobalId_] -= dmO2;
+ dMassSRCarrier[COGlobalId_] += dmCO;
+
+ const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+
+ // carrier sensible enthalpy exchange handled via change in mass
+
+ // Heat of reaction [J]
+//Info<<"Heat of reaction(Char)"<.
+
+Class
+ COxidationKD_CO
+
+Description
+ Kinetic/diffusion limited rate surface reaction model for coal parcels.
+ Limited to:
+
+ C(s) + Sb*O2 -> CO
+
+ where Sb is the stoichiometry of the reaction
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationKD_CO_H
+#define COxidationKD_CO_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template
+class COxidationKD_CO;
+
+/*---------------------------------------------------------------------------*\
+ Class COxidationKD_CO Declaration
+\*---------------------------------------------------------------------------*/
+
+template
+class COxidationKD_CO
+:
+ public SurfaceReactionModel
+{
+ // Private data
+
+ // Model constants
+
+ //- Stoichiometry of reaction
+ const scalar Sb_;
+
+ //- Mass diffusion limited rate constant, C1
+ const scalar C1_;
+
+ //- Kinetics limited rate pre-exponential constant, C2
+ const scalar C2_;
+
+ //- Kinetics limited rate activation energy
+ const scalar E_;
+
+
+ // Addressing
+
+ //- Cs positions in global/local lists
+ label CsLocalId_;
+
+ //- O2 position in global list
+ label O2GlobalId_;
+
+ //- CO2 positions in global list
+ label COGlobalId_;
+
+
+ // Local copies of thermo properties
+
+ //- Molecular weight of C [kg/kmol]
+ scalar WC_;
+
+ //- Molecular weight of O2 [kg/kmol]
+ scalar WO2_;
+
+ //- Formation enthalpy for CO2 [J/kg]
+ scalar HcCO_;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("COxidationKD_CO");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ COxidationKD_CO
+ (
+ const dictionary& dict,
+ CloudType& owner
+ );
+
+ //- Construct copy
+ COxidationKD_CO
+ (
+ const COxidationKD_CO& srm
+ );
+
+ //- Construct and return a clone
+ virtual autoPtr > clone() const
+ {
+ return autoPtr >
+ (
+ new COxidationKD_CO(*this)
+ );
+ }
+
+
+ //- Destructor
+ virtual ~COxidationKD_CO();
+
+
+ // Member Functions
+
+ //- Update surface reactions
+ virtual scalar calculate
+ (
+ const scalar dt,
+ const label cellI,
+ const scalar d,
+ const scalar T,
+ const scalar Tc,
+ const scalar pc,
+ const scalar rhoc,
+ const scalar mass,
+ const scalarField& YGas,
+ const scalarField& YLiquid,
+ const scalarField& YSolid,
+ const scalarField& YMixture,
+ const scalar N,
+ scalarField& dMassGas,
+ scalarField& dMassLiquid,
+ scalarField& dMassSolid,
+ scalarField& dMassSRCarrier
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "COxidationKD_CO.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 9cadf5511..f62ee8394 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -50,6 +50,12 @@ $(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
$(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
+/* thermo kinematic MPPIC parcel sub-models */
+THERMOMPPICPARCEL=$(DERIVEDPARCELS)/basicThermoKinematicMPPICParcel
+$(THERMOMPPICPARCEL)/defineBasicThermoKinematicMPPICParcel.C
+$(THERMOMPPICPARCEL)/makeBasicThermoKinematicMPPICParcelSubmodels.C
+
+
/* bolt-on models */
RADIATION=submodels/addOns/radiation
$(RADIATION)/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
index 702d2454b..f6cc7c950 100644
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@@ -70,12 +70,12 @@ Foam::MPPICCloud::MPPICCloud
const word& cloudName,
const volScalarField& rho,
const volVectorField& U,
- const volScalarField& mu,
const dimensionedVector& g,
+ const SLGThermo& thermo,
bool readFields
)
:
- CloudType(cloudName, rho, U, mu, g, false),
+ CloudType(cloudName, rho, U, g, thermo, false),
packingModel_(NULL),
dampingModel_(NULL),
isotropyModel_(NULL)
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
index 1d34b686e..c49125bd6 100644
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
@@ -43,6 +43,7 @@ SourceFiles
#include "fvMesh.H"
#include "volFields.H"
+#include "ThermoCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@@ -133,8 +134,8 @@ public:
const word& cloudName,
const volScalarField& rho,
const volVectorField& U,
- const volScalarField& mu,
const dimensionedVector& g,
+ const SLGThermo& thermo,
bool readFields = true
);
diff --git a/src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H b/src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H
new file mode 100644
index 000000000..dbf6b7c0a
--- /dev/null
+++ b/src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 .
+
+Class
+ Foam::basicThermoKinematicMPPICCloud
+
+Description
+ Cloud class to introduce Thermo kinematic MPPIC parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicThermoKinematicMPPICCloud_H
+#define basicThermoKinematicMPPICCloud_H
+
+#include "Cloud.H"
+#include "KinematicCloud.H"
+#include "MPPICCloud.H"
+#include "basicThermoKinematicMPPICParcel.H"
+
+#include "ThermoCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ typedef MPPICCloud
+ <
+ ThermoCloud
+ <
+ KinematicCloud
+ <
+ Cloud
+ <
+ basicThermoKinematicMPPICParcel
+ >
+ >
+ >
+ >
+ basicThermoKinematicMPPICCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H
new file mode 100644
index 000000000..fb2d9b5cf
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 .
+
+Class
+ Foam::basicThermoKinematicMPPICParcel
+
+Description
+ Definition of basic Thermo kinematic MPPIC parcel
+
+SourceFiles
+ basicThermoKinematicMPPICParcel.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicThermoKinematicMPPICParcel_H
+#define basicThermoKinematicMPPICParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "MPPICParcel.H"
+
+#include "ThermoParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+ typedef MPPICParcel < ThermoParcel < KinematicParcel < particle > > > basicThermoKinematicMPPICParcel;
+
+ template<>
+ inline bool contiguous()
+ {
+ return true;
+ }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C
new file mode 100644
index 000000000..f71ecafb9
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C
@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicThermoKinematicMPPICParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTemplateTypeNameAndDebug(basicThermoKinematicMPPICParcel, 0);
+ defineTemplateTypeNameAndDebug(Cloud, 0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C
new file mode 100644
index 000000000..771c6b821
--- /dev/null
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C
@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013-2015 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicThermoKinematicMPPICCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// Kinematic sub-models
+#include "makeParcelForces.H"
+#include "makeParcelDispersionModels.H"
+#include "makeParcelInjectionModels.H"
+#include "makeParcelPatchInteractionModels.H"
+#include "makeParcelStochasticCollisionModels.H"
+#include "makeParcelSurfaceFilmModels.H"
+
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
+// Thermo sub-models
+#include "makeParcelHeatTransferModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeParcelCloudFunctionObjects(basicThermoKinematicMPPICCloud);
+
+// Kinematic sub-models
+makeParcelForces(basicThermoKinematicMPPICCloud);
+
+makeParcelDispersionModels(basicThermoKinematicMPPICCloud);
+makeParcelInjectionModels(basicThermoKinematicMPPICCloud);
+makeParcelPatchInteractionModels(basicThermoKinematicMPPICCloud);
+makeParcelStochasticCollisionModels(basicThermoKinematicMPPICCloud);
+makeParcelSurfaceFilmModels(basicThermoKinematicMPPICCloud);
+
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicThermoKinematicMPPICCloud);
+makeMPPICParcelIsotropyModels(basicThermoKinematicMPPICCloud);
+makeMPPICParcelPackingModels(basicThermoKinematicMPPICCloud);
+
+// Thermo sub-models
+makeParcelHeatTransferModels(basicThermoKinematicMPPICCloud);
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/Make/files b/src/thermophysicalModels/radiation/Make/files
index 5f520b2dc..afca9e68e 100644
--- a/src/thermophysicalModels/radiation/Make/files
+++ b/src/thermophysicalModels/radiation/Make/files
@@ -25,6 +25,7 @@ submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmiss
submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
+submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
/* Soot model */
submodels/sootModel/sootModel/sootModel.C
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
new file mode 100644
index 000000000..86fe5fadc
--- /dev/null
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2013 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "wsggmAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "zeroGradientFvPatchFields.H"
+#include "basicSpecieMixture.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ namespace radiation
+ {
+ defineTypeNameAndDebug(wsggmAbsorptionEmission, 0);
+
+ addToRunTimeSelectionTable
+ (
+ absorptionEmissionModel,
+ wsggmAbsorptionEmission,
+ dictionary
+ );
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission
+(
+ const dictionary& dict,
+ const fvMesh& mesh
+)
+:
+ absorptionEmissionModel(dict, mesh),
+ coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+ thermo_(mesh.lookupObject("thermophysicalProperties")),
+ emissivityCoeffs_(coeffsDict_.lookup("emissivityCoeffs")),
+ fittingFactors_(coeffsDict_.lookup("fittingFactors")),
+ pathLength_(coeffsDict_.lookup("pathLength"))
+{
+ if (!isA(thermo_))
+ {
+ FatalErrorIn
+ (
+ "radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission"
+ "("
+ "const dictionary&, "
+ "const fvMesh&"
+ ")"
+ ) << "Model requires a multi-component thermo package"
+ << abort(FatalError);
+ }
+}
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::radiation::wsggmAbsorptionEmission::~wsggmAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::tmp
+Foam::radiation::wsggmAbsorptionEmission::aCont(const label bandI) const
+{
+
+ const basicSpecieMixture& mixture =
+ dynamic_cast(thermo_);
+
+ const volScalarField& T = thermo_.T();
+ const volScalarField& p = thermo_.p();
+
+ tmp ta
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "aCont" + name(bandI),
+ mesh().time().timeName(),
+ mesh(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh(),
+ dimensionedScalar("a", dimless/dimLength, 0.0),
+ zeroGradientFvPatchVectorField::typeName
+ )
+ );
+
+ scalarField& a = ta.ref().primitiveFieldRef();
+
+ label indexCO2, indexH2O;
+ scalar wCO2, wH2O;
+ forAll(mixture.Y(), specieI)
+ {
+ if(mixture.Y(specieI).name() == "CO2")
+ {
+ indexCO2 = specieI;
+ wCO2 = mixture.W(indexCO2);
+ }
+ else if(mixture.Y(specieI).name() == "H2O")
+ {
+ indexH2O = specieI;
+ wH2O = mixture.W(indexH2O);
+ }
+ }
+
+ forAll(a,celli)
+ {
+ scalar invWt = 0.0;
+ forAll(mixture.Y(), specieI)
+ {
+ invWt += mixture.Y(specieI)[celli]/mixture.W(specieI);
+ }
+ scalar meanWt = 1.0/invWt; //Mean molecular weight [kg/kmole]
+
+ // Sum of partial pressures of all absorbing gases [atm]
+ scalar Pcell_atm = paToAtm(meanWt*p[celli]*(mixture.Y(indexCO2)[celli]/wCO2 + mixture.Y(indexH2O)[celli]/wH2O));
+
+ // P*S at Eq.(1) of Ref. [atm*m]
+ scalar presPathLength = Pcell_atm*pathLength_.value();
+
+ // Limit cell temperature [K]
+ scalar Tcell = min(T[celli],2400.0);
+
+ scalar emissivity = 0.0;
+ forAll(emissivityCoeffs_, i) //k_i at Eq.(1) of Ref.
+ {
+ scalar weightingFactor = 0.0;
+ for(label j=0 ; j<4 ; j++)
+ {
+ //a_epsilon_i at Eq.(3) of Ref.
+ weightingFactor += fittingFactors_[i][j]*pow(Tcell, j);
+ }
+ //epsilon at Eq.(1) of Ref.
+ emissivity += weightingFactor * (1.0 - exp(-1.0*emissivityCoeffs_[i] * presPathLength));
+ }
+ emissivity = min(emissivity, 0.99999);
+
+ //Fluent theory guide(ver. 13) Eq.(5-106) [1/m]
+ a[celli] = max(1e-2, (-1.0)*log(1.0-emissivity)/pathLength_.value());
+ }
+
+ ta.ref().correctBoundaryConditions();
+
+ return ta;
+}
+
+
+Foam::tmp
+Foam::radiation::wsggmAbsorptionEmission::eCont(const label bandI) const
+{
+ return aCont(bandI);
+}
+
+
+Foam::tmp
+Foam::radiation::wsggmAbsorptionEmission::ECont(const label bandI) const
+{
+ tmp E
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "ECont" + name(bandI),
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+ )
+ );
+
+ return E;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H
new file mode 100644
index 000000000..d57f1ac1c
--- /dev/null
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 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 .
+
+Class
+ Foam::radiation::wsggmAbsorptionEmission
+
+Description
+ wsggmAbsorptionEmission radiation absorption and emission coefficients
+ for continuous phase
+
+ Ref. T.F Smith, Z.F. Shen and J.N. Friedman
+ Evaluation of Coefficients for the Weighted Sum of Gray Gases Model,
+ Transactions of the ASME, Vol. 104, pp. 602-608
+
+ Valid when Ptotal = 1atm, 0.001 < P*S < 10.0atm.m, 600 < T < 2400K
+ (where P*S = sum of partial pressures of absorbing gases
+ * mean beam length)
+
+ i.e. dictionary
+ \verbatim
+ wsggmAbsorptionEmissionCoeffs
+ {
+ pathLength pathLength [ 0 1 0 0 0 0 0] 0.251; //3.6V/S
+ // Pw/Pc = 2.0
+ emissivityCoeffs 3(0.4201 6.516 131.9);
+ fittingFactors
+ 3
+ (
+ 4(6.508e-1 -5.551e-4 3.029e-7 -5.353e-11)
+ 4(-0.2504e-1 6.112e-4 -3.882e-7 6.528e-11)
+ 4(2.718e-1 -3.118e-4 1.221e-7 -1.612e-11)
+ );
+ }
+ \endverbatim
+
+
+SourceFiles
+ wsggmAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wsggmAbsorptionEmission_H
+#define wsggmAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "fluidThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+ Class wsggmAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class wsggmAbsorptionEmission
+:
+ public absorptionEmissionModel
+{
+public:
+
+private:
+
+ // Private data
+
+ //- Absorption model dictionary
+ dictionary coeffsDict_;
+
+ //- SLG thermo package
+ const fluidThermo& thermo_;
+
+ //- absorption coefficitent for individual grey gases
+ scalarList emissivityCoeffs_;
+
+ //- fitting Factors
+ scalarListList fittingFactors_;
+
+ //characteristic path length [m] from dict
+ dimensionedScalar pathLength_;
+
+public:
+
+ //- Runtime type information
+ TypeName("wsggmAbsorptionEmission");
+
+
+ // Constructors
+
+ //- Construct from components
+ wsggmAbsorptionEmission(const dictionary& dict, const fvMesh& mesh);
+
+
+ //- Destructor
+ virtual ~wsggmAbsorptionEmission();
+
+
+ // Member Functions
+
+ // Access
+
+ // Absorption coefficient
+
+ //- Absorption coefficient for continuous phase
+ tmp aCont(const label bandI = 0) const;
+
+
+ // Emission coefficient
+
+ //- Emission coefficient for continuous phase
+ tmp eCont(const label bandI = 0) const;
+
+
+ // Emission contribution
+
+ //- Emission contribution for continuous phase
+ tmp ECont(const label bandI = 0) const;
+
+
+ // Member Functions
+
+ inline bool isGrey() const
+ {
+ return true;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //