applications/solvers/solidMechanics/viscoElasticSolidFoam/calculateSigma.H

   1 {
   2     Info << "\nCalculate total stress" << endl;
   3
   4     scalar t = runTime.value();
   5     scalar tNext = t + runTime.deltaT().value();
   6
   7     instantList Times = runTime.times();
   8
   9 //     Info << "Number of times: " << Times.size() << endl;
  10
  11     sigma = dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero);
  12     DSigmaCorr =
  13         dimensionedSymmTensor
  14         (
  15             "zero",
  16             dimForce/dimArea,
  17             symmTensor::zero
  18         );
  19
  20     for (label i=1; i<Times.size(); i++)
  21     {
  22         runTime.setTime(Times[i], i);
  23
  24         if(runTime.timeIndex() != i)
  25         {
  26             FatalErrorIn(args.executable())
  27                 << "The strain increment Depsilon must be stored for "
  28                     << "each calculated time step. "
  29                     << exit(FatalError);
  30         }
  31
  32         IOobject DepsilonHeader
  33         (
  34             "Depsilon",
  35             runTime.timeName(),
  36             mesh,
  37             IOobject::MUST_READ
  38         );
  39
  40         // Check Depsilon exists
  41         if (DepsilonHeader.headerOk())
  42         {
  43             volSymmTensorField DepsilonOld(DepsilonHeader, mesh);
  44
  45             scalar tau = runTime.value() - m*runTime.deltaT().value();
  46
  47             if(tau < 0)
  48             {
  49                 sigma += 2.0*rheology.mu(t)*DepsilonOld
  50                     + rheology.lambda(t)*(I*tr(DepsilonOld));
  51
  52                 DSigmaCorr += 2.0*rheology.mu(tNext)*DepsilonOld
  53                     + rheology.lambda(tNext)*(I*tr(DepsilonOld));
  54             }
  55             else
  56             {
  57                 sigma += 2.0*rheology.mu(t - tau)*DepsilonOld
  58                     + rheology.lambda(t - tau)*(I*tr(DepsilonOld));
  59
  60                 DSigmaCorr += 2.0*rheology.mu(tNext - tau)*DepsilonOld
  61                     + rheology.lambda(tNext - tau)*(I*tr(DepsilonOld));
  62             }
  63         }
  64         else
  65         {
  66             Info << "No Depsilon" << endl;
  67         }
  68     }
  69
  70     if(Times.size()>=2)
  71     {
  72         runTime++;
  73     }
  74
  75     scalar tau = runTime.value() - m*runTime.deltaT().value();
  76
  77     sigma += 2.0*rheology.mu(t - tau)*Depsilon
  78         + rheology.lambda(t - tau)*(I*tr(Depsilon));
  79
  80     DSigmaCorr += 2.0*rheology.mu(tNext - tau)*Depsilon
  81         + rheology.lambda(tNext - tau)*(I*tr(Depsilon));
  82
  83     DSigmaCorr -= sigma;
  84
  85     //Info << "Current time = " << runTime.timeName() << endl;
  86 }
  87