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