applications/solvers/solidMechanics/elasticPlasticNonLinTLSolidFoam/calculateThirdOrderDissipativeTerm.H

   1 // Calculate third order dissipative term as outlined by Demirdzic
   2 // This removes unphysical high frequency oscillations that may
   3 // appear in the solution. This term goes to zero on convergence
   4 // assuming a linear/quadratic displacement field, or the error
   5 // is second order
   6
   7 // volVectorField divThirdOrderTerm
   8 // (
   9 //  IOobject
  10 //  (
  11 //   "divThirdOrderTerm",
  12 //   runTime.timeName(),
  13 //   mesh,
  14 //   IOobject::NO_READ,
  15 //   IOobject::NO_WRITE
  16 //   ),
  17 //  mesh,
  18 //  dimensionedVector("zero", dimForce/dimVolume, vector::zero)
  19 //  );
  20
  21 // average gradDU of neighbouring cell centres
  22 // interpolation scheme should be midPoint
  23 surfaceTensorField averageGradDU("averageGradDU", fvc::interpolate(gradDU, "averageGradDU"));
  24
  25 // average face gradDU extrapolated from neighbouring cell centres
  26 surfaceTensorField extrapGradDU
  27 (
  28  IOobject
  29  (
  30   "extrapGradDU",
  31   runTime.timeName(),
  32   mesh,
  33   IOobject::NO_READ,
  34   IOobject::NO_WRITE
  35   ),
  36  mesh,
  37  dimensionedTensor("zero", dimless, tensor::zero)
  38  );
  39
  40 volVectorField gradGradDUcompXX = fvc::grad(gradDU.component(tensor::XX), "gradGradDU");
  41 volVectorField gradGradDUcompXY = fvc::grad(gradDU.component(tensor::XY), "gradGradDU");
  42 volVectorField gradGradDUcompXZ = fvc::grad(gradDU.component(tensor::XZ), "gradGradDU");
  43 volVectorField gradGradDUcompYX = fvc::grad(gradDU.component(tensor::YX), "gradGradDU");
  44 volVectorField gradGradDUcompYY = fvc::grad(gradDU.component(tensor::YY), "gradGradDU");
  45 volVectorField gradGradDUcompYZ = fvc::grad(gradDU.component(tensor::YZ), "gradGradDU");
  46 volVectorField gradGradDUcompZX = fvc::grad(gradDU.component(tensor::ZX), "gradGradDU");
  47 volVectorField gradGradDUcompZY = fvc::grad(gradDU.component(tensor::ZY), "gradGradDU");
  48 volVectorField gradGradDUcompZZ = fvc::grad(gradDU.component(tensor::ZZ), "gradGradDU");
  49
  50 // third order tensor components
  51 volScalarField gradGradDUXXX = gradGradDUcompXX.component(vector::X);
  52 volScalarField gradGradDUXXY = gradGradDUcompXY.component(vector::X);
  53 volScalarField gradGradDUXXZ = gradGradDUcompXZ.component(vector::X);
  54
  55 volScalarField gradGradDUXYX = gradGradDUcompYX.component(vector::X);
  56 volScalarField gradGradDUXYY = gradGradDUcompYY.component(vector::X);
  57 volScalarField gradGradDUXYZ = gradGradDUcompYZ.component(vector::X);
  58
  59 volScalarField gradGradDUXZX = gradGradDUcompZX.component(vector::X);
  60 volScalarField gradGradDUXZY = gradGradDUcompZY.component(vector::X);
  61 volScalarField gradGradDUXZZ = gradGradDUcompZZ.component(vector::X);
  62
  63 volScalarField gradGradDUYXX = gradGradDUcompXX.component(vector::Y);
  64 volScalarField gradGradDUYXY = gradGradDUcompXY.component(vector::Y);
  65 volScalarField gradGradDUYXZ = gradGradDUcompXZ.component(vector::Y);
  66
  67 volScalarField gradGradDUYYX = gradGradDUcompYX.component(vector::Y);
  68 volScalarField gradGradDUYYY = gradGradDUcompYY.component(vector::Y);
  69 volScalarField gradGradDUYYZ = gradGradDUcompYZ.component(vector::Y);
  70
  71 volScalarField gradGradDUYZX = gradGradDUcompZX.component(vector::Y);
  72 volScalarField gradGradDUYZY = gradGradDUcompZY.component(vector::Y);
  73 volScalarField gradGradDUYZZ = gradGradDUcompZZ.component(vector::Y);
  74
  75 volScalarField gradGradDUZXX = gradGradDUcompXX.component(vector::Z);
  76 volScalarField gradGradDUZXY = gradGradDUcompXY.component(vector::Z);
  77 volScalarField gradGradDUZXZ = gradGradDUcompXZ.component(vector::Z);
  78
  79 volScalarField gradGradDUZYX = gradGradDUcompYX.component(vector::Z);
  80 volScalarField gradGradDUZYY = gradGradDUcompYY.component(vector::Z);
  81 volScalarField gradGradDUZYZ = gradGradDUcompYZ.component(vector::Z);
  82
  83 volScalarField gradGradDUZZX = gradGradDUcompZX.component(vector::Z);
  84 volScalarField gradGradDUZZY = gradGradDUcompZY.component(vector::Z);
  85 volScalarField gradGradDUZZZ = gradGradDUcompZZ.component(vector::Z);
  86
  87 forAll(extrapGradDU.internalField(), facei)
  88 {
  89   const label own = mesh.owner()[facei];
  90   const label nei = mesh.neighbour()[facei];
  91   const vector deltaOwn = mesh.Cf()[facei] - mesh.C()[own];
  92   const vector deltaNei = mesh.Cf()[facei] - mesh.C()[nei];
  93   const tensor& gradDUOwn = gradDU.internalField()[own];
  94   const tensor& gradDUNei = gradDU.internalField()[nei];
  95
  96   // as there is there is no thirdOrderTensor class, we will
  97   // calculate (deltaOwn&gradGradDUOwn) out manually
  98   // tensor deltaOwnDotgradGradDUOwn = tensor::zero;
  99   // tensor deltaNeiDotgradGradDUNei = tensor::zero;
 100
 101   // deltaOwnDotgradGradDUOwn[tensor::XX] = deltaOwn & gradGradDUcompXX.internalField()[own];
 102   // deltaNeiDotgradGradDUNei[tensor::XX] = deltaNei & gradGradDUcompXX.internalField()[nei];
 103   // deltaOwnDotgradGradDUOwn[tensor::XY] = deltaOwn & gradGradDUcompXY.internalField()[own];
 104   // deltaNeiDotgradGradDUNei[tensor::XY] = deltaNei & gradGradDUcompXY.internalField()[nei];
 105   // deltaOwnDotgradGradDUOwn[tensor::XZ] = deltaOwn & gradGradDUcompXZ.internalField()[own];
 106   // deltaNeiDotgradGradDUNei[tensor::XZ] = deltaNei & gradGradDUcompXZ.internalField()[nei];
 107
 108   // deltaOwnDotgradGradDUOwn[tensor::YX] = deltaOwn & gradGradDUcompYX.internalField()[own];
 109   // deltaNeiDotgradGradDUNei[tensor::YX] = deltaNei & gradGradDUcompYX.internalField()[nei];
 110   // deltaOwnDotgradGradDUOwn[tensor::YY] = deltaOwn & gradGradDUcompYY.internalField()[own];
 111   // deltaNeiDotgradGradDUNei[tensor::YY] = deltaNei & gradGradDUcompYY.internalField()[nei];
 112   // deltaOwnDotgradGradDUOwn[tensor::YZ] = deltaOwn & gradGradDUcompYZ.internalField()[own];
 113   // deltaNeiDotgradGradDUNei[tensor::YZ] = deltaNei & gradGradDUcompYZ.internalField()[nei];
 114
 115   // deltaOwnDotgradGradDUOwn[tensor::ZX] = deltaOwn & gradGradDUcompZX.internalField()[own];
 116   // deltaNeiDotgradGradDUNei[tensor::ZX] = deltaNei & gradGradDUcompZX.internalField()[nei];
 117   // deltaOwnDotgradGradDUOwn[tensor::ZY] = deltaOwn & gradGradDUcompZY.internalField()[own];
 118   // deltaNeiDotgradGradDUNei[tensor::ZY] = deltaNei & gradGradDUcompZY.internalField()[nei];
 119   // deltaOwnDotgradGradDUOwn[tensor::ZZ] = deltaOwn & gradGradDUcompZZ.internalField()[own];
 120   // deltaNeiDotgradGradDUNei[tensor::ZZ] = deltaNei & gradGradDUcompZZ.internalField()[nei];
 121
 122   scalar gradGradDUXXXOwn = gradGradDUXXX.internalField()[own];
 123   scalar gradGradDUXXYOwn = gradGradDUXXY.internalField()[own];
 124   scalar gradGradDUXXZOwn = gradGradDUXXZ.internalField()[own];
 125
 126   scalar gradGradDUXYXOwn = gradGradDUXYX.internalField()[own];
 127   scalar gradGradDUXYYOwn = gradGradDUXYY.internalField()[own];
 128   scalar gradGradDUXYZOwn = gradGradDUXYZ.internalField()[own];
 129
 130   scalar gradGradDUXZXOwn = gradGradDUXZX.internalField()[own];
 131   scalar gradGradDUXZYOwn = gradGradDUXZY.internalField()[own];
 132   scalar gradGradDUXZZOwn = gradGradDUXZZ.internalField()[own];
 133
 134
 135   scalar gradGradDUYXXOwn = gradGradDUYXX.internalField()[own];
 136   scalar gradGradDUYXYOwn = gradGradDUYXY.internalField()[own];
 137   scalar gradGradDUYXZOwn = gradGradDUYXZ.internalField()[own];
 138
 139   scalar gradGradDUYYXOwn = gradGradDUYYX.internalField()[own];
 140   scalar gradGradDUYYYOwn = gradGradDUYYY.internalField()[own];
 141   scalar gradGradDUYYZOwn = gradGradDUYYZ.internalField()[own];
 142
 143   scalar gradGradDUYZXOwn = gradGradDUYZX.internalField()[own];
 144   scalar gradGradDUYZYOwn = gradGradDUYZY.internalField()[own];
 145   scalar gradGradDUYZZOwn = gradGradDUYZZ.internalField()[own];
 146
 147
 148   scalar gradGradDUZXXOwn = gradGradDUZXX.internalField()[own];
 149   scalar gradGradDUZXYOwn = gradGradDUZXY.internalField()[own];
 150   scalar gradGradDUZXZOwn = gradGradDUZXZ.internalField()[own];
 151
 152   scalar gradGradDUZYXOwn = gradGradDUZYX.internalField()[own];
 153   scalar gradGradDUZYYOwn = gradGradDUZYY.internalField()[own];
 154   scalar gradGradDUZYZOwn = gradGradDUZYZ.internalField()[own];
 155
 156   scalar gradGradDUZZXOwn = gradGradDUZZX.internalField()[own];
 157   scalar gradGradDUZZYOwn = gradGradDUZZY.internalField()[own];
 158   scalar gradGradDUZZZOwn = gradGradDUZZZ.internalField()[own];
 159
 160
 161   // nei
 162   scalar gradGradDUXXXNei = gradGradDUXXX.internalField()[nei];
 163   scalar gradGradDUXXYNei = gradGradDUXXY.internalField()[nei];
 164   scalar gradGradDUXXZNei = gradGradDUXXZ.internalField()[nei];
 165
 166   scalar gradGradDUXYXNei = gradGradDUXYX.internalField()[nei];
 167   scalar gradGradDUXYYNei = gradGradDUXYY.internalField()[nei];
 168   scalar gradGradDUXYZNei = gradGradDUXYZ.internalField()[nei];
 169
 170   scalar gradGradDUXZXNei = gradGradDUXZX.internalField()[nei];
 171   scalar gradGradDUXZYNei = gradGradDUXZY.internalField()[nei];
 172   scalar gradGradDUXZZNei = gradGradDUXZZ.internalField()[nei];
 173
 174
 175   scalar gradGradDUYXXNei = gradGradDUYXX.internalField()[nei];
 176   scalar gradGradDUYXYNei = gradGradDUYXY.internalField()[nei];
 177   scalar gradGradDUYXZNei = gradGradDUYXZ.internalField()[nei];
 178
 179   scalar gradGradDUYYXNei = gradGradDUYYX.internalField()[nei];
 180   scalar gradGradDUYYYNei = gradGradDUYYY.internalField()[nei];
 181   scalar gradGradDUYYZNei = gradGradDUYYZ.internalField()[nei];
 182
 183   scalar gradGradDUYZXNei = gradGradDUYZX.internalField()[nei];
 184   scalar gradGradDUYZYNei = gradGradDUYZY.internalField()[nei];
 185   scalar gradGradDUYZZNei = gradGradDUYZZ.internalField()[nei];
 186
 187
 188   scalar gradGradDUZXXNei = gradGradDUZXX.internalField()[nei];
 189   scalar gradGradDUZXYNei = gradGradDUZXY.internalField()[nei];
 190   scalar gradGradDUZXZNei = gradGradDUZXZ.internalField()[nei];
 191
 192   scalar gradGradDUZYXNei = gradGradDUZYX.internalField()[nei];
 193   scalar gradGradDUZYYNei = gradGradDUZYY.internalField()[nei];
 194   scalar gradGradDUZYZNei = gradGradDUZYZ.internalField()[nei];
 195
 196   scalar gradGradDUZZXNei = gradGradDUZZX.internalField()[nei];
 197   scalar gradGradDUZZYNei = gradGradDUZZY.internalField()[nei];
 198   scalar gradGradDUZZZNei = gradGradDUZZZ.internalField()[nei];
 199
 200   tensor deltaOwnDotgradGradDUOwn =
 201     tensor(
 202            deltaOwn.x()*gradGradDUXXXOwn + deltaOwn.y()*gradGradDUYXXOwn + deltaOwn.z()*gradGradDUZXXOwn,
 203            deltaOwn.x()*gradGradDUXXYOwn + deltaOwn.y()*gradGradDUYXYOwn + deltaOwn.z()*gradGradDUZXYOwn,
 204            deltaOwn.x()*gradGradDUXXZOwn + deltaOwn.y()*gradGradDUYXZOwn + deltaOwn.z()*gradGradDUZXZOwn,
 205
 206            deltaOwn.x()*gradGradDUXYXOwn + deltaOwn.y()*gradGradDUYYXOwn + deltaOwn.z()*gradGradDUZYXOwn,
 207            deltaOwn.x()*gradGradDUXYYOwn + deltaOwn.y()*gradGradDUYYYOwn + deltaOwn.z()*gradGradDUZYYOwn,
 208            deltaOwn.x()*gradGradDUXYZOwn + deltaOwn.y()*gradGradDUYYZOwn + deltaOwn.z()*gradGradDUZYZOwn,
 209
 210            deltaOwn.x()*gradGradDUXZXOwn + deltaOwn.y()*gradGradDUYZXOwn + deltaOwn.z()*gradGradDUZZXOwn,
 211            deltaOwn.x()*gradGradDUXZYOwn + deltaOwn.y()*gradGradDUYZYOwn + deltaOwn.z()*gradGradDUZZYOwn,
 212            deltaOwn.x()*gradGradDUXZZOwn + deltaOwn.y()*gradGradDUYZZOwn + deltaOwn.z()*gradGradDUZZZOwn
 213            );
 214
 215   tensor deltaNeiDotgradGradDUNei =
 216     tensor(
 217            deltaNei.x()*gradGradDUXXXNei + deltaNei.y()*gradGradDUYXXNei + deltaNei.z()*gradGradDUZXXNei,
 218            deltaNei.x()*gradGradDUXXYNei + deltaNei.y()*gradGradDUYXYNei + deltaNei.z()*gradGradDUZXYNei,
 219            deltaNei.x()*gradGradDUXXZNei + deltaNei.y()*gradGradDUYXZNei + deltaNei.z()*gradGradDUZXZNei,
 220
 221            deltaNei.x()*gradGradDUXYXNei + deltaNei.y()*gradGradDUYYXNei + deltaNei.z()*gradGradDUZYXNei,
 222            deltaNei.x()*gradGradDUXYYNei + deltaNei.y()*gradGradDUYYYNei + deltaNei.z()*gradGradDUZYYNei,
 223            deltaNei.x()*gradGradDUXYZNei + deltaNei.y()*gradGradDUYYZNei + deltaNei.z()*gradGradDUZYZNei,
 224
 225            deltaNei.x()*gradGradDUXZXNei + deltaNei.y()*gradGradDUYZXNei + deltaNei.z()*gradGradDUZZXNei,
 226            deltaNei.x()*gradGradDUXZYNei + deltaNei.y()*gradGradDUYZYNei + deltaNei.z()*gradGradDUZZYNei,
 227            deltaNei.x()*gradGradDUXZZNei + deltaNei.y()*gradGradDUYZZNei + deltaNei.z()*gradGradDUZZZNei
 228            );
 229
 230
 231   // get average of extrapolated values
 232   tensor extrapFaceGrad =
 233     0.5*
 234     (
 235      gradDUOwn + (deltaOwnDotgradGradDUOwn)
 236      +
 237      gradDUNei + (deltaNeiDotgradGradDUNei)
 238      );
 239
 240   extrapGradDU.internalField()[facei] = extrapFaceGrad;
 241
 242 }
 243
 244 // correction is zero on boundary
 245 forAll(extrapGradDU.boundaryField(), patchi)
 246 {
 247   extrapGradDU.boundaryField()[patchi] =
 248     averageGradDU.boundaryField()[patchi];
 249 }
 250
 251 // calculate thirdOrderTerm
 252 volVectorField divThirdOrderTerm (
 253                                   "thirdOrderTerm",
 254                                   fvc::div(
 255                                            (2*muf+lambdaf)*mesh.Sf()
 256                                            & (extrapGradDU - averageGradDU)
 257                                            )
 258                                   );
 259
 260 // if(runTime.outputTime())
 261 //   {
 262 //     divThirdOrderTerm.write();
 263 //     averageGradDU.write();
 264 //     extrapGradDU.write();
 265 //   }