ENH: Add local rho into reactingOneDim.C

This commit is contained in:
sergio
2012-08-07 16:33:33 +01:00
parent 043576f61f
commit 4604b692f8
3 changed files with 64 additions and 37 deletions

View File

@ -232,12 +232,22 @@ void reactingOneDim::solveContinuity()
Info<< "reactingOneDim::solveContinuity()" << endl;
}
solve
(
fvm::ddt(rho_)
==
- solidChemistry_->RRg()
);
if (moveMesh_)
{
const scalarField mass0 = rho_*regionMesh().V();
fvScalarMatrix rhoEqn
(
fvm::ddt(rho_)
==
- solidChemistry_->RRg()
);
rhoEqn.solve();
updateMesh(mass0);
}
}
@ -261,14 +271,15 @@ void reactingOneDim::solveSpeciesMass()
solidChemistry_->RRs(i)
);
if (moveMesh_)
if (regionMesh().moving())
{
surfaceScalarField phiRhoMesh
surfaceScalarField phiYiRhoMesh
(
fvc::interpolate(Yi*rho_)*regionMesh().phi()
);
YiEqn -= fvc::div(phiRhoMesh);
YiEqn += fvc::div(phiYiRhoMesh);
}
YiEqn.solve(regionMesh().solver("Yi"));
@ -303,14 +314,14 @@ void reactingOneDim::solveEnergy()
+ fvc::div(phiGas)
);
if (moveMesh_)
if (regionMesh().moving())
{
surfaceScalarField phiMesh
surfaceScalarField phihMesh
(
fvc::interpolate(rho_*h_)*regionMesh().phi()
);
hEqn -= fvc::div(phiMesh);
hEqn += fvc::div(phihMesh);
}
hEqn.relax();
@ -349,7 +360,18 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
pyrolysisModel(modelType, mesh),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solid()),
rho_(solidThermo_.rho()),
rho_
(
IOobject
(
"rho",
regionMesh().time().timeName(),
regionMesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
solidThermo_.rho()
),
Ys_(solidThermo_.composition().Y()),
h_(solidThermo_.he()),
primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
@ -449,7 +471,18 @@ reactingOneDim::reactingOneDim
pyrolysisModel(modelType, mesh, dict),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solid()),
rho_(solidThermo_.rho()),
rho_
(
IOobject
(
"rho",
regionMesh().time().timeName(),
regionMesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
solidThermo_.rho()
),
Ys_(solidThermo_.composition().Y()),
h_(solidThermo_.he()),
primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
@ -681,17 +714,15 @@ void reactingOneDim::evolveRegion()
{
Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
const scalarField mass0 = rho_*regionMesh().V();
solidChemistry_->solve
(
time().value() - time().deltaTValue(),
time().deltaTValue()
);
solveContinuity();
calculateMassTransfer();
updateMesh(mass0);
solveContinuity();
chemistrySh_ = solidChemistry_->Sh()();
@ -704,9 +735,9 @@ void reactingOneDim::evolveRegion()
solveEnergy();
}
calculateMassTransfer();
solidThermo_.correct();
rho_ = solidThermo_.rho();
}

View File

@ -84,14 +84,8 @@ protected:
// Reference to solid thermo properties
// //- Absorption coefficient [1/m]
// const volScalarField& kappaRad_;
//
// //- Thermal conductivity [W/m/K]
// const volScalarField& kappa_;
//- Density [kg/m3]
volScalarField& rho_;
volScalarField rho_;
//- List of solid components
PtrList<volScalarField>& Ys_;
@ -221,8 +215,8 @@ public:
//- Fields
//- Return density [kg/m3]
virtual const volScalarField& rho() const;
//- Return const density [Kg/m3]
const volScalarField& rho() const;
//- Return const temperature [K]
virtual const volScalarField& T() const;

View File

@ -117,6 +117,8 @@ void Foam::regionModels::regionModel1D::initialise()
}
boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI);
boundaryFaceFaces_.setSize(localPyrolysisFaceI);
boundaryFaceCells_.setSize(localPyrolysisFaceI);
surfaceScalarField& nMagSf = nMagSfPtr_();
@ -192,6 +194,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
label totalFaceId = 0;
forAll(intCoupledPatchIDs_, localPatchI)
{
label patchI = intCoupledPatchIDs_[localPatchI];
@ -200,8 +203,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
forAll(pp, patchFaceI)
{
const labelList& faces = boundaryFaceFaces_[patchFaceI];
const labelList& cells = boundaryFaceCells_[patchFaceI];
const labelList& faces = boundaryFaceFaces_[totalFaceId];
const labelList& cells = boundaryFaceCells_[totalFaceId];
const vector n = pp.faceNormals()[patchFaceI];
const vector sf = pp.faceAreas()[patchFaceI];
@ -231,7 +235,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
if
(
((nbrCf - (oldPoints[pointI] + newDelta)) & n)
mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
> minDelta
)
{
@ -242,19 +246,17 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
}
nbrCf = oldCf[i + 1] + localDelta;
}
// Modify boundary
const label bFaceI = boundaryFaceOppositeFace_[patchFaceI];
const label bFaceI = boundaryFaceOppositeFace_[totalFaceId];
const face f = regionMesh().faces()[bFaceI];
const label cellI = cells[cells.size() - 1];
newDelta += (deltaV[cellI]/mag(sf))*n;
forAll(f, pti)
{
const label pointI = f[pti];
if
(
((nbrCf - (oldPoints[pointI] + newDelta)) & n)
mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
> minDelta
)
{
@ -262,9 +264,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
cellMoveMap[cellI] = 1;
}
}
totalFaceId ++;
}
}
// Move points
regionMesh().movePoints(newPoints);