Calculation of traction force

Question

Hi,
I have calculated my shear stress 'sigma' using the pressure and velocity calculated from Navier-Stokes solver as follows,

n = FacetNormal(mesh)
W = TensorFunctionSpace(mesh, "Lagrange", 2)
sigma = project(-p*I + 2.0*nu*sym(grad(u)),W)
Traction = dot(sigma,n) #sigma*n

My problem is that I would like to know the Traction (force per unit area) as a vector Tx,Ty,Tz at all points in the mesh. How do I do this? when I try to print Traction[0] or Traction[1], I get ,

((({ A | A_{i_{10}, i_{11}} = I[i_{10}, i_{11}] * -1 * f_60[2] }) + ({ A | A_{i_{12}, i_{13}} = (sym(grad([f_60[0], f_60[1]])))[i_{12}, i_{13}] * 2.0 * f_29 })) . (n))[0]

instead of values.

Thanks

Answer 1

The normal vector and traction force are only well-defined on the boundary of the domain. Do you want to plot the vector on just the boundary?

Comments

Thank you for your reply, Kristian.I am using complex boundaries which are made of small particles clubbed together (a cluster of sand-particles) in my domain.Hence I wont be able to mark my domain(s) . I want to know the force exerted by the fluid in x,y,z directions on each particle.

After some online searching, I found,

scalar1=FunctionSpace(mesh,"CG",1)
normal_stresses = Function(scalar1)
v=TestFunction(scalar1)
ds = Measure("ds")
fx=(1/FacetArea(mesh))*v*Traction[0]*ds
fy=(1/FacetArea(mesh))*v*Traction[1]*ds
Traction_x=assemble(fx,tensor=normal_stresses.vector())

and now I am able to get some values for the force in all the cells (value 0 everywhere other than the boundaries). But I really do not understand why the traction is multiplied by surface ds and divided by FacetArea. Is this approach correct? Any hint is greatly appreciated.

Thanks

---

If they are spherical, you should be able to mark them with boolean functions.

If they are not touching, you could use the mesh topology to mark them.

I do not see a way around marking, when you want the force on every particle. Do you plan to move the particles?

---

There are many chunks of particles. Say, around 40 spheres immersed in the fluid, and each sphere made of around 100 particles. I can mark all the particles with the same domain number. The particles will be moved after a few time steps based on some external calculations.

Thanks

---

How do you represent the particles/spheres in FEniCS?

---

I apply Dirichlet boundary condition in the points where the particles exist.
pts=np.loadtxt("particlecoordinates.dat")
b_v = DirichletBC(V, Constant((0.0, 0.0,0.0)),Pinpoint(pts))

---

Do you not need vertices at the particle coordinates then?

If the particles are moving, should a non-zero velocity not be prescribed?

Has this numerical method been used in your field before? Has it been verified against references?

---

For now I am just marking the vertices nearest to the particle positions. It is an approximation.
The particles are moved only after a fixed time interval. This method has not yet been established. We are experimenting with our model and found that fenics suits our needs.

Thanks