diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C index 1160760bb6..dc1c09e744 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C @@ -71,7 +71,7 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readPISOControls.H" + #include "readPIMPLEControls.H" #include "readChemistryProperties.H" #include "readAdditionalSolutionControls.H" #include "readTimeControls.H" @@ -86,17 +86,27 @@ int main(int argc, char *argv[]) #include "timeScales.H" #include "rhoEqn.H" - #include "UEqn.H" - #include "YEqn.H" - #include "hsEqn.H" - // --- PISO loop - for (int corr=0; corrcorrect(); + turbulence->correct(); + + #include "UEqn.H" + #include "YEqn.H" + #include "hsEqn.H" + + // --- PISO loop + for (int corr=0; corr > mvConvection ( fv::convectionScheme::New @@ -10,7 +9,6 @@ tmp > mvConvection ) ); - if (solveSpecies) { label inertIndex = -1; diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H index 5ca55b84f6..cc3c74a535 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H @@ -73,12 +73,12 @@ dimensionedScalar rhoMax ( - mesh.solutionDict().subDict("PISO").lookup("rhoMax") + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax") ); dimensionedScalar rhoMin ( - mesh.solutionDict().subDict("PISO").lookup("rhoMin") + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin") ); Info<< "Creating turbulence model\n" << endl; @@ -116,11 +116,11 @@ dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0) ); - volScalarField invTauFlow + volScalarField invTau ( IOobject ( - "invTauFlow", + "invTau", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H index 9fb91f1d36..3e58534ea0 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H @@ -40,7 +40,20 @@ - fvm::laplacian(rho*rAU, p) ); - pEqn.solve(); + pEqn.solve + ( + mesh.solver + ( + p.select + ( + ( + finalIter + && corr == nCorr-1 + && nonOrth == nNonOrthCorr + ) + ) + ) + ); if (nonOrth == nNonOrthCorr) { diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H index f0763d8421..8136af1ad3 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H @@ -1,7 +1,9 @@ dictionary additional = mesh.solutionDict().subDict("additional"); -bool eWork = additional.lookupOrDefault("eWork", true); -bool hWork = additional.lookupOrDefault("hWork", true); +// pressure work term for enthalpy equation +bool pressureWork = additional.lookupOrDefault("pressureWork", true); +bool pressureWorkTimeDerivative = + additional.lookupOrDefault("pressureWorkTimeDerivative", true); // flag to activate solve transport for each specie (Y vector) bool solveSpecies = additional.lookupOrDefault("solveSpecies", true); diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H index 69d9a21731..c5dd4b7451 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H @@ -1,10 +1,13 @@ DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0); -if (eWork) +if (pressureWork) { - DpDt += -p*fvc::div(phi/fvc::interpolate(rho)); -} -if (hWork) -{ - DpDt += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p)); + surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho)); + + DpDt += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU); + + if (pressureWorkTimeDerivative) + { + DpDt += fvc::ddt(p); + } } diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H index b44271ab8c..2f8e036554 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H +++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H @@ -28,129 +28,40 @@ Info<< "Time scales min/max:" << endl; { // Cache old time scale field - tmp tinvTauFlow0 + tmp tinvTau0 ( new volScalarField ( IOobject ( - "invTauFlow0", + "invTau0", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE, false ), - invTauFlow + invTau ) ); - const volScalarField& invTauFlow0 = tinvTauFlow0(); + const volScalarField& invTau0 = tinvTau0(); // Flow time scale // ~~~~~~~~~~~~~~~ { - invTauFlow = + invTau = fvc::surfaceSum ( mag(phi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf()) ) /rho; - invTauFlow.max(1.0/maxDeltaT); + invTau.max(1.0/maxDeltaT); Info<< " Flow = " - << gMin(1/invTauFlow.internalField()) << ", " - << gMax(1/invTauFlow.internalField()) << endl; - } - - - // Mass source time scale - // ~~~~~~~~~~~~~~~~~~~~~~ - - { - scalarField tau - ( - runTime.deltaTValue()*mag(parcels.Srho() + massSource.SuTot()) - ); - - tau = alphaTauRho*rho/(tau + ROOTVSMALL); - - Info<< " Density = " - << min(maxDeltaT, gMin(tau)) << ", " - << min(maxDeltaT, gMax(tau)) << endl; - - invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); - } - - - // Momentum source time scale - // ~~~~~~~~~~~~~~~~~~~~~~~~~~ - - { -/* - // Method 1 - mag(U) limit using 'small' nominal velocity - scalarField tau - ( - runTime.deltaTValue() - *mag - ( - rho.dimensionedInternalField()*g - + parcels.UTrans()/(mesh.V()*runTime.deltaT()) - + momentumSource.Su() - ) - /rho - ); - - const scalar nomMagU(dimensionedScalar("1", dimVelocity, 1)); - tau = alphaTauU*(nomMagU + mag(U))/(tau + ROOTVSMALL); -*/ -/* - // Method 2 - based on fluxes and Co-like limit - volVectorField UEqnRhs - ( - IOobject - ( - "UEqnRhs", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh, - dimensionedVector("zero", dimDensity*dimAcceleration, vector::zero), - zeroGradientFvPatchVectorField::typeName - ); - - UEqnRhs.internalField() = - rho.dimensionedInternalField()*g - + parcels.UTrans()/(mesh.V()*runTime.deltaT()) - + momentumSource.Su(); - - surfaceScalarField phiSU - ( - "phiSU", - fvc::interpolate(runTime.deltaT()*UEqnRhs) & mesh.Sf() - ); - - scalarField tau - ( - alphaTauU*rho - /fvc::surfaceSum - ( - mag(phi + phiSU)*mesh.deltaCoeffs()/mesh.magSf() - + dimensionedScalar("SMALL", dimDensity/dimTime, ROOTVSMALL) - ) - ); - -*/ -/* - Info<< " Momentum = " << min(maxDeltaT, gMin(tau)) << ", " - << min(maxDeltaT, gMax(tau)) << endl; - - invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); -*/ + << gMin(1/invTau.internalField()) << ", " + << gMax(1/invTau.internalField()) << endl; } @@ -176,42 +87,7 @@ Info<< "Time scales min/max:" << endl; Info<< " Temperature = " << min(maxDeltaT, gMin(tau)) << ", " << min(maxDeltaT, gMax(tau)) << endl; - invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); - } - - - // Specie source time scale - // ~~~~~~~~~~~~~~~~~~~~~~~~ - - { - scalarField tau(mesh.nCells(), ROOTVGREAT); - forAll(Y, fieldI) - { - const volScalarField& Yi = Y[fieldI]; - const scalarField deltaYi - ( - runTime.deltaTValue() - *mag - ( - kappa*chemistry.RR(fieldI)() - + massSource.Su(fieldI) - + parcels.Srho(fieldI) - ) - /rho - ); - tau = - min - ( - tau, - alphaTauSpecie - /(deltaYi/(Yi + specieMaxUnbound) + ROOTVSMALL) - ); - } - - Info<< " Specie = " << min(maxDeltaT, gMin(tau)) << ", " - << min(maxDeltaT, gMax(tau)) << endl; - - invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); + invTau.internalField() = max(invTau.internalField(), 1/tau); } @@ -219,21 +95,21 @@ Info<< "Time scales min/max:" << endl; // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // - reduce as much as required for flow, but limit source contributions const dimensionedScalar deltaTRamp("deltaTRamp", dimless, 1/(1 + 0.2)); - invTauFlow = max(invTauFlow, invTauFlow0*deltaTRamp); - tinvTauFlow0.clear(); + invTau = max(invTau, invTau0*deltaTRamp); + tinvTau0.clear(); // Limit the largest time scale // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - invTauFlow.max(1/maxDeltaT); + invTau.max(1/maxDeltaT); // Spatially smooth the time scale field // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - fvc::smooth(invTauFlow, alphaTauSmooth); + fvc::smooth(invTau, alphaTauSmooth); - Info<< " Overall = " << min(1/invTauFlow).value() - << ", " << max(1/invTauFlow).value() << nl << endl; + Info<< " Overall = " << min(1/invTau).value() + << ", " << max(1/invTau).value() << nl << endl; } diff --git a/src/finiteVolume/cfdTools/general/include/addFinalIter.H b/src/finiteVolume/cfdTools/general/include/addFinalIter.H new file mode 100644 index 0000000000..5e4cac2d7b --- /dev/null +++ b/src/finiteVolume/cfdTools/general/include/addFinalIter.H @@ -0,0 +1,5 @@ +if (finalIter) +{ + mesh.data::add("finalIteration", true); +} + diff --git a/src/finiteVolume/cfdTools/general/include/removeFinalIter.H b/src/finiteVolume/cfdTools/general/include/removeFinalIter.H new file mode 100644 index 0000000000..db56e97557 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/include/removeFinalIter.H @@ -0,0 +1,4 @@ +if (finalIter) +{ + mesh.data::remove("finalIteration"); +}