applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "contactPatchPair.H"
  29 #include "contactProblem.H"
  30 #include "surfaceFields.H"
  31
  32 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  33
  34 // Construct from components
  35 Foam::contactPatchPair::contactPatchPair
  36 (
  37     const word& name,
  38     const contactProblem& cp,
  39     const word& masterPatchName,
  40     const word& slavePatchName,
  41     const dimensionedScalar& frictionCoeff,
  42     const scalar contactTol,
  43     const intersection::algorithm alg,
  44     const intersection::direction dir
  45 )
  46 :
  47     name_(name),
  48     cp_(cp),
  49     masterPatch_(masterPatchName, cp.mesh().boundaryMesh()),
  50     slavePatch_(slavePatchName, cp.mesh().boundaryMesh()),
  51     frictionCoeff_(frictionCoeff),
  52     contactTol_(contactTol),
  53     masterInterpolate_
  54     (
  55         cp.mesh().boundaryMesh()[masterPatch_.index()]
  56     ),
  57     slaveInterpolate_
  58     (
  59         cp.mesh().boundaryMesh()[slavePatch_.index()]
  60     ),
  61     masterToSlaveInterpolate_
  62     (
  63         cp.mesh().boundaryMesh()[masterPatch_.index()],   // from patch
  64         cp.mesh().boundaryMesh()[slavePatch_.index()],    // to patch
  65         alg,
  66         dir
  67     ),
  68     slaveToMasterInterpolate_
  69     (
  70         cp.mesh().boundaryMesh()[slavePatch_.index()],    // from patch
  71         cp.mesh().boundaryMesh()[masterPatch_.index()],   // to patch
  72         alg,
  73         dir
  74     )
  75 {}
  76
  77
  78 // Construct from dictionary
  79 Foam::contactPatchPair::contactPatchPair
  80 (
  81     const word& name,
  82     const contactProblem& cp,
  83     const dictionary& dict
  84 )
  85 :
  86     name_(name),
  87     cp_(cp),
  88     masterPatch_(dict.lookup("masterPatch"), cp.mesh().boundaryMesh()),
  89     slavePatch_(dict.lookup("slavePatch"), cp.mesh().boundaryMesh()),
  90     frictionCoeff_(dict.lookup("frictionCoeff")),
  91     contactTol_(readScalar(dict.lookup("contactTol"))),
  92     masterInterpolate_
  93     (
  94         cp.mesh().boundaryMesh()[masterPatch_.index()]
  95     ),
  96     slaveInterpolate_
  97     (
  98         cp.mesh().boundaryMesh()[slavePatch_.index()]
  99     ),
 100     masterToSlaveInterpolate_
 101     (
 102         cp.mesh().boundaryMesh()[masterPatch_.index()],    // from patch
 103         cp.mesh().boundaryMesh()[slavePatch_.index()],     // to patch
 104         intersection::algorithmNames_.read(dict.lookup("projectionAlgo")),
 105         intersection::directionNames_.read(dict.lookup("projectionDir"))
 106
 107     ),
 108     slaveToMasterInterpolate_
 109     (
 110         cp.mesh().boundaryMesh()[slavePatch_.index()],     // from patch
 111         cp.mesh().boundaryMesh()[masterPatch_.index()],    // to patch
 112         intersection::algorithmNames_.read(dict.lookup("projectionAlgo")),
 113         intersection::directionNames_.read(dict.lookup("projectionDir"))
 114
 115     )
 116 {}
 117
 118
 119 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 120
 121
 122 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 123
 124 Foam::tmp<Foam::scalarField>
 125 Foam::contactPatchPair::masterTouchFraction() const
 126 {
 127     // Get reference to displacement field and mesh
 128     const volVectorField& U = cp_.U();
 129     const fvMesh& mesh = cp_.mesh();
 130
 131     // Interpolate slave displacement into master vertices
 132     vectorField masterVertexU =
 133         slaveToMasterInterpolate_.pointInterpolate<vector>
 134         (
 135             slaveInterpolate_.faceToPointInterpolate
 136             (
 137                 U.boundaryField()[slavePatch_.index()]
 138             )
 139         );
 140
 141     const vectorField& projectionDir =
 142         mesh.boundaryMesh()[masterPatch_.index()].pointNormals();
 143
 144
 145     // Calculate master gap function
 146     scalarField vertexMasterGap =
 147     (
 148         (
 149             masterVertexU
 150           - masterInterpolate_.faceToPointInterpolate
 151             (
 152                 U.boundaryField()[masterPatch_.index()]
 153             )
 154         )
 155         & projectionDir
 156     ) + slaveToMasterInterpolate_.pointDistanceToIntersection() - contactTol_;
 157
 158     // Calculate area in contact
 159
 160     const faceList& masterPatchLocalFaces =
 161         mesh.boundaryMesh()[masterPatch_.index()].localFaces();
 162
 163     const pointField& masterPatchLocalPoints =
 164         mesh.boundaryMesh()[masterPatch_.index()].localPoints();
 165
 166     tmp<scalarField> ttouchFrac
 167     (
 168         new scalarField(masterPatchLocalFaces.size(), 0)
 169     );
 170     scalarField& touchFrac = ttouchFrac();
 171
 172     forAll (masterPatchLocalFaces, faceI)
 173     {
 174         touchFrac[faceI] =
 175             masterPatchLocalFaces[faceI].areaInContact
 176             (
 177                 masterPatchLocalPoints,
 178                 vertexMasterGap
 179             );
 180     }
 181
 182     return ttouchFrac;
 183 }
 184
 185
 186 Foam::tmp<Foam::scalarField>
 187 Foam::contactPatchPair::slaveTouchFraction() const
 188 {
 189     // Get reference to displacement field and mesh
 190     const volVectorField& U = cp_.U();
 191     const fvMesh& mesh = cp_.mesh();
 192
 193     // Interpolate master displacement into slave vertices
 194     vectorField slaveVertexU =
 195         masterToSlaveInterpolate_.pointInterpolate<vector>
 196         (
 197             masterInterpolate_.faceToPointInterpolate
 198             (
 199                 U.boundaryField()[masterPatch_.index()]
 200             )
 201         );
 202
 203     const vectorField& projectionDir =
 204         mesh.boundaryMesh()[slavePatch_.index()].pointNormals();
 205
 206
 207     // Calculate slave gap function
 208     scalarField vertexSlaveGap =
 209     (
 210         (
 211             slaveVertexU
 212           - slaveInterpolate_.faceToPointInterpolate
 213             (
 214                 U.boundaryField()[slavePatch_.index()]
 215             )
 216         )
 217         & projectionDir
 218     ) + masterToSlaveInterpolate_.pointDistanceToIntersection() - contactTol_;
 219
 220     // Calculate area in contact
 221
 222     const faceList& slavePatchLocalFaces =
 223         mesh.boundaryMesh()[slavePatch_.index()].localFaces();
 224
 225     const pointField& slavePatchLocalPoints =
 226         mesh.boundaryMesh()[slavePatch_.index()].localPoints();
 227
 228     tmp<scalarField> ttouchFrac
 229     (
 230         new scalarField(slavePatchLocalFaces.size(), 0)
 231     );
 232     scalarField& touchFrac = ttouchFrac();
 233
 234     forAll (slavePatchLocalFaces, faceI)
 235     {
 236         touchFrac[faceI] =
 237             slavePatchLocalFaces[faceI].areaInContact
 238             (
 239                 slavePatchLocalPoints,
 240                 vertexSlaveGap
 241             );
 242     }
 243
 244     return ttouchFrac;
 245 }
 246
 247
 248 void Foam::contactPatchPair::correct
 249 (
 250     const FieldField<Field, vector>& curTraction,
 251     FieldField<Field, vector>& newTraction,
 252     FieldField<Field, vector>& refValue,
 253     FieldField<Field, scalar>& valueFraction
 254 )
 255 {
 256     // Get reference to displacement field and mesh
 257     const volVectorField::GeometricBoundaryField& Upatches =
 258         cp_.U().boundaryField();
 259
 260     const fvMesh& mesh = cp_.mesh();
 261     const surfaceVectorField::GeometricBoundaryField& Apatches =
 262         mesh.Sf().boundaryField();
 263     const surfaceScalarField::GeometricBoundaryField& magApatches =
 264         mesh.magSf().boundaryField();
 265
 266     // Get patch indices
 267     const label masterIndex = masterPatch_.index();
 268     const label slaveIndex = slavePatch_.index();
 269
 270     // Calculate patch normals
 271     vectorField nMasterPatch = Apatches[masterIndex]/magApatches[masterIndex];
 272
 273     vectorField nSlavePatch = Apatches[slaveIndex]/magApatches[slaveIndex];
 274
 275
 276     // Calculate slave pressure and tangential force
 277
 278     scalarField slavePressure = -( nSlavePatch & curTraction[slaveIndex]);
 279
 280     // Enforce gradient condition on the master patch
 281
 282     // Calculate relative tangential velocity for master patch
 283     vectorField relUmaster =
 284         slaveToMasterInterpolate_.faceInterpolate<vector>
 285         (
 286             Upatches[slaveIndex]
 287         )
 288       - Upatches[masterIndex];
 289
 290     relUmaster -= nMasterPatch*(nMasterPatch & relUmaster);
 291     relUmaster /= mag(relUmaster) + VSMALL;
 292
 293     // Calculate tangential master traction
 294     scalarField magMasterTangential =
 295         Foam::mag((I - nMasterPatch*nMasterPatch) & curTraction[masterIndex]);
 296
 297     // Calculate master pressure
 298     scalarField masterPressure =
 299         max
 300         (
 301             slaveToMasterInterpolate_.faceInterpolate<scalar>
 302             (
 303                 slavePressure
 304             ),
 305             scalar(0)
 306         );
 307
 308     // Calculate master traction, using the positive part of
 309     // slave pressure and tangential fricton
 310     // Mind the signs: pressure = negative gradient (minus master normal)
 311     //                 friction = positive pressure
 312     newTraction[masterIndex] +=
 313         masterTouchFraction()*
 314         (
 315            -nMasterPatch*masterPressure
 316           + relUmaster*
 317             min
 318             (
 319                 frictionCoeff_.value()*masterPressure,
 320                 magMasterTangential
 321             )
 322         );
 323
 324     // Enforce direction mixed condition on the slave patch
 325
 326     // Calculate slave fraction.  Correct for negative pressure
 327     // (if the pressure is negative, the contact is released)
 328     //HJ, fiddle pos pressure!!!
 329     scalarField slaveFrac = slaveTouchFraction();
 330
 331     // Calculate slave displacement
 332     vectorField slaveVertexU =
 333         masterToSlaveInterpolate_.pointInterpolate<vector>
 334         (
 335             masterInterpolate_.faceToPointInterpolate
 336             (
 337                 Upatches[masterIndex]
 338             )
 339         );
 340
 341     const vectorField& projectionDir =
 342         mesh.boundaryMesh()[slaveIndex].pointNormals();
 343
 344     // Calculate slave displacement
 345     vectorField slaveDisp =
 346         slaveInterpolate_.pointToFaceInterpolate
 347         (
 348             slaveVertexU
 349           + masterToSlaveInterpolate_.pointDistanceToIntersection()
 350             *projectionDir
 351         );
 352
 353     // Accumulate normal of slave displacement
 354     refValue[slaveIndex] +=
 355         nSlavePatch*
 356         min
 357         (
 358             pos(slaveFrac)*
 359             (
 360                 (nSlavePatch & Upatches[slaveIndex])
 361               + slaveFrac*contactTol_
 362             ),
 363             (nSlavePatch & slaveDisp)
 364         );
 365
 366
 367     // Accumulate slave friction
 368
 369     // Calculate relative tangential velocity for slave patch
 370     vectorField relUslave =
 371         masterToSlaveInterpolate_.faceInterpolate<vector>
 372         (
 373             Upatches[masterIndex]
 374         )
 375       - Upatches[slaveIndex];
 376
 377     relUslave -= nSlavePatch*(nSlavePatch & relUslave);
 378     relUslave /= mag(relUslave) + VSMALL;
 379
 380     // Take out normal component out of slave traction and find the
 381     // magnitude of the tangential traction.
 382     scalarField magSlaveTangential =
 383         Foam::mag((I - nSlavePatch*nSlavePatch) & curTraction[slaveIndex]);
 384
 385     // Calculate slave traction
 386     newTraction[slaveIndex] +=
 387         slaveFrac*relUslave*
 388         min
 389         (
 390             frictionCoeff_.value()*max(slavePressure, scalar(0)),
 391             magSlaveTangential
 392         );
 393
 394     // Accumulate slave touch fraction
 395     valueFraction[slaveIndex] += slaveFrac;
 396
 397 /*
 398     Info << "slavePressure: " << slavePressure << nl
 399          << "slaveTouchFrac: " << slaveTouchFraction() << nl
 400 //          << "slaveFrac: " << slaveFrac << nl
 401          << "refValueSlave: " << refValue[slaveIndex].component(vector::Y) << nl
 402 //          << "slaveTraction: " << newTraction[slaveIndex] << nl
 403          << "masterTouchFrac: " << masterTouchFraction() << nl
 404 //          << "interpolated slave pressure: "
 405 //          << slaveToMasterInterpolate_.faceInterpolate<scalar>
 406 //             (
 407 //                 slavePressure
 408 //             )
 409 //          << nl
 410 //          << "masterTraction: "
 411 //          << newTraction[masterIndex].component(vector::Y)
 412          << endl;
 413 */
 414 }
 415
 416
 417 void Foam::contactPatchPair::writeDict(Ostream& os) const
 418 {
 419     os  << nl << name() << nl << token::BEGIN_BLOCK;
 420
 421     os  << "masterPatch " << masterPatch_.name() << token::END_STATEMENT << nl
 422         << "slavePatch " << slavePatch_.name() << token::END_STATEMENT << nl
 423         << "frictionCoeff " << frictionCoeff_ << token::END_STATEMENT << nl
 424         << "contactTol " << contactTol_ << token::END_STATEMENT << nl
 425         << "projectionAlgo "
 426         << intersection::algorithmNames_
 427                [masterToSlaveInterpolate_.projectionAlgo()]
 428         << token::END_STATEMENT << nl
 429         << "projectionDir "
 430         << intersection::directionNames_
 431                [masterToSlaveInterpolate_.projectionDir()]
 432         << token::END_STATEMENT << nl
 433         << token::END_BLOCK << endl;
 434 }
 435
 436
 437 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 438
 439
 440 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 441
 442
 443 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 444
 445
 446 // ************************************************************************* //