applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactProblem.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     Class describes a multiple body contact problem.  Each individual contact
  26     is described by a contactPatchPair.  contactProblem handles
  27     multiple contact updates and sets the boundary conditions on the
  28     displacement field.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "contactProblem.H"
  33 #include "fvMesh.H"
  34 #include "FieldFields.H"
  35 #include "directionMixedFvPatchFields.H"
  36 #include "surfaceFields.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  44
  45 defineTypeNameAndDebug(contactProblem, 0);
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // Read constructor given IOobject
  50 contactProblem::contactProblem
  51 (
  52     volVectorField& U,
  53     const volTensorField& gradU
  54 )
  55 :
  56     IOdictionary
  57     (
  58         IOobject
  59         (
  60             "contactProperties",
  61             U.time().constant(),
  62             U.db(),
  63             IOobject::MUST_READ,
  64             IOobject::NO_WRITE
  65         )
  66     ),
  67     contactPatchPairList(),
  68     U_(U),
  69     gradU_(gradU),
  70     urfValue_(readScalar(lookup("urfValue"))),
  71     urfTraction_(readScalar(lookup("urfTraction"))),
  72     urfFraction_(readScalar(lookup("urfFraction")))
  73 {
  74     // Read contactPatchPairList
  75     Istream& is = lookup("contacts");
  76
  77     PtrList<entry> contactEntries(is);
  78
  79     contactPatchPairList& contacts = *this;
  80
  81     contacts.setSize(contactEntries.size());
  82
  83     forAll(contacts, contactI)
  84     {
  85         contacts.set
  86         (
  87             contactI,
  88             new contactPatchPair
  89             (
  90                 contactEntries[contactI].keyword(),
  91                 *this,
  92                 contactEntries[contactI].dict()
  93             )
  94         );
  95     }
  96 }
  97
  98
  99 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 100
 101
 102 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 103
 104 void contactProblem::correct()
 105 {
 106     contactPatchPairList& contacts = *this;
 107
 108     // Create fields for accumulation
 109     volVectorField::GeometricBoundaryField& Upatches = U().boundaryField();
 110
 111     FieldField<Field, vector> curTraction(Upatches.size());
 112     FieldField<Field, vector> newTraction(Upatches.size());
 113     FieldField<Field, vector> refValue(Upatches.size());
 114     FieldField<Field, scalar> valueFraction(Upatches.size());
 115
 116     forAll (Upatches, patchI)
 117     {
 118         curTraction.set
 119         (
 120             patchI,
 121             new vectorField(Upatches[patchI].size(), vector::zero)
 122         );
 123
 124         newTraction.set
 125         (
 126             patchI,
 127             new vectorField(Upatches[patchI].size(), vector::zero)
 128         );
 129
 130         refValue.set
 131         (
 132             patchI,
 133             new vectorField(Upatches[patchI].size(), vector::zero)
 134         );
 135
 136         valueFraction.set
 137         (
 138             patchI,
 139             new scalarField(Upatches[patchI].size(), 0)
 140         );
 141     }
 142
 143     // Collect patches involved in contact
 144     boolList contactPatches(Upatches.size(), false);
 145
 146     forAll (contacts, contactI)
 147     {
 148         contactPatches[contacts[contactI].masterPatch().index()] = true;
 149         contactPatches[contacts[contactI].slavePatch().index()] = true;
 150     }
 151
 152     // Calculate the traction for all involved patches
 153
 154     // Collect fields
 155     const volTensorField::GeometricBoundaryField& gradUpatches =
 156         gradU().boundaryField();
 157
 158     const surfaceVectorField::GeometricBoundaryField& Apatches =
 159         mesh().Sf().boundaryField();
 160     const surfaceScalarField::GeometricBoundaryField& magApatches =
 161         mesh().magSf().boundaryField();
 162
 163     // Lookup mu and lambda form object registry
 164     const volScalarField& mu =
 165         mesh().objectRegistry::lookupObject<volScalarField>("mu");
 166
 167     const volScalarField::GeometricBoundaryField& muPatches =
 168         mu.boundaryField();
 169
 170     const volScalarField& lambda =
 171         mesh().objectRegistry::lookupObject<volScalarField>("lambda");
 172
 173     const volScalarField::GeometricBoundaryField& lambdaPatches =
 174         lambda.boundaryField();
 175
 176     forAll (Upatches, patchI)
 177     {
 178         if (contactPatches[patchI])
 179         {
 180             vectorField nPatch = Apatches[patchI]/magApatches[patchI];
 181
 182             curTraction[patchI] =
 183                 nPatch &
 184                 (
 185                     muPatches[patchI]*
 186                     (
 187                         gradUpatches[patchI]
 188                       + gradUpatches[patchI].T()
 189                     )
 190                   + I*(lambdaPatches[patchI]*tr(gradUpatches[patchI]))
 191                 );
 192         }
 193     }
 194
 195     // Accumulate contact data and active patches
 196     forAll (contacts, contactI)
 197     {
 198         contacts[contactI].correct
 199         (
 200             curTraction,
 201             newTraction,
 202             refValue,
 203             valueFraction
 204         );
 205     }
 206
 207     // Enforce accumulated contact onto the patches
 208     forAll (Upatches, patchI)
 209     {
 210         if (contactPatches[patchI])
 211         {
 212             // Cast the patch into direction mixed type
 213             directionMixedFvPatchVectorField& curUPatch =
 214                 refCast<directionMixedFvPatchVectorField>(Upatches[patchI]);
 215
 216             // Set the values using under-relaxation
 217             curUPatch.refValue() =
 218                 (1.0 - urfValue_)*curUPatch.refValue()
 219               + urfValue_*refValue[patchI];
 220
 221             // Calculate the gradient from under-relaxad accumulated traction
 222             vectorField nPatch = Apatches[patchI]/magApatches[patchI];
 223
 224             curUPatch.refGrad() =
 225             (
 226                 (1.0 - urfTraction_)*curTraction[patchI]
 227               + urfTraction_*newTraction[patchI]
 228
 229               - (nPatch &
 230                 (
 231                     muPatches[patchI]*gradUpatches[patchI].T()
 232                   - (
 233                         muPatches[patchI]
 234                       + lambdaPatches[patchI]
 235                     )*gradUpatches[patchI]
 236                 )
 237                 )
 238
 239               - nPatch*
 240                 (
 241                     lambdaPatches[patchI]*tr(gradUpatches[patchI])
 242                 )
 243
 244             )/(2.0*muPatches[patchI] + lambdaPatches[patchI]);
 245
 246             // Set the value fractions
 247             curUPatch.valueFraction() =
 248                 (1.0 - urfFraction_)*curUPatch.valueFraction()
 249               + I*urfFraction_*valueFraction[patchI];
 250         }
 251     }
 252 }
 253
 254
 255 tmp<volScalarField> contactProblem::contactArea() const
 256 {
 257     tmp<volScalarField> tca
 258     (
 259         new volScalarField
 260         (
 261             IOobject
 262             (
 263                 "contactArea",
 264                 U().time().timeName(),
 265                 U().db(),
 266                 IOobject::NO_READ,
 267                 IOobject::AUTO_WRITE
 268             ),
 269             mesh(),
 270             dimensionedScalar(0)
 271         )
 272     );
 273
 274     volScalarField& ca = tca();
 275
 276     // Set contact area boundary
 277     const contactPatchPairList& contacts = *this;
 278
 279     forAll (contacts, contactI)
 280     {
 281         // Get master contact
 282         ca.boundaryField()[contacts[contactI].masterPatch().index()] +=
 283             contacts[contactI].masterTouchFraction();
 284
 285         // Get slave contact
 286         ca.boundaryField()[contacts[contactI].slavePatch().index()] +=
 287             contacts[contactI].slaveTouchFraction();
 288     }
 289
 290     return tca;
 291 }
 292
 293
 294 // Return a list of contactPatchPair names
 295 wordList contactProblem::names() const
 296 {
 297     const contactPatchPairList& contacts = *this;
 298
 299     wordList t(contacts.size());
 300
 301     forAll (contacts, contactI)
 302     {
 303         t[contactI] = contacts[contactI].name();
 304     }
 305
 306     return t;
 307 }
 308
 309
 310 bool contactProblem::read()
 311 {
 312     if (regIOobject::read())
 313     {
 314         urfValue_ = readScalar(lookup("urfValue"));
 315         urfTraction_ = readScalar(lookup("urfTraction"));
 316         urfFraction_ = readScalar(lookup("urfFraction"));
 317
 318         // Decided not to re-read contactPatchPairList.  HJ, 10/Jul/2004
 319         return true;
 320     }
 321     else
 322     {
 323         return false;
 324     }
 325 }
 326
 327
 328 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 329
 330
 331 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 332
 333 } // End namespace Foam
 334
 335 // ************************************************************************* //