applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H

   1 nFacesToBreak = 0;
   2 nCoupledFacesToBreak = 0;
   3 {
   4     // Check internal faces
   5
   6     // scalarField effTraction =
   7     //   cohesiveZone.internalField() *
   8     //   mag(traction.internalField());
   9     scalarField normalTraction =
  10         cohesiveZone.internalField()*
  11         ( n.internalField() & traction.internalField() );
  12
  13     // only consider tensile tractions
  14     normalTraction = max(normalTraction, scalar(0));
  15     scalarField shearTraction =
  16         cohesiveZone.internalField()*
  17         mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
  18
  19     // the traction fraction is monitored to decide which faces to break:
  20     // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
  21
  22     const surfaceScalarField sigmaMax = rheology.cohLaw().sigmaMax();
  23     const surfaceScalarField tauMax = rheology.cohLaw().tauMax();
  24
  25     const scalarField& sigmaMaxI = sigmaMax.internalField();
  26     const scalarField& tauMaxI = tauMax.internalField();
  27
  28     //scalarField effTractionFraction = effTraction/sigmaMax;
  29     scalarField effTractionFraction(normalTraction.size(), 0.0);
  30
  31     if(cohesivePatchUPtr)
  32     {
  33         effTractionFraction =
  34             (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
  35           + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
  36     }
  37     else
  38     {
  39         // solidCohesiveFixedModeMix only uses sigmaMax
  40         effTractionFraction =
  41             (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
  42           + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
  43     }
  44
  45     maxEffTractionFraction = gMax(effTractionFraction);
  46
  47     SLList<label> facesToBreakList;
  48     SLList<scalar> facesToBreakEffTractionFractionList;
  49
  50     forAll(effTractionFraction, faceI)
  51     {
  52         if (effTractionFraction[faceI] > 1.0)
  53         {
  54             facesToBreakList.insert(faceI);
  55             facesToBreakEffTractionFractionList.insert
  56             (
  57                 effTractionFraction[faceI]
  58             );
  59         }
  60     }
  61
  62     labelList facesToBreak(facesToBreakList);
  63     List<scalar> facesToBreakEffTractionFraction
  64     (
  65         facesToBreakEffTractionFractionList
  66     );
  67
  68     nFacesToBreak = facesToBreak.size();
  69
  70     // Break only one face per topo change
  71     if (nFacesToBreak > 1)
  72     {
  73         nFacesToBreak = 1;
  74     }
  75
  76     // philipc - select face with maximum effective traction fraction
  77     label faceToBreakIndex = -1;
  78     scalar faceToBreakEffTractionFraction = 0;
  79     forAll(facesToBreakEffTractionFraction, faceI)
  80     {
  81         if
  82         (
  83             facesToBreakEffTractionFraction[faceI]
  84           > faceToBreakEffTractionFraction
  85         )
  86         {
  87             faceToBreakEffTractionFraction =
  88                 facesToBreakEffTractionFraction[faceI];
  89             faceToBreakIndex = facesToBreak[faceI];
  90         }
  91     }
  92
  93     scalar gMaxEffTractionFraction =
  94         returnReduce(faceToBreakEffTractionFraction, maxOp<scalar>());
  95
  96     if (Pstream::parRun())
  97     {
  98         bool procHasFaceToBreak = false;
  99         if (nFacesToBreak > 0)
 100         {
 101             if
 102             (
 103                 mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
 104               < SMALL
 105             )
 106             {
 107                 // philipc - Maximum traction fraction is on this processor
 108                 procHasFaceToBreak = true;
 109             }
 110         }
 111
 112         // Check if maximum is present on more then one processors
 113         label procID = Pstream::nProcs();
 114         if (procHasFaceToBreak)
 115         {
 116             procID = Pstream::myProcNo();
 117         }
 118
 119         label minProcID =
 120             returnReduce<label>(procID, minOp<label>());
 121
 122         if (procID != minProcID)
 123         {
 124             nFacesToBreak = 0;
 125         }
 126     }
 127
 128
 129     // Check coupled (processor) patches
 130
 131     SLList<label> coupledFacesToBreakList;
 132     SLList<scalar> coupledFacesToBreakEffTractionFractionList;
 133     forAll(mesh.boundary(), patchI)
 134     {
 135         if (mesh.boundary()[patchI].coupled())
 136         {
 137             // scalarField pEffTraction =
 138             //   cohesiveZone.boundaryField()[patchI]*
 139             //   mag(traction.boundaryField()[patchI]);
 140             // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
 141
 142             scalarField pNormalTraction =
 143                 cohesiveZone.boundaryField()[patchI]*
 144                 ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
 145
 146             // only consider tensile tractions
 147             pNormalTraction = max(pNormalTraction, scalar(0));
 148
 149             scalarField pShearTraction =
 150                 cohesiveZone.boundaryField()[patchI]*
 151                 mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
 152
 153             // the traction fraction is monitored to decide which faces to break:
 154             // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
 155             const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
 156             const scalarField& pTauMax = tauMax.boundaryField()[patchI];
 157
 158             scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
 159             if(cohesivePatchUPtr)
 160             {
 161                 pEffTractionFraction =
 162                     (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
 163                   + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
 164             }
 165             else
 166             {
 167                 // solidCohesiveFixedModeMix only uses sigmaMax
 168                 pEffTractionFraction =
 169                     (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
 170                   + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
 171             }
 172
 173             label start = mesh.boundaryMesh()[patchI].start();
 174
 175             forAll(pEffTractionFraction, faceI)
 176             {
 177                 if (pEffTractionFraction[faceI] > maxEffTractionFraction)
 178                 {
 179                     maxEffTractionFraction = pEffTractionFraction[faceI];
 180                 }
 181
 182                 if (pEffTractionFraction[faceI] > 1.0)
 183                 {
 184                     coupledFacesToBreakList.insert(start + faceI);
 185                     coupledFacesToBreakEffTractionFractionList.insert
 186                     (
 187                         pEffTractionFraction[faceI]
 188                     );
 189                 }
 190             }
 191         }
 192     }
 193
 194     labelList coupledFacesToBreak(coupledFacesToBreakList);
 195     List<scalar> coupledFacesToBreakEffTractionFraction
 196     (
 197         coupledFacesToBreakEffTractionFractionList
 198     );
 199
 200     nCoupledFacesToBreak = coupledFacesToBreak.size();
 201
 202     // Break only one face per topo change
 203     if (nCoupledFacesToBreak > 1)
 204     {
 205         nCoupledFacesToBreak = 1;
 206     }
 207
 208     // Select coupled face with maximum effective traction fraction
 209     label coupledFaceToBreakIndex = -1;
 210     scalar coupledFaceToBreakEffTractionFraction = 0;
 211     forAll(coupledFacesToBreakEffTractionFraction, faceI)
 212     {
 213         if
 214         (
 215             coupledFacesToBreakEffTractionFraction[faceI]
 216           > coupledFaceToBreakEffTractionFraction
 217         )
 218         {
 219             coupledFaceToBreakEffTractionFraction =
 220                 coupledFacesToBreakEffTractionFraction[faceI];
 221             coupledFaceToBreakIndex = coupledFacesToBreak[faceI];
 222         }
 223     }
 224
 225     scalar gMaxCoupledEffTractionFraction =
 226         returnReduce(coupledFaceToBreakEffTractionFraction, maxOp<scalar>());
 227
 228     if (Pstream::parRun())
 229     {
 230
 231         bool procHasCoupledFaceToBreak = false;
 232         if (nCoupledFacesToBreak > 0)
 233         {
 234             if
 235             (
 236                 mag(gMaxCoupledEffTractionFraction - coupledFaceToBreakEffTractionFraction)
 237               < SMALL
 238             )
 239             {
 240                 // Maximum traction fraction is on this processor
 241                 procHasCoupledFaceToBreak = true;
 242             }
 243         }
 244
 245         // Check if maximum is present on more then one processors
 246         label procID = Pstream::nProcs();
 247         if (procHasCoupledFaceToBreak)
 248         {
 249             procID = Pstream::myProcNo();
 250         }
 251
 252         label minProcID =
 253             returnReduce<label>(procID, minOp<label>());
 254
 255         if (procID != minProcID)
 256         {
 257             nCoupledFacesToBreak = 0;
 258         }
 259     }
 260
 261     if (gMaxCoupledEffTractionFraction > gMaxEffTractionFraction)
 262     {
 263         // Break coupled face
 264         nFacesToBreak = 0;
 265     }
 266     else
 267     {
 268         // Break internal face
 269         nCoupledFacesToBreak = 0;
 270     }
 271
 272     // Make sure that coupled faces are broken in pairs
 273     labelList ngbProc(Pstream::nProcs(), -1);
 274     labelList index(Pstream::nProcs(), -1);
 275     if (nCoupledFacesToBreak)
 276     {
 277         label patchID =
 278             mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
 279
 280         label start = mesh.boundaryMesh()[patchID].start();
 281         label localIndex = coupledFaceToBreakIndex - start;
 282
 283         const processorPolyPatch& procPatch =
 284             refCast<const processorPolyPatch>(mesh.boundaryMesh()[patchID]);
 285         label ngbProcNo = procPatch.neighbProcNo();
 286
 287         ngbProc[Pstream::myProcNo()] = ngbProcNo;
 288         index[Pstream::myProcNo()] = localIndex;
 289     }
 290
 291     if (returnReduce(nCoupledFacesToBreak, maxOp<label>()))
 292     {
 293         reduce(ngbProc, maxOp<labelList>());
 294         reduce(index, maxOp<labelList>());
 295
 296         for (label procI = 0; procI < Pstream::nProcs(); procI++)
 297         {
 298             if (procI != Pstream::myProcNo())
 299             {
 300                 if (ngbProc[procI] == Pstream::myProcNo())
 301                 {
 302                     forAll(mesh.boundaryMesh(), patchI)
 303                     {
 304                         if
 305                         (
 306                             mesh.boundaryMesh()[patchI].type()
 307                          == processorPolyPatch::typeName
 308                         )
 309                         {
 310                             const processorPolyPatch& procPatch =
 311                                 refCast<const processorPolyPatch>
 312                                 (
 313                                     mesh.boundaryMesh()[patchI]
 314                                 );
 315                             label ngbProcNo = procPatch.neighbProcNo();
 316
 317                             if (ngbProcNo == procI)
 318                             {
 319                                 label start =
 320                                     mesh.boundaryMesh()[patchI].start();
 321                                 coupledFaceToBreakIndex = start + index[procI];
 322                                 nCoupledFacesToBreak = 1;
 323                             }
 324                         }
 325                     }
 326                 }
 327             }
 328         }
 329     }
 330
 331
 332     vector faceToBreakTraction = vector::zero;
 333     vector faceToBreakNormal = vector::zero;
 334     scalar faceToBreakSigmaMax = 0.0;
 335     scalar faceToBreakTauMax = 0.0;
 336
 337     // Set faces to break
 338     if (nFacesToBreak > 0)
 339     {
 340         faceToBreakTraction = traction.internalField()[faceToBreakIndex];
 341         faceToBreakNormal = n.internalField()[faceToBreakIndex];
 342
 343         // Scale broken face traction
 344         faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex];
 345         faceToBreakTauMax = tauMaxI[faceToBreakIndex];
 346         scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
 347         normalTrac = max(normalTrac, 0.0);
 348         scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
 349         scalar scaleFactor = 1;
 350         if(cohesivePatchUPtr)
 351         {
 352             scaleFactor =
 353                 Foam::sqrt
 354                 (
 355                     1 /
 356                     (
 357                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 358                       + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
 359                     )
 360                 );
 361         }
 362         else
 363         {
 364             // solidCohesiveFixedModeMix only uses sigmaMax
 365             scaleFactor =
 366                 Foam::sqrt
 367                 (
 368                     1 /
 369                     (
 370                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 371                       + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
 372                     )
 373                 );
 374         }
 375
 376         faceToBreakTraction *= scaleFactor;
 377
 378         topoChange = true;
 379     }
 380     else if (nCoupledFacesToBreak > 0)
 381     {
 382         label patchID =
 383             mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
 384         label start = mesh.boundaryMesh()[patchID].start();
 385         label localIndex = coupledFaceToBreakIndex - start;
 386
 387         faceToBreakTraction = traction.boundaryField()[patchID][localIndex];
 388         faceToBreakNormal = n.boundaryField()[patchID][localIndex];
 389
 390         // Scale broken face traction
 391         faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
 392         faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
 393         scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
 394         normalTrac = max(normalTrac, 0.0);
 395         scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
 396         scalar scaleFactor = 1;
 397         if(cohesivePatchUPtr)
 398         {
 399             scaleFactor =
 400                 Foam::sqrt
 401                 (
 402                     1 /
 403                     (
 404                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 405                       + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
 406                     )
 407                 );
 408         }
 409         else
 410         {
 411             // solidCohesiveFixedModeMix only uses sigmaMax
 412             scaleFactor =
 413                 Foam::sqrt
 414                 (
 415                     1 /
 416                     (
 417                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 418                       + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
 419                     )
 420                 );
 421         }
 422
 423         faceToBreakTraction *= scaleFactor;
 424
 425         topoChange = true;
 426     }
 427
 428     reduce(topoChange, orOp<bool>());
 429
 430     labelList faceToBreak(nFacesToBreak, faceToBreakIndex);
 431     boolList faceToBreakFlip(nFacesToBreak, false);
 432     labelList coupledFaceToBreak
 433     (
 434         nCoupledFacesToBreak,
 435         coupledFaceToBreakIndex
 436     );
 437
 438     reduce(nFacesToBreak, maxOp<label>());
 439     reduce(nCoupledFacesToBreak, maxOp<label>());
 440
 441     if (nFacesToBreak || nCoupledFacesToBreak)
 442     {
 443         Pout << "Internal face to break: " << faceToBreak << endl;
 444         Pout << "Coupled face to break: " << coupledFaceToBreak << endl;
 445
 446         mesh.setBreak(faceToBreak, faceToBreakFlip, coupledFaceToBreak);
 447         mesh.update();
 448
 449         const labelList& faceMap = mesh.topoChangeMap().faceMap();
 450         label start = mesh.boundaryMesh()[cohesivePatchID].start();
 451
 452         mu = rheology.mu();
 453         lambda = rheology.lambda();
 454         muf = fvc::interpolate(mu);
 455         lambdaf = fvc::interpolate(lambda);
 456
 457         // we need to modify propertiess after cracking otherwise momentum equation is wrong
 458         // but solidInterface seems to hold some information about old mesh
 459         // so we will delete it and make another
 460         // we could probably add a public clearout function
 461         // create new solidInterface
 462         //Pout << "Creating new solidInterface" << endl;
 463         //delete solidInterfacePtr;
 464         //solidInterfacePtr = new solidInterface(mesh, rheology);
 465         // delete demand driven data as the mesh has changed
 466         if(rheology.solidInterfaceActive())
 467         {
 468             rheology.solInterface().clearOut();
 469             solidInterfacePtr->modifyProperties(muf, lambdaf);
 470         }
 471
 472         // Local crack displacement
 473         vectorField UpI =
 474             U.boundaryField()[cohesivePatchID].patchInternalField();
 475         vectorField oldUpI =
 476             U.oldTime().boundaryField()[cohesivePatchID].patchInternalField();
 477
 478         // Global crack displacement
 479         vectorField globalUpI = mesh.globalCrackField(UpI);
 480         vectorField globalOldUpI = mesh.globalCrackField(oldUpI);
 481
 482         // mu and lambda field on new crack faces must be updated
 483         scalarField muPI = mu.boundaryField()[cohesivePatchID].patchInternalField();
 484         scalarField lambdaPI = lambda.boundaryField()[cohesivePatchID].patchInternalField();
 485         scalarField globalMuPI = mesh.globalCrackField(muPI);
 486         scalarField globalLambdaPI = mesh.globalCrackField(lambdaPI);
 487
 488         // cohesivePatchU.size()
 489         int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
 490
 491         // Initialise U for new cohesive face
 492         const labelList& gcfa = mesh.globalCrackFaceAddressing();
 493         label globalIndex = mesh.localCrackStart();
 494         // for (label i=0; i<cohesivePatchU.size(); i++)
 495         for (label i=0; i<cohesivePatchSize; i++)
 496         {
 497             label oldFaceIndex = faceMap[start+i];
 498
 499             // If new face
 500             if (oldFaceIndex == faceToBreakIndex)
 501             {
 502                 U.boundaryField()[cohesivePatchID][i] =
 503                     0.5
 504                    *(
 505                        globalUpI[globalIndex]
 506                      + globalUpI[gcfa[globalIndex]]
 507                     );
 508                 U.oldTime().boundaryField()[cohesivePatchID][i] =
 509                     0.5
 510                    *(
 511                        globalOldUpI[globalIndex]
 512                      + globalOldUpI[gcfa[globalIndex]]
 513                     );
 514
 515                 // initialise mu and lambda on new faces
 516                 // set new face value to value of internal cell
 517                 muf.boundaryField()[cohesivePatchID][i] = globalMuPI[globalIndex];
 518                 lambdaf.boundaryField()[cohesivePatchID][i] = globalLambdaPI[globalIndex];
 519
 520                 globalIndex++;
 521             }
 522             else
 523             {
 524                 globalIndex++;
 525             }
 526         }
 527
 528         // we must calculate grad using interface
 529         // U at the interface has not been calculated yet as interface.correct()
 530         // has not been called yet
 531         // not really a problem as gradU is correct in second outer iteration
 532         // as long as this does not cause convergence problems for the first iterations.
 533         // we should be able to calculate the interface displacements without
 534         // having to call interface.correct()
 535         // todo: add calculateInterfaceU() function
 536         // interface grad uses Gauss, we need least squares
 537         //gradU = solidInterfacePtr->grad(U);
 538         gradU = fvc::grad(U); // leastSquaresSolidInterface grad scheme
 539         //snGradU = fvc::snGrad(U);
 540
 541 #       include "calculateTraction.H"
 542         //if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
 543
 544         // Initialise initiation traction for new cohesive patch face
 545         // for (label i=0; i<cohesivePatchU.size(); i++)
 546         for (label i=0; i<cohesivePatchSize; i++)
 547         {
 548             label oldFaceIndex = faceMap[start+i];
 549
 550             // If new face
 551             if
 552             (
 553                 (oldFaceIndex == faceToBreakIndex)
 554              || (oldFaceIndex == coupledFaceToBreakIndex)
 555             )
 556             {
 557                 vector n0 =
 558                     mesh.Sf().boundaryField()[cohesivePatchID][i]
 559                    /mesh.magSf().boundaryField()[cohesivePatchID][i];
 560                 //vector n1 = -n0;
 561
 562                 if ((n0 & faceToBreakNormal) > SMALL)
 563                 {
 564                     traction.boundaryField()[cohesivePatchID][i] =
 565                         faceToBreakTraction;
 566
 567                     traction.oldTime().boundaryField()[cohesivePatchID][i] =
 568                         faceToBreakTraction;
 569
 570                     if(cohesivePatchUPtr)
 571                     {
 572                         cohesivePatchUPtr->traction()[i] = faceToBreakTraction;
 573                     }
 574                     else
 575                     {
 576                         cohesivePatchUFixedModePtr->traction()[i] = faceToBreakTraction;
 577                         cohesivePatchUFixedModePtr->initiationTraction()[i] = faceToBreakTraction;
 578                     }
 579                 }
 580                 else
 581                 {
 582                     traction.boundaryField()[cohesivePatchID][i] =
 583                         -faceToBreakTraction;
 584                     traction.oldTime().boundaryField()[cohesivePatchID][i] =
 585                         -faceToBreakTraction;
 586
 587                     //cohesivePatchU.traction()[i] = -faceToBreakTraction;
 588                     if(cohesivePatchUPtr)
 589                     {
 590                         cohesivePatchUPtr->traction()[i] = -faceToBreakTraction;
 591                     }
 592                     else
 593                     {
 594                         cohesivePatchUFixedModePtr->traction()[i] =
 595                             -faceToBreakTraction;
 596                         cohesivePatchUFixedModePtr->initiationTraction()[i] =
 597                             -faceToBreakTraction;
 598                     }
 599                 }
 600             }
 601         }
 602
 603         // hmmnn we only need a reference for very small groups of cells
 604         // turn off for now
 605         //#       include "updateReference.H"
 606     }
 607 }