From 85e5abd3440c55831633d27b3707bc82061a43cc Mon Sep 17 00:00:00 2001 From: andy <a.heather@opencfd.co.uk> Date: Tue, 26 May 2009 12:01:39 +0100 Subject: [PATCH] updates + added porous zones --- .../trackedReactingParcelFoam/UEqn.H | 31 +++++++ .../trackedReactingParcelFoam/createFields.H | 80 ++++++------------- .../trackedReactingParcelFoam/pEqn.H | 54 ++++++++++--- .../trackedReactingParcelFoam.C | 56 ++++++------- 4 files changed, 129 insertions(+), 92 deletions(-) diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H index b34cd840019..c82307305b8 100644 --- a/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H +++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H @@ -10,6 +10,37 @@ UEqn.relax(); + tmp<volScalarField> trAU; + tmp<volTensorField> trTU; + + if (pressureImplicitPorosity) + { + tmp<volTensorField> tTU = tensor(I)*UEqn.A(); + pZones.addResistance(UEqn, tTU()); + trTU = inv(tTU()); + trTU().rename("rAU"); + + volVectorField gradp = fvc::grad(p); + + for (int UCorr=0; UCorr<nUCorr; UCorr++) + { + U = trTU() & (UEqn.H() - gradp); + } + U.correctBoundaryConditions(); + } + else + { + pZones.addResistance(UEqn); + + solve + ( + UEqn == -fvc::grad(p) + ); + + trAU = 1.0/UEqn.A(); + trAU().rename("rAU"); + } + if (momentumPredictor) { solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H index 079a4f2e1d1..b9cea162f4a 100644 --- a/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H +++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H @@ -1,4 +1,4 @@ - Info<< "Reading thermophysical properties\n" << endl; + Info<< "Reading thermophysical properties" << nl << endl; autoPtr<hCombustionThermo> thermo ( @@ -28,52 +28,7 @@ thermo->rho() ); -// lagrangian coal density field -/* volScalarField rholc - ( - IOobject - ( - "rholc", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0) - ); - -// lagrangian limestone density field - volScalarField rhols - ( - IOobject - ( - "rhols", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0) - ); - -// lagrangian total density field - volScalarField rhol - ( - IOobject - ( - "rhol", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0) - );*/ - - Info<< "\nReading field U\n" << endl; + Info<< "Reading field U" << nl << endl; volVectorField U ( IOobject @@ -103,7 +58,7 @@ dimensionedScalar("zero", dimless, 0.0) ); - Info<< "Creating turbulence model\n" << endl; + Info<< "Creating turbulence model" << nl << endl; autoPtr<compressible::turbulenceModel> turbulence ( compressible::turbulenceModel::New @@ -115,11 +70,11 @@ ) ); - Info<< "Creating field DpDt\n" << endl; + Info<< "Creating field DpDt" << nl << endl; volScalarField DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); - Info << "Constructing chemical mechanism" << endl; + Info << "Constructing chemical mechanism" << nl << endl; chemistryModel chemistry ( thermo(), @@ -134,8 +89,23 @@ } fields.add(h); - Info<< "Creating radiation model\n" << endl; - autoPtr<radiation::radiationModel> radiation - ( - radiation::radiationModel::New(T) - ); + Info<< "Creating porous zones" << nl << endl; + porousZones pZones(mesh); + Switch pressureImplicitPorosity(false); + + label nUCorr = 0; + if (pZones.size()) + { + // nUCorrectors for pressureImplicitPorosity + if (mesh.solutionDict().subDict("PISO").found("nUCorrectors")) + { + mesh.solutionDict().subDict("PISO").lookup("nUCorrectors") + >> nUCorr; + } + + if (nUCorr > 0) + { + pressureImplicitPorosity = true; + } + } + diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H index 3dee469329c..b26a6594aa1 100644 --- a/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H +++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H @@ -1,7 +1,13 @@ rho = thermo->rho(); -volScalarField rUA = 1.0/UEqn.A(); -U = rUA*UEqn.H(); +if (pressureImplicitPorosity) +{ + U = trTU()&UEqn.H(); +} +else +{ + U = trAU()*UEqn.H(); +} if (transonic) { @@ -11,17 +17,28 @@ if (transonic) fvc::interpolate(thermo->psi()) *( (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) +// + fvc::ddtPhiCorr(rUA, rho, U, phi) ) ); for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { + tmp<fvScalarMatrix> lapTerm; + + if (pressureImplicitPorosity) + { + lapTerm = fvm::laplacian(rho*trTU(), p); + } + else + { + lapTerm = fvm::laplacian(rho*trAU(), p); + } + fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rUA, p) + - lapTerm() == reactingParcels.Srho() ); @@ -37,19 +54,30 @@ if (transonic) else { phi = - fvc::interpolate(rho)* - ( + fvc::interpolate(rho) + *( (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) + + fvc::ddtPhiCorr(trAU(), rho, U, phi) ); for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { + tmp<fvScalarMatrix> lapTerm; + + if (pressureImplicitPorosity) + { + lapTerm = fvm::laplacian(rho*trTU(), p); + } + else + { + lapTerm = fvm::laplacian(rho*trAU(), p); + } + fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phi) - - fvm::laplacian(rho*rUA, p) + - lapTerm() == reactingParcels.Srho() ); @@ -66,7 +94,15 @@ else #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rUA*fvc::grad(p); +if (pressureImplicitPorosity) +{ + U -= trTU()&fvc::grad(p); +} +else +{ + U -= trAU()*fvc::grad(p); +} U.correctBoundaryConditions(); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); + diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C b/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C index 07838ea813e..77a67502dab 100644 --- a/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C +++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C @@ -36,60 +36,60 @@ Description #include "chemistrySolver.H" #include "ReactingCloudThermoTypes.H" #include "radiationModel.H" +#include "porousZones.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { -# include "setRootCase.H" - -# include "createTime.H" -# include "createMesh.H" -# include "readChemistryProperties.H" -# include "readEnvironmentalProperties.H" -# include "createFields.H" -# include "createClouds.H" -# include "readPISOControls.H" -# include "initContinuityErrs.H" -# include "readTimeControls.H" -# include "compressibleCourantNo.H" -# include "setInitialDeltaT.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "readChemistryProperties.H" + #include "readEnvironmentalProperties.H" + #include "createFields.H" + #include "createRadiationModel.H" + #include "createClouds.H" + #include "readPISOControls.H" + #include "initContinuityErrs.H" + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { -# include "readTimeControls.H" -# include "readPISOControls.H" -# include "compressibleCourantNo.H" -# include "setDeltaT.H" + #include "readTimeControls.H" + #include "readPISOControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; - Info << "Evolving reacting cloud" << endl; - reactingParcels.evolve(); reactingParcels.info(); -# include "chemistry.H" -# include "rhoEqn.H" + #include "chemistry.H" + #include "rhoEqn.H" // --- PIMPLE loop for (int ocorr=1; ocorr<=nOuterCorr; ocorr++) { -# include "UEqn.H" -# include "YEqn.H" + #include "UEqn.H" + #include "YEqn.H" // --- PISO loop for (int corr=1; corr<=nCorr; corr++) { -# include "hEqn.H" -# include "pEqn.H" + #include "hEqn.H" + #include "pEqn.H" } Info<< "T gas min/max = " << min(T).value() << ", " @@ -102,7 +102,7 @@ int main(int argc, char *argv[]) if (runTime.write()) { -# include "additionalOutput.H" + #include "additionalOutput.H" } Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" -- GitLab