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   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
 415         // we need to modify propertiess after cracking otherwise momentum equation is wrong
 416         // but solidInterface seems to hold some information about old mesh
 417         // so we will delete it and make another
 418         // we could probably add a public clearout function
 419         // create new solidInterface
 420         //Pout << "Creating new solidInterface" << endl;
 421         //delete solidInterfacePtr;
 422         //solidInterfacePtr = new solidInterface(mesh, rheology);
 423         // delete demand driven data as the mesh has changed
 424         if(rheology.solidInterfaceActive())
 425           {
 426             rheology.solInterface().clearOut();
 427             solidInterfacePtr->modifyProperties(muf, lambdaf);
 428           }
 429
 430         // Local crack displacement
 431         vectorField UpI =
 432             U.boundaryField()[cohesivePatchID].patchInternalField();
 433         vectorField oldUpI =
 434             U.oldTime().boundaryField()[cohesivePatchID].patchInternalField();
 435
 436         // Global crack displacement
 437         vectorField globalUpI = mesh.globalCrackField(UpI);
 438         vectorField globalOldUpI = mesh.globalCrackField(oldUpI);
 439
 440         // mu and lambda field on new crack faces must be updated
 441         scalarField muPI = mu.boundaryField()[cohesivePatchID].patchInternalField();
 442         scalarField lambdaPI = lambda.boundaryField()[cohesivePatchID].patchInternalField();
 443         scalarField globalMuPI = mesh.globalCrackField(muPI);
 444         scalarField globalLambdaPI = mesh.globalCrackField(lambdaPI);
 445
 446         // cohesivePatchU.size()
 447         int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
 448
 449         // Initialise U for new cohesive face
 450         const labelList& gcfa = mesh.globalCrackFaceAddressing();
 451         label globalIndex = mesh.localCrackStart();
 452         // for (label i=0; i<cohesivePatchU.size(); i++)
 453         for (label i=0; i<cohesivePatchSize; i++)
 454           {
 455             label oldFaceIndex = faceMap[start+i];
 456
 457             // If new face
 458             if (oldFaceIndex == faceToBreakIndex)
 459             {
 460                 U.boundaryField()[cohesivePatchID][i] =
 461                     0.5
 462                    *(
 463                        globalUpI[globalIndex]
 464                      + globalUpI[gcfa[globalIndex]]
 465                     );
 466                 U.oldTime().boundaryField()[cohesivePatchID][i] =
 467                     0.5
 468                    *(
 469                        globalOldUpI[globalIndex]
 470                      + globalOldUpI[gcfa[globalIndex]]
 471                     );
 472
 473                 // initialise mu and lambda on new faces
 474                 // set new face value to value of internal cell
 475                 muf.boundaryField()[cohesivePatchID][i] = globalMuPI[globalIndex];
 476                 lambdaf.boundaryField()[cohesivePatchID][i] = globalLambdaPI[globalIndex];
 477
 478                 globalIndex++;
 479             }
 480             else
 481             {
 482                 globalIndex++;
 483             }
 484         }
 485
 486         // we must calculate grad using interface
 487         // U at the interface has not been calculated yet as interface.correct()
 488         // has not been called yet
 489         // not really a problem as gradU is correct in second outer iteration
 490         // as long as this does not cause convergence problems for the first iterations.
 491         // we should be able to calculate the interface displacements without
 492         // having to call interface.correct()
 493         // todo: add calculateInterfaceU() function
 494         // interface grad uses Gauss, we need least squares
 495         //gradU = solidInterfacePtr->grad(U);
 496         gradU = fvc::grad(U); // leastSquaresSolidInterface grad scheme
 497         //snGradU = fvc::snGrad(U);
 498
 499 #       include "calculateTraction.H"
 500         //if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
 501
 502         // Initialise initiation traction for new cohesive patch face
 503         // for (label i=0; i<cohesivePatchU.size(); i++)
 504         for (label i=0; i<cohesivePatchSize; i++)
 505         {
 506             label oldFaceIndex = faceMap[start+i];
 507
 508             // If new face
 509             if
 510             (
 511                 (oldFaceIndex == faceToBreakIndex)
 512              || (oldFaceIndex == coupledFaceToBreakIndex)
 513             )
 514             {
 515                 vector n0 =
 516                     mesh.Sf().boundaryField()[cohesivePatchID][i]
 517                    /mesh.magSf().boundaryField()[cohesivePatchID][i];
 518                 //vector n1 = -n0;
 519
 520                 if ((n0&faceToBreakNormal) > SMALL)
 521                 {
 522                   traction.boundaryField()[cohesivePatchID][i] =
 523                     faceToBreakTraction;
 524
 525                   traction.oldTime().boundaryField()[cohesivePatchID][i] =
 526                     faceToBreakTraction;
 527
 528                   if(cohesivePatchUPtr)
 529                     {
 530                       cohesivePatchUPtr->traction()[i] = faceToBreakTraction;
 531                     }
 532                   else
 533                     {
 534                       cohesivePatchUFixedModePtr->traction()[i] = faceToBreakTraction;
 535                       cohesivePatchUFixedModePtr->initiationTraction()[i] = faceToBreakTraction;
 536                     }
 537                 }
 538                 else
 539                 {
 540                   traction.boundaryField()[cohesivePatchID][i] =
 541                     -faceToBreakTraction;
 542
 543                   traction.oldTime().boundaryField()[cohesivePatchID][i] =
 544                     -faceToBreakTraction;
 545
 546                   //cohesivePatchU.traction()[i] = -faceToBreakTraction;
 547                   if(cohesivePatchUPtr)
 548                     {
 549                       cohesivePatchUPtr->traction()[i] = -faceToBreakTraction;
 550                     }
 551                   else
 552                     {
 553                       cohesivePatchUFixedModePtr->traction()[i] = -faceToBreakTraction;
 554                       cohesivePatchUFixedModePtr->initiationTraction()[i] = -faceToBreakTraction;
 555                     }
 556                 }
 557             }
 558         }
 559
 560         // hmmnn we only need a reference for very small groups of cells
 561         // turn off for now
 562         //#       include "updateReference.H"
 563     }
 564 }