applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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     icoDyMFoam
  27
  28 Description
  29     Transient solver for incompressible, laminar flow of Newtonian fluids
  30     with dynamic mesh.
  31
  32 Author
  33     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "fvCFD.H"
  38 #include "dynamicFvMesh.H"
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 int main(int argc, char *argv[])
  43 {
  44 #   include "setRootCase.H"
  45 #   include "createTime.H"
  46 #   include "createDynamicFvMesh.H"
  47 #   include "initContinuityErrs.H"
  48 #   include "initTotalVolume.H"
  49 #   include "createFields.H"
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53     Info<< "\nStarting time loop\n" << endl;
  54
  55     while (runTime.run())
  56     {
  57 #       include "readControls.H"
  58 #       include "checkTotalVolume.H"
  59
  60         // Make the fluxes absolute
  61         fvc::makeAbsolute(phi, U);
  62
  63 #       include "CourantNo.H"
  64 #       include "setDeltaT.H"
  65
  66         runTime++;
  67
  68         Info<< "Time = " << runTime.timeName() << nl << endl;
  69
  70         bool meshChanged = mesh.update();
  71         reduce(meshChanged, orOp<bool>());
  72
  73 #       include "volContinuity.H"
  74
  75         if (correctPhi && (mesh.moving() || meshChanged))
  76         {
  77             // Fluxes will be corrected to absolute velocity
  78             // HJ, 6/Feb/2009
  79 #           include "correctPhi.H"
  80         }
  81
  82         // Make the fluxes relative to the mesh motion
  83         fvc::makeRelative(phi, U);
  84
  85         if (mesh.moving() && checkMeshCourantNo)
  86         {
  87 #           include "meshCourantNo.H"
  88         }
  89
  90 #       include "UEqn.H"
  91
  92         // --- PISO loop
  93
  94         for (int corr = 0; corr < nCorr; corr++)
  95         {
  96             rAU = 1.0/UEqn.A();
  97
  98             U = rAU*UEqn.H();
  99             phi = (fvc::interpolate(U) & mesh.Sf());
 100               //+ fvc::ddtPhiCorr(rAU, U, phi);
 101
 102             adjustPhi(phi, U, p);
 103
 104             for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 105             {
 106                 fvScalarMatrix pEqn
 107                 (
 108                     fvm::laplacian(rAU, p) == fvc::div(phi)
 109                 );
 110
 111                 pEqn.setReference(pRefCell, pRefValue);
 112
 113                 if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
 114                 {
 115                     pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
 116                 }
 117                 else
 118                 {
 119                     pEqn.solve(mesh.solutionDict().solver(p.name()));
 120                 }
 121
 122                 if (nonOrth == nNonOrthCorr)
 123                 {
 124                     phi -= pEqn.flux();
 125                 }
 126             }
 127
 128 #           include "continuityErrs.H"
 129
 130             // Make the fluxes relative to the mesh motion
 131             fvc::makeRelative(phi, U);
 132
 133             U -= rAU*fvc::grad(p);
 134             U.correctBoundaryConditions();
 135         }
 136
 137         runTime.write();
 138
 139         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 140             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 141             << nl << endl;
 142     }
 143
 144     Info<< "End\n" << endl;
 145
 146     return(0);
 147 }
 148
 149
 150 // ************************************************************************* //