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 (C) 2011 OpenFOAM Foundation
   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     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 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
  40 {
  41     bool writeResults = !args.optionFound("noWrite");
  42
  43     IOobject Uheader
  44     (
  45         "U",
  46         runTime.timeName(),
  47         mesh,
  48         IOobject::MUST_READ
  49     );
  50
  51     IOobject Theader
  52     (
  53         "T",
  54         runTime.timeName(),
  55         mesh,
  56         IOobject::MUST_READ
  57     );
  58
  59     // Check U and T exists
  60     if (Uheader.headerOk() && Theader.headerOk())
  61     {
  62         autoPtr<volScalarField> MachPtr;
  63
  64         volVectorField U(Uheader, mesh);
  65
  66         if
  67         (
  68             IOobject
  69             (
  70                 "thermophysicalProperties",
  71                 runTime.constant(),
  72                 mesh
  73             ).headerOk()
  74         )
  75         {
  76             // thermophysical Mach
  77             autoPtr<basicPsiThermo> thermo
  78             (
  79                 basicPsiThermo::New(mesh)
  80             );
  81
  82             volScalarField Cp(thermo->Cp());
  83             volScalarField Cv(thermo->Cv());
  84
  85             MachPtr.set
  86             (
  87                 new volScalarField
  88                 (
  89                     IOobject
  90                     (
  91                         "Ma",
  92                         runTime.timeName(),
  93                         mesh
  94                     ),
  95                     mag(U)/(sqrt((Cp/Cv)*(Cp - Cv)*thermo->T()))
  96                 )
  97             );
  98         }
  99         else
 100         {
 101             // thermodynamic Mach
 102             IOdictionary thermoProps
 103             (
 104                 IOobject
 105                 (
 106                     "thermodynamicProperties",
 107                     runTime.constant(),
 108                     mesh,
 109                     IOobject::MUST_READ_IF_MODIFIED,
 110                     IOobject::NO_WRITE
 111                 )
 112             );
 113
 114             dimensionedScalar R(thermoProps.lookup("R"));
 115             dimensionedScalar Cv(thermoProps.lookup("Cv"));
 116
 117             volScalarField T(Theader, mesh);
 118
 119             MachPtr.set
 120             (
 121                 new volScalarField
 122                 (
 123                     IOobject
 124                     (
 125                         "Ma",
 126                         runTime.timeName(),
 127                         mesh
 128                     ),
 129                     mag(U)/(sqrt(((Cv + R)/Cv)*R*T))
 130                 )
 131             );
 132         }
 133
 134         Info<< "Mach max : " << max(MachPtr()).value() << endl;
 135
 136         if (writeResults)
 137         {
 138             MachPtr().write();
 139         }
 140     }
 141     else
 142     {
 143         Info<< "    Missing U or T" << endl;
 144     }
 145
 146     Info<< "\nEnd\n" << endl;
 147 }
 148
 149
 150 // ************************************************************************* //