applications/solvers/solidMechanics/elasticSolidFoam/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 gradU of neighbouring cell centres
  22 // interpolation scheme should be midPoint
  23 surfaceTensorField averageGradU("averageGradU", fvc::interpolate(gradU, "averageGradU"));
  24
  25 // average face gradU extrapolated from neighbouring cell centres
  26 surfaceTensorField extrapGradU
  27 (
  28  IOobject
  29  (
  30   "extrapGradU",
  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 gradGradUcompXX = fvc::grad(gradU.component(tensor::XX), "gradGradU");
  41 volVectorField gradGradUcompXY = fvc::grad(gradU.component(tensor::XY), "gradGradU");
  42 volVectorField gradGradUcompXZ = fvc::grad(gradU.component(tensor::XZ), "gradGradU");
  43 volVectorField gradGradUcompYX = fvc::grad(gradU.component(tensor::YX), "gradGradU");
  44 volVectorField gradGradUcompYY = fvc::grad(gradU.component(tensor::YY), "gradGradU");
  45 volVectorField gradGradUcompYZ = fvc::grad(gradU.component(tensor::YZ), "gradGradU");
  46 volVectorField gradGradUcompZX = fvc::grad(gradU.component(tensor::ZX), "gradGradU");
  47 volVectorField gradGradUcompZY = fvc::grad(gradU.component(tensor::ZY), "gradGradU");
  48 volVectorField gradGradUcompZZ = fvc::grad(gradU.component(tensor::ZZ), "gradGradU");
  49
  50 // third order tensor components
  51 volScalarField gradGradUXXX = gradGradUcompXX.component(vector::X);
  52 volScalarField gradGradUXXY = gradGradUcompXY.component(vector::X);
  53 volScalarField gradGradUXXZ = gradGradUcompXZ.component(vector::X);
  54
  55 volScalarField gradGradUXYX = gradGradUcompYX.component(vector::X);
  56 volScalarField gradGradUXYY = gradGradUcompYY.component(vector::X);
  57 volScalarField gradGradUXYZ = gradGradUcompYZ.component(vector::X);
  58
  59 volScalarField gradGradUXZX = gradGradUcompZX.component(vector::X);
  60 volScalarField gradGradUXZY = gradGradUcompZY.component(vector::X);
  61 volScalarField gradGradUXZZ = gradGradUcompZZ.component(vector::X);
  62
  63 volScalarField gradGradUYXX = gradGradUcompXX.component(vector::Y);
  64 volScalarField gradGradUYXY = gradGradUcompXY.component(vector::Y);
  65 volScalarField gradGradUYXZ = gradGradUcompXZ.component(vector::Y);
  66
  67 volScalarField gradGradUYYX = gradGradUcompYX.component(vector::Y);
  68 volScalarField gradGradUYYY = gradGradUcompYY.component(vector::Y);
  69 volScalarField gradGradUYYZ = gradGradUcompYZ.component(vector::Y);
  70
  71 volScalarField gradGradUYZX = gradGradUcompZX.component(vector::Y);
  72 volScalarField gradGradUYZY = gradGradUcompZY.component(vector::Y);
  73 volScalarField gradGradUYZZ = gradGradUcompZZ.component(vector::Y);
  74
  75 volScalarField gradGradUZXX = gradGradUcompXX.component(vector::Z);
  76 volScalarField gradGradUZXY = gradGradUcompXY.component(vector::Z);
  77 volScalarField gradGradUZXZ = gradGradUcompXZ.component(vector::Z);
  78
  79 volScalarField gradGradUZYX = gradGradUcompYX.component(vector::Z);
  80 volScalarField gradGradUZYY = gradGradUcompYY.component(vector::Z);
  81 volScalarField gradGradUZYZ = gradGradUcompYZ.component(vector::Z);
  82
  83 volScalarField gradGradUZZX = gradGradUcompZX.component(vector::Z);
  84 volScalarField gradGradUZZY = gradGradUcompZY.component(vector::Z);
  85 volScalarField gradGradUZZZ = gradGradUcompZZ.component(vector::Z);
  86
  87 forAll(extrapGradU.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& gradUOwn = gradU.internalField()[own];
  94   const tensor& gradUNei = gradU.internalField()[nei];
  95
  96   // as there is there is no thirdOrderTensor class, we will
  97   // calculate (deltaOwn&gradGradUOwn) out manually
  98   // tensor deltaOwnDotgradGradUOwn = tensor::zero;
  99   // tensor deltaNeiDotgradGradUNei = tensor::zero;
 100
 101   // deltaOwnDotgradGradUOwn[tensor::XX] = deltaOwn & gradGradUcompXX.internalField()[own];
 102   // deltaNeiDotgradGradUNei[tensor::XX] = deltaNei & gradGradUcompXX.internalField()[nei];
 103   // deltaOwnDotgradGradUOwn[tensor::XY] = deltaOwn & gradGradUcompXY.internalField()[own];
 104   // deltaNeiDotgradGradUNei[tensor::XY] = deltaNei & gradGradUcompXY.internalField()[nei];
 105   // deltaOwnDotgradGradUOwn[tensor::XZ] = deltaOwn & gradGradUcompXZ.internalField()[own];
 106   // deltaNeiDotgradGradUNei[tensor::XZ] = deltaNei & gradGradUcompXZ.internalField()[nei];
 107
 108   // deltaOwnDotgradGradUOwn[tensor::YX] = deltaOwn & gradGradUcompYX.internalField()[own];
 109   // deltaNeiDotgradGradUNei[tensor::YX] = deltaNei & gradGradUcompYX.internalField()[nei];
 110   // deltaOwnDotgradGradUOwn[tensor::YY] = deltaOwn & gradGradUcompYY.internalField()[own];
 111   // deltaNeiDotgradGradUNei[tensor::YY] = deltaNei & gradGradUcompYY.internalField()[nei];
 112   // deltaOwnDotgradGradUOwn[tensor::YZ] = deltaOwn & gradGradUcompYZ.internalField()[own];
 113   // deltaNeiDotgradGradUNei[tensor::YZ] = deltaNei & gradGradUcompYZ.internalField()[nei];
 114
 115   // deltaOwnDotgradGradUOwn[tensor::ZX] = deltaOwn & gradGradUcompZX.internalField()[own];
 116   // deltaNeiDotgradGradUNei[tensor::ZX] = deltaNei & gradGradUcompZX.internalField()[nei];
 117   // deltaOwnDotgradGradUOwn[tensor::ZY] = deltaOwn & gradGradUcompZY.internalField()[own];
 118   // deltaNeiDotgradGradUNei[tensor::ZY] = deltaNei & gradGradUcompZY.internalField()[nei];
 119   // deltaOwnDotgradGradUOwn[tensor::ZZ] = deltaOwn & gradGradUcompZZ.internalField()[own];
 120   // deltaNeiDotgradGradUNei[tensor::ZZ] = deltaNei & gradGradUcompZZ.internalField()[nei];
 121
 122   scalar gradGradUXXXOwn = gradGradUXXX.internalField()[own];
 123   scalar gradGradUXXYOwn = gradGradUXXY.internalField()[own];
 124   scalar gradGradUXXZOwn = gradGradUXXZ.internalField()[own];
 125
 126   scalar gradGradUXYXOwn = gradGradUXYX.internalField()[own];
 127   scalar gradGradUXYYOwn = gradGradUXYY.internalField()[own];
 128   scalar gradGradUXYZOwn = gradGradUXYZ.internalField()[own];
 129
 130   scalar gradGradUXZXOwn = gradGradUXZX.internalField()[own];
 131   scalar gradGradUXZYOwn = gradGradUXZY.internalField()[own];
 132   scalar gradGradUXZZOwn = gradGradUXZZ.internalField()[own];
 133
 134
 135   scalar gradGradUYXXOwn = gradGradUYXX.internalField()[own];
 136   scalar gradGradUYXYOwn = gradGradUYXY.internalField()[own];
 137   scalar gradGradUYXZOwn = gradGradUYXZ.internalField()[own];
 138
 139   scalar gradGradUYYXOwn = gradGradUYYX.internalField()[own];
 140   scalar gradGradUYYYOwn = gradGradUYYY.internalField()[own];
 141   scalar gradGradUYYZOwn = gradGradUYYZ.internalField()[own];
 142
 143   scalar gradGradUYZXOwn = gradGradUYZX.internalField()[own];
 144   scalar gradGradUYZYOwn = gradGradUYZY.internalField()[own];
 145   scalar gradGradUYZZOwn = gradGradUYZZ.internalField()[own];
 146
 147
 148   scalar gradGradUZXXOwn = gradGradUZXX.internalField()[own];
 149   scalar gradGradUZXYOwn = gradGradUZXY.internalField()[own];
 150   scalar gradGradUZXZOwn = gradGradUZXZ.internalField()[own];
 151
 152   scalar gradGradUZYXOwn = gradGradUZYX.internalField()[own];
 153   scalar gradGradUZYYOwn = gradGradUZYY.internalField()[own];
 154   scalar gradGradUZYZOwn = gradGradUZYZ.internalField()[own];
 155
 156   scalar gradGradUZZXOwn = gradGradUZZX.internalField()[own];
 157   scalar gradGradUZZYOwn = gradGradUZZY.internalField()[own];
 158   scalar gradGradUZZZOwn = gradGradUZZZ.internalField()[own];
 159
 160
 161   // nei
 162   scalar gradGradUXXXNei = gradGradUXXX.internalField()[nei];
 163   scalar gradGradUXXYNei = gradGradUXXY.internalField()[nei];
 164   scalar gradGradUXXZNei = gradGradUXXZ.internalField()[nei];
 165
 166   scalar gradGradUXYXNei = gradGradUXYX.internalField()[nei];
 167   scalar gradGradUXYYNei = gradGradUXYY.internalField()[nei];
 168   scalar gradGradUXYZNei = gradGradUXYZ.internalField()[nei];
 169
 170   scalar gradGradUXZXNei = gradGradUXZX.internalField()[nei];
 171   scalar gradGradUXZYNei = gradGradUXZY.internalField()[nei];
 172   scalar gradGradUXZZNei = gradGradUXZZ.internalField()[nei];
 173
 174
 175   scalar gradGradUYXXNei = gradGradUYXX.internalField()[nei];
 176   scalar gradGradUYXYNei = gradGradUYXY.internalField()[nei];
 177   scalar gradGradUYXZNei = gradGradUYXZ.internalField()[nei];
 178
 179   scalar gradGradUYYXNei = gradGradUYYX.internalField()[nei];
 180   scalar gradGradUYYYNei = gradGradUYYY.internalField()[nei];
 181   scalar gradGradUYYZNei = gradGradUYYZ.internalField()[nei];
 182
 183   scalar gradGradUYZXNei = gradGradUYZX.internalField()[nei];
 184   scalar gradGradUYZYNei = gradGradUYZY.internalField()[nei];
 185   scalar gradGradUYZZNei = gradGradUYZZ.internalField()[nei];
 186
 187
 188   scalar gradGradUZXXNei = gradGradUZXX.internalField()[nei];
 189   scalar gradGradUZXYNei = gradGradUZXY.internalField()[nei];
 190   scalar gradGradUZXZNei = gradGradUZXZ.internalField()[nei];
 191
 192   scalar gradGradUZYXNei = gradGradUZYX.internalField()[nei];
 193   scalar gradGradUZYYNei = gradGradUZYY.internalField()[nei];
 194   scalar gradGradUZYZNei = gradGradUZYZ.internalField()[nei];
 195
 196   scalar gradGradUZZXNei = gradGradUZZX.internalField()[nei];
 197   scalar gradGradUZZYNei = gradGradUZZY.internalField()[nei];
 198   scalar gradGradUZZZNei = gradGradUZZZ.internalField()[nei];
 199
 200   tensor deltaOwnDotgradGradUOwn =
 201     tensor(
 202        deltaOwn.x()*gradGradUXXXOwn + deltaOwn.y()*gradGradUYXXOwn + deltaOwn.z()*gradGradUZXXOwn,
 203        deltaOwn.x()*gradGradUXXYOwn + deltaOwn.y()*gradGradUYXYOwn + deltaOwn.z()*gradGradUZXYOwn,
 204        deltaOwn.x()*gradGradUXXZOwn + deltaOwn.y()*gradGradUYXZOwn + deltaOwn.z()*gradGradUZXZOwn,
 205
 206        deltaOwn.x()*gradGradUXYXOwn + deltaOwn.y()*gradGradUYYXOwn + deltaOwn.z()*gradGradUZYXOwn,
 207        deltaOwn.x()*gradGradUXYYOwn + deltaOwn.y()*gradGradUYYYOwn + deltaOwn.z()*gradGradUZYYOwn,
 208        deltaOwn.x()*gradGradUXYZOwn + deltaOwn.y()*gradGradUYYZOwn + deltaOwn.z()*gradGradUZYZOwn,
 209
 210        deltaOwn.x()*gradGradUXZXOwn + deltaOwn.y()*gradGradUYZXOwn + deltaOwn.z()*gradGradUZZXOwn,
 211        deltaOwn.x()*gradGradUXZYOwn + deltaOwn.y()*gradGradUYZYOwn + deltaOwn.z()*gradGradUZZYOwn,
 212        deltaOwn.x()*gradGradUXZZOwn + deltaOwn.y()*gradGradUYZZOwn + deltaOwn.z()*gradGradUZZZOwn
 213        );
 214
 215   tensor deltaNeiDotgradGradUNei =
 216     tensor(
 217        deltaNei.x()*gradGradUXXXNei + deltaNei.y()*gradGradUYXXNei + deltaNei.z()*gradGradUZXXNei,
 218        deltaNei.x()*gradGradUXXYNei + deltaNei.y()*gradGradUYXYNei + deltaNei.z()*gradGradUZXYNei,
 219        deltaNei.x()*gradGradUXXZNei + deltaNei.y()*gradGradUYXZNei + deltaNei.z()*gradGradUZXZNei,
 220
 221        deltaNei.x()*gradGradUXYXNei + deltaNei.y()*gradGradUYYXNei + deltaNei.z()*gradGradUZYXNei,
 222        deltaNei.x()*gradGradUXYYNei + deltaNei.y()*gradGradUYYYNei + deltaNei.z()*gradGradUZYYNei,
 223        deltaNei.x()*gradGradUXYZNei + deltaNei.y()*gradGradUYYZNei + deltaNei.z()*gradGradUZYZNei,
 224
 225        deltaNei.x()*gradGradUXZXNei + deltaNei.y()*gradGradUYZXNei + deltaNei.z()*gradGradUZZXNei,
 226        deltaNei.x()*gradGradUXZYNei + deltaNei.y()*gradGradUYZYNei + deltaNei.z()*gradGradUZZYNei,
 227        deltaNei.x()*gradGradUXZZNei + deltaNei.y()*gradGradUYZZNei + deltaNei.z()*gradGradUZZZNei
 228        );
 229
 230
 231   // get average of extrapolated values
 232   tensor extrapFaceGrad =
 233     0.5*
 234     (
 235      gradUOwn + (deltaOwnDotgradGradUOwn)
 236      +
 237      gradUNei + (deltaNeiDotgradGradUNei)
 238      );
 239
 240   extrapGradU.internalField()[facei] = extrapFaceGrad;
 241
 242 }
 243
 244 // correction is zero on boundary
 245 forAll(extrapGradU.boundaryField(), patchi)
 246 {
 247   extrapGradU.boundaryField()[patchi] =
 248     averageGradU.boundaryField()[patchi];
 249 }
 250
 251 // calculate thirdOrderTerm
 252 volVectorField divThirdOrderTerm (
 253                   "thirdOrderTerm",
 254                   fvc::div(
 255                        (2*muf+lambdaf)*mesh.Sf()
 256                        & (extrapGradU - averageGradU)
 257                        )
 258                   );
 259
 260 // if(runTime.outputTime())
 261 //   {
 262 //     divThirdOrderTerm.write();
 263 //     averageGradU.write();
 264 //     extrapGradU.write();
 265 //   }