<div dir="ltr">Dear Neks, <div><br></div><div>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:</div><div><br></div><div><div>      scale = 2.  ! Cd = F/(.5 rho U^2 ) = 2*F</div><div>      if (mod(istep,10).eq.0) call torque_calc(scale,x0,.true.,.false.)</div></div><div><br></div><div>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".</div><div><br></div><div>What I've been doing is something similar to the following:</div><div><br></div><div><a href="https://lists.mcs.anl.gov/mailman/htdig/nek5000-users/2011-September/001473.html">https://lists.mcs.anl.gov/mailman/htdig/nek5000-users/2011-September/001473.html</a><br></div><div><br></div><div>but it can't get the right Cd value.</div><div><br></div><div><br></div><div><br></div><div>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.</div><div><br></div><div>My set object subroutine is:</div><div><br></div><div><div>      subroutine set_obj  ! define objects for surface integrals</div><div>c</div><div>      include 'SIZE'</div><div>      include 'TOTAL'</div><div><br></div><div>      integer e,f,eg</div><div>      real xx,yy,nx,ny,rn</div><div><br></div><div>      nobj = 1</div><div>      iobj = 0</div><div>      do ii=nhis+1,nhis+nobj</div><div>         iobj = iobj+1</div><div>         hcode(10,ii) = 'I'</div><div>         hcode( 1,ii) = 'F'</div><div>         hcode( 2,ii) = 'F'</div><div>         hcode( 3,ii) = 'F'</div><div>         lochis(1,ii) = iobj</div><div>      enddo</div><div>      nhis = nhis + nobj</div><div><br></div><div>      if (maxobj.lt.nobj) call exitti('increase maxobj in SIZE$',nobj)</div><div><br></div><div>      nxyz  = nx1*ny1*nz1</div><div>      nface = 2*ndim</div><div><br></div><div>      do e=1,nelv</div><div>       do ix=1,nx1</div><div>        do iy=1,ny1</div><div>        do iz=1,nz1</div><div>         xx=xm1(ix,iy,iz,e)</div><div>         yy=ym1(ix,iy,iz,e)</div><div>         if ((sqrt(xx*xx+yy*yy)).lt.(0.51)) then<span class="gmail-Apple-tab-span" style="white-space:pre">              </span></div><div>          do f=1,nface</div><div>           if (cbc(f,e,1).eq.'W  ') then</div><div>            iobj  = 1</div><div>            if (iobj.gt.0) then</div><div>             nmember(iobj) = nmember(iobj) + 1</div><div>             mem = nmember(iobj)</div><div>             eg  = lglel(e)</div><div>             object(iobj,mem,1) = eg</div><div>             object(iobj,mem,2) = f</div><div>            endif</div><div>           endif</div><div>          enddo</div><div>         endif</div><div>        enddo</div><div>        enddo</div><div>       enddo</div><div>      enddo</div><div><br></div><div><br></div><div>c     write(6,*) 'number',(nmember(k),k=1,4)</div><div>c</div><div>      return</div><div>      end</div></div><div><br></div><div><br></div><div><br></div><div>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.</div><div><br></div><div>What am I doing wrong?</div><div><br></div><div>Thanks in advace,</div><div><br></div><div>Juan Pablo.</div></div>