% PUCKS_ON_ICE.M % % Pucks_on_ice examines the effect of an impulsively-imposed body % force upon a series of "hockey pucks" of uniform mass that are % distributed about the surface of a frozen ocean. The pucks are % initially at rest. The body force, which pushes the pucks to the % northeast, acts for 3 days and then stops. The integration of % the ensuing motion continues for period of time after the forcing % has ceased. % % This is an "f-plane" calculation in two dimensions, in which the % Coriolis force is set to a fixed latitude (30 degrees north latitude % in this example). There are a number of "forces" in this problem, % in addition to the body force. First, Rayleigh friction is weak and % opposed to the motion of the particles. If the particle velocity is u, % then the frictional force per unit mass is expressed as -r*u. Second, % because the frozen ocean surface is not flat, gravity (g, directed % downwards) accelerates the pucks down the sloping surface. In % the x-direction this acceleration, or force per unit mass, is given % by -g*(dh/dx), where dh is the change in surface height over a % distance, dx. In this example the inertial period, 2\pi/f, is 24 hours. % The equations used are as follows: % % dx/dt = u % dy/dt = v % du/dt - fv = -ru - g(dh/dx) + Fx % dv/dt + fu = -rv - g(dh/dy) + Fy, % % where (x,y) are the eastward and northward particle position, (u,v) % are the eastward and northward particle velocity, f is the Coriolis % parameter, and (Fx,Fy) are the eastward and northward components of % the external body force. % % The frozen sea surface has a "bump" that is Gaussian in shape, % centered at (x_m,y_m), with an amplitude, h_0, and a length scale % defined by (lx,ly). The sea surface height is given by the following: % % h(x,y) = h_0*exp(-((x-x_m)/lx)^2 -((y-y_m)/ly)^2), % where h_0 = 5 cm, x_m = y_m = 400 km, and lx = ly = 30 km. % % This example illustrates several important phenomena: % (1) Particles move at right angles to the surface force (Ekman layers). % (2) Particles oscillate with an inertial period. % (3) Particles near the sea surface bump get "trapped" and continue to % move long after those far away from the bump have come to rest. % This illustrates the "geostrophic" balance. % % All units used in the calculation are given in MKS, but are % converted to km and hours for display. In Figure 1 the trajectories % of 30 particles are plotted with an 'x' every 24 hours. In Figure 2 % the trajectories of two particular particles are shown (A & B) relative % to where they started out. % % Terry Joyce, July 2002, modified for low friction limit Sept. 2003 %---------------------------------------------------------------------------- % DISPLAY TEXT clear close all display([' Pucks on ice']); str2(1) = {'Pucks on ice, Terry Joyce, July 2002:' }; str2(2) = {' '}; str2(3) = {'Pucks\_on\_ice examines the effect of an impulsively-imposed body'}; str2(4) = {'force upon a series of "hockey pucks" of uniform mass that are'}; str2(5) = {'distributed about the surface of a frozen ocean. The pucks are'}; str2(6) = {'initially at rest. The body force, which pushes the pucks to the'}; str2(7) = {'northeast, acts for 3 days and then stops. The integration of'}; str2(8) = {'the ensuing motion continues for period of time after the forcing'}; str2(9) = {'has ceased.'}; str2(10) = {''}; str2(11) = {'This is an "f-plane" calculation in two dimensions, in which the'}; str2(12) = {'Coriolis force is set to a fixed latitude (30 degrees north latitude'}; str2(13) = {'in this example). There are a number of "forces" in this problem,'}; str2(14) = {'in addition to the body force. First, Rayleigh friction is weak and'}; str2(15) = {'opposed to the motion of the particles. If the particle velocity is u,'}; str2(16) = {'then the frictional force per unit mass is expressed as -r*u. Second,'}; str2(17) = {'because the frozen ocean surface is not flat, gravity (g, directed'}; str2(18) = {'downwards) accelerates the pucks down the sloping surface. In'}; str2(19) = {'the x-direction this acceleration, or force per unit mass, is given'}; str2(20) = {'by -g*(dh/dx), where dh is the change in surface height over a'}; str2(21) = {'distance, dx. In this example the inertial period, 2\pi/f, is 24 hours.'}; str2(22) = {'The equations used are as follows:'}; str2(23) = {''}; str2(24) = {'dx/dt = u'}; str2(25) = {'dy/dt = v'}; str2(26) = {'du/dt - fv = -ru - g(dh/dx) + Fx'}; str2(27) = {'dv/dt + fu = -rv - g(dh/dy) + Fy,'}; str2(28) = {''}; str2(29) = {'where (x,y) are the eastward and northward particle position, (u,v)'}; str2(30) = {'are the eastward and northward particle velocity, f is the Coriolis'}; str2(31) = {'parameter, and (Fx,Fy) are the eastward and northward components of'}; str2(32) = {'the external body force.'}; str2(33) = {''}; str2(34) = {'The frozen sea surface has a "bump" that is Gaussian in shape, '}; str2(35) = {'centered at (x_m,y_m), with an amplitude, h_0, and a length scale'}; str2(36) = {'defined by (lx,ly). The sea surface height is given by the following:'}; str2(37) = {''}; str2(38) = {'h(x,y) = h_0*exp(-((x-x_m)/lx)^2 -((y-y_m)/ly)^2),'}; str2(39) = {'where h_0 = 5 cm, x_m = y_m = 400 km, and lx = ly = 30 km.'}; str2(40) = {''}; str2(41) = {'This example illustrates several important phenomena:'}; str2(42) = {'(1) Particles move at right angles to the surface force (Ekman layers).'}; str2(43) = {'(2) Particles oscillate with an inertial period.'}; str2(44) = {'(3) Particles near the sea surface bump get "trapped" and continue to'}; str2(45) = {' move long after those far away from the bump have come to rest.'}; str2(46) = {' This illustrates the "geostrophic" balance.'}; str2(47) = {''}; str2(48) = {'*** Hit any key to begin...& be patient! ***'}; hf3 = figure(5); clf set(hf3,'Position',[40 40 700 900]) set(gca,'Visible','off') text(-0.1, 0.50, str2,'FontSize',10,'HorizontalAlignment','left') pause %---------------------------------------------------------------------------- % START %---------------------------------------------------------------------------- % Define parameters count=0; % number of particles t1=0:1000; % time in hours t=3600*t1; % time in seconds rfact=.05; % set the coefficient for Rayleigh friction %-----------------> here are some key paprameters that can be changed h0 = .05; % set topography here lat = 30; % set latitude here %---------------------------------------------------------------------------- % Set initial positions (x0,y0) in [m] and velocities (u0,v0) in [m/s] of 30 % particles, then integrate the equations governing their motion (using ode45), % and plot each particle trajectory more off figure(1), clf for i=1:5 for j=1:6 count =count+1; display([' calculating trajectory for particle # ' num2str(count),' of 30']); x0(i,j)=(300+50*(i-1))*10^3; u0(i,j)=0; y0(i,j)=(300+40*(j-1))*10^3; v0(i,j)=0; % The function 'ode45' integrates the differential equations in 'deriv', given the % integration time, t, and the initial position (x0,y0) and velocity (u0,v0) of % the particle, and the topography amplitude, h0 and the reference latitude. % The last two don't change, but putting them here passes them to the function. [T,XY]=ode45('deriv',t,[x0(i,j),y0(i,j),u0(i,j),v0(i,j),h0,lat]); % Plot particle trajectory over time t plot(XY(:,1)/10^3,XY(:,2)/10^3); hold on h=plot(XY(1,1)/10^3,XY(1,2)/10^3,'ro'); set(h,'linewidth',2) figure(1) plot(XY(1:24:length(T),1)/10^3,XY(1:24:length(T),2)/10^3,'bx') axis([250 550 250 550]) axis square % Save trajectories of 2 particles for plotting in Figure 2 if(i==1 & j==1) XYA=XY; end if(i==3 & j==3) XYB=XY; end end end % Put labels on the plot xlabel('East distance [km]') ylabel('North distance [km]') text(x0(1,1)/10^3+5,y0(1,1)/10^3+5,'A') text(x0(3,3)/10^3+5,y0(3,3)/10^3+5,'B') % plot force vector quiver(260,510,35,35,0,'k') text(255,502,'F') title('Pucks on ice, Fig. 1', 'fontsize', 14) %---------------------------------------------------------------------------- % Plot close-up of particles A and B in Figure 2 figure(2), clf axis([-70 70 -70 70]) xlabel('East distance from initial position [km]') ylabel('North distance from initial position [km]') axis square hold on title('Pucks on ice, Fig. 2', 'fontsize', 14) for j=1:floor(length(T)/24) plot((XYA(24*j-23,1)-XYA(1,1))/10^3,(XYA(24*j-23,2)-XYA(1,2))/10^3,'bx') plot((XYB(24*j-23,1)-XYB(1,1))/10^3,(XYB(24*j-23,2)-XYB(1,2))/10^3,'rx') plot((XYA(24*j-23:24*j,1)-XYA(1,1))/10^3,(XYA(24*j-23:24*j,2)-XYA(1,2))/10^3,'b-'); plot((XYB(24*j-23:24*j,1)-XYB(1,1))/10^3,(XYB(24*j-23:24*j,2)-XYB(1,2))/10^3,'r-'); pause(0.2) end legend('particle A','particle B',4) quiver(0,0,20,20,0,'k') text(22,22,'F') %---------------------------------------------------------------------------- % Define and plot sea surface height/topography in the domain % Define the topography of the sea surface to be a symmetric Gaussian "bump" % All quantities have units of [m] xtop=min(x0(:,1)-10^5):5*10^3:max(x0(:,1)+10^5); ytop=xtop; [xx,yy]=meshgrid(xtop,ytop); % define x,y points in the domain lx=30*10^3; % x length scale of bump ly=30*10^3; % y length scale of bump xm=400*10^3; % x location of the maximum height of bump ym=400*10^3; % y location of the maximum height of bump bump=h0*exp(-((xx-xm)/lx).^2 - ((yy-ym)/ly).^2); hpuck=h0*exp(-((x0-xm)/lx).^2 - ((y0-ym)/ly).^2); figure(3), clf colormap(cool) surf(xx/10^3,yy/10^3,bump) shading interp hold on if(h0>0) zmin=0; zmax=2*h0; else zmin=2*h0; zmax=0; end plot3(x0/10^3,y0/10^3,hpuck,'k*') set(gca,'fontsize',14) axis([200 600 200 600 zmin-abs(h0/10) zmax+abs(h0/10)]) view(-40,20); %label particles A, B again za=h0*exp(-((x0(1,1)-xm)/lx).^2 - ((y0(1,1)-ym)/ly).^2)+abs(h0/10); zb=h0*exp(-((x0(3,3)-xm)/lx).^2 - ((y0(3,3)-ym)/ly).^2)+h0/10; text(x0(1,1)/10^3+10,y0(1,1)/10^3+5,za,'A','fontsize',14) text(x0(3,3)/10^3-10,y0(3,3)/10^3-5,zb,'B','fontsize',14) %other stuff xlabel('East distance [km]') ylabel('North distance [km]') zlabel('Height [m]') title('Pucks on ice, Fig. 3') %---------------------------------------------------------------------------- % Text explaining figures str1(1) = {'Pucks\_on\_ice, discussion & homework problem:' }; str1(2) = {''}; str1(3) = {'Figure 1 shows the trajectories of 30 particles that were initially'}; str1(4) = {'at rest. The "wind" body force, F, is applied for 3 days in the '}; str1(5) = {'direction indicated by the vector. Red circles indicate the initial'}; str1(6) = {'position of each particle, and x''s mark the particle position after each'}; str1(7) = {'day. Many of the particles stop their translation at right angles to the'}; str1(8) = {'wind after the forcing ceases, and their motion decays in inertial circles'}; str1(9) = {'at a rate given by the friction. Some of the particles continue to move'}; str1(10) = {'in a clockwise fashion around the topographic bump located at'}; str1(11) = {'(400,400) km in the plot. '}; str1(12) = {''}; str1(13) = {'Use the "zoom" feature to observe the trajectories of different particles'}; str1(14) = {'both during the forcing and after the forcing has ceased.'}; str1(15) = {''}; str1(16) = {'The trajectories of the particles labelled "A" and "B" are plotted in Figure 2.'}; str1(17) = {'To see these plotted again, run draw\_fig2.m.'}; str1(18) = {'Compare the difference in the motion of the two particles relative to their'}; str1(19) = {'initial positions. Particle A is far from the topography and does not feel any'}; str1(20) = {'"force" down the slope. Particle B, however, feels the slope, begins to'}; str1(21) = {'slide down it and is deflected to the right by the Coriolis force. It'}; str1(22) = {'continues to move long after the forcing stops because the sloping surface'}; str1(23) = {'remains. Except for friction (which causes the outward spiral), it is in'}; str1(24) = {'"geostrophic balance".'}; str1(25) = {''}; str1(26) = {'Figure 3 illustrates the shape of the frozen sea surface, and the initial'}; str1(27) = {'positions of the particles.'}; str1(28) = {''}; str1(29) = {'Adjust the following parameters and observe how the particle motion changes:'}; str1(30) = {'(1) The latitude, to observe how the period of the inertial oscillations changes'}; str1(31) = {'(2) The sign of h0 in "pucks\_on\_ice.m", to change the "bump" into a "dimple"'}; str1(32) = {'(3) The hemisphere, by choosing a negative latitude in "pucks\_on\_ice.m"'}; str1(33) = {''}; str1(34) = {'For each case explain what happens & show plots of the results.'}; str1(35) = {''}; str1(36) = {'Terry Joyce, 12.808 study problem, July 2002'}; hf3=figure(6); clf set(hf3,'Position',[40 40 700 700]) set(gca,'Visible','off') text(-0.1, 0.50, str1,'FontSize',10,'HorizontalAlignment','left')