applications/utilities/postProcessing/velocityField/Mach/Mach.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     Mach
  27
  28 Description
  29     Calculates and optionally writes the local Mach number from the velocity
  30     field U at each time.
  31
  32     The -nowrite option just outputs the max value without writing the field.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "calc.H"
  37 #include "basicPsiThermo.H"
  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     Info<< "\nEnd\n" << endl;
 140 }
 141
 142
 143 // ************************************************************************* //