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 = tensor
 201     (
 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 = tensor
 216     (
 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 (
 254     "thirdOrderTerm",
 255     fvc::div
 256     (
 257         (2*muf+lambdaf)*mesh.Sf() & (extrapGradDU - averageGradDU)
 258     )
 259 );
 260
 261 // if(runTime.outputTime())
 262 // {
 263 //     divThirdOrderTerm.write();
 264 //     averageGradDU.write();
 265 //     extrapGradDU.write();
 266 // }