% Blasius,m Last modified 8/20/2008 % This script animates the solution to the Balsius' model for a boundary % layer flow over a flat plate. This model assumes steady flow, constant % density and viscosity. First one must solve the Blasius differential % equation f'''+0.5*f*f''= 0 with f(0)=0, f'(0)=0 AND f'(infinity)=0 where % the prime denotes differentiation wuth respect to eta=y*(U*rho/mu)^.5. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); eta=linspace(0,10,2001); deta=.005; z1=0; z2=0; z3=.33193;%input('The second derivative at eta = 0 is'); % Found by trial and error assuming various values until the final % condition was satisfied. X=[z1;z2;z3]; for i=1:2000 z1=z1+z2*deta; z2=z2+z3*deta; z3=z3-0.5*z1*z3*deta; X=[X [z1;z2;z3]]; end figure(1);clf; subplot(3,1,1) plot(eta,X(1,:)) xlabel('Variable, \eta') ylabel('Solution, f(\eta)') subplot(3,1,2) plot(eta,X(2,:)) xlabel('Variable, \eta') ylabel('Solution, df(\eta)/d\eta') subplot(3,1,3) plot(eta,X(3,:)) xlabel('Variable, \eta') ylabel('Solution, d^2f(\eta)/d\eta^2') F=eta.*X(2,:)-X(1,:); hold off text(5,1.5,'Press Enter to Continue') pause figure(2);clf; plot(eta,F) xlabel('Variable, \eta') ylabel('Solution, \etadf(\eta)/d\eta-f(\eta)') hold off text(4,.8,'Press Enter to Continue') pause y=linspace(0,.003,51); x=[0 .05 .1 .15 .2 .25 .3 .35 .4 .45 .5]*.1; figure(3);clf; axis([0 1 -.0001 .003]) hold on plot([0 1],[0 0]) hold on u=zeros(51,11); eta3=[]; % (U*rho/mu)=1e5:Corresponds to water at 20 C with U=0.1 m/s for i=2:51 % y loop for j=2:11 % x loop eta1=sqrt(1e5/(x(1,j)))*y(1,i); if y(1,i)==0 fprime=0; fact=0; else fprime=interp1(eta,X(2,:),eta1); fact=interp1(eta,F,eta1); end u(i,j)=fprime; v(i,j)=sqrt(1/(1e5*x(1,j)))*fact; eta3(i,j)=eta1; end end u(:,1)=ones(51,1); v(:,1)=zeros(51,1); v(1,:)=zeros(1,11); for i=39:51 u(i,2)=1; end for j=1:11 plot(u(:,j),y) hold on end text(.8,.6e-3,'x = 0.005'); text(.7,.7e-3,'0.01'); text(.38,1e-3,'0.05'); box on xlabel('Dimensionless Velocity, u(x,y)/U') ylabel('Distance Above the Surface, y') text(0.2,2e-3,['\rhoU/\mu = ' num2str(1e5) ' m^-^1']) hold off text(.2,2.5e-3,'Press Enter to Continue') pause figure(4);clf; axis([-.01 .06 -.0002 .006]) hold on patch([0 .005 .062 .06 0],[0 -.0002 -.0002 0 0],'r') box on plot([0 .004 .004 0 0],[0 0 3e-3 3e-3 0]) patch([.004 .003 .003 .004],[3e-3 (3e-3)-.0001 (3e-3)+.0001 3e-3],'b') hold on for i=2:11 plot(x(1,i)*ones(1,51)+.004*u(:,i)',y(1,:)) xtip=x(1,i)+.004*u(51,i); ytip=y(1,51); plot([x(1,i) x(1,i) xtip],[0 ytip ytip]) patch([xtip xtip-.001 xtip-.001 xtip ],[ytip ytip-.0001 ytip+.0001 ytip],'b') hold on end x3=linspace(0,.06,201); delta3=5*sqrt(x3/1e5); plot(x3,delta3,'m') hold on plot([x3(1,150) .033],[delta3(1,150) 3.5e-3],'m'); text(.004,3.5e-3, 'Boundary Layer Thickness') xlabel('Distance along the plate x, m') ylabel('Distance above the plate y, m') text(-.005,4.8e-3,'U\rho/\mu = 100000 m^-^1, Water at U = 0.1 m/s and 20 C') text(-.005,5.5e-3,'Horizontal Velocity Profiles vs Distance along the Plate') hold off text(-.005, 4e-3,'Press Enter to Continue') pause dt=.00005; x1=-5*ones(8,1); beginningx=-5; y1=[.5e-5 .00001 .00002 .00003 .00004 .00005 .00006 .00007]'*2; hold on v1=zeros(8,1); u1=ones(8,1); figure(5);clf;%Do animation axis([beginningx 200 -2 80]); hold on patch([0 4 200 200 0],[0 -2 -2 0 0],'r') x1new=.000001*ones(8,1); y1new=y1; Udt=.1*dt; plot([x1 x1new]'/Udt,[y1 y1new]'/Udt, 'k'); box on xlabel('Dimensionless Location, x/U\Deltat') ylabel('Dimensionless Location, y/U\Deltat') text(10,70 ,['\rhoU/\mu = ' num2str(1e5) ' m^-^1, U = 0.1 m/s']) T=text(100,5,'Press Enter to Animate'); text(10,65,'\Deltat = 0.00005 s') pause set(T,'String',' '); x1=x1new; y1=y1new; for t=1:200%time loop for i=1:8%streamline loop eta3=sqrt(1e5/(x1(i,1)))*y1(i,1); if eta3>=10 fprime=1; F1=1.7233; else fprime=interp1(eta, X(2,:),eta3); F1=interp1(eta,F,eta3); end u1(i,1)=fprime; v1(i,1)=sqrt(1/(4e5*x1(i,1)))*F1; x1new(i,1)=x1(i,1)+dt*u1(i,1)*.1; y1new(i,1)=y1(i,1)+dt*v1(i,1)*.1; hold on end plot([x1 x1new]'/Udt,[y1 y1new]'/Udt, 'k'); x1=x1new; y1=y1new; pause(.01) end x4=linspace(0,200,201); delta4=5*sqrt(x4/.50); plot(x4,delta4,'m','Linewidth',2) text(100,62,'Boundary Layer Thickness') plot([x4(1,90),97],[delta4(1,90),63],'m','LineWidth',2); hold off set(T,'String','Press Enter to Continue') pause x1=-5*ones(8,1); beginningx=-5; y1=[.5e-5 .00001 .00002 .00003 .00004 .00005 .00006 .00007]'*5; v1=zeros(8,1); u1=ones(8,1); figure(6);clf;%Do animation axis([beginningx 200 -3 120]); hold on patch([0 4 200 200 0],[0 -3 -3 0 0],'r') x1new=.000001*ones(8,1); y1new=y1; plot([x1 x1new]'/Udt,[y1 y1new]'/Udt, 'k'); box on xlabel('Dimensionless Location, x/U\Deltat') ylabel('Dimensionless Location, y/U\Deltat') text(20,110 ,['\rhoU/\mu = ' num2str(1e5) ' m^-^1, U = 0.1 m/s']) %text(20,110 ,['\rhoU/\mu = ' num2str(1e5) ' m^-^1']) T=text(120,10,'Press Enter to Animate'); text(20,100,'\Deltat = 0.00005 s') pause set(T,'String',' '); x1=x1new; y1=y1new; u8=[]; for t=1:200%time loop for i=1:8%streamline loop eta3=sqrt(1e5/(x1(i,1)))*y1(i,1); if eta3>=10 fprime=1; F1=1.7233; else fprime=interp1(eta, X(2,:),eta3); F1=interp1(eta,F,eta3); end u1(i,1)=fprime;%actuallu u/U v1(i,1)=sqrt(1/(4e5*x1(i,1)))*F1;%actually v/U x1new(i,1)=x1(i,1)+dt*u1(i,1)*.1; y1new(i,1)=y1(i,1)+dt*v1(i,1)*.1; hold on end plot([x1 x1new]'/Udt,[y1 y1new]'/Udt, 'k'); hold on x1=x1new; y1=y1new; pause(.01) end plot(x4,delta4,'m','Linewidth',2) plot([x4(1,90),99],[delta4(1,90),19],'m','Linewidth',2); text(100,18,'Boundary Layer Thickness') set(T,'String','Press Enter to End') hold off pause