src/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "supersonicFreestreamFvPatchVectorField.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "fvPatchFieldMapper.H"
  30 #include "volFields.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  38
  39 supersonicFreestreamFvPatchVectorField::supersonicFreestreamFvPatchVectorField
  40 (
  41     const fvPatch& p,
  42     const DimensionedField<vector, volMesh>& iF
  43 )
  44 :
  45     mixedFvPatchVectorField(p, iF),
  46     UInf_(vector::zero),
  47     pInf_(0),
  48     TInf_(0),
  49     gamma_(0)
  50 {
  51     refValue() = patchInternalField();
  52     refGrad() = vector::zero;
  53     valueFraction() = 1;
  54 }
  55
  56
  57 supersonicFreestreamFvPatchVectorField::supersonicFreestreamFvPatchVectorField
  58 (
  59     const supersonicFreestreamFvPatchVectorField& ptf,
  60     const fvPatch& p,
  61     const DimensionedField<vector, volMesh>& iF,
  62     const fvPatchFieldMapper& mapper
  63 )
  64 :
  65     mixedFvPatchVectorField(ptf, p, iF, mapper),
  66     UInf_(ptf.UInf_),
  67     pInf_(ptf.pInf_),
  68     TInf_(ptf.TInf_),
  69     gamma_(ptf.gamma_)
  70 {}
  71
  72
  73 supersonicFreestreamFvPatchVectorField::supersonicFreestreamFvPatchVectorField
  74 (
  75     const fvPatch& p,
  76     const DimensionedField<vector, volMesh>& iF,
  77     const dictionary& dict
  78 )
  79 :
  80     mixedFvPatchVectorField(p, iF),
  81     UInf_(dict.lookup("UInf")),
  82     pInf_(readScalar(dict.lookup("pInf"))),
  83     TInf_(readScalar(dict.lookup("TInf"))),
  84     gamma_(readScalar(dict.lookup("gamma")))
  85 {
  86     if (dict.found("value"))
  87     {
  88         fvPatchField<vector>::operator=
  89         (
  90             vectorField("value", dict, p.size())
  91         );
  92     }
  93     else
  94     {
  95         fvPatchField<vector>::operator=(patchInternalField());
  96     }
  97
  98     refValue() = *this;
  99     refGrad() = vector::zero;
 100     valueFraction() = 1;
 101
 102     if (pInf_ < SMALL)
 103     {
 104         FatalIOErrorIn
 105         (
 106             "supersonicFreestreamFvPatchVectorField::"
 107             "supersonicFreestreamFvPatchVectorField"
 108             "(const fvPatch&, const vectorField&, const dictionary&)",
 109             dict
 110         )   << "    unphysical pInf specified (pInf <= 0.0)"
 111             << "\n    on patch " << this->patch().name()
 112             << " of field " << this->dimensionedInternalField().name()
 113             << " in file " << this->dimensionedInternalField().objectPath()
 114             << exit(FatalIOError);
 115     }
 116 }
 117
 118
 119 supersonicFreestreamFvPatchVectorField::supersonicFreestreamFvPatchVectorField
 120 (
 121     const supersonicFreestreamFvPatchVectorField& sfspvf
 122 )
 123 :
 124     mixedFvPatchVectorField(sfspvf),
 125     UInf_(sfspvf.UInf_),
 126     pInf_(sfspvf.pInf_),
 127     TInf_(sfspvf.TInf_),
 128     gamma_(sfspvf.gamma_)
 129 {}
 130
 131
 132 supersonicFreestreamFvPatchVectorField::supersonicFreestreamFvPatchVectorField
 133 (
 134     const supersonicFreestreamFvPatchVectorField& sfspvf,
 135     const DimensionedField<vector, volMesh>& iF
 136 )
 137 :
 138     mixedFvPatchVectorField(sfspvf, iF),
 139     UInf_(sfspvf.UInf_),
 140     pInf_(sfspvf.pInf_),
 141     TInf_(sfspvf.TInf_),
 142     gamma_(sfspvf.gamma_)
 143 {}
 144
 145
 146 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 147
 148 void supersonicFreestreamFvPatchVectorField::updateCoeffs()
 149 {
 150     if (!size() || updated())
 151     {
 152         return;
 153     }
 154
 155     const fvPatchField<scalar>& pT =
 156         patch().lookupPatchField<volScalarField, scalar>("T");
 157
 158     const fvPatchField<scalar>& pp =
 159         patch().lookupPatchField<volScalarField, scalar>("p");
 160
 161     const fvPatchField<scalar>& ppsi =
 162         patch().lookupPatchField<volScalarField, scalar>("psi");
 163
 164     // Need R of the free-stream flow.  Assume R is independent of location
 165     // along patch so use face 0
 166     scalar R = 1.0/(ppsi[0]*pT[0]);
 167
 168     scalar MachInf = mag(UInf_)/sqrt(gamma_*R*TInf_);
 169
 170     if (MachInf < 1.0)
 171     {
 172         FatalErrorIn
 173         (
 174             "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
 175         )   << "    MachInf < 1.0, free stream must be supersonic"
 176             << "\n    on patch " << this->patch().name()
 177             << " of field " << this->dimensionedInternalField().name()
 178             << " in file " << this->dimensionedInternalField().objectPath()
 179             << exit(FatalError);
 180     }
 181
 182     vectorField& Up = refValue();
 183     valueFraction() = 1;
 184
 185     // get the near patch internal cell values
 186     vectorField U = patchInternalField();
 187
 188
 189     // Find the component of U normal to the free-stream flow and in the
 190     // plane of the free-stream flow and the patch normal
 191
 192     // Direction of the free-stream flow
 193     vector UInfHat = UInf_/mag(UInf_);
 194
 195     // Normal to the plane defined by the free-stream and the patch normal
 196     vectorField nnInfHat = UInfHat ^ patch().nf();
 197
 198     //vectorField UnInf = U - nnInfHat*(nnInfHat & U);
 199     //vectorField Un = UnInf - UInfHat*(UInfHat & UnInf);
 200     //vectorField nHatInf =
 201     //    (Un/(mag(Un) + SMALL)) * sign(patch().nf() & Un);
 202
 203     // Normal to the free-stream in the plane defined by the free-stream
 204     // and the patch normal
 205     vectorField nHatInf = nnInfHat ^ UInfHat;
 206
 207     // Component of U normal to the free-stream in the plane defined by the
 208     // free-stream and the patch normal
 209     vectorField Un = nHatInf*(nHatInf & U);
 210
 211     // The tangential component is
 212     vectorField Ut = U - Un;
 213
 214     // Calculate the Prandtl-Meyer function of the free-stream
 215     scalar nuMachInf =
 216         sqrt((gamma_ + 1)/(gamma_ - 1))
 217        *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(MachInf) - 1)))
 218       - atan(sqr(MachInf) - 1);
 219
 220
 221     // Set the patch boundary condition based on the Mach number and direction
 222     // of the flow dictated by the boundary/free-stream pressure difference
 223
 224     forAll(Up, facei)
 225     {
 226         if (pp[facei] >= pInf_) // If outflow
 227         {
 228             // Assume supersonic outflow and calculate the boundary velocity
 229             // according to ???
 230
 231             scalar fpp =
 232                 sqrt(sqr(MachInf) - 1)
 233                /(gamma_*sqr(MachInf))*mag(Ut[facei])*log(pp[facei]/pInf_);
 234
 235             Up[facei] = Ut[facei] + fpp*nHatInf[facei];
 236
 237             // Calculate the Mach number of the boundary velocity
 238             scalar Mach = mag(Up[facei])/sqrt(gamma_/ppsi[facei]);
 239
 240             if (Mach <= 1) // If subsonic
 241             {
 242                 // Zero-gradient subsonic outflow
 243
 244                 Up[facei] = U[facei];
 245                 valueFraction()[facei] = 0;
 246             }
 247         }
 248         else // if inflow
 249         {
 250             // Calculate the Mach number of the boundary velocity
 251             // from the boundary pressure assuming constant total pressure
 252             // expansion from the free-stream
 253             scalar Mach =
 254                 sqrt
 255                 (
 256                     (2/(gamma_ - 1))*(1 + ((gamma_ - 1)/2)*sqr(MachInf))
 257                    *pow(pp[facei]/pInf_, (1 - gamma_)/gamma_)
 258                   - 2/(gamma_ - 1)
 259                 );
 260
 261             if (Mach > 1) // If supersonic
 262             {
 263                 // Supersonic inflow is assumed to occur according to the
 264                 // Prandtl-Meyer expansion process
 265
 266                 scalar nuMachf =
 267                     sqrt((gamma_ + 1)/(gamma_ - 1))
 268                    *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(Mach) - 1)))
 269                   - atan(sqr(Mach) - 1);
 270
 271                 scalar fpp = (nuMachInf - nuMachf)*mag(Ut[facei]);
 272
 273                 Up[facei] = Ut[facei] + fpp*nHatInf[facei];
 274             }
 275             else // If subsonic
 276             {
 277                 FatalErrorIn
 278                 (
 279                     "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
 280                 )   << "unphysical subsonic inflow has been generated"
 281                     << "\n    on patch " << this->patch().name()
 282                     << " of field " << this->dimensionedInternalField().name()
 283                     << " in file "
 284                     << this->dimensionedInternalField().objectPath()
 285                     << exit(FatalError);
 286             }
 287         }
 288     }
 289
 290     mixedFvPatchVectorField::updateCoeffs();
 291 }
 292
 293
 294 void supersonicFreestreamFvPatchVectorField::write(Ostream& os) const
 295 {
 296     fvPatchVectorField::write(os);
 297     os.writeKeyword("UInf") << UInf_ << token::END_STATEMENT << nl;
 298     os.writeKeyword("pInf") << pInf_ << token::END_STATEMENT << nl;
 299     os.writeKeyword("TInf") << TInf_ << token::END_STATEMENT << nl;
 300     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
 301     writeEntry("value", os);
 302 }
 303
 304
 305 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 306
 307 makePatchTypeField
 308 (
 309     fvPatchVectorField,
 310     supersonicFreestreamFvPatchVectorField
 311 );
 312
 313 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 314
 315 } // End namespace Foam
 316
 317 // ************************************************************************* //