Integrating fluid-stress over a surface to get forces (drag)

Question

Hello,

I was wondering what is the correct way of going from flow velocity and pressure to drag and lift forces using FEniCS code. I am trying this, right now, on the typical 2D cylinder in a channel flow simulation. I see the Von Karman vortex street forming very well, and now I am interested in the drag and lift.

I have done the usual marking of boundaries, and have extracted the boundary for the circular section (for the cylinder) as:

domainBoundaries = FacetFunction("size_t", mesh)
domainBoundaries.set_all(0)
.... marking other walls ...
cylinderwall.mark(domainBoundaries,5)
ds = Measure("ds")[domainBoundaries]
GammaP = ds(5)

and I have the following function defined to calculate the total drag and lift forces:

def integrateFluidStress(a_V, a_P, a_Nu, a_Normal, a_GammaP, a_Dim = 2):

    eps   = 0.5*(nabla_grad(a_V) + nabla_grad(a_V).T)
    sig   = -a_P*Identity(len(a_V)) + 2.0*a_Nu*eps

    traction  = dot(sig, a_Normal)

    if a_Dim == 2:

       forceX  = traction[0]*a_GammaP
       forceY  = traction[1]*a_GammaP
       fX      = assemble(forceX)
       fY      = assemble(forceY)

      return (fX, fY)

    elif a_Dim == 3:

       forceX  = traction[0]*a_GammaP
       forceY  = traction[1]*a_GammaP
       forceZ  = traction[2]*a_GammaP
       fX      = assemble(forceX)
       fY      = assemble(forceY)
       fZ      = assemble(forceZ)

      return (fX, fY, fZ)

I then calculate the forces by using the following:

(fX,fY)     = integrateFluidStress(v, p, Nu, n, GammaP)

where v and p are the solutions obtained from the incompressible Navier Stokes. However, after doing this, I do not get the drag and lift forces to match available literature (in fact the mean drag is coming out to be 0), even though the vortex patterns match.

Is there something that I am doing wrong while calculating the integrals? I am also looking for perhaps other (more efficient?) ways of doing this fore calculation, if any.

Any help will be much appreciated.

Comments

Do you have a complete example? Your code looks correct, but a hands on example might be good to find potential hidden flaws.

---

I agree that the code looks correct. Have you tried some basic debugging like plotting the traction before the integrals are calculated?

---

Hello Christian Waluga and Gregor Mitscha-Eibi - Thanks for your response.
I figured out the issue, it had to do with correctly updating the surface normal to the cylinder in my code. There were multiple surface normals, and I was not updating the correct normals to be used for force calculation.

The code now works correctly (i think ! unless any other bug shows up :) )

Thanks