% gwavemovie; make a movie of wave fluid velocity % A pmovie of a progating gravity wave. wavelength = % 100 m; you must input a water depth. The surface % displacement is shown as a red line; its amplitude % is not shown to the scale of the water depth or of the % particle velocity (i.e, amplitude is arbitrary). % % Jim Price, January 2001 % clear clear memory close all path(path, 'G:\educ\rotation') % Set some default graphics things. set(0,'DefaultTextFontSize',14) set(0,'DefaultAxesFontSize',14) set(0,'DefaultAxesLineWidth',1.6) set(0,'DefaultLineLineWidth',1.4) g = 9.8; denw = 1025.; sig = 7.3e-2; % surface tension, N/m a = 3.; lambda = 100. ; % wavelength, m k = 2*pi./lambda; H = input(' Enter the water depth (10 - 5000 m) '); omega = sqrt((g*k + (sig/denw)*(k^3))*tanh(k*H)); T = 2*pi/omega; c = sqrt((g/k + (sig/denw)*k)*tanh(k*H)); x = 0:5:100; nx = max(size(x)); z = 0:-10:-70; nz = max(size(z)); nt = 125; xp = zeros(4,nt+1); zp = xp; xp(:,1) = [40. 40. 40. 40.]; zp(:,1) = [-5. -15. -25. -35.]; dt = 0.3; for jt = 1:nt time = jt*dt; pht = omega*time; eta = a*cos(x*k - pht); for jx=1:nx for jz=1:nz u(jz,jx) = a*omega*(cosh(k*(z(jz) + H))/sinh(k*H))*cos(k*x(jx)-pht); w(jz,jx) = a*omega*(sinh(k*(z(jz) + H))/sinh(k*H))*sin(k*x(jx)-pht); if -z(jz) > H u(jz,jx) = 0.; w(jz, jx) = 0.; end xq(jz,jx) = x(jx); zq(jz,jx) = z(jz); end end % track a particle for m=1:4 xph(m) = xp(m,jt) + 0.5*dt*a*omega*(cosh(k*(zp(m,jt) + H))/sinh(k*H))*cos(k*xp(m,jt)-pht); zph(m) = zp(m,jt) + 0.5*dt*a*omega*(sinh(k*(zp(m,jt) + H))/sinh(k*H))*sin(k*xp(m,jt)-pht); xp(m,jt+1) = xp(m,jt) + dt*a*omega*(cosh(k*(zph(m) + H))/sinh(k*H))*cos(k*xph(m)-(pht+dt/2)); zp(m,jt+1) = zp(m,jt) + dt*a*omega*(sinh(k*(zph(m) + H))/sinh(k*H))*sin(k*xph(m)-(pht+dt/2)); end figure(1) clf reset set(gcf,'position', [100 100 600 500]) quiver(xq, zq, u, w, 0) axis([0 80 -60 20]) ylabel('depth, m') xlabel('x distance, m') hold on plot(x, eta, 'r') hold on; plot([0 150], [0 0], ':') hstr = num2str(H); kh = floor(k*H*10)/10; khstr = num2str(kh); dep = ['\lambda = 100 m, h = ', hstr, ' m, kh = ', khstr ]; text(15, -55, dep) timtex = ['time = ', num2str(time), ' sec']; text(10, 10, timtex) title('Eulerian velocity and particle trajectories') for m=1:4 hold on plot([xp(m,jt+1)], [zp(m,jt+1)], 'g*') hold on plot(xp(m,1:jt), zp(m,1:jt), 'r') end M(jt) = getframe(gcf, [0 0 600 500]); end % print -depsc2 gwaveposter.eps % mpgwrite(M, colorbar, 'gwavemovie.mpg')