From bab72d5139d85d21768f4c0ab40345ca27623b7c Mon Sep 17 00:00:00 2001 From: Johan Roenby Date: Tue, 12 Jun 2018 23:14:10 +0200 Subject: [PATCH] ENH: multiple updates to interIsoFoam related code - Updated tutorial headers - Added copyright note to isoAdvector src - Removed outcommented code lines in interIsoFoam solver - Removed all LTS from interIsoFoam since this is not currently supported - Confirmed that discInConstantFlow gives identical results with N subCylces and time step N*dt - Confirmed that this also holds when nOuterCorrectors > 1. --- .../multiphase/interIsoFoam/alphaControls.H | 26 ---- .../multiphase/interIsoFoam/alphaEqn.H | 99 ++++++------- .../interIsoFoam/alphaEqnSubCycle.H | 8 -- .../multiphase/interIsoFoam/alphaSuSp.H | 3 - .../multiphase/interIsoFoam/createFields.H | 20 --- .../multiphase/interIsoFoam/interIsoFoam.C | 28 +--- .../multiphase/interIsoFoam/setRDeltaT.H | 136 ------------------ .../isoAdvection/isoAdvection/isoAdvection.C | 1 + .../isoAdvection/isoAdvection/isoAdvection.H | 1 + .../isoAdvection/isoCutCell/isoCutCell.C | 1 + .../isoAdvection/isoCutCell/isoCutCell.H | 1 + .../isoAdvection/isoCutFace/isoCutFace.C | 1 + .../isoAdvection/isoCutFace/isoCutFace.H | 1 + .../damBreakWithObstacle/0.orig/U | 4 +- .../damBreakWithObstacle/0.orig/alpha.water | 4 +- .../damBreakWithObstacle/0.orig/p_rgh | 4 +- .../constant/dynamicMeshDict | 4 +- .../damBreakWithObstacle/constant/g | 4 +- .../constant/transportProperties | 4 +- .../constant/turbulenceProperties | 4 +- .../damBreakWithObstacle/system/blockMeshDict | 4 +- .../damBreakWithObstacle/system/controlDict | 4 +- .../system/decomposeParDict | 4 +- .../damBreakWithObstacle/system/fvSchemes | 4 +- .../damBreakWithObstacle/system/fvSolution | 4 +- .../damBreakWithObstacle/system/setFieldsDict | 4 +- .../damBreakWithObstacle/system/topoSetDict | 4 +- .../discInConstantFlowCyclicBCs/0.orig/U | 3 +- .../constant/dynamicMeshDict | 3 +- .../discInConstantFlowCyclicBCs/constant/g | 3 +- .../constant/transportProperties | 3 +- .../constant/turbulenceProperties | 3 +- .../system/blockMeshDict | 3 +- .../system/controlDict | 3 +- .../system/decomposeParDict | 3 +- .../system/fvSchemes | 3 +- .../system/fvSolution | 3 +- .../system/setAlphaFieldDict | 3 +- 38 files changed, 112 insertions(+), 303 deletions(-) delete mode 100644 applications/solvers/multiphase/interIsoFoam/alphaSuSp.H delete mode 100644 applications/solvers/multiphase/interIsoFoam/setRDeltaT.H diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H index ebe73c3c2c..db77d94af4 100644 --- a/applications/solvers/multiphase/interIsoFoam/alphaControls.H +++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H @@ -1,29 +1,3 @@ const dictionary& alphaControls = mesh.solverDict(alpha1.name()); -//label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr"))); - label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); - -/* -bool MULESCorr(alphaControls.lookupOrDefault("MULESCorr", false)); - -// Apply the compression correction from the previous iteration -// Improves efficiency for steady-simulations but can only be applied -// once the alpha field is reasonably steady, i.e. fully developed -bool alphaApplyPrevCorr -( - alphaControls.lookupOrDefault("alphaApplyPrevCorr", false) -); - -// Isotropic compression coefficient -scalar icAlpha -( - alphaControls.lookupOrDefault("icAlpha", 0) -); - -// Shear compression coefficient -scalar scAlpha -( - alphaControls.lookupOrDefault("scAlpha", 0) -); -*/ diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H index 1bda9eef18..a73391a996 100644 --- a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H @@ -1,50 +1,53 @@ +{ + // Note for AMR implementation: + // At this point we have just entered the new time step, + // the mesh has been refined and the alpha, phi and U contain + // the field values at the beginning of this time step mapped + // to the new mesh. + + // The above assumes that we are in firstIter() of the outer + // corrector loop. If we are in any subsequent iter of this loop + // the alpha1, U and phi will be overwritten with the new time step + // values but still on the same mesh. + + + if (pimple.firstIter()) { - // Note for AMR implementation: - // At this point we have just entered the new time step, - // the mesh has been refined and the alpha, phi and U contain - // the field values at the beginning of this time step mapped - // to the new mesh. - - // The above assumes that we are in firstIter() of the outer - // corrector loop. If we are in any subsequent iter of this loop - // the alpha1, U and phi will be overwritten with the new time step - // values but still on the same mesh. - - - if (pimple.firstIter()) - { - // Note: This assumes moveMeshOuterCorrectors is false - alpha10 = alpha1; - U0 = U; - phi0 = phi; - advector.advect(); - - #include "rhofs.H" - rhoPhi = advector.getRhoPhi(rho1f, rho2f); - } - else - { - alpha1 = alpha10; - // Temporarily setting U and phi to average of old and new value - // Note: Illegal additions if mesh is topoChanging - // (e.g. if moveMeshOuterCorrectors and AMR) - U = 0.5*U0 + 0.5*U; - phi = 0.5*phi0 + 0.5*phi; - isoAdvection advector(alpha1, phi, U); - advector.advect(); - #include "rhofs.H" - rhoPhi = advector.getRhoPhi(rho1f, rho2f); - // Resetting U and phi to the new value - U = 2.0*U - U0; - phi = 2.0*phi - phi0; - } - - alpha2 = 1.0 - alpha1; - mixture.correct(); - + // Note: This assumes moveMeshOuterCorrectors is false + // Saving field values before advection in first PIMPLE iteration + alpha10 = alpha1; + U0 = U; + phi0 = phi; } - Info<< "Phase-1 volume fraction = " - << alpha1.weightedAverage(mesh.Vsc()).value() - << " Min(" << alpha1.name() << ") = " << min(alpha1).value() - << " Max(" << alpha1.name() << ") = " << max(alpha1).value() - << endl; + else + { + // Resetting alpha1 to value before advection in first PIMPLE iteration + alpha1 = alpha10; + // Temporarily setting U and phi to average of old and new value + // Note: Illegal additions if mesh is topoChanging + // (e.g. if moveMeshOuterCorrectors and AMR) + U = 0.5*U0 + 0.5*U; + phi = 0.5*phi0 + 0.5*phi; + } + + advector.advect(); + #include "rhofs.H" + rhoPhi = advector.getRhoPhi(rho1f, rho2f); + + if (!pimple.firstIter()) + { + // Resetting U and phi to the new value + U = 2.0*U - U0; + phi = 2.0*phi - phi0; + } + + alpha2 = 1.0 - alpha1; + mixture.correct(); + +} + +Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(" << alpha1.name() << ") = " << min(alpha1).value() + << " Max(" << alpha1.name() << ") = " << max(alpha1).value() + << endl; diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H index c53d6e106b..3c8ba18f2b 100644 --- a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H @@ -13,14 +13,6 @@ if (nAlphaSubCycles > 1) dimensionedScalar(rhoPhi.dimensions(), Zero) ); - tmp trSubDeltaT; - - if (LTS) - { - trSubDeltaT = - fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles); - } - for ( subCycle alphaSubCycle(alpha1, nAlphaSubCycles); diff --git a/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H b/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H deleted file mode 100644 index 0b8e05e187..0000000000 --- a/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H +++ /dev/null @@ -1,3 +0,0 @@ -zeroField Su; -zeroField Sp; -zeroField divU; diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H index 02269e379a..76b9938107 100644 --- a/applications/solvers/multiphase/interIsoFoam/createFields.H +++ b/applications/solvers/multiphase/interIsoFoam/createFields.H @@ -121,23 +121,6 @@ if (p_rgh.needReference()) mesh.setFluxRequired(p_rgh.name()); mesh.setFluxRequired(alpha1.name()); -/* -// MULES compressed flux is registered in case scalarTransport FO needs it. -surfaceScalarField alphaPhiUn -( - IOobject - ( - "alphaPhiUn", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar(phi.dimensions(), Zero) -); -*/ - #include "createMRF.H" #include "createFvOptions.H" @@ -155,7 +138,6 @@ surfaceScalarField phi0 phi ); - volVectorField U0 ( IOobject @@ -169,7 +151,6 @@ volVectorField U0 U ); - volScalarField alpha10 ( IOobject @@ -183,5 +164,4 @@ volScalarField alpha10 alpha1 ); - isoAdvection advector(alpha1, phi, U); diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C index 024e502633..1593cae755 100644 --- a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C +++ b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C @@ -75,17 +75,13 @@ int main(int argc, char *argv[]) #include "initContinuityErrs.H" #include "createDyMControls.H" #include "createFields.H" -// #include "createAlphaFluxes.H" #include "initCorrectPhi.H" #include "createUfIfPresent.H" turbulence->validate(); - if (!LTS) - { - #include "CourantNo.H" - #include "setInitialDeltaT.H" - } + #include "CourantNo.H" + #include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -93,17 +89,9 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readDyMControls.H" - - if (LTS) - { - #include "setRDeltaT.H" - } - else - { - #include "CourantNo.H" - #include "alphaCourantNo.H" - #include "setDeltaT.H" - } + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" runTime++; @@ -118,12 +106,6 @@ int main(int argc, char *argv[]) if (mesh.changing()) { - // Do not apply previous time-step mesh compression flux - // if the mesh topology changed - if (mesh.topoChanging()) - { -// talphaPhi1Corr0.clear(); - } gh = (g & mesh.C()) - ghRef; ghf = (g & mesh.Cf()) - ghRef; diff --git a/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H deleted file mode 100644 index 26b237e795..0000000000 --- a/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H +++ /dev/null @@ -1,136 +0,0 @@ -{ - volScalarField& rDeltaT = trDeltaT.ref(); - - const dictionary& pimpleDict = pimple.dict(); - - scalar maxCo - ( - pimpleDict.lookupOrDefault("maxCo", 0.9) - ); - - scalar maxAlphaCo - ( - pimpleDict.lookupOrDefault("maxAlphaCo", 0.2) - ); - - scalar rDeltaTSmoothingCoeff - ( - pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.1) - ); - - label nAlphaSpreadIter - ( - pimpleDict.lookupOrDefault