* Program to simulate "standing" wave response to obstacle at * bottom of steady-flowing stream of water in trough. * Formulas are based on Lighthill's "Waves in Fluids", chap. 3 * In this program the obstacle is a "cosine bump" of length L, i.e., * half of a single period of a cosine wave of wavelength 2*L, * centered on x=0. The maximum height of the bump is a. * The depth z of the bottom is therefore * | -h if |x| > L/2 * z = | * | a*cos(pi*x/L) if |x| <= L/2 program stand1 implicit none integer iv,ivlo,ivhi,niter,icol(50) double precision h,l,a,xmin,xmax,zmin,zmax,g,pi,vmax, + vtry(50),v,f,alpha,diff,k0,kh,kl,zeta, + x1,y1,y2,y3,y4,yoff,xv,yv,gamma,factor * Define nominal depth of trough (in feet) h = 2.d0 * Define length of "cosine bump" (in feet) l = 0.5d0 * Define maximum height of "cosine bump" (in feet) a = 0.4d0 * Define distances upstream and downstream to plot xmin = -2.d0 xmax = +6.d0 * Define minimum and maximum z values of plot zmin = - h * 1.2d0 zmax = + h * 1.2d0 * Define constants g = 32.d0 pi = acos(-1.d0) * Find maximum allowed velocity for which standing waves * can be created vmax = sqrt(g*h) * Define the stream velocity values to try (as fractions of vmax) ivlo = 1 ivhi = 5 vtry(1) = 0.30d0 vtry(2) = 0.50d0 vtry(3) = 0.75d0 vtry(4) = 0.90d0 vtry(5) = 0.95d0 * vtry(6) = 0.99d0 * Select colors for plotting icol(1) = 1 icol(2) = 2 icol(3) = 3 icol(4) = 4 icol(5) = 6 * icol(6) = 6 * Initialize plot file x1 = xmax * 0.4d0 y1 = zmax * 0.9d0 yoff = zmax*0.15d0 y2 = y1 - yoff y3 = y1 - 2.d0*yoff y4 = y1 - 3.d0*yoff open(1,file='stand1.kumac',status='new') write(1,5) xmin,xmax,zmin,zmax, + pi,a,l,h, + xmax,xmin, + -l/2.d0,l/2.d0, + xmin,-h,-l/2.d0,-h, + l/2.d0,-h,xmax,-h, + x1,y1,h, + x1,y2,vmax, + x1,y3,l, + x1,y4,a 5 format('fort/file 66 stand1.ps'/ + 'meta 66 -111'/ + 'zone 1 1'/ + 'null',4(1x,f10.6)/ + 'igset plci 1'/ + 'igset txci 1'/ + 'igset chhe 0.35'/ + 'nxd = 100'/ + 'nxu = 100'/ + 'nxb = 50'/ + 'pi = ',f10.8/ + 'a = ',f10.6/ + 'l = ',f10.6/ + 'h = ',f10.6/ + 'sigma xd = array([nxd],0#',f10.6,')'/ + 'sigma xu = array([nxu],',f10.6,'#0)'/ + 'sigma xb = array([nxb],',f10.6,'#',f10.6,')'/ + 'sigma bump = -[h] + [a]*cos([pi]*xb/[l])'/ + 'graph [nxb] xb bump l'/ + 2('line',4(1x,f10.6)/)/ + 'itx ',f10.6,1x,f10.6,' Depth = ',f4.2,' ft'/ + 'itx ',f10.6,1x,f10.6,' Max V = ',f5.2,' ft/s ("5#)'/ + 'itx ',f10.6,1x,f10.6,' Bump length = ',f4.2,' ft'/ + 'itx ',f10.6,1x,f10.6,' Bump height = ',f4.2,' ft') * Loop over trial values of the stream velocity do 100 iv = ivlo, ivhi v = vtry(iv)*vmax * Determine the gravity wave number k0 consistent with the depth * and current trial value of v * * Must numerically solve v**2/gh = tanh(k0*h)/(k0*h) * * or f = tanh(alpha)/alpha where f = vtry(iv)**2 f = vtry(iv)**2 alpha = 1.d0/f niter = 0 10 niter = niter + 1 diff = f*alpha - tanh(alpha) if(abs(diff).lt.1.d-10)then k0 = alpha/h write(6,15) v,vmax,f,alpha,niter 15 format(' For (v=',f10.6,'/vmax=',f10.6,')**2=',f8.6, + ', k0*h=',f10.6,'(niter=',i2,')') goto 20 else alpha = alpha - diff/(f-1./cosh(alpha)**2) endif if(niter.gt.20)then write(6,*)'Error: Failed to converge on alpha' write(6,*)' Trying to solve for f=',f stop endif goto 10 * Compute amplitude of downstream standing wave 20 kl = k0*l kh = k0*h zeta = -16.d0*pi*a * kl*cos(kl/2.d0)/(pi**2-kl**2) + *sinh(kh)/(sinh(2.d0*kh)-2.d0*kh) if(abs(zeta).gt.h)then write(6,*)'Warning: Zeta=',zeta,' skipping this V' goto 100 endif write(6,22) kl,zeta 22 format(' -> k0*L=',e14.6,' zeta=',e14.6) * Write downstream sine-wave info to plot file xv = xmin * 0.5d0 yv = y1 - yoff*(iv-1) write(1,25) zeta,k0,icol(iv),icol(iv),xv,yv,vtry(iv) 25 format('zeta = ',e14.6/ + 'k0 = ',e14.6/ + 'sigma zd = [zeta]*sin([k0]*xd)'/ + 'igset plci ',i1/ + 'graph [nxd] xd zd l'/ + 'igset txci ',i1/ + 'itx ',f10.6,1x,f10.6,' V/V?MAX! = ',f4.2) * Determine the tension wave number k0 consistent with the * trial value of v, using the approximate formula: * v = sqrt(k0*T/rho) * where T is the surface tension and rho is water density: * gamma = T/rho = 0.0026 ft**3/s**2 gamma = 2.6d-3 k0 = v**2/gamma * Compute amplitude of upstream standing wave kl = k0*l kh = k0*h if(kh.lt.50.d0)then factor = sinh(kh)/(sinh(2.d0*kh)-2.d0*kh) else factor = 1.d0/(1.d0-2.d0*kh) endif zeta = +16.d0*pi*a * kl*cos(kl/2.d0)/(pi**2-kl**2) + *factor write(6,22) kl,zeta * Write upstream sine-wave info to plot file write(1,30) zeta,k0 30 format('zeta = ',e14.6/ + 'k0 = ',e14.6/ + 'sigma zu = [zeta]*sin([k0]*xu)'/ + 'graph [nxu] xu zu l') 100 continue stop end