src/vofDynamicMesh/floatingBody/floatingBody.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 Author
  26     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "volFields.H"
  31 #include "surfaceFields.H"
  32 #include "floatingBody.H"
  33 #include "transform.H"
  34 #include "transformField.H"
  35 #include "fixedValueTetPolyPatchFields.H"
  36
  37 #include "incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H"
  38 #include "incompressible/RAS/RASModel/RASModel.H"
  39 #include "incompressible/LES/LESModel/LESModel.H"
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::floatingBody::check() const
  44 {
  45     forAll (hullPatches_, patchI)
  46     {
  47         if (!hullPatches_[patchI].active())
  48         {
  49             FatalErrorIn("void floatingBody::check() const")
  50                 << "Hull patch named " << hullPatches_[patchI].name()
  51                 << " not found.  Available patches: "
  52                 << mesh_.boundaryMesh().names()
  53                 << abort(FatalError);
  54         }
  55     }
  56 }
  57
  58
  59 Foam::tmp<Foam::volSymmTensorField> Foam::floatingBody::devRhoReff() const
  60 {
  61     if (mesh_.foundObject<incompressible::RASModel>("RASProperties"))
  62     {
  63         const incompressible::RASModel& ras =
  64             mesh_.lookupObject<incompressible::RASModel>("RASProperties");
  65
  66         return ras.devReff();
  67     }
  68     else if (mesh_.foundObject<incompressible::LESModel>("LESProperties"))
  69     {
  70         const incompressible::LESModel& les =
  71             mesh_.lookupObject<incompressible::LESModel>("LESProperties");
  72
  73         return les.devBeff();
  74     }
  75     else if (mesh_.foundObject<twoPhaseMixture>("transportProperties"))
  76     {
  77         const twoPhaseMixture& twoPhaseProperties =
  78             mesh_.lookupObject<twoPhaseMixture>("transportProperties");
  79
  80         const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
  81
  82         return twoPhaseProperties.nu()*dev(twoSymm(fvc::grad(U)));
  83     }
  84     else
  85     {
  86         FatalErrorIn("floatingBody::devRhoReff()")
  87             << "No valid model for viscous stress calculation."
  88             << exit(FatalError);
  89
  90         return volSymmTensorField::null();
  91     }
  92 }
  93
  94
  95 Foam::vector Foam::floatingBody::dragForce() const
  96 {
  97     // Note: Return zero if there is no flow solution
  98     // HR, 6/Jan/2010
  99     if
 100     (
 101         !mesh_.foundObject<volVectorField>("U")
 102         || !mesh_.foundObject<volScalarField>("p")
 103     )
 104     {
 105         WarningIn
 106         (
 107             "floatingBody::dragForce() const"
 108         )   << "No flow solution found. Returning zero."
 109             << endl;
 110
 111         return vector::zero;
 112     }
 113
 114     const volScalarField& alpha = mesh_.lookupObject<volScalarField>("alpha1");
 115     const volScalarField& p = mesh_.lookupObject<volScalarField>("p");
 116     const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
 117
 118     const surfaceVectorField::GeometricBoundaryField& Sfb =
 119         mesh_.Sf().boundaryField();
 120
 121     tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
 122     const volSymmTensorField::GeometricBoundaryField& devRhoReffb
 123         = tdevRhoReff().boundaryField();
 124
 125     vector pressureForce = vector::zero;
 126     vector viscousForce = vector::zero;
 127
 128     forAll (hullPatches_, patchI)
 129     {
 130         const label curPatch = hullPatches_[patchI].index();
 131
 132         // Note: calculating forces on the wetted surface only
 133         // HJ, 9/Mar/2009
 134
 135         pressureForce +=
 136             gSum
 137             (
 138                 alpha.boundaryField()[curPatch]*
 139                 p.boundaryField()[curPatch]*Sfb[curPatch]
 140             );
 141
 142         viscousForce +=
 143             gSum
 144             (
 145                 alpha.boundaryField()[curPatch]*
 146                 (Sfb[curPatch] & devRhoReffb[curPatch])
 147             );
 148     }
 149
 150     vector totalForce = pressureForce + viscousForce;
 151
 152     Info<< "Pressure force = " << pressureForce << nl
 153         << "Viscous force = " << viscousForce << nl
 154         << "Total force = " << totalForce << endl;
 155
 156     return totalForce;
 157 }
 158
 159
 160 Foam::vector Foam::floatingBody::dragMoment() const
 161 {
 162     // Note: Return zero if there is no flow solution
 163     // HR, 6/Jan/2010
 164     if
 165     (
 166         !mesh_.foundObject<volVectorField>("U")
 167         || !mesh_.foundObject<volScalarField>("p")
 168     )
 169     {
 170         WarningIn
 171         (
 172             "floatingBody::dragForce() const"
 173         )   << "No flow solution found. Returning zero."
 174             << endl;
 175
 176         return vector::zero;
 177     }
 178
 179     const volScalarField& alpha = mesh_.lookupObject<volScalarField>("alpha1");
 180     const volScalarField& p = mesh_.lookupObject<volScalarField>("p");
 181     const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
 182
 183     const surfaceVectorField::GeometricBoundaryField& Sfb =
 184         mesh_.Sf().boundaryField();
 185
 186     tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
 187     const volSymmTensorField::GeometricBoundaryField& devRhoReffb
 188         = tdevRhoReff().boundaryField();
 189
 190     vector pressureMoment = vector::zero;
 191     vector viscousMoment = vector::zero;
 192
 193     forAll (hullPatches_, patchI)
 194     {
 195         // Get index of current patch
 196         const label curPatch = hullPatches_[patchI].index();
 197
 198         // Note: calculating forces on the wetted surface only
 199         // HJ, 9/Mar/2009
 200
 201         pressureMoment +=
 202             gSum
 203             (
 204                 (mesh_.Cf().boundaryField()[curPatch] - sixDOF_.X().value())
 205               ^ (
 206                     alpha.boundaryField()[curPatch]*
 207                     p.boundaryField()[curPatch]*Sfb[curPatch]
 208                 )
 209             );
 210
 211         viscousMoment +=
 212             gSum
 213             (
 214                 (mesh_.Cf().boundaryField()[curPatch] - sixDOF_.X().value())
 215               ^ (
 216                     alpha.boundaryField()[curPatch]*
 217                     (Sfb[curPatch] & devRhoReffb[curPatch])
 218                 )
 219             );
 220     }
 221
 222     vector totalMoment = pressureMoment + viscousMoment;
 223
 224     Info<< "Pressure moment = " << pressureMoment << nl
 225         << "Viscous moment = " << viscousMoment << nl
 226         << "Total moment = " << totalMoment << endl;
 227
 228     return totalMoment;
 229 }
 230
 231
 232 void Foam::floatingBody::fixTranslation(vector& v) const
 233 {
 234     if (fixedSurge_)
 235     {
 236         v.x() = 0;
 237     }
 238
 239     if (fixedSway_)
 240     {
 241         v.y() = 0;
 242     }
 243
 244     if (fixedHeave_)
 245     {
 246         v.z() = 0;
 247     }
 248 }
 249
 250
 251 void Foam::floatingBody::fixRotation(vector& rot) const
 252 {
 253     if (fixedRoll_)
 254     {
 255         rot.x() = 0;
 256     }
 257
 258     if (fixedPitch_)
 259     {
 260         rot.y() = 0;
 261     }
 262
 263     if (fixedYaw_)
 264     {
 265         rot.z() = 0;
 266     }
 267
 268 }
 269
 270
 271 Foam::vector Foam::floatingBody::rotCentreVelocity() const
 272 {
 273     return
 274     (
 275         // New position of rotation centroid in absolute coordinate system
 276         transform(sixDOF_.toAbsolute(), centreOfSlider_)
 277         // Motion of centre of gravity
 278       + sixDOF_.X().value()
 279         // Old position of rotation centroid
 280       - rotCentre_
 281     )/mesh_.time().deltaT().value();
 282 }
 283
 284
 285 Foam::tmp<Foam::vectorField> Foam::floatingBody::boundaryVelocity
 286 (
 287     const vectorField& boundaryPoints
 288 ) const
 289 {
 290     return
 291         // Translation of centroid
 292         rotCentreVelocity()
 293         // Rotation around centroid
 294       + (
 295             // New points,scaled around centroid
 296             transform
 297             (
 298                 sixDOF_.rotIncrementTensor().T(),
 299                 boundaryPoints - rotCentre_
 300             )
 301             // Old points, scaled around centroid
 302           - (boundaryPoints - rotCentre_)
 303         )/mesh_.time().deltaT().value();
 304 }
 305
 306
 307 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 308
 309 Foam::floatingBody::floatingBody
 310 (
 311     const word& name,
 312     const fvMesh& mesh,
 313     const dictionary& dict
 314 )
 315 :
 316     name_(name),
 317     mesh_(mesh),
 318
 319     hullPatches_(),
 320     hullSlider_(word(dict.lookup("hullSlider")), mesh_.boundaryMesh()),
 321     hullSliderZone_(word(dict.lookup("hullSliderZone")), mesh_.faceZones()),
 322     fixedSlider_(word(dict.lookup("fixedSlider")), mesh_.boundaryMesh()),
 323     fixedSliderZone_(word(dict.lookup("fixedSliderZone")), mesh_.faceZones()),
 324     centreOfSlider_(dict.lookup("centreOfSlider")),
 325
 326     addPropulsionForce_(dict.lookup("addPropulsionForce")),
 327     propulsionDirection_(dict.lookup("propulsionDirection")),
 328     propulsionForceCentroid_(dict.lookup("propulsionForceCentroid")),
 329
 330     sixDOF_
 331     (
 332         IOobject
 333         (
 334             name,
 335             mesh_.time().timeName(),
 336             "uniform",
 337             mesh_,
 338             IOobject::MUST_READ,
 339             IOobject::AUTO_WRITE
 340         )
 341     ),
 342     solver_
 343     (
 344         ODESolver::New
 345         (
 346             dict.lookup("odeSolver"),
 347             sixDOF_
 348         )
 349     ),
 350     epsilon_(readScalar(dict.lookup("epsilon"))),
 351
 352     fixedSurge_(dict.lookup("fixedSurge")),
 353     fixedSway_(dict.lookup("fixedSway")),
 354     fixedHeave_(dict.lookup("fixedHeave")),
 355     fixedRoll_(dict.lookup("fixedRoll")),
 356     fixedPitch_(dict.lookup("fixedPitch")),
 357     fixedYaw_(dict.lookup("fixedYaw")),
 358
 359     steadyState_(dict.lookup("steadyState")),
 360     curTimeIndex_(-1),
 361     rotCentre_(sixDOF_.X().value()),
 362     oldForce_(vector::zero),
 363     oldMoment_(vector::zero),
 364     // Dealing with writing a motion file in parallel.  HJ, 14/Aug/2008
 365     of_(mesh_.time().caseConstant()/".."/name_ + ".dat")
 366 {
 367     // Read hull patch names
 368     wordList hullPatchNames(dict.lookup("hullPatches"));
 369
 370     hullPatches_.setSize(hullPatchNames.size());
 371
 372     forAll (hullPatchNames, patchI)
 373     {
 374         hullPatches_.set
 375         (
 376             patchI,
 377             new polyPatchID
 378             (
 379                 hullPatchNames[patchI],
 380                 mesh_.boundaryMesh()
 381             )
 382         );
 383     }
 384
 385     // Rescale propulsion direction vector
 386     if (mag(propulsionDirection_) < SMALL)
 387     {
 388         FatalErrorIn
 389         (
 390             "floatingBody::floatingBody(const word& name, const fvMesh& mesh,"
 391             "const dictionary& dict)"
 392         )   << " Incorrect propulsion direction: " << propulsionDirection_
 393             << abort(FatalError);
 394     }
 395     else
 396     {
 397         propulsionDirection_ /= mag(propulsionDirection_);
 398     }
 399
 400     // Block (initial) motion components as requested
 401
 402     // Translation
 403     fixTranslation(sixDOF_.U().value());
 404
 405     // Rotation
 406     fixRotation(sixDOF_.omega().value());
 407
 408     check();
 409 }
 410
 411
 412 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 413
 414 void Foam::floatingBody::setMotion(tetPointVectorField& motionU)
 415 {
 416     // Calculate forces and moments
 417     vector curForce = dragForce();
 418     vector curMoment = dragMoment();
 419
 420     if (addPropulsionForce_)
 421     {
 422         // Here, curForce = drag.  HJ, 11/Mar/2009
 423
 424         const vector propulsionForce =
 425             - propulsionDirection_*(propulsionDirection_ & curForce);
 426
 427         const vector propulsionMoment =
 428             (
 429                 (propulsionForceCentroid_ - sixDOF_.X().value())
 430               ^ propulsionForce
 431             );
 432
 433         Info<< "Propulsion force = " << propulsionForce
 434             << " Propulsion moment = " << propulsionMoment << endl;
 435
 436         curForce += propulsionForce;
 437         curMoment += propulsionMoment;
 438     }
 439
 440     const uniformDimensionedVectorField& g =
 441         mesh_.lookupObject<uniformDimensionedVectorField>("g");
 442
 443     // Grab initial force for the first time-step
 444     if (curTimeIndex_ == -1)
 445     {
 446         oldForce_ = curForce;
 447         oldMoment_ = curMoment;
 448         curTimeIndex_ = mesh_.time().timeIndex();
 449     }
 450
 451     // Set force and moment
 452     sixDOF_.force().value() =
 453         0.5*(curForce + oldForce_)
 454       + sixDOF_.mass().value()*g.value();
 455
 456     sixDOF_.moment().value() =
 457         0.5*(curMoment + oldMoment_);
 458
 459     // Manage multiple calls to mesh motion within a single time-step
 460     if (curTimeIndex_ < mesh_.time().timeIndex())
 461     {
 462         oldForce_ = curForce;
 463         oldMoment_ = curMoment;
 464         curTimeIndex_ = mesh_.time().timeIndex();
 465     }
 466
 467     // Manage steady-state: reset velocities to zero
 468     if (steadyState_)
 469     {
 470         sixDOF_.U().value() = vector::zero;
 471         sixDOF_.omega().value() = vector::zero;
 472     }
 473     else
 474     {
 475         // Block motion components as requested
 476
 477         // Translation
 478         fixTranslation(sixDOF_.force().value());
 479
 480         // Rotation
 481         fixRotation(sixDOF_.moment().value());
 482     }
 483
 484     // Grab centre of rotation before the ODE solution
 485     rotCentre_ = sixDOF_.X().value()
 486         + transform(sixDOF_.toAbsolute(), centreOfSlider_);
 487
 488     // Solve for motion
 489     solver_->solve
 490     (
 491         mesh_.time().value(),
 492         mesh_.time().value() +  mesh_.time().deltaT().value(),
 493         epsilon_,
 494         mesh_.time().deltaT().value()
 495     );
 496
 497     Info<< "body: " << name() << nl
 498         << "f = " << sixDOF_.force().value() << tab
 499         << "mg = " << sixDOF_.mass().value()*g.value() << tab
 500         << "M = " << sixDOF_.moment().value() << nl
 501         << "x = " << sixDOF_.X().value() << tab
 502         << "u = " << sixDOF_.Uaverage().value() << nl
 503         << "rot vector = " << sixDOF_.rotVector() << tab
 504         << "rot angle = " << sixDOF_.rotAngle().value() << tab
 505         << "omega = " << sixDOF_.omegaAverageAbsolute().value() << endl;
 506
 507     // Set boundary motion on the hull patches
 508     forAll (hullPatches_, patchI)
 509     {
 510         fixedValueTetPolyPatchVectorField& motionUHull =
 511             refCast<fixedValueTetPolyPatchVectorField>
 512             (
 513                 motionU.boundaryField()[hullPatches_[patchI].index()]
 514             );
 515
 516         // Set motion based on current hull geometry
 517         motionUHull ==
 518             boundaryVelocity
 519             (
 520                 motionUHull.patch().localPoints()
 521             );
 522     }
 523
 524     // Set motion on slider surfaces
 525     if (hullSlider_.active())
 526     {
 527         fixedValueTetPolyPatchVectorField& motionUHullSlider =
 528             refCast<fixedValueTetPolyPatchVectorField>
 529             (
 530                 motionU.boundaryField()[hullSlider_.index()]
 531             );
 532
 533         motionUHullSlider ==
 534             boundaryVelocity
 535             (
 536                 motionUHullSlider.patch().localPoints()
 537             );
 538     }
 539
 540     if (fixedSlider_.active())
 541     {
 542         fixedValueTetPolyPatchVectorField& motionUFixedSlider =
 543             refCast<fixedValueTetPolyPatchVectorField>
 544             (
 545                 motionU.boundaryField()[fixedSlider_.index()]
 546             );
 547
 548         motionUFixedSlider == rotCentreVelocity();
 549     }
 550
 551
 552     // Report position and rotation.  Write only if master
 553     if (Pstream::master())
 554     {
 555         of_ << mesh_.time().value() << tab;
 556
 557         // Translation
 558         if (!fixedSurge_)
 559         {
 560             of_ << sixDOF_.X().value().x() << tab;
 561         }
 562
 563         if (!fixedSway_)
 564         {
 565             of_ << sixDOF_.X().value().y() << tab;
 566         }
 567
 568         if (!fixedHeave_)
 569         {
 570             of_ << sixDOF_.X().value().z() << tab;
 571         }
 572
 573         // Rotation
 574         if (!fixedRoll_ || !fixedPitch_ || !fixedYaw_)
 575         {
 576             of_ << sixDOF_.rotAngle().value() << tab;
 577         }
 578
 579         of_ << endl;
 580     }
 581 }
 582
 583
 584 void Foam::floatingBody::correctMotion(vectorField& newAllPoints) const
 585 {
 586     // Adjust motion of slider zone points
 587     if (hullSliderZone_.active())
 588     {
 589         // Hull slider zone takes translation and rotation
 590         const pointField newP =
 591             boundaryVelocity
 592             (
 593                 mesh_.faceZones()[hullSliderZone_.index()]().localPoints()
 594             );
 595
 596         const labelList& addr =
 597             mesh_.faceZones()[hullSliderZone_.index()]().meshPoints();
 598
 599         forAll (addr, i)
 600         {
 601             if (addr[i] >= mesh_.nPoints())
 602             {
 603                 newAllPoints[addr[i]] += newP[i]*mesh_.time().deltaT().value();
 604             }
 605         }
 606     }
 607
 608     if (fixedSliderZone_.active())
 609     {
 610         // Fixed slider zone takes translation only
 611         const vector rcv = rotCentreVelocity();
 612
 613         const labelList& addr =
 614             mesh_.faceZones()[fixedSliderZone_.index()]().meshPoints();
 615
 616         forAll (addr, i)
 617         {
 618             if (addr[i] >= mesh_.nPoints())
 619             {
 620                 newAllPoints[addr[i]] += rcv*mesh_.time().deltaT().value();
 621             }
 622         }
 623     }
 624 }
 625
 626
 627 // ************************************************************************* //