diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwclean b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwclean
new file mode 100755
index 0000000000..c724941863
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+set -x
+
+wclean libso phasesSystem
+wclean libso meltingEvaporationModel
+wclean libso CompressibleMultiPhaseTurbulenceModels
+wclean libso laserDTRM
+wclean
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwmake b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwmake
new file mode 100755
index 0000000000..60693c27bc
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwmake
@@ -0,0 +1,13 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+set -x
+
+wmakeLnInclude meltingEvaporationModel
+wmake libso phasesSystem
+wmake libso meltingEvaporationModel
+wmake libso CompressibleMultiPhaseTurbulenceModels
+wmake libso laserDTRM
+wmake
+
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/CompressibleMultiPhaseTurbulenceModels.C b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/CompressibleMultiPhaseTurbulenceModels.C
new file mode 100644
index 0000000000..ab4f06d8ef
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/CompressibleMultiPhaseTurbulenceModels.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 OpenCFD Ltd
+ \\/ 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 "CompressibleTurbulenceModel.H"
+#include "compressibleTurbulenceModel.H"
+#include "multiphaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
+
+#include "ThermalDiffusivity.H"
+
+#include "laminarModel.H"
+#include "RASModel.H"
+#include "LESModel.H"
+
+makeBaseTurbulenceModel
+(
+ geometricOneField,
+ volScalarField,
+ compressibleTurbulenceModel,
+ CompressibleTurbulenceModel,
+ ThermalDiffusivity,
+ multiphaseSystem
+);
+
+#define makeLaminarModel(Type) \
+ makeTemplatedLaminarModel \
+ (multiphaseSystemCompressibleTurbulenceModel, laminar, Type)
+
+#define makeRASModel(Type) \
+ makeTemplatedTurbulenceModel \
+ (multiphaseSystemCompressibleTurbulenceModel, RAS, Type)
+
+#define makeLESModel(Type) \
+ makeTemplatedTurbulenceModel \
+ (multiphaseSystemCompressibleTurbulenceModel, LES, Type)
+
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
+
+#include "Smagorinsky.H"
+makeLESModel(Smagorinsky);
+
+#include "kEqn.H"
+makeLESModel(kEqn);
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/files b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/files
new file mode 100644
index 0000000000..c85d7222d2
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/files
@@ -0,0 +1,3 @@
+CompressibleMultiPhaseTurbulenceModels.C
+
+LIB = $(FOAM_LIBBIN)/libCompressibleMultiPhaseTurbulenceModels
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/options b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/options
new file mode 100644
index 0000000000..01969d251b
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/CompressibleMultiPhaseTurbulenceModels/Make/options
@@ -0,0 +1,11 @@
+EXE_INC = \
+ -I../phasesSystem/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/transportModels \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/files b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/files
new file mode 100644
index 0000000000..1e257d9c1a
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/files
@@ -0,0 +1,3 @@
+icoReactingMultiphaseInterFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/icoReactingMultiphaseInterFoam
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/options b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/options
new file mode 100644
index 0000000000..37c05c04aa
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/options
@@ -0,0 +1,34 @@
+
+EXE_INC = \
+ -I./phasesSystem/lnInclude \
+ -I./CompressibleMultiPhaseTurbulenceModels/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/fvOptions/lnInclude\
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude
+
+EXE_LIBS = \
+ -lturbulenceModels \
+ -lcompressibleTurbulenceModels \
+ -lcompressibleTransportModels \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lfvOptions \
+ -lsampling \
+ -lradiationModels \
+ -lfluidThermophysicalModels \
+ -lIncompressibleMultiphaseSystems \
+ -lCompressibleMultiPhaseTurbulenceModels \
+ -lmeltingEvaporationModels \
+ -lsolidThermo \
+ -lsolidSpecie \
+ -ltwoPhaseProperties
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/TEqn.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/TEqn.H
new file mode 100644
index 0000000000..674d03d751
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/TEqn.H
@@ -0,0 +1,102 @@
+{
+ radiation->correct();
+
+ tmp texpSource
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "texpSource",
+ mesh.time().timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("zero", dimTemperature/dimTime, 0),
+ zeroGradientFvPatchScalarField::typeName
+ )
+ );
+ volScalarField& expSource = texpSource.ref();
+
+ tmp tkappaEff
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "kappaEff",
+ mesh.time().timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("zero", sqr(dimLength)/dimTime, 0),
+ zeroGradientFvPatchScalarField::typeName
+ )
+ );
+
+ volScalarField& kappaEff = tkappaEff.ref();
+
+ //const surfaceScalarField rhoTempPhi("phi", rhoPhi/fvc::interpolate(rho));
+
+ const surfaceScalarField rhoTempPhi("phi", fluid.phi());
+
+ const volScalarField divU(fvc::div(phi));
+
+ forAllIter(UPtrList, fluid.phases(), iter)
+ {
+ phaseModel& phase = iter();
+ const volScalarField& alpha = phase;
+ const volScalarField DDtAlpha(fvc::DDt(phi, alpha));
+
+ const volScalarField invCpRho(1.0/phase.rho()/phase.Cp());
+
+ if (fluid.dpdt())
+ {
+ const volScalarField ddtp(fvc::ddt(p));
+ expSource += (DDtAlpha*p + alpha*(p*divU + ddtp))*invCpRho;
+ }
+ else
+ {
+ expSource += (DDtAlpha*p + alpha*(p*divU))*invCpRho;
+ }
+
+ kappaEff += alpha*phase.kappa()*invCpRho;
+
+ DebugVar(max(alpha*phase.kappa()));
+ DebugVar(max(alpha*phase.Cp()));
+ }
+
+ kappaEff += turbulence->nut()/fluid.Prt();
+
+ if (mesh.time().outputTime())
+ {
+ expSource.write();
+ kappaEff.write();
+ }
+
+ //dimensionedScalar S("S", dimEnergy/dimVolume/dimTime, 1.225e8);
+
+ fvScalarMatrix TEqn
+ (
+ fvm::ddt(T)
+ + fvm::div(rhoTempPhi, T)
+ - fvm::laplacian(kappaEff, T, "laplacian(kappa,T)")
+ ==
+ fluid.heatTransfer(T)
+ + expSource
+ + radiation->ST(fluid.Cp()*rho, T)
+// + S/Cp/rho
+ );
+
+ TEqn.relax();
+ TEqn.solve();
+
+ fluid.correct();
+
+ Info<< "min/max(T) = "
+ << min(T).value() << ", " << max(T).value() << endl;
+}
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/UEqn.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/UEqn.H
new file mode 100644
index 0000000000..7f1f0e4423
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/UEqn.H
@@ -0,0 +1,34 @@
+ fvVectorMatrix UEqn
+ (
+ fvm::ddt(rho, U)
+ + fvm::div(rhoPhi, U)
+ + turbulence->divDevRhoReff(U)
+ ==
+ fvOptions(rho, U)
+ );
+
+ UEqn.relax();
+
+ fluid.addInterfacePorosity(UEqn);
+
+ if (pimple.momentumPredictor())
+ {
+
+
+ solve
+ (
+ UEqn
+ ==
+ fvc::reconstruct
+ (
+ (
+ fluid.surfaceTensionForce()
+ - ghf*fvc::snGrad(rho)
+ - fvc::snGrad(p_rgh)
+ ) * mesh.magSf()
+ )
+ );
+
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
+ }
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/YEqns.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/YEqns.H
new file mode 100644
index 0000000000..8be35a2f41
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/YEqns.H
@@ -0,0 +1,76 @@
+{
+ // Semi-implicit mass transfer for species
+ // Initilize dmdt for alpha Eq's for mass transfers driven by species
+ autoPtr
+ massTransferPtr(fluid.massTransfer(T));
+
+ //phaseSystem::massTransferTable& massTransfer(massTransferPtr());
+
+ forAllIter(UPtrList, fluid.phases(), iter)
+ {
+ phaseModel& phase = iter();
+ PtrList& Y = phase.Y();
+ //const surfaceScalarField& alphaPhi = phase.alphaPhi();
+
+ if (!Y.empty())
+ {
+ const volScalarField& alpha = phase;
+
+ label inertIndex = -1;
+ volScalarField Yt(0.0*Y[0]);
+
+ forAll(Y, i)
+ {
+ tmp YiEqn(phase.YiEqn(Y[i]));
+
+ if (YiEqn.valid())
+ {
+ YiEqn.ref() =
+ (
+ YiEqn()
+ - fvm::laplacian
+ (
+ alpha*turbulence->nuEff(),
+ Y[i]
+ )
+ ==
+ // (*massTransfer[Y[i].name()])(/phase.rho()
+ fvc::ddt(alpha)
+ *pos
+ (
+ fluid.dmdtYi(Y[i].name())
+ - dimensionedScalar("zero", dimDensity/dimTime, 1e-3)
+ )
+ //fluid.dmdtYi(Y[i].name())/phase.rho()
+ //explicit mass sources (P or T)
+ );
+
+ YiEqn->relax();
+ YiEqn->solve(mesh.solver("Yi"));
+ Y[i].max(0.0);
+ Y[i].min(1.0);
+ Yt += Y[i];
+ }
+ else
+ {
+ inertIndex = i;
+ }
+
+// if (mesh.time().outputTime())
+// {
+// volScalarField dmdtYi("dmdtYi", pos(fluid.dmdtYi(Y[i].name())));
+// dmdtYi.write();
+// }
+
+ Info << "Min/Max : " << min(Y[i]) << " " << max(Y[i]) << endl;
+ Info<< "Max dmdtYi : "
+ << max(fluid.dmdtYi(Y[i].name())().internalField()) << endl;
+ Info<< "Min dmdtYi : "
+ << min(fluid.dmdtYi(Y[i].name())().internalField()) << endl;
+ }
+ Y[inertIndex] = scalar(1) - Yt;
+ Y[inertIndex].max(0.0);
+ }
+ }
+}
+
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/alphaCourantNo.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/alphaCourantNo.H
new file mode 100644
index 0000000000..b9da86e23c
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/alphaCourantNo.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2014 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
+ alphaCourantNo
+
+Description
+ Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+ readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+ scalarField sumPhi
+ (
+ fluid.nearInterface()().internalField()
+ *fvc::surfaceSum(mag(phi))().internalField()
+ );
+
+ alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+ meanAlphaCoNum =
+ 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+ << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/createFields.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/createFields.H
new file mode 100644
index 0000000000..893608fc95
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/createFields.H
@@ -0,0 +1,137 @@
+ Info<< "Reading field p_rgh\n" << endl;
+ volScalarField p_rgh
+ (
+ IOobject
+ (
+ "p_rgh",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+ Info<< "Reading field U\n" << endl;
+ volVectorField U
+ (
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+
+ Info<< "Calculating field g.h\n" << endl;
+ volScalarField gh("gh", g & mesh.C());
+ surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+ volScalarField p
+ (
+ IOobject
+ (
+ "p",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ p_rgh //+ rho*gh
+ );
+
+ Info<< "Creating multiphaseSystem\n" << endl;
+ autoPtr fluidPtr = multiphaseSystem::New(mesh);
+
+ multiphaseSystem& fluid = fluidPtr();
+
+ //volScalarField& e = fluid.he();
+ volScalarField& T = fluid.T();
+
+ // Need to store rho for ddt(rho, U)
+ volScalarField rho
+ (
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ fluid.rho()
+ );
+ rho.oldTime();
+
+ // Update p using fluid.rho()
+ p = p_rgh + rho*gh;
+
+ label pRefCell = 0;
+ scalar pRefValue = 0.0;
+ setRefCell
+ (
+ p,
+ p_rgh,
+ pimple.dict(),
+ pRefCell,
+ pRefValue
+ );
+
+ if (p_rgh.needReference())
+ {
+ p += dimensionedScalar
+ (
+ "p",
+ p.dimensions(),
+ pRefValue - getRefCellValue(p, pRefCell)
+ );
+ p_rgh = p - rho*gh;
+ }
+
+ // Total volumetric flux
+ surfaceScalarField& phi = fluid.phi();
+
+ // Mass flux
+ surfaceScalarField& rhoPhi = fluid.rhoPhi();
+
+ // Construct incompressible turbulence model
+ autoPtr > turbulence
+ (
+ CompressibleTurbulenceModel::New
+ (
+ rho,
+ U,
+ rhoPhi,
+ fluid
+ )
+ );
+
+ // Creating radiation model
+ autoPtr radiation
+ (
+ radiation::radiationModel::New(T)
+ );
+/*
+ Info<< "Calculating field kappaEff\n" << endl;
+ volScalarField kappaEff
+ (
+ IOobject
+ (
+ "kappaEff",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ fluid.kappa()
+ );
+
+ kappaEff.correctBoundaryConditions();
+*/
+
+ Info<< "Creating field kinetic energy K\n" << endl;
+ volScalarField K("K", 0.5*magSqr(U));
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/icoReactingMultiphaseInterFoam.C b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/icoReactingMultiphaseInterFoam.C
new file mode 100644
index 0000000000..a54ec210c0
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/icoReactingMultiphaseInterFoam.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
+ \\/ 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
+ icoReactingMultiphaseInterFoam
+
+Group
+ grpMultiphaseSolvers
+
+Description
+ Solver for n incompressible, non-isothermal immiscible fluids with
+ phase-change (evaporation-condensation). Uses a VOF (volume of fluid)
+ phase-fraction based interface capturing approach.
+
+ The momentum, energy and other fluid properties are of the "mixture" and a
+ single momentum equation is solved.
+
+ Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "multiphaseSystem.H"
+#include "turbulentFluidThermoModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+#include "radiationModel.H"
+#include "HashPtrTable.H"
+#include "fvcDDt.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createMesh.H"
+
+ pimpleControl pimple(mesh);
+
+ #include "readGravitationalAcceleration.H"
+ #include "createFields.H"
+ #include "createFvOptions.H"
+ #include "createTimeControls.H"
+ #include "CourantNo.H"
+ #include "alphaCourantNo.H"
+ #include "setInitialDeltaT.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+ #include "readTimeControls.H"
+ #include "CourantNo.H"
+ #include "alphaCourantNo.H"
+ #include "setDeltaT.H"
+
+ runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ fluid.solve();
+ rho = fluid.rho();
+
+ // --- Pressure-velocity PIMPLE corrector loop
+ while (pimple.loop())
+ {
+ solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+ #include "UEqn.H"
+ #include "YEqns.H"
+ #include "TEqn.H"
+ // --- Pressure corrector loop
+ while (pimple.correct())
+ {
+ #include "pEqn.H"
+ }
+
+ if (pimple.turbCorr())
+ {
+ turbulence->correct();
+ }
+ }
+
+ rho = fluid.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/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C
new file mode 100644
index 0000000000..54df98a54e
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C
@@ -0,0 +1,256 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd
+ \\/ 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 "DTRMParticle.H"
+#include "constants.H"
+#include "physicoChemicalConstants.H"
+
+using namespace Foam::constant;
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::DTRMParticle::DTRMParticle
+(
+ const polyMesh& mesh,
+ const vector& position,
+ const vector& targetPosition,
+ const scalar I,
+ const label cellI,
+ const scalar dA,
+ const label reflectedId,
+ const scalar Imin,
+ bool doCellFacePt
+)
+:
+ particle(mesh, position, cellI, doCellFacePt),
+ p0_(position),
+ p1_(targetPosition),
+ I0_(I),
+ I_(I),
+ dA_(dA),
+ reflectedId_(reflectedId),
+ Imin_(Imin)
+{}
+
+
+Foam::DTRMParticle::DTRMParticle(const DTRMParticle& p)
+:
+ particle(p),
+ p0_(p.p0_),
+ p1_(p.p1_),
+ I0_(p.I0_),
+ I_(p.I_),
+ dA_(p.dA_),
+ reflectedId_(p.reflectedId_),
+ Imin_(p.Imin_)
+{}
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::DTRMParticle::move
+(
+ trackingData& td,
+ const scalar trackTime
+)
+{
+ td.switchProcessor = false;
+ td.keepParticle = true;
+
+ const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();
+
+ while (td.keepParticle && !td.switchProcessor)
+ {
+ point p0 = position();
+ label cell0 = cell();
+
+ scalar dt = trackToFace(p1_, td);
+
+ // Consider the cell between f0(start of tracking) and f1
+ label celli = cell();
+
+ const vector dsv = position() - p0;
+ const scalar ds = mag(dsv);
+
+ // Boltzman constant
+ const scalar sigma = physicoChemical::sigma.value();
+ if
+ (
+ (!td.relfectedCells()[celli] > 0 && reflectedId_ == 0)
+ || reflectedId_ > 0
+ )
+ {
+ scalar a = td.aInterp().interpolate(position(), cell0);
+ scalar e = td.eInterp().interpolate(position(), cell0);
+ scalar E = td.EInterp().interpolate(position(), cell0);
+ scalar T = td.TInterp().interpolate(position(), cell0);
+
+ const scalar I1 =
+ (
+ I_
+ + ds*(e*sigma*pow4(T)/mathematical::pi + E)
+ ) / (1 + ds*a);
+
+ td.Q(cell0) += (I_ - I1)*dA_;
+
+ I_ = I1;
+
+ if ((I_ <= 0.01*I0_) || (I_ < Imin_))
+ {
+ break;
+ }
+ }
+ else
+ {
+ scalar rho(0);
+ // Create a new reflected particle when the particles is not
+ // transmissive and larger than an absolute I
+ if (reflectedId_ == 0 && I_ > Imin_)
+ {
+ vector pDir = dsv/ds;
+
+ cellPointWeight cpw(mesh_, position(), celli, face());
+ //vector nHat = td.nHatCells()[celli];
+ vector nHat = td.nHatInterp().interpolate(cpw);
+
+ nHat /= mag(nHat);
+ scalar cosTheta(-pDir & nHat);
+ // Only new incoming rays
+ if (cosTheta > SMALL)
+ {
+ vector newDir = td.reflection().R(pDir, nHat);
+
+ //scalar theta = acos(-pDir & nHat);
+
+ // reflectivity
+ rho = min(max(td.reflection().rho(cosTheta), 0.0), 0.98);
+
+ scalar delaM = sqrt(mesh_.cellVolumes()[cell0]);
+
+ DTRMParticle* pPtr = new DTRMParticle
+ (
+ mesh_,
+ position() - pDir*0.1*delaM,
+ position() + newDir*mesh_.bounds().mag(),
+ I_*rho,
+ cell0,
+ dA_,
+ reflectedId_,
+ Imin_,
+ true
+ );
+ // Add to cloud
+ td.cloud().addParticle(pPtr);
+ }
+ }
+
+ reflectedId_++;
+
+ const point p0 = position();
+
+ // Try to locate this particle across the reflecting surface in
+ // a pure phase face
+ scalar dt = trackToFace(p1_, td);
+ const scalar ds = mag(position() - p0);
+
+ scalar a = td.aInterp().interpolate(position(), celli);
+ scalar e = td.eInterp().interpolate(position(), celli);
+ scalar E = td.EInterp().interpolate(position(), celli);
+ scalar T = td.TInterp().interpolate(position(), celli);
+
+ // Left intensity after reflection
+ const scalar Itran = I_*(1.0 - rho);
+ const scalar I1 =
+ (
+ Itran
+ + ds*(e*sigma*pow4(T)/mathematical::pi + E)
+ ) / (1 + ds*a);
+
+ td.Q(celli) += (Itran - I1)*dA_;
+
+ I_ = I1;
+
+ if (I_ <= 0.01*I0_ || I_ < Imin_)
+ {
+ break;
+ }
+ }
+
+ if (onBoundary() && td.keepParticle)
+ {
+ if (isA(pbMesh[patch(face())]))
+ {
+ td.switchProcessor = true;
+ }
+ }
+ }
+
+ return td.keepParticle;
+}
+
+
+bool Foam::DTRMParticle::hitPatch
+(
+ const polyPatch&,
+ particle::TrackingData>&,
+ const label,
+ const scalar,
+ const tetIndices&
+)
+{
+ return false;
+}
+
+
+void Foam::DTRMParticle::hitProcessorPatch
+(
+ const processorPolyPatch&,
+ particle::TrackingData>& td
+)
+{
+ td.switchProcessor = true;
+}
+
+
+void Foam::DTRMParticle::hitWallPatch
+(
+ const wallPolyPatch& wpp,
+ particle::TrackingData>& td,
+ const tetIndices& tetIs
+)
+{
+ td.keepParticle = false;
+}
+
+
+void Foam::DTRMParticle::hitPatch
+(
+ const polyPatch&,
+ particle::TrackingData>& td
+)
+{
+ td.keepParticle = false;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H
new file mode 100644
index 0000000000..e99925b2cf
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H
@@ -0,0 +1,347 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd
+ \\/ 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::DTRMParticle
+
+Description
+ Discrete Transfer Radiation Model (DTRM) particle
+
+SourceFiles
+ DTRMParticle.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DTRMParticle_H
+#define DTRMParticle_H
+
+#include "particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+#include "interpolationCell.H"
+#include "volFieldsFwd.H"
+#include "reflectionModel.H"
+#include "interpolationCellPoint.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class DTRMParticle;
+using namespace Foam::radiation;
+
+// Forward declaration of friend functions
+Ostream& operator<<(Ostream&, const DTRMParticle&);
+
+/*---------------------------------------------------------------------------*\
+ Class DTRMParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class DTRMParticle
+:
+ public particle
+{
+private:
+
+ // Private data
+
+ //- Size in bytes of the fields
+ static const std::size_t sizeofFields_;
+
+ //- Initial position
+ point p0_;
+
+ //- Target position
+ point p1_;
+
+ //- Initial radiation intensity [W/m2]
+ scalar I0_;
+
+ //- Radiation intensity [W/m2]
+ scalar I_;
+
+ //- Area of radiation
+ scalar dA_;
+
+ //- Reflected index
+ label reflectedId_;
+
+ //- Minimum radiation intensity to which the particle is tracked [W/m2]
+ scalar Imin_;
+
+
+public:
+
+ friend class Cloud;
+
+ //- Class used to pass tracking data to the trackToFace function
+ class trackingData
+ :
+ public particle::TrackingData>
+ {
+ // Interpolators for continuous phase fields
+
+ const interpolationCell& aInterp_;
+ const interpolationCell& eInterp_;
+ const interpolationCell& EInterp_;
+ const interpolationCell& TInterp_;
+
+ const interpolationCellPoint& nHatInterp_;
+
+ const labelField& relfectedCells_;
+ const reflectionModel& reflection_;
+
+ //- Heat source term
+ volScalarField& Q_;
+
+
+ public:
+
+ // Constructors
+
+ inline trackingData
+ (
+ Cloud& spc,
+ const interpolationCell& aInterp,
+ const interpolationCell& eInterp,
+ const interpolationCell& EInterp,
+ const interpolationCell& TInterp,
+ const interpolationCellPoint& nHatInterp,
+ const labelField&,
+ const reflectionModel&,
+ volScalarField& Q
+ );
+
+ // Public data
+
+
+
+
+ // Member functions
+
+ inline const interpolationCell& aInterp() const;
+ inline const interpolationCell& eInterp() const;
+ inline const interpolationCell& EInterp() const;
+ inline const interpolationCell& TInterp() const;
+ inline const interpolationCellPoint& nHatInterp() const;
+ inline const labelField& relfectedCells() const;
+ inline const reflectionModel& reflection() const;
+
+ inline scalar& Q(label celli);
+
+ };
+
+ // Static data members
+
+ //- String representation of properties
+ AddToPropertyList
+ (
+ particle,
+ " p0"
+ + " p1"
+ + " I0"
+ + " I"
+ + " dA";
+ );
+
+ //- String representation of property types
+ AddToPropertyTypes
+ (
+ particle,
+ "{point"
+ + " point"
+ + " scalar"
+ + " scalar}"
+ );
+
+
+ // Constructors
+
+
+ //- Construct from components, with searching for tetFace and
+ // tetPt unless disabled by doCellFacePt = false.
+ DTRMParticle
+ (
+ const polyMesh& mesh,
+ const vector& position,
+ const vector& targetPosition,
+ const scalar I,
+ const label cellI,
+ const scalar dA,
+ const label reflectedId,
+ const scalar Imin,
+ bool doCellFacePt = true
+ );
+
+ //- Construct from Istream
+ DTRMParticle
+ (
+ const polyMesh& mesh,
+ Istream& is,
+ bool readFields = true
+ );
+
+ //- Construct as copy
+ DTRMParticle(const DTRMParticle& p);
+
+ //- Construct and return a clone
+ virtual autoPtr clone() const
+ {
+ return autoPtr(new DTRMParticle(*this));
+ }
+
+
+ //- Factory class to read-construct particles used for
+ // parallel transfer
+ class iNew
+ {
+ const polyMesh& mesh_;
+
+ public:
+
+ iNew(const polyMesh& mesh)
+ :
+ mesh_(mesh)
+ {}
+
+ autoPtr operator()(Istream& is) const
+ {
+ return autoPtr
+ (
+ new DTRMParticle(mesh_, is, true)
+ );
+ }
+ };
+
+
+ // Access
+
+ //- Return const access to the initial position
+ inline const point& p0() const;
+
+ //- Return const access to the target position
+ inline const point& p1() const;
+
+ //- Return const access to the initial intensity
+ inline scalar I0() const;
+
+ //- Return const access to the current intensity
+ inline scalar I() const;
+
+ //- Return const access dA
+ inline scalar dA() const;
+
+
+ // Edit
+
+ //- Return access to the initial position
+ inline point& p0();
+
+ //- Return access to the target position
+ inline point& p1();
+
+ //- Return access to the initial intensity
+ inline scalar& I0();
+
+ //- Return access to the current intensity
+ inline scalar& I();
+
+ //- Return access to dA
+ inline scalar& dA();
+
+ //- Return access to reflectedId
+ inline label& reflectedId();
+
+ //- Return access to initial tet face
+ //inline label& tetFace0();
+
+ //- Return access to initial tet point
+ //inline label& tetPt0();
+
+ //- Return access to initial proc Id
+ //inline label& origProc0();
+
+
+ // Tracking
+
+ //- Move
+ bool move(trackingData&, const scalar);
+
+
+ // Member Operators
+
+ //- Overridable function to handle the particle hitting a patch
+ // Executed before other patch-hitting functions
+ bool hitPatch
+ (
+ const polyPatch&,
+ particle::TrackingData>& td,
+ const label patchi,
+ const scalar trackFraction,
+ const tetIndices& tetIs
+ );
+
+ //- Overridable function to handle the particle hitting a processorPatch
+ void hitProcessorPatch
+ (
+ const processorPolyPatch&,
+ particle::TrackingData>& td
+ );
+
+ //- Overridable function to handle the particle hitting a wallPatch
+ void hitWallPatch
+ (
+ const wallPolyPatch&,
+ particle::TrackingData>& td,
+ const tetIndices&
+ );
+
+ //- Overridable function to handle the particle hitting a polyPatch
+ void hitPatch
+ (
+ const polyPatch&,
+ particle::TrackingData>& td
+ );
+
+
+ // Ostream Operator
+
+ friend Ostream& operator<<(Ostream& os, const DTRMParticle& p);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DTRMParticleI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H
new file mode 100644
index 0000000000..5aea5a45a0
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H
@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 OpenCFD Ltd
+ \\/ 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 .
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+inline Foam::DTRMParticle::trackingData::trackingData
+(
+ Cloud& spc,
+ const interpolationCell& aInterp,
+ const interpolationCell& eInterp,
+ const interpolationCell& EInterp,
+ const interpolationCell& TInterp,
+ const interpolationCellPoint& nHatInterp,
+ const labelField& relfectedCell,
+ const reflectionModel& reflection,
+ volScalarField& Q
+)
+:
+ particle::TrackingData>(spc),
+ aInterp_(aInterp),
+ eInterp_(eInterp),
+ EInterp_(EInterp),
+ TInterp_(TInterp),
+ nHatInterp_(nHatInterp),
+ relfectedCells_(relfectedCell),
+ reflection_(reflection),
+ Q_(Q)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+inline const Foam::interpolationCell&
+Foam::DTRMParticle::trackingData::aInterp() const
+{
+ return aInterp_;
+}
+
+
+inline const Foam::interpolationCell&
+Foam::DTRMParticle::trackingData::eInterp() const
+{
+ return eInterp_;
+}
+
+
+inline const Foam::interpolationCell&
+Foam::DTRMParticle::trackingData::EInterp() const
+{
+ return EInterp_;
+}
+
+
+inline const Foam::interpolationCell&
+Foam::DTRMParticle::trackingData::TInterp() const
+{
+ return TInterp_;
+}
+
+inline const Foam::interpolationCellPoint&
+Foam::DTRMParticle::trackingData::nHatInterp() const
+{
+ return nHatInterp_;
+}
+
+inline const Foam::labelField&
+Foam::DTRMParticle::trackingData::relfectedCells() const
+{
+ return relfectedCells_;
+}
+
+
+inline const Foam::reflectionModel&
+Foam::DTRMParticle::trackingData::reflection() const
+{
+ return reflection_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::trackingData::Q(label celli)
+{
+ return Q_[celli];
+}
+
+
+inline const Foam::point& Foam::DTRMParticle::p0() const
+{
+ return p0_;
+}
+
+
+inline const Foam::point& Foam::DTRMParticle::p1() const
+{
+ return p1_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I0() const
+{
+ return I0_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I() const
+{
+ return I_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::dA() const
+{
+ return dA_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::dA()
+{
+ return dA_;
+}
+
+
+inline Foam::point& Foam::DTRMParticle::p0()
+{
+ return p0_;
+}
+
+
+inline Foam::label& Foam::DTRMParticle::reflectedId()
+{
+ return reflectedId_;
+}
+
+
+inline Foam::point& Foam::DTRMParticle::p1()
+{
+ return p1_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I0()
+{
+ return I0_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I()
+{
+ return I_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C
new file mode 100644
index 0000000000..6e0af2a56b
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd
+ \\/ 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 "DTRMParticle.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::string Foam::DTRMParticle::propertyList_ =
+ Foam::DTRMParticle::propertyList();
+
+const std::size_t Foam::DTRMParticle::sizeofFields_
+(
+ sizeof(DTRMParticle) - sizeof(particle)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::DTRMParticle::DTRMParticle
+(
+ const polyMesh& mesh,
+ Istream& is,
+ bool readFields
+)
+:
+ particle(mesh, is, readFields),
+ p0_(position_),
+ p1_(point::zero),
+ I0_(0),
+ I_(0),
+ dA_(0)
+{
+ if (readFields)
+ {
+ if (is.format() == IOstream::ASCII)
+ {
+ is >> p0_ >> p1_ >> I0_ >> I_ >> dA_;
+ }
+ else
+ {
+ is.read(reinterpret_cast(&p0_), sizeofFields_);
+ }
+ }
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const DTRMParticle& p)
+{
+ if (os.format() == IOstream::ASCII)
+ {
+ os << static_cast(p)
+ << token::SPACE << p.p0_
+ << token::SPACE << p.p1_
+ << token::SPACE << p.I0_
+ << token::SPACE << p.I_
+ << token::SPACE << p.dA_;
+ }
+ else
+ {
+ os << static_cast(p);
+ os.write
+ (
+ reinterpret_cast(&p.p0_),
+ DTRMParticle::sizeofFields_
+ );
+ }
+
+ // Check state of Ostream
+ os.check("Ostream& operator<<(Ostream&, const DTRMParticle&)");
+
+ return os;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/files b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/files
new file mode 100644
index 0000000000..04d0fb42fc
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/files
@@ -0,0 +1,13 @@
+laserDTRM.C
+
+DTRMParticle/DTRMParticle.C
+DTRMParticle/DTRMParticleIO.C
+localDensityAbsorptionEmission/localDensityAbsorptionEmission.C
+reflectionModel/reflectionModel/reflectionModel.C
+reflectionModel/reflectionModel/reflectionModelNew.C
+reflectionModel/noReflection/noReflection.C
+reflectionModel/FresnelLaser/FresnelLaser.C
+reflectionModel/Fresnel/Fresnel.C
+
+
+LIB = $(FOAM_LIBBIN)/liblaserDTRM
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/options b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/options
new file mode 100644
index 0000000000..1a31e07a96
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/options
@@ -0,0 +1,11 @@
+EXE_INC = \
+ -DFULLDEBUG -g -O0 \
+ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+ -lradiationModels \
+ -llagrangian \
+ -lfiniteVolume
diff --git a/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/laserDTRM.C b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/laserDTRM.C
new file mode 100644
index 0000000000..c89e0ad3e6
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/laserDTRM.C
@@ -0,0 +1,753 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 OpenCFD Ltd
+ \\/ 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 "laserDTRM.H"
+#include "fvmLaplacian.H"
+#include "fvmSup.H"
+#include "absorptionEmissionModel.H"
+#include "scatterModel.H"
+#include "constants.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "interpolationCell.H"
+#include "interpolationCellPoint.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ namespace radiation
+ {
+ defineTypeNameAndDebug(laserDTRM, 0);
+ addToRadiationRunTimeSelectionTables(laserDTRM);
+ }
+
+ defineTemplateTypeNameAndDebugWithName
+ (
+ Cloud,
+ "DTRMCloud",
+ 0
+ );
+
+ namespace radiation
+ {
+ template<>
+ const char* Foam::NamedEnum
+ <
+ Foam::radiation::laserDTRM::powerDistributionMode,
+ 3
+ >::names[] =
+ {
+ "Gaussian",
+ "manual",
+ "uniform"
+ };
+ }
+}
+
+const Foam::NamedEnum
+<
+ Foam::radiation::laserDTRM::powerDistributionMode,
+ 3
+> Foam::radiation::laserDTRM::powerDistypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::scalar Foam::radiation::laserDTRM::calculateIp(scalar r, scalar theta)
+{
+ const scalar t = mesh_.time().value();
+ switch(mode_)
+ {
+ case pdGaussian:
+ {
+ scalar I0 =
+ laserPower_->value(t)/(mathematical::twoPi*sqr(sigma_));
+
+ return(I0*exp(-sqr(r)/2.0/sqr(sigma_)));
+
+ break;
+ }
+ case pdManual:
+ {
+ return(laserPower_->value(t)*powerDistribution_()(theta, r));
+ break;
+ }
+ case pdUniform:
+ {
+ return
+ (
+ laserPower_->value(t)/(mathematical::pi*sqr(focalLaserRadius_))
+ );
+ }
+ default:
+ {
+ FatalErrorInFunction
+ << "Unhandled type " << powerDistypeNames_
+ << abort(FatalError);
+ return(0);
+ }
+ }
+}
+
+
+Foam::tmp Foam::radiation::laserDTRM::nHatfv
+(
+ const volScalarField& alpha1,
+ const volScalarField& alpha2
+) const
+{
+ const dimensionedScalar deltaN
+ (
+ "deltaN",
+ 1e-8/pow(average(mesh_.V()), 1.0/3.0)
+ );
+
+ const dimensionedScalar minAlpha
+ (
+ "minAlpha", dimless, 1e-3
+ );
+
+ volVectorField gradAlphaf
+ (
+ "gradAlphaf",
+ (alpha2 + minAlpha)*fvc::grad(alpha1)
+ - (alpha1 + minAlpha)*fvc::grad(alpha2)
+ );
+
+ // Face unit interface normal
+ return gradAlphaf/(mag(gradAlphaf)+ deltaN);
+}
+
+
+Foam::tmpFoam::radiation::laserDTRM::nearInterface
+(
+ const volScalarField& alpha1,
+ const volScalarField& alpha2
+) const
+{
+ return
+ pos(alpha1 - 0.1)*pos(0.9 - alpha1)
+ * pos(alpha2 - 0.1)*pos(0.9 - alpha2);
+}
+
+
+void Foam::radiation::laserDTRM::initialise()
+{
+ // Initialise the DTRM particles
+ DTRMCloud_.clear();
+
+ const scalar t = mesh_.time().value();
+ const vector lPosition = focalLaserPosition_->value(t);
+ vector lDir = laserDirection_->value(t);
+ lDir /= mag(lDir);
+
+ if (debug)
+ {
+ Info << "Laser position : " << lPosition << endl;
+ Info << "Laser direction : " << lDir << endl;
+ }
+
+ // Find a vector on the area plane. Normal to laser direction
+ vector rArea = vector::zero;
+ scalar magr = 0.0;
+ cachedRandom rnd(1234, -1);
+ while (magr < VSMALL)
+ {
+ vector v = rnd.sample01();
+ rArea = v - (v & lDir)*lDir;
+ magr = mag(rArea);
+ }
+ rArea /= mag(rArea);
+
+ scalar dr = focalLaserRadius_/ndr_;
+ scalar dTheta = mathematical::twoPi/ndTheta_;
+
+ nParticles_ = ndr_*ndTheta_;
+
+ switch(mode_)
+ {
+ case pdGaussian:
+ {
+ sigma_ = readScalar(lookup("sigma"));
+ break;
+ }
+ case pdManual:
+ {
+ powerDistribution_.reset
+ (
+ new interpolation2DTable(*this)
+ );
+
+ // Check dimensions ndr and ndTheta
+// if
+// (
+// (powerDistribution_->size() != ndTheta_)
+// || (powerDistribution_().first().second().size() != ndr_)
+// )
+// {
+// FatalErrorInFunction
+// << " The table dimensions should correspond with ndTheta "
+// << " and ndr "
+// << exit(FatalError);
+// }
+
+ break;
+ }
+ case pdUniform:
+ {
+ break;
+ }
+ }
+
+ // Target position
+ point p1 = vector::zero;
+
+ // Seed DTRM particles
+ // TODO: currently only applicable to 3-D cases
+ point p0 = lPosition;
+ scalar power(0.0);
+ scalar area(0.0);
+ p1 = p0;
+ if (mesh_.nGeometricD() == 3)
+ {
+ //scalar r0 = dr/2.0;
+ //scalar r1Max0 = drMax/2.0;
+
+ for (label ri = 0; ri < ndr_; ri++)
+ {
+ scalar r1 = SMALL + dr*ri;
+
+ scalar r2 = r1 + dr;
+
+ scalar rP = ((r1 + r2)/2);
+
+ // local radius on disk
+ vector localR = ((r1 + r2)/2)*rArea;
+
+ // local final radius on disk
+ vector finalR = rP*rArea;
+
+ scalar theta0 = 0.0;//dTheta/2.0;
+ for (label thetai = 0; thetai < ndTheta_; thetai++)
+ {
+ scalar theta1 = theta0 + SMALL + dTheta*thetai;
+
+ scalar theta2 = theta1 + dTheta;
+
+ scalar thetaP = (theta1 + theta2)/2.0;
+
+ quaternion Q(lDir, thetaP);
+
+ // Initial position on disk
+ vector initialPos = (Q.R() & localR);
+
+ // Final position on disk
+ vector finalPos = (Q.R() & finalR);
+
+ // Initial position
+ p0 = lPosition + initialPos;
+
+ // calculate target point using new deviation rl
+ p1 = lPosition + finalPos + (0.5*maxTrackLength_*lDir);
+
+ //scalar p = magSqr(p0 - lPosition);
+
+ scalar Ip = calculateIp(rP, thetaP);
+
+ scalar dAi = (sqr(r2) - sqr(r1))*(theta2 - theta1)/2.0;
+
+ power += Ip*dAi;
+ area += dAi;
+
+ label cellI = mesh_.findCell(p0);
+
+ if (cellI != -1)
+ {
+ // Create a new particle
+ DTRMParticle* pPtr = new DTRMParticle
+ (mesh_, p0, p1, Ip, cellI, dAi, 0 , 0.01*Ip, true);
+
+ // Add to cloud
+ DTRMCloud_.addParticle(pPtr);
+ }
+
+ if (returnReduce(cellI, maxOp