src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.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 "advectiveFvPatchField.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "fvPatchFieldMapper.H"
  30 #include "volFields.H"
  31 #include "EulerDdtScheme.H"
  32 #include "CrankNicholsonDdtScheme.H"
  33 #include "backwardDdtScheme.H"
  34 #include "steadyStateDdtScheme.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  42
  43 template<class Type>
  44 advectiveFvPatchField<Type>::advectiveFvPatchField
  45 (
  46     const fvPatch& p,
  47     const DimensionedField<Type, volMesh>& iF
  48 )
  49 :
  50     mixedFvPatchField<Type>(p, iF),
  51     phiName_("phi"),
  52     rhoName_("rho"),
  53     fieldInf_(pTraits<Type>::zero),
  54     lInf_(0.0),
  55     inletOutlet_(false),
  56     correctSupercritical_(false)
  57 {
  58     this->refValue() = pTraits<Type>::zero;
  59     this->refGrad() = pTraits<Type>::zero;
  60     this->valueFraction() = 0.0;
  61 }
  62
  63
  64 template<class Type>
  65 advectiveFvPatchField<Type>::advectiveFvPatchField
  66 (
  67     const fvPatch& p,
  68     const DimensionedField<Type, volMesh>& iF,
  69     const dictionary& dict
  70 )
  71 :
  72     mixedFvPatchField<Type>(p, iF),
  73     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
  74     rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
  75     fieldInf_(pTraits<Type>::zero),
  76     lInf_(0.0),
  77     inletOutlet_(dict.lookup("inletOutlet")),
  78     correctSupercritical_(dict.lookup("correctSupercritical"))
  79 {
  80     if (dict.found("value"))
  81     {
  82         fvPatchField<Type>::operator=
  83         (
  84             Field<Type>("value", dict, p.size())
  85         );
  86     }
  87     else
  88     {
  89         fvPatchField<Type>::operator=(this->patchInternalField());
  90     }
  91
  92     this->refValue() = *this;
  93     this->refGrad() = pTraits<Type>::zero;
  94     this->valueFraction() = 0.0;
  95
  96     if (dict.readIfPresent("lInf", lInf_))
  97     {
  98         dict.lookup("fieldInf") >> fieldInf_;
  99
 100         if (lInf_ < 0.0)
 101         {
 102             FatalIOErrorIn
 103             (
 104                 "advectiveFvPatchField<Type>::"
 105                 "advectiveFvPatchField"
 106                 "(const fvPatch&, const Field<Type>&, const dictionary&)",
 107                 dict
 108             )   << "unphysical lInf specified (lInf < 0)\n"
 109                 << "    on patch " << this->patch().name()
 110                 << " of field " << this->dimensionedInternalField().name()
 111                 << " in file " << this->dimensionedInternalField().objectPath()
 112                 << exit(FatalIOError);
 113         }
 114     }
 115 }
 116
 117
 118 template<class Type>
 119 advectiveFvPatchField<Type>::advectiveFvPatchField
 120 (
 121     const advectiveFvPatchField& ptf,
 122     const fvPatch& p,
 123     const DimensionedField<Type, volMesh>& iF,
 124     const fvPatchFieldMapper& mapper
 125 )
 126 :
 127     mixedFvPatchField<Type>(ptf, p, iF, mapper),
 128     phiName_(ptf.phiName_),
 129     rhoName_(ptf.rhoName_),
 130     fieldInf_(ptf.fieldInf_),
 131     lInf_(ptf.lInf_),
 132     inletOutlet_(ptf.inletOutlet_),
 133     correctSupercritical_(ptf.correctSupercritical_)
 134 {}
 135
 136
 137 template<class Type>
 138 advectiveFvPatchField<Type>::advectiveFvPatchField
 139 (
 140     const advectiveFvPatchField& ptpsf
 141 )
 142 :
 143     mixedFvPatchField<Type>(ptpsf),
 144     phiName_(ptpsf.phiName_),
 145     rhoName_(ptpsf.rhoName_),
 146     fieldInf_(ptpsf.fieldInf_),
 147     lInf_(ptpsf.lInf_),
 148     inletOutlet_(ptpsf.inletOutlet_),
 149     correctSupercritical_(ptpsf.correctSupercritical_)
 150 {}
 151
 152
 153 template<class Type>
 154 advectiveFvPatchField<Type>::advectiveFvPatchField
 155 (
 156     const advectiveFvPatchField& ptpsf,
 157     const DimensionedField<Type, volMesh>& iF
 158 )
 159 :
 160     mixedFvPatchField<Type>(ptpsf, iF),
 161     phiName_(ptpsf.phiName_),
 162     rhoName_(ptpsf.rhoName_),
 163     fieldInf_(ptpsf.fieldInf_),
 164     lInf_(ptpsf.lInf_),
 165     inletOutlet_(ptpsf.inletOutlet_),
 166     correctSupercritical_(ptpsf.correctSupercritical_)
 167 {}
 168
 169
 170 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 171
 172 template<class Type>
 173 tmp<scalarField> advectiveFvPatchField<Type>::advectionSpeed() const
 174 {
 175     const surfaceScalarField& phi =
 176         this->db().objectRegistry::lookupObject<surfaceScalarField>(phiName_);
 177
 178     fvsPatchField<scalar> phip = this->patch().lookupPatchField
 179     (
 180         phiName_,
 181         reinterpret_cast<const surfaceScalarField*>(0),
 182         reinterpret_cast<const scalar*>(0)
 183     );
 184
 185     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
 186     {
 187         const fvPatchScalarField& rhop = this->patch().lookupPatchField
 188         (
 189             rhoName_,
 190             reinterpret_cast<const volScalarField*>(0),
 191             reinterpret_cast<const scalar*>(0)
 192         );
 193
 194         return phip/(rhop*this->patch().magSf());
 195     }
 196     else
 197     {
 198         return phip/this->patch().magSf();
 199     }
 200 }
 201
 202
 203 template<class Type>
 204 tmp<scalarField> advectiveFvPatchField<Type>::supercritical() const
 205 {
 206     // In base class, the condition is never supercritical
 207     return tmp<scalarField>(new scalarField(this->size(), scalar(0)));
 208 }
 209
 210
 211 template<class Type>
 212 void advectiveFvPatchField<Type>::updateCoeffs()
 213 {
 214     if (this->updated())
 215     {
 216         return;
 217     }
 218
 219     const GeometricField<Type, fvPatchField, volMesh>& field =
 220         this->db().objectRegistry::
 221         lookupObject<GeometricField<Type, fvPatchField, volMesh> >
 222         (
 223             this->dimensionedInternalField().name()
 224         );
 225
 226     word ddtScheme
 227     (
 228         this->dimensionedInternalField().mesh().ddtScheme(field.name())
 229     );
 230
 231     // Calculate the advection speed of the field wave
 232     // If the wave is incoming set the speed to 0.
 233     scalarField w = Foam::max(advectionSpeed(), scalar(0));
 234
 235     // Get time-step
 236     scalarField deltaT(this->size(), this->db().time().deltaT().value());
 237
 238     // For steady-state, calculate formal deltaT from the
 239     // under-relaxation factor
 240     if (ddtScheme == fv::steadyStateDdtScheme<Type>::typeName)
 241     {
 242         // Get under-relaxation factor
 243 //         scalar urf =
 244 //             this->dimensionedInternalField().mesh().relaxationFactor
 245 //             (
 246 //                 field.name()
 247 //             );
 248
 249 //         deltaT = urf/(1 - urf)*1/(w*this->patch().deltaCoeffs() + SMALL);
 250
 251         // Set delta t for Co = 1
 252         deltaT = 1/(w*this->patch().deltaCoeffs() + SMALL);
 253     }
 254
 255     // Calculate the field wave coefficient alpha (See notes)
 256     scalarField alpha = w*deltaT*this->patch().deltaCoeffs();
 257
 258     label patchi = this->patch().index();
 259
 260     // Non-reflecting outflow boundary
 261     // If lInf_ defined setup relaxation to the value fieldInf_.
 262     // HJ, 25/Sep/2009
 263     if (lInf_ > SMALL)
 264     {
 265         // Calculate the field relaxation coefficient k (See notes)
 266         scalarField k = w*deltaT/lInf_;
 267
 268         if
 269         (
 270             ddtScheme == fv::EulerDdtScheme<Type>::typeName
 271          || ddtScheme == fv::CrankNicholsonDdtScheme<Type>::typeName
 272         )
 273         {
 274             this->refValue() =
 275             (
 276                 field.oldTime().boundaryField()[patchi] + k*fieldInf_
 277             )/(1.0 + k);
 278
 279             this->valueFraction() = (1.0 + k)/(1.0 + alpha + k);
 280         }
 281         else if (ddtScheme == fv::backwardDdtScheme<Type>::typeName)
 282         {
 283             this->refValue() =
 284             (
 285                 2.0*field.oldTime().boundaryField()[patchi]
 286               - 0.5*field.oldTime().oldTime().boundaryField()[patchi]
 287               + k*fieldInf_
 288             )/(1.5 + k);
 289
 290             this->valueFraction() = (1.5 + k)/(1.5 + alpha + k);
 291         }
 292         else if (ddtScheme == fv::steadyStateDdtScheme<Type>::typeName)
 293         {
 294             this->refValue() =
 295             (
 296                 field.prevIter().boundaryField()[patchi] + k*fieldInf_
 297             )/(1.0 + k);
 298
 299             this->valueFraction() = (1.0 + k)/(1.0 + alpha + k);
 300         }
 301         else
 302         {
 303             FatalErrorIn
 304             (
 305                 "advectiveFvPatchField<Type>::updateCoeffs()"
 306             )   << "    Unsupported temporal differencing scheme : "
 307                 << ddtScheme
 308                 << "\n    on patch " << this->patch().name()
 309                 << " of field " << this->dimensionedInternalField().name()
 310                 << " in file " << this->dimensionedInternalField().objectPath()
 311                 << exit(FatalError);
 312         }
 313     }
 314     else
 315     {
 316         if
 317         (
 318             ddtScheme == fv::EulerDdtScheme<Type>::typeName
 319          || ddtScheme == fv::CrankNicholsonDdtScheme<Type>::typeName
 320         )
 321         {
 322             this->refValue() = field.oldTime().boundaryField()[patchi];
 323
 324             this->valueFraction() = 1.0/(1.0 + alpha);
 325         }
 326         else if (ddtScheme == fv::backwardDdtScheme<Type>::typeName)
 327         {
 328             this->refValue() =
 329             (
 330                 2.0*field.oldTime().boundaryField()[patchi]
 331               - 0.5*field.oldTime().oldTime().boundaryField()[patchi]
 332             )/1.5;
 333
 334             this->valueFraction() = 1.5/(1.5 + alpha);
 335         }
 336         else if (ddtScheme == fv::steadyStateDdtScheme<Type>::typeName)
 337         {
 338             this->refValue() = field.prevIter().boundaryField()[patchi];
 339
 340             this->valueFraction() = 1.0/(1.0 + alpha);
 341         }
 342         else
 343         {
 344             FatalErrorIn
 345             (
 346                 "advectiveFvPatchField<Type>::updateCoeffs()"
 347             )   << "    Unsupported temporal differencing scheme : "
 348                 << ddtScheme
 349                 << "\n    on patch " << this->patch().name()
 350                 << " of field " << this->dimensionedInternalField().name()
 351                 << " in file " << this->dimensionedInternalField().objectPath()
 352                 << exit(FatalError);
 353         }
 354     }
 355
 356     // Get access to flux field
 357     fvsPatchField<scalar> phip = this->patch().lookupPatchField
 358     (
 359         phiName_,
 360         reinterpret_cast<const surfaceScalarField*>(NULL),
 361         reinterpret_cast<const scalar*>(NULL)
 362     );
 363
 364     // Treatment for supercritical inlet or outlet.  HJ, 28/Oct/2009
 365     // If the flow is faster than critical speed, the boundary condition
 366     // needs to behave as zero gradient at outflow and fixed value at inflow.
 367     // This is the case when eg.  a wave-transmissive pressure outlet
 368     // becomes supersonic.  Face selection for supercritical outlet is
 369     // done through a supercritical() virtual function: return 1 if
 370     // flow is supercritical and zero otherwise
 371     if (correctSupercritical_)
 372     {
 373         this->valueFraction() =
 374         (
 375             // If the flow is going out
 376             // - set zero gradient for supercritical; otherwise use current
 377             pos(phip)*
 378             (
 379                 (scalar(1) - this->supercritical())
 380               + this->supercritical()*this->valueFraction()
 381             )
 382             // If the flow is going in
 383             // - set fixed value for supercritical; otherwise use current
 384           + neg(phip)*
 385             (
 386                 this->supercritical()
 387               + (scalar(1) - this->supercritical())*this->valueFraction()
 388             )
 389         );
 390
 391         this->refValue() =
 392             // If the flow is going out
 393             // - for supercritical, the value does not matter.  use fieldInf_
 394             // - for subcritical, use current value
 395             pos(phip)*
 396             (
 397                 this->supercritical()*fieldInf_
 398               + (scalar(1) - this->supercritical())*this->refValue()
 399             )
 400             // If the flow is going in
 401             // - for supercritical use fieldInf_
 402           + neg(phip)*
 403             (
 404                 this->supercritical()*fieldInf_
 405               + (scalar(1) - this->supercritical())*this->refValue()
 406             );
 407     }
 408
 409     // Inlet-outlet treatment at the boundary.  HJ, 28/Oct/2009
 410     // If the flux is outgoing, the normal condition is used
 411     // if the flux is incoming, incoming value is set to fieldInf_
 412     // (value = fieldInf, valueFraction = 1)
 413     if (inletOutlet_)
 414     {
 415         this->refValue() = pos(phip)*this->refValue() + neg(phip)*fieldInf_;
 416         this->refGrad() = pTraits<Type>::zero;
 417         this->valueFraction() = pos(phip)*this->valueFraction() + neg(phip);
 418     }
 419
 420     mixedFvPatchField<Type>::updateCoeffs();
 421 }
 422
 423
 424 template<class Type>
 425 void advectiveFvPatchField<Type>::write(Ostream& os) const
 426 {
 427     fvPatchField<Type>::write(os);
 428
 429     if (phiName_ != "phi")
 430     {
 431         os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
 432     }
 433     if (rhoName_ != "rho")
 434     {
 435         os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
 436     }
 437
 438     if (lInf_ > SMALL)
 439     {
 440         os.writeKeyword("fieldInf") << fieldInf_
 441             << token::END_STATEMENT << nl;
 442         os.writeKeyword("lInf") << lInf_
 443         << token::END_STATEMENT << endl;
 444     }
 445
 446     os.writeKeyword("inletOutlet") << inletOutlet_
 447         << token::END_STATEMENT << nl;
 448
 449     os.writeKeyword("correctSupercritical") << correctSupercritical_
 450             << token::END_STATEMENT << nl;
 451
 452     this->writeEntry("value", os);
 453 }
 454
 455
 456 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 457
 458 } // End namespace Foam
 459
 460 // ************************************************************************* //