From db568490b715c7370f56f3b6adb34e7983b9b2a5 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 11:23:41 +0200 Subject: [PATCH 1/7] incompressible/icoFoam, solver + tutorials --- .../solvers/incompressible/icoFoam/icoFoam.C | 37 +++++++++++++------ .../icoFoam/elbow/system/fvSchemes | 2 +- .../0_orig/U | 6 --- .../twoBlocksMixingPlane_dirY_spanZ/0_orig/U | 6 --- .../twoBlocksMixingPlane_dirZ_spanY/0_orig/U | 6 --- 5 files changed, 26 insertions(+), 31 deletions(-) diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C index 164f5f6e9..5073f4e61 100644 --- a/applications/solvers/incompressible/icoFoam/icoFoam.C +++ b/applications/solvers/incompressible/icoFoam/icoFoam.C @@ -26,6 +26,10 @@ Application Description Transient solver for incompressible, laminar flow of Newtonian fluids. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved \*---------------------------------------------------------------------------*/ @@ -53,32 +57,36 @@ int main(int argc, char *argv[]) # include "readPISOControls.H" # include "CourantNo.H" - fvVectorMatrix UEqn + // Convection-diffusion matrix + fvVectorMatrix HUEqn ( - fvm::ddt(U) - + fvm::div(phi, U) + fvm::div(phi, U) - fvm::laplacian(nu, U) ); - solve(UEqn == -fvc::grad(p)); + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + solve(ddtUEqn + HUEqn == -fvc::grad(p)); + + // Prepare clean Ap without time derivative contribution + // HJ, 26/Oct/2015 + volScalarField aU = HUEqn.A(); // --- PISO loop for (int corr = 0; corr < nCorr; corr++) { - volScalarField rUA = 1.0/UEqn.A(); - - U = rUA*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, U, phi); + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); adjustPhi(phi, U, p); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( - fvm::laplacian(rUA, p) == fvc::div(phi) + fvm::laplacian(1/aU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); @@ -92,7 +100,12 @@ int main(int argc, char *argv[]) # include "continuityErrs.H" - U -= rUA*fvc::grad(p); + // Note: cannot call H(U) here because the velocity is not complete + // HJ, 22/Jan/2016 + U = 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ); U.correctBoundaryConditions(); } diff --git a/tutorials/incompressible/icoFoam/elbow/system/fvSchemes b/tutorials/incompressible/icoFoam/elbow/system/fvSchemes index 0cd903c76..7c96ea2c4 100644 --- a/tutorials/incompressible/icoFoam/elbow/system/fvSchemes +++ b/tutorials/incompressible/icoFoam/elbow/system/fvSchemes @@ -28,7 +28,7 @@ gradSchemes divSchemes { default none; - div(phi,U) Gauss limitedLinearV 1; + div(phi,U) Gauss linearUpwind leastSquares; } laplacianSchemes diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U index c9b0cd956..e5afe4b5d 100644 --- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U +++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlaneMismatch_dirY_spanZ/0_orig/U @@ -54,12 +54,6 @@ boundaryField type mixingPlane; } - top - { - type fixedValue; - value uniform (0 0 0); - } - frontAndBack { type empty; diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U index c9b0cd956..e5afe4b5d 100644 --- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U +++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirY_spanZ/0_orig/U @@ -54,12 +54,6 @@ boundaryField type mixingPlane; } - top - { - type fixedValue; - value uniform (0 0 0); - } - frontAndBack { type empty; diff --git a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U index c9b0cd956..e5afe4b5d 100644 --- a/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U +++ b/tutorials/incompressible/icoFoam/mixingPlane/twoBlocksMixingPlane_dirZ_spanY/0_orig/U @@ -54,12 +54,6 @@ boundaryField type mixingPlane; } - top - { - type fixedValue; - value uniform (0 0 0); - } - frontAndBack { type empty; From 1670cebef1fcc611b59ff5b95959f967ddc9471d Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 12:55:56 +0200 Subject: [PATCH 2/7] incompressible/nonNewtonianIcoFoam, solver + tutorials --- .../nonNewtonianIcoFoam/nonNewtonianIcoFoam.C | 39 ++++++++++++------- .../offsetCylinder/system/fvSchemes | 2 +- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C index 6420661d7..624c47c6c 100644 --- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C +++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C @@ -26,6 +26,10 @@ Application Description Transient solver for incompressible, laminar flow of non-Newtonian fluids. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved \*---------------------------------------------------------------------------*/ @@ -39,7 +43,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" -# include "createMeshNoClear.H" +# include "createMesh.H" # include "createFields.H" # include "initContinuityErrs.H" @@ -56,32 +60,36 @@ int main(int argc, char *argv[]) fluid.correct(); - fvVectorMatrix UEqn + // Convection-diffusion matrix + fvVectorMatrix HUEqn ( - fvm::ddt(U) - + fvm::div(phi, U) + fvm::div(phi, U) - fvm::laplacian(fluid.nu(), U) ); - solve(UEqn == -fvc::grad(p)); + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + solve(ddtUEqn + HUEqn == -fvc::grad(p)); + + // Prepare clean Ap without time derivative contribution + // HJ, 26/Oct/2015 + volScalarField aU = HUEqn.A(); // --- PISO loop for (int corr = 0; corr < nCorr; corr++) { - volScalarField rUA = 1.0/UEqn.A(); - - U = rUA*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, U, phi); + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); adjustPhi(phi, U, p); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( - fvm::laplacian(rUA, p) == fvc::div(phi) + fvm::laplacian(1/aU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); @@ -95,7 +103,12 @@ int main(int argc, char *argv[]) # include "continuityErrs.H" - U -= rUA*fvc::grad(p); + // Note: cannot call H(U) here because the velocity is not complete + // HJ, 22/Jan/2016 + U = 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ); U.correctBoundaryConditions(); } diff --git a/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/fvSchemes b/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/fvSchemes index 1a04dfafd..5e09611e2 100644 --- a/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/fvSchemes +++ b/tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder/system/fvSchemes @@ -28,7 +28,7 @@ gradSchemes divSchemes { default none; - div(phi,U) Gauss linear; + div(phi,U) Gauss linearUpwind leastSquares; } laplacianSchemes From edb598687a28044970dc9476c52633a99ea7f17b Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 13:53:56 +0200 Subject: [PATCH 3/7] incompressible/channelFoam, solver + tutorials --- .../incompressible/channelFoam/channelFoam.C | 43 ++++++++++++------- .../channelFoam/channel395/system/fvSchemes | 2 +- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/applications/solvers/incompressible/channelFoam/channelFoam.C b/applications/solvers/incompressible/channelFoam/channelFoam.C index ba2b77038..40730616f 100644 --- a/applications/solvers/incompressible/channelFoam/channelFoam.C +++ b/applications/solvers/incompressible/channelFoam/channelFoam.C @@ -26,6 +26,10 @@ Application Description Incompressible LES solver for flow in a channel. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved \*---------------------------------------------------------------------------*/ @@ -60,38 +64,41 @@ int main(int argc, char *argv[]) sgsModel->correct(); - fvVectorMatrix UEqn + // Convection-diffusion matrix + fvVectorMatrix HUEqn ( - fvm::ddt(U) - + fvm::div(phi, U) + fvm::div(phi, U) + sgsModel->divDevBeff(U) == flowDirection*gradP ); + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + if (momentumPredictor) { - solve(UEqn == -fvc::grad(p)); + solve(ddtUEqn + HUEqn == -fvc::grad(p)); } + // Prepare clean Ap without time derivative contribution + // HJ, 26/Oct/2015 + volScalarField aU = HUEqn.A(); // --- PISO loop - volScalarField rUA = 1.0/UEqn.A(); - for (int corr = 0; corr < nCorr; corr++) { - U = rUA*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, U, phi); + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); adjustPhi(phi, U, p); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( - fvm::laplacian(rUA, p) == fvc::div(phi) + fvm::laplacian(1/aU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); @@ -111,13 +118,17 @@ int main(int argc, char *argv[]) } } - #include "continuityErrs.H" +# include "continuityErrs.H" - U -= rUA*fvc::grad(p); + // Note: cannot call H(U) here because the velocity is not complete + // HJ, 22/Jan/2016 + U = 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ); U.correctBoundaryConditions(); } - // Correct driving force for a constant mass flow rate // Extract the velocity in the flow direction @@ -127,9 +138,9 @@ int main(int argc, char *argv[]) // Calculate the pressure gradient increment needed to // adjust the average flow-rate to the correct value dimensionedScalar gragPplus = - (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V()); + (magUbar - magUbarStar)*aU.weightedAverage(mesh.V()); - U += flowDirection*rUA*gragPplus; + U += gragPplus/aU*flowDirection; gradP += gragPplus; diff --git a/tutorials/incompressible/channelFoam/channel395/system/fvSchemes b/tutorials/incompressible/channelFoam/channel395/system/fvSchemes index b8b2e5c55..1e1237665 100644 --- a/tutorials/incompressible/channelFoam/channel395/system/fvSchemes +++ b/tutorials/incompressible/channelFoam/channel395/system/fvSchemes @@ -29,7 +29,7 @@ gradSchemes divSchemes { default none; - div(phi,U) Gauss linear; + div(phi,U) Gauss linearUpwind Gauss linear; div(phi,k) Gauss limitedLinear 1; div(phi,B) Gauss limitedLinear 1; div(B) Gauss linear; From 74b90dd250a408e44bde52a1e394da0f37011805 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 14:06:15 +0200 Subject: [PATCH 4/7] incompressible/pimpleDyMFoam, solver + tutorials --- .../incompressible/pimpleDyMFoam/UEqn.H | 32 +++++---- .../incompressible/pimpleDyMFoam/correctPhi.H | 2 +- .../pimpleDyMFoam/createFields.H | 8 +-- .../incompressible/pimpleDyMFoam/pEqn.H | 70 +++++++++++++++++++ .../pimpleDyMFoam/pimpleDyMFoam.C | 61 ++-------------- .../pimpleDyMFoam/movingCone/system/fvSchemes | 6 +- .../movingCone/system/fvSolution | 6 ++ .../movingCylinders/system/fvSchemes | 5 +- 8 files changed, 112 insertions(+), 78 deletions(-) create mode 100644 applications/solvers/incompressible/pimpleDyMFoam/pEqn.H diff --git a/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H index 6aa608d3a..0fe0a6284 100644 --- a/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H +++ b/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H @@ -1,27 +1,33 @@ - fvVectorMatrix UEqn + // Convection-diffusion matrix + fvVectorMatrix HUEqn ( - fvm::ddt(U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevReff(U) ); + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + // Get under-relaxation factor + scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); + if (oCorr == nOuterCorr - 1) { if (mesh.solutionDict().relax("UFinal")) { - UEqn.relax(mesh.solutionDict().relaxationFactor("UFinal")); + UUrf = mesh.solutionDict().relaxationFactor("UFinal"); } else { - UEqn.relax(1); + UUrf = 1; } } - else - { - UEqn.relax(); - } - if (momentumPredictor) - { - solve(UEqn == -fvc::grad(p)); - } + // Solve momentum predictor + solve + ( + ddtUEqn + + relax(HUEqn, UUrf) + == + - fvc::grad(p) + ); diff --git a/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H index ce45e1bd3..60b36ceef 100644 --- a/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H +++ b/applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H @@ -39,7 +39,7 @@ { fvScalarMatrix pcorrEqn ( - fvm::laplacian(rAU, pcorr) == fvc::div(phi) + fvm::laplacian(1/aU, pcorr) == fvc::div(phi) ); pcorrEqn.setReference(pRefCell, pRefValue); diff --git a/applications/solvers/incompressible/pimpleDyMFoam/createFields.H b/applications/solvers/incompressible/pimpleDyMFoam/createFields.H index bc798ba98..1c04f5dbe 100644 --- a/applications/solvers/incompressible/pimpleDyMFoam/createFields.H +++ b/applications/solvers/incompressible/pimpleDyMFoam/createFields.H @@ -41,18 +41,18 @@ incompressible::turbulenceModel::New(U, phi, laminarTransport) ); - Info<< "Reading field rAU if present\n" << endl; - volScalarField rAU + Info<< "Reading field aU if present\n" << endl; + volScalarField aU ( IOobject ( - "rAU", + "aU", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), mesh, - runTime.deltaT(), + 1/runTime.deltaT(), zeroGradientFvPatchScalarField::typeName ); diff --git a/applications/solvers/incompressible/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleDyMFoam/pEqn.H new file mode 100644 index 000000000..b7b006fde --- /dev/null +++ b/applications/solvers/incompressible/pimpleDyMFoam/pEqn.H @@ -0,0 +1,70 @@ +{ + p.boundaryField().updateCoeffs(); + + // Prepare clean Ap without time derivative contribution and + // without contribution from under-relaxation + // HJ, 26/Oct/2015 + aU = HUEqn.A(); + + // Store velocity under-relaxation point before using U for + // the flux precursor + U.storePrevIter(); + + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); + + // Global flux balance + adjustPhi(phi, U, p); + + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + + if + ( + corr == nCorr - 1 + && nonOrth == nNonOrthCorr + ) + { + pEqn.solve + ( + mesh.solutionDict().solver(p.name() + "Final") + ); + } + else + { + pEqn.solve(mesh.solutionDict().solver(p.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + } + + // Explicitly relax pressure for momentum corrector + if (oCorr != nOuterCorr - 1) + { + p.relax(); + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "movingMeshContinuityErrs.H" + + U = UUrf* + ( + 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ) + ) + + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C index 6c24eed98..f3ee49d3e 100644 --- a/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C +++ b/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C @@ -30,6 +30,11 @@ Description Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" @@ -106,61 +111,7 @@ int main(int argc, char *argv[]) // --- PISO loop for (int corr = 0; corr < nCorr; corr++) { - rAU = 1.0/UEqn.A(); - - U = rAU*UEqn.H(); - phi = (fvc::interpolate(U) & mesh.Sf()); - // ddtPhiCorr does not work. HJ, 20/Nov/2013 - - adjustPhi(phi, U, p); - - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(rAU, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - - if - ( -// oCorr == nOuterCorr - 1 - corr == nCorr - 1 - && nonOrth == nNonOrthCorr - ) - { - pEqn.solve - ( - mesh.solutionDict().solver(p.name() + "Final") - ); - } - else - { - pEqn.solve(mesh.solutionDict().solver(p.name())); - } - - if (nonOrth == nNonOrthCorr) - { - phi -= pEqn.flux(); - } - } - -# include "continuityErrs.H" - - // Explicitly relax pressure for momentum corrector - if (oCorr != nOuterCorr - 1) - { - p.relax(); - } - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); - -# include "movingMeshContinuityErrs.H" - - U -= rAU*fvc::grad(p); - U.correctBoundaryConditions(); +# include "pEqn.H" } turbulence->correct(); diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes index b55d4e0a1..650771f63 100644 --- a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes +++ b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes @@ -29,7 +29,7 @@ gradSchemes divSchemes { default none; - div(phi,U) Gauss linear; + div(phi,U) Gauss linearUpwind Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; } @@ -37,8 +37,8 @@ laplacianSchemes { default none; laplacian(nu,U) Gauss linear corrected; - laplacian(rAU,pcorr) Gauss linear corrected; - laplacian(rAU,p) Gauss linear corrected; + laplacian((1|aU),pcorr) Gauss linear corrected; + laplacian((1|aU),p) Gauss linear corrected; laplacian(diffusivity,cellMotionU) Gauss linear uncorrected; laplacian(nuEff,U) Gauss linear uncorrected; } diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution index e610ff1b1..3b88c2f9b 100644 --- a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution +++ b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution @@ -65,4 +65,10 @@ PIMPLE nNonOrthogonalCorrectors 0; } +relaxationFactors +{ + U 1; + UFinal 1; +} + // ************************************************************************* // diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCylinders/system/fvSchemes b/tutorials/incompressible/pimpleDyMFoam/movingCylinders/system/fvSchemes index 0ac0eb72a..97781f0b3 100644 --- a/tutorials/incompressible/pimpleDyMFoam/movingCylinders/system/fvSchemes +++ b/tutorials/incompressible/pimpleDyMFoam/movingCylinders/system/fvSchemes @@ -41,14 +41,15 @@ laplacianSchemes laplacian(1,p) Gauss linear limited 0.5; - laplacian((1|A(U)),p) Gauss linear limited 0.5; + laplacian((1|aU),pcorr) Gauss linear limited 0.5; + laplacian((1|aU),p) Gauss linear limited 0.5; } interpolationSchemes { default linear; interpolate(HbyA) linear; - interpolate(1|A) linear; + interpolate(1|aU) linear; } snGradSchemes From 3a0021239de4ab3ac291e7b015509e55f849f071 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 15:08:06 +0200 Subject: [PATCH 5/7] incompressible/pimpleFoam, solver + tutorials --- .../solvers/incompressible/pimpleFoam/UEqn.H | 41 ++++---- .../incompressible/pimpleFoam/createFields.H | 16 ++++ .../solvers/incompressible/pimpleFoam/pEqn.H | 96 +++++++++++-------- .../incompressible/pimpleFoam/pimpleFoam.C | 5 + .../pimpleFoam/t-junction/system/fvSchemes | 2 +- 5 files changed, 96 insertions(+), 64 deletions(-) diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H index 64ffbea02..f6c83eee3 100644 --- a/applications/solvers/incompressible/pimpleFoam/UEqn.H +++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H @@ -1,43 +1,40 @@ -// Solve the Momentum equation - -tmp UEqn +// Convection-diffusion matrx +fvVectorMatrix HUEqn ( - fvm::ddt(U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevReff(U) ); +// Time derivative matrix +fvVectorMatrix ddtUEqn(fvm::ddt(U)); + +// Get under-relaxation factor +scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); + if (oCorr == nOuterCorr - 1) { if (mesh.solutionDict().relax("UFinal")) { - UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal")); + UUrf = mesh.solutionDict().relaxationFactor("UFinal"); } else { - UEqn().relax(1); + UUrf = 1; } } -else -{ - UEqn().relax(); -} - -volScalarField rUA = 1.0/UEqn().A(); if (momentumPredictor) { - if (oCorr == nOuterCorr-1) - { - solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal")); - } - else - { - solve(UEqn() == -fvc::grad(p)); - } + solve + ( + ddtUEqn + + relax(HUEqn, UUrf) + == + - fvc::grad(p) + ); } else { - U = rUA*(UEqn().H() - fvc::grad(p)); + U = (ddtUEqn.H() + HUEqn.H() - fvc::grad(p))/(HUEqn.A() + ddtUEqn.A()); U.correctBoundaryConditions(); } diff --git a/applications/solvers/incompressible/pimpleFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/createFields.H index e4127150c..7a8fdfb8a 100644 --- a/applications/solvers/incompressible/pimpleFoam/createFields.H +++ b/applications/solvers/incompressible/pimpleFoam/createFields.H @@ -40,3 +40,19 @@ autoPtr turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); + +Info<< "Reading field aU if present\n" << endl; +volScalarField aU +( + IOobject + ( + "aU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + 1/runTime.deltaT(), + zeroGradientFvPatchScalarField::typeName +); diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H index 3a750deb8..72143f6a9 100644 --- a/applications/solvers/incompressible/pimpleFoam/pEqn.H +++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H @@ -1,53 +1,67 @@ -U = rUA*UEqn().H(); - -if (nCorr <= 1) { - UEqn.clear(); -} + p.boundaryField().updateCoeffs(); -phi = (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, U, phi); + // Prepare clean Ap without time derivative contribution and + // without contribution from under-relaxation + // HJ, 26/Oct/2015 + aU = HUEqn.A(); -adjustPhi(phi, U, p); + // Store velocity under-relaxation point before using U for the flux + // precursor + U.storePrevIter(); -// Non-orthogonal pressure corrector loop -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) -{ - // Pressure corrector - fvScalarMatrix pEqn - ( - fvm::laplacian(rUA, p) == fvc::div(phi) - ); + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); - pEqn.setReference(pRefCell, pRefValue); + // Global flux balance + adjustPhi(phi, U, p); - if - ( - oCorr == nOuterCorr - 1 - && corr == nCorr - 1 - && nonOrth == nNonOrthCorr - ) + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { - pEqn.solve(mesh.solutionDict().solver("pFinal")); - } - else - { - pEqn.solve(); + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + + if + ( + corr == nCorr - 1 + && nonOrth == nNonOrthCorr + ) + { + pEqn.solve + ( + mesh.solutionDict().solver(p.name() + "Final") + ); + } + else + { + pEqn.solve(mesh.solutionDict().solver(p.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } } - if (nonOrth == nNonOrthCorr) + // Explicitly relax pressure for momentum corrector + if (oCorr != nOuterCorr - 1) { - phi -= pEqn.flux(); + p.relax(); } + +# include "movingMeshContinuityErrs.H" + + U = UUrf* + ( + 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ) + ) + + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); } - -#include "continuityErrs.H" - -// Explicitly relax pressure for momentum corrector except for last corrector -if (oCorr != nOuterCorr-1) -{ - p.relax(); -} - -U -= rUA*fvc::grad(p); -U.correctBoundaryConditions(); diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C index 30d03c55d..aadcc861c 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C +++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C @@ -30,6 +30,11 @@ Description Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" diff --git a/tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes b/tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes index ca497019e..b81a47b4e 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes +++ b/tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes @@ -43,7 +43,7 @@ laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear corrected; - laplacian((1|A(U)),p) Gauss linear corrected; + laplacian((1|aU),p) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; From fb2d756eba4e343017bbdcff278f90a1cff81175 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 7 Apr 2016 15:39:09 +0200 Subject: [PATCH 6/7] incompressible/simpleFoam, solver + tutorials --- .../solvers/incompressible/simpleFoam/UEqn.H | 10 ++++--- .../solvers/incompressible/simpleFoam/pEqn.H | 27 ++++++++++++++----- .../incompressible/simpleFoam/simpleFoam.C | 4 +++ .../simpleFoam/motorBike/system/fvSchemes | 7 +---- .../simpleFoam/motorBike/system/fvSolution | 2 +- 5 files changed, 34 insertions(+), 16 deletions(-) diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H index 62c9f5536..a006813c4 100644 --- a/applications/solvers/incompressible/simpleFoam/UEqn.H +++ b/applications/solvers/incompressible/simpleFoam/UEqn.H @@ -1,16 +1,20 @@ // Solve the momentum equation - tmp UEqn + tmp HUEqn ( fvm::div(phi, U) + turbulence->divDevReff(U) ); - UEqn().relax(); + // Get under-relaxation factor + const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); + // Momentum solution eqnResidual = solve ( - UEqn() == -fvc::grad(p) + relax(HUEqn(), UUrf) + == + -fvc::grad(p) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H index 42ebfc225..021350f09 100644 --- a/applications/solvers/incompressible/simpleFoam/pEqn.H +++ b/applications/solvers/incompressible/simpleFoam/pEqn.H @@ -1,17 +1,30 @@ p.boundaryField().updateCoeffs(); - volScalarField AU = UEqn().A(); - U = UEqn().H()/AU; - UEqn.clear(); + // Prepare clean 1/Ap without contribution from under-relaxation + // HJ, 26/Oct/2015 + volScalarField rUA + ( + "(1|A(U))", + 1/HUEqn().A() + ); + + // Store velocity under-relaxation point before using U for + // the flux precursor + U.storePrevIter(); + + U = rUA*HUEqn().H(); + HUEqn.clear(); phi = fvc::interpolate(U) & mesh.Sf(); + + // Global flux balance adjustPhi(phi, U, p); // Non-orthogonal pressure corrector loop - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( - fvm::laplacian(1.0/AU, p) == fvc::div(phi) + fvm::laplacian(rUA, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); @@ -39,6 +52,8 @@ p.relax(); // Momentum corrector - U -= fvc::grad(p)/AU; + // Note: since under-relaxation does not change aU, H/a in U can be + // re-used. HJ, 22/Jan/2016 + U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter(); U.correctBoundaryConditions(); diff --git a/applications/solvers/incompressible/simpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/simpleFoam.C index 5e80c10f7..5a6e7ef7e 100644 --- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C +++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C @@ -26,6 +26,10 @@ Application Description Steady-state solver for incompressible, turbulent flow + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes index c9329c21e..514e3240f 100644 --- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes +++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes @@ -21,10 +21,7 @@ ddtSchemes gradSchemes { - default Gauss linear; - grad(p) Gauss linear; - grad(U) Gauss linear; -// grad(U) cellLimited Gauss linear 1; + default cellLimited leastSquares 1; } divSchemes @@ -39,8 +36,6 @@ divSchemes laplacianSchemes { default Gauss linear corrected; -// default Gauss linear limited 0.5; -// default Gauss linear limited 0.333; } interpolationSchemes diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution index 61742eba9..928d503a6 100644 --- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution +++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution @@ -60,7 +60,7 @@ solvers SIMPLE { - nNonOrthogonalCorrectors 0; + nNonOrthogonalCorrectors 2; } relaxationFactors From 59b857165abbcb6b3923ab9ad3cb78423c34a61c Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Fri, 8 Apr 2016 15:19:19 +0200 Subject: [PATCH 7/7] incompressible/simpleSRFFoam, solver + tutorials --- .../incompressible/simpleSRFFoam/UEqn.H | 14 +++++ .../simpleSRFFoam/createFields.H | 1 - .../incompressible/simpleSRFFoam/pEqn.H | 48 +++++++++++++++++ .../simpleSRFFoam/simpleSRFFoam.C | 54 +++---------------- .../simpleSRFFoam/mixer/system/fvSolution | 6 +-- 5 files changed, 72 insertions(+), 51 deletions(-) create mode 100644 applications/solvers/incompressible/simpleSRFFoam/UEqn.H create mode 100644 applications/solvers/incompressible/simpleSRFFoam/pEqn.H diff --git a/applications/solvers/incompressible/simpleSRFFoam/UEqn.H b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H new file mode 100644 index 000000000..070e768d8 --- /dev/null +++ b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H @@ -0,0 +1,14 @@ + // Solve the momentum equation + + tmp HUrelEqn + ( + fvm::div(phi, Urel) + + turbulence->divDevReff(Urel) + + SRF->Su() + ); + + // Get under-relaxation factor + const scalar UUrf = mesh.solutionDict().relaxationFactor(Urel.name()); + + // Momentum solution + solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p)); diff --git a/applications/solvers/incompressible/simpleSRFFoam/createFields.H b/applications/solvers/incompressible/simpleSRFFoam/createFields.H index 864b1dd93..681f75cb9 100644 --- a/applications/solvers/incompressible/simpleSRFFoam/createFields.H +++ b/applications/solvers/incompressible/simpleSRFFoam/createFields.H @@ -71,4 +71,3 @@ ), Urel + SRF->U() ); - diff --git a/applications/solvers/incompressible/simpleSRFFoam/pEqn.H b/applications/solvers/incompressible/simpleSRFFoam/pEqn.H new file mode 100644 index 000000000..04e6813d2 --- /dev/null +++ b/applications/solvers/incompressible/simpleSRFFoam/pEqn.H @@ -0,0 +1,48 @@ + p.boundaryField().updateCoeffs(); + + // Prepare clean 1/Ap without contribution from under-relaxation + // HJ, 26/Oct/2015 + volScalarField rUA + ( + "(1|A(Urel))", + 1/HUrelEqn().A() + ); + + // Store velocity under-relaxation point before using U for the flux + // precursor + Urel.storePrevIter(); + + Urel = rUA*HUrelEqn().H(); + HUrelEqn.clear(); + phi = fvc::interpolate(Urel) & mesh.Sf(); + + // Global flux balance + adjustPhi(phi, Urel, p); + + // Non-orthogonal pressure corrector loop + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rUA, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi -= pEqn.flux(); + } + } + +# include "continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + // Momentum corrector + // Note: since under-relaxation does not change aU, H/a in U can be + // re-used. HJ, 22/Jan/2016 + Urel = UUrf*(Urel - rUA*fvc::grad(p)) + (1 - UUrf)*Urel.prevIter(); + Urel.correctBoundaryConditions(); diff --git a/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C b/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C index 528c913b1..4e844b1d1 100644 --- a/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C +++ b/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C @@ -25,8 +25,12 @@ Application simpleSRFFoam Description - Steady-state solver for incompressible, turbulent flow of non-Newtonian + Steady-state solver for incompressible, turbulent flow of Newtonian fluids with single rotating frame. + Consistent formulation without time-step and relaxation dependence by Jasak + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved \*---------------------------------------------------------------------------*/ @@ -47,8 +51,6 @@ int main(int argc, char *argv[]) # include "createFields.H" # include "initContinuityErrs.H" - //mesh.clearPrimitives(); - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -63,50 +65,8 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { - // Momentum predictor - tmp UrelEqn - ( - fvm::div(phi, Urel) - + turbulence->divDevReff(Urel) - + SRF->Su() - ); - - UrelEqn().relax(); - - solve(UrelEqn() == -fvc::grad(p)); - - p.boundaryField().updateCoeffs(); - volScalarField AUrel = UrelEqn().A(); - Urel = UrelEqn().H()/AUrel; - UrelEqn.clear(); - phi = fvc::interpolate(Urel) & mesh.Sf(); - adjustPhi(phi, Urel, p); - - // Non-orthogonal pressure corrector loop - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(1.0/AUrel, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - pEqn.solve(); - - if (nonOrth == nNonOrthCorr) - { - phi -= pEqn.flux(); - } - } - -# include "continuityErrs.H" - - // Explicitly relax pressure for momentum corrector - p.relax(); - - // Momentum corrector - Urel -= fvc::grad(p)/AUrel; - Urel.correctBoundaryConditions(); +# include "UEqn.H" +# include "pEqn.H" } turbulence->correct(); diff --git a/tutorials/incompressible/simpleSRFFoam/mixer/system/fvSolution b/tutorials/incompressible/simpleSRFFoam/mixer/system/fvSolution index f60765a81..0a7837c20 100644 --- a/tutorials/incompressible/simpleSRFFoam/mixer/system/fvSolution +++ b/tutorials/incompressible/simpleSRFFoam/mixer/system/fvSolution @@ -25,8 +25,8 @@ solvers }; Urel { - solver PBiCG; - preconditioner DILU; + solver smoothSolver; + smoother GaussSeidel; tolerance 1e-05; relTol 0.1; }; @@ -75,7 +75,7 @@ SIMPLE relaxationFactors { p 0.3; - Urel 0.7; + Urel 0.5; k 0.7; epsilon 0.7; omega 0.7;