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

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

Dear Paul,

Thank you for the reply. I tried to implement the code in the usrchk
subroutine but found that the values of dgtq and dg are zeros. I tried to
identify the wall boundary conditions and found the corresponding the
elements , faces and used the code snippet. I am sure I am missing something
here.. Any thoughts ?

However the way you said with the postx works. It did give me a profile.out
and I am currently working on a way to post-process it. I have close to 25k
elements and probably a few hundred elements near the blade and I am not
quite sure how to grab the coordinates to enter into the prof.in file.
Ideally I would like to get the yplus/z plus for the first couple of
elements or so. I tried the same way of grabbing the elements and faces near
the area of interest, but I am clueless as to how I access the co-ordinate
data of those elements.

I had another question regarding span-wise averaging of shear stress, Once I
calculate the shear stress using the above routine, could I use the
z_average() routine in navier5.f to do the span-wise averaging ? I was
looking at the turbchannel.usr where the wall shear stress and averaging
process is done but did not make a huge progress..

<code>
      subroutine usrchk
      include 'SIZE'
      include 'TOTAL'
      include 'ZPER'

      integer e,f,pf
      real x0(3)
      character*3 cb_surf
      real    n1,n2,n3
      real dg(3,2)
      real dgtq(3,4)
      real xm0 (lx1,ly1,lz1,lelt)
      real ym0 (lx1,ly1,lz1,lelt)
      real zm0 (lx1,ly1,lz1,lelt)
      real sij (lx1,ly1,lz1,3*ldim-3,lelv)
      real pm1 (lx1,ly1,lz1,lelv)
      real visc(lx1,ly1,lz1,lelv)

       n=nx1*ny1*nz1*nelv

      call rzero(dgtq,12)

      nface = 2 * ndim
      cb_surf = 'W '
      do e = 1,nelv
         flag = 1
         do f = 1,nface
      if(cbc(f,e,1).eq.cb_surf.and.x.gt.-25E-02 .and. x.lt.125E-02)then
      flag = 0                  ! Wall Element
      pf     = eface1(f)         ! convert from preproc. notation

      js1    = skpdat(1,pf)
      jf1    = skpdat(2,pf)
      jskip1 = skpdat(3,pf)
      js2    = skpdat(4,pf)
      jf2    = skpdat(5,pf)
      jskip2 = skpdat(6,pf)

        if (if3d.or.ifaxis) then
       i = 0
       a = 0
       do j2=js2,jf2,jskip2
       do j1=js1,jf1,jskip1
         i = i+1
         n1 = unx(i,1,f,e)*area(i,1,f,e)
         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)
         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)
c
         r1 = xm0(j1,j2,1,e)
         r2 = ym0(j1,j2,1,e)
         r3 = zm0(j1,j2,1,e)
c
         do l=1,2
         do k=1,3
            dgtq(k,l) = dgtq(k,l) + dg(k,l)
         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
        endif
        endif
        enddo
        enddo

      return
      end

</code>


Thanks and Regards
Shriram
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20100414/723fc290/attachment.html>