src/vofDynamicMesh/lnInclude/threeDofMotion.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "threeDofMotion.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "motionSolver.H"
  30 #include "volFields.H"
  31 #include "surfaceFields.H"
  32 #include "uniformDimensionedFields.H"
  33
  34 #include "tetFemMatrices.H"
  35 #include "tetPointFields.H"
  36 #include "faceTetPolyPatch.H"
  37 #include "fixedValueTetPolyPatchFields.H"
  38
  39 #include "incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H"
  40 #include "incompressible/RAS/RASModel/RASModel.H"
  41 #include "incompressible/LES/LESModel/LESModel.H"
  42
  43
  44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48     defineTypeNameAndDebug(threeDofMotion, 0);
  49
  50     addToRunTimeSelectionTable(dynamicFvMesh, threeDofMotion, IOobject);
  51 } // End namespace Foam
  52
  53
  54 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  55
  56 Foam::tmp<Foam::volSymmTensorField> Foam::threeDofMotion::devRhoReff() const
  57 {
  58     if (this->foundObject<incompressible::RASModel>("RASProperties"))
  59     {
  60         const incompressible::RASModel& ras =
  61             this->lookupObject<incompressible::RASModel>("RASProperties");
  62
  63         return ras.devReff();
  64     }
  65     else if (this->foundObject<incompressible::LESModel>("LESProperties"))
  66     {
  67         const incompressible::LESModel& les =
  68             this->lookupObject<incompressible::LESModel>("LESProperties");
  69
  70         return les.devBeff();
  71     }
  72     else if (this->foundObject<twoPhaseMixture>("transportProperties"))
  73     {
  74         const twoPhaseMixture& twoPhaseProperties =
  75             this->lookupObject<twoPhaseMixture>("transportProperties");
  76
  77         const volVectorField& U = this->lookupObject<volVectorField>("U");
  78
  79         return twoPhaseProperties.nu()*dev(twoSymm(fvc::grad(U)));
  80     }
  81     else
  82     {
  83         FatalErrorIn("floatingBody::devRhoReff()")
  84             << "No valid model for viscous stress calculation."
  85             << exit(FatalError);
  86
  87         return volSymmTensorField::null();
  88     }
  89 }
  90
  91
  92 Foam::vector Foam::threeDofMotion::hullForce() const
  93 {
  94     const fvMesh& mesh = *this;
  95
  96     const volScalarField& alpha = mesh.lookupObject<volScalarField>("alpha1");
  97     const volScalarField& p = mesh.lookupObject<volScalarField>("p");
  98     const volVectorField& U = mesh.lookupObject<volVectorField>("U");
  99
 100     tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
 101     const volSymmTensorField::GeometricBoundaryField& devRhoReffb
 102         = tdevRhoReff().boundaryField();
 103
 104     // Get index of current patch
 105     const label curPatch = hullPatch_.index();
 106
 107     vector pressureForce =
 108         gSum
 109         (
 110             alpha.boundaryField()[curPatch]*
 111             p.boundaryField()[curPatch]*
 112             mesh.Sf().boundaryField()[curPatch]
 113         );
 114
 115     vector viscousForce =
 116         gSum
 117         (
 118             alpha.boundaryField()[curPatch]*
 119             (mesh.Sf().boundaryField()[curPatch] & devRhoReffb[curPatch])
 120         );
 121
 122     vector totalForce = pressureForce + viscousForce;
 123
 124     Info<< "Pressure force = " << pressureForce << nl
 125         << "Viscous force = " << viscousForce << nl
 126         << "Total force = " << totalForce << endl;
 127
 128     return totalForce;
 129 }
 130
 131
 132 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 133
 134 Foam::threeDofMotion::threeDofMotion(const IOobject& io)
 135 :
 136     dynamicFvMesh(io),
 137     dict_
 138     (
 139         IOdictionary
 140         (
 141             IOobject
 142             (
 143                 "dynamicMeshDict",
 144                 time().constant(),
 145                 *this,
 146                 IOobject::MUST_READ,
 147                 IOobject::NO_WRITE
 148             )
 149         ).subDict(typeName + "Coeffs")
 150     ),
 151     equation_
 152     (
 153         IOobject
 154         (
 155             dict_.lookup("name"),
 156             time().timeName(),
 157             *this,
 158             IOobject::MUST_READ,
 159             IOobject::AUTO_WRITE
 160         )
 161     ),
 162     solver_
 163     (
 164         ODESolver::New
 165         (
 166             dict_.lookup("solver"),
 167             equation_
 168         )
 169     ),
 170     epsilon_(readScalar(dict_.lookup("epsilon"))),
 171     motionPtr_(motionSolver::New(*this)),
 172     hullPatch_(word(dict_.lookup("hullPatchName")), boundaryMesh()),
 173     curTimeIndex_(-1),
 174     oldForce_(vector::zero)
 175 {
 176     // Check that the hull patch has been found
 177     if (!hullPatch_.active())
 178     {
 179         FatalErrorIn("threeDofMotion::threeDofMotion(const IOobject& io)")
 180             << "Cannot find hull patch "<< hullPatch_.name() << nl
 181             << "Valid patches are: " << boundaryMesh().names()
 182             << abort(FatalError);
 183     }
 184 }
 185
 186
 187 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 188
 189 Foam::threeDofMotion::~threeDofMotion()
 190 {}
 191
 192
 193 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 194
 195 bool Foam::threeDofMotion::update()
 196 {
 197     // Grab force
 198     vector curForce = hullForce();
 199
 200     const uniformDimensionedVectorField& g =
 201         this->lookupObject<uniformDimensionedVectorField>("g");
 202
 203     equation_.force().value() = 0.5*(curForce + oldForce_)
 204         + equation_.mass().value()*g.value();
 205
 206     if (curTimeIndex_ == -1 || curTimeIndex_ < time().timeIndex())
 207     {
 208         Info << "Storing current force as old force" << endl;
 209         oldForce_ = curForce;
 210         curTimeIndex_ = time().timeIndex();
 211     }
 212
 213     solver_->solve
 214     (
 215         time().value(),
 216         time().value() + time().deltaT().value(),
 217         epsilon_,
 218         time().deltaT().value()
 219     );
 220
 221     // Set the motion onto the motion patch
 222     tetPointVectorField& motionU =
 223         const_cast<tetPointVectorField&>
 224         (
 225             this->objectRegistry::lookupObject<tetPointVectorField>("motionU")
 226         );
 227
 228     fixedValueTetPolyPatchVectorField& motionUHullPatch =
 229         refCast<fixedValueTetPolyPatchVectorField>
 230         (
 231             motionU.boundaryField()[hullPatch_.index()]
 232         );
 233
 234     Info << "Boundary velocity = " << equation_.U().value() << endl;
 235
 236     motionUHullPatch == equation_.U().value();
 237
 238     fvMesh::movePoints(motionPtr_->newPoints());
 239
 240     // Mesh motion only - return false
 241     return false;
 242 }
 243
 244
 245 // ************************************************************************* //