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 }