applications/solvers/surfaceTracking/freeSurface/functionObjects/bubbleHistory/bubbleHistory.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 Author
  25     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  26
  27 \*----------------------------------------------------------------------------*/
  28
  29 #include "bubbleHistory.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "volFields.H"
  32 #include "boundBox.H"
  33 #include "freeSurface.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(bubbleHistory, 0);
  40
  41     addToRunTimeSelectionTable
  42     (
  43         functionObject,
  44         bubbleHistory,
  45         dictionary
  46     );
  47 }
  48
  49
  50 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  51
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 Foam::bubbleHistory::bubbleHistory
  56 (
  57     const word& name,
  58     const Time& t,
  59     const dictionary& dict
  60 )
  61 :
  62     functionObject(name),
  63     name_(name),
  64     time_(t),
  65     regionName_(polyMesh::defaultRegion),
  66     V0_(SMALL),
  67     historyFilePtr_(NULL)
  68 {
  69     if (dict.found("region"))
  70     {
  71         dict.lookup("region") >> regionName_;
  72     }
  73
  74     if (Pstream::master())
  75     {
  76         if (!time_.processorCase())
  77         {
  78             mkDir
  79             (
  80                 time_.path()
  81                /"history"
  82                /time_.timeName()
  83             );
  84
  85             historyFilePtr_ =
  86                 new OFstream
  87                 (
  88                     time_.path()
  89                    /"history"
  90                    /time_.timeName()
  91                    /"history.dat"
  92                 );
  93         }
  94         else
  95         {
  96             mkDir
  97             (
  98                 time_.path()/".."/"history"
  99                /time_.timeName()
 100             );
 101
 102             historyFilePtr_ =
 103                 new OFstream
 104                 (
 105                     time_.path()/".."
 106                    /"history"
 107                    /time_.timeName()
 108                    /"history.dat"
 109                 );
 110         }
 111
 112         (*historyFilePtr_)
 113             << "Time" << tab
 114                 << "Cx" << tab
 115                 << "Cy" << tab
 116                 << "Cz" << tab
 117                 << "Ux" << tab
 118                 << "Uy" << tab
 119                 << "Uz" << tab
 120                 << "ax" << tab
 121                 << "ay" << tab
 122                 << "az" << tab
 123                 << "Fx" << tab
 124                 << "Fy" << tab
 125                 << "Fz" << tab
 126                 << "D" << tab
 127                 << "Lx" << tab
 128                 << "Ly" << tab
 129                 << "Lz" << tab
 130                 << "CD" << tab
 131                 << "CVM" << tab
 132                 << "RE" << tab
 133                 << "WE" << tab
 134                 << "V/V0" << tab
 135                 << "A" << tab
 136                 << "B" << tab
 137                 << "C" << endl;
 138     }
 139
 140     const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName_);
 141
 142     const volScalarField& fluidIndicator =
 143         mesh.lookupObject<volScalarField>("fluidIndicator");
 144
 145     V0_ = gSum((1 - fluidIndicator.internalField())*mesh.V().field());
 146
 147     const freeSurface& fs =
 148         mesh.lookupObject<freeSurface>("freeSurfaceProperties");
 149
 150     if (!fs.twoFluids())
 151     {
 152         V0_ =
 153           - gSum
 154             (
 155                 mesh.Cf().boundaryField()[fs.aPatchID()]
 156               & mesh.Sf().boundaryField()[fs.aPatchID()]
 157             )/3;
 158
 159         if (mesh.nGeometricD() != 3)
 160         {
 161             FatalErrorIn("bubbleHistory::")
 162                 << "One-fluid bubble centroid calc "
 163                     << "is not implemented for 2d mesh"
 164                     << abort(FatalError);
 165         }
 166     }
 167 }
 168
 169
 170 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 171
 172
 173 bool Foam::bubbleHistory::start()
 174 {
 175     const fvMesh& mesh =
 176         time_.lookupObject<fvMesh>(regionName_);
 177
 178     const volScalarField& fluidIndicator =
 179         mesh.lookupObject<volScalarField>("fluidIndicator");
 180
 181     const freeSurface& fs =
 182         mesh.lookupObject<freeSurface>("freeSurfaceProperties");
 183
 184     scalar V = gSum((1 - fluidIndicator.internalField())*mesh.V().field());
 185
 186     if (!fs.twoFluids())
 187     {
 188         V =
 189           - gSum
 190             (
 191                 mesh.Cf().boundaryField()[fs.aPatchID()]
 192                 & mesh.Sf().boundaryField()[fs.aPatchID()]
 193             )/3;
 194
 195         if (mesh.nGeometricD() != 3)
 196         {
 197             FatalErrorIn("bubbleHistory::")
 198                 << "One-fluid bubble centroid calc "
 199                     << "is not implemented for 2d mesh"
 200                     << abort(FatalError);
 201         }
 202     }
 203
 204     const IOdictionary& mrf =
 205         mesh.lookupObject<IOdictionary>("movingReferenceFrame");
 206
 207     dimensionedVector C(mrf.lookup("XF"));
 208     dimensionedVector U(mrf.lookup("UF"));
 209     dimensionedVector a(mrf.lookup("aF"));
 210
 211     vector F = fs.totalViscousForce() + fs.totalPressureForce();
 212
 213     vector dragDir;
 214
 215     if(mag(U.value()) > SMALL)
 216     {
 217         dragDir = -U.value()/mag(U.value());
 218     }
 219     else
 220     {
 221         dragDir = fs.g().value()/(mag(fs.g().value()) + SMALL);
 222     }
 223
 224     scalar dragForce = (dragDir&F);
 225
 226     vector liftForce =  transform(I - dragDir*dragDir, F);
 227
 228     scalar Deq = 2*pow(3*V/(4*M_PI), 1.0/3.0);
 229
 230     scalar Ug = -(U.value()&(fs.g().value()/(mag(fs.g().value()) + SMALL)));
 231
 232     scalar CD =
 233         (4.0/3.0)*(fs.rhoFluidA().value() - fs.rhoFluidB().value())
 234        *mag(fs.g().value())*Deq
 235        /(fs.rhoFluidA().value()*mag(U.value())*Ug + SMALL);
 236
 237     scalar ag = -(a.value()&(fs.g().value()/(mag(fs.g().value()) + SMALL)));
 238
 239     scalar CVM =
 240         (fs.rhoFluidA().value() - fs.rhoFluidB().value())*mag(fs.g().value())
 241        /(fs.rhoFluidA().value()*ag + SMALL)
 242       - (fs.rhoFluidB().value()/(fs.rhoFluidA().value() + SMALL));
 243
 244     scalar RE = Ug*Deq*fs.rhoFluidA().value()/(fs.muFluidA().value() + SMALL);
 245
 246     scalar WE = fs.rhoFluidA().value()*sqr(Ug)*Deq
 247        /(fs.cleanInterfaceSurfTension().value() + SMALL);
 248
 249     boundBox box(mesh.C().boundaryField()[fs.aPatchID()]);
 250
 251     if (Pstream::master())
 252     {
 253         historyFilePtr_->precision(12);
 254
 255         (*historyFilePtr_) << time_.value() << tab
 256             << C.value().x() << tab
 257             << C.value().y() << tab
 258             << C.value().z() << tab
 259             << U.value().x() << tab
 260             << U.value().y() << tab
 261             << U.value().z() << tab
 262             << a.value().x() << tab
 263             << a.value().y() << tab
 264             << a.value().z() << tab
 265             << F.x() << tab
 266             << F.y() << tab
 267             << F.z() << tab
 268             << dragForce << tab
 269             << liftForce.x() << tab
 270             << liftForce.y() << tab
 271             << liftForce.z() << tab
 272             << CD << tab
 273             << CVM << tab
 274             << RE << tab
 275             << WE << tab
 276             << mag(1.0 - V/V0_) << tab
 277             << (box.max().x()-box.min().x())/2 << tab
 278             << (box.max().y()-box.min().y())/2 << tab
 279             << (box.max().z()-box.min().z())/2 << endl;
 280
 281         return true;
 282     }
 283
 284     return false;
 285 }
 286
 287
 288 bool Foam::bubbleHistory::execute()
 289 {
 290     const fvMesh& mesh =
 291         time_.lookupObject<fvMesh>(regionName_);
 292
 293     const volScalarField& fluidIndicator =
 294         mesh.lookupObject<volScalarField>("fluidIndicator");
 295
 296     const freeSurface& fs =
 297         mesh.lookupObject<freeSurface>("freeSurfaceProperties");
 298
 299     scalar V = gSum((1 - fluidIndicator.internalField())*mesh.V().field());
 300
 301     if (!fs.twoFluids())
 302     {
 303         V =
 304           - gSum
 305             (
 306                 mesh.Cf().boundaryField()[fs.aPatchID()]
 307                 & mesh.Sf().boundaryField()[fs.aPatchID()]
 308             )/3;
 309
 310         if (mesh.nGeometricD() != 3)
 311         {
 312             FatalErrorIn("bubbleHistory")
 313                 << "One-fluid bubble centroid calc "
 314                     << "is not implemented for 2d mesh"
 315                     << abort(FatalError);
 316         }
 317     }
 318
 319     const IOdictionary& mrf =
 320         mesh.lookupObject<IOdictionary>("movingReferenceFrame");
 321
 322     dimensionedVector C(mrf.lookup("XF"));
 323     dimensionedVector U(mrf.lookup("UF"));
 324     dimensionedVector a(mrf.lookup("aF"));
 325
 326
 327     vector F = fs.totalViscousForce() + fs.totalPressureForce();
 328
 329     vector dragDir;
 330
 331     if(mag(U.value()) > SMALL)
 332     {
 333         dragDir = -U.value()/mag(U.value());
 334     }
 335     else
 336     {
 337         dragDir = fs.g().value()/(mag(fs.g().value()) + SMALL);
 338     }
 339
 340     scalar dragForce = (dragDir&F);
 341
 342     vector liftForce =  transform(I - dragDir*dragDir, F);
 343
 344     scalar Deq = pow(3*V/(4*M_PI), 1.0/3.0);
 345
 346     scalar Ug = -(U.value()&(fs.g().value()/(mag(fs.g().value()) + SMALL)));
 347
 348     scalar CD =
 349         (4.0/3.0)*(fs.rhoFluidA().value() - fs.rhoFluidB().value())
 350        *mag(fs.g().value())*Deq
 351        /(fs.rhoFluidA().value()*mag(U.value())*Ug + SMALL);
 352
 353     scalar ag = -(a.value()&(fs.g().value()/(mag(fs.g().value()) + SMALL)));
 354
 355     scalar CVM =
 356         (fs.rhoFluidA().value() - fs.rhoFluidB().value())*mag(fs.g().value())
 357        /(fs.rhoFluidA().value()*ag + SMALL)
 358       - (fs.rhoFluidB().value()/(fs.rhoFluidA().value() + SMALL));
 359
 360     scalar RE = Ug*Deq*fs.rhoFluidA().value()/(fs.muFluidA().value() + SMALL);
 361
 362     scalar WE = fs.rhoFluidA().value()*sqr(Ug)*Deq
 363        /(fs.cleanInterfaceSurfTension().value() + SMALL);
 364
 365     boundBox box(mesh.C().boundaryField()[fs.aPatchID()]);
 366
 367     if (Pstream::master())
 368     {
 369         historyFilePtr_->precision(12);
 370
 371         (*historyFilePtr_) << time_.value() << tab
 372             << C.value().x() << tab
 373             << C.value().y() << tab
 374             << C.value().z() << tab
 375             << U.value().x() << tab
 376             << U.value().y() << tab
 377             << U.value().z() << tab
 378             << a.value().x() << tab
 379             << a.value().y() << tab
 380             << a.value().z() << tab
 381             << F.x() << tab
 382             << F.y() << tab
 383             << F.z() << tab
 384             << dragForce << tab
 385             << liftForce.x() << tab
 386             << liftForce.y() << tab
 387             << liftForce.z() << tab
 388             << CD << tab
 389             << CVM << tab
 390             << RE << tab
 391             << WE << tab
 392             << mag(1.0 - V/V0_) << tab
 393             << (box.max().x()-box.min().x())/2 << tab
 394             << (box.max().y()-box.min().y())/2 << tab
 395             << (box.max().z()-box.min().z())/2 << endl;
 396
 397         return true;
 398     }
 399
 400     return false;
 401 }
 402
 403
 404 bool Foam::bubbleHistory::read(const dictionary& dict)
 405 {
 406     if (dict.found("region"))
 407     {
 408         dict.lookup("region") >> regionName_;
 409     }
 410
 411     return true;
 412 }
 413
 414 // ************************************************************************* //