applications/solvers/incompressible/pisoFoam/pisoFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Application
  26     pisoFoam
  27
  28 Description
  29     Transient solver for incompressible flow.
  30
  31     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "fvCFD.H"
  36 #include "singlePhaseTransportModel.H"
  37 #include "turbulenceModel.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 int main(int argc, char *argv[])
  42 {
  43     #include "setRootCase.H"
  44
  45     #include "createTime.H"
  46     #include "createMesh.H"
  47     #include "createFields.H"
  48     #include "initContinuityErrs.H"
  49
  50     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52     Info<< "\nStarting time loop\n" << endl;
  53
  54     while (runTime.loop())
  55     {
  56         Info<< "Time = " << runTime.timeName() << nl << endl;
  57
  58         #include "readPISOControls.H"
  59         #include "CourantNo.H"
  60
  61         // Pressure-velocity PISO corrector
  62         {
  63             // Momentum predictor
  64
  65             fvVectorMatrix UEqn
  66             (
  67                 fvm::ddt(U)
  68               + fvm::div(phi, U)
  69               + turbulence->divDevReff(U)
  70             );
  71
  72             UEqn.relax();
  73
  74             if (momentumPredictor)
  75             {
  76                 solve(UEqn == -fvc::grad(p));
  77             }
  78
  79             // --- PISO loop
  80
  81             for (int corr=0; corr<nCorr; corr++)
  82             {
  83                 volScalarField rUA = 1.0/UEqn.A();
  84
  85                 U = rUA*UEqn.H();
  86                 phi = (fvc::interpolate(U) & mesh.Sf())
  87                     + fvc::ddtPhiCorr(rUA, U, phi);
  88
  89                 adjustPhi(phi, U, p);
  90
  91                 // Non-orthogonal pressure corrector loop
  92                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  93                 {
  94                     // Pressure corrector
  95
  96                     fvScalarMatrix pEqn
  97                     (
  98                         fvm::laplacian(rUA, p) == fvc::div(phi)
  99                     );
 100
 101                     pEqn.setReference(pRefCell, pRefValue);
 102
 103                     if
 104                     (
 105                         corr == nCorr-1
 106                      && nonOrth == nNonOrthCorr
 107                     )
 108                     {
 109                         pEqn.solve(mesh.solver("pFinal"));
 110                     }
 111                     else
 112                     {
 113                         pEqn.solve();
 114                     }
 115
 116                     if (nonOrth == nNonOrthCorr)
 117                     {
 118                         phi -= pEqn.flux();
 119                     }
 120                 }
 121
 122                 #include "continuityErrs.H"
 123
 124                 U -= rUA*fvc::grad(p);
 125                 U.correctBoundaryConditions();
 126             }
 127         }
 128
 129         turbulence->correct();
 130
 131         runTime.write();
 132
 133         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 134             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 135             << nl << endl;
 136     }
 137
 138     Info<< "End\n" << endl;
 139
 140     return 0;
 141 }
 142
 143
 144 // ************************************************************************* //