src/engine/engineTopoChangerMesh/simpleTwoStroke/addSimpleTwoStrokeModifiers.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "simpleTwoStroke.H"
  27 #include "slidingInterface.H"
  28 #include "layerAdditionRemoval.H"
  29 #include "surfaceFields.H"
  30 #include "regionSplit.H"
  31
  32 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  33
  34 void Foam::simpleTwoStroke::addZonesAndModifiers()
  35 {
  36     // Add the zones and mesh modifiers to operate piston motion
  37
  38     if
  39     (
  40         pointZones().size() > 0
  41      || faceZones().size() > 0
  42      || cellZones().size() > 0
  43     )
  44     {
  45         Info<< "Time = " << engTime().theta() << endl;
  46         Info<< "void simpleTwoStroke::addZonesAndModifiers() : "
  47             << "Zones and modifiers already present.  Skipping."
  48             << endl;
  49
  50         if (topoChanger_.size() == 0)
  51         {
  52             FatalErrorIn
  53             (
  54                 "void simpleTwoStroke::addZonesAndModifiers()"
  55             )   << "Mesh modifiers not read properly"
  56                 << abort(FatalError);
  57         }
  58
  59         setVirtualPistonPosition();
  60         checkAndCalculate();
  61
  62         return;
  63
  64
  65     }
  66
  67     Info << "checkAndCalculate()" << endl;
  68     checkAndCalculate();
  69
  70     Info<< "Time = " << engTime().theta() << endl
  71         << "Adding zones to the engine mesh" << endl;
  72
  73
  74     //fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
  75     //pz = 2: piston points, cutPointZone
  76     //cz = 1: moving mask
  77
  78     List<pointZone*> pz(3);
  79     List<faceZone*> fz(4);
  80     List<cellZone*> cz(1);
  81
  82     label nPointZones = 0;
  83     label nFaceZones = 0;
  84     label nCellZones = 0;
  85
  86     // Add the piston zone
  87     if (piston().patchID().active())
  88     {
  89
  90         // Piston position
  91
  92         Info << "Adding face zone for piston layer addition/removal" << endl;
  93
  94         label pistonPatchID = piston().patchID().index();
  95
  96         scalar zPist =
  97             max(boundary()[pistonPatchID].patch().localPoints()).z();
  98
  99         scalar zPistV = zPist + offSet();
 100
 101         labelList zone1(faceCentres().size());
 102         boolList flipZone1(faceCentres().size(), false);
 103         label nZoneFaces1 = 0;
 104
 105         bool foundAtLeastOne = false;
 106         scalar zHigher = GREAT;
 107         scalar dh = GREAT;
 108         scalar dl = GREAT;
 109
 110         forAll (faceCentres(), faceI)
 111         {
 112             // The points have to be in the cylinder and not in the ports....
 113
 114             scalar zc = faceCentres()[faceI].z();
 115
 116             scalar xc = faceCentres()[faceI].x();
 117             scalar yc = faceCentres()[faceI].y();
 118
 119             vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
 120             scalar dd = n & vector(0,0,1);
 121
 122             if(sqrt(sqr(xc)+sqr(yc)) <  0.5 * engTime().bore().value())
 123             {
 124                 if (dd > 0.1)
 125                 {
 126                     if (zPistV - zc > 0 && zPistV - zc < dl)
 127                     {
 128                         dl = zPistV - zc;
 129                     }
 130
 131                     if (zc - zPistV > 0 && zc - zPistV < dh)
 132                     {
 133                         zHigher = zc;
 134                         dh = zc - zHigher;
 135                     }
 136
 137                     if
 138                     (
 139                         zc > zPistV - delta()
 140                         && zc < zPistV + delta()
 141                     )
 142                     {
 143                         foundAtLeastOne = true;
 144                         if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
 145                         {
 146                             flipZone1[nZoneFaces1] = true;
 147                         }
 148
 149                         zone1[nZoneFaces1] = faceI;
 150                         nZoneFaces1++;
 151                     }
 152                 }
 153             }
 154         }
 155
 156         // if no cut was found use the layer above
 157         if (!foundAtLeastOne)
 158         {
 159             zPistV = zHigher;
 160
 161             forAll (faceCentres(), faceI)
 162             {
 163                 scalar zc = faceCentres()[faceI].z();
 164
 165                 scalar xc = faceCentres()[faceI].x();
 166                 scalar yc = faceCentres()[faceI].y();
 167
 168                 vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
 169                 scalar dd = n & vector(0,0,1);
 170
 171                 if(sqrt(sqr(xc)+sqr(yc)) <  0.5 * engTime().bore().value())
 172                 {
 173                     if (dd > 0.1)
 174                     {
 175
 176                         if
 177                         (
 178                             zc > zPistV - delta()
 179                             && zc < zPistV + delta()
 180                         )
 181                         {
 182                             if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
 183                             {
 184                                 flipZone1[nZoneFaces1] = true;
 185                             }
 186
 187                             zone1[nZoneFaces1] = faceI;
 188                             nZoneFaces1++;
 189                         }
 190                     }
 191                 }
 192
 193             }
 194         }
 195
 196         zone1.setSize(nZoneFaces1);
 197         flipZone1.setSize(nZoneFaces1);
 198
 199         fz[nFaceZones] =
 200             new faceZone
 201             (
 202                 "pistonLayerFaces",
 203                 zone1,
 204                 flipZone1,
 205                 nFaceZones,
 206                 faceZones()
 207             );
 208
 209         nFaceZones++;
 210
 211         // Points which don't move (= cylinder head)
 212         DynamicList<label> headPoints(nPoints() / 10);
 213
 214         label nHeadPoints = 0;
 215         forAll (points(), pointI)
 216         {
 217             scalar zCoord = points()[pointI].z();
 218
 219             if (zCoord > deckHeight() - delta())
 220             {
 221                 headPoints.append(pointI);
 222                 nHeadPoints++;
 223             }
 224         }
 225
 226         Info << "Number of head points = " << nHeadPoints << endl;
 227         pz[nPointZones] =
 228             new pointZone
 229             (
 230                 "headPoints",
 231                 headPoints.shrink(),
 232                 nPointZones,
 233                 pointZones()
 234             );
 235
 236         nPointZones++;
 237
 238     }
 239
 240     //  Sliding interface for scavenging ports
 241
 242     if (foundScavPorts())
 243     {
 244         // Inner slider
 245
 246         const polyPatch& innerScav =
 247             boundaryMesh()[boundaryMesh().findPatchID(scavInCylPatchName_)];
 248
 249         labelList isf(innerScav.size());
 250
 251         forAll (isf, i)
 252         {
 253             isf[i] = innerScav.start() + i;
 254         }
 255
 256         fz[nFaceZones] = new faceZone
 257         (
 258             scavInCylPatchName_ + "Zone",
 259             isf,
 260             boolList(innerScav.size(), false),
 261             nFaceZones,
 262             faceZones()
 263         );
 264
 265         nFaceZones++;
 266
 267         // Outer slider
 268
 269         const polyPatch& outerScav =
 270             boundaryMesh()[boundaryMesh().findPatchID(scavInPortPatchName_)];
 271
 272         labelList osf(outerScav.size());
 273
 274         forAll (osf, i)
 275         {
 276             osf[i] = outerScav.start() + i;
 277         }
 278
 279         fz[nFaceZones] = new faceZone
 280         (
 281             scavInPortPatchName_ + "Zone",
 282             osf,
 283             boolList(outerScav.size(), false),
 284             nFaceZones,
 285             faceZones()
 286         );
 287
 288         nFaceZones++;
 289
 290         fz[nFaceZones] = new faceZone
 291         (
 292             "cutFaceZone",
 293             labelList(0),
 294             boolList(0, false),
 295             nFaceZones,
 296             faceZones()
 297         );
 298
 299         nFaceZones++;
 300
 301         pz[nPointZones] = new pointZone
 302         (
 303             "cutPointZone",
 304             labelList(0),
 305             nPointZones,
 306             pointZones()
 307         );
 308
 309         nPointZones++;
 310     }
 311
 312
 313     {
 314         labelList movingCells(nCells());
 315         label nMovingCells = 0;
 316
 317         scalar pistonX = piston().cs().origin().x();
 318         scalar pistonY = piston().cs().origin().y();
 319
 320         forAll(cellCentres(),cellI)
 321         {
 322             const vector& v = cellCentres()[cellI];
 323
 324             if
 325             (
 326                 sqrt(sqr(v.x()-pistonX)+sqr(v.y()-pistonY))
 327              < 0.5*engTime().bore().value()
 328             )
 329             {
 330                 movingCells[nMovingCells] = cellI;
 331                 nMovingCells++;
 332             }
 333         }
 334
 335         movingCells.setSize(nMovingCells);
 336         Info<< "Number of cells in the moving region: "
 337             << nMovingCells << endl;
 338
 339         cz[nCellZones] = new cellZone
 340         (
 341             "movingCells",
 342             movingCells,
 343             nCellZones,
 344             cellZones()
 345         );
 346
 347         nCellZones++;
 348     }
 349
 350
 351     Info<< "Adding " << nPointZones << " point, "
 352         << nFaceZones << " face zones and "
 353         << nCellZones << " cell zones" << endl;
 354
 355     pz.setSize(nPointZones);
 356     fz.setSize(nFaceZones);
 357     cz.setSize(nCellZones);
 358     addZones(pz, fz, cz);
 359
 360     List<polyMeshModifier*> tm(2);
 361     label nMods = 0;
 362
 363     // Add piston layer addition
 364     if (piston().patchID().active())
 365     {
 366
 367         topoChanger_.setSize(topoChanger_.size() + 1);
 368
 369         topoChanger_.set
 370         (
 371             nMods,
 372             new layerAdditionRemoval
 373             (
 374                 "pistonLayer",
 375                 nMods,
 376                 topoChanger_,
 377                 "pistonLayerFaces",
 378                 piston().minLayer(),
 379                 piston().maxLayer()
 380             )
 381         );
 382         nMods++;
 383     }
 384
 385
 386     if(foundScavPorts())
 387     {
 388         topoChanger_.setSize(topoChanger_.size() + 1);
 389
 390         topoChanger_.set
 391         (
 392             nMods,
 393             new slidingInterface
 394             (
 395                 "valveSlider",
 396                 nMods,
 397                 topoChanger_,
 398                 scavInCylPatchName_ + "Zone",
 399                 scavInPortPatchName_ + "Zone",
 400                 "cutPointZone",
 401                 "cutFaceZone",
 402                 scavInCylPatchName_,
 403                 scavInPortPatchName_,
 404                 slidingInterface::INTEGRAL,
 405 //              slidingInterface::PARTIAL,
 406                 true,                          // Attach-detach action
 407                 intersection::VISIBLE         // Projection algorithm
 408             )
 409         );
 410         nMods++;
 411     }
 412
 413     Info << "Adding " << nMods << " topology modifiers" << endl;
 414
 415     // Calculating the virtual piston position
 416
 417     setVirtualPistonPosition();
 418
 419     topoChanger_.setSize(nMods);
 420
 421     // Write mesh modifiers
 422     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
 423     topoChanger_.write();
 424     write();
 425
 426     Info << "virtualPistonPosition = " << virtualPistonPosition() << endl;
 427     Info << "piston position = " << pistonPosition() << endl;
 428 }
 429
 430
 431 // ************************************************************************* //