|
|
| Línea 57: |
Línea 57: |
| | | | |
| | Once the system has been formulated , we start to solve it. | | 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);1
| |
| − | 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 provokes 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]]
| |
First, we must assume that the fluid velocity before reaching the plate is constant, as in remote areas after passing the plate.
