% Solve a 1-D Couette flow numerically. % The IC is a state of rest. The fluid is forced % by the BC at the top of the column (z=0), either % a constant velocity, or, % a constant stress, determined by stresBC. % The lower boundary condition is either % u = 0 (as if a fixed wall and no-slip), % or dudz = 0 (as if free slip). % % There has been no effort to make this efficient or % the least bit sophisticated. The plotting is mostly % done in the dimensional variables, mainly because % the BC options make it hard to use non-d variables. % % Jim Price, Oct 99; Modified in Dec 01. % % This is public domain. clear close all set(0,'DefaultLineLineWidth',1.2) set(0,'DefaultTextFontSize',13) set(0,'DefaultAxesLineWidth',1.1) set(0,'DefaultAxesFontSize',13) % the text can go to about here > X str2(1) = {'Couette:' }; str2(2) = {' '}; str2(3) = {'Solve for diffusive flow in a channel, sometimes '}; str2(4) = {'called Couette flow, by numerical methods. '}; str2(5) = {'The only input is required to specify the upper'}; str2(6) = {'and lower boundary conditions. The upper BC '}; str2(7) = {'drives the flow by either an imposed speed or'}; str2(8) = {'an imposed stress; the lower BC is either '}; str2(9) = {'no-slip or free-slip. Solutions are compared to'}; str2(10) = {'a similarity solution that assumes an infinite'}; str2(11) = {'domain, and various energy and vorticity '}; str2(12) = {'budgets are evaluated. The dimensional values'}; str2(13) = {'are consistent with a turbulent upper ocean.'}; str2(14) = {''}; str2(15) = {'This code is public domain.'}; str2(16) = {''}; str2(17) = {'Jim Price, March, 2000.'}; hf3 = figure(10); clf set(hf3,'Position',[100 100 600 600]) set(gca,'Visible','off') text(0, 0.50, str2,'FontSize',13,'HorizontalAlignment','left') % set some parameters % here we use dimensional variables, with scales % like those of a 'turbulent' upper ocean dz = 1.; % the grid interval, m L = 50. % the grid extent (length), m z = 0:-dz:-L ; % the grid defined s = 0*z; % the dependent variable, think of it as velocity [nn ns] = size(s) delsq = s; nstep = 3000; % the number of time steps dt = 25.; % the time step, sec K = 1.e-2; % the diffusivity (MKS, appropriate for the upper ocean) w = dt*K/dz^2 % must be less then 0.5 for numerical stability rho = 1025.; % nominal constant density of water % define the type of BC: Uo = 0.1; % an imposed speed upper BC tauBC = 0.1/1025; % imposed constant stress/density as upper BC stresBC = menu('Choose the upper BC:', 'Imposed speed', ... 'Imposed stress') stresBC = stresBC - 1; freeslip = menu('Choose the lower BC:', 'No slip', 'Free slip') freeslip = freeslip - 1; figure(1) clf reset km = 1:ns-2; kp = km + 2; kmid = 2:ns-1; % begin time-stepping; the 'scheme' used here is forward in time % and centered in space - very crude, and it works pretty well! for j=1:nstep % estimate the Laplacian delsq(kmid) = s(km) + s(kp) - 2.*s(kmid); delsq(1) = 0.; delsq(ns) = 0.; % step forward s = s + w*delsq; % now set the upper boundary condition if stresBC == 0 s(1) = Uo; else s(1) = s(2) + dz*tauBC/K; end % the lower boundary condition if freeslip == 0 s(ns) = 0.; else s(ns) = s(ns-1); end % compute and store some diagnostics time(j) = (j-1)*dt; topstress(j) = -K*(s(2) - s(1))/dz; botstress(j) = -K*(s(ns) - s(ns-1))/dz; trans(j) = (sum(s) - s(1))*dz; % do a volume integrals over these limits nv1 = round(20/dz); nv2 = round(30/dz); topwork(j) = -s(nv1)*K*(s(nv1+1) - s(nv1-1))/(2*dz); botwork(j) = s(nv2)*K*(s(nv2+1) - s(nv2-1))/(2*dz); ds = 0.; kes = 0.; for m=nv1:nv2 ds = ds - dz*K*(((s(m+1)-s(m))/dz)^2); kes = kes + 0.5*dz*(s(m)^2); end diss(j) = ds; ke(j) = kes; efold = 1 - exp(-1); ssum = cumsum((s+1.e-4)); ze(j) = -interp1(ssum,z,efold*ssum(ns), 'linear'); zK(j) = sqrt(0.5*K*time(j)); % vorticity budget over this volume avgvort(j) = s(nv1) - s(nv2); topvflux(j) = -K*delsq(nv1)/dz^2; botvflux(j) = -K*delsq(nv2)/dz^2; netvflux(j) = botvflux(j) - topvflux(j); % plot profiles every nplot steps nplot = 100; if mod(j,nplot) == nplot - 1 hold on znd = z/sqrt(K*time(j)); Und = Uo; if stresBC == 1 Und = tauBC*(tauBC/K); end figure(1) subplot(2,1,1) hold on plot(s, z) subplot(2,1,2) hold on if stresBC == 0 plot(s/Und, znd) else for m=1:ns-1 shear(m) = -(s(m+1) - s(m))/dz; zshear(m) = 0.5*(z(m+1) + z(m)); end plot(shear/(tauBC/K), zshear/sqrt(K*time(j))) end % some things done once to figure 1: if j <= 2*nplot % evaluate and plot the similarity solution, even if it % may not be appropriate! deta = 0.01; eta = 0:deta:5; e2 = exp(-0.25*(eta.^2)); integ = cumsum(e2)*deta; usim = (1 - integ/sqrt(pi)); plot(fliplr(usim), fliplr(-eta),'r.', 'LineWidth',1.6) axis([-0.2 1. -6 0]) xxx = plot([-99 -90], [-93 -90], 'b', [-99 -97], [-93 -99], 'r', 'LineWidth',2.0) legend(xxx, 'numerical soln', 'similarity soln',4) % add scales and labels subplot(2,1,1) xlabel('U, m/s') ylabel('Z, m') delt = nplot*dt/(L^2/K); tit = strcat('Couette flow, numerically, dt = ',num2str(delt), ... ' (non-d)') title(tit) subplot(2,1,2) if stresBC == 0 xlabel('current, U/U_o') else xlabel('shear, \partialU/\partialz/(U_{*}^{2}/K)') end ylabel('Z/sqrt(K*t)') end % if on things done once to fig 1 drawnow % flush the print buffer end % if on whether to plot this time step end % loop on j for the time step % time-integrate the top and bottom stress to compare with the % transport totstr = cumsum(topstress-botstress)*dt; figure(4) clf reset subplot(3,1,1) plot(time, topstress, time, botstress,':') %axis([ 0 4e6 0 2e-5]) ylabel('kin stress, m*m/s*s') legend('top stress', 'bottom stress') title('overall momentum budget') subplot(3,1,2) plot(time, totstr, time, trans, ':') ylabel('stress*t m*m/s') legend('integrated stress', 'transport') subplot(3,1,3) plot(time, ze, time, zK,':') xlabel('time, s') ylabel('bl thickness, m') legend('e-folding', 'sqrt(0.5*K*t)') % energy budget figure(6) clf reset dkedt = diff(ke)/dt; tdkedt = cumsum(diff(time)); plot(time, topwork,':', time, botwork,'-.', ... time, diss, '--', tdkedt, dkedt) legend('top work rate', 'bot work rate', 'dissipation rate','dKE/dt') ylabel('Work rate, m*m/s*s*s') xlabel('time, s') title('energy budget, 20-30 m') % vorticity budget figure(7) clf reset dvortdt = diff(avgvort)/dt; timedv = cumsum(diff(time)); plot(time, netvflux, '-', time, botvflux, '--', time, topvflux,':', ... timedv, dvortdt, '-.') xlabel('time, s') ylabel('vorticity fluxes') legend('-div(flux)', 'bottom v flux', 'top v flux','dvort/dt') title('Vorticity budget, Z = 20-30 m')