src/dynamicMesh/dynamicFvMesh/dynamicBodyFvMesh/dynamicBodyFvMesh.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 "dynamicBodyFvMesh.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "motionSolver.H"
  30 #include "volFields.H"
  31 #include "mathematicalConstants.H"
  32 #include "tetMotionSolver.H"
  33 #include "laplaceTetMotionSolver.H"
  34 #include "fixedValueTetPolyPatchFields.H"
  35 #include "transformField.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41     defineTypeNameAndDebug(dynamicBodyFvMesh, 0);
  42     addToRunTimeSelectionTable(dynamicFvMesh, dynamicBodyFvMesh, IOobject);
  43 }
  44
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 Foam::dynamicBodyFvMesh::dynamicBodyFvMesh(const IOobject& io)
  49 :
  50     dynamicFvMesh(io),
  51     dynamicMeshCoeffs_
  52     (
  53         IOdictionary
  54         (
  55             IOobject
  56             (
  57                 "dynamicMeshDict",
  58                 io.time().constant(),
  59                 *this,
  60                 IOobject::MUST_READ,
  61                 IOobject::NO_WRITE
  62             )
  63         ).subDict(typeName + "Coeffs")
  64     ),
  65     motionPtr_(motionSolver::New(*this)),
  66     bodyPatchName_
  67     (
  68         dynamicMeshCoeffs_.lookup("bodyPatchName")
  69     ),
  70     bodyPatchID_(-1),
  71     translationDirection_
  72     (
  73         dynamicMeshCoeffs_.lookup("translationDirection")
  74     ),
  75     translationAmplitude_
  76     (
  77         readScalar(dynamicMeshCoeffs_.lookup("translationAmplitude"))
  78     ),
  79     translationFrequency_
  80     (
  81         readScalar(dynamicMeshCoeffs_.lookup("translationFrequency"))
  82     ),
  83     initialRotationOrigin_
  84     (
  85         dynamicMeshCoeffs_.lookup("initialRotationOrigin")
  86     ),
  87     rotationAxis_
  88     (
  89         dynamicMeshCoeffs_.lookup("rotationAxis")
  90     ),
  91     rotationAmplitude_
  92     (
  93         readScalar(dynamicMeshCoeffs_.lookup("rotationAmplitude"))
  94     ),
  95     rotationFrequency_
  96     (
  97         readScalar(dynamicMeshCoeffs_.lookup("rotationFrequency"))
  98     )
  99 {
 100     bodyPatchID_ = boundaryMesh().findPatchID(bodyPatchName_);
 101
 102     if(bodyPatchID_<0)
 103     {
 104         FatalErrorIn
 105         (
 106             "dynamicBodyFvMesh::dynamicBodyFvMesh(const IOobject& io)"
 107         )
 108             << "Can't find patch: " << bodyPatchName_
 109                 << exit(FatalError);
 110     }
 111
 112     translationDirection_ /= mag(translationDirection_) + SMALL;
 113
 114     rotationAxis_ /= mag(rotationAxis_) + SMALL;
 115 }
 116
 117
 118 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 119
 120 Foam::dynamicBodyFvMesh::~dynamicBodyFvMesh()
 121 {}
 122
 123
 124 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 125
 126 bool Foam::dynamicBodyFvMesh::update()
 127 {
 128     scalar curTime = time().value();
 129     scalar oldTime = curTime - time().deltaT().value();
 130
 131     if
 132     (
 133         motionPtr_->type()
 134      == laplaceTetMotionSolver::typeName
 135     )
 136     {
 137         tetMotionSolver& mSolver =
 138             dynamic_cast<tetMotionSolver&>
 139             (
 140                 motionPtr_()
 141             );
 142
 143         vector trans =
 144             translationAmplitude_
 145            *(
 146                 sin(2*mathematicalConstant::pi*translationFrequency_*curTime)
 147               - sin(2*mathematicalConstant::pi*translationFrequency_*oldTime)
 148             )
 149            *translationDirection_;
 150
 151
 152         scalar rotAngle =
 153             rotationAmplitude_
 154            *(
 155                 sin(2*mathematicalConstant::pi*rotationFrequency_*curTime)
 156               - sin(2*mathematicalConstant::pi*rotationFrequency_*oldTime)
 157             );
 158
 159         vector curRotationOrigin =
 160             initialRotationOrigin_
 161           + translationDirection_
 162            *translationAmplitude_
 163            *sin(2*mathematicalConstant::pi*translationFrequency_*oldTime);
 164
 165         const pointField& oldPoints =
 166             mSolver.tetMesh().boundary()[bodyPatchID_].localPoints();
 167
 168         vector r0(1, 1, 1);
 169         r0 -= rotationAxis_*(rotationAxis_ & r0);
 170         r0 /= mag(r0);
 171
 172         // http://mathworld.wolfram.com/RotationFormula.html
 173         vector r1 =
 174             r0*cos(rotAngle)
 175           + rotationAxis_*(rotationAxis_ & r0)*(1 - cos(rotAngle))
 176           + (r0 ^ rotationAxis_)*sin(rotAngle);
 177
 178         tensor T = rotationTensor(r0, r1);
 179
 180         vectorField rot =
 181             transform(T, oldPoints - curRotationOrigin)
 182           + curRotationOrigin
 183           - oldPoints;
 184
 185
 186         tetPointVectorField& motionU = mSolver.motionU();
 187
 188         if
 189         (
 190             motionU.boundaryField()[bodyPatchID_].type()
 191          == fixedValueTetPolyPatchVectorField::typeName
 192         )
 193         {
 194             fixedValueTetPolyPatchVectorField& motionUBodyPatch =
 195                 refCast<fixedValueTetPolyPatchVectorField>
 196                 (
 197                     motionU.boundaryField()[bodyPatchID_]
 198                 );
 199
 200             motionUBodyPatch == (trans + rot)/time().deltaT().value();
 201         }
 202         else
 203         {
 204             FatalErrorIn("dynamicBodyFvMesh::update()")
 205                 << "Bounary condition on " << motionU.name()
 206                     <<  " for " << bodyPatchName_ << " patch is "
 207                     << motionU.boundaryField()[bodyPatchID_].type()
 208                     << ", instead "
 209                     << fixedValueTetPolyPatchVectorField::typeName
 210                     << exit(FatalError);
 211         }
 212     }
 213     else
 214     {
 215         FatalErrorIn("dynamicBodyFvMesh::update()")
 216             << "Selected mesh motion solver is "
 217                 << motionPtr_->type()
 218                 << ", instead "
 219                 << tetMotionSolver::typeName
 220                 << exit(FatalError);
 221     }
 222
 223     fvMesh::movePoints(motionPtr_->newPoints());
 224
 225     // Mesh motion only - return false
 226     return false;
 227 }
 228
 229
 230 // ************************************************************************* //