src/engine/engineTopoChangerMesh/engineValveSliding/engineValveSlidingMove.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2005-2010 Tommaso Lucchini
   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 Class
  26     verticalValvesGambit
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "engineValveSliding.H"
  31 #include "slidingInterface.H"
  32 #include "layerAdditionRemoval.H"
  33 #include "attachDetach.H"
  34 #include "surfaceFields.H"
  35 #include "regionSplit.H"
  36 #include "componentMixedTetPolyPatchVectorField.H"
  37 #include "mapPolyMesh.H"
  38
  39 #include "tetPolyMesh.H"
  40 #include "tetPointFields.H"
  41 #include "elementFields.H"
  42 #include "fixedValueTetPolyPatchFields.H"
  43 #include "slipTetPolyPatchFields.H"
  44
  45 #include "tetFem.H"
  46
  47 #include "symmetryFvPatch.H"
  48 #include "wedgeFvPatch.H"
  49 #include "emptyFvPatch.H"
  50 #include "zeroGradientTetPolyPatchFields.H"
  51 #include "tetDecompositionMotionSolver.H"
  52
  53 #include "fixedValueTetPolyPatchFields.H"
  54 #include "mixedTetPolyPatchFields.H"
  55 #include "slipTetPolyPatchFields.H"
  56 #include "zeroGradientTetPolyPatchFields.H"
  57
  58 #include "zeroGradientFvPatchFields.H"
  59 #include "fvCFD.H"
  60
  61 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  62
  63 void Foam::engineValveSliding::makeLayersLive()
  64 {
  65     // Enable layering
  66     forAll (topoChanger_, modI)
  67     {
  68         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  69         {
  70             topoChanger_[modI].enable();
  71         }
  72         else if (isA<slidingInterface>(topoChanger_[modI]))
  73         {
  74             topoChanger_[modI].disable();
  75         }
  76         else if (isA<attachDetach>(topoChanger_[modI]))
  77         {
  78             topoChanger_[modI].enable();
  79         }
  80         else
  81         {
  82             FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
  83                 << "Don't know what to do with mesh modifier "
  84                 << modI << " of type " << topoChanger_[modI].type()
  85                 << abort(FatalError);
  86         }
  87     }
  88 }
  89
  90 void Foam::engineValveSliding::makeSlidersLive()
  91 {
  92     // Enable sliding interface
  93     forAll (topoChanger_, modI)
  94     {
  95         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  96         {
  97             topoChanger_[modI].disable();
  98         }
  99         else if (isA<slidingInterface>(topoChanger_[modI]))
 100         {
 101             topoChanger_[modI].enable();
 102         }
 103         else if (isA<attachDetach>(topoChanger_[modI]))
 104         {
 105             topoChanger_[modI].enable();
 106         }
 107         else
 108         {
 109             FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
 110                 << "Don't know what to do with mesh modifier "
 111                 << modI << " of type " << topoChanger_[modI].type()
 112                 << abort(FatalError);
 113         }
 114     }
 115
 116 }
 117
 118 bool Foam::engineValveSliding::attached() const
 119 {
 120     const polyTopoChanger& morphs = topoChanger_;
 121
 122     bool result = false;
 123
 124     forAll (morphs, modI)
 125     {
 126         if (typeid(morphs[modI]) == typeid(slidingInterface))
 127         {
 128             result =
 129                 result
 130              || refCast<const slidingInterface>(morphs[modI]).attached();
 131         }
 132     }
 133
 134     // Check thal all sliders are in sync (debug only)
 135     forAll (morphs, modI)
 136     {
 137         if (typeid(morphs[modI]) == typeid(slidingInterface))
 138         {
 139             if
 140             (
 141                 result
 142              != refCast<const slidingInterface>(morphs[modI]).attached()
 143             )
 144             {
 145                 FatalErrorIn("bool movingSquaresTM::attached() const")
 146                     << "Slider " << modI << " named " << morphs[modI].name()
 147                     << " out of sync: Should be" << result
 148                     << abort(FatalError);
 149             }
 150         }
 151     }
 152
 153     return result;
 154 }
 155
 156
 157 void Foam::engineValveSliding::valveDetach()
 158 {
 159     // Enable sliding interface
 160     forAll (topoChanger_, modI)
 161     {
 162         if (isA<attachDetach>(topoChanger_[modI]))
 163         {
 164             const attachDetach& ad =
 165                 refCast<const attachDetach>(topoChanger_[modI]);
 166
 167             const word masterName = ad.masterPatchID().name();
 168
 169             // Find the valve with that name
 170             label valveIndex = -1;
 171
 172             forAll (valves_, valveI)
 173             {
 174                 if
 175                 (
 176                     valves_[valveI].detachInCylinderPatchID().name()
 177                  == masterName
 178                 )
 179                 {
 180                     valveIndex = valveI;
 181                     break;
 182                 }
 183             }
 184
 185             if (valveIndex < 0)
 186             {
 187                 FatalErrorIn
 188                 (
 189                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 190                 )   << "Cannot match patch for attach/detach " << modI
 191                     << abort(FatalError);
 192             }
 193
 194             if (debug)
 195             {
 196                 Info<< " valveI: " << valveIndex << " attached: "
 197                     << ad.attached()
 198                     << " valve open: " << valves_[valveIndex].isOpen()
 199                     << endl;
 200             }
 201
 202             ad.setDetach();
 203         }
 204     }
 205 }
 206
 207 void Foam::engineValveSliding::valveAttach()
 208 {
 209     // Enable sliding interface
 210     forAll (topoChanger_, modI)
 211     {
 212         if (isA<attachDetach>(topoChanger_[modI]))
 213         {
 214             const attachDetach& ad =
 215                 refCast<const attachDetach>(topoChanger_[modI]);
 216
 217             const word masterName = ad.masterPatchID().name();
 218
 219             // Find the valve with that name
 220             label valveIndex = -1;
 221
 222             forAll (valves_, valveI)
 223             {
 224                 if
 225                 (
 226                     valves_[valveI].detachInCylinderPatchID().name()
 227                  == masterName
 228                 )
 229                 {
 230                     valveIndex = valveI;
 231                     break;
 232                 }
 233             }
 234
 235             if (valveIndex < 0)
 236             {
 237                 FatalErrorIn
 238                 (
 239                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 240                 )   << "Cannot match patch for attach/detach " << modI
 241                     << abort(FatalError);
 242             }
 243
 244             if (debug)
 245             {
 246                 Info<< " valveI: " << valveIndex << " attached: "
 247                     << ad.attached()
 248                     << " valve open: " << valves_[valveIndex].isOpen()
 249                     << endl;
 250             }
 251
 252             ad.setAttach();
 253         }
 254     }
 255 }
 256
 257
 258 void Foam::engineValveSliding::prepareValveDetach()
 259 {
 260     // Enable sliding interface
 261     forAll (topoChanger_, modI)
 262     {
 263         if (isA<attachDetach>(topoChanger_[modI]))
 264         {
 265             const attachDetach& ad =
 266                 refCast<const attachDetach>(topoChanger_[modI]);
 267
 268             const word masterName = ad.masterPatchID().name();
 269
 270             // Find the valve with that name
 271             label valveIndex = -1;
 272
 273             forAll (valves_, valveI)
 274             {
 275                 if
 276                 (
 277                     valves_[valveI].detachInCylinderPatchID().name()
 278                  == masterName
 279                 )
 280                 {
 281                     valveIndex = valveI;
 282                     break;
 283                 }
 284             }
 285
 286             if (valveIndex < 0)
 287             {
 288                 FatalErrorIn
 289                 (
 290                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 291                 )   << "Cannot match patch for attach/detach " << modI
 292                     << abort(FatalError);
 293             }
 294
 295             if (debug)
 296             {
 297                 Info<< " valveI: " << valveIndex << " attached: "
 298                     << ad.attached()
 299                     << " valve open: " << valves_[valveIndex].isOpen()
 300                     << endl;
 301             }
 302
 303             if (valves_[valveIndex].isOpen())
 304             {
 305                 ad.setAttach();
 306             }
 307             else
 308             {
 309                 ad.setDetach();
 310             }
 311         }
 312     }
 313 }
 314
 315
 316 bool Foam::engineValveSliding::update()
 317 {
 318
 319     tetDecompositionMotionSolver& mSolver =
 320         refCast<tetDecompositionMotionSolver>(msPtr_());
 321
 322     // Detaching the interfacethobois2DSlidingDeform
 323     if (attached())
 324     {
 325         Info << "Decoupling sliding interfaces" << endl;
 326         makeSlidersLive();
 327         valveDetach();
 328         autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 329
 330         Info << "sliding interfaces successfully decoupled!!!" << endl;
 331         if (topoChangeMap1->morphing())
 332         {
 333             mSolver.updateMesh(topoChangeMap1());
 334         }
 335     }
 336     else
 337     {
 338         Info << "Sliding interfaces decoupled" << endl;
 339         valveDetach();
 340     }
 341
 342     Info << "Executing layer action" << endl;
 343
 344     scalar deltaZ = engTime().pistonDisplacement().value();
 345
 346     // Piston Layering
 347     makeLayersLive();
 348
 349     {
 350         forAll(valves_, valveI)
 351         {
 352             scalar valveDisplacement = valves_[valveI].curVelocity()*
 353                 valves_[valveI].cs().axis().z()*engTime().deltaT().value();
 354
 355             Info << "valveDisplacement = " << valveDisplacement << nl
 356                 << "valve velocity =" << valves_[valveI].curVelocity() << nl
 357                 << "valve lift =" << valves_[valveI].curLift() << nl
 358                 << "valve deformation lift ="
 359                 << valves_[valveI].deformationLift() << endl;
 360         }
 361
 362     }
 363
 364     forAll(valves_,valveI)
 365     {
 366         if(valves_[valveI].curLift() >= valves_[valveI].deformationLift())
 367         {
 368             if(valves_[valveI].poppetPatchID().active())
 369             {
 370                 // Find valve layer mesh modifier
 371                 const label valveLayerID =
 372                     topoChanger_.findModifierID
 373                     (
 374                         "valvePoppetLayer"
 375                       + Foam::name(valveI + 1)
 376                     );
 377
 378                 if (valveLayerID < 0)
 379                 {
 380                         FatalErrorIn("void engineFvMesh::moveAndMorph()")
 381                            << "valve modifier not found."
 382                            << abort(FatalError);
 383                 }
 384
 385                 topoChanger_[valveLayerID].enable();
 386             }
 387         }
 388         else
 389         {
 390             if(valves_[valveI].poppetPatchID().active())
 391             {
 392                 // Find valve layer mesh modifier
 393                 const label valveLayerID =
 394                     topoChanger_.findModifierID
 395                     (
 396                         "valvePoppetLayer"
 397                       + Foam::name(valveI + 1)
 398                     );
 399
 400                 if (valveLayerID < 0)
 401                 {
 402                         FatalErrorIn("void engineFvMesh::moveAndMorph()")
 403                            << "valve modifier not found."
 404                            << abort(FatalError);
 405                 }
 406
 407                 topoChanger_[valveLayerID].disable();
 408             }
 409         }
 410     }
 411
 412     pointField newPoints = points();
 413
 414     autoPtr<mapPolyMesh> topoChangeMap2 = topoChanger_.changeMesh();
 415
 416     // Changing topology by hand
 417     if (topoChangeMap2->morphing())
 418     {
 419         mSolver.updateMesh(topoChangeMap2());
 420
 421         if (topoChangeMap2->hasMotionPoints())
 422         {
 423             movePoints(topoChangeMap2->preMotionPoints());
 424             newPoints = topoChangeMap2->preMotionPoints();
 425         }
 426         setV0();
 427         resetMotion();
 428     }
 429
 430     // Work array for new points position.
 431
 432 // NEW !
 433
 434     // Reset the position of layered interfaces
 435
 436     {
 437
 438
 439 #        include "setMotionBoundaryConditionEngineValveSliding.H"
 440
 441         DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
 442         DynamicList<vector> constrainedVelocity
 443         (
 444             mSolver.curPoints()().size()/100
 445         );
 446
 447 #       include "setEngineValveSlidingConstraint.H"
 448
 449         forAll (constrainedPoints, i)
 450         {
 451             mSolver.setConstraint
 452             (
 453                 constrainedPoints[i],
 454                 constrainedVelocity[i]
 455             );
 456         }
 457
 458 //        Info << mSolver.motionU() << endl;
 459
 460         Info << "pre solve" << endl;
 461
 462         mSolver.solve();
 463
 464         newPoints = mSolver.curPoints();
 465         movePoints(newPoints);
 466
 467         Info << "max mesh phi, 1 = " << max(phi()) << endl;
 468         Info << "min mesh phi, 1 = " << min(phi()) << endl;
 469
 470         setVirtualPositions();
 471         mSolver.clearConstraints();
 472
 473
 474
 475 //      layering
 476
 477 #       include "moveValvePointsEngineValveSliding.H"
 478         movePoints(newPoints);
 479 //        newPoints = points();
 480         newPoints = allPoints();
 481         setVirtualPositions();
 482
 483         Info << "max mesh phi, 2 = " << max(phi()) << endl;
 484         Info << "min mesh phi, 2 = " << min(phi()) << endl;
 485
 486     }
 487
 488     Info << ", deckHeight = " << deckHeight()
 489         << ", pistonPosition = " << pistonPosition()    << endl;
 490
 491     pistonPosition() += deltaZ;
 492
 493     Info << "max V()-V0() before = " << max(mag(V() - V0())) << endl;
 494     Info << "max V() before = " << max(mag(V())) << endl;
 495
 496     Info<< "BEFORE mag(points - oldPoints)/deltaT = "
 497         << max(mag(points() - oldPoints()))/engTime().deltaT().value() << endl;
 498
 499 /*
 500     forAll(points(), pp)
 501     {
 502         Info<< "OLD " <<  "point " << pp << " " << points()[pp]
 503             << " old point = " << oldPoints()[pp] << mag(points()[pp] -
 504         oldPoints()[pp])/engTime().deltaT().value() << endl;
 505     }
 506
 507 */
 508
 509 //*/ //Tommaso, 23/5/2008
 510
 511     {
 512         // Grab old points to correct the motion
 513         pointField oldPointsNew = oldAllPoints();
 514 //        pointField oldPointsNew = oldPoints();
 515
 516         // Attach the interface
 517         Info << "Coupling sliding interfaces" << endl;
 518         makeSlidersLive();
 519
 520         // Changing topology by hand
 521         autoPtr<mapPolyMesh> topoChangeMap4 = topoChanger_.changeMesh();
 522
 523         // Changing topology by hand
 524         if(topoChangeMap4->morphing())
 525         {
 526             mSolver.updateMesh(topoChangeMap4());
 527
 528             if (topoChangeMap4->hasMotionPoints())
 529             {
 530 //                 Info<< "Topology change; executing pre-motion "
 531 //                     << "after sliding attach" << endl;
 532
 533 //                 Info<< "topoChangeMap4->preMotionPoints().size() = "
 534 //                     << topoChangeMap4->preMotionPoints().size() << endl;
 535 //                 Info << "allPoints.size() = " << allPoints().size() << endl;
 536 //                 Info << "points.size() = " << points().size() << endl;
 537
 538                 movePoints(topoChangeMap4->preMotionPoints());
 539 //            newPoints = topoChangeMap4->preMotionPoints();
 540 //            newPoints = allPoints();
 541                 newPoints = allPoints();
 542             }
 543         }
 544
 545         surfaceScalarField correctedMeshPhi = phi();
 546
 547
 548         Info << "max V()-V0() after = " << max(mag(V() - V0())) << endl;
 549         Info << "max V() before = " << max(mag(V())) << endl;
 550
 551         {
 552             Info << "AFTER mag(points - oldPoints)/deltaT = "
 553                 << max(mag(points() - oldPoints()))/engTime().deltaT().value()
 554                 << endl;
 555             Info << "max V()-V0() after movePoints?!? = "
 556                 << max(mag(V() - V0())) << endl;
 557
 558             pointField mappedOldPointsNew(allPoints().size());
 559
 560             mappedOldPointsNew.map(oldPointsNew, topoChangeMap4->pointMap());
 561
 562             forAll(valves(), valveI)
 563             {
 564                 const labelList& cutPointsAddressing =
 565                     pointZones()
 566                     [
 567                         pointZones().findZoneID
 568                         (
 569                             "cutPointsV"
 570                             + Foam::name(valveI + 1)
 571                         )
 572                     ];
 573
 574                 forAll(cutPointsAddressing, i)
 575                 {
 576                     mappedOldPointsNew[cutPointsAddressing[i]] =
 577                         newPoints[cutPointsAddressing[i]];
 578                 }
 579             }
 580
 581             pointField newPoints = allPoints();
 582
 583             movePoints(mappedOldPointsNew);
 584
 585             resetMotion();
 586             setV0();
 587             movePoints(newPoints);
 588
 589             Info<< "max mesh phi, sliding mesh attached = "
 590                 << max(mag(phi().internalField())) << endl;
 591             Info<< "max mesh phi, sliding mesh attached = "
 592                 << min(mag(phi().internalField())) << endl;
 593
 594             Info<< "AFTER 2 mag(points - oldPoints)/deltaT = "
 595                 << max(mag(points() - oldPoints()))/engTime().deltaT().value()
 596                 << endl;
 597             Info << "AFTER 2 max V()-V0() after = "
 598                 << max(mag(V() - V0())) << endl;
 599
 600         }
 601     }
 602
 603     return true;
 604 }
 605
 606