src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.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 "solidCohesiveFvPatchVectorField.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "volFields.H"
  29 #include "constitutiveModel.H"
  30 #include "regionSplit.H"
  31 #include "crackerFvMesh.H"
  32 #include "tractionBoundaryGradient.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  40
  41 solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField
  42 (
  43     const fvPatch& p,
  44     const DimensionedField<vector, volMesh>& iF
  45 )
  46 :
  47     directionMixedFvPatchVectorField(p, iF),
  48     fieldName_("undefined"),
  49     relaxationFactor_(1.0),
  50     traction_(p.size(), vector::zero),
  51     //minUnloadingSeparationDistance_(0.0),
  52     contact_(false),
  53     cracked_(0,false),
  54     curTractionN_(0, 0.0),
  55     oldTractionN_(0, 0.0),
  56     curTractionS_(0, 0.0),
  57     oldTractionS_(0, 0.0),
  58     deltaN_(0, 0.0),
  59     oldDeltaN_(0, 0.0),
  60     deltaS_(0, 0.0),
  61     oldDeltaS_(0, 0.0),
  62     unloadingDeltaEff_(0, 0.0),
  63     currentGI_(0, 0.0),
  64     oldGI_(0, 0.0),
  65     currentGII_(0, 0.0),
  66     oldGII_(0, 0.0),
  67     curTimeIndex_(-1),
  68     penaltyFactorPtr_(NULL),
  69     penaltyScale_(0.0),
  70     frictionCoeff_(0.0),
  71     explicitSeparationDistance_(false),
  72     curDeltaN_(0, 0.0),
  73     curDeltaS_(0, 0.0),
  74     updateGlobalPatchMaterials_(true),
  75     globalPatchMaterials_(0, 0.0),
  76     nonLinear_(nonLinearGeometry::OFF),
  77     orthotropic_(false)
  78 {}
  79
  80
  81 solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField
  82 (
  83     const fvPatch& p,
  84     const DimensionedField<vector, volMesh>& iF,
  85     const dictionary& dict
  86 )
  87 :
  88     directionMixedFvPatchVectorField(p, iF),
  89     fieldName_(dimensionedInternalField().name()),
  90     relaxationFactor_(readScalar(dict.lookup("relaxationFactor"))),
  91     traction_(p.size(), vector::zero),
  92     contact_(dict.lookup("contact")),
  93     cracked_(p.size(),false),
  94     curTractionN_(p.size(), 0.0),
  95     oldTractionN_(p.size(), 0.0),
  96     curTractionS_(p.size(), 0.0),
  97     oldTractionS_(p.size(), 0.0),
  98     deltaN_(p.size(), 0.0),
  99     oldDeltaN_(p.size(), 0.0),
 100     deltaS_(p.size(), 0.0),
 101     oldDeltaS_(p.size(), 0.0),
 102     unloadingDeltaEff_(p.size(), 0.0),
 103     currentGI_(p.size(), 0.0),
 104     oldGI_(p.size(), 0.0),
 105     currentGII_(p.size(), 0.0),
 106     oldGII_(p.size(), 0.0),
 107     curTimeIndex_(-1),
 108     penaltyFactorPtr_(NULL),
 109     penaltyScale_(readScalar(dict.lookup("penaltyScale"))),
 110     frictionCoeff_(readScalar(dict.lookup("frictionCoeff"))),
 111     explicitSeparationDistance_(dict.lookup("explicitSeparationDistance")),
 112     curDeltaN_(p.size(), 0.0),
 113     curDeltaS_(p.size(), 0.0),
 114     updateGlobalPatchMaterials_(true),
 115     globalPatchMaterials_(0, 0.0),
 116     nonLinear_(nonLinearGeometry::OFF),
 117     orthotropic_(false)
 118 {
 119     Info<< "Creating solidCohesive patch" << nl
 120         << "\tOnly Dugdale law currently available!" << endl;
 121
 122     // Check if traction boundary is for non linear solver
 123     if (dict.found("nonLinear"))
 124     {
 125         nonLinear_ = nonLinearGeometry::nonLinearNames_.read
 126         (
 127             dict.lookup("nonLinear")
 128         );
 129
 130         if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN)
 131         {
 132             Info << "\tnonLinear set to updated Lagrangian"
 133                 << endl;
 134         }
 135         else if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN)
 136         {
 137             Info << "\tnonLinear set to total Lagrangian"
 138                 << endl;
 139         }
 140     }
 141
 142     if (dict.found("orthotropic"))
 143     {
 144         orthotropic_ = Switch(dict.lookup("orthotropic"));
 145         Info << "\t\torthotropic set to " << orthotropic_ << endl;
 146     }
 147
 148     // if (minUnloadingSeparationDistance_ < 0.001)
 149     // {
 150     //     minUnloadingSeparationDistance_ = 0.001;
 151     // }
 152
 153     // if (minUnloadingSeparationDistance_ > 1.0)
 154     // {
 155     //     minUnloadingSeparationDistance_ = 1.0;
 156     // }
 157
 158     // Info << "minUnloadingSeparationDistance: "
 159     //     << minUnloadingSeparationDistance_ << endl;
 160
 161     if (dict.found("refValue"))
 162     {
 163         this->refValue() = vectorField("refValue", dict, p.size());
 164     }
 165     else
 166     {
 167         this->refValue() = vector::zero;
 168     }
 169
 170     if (dict.found("refGradient"))
 171     {
 172         this->refGrad() = vectorField("refGradient", dict, p.size());
 173     }
 174     else
 175     {
 176         this->refGrad() = vector::zero;
 177     }
 178
 179     if (dict.found("valueFraction"))
 180     {
 181         this->valueFraction() =
 182             symmTensorField("valueFraction", dict, p.size());
 183     }
 184     else
 185     {
 186         this->valueFraction() = symmTensor::zero;
 187     }
 188
 189     if (dict.found("value"))
 190     {
 191         vectorField::operator=(vectorField("value", dict, p.size()));
 192     }
 193     else
 194     {
 195         vectorField normalValue = transform(valueFraction(), refValue());
 196
 197         vectorField gradValue =
 198             this->patchInternalField() + refGrad()/this->patch().deltaCoeffs();
 199
 200         vectorField transformGradValue =
 201             transform(I - valueFraction(), gradValue);
 202
 203         vectorField::operator=(normalValue + transformGradValue);
 204     }
 205
 206     if (dict.found("traction"))
 207     {
 208         traction_ =
 209             vectorField("traction", dict, p.size());
 210     }
 211
 212     // if (dict.found("unloadingSeparationDistance"))
 213     // {
 214     //     unloadingSeparationDistance_ =
 215     //         vectorField("unloadingSeparationDistance", dict, p.size());
 216     // }
 217
 218     if (dict.found("cracked"))
 219     {
 220 //      cracked_ = Field<bool>("cracked", dict, p.size());
 221         cracked_ = Field<scalar>("cracked", dict, p.size());
 222     }
 223
 224     if (dict.found("curTractionN"))
 225     {
 226         curTractionN_ = scalarField("curTractionN", dict, p.size());
 227     }
 228     if (dict.found("curTractionS"))
 229     {
 230         curTractionS_ = scalarField("curTractionS", dict, p.size());
 231     }
 232     if (dict.found("oldTractionN"))
 233     {
 234         oldTractionN_ = scalarField("oldTractionN", dict, p.size());
 235     }
 236     if (dict.found("oldTractionS"))
 237     {
 238         oldTractionS_ = scalarField("oldTractionS", dict, p.size());
 239     }
 240
 241     if (dict.found("deltaN"))
 242     {
 243         deltaN_ = scalarField("deltaN", dict, p.size());
 244     }
 245     if (dict.found("deltaS"))
 246     {
 247         deltaS_ = scalarField("deltaS", dict, p.size());
 248     }
 249     if (dict.found("oldDeltaN"))
 250     {
 251         oldDeltaN_ = scalarField("oldDeltaN", dict, p.size());
 252     }
 253     if (dict.found("oldDeltaS"))
 254     {
 255         oldDeltaS_ = scalarField("oldDeltaS", dict, p.size());
 256     }
 257     if (dict.found("curDeltaN"))
 258     {
 259         curDeltaN_ = scalarField("curDeltaN", dict, p.size());
 260     }
 261     if (dict.found("curDeltaS"))
 262     {
 263         curDeltaS_ = scalarField("curDeltaS", dict, p.size());
 264     }
 265     if (dict.found("unloadingDeltaEff"))
 266     {
 267         unloadingDeltaEff_ = scalarField("unloadingDeltaEff", dict, p.size());
 268     }
 269
 270     if (dict.found("currentGI"))
 271     {
 272         currentGI_ = scalarField("currentGI", dict, p.size());
 273     }
 274     if (dict.found("currentGII"))
 275     {
 276         currentGII_ = scalarField("currentGII", dict, p.size());
 277     }
 278     if (dict.found("oldGI"))
 279     {
 280         oldGI_ = scalarField("oldGI", dict, p.size());
 281     }
 282     if (dict.found("oldGII"))
 283     {
 284         oldGII_ = scalarField("oldGII", dict, p.size());
 285     }
 286 }
 287
 288
 289 solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField
 290 (
 291     const solidCohesiveFvPatchVectorField& cpf
 292 )
 293 :
 294     directionMixedFvPatchVectorField(cpf),
 295     fieldName_(cpf.fieldName_),
 296     relaxationFactor_(cpf.relaxationFactor_),
 297     traction_(cpf.traction_),
 298     //minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 299     contact_(cpf.contact_),
 300     cracked_(cpf.cracked_),
 301     curTractionN_(cpf.curTractionN_),
 302     oldTractionN_(cpf.oldTractionN_),
 303     curTractionS_(cpf.curTractionS_),
 304     oldTractionS_(cpf.oldTractionS_),
 305     deltaN_(cpf.deltaN_),
 306     oldDeltaN_(cpf.oldDeltaN_),
 307     deltaS_(cpf.deltaS_),
 308     oldDeltaS_(cpf.oldDeltaS_),
 309     unloadingDeltaEff_(cpf.unloadingDeltaEff_),
 310     currentGI_(cpf.currentGI_),
 311     oldGI_(cpf.oldGI_),
 312     currentGII_(cpf.currentGII_),
 313     oldGII_(cpf.oldGII_),
 314     curTimeIndex_(-1),
 315     penaltyFactorPtr_(cpf.penaltyFactorPtr_),
 316     penaltyScale_(cpf.penaltyScale_),
 317     frictionCoeff_(cpf.frictionCoeff_),
 318     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 319     curDeltaN_(cpf.curDeltaN_),
 320     curDeltaS_(cpf.curDeltaS_),
 321     updateGlobalPatchMaterials_(cpf.updateGlobalPatchMaterials_),
 322     globalPatchMaterials_(cpf.globalPatchMaterials_),
 323     nonLinear_(cpf.nonLinear_),
 324     orthotropic_(cpf.orthotropic_)
 325 {}
 326
 327
 328 solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField
 329 (
 330     const solidCohesiveFvPatchVectorField& cpf,
 331     const fvPatch& p,
 332     const DimensionedField<vector, volMesh>& iF,
 333     const fvPatchFieldMapper& mapper
 334 )
 335 :
 336     directionMixedFvPatchVectorField(cpf, p, iF, mapper),
 337     fieldName_(cpf.fieldName_),
 338     relaxationFactor_(cpf.relaxationFactor_),
 339     traction_(cpf.traction_, mapper),
 340     //minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 341     contact_(cpf.contact_),
 342     cracked_(cpf.cracked_),
 343     curTractionN_(cpf.curTractionN_),
 344     oldTractionN_(cpf.oldTractionN_),
 345     curTractionS_(cpf.curTractionS_),
 346     oldTractionS_(cpf.oldTractionS_),
 347     deltaN_(cpf.deltaN_),
 348     oldDeltaN_(cpf.oldDeltaN_),
 349     deltaS_(cpf.deltaS_),
 350     oldDeltaS_(cpf.oldDeltaS_),
 351     unloadingDeltaEff_(cpf.unloadingDeltaEff_),
 352     currentGI_(cpf.currentGI_),
 353     oldGI_(cpf.oldGI_),
 354     currentGII_(cpf.currentGII_),
 355     oldGII_(cpf.oldGII_),
 356     curTimeIndex_(-1),
 357     penaltyFactorPtr_(cpf.penaltyFactorPtr_),
 358     penaltyScale_(cpf.penaltyScale_),
 359     frictionCoeff_(cpf.frictionCoeff_),
 360     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 361     curDeltaN_(cpf.curDeltaN_),
 362     curDeltaS_(cpf.curDeltaS_),
 363     updateGlobalPatchMaterials_(cpf.updateGlobalPatchMaterials_),
 364     globalPatchMaterials_(cpf.globalPatchMaterials_),
 365     nonLinear_(cpf.nonLinear_),
 366     orthotropic_(cpf.orthotropic_)
 367 {}
 368
 369
 370 solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField
 371 (
 372     const solidCohesiveFvPatchVectorField& cpf,
 373     const DimensionedField<vector, volMesh>& iF
 374 )
 375 :
 376     directionMixedFvPatchVectorField(cpf, iF),
 377     fieldName_(cpf.fieldName_),
 378     relaxationFactor_(cpf.relaxationFactor_),
 379     traction_(cpf.traction_),
 380     //minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 381     contact_(cpf.contact_),
 382     cracked_(cpf.cracked_),
 383     curTractionN_(cpf.curTractionN_),
 384     oldTractionN_(cpf.oldTractionN_),
 385     curTractionS_(cpf.curTractionS_),
 386     oldTractionS_(cpf.oldTractionS_),
 387     deltaN_(cpf.deltaN_),
 388     oldDeltaN_(cpf.oldDeltaN_),
 389     deltaS_(cpf.deltaS_),
 390     oldDeltaS_(cpf.oldDeltaS_),
 391     unloadingDeltaEff_(cpf.unloadingDeltaEff_),
 392     currentGI_(cpf.currentGI_),
 393     oldGI_(cpf.oldGI_),
 394     currentGII_(cpf.currentGII_),
 395     oldGII_(cpf.oldGII_),
 396     curTimeIndex_(-1),
 397     penaltyFactorPtr_(cpf.penaltyFactorPtr_),
 398     penaltyScale_(cpf.penaltyScale_),
 399     frictionCoeff_(cpf.frictionCoeff_),
 400     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 401     curDeltaN_(cpf.curDeltaN_),
 402     curDeltaS_(cpf.curDeltaS_),
 403     updateGlobalPatchMaterials_(cpf.updateGlobalPatchMaterials_),
 404     globalPatchMaterials_(cpf.globalPatchMaterials_),
 405     nonLinear_(cpf.nonLinear_),
 406     orthotropic_(cpf.orthotropic_)
 407 {}
 408
 409
 410 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 411
 412 tmp<scalarField> solidCohesiveFvPatchVectorField::crackingAndDamage() const
 413 {
 414     tmp<scalarField> tcrackingAndDamage
 415     (
 416         new scalarField(size(), 1.0)
 417     );
 418     scalarField& cad = tcrackingAndDamage();
 419
 420     forAll(cad, facei)
 421     {
 422         if (cracked_[facei])
 423         {
 424             cad[facei] = 2.0;
 425         }
 426     }
 427
 428     return tcrackingAndDamage;
 429 }
 430
 431
 432 tmp<scalarField> solidCohesiveFvPatchVectorField::GI() const
 433 {
 434     tmp<scalarField> tGI
 435     (
 436         new scalarField(size(), 0.0)
 437     );
 438     scalarField& GI = tGI();
 439
 440     forAll(GI, facei)
 441     {
 442         GI[facei] = currentGI_[facei];
 443     }
 444
 445     return tGI;
 446 }
 447
 448
 449 tmp<scalarField> solidCohesiveFvPatchVectorField::GII() const
 450 {
 451     tmp<scalarField> tGII
 452     (
 453         new scalarField(size(), 0.0)
 454     );
 455     scalarField& GII = tGII();
 456
 457     forAll(GII, facei)
 458     {
 459         GII[facei] = currentGII_[facei];
 460     }
 461
 462     return tGII;
 463 }
 464
 465
 466 bool solidCohesiveFvPatchVectorField::cracking()
 467 {
 468     const fvMesh& mesh = patch().boundaryMesh().mesh();
 469     const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh);
 470
 471     // global fields
 472     Field<scalar> globalCracked =
 473         crackerMesh.globalCrackField(cracked_);
 474
 475     label sumDamaged = 0;
 476     label sumCracked = 0;
 477
 478     forAll(globalCracked, facei)
 479     {
 480         if (globalCracked[facei] > 0.0)
 481         {
 482             sumCracked++;
 483         }
 484         else
 485         {
 486             sumDamaged++;
 487         }
 488     }
 489     Info<< "\t\tThere are " << sumDamaged << " damaged faces" << nl
 490         << "\t\tThere are " << sumCracked << " cracked faces" << endl;
 491
 492     return false;
 493 }
 494
 495
 496 void solidCohesiveFvPatchVectorField::autoMap
 497 (
 498     const fvPatchFieldMapper& m
 499 )
 500 {
 501     directionMixedFvPatchVectorField::autoMap(m);
 502
 503     traction_.autoMap(m);
 504     cracked_.autoMap(m);
 505     curTractionN_.autoMap(m);
 506     oldTractionN_.autoMap(m);
 507     curTractionS_.autoMap(m);
 508     oldTractionS_.autoMap(m);
 509     deltaN_.autoMap(m);
 510     oldDeltaN_.autoMap(m);
 511     deltaS_.autoMap(m);
 512     oldDeltaS_.autoMap(m);
 513     unloadingDeltaEff_.autoMap(m);
 514     currentGI_.autoMap(m);
 515     oldGI_.autoMap(m);
 516     currentGII_.autoMap(m);
 517     oldGII_.autoMap(m);
 518     curDeltaN_.autoMap(m);
 519     curDeltaS_.autoMap(m);
 520     globalPatchMaterials_.autoMap(m);
 521
 522     // force this field to be reset next time updateCoeffs is called
 523     updateGlobalPatchMaterials_ = true;
 524
 525     // Must use primitive mesh data because fvMesh data is packed in fields
 526     // and cannot be accessed during mapping.  HJ, 12/Dec/2008
 527     vectorField n = patch().patch().faceNormals();
 528
 529     const labelList& addressing = m.directAddressing();
 530
 531     label nNewFaces = m.size() - m.sizeBeforeMapping();
 532
 533     // philipc:
 534     // when a new face is created in the crack, we need to set
 535     // its traction, but how do we do it:
 536     // we calculate the traciton in the solver (updatedCrack.H)
 537     // but we cannot access that field here during mapping
 538     // So in the solver we will have to cast this patch and update
 539     // the traction_ field after mapping and call updateCoeffs.
 540     // For here, we will just set traction to zero.
 541
 542     if ( (patch().size() == 1) && (nNewFaces == 1) )
 543     {
 544         label i=0;
 545         this->valueFraction()[i] = symmTensor::zero;
 546
 547         traction_[i] = vector::zero;
 548         cracked_[i] = false;
 549         curTractionN_[i] = 0.0;
 550         oldTractionN_[i] = 0.0;
 551         curTractionS_[i] = 0.0;
 552         oldTractionS_[i] = 0.0;
 553         deltaN_[i] = 0.0;
 554         oldDeltaN_[i] = 0.0;
 555         deltaS_[i] = 0.0;
 556         oldDeltaS_[i] = 0.0;
 557         curDeltaN_[i] = 0.0;
 558         curDeltaS_[i] = 0.0;
 559         unloadingDeltaEff_[i] = 0.0;
 560         currentGI_[i] = 0.0;
 561         oldGI_[i] = 0.0;
 562         currentGII_[i] = 0.0;
 563         oldGII_[i] = 0.0;
 564     }
 565     else if ( (patch().size() == 2) && (nNewFaces == 1) )
 566     {
 567         label i=1;
 568         this->valueFraction()[i] = symmTensor::zero;
 569         traction_[i] = vector::zero;
 570         cracked_[i] = false;
 571
 572         curTractionN_[i] = 0.0;
 573         oldTractionN_[i] = 0.0;
 574         curTractionS_[i] = 0.0;
 575         oldTractionS_[i] = 0.0;
 576         deltaN_[i] = 0.0;
 577         oldDeltaN_[i] = 0.0;
 578         deltaS_[i] = 0.0;
 579         oldDeltaS_[i] = 0.0;
 580         curDeltaN_[i] = 0.0;
 581         curDeltaS_[i] = 0.0;
 582         unloadingDeltaEff_[i] = 0.0;
 583         currentGI_[i] = 0.0;
 584         oldGI_[i] = 0.0;
 585         currentGII_[i] = 0.0;
 586         oldGII_[i] = 0.0;
 587     }
 588     else if ( (patch().size()==2) && (nNewFaces == 2) )
 589     {
 590         label i=0;
 591         this->valueFraction()[i] = symmTensor::zero;
 592         traction_[i] = vector::zero;
 593         i=1;
 594         this->valueFraction()[i] = symmTensor::zero;
 595         traction_[i] = vector::zero;
 596
 597         cracked_[i] = false;
 598         curTractionN_[i] = 0.0;
 599         oldTractionN_[i] = 0.0;
 600         curTractionS_[i] = 0.0;
 601         oldTractionS_[i] = 0.0;
 602         deltaN_[i] = 0.0;
 603         oldDeltaN_[i] = 0.0;
 604         deltaS_[i] = 0.0;
 605         oldDeltaS_[i] = 0.0;
 606         curDeltaN_[i] = 0.0;
 607         curDeltaS_[i] = 0.0;
 608         unloadingDeltaEff_[i] = 0.0;
 609         currentGI_[i] = 0.0;
 610         oldGI_[i] = 0.0;
 611         currentGII_[i] = 0.0;
 612         oldGII_[i] = 0.0;
 613     }
 614     else
 615     {
 616         for (label i = 1; i < patch().size(); i++)
 617         {
 618             if (addressing[i] == 0)
 619             {
 620                 this->valueFraction()[i] = symmTensor::zero;
 621                 traction_[i] = vector::zero;
 622
 623                 cracked_[i] = false;
 624                 curTractionN_[i] = 0.0;
 625                 oldTractionN_[i] = 0.0;
 626                 curTractionS_[i] = 0.0;
 627                 oldTractionS_[i] = 0.0;
 628                 deltaN_[i] = 0.0;
 629                 oldDeltaN_[i] = 0.0;
 630                 deltaS_[i] = 0.0;
 631                 oldDeltaS_[i] = 0.0;
 632                 curDeltaN_[i] = 0.0;
 633                 curDeltaS_[i] = 0.0;
 634                 unloadingDeltaEff_[i] = 0.0;
 635                 currentGI_[i] = 0.0;
 636                 oldGI_[i] = 0.0;
 637                 currentGII_[i] = 0.0;
 638                 oldGII_[i] = 0.0;
 639             }
 640         }
 641     }
 642 }
 643
 644
 645 // Reverse-map the given fvPatchField onto this fvPatchField
 646 void solidCohesiveFvPatchVectorField::rmap
 647 (
 648     const fvPatchVectorField& ptf,
 649     const labelList& addr
 650 )
 651 {
 652     directionMixedFvPatchVectorField::rmap(ptf, addr);
 653
 654     const solidCohesiveFvPatchVectorField& dmptf =
 655         refCast<const solidCohesiveFvPatchVectorField>(ptf);
 656
 657     relaxationFactor_ = dmptf.relaxationFactor_;
 658     traction_.rmap(dmptf.traction_, addr);
 659     contact_ = dmptf.contact_;
 660     cracked_ = dmptf.cracked_;
 661     curTractionN_ = dmptf.curTractionN_;
 662     oldTractionN_ = dmptf.oldTractionN_;
 663     curTractionS_ = dmptf.curTractionS_;
 664     oldTractionS_ = dmptf.oldTractionS_;
 665     deltaN_ = dmptf.deltaN_;
 666     oldDeltaN_ = dmptf.oldDeltaN_;
 667     deltaS_ = dmptf.deltaS_;
 668     oldDeltaS_ = dmptf.oldDeltaS_;
 669     curDeltaN_ = dmptf.curDeltaN_;
 670     curDeltaS_ = dmptf.curDeltaS_;
 671     globalPatchMaterials_ = dmptf.globalPatchMaterials_;
 672     unloadingDeltaEff_ = dmptf.unloadingDeltaEff_;
 673     currentGI_ = dmptf.currentGI_;
 674     currentGII_ = dmptf.currentGII_;
 675 }
 676
 677
 678 // Update the coefficients associated with the patch field
 679 void solidCohesiveFvPatchVectorField::updateCoeffs()
 680 {
 681     if (updated())
 682     {
 683         return;
 684     }
 685
 686     vectorField n = patch().nf();
 687
 688     // lookup cohesive properties from rheology
 689     const constitutiveModel& rheology =
 690         this->db().objectRegistry::lookupObject<constitutiveModel>
 691         (
 692             "rheologyProperties"
 693         );
 694
 695     label patchID = patch().index();
 696     const scalarField sigmaMax =
 697         rheology.cohLaw().sigmaMax()().boundaryField()[patchID];
 698     const scalarField tauMax =
 699         rheology.cohLaw().tauMax()().boundaryField()[patchID];
 700     const scalarField GIc =
 701         rheology.cohLaw().GIc()().boundaryField()[patchID];
 702     const scalarField GIIc =
 703         rheology.cohLaw().GIIc()().boundaryField()[patchID];
 704
 705     // update old values
 706     if (curTimeIndex_ != this->db().time().timeIndex())
 707     {
 708         // we force the penalty factor to be calculated here for the first time
 709         // as all processors must call this at the same time
 710         if (contact_ && patch().size() > 0)
 711         {
 712             // force calculation of penalty factor here
 713             penaltyFactor();
 714         }
 715
 716         oldDeltaN_ = deltaN_;
 717         oldDeltaS_ = deltaS_;
 718         // curDelta is always latest delta
 719         // whereas delta is only latest when explicitDist is false
 720         //oldDeltaN_ = curDeltaN_;
 721         //oldDeltaS_ = curDeltaS_;
 722         oldTractionN_ = curTractionN_;
 723         oldTractionS_ = curTractionS_;
 724         oldGI_ = currentGI_;
 725         oldGII_ = currentGII_;
 726         unloadingDeltaEff_ =
 727             max
 728             (
 729                 unloadingDeltaEff_,
 730                 Foam::sqrt
 731                 (
 732                     sqr(max(deltaN_, scalar(0))) + sqr(deltaS_)
 733                 )
 734             );
 735
 736         curTimeIndex_ = this->db().time().timeIndex();
 737     }
 738
 739
 740     // Get face cells regions
 741     const unallocLabelList& faceCells = patch().faceCells();
 742     const fvMesh& mesh = patch().boundaryMesh().mesh();
 743
 744     if (!isA<crackerFvMesh>(mesh))
 745     {
 746         FatalErrorIn
 747         (
 748             "void solidCohesiveFvPatchVectorField::updateCoeffs() const"
 749         )   << "Mesh should be of type: " << crackerFvMesh::typeName
 750             << abort(FatalError);
 751     }
 752
 753     const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh);
 754
 755     // global patch material field is needed for multimaterial cohesive laws
 756     if (updateGlobalPatchMaterials_)
 757     {
 758         updateGlobalPatchMaterials_ = false;
 759
 760         if (mesh.objectRegistry::foundObject<volScalarField>("materials"))
 761         {
 762             scalarField localPatchMaterials =
 763                 patch().lookupPatchField<volScalarField, scalar>
 764                 ("materials").patchInternalField();
 765             scalarField newGlobalPatchMaterials =
 766                 crackerMesh.globalCrackField(localPatchMaterials);
 767
 768             globalPatchMaterials_.setSize(newGlobalPatchMaterials.size());
 769             globalPatchMaterials_ = newGlobalPatchMaterials;
 770         }
 771     }
 772
 773     const regionSplit& regions = crackerMesh.regions();
 774
 775     labelField faceCellRegion(size(), -1);
 776
 777     forAll(faceCellRegion, faceI)
 778     {
 779         faceCellRegion[faceI] = regions[faceCells[faceI]];
 780     }
 781
 782     labelField globalFaceCellRegion =
 783         crackerMesh.globalCrackField(faceCellRegion);
 784
 785     // Patch displacement
 786     vectorField UPatch = *this;
 787     if (fieldName_ == "DU")
 788     {
 789         UPatch += patch().lookupPatchField<volVectorField, vector>("U");
 790     }
 791
 792     // Global displacement
 793     vectorField globalUPatch = crackerMesh.globalCrackField(UPatch);
 794
 795
 796     // Update separation distance
 797     vectorField delta(patch().size(), vector::zero);
 798     const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 799     label globalIndex = crackerMesh.localCrackStart();
 800     for (label i = 0; i < patch().size(); i++)
 801     {
 802         // newSeparationDistance[i] =
 803         delta[i] =
 804             globalUPatch[gcfa[globalIndex]]
 805           - globalUPatch[globalIndex];
 806
 807         globalIndex++;
 808     }
 809
 810     //globalIndex = crackerMesh.localCrackStart();
 811     for (label i = 0; i < patch().size(); i++)
 812     {
 813         // update deltas
 814         curDeltaN_[i] = n[i] & delta[i];
 815         curDeltaS_[i] = mag( (I - sqr(n[i])) & delta[i]); // shearing
 816         if (explicitSeparationDistance_)
 817         {
 818             deltaN_[i] = oldDeltaN_[i];
 819             deltaS_[i] = oldDeltaS_[i];
 820
 821             if (mag(deltaN_[i]) < SMALL)
 822             {
 823                 deltaN_[i] = 2*SMALL;
 824             }
 825             FatalError
 826                 << "explicitSeparationDistance_ is true!" << abort(FatalError);
 827         }
 828         else
 829         {
 830             // relax
 831             deltaN_[i] =
 832                 relaxationFactor_*curDeltaN_[i]
 833               + (1.0 - relaxationFactor_)*deltaN_[i];
 834
 835             deltaS_[i] =
 836                 relaxationFactor_*curDeltaS_[i]
 837               + (1.0 - relaxationFactor_)*deltaS_[i];
 838         }
 839
 840         // update tractions
 841         curTractionN_[i] = n[i] & traction_[i];
 842         curTractionS_[i] = mag( (I-sqr(n[i])) & traction_[i] );
 843
 844         // effective delta - only positive deltaN is considered
 845         const scalar deltaEff =
 846             Foam::sqrt(max(deltaN_[i], 0.0)*max(deltaN_[i], 0.0)
 847           + deltaS_[i]*deltaS_[i]);
 848
 849         // update energies
 850         // stop calculating after cracking for convergence
 851         // because crack might jump in and out of damaged/failed
 852         // only update energies if there is loading
 853         if ( !cracked_[i] && (deltaEff > (unloadingDeltaEff_[i]-SMALL)) )
 854         {
 855             // if the average normal stress is tensile
 856             if ((curTractionN_[i]+oldTractionN_[i]) > 0.0)
 857             {
 858                 // trapezoidal rule
 859                 currentGI_[i] =
 860                     oldGI_[i]
 861                   + ((0.5*(curTractionN_[i]+oldTractionN_[i]))*
 862                     (deltaN_[i]-oldDeltaN_[i]));
 863             }
 864             else
 865             {
 866                 // no modeI energy lost if in compression
 867                 currentGI_[i] = oldGI_[i];
 868             }
 869
 870             // mode II - trapezoidal rule
 871             currentGII_[i] = oldGII_[i]
 872               + ((0.5*(curTractionS_[i]+oldTractionS_[i]))*
 873                 (deltaS_[i]-oldDeltaS_[i]));
 874         }
 875         //scalar currentG = currentGI_[i] + currentGII_[i];
 876
 877         // if
 878         // (
 879         //  ( globalFaceCellRegion[globalIndex]
 880         //    == globalFaceCellRegion[gcfa[globalIndex]] )
 881         //  ||
 882         //  !cracked_[i]
 883         //  )
 884
 885         // new tractions to be calculated
 886         scalar curNormalTraction = 0;
 887         vector curTangentialTraction = vector::zero;
 888
 889         // loading
 890         // if the effective delta is greater than unloadingDeltaEff then
 891         // there is loading
 892         // unloadingDeltaEff is the maximum previsouly reached deltaEff
 893         if (deltaEff > (unloadingDeltaEff_[i]-SMALL) )
 894         {
 895             // at the moment loading and unloading are the same
 896             // if total energy is dissipated, then fully crack face
 897             //if ( currentG > GIc[i] || cracked_[i] ) //Gc )
 898             // propagation
 899             if ( ((currentGI_[i]/GIc[i]) + (currentGII_[i]/GIIc[i])) >= 1 )
 900             {
 901                 //Pout << "GIc[i] is " << GIc[i] << ", curG is " << currentG << endl;
 902                 if (!cracked_[i])
 903                 {
 904                     Pout << "Face " << i << " is fully cracked" << endl;
 905                 }
 906
 907                 cracked_[i] = true;
 908
 909
 910                 // failed faces might come in to contact so we need to deal with them
 911                 // here we could use the face delta values but that assumes that the
 912                 // opposing crack faces are aligned direclty opposite one another, which
 913                 // in general they are not. So we actually need to calculate the actual
 914                 // face penetration distances and then a standard penalty approach
 915                 // contact is probably fine. We will use simple assumption for now which
 916                 // is fine is the relative displacement
 917                 // of the opposing faces is not too large
 918                 // To-do: we must calculate actual distances
 919                 curNormalTraction = 0.0;
 920                 curTangentialTraction = vector::zero;
 921                 if ( contact_ && deltaN_[i] <= 0.0 )
 922                 {
 923                     curNormalTraction = deltaN_[i]*penaltyFactor();
 924                     //Info << "penaltyFactor() is " << penaltyFactor() << endl;
 925
 926                     // friction
 927                     vector sDir = (I - sqr(n[i]))&delta[i];
 928                     sDir /= mag(sDir+vector(SMALL,SMALL,SMALL));
 929                     scalar slip = mag(deltaS_[gcfa[i]] - deltaS_[i]);
 930                     scalar fricMag =
 931                         min
 932                         (
 933                             slip*penaltyFactor(),
 934                             frictionCoeff_*curNormalTraction // stick/slip
 935                         );
 936                     curTangentialTraction = fricMag*sDir;
 937                 }
 938
 939                 // relax tractions
 940                 curNormalTraction =
 941                     relaxationFactor_*curNormalTraction
 942                   + (1-relaxationFactor_)*(n[i] & traction_[i]);
 943
 944                 curTangentialTraction =
 945                     relaxationFactor_*curTangentialTraction
 946                   + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]);
 947
 948                 traction_[i] = curNormalTraction*n[i] + curTangentialTraction;
 949             }
 950             // damging face with positive normal delta
 951             else if ( deltaN_[i] > 0.0 )
 952             {
 953                 if (cracked_[i])
 954                 {
 955                     Pout << "Face " << i << " is un-cracked" << endl;
 956                 }
 957
 958                 cracked_[i] = false;
 959
 960                 // set traction in a fixed point iteration manner to force
 961                 // (tN/sigmaMax)^2 + (tS/tauMax)^2 = 1
 962                 // while also allowing varying mode mix by assuming colinear traction
 963                 // and delta
 964                 // Dugdale version
 965                 // scalar tN =
 966                 //   sigmaMax[i] * deltaN_[i] /
 967                 //   (SMALL + ::sqrt( (deltaN_[i]*deltaN_[i])
 968                 // + (deltaS_[i]*deltaS_[i])*(sigmaMax[i]*sigmaMax[i]
 969                 // /(tauMax[i]*tauMax[i])) ));
 970                 // scalar tS =
 971                 //   tauMax[i] * deltaS_[i] /
 972                 //   (SMALL + ::sqrt( (deltaS_[i]*deltaS_[i])
 973                 // + (deltaN_[i]*deltaN_[i])*(tauMax[i]*tauMax[i]
 974                 // /(sigmaMax[i]*sigmaMax[i])) ));
 975
 976                 // Given current deltaN and deltaS, the cohesive law
 977                 // (Dugdale, linear, etc)gives back current traction:
 978                 scalar tN = n[i] & traction_[i];
 979                 scalar tS = mag( (I - sqr(n[i])) & traction_[i] );
 980                 // Update tN and tS
 981                 rheology.cohLaw().damageTractions
 982                 (
 983                     tN,
 984                     tS,
 985                     deltaN_[i],
 986                     deltaS_[i],
 987                     currentGI_[i],
 988                     currentGII_[i],
 989                     i,
 990                     globalPatchMaterials_
 991                 );
 992
 993                 // shear delta direction
 994                 vector sDir = (I - sqr(n[i]))&delta[i];
 995                 sDir /= mag(sDir+vector(SMALL,SMALL,SMALL));
 996
 997                 // relax tractions
 998                 curNormalTraction =
 999                     relaxationFactor_*tN
1000                   + (1-relaxationFactor_)*(n[i] & traction_[i]);
1001
1002                 curTangentialTraction =
1003                     relaxationFactor_*(tS*sDir)
1004                   + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]);
1005
1006                 traction_[i] = curNormalTraction*n[i] + curTangentialTraction;
1007             }
1008             // damaging faces with negative normal delta
1009             else
1010             {
1011                 if (cracked_[i])
1012                 {
1013                     Pout << "Face " << i << " is un-cracked" << endl;
1014                 }
1015
1016                 cracked_[i] = false;
1017                 //Pout << "Contact and shearing face " << i << endl;
1018
1019                 //Pout << "Face " << i << " is damaging but in contact" << endl;
1020                 // set shear traction from cohesive law
1021                 // set normal traction to contact condition
1022                 scalar tS = tauMax[i];
1023
1024                 vector sDir = (I - sqr(n[i])) & delta[i];
1025                 sDir /= mag(sDir+vector(SMALL,SMALL,SMALL));
1026
1027                 // Simple penalty condition
1028                 scalar penaltyFac = 0.0;
1029                 if (contact_)
1030                 {
1031                     penaltyFac = penaltyFactor();
1032                 }
1033                 scalar contactTN = deltaN_[i]*penaltyFac;
1034
1035                 // relax tractions
1036                 curNormalTraction =
1037                     relaxationFactor_*contactTN
1038                   + (1-relaxationFactor_)*(n[i] & traction_[i]);
1039
1040                 curTangentialTraction =
1041                     relaxationFactor_*(tS*sDir)
1042                   + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]);
1043
1044                 traction_[i] = curNormalTraction*n[i] + curTangentialTraction;
1045             }
1046         }
1047         else
1048         {
1049             // unloading
1050             // as the current effective delta is less than the old time
1051             // effective delta
1052
1053             // We have two choice for unloading:
1054             // (a) ductile approach
1055             //         unload with initial stiffness (which is infinite)
1056             //         similar to ductilve metal stress-strain curve
1057             //         as we use infinite intial stiffness, this means to elastic
1058             //         energy is recovered
1059             //         and we immediately start loading in opposite direction
1060             // (b) brittle approach
1061             //         unload straight back to the origin
1062             //         this means we recover elastic energy
1063             //
1064             // approach (b) is numerically "nicer"; however, for Dugdale cohesive
1065             // zone this implys that only half the energy is dissipated just before
1066             // failure and then the other half is dissipated at failure. This may
1067             // not be an issue in practive but it does not really make sense.
1068             // Obviously it is fine for a linear cohesive zone as the energy
1069             // is smoothly dissipated up to failure. For now, we will implement
1070             // approach (b), but this requires more thinking...
1071
1072             // reduce traction linearly with the reduction in delta
1073             scalar scaleFactor = deltaEff/(unloadingDeltaEff_[i]);
1074             traction_[i] =
1075                 relaxationFactor_*scaleFactor*traction_[i]
1076               + (1 - relaxationFactor_)*traction_[i];
1077         }
1078     }
1079
1080     bool incremental(fieldName_ == "DU");
1081
1082     this->refGrad() = tractionBoundaryGradient::snGrad
1083     (
1084         traction_,
1085         scalarField(traction_.size(), 0.0),
1086         fieldName_,
1087         "U",
1088         patch(),
1089         orthotropic_,
1090         nonLinear_,
1091         incremental
1092     );
1093
1094     directionMixedFvPatchVectorField::updateCoeffs();
1095 }
1096
1097
1098 void solidCohesiveFvPatchVectorField::calcPenaltyFactor()
1099 {
1100     // calculate penalty factor similar to standardPenalty contact model
1101     // approx penaltyFactor from mechanical properties
1102     // this can then be scaled using the penaltyScale
1103     // to-do: write equivalent for orthotropic
1104     if (orthotropic_)
1105     {
1106         FatalErrorIn
1107         (
1108             "void solidCohesiveFvPatchVectorField::calcPenaltyFactor()"
1109         )   << "solidCohesiveFvPatchVectorField::calcPenaltyFactor()"
1110             << " has yet to be written for orthotropic"
1111             << exit(FatalError);
1112     }
1113
1114     const label patchID = patch().index();
1115     const fvMesh& mesh = patch().boundaryMesh().mesh();
1116
1117     const scalarField& mu =
1118         mesh.objectRegistry::lookupObject<volScalarField>
1119         ("mu").boundaryField()[patchID];
1120
1121     const scalarField& lambda =
1122         mesh.objectRegistry::lookupObject<volScalarField>
1123         ("lambda").boundaryField()[patchID];
1124
1125     // avarage contact patch bulk modulus
1126     scalar bulkModulus = gAverage(lambda + (2.0/3.0)*mu);
1127
1128     // average contact patch face area
1129     scalar faceArea = gAverage(mesh.magSf().boundaryField()[patchID]);
1130
1131     // average contact patch cell volume
1132     scalar cellVolume = 0.0;
1133
1134     const volScalarField::DimensionedInternalField & V = mesh.V();
1135     {
1136         const unallocLabelList& faceCells =
1137             mesh.boundary()[patchID].faceCells();
1138
1139         forAll(mesh.boundary()[patchID], facei)
1140         {
1141             cellVolume += V[faceCells[facei]];
1142         }
1143     }
1144     cellVolume /= patch().size();
1145
1146     // approximate penalty factor based on Hallquist et al.
1147     // we approximate penalty factor for traction instead of force
1148     penaltyFactorPtr_ =
1149         new scalar(penaltyScale_*bulkModulus*faceArea/cellVolume);
1150 }
1151
1152
1153 // Write
1154 void solidCohesiveFvPatchVectorField::write(Ostream& os) const
1155 {
1156     directionMixedFvPatchVectorField::write(os);
1157     traction_.writeEntry("traction", os);
1158     os.writeKeyword("relaxationFactor")
1159         << relaxationFactor_ << token::END_STATEMENT << nl;
1160
1161     cracked_.writeEntry("cracked", os);
1162     curTractionN_.writeEntry("curTractionN", os);
1163     curTractionS_.writeEntry("curTractionS", os);
1164     oldTractionN_.writeEntry("oldTractionN", os);
1165     oldTractionS_.writeEntry("oldTractionS", os);
1166     deltaN_.writeEntry("deltaN", os);
1167     deltaS_.writeEntry("deltaS", os);
1168     oldDeltaN_.writeEntry("oldDeltaN", os);
1169     oldDeltaS_.writeEntry("oldDeltaS", os);
1170     unloadingDeltaEff_.writeEntry("unloadingDeltaEff", os);
1171     currentGI_.writeEntry("currentGI", os);
1172     currentGII_.writeEntry("currentGII", os);
1173     oldGI_.writeEntry("oldGI", os);
1174     oldGII_.writeEntry("oldGII", os);
1175
1176     os.writeKeyword("curTimeIndex")
1177         << curTimeIndex_ << token::END_STATEMENT << nl;
1178     os.writeKeyword("contact") << contact_<< token::END_STATEMENT << nl;
1179     os.writeKeyword("penaltyScale")
1180         << penaltyScale_ << token::END_STATEMENT << nl;
1181     os.writeKeyword("nonLinear")
1182         << nonLinearGeometry::nonLinearNames_[nonLinear_]
1183         << token::END_STATEMENT << nl;
1184     os.writeKeyword("orthotropic")
1185         << orthotropic_ << token::END_STATEMENT << nl;
1186 }
1187
1188
1189 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1190
1191 makePatchTypeField(fvPatchVectorField, solidCohesiveFvPatchVectorField);
1192
1193 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1194
1195 } // End namespace Foam
1196
1197 // ************************************************************************* //