Revisión del 11:01 7 mar 2014

In this numerical project we have studied what happens when we introduce a flat plate in a laminar fluid whose speed and viscosity are constant.

1 Blasius equation

First, we must assume that the fluid velocity before reaching the plate is constant, as in remote areas after passing the plate.

In our case we take this constant as [math]2[/math], such that

[math]\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 [/math]

Laminar fluid with velocity [math] u_0 [/math]

Then, we must define the fluid stream function

[math]\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)[/math]

Where we take the viscosity [math] \nu [/math] as a unit value and

[math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math]

[math] f(\eta) [/math] Satisfies the Blasius equation, and therefore we will raise the initial value problem associated with this equation with the following initial conditions

[math] \begin{cases} f’’’(\eta)+\frac{1}{2}f(\eta)f’’(\eta)=0 ; \\ f(0)=f’(0)=0, \lim_{\eta \to \infty}f’(\eta)= 1 ; \end{cases} [/math]

However, on time of programming we cannot introduce a conditional limit, so we have to replace them by the condition [math] f’’(\eta)=k [/math], and vary the values of k to find the one that satisfies the limit.

We cannot solve the differential equation like such, to apply the numerical methods, we need to pass it to a system of equations:

[math] f(\eta)=y_1,f’(\eta)=y_2,f’’(\eta)=y_3\\ \begin{cases} y_1’=y_2;\\ y_2’=y_3;\\ y_3’=-\frac{1}{2}y_1y_3;\\ y_1(0)=y_2(0)=0; y_3(0)=k \end{cases} [/math]

Once the system has been formulated , we start to solve it.

1.1 Resolution with the modified Euler method

Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of [math]k[/math] , with [math] k \in \mbox{(0,1;1)}[/math] and [math]dk=0,01[/math] ,with [math] \eta \in \mbox{(0,20)}[/math] and with [math]h=0,05[/math] .

%Resolution of Blasius equation(with modified Euler)
clear all
%Initial conditions
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);%We create the matrix F2 where we will store the different
%solutions of f2 for each value of k
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
    A=[0 1 0;0 0 1;(-y(3)/2) 0 0]; %To simplify and solve using matrices, we create
 %the matrix A in the loop with different values of f3
    z=y+h*A*y;
    y=y+(h/2)*(A*y+A*z);
    y1(n+1)=y(1);
    y2(n+1)=y(2);
    y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k
F2(num,:)=y2;
end
k1=[0.1:0.01:1]; %Vector to represent the values of f2 in 20
f20=F2(:,401);
f20=f20';
o=ones(1,91);%We represented f = 1 for better viewing
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r') 
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

Graph of [math]f’(20)[/math] for each [math]k[/math]

As noted in the graph the value for which the function is closer to [math]1[/math] is [math]k=0,33[/math].

1.2 Resolution with 4th order Runge-Kutta method

Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of [math]k[/math] , with [math] k \in \mbox{(0,1;1)}[/math] and [math]dk=0,01[/math] ,with [math] \eta \in \mbox{(0,20)}[/math] and with [math]h=0,05[/math] .

% Resolution of Blasius  equation(with Runge-Kutta)
clear all
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
    A=[0 1 0;0 0 1;(-y(3)/2) 0 0];
    k1=A*y;
    k2=A*(y+(h/2)*k1);
    k3=A*(y+(h/2)*k2);
    k4=A*(y+h*k3);
    y=y+(h/6)*(k1+2*k2+2*k3+k4);
    y1(n+1)=y(1);
    y2(n+1)=y(2);
    y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k 
F2(num,:)=y2;
end
t=[t0:h:tN];
k1=[0.1:0.01:1];
f20=F2(:,401);
f20=f20';
o=ones(1,91);
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r')
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

Graph of [math]f’(20)[/math] for each [math]k[/math]

1.3 Resolution with Euler method

Then is exposed the Matlab code that numerically that solves the Blasius equation for different values of [math]k[/math] , with [math] k \in \mbox{(0,1;1)}[/math] and [math]dk=0,01[/math] ,with [math] \eta \in \mbox{(0,20)}[/math] and with [math]h=0,05[/math] .

%Resolution of Blasius  equation(with Euler)

clear all
t0=0;
tN=20;
h=0.05;
N=(tN-t0)/h;
F2=zeros(91,401);
for k=0.1:0.01:1
y=[0;0;k];
y1=y(1);
y2=y(2);
y3=y(3);
for n=1:N
    A=[0 1 0;0 0 1;(-y(3)/2) 0 0];
    y=y+h*A*y;
    y1(n+1)=y(1);
    y2(n+1)=y(2);
    y3(n+1)=y(3);
end
t=[t0:h:tN];
num=int8(100*(k-0.1+0.01));
%F2 has as rows approximations of y2 for the different values of k 
F2(num,:)=y2;
end
t=[t0:h:tN];
k1=[0.1:0.01:1];
f20=F2(:,401);
f20=f20';
o=ones(1,91);
figure(1)
hold on
plot(k1,f20,'+')
plot(k1,o,'r')
xlabel('k')
ylabel('f´(20)')
legend('f´(20)','y=1')
hold off

Graph of [math]f’(20)[/math] for each [math]k[/math]

1.4 Conclusion

Comparing the graphs, we see that the difference between modified Euler and 4th order Runge Kutta methods is minimal and the value for each parameter is [math]k=0,33[/math] , on the other hand, by using the Euler method (less accurate than the above) the value of the parameter is [math]k=0,32[/math], although graphically this difference is hardly seen.

1.5 Graph of [math]f´(\eta)[/math]

The graph of [math]f'(\eta)[/math]has been realized in [math]\eta \in \mbox {(0,20)}[/math] for the value of parameter [math]k[/math] obtained in the modified Euler method, ie [math]k=0.33[/math].

To get this graph we add into the modified Euler method, the following MATLAB code:

%We plot f2 for k = 0.33 which is in row 24 of the matrix
f2=F2(24,:);
figure(2)
plot(t,f2,'*')
xlabel('\eta')
ylabel('f´(\eta)')

Graph of [math]f´(\eta)[/math] for [math]k=0,33[/math]

As can be seen, if we run the full program and see the vector [math]f2[/math] the value of [math]\eta_0[/math] for which [math]\ \vert f'(\eta)-1 \vert \lt 0,01[/math] , if [math]\eta\gt\eta _0[/math], is [math]\eta_0\ge5,95[/math] for [math] \eta \in \mbox{(0,20)}[/math].

2 Horizontal velocity of the fluid

Once we have numerically calculated the [math]f( \eta)[/math] we proceed to calculate the horizontal component of the fluid velocity [math] u_1 [/math]:

[math]\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})[/math] Thus, as defined above:

[math]u_1=\frac{\partial \psi}{\partial y}=\frac{\partial }{\partial y}(\sqrt[]{\nu u_0 x} f(\eta))= \sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial y}=\sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial \eta} \frac{\partial \eta}{\partial y}=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)[/math]

To translate this result graphically we calculate [math]u_1(x_k,y)[/math] with the Modified Euler method where [math]x_k=0.05,0.2,0.4,0.6,0.8[/math] and [math]y \in \mbox{(0,3)}[/math] with [math]h=0.01[/math].

%We calculate u1 with the different values of xk(with modified Euler)
clear all
xk=[0.05,0.2,0.4,0.6,0.8];
nu=1; u0=2;
y0=0; yN=3; hy=0.01;
N=(yN-y0)/hy;
y=y0:hy:yN;
for m=1:5
%We define eta ('t') for each value of  xk, each one 
%corresponding to a row, and with 'y' in (0,3)
t(m,:)=sqrt(u0/(nu*xk(m)))*y;
h=sqrt(u0/(nu*xk(m)))*0.01;
f0=[0;0;0.33];
f=[f0(1);f0(2);f0(3)];
for n=1:N
    A=[0 1 0;0 0 1;(-f0(3)/2) 0 0];
    z=f0+h*A*f0;
    f0=f0+(h/2)*(A*f0+A*z);
    f(:,n+1)=[f0(1);f0(2);f0(3)];
end
Y(m,:)=f(2,:);
end
F=u0*Y;
hold on
plot(y,F(1,:),'k')
plot(y,F(2,:))
plot(y,F(3,:),'r')
plot(y,F(4,:),'m')
plot(y,F(5,:),'g')
legend('u_{1} for x_{k}=0.05','u_{1} for x_{k}=0.2','u_{1} for x_{k}=0.4','u_{1} for x_{k}=0.6','u_{1} for x_{k}=0.8','location','best')
xlabel('y')
ylabel('u_{1}')
hold off

In this picture are shown [math]5[/math] graphs, each one corresponds to a different value of [math]k[/math].

According to this graph we can appreciate that the fluid, when is moving along the [math]x[/math] axis, has to achieve higher height to get limit velocity [math]u_0[/math], i.e., when the fluid moves, it must be getting over the plate to offset the perturbation that the plate causes to it ( to a higher value of [math]y[/math] is greater the transition zone between zero velocity of the plate and the limit velocity [math]u_0[/math] with which the fluid initially starts).

3 Laminar boundary layer of the fluid

According to the above findings, it can be deduced that there is for each value of [math]x[/math] a limit value [math]y[/math] from which the fluid velocity becomes constant speed again, with the same value that it had initially before reaching the area of the plate. Obviously the value of the boundary layer will be related to the value [math]\eta_0[/math], calculated above, for which the function [math]f’(\eta)[/math] will become almost constant. The relationship is expressed as follows: [math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math]

If [math]\eta=\eta_0[/math], then

Therefore we can interpret this function [math]g(x)[/math] as the fluid boundary layer. Simple Matlab code is exposed for its representation:

%plot function g(x)
x=[0:0.05:10];
eta0=5.95;
y=eta0*(x/2).^(1/2);
plot(x,y,'r')
xlabel('x')
ylabel('g(x)')
legend('g(x), interpreted as boundary layer','location','best')

Graph of [math]g(x)[/math]

In view of the graph, the findings are similar to previous ones. As the value [math]x[/math] increases, you need a higher value of [math]y[/math] for the fluid velocity stabilizes and becomes the speed that we had initially.