Boundary layer in laminar fluids (Grupo 1B)

Trabajo realizado por estudiantes
Título	Boundary layer in laminar fluids. Grupo 1-B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Sandro Andrés Martínez David Ayala Díez Claudia Cózar Coarasa Lorena de la Fuente Sanz Marino Rivera Muñoz José Manuel Torres Serrano
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

In this numerical project we have studied what happens when we introduce a flat plate in a laminar fluid whose speed and viscosity are constant.

1 Condición de frontera tipo Neumann

Ahora vamos a considerar una condición de frontera diferente en el extremo derecho. En vez de suponer una temperatura constante en ese extremo como anteriormente, vamos a colocar en él una pieza aislante. Este aislante provoca que no haya pérdida de calor por el extremo derecho, es decir, que el flujo de temperatura sea nulo. Esta condición es de tipo Neumann, a diferencia de las anteriores que eran de tipo Dirichlet. Así, mantenemos la condición en el extremo izquierdo y aplicamos la nueva en el extremo derecho, que resulta

[math]-ku_x(3,t)=0 \rightarrow u_x(3,t)=0[/math]

En esta situación, la temperatura de la varilla viene dada por el siguiente problema

[math] \begin{cases} u_t-u_{xx}=0, \qquad x\epsilon(0,3), t\gt0\\ u(0,t)=0; u_x(3,t)=0, \qquad t\gt0\\ u(x,0)=u_0(x), \qquad x\epsilon[0,3] \end{cases} [/math]

Por lo tanto, ahora en el estado estacionario para tiempos grandes la función [math]u(x,t)[/math] que modeliza la temperatura de la varilla es solución del siguiente problema de contorno (lo llamamos así pues la ecuación diferencial depende únicamente de [math]x[/math] en el estado estacionario)

[math] \begin{cases} u_{xx}=0, \qquad x\epsilon(0,3), t\gt0\\ u(0,t)=0, u_x(3,t)=0, \qquad t\gt0 \end{cases} [/math]

Resolvemos la ecuación diferencial obteniendo el mismo resultado que anteriormente

[math]u_{xx}=0 \rightarrow u(x,t)=C_1(t)x+C_2(t)[/math]

Aplicando las condiciones de frontera: [math]C_1(t)=C_2(t)=0 \rightarrow u(x,t)=0[/math]

Este resultado refleja que al cabo de un tiempo lo suficientemente grande como para considerar un estado estacionario en la varilla, esta adquiere una temperatura uniforme nula. El comportamiento de la varilla es coherente con las condiciones de frontera, pues la temperatura final de la varilla es la que se mantiene constante en el extremo izquierdo, que permite flujo de temperatura.

1.1 Método de diferencias finitas

En la siguiente imagen se muestra una aproximación del problema mediante el método de diferencias finitas. En ella se puede apreciar como para un valor elevado del tiempo la temperatura en la varilla se puede considerar constante y uniforme en toda ella, alcanzando el valor estacionario [math]u(x,t)=0[/math].

(IMAGEN)

En concreto, a partir de un tiempo [math]t=…[/math] podemos considerar que la temperatura alcanza dicho valor estacionario con un error del 5%, es decir, en ese instante la diferencia entre la distribución térmica calculada y la estacionaria es… (ESTA PARTE HAY QUE COMPLETARLA CUANDO PREGUNTEMOS A CARLOS A QUÉ SE REFIERE EL 5%).

A continuación se refleja el código de Matlab que aproxima la temperatura de la varilla por diferencias finitas usando el método de Euler implícito por ser el que aporta una mejor aproximación.

1.2 Método de Fourier

Planteamos ahora el mismo problema usando el método de Fourier. Por tanto, buscamos soluciones de la forma [math]u(x,t)=\varphi(x)T(t)[/math], donde [math]\varphi(x)[/math] satisface el siguiente problema de autovalores:

[math] \begin{cases} \varphi’’(x)+\lambda\varphi(x)=0 ;\\ \varphi(0)=0; \varphi’(3)=0 \end{cases} [/math]

La solución de la ecuación diferencial es [math]\varphi(x)=a\cos(\sqrt{\lambda}x)+b\sin(\sqrt{\lambda}x)[/math]

Aplicando las condiciones de frontera obtenemos los autovalores [math]\mu_k[/math] y autofunciones [math]\varphi_k(x)[/math] del problema

[math]\varphi(0)=0 \rightarrow a=0[/math]; [math]\varphi’(3)=0 \rightarrow \sqrt{\lambda}b\cos(\sqrt{\lambda}\,3)=0 \rightarrow \sqrt{\lambda}\,3=(k-{1\over2})\pi \rightarrow \lambda=\mu_k=(k-{1\over2})^2{\pi^2\over9}[/math];

Así [math]\boldsymbol{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x} \qquad k=1,2,3…N[/math]

Si aplicamos que [math]u_k(x,t)=\varphi_k(x)T_k(t)[/math] satisface la ecuación diferencial [math]u_t-u{xx}=0[/math], obtenemos la ecuación diferencial que determina [math]T_k(t)[/math]

[math]u_t-u{xx}=0 \rightarrow \varphi_k(x)T_k’(t)-\varphi_k’’(x)T_k(t)=\varphi_k(x)T_k’(t)-(-\mu_k)\varphi_k(x)T_k(t)=0 \rightarrow T_k’(t)+\mu_kT_k(t)=0[/math]

La solución de esta ecuación diferencial es [math]T_k(t)=C_ke^{-\mu_kt}= C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}[/math] y por tanto [math]u_k(x,t)= \varphi_k(x)T_k(t)=C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}\sin(k-{1\over2}){\pi\over3}x[/math]

Si expresamos la solución del problema como [math]u(x,t)=\sum_{k=1}^Nu_k(x,t)= \sum_{k=1}^N C_ke^{-(k-{1\over2})^2{\pi^2\over9}t}\sin(k-{1\over2}){\pi\over3}x[/math] y hacemos que satisfaga la condición inicial obtenemos [math]u(x,0)=\sum_{k=1}^NC_k\sin(k-{1\over2}){\pi\over3}x[/math]. De esta forma, por unicidad de los coeficientes de Fourier, los coeficientes [math]C_k[/math] coinciden con los de la serie de Fourier con respecto a las autofunciones [math]\varphi_k(x)[/math] de la función a trozos que determina la condición inicial (expresada al principio del artículo).

El problema se limita así al cálculo de dichos coeficientes, de acuerdo con la expresión

[math]C_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}[/math]

El grado de exactitud de la aproximación con este método depende del número de elementos de la serie de Fourier, es decir, del valor de [math]k[/math]. Estudiamos la temperatura de la varilla tomando los valores [math]k=1,3,5,10,20[/math], como se puede ver en la siguiente imagen. En ella se aprecia, sobre todo en la condición inicial, que a medida que se aumente el valor de [math]k[/math] la función aproximada se acerca más a la real.

Estos resultados se pueden comparar mejor tomando un instante fijo [math]t=0.5[/math] y representando en él cada función en una misma gráfica, como podemos ver a continuación.

First, we must assume that the fluid velocity before reaching the plate is constant, as in remote areas after passing the plate.

In our case we take this constant as [math]2[/math], such that

[math]\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 [/math]

Laminar fluid with velocity [math] u_0 [/math]

Then, we must define the fluid stream function

[math]\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)[/math]

Where we take the viscosity [math] \nu [/math] as a unit value and

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

[math] f(\eta) [/math] Satisfies the Blasius equation, and therefore we will raise the initial value problem associated with this equation with the following initial conditions

[math] \begin{cases} f’’’(\eta)+\frac{1}{2}f(\eta)f’’(\eta)=0 ; \\ f(0)=f’(0)=0, \lim_{\eta \to \infty}f’(\eta)= 1 ; \end{cases} [/math]

However, on time of programming we cannot introduce a conditional limit, so we have to replace them by the condition [math] f’’(0)=k [/math], and vary the values of k to find the one that satisfies the limit.

We cannot solve the differential equation like such, to apply the numerical methods, we need to pass it to a system of equations:

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

Once the system has been formulated , we start to solve it.

1.3 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] .

%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);
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

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].

1.4 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] .

% 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

Graph of [math]f’(20)[/math] for each [math]k[/math]

1.5 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] .

%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

Graph of [math]f’(20)[/math] for each [math]k[/math]

1.6 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.

1.7 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:

%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)')

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 \lt 0,01[/math] , if [math]\eta\gt\eta _0[/math], is [math]\eta_0\ge5,95[/math] for [math] \eta \in \mbox{(0,20)}[/math].

2 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].

%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

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 causes 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).

3 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

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

Therefore we can interpret this function [math]g(x)[/math] as the fluid boundary layer. Simple Matlab code is exposed for its representation:

%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')

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.