applications/solvers/solidMechanics/icoFsiElasticNonLinULSolidFoam/setInterfaceForce.H

   1 {
   2     Info << "Setting traction on solid patch" << endl;
   3
   4 //     vectorField fluidPatchTraction =
   5 //        -rhoFluid.value()*nu.value()
   6 //        *U.boundaryField()[fluidPatchID].snGrad()
   7 //       + rhoFluid.value()*p.boundaryField()[fluidPatchID]
   8 //        *mesh.boundary()[fluidPatchID].nf();
   9
  10     vectorField fluidPatchTraction =
  11        -rhoFluid.value()*nu.value()
  12        *U.boundaryField()[fluidPatchID].snGrad();
  13
  14     scalarField fluidPatchPressure =
  15         rhoFluid.value()*p.boundaryField()[fluidPatchID];
  16
  17     vectorField fluidZoneTraction
  18         (
  19             mesh.faceZones()[fluidZoneID].size(),
  20             vector::zero
  21         );
  22
  23     const label fluidPatchStart =
  24         mesh.boundaryMesh()[fluidPatchID].start();
  25
  26     forAll(fluidPatchTraction, i)
  27     {
  28         fluidZoneTraction
  29         [
  30             mesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
  31         ] =
  32             fluidPatchTraction[i];
  33     }
  34
  35     // Parallel data exchange: collect pressure field on all processors
  36     reduce(fluidZoneTraction, sumOp<vectorField>());
  37
  38
  39     scalarField fluidZonePressure
  40         (
  41             mesh.faceZones()[fluidZoneID].size(),
  42             0.0
  43         );
  44
  45     forAll(fluidPatchPressure, i)
  46     {
  47         fluidZonePressure
  48         [
  49             mesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
  50         ] =
  51             fluidPatchPressure[i];
  52     }
  53
  54     // Parallel data exchange: collect pressure field on all processors
  55     reduce(fluidZonePressure, sumOp<scalarField>());
  56
  57     vectorField solidZoneTraction =
  58         interpolatorFluidSolid.faceInterpolate
  59         (
  60             fluidZoneTraction
  61         );
  62
  63     scalarField solidZonePressure =
  64         interpolatorFluidSolid.faceInterpolate
  65         (
  66             fluidZonePressure
  67         );
  68
  69     const label solidPatchStart =
  70         stressMesh.boundaryMesh()[solidPatchID].start();
  71
  72     forAll(solidPatchTraction, i)
  73     {
  74         solidPatchTraction[i] =
  75             solidZoneTraction
  76             [
  77                 stressMesh.faceZones()[solidZoneID]
  78                .whichFace(solidPatchStart + i)
  79             ];
  80     }
  81
  82     forAll(solidPatchPressure, i)
  83     {
  84         solidPatchPressure[i] =
  85             solidZonePressure
  86             [
  87                 stressMesh.faceZones()[solidZoneID]
  88                .whichFace(solidPatchStart + i)
  89             ];
  90     }
  91
  92     if (fsi)
  93     {
  94         tForce.traction() = solidPatchTraction;
  95         tForce.pressure() = solidPatchPressure;
  96     }
  97
  98     vector totalTractionForce =
  99         gSum
 100         (
 101             solidPatchTraction
 102            *stressMesh.magSf().boundaryField()[solidPatchID]
 103         );
 104
 105     Info << "Total traction force = "
 106         << totalTractionForce << endl;
 107 }