[Nek5000-users] About the iobj in the drag calculation

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

Hi,Paul:
     I'm confused about part of the subroutine userchk in turbChannel.usr
     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)
        atime = 0.
        timel = time
        call planar_average_s(yy     ,ym1 ,w1,w2)
        icalld = 1
      endif

      dtime = time - timel
      atime = atime + dtime

      if (atime.ne.0. .and. dtime.ne.0.) then
        beta      = dtime/atime
        alpha     = 1.-beta
        ifverbose = .false.

        call avg1(uavg,vx   ,alpha,beta,n,'uavg',ifverbose)
        call avg2(urms,vx   ,alpha,beta,n,'urms',ifverbose)
        call avg2(vrms,vy   ,alpha,beta,n,'vrms',ifverbose)
        call avg2(wrms,vz   ,alpha,beta,n,'wrms',ifverbose)
        call avg3(uvms,vx,vy,alpha,beta,n,'uvmm',ifverbose)

        dragx_avg = alpha*dragx_avg + beta*0.5*(dragx(1)+dragx(2))

        ! averaging over statistical homogeneous directions (r-t)
        call planar_average_s(uavg_pl,uavg,w1,w2)
        call planar_average_s(urms_pl,urms,w1,w2)
        call planar_average_s(vrms_pl,vrms,w1,w2)
        call planar_average_s(wrms_pl,wrms,w1,w2)
        call planar_average_s(uvms_pl,uvms,w1,w2)

        ! average over half the channel height
        m = ny1*nely
        do i=1,ny1*nely/2
           uavg_pl(i) = 0.5 * (uavg_pl(i) + uavg_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
      endif

      tw     = dragx_avg/A_w + 1.e-50
      u_tau  = sqrt(tw/rho)
      Re_tau = u_tau*delta/dnu

      diff   = (tw-ffx_new)/tw
      if(nid.eq.0) write(6,1) istep,time,e2,ffx_new,diff
    1              format(i6,1p4e17.9,' err2')

      ! write statistics to file
      iostep_avg = param(68)
      if(nid.eq.0 .and. istep.gt.0 .and.
     &   mod(istep,iostep_avg).eq.0) then
            write(6,*) 'Dumping statistics ...'
            open(unit=56,file='reystresses.dat')
            write(56,'(A,1pe14.7)') '#time = ', time
            write(56,'(A)')
     &      '#  y     y+     R_uu     R_vv     R_ww     R_uv'
            open(unit=57,file='means.dat')
            write(57,'(A,1pe14.7)') '#time = ', time
            write(57,'(A)')
     &      '#  y     y+    Umean'

            m = ny1*nely/2
            do i=1,m
                write(56,3) yy(i)+1
     &                     ,(yy(i)+1)*Re_tau
     &                     , (urms_pl(i)-(uavg_pl(i))**2)/u_tau**2
     &                     , vrms_pl(i)/u_tau**2
     &                     , wrms_pl(i)/u_tau**2
     &                     , uvms_pl(i)/u_tau**2
                write(57,3)  yy(i) + 1.
     &                     , (yy(i)+1.)*Re_tau
     &                     , uavg_pl(i)/u_tau
    3           format(1p15e17.9)
            enddo

            close(56)
            close(57)
      endif

      return
      end

I look into the code of navier5.f and find that avg1, avg2 returns the time average of vx and vx**2 in this example, so if we want to calculate the Reynolds stress, the formulation should be :
     R_uu=urms-uavg**2
also, if normalized,
     R_uu=(urms-uavg**2)/u_tau**2
so are the other components
I think this is the right way instead of that in the example









At 2014-05-29 06:54:10, nek5000-users at lists.mcs.anl.gov wrote:
>
>Hi Xianbei,
>
>You can define more objects by picking your
>own surface discriminators (which will 
>be dependent on your particular application).
>
>If you need more than 4 objects, you can increase maxobj
>in SIZE.
>
>Typically, dragx(0) would be the sum of the drag on all
>the objects --- this is sometimes useful.
>
>If you have (say) a 4-sided channel, you can just define
>one object instead of 4, and set it to be the collection
>of all faces for which cbc(f,e,1)='W  '.
>
>Paul
>
>
>
>On Thu, 29 May 2014, nek5000-users at lists.mcs.anl.gov wrote:
>
>> Hi,all:
>     I read Paul's reply in this mail https://lists.mcs.anl.gov/mailman/htdig/nek5000-users/2014-May/002756.html
>     And got to know something about objects. As you say, the object is a collection of faces, here , I have a question about the iboj, how could I know which face it represent with iboj? As in turbChannel.usr,
>      dragx_avg = alpha*dragx_avg + beta*0.5*(dragx(1)+dragx(2))
>      I think it means to do the time-average of the average dragx on planey1 and planey2, as
>      0.5*(dragx(1)+dragx(2))
>  indicates, while how to calculate the dragx on the planex1, planex2, planez1 and planez2? for the maxobj in this example is only 4 and dragx is defined as
>      dragx(0:maxobj)
>
>Thank you
>Xianbei
>_______________________________________________
>Nek5000-users mailing list
>Nek5000-users at lists.mcs.anl.gov
>https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20140529/677884d7/attachment-0001.html>