src/solidModels/contactModels/frictionContactModels/iterativePenaltyFriction/iterativePenaltyFriction.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 "iterativePenaltyFriction.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "zeroGradientFvPatchFields.H"
  29 #include "PatchToPatchInterpolationTemplate.H"
  30 #include "ggiInterpolation.H"
  31 #include "constitutiveModel.H"
  32 #include "fvc.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(iterativePenaltyFriction, 0);
  39     addToRunTimeSelectionTable
  40     (frictionContactModel, iterativePenaltyFriction, dictionary);
  41 }
  42
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // Construct from dictionary
  50 Foam::iterativePenaltyFriction::iterativePenaltyFriction
  51 (
  52     const word& name,
  53     const fvPatch& patch,
  54     const dictionary& dict,
  55     const label masterPatchID,
  56     const label slavePatchID,
  57     const label masterFaceZoneID,
  58     const label slaveFaceZoneID
  59 )
  60 :
  61   frictionContactModel
  62   (
  63       name,
  64       patch,
  65       dict,
  66       masterPatchID,
  67       slavePatchID,
  68       masterFaceZoneID,
  69       slaveFaceZoneID
  70   ),
  71   frictionContactModelDict_(dict.subDict(name+"FrictionModelDict")),
  72   frictionLawPtr_(NULL),
  73   mesh_(patch.boundaryMesh().mesh()),
  74   slaveDisp_(mesh().boundaryMesh()[slavePatchID].size(), vector::zero),
  75   slaveTraction_(mesh().boundaryMesh()[slavePatchID].size(), vector::zero),
  76   slaveValueFrac_(mesh().boundaryMesh()[slavePatchID].size(), symmTensor::zero),
  77   frictionPenaltyFactorPtr_(NULL),
  78   frictionPenaltyScale_
  79   (readScalar(frictionContactModelDict_.lookup("penaltyScale"))),
  80   contactIterNum_(0),
  81   infoFreq_(readInt(frictionContactModelDict_.lookup("infoFrequency"))),
  82   oscillationCorr_(frictionContactModelDict_.lookup("oscillationCorrection")),
  83   oscillationCorrFac_
  84   (readScalar(frictionContactModelDict_.lookup("oscillationCorrectionFactor"))),
  85   contactFilePtr_(NULL)
  86 {
  87   // create friction law
  88   frictionLawPtr_ = frictionLaw::New(
  89                      frictionContactModelDict_.lookup("frictionLaw"),
  90                      frictionContactModelDict_
  91                      ).ptr();
  92
  93   // master proc open contact info file
  94   if (Pstream::master())
  95     {
  96       word masterName = mesh_.boundary()[masterPatchID].name();
  97       word slaveName = mesh_.boundary()[slavePatchID].name();
  98       contactFilePtr_ =
  99           new OFstream
 100           (fileName("frictionContact_"+masterName+"_"+slaveName+".txt"));
 101       OFstream& contactFile = *contactFilePtr_;
 102       int width = 20;
 103       contactFile << "time";
 104       contactFile.width(width);
 105       contactFile << "iterNum";
 106       contactFile.width(width);
 107       contactFile << "penaltyScale";
 108       contactFile.width(width);
 109       contactFile << "slipFaces";
 110       contactFile.width(width);
 111       contactFile << "stickFaces";
 112       contactFile.width(width);
 113       contactFile << "maxMagSlaveTraction" << endl;
 114     }
 115 }
 116
 117
 118 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 119
 120   void Foam::iterativePenaltyFriction::correct
 121   (
 122    const vectorField& slavePressure,
 123    const PrimitivePatch<face, List, pointField>& masterFaceZonePatch,
 124    const PrimitivePatch<face, List, pointField>& slaveFaceZonePatch,
 125    const intersection::algorithm alg,
 126    const intersection::direction dir,
 127    const word interpolationMethod,
 128    const word fieldName,
 129    const Switch orthotropic,
 130    const word nonLinear,
 131    const vectorField& slaveFaceNormals
 132    )
 133   {
 134     const fvMesh& mesh = mesh_;
 135     const label slavePatchIndex = slavePatchID();
 136     const label masterPatchIndex = masterPatchID();
 137     contactIterNum_++;
 138
 139     // we have local masterDU and we want to interpolate it to the slave
 140     // to get local masterDUInterpToSlave (i.e. masterDU interpolated to
 141     // the slave)
 142     // so the method is:
 143     // create global masterDU field
 144     // interpolate global masterDU from master global face zone to slave
 145     // global zone
 146     // then find local masterDUInterpToSlave from the global interpolated field
 147
 148     vectorField masterDUInterpToSlave
 149         (mesh.boundaryMesh()[slavePatchIndex].size(), vector::zero);
 150
 151     // global master DU
 152     vectorField globalMasterDU(masterFaceZonePatch.size(), vector::zero);
 153
 154     // lookup current displacement field
 155     const volVectorField& dispField =
 156       mesh.objectRegistry::lookupObject<volVectorField>(fieldName);
 157
 158     // local master and slave DU increment
 159     vectorField masterDU = dispField.boundaryField()[masterPatchIndex];
 160     vectorField slaveDU = dispField.boundaryField()[slavePatchIndex];
 161
 162     if (fieldName == "U")
 163       {
 164     // lookup old U
 165     const volVectorField& dispOldField =
 166       mesh.objectRegistry::lookupObject<volVectorField>(fieldName+"_0");
 167
 168     // subtract old U
 169     masterDU -= dispOldField.boundaryField()[masterPatchIndex];
 170     slaveDU -= dispOldField.boundaryField()[slavePatchIndex];
 171       }
 172     else if (fieldName != "DU")
 173       {
 174     FatalError << "iterativePenaltyFunction::correct()\n"
 175       " The displacement field must be called U or DU"
 176            << exit(FatalError);
 177       }
 178
 179     // put local masterDU into globalMasterDU
 180     const label masterPatchStart
 181       = mesh.boundaryMesh()[masterPatchIndex].start();
 182     forAll(masterDU, i)
 183     {
 184         globalMasterDU
 185             [
 186                 mesh.faceZones()[masterFaceZoneID()
 187                     ].whichFace(masterPatchStart + i)] =
 188             masterDU[i];
 189     }
 190     //- exchange parallel data
 191     // sum because each face is only on one proc
 192     reduce(globalMasterDU, sumOp<vectorField>());
 193
 194     // globalMasterDU is interpolated to the slave
 195     vectorField globalMasterDUInterpToSlave
 196         (slaveFaceZonePatch.size(), vector::zero);
 197
 198     // interpolate DU from master to slave using inverse distance or ggi
 199     if (interpolationMethod == "patchToPatch")
 200       {
 201     PatchToPatchInterpolation<
 202       PrimitivePatch<
 203           face, List, pointField
 204           >, PrimitivePatch<face, List, pointField>
 205       > masterToSlavePatchToPatchInterpolator
 206       (
 207        masterFaceZonePatch, // from zone
 208        slaveFaceZonePatch, // to zone
 209        alg,
 210        dir
 211        );
 212     globalMasterDUInterpToSlave =
 213       masterToSlavePatchToPatchInterpolator.faceInterpolate<vector>
 214       (
 215        globalMasterDU
 216        );
 217       }
 218     else if (interpolationMethod == "ggi")
 219       {
 220     GGIInterpolation<
 221       PrimitivePatch<
 222           face, List, pointField
 223           >, PrimitivePatch< face, List, pointField >
 224       > masterToSlaveGgiInterpolator
 225       (
 226        masterFaceZonePatch, // master zone
 227        slaveFaceZonePatch, // slave zone
 228        tensorField(0),
 229        tensorField(0),
 230        vectorField(0),
 231        0.0,
 232        0.0,
 233        true,
 234        ggiInterpolation::AABB
 235        );
 236     globalMasterDUInterpToSlave =
 237       masterToSlaveGgiInterpolator.masterToSlave
 238       (
 239        globalMasterDU
 240        );
 241       }
 242     else
 243       {
 244     FatalError << "iterativePenaltyFunction::correct()\n"
 245       "interpolationMethod " << interpolationMethod << " not known\n"
 246       "interpolationMethod must be patchToPatch or ggi"
 247            << exit(FatalError);
 248       }
 249
 250     // now put global back into local
 251     const label slavePatchStart
 252       = mesh.boundaryMesh()[slavePatchIndex].start();
 253
 254     forAll(masterDUInterpToSlave, i)
 255       {
 256     masterDUInterpToSlave[i] =
 257       globalMasterDUInterpToSlave
 258       [
 259        mesh.faceZones()[slaveFaceZoneID()].whichFace(slavePatchStart + i)
 260        ];
 261       }
 262
 263     // Now masterDUInterpToSlave should be masterDU interpolated to the slave
 264
 265     // calculate slave shear traction increments
 266     const scalarField magSlavePressure = mag(slavePressure);
 267     label numSlipFaces = 0;
 268     label numStickFaces = 0;
 269     scalarField& stickSlip = stickSlipFaces();
 270     scalar frictionPenaltyFac = frictionPenaltyFactor();
 271     scalar maxMagSlip = 0.0; // of stick faces
 272     //scalar slipTol = 5e-8;
 273
 274     // I am playing with things at the moment
 275     FatalError << "iterativePenaltyFriction can not be used at the moment"
 276            << " as I am changing a few things around: philipc"
 277            << exit(FatalError);
 278
 279     forAll(magSlavePressure, faceI)
 280       {
 281     // there can only be a frictional tangential force when there is
 282     // a positive pressure
 283     if (magSlavePressure[faceI] > 1e5) //SMALL)
 284       {
 285         // compute slip
 286         //- we need the difference of DU between the master and slave
 287         vector slip =
 288           slaveDU[faceI] - masterDUInterpToSlave[faceI];
 289         //- the shear traction direction is got by removing the normal
 290         // component of the DU
 291         //- (I - sqr(n)) removes the normal
 292         //- sqr(n) would remove the shear
 293         slip = (I - sqr(slaveFaceNormals[faceI])) & slip;
 294
 295         // we define stick faces as those with a slip less than slipTol
 296         // if (mag(slip) < slipTol)
 297         //   slip = vector::zero;
 298         // slip -= slipTol*slip/mag(slip);
 299
 300         // traction to cause slipping
 301         //scalar slipTrac =
 302         //   frictionLawPtr_->slipTraction(magSlavePressure[faceI]);
 303         vector& slaveTrac = slaveTraction_[faceI];
 304
 305         // if mag(slaveTrac) is greater than slipTrac and
 306         // it acts in the opposite direction to slip
 307         // then the face is a slip face
 308         //scalar slipFunc = mag(slaveTrac) - 0.99*slipTrac;
 309
 310         // if (slipFunc > 0) // slip face
 311         //   {
 312         //      // slave traction acts in opposite direction to slip
 313         //      // slowly bring traction to the slip traction
 314         //      //slaveTrac +=
 315         // frictionPenaltyScale_*( (-slipTrac * (slip/mag(slip))) - slaveTrac );
 316         //      slaveTrac = -slipTrac * slip / mag(slip);
 317         //      numSlipFaces++;
 318         //      stickSlip[faceI] = 1;
 319         //   }
 320         // else // stick face
 321           {
 322         // we iteratively increase the traction until
 323         // the slip is zero (ie less than the slipTol)
 324         // when slip is zero, the slave trac remains unchanged
 325         // if (mag(slip) < slipTol)
 326         //   {
 327         //     slip = vector::zero;
 328         //   }
 329         // else
 330         //   {
 331         //     slip -= slipTol*slip/mag(slip);
 332         //   }
 333         // if (mag(slip) > slipTol)
 334         //   {
 335         //     slaveTrac -= frictionPenaltyFac * slip;
 336         //   }
 337
 338         // pure penalty
 339         slaveTrac = -frictionPenaltyFac * slip;
 340         numStickFaces++;
 341         //stickSlip[faceI] = 2;
 342         stickSlip[faceI] = slip.x();
 343         maxMagSlip = max(maxMagSlip, mag(slip));
 344           }
 345
 346
 347         //maxMagSlip = max(maxMagSlip, max(mag(Slip)));
 348
 349         //- reduce the traction
 350         //scalar slipTrac = frictionCoeff_*magSlavePressure[faceI];
 351         // scalar slipTrac =
 352           // frictionLawPtr_->slipTraction(magSlavePressure[faceI]);
 353         // vector& slaveTrac = slaveTraction_[faceI];
 354         // slaveTrac -=
 355         //   frictionPenaltyFac * slip;
 356
 357         //- limit traction
 358         // if (mag(slaveTrac) > slipTrac)
 359         // if (mag(slaveTrac) > 0.99*slipTrac)
 360         //   {
 361         //      //vector dir = slaveTrac / mag(slaveTrac);
 362         //      vector dir = -slip/mag(slip);
 363         //      // slaveTrac = slipTrac * dir;
 364         //      slaveTrac = slipTrac * dir;
 365         //      numSlipFaces++;
 366         //      stickSlip[faceI] = 1;
 367         //   }
 368         // else if (mag(slip) < SMALL)
 369         //   {
 370         //      numStickFaces++;
 371         //      stickSlip[faceI] = 2;
 372         //   }
 373         // else if ((slaveTrac & slip) > 0.0) // acting in same direction
 374         //   {
 375         //      slaveTrac = vector::zero;
 376         //      numSlipFaces++;
 377         //      stickSlip[faceI] = 1;
 378         //   }
 379         // else
 380         //   {
 381         //      numStickFaces++;
 382         //      stickSlip[faceI] = 2;
 383         //   }
 384       }
 385     // no friction if pressure is negative or zero
 386     else
 387       {
 388         slaveTraction_[faceI] = vector::zero;
 389         //slaveTraction_[faceI] -= frictionPenaltyScale_*slaveTraction_[faceI];
 390         stickSlip[faceI] = 0;
 391       }
 392       }
 393
 394     if (oscillationCorr_)
 395       {
 396     correctOscillations(slaveTraction_, slaveFaceZonePatch);
 397       }
 398
 399     scalar maxMagSlaveTraction = gMax(mag(slaveTraction_));
 400     reduce(numSlipFaces, sumOp<int>());
 401     reduce(numStickFaces, sumOp<int>());
 402
 403     // master writes to contact info file
 404     if (Pstream::master() && (contactIterNum_ %  infoFreq_ == 0))
 405       {
 406           OFstream& contactFile = *contactFilePtr_;
 407           int width = 20;
 408           contactFile << mesh.time().value();
 409           contactFile.width(width);
 410           contactFile << contactIterNum_;
 411           contactFile.width(width);
 412           contactFile << frictionPenaltyScale_;
 413           contactFile.width(width);
 414           contactFile << numSlipFaces;
 415           contactFile.width(width);
 416           contactFile << numStickFaces;
 417           contactFile.width(width);
 418           contactFile << maxMagSlaveTraction;
 419           contactFile.width(width);
 420           contactFile << maxMagSlip;
 421           contactFile << endl;
 422       }
 423   }
 424
 425
 426   void Foam::iterativePenaltyFriction::correctOscillations
 427   (
 428    vectorField& slaveTraction,
 429    const PrimitivePatch<face, List, pointField>& slaveFaceZonePatch
 430    )
 431   {
 432     // oscillations sometimes appear in contact shear traction
 433     // so we will try to limit them here
 434     // we will weight the current face slaveTraction with the average of the
 435     // neighbours using the weight oscillationCorrectionFactor_
 436     // we will just change the magnitude of slaveTraction and leave the
 437       // direction
 438
 439     //Pout << "Applying contact shear oscillation correction..." << flush;
 440
 441     const fvMesh& mesh = mesh_;
 442     const label slavePatchIndex = slavePatchID();
 443     const labelListList& faceFaces = slaveFaceZonePatch.faceFaces();
 444
 445     // create global slaveTraction
 446     vectorField globalSlaveTraction(slaveFaceZonePatch.size(), vector::zero);
 447     const label slavePatchStart
 448       = mesh.boundaryMesh()[slavePatchIndex].start();
 449     forAll(slaveTraction, i)
 450     {
 451         globalSlaveTraction
 452             [
 453                 mesh.faceZones()[slaveFaceZoneID()
 454                     ].whichFace(slavePatchStart + i)] =
 455             slaveTraction[i];
 456       }
 457     // sum because each face is only on one proc
 458     reduce(globalSlaveTraction, sumOp<vectorField>());
 459
 460     //scalarField globalSlaveMagTraction = mag(globalSlaveTraction);
 461
 462     // smooth mag of slaveTraction with face face disps
 463     forAll(faceFaces, facei)
 464       {
 465     //label numFaceFaces = faceFaces[facei].size();
 466
 467     // only smooth faces in contact and don't smooth
 468     // end/corner faces
 469     // if (globalSlaveMagTraction[facei] > SMALL)
 470     if (mag(globalSlaveTraction[facei]) > SMALL)
 471       //&&
 472       //numFaceFaces > 1) // don't smooth end/corner faces
 473       {
 474         //scalar avTracMag = 0.0;
 475         vector avTrac = vector::zero;
 476         int numNei = 0;
 477         forAll(faceFaces[facei], ffi)
 478           {
 479         label faceFace = faceFaces[facei][ffi];
 480
 481         //avTracMag += globalSlaveMagTraction[faceFace];
 482         avTrac += globalSlaveTraction[faceFace];
 483         numNei++;
 484           }
 485
 486         // avTracMag /= numNei;
 487         avTrac /= numNei;
 488
 489         // if (numFaceFaces == 1)
 490         //   {
 491         //      // for corner/end faces, decrease the weight of the neighbours
 492         //      avTracMag += globalSlaveTraction[facei];
 493         //      avTracMag /= 2;
 494         //   }
 495
 496         // weighted-average with face-faces
 497         //vector dir = globalSlaveTraction[facei];
 498         //dir /= mag(dir);
 499         globalSlaveTraction[facei] =
 500           oscillationCorrFac_*globalSlaveTraction[facei]
 501             + (1.0-oscillationCorrFac_)*avTrac;
 502       }
 503       }
 504
 505     // convert global back to local
 506     forAll(slaveTraction, facei)
 507       {
 508     slaveTraction[facei] =
 509       globalSlaveTraction
 510       [
 511        mesh.faceZones()[slaveFaceZoneID()].whichFace(slavePatchStart + facei)
 512        ];
 513       }
 514
 515     //Pout << "\tdone" << endl;
 516   }
 517
 518
 519   void Foam::iterativePenaltyFriction::calcFrictionPenaltyFactor()
 520   {
 521     // set penalty factor using a similar method to the normal
 522     // contact where we approx penaltyFactor from mechanical properties
 523     // this can then be scaled using the penaltyScale
 524     const label masterPatchIndex =  masterPatchID();
 525     const label slavePatchIndex =  slavePatchID();
 526     const constitutiveModel& rheology =
 527       mesh_.objectRegistry::lookupObject<constitutiveModel>
 528         ("rheologyProperties");
 529     scalarField masterMu =
 530       rheology.mu()().boundaryField()[masterPatchIndex];
 531     scalarField slaveMu =
 532       rheology.mu()().boundaryField()[slavePatchIndex];
 533
 534     // avarage contact patch shear modulus
 535     scalar shearModulus = 0.5*(gAverage(masterMu)+gAverage(slaveMu));
 536
 537     // average contact patch face area
 538     scalar masterMagSf =
 539         gAverage(mesh_.magSf().boundaryField()[masterPatchIndex]);
 540     scalar slaveMagSf =
 541         gAverage(mesh_.magSf().boundaryField()[slavePatchIndex]);
 542     scalar faceArea = 0.5*(masterMagSf+slaveMagSf);
 543
 544     // average contact patch cell volume
 545     scalarField masterV(mesh_.boundary()[masterPatchIndex].size(), 0.0);
 546     scalarField slaveV(mesh_.boundary()[slavePatchIndex].size(), 0.0);
 547     const volScalarField::DimensionedInternalField & V = mesh_.V();
 548     {
 549       const unallocLabelList& faceCells =
 550           mesh_.boundary()[masterPatchIndex].faceCells();
 551       forAll(mesh_.boundary()[masterPatchIndex], facei)
 552     {
 553       masterV[facei] = V[faceCells[facei]];
 554     }
 555     }
 556     {
 557       const unallocLabelList& faceCells =
 558           mesh_.boundary()[slavePatchIndex].faceCells();
 559       forAll(mesh_.boundary()[slavePatchIndex], facei)
 560     {
 561       slaveV[facei] = V[faceCells[facei]];
 562     }
 563     }
 564     scalar cellVolume = 0.5*(gAverage(masterV)+gAverage(slaveV));
 565
 566     // approximate penalty factor based on Hallquist et al.
 567     // we approximate penalty factor for traction instead of force
 568     frictionPenaltyFactorPtr_ =
 569         new scalar(frictionPenaltyScale_*shearModulus*faceArea/cellVolume);
 570   }
 571
 572   void Foam::iterativePenaltyFriction::writeDict(Ostream& os) const
 573   {
 574     word keyword(name()+"FrictionModelDict");
 575     os.writeKeyword(keyword)
 576       << frictionContactModelDict_;
 577   }
 578
 579
 580 // ************************************************************************* //