c----------------------------------------------------------------------- c c c OSCILLATING CYLINDER EXAMPLE c ============================ c c c This case illustrates the Nek5000 ALE formulation for c flow past an oscillating cylinder in 2D. c c (Note: The case "ocyl" is the first of two examples to c be posted. An improved version with faster computation c of mesh velocity will be posted soon.) c c Frequency and amplitude of cylinder motion are set in the c .rea file: c c 0.200000 p033 user: frequency of oscillation c 0.500000 p034 user: amplitude of oscillation c c Note that the amplitude cannot be too large or the c mesh will become entangled and the simulation will stop. c c c The ALE formulation in this case solves an elasticity c problem at every timestep in order compute the mesh c velocity. While robust, this approach is slow and not c recommended for production runs. A better approach is c to compute the mesh velocity for a single periodic cycle c and to then use a short Fourier exapansion based on a c small set of time points over the full cycle to reconstruct c the desired mesh velocity at any future time point. c Such an approach would avoid repetitious computation of c the mesh velocity, which for this example is far more expensive c than computing the velocity/pressure for the fluid. There c would be no loss in accuracy because the ALE formulation c allows for arbitrary mesh velocities, provided that they c meet the kinematic condition that the normal mesh velocity c component coincides with the normal fluid c velocity component on the boundary. c c c Prescribed boundary conditions are given in Cartesian coordinates. c The prescribed surface motion is indicated by the "mv " bc on your c moving surface (i.e., the cylinder surface). c c To use mesh motion, you must also have c c T ifstrs c T ifmvbd c c in the .rea file and SIZE must have the following settings: c c parameter (lx2=lx1-2) c parameter (ly2=ly1-2) c parameter (lz2=lz1 ) c parameter (lx1m=lx1,ly1m=ly1,lz1m=lz1) c c c An example of the result is shown in ocyl.gif. c c c c c----------------------------------------------------------------------- subroutine uservp (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,eg,f udiff = 0 utrans = 0 return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,eg,f ffx = 0. ffy = 0. ffz = 0. return end c----------------------------------------------------------------------- subroutine userq (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,eg,f qvol = 0.0 source = 0.0 return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,f,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,eg,f real ty common /c_mybc/ umeshx(lx1,ly1,lz1,lelt) $ , umeshy(lx1,ly1,lz1,lelt) $ , umeshz(lx1,ly1,lz1,lelt) if (x .lt. 0.2) then ! At inflow ux=1 uy=0. uz=0. else e = gllel(eg) ux=umeshx(ix,iy,iz,e) uy=umeshy(ix,iy,iz,e) uz=umeshz(ix,iy,iz,e) endif return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,eg,f real ty ty = (y+12.0)/24 ux=1.0*ty/100 uy=0.0 uz=0.0 return end c----------------------------------------------------------------------- subroutine usrdat include 'SIZE' include 'TOTAL' if (param(33) .eq. 0.0) param(33) = 1.0 if (param(34) .eq. 0.0) param(34) = 1.0e-2 return end c----------------------------------------------------------------------- subroutine usrdat3 return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' param(59) = 1 ! all elements deformed param(66) = 4 param(67) = 4 ifxyo = .true. ifusermv = .false. ! do not define our own mesh velocity return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' include 'RESTART' real x0(3) save x0 data x0 /3*0/ parameter (lt=lx1*ly1*lz1*lelv) common /scrns/ vort(lt,3), w1(lt), w2(lt) n = nx1*ny1*nz1*nelv call comp_vort3(vort , w1, w2, vx, vy, vz) call copy (T,vort,n) ! Vorticity --> T ifto = .true. ! Dump vorticity as T if (istep.eq.0) call set_obj ! define objects for surface integrals scale = 2. ! Cd = F/(.5 rho U^2 ) = 2*F c if (mod(istep,10).eq.0) c call torque_calc(scale,x0,.true.,.true.) ifusermv = .true. if (ifusermv) call my_meshv ! Compute our own mesh velocity return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine my_meshv include 'SIZE' include 'TOTAL' parameter (lt=lx1*ly1*lz1*lelv) common /cmeshv/ basev(lt) common /scrns/ h1(lt),h2(lt),rhs(lt),msk(lt),tmp(lt) real msk real torq,alpha_i,alpha_dot_i,alpha_f,alpha_dot_f real resultant(nx1*ny1*nz1*nelv),costheta(nx1*ny1*nz1*nelv) real sintheta(nx1*ny1*nz1*nelv),vel(nx1*ny1*nz1*nelv) real ucx(nx1*ny1*nz1*nelv),ucy(nx1*ny1*nz1*nelv) real ucz(nx1*ny1*nz1*nelv),xx(nx1*ny1*nz1*nelv) real yy(nx1*ny1*nz1*nelv) common /ctorq/ dragx(0:maxobj),dragpx(0:maxobj),dragvx(0:maxobj) $ , dragy(0:maxobj),dragpy(0:maxobj),dragvy(0:maxobj) $ , dragz(0:maxobj),dragpz(0:maxobj),dragvz(0:maxobj) c $ , torqx(0:maxobj),torqpx(0:maxobj),torqvx(0:maxobj) $ , torqy(-1:maxobj),torqpy(0:maxobj),torqvy(0:maxobj) $ , torqz(0:maxobj),torqpz(0:maxobj),torqvz(0:maxobj) c $ , dpdx_mean,dpdy_mean,dpdz_mean $ , dgtq(3,4) common /c_mybc/ umeshx(lx1,ly1,lz1,lelt) $ , umeshy(lx1,ly1,lz1,lelt) $ , umeshz(lx1,ly1,lz1,lelt) if (istep.eq.0) call my_base_meshv(basev,h1,h2,rhs,msk,tmp) aa = 0.01 n = nx1*ny1*nz1*nelv xo = 0.44 zo = 0 yo = 0.53 c both method are working for ucx and ucy do i = 1,n xx(i) =xm1(i,1,1,1)-xo yy(i) = ym1(i,1,1,1)-yo resultant(i) = sqrt( ( (xm1(i,1,1,1) - xo) $ * (xm1(i,1,1,1) - xo) ) $ + ( (ym1(i,1,1,1) - yo) * (ym1(i,1,1,1) - yo) ) ) costheta(i) = xx(i) / resultant(i) sintheta(i) = yy(i) / resultant(i) vel(i) = aa * resultant(i) ucy(i) = - vel(i) * costheta(i) ucz(i) = 0.0 ucx(i) = vel(i) * sintheta(i) c ucy(i) = - xx(i) c ucx(i) = yy(i) enddo c enddo n = nx1*ny1*nz1*nelv do i=1,n ! Translational velocity wx(i,1,1,1) = basev(i)*ucx(i) ! component. wy(i,1,1,1) = basev(i)*ucy(i) wz(i,1,1,1) = basev(i)*ucz(i) umeshx(i,1,1,1) = wx(i,1,1,1) umeshy(i,1,1,1) = wy(i,1,1,1) umeshz(i,1,1,1) = wz(i,1,1,1) enddo nface = 4 n = nx1*ny1*nz1*nelv do e=1,n do f=1,nface if (cbc(f,e,1).eq.'mv ') then open (unit=10,file="rectangle.dat") write (10,7020) e,f,xm1(1,1,f,e),ym1(1,1,f,e) 7020 format(7f20.8) endif enddo enddo return end c----------------------------------------------------------------------- subroutine my_base_meshv(basev,h1,h2,rhs,msk,tmp) include 'SIZE' include 'TOTAL' real basev(lx1,ly1,lz1,lelt),tmp(lx1,ly1,lz1,lelt) real h1(1),h2(1),rhs(1),msk(1) real m1,denom integer e,f n = nx1*ny1*nz1*nelv call rone (h1 ,n) call rzero(h2 ,n) call rone (msk,n) call rzero(tmp,n) delta = 0.5 do i=1,n rr = h1(i) arg = -rr/(delta**2) h1(i) = 1. + 9.0*exp(arg) enddo m1 = -1. z0 = 0. nface = 2*ndim do e=1,nelv do f=1,nface c Set Dirichlet for mesh velocity on all boundaries if (cbc(f,e,1).ne.'E ') call facev(msk,e,f,z0,nx1,ny1,nz1) c Inhomogeneous Dirichlet on cylinder only if (cbc(f,e,1).eq.'mv ') then call fcaver(xavg,xm1,e,f) if (xavg.gt.-1.2) call facev(tmp,e,f,m1,nx1,ny1,nz1) endif enddo enddo call axhelm (rhs,tmp,h1,h2,1,1) tol = 1.e-4 imsh = 1 ifield = 1 call hmholtz('mshv',basev,rhs,h1,h2,msk,vmult,imsh,tol,200,1) call sub2(basev,tmp,n) c call outpost(basev,basev,basev,basev,basev,' ') c stop return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine set_obj ! define objects for surface integrals c include 'SIZE' include 'TOTAL' integer e,f,eg 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 c write(6,*) 'Number of history=', nhis nhis = nhis + nobj c write(6,*) 'Number of history+obj=', nhis if (maxobj.lt.nobj) call exitti('increase maxobj in SIZE$',nobj) nxyz = nx1*ny1*nz1 nface = 2*ndim c write(6,*) 'Number of elements',nelv c write(6,*) 'Number of faces', nface do e=1,nelv do f=1,nface if (cbc(f,e,1).eq.'mv ') then iobj = 1 if (iobj.gt.0) then c write(6,*) 'DEBUG iobj nmember(iobj)', iobj, nmember(iobj) nmember(iobj) = nmember(iobj) + 1 c write(6,*) 'DEBUG nmember(iobj)', nmember(iobj) mem = nmember(iobj) c write(6,*) 'DEBUG e,f,mem, iobj, mem', e,f,mem, iobj eg = lglel(e) object(iobj,mem,1) = eg object(iobj,mem,2) = f c write(6,1) iobj,mem,f,eg,e,nid,' OBJ' c 1 format(6i9,a4) endif endif enddo enddo c write(6,*) 'number',(nmember(k),k=1,4) c return end c-----------------------------------------------------------------------