[Nek5000-users] Calculating wall shear stress and dimensionless wall distances

[nek5000-users at lists.mcs.anl.gov]

One option is to simply compute du_i/dx_j and then contract this
with the unit normal on the walls in question.

You could then also dump out the x,y,z, tau_w at all the points
on the walls of interest...

The one unfortunate thing here is that we don't have  a great
mechanism for interrogating/visualizing surface structure data.
(I usually resort to matlab or gnuplot at that point...)

Note also that you can define a quantity, let's call it u_tan,
for all elements that have a "W  " (wall) boundary condition.

This would be a volumetric quantity for only those elements having
W bc and would be defined as:

     U - (U.nhat)/||U|| U

(which I _think_ is right, but should be checked...).  Differentiating
this quantity should give the shear stress at the wall.  Another trick
that works in limited cases (i.e., fixed surface) is vort.n_hat.

After searching around a bit, probably the best example of how
to compute tau_w comes from  navier5.f in the /nek directory,
subruoutine drgtrq().  Here is an annotated code snippet:


c
       if (if3d.or.ifaxis) then
        i = 0
        a = 0
        do j2=js2,jf2,jskip2   ! SCAN POINTS ON FACE f of ELEMENT e
        do j1=js1,jf1,jskip1
          i = i+1
          n1 = unx(i,1,f,e)*area(i,1,f,e)  ! unit normal x surface Jacobian
          n2 = uny(i,1,f,e)*area(i,1,f,e)
          n3 = unz(i,1,f,e)*area(i,1,f,e)
          a  = a +          area(i,1,f,e)
c
          v  = visc(j1,j2,1,e)
c
          s11 = sij(j1,j2,1,1,e)           ! Sij at point of interest
          s21 = sij(j1,j2,1,4,e)
          s31 = sij(j1,j2,1,6,e)
c
          s12 = sij(j1,j2,1,4,e)
          s22 = sij(j1,j2,1,2,e)
          s32 = sij(j1,j2,1,5,e)
c
          s13 = sij(j1,j2,1,6,e)
          s23 = sij(j1,j2,1,5,e)
          s33 = sij(j1,j2,1,3,e)
c
          dg(1,1) = pm1(j1,j2,1,e)*n1     ! pressure drag
          dg(2,1) = pm1(j1,j2,1,e)*n2
          dg(3,1) = pm1(j1,j2,1,e)*n3
c
          dg(1,2) = -v*(s11*n1 + s12*n2 + s13*n3) ! viscous drag
          dg(2,2) = -v*(s21*n1 + s22*n2 + s23*n3)
          dg(3,2) = -v*(s31*n1 + s32*n2 + s33*n3) ! Sij contracted with n_j
c
          r1 = xm0(j1,j2,1,e)                     ! Position relative
          r2 = ym0(j1,j2,1,e)                     ! to where torque is desired
          r3 = zm0(j1,j2,1,e)
c
          do l=1,2
          do k=1,3
             dgtq(k,l) = dgtq(k,l) + dg(k,l)              ! Drag
          enddo
          enddo
c
          dgtq(1,3) = dgtq(1,3) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure
          dgtq(2,3) = dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque
          dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))
c
          dgtq(1,4) = dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous
          dgtq(2,4) = dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque
          dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(1,2))
        enddo
        enddo


-- Paul



On Mon, 12 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:

> Hello,
>
> Background : I am trying to perform a LES in a turbine blade where the
> stream-wise, pitch-wise and span-wise direction is along X,Y and Z axis.
>
> Question : I would like to calculate the yplus and zplus values for the
> nodes close to the blade. For calculating the dimensionless wall distances,
> one would require the velocity gradients in the respective direction or the
> wall shear stress . I would like to know if there are any subroutines
> available in nek that would help calculate the wall normal distances or the
> velocity gradients ?
>
> Any help is greatly appreciated. Thanks
>
> Regards
> Shriram Jagannathan
>