From fcab778f57724ee90ef8bfe490749746fa70f206 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Wed, 25 Jan 2023 13:53:00 +0000 Subject: [PATCH] solvers: Adjust time step even if Courant number is zero This change means that even if the Courant number is zero, the time step is adjusted based on maximum time step settings and/or constraints specified by active fvModels. If none of these additional constraints are present then adjustment is deactivated. --- applications/solvers/foamMultiRun/setDeltaT.C | 15 ++++---- applications/solvers/foamRun/setDeltaT.C | 10 +++--- .../solvers/modules/VoFSolver/VoFSolver.C | 20 ++++------- .../solvers/modules/VoFSolver/setRDeltaT.C | 29 +++++----------- .../solvers/modules/fluidSolver/fluidSolver.C | 13 ++++--- .../modules/incompressibleFluid/setRDeltaT.C | 31 ++++++----------- .../modules/isothermalFluid/setRDeltaT.C | 34 ++++++------------- .../modules/multicomponentFluid/setRDeltaT.C | 13 +++---- .../multiphaseEuler/setRDeltaT.C | 31 ++++++----------- .../solvers/modules/shockFluid/setRDeltaT.C | 16 +++------ applications/solvers/modules/solid/solid.C | 13 ++++--- .../cfdTools/general/fvModels/fvModel.C | 4 +-- .../cfdTools/general/fvModels/fvModels.C | 4 +-- 13 files changed, 85 insertions(+), 148 deletions(-) diff --git a/applications/solvers/foamMultiRun/setDeltaT.C b/applications/solvers/foamMultiRun/setDeltaT.C index 3d2fd40d5b..ff1a79e684 100644 --- a/applications/solvers/foamMultiRun/setDeltaT.C +++ b/applications/solvers/foamMultiRun/setDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,17 +35,19 @@ void Foam::setDeltaT(Time& runTime, const PtrList& solvers) && runTime.controlDict().lookupOrDefault("adjustTimeStep", false) ) { - scalar deltaT = great; + bool transient = false; + scalar deltaT = vGreat; forAll(solvers, i) { if (solvers[i].transient()) { + transient = true; deltaT = min(deltaT, solvers[i].maxDeltaT()); } } - if (deltaT != great) + if (transient && deltaT < rootVGreat) { runTime.setDeltaT(min(runTime.deltaTValue(), deltaT)); } @@ -58,9 +60,8 @@ void Foam::adjustDeltaT(Time& runTime, const PtrList& solvers) // Update the time-step limited by the solvers maxDeltaT if (runTime.controlDict().lookupOrDefault("adjustTimeStep", false)) { - scalar deltaT = 1.2*runTime.deltaTValue(); - bool transient = false; + scalar deltaT = vGreat; forAll(solvers, i) { @@ -71,9 +72,9 @@ void Foam::adjustDeltaT(Time& runTime, const PtrList& solvers) } } - if (transient) + if (transient && deltaT < rootVGreat) { - runTime.setDeltaT(deltaT); + runTime.setDeltaT(min(1.2*runTime.deltaTValue(), deltaT)); Info<< "deltaT = " << runTime.deltaTValue() << endl; } } diff --git a/applications/solvers/foamRun/setDeltaT.C b/applications/solvers/foamRun/setDeltaT.C index bae4af87f3..7960fe8fb0 100644 --- a/applications/solvers/foamRun/setDeltaT.C +++ b/applications/solvers/foamRun/setDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,8 +32,9 @@ void Foam::setDeltaT(Time& runTime, const solver& solver) if ( runTime.timeIndex() == 0 - && solver.transient() && runTime.controlDict().lookupOrDefault("adjustTimeStep", false) + && solver.transient() + && solver.maxDeltaT() < rootVGreat ) { runTime.setDeltaT(min(runTime.deltaTValue(), solver.maxDeltaT())); @@ -46,8 +47,9 @@ void Foam::adjustDeltaT(Time& runTime, const solver& solver) // Update the time-step limited by the solver maxDeltaT if ( - solver.transient() - && runTime.controlDict().lookupOrDefault("adjustTimeStep", false) + runTime.controlDict().lookupOrDefault("adjustTimeStep", false) + && solver.transient() + && solver.maxDeltaT() < rootVGreat ) { runTime.setDeltaT(min(1.2*runTime.deltaTValue(), solver.maxDeltaT())); diff --git a/applications/solvers/modules/VoFSolver/VoFSolver.C b/applications/solvers/modules/VoFSolver/VoFSolver.C index da491b1469..9f7d46052a 100644 --- a/applications/solvers/modules/VoFSolver/VoFSolver.C +++ b/applications/solvers/modules/VoFSolver/VoFSolver.C @@ -172,25 +172,17 @@ Foam::solvers::VoFSolver::~VoFSolver() Foam::scalar Foam::solvers::VoFSolver::maxDeltaT() const { - const scalar maxAlphaCo - ( - runTime.controlDict().lookup("maxAlphaCo") - ); + const scalar maxAlphaCo = + runTime.controlDict().lookup("maxAlphaCo"); - const scalar deltaT = fluidSolver::maxDeltaT(); + scalar deltaT = fluidSolver::maxDeltaT(); if (alphaCoNum > small) { - return min - ( - deltaT, - maxAlphaCo/(alphaCoNum + small)*runTime.deltaTValue() - ); - } - else - { - return deltaT; + deltaT = min(deltaT, maxAlphaCo/alphaCoNum*runTime.deltaTValue()); } + + return deltaT; } diff --git a/applications/solvers/modules/VoFSolver/setRDeltaT.C b/applications/solvers/modules/VoFSolver/setRDeltaT.C index fb87416236..3b557e37b0 100644 --- a/applications/solvers/modules/VoFSolver/setRDeltaT.C +++ b/applications/solvers/modules/VoFSolver/setRDeltaT.C @@ -39,30 +39,19 @@ void Foam::solvers::VoFSolver::setRDeltaT() pimpleDict.lookupOrDefault("maxCo", 0.9) ); - const scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - - const scalar minDeltaT - ( - pimpleDict.lookupOrDefault("minDeltaT", small) - ); - const volScalarField rDeltaT0("rDeltaT0", rDeltaT); // Set the reciprocal time-step from the local Courant number // and maximum and minimum time-steps - rDeltaT.ref() = min - ( - 1/dimensionedScalar(dimTime, minDeltaT), - max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(mag(phi))()() - /((2*maxCo)*mesh.V()) - ) - ); + rDeltaT.ref() = fvc::surfaceSum(mag(phi))()()/((2*maxCo)*mesh.V()); + if (pimpleDict.found("maxDeltaT")) + { + rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); + } + if (pimpleDict.found("minDeltaT")) + { + rDeltaT.min(1/pimpleDict.lookup("minDeltaT")); + } Info<< "Flow time scale min/max = " << gMin(1/rDeltaT.primitiveField()) diff --git a/applications/solvers/modules/fluidSolver/fluidSolver.C b/applications/solvers/modules/fluidSolver/fluidSolver.C index 3c3b2577e4..ed9c49a699 100644 --- a/applications/solvers/modules/fluidSolver/fluidSolver.C +++ b/applications/solvers/modules/fluidSolver/fluidSolver.C @@ -48,7 +48,7 @@ void Foam::solvers::fluidSolver::readControls() runTime.controlDict().lookupOrDefault("maxCo", 1.0); maxDeltaT_ = - runTime.controlDict().lookupOrDefault("maxDeltaT", great); + runTime.controlDict().lookupOrDefault("maxDeltaT", vGreat); correctPhi = pimple.dict().lookupOrDefault ( @@ -205,15 +205,14 @@ Foam::solvers::fluidSolver::~fluidSolver() Foam::scalar Foam::solvers::fluidSolver::maxDeltaT() const { + scalar deltaT = min(fvModels().maxDeltaT(), maxDeltaT_); + if (CoNum > small) { - const scalar deltaT = maxCo*runTime.deltaTValue()/CoNum; - return min(min(deltaT, fvModels().maxDeltaT()), maxDeltaT_); - } - else - { - return runTime.deltaTValue(); + deltaT = min(deltaT, maxCo/CoNum*runTime.deltaTValue()); } + + return deltaT; } diff --git a/applications/solvers/modules/incompressibleFluid/setRDeltaT.C b/applications/solvers/modules/incompressibleFluid/setRDeltaT.C index f0289ea3ce..ae09490dda 100644 --- a/applications/solvers/modules/incompressibleFluid/setRDeltaT.C +++ b/applications/solvers/modules/incompressibleFluid/setRDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -49,30 +49,19 @@ void Foam::solvers::incompressibleFluid::setRDeltaT() pimpleDict.lookupOrDefault("rDeltaTDampingCoeff", 1.0) ); - const scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - - const scalar minDeltaT - ( - pimpleDict.lookupOrDefault("minDeltaT", small) - ); - const volScalarField rDeltaT0("rDeltaT0", rDeltaT); // Set the reciprocal time-step from the local Courant number // and maximum and minimum time-steps - rDeltaT.ref() = min - ( - 1/dimensionedScalar(dimTime, minDeltaT), - max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(mag(phi))()() - /((2*maxCo)*mesh.V()) - ) - ); + rDeltaT.ref() = fvc::surfaceSum(mag(phi))()()/((2*maxCo)*mesh.V()); + if (pimpleDict.found("maxDeltaT")) + { + rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); + } + if (pimpleDict.found("minDeltaT")) + { + rDeltaT.min(1/pimpleDict.lookup("minDeltaT")); + } // Update the boundary values of the reciprocal time-step rDeltaT.correctBoundaryConditions(); diff --git a/applications/solvers/modules/isothermalFluid/setRDeltaT.C b/applications/solvers/modules/isothermalFluid/setRDeltaT.C index 7998fedd0f..e6de67e232 100644 --- a/applications/solvers/modules/isothermalFluid/setRDeltaT.C +++ b/applications/solvers/modules/isothermalFluid/setRDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,30 +47,19 @@ void Foam::solvers::isothermalFluid::setRDeltaT() pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.02) ); - const scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - - const scalar minDeltaT - ( - pimpleDict.lookupOrDefault("minDeltaT", small) - ); - const volScalarField rDeltaT0("rDeltaT0", rDeltaT); // Set the reciprocal time-step from the local Courant number // and maximum and minimum time-steps - rDeltaT.ref() = min - ( - 1/dimensionedScalar(dimTime, minDeltaT), - max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(mag(phi))()() - /((2*maxCo)*mesh.V()*rho()) - ) - ); + rDeltaT.ref() = fvc::surfaceSum(mag(phi))()()/((2*maxCo)*mesh.V()*rho()); + if (pimpleDict.found("maxDeltaT")) + { + rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); + } + if (pimpleDict.found("minDeltaT")) + { + rDeltaT.min(1/pimpleDict.lookup("minDeltaT")); + } if (pimple.transonic()) { @@ -83,8 +72,7 @@ void Foam::solvers::isothermalFluid::setRDeltaT() rDeltaT.ref() = max ( rDeltaT(), - fvc::surfaceSum(mag(phid))()() - /((2*maxCo)*mesh.V()*psi()) + fvc::surfaceSum(mag(phid))()()/((2*maxCo)*mesh.V()*psi()) ); } diff --git a/applications/solvers/modules/multicomponentFluid/setRDeltaT.C b/applications/solvers/modules/multicomponentFluid/setRDeltaT.C index f588fb4a67..ef1c640738 100644 --- a/applications/solvers/modules/multicomponentFluid/setRDeltaT.C +++ b/applications/solvers/modules/multicomponentFluid/setRDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,19 +44,14 @@ void Foam::solvers::multicomponentFluid::setRDeltaT() // Maximum flow Courant number const scalar maxCo(pimpleDict.lookup("maxCo")); + // Set the reciprocal time-step from the local Courant number + // and maximum and minimum time-steps rDeltaT.ref() = - ( - fvc::surfaceSum(mag(phi))()() - /((2*maxCo)*mesh.V()*rho()) - ); - - // Limit the largest time step + fvc::surfaceSum(mag(phi))()()/((2*maxCo)*mesh.V()*rho()); if (pimpleDict.found("maxDeltaT")) { rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); } - - // Limit the smallest time step if (pimpleDict.found("minDeltaT")) { rDeltaT.min(1/pimpleDict.lookup("minDeltaT")); diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/setRDeltaT.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/setRDeltaT.C index 99bbeba635..f2c290bfa3 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/setRDeltaT.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/setRDeltaT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -40,16 +40,6 @@ void Foam::solvers::multiphaseEuler::setRDeltaT() pimpleDict.lookupOrDefault("maxCo", 0.2) ); - const scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - - const scalar minDeltaT - ( - pimpleDict.lookupOrDefault("minDeltaT", small) - ); - const scalar rDeltaTSmoothingCoeff ( pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.02) @@ -64,16 +54,15 @@ void Foam::solvers::multiphaseEuler::setRDeltaT() // Set the reciprocal time-step from the local Courant number // and maximum and minimum time-steps - rDeltaT.ref() = min - ( - 1/dimensionedScalar(dimTime, minDeltaT), - max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(maxPhi)()() - /((2*maxCo)*mesh.V()) - ) - ); + rDeltaT.ref() = fvc::surfaceSum(maxPhi)()()/((2*maxCo)*mesh.V()); + if (pimpleDict.found("maxDeltaT")) + { + rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); + } + if (pimpleDict.found("minDeltaT")) + { + rDeltaT.min(1/pimpleDict.lookup("minDeltaT")); + } // Update the boundary values of the reciprocal time-step rDeltaT.correctBoundaryConditions(); diff --git a/applications/solvers/modules/shockFluid/setRDeltaT.C b/applications/solvers/modules/shockFluid/setRDeltaT.C index 5e1902f807..d0d26a7a73 100644 --- a/applications/solvers/modules/shockFluid/setRDeltaT.C +++ b/applications/solvers/modules/shockFluid/setRDeltaT.C @@ -48,18 +48,12 @@ void Foam::solvers::shockFluid::setRDeltaT(const surfaceScalarField& amaxSf) ) ); - const scalar maxDeltaT - ( - pimpleDict.lookupOrDefault("maxDeltaT", great) - ); - // Set the reciprocal time-step from the local Courant number - rDeltaT.ref() = max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(amaxSf)()() - /((2*maxCo)*mesh.V()) - ); + rDeltaT.ref() = fvc::surfaceSum(amaxSf)()()/((2*maxCo)*mesh.V()); + if (pimpleDict.found("maxDeltaT")) + { + rDeltaT.max(1/pimpleDict.lookup("maxDeltaT")); + } // Update the boundary values of the reciprocal time-step rDeltaT.correctBoundaryConditions(); diff --git a/applications/solvers/modules/solid/solid.C b/applications/solvers/modules/solid/solid.C index 2549b1a1f2..4ad044677a 100644 --- a/applications/solvers/modules/solid/solid.C +++ b/applications/solvers/modules/solid/solid.C @@ -47,7 +47,7 @@ void Foam::solvers::solid::readControls() runTime.controlDict().lookupOrDefault("maxDi", 1.0); maxDeltaT_ = - runTime.controlDict().lookupOrDefault("maxDeltaT", great); + runTime.controlDict().lookupOrDefault("maxDeltaT", vGreat); } @@ -131,15 +131,14 @@ Foam::solvers::solid::~solid() Foam::scalar Foam::solvers::solid::maxDeltaT() const { + scalar deltaT = min(fvModels().maxDeltaT(), maxDeltaT_); + if (DiNum > small) { - const scalar deltaT = maxDi*runTime.deltaTValue()/DiNum; - return min(min(deltaT, fvModels().maxDeltaT()), maxDeltaT_); - } - else - { - return maxDeltaT_; + deltaT = min(deltaT, maxDi/DiNum*runTime.deltaTValue()); } + + return deltaT; } diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModel.C b/src/finiteVolume/cfdTools/general/fvModels/fvModel.C index 467dfda744..d80732487e 100644 --- a/src/finiteVolume/cfdTools/general/fvModels/fvModel.C +++ b/src/finiteVolume/cfdTools/general/fvModels/fvModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2021-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2021-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -148,7 +148,7 @@ bool Foam::fvModel::addsSupToField(const word& fieldName) const Foam::scalar Foam::fvModel::maxDeltaT() const { - return great; + return vGreat; } diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModels.C b/src/finiteVolume/cfdTools/general/fvModels/fvModels.C index 436951decc..04f3e1512d 100644 --- a/src/finiteVolume/cfdTools/general/fvModels/fvModels.C +++ b/src/finiteVolume/cfdTools/general/fvModels/fvModels.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2021-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2021-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -241,7 +241,7 @@ Foam::scalar Foam::fvModels::maxDeltaT() const { const PtrListDictionary& modelList(*this); - scalar maxDeltaT = great; + scalar maxDeltaT = vGreat; forAll(modelList, i) {