Capa límite de un fluido laminar sobre una placa plana - (Grupo 17B)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Capa límite de fluido laminar. Grupo b-17
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Eduardo Pereiras Navas

Álvaro Ramírez Fernández de la Puente

Jesús Infestas Robles

Miguel Fandiño Álvarez

Rachel L’Hostis

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


Se estudiará en lo siguiente el comportamiento de un fluido laminar cuando se introduce una placa, según su sección transversal. Para poder conocer que sucede, debemos concretar la ecuación de corriente, que a su vez depende de una función que es solución de una EDO de tercer grado. La ecuación del campo de velocidades del fluido es: [math]\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})[/math] Se supone que las velocidades toman valores constantes en las capas alejadas de la placa, donde [math]\overrightarrow{u} = u_0\cdot\overrightarrow{i}[/math]

La corriente de fluido [math]\psi[/math] viene dada por la función: [math]\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)[/math] Donde [math] f(\eta) [/math] con [math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math] satisface la Ecuación de Blasius: [math]y''' + \frac{1}{2}y y'' =0[/math]

Supondemos la sección de una paca plana definida en [math] (0,x) [/math], con [math]x\gt0[/math] e [math]y\gt0[/math] La complejidad del ejercicio viene debido a resolver numéricamente esta función, llamada Ecuación de Blasius, donde tomaremos como valores iniciales: [math]f(0)=f(0)'=0; f(0)''=k[/math], [math] f'\rightarrow 1[/math] cuando [math]\eta\rightarrow \infty[/math] , [math] u_0 = 2 [/math], la viscosidad [math] \nu [/math] toma valor 1. Donde la [math]k=0.1,0.11,0.12,…1[/math]

1 Solución numérica de la ecuación de blasius

Para la solución de la ecuación diferencial utilizaremos las condiciones iniciales anteriormente dadas.

Fluido laminar sobre sobre la sección de una placa


[math] y(\nu) [/math] satisface la ecuación de Blasius, que para resolverla numéricamente se tiene que transformar en un sistema de tres ecuaciones diferenciales

[math] y(\eta)=y_1,y’(\eta)=y_2,y’’(\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]


1.1 Resolución del Sistema con el Método de Euler Modificado

A continuación se añade la solución numérica que nos ha dado el programa Matlab 

%primero definimos las variables e inicializamos las matrices
k=[0.1:0.01:1];
l=length(k);
t0=0;
tf=20;
h=0.05;
N=(tf-t0)/h;
t=[t0:h:tf];
Y1=zeros(l,N+1);
Y2=zeros(l,N+1);
Y3=zeros(l,N+1);
yy=zeros(3,l);
%comenzamos el bucle para todas las k de los valores iniciales
%dentro de este, ponemos otro bucle que los resuelva la ecuacion diferencial. Metodo de Euler
for j=1:l
  y0=[0; 0; k(j)];
  yy=y0;
  for n=1:N
    A=[0 ,1, 0;0 0 1; 0 0 -0.5*yy(1,n)];
    K1=A*yy(:,n);
    K2=A*[yy(:,n)+h*K1];
    yy(:,n+1)=yy(:,n)+0.5*h*(K1+K2);
  end
  Y1(j,:)=yy(1,:);
  Y2(j,:)=yy(2,:);
  Y3(j,:)=yy(3,:);
end
%las soluciones y1, y2, y3, se guardan en matrices donde cada fila son los valores (y1, y2, o y3) para una k concreta.
%Buscaremos para todos los valores de k los diferentes valores de f(20)
  fy=Y2(:,l);
  plot(fy,k,'*')
%¿Cual es el valor del parametro para el cual f’(20) se acerca más a 1?
%Para obetener el valor más proximo a uno, utilizamos un bucle que nos de la posición del vector k tal que la diferencia entre f'(y)-1 sea mínima
for m=1:l
ll=m;
q=abs(1-fy(m));
if min(abs(1-fy))==q
break
end
end
display(ll)
kproxima=k(ll)
Gráfica [math]f’(20)[/math] para cada [math]k[/math]

El valor de kpoxima, que será el [math]k[/math] tal que la distancia a [math]f'(20)=1[/math]sea mínima será [math]k=0.31[/math], que en el vector [math]k[/math] estará en la posición 22.

1.2 Resolución mediante Runge-Kutta 4

Procederemos de la misma forma, cambiando aquellos puntos característicos de Runge-Kutta 4 

k=[0.1:0.01:1];
l=length(k); 
t0=0;
tf=20;
h=0.05;
N=(tf-t0)/h;
t=[t0:h:tf];
Y1=zeros(l,N+1);
Y2=zeros(l,N+1);
Y3=zeros(l,N+1);
yy=zeros(3,l);
%comenzamos el bucle para todas las k de los valores iniciales
%dentro de este, ponemos otro bucle que los resuelva la ecuacion diferencial. Metodo de Runge-cuta
for j=1:l
  y0=[0; 0; k(j)];
  yy(:,1)=y0;
  for n=1:N
    A=[0 ,1, 0;0 0 1; 0 0 -0.5*yy(1,n)];
    K1=A*yy(:,n);
    K2=A*[yy(:,n)+0.5*h*K1];
    K3=A*[yy(:,n)+0.5*h*K2];
    K4=A*[yy(:,n)+h*K3];
    yy(:,n+1)=yy(:,n)+0.5*h*(K1+2*K2+2*K3+K4);
  end
  Y1(j,:)=yy(1,:);
  Y2(j,:)=yy(2,:);
  Y3(j,:)=yy(3,:);
end
%las soluciones y1, y2, y3, se guardan en matrices donde cada fila son los valores (y1, y2, o y3) para una k concreta.
%gráficaf'(20) para cada valor del parametro
  fy=Y2(:,21);
  plot(fy,k,'*')
%¿Cual es el valor del parámetro para el cual f0(20) se acerca mas a 1?
  m=max(fy)


1.3 Resolución mediante el Método de Euler

De la misma forma aplicamos: 

k=[0.1:0.01:1];
l=length(k);
t0=0;
tf=20;
h=0.05;
N=(tf-t0)/h;
t=[t0:h:tf];
Y1=zeros(l,N+1);
Y2=zeros(l,N+1);
Y3=zeros(l,N+1);
yy=zeros(3,l);
%comenzamos el bucle para todas las k de los valores iniciales
%dentro de este, ponemos otro bucle que los resuelva la ecuacion diferencial. Metodo de Euler
for j=1:l
  y0=[0; 0; k(j)];
  yy(:,1)=y0;
  for n=1:N
    A=[0 ,1, 0;0 0 1; 0 0 -0.5*yy(1,n)];
    yy(:,n+1)=yy(:,n)+h*A*yy(:,n);
  end
  Y1(j,:)=yy(1,:);
  Y2(j,:)=yy(2,:);
  Y3(j,:)=yy(3,:);
end
  fy=Y2(:,21);
  plot(fy,k,'*')


1.4 Conclusiones

Utilizando el método de Euler modificado y el de Runge Kutta 4, obtenemos valores de k muy similares, que además son muy próximos al valor real por la precisión de estos métodos. Por otra parte en el de Euler sale otro valor diferente alejado.

2 Gráfica [math]f´(\eta)[/math] para un valor [math]k[/math] determinado

Nuestro valor será [math]k=0.31[/math]: 

%dibujamos para k=0.31
yy2=Y2(22,:);
plot(t,yy2,'*')
xlabel('\eta')
ylabel('y´(\eta)')
Gráfica [math]f´(\eta)[/math] para [math]k=0.31[/math]

Para encontrar el valor de [math]\eta_0[/math] por el cual [math]\ \vert f'(\eta)-1 \vert \lt 0,01[/math] , si [math]\eta\gt\eta _0[/math], es [math]\eta_0\ge4.5[/math] para [math] \eta \in \mbox{(0,20)}[/math]. 

g=abs(1-fy);
w=1;
for s=1:l
 if g(s)<0.99
  u(w)=g(s);
  w=w+1;
  end
end

tam=length(u);
ultimo=u(tam);
primero=u(1);

for q=1:l
  if g(q)=primero
    inicio=q;
    end
  if g(q)=ultimo
    final=q;
    end
  end
eta0=t(final)


3 Componente horizontal de la velocidad

En este apartado, vamos a calcular analiticamente la componente horizontal del campo de velocidades del fluido, u1, en funcion de f(n).

Vamos a calcular dicha componente derivando y utilizando la regla de la cadena: [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} [/math] Sabiendo que la función [math] \eta(x,y) es: \eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math] Nos queda la ecuacion diferencial: [math]u_1=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)[/math] Ahora procedemos a resolverla por el metodo Runge-Kutta 4 para los valores [math]y \in (0,3)[/math] y [math]x \in \{0.05, 0.2, 0.4, 0.6, 0.8\}[/math]. Este programa nos permite sacar la grafica de los valores de [math]u_1(x_k,y)[/math] para dichos valores. 

x=[0.05,0.2,0.4,0.6,0.8]; l=length(x);
v=1; u0=2; k=0.31;
y0=0; yf=3; h=0.01;
N=(yf-y0)/h;
y=y0:h:yf;
f0=[0;0;k];
f1=f0;
f=zeros(3,N+1);
f(:,1)=f0;
for i=1:l
	t(1,:)=sqrt(u0/(v*x(i)))*y;
	H=sqrt(u0/(v*x(i)))*h;
	for n=1:N
		k1=[f1(2);f1(3);-1/2*f1(1)*f1(3)];
		t1=t(n)+H/2; 
   		ff=f1+1/2*H*k1; 
   		k2=[ff(2);ff(3);-1/2*ff(1)*ff(3)];
   		ff=f1+1/2*H*k2;
        	k3=[ff(2); ff(3); -1/2*ff(1)*ff(3)];
        	t1=t(n)+H;
   		ff=f1+H*k3; 
   		k4=[ff(2); ff(3); -1/2*ff(1)*ff(3)];
   		f1=f1+H*(k1+2*k2+2*k3+k4)/6;
   		f(:,n+1)=f1;
	end
	F(i,:)=f(2,:);
end
u1=u0.*F;
hold on
plot(y,u1(1,:))
plot(y,u1(2,:),'m')
plot(y,u1(3,:),'b')
plot(y,u1(4,:),'r')
plot(y,u1(5,:),'g')
legend('x=0.05','x=0.2','x=0.4','x=0.6','x=0.8')
xlabel('y')
ylabel('u1')
title ('Componente u1 de la velocidad')
hold off
Gráfica para cada [math]k[/math]

3.1 Conclusión

Si nos fijamos en la gráfica, observamos que el movimiento del fluido tiene una velocidad limite mas estable cuando tiene una mayor altura respecto la placa plana o eje de abscisas. Sin embargo, para los primeros valores de [math]x_k[/math] la gráfica tiene la forma parecida a una recta y para los mayores es muy próxima a una parábola.

4 Borde de la capa límite del fluido

Utilizando los resultados obtenidos y las condiciones iniciales, deducimos que existen unos valores de posición[math]x[/math] e [math]y[/math] por la cual la velocidad toma un valor estacionario. El valor de la capa este será la constante [math]\eta_0[/math], calculada anteriormente [math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math]

Si [math]\eta=\eta_0[/math],entonces

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

Podemos dibujar esta gráfica (ya que [math]y[/math] no toma un valor constante) mediante: 

x=[0:0.02:10];
eta0=4.5;
y=eta0*(x/2).^(1/2);
plot(x,y,'*')
xlabel('x')
ylabel('g(x)')
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Capa_límite_de_un_fluido_laminar_sobre_una_placa_plana_-_(Grupo_17B)&oldid=10659»