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