Commit 677564e9 authored by Andrew Heather's avatar Andrew Heather
Browse files

Merge branch 'integration-foundation' into 'develop'

Integration openfoam.org

See merge request !144
parents 15694bdc 884f2f86
const volScalarField& psi = thermo.psi();
const volScalarField& T = thermo.T();
filmModelType& surfaceFilm = tsurfaceFilm();
regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
const label inertIndex(composition.species()[inertSpecie]);
Info<< "\nConstructing surface film model" << endl;
typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
(
regionModels::surfaceFilmModel::New(mesh, g)
);
......@@ -9,7 +9,7 @@ IOobject phiBHeader
surfaceScalarField* phiBPtr = nullptr;
if (phiBHeader.typeHeaderOk<surfaceScalarField>(true))
if (phiBHeader.typeHeaderOk<surfaceScalarField>(true))
{
Info<< "Reading face flux ";
......
......@@ -11,7 +11,7 @@ IOobject turbulencePropertiesHeader
false
);
if (turbulencePropertiesHeader.typeHeaderOk<IOdictionary>(false))
if (turbulencePropertiesHeader.typeHeaderOk<IOdictionary>(true))
{
autoPtr<compressible::turbulenceModel> turbulence
(
......
......@@ -77,6 +77,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
// Store the particle positions
kinematicCloud.storeGlobalPositions();
mesh.update();
// Calculate absolute flux from the mapped surface velocity
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -68,6 +68,8 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
kinematicCloud.storeGlobalPositions();
mesh.update();
U.correctBoundaryConditions();
......
reactingParcelFilmFoam.C
EXE = $(FOAM_APPBIN)/reactingParcelFilmFoam
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
parcels.SU(U)
+ fvOptions(rho, U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
Info<< "\nConstructing reacting cloud" << endl;
basicReactingCloud parcels
(
"reactingCloud1",
rho,
U,
g,
slgThermo
);
Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::psiCombustionModel> combustion
(
combustionModels::psiCombustionModel::New(mesh)
);
psiReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo);
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
const word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.species().found(inertSpecie))
{
FatalIOErrorIn(args.executable().c_str(), thermo)
<< "Inert specie " << inertSpecie << " not found in available species "
<< composition.species()
<< exit(FatalIOError);
}
Info<< "Creating field rho\n" << endl;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
// Set the turbulence into the combustion model
combustion->setTurbulence(turbulence());
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
mesh.setFluxRequired(p_rgh.name());
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
IOdictionary additionalControlsDict
(
IOobject
(
"additionalControls",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Switch solvePrimaryRegion
(
additionalControlsDict.lookup("solvePrimaryRegion")
);
volScalarField Qdot
(
IOobject
(
"Qdot",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
);
#include "createDpdt.H"
#include "createK.H"
#include "createMRF.H"
#include "createClouds.H"
#include "createRadiationModel.H"
#include "createSurfaceFilmModel.H"
Info<< "\nConstructing surface film model" << endl;
typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- fvm::laplacian(rhorAUf, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
+ fvOptions(psi, p_rgh, rho.name())
);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
p = p_rgh + rho*gh;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
......@@ -19,6 +19,7 @@
==
rho*(U&g)
+ parcels.Sh(he)
+ surfaceFilm.Sh()
+ radiation->Sh(thermo, he)
+ Qdot
+ fvOptions(rho, he)
......@@ -35,6 +36,6 @@
thermo.correct();
radiation->correct();
Info<< "T gas min/max " << min(T).value() << ", "
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}
EXE_INC = \
-I. \
-I../reactingParcelFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/sampling/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 \
......@@ -17,33 +16,33 @@ EXE_INC = \
-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)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(FOAM_SOLVERS)/combustion/reactingFoam
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lsampling \
-lmeshTools \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-llagrangian \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lspecie \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lthermophysicalProperties \
-lreactionThermophysicalModels \
-lSLGThermo \
-lchemistryModel \
-lradiationModels \
-lODE \
-lregionModels \
-lradiationModels \
-lsurfaceFilmModels \
-lcombustionModels \
-lfvOptions \
-lsampling
-lsurfaceFilmDerivedFvPatchFields \
-llagrangian \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lODE \
-lcombustionModels
......@@ -6,8 +6,7 @@
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
rho()*g
+ parcels.SU(U)
parcels.SU(U)
+ fvOptions(rho, U)
);
......@@ -17,7 +16,18 @@
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvOptions.correct(U);
K = 0.5*magSqr(U);
......
......@@ -9,6 +9,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
)
);
{
combustion->correct();
Qdot = combustion->Qdot();
......@@ -24,11 +25,12 @@ tmp<fv::convectionScheme<scalar>> mvConvection
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
- fvm::laplacian(turbulence->alphaEff(), Yi)
==
parcels.SYi(i, Yi)
+ combustion->R(Yi)
+ fvOptions(rho, Yi)
+ combustion->R(Yi)
+ surfaceFilm.Srho(i)
);
YEqn.relax();
......
const label inertIndex(composition.species()[inertSpecie]);
const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi();
const label inertIndex(composition.species()[inertSpecie]);
regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
#include "createRDeltaT.H"
#include "readGravitationalAcceleration.H"
Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::rhoCombustionModel> combustion
......@@ -26,8 +24,7 @@ if (!composition.species().found(inertSpecie))
<< exit(FatalIOError);
}
volScalarField& p = thermo.p();
Info<< "Creating field rho\n" << endl;
volScalarField rho
(
IOobject
......@@ -41,6 +38,8 @@ volScalarField rho
thermo.rho()
);
volScalarField& p = thermo.p();
Info<< "\nReading field U\n" << endl;
volVectorField U
(
......@@ -57,30 +56,6 @@ volVectorField U
#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<compressible::turbulenceModel> turbulence
(
......@@ -96,6 +71,31 @@ autoPtr<compressible::turbulenceModel> turbulence
// Set the turbulence into the combustion model
combustion->setTurbulence(turbulence());
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
pressureControl pressureControl(p, rho, pimple.dict(), false);
mesh.setFluxRequired(p_rgh.name());
Info<< "Creating multi-variate interpolation scheme\n" << endl;
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
......@@ -105,6 +105,11 @@ forAll(Y, i)
}
fields.add(thermo.he());
Switch solvePrimaryRegion
(
pimple.dict().lookupOrDefault<Switch>("solvePrimaryRegion", true)
);
volScalarField Qdot
(
IOobject
......@@ -126,3 +131,4 @@ volScalarField Qdot
#include "createMRF.H"
#include "createRadiationModel.H"
#include "createClouds.H"
#include "createSurfaceFilmModel.H"
Info<< "\nConstructing surface film model" << endl;
autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
(
regionModels::surfaceFilmModel::New(mesh, g)
);
rho = thermo.rho();
if (!pimple.SIMPLErho())
{
rho = thermo.rho();
}
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
......@@ -7,6 +10,9 @@ const volScalarField psip0(psi*p);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
......@@ -14,58 +20,68 @@ surfaceScalarField phiHbyA
fvc::flux(rho*HbyA)