applications/utilities/postProcessing/velocityField/Mach/Mach.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     Mach
  26
  27 Description
  28     Calculates and optionally writes the local Mach number from the velocity
  29     field U at each time.
  30
  31     The -nowrite option just outputs the max value without writing the field.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "calc.H"
  36 #include "basicPsiThermo.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
  41 {
  42     bool writeResults = !args.optionFound("noWrite");
  43
  44     IOobject Uheader
  45     (
  46         "U",
  47         runTime.timeName(),
  48         mesh,
  49         IOobject::MUST_READ
  50     );
  51
  52     IOobject Theader
  53     (
  54         "T",
  55         runTime.timeName(),
  56         mesh,
  57         IOobject::MUST_READ
  58     );
  59
  60     // Check U and T exists
  61     if (Uheader.headerOk() && Theader.headerOk())
  62     {
  63         autoPtr<volScalarField> MachPtr;
  64
  65         volVectorField U(Uheader, mesh);
  66
  67         if (isFile(runTime.constantPath()/"thermophysicalProperties"))
  68         {
  69             // thermophysical Mach
  70             autoPtr<basicPsiThermo> thermo
  71             (
  72                 basicPsiThermo::New(mesh)
  73             );
  74
  75             volScalarField Cp = thermo->Cp();
  76             volScalarField Cv = thermo->Cv();
  77
  78             MachPtr.set
  79             (
  80                 new volScalarField
  81                 (
  82                     IOobject
  83                     (
  84                         "Ma",
  85                         runTime.timeName(),
  86                         mesh
  87                     ),
  88                     mag(U)/(sqrt((Cp/Cv)*(Cp - Cv)*thermo->T()))
  89                 )
  90             );
  91         }
  92         else
  93         {
  94             // thermodynamic Mach
  95             IOdictionary thermoProps
  96             (
  97                 IOobject
  98                 (
  99                     "thermodynamicProperties",
 100                     runTime.constant(),
 101                     mesh,
 102                     IOobject::MUST_READ,
 103                     IOobject::NO_WRITE
 104                 )
 105             );
 106
 107             dimensionedScalar R(thermoProps.lookup("R"));
 108             dimensionedScalar Cv(thermoProps.lookup("Cv"));
 109
 110             volScalarField T(Theader, mesh);
 111
 112             MachPtr.set
 113             (
 114                 new volScalarField
 115                 (
 116                     IOobject
 117                     (
 118                         "Ma",
 119                         runTime.timeName(),
 120                         mesh
 121                     ),
 122                     mag(U)/(sqrt(((Cv + R)/Cv)*R*T))
 123                 )
 124             );
 125         }
 126
 127         Info<< "Mach max : " << max(MachPtr()).value() << endl;
 128
 129         if (writeResults)
 130         {
 131             MachPtr().write();
 132         }
 133     }
 134     else
 135     {
 136         Info<< "    Missing U or T" << endl;
 137     }
 138 }
 139
 140
 141 // ************************************************************************* //