[Nek5000-users] Problem with set_object

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

Dear Neks,

I'm trying to calculate the drag on a cylinder surface in a flow. I started
from "ext_cyl", "TurbChannel" and some threads on this list to understand
how to set objects. I'm using torque_calc in user check:

      scale = 2.  ! Cd = F/(.5 rho U^2 ) = 2*F
      if (mod(istep,10).eq.0) call torque_calc(scale,x0,.true.,.false.)

The problem is that the cylinder is rotating, so its boundary conditions
are "v  ", and the values are set to be the tangential velocity that I
need, for ux and uy (2D case), so I can't follow exactly the procedure en
"ext_cyl".

What I've been doing is something similar to the following:

https://lists.mcs.anl.gov/mailman/htdig/nek5000-users/2011-September/001473.html

but it can't get the right Cd value.



I'm testing my .usr file with "exy_cyl" case, trying to obtain the same
result but it doesn't work. For example, with Re=20, I should get a Cd=2.0
aprox., but I'm getting 25.0 aprox.

My set object subroutine is:

      subroutine set_obj  ! define objects for surface integrals
c
      include 'SIZE'
      include 'TOTAL'

      integer e,f,eg
      real xx,yy,nx,ny,rn

      nobj = 1
      iobj = 0
      do ii=nhis+1,nhis+nobj
         iobj = iobj+1
         hcode(10,ii) = 'I'
         hcode( 1,ii) = 'F'
         hcode( 2,ii) = 'F'
         hcode( 3,ii) = 'F'
         lochis(1,ii) = iobj
      enddo
      nhis = nhis + nobj

      if (maxobj.lt.nobj) call exitti('increase maxobj in SIZE$',nobj)

      nxyz  = nx1*ny1*nz1
      nface = 2*ndim

      do e=1,nelv
       do ix=1,nx1
        do iy=1,ny1
        do iz=1,nz1
         xx=xm1(ix,iy,iz,e)
         yy=ym1(ix,iy,iz,e)
         if ((sqrt(xx*xx+yy*yy)).lt.(0.51)) then
          do f=1,nface
           if (cbc(f,e,1).eq.'W  ') then
            iobj  = 1
            if (iobj.gt.0) then
             nmember(iobj) = nmember(iobj) + 1
             mem = nmember(iobj)
             eg  = lglel(e)
             object(iobj,mem,1) = eg
             object(iobj,mem,2) = f
            endif
           endif
          enddo
         endif
        enddo
        enddo
       enddo
      enddo


c     write(6,*) 'number',(nmember(k),k=1,4)
c
      return
      end



So what I'm trying to do is to capture the cylinder (radius=0.5) on a 0.51
radius circle, and to identify the faces with "W  " B.C. I think it should
work and give me the same result as the one obtained with the original
ext_cyl.usr.

What am I doing wrong?

Thanks in advace,

Juan Pablo.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20170609/4fe3b392/attachment.html>