Diferencia entre revisiones de «Usuario:Grupo 4»

De MateWiki
Saltar a: navegación, buscar
(→‎Resolución numérica del problema.)
(→‎Representación gráfica e interpretación del borde de la capa limite)
 
(No se muestran 106 ediciones intermedias del mismo usuario)
Línea 2: Línea 2:
  
 
==Introducción.==
 
==Introducción.==
−
En el trabajo que nos ocupa se va a estudiar la capa límite de un fluido laminar sobre una placa plana. todo esto se va a efectuar mediante métodos numéricos ayudados con el programa MATLAB.
+
En el trabajo que nos ocupa se va a estudiar la capa límite de un fluido laminar sobre una placa plana. Todo esto se va a efectuar mediante métodos numéricos ayudados con el programa MATLAB.
 +
[[Archivo:1515.png|marco|derecha]]
 +
 
 
Para resolver el problema tomaremos que la placa ocupa la semirrecta "y=0" estudiando lo que ocurre por encima de la misma, es decir, en los puntos "(x,y) con x>0 y y>0". Antes de llegar a la placa, el fluido mantiene una velocidad constante:
 
Para resolver el problema tomaremos que la placa ocupa la semirrecta "y=0" estudiando lo que ocurre por encima de la misma, es decir, en los puntos "(x,y) con x>0 y y>0". Antes de llegar a la placa, el fluido mantiene una velocidad constante:
−
con "u<sub>0</sub>=2",al igual que una vez superada. En cambio en los puntos de la placa la velocidad será nula existiendo también una zona de transición. En esa zona la corriente del fluido "ψ(x,y)" viene dada por la función:
 
  
−
siendo "η":  
+
[[Archivo:333.png|sin marco|centro]]
 +
 
 +
con "u<sub>0</sub>=2",al igual que una vez superada. En cambio en los puntos de la placa la velocidad será nula existiendo también una zona de transición. En esa zona la corriente del fluido "ψ(x,y)" viene dada por la función "ψ" y siendo la variable de similaridad "η":  
 +
 
 +
[[Archivo:33.png|sin marco|centro]]
  
 
donde la viscosidad del fluido "ν" será (ν=1) y "f" satisface la ecuación de Blasius:
 
donde la viscosidad del fluido "ν" será (ν=1) y "f" satisface la ecuación de Blasius:
 +
 +
[[Archivo:3223.png|sin marco|centro]]
  
 
Las condiciones iniciales serán: f(0)=f'(0)=0 y:
 
Las condiciones iniciales serán: f(0)=f'(0)=0 y:
 +
 +
[[Archivo:222.png|sin marco|centro]]
  
 
El campo de velocidades del fluido viene dado por:
 
El campo de velocidades del fluido viene dado por:
 +
 +
[[Archivo:2222.png|sin marco|centro]]
  
 
A continuación resolveremos la Ecuación de Blasius con las condiciones iniciales dadas mediante el uso de métodos numéricos buscando a su vez el valor inicial óptimo de la segunda derivada del factor de forma "k".
 
A continuación resolveremos la Ecuación de Blasius con las condiciones iniciales dadas mediante el uso de métodos numéricos buscando a su vez el valor inicial óptimo de la segunda derivada del factor de forma "k".
  
 
==Resolución numérica del problema.==
 
==Resolución numérica del problema.==
 +
 +
[[Archivo:5555.png|centro]]
  
 
===Resolución mediante Euler modificado.===
 
===Resolución mediante Euler modificado.===
 +
 +
En este primer apartado se va a resolver la ecuación de Blasius, en función de las siguientes condiciones iniciales:
 +
[[Archivo:33333.png|sin marco|centro]]
 +
 +
donde "k" toma valores de 0,1 a 1 con un intervalo de 0,01; mediante el método numérico de Euler modificado siendo "h=0,05".
 +
<br />
 +
 +
 +
[[Archivo:EULERMODIFICADO.png|marco|derecha]]
 +
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%En este caso la variable "t" es la variable de similaridad
−
%Datos iniciales
+
 
%Numero de ecuaciones
 
%Numero de ecuaciones
 
ne=3;
 
ne=3;
−
%Tiempo inicial
+
%Se define el intervalo temporal
 
t0=0;
 
t0=0;
−
%Tiempo final
 
 
tn=20;
 
tn=20;
 
%Tamaño del salto
 
%Tamaño del salto
Línea 35: Línea 56:
 
%Determinacion del vector t
 
%Determinacion del vector t
 
t=linspace(t0,tn,N+1);
 
t=linspace(t0,tn,N+1);
−
%Vector k, valor inicial de f(3)
+
%Se establece el vector K
 
k0=0.1;
 
k0=0.1;
 
kf=1;
 
kf=1;
 
k=k0:0.01:kf;
 
k=k0:0.01:kf;
 
G=length(k);
 
G=length(k);
−
%Sistema de ecuaciones
+
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
f=zeros(ne,N+1);
 
f=zeros(ne,N+1);
−
%Resolucion
+
%Se resuelve el problema mediante el método de Euler modificado aplicado a sistemas
 
for j=1:G
 
for j=1:G
 
   f0=[0 0 k(j)];
 
   f0=[0 0 k(j)];
Línea 55: Línea 76:
 
   Y(j)=f(2,end);
 
   Y(j)=f(2,end);
 
end
 
end
 +
%Representación gráfica de la primera derivada del factor de forma para n=20
 +
hold on
 
plot(k,Y)
 
plot(k,Y)
 +
legend('f´(20)')
 +
xlabel('k')
 +
ylabel('f´(20)')
 +
title('Representación de f(20)´ para los posibles k')
 +
hold off
 +
%Valor de k para el cual f'(20) se acerca mas a 1
 
for q=1:G
 
for q=1:G
 
   a(q)=abs(Y(q)-1);
 
   a(q)=abs(Y(q)-1);
Línea 68: Línea 97:
 
}}
 
}}
  
−
[[Archivo:EULERMODIFICADO.png|marco|derecha]]
+
El parámetro de k para el cual f'(20) se acerca mas a 1 es k=0.33
  
 
===Resolución mediante Runge-Kutta de orden 4.===
 
===Resolución mediante Runge-Kutta de orden 4.===
  
 +
En este segundo apartado se resuelve la ecuación de Blasius con las mismas condiciones iniciales, pero en este caso utilizando el método de Runge-Kutta.
 +
<br />
 +
 +
[[Archivo:RUNGEKUTTA2.png|marco|derecha]]
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%En este caso la variable "t" es la variable de similaridad
−
%Datos iniciales
+
 
%Numero de ecuaciones
 
%Numero de ecuaciones
 
ne=3;
 
ne=3;
−
%Tiempo inicial
+
%Se define el intervalo temporal
 
t0=0;
 
t0=0;
−
%Tiempo final
 
 
tn=20;
 
tn=20;
 
%Tamaño del salto
 
%Tamaño del salto
Línea 88: Línea 119:
 
%Determinacion del vector t
 
%Determinacion del vector t
 
t=linspace(t0,tn,N+1);
 
t=linspace(t0,tn,N+1);
−
%Vector k, valor inicial de f(3)
+
%Se establece el vector K
 
k0=0.1;
 
k0=0.1;
 
kf=1;
 
kf=1;
 
k=k0:0.01:kf;
 
k=k0:0.01:kf;
 
G=length(k);
 
G=length(k);
−
%Sistema de ecuaciones
+
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
f=zeros(ne,N+1);
 
f=zeros(ne,N+1);
−
%Resolucion
+
%Se resuelve el problema mediante el método de Runge-Kutta de orden 4 aplicado a sistemas
 
for j=1:G
 
for j=1:G
 
   f0=[0 0 k(j)];
 
   f0=[0 0 k(j)];
 
   f(:,1)=f0;
 
   f(:,1)=f0;
 
for i=1:N
 
for i=1:N
 +
      %Definicion de los parametros K1,K2,K3 y K4 necesarios para la resolución
 
       K1=w(t(i),f(:,i));
 
       K1=w(t(i),f(:,i));
 
       t(i)=t(i)+h/2;
 
       t(i)=t(i)+h/2;
Línea 113: Línea 145:
 
   Y(j)=f(2,end);
 
   Y(j)=f(2,end);
 
end
 
end
−
plot(k,Y)
+
%Representación gráfica de la primera derivada del factor de forma para n=20
 +
hold on
 +
plot(k,Y,'r')
 +
legend('f´(20)')
 +
xlabel('k')
 +
ylabel('f´(20)')
 +
title('Representación de f(20)´ para los posibles k')
 +
hold off
 +
%Valor de k para el cual f'(20) se acerca mas a 1
 
for q=1:G
 
for q=1:G
 
     a(q)=abs(Y(q)-1);
 
     a(q)=abs(Y(q)-1);
Línea 125: Línea 165:
 
end
 
end
 
}}
 
}}
 +
 +
El parámetro de k para el cual f'(20) se acerca mas a 1 es k=0.33
  
 
===Resolución mediante Euler.===
 
===Resolución mediante Euler.===
  
 +
En este tercer apartado se resuelve la ecuación de Blasius en las mismas condiciones, pero en este caso mediante el método numérico de Euler.
 +
<br />
 +
[[Archivo:EULER33.png|marco|derecha]]
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%En este caso la variable "t" es la variable de similaridad
−
%Datos iniciales
+
 
%Numero de ecuaciones
 
%Numero de ecuaciones
 
ne=3;
 
ne=3;
−
%Tiempo inicial
+
%Se define el intervalo temporal
 
t0=0;
 
t0=0;
−
%Tiempo final
 
 
tn=20;
 
tn=20;
 
%Tamaño del salto
 
%Tamaño del salto
Línea 144: Línea 187:
 
%Determinacion del vector t
 
%Determinacion del vector t
 
t=linspace(t0,tn,N+1);
 
t=linspace(t0,tn,N+1);
−
%Vector k, valor inicial de f(3)
+
%Se establece el vector K
 
k0=0.1;
 
k0=0.1;
 
kf=1;
 
kf=1;
 
k=k0:0.01:kf;
 
k=k0:0.01:kf;
 
G=length(k);
 
G=length(k);
−
%Sistema de ecuaciones
+
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
f=zeros(ne,N+1);
 
f=zeros(ne,N+1);
−
%Resolucion
+
%Se resuelve el problema mediante el método de Euler aplicado a sistemas
 
for j=1:G
 
for j=1:G
 
   f0=[0 0 k(j)];
 
   f0=[0 0 k(j)];
Línea 161: Línea 204:
 
   Y(j)=f(2,end);
 
   Y(j)=f(2,end);
 
end
 
end
−
plot(k,Y)
+
%Representación gráfica de la primera derivada del factor de forma para n=20
 +
hold on
 +
plot(k,Y,'g')
 +
legend('f´(20)')
 +
xlabel('k')
 +
ylabel('f´(20)')
 +
title('Representación de f(20)´ para los posibles k')
 +
hold off
 +
%Valor de k para el cual f'(20) se acerca mas a 1
 
for q=1:G
 
for q=1:G
 
   a(q)=abs(Y(q)-1);
 
   a(q)=abs(Y(q)-1);
Línea 175: Línea 226:
 
end
 
end
 
}}
 
}}
 +
 +
Parametro de k para el cual f'(20) se acerca mas a 1 es k=0.32
  
 
===Comparativa y elección de k===
 
===Comparativa y elección de k===
 +
Analizando los resultados de los tres métodos utilizados para el estudio de la ecuación de Blasius, y fijándonos en las gráficas y en el valor de K obtenido de cada uno de ellos para que f'(20) sea lo más próximo a la unidad, podemos afirmar que la diferencia entre ellos es muy pequeña. El valor obtenido para los métodos Runge-kutta y Euler modificado es de 0,33 y para el método de Euler es de 0,32. Esto último se debe a una mayor precisión local de los dos primeros métodos numéricos y un mayor error en el de Euler siendo en este caso prácticamente despreciable debido a que "h" (valor de salto) es muy pequeño.<br />
 +
 +
El método que vamos a elegir para determinar el valor de k es ek de Runge-kutta de orden 4, ya que el error de paso es de orden O(h<sup>5</sup>) mientras que el error total acumulado es O(h<sup>4</sup>), razón por la cual es un método muy preciso y además el más usado en los desarrollos computacionales. El inconveniente a la hora de usar este método es una mayor cantidad de evaluaciones de la función, resultando un mayor tiempo de cálculo y una mayor solicitación de la computadora.
  
 
==Interpretación de la derivada de la función de forma f(η)==
 
==Interpretación de la derivada de la función de forma f(η)==
 +
En este caso, vamos a estudiar el valor de la variable de similaridad "η" a partir del cual la diferencia entre cualquier valor de la primera derivada de la función de forma y 1 es menor de 0,01. El método numérico escogido para la resolución de este problema es el de Runge-kutta debido a las propiedades mencionadas en el anterior apartado y la k obtenida es 0,33. Como se puede apreciar en el gráfico, el valor de η para el cual la función f'(η) se acerca a menos de 0,01 de 1, que es su límite cuando la variable de similaridad tiende a infinito, es de 5,225 en nuestro caso.
 +
<br />
 +
Realizando otra interpretación en la búsqueda de este valor podemos encontrarnos con el concepto de grosor de capa, siendo éste el valor para el cual la primera derivada del factor de forma llega al 99% de su límite.
 +
<br />
 +
[[Archivo:APARTADO3.png|marco|derecha]]
  
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%En este caso la variable "t" es la variable de similaridad
−
%Datos iniciales
+
 
%Numero de ecuaciones
 
%Numero de ecuaciones
 
ne=3;
 
ne=3;
−
%Tiempo inicial
+
%Se define el intervalo temporal
 
t0=0;
 
t0=0;
−
%Tiempo final
 
 
tn=20;
 
tn=20;
 
%Tamaño del salto
 
%Tamaño del salto
Línea 196: Línea 255:
 
%Determinacion del vector t
 
%Determinacion del vector t
 
t=linspace(t0,tn,N+1);
 
t=linspace(t0,tn,N+1);
−
%Vector k, valor inicial de f(3)
+
%Determinacion del vector t
 
k0=0.1;
 
k0=0.1;
 
kf=1;
 
kf=1;
 
k=k0:0.01:kf;
 
k=k0:0.01:kf;
 
G=length(k);
 
G=length(k);
−
%Sistema de ecuaciones
+
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
f=zeros(ne,N+1);
 
f=zeros(ne,N+1);
−
%Resolucion
+
%Se resuelve el problema mediante el método de Runge-Kutta de orden 4 aplicado a sistemas para k=0.33
−
  f0=[0 0 k(24)];
+
f0=[0 0 k(24)];
−
  f(:,1)=f0;
+
f(:,1)=f0;
−
  for i=1:N
+
for i=1:N
−
  K1=w(t(i),f(:,i));
+
  K1=w(t(i),f(:,i));
−
  t(i)=t(i)+h/2;
+
  t(i)=t(i)+h/2;
−
  z(:,i)=f(:,i)+K1*h/2;
+
  z(:,i)=f(:,i)+K1*h/2;
−
  K2=w(t(i),z(:,i));
+
  K2=w(t(i),z(:,i));
−
  v(:,i)=f(:,i)+K2*h/2;
+
  v(:,i)=f(:,i)+K2*h/2;
−
  K3=w(t(i),v(:,i));
+
  K3=w(t(i),v(:,i));
−
  c(:,i)=f(:,i)+K3*h;
+
  c(:,i)=f(:,i)+K3*h;
−
  K4=w(t(i)-h/2,c(:,i));
+
  K4=w(t(i)-h/2,c(:,i));
−
  f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);
+
  f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);
−
 
+
end
−
  end
+
O=ones(1,N+1);
−
O=ones(1,N+1);
+
%Representación de  f´ para k=0,33
−
    plot(t,f(2,:))
+
hold on
−
  for j=1:N+1   
+
plot(t,f(2,:))
 +
legend('f´(n)')
 +
xlabel('n')
 +
ylabel('f´(n)')
 +
title('Representación de f´ para k=0,33')
 +
hold off
 +
%Determinación de n para que la diferencia entre f' y 1 sea menor de 0.01
 +
for j=1:N+1   
 
   if abs(f(2,j)-O)<0.01
 
   if abs(f(2,j)-O)<0.01
 
       t(j)
 
       t(j)
−
      j
 
 
       break
 
       break
 
   end
 
   end
 
end
 
end
 +
for m=1:N+1
 +
    b(m,:)=abs(f(2,m)-O);
 +
end
 +
%Representacion de la diferencia en funcion de n
 +
figure(2)
 +
hold on
 +
plot(t,b)
 +
xlabel('n')
 +
ylabel('Modulo de la diferencia')
 +
title('Representacion de la diferencia en funcion de n')
 +
hold off
 +
 
}}
 
}}
  
 
==Estudio de la velocidad en función del factor de forma  f(η)==
 
==Estudio de la velocidad en función del factor de forma  f(η)==
  
−
===Resolución numérica y representación gráfica de la componente U<sub>1</sub> del campo de velocidades===
+
<br />
  
 +
[[Archivo:2323.png|centro]]
 +
 +
 +
===Resolución numérica y representación gráfica de la componente U<sub>1</sub> del campo de velocidades===
 +
Sabiendo que la primera componente u<sub>1</sub> del campo de velocidades del fluido es proporcional a la primera derivada del factor de forma, siendo el factor de proporcionalidad la velocidad constante u<sub>0</sub>=2, podemos en este caso dar utilidad al método numérico de Euler modificado usado en el apartado 2.1 para la representación gráfica de u<sub>1</sub>. Elegimos este método porque el resultado obtenido para la k es el mismo que en el de Runge-kutta, es decir nos aporta buena precisión, a la vez que reduce carga de trabajo a la computadora.
 +
<br />
 +
[[Archivo:APARTADO4.png|marco|derecha]]
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%En este caso la variable "t" es la variable de similaridad
−
%Datos iniciales
+
%Numero de ecuacionesne=3;
−
ne=3;
+
 
h=0.05;
 
h=0.05;
 +
%Vectores x e y
 
x=[0.05 0.2 0.4 0.6 0.8];
 
x=[0.05 0.2 0.4 0.6 0.8];
 
y0=0;
 
y0=0;
Línea 246: Línea 330:
 
P=length(y);
 
P=length(y);
 
N=round((yn-y0)/h);
 
N=round((yn-y0)/h);
 +
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
 
f=zeros(ne,N+1);
 
f=zeros(ne,N+1);
 +
%Se resuelve el problema mediante el método de Euler modificado aplicado a sistemas
 
for j=1:length(x)
 
for j=1:length(x)
 
   H2=h*sqrt(2/x(j));
 
   H2=h*sqrt(2/x(j));
Línea 259: Línea 345:
 
   end
 
   end
 
   Z(j,:)=f(2,:);
 
   Z(j,:)=f(2,:);
 +
  %Representación de la primera variable del campo de velocidades para las diferentes x
 
   hold on  
 
   hold on  
 
   plot(y,2*Z(j,:))
 
   plot(y,2*Z(j,:))
−
  hold off
 
 
end
 
end
 +
legend('x=0.05','x=0.2','x=0.4','x=0.6','x=0.8')
 +
xlabel('y')
 +
ylabel('U1')
 +
title ('Representacion de U1')
 +
hold off
 +
 
}}
 
}}
  
 
===Observaciones===
 
===Observaciones===
 +
 +
Interpretando las diferentes gráficas para los valores de x indicados, se observa que, cuanto mayor es el valor de la componente horizontal, mayor es la altura que necesita el fluido para alcanzar la velocidad constante u<sub>0</sub>=2, lo que nos sugiere la idea de que mayor es el grosor de la capa límite.
  
 
==Representación gráfica e interpretación del borde de la capa limite==
 
==Representación gráfica e interpretación del borde de la capa limite==
  
 +
En esta última ilustración se observa gráficamente el borde de la capa límite exactamente igual como si se viera en un laboratorio realizando el experimento físicamente con los datos indicados, dando a entender el comportamiento real del fluido en esta situación. Por debajo del borde el fluido comienza a reducir su velocidad hasta llegar a ser nula justo en la superficie de la placa, análogamente, a medida que nos alejamos del origen de ella, el fluido tiende a volver a ganar la velocidad constante anterior al contacto.
 +
<br /><br />
 +
 +
[[Archivo:Ejemploreal1.png|sinmarco|centro|460px]]
 +
 +
 +
[[Archivo:APARTADO5.png|sinmarco|derecha|505px]]
 
{{matlab|codigo=
 
{{matlab|codigo=
−
clc
+
%Se definen los datos del problema
−
clear all
+
%Vector x (valor de la longitud de la chapa)
 
x0=0;
 
x0=0;
 
xn=10;
 
xn=10;
 
h=0.05;
 
h=0.05;
 
x=x0:h:xn;
 
x=x0:h:xn;
−
eta=5.225;
+
%Valor de la variable de similaridad en el borde de la capa
 +
n=5.225;
 
for i=1:length(x)
 
for i=1:length(x)
−
     y(i)=eta*sqrt(x(i)/2);
+
     y(i)=n*sqrt(x(i)/2);
 
end
 
end
 +
%Representación de la capa de transición
 +
hold on
 
plot(x,y);
 
plot(x,y);
 +
legend('g(x)')
 +
xlabel('X')
 +
ylabel('Y')
 +
title ('Borde de la capa de transición')
 +
hold off
 +
 
}}
 
}}
  

Revisión actual del 03:12 30 abr 2016

Trabajo realizado por estudiantes
Título Fluido por encima de una placa plana. Grupo 4
Asignatura Ecuaciones Diferenciales
Curso Curso 2015-16
Autores Ruben Martos López,
Guillermo Megino León,
Silviu Popa
Alejandro Sistac Ara,
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción.

En el trabajo que nos ocupa se va a estudiar la capa límite de un fluido laminar sobre una placa plana. Todo esto se va a efectuar mediante métodos numéricos ayudados con el programa MATLAB.

derecha

Para resolver el problema tomaremos que la placa ocupa la semirrecta "y=0" estudiando lo que ocurre por encima de la misma, es decir, en los puntos "(x,y) con x>0 y y>0". Antes de llegar a la placa, el fluido mantiene una velocidad constante:

centro

con "u0=2",al igual que una vez superada. En cambio en los puntos de la placa la velocidad será nula existiendo también una zona de transición. En esa zona la corriente del fluido "ψ(x,y)" viene dada por la función "ψ" y siendo la variable de similaridad "η":

centro

donde la viscosidad del fluido "ν" será (ν=1) y "f" satisface la ecuación de Blasius:

centro

Las condiciones iniciales serán: f(0)=f'(0)=0 y:

centro

El campo de velocidades del fluido viene dado por:

centro

A continuación resolveremos la Ecuación de Blasius con las condiciones iniciales dadas mediante el uso de métodos numéricos buscando a su vez el valor inicial óptimo de la segunda derivada del factor de forma "k".

2 Resolución numérica del problema.

centro

2.1 Resolución mediante Euler modificado.

En este primer apartado se va a resolver la ecuación de Blasius, en función de las siguientes condiciones iniciales: centro

donde "k" toma valores de 0,1 a 1 con un intervalo de 0,01; mediante el método numérico de Euler modificado siendo "h=0,05".


derecha
%Se definen los datos del problema
%En este caso la variable "t" es la variable de similaridad
%Numero de ecuaciones
ne=3;
%Se define el intervalo temporal
t0=0;
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Se establece el vector K
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Se resuelve el problema mediante el método de Euler modificado aplicado a sistemas
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      K1=w(t(i),f(:,i));
      t(i)=t(i)+h/2;
      z(:,i)=f(:,i)+K1*h/2;
      f(:,i+1)=f(:,i)+h*w(t(i),z(:,i));
end
   Y(j)=f(2,end);
end
%Representación gráfica de la primera derivada del factor de forma para n=20
hold on
plot(k,Y)
legend('f´(20)')
xlabel('k')
ylabel('f´(20)')
title('Representación de f(20)´ para los posibles k')
hold off
%Valor de k para el cual f'(20) se acerca mas a 1
for q=1:G
   a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
      k(q)
      break
end


El parámetro de k para el cual f'(20) se acerca mas a 1 es k=0.33

2.2 Resolución mediante Runge-Kutta de orden 4.

En este segundo apartado se resuelve la ecuación de Blasius con las mismas condiciones iniciales, pero en este caso utilizando el método de Runge-Kutta.

derecha
%Se definen los datos del problema
%En este caso la variable "t" es la variable de similaridad
%Numero de ecuaciones
ne=3;
%Se define el intervalo temporal
t0=0;
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Se establece el vector K
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Se resuelve el problema mediante el método de Runge-Kutta de orden 4 aplicado a sistemas
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      %Definicion de los parametros K1,K2,K3 y K4 necesarios para la resolución
      K1=w(t(i),f(:,i));
      t(i)=t(i)+h/2;
      z(:,i)=f(:,i)+K1*h/2;
      K2=w(t(i),z(:,i));
      v(:,i)=f(:,i)+K2*h/2;
      K3=w(t(i),v(:,i));
      c(:,i)=f(:,i)+K3*h;
      K4=w(t(i)-h/2,c(:,i));
      f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);
end
   Y(j)=f(2,end);
end
%Representación gráfica de la primera derivada del factor de forma para n=20
hold on
plot(k,Y,'r')
legend('f´(20)')
xlabel('k')
ylabel('f´(20)')
title('Representación de f(20)´ para los posibles k')
hold off
%Valor de k para el cual f'(20) se acerca mas a 1
for q=1:G
    a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
      k(q)
      break
end


El parámetro de k para el cual f'(20) se acerca mas a 1 es k=0.33

2.3 Resolución mediante Euler.

En este tercer apartado se resuelve la ecuación de Blasius en las mismas condiciones, pero en este caso mediante el método numérico de Euler.

derecha
%Se definen los datos del problema
%En este caso la variable "t" es la variable de similaridad
%Numero de ecuaciones
ne=3;
%Se define el intervalo temporal
t0=0;
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Se establece el vector K
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Se resuelve el problema mediante el método de Euler aplicado a sistemas
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      f(:,i+1)=f(:,i)+h*w(t(i),f(:,i));
end
   Y(j)=f(2,end);
end
%Representación gráfica de la primera derivada del factor de forma para n=20
hold on
plot(k,Y,'g')
legend('f´(20)')
xlabel('k')
ylabel('f´(20)')
title('Representación de f(20)´ para los posibles k')
hold off
%Valor de k para el cual f'(20) se acerca mas a 1
for q=1:G
   a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
   k(q)
   q
   break
end
end


Parametro de k para el cual f'(20) se acerca mas a 1 es k=0.32

2.4 Comparativa y elección de k

Analizando los resultados de los tres métodos utilizados para el estudio de la ecuación de Blasius, y fijándonos en las gráficas y en el valor de K obtenido de cada uno de ellos para que f'(20) sea lo más próximo a la unidad, podemos afirmar que la diferencia entre ellos es muy pequeña. El valor obtenido para los métodos Runge-kutta y Euler modificado es de 0,33 y para el método de Euler es de 0,32. Esto último se debe a una mayor precisión local de los dos primeros métodos numéricos y un mayor error en el de Euler siendo en este caso prácticamente despreciable debido a que "h" (valor de salto) es muy pequeño.

El método que vamos a elegir para determinar el valor de k es ek de Runge-kutta de orden 4, ya que el error de paso es de orden O(h5) mientras que el error total acumulado es O(h4), razón por la cual es un método muy preciso y además el más usado en los desarrollos computacionales. El inconveniente a la hora de usar este método es una mayor cantidad de evaluaciones de la función, resultando un mayor tiempo de cálculo y una mayor solicitación de la computadora.

3 Interpretación de la derivada de la función de forma f(η)

En este caso, vamos a estudiar el valor de la variable de similaridad "η" a partir del cual la diferencia entre cualquier valor de la primera derivada de la función de forma y 1 es menor de 0,01. El método numérico escogido para la resolución de este problema es el de Runge-kutta debido a las propiedades mencionadas en el anterior apartado y la k obtenida es 0,33. Como se puede apreciar en el gráfico, el valor de η para el cual la función f'(η) se acerca a menos de 0,01 de 1, que es su límite cuando la variable de similaridad tiende a infinito, es de 5,225 en nuestro caso.
Realizando otra interpretación en la búsqueda de este valor podemos encontrarnos con el concepto de grosor de capa, siendo éste el valor para el cual la primera derivada del factor de forma llega al 99% de su límite.

derecha
%Se definen los datos del problema
%En este caso la variable "t" es la variable de similaridad
%Numero de ecuaciones
ne=3;
%Se define el intervalo temporal
t0=0;
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Determinacion del vector t
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Se resuelve el problema mediante el método de Runge-Kutta de orden 4 aplicado a sistemas para k=0.33
f0=[0 0 k(24)];
f(:,1)=f0;
for i=1:N
   K1=w(t(i),f(:,i));
   t(i)=t(i)+h/2;
   z(:,i)=f(:,i)+K1*h/2;
   K2=w(t(i),z(:,i));
   v(:,i)=f(:,i)+K2*h/2;
   K3=w(t(i),v(:,i));
   c(:,i)=f(:,i)+K3*h;
   K4=w(t(i)-h/2,c(:,i));
   f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);
end
O=ones(1,N+1);
%Representación de  f´ para k=0,33
hold on
plot(t,f(2,:))
legend('f´(n)')
xlabel('n')
ylabel('f´(n)')
title('Representación de  f´ para k=0,33')
hold off
%Determinación de n para que la diferencia entre f' y 1 sea menor de 0.01
for j=1:N+1   
   if abs(f(2,j)-O)<0.01
       t(j)
       break
   end
end
for m=1:N+1
    b(m,:)=abs(f(2,m)-O);
end
%Representacion de la diferencia en funcion de n
figure(2)
hold on
plot(t,b)
xlabel('n')
ylabel('Modulo de la diferencia')
title('Representacion de la diferencia en funcion de n')
hold off


4 Estudio de la velocidad en función del factor de forma f(η)


centro


4.1 Resolución numérica y representación gráfica de la componente U1 del campo de velocidades

Sabiendo que la primera componente u1 del campo de velocidades del fluido es proporcional a la primera derivada del factor de forma, siendo el factor de proporcionalidad la velocidad constante u0=2, podemos en este caso dar utilidad al método numérico de Euler modificado usado en el apartado 2.1 para la representación gráfica de u1. Elegimos este método porque el resultado obtenido para la k es el mismo que en el de Runge-kutta, es decir nos aporta buena precisión, a la vez que reduce carga de trabajo a la computadora.

derecha
%Se definen los datos del problema
%En este caso la variable "t" es la variable de similaridad
%Numero de ecuacionesne=3;
h=0.05;
%Vectores x e y 
x=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yn=3;
y=y0:h:yn;
P=length(y);
N=round((yn-y0)/h);
%Definicion de la función w siendo esta la parametrizacion del sistema en forma general
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Se resuelve el problema mediante el método de Euler modificado aplicado a sistemas
for j=1:length(x)
   H2=h*sqrt(2/x(j));
   f0=[0 0 0.33];
   f(:,1)=f0;
   for i=1:N
   t(i)=y(i)*sqrt(2/x(j))+H2/2;
   K1=w(t(i),f(:,i));
   q(:,i)=f(:,i)+K1*H2/2;
   f(:,i+1)=f(:,i)+H2*w(t(i),q(:,i));
   end
   Z(j,:)=f(2,:);
   %Representación de la primera variable del campo de velocidades para las diferentes x
   hold on 
   plot(y,2*Z(j,:))
end
legend('x=0.05','x=0.2','x=0.4','x=0.6','x=0.8')
xlabel('y')
ylabel('U1')
title ('Representacion de U1')
hold off


4.2 Observaciones

Interpretando las diferentes gráficas para los valores de x indicados, se observa que, cuanto mayor es el valor de la componente horizontal, mayor es la altura que necesita el fluido para alcanzar la velocidad constante u0=2, lo que nos sugiere la idea de que mayor es el grosor de la capa límite.

5 Representación gráfica e interpretación del borde de la capa limite

En esta última ilustración se observa gráficamente el borde de la capa límite exactamente igual como si se viera en un laboratorio realizando el experimento físicamente con los datos indicados, dando a entender el comportamiento real del fluido en esta situación. Por debajo del borde el fluido comienza a reducir su velocidad hasta llegar a ser nula justo en la superficie de la placa, análogamente, a medida que nos alejamos del origen de ella, el fluido tiende a volver a ganar la velocidad constante anterior al contacto.

centro


derecha 

%Se definen los datos del problema
%Vector x (valor de la longitud de la chapa)
x0=0;
xn=10;
h=0.05;
x=x0:h:xn;
%Valor de la variable de similaridad en el borde de la capa
n=5.225;
for i=1:length(x)
    y(i)=n*sqrt(x(i)/2);
end
%Representación de la capa de transición
hold on
plot(x,y);
legend('g(x)')
xlabel('X')
ylabel('Y')
title ('Borde de la capa de transición')
hold off
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Usuario:Grupo_4&oldid=34719»