Archivo:File3.jpg

2014-05-05T09:34:49Z

Jose manuel torres serrano:

Archivo:File2.jpg

2014-05-05T09:33:48Z

Jose manuel torres serrano:

Archivo:File1.jpg

2014-05-05T09:31:56Z

Jose manuel torres serrano:

Archivo:Grafica Runge Kutta.jpeg

2014-05-05T09:14:10Z

Jose manuel torres serrano:

Archivo:Superficie Runge Kutta.jpg

2014-05-05T09:13:10Z

Jose manuel torres serrano:

Archivo:Euler1.3.jpg

2014-05-05T09:12:38Z

Jose manuel torres serrano:

Archivo:Euler1.2.jpg

2014-05-05T09:12:25Z

Jose manuel torres serrano:

Archivo:Grafica Euler implícito.jpg

2014-05-05T08:31:56Z

Jose manuel torres serrano:

Archivo:Superficie Euler impli.jpeg

2014-05-05T08:31:15Z

Jose manuel torres serrano:

Archivo:Figure1.jpg

2014-05-04T22:43:21Z

Jose manuel torres serrano:

Boundary layer in laminar fluids (Grupo 1B)

2014-03-07T00:51:06Z

Jose manuel torres serrano: /* Horizontal velocity of the fluid */

{{ TrabajoED |Boundary layer in laminar fluids. Grupo 1-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandro Andrés Martínez

David Ayala Díez

Claudia Cózar Coarasa

Lorena de la Fuente Sanz

Marino Rivera Muñoz

José Manuel Torres Serrano }}

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.

= 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>
[[Archivo:Flujo.png|thumb|200px|left|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(\nu) </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.

== 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> .

{{matlab|codigo=
%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

}}

{|
|-
|[[Archivo:Graficaf'(20)M.jpg|thumb|500px|left|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>.

==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> .

{{matlab|codigo=
% 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

}}
{|
|-
| [[Archivo: Graficaf'(20)RK.jpg|thumb|500px|left|Graph of <math>f’(20)</math> for each <math>k</math>]]
|}

==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> .

{{matlab|codigo=
%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

}}
{|
|-
| [[Archivo: Graficaf'(20).jpg|thumb|500px|left|Graph of <math>f’(20)</math> for each <math>k</math>]]
|}
== 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.

==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:

{{matlab|codigo=

%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)')
}}

{|

|-

| [[Archivo: Graficaf'eta.jpg|thumb|800px|left|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 < 0,01</math> , if <math>\eta>\eta _0</math>, is <math>\eta_0\ge5,95</math> for <math> \eta \in \mbox{(0,20)}</math>.

=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>.

{{matlab|codigo=
%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
}}

[[Archivo:Ap4a.png|thumb|800px|centre|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).

= 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

<center><math>\eta_0=y \sqrt[]{ \frac{2}{x}}; y=\frac{\eta_0\sqrt[]{x}}{\sqrt[]{2}}=g(x);</math></center>

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

{{matlab|codigo=
%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')

}}

{|
|-
| [[Archivo: Figureg(x).jpg|thumb|500px|left|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.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

MateWiki - Contribuciones del usuario [es]

Archivo:File3.jpg

Archivo:File2.jpg

Archivo:File1.jpg

Archivo:Grafica Runge Kutta.jpeg

Archivo:Superficie Runge Kutta.jpg

Archivo:Euler1.3.jpg

Archivo:Euler1.2.jpg

Archivo:Grafica Euler implícito.jpg

Archivo:Superficie Euler impli.jpeg

Archivo:Figure1.jpg

Boundary layer in laminar fluids (Grupo 1B)