src/engine/simpleEngineTopoFvMesh/simpleEngineTopoFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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 "simpleEngineTopoFvMesh.H"
  27 #include "slidingInterface.H"
  28 #include "layerAdditionRemoval.H"
  29 #include "attachDetach.H"
  30 #include "componentMixedTetPolyPatchVectorField.H"
  31 #include "mapPolyMesh.H"
  32 #include "polyTopoChange.H"
  33 #include "tetMotionSolver.H"
  34 #include "volMesh.H"
  35 #include "addToRunTimeSelectionTable.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41     defineTypeNameAndDebug(simpleEngineTopoFvMesh, 0);
  42
  43     addToRunTimeSelectionTable
  44     (
  45         engineTopoChangerMesh,
  46         simpleEngineTopoFvMesh,
  47         IOobject
  48     );
  49 }
  50
  51
  52 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  53
  54 void Foam::simpleEngineTopoFvMesh::makeLayersLive()
  55 {
  56     // Enable layering
  57     forAll (topoChanger_, modI)
  58     {
  59         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  60         {
  61             if (debug)
  62             {
  63                 Info<< "Enabling layer modifier "
  64                     << topoChanger_[modI].name() << endl;
  65             }
  66
  67             topoChanger_[modI].enable();
  68         }
  69         else if (isA<slidingInterface>(topoChanger_[modI]))
  70         {
  71             if (debug)
  72             {
  73                 Info<< "Disabling slider modifier "
  74                     << topoChanger_[modI].name() << endl;
  75             }
  76
  77             topoChanger_[modI].disable();
  78         }
  79         else if (isA<attachDetach>(topoChanger_[modI]))
  80         {
  81             topoChanger_[modI].enable();
  82         }
  83         else
  84         {
  85             FatalErrorIn("void Foam::simpleEngineTopoFvMesh::makeLayersLive()")
  86                 << "Don't know what to do with mesh modifier "
  87                 << modI << " of type " << topoChanger_[modI].type()
  88                 << abort(FatalError);
  89         }
  90     }
  91 }
  92
  93
  94 void Foam::simpleEngineTopoFvMesh::makeSlidersLive()
  95 {
  96     // Enable sliding interface
  97     forAll (topoChanger_, modI)
  98     {
  99         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
 100         {
 101             if (debug)
 102             {
 103                 Info<< "Disabling layer modifier "
 104                     << topoChanger_[modI].name() << endl;
 105             }
 106
 107             topoChanger_[modI].disable();
 108         }
 109         else if (isA<slidingInterface>(topoChanger_[modI]))
 110         {
 111             if (debug)
 112             {
 113                 Info<< "Enabling slider modifier "
 114                     << topoChanger_[modI].name() << endl;
 115             }
 116
 117             topoChanger_[modI].enable();
 118         }
 119         else if (isA<attachDetach>(topoChanger_[modI]))
 120         {
 121             topoChanger_[modI].enable();
 122         }
 123         else
 124         {
 125             FatalErrorIn("void Foam::simpleEngineTopoFvMesh::makeSlidersLive()")
 126                 << "Don't know what to do with mesh modifier "
 127                 << modI << " of type " << topoChanger_[modI].type()
 128                 << abort(FatalError);
 129         }
 130     }
 131 }
 132
 133
 134 void Foam::simpleEngineTopoFvMesh::prepareValveDetach()
 135 {
 136     // Enable sliding interface
 137     forAll (topoChanger_, modI)
 138     {
 139         if (isA<attachDetach>(topoChanger_[modI]))
 140         {
 141             const attachDetach& ad =
 142                 refCast<const attachDetach>(topoChanger_[modI]);
 143
 144             const word masterName = ad.masterPatchID().name();
 145
 146             // Find the valve with that name
 147             label valveIndex = -1;
 148
 149             forAll (valves_, valveI)
 150             {
 151                 if
 152                 (
 153                     valves_[valveI].detachInCylinderPatchID().name()
 154                  == masterName
 155                 )
 156                 {
 157                     valveIndex = valveI;
 158                     break;
 159                 }
 160             }
 161
 162             if (valveIndex < 0)
 163             {
 164                 FatalErrorIn
 165                 (
 166                     "void Foam::simpleEngineTopoFvMesh::prepareValveDetach()"
 167                 )   << "Cannot match patch for attach/detach " << modI
 168                     << abort(FatalError);
 169             }
 170
 171             if (debug)
 172             {
 173                 Info<< " valveI: " << valveIndex << " attached: "
 174                     << ad.attached()
 175                     << " valve open: " << valves_[valveIndex].isOpen()
 176                     << endl;
 177             }
 178
 179             if (valves_[valveIndex].isOpen())
 180             {
 181                 ad.setAttach();
 182             }
 183             else
 184             {
 185                 ad.setDetach();
 186             }
 187         }
 188     }
 189 }
 190
 191
 192 bool Foam::simpleEngineTopoFvMesh::attached() const
 193 {
 194     bool result = false;
 195
 196     forAll (topoChanger_, modI)
 197     {
 198         if (isA<slidingInterface>(topoChanger_[modI]))
 199         {
 200             result =
 201                 result
 202              || refCast<const slidingInterface>(topoChanger_[modI]).attached();
 203         }
 204     }
 205
 206     // Check thal all sliders are in sync (debug only)
 207     forAll (topoChanger_, modI)
 208     {
 209         if (isA<slidingInterface>(topoChanger_[modI]))
 210         {
 211             if
 212             (
 213                 result
 214              != refCast<const slidingInterface>(topoChanger_[modI]).attached()
 215             )
 216             {
 217                 FatalErrorIn("bool simpleEngineTopoFvMesh::attached() const")
 218                     << "Slider " << modI << " named "
 219                     << topoChanger_[modI].name()
 220                     << " out of sync: Should be" << result
 221                     << abort(FatalError);
 222             }
 223         }
 224     }
 225
 226     if (debug)
 227     {
 228         if (result)
 229         {
 230             Info << "simpleEngineTopoFvMesh is attached" << endl;
 231         }
 232         else
 233         {
 234             Info << "simpleEngineTopoFvMesh is detached" << endl;
 235         }
 236     }
 237
 238     return result;
 239 }
 240
 241
 242 void Foam::simpleEngineTopoFvMesh::setBoundaryMotion()
 243 {
 244     // Set the boundary conditions on motion field in order to solve
 245     // for motion.  Deformation only happens within the cylinder and
 246     // not in ports - the motion of valve top is set to zero.  Correct
 247     // using setBoundaryPosition()
 248     // HJ, 10/Jun/2004  Reconsider
 249     if (debug)
 250     {
 251         Info << "Setting boundary motion" << endl;
 252     }
 253
 254     tetMotionSolver& mSolver =
 255         refCast<tetMotionSolver>(msPtr_());
 256
 257     tetPointVectorField& motionU = mSolver.motionU();
 258
 259     // Set valve velocity
 260     forAll (valves_, valveI)
 261     {
 262         // If valve is present in geometry, set the motion
 263         if (valves_[valveI].bottomPatchID().active())
 264         {
 265             vector valveVel =
 266                 valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
 267
 268             // Bottom of the valve moves with given velocity
 269             motionU.boundaryField()[valves_[valveI].bottomPatchID().index()] ==
 270                 valveVel;
 271
 272             if (debug)
 273             {
 274                 Info<< "Valve " << valveI << " lift: "
 275                     << valves_[valveI].curLift()
 276                     << " velocity: " << valves_[valveI].curVelocity()
 277                     << endl;
 278             }
 279         }
 280
 281         if (valves_[valveI].poppetPatchID().active())
 282         {
 283             // Top of the valve does not move
 284             motionU.boundaryField()[valves_[valveI].poppetPatchID().index()] ==
 285                 vector::zero;
 286         }
 287
 288         if (valves_[valveI].curtainInCylinderPatchID().active())
 289         {
 290 //             label cicPatchIndex =
 291 //                 valves_[valveI].curtainInCylinderPatchID().index();
 292
 293 //             componentMixedTetPolyPatchVectorField& pf =
 294 //                 refCast<componentMixedTetPolyPatchVectorField>
 295 //                 (
 296 //                     motionU.boundaryField()[cicPatchIndex]
 297 //                 );
 298
 299 //             if (valves_[valveI].isOpen())
 300 //             {
 301 //                 // Get valve coordinate system
 302 //                 const coordinateSystem& vcs = valves_[valveI].cs();
 303 //                 const scalar r = 0.5*valves_[valveI].diameter();
 304
 305 //                 // Get local points in the patch
 306 //                 const pointField& cpGlobal =
 307 //                 motionU.boundaryField()
 308 //                     [cicPatchIndex].patchMesh().localPoints();
 309
 310 //                 pointField cpLocal(vcs.toLocal(cpGlobal));
 311 //                 scalarField mask =
 312 //                     pos
 313 //                     (
 314 //                         cpLocal.component(vector::Z)
 315 //                       + valves_[valveI].curLift()
 316 //                       - valvePosTol_
 317 //                     );
 318
 319 //                 // Calculate motion of valve centre
 320 //                 cpLocal.replace
 321 //                 (
 322 //                     vector::X,
 323 //                     mask*r + (1.0 - mask)*0.0
 324 //                 );
 325
 326 //                 pf.refValue() =
 327 //                 (
 328 //                     mask*(vcs.toGlobal(cpLocal) - cpGlobal)/
 329 //                     engineTime_.deltaT().value()
 330 //                   + (1.0 - mask)*
 331 //                     vector(vcs.axis().x()*valves_[valveI].curVelocity(), 0, 0)
 332 //                 );
 333
 334 //                 pf.valueFraction() =
 335 //                     mask*vector::one + (1.0 - mask)*vector(1, 0, 0);
 336 //             }
 337 //             else
 338 //             {
 339 //                 pf.refValue() = vector::zero;
 340 //                 pf.valueFraction() = vector::one;
 341 //             }
 342         }
 343     }
 344
 345     // Set piston velocity
 346     if (piston().patchID().active())
 347     {
 348         vector pistonVel =
 349             piston().cs().axis()*engineTime_.pistonSpeed().value();
 350
 351         if (debug)
 352         {
 353             Info<< "Piston velocity: " << pistonVel;
 354         }
 355
 356         componentMixedTetPolyPatchVectorField& pp =
 357             refCast<componentMixedTetPolyPatchVectorField>
 358             (
 359                 motionU.boundaryField()[piston().patchID().index()]
 360             );
 361
 362         if (deformation())
 363         {
 364             if (debug)
 365             {
 366                 Info << " deformation"  << endl;
 367             }
 368
 369             pp.refValue() = pistonVel;
 370         }
 371         else
 372         {
 373             if (debug)
 374             {
 375                 Info << " layering" << endl;
 376             }
 377
 378             pp.refValue() = vector::zero;
 379         }
 380     }
 381 }
 382
 383
 384 void Foam::simpleEngineTopoFvMesh::setBoundaryPosition()
 385 {
 386     // Set the boundary position for layer modifiers
 387     if (debug)
 388     {
 389         Info << "Setting boundary position" << endl;
 390     }
 391
 392     tetMotionSolver& mSolver =
 393         refCast<tetMotionSolver>(msPtr_());
 394
 395     tetPointVectorField& motionU = mSolver.motionU();
 396
 397     // Set valve velocity
 398     forAll (valves_, valveI)
 399     {
 400         vector valveVel =
 401             valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
 402
 403         // If valve is present in geometry, set the motion
 404         if (valves_[valveI].bottomPatchID().active())
 405         {
 406             // Bottom of the valve moves with given velocity
 407             motionU.boundaryField()[valves_[valveI].bottomPatchID().index()] ==
 408                 valveVel;
 409
 410             if (debug)
 411             {
 412                 Info<< "Valve " << valveI << " lift: "
 413                     << valves_[valveI].curLift()
 414                     << " velocity: " << valves_[valveI].curVelocity()
 415                     << endl;
 416             }
 417
 418         }
 419
 420         if (valves_[valveI].poppetPatchID().active())
 421         {
 422             // Top of the valve does not move
 423             motionU.boundaryField()[valves_[valveI].poppetPatchID().index()] ==
 424                 valveVel;
 425         }
 426     }
 427
 428     // Set piston velocity
 429     if (piston().patchID().active())
 430     {
 431         vector pistonVel =
 432             piston().cs().axis()*engineTime_.pistonSpeed().value();
 433
 434         componentMixedTetPolyPatchVectorField& pp =
 435             refCast<componentMixedTetPolyPatchVectorField>
 436             (
 437                 motionU.boundaryField()[piston().patchID().index()]
 438             );
 439
 440         if (debug)
 441         {
 442             Info<< "Piston velocity: " << pistonVel << endl;
 443         }
 444
 445         pp.refValue() = pistonVel;
 446     }
 447
 448     motionU.correctBoundaryConditions();
 449 }
 450
 451
 452 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 453
 454 // Construct from components
 455 Foam::simpleEngineTopoFvMesh::simpleEngineTopoFvMesh
 456 (
 457     const IOobject& io
 458 )
 459 :
 460     engineTopoChangerMesh(io),
 461     valves_(*this, engineTime_.engineDict().lookup("valves")),
 462     piston_(*this, engineTime_.engineDict().subDict("piston")),
 463     msPtr_(motionSolver::New(*this)),
 464     deformSwitch_(readScalar(engineTime_.engineDict().lookup("deformAngle"))),
 465     valvePosTol_(readScalar(engineTime_.engineDict().lookup("valvePosTol")))
 466 {
 467     // Add zones and modifiers if not already there.
 468     addZonesAndModifiers();
 469 }
 470
 471
 472 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 473
 474
 475 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 476
 477 bool Foam::simpleEngineTopoFvMesh::update()
 478 {
 479     tetMotionSolver& mSolver =
 480         refCast<tetMotionSolver>(msPtr_());
 481
 482     // Detaching the interface
 483     if (attached())
 484     {
 485         if (debug)
 486         {
 487             Info << "Decoupling sliding interfaces" << endl;
 488         }
 489
 490         makeSlidersLive();
 491
 492         // Changing topology by hand
 493         autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 494
 495         if (topoChangeMap1->morphing())
 496         {
 497             mSolver.updateMesh(topoChangeMap1());
 498         }
 499     }
 500     else
 501     {
 502         if (debug)
 503         {
 504             Info << "Sliding interfaces decoupled" << endl;
 505         }
 506     }
 507
 508     // Perform layer action and mesh motion
 509     makeLayersLive();
 510
 511     if (debug)
 512     {
 513         Info << "Executing layer action" << endl;
 514     }
 515
 516     // Find piston mesh modifier
 517     const label pistonLayerID =
 518         topoChanger_.findModifierID("pistonLayer");
 519
 520     if (pistonLayerID < 0)
 521     {
 522         FatalErrorIn("void simpleEngineTopoFvMesh::moveAndMorph()")
 523             << "Piston modifier not found."
 524             << abort(FatalError);
 525     }
 526
 527     if (deformation())
 528     {
 529         // Dectivate piston layer
 530         if (debug)
 531         {
 532             Info << "Disabling piston layer (deformation)"<< endl;
 533         }
 534
 535         topoChanger_[pistonLayerID].disable();
 536     }
 537     else
 538     {
 539         // Activate piston layer
 540         if (debug)
 541         {
 542             Info << "Enabling piston layer (deformation)"<< endl;
 543         }
 544
 545         topoChanger_[pistonLayerID].enable();
 546     }
 547
 548     // Changing topology by hand
 549     {
 550         autoPtr<mapPolyMesh> topoChangeMap2 = topoChanger_.changeMesh();
 551
 552         if (topoChangeMap2->morphing())
 553         {
 554             mSolver.updateMesh(topoChangeMap2());
 555
 556             if (debug)
 557             {
 558                 Info << "Topology change; executing pre-motion" << endl;
 559             }
 560
 561             movePoints(topoChangeMap2->preMotionPoints());
 562             setV0();
 563             resetMotion();
 564
 565         }
 566     }
 567
 568     if (deformation())
 569     {
 570         if (debug)
 571         {
 572             Info << "Mesh deformation mode" << endl;
 573         }
 574
 575         setBoundaryMotion();
 576
 577         // Solve for motion
 578         mSolver.solve();
 579
 580         // Dectivate piston layer
 581         if (debug)
 582         {
 583             Info << "Disabling piston layer (topo 2)"<< endl;
 584         }
 585
 586         topoChanger_[pistonLayerID].disable();
 587     }
 588     else
 589     {
 590         if (debug)
 591         {
 592             Info << "Piston layering mode" << endl;
 593         }
 594
 595         bool tiltedValves = true;
 596
 597         if (tiltedValves)
 598         {
 599             setBoundaryMotion();
 600
 601             // Solve for motion
 602             mSolver.solve();
 603         }
 604
 605         // Blocking vertical motion
 606         mSolver.motionU().internalField().replace(vector::Z, 0);
 607
 608         // Activate piston layer
 609         if (debug)
 610         {
 611             Info << "Enabling piston layer (topo 2)"<< endl;
 612         }
 613
 614         topoChanger_[pistonLayerID].enable();
 615     }
 616
 617     // Reset the position of layered interfaces
 618     setBoundaryPosition();
 619
 620     movePoints(mSolver.curPoints());
 621
 622     // Attach the interface
 623     if (debug)
 624     {
 625         Info << "Coupling sliding interfaces" << endl;
 626     }
 627
 628     makeSlidersLive();
 629     prepareValveDetach();
 630
 631     // Changing topology by hand
 632     {
 633         // Grab old points to correct the motion
 634         pointField oldPointsNew = oldAllPoints();
 635
 636         autoPtr<mapPolyMesh> topoChangeMap3 = topoChanger_.changeMesh();
 637
 638         if (debug)
 639         {
 640             Info << "Moving points post slider attach" << endl;
 641         }
 642
 643         if (topoChangeMap3->morphing())
 644         {
 645             mSolver.updateMesh(topoChangeMap3());
 646
 647             if (debug)
 648             {
 649                 Info << "Moving points post slider attach" << endl;
 650             }
 651
 652             pointField newPoints = allPoints();
 653             pointField mappedOldPointsNew(newPoints.size());
 654
 655             mappedOldPointsNew.map(oldPointsNew, topoChangeMap3->pointMap());
 656
 657             // Solve the correct mesh motion to make sure motion fluxes
 658             // are solved for and not mapped
 659             movePoints(mappedOldPointsNew);
 660             resetMotion();
 661             setV0();
 662             movePoints(newPoints);
 663         }
 664     }
 665
 666     if (debug)
 667     {
 668         Info << "Sliding interfaces coupled: " << attached() << endl;
 669     }
 670
 671     return true;
 672 }
 673
 674
 675 // ************************************************************************* //