From 74b90dd250a408e44bde52a1e394da0f37011805 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 5 Apr 2016 14:06:15 +0200 Subject: [PATCH] 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