From 5d4584a24c9ef3394d823178145d1c7cf5e63c19 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 28 Feb 2018 12:21:37 +0100 Subject: [PATCH 1/2] Consistency fix Ratio of old volumes to current volumes needs to be taken into account when creating time-step consistent face flux in ddt schemes. --- .../CoEulerDdtScheme/CoEulerDdtScheme.C | 30 ++++++++- .../CrankNicolsonDdtScheme.C | 39 +++++++++--- .../EulerDdtScheme/EulerDdtScheme.C | 30 ++++++++- .../ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C | 30 ++++++++- .../backwardDdtScheme/backwardDdtScheme.C | 62 ++++++++++++++++--- .../steadyInertialDdtScheme.C | 30 ++++++++- 6 files changed, 202 insertions(+), 19 deletions(-) diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C index c75f9dfd8..2a54fbd6f 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C @@ -660,7 +660,35 @@ CoEulerDdtScheme::fvcDdtConsistentPhiCorr const surfaceScalarField& rAUf ) { - return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT(); + tmp toldTimeFlux = + (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT(); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes + volScalarField V0ByV + ( + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct the flux with interpolated volume ratio + toldTimeFlux() *= fvc::interpolate(V0ByV); + } + + return toldTimeFlux; } diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C index e632f8119..dd58b3996 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C @@ -1201,16 +1201,39 @@ CrankNicolsonDdtScheme::fvcDdtConsistentPhiCorr - offCentre_(faceUDdt0()) ); } + + // Calculate old time flux + fluxFieldType oldTimeFlux = + rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & faceU.oldTime()); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes + volScalarField V0ByV + ( + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct the flux with interpolated volume ratio + oldTimeFlux *= fvc::interpolate(V0ByV); + } return - rAUf* - ( - mesh().Sf() - & ( - rDtCoef_(faceUDdt0)*faceU.oldTime() - + offCentre_(faceUDdt0()) - ) - ); + oldTimeFlux + + rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & offCentre_(faceUDdt0())); } diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C index 96ade73df..9f6187fa3 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C @@ -522,7 +522,35 @@ EulerDdtScheme::fvcDdtConsistentPhiCorr const surfaceScalarField& rAUf ) { - return (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT(); + tmp toldTimeFlux = + (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT(); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes + volScalarField V0ByV + ( + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct the flux with interpolated volume ratio + toldTimeFlux() *= fvc::interpolate(V0ByV); + } + + return toldTimeFlux; } diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C index e33ff3956..0452c17c6 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C @@ -665,7 +665,35 @@ SLTSDdtScheme::fvcDdtConsistentPhiCorr const surfaceScalarField& rAUf ) { - return (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT()); + tmp toldTimeFlux = + (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT()); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes + volScalarField V0ByV + ( + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct the flux with interpolated volume ratio + toldTimeFlux() *= fvc::interpolate(V0ByV); + } + + return toldTimeFlux; } diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C index 097d4eaf9..a45577e51 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C @@ -732,18 +732,66 @@ backwardDdtScheme::fvcDdtConsistentPhiCorr const scalar rDeltaT = 1.0/deltaT; + // Note: minus sign in gamma coefficient so we can simply add the fluxes + // together at the end const dimensionedScalar beta("beta", dimless/dimTime, coefft0*rDeltaT); const dimensionedScalar gamma("gamma", dimless/dimTime, -coefft00*rDeltaT); - return - rAUf* + // Calculate old and old-old flux contributions + fluxFieldType oldTimeFlux = + beta*rAUf*(mesh().Sf() & faceU.oldTime()); + fluxFieldType oldOldTimeFlux = + gamma*rAUf*(mesh().Sf() & faceU.oldTime().oldTime()); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes for old flux contribution + volScalarField V0ByV ( - mesh().Sf() - & ( - beta*faceU.oldTime() - + gamma*faceU.oldTime().oldTime() - ) + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct old time flux contribution + oldTimeFlux *= fvc::interpolate(V0ByV); + + + // Also need to take into account the ratio between old-old and current + // cell volumes for old-old time flux contribution + volScalarField V00ByV + ( + IOobject + ( + "V00ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V00ByV.internalField() = mesh().V00()/mesh().V(); + V00ByV.correctBoundaryConditions(); + + // Correct old-old time flux contribution + oldOldTimeFlux *= fvc::interpolate(V00ByV); + } + + return oldTimeFlux + oldOldTimeFlux; } diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C index ca59ca10e..93d27f608 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C @@ -667,7 +667,35 @@ steadyInertialDdtScheme::fvcDdtConsistentPhiCorr const surfaceScalarField& rAUf ) { - return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT(); + tmp toldTimeFlux = + (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT(); + + if (mesh().moving()) + { + // Mesh is moving, need to take into account the ratio between old and + // current cell volumes + volScalarField V0ByV + ( + IOobject + ( + "V0ByV", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + V0ByV.internalField() = mesh().V0()/mesh().V(); + V0ByV.correctBoundaryConditions(); + + // Correct the flux with interpolated volume ratio + toldTimeFlux() *= fvc::interpolate(V0ByV); + } + + return toldTimeFlux; } From 0e0debd52941268767ef13edbb515d2a277fe4ca Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 28 Feb 2018 12:22:57 +0100 Subject: [PATCH 2/2] Updates to verificationSuite cases --- .../incompressible/icoFoam/consistencyTest/cavity/Allrun | 2 +- .../icoFoam/consistencyTest/cavity/explanationAndReference | 2 +- .../incompressible/simpleFoam/consistencyTest/cavity/Allrun | 2 +- .../simpleFoam/consistencyTest/cavity/explanationAndReference | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/Allrun b/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/Allrun index ef103931e..e24fa1266 100755 --- a/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/Allrun +++ b/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/Allrun @@ -18,4 +18,4 @@ done # Print out the converged pressure for all time steps for visual check whether # the solution does not depend on the time step -tail -n 1 */probes/0/p +tail -n 1 */postProcessing/probes/0/p diff --git a/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/explanationAndReference b/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/explanationAndReference index f4c898a5d..61a210994 100644 --- a/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/explanationAndReference +++ b/validationAndVerificationSuite/verification/incompressible/icoFoam/consistencyTest/cavity/explanationAndReference @@ -5,6 +5,6 @@ simulations are performed with four different time steps spanning four orders of magnitude: 0.01, 0.001, 0.0001 and 0.00001 s. Tolerances for all equations are very small, yielding extremely small (O(1e-11)) differences in converged pressure field. If the differences are larger - it means that the converged -solution is not independent to time step size. +solution depends on the time step size. Author: Vuko Vukcevic, vuko.vukcevic@fsb.hr diff --git a/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/Allrun b/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/Allrun index 145805799..44fc5d665 100755 --- a/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/Allrun +++ b/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/Allrun @@ -18,4 +18,4 @@ done # Print out the converged pressure for all relaxation factors for visual check # whether the solution does not depend on the under-relaxation factors -tail -n 1 */probes/0/p +tail -n 1 */postProcessing/probes/0/p diff --git a/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/explanationAndReference b/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/explanationAndReference index ddff443a2..5aacb29e3 100644 --- a/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/explanationAndReference +++ b/validationAndVerificationSuite/verification/incompressible/simpleFoam/consistencyTest/cavity/explanationAndReference @@ -9,6 +9,6 @@ simulations are performed with five different under-relaxation pairs: 5. alphaU = 0.4, alphap = 0.6 Tolerances for all equations are very small, yielding extremely small (O(1e-11)) differences in converged pressure field. If the differences are larger - it -means that the converged solution is not independent to relaxation factors. +means that the converged solution is dependend on relaxation factors. Author: Vuko Vukcevic, vuko.vukcevic@fsb.hr