[philippose]

Hi,
I just happened to browse through this dialogue regarding calculating wall shear stresses.

In the simple case, if you take laminar flow, the wall shear stress is calculated using the equation:

Tau (shear stress) = mu * (dU/dy)

In this equation, "mu" is the dynamic viscosity of the medium, and "(dU/dy)" is the change in velocity "U" with change in "y", which is the distance from the wall in a direction normal to the wall.

Now, in the case of incompressible solvers in openFOAM, the whole system is normalised with respect to the fluid density "rho".... so in this case, the dynamic viscosity (mu) is given by:

mu = nu * rho

Where "nu" is the kinematic viscosity, defined in the "transportProperties" dictionary, and "rho" would depend on the fluid you are working with.

Now, to implement this in openFOAM, there exists a class member called "snGrad". This directly gives you the gradient of velocity (dU/dy) normal to each face of the wall patches.

So, you would use:

wallShear = nu * rho * mesh.boundaryField()[patchID].snGrad()

This directly gives you the wall shear stress (assuming the "patchID" refers to a wall patch ofcourse). Here, you can read "nu" from the transportProperties dictionary, and you will have to define "rho" somewhere.

When the case is not laminar, but turbulent, then things look a bit different, which is why you may have gotten confused when you looked through the wallShearStress source-code.

In this case, the Reynolds Stress Tensor gives you the shear stress, but not necessarily normal to the face, and you need to resolve the tensor along the normal to the face in order to obtain the wall shear stress.

This is because in the case of turbulent simulations, the viscosity is not a constant, and is an effective value which is a result of the turbulence model, and the Reynolds Stress Tensor.

So, to calculate the wall shear stress in the case of a turbulent simulation, you need to do something like:

wallShear = (-mesh.Sf().boundaryField()[patchID]/(mag(mesh.Sf().boundaryField()[patchID])) ) & (turbulence->R()().boundaryField()[patchID])

The above equation will give you the wall shear stress with the velocity gradient resolved in the direction normal to each face in the patch.

Have a nice day!

Philippose