applications/solvers/solidMechanics/elasticAcpSolidFoam/elasticAcpSolidFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM 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 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Application
  26     elasticAcpSolidFoam
  27
  28 Description
  29     Arbitrary crack propagation solver.
  30     Cracks may propagate along any mesh internal face.
  31
  32     Please cite:
  33     Carolan D, Tuković Z, Murphy N, Ivankovic A, Arbitrary crack propagation
  34     in multi-phase materials using the finite volume method, Computational
  35     Materials Science, 2013, http://dx.doi.org/10.1016/j.commatsci.2012.11.049.
  36
  37 Author
  38     Zeljko Tukovic, FSB Zagreb
  39     Declan Carolan UCD
  40     Philip Cardiff UCD
  41
  42 \*---------------------------------------------------------------------------*/
  43
  44 #include "fvCFD.H"
  45 #include "constitutiveModel.H"
  46 //#include "componentReferenceList.H"
  47 #include "crackerFvMesh.H"
  48 #include "processorPolyPatch.H"
  49 #include "SortableList.H"
  50 #include "solidInterface.H"
  51 #include "solidCohesiveFvPatchVectorField.H"
  52 #include "solidCohesiveFixedModeMixFvPatchVectorField.H"
  53
  54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  55
  56 int main(int argc, char *argv[])
  57 {
  58 # include "setRootCase.H"
  59 # include "createTime.H"
  60 # include "createCrackerMesh.H"
  61 # include "createFields.H"
  62 # include "createCrack.H"
  63 //# include "createReference.H"
  64 # include "createHistory.H"
  65 # include "readDivSigmaExpMethod.H"
  66 # include "createSolidInterfaceNoModify.H"
  67
  68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  69
  70   Info<< "\nStarting time loop\n" << endl;
  71
  72   lduMatrix::debug = 0;
  73
  74   scalar maxEffTractionFraction = 0;
  75
  76   // time rates for predictor
  77   volTensorField gradV = fvc::ddt(gradU);
  78   surfaceVectorField snGradV =
  79     (snGradU - snGradU.oldTime())/runTime.deltaT();
  80
  81   //# include "initialiseSolution.H"
  82
  83   while (runTime.run())
  84   {
  85 #     include "readSolidMechanicsControls.H"
  86 #     include "setDeltaT.H"
  87
  88       runTime++;
  89
  90       Info<< "\nTime: " << runTime.timeName() << " s\n" << endl;
  91
  92       volScalarField rho = rheology.rho();
  93       volScalarField mu = rheology.mu();
  94       volScalarField lambda = rheology.lambda();
  95       surfaceScalarField muf = fvc::interpolate(mu);
  96       surfaceScalarField lambdaf = fvc::interpolate(lambda);
  97
  98       if (solidInterfaceCorr)
  99       {
 100           solidInterfacePtr->modifyProperties(muf, lambdaf);
 101       }
 102
 103       //#     include "waveCourantNo.H"
 104
 105       int iCorr = 0;
 106       lduMatrix::solverPerformance solverPerf;
 107       scalar initialResidual = 0;
 108       scalar relativeResidual = 1;
 109       //scalar forceResidual = 1;
 110       label nFacesToBreak = 0;
 111       label nCoupledFacesToBreak = 0;
 112       bool topoChange = false;
 113
 114       // Predictor step using time rates
 115       if (predictor)
 116         {
 117             Info<< "Predicting U, gradU and snGradU using velocity"
 118                 << endl;
 119             U += V*runTime.deltaT();
 120             gradU += gradV*runTime.deltaT();
 121             snGradU += snGradV*runTime.deltaT();
 122         }
 123
 124       do
 125       {
 126           surfaceVectorField n = mesh.Sf()/mesh.magSf();
 127           do
 128           {
 129               U.storePrevIter();
 130
 131 #             include "calculateDivSigmaExp.H"
 132
 133               fvVectorMatrix UEqn
 134                   (
 135                       rho*fvm::d2dt2(U)
 136                       ==
 137                       fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
 138                       + divSigmaExp
 139                       );
 140
 141               //#             include "setReference.H"
 142
 143               if (solidInterfaceCorr)
 144               {
 145                   solidInterfacePtr->correct(UEqn);
 146               }
 147
 148               if (relaxEqn)
 149               {
 150                   UEqn.relax();
 151               }
 152
 153               solverPerf = UEqn.solve();
 154
 155               if (aitkenRelax)
 156               {
 157 #                 include "aitkenRelaxation.H"
 158               }
 159               else
 160               {
 161                   U.relax();
 162               }
 163
 164               if (iCorr == 0)
 165               {
 166                   initialResidual = solverPerf.initialResidual();
 167                   aitkenInitialRes = gMax(mag(U.internalField()));
 168               }
 169
 170               //gradU = solidInterfacePtr->grad(U);
 171               // use leastSquaresSolidInterface grad scheme
 172               gradU = fvc::grad(U);
 173
 174 #             include "calculateRelativeResidual.H"
 175               //#             include "calculateForceResidual.H"
 176
 177               if (iCorr % infoFrequency == 0)
 178               {
 179                   Info << "\tTime " << runTime.value()
 180                        << ", Corr " << iCorr
 181                        << ", Solving for " << U.name()
 182                        << " using " << solverPerf.solverName()
 183                        << ", res = " << solverPerf.initialResidual()
 184                        << ", rel res = " << relativeResidual;
 185                   if (aitkenRelax)
 186                   {
 187                       Info << ", aitken = " << aitkenTheta;
 188                   }
 189                   Info << ", inner iters " << solverPerf.nIterations() << endl;
 190               }
 191           }
 192           while
 193               (
 194                   //iCorr++ == 0
 195                   iCorr++ < 2
 196                   ||
 197                   (
 198                       solverPerf.initialResidual() > convergenceTolerance
 199                       //relativeResidual > convergenceTolerance
 200                       &&
 201                       iCorr < nCorr
 202                       )
 203                   );
 204
 205           Info << "Solving for " << U.name() << " using "
 206                << solverPerf.solverName()
 207                << ", Initial residual = " << initialResidual
 208                << ", Final residual = " << solverPerf.initialResidual()
 209                << ", No outer iterations " << iCorr
 210                << ", Relative residual " << relativeResidual << endl;
 211
 212 #         include "calculateTraction.H"
 213 #         include "updateCrack.H"
 214
 215           Info<< "Max effective traction fraction: "
 216               << maxEffTractionFraction << endl;
 217
 218           // reset counter if faces want to crack
 219           if ((nFacesToBreak > 0)  || (nCoupledFacesToBreak > 0)) iCorr = 0;
 220       }
 221       while( (nFacesToBreak > 0)  || (nCoupledFacesToBreak > 0));
 222
 223       if (cohesivePatchUPtr)
 224       {
 225           if (returnReduce(cohesivePatchUPtr->size(), sumOp<label>()))
 226           {
 227               cohesivePatchUPtr->cracking();
 228           }
 229       }
 230       else
 231       {
 232           if
 233               (
 234                   returnReduce
 235                   (
 236                       cohesivePatchUFixedModePtr->size(),
 237                       sumOp<label>()
 238                       )
 239                   )
 240           {
 241               Pout << "Number of faces in crack: "
 242                << cohesivePatchUFixedModePtr->size() << endl;
 243               cohesivePatchUFixedModePtr->relativeSeparationDistance();
 244           }
 245       }
 246
 247       // update time rates for predictor
 248       if (predictor)
 249       {
 250           V = fvc::ddt(U);
 251           gradV = fvc::ddt(gradU);
 252           snGradV = (snGradU - snGradU.oldTime())/runTime.deltaT();
 253       }
 254
 255 #     include "calculateEpsilonSigma.H"
 256 #     include "writeFields.H"
 257 #     include "writeHistory.H"
 258
 259       Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 260       << "  ClockTime = " << runTime.elapsedClockTime() << " s\n\n"
 261       << endl;
 262     }
 263
 264   Info<< "End\n" << endl;
 265
 266   return(0);
 267 }
 268
 269
 270 // ************************************************************************* //