[Nek5000-users] Help: compute dissipation

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

Hi NEK users,

I am trying to calculate dissipation for fully developed turbulent channel flow. I tried the code below, but result I got from this totally disagrees with Vreman’s 2014 paper. x-axis is streamwise, y-axis is wall-normal and z-axis is spanwise. Any errors in the code? Any thoughts?

     real diss_sqr(lx1,ly1,lz1, lelt)
      real s11
      real s12
      real s13
      real s22
      real s23
      real s33


  common /plane/  uavg_pl(ly1*lely)
     $             ,  vavg_pl(ly1*lely)
     $             ,  wavg_pl(ly1*lely)
     $             ,  urms_pl(ly1*lely)
     $             ,  vrms_pl(ly1*lely)
     $             ,  wrms_pl(ly1*lely)
     $             ,  uvms_pl(ly1*lely)
     $             ,  diss_sqr_pl(ly1*lely)
     $             ,  yy(ly1*lely)
     $             ,  w1(ly1*lely),w2(ly1*lely)
     $             ,  ffx_avg, dragx_avg


parameter (lr=lx1*ly1*lz1)
      common /scruz/  ur(lr),us(lr),ut(lr)
     $             ,  vr(lr),vs(lr),vt(lr)
     $             ,  wr(lr),ws(lr),wt(lr)


common /avg/    uavg(lx1,ly1,lz1,lelv)
     &             ,  vavg(lx1,ly1,lz1,lelv)
     &             ,  wavg(lx1,ly1,lz1,lelv)
     &             ,  urms(lx1,ly1,lz1,lelv)
     &             ,  vrms(lx1,ly1,lz1,lelv)
     &             ,  wrms(lx1,ly1,lz1,lelv)
     &             ,  uvms(lx1,ly1,lz1,lelv)
     &             ,  dissavg(lx1,ly1,lz1,lelv)


common /scrns/  sij (lx1*ly1*lz1,6,lelv)

if(icalld.eq.0) then
        call rzero(uavg,n)
        call rzero(urms,n)
        call rzero(vrms,n)
        call rzero(wrms,n)
        call rzero(uvms,n)
        call rzero(dissavg,n)
        call rzero(wavg,n)
        call rzero(vavg,n)
        call rzero(diss_sqr,n)


call comp_sij(sij,6,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)

       do e=1,nelt
          do i = 1,nx1*ny1*nz1
            s11=sij(i,1,e)
            s12=sij(i,4,e)
            s22=sij(i,2,e)
            s13=sij(i,6,e)
            s23=sij(i,5,e)
            s33=sij(i,3,e)

          diss_sqr(i,1,1,e) = (s11*s11) +
     $                       (s22*s22) +
     $                       (s33*s33) +
     $                       (s12*s12)*2 + (s23*s23)*2 +
     $                       (s13*s13)*2

         enddo
         enddo

call avg1(dissavg,diss_sqr,alpha,beta,n,'diss',ifverbose)

call planar_average_s(diss_sqr_pl,dissavg,w1,w2)

do i=1,ny1*nely/2
           uavg_pl(i) = 0.5 * (uavg_pl(i) + uavg_pl(m-i+1))
           vavg_pl(i) = 0.5 * (vavg_pl(i) + vavg_pl(m-i+1))
           wavg_pl(i) = 0.5 * (wavg_pl(i) + wavg_pl(m-i+1))
           diss_sqr_pl(i) = 0.5 * (diss_sqr_pl(i) + diss_sqr_pl(m-i+1))
           urms_pl(i) = 0.5 * (urms_pl(i) + urms_pl(m-i+1))
           vrms_pl(i) = 0.5 * (vrms_pl(i) + vrms_pl(m-i+1))
           wrms_pl(i) = 0.5 * (wrms_pl(i) + wrms_pl(m-i+1))

       enddo



Many thanks.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20170623/331b477e/attachment-0001.html>