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

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

For the postx approach, you don't need points associated w/ the
mesh, just points that are on the surface of your geometry.

(Sorry for the short response, but I have to leave the office
right now...)

Paul




On Wed, 14 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:

> 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
>