applications/solvers/incompressible/channelFoam/channelFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     channelFoam
  26
  27 Description
  28     Incompressible LES solver for flow in a channel.
  29
  30 Usage
  31     - channelFoam [OPTION]
  32
  33     @param -case \<dir\> \n
  34     Specify the case directory
  35
  36     @param -parallel \n
  37     Run the case in parallel
  38
  39     @param -help \n
  40     Display short usage message
  41
  42     @param -doc \n
  43     Display Doxygen documentation page
  44
  45     @param -srcDoc \n
  46     Display source code
  47
  48 \*---------------------------------------------------------------------------*/
  49
  50 #include <finiteVolume/fvCFD.H>
  51 #include <incompressibleTransportModels/singlePhaseTransportModel.H>
  52 #include <incompressibleLESModels/LESModel.H>
  53 #include <OpenFOAM/IFstream.H>
  54 #include <OpenFOAM/OFstream.H>
  55 #include <OpenFOAM/Random.H>
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59 int main(int argc, char *argv[])
  60 {
  61     #include <OpenFOAM/setRootCase.H>
  62     #include <OpenFOAM/createTime.H>
  63     #include <OpenFOAM/createMesh.H>
  64     #include "readTransportProperties.H"
  65     #include "createFields.H"
  66     #include <finiteVolume/initContinuityErrs.H>
  67     #include "createGradP.H"
  68
  69     Info<< "\nStarting time loop\n" << endl;
  70
  71     while (runTime.loop())
  72     {
  73         Info<< "Time = " << runTime.timeName() << nl << endl;
  74
  75         #include <finiteVolume/readPISOControls.H>
  76
  77         #include <finiteVolume/CourantNo.H>
  78
  79         sgsModel->correct();
  80
  81         fvVectorMatrix UEqn
  82         (
  83             fvm::ddt(U)
  84           + fvm::div(phi, U)
  85           + sgsModel->divDevBeff(U)
  86          ==
  87             flowDirection*gradP
  88         );
  89
  90         if (momentumPredictor)
  91         {
  92             solve(UEqn == -fvc::grad(p));
  93         }
  94
  95
  96         // --- PISO loop
  97
  98         volScalarField rUA = 1.0/UEqn.A();
  99
 100         for (int corr=0; corr<nCorr; corr++)
 101         {
 102             U = rUA*UEqn.H();
 103             phi = (fvc::interpolate(U) & mesh.Sf())
 104                 + fvc::ddtPhiCorr(rUA, U, phi);
 105
 106             adjustPhi(phi, U, p);
 107
 108             for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 109             {
 110                 fvScalarMatrix pEqn
 111                 (
 112                     fvm::laplacian(rUA, p) == fvc::div(phi)
 113                 );
 114
 115                 pEqn.setReference(pRefCell, pRefValue);
 116
 117                 if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
 118                 {
 119                     pEqn.solve(mesh.solver(p.name() + "Final"));
 120                 }
 121                 else
 122                 {
 123                     pEqn.solve(mesh.solver(p.name()));
 124                 }
 125
 126                 if (nonOrth == nNonOrthCorr)
 127                 {
 128                     phi -= pEqn.flux();
 129                 }
 130             }
 131
 132             #include <finiteVolume/continuityErrs.H>
 133
 134             U -= rUA*fvc::grad(p);
 135             U.correctBoundaryConditions();
 136         }
 137
 138
 139         // Correct driving force for a constant mass flow rate
 140
 141         // Extract the velocity in the flow direction
 142         dimensionedScalar magUbarStar =
 143             (flowDirection & U)().weightedAverage(mesh.V());
 144
 145         // Calculate the pressure gradient increment needed to
 146         // adjust the average flow-rate to the correct value
 147         dimensionedScalar gragPplus =
 148             (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());
 149
 150         U += flowDirection*rUA*gragPplus;
 151
 152         gradP += gragPplus;
 153
 154         Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
 155             << "pressure gradient = " << gradP.value() << endl;
 156
 157         runTime.write();
 158
 159         #include "writeGradP.H"
 160
 161         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 162             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 163             << nl << endl;
 164     }
 165
 166     Info<< "End\n" << endl;
 167
 168     return 0;
 169 }
 170
 171
 172 // ************************ vim: set sw=4 sts=4 et: ************************ //