src/solidModels/constitutiveModel/cohesiveLaws/linearCohesiveLaw/linearCohesiveLaw.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 "linearCohesiveLaw.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "zeroGradientFvPatchFields.H"
  29 #include "fvc.H"
  30 #include "cohesiveFvPatch.H"
  31 #include "cohesiveLaw.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     defineTypeNameAndDebug(linearCohesiveLaw, 0);
  38     addToRunTimeSelectionTable(cohesiveLaw, linearCohesiveLaw, dictionary);
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 // Construct from dictionary
  47 Foam::linearCohesiveLaw::linearCohesiveLaw
  48 (
  49     const word& name,
  50     const volSymmTensorField& sigma,
  51     const dictionary& dict
  52 )
  53 :
  54   cohesiveLaw(name, sigma, dict),
  55   GIc_(dict.lookup("GIc")),
  56   GIIc_(dict.lookup("GIIc")),
  57   sigmaMax_(dict.lookup("sigmaMax")),
  58   tauMax_(dict.lookup("tauMax"))
  59 {}
  60
  61
  62 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  63
  64 Foam::linearCohesiveLaw::~linearCohesiveLaw()
  65 {}
  66
  67
  68 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  69
  70 Foam::tmp<Foam::volScalarField> Foam::linearCohesiveLaw::materials() const
  71 {
  72   notImplemented(type() + "::materials()");
  73
  74     return tmp<volScalarField>
  75     (
  76         new volScalarField
  77         (
  78             IOobject
  79             (
  80                 "materials",
  81                 mesh().time().timeName(),
  82                 mesh(),
  83                 IOobject::NO_READ,
  84                 IOobject::NO_WRITE
  85             ),
  86             mesh(),
  87             dimensionedScalar("zero", dimless, 0.0)
  88         )
  89     );
  90 }
  91
  92 Foam::tmp<Foam::surfaceScalarField> Foam::linearCohesiveLaw::sigmaMax() const
  93 {
  94     return tmp<surfaceScalarField>
  95     (
  96         new surfaceScalarField
  97         (
  98             IOobject
  99             (
 100                 "sigmaMax",
 101                 mesh().time().timeName(),
 102                 mesh(),
 103                 IOobject::NO_READ,
 104                 IOobject::NO_WRITE
 105             ),
 106             mesh(),
 107             sigmaMax_
 108         )
 109     );
 110 }
 111
 112 Foam::tmp<Foam::surfaceScalarField> Foam::linearCohesiveLaw::tauMax() const
 113 {
 114     return tmp<surfaceScalarField>
 115     (
 116         new surfaceScalarField
 117         (
 118             IOobject
 119             (
 120                 "tauMax",
 121                 mesh().time().timeName(),
 122                 mesh(),
 123                 IOobject::NO_READ,
 124                 IOobject::NO_WRITE
 125             ),
 126             mesh(),
 127             tauMax_
 128         )
 129     );
 130 }
 131
 132 Foam::tmp<Foam::surfaceScalarField> Foam::linearCohesiveLaw::GIc() const
 133 {
 134     return tmp<surfaceScalarField>
 135     (
 136         new surfaceScalarField
 137         (
 138             IOobject
 139             (
 140                 "GIc",
 141                 mesh().time().timeName(),
 142                 mesh(),
 143                 IOobject::NO_READ,
 144                 IOobject::NO_WRITE
 145             ),
 146             mesh(),
 147             GIc_
 148         )
 149     );
 150 }
 151
 152 Foam::tmp<Foam::surfaceScalarField> Foam::linearCohesiveLaw::GIIc() const
 153 {
 154     return tmp<surfaceScalarField>
 155     (
 156         new surfaceScalarField
 157         (
 158             IOobject
 159             (
 160                 "GIIc",
 161                 mesh().time().timeName(),
 162                 mesh(),
 163                 IOobject::NO_READ,
 164                 IOobject::NO_WRITE
 165             ),
 166             mesh(),
 167             GIIc_
 168         )
 169     );
 170 }
 171
 172
 173 void Foam::linearCohesiveLaw::damageTractions
 174 (
 175  scalar& tN,
 176  scalar& tS,
 177  const scalar deltaN,
 178  const scalar deltaS,
 179  const scalar GI,
 180  const scalar GII,
 181  const label faceID,
 182  const scalarField& globalPatchMaterials
 183  ) const
 184 {
 185   // Step 1
 186   // Calculate apparent initiation normal traction tNia
 187   // tNia is the initiation traction assuming that the ratio tN/tS
 188   // has stayed constant at the current value i.e.
 189   // tNia = B * tN
 190   // tSia = B * tS
 191   // where we term B the initiation traction constant.
 192
 193   // Calculate current B
 194   const scalar B = 1.0 /
 195     ::sqrt(
 196        (tN/sigmaMax_.value())*(tN/sigmaMax_.value())
 197        + (tS/tauMax_.value())*(tS/tauMax_.value())
 198        );
 199
 200   // Calculate apparent initiation tractions
 201   const scalar tNia = B * tN;
 202   const scalar tSia = B * tS;
 203
 204   // Step 2
 205   // Calculate apparent critical delta deltaCa
 206   // deltaCa would be the deltaC if the mode mix stayed constant at
 207   // the current value
 208   // GIca = C * GII
 209   // GIIca = C * GII
 210   // where we term C the mode-mix constant
 211
 212   // Calculate C
 213   const scalar C = 1.0 /
 214     ( (GI/GIc_.value()) + (GII/GIIc_.value()) );
 215
 216   // Calculate apparent critical energies
 217   const scalar GIca = C * GI;
 218   const scalar GIIca = C * GII;
 219
 220   // Calculate apparent critical deltas
 221   const scalar deltaNca = 2.0 * GIca / tNia;
 222   const scalar deltaSca = 2.0 * GIIca / tSia;
 223
 224   // Step 3
 225   // Set current tractions linearly decreasing from
 226   // tia to deltaCa
 227   // When the face is close to propagation the tN and
 228   // tS will get close to zero so we will limit
 229   // them to be positive as they are magnitudes
 230   tN = max( 0.0, tNia * ( 1.0 - (deltaN/deltaNca) ) );
 231   tS = max( 0.0, tSia * ( 1.0 - (deltaS/deltaSca) ) );
 232
 233   // Step 4
 234   // Hold the effective traction traction and
 235   // allow the mix of tractions to vary based on:
 236   // (tN/tS) = a * (deltaN/deltaS)
 237   // where a=1 is assumed here.
 238   // This essentially allows the mode-mixity
 239   // to vary depending on the deltas
 240   tN = (
 241     (sigmaMax_.value()/B) * deltaN /
 242     (SMALL + ::sqrt(
 243             (deltaN*deltaN)
 244             + (deltaS*deltaS)*(sigmaMax_.value()*sigmaMax_.value()
 245                                /(tauMax_.value()*tauMax_.value()))
 246             ))
 247     );
 248   tS = (
 249     (tauMax_.value()/B) * deltaS /
 250     (SMALL + ::sqrt(
 251             (deltaS*deltaS)
 252             + (deltaN*deltaN)*(tauMax_.value()*tauMax_.value()
 253                                /(sigmaMax_.value()*sigmaMax_.value()))
 254             ))
 255     );
 256 }
 257
 258
 259 Foam::tmp<Foam::surfaceVectorField>
 260 Foam::linearCohesiveLaw::interfaceTraction
 261 (
 262  surfaceVectorField n,
 263  volVectorField U,
 264  volTensorField gradU,
 265  volScalarField mu,
 266  volScalarField lambda
 267  ) const
 268 {
 269   notImplemented(type() + "::interfaceTraction()");
 270
 271     tmp<surfaceVectorField> tresult
 272     (
 273         new surfaceVectorField
 274         (
 275             IOobject
 276             (
 277             "interfaceTraction",
 278             mesh().time().timeName(),
 279             mesh(),
 280             IOobject::NO_READ,
 281             IOobject::NO_WRITE
 282         ),
 283             mesh(),
 284             dimensionedVector("zero", dimForce/dimArea, vector(0, 0, 0))
 285     )
 286     );
 287
 288     return tresult;
 289 }
 290
 291 // ************************************************************************* //