applications/solvers/solidMechanics/viscoElasticSolidFoam/calcInitialState.H

   1
   2     {
   3         Info<< "Time: " << runTime.timeName() << nl << endl;
   4
   5 #       include "readSolidMechanicsControls.H"
   6
   7         volScalarField mu = rheology.mu(0);
   8         volScalarField lambda = rheology.lambda(0);
   9
  10         Info << "mu = " << average(mu.internalField()) << endl;
  11         Info << "lambda = " << average(lambda.internalField()) << endl;
  12
  13         int iCorr = 0;
  14         lduMatrix::solverPerformance solverPerf;
  15         scalar initialResidual = 0;
  16         scalar err = GREAT;
  17
  18         do
  19         {
  20             DU.storePrevIter();
  21
  22             fvVectorMatrix DUEqn
  23             (
  24                 fvm::laplacian(2*mu+lambda, DU, "laplacian(DDU,DU)")
  25               + fvc::div
  26                 (
  27                     mu*gradDU.T()
  28                   + lambda*(I*tr(gradDU))
  29                   - (mu + lambda)*gradDU,
  30                     "div(sigma)"
  31                 )
  32             );
  33
  34             solverPerf = DUEqn.solve();
  35
  36             DU.relax();
  37
  38             if(iCorr == 0)
  39             {
  40                 initialResidual = solverPerf.initialResidual();
  41             }
  42
  43             gradDU = fvc::grad(DU);
  44         }
  45         while
  46         (
  47             solverPerf.initialResidual() > convergenceTolerance
  48          && ++iCorr < nCorr
  49         );
  50
  51         Info << "Solving for " << DU.name() << " using "
  52             << solverPerf.solverName() << " solver"
  53             << ", Initial residula = " << initialResidual
  54             << ", Final residual = " << solverPerf.initialResidual()
  55             << ", No outer iterations " << iCorr
  56             << ", Relative error: " << err << endl;
  57
  58         U += DU;
  59
  60 #       include "calculateDSigma.H"
  61
  62         sigma += DSigma;
  63
  64         {
  65             DSigmaCorr =
  66                 dimensionedSymmTensor
  67                 (
  68                     "zero",
  69                     dimForce/dimArea,
  70                     symmTensor::zero
  71                 );
  72             scalar t = runTime.value();
  73             scalar tNext = t + runTime.deltaT().value();
  74             DSigmaCorr += 2.0*rheology.mu(tNext)*Depsilon
  75                 + rheology.lambda(tNext)*(I*tr(Depsilon));
  76             DSigmaCorr -= sigma;
  77         }
  78
  79 #       include "calculateStress.H"
  80
  81 #       include "writeHistory.H"
  82
  83         Depsilon.write();
  84
  85         Info<< "ExecutionTime = "
  86             << runTime.elapsedCpuTime()
  87             << " s\n\n" << endl;
  88     }