src/engine/engineTopoChangerMesh/layerARGambit/layerARGambit.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 "layerARGambit.H"
  28 #include "engineTime.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "pointField.H"
  31 #include "mapPolyMesh.H"
  32 #include "polyTopoChange.H"
  33 #include "addToRunTimeSelectionTable.H"
  34 #include "GeometricField.H"
  35 #include "volMesh.H"
  36 #include "fvPatchField.H"
  37 #include "volPointInterpolation.H"
  38 #include "fvcMeshPhi.H"
  39 #include "fvcVolumeIntegrate.H"
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44     defineTypeNameAndDebug(layerARGambit, 0);
  45     addToRunTimeSelectionTable(engineTopoChangerMesh, layerARGambit, IOobject);
  46 }
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 void Foam::layerARGambit::makeLayersLive()
  52 {
  53     const polyTopoChanger& topoChanges = topoChanger_;
  54
  55     // Enable layering
  56     forAll (topoChanges, modI)
  57     {
  58         if (typeid(topoChanges[modI]) == typeid(layerAdditionRemoval))
  59         {
  60             topoChanges[modI].enable();
  61         }
  62         else
  63         {
  64             FatalErrorIn("void Foam::layerARGambit::makeLayersLive()")
  65                 << "Don't know what to do with mesh modifier "
  66                 << modI << " of type " << topoChanges[modI].type()
  67                 << abort(FatalError);
  68         }
  69     }
  70 }
  71
  72
  73 void Foam::layerARGambit::checkAndCalculate()
  74 {
  75     label pistonIndex = -1;
  76     bool foundPiston = false;
  77
  78     label linerIndex = -1;
  79     bool foundLiner = false;
  80
  81     label cylinderHeadIndex = -1;
  82     bool foundCylinderHead = false;
  83
  84     forAll(boundary(), i)
  85     {
  86         Info << boundary()[i].name() << endl;
  87         if (boundary()[i].name() == "piston")
  88         {
  89             pistonIndex = i;
  90             foundPiston = true;
  91         }
  92         else if (boundary()[i].name() == "liner")
  93         {
  94             linerIndex = i;
  95             foundLiner = true;
  96         }
  97         else if (boundary()[i].name() == "cylinderHead")
  98         {
  99             cylinderHeadIndex = i;
 100             foundCylinderHead = true;
 101         }
 102     }
 103
 104     reduce(foundPiston, orOp<bool>());
 105     reduce(foundLiner, orOp<bool>());
 106     reduce(foundCylinderHead, orOp<bool>());
 107
 108     if (!foundPiston)
 109     {
 110         FatalErrorIn("Foam::layerARGambit::checkAndCalculate()")
 111             << " : cannot find piston patch"
 112             << abort(FatalError);
 113     }
 114
 115     if (!foundLiner)
 116     {
 117         FatalErrorIn("Foam::layerARGambit::checkAndCalculate()")
 118             << " : cannot find liner patch"
 119             << abort(FatalError);
 120     }
 121
 122     if (!foundCylinderHead)
 123     {
 124         FatalErrorIn("Foam::layerARGambit::checkAndCalculate()")
 125             << " : cannot find cylinderHead patch"
 126             << exit(FatalError);
 127     }
 128
 129     {
 130         if (linerIndex != -1)
 131         {
 132             pistonPosition() =
 133                 max(boundary()[pistonIndex].patch().localPoints()).z();
 134         }
 135         reduce(pistonPosition(), minOp<scalar>());
 136
 137         if (cylinderHeadIndex != -1)
 138         {
 139             deckHeight() = min
 140             (
 141                 boundary()[cylinderHeadIndex].patch().localPoints()
 142             ).z();
 143         }
 144         reduce(deckHeight(), minOp<scalar>());
 145
 146         Info<< "deckHeight: " << deckHeight() << nl
 147             << "piston position: " << pistonPosition() << endl;
 148     }
 149 }
 150
 151 void Foam::layerARGambit::setVirtualPistonPosition()
 152 {
 153
 154     label pistonFaceIndex = faceZones().findZoneID("pistonLayerFaces");
 155
 156     bool foundPistonFace = (pistonFaceIndex != -1);
 157
 158     Info << "piston face index = " << pistonFaceIndex << endl;
 159
 160     if(!foundPistonFace)
 161     {
 162         FatalErrorIn("Foam::layerARGambit::setVirtualPistonPosition()")
 163             << " : cannot find the pistonLayerFaces"
 164             << exit(FatalError);
 165     }
 166
 167     const labelList& pistonFaces = faceZones()[pistonFaceIndex];
 168     forAll(pistonFaces, i)
 169     {
 170         const face& f = faces()[pistonFaces[i]];
 171
 172         // should loop over facepoints...
 173         forAll(f, j)
 174         {
 175             virtualPistonPosition() =
 176                 Foam::max(virtualPistonPosition(), points()[f[j]].z());
 177         }
 178     }
 179
 180     reduce(virtualPistonPosition(), maxOp<scalar>());
 181
 182 }
 183
 184
 185 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 186
 187 // Construct from components
 188 Foam::layerARGambit::layerARGambit
 189 (
 190     const IOobject& io
 191 )
 192 :
 193     engineTopoChangerMesh(io),
 194     piston_(*this, engTime().engineDict().subDict("piston")),
 195     deformSwitch_(readScalar(engTime().engineDict().lookup("deformAngle"))),
 196     delta_(readScalar(engTime().engineDict().lookup("delta"))),
 197     offSet_(readScalar(engTime().engineDict().lookup("offSet"))),
 198     pistonPosition_(-GREAT),
 199     virtualPistonPosition_(-GREAT),
 200     deckHeight_(GREAT)
 201 {
 202     // Add zones and modifiers if not already there.
 203     addZonesAndModifiers();
 204 }
 205
 206
 207 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 208
 209 Foam::layerARGambit::~layerARGambit()
 210 {}
 211
 212
 213 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 214
 215 bool Foam::layerARGambit::update()
 216 {
 217     scalar deltaZ = engTime().pistonDisplacement().value();
 218     Info<< "deltaZ = " << deltaZ << " Piston at:" << pistonPosition()
 219         << endl;
 220     virtualPistonPosition() += deltaZ;
 221
 222     // Find piston mesh modifier
 223     const label pistonLayerID =
 224         topoChanger_.findModifierID("pistonLayer");
 225
 226     Info << "pistonLayerID: " << pistonLayerID << endl;
 227
 228     const layerAdditionRemoval& pistonLayers =
 229         dynamicCast<const layerAdditionRemoval>
 230         (
 231             topoChanger_[pistonLayerID]
 232         );
 233
 234     bool realDeformation = deformation();
 235
 236     if
 237     (
 238         virtualPistonPosition() + deltaZ
 239      > deckHeight()-engTime().clearance().value() - SMALL
 240     )
 241     {
 242         realDeformation = true;
 243     }
 244
 245     if (realDeformation)
 246     {
 247         // Disable layer addition
 248         Info << "Mesh deformation mode" << endl;
 249         topoChanger_[pistonLayerID].disable();
 250     }
 251     else
 252     {
 253         // Enable layer addition
 254         Info << "Piston layering mode" << endl;
 255         topoChanger_[pistonLayerID].enable();
 256     }
 257
 258     scalar minLayerThickness = pistonLayers.minLayerThickness();
 259
 260     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 261
 262     bool meshChanged = topoChangeMap->morphing();
 263
 264     pointField newPoints = allPoints();
 265
 266     if (meshChanged)
 267     {
 268         if (topoChangeMap().hasMotionPoints())
 269         {
 270             movePoints(topoChangeMap().preMotionPoints());
 271             newPoints = topoChangeMap().preMotionPoints();
 272         }
 273         setV0();
 274         resetMotion();
 275     }
 276
 277     Info << "virtualPistonPosition = " << virtualPistonPosition()
 278     << ", deckHeight = " << deckHeight() << endl;
 279
 280     // Mesh in three parts:
 281     // - pistonPoints - move with deltaZ
 282     // - headPoints - do not move
 283
 284     const pointZoneMesh& pZones = pointZones();
 285     label headPtsIndex = pZones.findZoneID("headPoints");
 286     label pistonPtsIndex = pZones.findZoneID("pistonPoints");
 287     const labelList& pistonPoints = pZones[pistonPtsIndex];
 288     const labelList& headPoints = pZones[headPtsIndex];
 289
 290
 291     // Whether point displacement is by scaling
 292     boolList scaleDisp(nPoints(), true);
 293     label nScaled = nPoints();
 294     List<bool> pistonPoint(newPoints.size(), false);
 295     List<bool> headPoint(newPoints.size(), false);
 296
 297     forAll(pistonPoints, i)
 298     {
 299         label pointI = pistonPoints[i];
 300         pistonPoint[pointI] = true;
 301         point& p = newPoints[pointI];
 302
 303         if (p.z() < pistonPosition() - 1.0e-6)
 304         {
 305             scaleDisp[pointI] = false;
 306             nScaled--;
 307         }
 308     }
 309
 310     forAll(headPoints, i)
 311     {
 312         headPoint[headPoints[i]] = true;
 313         scaleDisp[headPoints[i]] = false;
 314         nScaled--;
 315     }
 316
 317     if (realDeformation)
 318     {
 319
 320         forAll(scaleDisp, pointI)
 321         {
 322             point& p = newPoints[pointI];
 323
 324             if (scaleDisp[pointI])
 325             {
 326                 p.z() +=
 327                     deltaZ
 328                   * (deckHeight() - p.z())/(deckHeight() - pistonPosition());
 329             }
 330             else
 331             {
 332                 if (!headPoint[pointI])
 333                 {
 334                     p.z() += deltaZ;
 335                 }
 336             }
 337         }
 338     }
 339     else
 340     {
 341
 342         // Always move piston
 343         scalar pistonTopZ = -GREAT;
 344         forAll(pistonPoints, i)
 345         {
 346             point& p = newPoints[pistonPoints[i]];
 347             p.z() += deltaZ;
 348             pistonTopZ = max(pistonTopZ, p.z());
 349         }
 350
 351
 352         // NN! fix. only needed for compression
 353         if (deltaZ > 0.0)
 354         {
 355             // check if piston-points have moved beyond the layer above
 356             forAll(newPoints, i)
 357             {
 358                 if (!pistonPoint[i])
 359                 {
 360                     if (virtualPistonPosition() > newPoints[i].z())
 361                     {
 362                         newPoints[i].z() = pistonTopZ + 0.9*minLayerThickness;
 363                     }
 364                 }
 365             }
 366         }
 367     }
 368
 369     movePoints(newPoints);
 370
 371     pistonPosition() += deltaZ;
 372     scalar pistonSpeed = deltaZ/engTime().deltaT().value();
 373
 374     Info<< "clearance: " << deckHeight() - pistonPosition() << nl
 375         << "Piston speed = " << pistonSpeed << " m/s" << endl;
 376
 377     Info << "Total cylinder volume at CA " << engTime().timeName() << " = " <<
 378     sum(V()) << endl;
 379
 380     return meshChanged;
 381
 382 }
 383
 384 void Foam::layerARGambit::setBoundaryVelocity(volVectorField& U)
 385 {
 386 // Does nothing, using the movingWallVelocity boundary condition for U in the piston patch...
 387
 388 //    vector pistonVel = piston().cs().axis()*engTime().pistonSpeed().value();
 389
 390     //mean piston velocityy
 391 /*
 392     vector pistonVel = 0.5 * piston().cs().axis()*
 393                             dimensionedScalar
 394                             (
 395                                 "meanPistonSpeed",
 396                                 dimLength/dimTime,
 397                                 2.0*engTime().stroke().value()*engTime().rpm().value()/60.0
 398                             ).value();
 399 */
 400
 401 //    U.boundaryField()[piston().patchID().index()] = pistonVel;
 402 }
 403
 404 // ************************************************************************* //