src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "activePressureForceBaffleVelocityFvPatchVectorField.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "volFields.H"
  29 #include "surfaceFields.H"
  30 #include "cyclicFvPatch.H"
  31
  32 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  33
  34 Foam::activePressureForceBaffleVelocityFvPatchVectorField::
  35 activePressureForceBaffleVelocityFvPatchVectorField
  36 (
  37     const fvPatch& p,
  38     const DimensionedField<vector, volMesh>& iF
  39 )
  40 :
  41     fixedValueFvPatchVectorField(p, iF),
  42     pName_("p"),
  43     cyclicPatchName_(),
  44     cyclicPatchLabel_(-1),
  45     orientation_(1),
  46     initWallSf_(0),
  47     initCyclicSf_(0),
  48     nbrCyclicSf_(0),
  49     openFraction_(0),
  50     openingTime_(0),
  51     maxOpenFractionDelta_(0),
  52     curTimeIndex_(-1),
  53     minThresholdValue_(0),
  54     fBased_(1),
  55     baffleActivated_(0)
  56 {}
  57
  58
  59 Foam::activePressureForceBaffleVelocityFvPatchVectorField::
  60 activePressureForceBaffleVelocityFvPatchVectorField
  61 (
  62     const activePressureForceBaffleVelocityFvPatchVectorField& ptf,
  63     const fvPatch& p,
  64     const DimensionedField<vector, volMesh>& iF,
  65     const fvPatchFieldMapper& mapper
  66 )
  67 :
  68     fixedValueFvPatchVectorField(ptf, p, iF, mapper),
  69     pName_(ptf.pName_),
  70     cyclicPatchName_(ptf.cyclicPatchName_),
  71     cyclicPatchLabel_(ptf.cyclicPatchLabel_),
  72     orientation_(ptf.orientation_),
  73     initWallSf_(ptf.initWallSf_),
  74     initCyclicSf_(ptf.initCyclicSf_),
  75     nbrCyclicSf_(ptf.nbrCyclicSf_),
  76     openFraction_(ptf.openFraction_),
  77     openingTime_(ptf.openingTime_),
  78     maxOpenFractionDelta_(ptf.maxOpenFractionDelta_),
  79     curTimeIndex_(-1),
  80     minThresholdValue_(ptf.minThresholdValue_),
  81     fBased_(ptf.fBased_),
  82     baffleActivated_(ptf.baffleActivated_)
  83 {}
  84
  85
  86 Foam::activePressureForceBaffleVelocityFvPatchVectorField::
  87 activePressureForceBaffleVelocityFvPatchVectorField
  88 (
  89     const fvPatch& p,
  90     const DimensionedField<vector, volMesh>& iF,
  91     const dictionary& dict
  92 )
  93 :
  94     fixedValueFvPatchVectorField(p, iF),
  95     pName_("p"),
  96     cyclicPatchName_(dict.lookup("cyclicPatch")),
  97     cyclicPatchLabel_(p.patch().boundaryMesh().findPatchID(cyclicPatchName_)),
  98     orientation_(readLabel(dict.lookup("orientation"))),
  99     initWallSf_(0),
 100     initCyclicSf_(0),
 101     nbrCyclicSf_(0),
 102     openFraction_(readScalar(dict.lookup("openFraction"))),
 103     openingTime_(readScalar(dict.lookup("openingTime"))),
 104     maxOpenFractionDelta_(readScalar(dict.lookup("maxOpenFractionDelta"))),
 105     curTimeIndex_(-1),
 106     minThresholdValue_(readScalar(dict.lookup("minThresholdValue"))),
 107     fBased_(readBool(dict.lookup("forceBased"))),
 108     baffleActivated_(0)
 109 {
 110     fvPatchVectorField::operator=(vector::zero);
 111
 112     if (p.size() > 0)
 113     {
 114         initWallSf_ = p.Sf();
 115         initCyclicSf_ = p.boundaryMesh()[cyclicPatchLabel_].Sf();
 116         nbrCyclicSf_ =  refCast<const cyclicFvPatch>
 117         (
 118             p.boundaryMesh()[cyclicPatchLabel_]
 119         ).neighbFvPatch().Sf();
 120     }
 121
 122     if (dict.found("p"))
 123     {
 124         dict.lookup("p") >> pName_;
 125     }
 126 }
 127
 128
 129 Foam::activePressureForceBaffleVelocityFvPatchVectorField::
 130 activePressureForceBaffleVelocityFvPatchVectorField
 131 (
 132     const activePressureForceBaffleVelocityFvPatchVectorField& ptf
 133 )
 134 :
 135     fixedValueFvPatchVectorField(ptf),
 136     pName_(ptf.pName_),
 137     cyclicPatchName_(ptf.cyclicPatchName_),
 138     cyclicPatchLabel_(ptf.cyclicPatchLabel_),
 139     orientation_(ptf.orientation_),
 140     initWallSf_(ptf.initWallSf_),
 141     initCyclicSf_(ptf.initCyclicSf_),
 142     nbrCyclicSf_(ptf.nbrCyclicSf_),
 143     openFraction_(ptf.openFraction_),
 144     openingTime_(ptf.openingTime_),
 145     maxOpenFractionDelta_(ptf.maxOpenFractionDelta_),
 146     curTimeIndex_(-1),
 147     minThresholdValue_(ptf.minThresholdValue_),
 148     fBased_(ptf.fBased_),
 149     baffleActivated_(ptf.baffleActivated_)
 150 {}
 151
 152
 153 Foam::activePressureForceBaffleVelocityFvPatchVectorField::
 154 activePressureForceBaffleVelocityFvPatchVectorField
 155 (
 156     const activePressureForceBaffleVelocityFvPatchVectorField& ptf,
 157     const DimensionedField<vector, volMesh>& iF
 158 )
 159 :
 160     fixedValueFvPatchVectorField(ptf, iF),
 161     pName_(ptf.pName_),
 162     cyclicPatchName_(ptf.cyclicPatchName_),
 163     cyclicPatchLabel_(ptf.cyclicPatchLabel_),
 164     orientation_(ptf.orientation_),
 165     initWallSf_(ptf.initWallSf_),
 166     initCyclicSf_(ptf.initCyclicSf_),
 167     nbrCyclicSf_(ptf.nbrCyclicSf_),
 168     openFraction_(ptf.openFraction_),
 169     openingTime_(ptf.openingTime_),
 170     maxOpenFractionDelta_(ptf.maxOpenFractionDelta_),
 171     curTimeIndex_(-1),
 172     minThresholdValue_(ptf.minThresholdValue_),
 173     fBased_(ptf.fBased_),
 174     baffleActivated_(ptf.baffleActivated_)
 175 {}
 176
 177
 178 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 179
 180 void Foam::activePressureForceBaffleVelocityFvPatchVectorField::autoMap
 181 (
 182     const fvPatchFieldMapper& m
 183 )
 184 {
 185     fixedValueFvPatchVectorField::autoMap(m);
 186
 187     //- Note: cannot map field from cyclic patch anyway so just recalculate
 188     //  Areas should be consistent when doing autoMap except in case of
 189     //  topo changes.
 190     //- Note: we don't want to use Sf here since triggers rebuilding of
 191     //  fvMesh::S() which will give problems when mapped (since already
 192     //  on new mesh)
 193     forAll (patch().boundaryMesh().mesh().faceAreas(), i)
 194     {
 195         if (mag(patch().boundaryMesh().mesh().faceAreas()[i]) == 0)
 196         {
 197             Info << "faceArea[active] "<< i << endl;
 198         }
 199     }
 200     if (patch().size() > 0)
 201     {
 202         const vectorField& areas = patch().boundaryMesh().mesh().faceAreas();
 203         initWallSf_ = patch().patchSlice(areas);
 204         initCyclicSf_ = patch().boundaryMesh()
 205         [
 206             cyclicPatchLabel_
 207         ].patchSlice(areas);
 208         nbrCyclicSf_ = refCast<const cyclicFvPatch>
 209         (
 210             patch().boundaryMesh()
 211             [
 212                 cyclicPatchLabel_
 213             ]
 214         ).neighbFvPatch().patch().patchSlice(areas);
 215     }
 216 }
 217
 218 void Foam::activePressureForceBaffleVelocityFvPatchVectorField::rmap
 219 (
 220     const fvPatchVectorField& ptf,
 221     const labelList& addr
 222 )
 223 {
 224     fixedValueFvPatchVectorField::rmap(ptf, addr);
 225
 226     // See autoMap.
 227     const vectorField& areas = patch().boundaryMesh().mesh().faceAreas();
 228     initWallSf_ = patch().patchSlice(areas);
 229     initCyclicSf_ = patch().boundaryMesh()
 230     [
 231         cyclicPatchLabel_
 232     ].patchSlice(areas);
 233     nbrCyclicSf_ = refCast<const cyclicFvPatch>
 234     (
 235         patch().boundaryMesh()
 236         [
 237             cyclicPatchLabel_
 238         ]
 239     ).neighbFvPatch().patch().patchSlice(areas);
 240 }
 241
 242
 243 void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
 244 {
 245     if (updated())
 246     {
 247         return;
 248     }
 249     // Execute the change to the openFraction only once per time-step
 250     if (curTimeIndex_ != this->db().time().timeIndex())
 251     {
 252         const volScalarField& p = db().lookupObject<volScalarField>
 253         (
 254             pName_
 255         );
 256
 257         const fvPatch& cyclicPatch = patch().boundaryMesh()[cyclicPatchLabel_];
 258         const labelList& cyclicFaceCells = cyclicPatch.patch().faceCells();
 259         const fvPatch& nbrPatch = refCast<const cyclicFvPatch>
 260         (
 261             cyclicPatch
 262         ).neighbFvPatch();
 263
 264         const labelList& nbrFaceCells = nbrPatch.patch().faceCells();
 265
 266         scalar valueDiff = 0;
 267
 268         if (fBased_)
 269         {
 270              // Add this side
 271             forAll(cyclicFaceCells, facei)
 272             {
 273                 valueDiff +=p[cyclicFaceCells[facei]]*mag(initCyclicSf_[facei]);
 274             }
 275
 276             // Remove other side
 277             forAll(nbrFaceCells, facei)
 278             {
 279                 valueDiff -=p[nbrFaceCells[facei]]*mag(initCyclicSf_[facei]);
 280             }
 281         }
 282         else //pressure based
 283         {
 284             forAll(cyclicFaceCells, facei)
 285             {
 286                 valueDiff += p[cyclicFaceCells[facei]];
 287             }
 288
 289             forAll(nbrFaceCells, facei)
 290             {
 291                 valueDiff -= p[nbrFaceCells[facei]];
 292             }
 293         }
 294
 295         if ((mag(valueDiff) > mag(minThresholdValue_) || baffleActivated_))
 296         {
 297             openFraction_ =
 298                 max(
 299                     min(
 300                         openFraction_
 301                       + max
 302                         (
 303                           this->db().time().deltaT().value()/openingTime_,
 304                           maxOpenFractionDelta_
 305                         )*(orientation_),
 306                         1 - 1e-6
 307                         ),
 308                     1e-6
 309                     );
 310
 311              baffleActivated_ = true;
 312         }
 313         else
 314         {
 315             openFraction_ = max(min(1 - 1e-6, openFraction_), 1e-6);
 316         }
 317
 318         Info<< "Open fraction = " << openFraction_ << endl;
 319         Info<< "Pressure difference = " << valueDiff << endl;
 320
 321         vectorField::subField Sfw = patch().patch().faceAreas();
 322         vectorField newSfw((1 - openFraction_)*initWallSf_);
 323         forAll(Sfw, facei)
 324         {
 325             Sfw[facei] = newSfw[facei];
 326         }
 327         const_cast<scalarField&>(patch().magSf()) = mag(patch().Sf());
 328
 329         // Update owner side of cyclic
 330         const_cast<vectorField&>(cyclicPatch.Sf()) =
 331             openFraction_*initCyclicSf_;
 332
 333         const_cast<scalarField&>(cyclicPatch.magSf()) =
 334             mag(cyclicPatch.Sf());
 335
 336         // Update neighbour side of cyclic
 337         const_cast<vectorField&>(nbrPatch.Sf()) =
 338             openFraction_*nbrCyclicSf_;
 339
 340         const_cast<scalarField&>(nbrPatch.magSf()) =
 341             mag(nbrPatch.Sf());
 342
 343         curTimeIndex_ = this->db().time().timeIndex();
 344     }
 345
 346     fixedValueFvPatchVectorField::updateCoeffs();
 347 }
 348
 349
 350 void Foam::activePressureForceBaffleVelocityFvPatchVectorField::
 351 write(Ostream& os) const
 352 {
 353     fvPatchVectorField::write(os);
 354     os.writeKeyword("cyclicPatch")
 355         << cyclicPatchName_ << token::END_STATEMENT << nl;
 356     os.writeKeyword("orientation")
 357         << orientation_ << token::END_STATEMENT << nl;
 358     os.writeKeyword("openingTime")
 359         << openingTime_ << token::END_STATEMENT << nl;
 360     os.writeKeyword("maxOpenFractionDelta")
 361         << maxOpenFractionDelta_ << token::END_STATEMENT << nl;
 362     os.writeKeyword("openFraction")
 363         << openFraction_ << token::END_STATEMENT << nl;
 364     os.writeKeyword("p")
 365         << pName_ << token::END_STATEMENT << nl;
 366     os.writeKeyword("minThresholdValue")
 367         << minThresholdValue_ << token::END_STATEMENT << nl;
 368     os.writeKeyword("forceBased")
 369         << fBased_ << token::END_STATEMENT << nl;
 370     writeEntry("value", os);
 371 }
 372
 373
 374 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 375
 376 namespace Foam
 377 {
 378     makePatchTypeField
 379     (
 380         fvPatchVectorField,
 381         activePressureForceBaffleVelocityFvPatchVectorField
 382     );
 383 }
 384
 385 // ************************************************************************* //