% 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,:);
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)')
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'])
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')
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);
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')
pause