applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/elasticIncrAcpSolidFoam.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 Application
  25     elasticIncrAcpSolidFoam
  26
  27 Description
  28     Incremental form of elasticAcpSolidFoam
  29     arbitrary crack propagation solver
  30
  31 Author
  32     Zeljko Tukovic, FSB Zagreb
  33     Declan Carolan UCD
  34     Philip Cardiff UCD
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "fvCFD.H"
  39 #include "constitutiveModel.H"
  40 //#include "componentReferenceList.H"
  41 #include "crackerFvMesh.H"
  42 #include "processorPolyPatch.H"
  43 #include "SortableList.H"
  44 #include "solidInterface.H"
  45 #include "solidCohesiveFvPatchVectorField.H"
  46 #include "solidCohesiveFixedModeMixFvPatchVectorField.H"
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 int main(int argc, char *argv[])
  51 {
  52 #   include "setRootCase.H"
  53 #   include "createTime.H"
  54 #   include "createCrackerMesh.H"
  55 #   include "createFields.H"
  56 #   include "createCrack.H"
  57 //#   include "createReference.H"
  58 #   include "createHistory.H"
  59 #   include "readDivDSigmaExpMethod.H"
  60 #   include "createSolidInterfaceIncrNoModify.H"
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64     Info<< "\nStarting time loop\n" << endl;
  65
  66     lduMatrix::debug = 0;
  67
  68     scalar maxEffTractionFraction = 0;
  69
  70     while (runTime.run())
  71     {
  72 #       include "readSolidMechanicsControls.H"
  73 #       include "setDeltaT.H"
  74
  75         runTime++;
  76
  77         Info<< "\nTime = " << runTime.timeName() << " s\n" << endl;
  78
  79         volScalarField rho = rheology.rho();
  80         volScalarField mu = rheology.mu();
  81         volScalarField lambda = rheology.lambda();
  82         surfaceScalarField muf = fvc::interpolate(mu);
  83         surfaceScalarField lambdaf = fvc::interpolate(lambda);
  84
  85         solidInterfacePtr->modifyProperties(muf, lambdaf);
  86         //#     include "waveCourantNo.H"
  87
  88         int iCorr = 0;
  89         lduMatrix::solverPerformance solverPerf;
  90         scalar initialResidual = 0;
  91         scalar relativeResidual = 1;
  92         //scalar forceResidual = 1;
  93         label nFacesToBreak = 0;
  94         label nCoupledFacesToBreak = 0;
  95         bool topoChange = false;
  96
  97         // DU from the previous timestep is usually a good guess
  98         // for the next timestep, but it can cause faces to prematurely
  99         // crack.
 100         // so I will reduce DU here to stop this happening
 101         if (!predictor)
 102         {
 103             DU *= 0.0;
 104         }
 105
 106         do
 107         {
 108             surfaceVectorField n = mesh.Sf()/mesh.magSf();
 109             do
 110             {
 111                 DU.storePrevIter();
 112
 113 #               include "calculateDivDSigmaExp.H"
 114
 115                 fvVectorMatrix DUEqn
 116                 (
 117                     rho*fvm::d2dt2(DU)
 118                     ==
 119                     fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
 120                   + divDSigmaExp
 121                 );
 122
 123                 //#             include "setReference.H"
 124
 125                 if(solidInterfacePtr)
 126                 {
 127                     solidInterfacePtr->correct(DUEqn);
 128                 }
 129
 130                 //DUEqn.relax();
 131
 132                 solverPerf = DUEqn.solve();
 133
 134                 if (aitkenRelax)
 135                 {
 136 #                   include "aitkenRelaxation.H"
 137                 }
 138                 else
 139                 {
 140                     DU.relax();
 141                 }
 142
 143                 if (iCorr == 0)
 144                 {
 145                     initialResidual = solverPerf.initialResidual();
 146                     aitkenInitialRes = gMax(mag(DU.internalField()));
 147                 }
 148
 149                 //gradDU = solidInterfacePtr->grad(DU);
 150                 // use leastSquaresSolidInterface grad scheme
 151                 gradDU = fvc::grad(DU);
 152
 153 #               include "calculateRelativeResidual.H"
 154
 155                 if (iCorr % infoFrequency == 0)
 156                 {
 157                     Info<< "\tTime " << runTime.value()
 158                         << ", Corr " << iCorr
 159                         << ", Solving for " << DU.name()
 160                         << " using " << solverPerf.solverName()
 161                         << ", res = " << solverPerf.initialResidual()
 162                         << ", rel res = " << relativeResidual;
 163                     if (aitkenRelax)
 164                     {
 165                         Info << ", aitken = " << aitkenTheta;
 166                     }
 167                     Info << ", inner iters " << solverPerf.nIterations() << endl;
 168                 }
 169             }
 170             while
 171             (
 172                 //iCorr++ == 0
 173                 iCorr++ < 10
 174                 ||
 175                 (
 176                     //solverPerf.initialResidual() > convergenceTolerance
 177                     relativeResidual > convergenceTolerance
 178                     &&
 179                     iCorr < nCorr
 180                 )
 181             );
 182
 183             Info<< "Solving for " << DU.name() << " using "
 184                 << solverPerf.solverName()
 185                 << ", Initial residual = " << initialResidual
 186                 << ", Final residual = " << solverPerf.initialResidual()
 187                 << ", No outer iterations " << iCorr
 188                 << ", Relative residual " << relativeResidual << endl;
 189
 190 #           include "calculateTraction.H"
 191 #           include "updateCrack.H"
 192
 193             Info<< "Max effective traction fraction: "
 194                 << maxEffTractionFraction << endl;
 195
 196             // reset counter if faces want to crack
 197             if ((nFacesToBreak > 0)  || (nCoupledFacesToBreak > 0)) iCorr = 0;
 198         }
 199         while( (nFacesToBreak > 0)  || (nCoupledFacesToBreak > 0));
 200
 201         if (cohesivePatchDUPtr)
 202         {
 203             if (returnReduce(cohesivePatchDUPtr->size(), sumOp<label>()))
 204             {
 205                 cohesivePatchDUPtr->cracking();
 206             }
 207         }
 208         else
 209         {
 210             if
 211             (
 212                 returnReduce
 213                 (
 214                     cohesivePatchDUFixedModePtr->size(),
 215                     sumOp<label>()
 216                 )
 217             )
 218             {
 219                 Pout << "Number of faces in crack: "
 220                     << cohesivePatchDUFixedModePtr->size() << endl;
 221                 cohesivePatchDUFixedModePtr->relativeSeparationDistance();
 222             }
 223         }
 224
 225 #       include "calculateDEpsilonDSigma.H"
 226
 227         // update total quantities
 228         U += DU;
 229         epsilon += DEpsilon;
 230         sigma += DSigma;
 231
 232 #       include "writeFields.H"
 233 #       include "writeHistory.H"
 234
 235         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 236             << "  ClockTime = " << runTime.elapsedClockTime() << " s\n\n"
 237             << endl;
 238     }
 239
 240     Info<< "End\n" << endl;
 241
 242     return(0);
 243 }
 244
 245
 246 // ************************************************************************* //