Usuario:Grupo 4

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. 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₀=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 "η":

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

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

El campo de velocidades del fluido viene dado por:

A continuación resolveremos la Ecuación de Blasius con las condiciones iniciales dadas mediante el uso de métodos numéricos

2 Resolución numérica del problema.

2.1 Resolución mediante Euler modificado.

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
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);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
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
plot(k,Y)
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

izquierda

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

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
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);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
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;
      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
plot(k,Y)
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

2.3 Resolución mediante Euler.

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
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);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
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
plot(k,Y)
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

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

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
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);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
   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);
    plot(t,f(2,:))
 for j=1:N+1   
   if abs(f(2,j)-O)<0.01
       t(j)
       j
       break
   end
end

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

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

clc
clear all
 %Datos iniciales
ne=3;
h=0.05;
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);
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
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,:);
   hold on 
   plot(y,2*Z(j,:))
   hold off
end

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

clc
clear all
x0=0;
xn=10;
h=0.05;
x=x0:h:xn;
eta=5.225;
for i=1:length(x)
    y(i)=eta*sqrt(x(i)/2);
end
plot(x,y);