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

Trabajo realizado por estudiantes
Título	Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	María García Fernández José Manuel Corbí Garrido Estefanía Jiménez Ocampo Diego García Vaquero José Francisco Aguilera Corral
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción y enunciado del problema

En este artículo, se desarrollará el estudio de la capa límite de un fluido laminar sobre una placa plana semi-infinita mediante métodos numéricos programados en lenguaje M. De esta forma el lector podrá reproducir (si así lo desea) los algoritmos y códigos aquí escritos en MATLAB u otras alternativas de software libre como Octave.

Fluido laminar por encima de una placa plana

El tratamiento de este problema mediante métodos numéricos requiere una definición precisa de las condiciones iniciales y de contorno. En nuestro caso definiremos la semi-placa en el plano cartesiano como la semirecta [math]y=0[/math], con [math]x∈(0,∞)[/math]; y estudiaremos lo que pasa por encima de ella [math](x\gt0, y\gt0)[/math]. El fluido laminar que se encuentra con la placa, antes de llegar a ella tiene una velocidad constante de [math]\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 [/math], iniciándose una fase de transición en su encuentro con la placa (la velocidad es nula en los puntos de la placa), y lejos del punto de encuentro, una vez finalizada la transición, volviendo a fluir a la velocidad constante inicial. El campo de velocidades del fluido se define como:: [math]\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})[/math] donde [math]\psi[/math] es la corriente del fluido:: [math]\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)[/math] y [math]ν[/math] es su viscosidad (que supondremos [math]ν = 1[/math]), [math] f(\eta) [/math] con [math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math] satisface la ecuación de Blasius::

[math]f''' + \frac{1}{2}f f'' =0[/math]

En los siguientes apartados resolveremos numéricamente esta ecuación diferencial ordinaria no lineal, para unas condiciones iniciales:: [math]f(0)=f(0)'=0; f(0)''=k[/math] [math] f'\rightarrow 1[/math] cuando [math]\eta\rightarrow \infty[/math] tomando [math]k[/math] valores [math] k \in \mbox{(0,1;1)}[/math] y [math]dk=0,01[/math] A continuación estudiaremos la resolución de la ecuación mediante diferentes métodos numéricos, y estudiaremos el [math]k[/math] "óptimo" para las condiciones iniciales enunciadas.

2 Resolución numérica de la ecuación de Blasius

Para poder resolver esta ecuación independientemente del posterior procedimiento numérico que vayamos a aplicar, primero hemos tenido que convertirla en un sistema, aplicando el método analítico que detallamos a continuación::

[math]y_{1}=f(x)[/math]:

[math]y_{2}=f'(x)[/math]:

[math]y_{3}=f''(x)[/math]

Así [math]y_{3}'=f'''(x)[/math], si despejamos de nuestra ecuación diferencial obtenemos que::

[math]y_{2}=y_{1}'[/math]:

[math]y_{3}=y_{2}'[/math]:

[math]y_{3}'=-\frac{1}{2}y_{1}y_{3}[/math]

Como resultado de esta transformación obtenemos un sistema para el cual podemos obtener su solución numéricamente.

2.1 Mediante el método de Euler modificado

A continuación se muestra el código en lenguaje matlab que proporciona las soluciones numéricas a la ecuación de Blasius para los diferentes valores de k (recordemos que tomaba valores de 0,1 hasta 1 con un paso de 0,01) en un intervalo [math] \eta \in \mbox{(0,20)}[/math]. Hemos tomado una [math]h=0,5[/math].

Implementando el método de Euler en este sistema:

%% Método de Euler modificado
 
% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
%% Datos del problema y discretización
x0=0; % las x son las eta de las que depende f
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta, k y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K  
    f0=[0 0 k(j)]';
    %% Inicialización
    f(:,1)=f0;
    ff=f0; % Valor de la aprox. en cada x
    %% Iteraciones y(x_n)-->y(x_n+1)
    for n=1:N
        % calculamos k1=f(x_n,y_n)
        k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
        % Calculamos k2=f(x_n+h/2,y_n+1/2k1*h)
        xp=x(n)+h; %argumento de x
        fp=ff+h*k1; %argumento y 
        k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
         
        ff=ff+h/2*(k1+k2);
        f(:,n+1)=ff;
 
    end
    laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solución
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, siendo k el que mejor aproxima f´(20) a 1')
corte=1;
figure (2)
hold on 
plot(k,laf2paraescogerk)
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')
hold off
%Para calcular la k de forma exacta utilizamos este bucle que nos dice con qué k la diferencia en valor absoluto entre la f´(n) y uno es menor
for n=1:K
     l=n;
     p=abs(1-f2(n));
 if min(abs(1-f2))==p 
     
   break
 end
end
kfinal=k(l)

f,f´,f´´ con k=0,33 para intervalo de eta de 0 a 20

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es [math]k=0.330[/math]

2.2 Mediante el método de Runge-Kutta de cuarto orden

La siguiente secuencia de código da solución a la ecuación de Blasius mediante el método de Runge-Kutta, siendo aplicado al sistema de ecuaciones diferenciales obtenido en el primer apartado mediante procedimientos analíticos:

% Resolver sistema  con un método RK4
%% Método de Runge-Kutta orden 4

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% La variable eta recibe el nombre x
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K  
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada tiempo

for n=1:N
    % calculamos k1=f(x_n,f_n)
    k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
    % Calculamos k2=f(x_n+h/2,f_n+1/2k1*h)
    xp=x(n)+h/2; %argumento de eta
    fp=ff+1/2*h*k1; %argumento f 
    k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
    %Calculamos k3=f(x_n+h/2,f_n+1/2k2*h)
    %mismo argumento eta que antes
    fp=ff+1/2*h*k2;%argumento y
    k3=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
    %calculamos k4=f(x_n+h/2,f_n+1/2k3*h)
    xp=x(n)+h;%argumento eta
    fp=ff+h*k3; %argumento f
    k4=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
    ff=ff+h/6*(k1+2*k2+2*k3+k4);
    f(:,n+1)=ff;
  end
 laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f1(x)','f2(x)','f3(x)')
corte=1;
figure (2)
plot(k,laf2paraescogerk)
hold on 
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
     l=n;
     p=abs(1-f2(n));
 if min(abs(1-f2))==p 
     
   break
 end
end

kfinal=k(l)

Según este método la k que mejor se ajusta a las condiciones iniciales es k=0.330, como podemos apreciar en el gráfico.

2.3 Mediante el método de Euler

Al igual que en los otros apartados resolvemos el sistema asociado a la ecuación de Blasius numéricamente, en este caso mediante el método de Euler:

%% Método de Euler

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% x es la variable eta del problema
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K  
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada valor de eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N
    
   y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
   ff=ff+h*y;
   f(:,n+1)=ff;
   
end
 laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
clf
corte=1;
figure (1)
plot(k,laf2paraescogerk)
hold on 
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')


for n=1:K
     l=n;
     p=abs(1-f2(n));
 if min(abs(1-f2))==p 
     
   break
 end
end

kfinal=k(l)
figure(2)
%ya con la k=0,33

x0=0;
xf=20;
k=0.32   
f0=[0 0 k]';
%% Datos discretizacion
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
f=zeros(3,N+1);
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
   ff=ff+h*y;
    f(:,n+1)=ff;
   
end
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, con Euler f´(20) a 1, ya con la k escogida')

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es k=0.320

2.4 Comparativa entre los tres métodos, elección de k

La diferencia entre los tres métodos numéricos empleados para resolver la ecuacion de Blasius (Euler modificado, Runge-Kutta y Euler), es muy pequeña. La forma de las gráficas es completamente similar y el valor que se obtiene para que f'(20) sea lo más proxima a uno es en el caso de Euler modificado y Runge-Kutta igual a 0.33 y en Euler es 0.32. Por la tanto la diferencia entre los tres métodos es pequeña.

De hecho a simple vista no se aprecia la diferencia.

De cara a obtener una mejor precisión con un paso "grande" es aconsejable el método de Runge-Kutta, pues con el método de Euler se comete un mayor error local frente al método de Runge-Kutta, aunque la ejecución de este último supone una mayor carga de trabajo para el ordenador. Para compensar este mayor error local en el método de Euler (diferencia entre la solución exacta para un punto y la estimada) se puede reducir el paso, aunque con ello el programa se ejecutará mas lentamente. En nuestro caso, para el paso considerado en los códigos, será más fiable el método de Runge-Kutta de cara a la obtención del valor de k, por ello, en lo sucesivo, consideraremos k=0.33

3 Representación e interpretación de [math]f’(\eta)[/math]

Vamos a estudiar a partir de qué valor de [math]\eta[/math] la gráfica difiere 0,01 o menos de 1.

Para esto se puede escoger cualquiera de los métodos numéricos anteriormente utilizados. Y después, para que nos muestre la gráfica de [math]f’(\eta)[/math], se escriben las siguientes líneas de código:

plot(eta,f(2,:),'g')
legend('f2(eta)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de  f´ para k=0,33')

Observamos como tiende rápidamente al valor 1, por lo que confirma que la k escogida en los apartados anteriores es correcta.

Ahora para encontrar a partir de qué valor de [math]\eta[/math] la gráfica difiere 0,01 o menos de 1; introducimos la siguiente secuencia de código:

figure(2)
gg=abs(f(2,:)-1);
for n=1:N
     l=n;
 if gg(n)<0.01
   break
 end
end
nu0=l-1
loquequermos=eta0*h

plot(eta0,gg)

Visualmente vemos que el valor buscado está más o menos entre 5 y 6. Con la condición del código de que la diferencia en valor absoluto sea menor que 0.01 llegamos a que el valor que buscamos concretamente es 5,2.

4 Campo de velocidades del fluido

Definimos el campo de velocidades del fluido [math]\vec{u}[/math] como::

[math]\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})[/math]

4.1 Cálculo y representación de [math]u_1[/math]

A continuación vamos a calcular analíticamente la componente de [math] u_1 [/math] de la velocidad del fluido en función de [math]f( \eta)[/math]

Para ello, desarrollamos la componente hasta la obtención de su expresión final::

[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]:

Definimos la función [math]\eta(x,y)[/math] como::

[math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math]:

[math]u_1=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)[/math]

Resolviendo esta ecuación para un [math]y \in (0,3)[/math] y unos valores de [math]x \in \{0.05, 0.2, 0.4, 0.6, 0.8\}[/math] mediante el siguiente código, que nos permite calcular las gráficas para cada una de las [math]x_{i}[/math] particulares:

% Componente u1 de la velocidad
% Introducimos las variables
clear all
u0=2;
nu=1;
xk=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yf=3;
h=0.05;
y=y0:h:yf;
N=(yf-y0)/h;
% Iteraciones
for i=1:5
eta(i,:)=y*sqrt(u0/(nu*xk(i)));
% H es el paso para resolver el sistema
H=0.05*sqrt(u0/(nu*xk(i)));
k=0.33;
f0=[0 0 k]';
ff=f0;
for n=1:N
    A=[0 1 0;0 0 1;(-ff(3)/2) 0 0];
    j=ff+H*A*ff;
    ff=ff+(H/2)*(A*ff+A*j);
    f(:,n+1)=[ff(1);ff(2);ff(3)];
end
G(i,:)=f(2,:);
end
u1=u0*G;

%% Dibujamos la solución
hold on
plot(y,u1)
legend('para x=0.05','para x=0.2','para x=0.4','para x=0.6','para x=0.8')
xlabel('y')
ylabel('u1')
title ('Componente u1 de la velocidad')
hold off

4.2 Conclusión

Al representar únicamente la componente [math]u_1[/math] del campo de velocidades, no tenemos una percepción global de dicho campo vectorial, si no únicamente de la variación de la corriente con respecto a la altura, es decir de las envolventes (tangentes) a esta componente de la velocidad [math]u_1[/math](definida como la derivada parcial de la corriente con respecto a y), que resultan ser una familia de líneas de corriente.

Estas líneas de corriente, cuya tangente en un punto cualquiera tiene la dirección de la velocidad, si tenemos en cuenta que nuestro fluido es laminar, realmente representarán las trayectorias descritas por las partículas del fluido. Observamos como estas líneas de corriente tienden a estabilizarse a medida que aumenta la altura respecto a la placa, dicho de otra forma, a medida que se alejan de ella. Análogamente si representásemos las líneas de corriente para la componente [math]u_2[/math] observaríamos como se estabilizan a medida que aumenta la distancia al origen de la placa (lugar en el que comienza una fase de transición).

Las partículas del fluido no hacen otra cosa si no evitar el "obstáculo" que hemos colocado, la placa, tratando de modificar lo menos posible su trayectoria rectilínea a velocidad constante [math]u_0[/math]. Por ello, a medida que nos alejamos del obstáculo (la placa), bien según las abscisas o las ordenadas, las partículas del fluido tienden a recuperar su trayectoria rectilínea con velocidad constante, mientras que en la zona de transición las líneas de corriente tienden a "envolverlo".

Si el lector lo desea, puede profundizar en el tema consultando los siguientes enlaces: [1] Artículo elaborado por el grupo 4A para la asignatura Teoría de Campos durante el curso 2013-2014 [2] [3]

5 Capa límite, transición entre la velocidad nula al comienzo de la placa y lejos de la misma

Hemos comprobado que la altura ([math]y[/math]) a la que el fluido alcanza la velocidad constante ([math]u_0=2[/math]) depende de la distancia recorrida sobre la placa ([math]x[/math]). Ahora vamos a dibujar la función [math]g(x)[/math], donde [math]g(x)[/math] es el valor [math]y[/math] para el cual [math]\eta=\eta_0[/math], para observar como es ese cambio de [math]y[/math] en función de la [math]x[/math]. La relación entre estas variables viene dada por la expresión: [math]\eta = y \sqrt[]{ \frac{u_0}{\nu x}}[/math] Despejamos la [math]y[/math] y suponemos [math]u_0=2[/math], [math]\nu=1[/math] y [math]\eta=\eta_0=5.2[/math]

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

Representamos esta función en Matlab.

% Dibujamos la función g(x), y en función eta0(cte.) y  de x.
% Introducimos variables
eta0=5.2;
u0=2;
v=1;
x=0:0.05:10;
% Calculamos y
y=eta0*sqrt(v*x/u0);
% Dibujamos la gráfica
plot(x,y)
legend('g(x)')
xlabel('x')
ylabel('y')
title('Capa límite de transición de la velocidad')

Interpretamos la función [math]g(x)[/math] como el borde de la capa límite que hace la transición entre la velocidad nula en la placa y la velocidad lejos de la misma del fluido. Se confirma lo que ya podíamos ver en el apartado anterior. A medida que el fluido se aleja del origen de la placa ([math]x=0[/math]), es decir, a medida que aumenta el valor de la variable [math]x[/math], el fluido alcanza la velocidad inicial constante a una mayor altura ([math]y[/math]).