applications/solvers/solidMechanics/icoFsiElasticNonLinULSolidFoam/calculateLiftAndDrag.H

   1 {
   2 //- calculate lift and drag for plate and cylinder due to pressure (ignore shear stress)
   3 label plateID = mesh.boundaryMesh().findPatchID("plate");
   4 label cylinderID = mesh.boundaryMesh().findPatchID("cylinder");
   5
   6 if(plateID == -1 || cylinderID == -1)
   7   {
   8     FatalError << "\n Cannot find the plate patch or the cylinder"
   9                << " patch to calculate lift and drag!"
  10                << exit(FatalError);
  11   }
  12
  13 scalar lift = 0;
  14 scalar drag = 0;
  15
  16  const vectorField& Sfp = mesh.boundary()[plateID].Sf();
  17 forAll(p.boundaryField()[plateID], facei)
  18 {
  19   vector faceForce = p.boundaryField()[plateID][facei] * Sfp[facei];
  20   lift += vector(0,1,0) & faceForce;
  21   drag += vector(1,0,0) & faceForce;
  22 }
  23
  24  const vectorField& Sfc = mesh.boundary()[cylinderID].Sf();
  25 forAll(p.boundaryField()[cylinderID], facei)
  26 {
  27   vector faceForce = p.boundaryField()[cylinderID][facei] * Sfc[facei];
  28   lift += vector(0,1,0) & faceForce;
  29   drag += vector(1,0,0) & faceForce;
  30 }
  31
  32  scalar width = 0.050668;
  33
  34 Info << "Total lift on the cylinder and plate boundaries is " << lift << " N, per unit width is " << (lift/width) << " N\n"
  35      << "Total drag on the cylinder and plate boundaries is " << drag << " N, per unit width is " << (drag/width) << " N\n"
  36      << endl;
  37 }