applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/elasticOrthoAcpSolidFoam.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     elasticOrthoAcpSolidFoam
  26
  27 Description
  28     Arbitrary crack propagation (ACP) solver
  29     allowing orthotropic material properties.
  30
  31     Please cite:
  32     Cardiff P, Karac A & Ivankovic A, A Large Strain Finite Volume Method for
  33     Orthotropic Bodies with General Material Orientations, Computer Methods
  34     in Applied Mechanics & Engineering, Sep 2013,
  35     http://dx.doi.org/10.1016/j.cma.2013.09.008.
  36
  37     Carolan D, Tuković Z, Murphy N, Ivankovic A, Arbitrary crack propagation
  38     in multi-phase materials using the finite volume method, Computational
  39     Materials Science, 2013, http://dx.doi.org/10.1016/j.commatsci.2012.11.049.
  40
  41 Author
  42     Philip Cardiff UCD
  43     ACP by Tukovic FSB and Carolan UCD
  44
  45 \*---------------------------------------------------------------------------*/
  46
  47 #include "fvCFD.H"
  48 #include "constitutiveModel.H"
  49 //#include "componentReferenceList.H"
  50 #include "crackerFvMesh.H"
  51 #include "processorPolyPatch.H"
  52 #include "SortableList.H"
  53 #include "solidInterface.H"
  54 #include "solidCohesiveFvPatchVectorField.H"
  55 #include "solidCohesiveFixedModeMixFvPatchVectorField.H"
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59 int main(int argc, char *argv[])
  60 {
  61 #   include "setRootCase.H"
  62 #   include "createTime.H"
  63 #   include "createCrackerMesh.H"
  64 #   include "createFields.H"
  65 #   include "createCrack.H"
  66 //# include "createReference.H"
  67 #   include "createHistory.H"
  68 #   include "readDivSigmaExpMethod.H"
  69 #   include "createSolidInterfaceNoModify.H"
  70
  71 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  72
  73     Info<< "\nStarting time loop\n" << endl;
  74
  75     lduMatrix::debug = 0;
  76
  77     scalar maxEffTractionFraction = 0;
  78
  79     // time rates for predictor
  80     volTensorField gradV = fvc::ddt(gradU);
  81     surfaceVectorField snGradV =
  82         (snGradU - snGradU.oldTime())/runTime.deltaT();
  83
  84     //# include "initialiseSolution.H"
  85
  86     while (runTime.run())
  87     {
  88 #       include "readSolidMechanicsControls.H"
  89 #       include "setDeltaT.H"
  90
  91         runTime++;
  92
  93         Info<< "\nTime = " << runTime.timeName() << " s\n" << endl;
  94
  95         volScalarField rho = rheology.rho();
  96         volDiagTensorField K = rheology.K();
  97         surfaceDiagTensorField Kf = fvc::interpolate(K, "K");
  98         volSymmTensor4thOrderField C = rheology.C();
  99         surfaceSymmTensor4thOrderField Cf = fvc::interpolate(C, "C");
 100
 101         solidInterfacePtr->modifyProperties(Cf, Kf);
 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         //bool noMoreCracks = false;
 115
 116         // Predictor step using time rates
 117         if (predictor)
 118         {
 119             Info<< "Predicting U, gradU and snGradU using velocity"
 120                 << endl;
 121             U += V*runTime.deltaT();
 122             gradU += gradV*runTime.deltaT();
 123             snGradU += snGradV*runTime.deltaT();
 124         }
 125
 126         do
 127         {
 128             surfaceVectorField n = mesh.Sf()/mesh.magSf();
 129             do
 130             {
 131                 U.storePrevIter();
 132
 133 #               include "calculateDivSigmaExp.H"
 134
 135                 fvVectorMatrix UEqn
 136                 (
 137                     rho*fvm::d2dt2(U)
 138                  ==
 139                     fvm::laplacian(Kf, U, "laplacian(K,U)")
 140                   + divSigmaExp
 141                 );
 142
 143 //#               include "setReference.H"
 144
 145                 if(solidInterfacePtr)
 146                 {
 147                     solidInterfacePtr->correct(UEqn);
 148                 }
 149
 150                 if (relaxEqn)
 151                 {
 152                     UEqn.relax();
 153                 }
 154
 155                 solverPerf = UEqn.solve();
 156
 157                 if (aitkenRelax)
 158                 {
 159 #                   include "aitkenRelaxation.H"
 160                 }
 161                 else
 162                 {
 163                     U.relax();
 164                 }
 165
 166                 if (iCorr == 0)
 167                 {
 168                     initialResidual = solverPerf.initialResidual();
 169                     aitkenInitialRes = gMax(mag(U.internalField()));
 170                 }
 171
 172                 //gradU = solidInterfacePtr->grad(U);
 173                 // use leastSquaresSolidInterface grad scheme
 174                 gradU = fvc::grad(U);
 175
 176 #               include "calculateRelativeResidual.H"
 177
 178                 if (iCorr % infoFrequency == 0)
 179                 {
 180                     Info<< "\tTime " << runTime.value()
 181                         << ", Corr " << iCorr
 182                         << ", Solving for " << U.name()
 183                         << " using " << solverPerf.solverName()
 184                         << ", res = " << solverPerf.initialResidual()
 185                         << ", rel res = " << relativeResidual;
 186                     if (aitkenRelax)
 187                     {
 188                         Info << ", aitken = " << aitkenTheta;
 189                     }
 190                     Info << ", inner iters " << solverPerf.nIterations() << endl;
 191                 }
 192             }
 193             while
 194             (
 195                 //iCorr++ == 0
 196                 iCorr++ < 10
 197                 ||
 198                 (
 199                     //solverPerf.initialResidual() > convergenceTolerance
 200                     relativeResidual > convergenceTolerance
 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 // ************************************************************************* //