src/postProcessing/functionObjects/forces/forces/forces.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 "forces.H"
  28 #include "volFields.H"
  29 #include "dictionary.H"
  30 #include "Time.H"
  31
  32 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
  33 #include "incompressible/RAS/RASModel/RASModel.H"
  34 #include "incompressible/LES/LESModel/LESModel.H"
  35
  36 #include "basicThermo.H"
  37 #include "compressible/RAS/RASModel/RASModel.H"
  38 #include "compressible/LES/LESModel/LESModel.H"
  39
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45     defineTypeNameAndDebug(forces, 0);
  46 }
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
  51 {
  52     if (obr_.foundObject<compressible::RASModel>("RASProperties"))
  53     {
  54         const compressible::RASModel& ras
  55             = obr_.lookupObject<compressible::RASModel>("RASProperties");
  56
  57         return ras.devRhoReff();
  58     }
  59     else if (obr_.foundObject<incompressible::RASModel>("RASProperties"))
  60     {
  61         const incompressible::RASModel& ras
  62             = obr_.lookupObject<incompressible::RASModel>("RASProperties");
  63
  64         return rho()*ras.devReff();
  65     }
  66     else if (obr_.foundObject<compressible::LESModel>("LESProperties"))
  67     {
  68         const compressible::LESModel& les =
  69         obr_.lookupObject<compressible::LESModel>("LESProperties");
  70
  71         return les.devRhoBeff();
  72     }
  73     else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
  74     {
  75         const incompressible::LESModel& les
  76             = obr_.lookupObject<incompressible::LESModel>("LESProperties");
  77
  78         return rho()*les.devBeff();
  79     }
  80     else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
  81     {
  82         const basicThermo& thermo =
  83              obr_.lookupObject<basicThermo>("thermophysicalProperties");
  84
  85         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
  86
  87         return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
  88     }
  89     else if
  90     (
  91         obr_.foundObject<singlePhaseTransportModel>("transportProperties")
  92     )
  93     {
  94         const singlePhaseTransportModel& laminarT =
  95             obr_.lookupObject<singlePhaseTransportModel>
  96             ("transportProperties");
  97
  98         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
  99
 100         return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
 101     }
 102     else if (obr_.foundObject<dictionary>("transportProperties"))
 103     {
 104         const dictionary& transportProperties =
 105              obr_.lookupObject<dictionary>("transportProperties");
 106
 107         dimensionedScalar nu(transportProperties.lookup("nu"));
 108
 109         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 110
 111         return -rho()*nu*dev(twoSymm(fvc::grad(U)));
 112     }
 113     else
 114     {
 115         FatalErrorIn("forces::devRhoReff()")
 116             << "No valid model for viscous stress calculation."
 117             << exit(FatalError);
 118
 119         return volSymmTensorField::null();
 120     }
 121 }
 122
 123
 124 Foam::tmp<Foam::volScalarField> Foam::forces::rho() const
 125 {
 126     if (rhoName_ == "rhoInf")
 127     {
 128         const fvMesh& mesh = refCast<const fvMesh>(obr_);
 129
 130         return tmp<volScalarField>
 131         (
 132             new volScalarField
 133             (
 134                 IOobject
 135                 (
 136                     "rho",
 137                     mesh.time().timeName(),
 138                     mesh
 139                 ),
 140                 mesh,
 141                 dimensionedScalar("rho", dimDensity, rhoRef_)
 142             )
 143         );
 144     }
 145     else
 146     {
 147         return(obr_.lookupObject<volScalarField>(rhoName_));
 148     }
 149 }
 150
 151
 152 Foam::scalar Foam::forces::rho(const volScalarField& p) const
 153 {
 154     if (p.dimensions() == dimPressure)
 155     {
 156         return 1.0;
 157     }
 158     else
 159     {
 160         if (rhoName_ != "rhoInf")
 161         {
 162             FatalErrorIn("forces::rho(const volScalarField& p)")
 163                 << "Dynamic pressure is expected but kinematic is provided."
 164                 << exit(FatalError);
 165         }
 166
 167         return rhoRef_;
 168     }
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 173
 174 Foam::forces::forces
 175 (
 176     const word& name,
 177     const objectRegistry& obr,
 178     const dictionary& dict,
 179     const bool loadFromFiles
 180 )
 181 :
 182     name_(name),
 183     obr_(obr),
 184     active_(true),
 185     log_(false),
 186     patchSet_(),
 187     pName_(word::null),
 188     UName_(word::null),
 189     rhoName_(word::null),
 190     directForceDensity_(false),
 191     fDName_(""),
 192     rhoRef_(VGREAT),
 193     pRef_(0),
 194     CofR_(vector::zero),
 195     forcesFilePtr_(NULL)
 196 {
 197     // Check if the available mesh is an fvMesh otherise deactivate
 198     if (!isA<fvMesh>(obr_))
 199     {
 200         active_ = false;
 201         WarningIn
 202         (
 203             "Foam::forces::forces"
 204             "("
 205                 "const word&, "
 206                 "const objectRegistry&, "
 207                 "const dictionary&, "
 208                 "const bool"
 209             ")"
 210         )   << "No fvMesh available, deactivating."
 211             << endl;
 212     }
 213
 214     read(dict);
 215 }
 216
 217
 218 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 219
 220 Foam::forces::~forces()
 221 {}
 222
 223
 224 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 225
 226 void Foam::forces::read(const dictionary& dict)
 227 {
 228     if (active_)
 229     {
 230         log_ = dict.lookupOrDefault<Switch>("log", false);
 231
 232         const fvMesh& mesh = refCast<const fvMesh>(obr_);
 233
 234         patchSet_ =
 235             mesh.boundaryMesh().patchSet(wordList(dict.lookup("patches")));
 236
 237         dict.readIfPresent("directForceDensity", directForceDensity_);
 238
 239         if (directForceDensity_)
 240         {
 241             // Optional entry for fDName
 242             fDName_ = dict.lookupOrDefault<word>("fDName", "fD");
 243
 244             // Check whether fDName exists, if not deactivate forces
 245             if
 246             (
 247                 !obr_.foundObject<volVectorField>(fDName_)
 248             )
 249             {
 250                 active_ = false;
 251                 WarningIn("void forces::read(const dictionary& dict)")
 252                 << "Could not find " << fDName_ << " in database." << nl
 253                     << "    De-activating forces."
 254                     << endl;
 255             }
 256         }
 257         else
 258         {
 259             // Optional entries U and p
 260             pName_ = dict.lookupOrDefault<word>("pName", "p");
 261             UName_ = dict.lookupOrDefault<word>("UName", "U");
 262             rhoName_ = dict.lookupOrDefault<word>("rhoName", "rho");
 263
 264             // Check whether UName, pName and rhoName exists,
 265             // if not deactivate forces
 266             if
 267             (
 268                 !obr_.foundObject<volVectorField>(UName_)
 269              || !obr_.foundObject<volScalarField>(pName_)
 270              || (
 271                     rhoName_ != "rhoInf"
 272                  && !obr_.foundObject<volScalarField>(rhoName_)
 273                 )
 274             )
 275             {
 276                 active_ = false;
 277
 278                 WarningIn("void forces::read(const dictionary& dict)")
 279                     << "Could not find " << UName_ << ", " << pName_;
 280
 281                 if (rhoName_ != "rhoInf")
 282                 {
 283                     Info<< " or " << rhoName_;
 284                 }
 285
 286                 Info<< " in database." << nl << "    De-activating forces."
 287                     << endl;
 288             }
 289
 290             // Reference density needed for incompressible calculations
 291             rhoRef_ = readScalar(dict.lookup("rhoInf"));
 292
 293             // Reference pressure, 0 by default
 294             pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
 295         }
 296
 297         // Centre of rotation for moment calculations
 298         CofR_ = dict.lookup("CofR");
 299     }
 300 }
 301
 302
 303 void Foam::forces::makeFile()
 304 {
 305     // Create the forces file if not already created
 306     if (forcesFilePtr_.empty())
 307     {
 308         if (debug)
 309         {
 310             Info<< "Creating forces file." << endl;
 311         }
 312
 313         // File update
 314         if (Pstream::master())
 315         {
 316             fileName forcesDir;
 317             word startTimeName =
 318                 obr_.time().timeName(obr_.time().startTime().value());
 319
 320             if (Pstream::parRun())
 321             {
 322                 // Put in undecomposed case (Note: gives problems for
 323                 // distributed data running)
 324                 forcesDir = obr_.time().path()/".."/name_/startTimeName;
 325             }
 326             else
 327             {
 328                 forcesDir = obr_.time().path()/name_/startTimeName;
 329             }
 330
 331             // Create directory if does not exist.
 332             mkDir(forcesDir);
 333
 334             // Open new file at start up
 335             forcesFilePtr_.reset(new OFstream(forcesDir/(type() + ".dat")));
 336
 337             // Add headers to output data
 338             writeFileHeader();
 339         }
 340     }
 341 }
 342
 343
 344 void Foam::forces::writeFileHeader()
 345 {
 346     if (forcesFilePtr_.valid())
 347     {
 348         forcesFilePtr_()
 349             << "# Time" << tab
 350             << "forces(pressure, viscous) moment(pressure, viscous)"
 351             << endl;
 352     }
 353 }
 354
 355
 356 void Foam::forces::execute()
 357 {
 358     // Do nothing - only valid on write
 359 }
 360
 361
 362 void Foam::forces::end()
 363 {
 364     // Do nothing - only valid on write
 365 }
 366
 367
 368 void Foam::forces::write()
 369 {
 370     if (active_)
 371     {
 372         // Create the forces file if not already created
 373         makeFile();
 374
 375         forcesMoments fm = calcForcesMoment();
 376
 377         if (Pstream::master())
 378         {
 379             forcesFilePtr_() << obr_.time().value() << tab << fm << endl;
 380
 381             if (log_)
 382             {
 383                 Info<< "forces output:" << nl
 384                     << "    forces(pressure, viscous)" << fm.first() << nl
 385                     << "    moment(pressure, viscous)" << fm.second() << nl
 386                     << endl;
 387             }
 388         }
 389     }
 390 }
 391
 392
 393 Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 394 {
 395     forcesMoments fm
 396     (
 397         pressureViscous(vector::zero, vector::zero),
 398         pressureViscous(vector::zero, vector::zero)
 399     );
 400
 401     if (directForceDensity_)
 402     {
 403         const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
 404
 405         const fvMesh& mesh = fD.mesh();
 406
 407         const surfaceVectorField::GeometricBoundaryField& Sfb =
 408             mesh.Sf().boundaryField();
 409
 410         forAllConstIter(labelHashSet, patchSet_, iter)
 411         {
 412             label patchi = iter.key();
 413
 414             vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
 415
 416             scalarField sA = mag(Sfb[patchi]);
 417
 418             // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity)
 419             vectorField fN =
 420                 Sfb[patchi]/sA
 421                *(
 422                     Sfb[patchi] & fD.boundaryField()[patchi]
 423                 );
 424
 425             fm.first().first() += sum(fN);
 426             fm.second().first() += sum(Md ^ fN);
 427
 428             // Tangential force (total force minus normal fN)
 429             vectorField fT = sA*fD.boundaryField()[patchi] - fN;
 430
 431             fm.first().second() += sum(fT);
 432             fm.second().second() += sum(Md ^ fT);
 433         }
 434     }
 435     else
 436     {
 437         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 438         const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
 439
 440         const fvMesh& mesh = U.mesh();
 441
 442         const surfaceVectorField::GeometricBoundaryField& Sfb =
 443             mesh.Sf().boundaryField();
 444
 445         tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
 446         const volSymmTensorField::GeometricBoundaryField& devRhoReffb
 447             = tdevRhoReff().boundaryField();
 448
 449         // Scale pRef by density for incompressible simulations
 450         scalar pRef = pRef_/rho(p);
 451
 452         forAllConstIter(labelHashSet, patchSet_, iter)
 453         {
 454             label patchi = iter.key();
 455
 456             vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
 457
 458             vectorField pf = Sfb[patchi]*(p.boundaryField()[patchi] - pRef);
 459
 460             fm.first().first() += rho(p)*sum(pf);
 461             fm.second().first() += rho(p)*sum(Md ^ pf);
 462
 463             vectorField vf = Sfb[patchi] & devRhoReffb[patchi];
 464
 465             fm.first().second() += sum(vf);
 466             fm.second().second() += sum(Md ^ vf);
 467         }
 468     }
 469
 470     reduce(fm, sumOp());
 471
 472     return fm;
 473 }
 474
 475
 476 // ************************************************************************* //