Commit dccad820 authored by sergio's avatar sergio
Browse files

ENH: Modification of solid thermo, thermo baffles, pyrolysis and

tutorial, solvers solve for h in the solid
parent e36f6bf7
......@@ -43,7 +43,7 @@ Description
#include "regionProperties.H"
#include "compressibleCourantNo.H"
#include "solidRegionDiffNo.H"
#include "basicSolidThermo.H"
#include "solidThermo.H"
#include "radiationModel.H"
#include "porousZones.H"
#include "IObasicSourceList.H"
......
......@@ -34,7 +34,7 @@ Description
#include "turbulenceModel.H"
#include "fixedGradientFvPatchFields.H"
#include "regionProperties.H"
#include "basicSolidThermo.H"
#include "solidThermo.H"
#include "radiationModel.H"
#include "porousZones.H"
#include "IObasicSourceList.H"
......
// Initialise solid field pointer lists
PtrList<basicSolidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
PtrList<volScalarField> betavSolid(porousSolidRegions.size());
......@@ -14,7 +14,7 @@
porousSolidThermos.set
(
i,
basicSolidThermo::New(porousSolidRegions[i])
solidThermo::New(porousSolidRegions[i])
);
Info<< " Adding sources\n" << endl;
solidHeatSources.set
......
const fvMesh& mesh = porousSolidRegions[i];
basicSolidThermo& thermo = porousSolidThermos[i];
solidThermo& thermo = porousSolidThermos[i];
const volScalarField& betav = betavSolid[i];
tmp<volScalarField> trho = thermo.rho();
......@@ -15,9 +15,9 @@
const volScalarField& kappa = tkappa();
//const volSymmTensorField& K = tK();
tmp<volScalarField> trhoCp = cp*rho;
const volScalarField& rhoCp = trhoCp();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& T = thermo.T();
volScalarField& h = thermo.he();
IObasicSourceList& sources = solidHeatSources[i];
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> TEqn
tmp<fvScalarMatrix> hEqn
(
- fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
+ sources(rhoCp, T)
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
TEqn().relax();
TEqn().solve();
hEqn().relax();
hEqn().solve();
}
Info<< "Min/max T:" << min(T).value() << ' ' << max(T).value() << endl;
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
// Initialise solid field pointer lists
PtrList<basicSolidThermo> thermos(solidRegions.size());
PtrList<solidThermo> thermos(solidRegions.size());
// Populate solid field pointer lists
forAll(solidRegions, i)
......@@ -11,6 +11,6 @@
thermos.set
(
i,
basicSolidThermo::New(solidRegions[i])
solidThermo::New(solidRegions[i])
);
}
fvMesh& mesh = solidRegions[i];
basicSolidThermo& thermo = thermos[i];
solidThermo& thermo = thermos[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
......@@ -11,4 +11,7 @@
//tmp<volSymmTensorField> tkappa = thermo.directionalkappa();
const volScalarField& kappa = tkappa();
volScalarField& T = thermo.T();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix tEqn
fvScalarMatrix hEqn
(
-fvm::laplacian(kappa, T)
-fvm::laplacian(alpha, h)
);
tEqn.relax();
tEqn.solve();
hEqn.relax();
hEqn.solve();
}
Info<< "Min/max T:" << min(T).value() << ' '
<< max(T).value() << endl;
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
// Initialise solid field pointer lists
PtrList<basicSolidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
PtrList<volScalarField> betavSolid(porousSolidRegions.size());
......@@ -14,7 +14,7 @@
porousSolidThermos.set
(
i,
basicSolidThermo::New(porousSolidRegions[i])
solidThermo::New(porousSolidRegions[i])
);
Info<< " Adding sources\n" << endl;
solidHeatSources.set
......
fvMesh& mesh = porousSolidRegions[i];
basicSolidThermo& thermo = porousSolidThermos[i];
solidThermo& thermo = porousSolidThermos[i];
const volScalarField& betav = betavSolid[i];
tmp<volScalarField> trho = thermo.rho();
......@@ -15,9 +15,12 @@
const volScalarField& kappa = tkappa();
//const volSymmTensorField& K = tK();
tmp<volScalarField> trhoCp = cp*rho;
const volScalarField& rhoCp = trhoCp();
//tmp<volScalarField> trhoCp = cp*rho;
//const volScalarField& rhoCp = trhoCp();
volScalarField& T = thermo.T();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();
IObasicSourceList& sources = solidHeatSources[i];
......@@ -6,22 +6,22 @@ if (finalIter)
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> TEqn
tmp<fvScalarMatrix> hEqn
(
fvm::ddt(betav*rho*cp, T)
- fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
+ sources(rhoCp, T)
fvm::ddt(betav*rho, h)
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
TEqn().relax();
TEqn().solve(mesh.solver(T.select(finalIter)));
hEqn().relax();
hEqn().solve(mesh.solver(h.select(finalIter)));
}
Info<< "Min/max T:" << min(T).value() << ' ' << max(T).value() << endl;
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
if (finalIter)
{
mesh.data::remove("finalIteration");
......
// Initialise solid field pointer lists
PtrList<basicSolidThermo> thermos(solidRegions.size());
PtrList<solidThermo> thermos(solidRegions.size());
// Populate solid field pointer lists
forAll(solidRegions, i)
......@@ -8,5 +8,5 @@
<< solidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
thermos.set(i, basicSolidThermo::New(solidRegions[i]));
thermos.set(i, solidThermo::New(solidRegions[i]));
}
fvMesh& mesh = solidRegions[i];
basicSolidThermo& thermo = thermos[i];
solidThermo& thermo = thermos[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
......@@ -7,9 +7,12 @@
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
tmp<volScalarField> tkappa = thermo.kappa();
const volScalarField& kappa = tkappa();
//tmp<volSymmTensorField> tkappa = thermo.directionalKappa();
//const volSymmTensorField& kappa = tkappa();
volScalarField& T = thermo.T();
volScalarField& h = thermo.he();
......@@ -6,20 +6,20 @@ if (finalIter)
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> TEqn
tmp<fvScalarMatrix> hEqn
(
fvm::ddt(rho*cp, T)
- fvm::laplacian(kappa, T)
fvm::ddt(rho, h)
- fvm::laplacian(alpha, h)
);
TEqn().relax();
TEqn().solve(mesh.solver(T.select(finalIter)));
hEqn().relax();
hEqn().solve(mesh.solver(h.select(finalIter)));
}
Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
if (finalIter)
{
mesh.data::remove("finalIteration");
......
......@@ -51,7 +51,7 @@ void noPyrolysis::constructThermoChemistry()
solidChemistryModel::New(regionMesh()).ptr()
);
solidThermo_.reset(&solidChemistry_->solidThermo());
solidThermo_.reset(&solidChemistry_->solid());
}
bool noPyrolysis::read()
......@@ -151,13 +151,13 @@ const tmp<volScalarField> noPyrolysis::Cp() const
}
const volScalarField& noPyrolysis::kappaRad() const
tmp<volScalarField> noPyrolysis::kappaRad() const
{
return (solidThermo_->kappaRad());
}
const volScalarField& noPyrolysis::kappa() const
tmp<volScalarField> noPyrolysis::kappa() const
{
return (solidThermo_->kappa());
}
......
......@@ -37,6 +37,7 @@ SourceFiles
#include "pyrolysisModel.H"
#include "volFieldsFwd.H"
#include "solidChemistryModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -83,7 +84,7 @@ protected:
autoPtr<solidChemistryModel> solidChemistry_;
//- Reference to solid thermo
autoPtr<basicSolidThermo> solidThermo_;
autoPtr<solidReactionThermo> solidThermo_;
public:
......@@ -124,10 +125,10 @@ public:
virtual const tmp<volScalarField> Cp() const;
//- Return the region absorptivity [1/m]
virtual const volScalarField& kappaRad() const;
virtual tmp<volScalarField> kappaRad() const;
//- Return the region thermal conductivity [W/m/k]
virtual const volScalarField& kappa() const;
virtual tmp<volScalarField> kappa() const;
//- Return the total gas mass flux to primary region [kg/m2/s]
virtual const surfaceScalarField& phiGas() const;
......
......@@ -37,8 +37,6 @@ SourceFiles
#include "runTimeSelectionTables.H"
#include "volFieldsFwd.H"
#include "solidChemistryModel.H"
#include "basicSolidThermo.H"
#include "regionModel1D.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -218,10 +216,10 @@ public:
virtual const tmp<volScalarField> Cp() const = 0;
//- Return the region absorptivity [1/m]
virtual const volScalarField& kappaRad() const = 0;
virtual tmp<volScalarField> kappaRad() const = 0;
//- Return the region thermal conductivity [W/m/k]
virtual const volScalarField& kappa() const = 0;
virtual tmp<volScalarField> kappa() const = 0;
//- Return the total gas mass flux to primary region [kg/m2/s]
virtual const surfaceScalarField& phiGas() const = 0;
......
......@@ -54,7 +54,7 @@ void reactingOneDim::readReactingOneDimControls()
{
const dictionary& solution = this->solution().subDict("SIMPLE");
solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
time_.controlDict().lookup("maxDi") >> maxDiff_;
time().controlDict().lookup("maxDi") >> maxDiff_;
coeffs().lookup("radFluxName") >> primaryRadFluxName_;
coeffs().lookup("minimumDelta") >> minimumDelta_;
......@@ -107,6 +107,8 @@ void reactingOneDim::updateQr()
Qrp = max(Qrp, scalar(0.0));
}
const volScalarField kappaRad_(kappaRad());
// Propagate Qr through 1-D regions
forAll(intCoupledPatchIDs_, i)
{
......@@ -147,18 +149,17 @@ void reactingOneDim::updatePhiGas()
forAll(gasTable, gasI)
{
tmp<volScalarField> tHsiGas =
solidChemistry_->gasHs(p, T_, gasI);
solidChemistry_->gasHs(solidThermo_.p(), solidThermo_.T(), gasI);
tmp<volScalarField> tRRiGas = solidChemistry_->RRg(gasI);
const volScalarField& HsiGas = tHsiGas();
const volScalarField& RRiGas = tRRiGas();
const surfaceScalarField HsiGasf(fvc::interpolate(HsiGas));
const surfaceScalarField RRiGasf(fvc::interpolate(RRiGas));
forAll(intCoupledPatchIDs_, i)
{
const label patchI = intCoupledPatchIDs_[i];
const scalarField& phiGasp = phiHsGas_.boundaryField()[patchI];
forAll(phiGasp, faceI)
......@@ -286,16 +287,16 @@ void reactingOneDim::solveEnergy()
Info<< "reactingOneDim::solveEnergy()" << endl;
}
const volScalarField rhoCp(rho_*solidThermo_.Cp());
tmp<volScalarField> alpha(solidThermo_.alpha());
const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
fvScalarMatrix TEqn
fvScalarMatrix hEqn
(
fvm::ddt(rhoCp, T_)
- fvm::laplacian(kappa_, T_)
fvm::ddt(rho_, h_)
- fvm::laplacian(alpha, h_)
==
chemistrySh_
+ fvc::div(phiQr)
......@@ -306,17 +307,17 @@ void reactingOneDim::solveEnergy()
{
surfaceScalarField phiMesh
(
fvc::interpolate(rhoCp*T_)*regionMesh().phi()
fvc::interpolate(rho_*h_)*regionMesh().phi()
);
TEqn -= fvc::div(phiMesh);
hEqn -= fvc::div(phiMesh);
}
TEqn.relax();
TEqn.solve();
hEqn.relax();
hEqn.solve();
Info<< "pyrolysis min/max(T) = " << min(T_).value() << ", "
<< max(T_).value() << endl;
Info<< "pyrolysis min/max(T) = " << min(solidThermo_.T()) << ", "
<< max(solidThermo_.T()) << endl;
}
......@@ -347,12 +348,10 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
:
pyrolysisModel(modelType, mesh),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solidThermo()),
kappaRad_(solidThermo_.kappaRad()),
kappa_(solidThermo_.kappa()),
rho_(solidThermo_.rho()),
solidThermo_(solidChemistry_->solid()),
rho_(solidThermo_.rhos()),
Ys_(solidThermo_.composition().Y()),
T_(solidThermo_.T()),
h_(solidThermo_.he()),
primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
nNonOrthCorr_(-1),
maxDiff_(10),
......@@ -449,12 +448,10 @@ reactingOneDim::reactingOneDim
:
pyrolysisModel(modelType, mesh, dict),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solidThermo()),
kappaRad_(solidThermo_.kappaRad()),
kappa_(solidThermo_.kappa()),
rho_(solidThermo_.rho()),
solidThermo_(solidChemistry_->solid()),
rho_(solidThermo_.rhos()),
Ys_(solidThermo_.composition().Y()),
T_(solidThermo_.T()),
h_(solidThermo_.he()),
primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
nNonOrthCorr_(-1),
maxDiff_(10),
......@@ -585,11 +582,11 @@ scalar reactingOneDim::solidRegionDiffNo() const
surfaceScalarField KrhoCpbyDelta
(
regionMesh().surfaceInterpolation::deltaCoeffs()
* fvc::interpolate(kappa_)
* fvc::interpolate(kappa())
/ fvc::interpolate(Cp()*rho_)
);
DiNum = max(KrhoCpbyDelta.internalField())*time_.deltaTValue();
DiNum = max(KrhoCpbyDelta.internalField())*time().deltaTValue();
}
return DiNum;
......@@ -610,7 +607,7 @@ const volScalarField& reactingOneDim::rho() const
const volScalarField& reactingOneDim::T() const
{
return T_;
return solidThermo_.T();
}
......@@ -620,15 +617,15 @@ const tmp<volScalarField> reactingOneDim::Cp() const
}
const volScalarField& reactingOneDim::kappaRad() const
tmp<volScalarField> reactingOneDim::kappaRad() const
{
return kappaRad_;
return solidThermo_.kappaRad();
}
const volScalarField& reactingOneDim::kappa() const
tmp<volScalarField> reactingOneDim::kappa() const
{
return kappa_;
return solidThermo_.kappa();
}
......@@ -643,7 +640,7 @@ void reactingOneDim::preEvolveRegion()
pyrolysisModel::preEvolveRegion();
// Initialise all cells as able to react
forAll(T_, cellI)
forAll(h_, cellI)
{
solidChemistry_->setCellReacting(cellI, true);
}
......
......@@ -36,6 +36,7 @@ SourceFiles
#define reactingOneDim_H
#include "pyrolysisModel.H"
#include "solidChemistryModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -78,16 +79,16 @@ protected:
autoPtr<solidChemistryModel> solidChemistry_;
//- Reference to solid thermo
basicSolidThermo& solidThermo_;
solidReactionThermo& solidThermo_;
// Reference to solid thermo properties
//- Absorption coefficient [1/m]
const volScalarField& kappaRad_;
//- Thermal conductivity [W/m/K]
const volScalarField& kappa_;
// //- Absorption coefficient [1/m]
// const volScalarField& kappaRad_;
//
// //- Thermal conductivity [W/m/K]
// const volScalarField& kappa_;
//- Density [kg/m3]
volScalarField& rho_;
......@@ -96,7 +97,7 @@ protected:
PtrList<volScalarField>& Ys_;
// Non-const access to temperature
volScalarField& T_;
volScalarField& h_;
//- Name of the radiative flux in the primary region
......@@ -230,10 +231,10 @@ public:
virtual const tmp<volScalarField> Cp() const;
//- Return the region absorptivity [1/m]
virtual const volScalarField& kappaRad() const;
virtual tmp<volScalarField> kappaRad() const;
//- Return the region thermal conductivity [W/m/k]
virtual const volScalarField& kappa() const;
virtual tmp<volScalarField> kappa() const;
//- Return the total gas mass flux to primary region [kg/m2/s]
virtual const surfaceScalarField& phiGas() const;
......
......@@ -222,9 +222,8 @@ void temperatureThermoBaffleFvPatchScalarField::write(Ostream& os) const
os.writeKeyword("thermoType") << solidThermoType_
<< token::END_STATEMENT << nl;
os