Flujo de Poiseuille (Grupo 23)

Trabajo realizado por estudiantes
Título	Flujo de Poiseuille (Grupo 23)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Ana Gastañaga Solana Jaime Casanova Navas Jorge Muñoz Jiménez Daniel Galarza Polo Óscar García Caballero
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

La ley de Poiseuille es la ley que permite determinar el flujo laminar estacionario de un líquido incompresible. En este caso, vamos a considerar el flujo de una tubería cilíndrica de radio 2, suponiendo que está centrada en el eje OZ.

Conociendo la función por la que viene dada la velocidad [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}},[/math] y la presión [math]p\left(x,y\right)[/math] que viene dada por [math]p\left(x,y\right)=p_{1}+\left (p_{2}-p_{1}\right)(z-1)/{2}[/math] representaremos los gráficos respectivos al comportamiento de la velocidad y la presión de la tubería.

Además veremos qué puntos del campo tienen mayor rotacional y observaremos también los máximos de la temperatura en la tubería, dada por el campo: [math]T\left (\rho,\theta,z\right)=log\left(1+\rho\right)e^\left(-\left(z-2\right)^{2}\right)[/math].

Por último calcularemos el caudal que pasa por la sección transversal de la tubería.

Destacar que hemos utilizado el programa Matlab para la representación del mallado y de los gráficos que aparecen en los distintos apartados del trabajo.

centro

2 Mallado de la sección transversal de la tubería

El mallado dibujado en dimensión 2 representa la sección transversal de la tubería cilíndrica de radio 2 [math] x_{1} = 0 [/math], fijando los ejes en la región [math] \left ( \rho,z \right )\epsilon \left [ 0,3 \right ]\times \left [ 0,10 \right ]. [/math]

centro

x=0:0.1:3;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y)                                                    
hold on
mesh(xx,yy,0*xx)
axis([0,3,0,10])
view(2);
title('Mallado de la sección de la tubería');
hold off

3 Ecuación Navier-Stokes

La velocidad esta definida por el campo vectorial [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}[/math] y su presión por el camp escalar [math]p\left(x,y\right)=p_{1}+\left (p_{2}-p_{1}\right)(z-1) [/math] siendo [math] p_{1} [/math] es la presión en [math] z=1 [/math], [math] p_{2} [/math] la presión en [math] z=3 [/math] y [math] \mu [/math] su coeficiente de viscosidad.

Ecuacion de Navier-Stokes: [math] \left ( \vec{u}\cdot \triangledown \right )\vec{u}+\triangledown p=\mu\Delta\vec{u}, [/math]

Si [math] \left ( \vec{u},\rho \right ) [/math] satisface la ecuación Navier-Stokes y despreciamos el termino convectivo de la misma llegamos a la siguiente ecuacion diferencial

[math] \frac{1}{\rho} \frac{\partial }{\partial \rho }\left ( \rho \frac{\partial f\left ( \rho \right )}{\partial \rho } \right )=\frac{p_{2}-p_{1}}{\mu} [/math].

Mutiplicando por [math] \rho [/math] e integramos dos veces obtenemos que

[math] f\left(\rho \right )=\frac{\left ( p_{2}-p_{1} \right )}{4\mu}\rho^{2} + \rho*c + k[/math].

Para encontrar la solución exacta suponemos que para [math] \rho=2 [/math] y para [math] \rho=2 [/math], [math] f\left ( \rho \right )=0 [/math],obteniendo:

[math]f\left ( \rho \right )=\frac{1}{\mu }\left [ \frac{p2-p1}{4 }\rho ^{2} +p1-p2\right ][/math]

Ademas, al ser la funcion de velocidad [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}[/math] y [math] f\left ( \rho \right ) [/math] no depender de [math] z [/math] su divergente es nulo lo que nos indica que el fluido es incompresible.

4 Representación del campo de presiones y velocidades

Suponiendo que [math] p_{1}=4, p_{2}=1 [/math] y [math]\mu=1 [/math], primeramente calculamos el campo de presiones y de velocidades.

- Calculamos el campo de presiones introduciendo los datos proporcionados en la ecuación de la presión [math]p\left ( x,y \right ) [/math]: [math]p\left ( x,y \right )=p_{1}+\left ( p_{2}-p_{1} \right )\left ( z-1 \right ){/2}=4+\left ( 1-4 \right )\left ( z-1 \right ){/2}={4}-\frac{3z}{2}+\frac{3}{2}=\frac{-3z}{2}+\frac{11}{2} [/math]

- Calculamos ahora el campo de velocidades utilizando la ecuación de la velocidad dada: [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}} [/math], siendo [math] f\left(\rho\right)=\frac{1}{\mu}\left(\frac{\left(p_{2}-p_{1}\right)}{4}\rho^{2}+p_{1}-p_{2}\right)\vec{e_{z}}[/math],(como hemos calculado previamente mediante la ecuación de Navier-Stockes). Sustituyendo los datos en la ecuación obtenemos: [math]f\left(\rho\right)=\left(\frac{-3}{4}\rho^{2}+{3}\right)\vec{e_{z}}[/math]

Una vez calculados procedemos a representarlos:

Campo de presiones: como podemos observar en el gráfico la relación entre el incremento de la altura y el incremento de presión es lineal. Debido a esta condición a medida que aumenta o disminuye la altura lo hace de igual forma la presión.

centro

z=0:0.05:10;
f=(-3*z+3)/2;
plot(f,z)
xlabel('Incremento de presión');
ylabel('Incremento de altura');
title('Representación del campo de presiones del fluido');

Campo de velocidades: el campo de velocidades del fluido se trata de un campo vectorial.

centro

x=0:0.1:2;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
uy=(-3/4).*xx.^2+3;
ux=0.*yy;
hold on
quiver(xx,yy,ux,uy)
axis([0,3,0,10])
hold off
view(2)
title('Representación del campo de velocidades del fluido')

5 Líneas de Corriente del Campo

Para dibujar las líneas de corriente primero debemos calcular el campo [math]\vec{v}[/math] que es ortogonal a [math]\vec{u}[/math] en cada punto. [math]\vec{v}=\vec{e_{\theta}}\times\vec{u}=\begin{vmatrix}\vec{e_{\rho}} & \vec{e_{\theta}} & \vec{e_{z}} \\ 0 & 1 & 0\\ 0 & 0 & f\left (\rho\right)\end{vmatrix}=f\left (\rho\right)\overrightarrow{e_{\rho}}\rightarrow \overrightarrow{v}=f\left (\rho\right)\overrightarrow{e_{\rho}}[/math].

Sustituyendo [math]f\left(\rho\right)[/math] en la ecuación se obtiene:

[math]\overrightarrow{v}=\frac{1}{\mu}\left(\frac{p_{2}-p_{1}}{4}\rho^{2} + p_{1}-p_{2}\right)\overrightarrow{e_{\rho}}[/math].

Como [math]\vec{v} [/math] tiene un potencial escalar [math]\psi[/math], calculamos [math]\psi[/math] sabiendo que [math]\bigtriangledown\psi=\vec{v}[/math], la cual se conoce como función de corriente de [math]\vec{u}[/math].

[math]\frac{d\psi}{d\rho}=\frac{p_{2}-p_{1}}{4\mu}\rho ^{2} +\frac{p_{1}-p_{2}}{\mu}\ [/math].

[math]\psi=\frac{p_{2}-p_{1}}{12\mu}\rho^{3} + \frac{p_{1}-p_{2}}{\mu}\rho[/math].

Sustituyendo los valores de [math]p_{1}[/math], [math]p_{2}[/math] y [math]\mu[/math] obtenemos:

[math]\psi=-\frac{1}{4}\rho^{3} + {3}\rho[/math]

Una vez calculado [math]\psi[/math] dibujamos las líneas de corriente del campo [math]\psi=cte[/math]:

centro

rho=0:0.05:2;
z=0:0.05:10;
[R,Z]=meshgrid(rho,z);
F=-(1/4).*(R.^3)+3.*R;
contour(R,Z,F);
colorbar
axis([0,3,0,10]);
title('Representación de las líneas de corriente')

6 Velocidad Máxima del Fluido

En la gráfica podemos observar el comportamiento del módulo de la velocidad. Vemos también como el máximo de la velocidad se alcanza en [math]\rho={0}[/math], en el eje de la tubería tomando un valor de [math]{3}[/math] y que disminuye según nos acercamos a las paredes de la tubería donde la velocidad se hace nula.

centro

rho=0:0.05:2;
f=abs((-3/4)*rho.^2+3);
plot(rho,f);
title('Comportamiento del módulo de la velocidad')
xlabel('Radio')
ylabel('Velocidad')
axis([0,3,0,5])

7 Rotacional del Campo

Calculamos el rotacional de el campo de velocidades.

[math]\bigtriangledown\times\overrightarrow{u}\left ( \rho ,\theta ,z \right )= \frac{1}{\rho }\begin{vmatrix} \overrightarrow{e_{\rho }} &\rho \overrightarrow{e_{\theta }} & \overrightarrow{e_{z}}\\ \frac{d}{d\rho }& \frac{d}{d\theta } &\frac{d}{dz} \\ 0 & 0 &f\left ( \rho \right ) \end{vmatrix}=\frac{1}{2}\frac{p2-p1}{\mu }\rho \overrightarrow{e_{\theta }}[/math]

Si sustituimos los valores anteriores llegamos a que [math]\bigtriangledown\times\overrightarrow{u}\left ( \rho ,\theta ,z \right )=\frac{-3\rho}2[/math] Podemos observar que los puntos con mayor rotacional se encuentran en los extremos, al depender solo de [math]\rho[/math].

centro

x=0:0.1:2;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs((-3/2).*xx);
hold on
surf(xx,yy,rot);
colorbar;
axis([0,2.5,0,10]);
title('Rotacional de  u');
hold off

8 Temperatura del fluido

La temperatura de nuestro fluido viene dada por el campo escalar: [math]{T}(\rho,\theta,z)= log (1 + \rho) e^{-(z-2)^2}[/math]

El campo de temperatura se observara a mayor detalle en la siguiente gráfica, realizada en Matlab

izquierda

centro

x=0:0.01:2;
y=0:0.05:10;
[X,Y]=meshgrid(x,y);
figure (1)
p=log(1+X).*exp(-(Y-2).^2);
pcolor(X,Y,p);
shading flat
hold on
title('Campo de temperaturas')
axis([0,3,0,10]);
colorbar 
hold off
figure(2)
hold on
contour(X,Y,p,'k'); 
title('Curvas de nivel')
hold off
axis([0,2,0,10]);

9 Gradiente de la Temperatura

Para estudiar la variación de temperatura, analizamos el gradiente que indica la dirección de crecimiento de la función temperatura en cada uno de los puntos. Con este gráfico creado, podemos observar que las curvas de nivel son ortogonales al gradiente de la temperatura. Esto se debe a que el gradiente es el responsable de la dirección de crecimiento y decrecimiento del campo.

Como podemos ver el campo escalar viene dado en coordenadas cilíndricas por lo que que calculo el gradiente con la siguiente fórmula:

[math]\bigtriangledown{T}(\rho,\theta,z)=\frac{df}{dp}\cdot(e\rho)+ \frac{1}{p} \cdot\frac{df}{d\theta} \cdot(e\theta)+ \frac{df}{dz}\cdot(ez)[/math]

Finalmente, el gradiente queda expresado en el siguiente campo:

[math]\bigtriangledown{T}(\rho,\theta,z)=\frac{e^{-(z-2)^2}}{1+\rho} -2 \cdot log(\rho + 1) \cdot (z-2) \cdot{e^{-(z-2)}}[/math]

centro

x=0:0.05:2;
y=0:0.05:10;
[X,Y]=meshgrid(x,y);
figure (1)
p=log(1+X).*exp(-(Y-2).^2);
[pX,pY]=gradient(p);
hold on
quiver(X,Y,pX,pY)
contour(X,Y,p,'k');
axis([0,2,0,10]);
title('Gradiente temperatura y curvas de nivel');
shading flat
grid on
hold off

10 Caudal

Para poder calcular el caudal que llevara la tubería, se estudia el volumen de fluido que pasa a través de la tubería. La integral de la superficie es la siguiente:

[math]\int_{S}\vec{v} \cdot d\vec{S}[/math]

La [math]\vec{e_{v}}[/math] es el campo de velocidades, este campo es el calculado anteriormente.

[math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}} [/math], siendo [math] f\left(\rho\right)=\frac{1}{\mu}\left(\frac{\left(p_{2}-p_{1}\right)}{4}\rho^{2}+p_{1}-p_{2}\right)\vec{e_{z}}[/math]

Sustituimos los valores [math] p_{1} = 4, p_{2} = 1, {\mu} = 1[/math]

[math]Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left(\frac{1}{\mu}\right) \left(\frac{\left(p_{2}-p_{1}\right)}{4}\rho^{2}+p_{1}-p_{2}\right)\vec{e_{z}} dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{1}{4}\frac{p2-p1}{\mu }\rho ^{2}+\frac{p1-p2}{\mu } \right )d_{\rho }d_{\theta }= =\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{3x{^3}}{4}+ 3x\right)d_{\rho }d_{\theta }=...= {6\Pi }[/math]

[math]= {6\pi} =[/math] 18,85(m/s)