Flujo de Poiseuille (Grupo 27-A)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Flujo de Poiseuille (Grupo 27-A) |
| Asignatura | Teoría de Campos |
| Curso | 2022-23 |
| Autores | Pablo Díaz Paz-Albo Gonzalo de la Flor Fernández Jeremy García Herrera Juan Carlos Santos Expósito Emilio Villegas Maroto |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 INTRODUCCIÓN
La ley de Poiseuille o ley de Hagen-Poiseuille es una ley que permite determinar el flujo laminar estacionario de un líquido incompresible y uniformemente viscoso (fluido newtoniano) a través de un tubo cilíndrico de sección circular constante.
Para la realización de este trabajo hemos utilizado el programa informático Matlab, para poder representar los resultados de una manera visual y ayudar al usuario a comprender mejor las interpretaciones. A continuación se muestra a los dos autores de la ley, Jean Léonard Marie Poiseuille y Gotthilf Heinrich Ludwig Hagen.
2 SUPERFICIE DE TRABAJO
Vamos a considerar el flujo de un fluido incompresible en una tubería cilíndrica de radio 2 que supondremos centrada en el eje OZ. Los ejes [math]\rho [/math], y z se han definido en los intervalos [math][0,3]×[0,10][/math] respectivamente.
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
mesh(xx,yy,0*yy);
axis([0,5,0,12]);
view(2);
3 ECUACIÓN DE NAVIER-STOKES ESTACIONARIA
La ecuación de Navier-Stokes nos permite modelizar el movimiento de un fluido newtoniano. Esto nos resultará muy útil para conocer el comportamiento de los fluidos en instalaciones hidráulicas como tuberias, canales, etc... Lo primero que tenemos que demostrar es que el campo de presiones y velocidades al que esta sometido el fluido cumple con la ecuación de Navier-Stokes. Para ello se debe partir de que el fluido es incompresible ya que esta ecuación determina el comportamiento de fluidos Newtonianos.
[math]\overrightarrow{u}\left ( \rho ,\theta ,z \right )=f\left ( \rho \right )\overrightarrow{ez}[/math] y su presión [math]p\left ( x,y \right )=p1+\left ( p2-p1 \right )\left ( z-1 \right )[/math], donde [math]\mu[/math] es el coeficiente de viscosidad, p1 es la presión en los puntos z = 1, p2 la presión en los puntos z = 2. Para entender la formula es importante conocer a que se refiere cada termino, pues [math]\mu[/math] es el coeficiente de viscosidad, p1 representa la presión para valores de z=1, y por consiguiente p2 lo hará para los valores de z=2. El primer térmido de la ecuación será: [math](\vec{u} \cdot \nabla) \vec{u})[/math] Si sustituimos por nuestro campo y lo desarrollamos veremos que nos queda el vector nulo por lo que el primer término de la ecuación de Navier-Stokes no interviene.
[math] \frac{1}{\rho } \frac{\partial }{\partial \rho } \left ( \rho \frac{\partial f(\rho )}{\partial \rho } \right ) = \frac{p_{2}-p_{1}}{\mu }[/math]
Esta última ecuación es una derivada segunda del campo requerido, por lo que debemos integrar dos veces, una vez múltiplicado por [math] \rho [/math]
[math]f\left ( \rho \right )=\frac{p2-p1}{4\mu}\rho ^{2} +Cln\left (\rho \right )+K[/math] Para continuar debemos comprender el sistema, y saber que en [math] \rho=2 [/math] la ecuación será 0 y para [math] \rho=2 [/math] no tenderá a infinito, lo que nos deja lo siguiente:
[math]f\left (2\right )=0[/math], [math]K=\frac{p1-p2}{\mu}[/math]
[math]f\left ( 0 \right )=0[/math], [math]C=0[/math]
[math]f\left ( \rho \right )=\frac{1}{\mu }\left [ \frac{p2-p1}{4 }\rho ^{2} +p1-p2\right ][/math]
Como ya se comentó al principio del apartado, el fluido es incompresible, para comprobarlo calcularemos la divergencia del campo, y al darnos 0 podemos concluir que la hipótesis inicial es valida, pues en caso de ser el resultado negativo significaría que se contrae y lo contrario en caso de resultado positivo.
[math]\bigtriangledown \cdot \overrightarrow{u}=0=\frac{1}{\rho } \left\{ \frac{\partial }{\partial \rho }\left ( \rho u_{\rho } \right )+ \frac{\partial }{\partial \theta }(u_{\theta })+ \frac{\partial }{\partial z}(\rho u_{z})\right\}[/math]
4 CAMPO DE VELOCIDADES
Para la representación del campo se han tomado los siguientes valores:
Sustituyendo en la expresión del campo obtenemos:
Para su representación se ha introducido en Matlab, el siguiente código:
[Y,Z] = meshgrid([0:0.1:1],[0:0.25:10]);
uy=inline('y.*0','y','z');
uz=inline('((-y.^2)./4)+1','y','z');
U=uy(Y,Z);
V=uz(Y,Z);
axis([0,5,0,12]);
quiver(Y,Z,U,V);
Según la Ley de Poiseuille la velocidad en las paredes o lateral del tubo debería ser nula, tal y como se muestra en la segunda imagen. En nuestro caso, la velocidad en el lateral [math] \rho =2 [/math] es nula, y la velocidad máxima se encuentra en el extremo opuesto, es decir, [math] \rho =0 [/math], como podemos observar en el gráfico realizado con Matlab.
4.1 CÁLCULO DE LA VELOCIDAD MÁXIMA DEL FLUIDO
Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:
donde:
Si igualamos a 0, [math] \frac{-\rho}{2}=0[/math] obtenemos que el punto donde se produce la máxima velocidad es [math] \rho =0 [/math], como hemos podido observar en el gráfico anterior del campo de velocidades.
5 CAMPO DE PRESIONES
El campo de presiones al que está sometido nuestro fluido es un campo escalar, con las siguiente expresión:
Para representarlo gráficamente se consideran nuevamente los valores: [math]p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1[/math] .Obteniendo, por tanto:
Implementado el código siguiente en Octave:
y=0:0.05:2;
z=0:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=3-Z;
surf(Y,Z,p)
view(2)
axis([0,5,0,12])
Como se puede apreciar en el gráfico a medida que el fluido avanza por el canal la presión va disminuyendo. Las presiones más altas están representadas con colores cálidos, y las presiones más bajas con colores fríos. Esto se debe a que para mantener el caudal de un fluido viscoso estable debe mantenerse una diferencia de presiones entre las paredes del canal, esta diferencia de presión es necesaria debida a la fuerza de arrastre o frenada que ejerce el canal sobre la capa de fluido en contacto con él. Estas fuerzas de arrastre o de frenado se denominan fuerzas viscosas. El resultado, hace que la velocidad del fluido no sea constante a lo largo del canal siendo mayor cerca del centro y menor cerca de las paredes.
5.1 LÍNEAS DE CORRIENTE
Vamos a dibujar las líneas de corriente del campo [math]\overrightarrow{u}[/math], es decir, las líneas que son tangenes a [math]\overrightarrow{u}[/math] en cada apunto, es necesario calcular el campo [math]\overrightarrow{v}[/math] que en cada punto es ortogonal a [math]\overrightarrow{u}[/math], Siendo [math]\overrightarrow{v}=\overrightarrow{e_{\theta }}\cdot \overrightarrow{u}[/math]. El campo [math]\overrightarrow{v}[/math] es irrotacional por ser nula la divergencia de [math]\overrightarrow{u}[/math] y ,además, tiene un potencial escalar [math]\psi [/math],([math]\bigtriangledown \psi =\overrightarrow{v}[/math]), que se conoce como función de corriende de [math]\overrightarrow{u}[/math].
Sustituyendo [math]f\left ( \rho \right )[/math] obtenemos:
Para calcular [math]\psi [/math], obtenemos su gradiente e igualmaos al campo [math]\overrightarrow{v}[/math], [math]\bigtriangledown \psi =\overrightarrow{v}[/math].
[math]\bigtriangledown \psi=\frac{d\psi }{d\rho }\overrightarrow{e_{\rho }}+\frac{1}{\rho }\frac{d\psi }{d\theta }\overrightarrow{e_{\theta }}+\frac{d\psi }{dz}\overrightarrow{e_{z}}=\overrightarrow{v}[/math]. Por ende: [math]\frac{d\psi }{d\rho }=\frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\ [/math].
Tras integrar nos queda finalmente que [math]\psi=\frac{p2-p1}{24\mu }\rho ^{2} +\frac{p1-p2}{\mu }\rho [/math].
Sustituimos en la anterior función Potencial y como resultado final:Sustituyendo en la anterior función los valores asignados, el resultado de la función potencial será: [math] \psi= -\frac{1}{24} \rho ^{2} + \rho[/math]. Con ayuda del programa Matlab se muestran las líneas de corriente de la siguiente manera:
x=0:0.2:2;
y=0:0.2:10;
[X,Y]=meshgrid(x,y);
lineas=(-1./24).*X.^2+X;
contour(X,Y,lineas);
axis([0,3,0,10]);
colorbar
title('Líneas de corriente');
6 ROTACIONAL
Para la realización del cálculo del rotacional del campo, utilizaremos la siguiente fórmula:
Sustituyendo en la fórmula:
Obtenemos que [math]\nabla\times\vec{u}=(\frac{-3\rho}{4}+\frac{1}{\rho})\frac{(p_1-p_2)}{μ}\vec{e_\theta}[/math], y dando los siguientes valores:
Tenemos que [math]\nabla\times\vec u=(\frac{-3\rho}{2}+\frac{2}{\rho})\vec{e_\theta}[/math], siendo el módulo del rotacional,[math]\left | \nabla\times\vec u \right |= (\frac{-3\rho}{2}+\frac{2}{\rho}) [/math], como se observa en la expresión, depende del valor de [math]\rho[/math]. Tomando los valores mínimos en los extremos del intervalo en el que está definido [math]\rho[/math], [math] [0,2][/math].
Para visualizar esto de manera gráfica, nuevamente a Octave,
x=0:0.1:2;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs(((-y.^2)./4)+1);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])
El rotacional muestra la tendencia del campo a inducir rotación alrededor de un punto. Por otra parte, el rotacional caracteriza la rotación de un fluido, por lo que en los extremos del canal el efecto de giro será máximo, particularizando en [math]\rho=2[/math] tendrá el valor máximo en el sentido positivo de [math]\vec{e_z}[/math]. En el gráfico los puntos con mayor tendencia a la rotación aparecen representados con colores cálidos, amarillos hacia verde azulado, y los puntos con menor tendencia a la rotación aparecen representados con colores fríos, azules. Como se puede apreciar en el gráfico aparece la tendencia del rotacional en torno a las paredes del canal [math]\rho=2[/math].
7 CAMPO DE TEMPERATURAS
La temperatura a la que está sometido el fluido viene determinada por el campo escalar:
Como se observa, la expresión del campo está en coordenadas porlares. Teniendo en cuenta que las relación para pasar de coordenadas polares a cartesianas es:
- [math]x=ρcosθ[/math]
- [math]y=ρsenθ[/math]
Pero como nosotros estamos en los ejes [math]\{\vec{j} {,} \vec{k} \}[/math], usaremos:
- [math]y=ρcosθ[/math]
- [math]z=ρsenθ[/math]
Podemos expresar el campo en coordenadas cartesianas, que es el sistema en el cual estamos trabajando, resultando:
REPRESENTACIÓN DEL CAMPO DE TEMPERATURAS
El campo de temperaturas es un campo escalar que representa la distribución espacial de la temperatura asociando un valor a cada punto del espacio. Depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión [math] T(y{,}z)=1+(y^2+z^2-1/2) e^{-((z-1)^2)}[/math], por lo tanto no depende del tiempo, se dice que es estacionario, sino exclusivamente de las componentes espaciales [math](y{,}z)[/math]
Para la representación se ha implementado el siguiente código en Octave:
y=0:0.05:8;
z=-2:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
colorbar
figure(2)
contour(Y,Z,p,10,'k');
grid on
axis([0,8,-1,10]);
En el gráfico se puede observar que la temperatura es mayor cuándo se representa con colores cálidos, esto sucede en los valores próximos a [math]z=1[/math].
REPRESENTACIÓN DE LAS CURVAS DE NIVEL
El código utilizado para la representación de este gráfico se incluyó con el de la representación del campo de temperaturas.
En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica estando más próximas entre sí cuando la temperatura aumenta, en torno a los puntos [math]z=1[/math] y más alejadas cuando la temperatura disminuya.
La representaciónes anteriores permite determinar que en el punto [math](y,z)= (8,1)[/math] la temperatura es máxima.
7.1 GRADIENTE DEL CAMPO DE TEMPERATURAS
El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente [math]\nabla T[/math], en nuestro caso:
Como observación, el gradiente analíticamente se calcula como [math]\frac{\partial T(y,z)}{\partial y}[/math] siendo ésta la compontente [math]\vec{j}[/math] del campo vectorial gradiente y [math]\frac{\partial T(y,z)}{\partial z }[/math] la componente [math]\vec{k}[/math].
REPRESENTACIÓN GRÁFICA DEL GRADIENTE
El gráfico del gradiente de la temperatura, se ha representado utilizando el siguiente código:
y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+((Z.^2).*exp(-((Y.^2)+(Z.^2)-1/2).^2));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
axis([0,8,-1,2]);
shading flat
grid on
hold off
Código Octave utilizado:
y=0:0.05:8;
z=0:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
contour(Y,Z,p,'k');
axis([0,8,-1,10]);
shading flat
grid on
hold off
Con este gráfico afirmamos que las curvas de nivel son ortogonales al gradiente de la temperatura.
Para estudiar la variación de temperatura a lo largo del canal analizamos el gradiente que indica la direccion de crecimiento de la función temperatura en cada uno de los puntos. Como ya se ha mencionado antes, junto al gradiente aparecen las curvas de nivel de la temperatura. Se puede comprobar que en torno al punto [math]z=1[/math] el módulo del gradiente aumenta, lo que indica que la temperatura tendrá una variación más rápida en esta zona. Se aprecia claramente como el gradiente es ortogonal a las curvas de nivel, esto se debe a que el gradiente indica la dirección de crecimiento y sentido de los valores del campo en cada punto y las curvas de nivel son el lugar geométrico de los puntos equipotenciales.
8 CAUDAL
El caudal es el volumen de un fluido que pasa a través de la tubería estudiada por unidad de tiempo y se calcula mediante la integral de superficie siguiente:
donde [math]\vec{v}[/math] es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.
En nuestro caso, el campo de velocidades es:
Esa expresión resulta de sustituir los valores siguientes en la expresión incial del campo:
La profundidad del caudal es de 1 metro, por lo que tras parametrizar realizamos el calculo del caudal.
CÁLCULO DEL CAUDAL
[math]Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left ( \frac{1}{4}\frac{\rho ^{2}}{\mu }+\frac{\rho}{\mu } \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}} dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{1}{4}\frac{\rho ^{2}}{\mu }+\frac{\rho}{\mu } \right )d_{\rho }d_{\theta }=11.52\left ( m/s \right )[/math].
