Flujo de Poiseuille (GRUPO 10A)

Trabajo realizado por estudiantes
Título	Flujo de Poiseuille (GRUPO 10A)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Alba Piedad Prats Moreno Carlos Muñoz González Carla De Juan Merchán Rodrigo Prado Fornos Miguel Vela Gonçalves Cerejeira
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

La ley de Poiseuille se utiliza para describir el flujo estacionario y laminar de un líquido incompresible; es decir, la densidad del líquido solo cambia si se le aplica una presión determinada.

En el estudio de esta ley nos enfocamos en el flujo de un líquido incompresible a través de una tubería cilíndrica cuyo radio es 2 centrada en el eje OZ.

La magnitud de este flujo viene determinada por el gradiente de presión y el radio de la propia tubería, teniendo en cuenta la viscosidad del líquido y la longitud de la tubería.

Para el desarrollo de este estudio hemos utilizado los programas Octave y Matlab y hemos trabajado en coordenadas cilíndricas.

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

Representación del mallado de dimensión 2, de la sección longitudinal del eje x=0. Consideraremos la región encerrada en las coordenadas (rho,z) = [0,3] x [0,10].

Mallado de la sección

%1. Definimos los ejes
rho=0:0.2:2; 
z=0:0.2:10; 
%2. Definimos mallado en 2 dimensiones
[xx,yy]=meshgrid(rho,z); %Mallado XY.
hold on
mesh(xx,yy,0.*xx); %Representamos la tubería
axis([0,3,0,10]); %Rango de los ejes
xlabel('ρ') ;
ylabel('z') ;
view(2);
title ('Mallado de la sección');
hold off

3 Ecuación de Navier-Stokes Estacionaria

La ecuación de Navier-Stokes describe como se mueve un fluido newtoniano. El fluido newtoniano se caracteriza porque su deformación es proporcional al esfuerzo cortante (fuerza tangencial).

Esta herramienta es esencial para comprender como se comportan los fluidos en sistemas hidráulicos, como tuberías, canales... Antes de sumergirnos en la demonstración es crucial establecer que el fluido es incompresible ya que esta ecuación rige el comportamiento de los fluidos newtonianos.

Vamos a definir nuestra fuerza según su presión y superficie, relacionándolo con la viscosidad del medio.

En el enunciado nos proporcionan la siguiente igualdad:

[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)/2 \right )[/math]

Las variables que compreenden la ecuación de Navier-Stokes son las siguiente:

- [math]\mu[/math] es el coeficiente de viscosidad.

- p1 representa la presión para valores de z=1.

- p2 lo hará para los valores de z=3.

El primer térmido de la ecuación será:

[math](\vec{u} \cdot \nabla) \vec{u})[/math]

Si desarrollamos este primer término, como el fluido es incompresible, no hay fuerzas externas (fuentes ni sumideros) por lo que la divergencia va a ser nula. Consequentemente, podemos despreciar el primer término de la ecuación Navier-Stokes obteniendo:

[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 ecuación representa la derivada segunda de f( [math] \rho [/math]). Para obtener el valor de [math] \rho [/math] podemos hacer una integral doble y a continuación multiplicar por [math] \rho [/math] ya que se trata de un gradiente en coordenadas cilíndricas:

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

Tenemos ahora que encontrar una solución particulas para las constantes de integración vamos a tener en cuenta el siguiente supuesto:

La velocidad del fluido converge a 0 en las paredes de la tubería. Consequentemente, mayor será la velocidad del fluido cuando nos acercamos al centro de la tibería. Como el radio de la tubería es 2, el valor de la velocidad en los extremos es nula.

Es decir, [math]f\left (2\right )=0[/math], [math]f\left ( 0 \right )=0[/math],

Vamos a proponer un sistemas de dos ecuaciones y obtener las constantes de integración:

[math]K=\frac{p1-p2}{\mu}[/math]

[math]C=0[/math]

Finalmente obtenemos la expresión f([math] \rho [/math]):

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

Para comprobar la condición de incomprensibilidad (el agua siempre ocupa el mismo volumen), vamos a comprobar que la divergencia del campo es nula.

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

En conclusión, como f(p) no depende de "z", la divergencia es nula.

4 Campo de presiones y campo de velocidades

Vamos a dar como dato: [math] p1=4 [/math], [math] p2=1 [/math] y [math] \mu=1 [/math].

4.1 Campo de velocidades

Anteriormente hemos obtenido la función de velocidad. Vamos a proceder a sustituir los valores dados:´

[math]\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4} +3)\vec{e_z}[/math]

Como podemos comprobar, a medida que nos acercamos al borde de la tubería, entra menos velocidad mientras que en el centro, hay mayor flujo de concentración de velocidades.

Campo de velocidades.

x=0:0.2:2; %Intervalo de 'ρ'
z=0:0.2:10; %Intervalo de 'z'
[X,Z]=meshgrid(x,z);
ux=(-3./4).*X.^2+3; %Dar valor a 'u'
uz=0.*zz;
hold on
quiver(X,Z,ux,uz)  %Representación
axis([0,5,0,12])
xlabel('ρ');
ylabel('z');
hold off
view(2)
title('Campo de velocidades')

4.2 Campo de presiones

La expresión del campo de presiones viene dada como:

[math] p(x,y,z)=p_1+(p_2-p_1)(z-1)[/math]

[math]p_1=4 \ {,} \ p_2=1 \ {,} \ μ=1[/math]

Obtenemos por tanto:

[math] p(x,y)=4-(3z+3)=1-3z[/math]

Como podemos comprobar, la presión NO aumenta con el radio, sino con el parámetro "z".

Campo de presiones

x=0:0.05:2; %Rango de 'ρ'
z=0:0.05:10; &Rango de 'z'
[X,Z]=meshgrid(x,z); 
figure (1)
p=1-3*Z; %Definimos rho
surf(X,Z,p)
view(2)
axis([0,5,0,12])
colorbar
xlabel('ρ')
ylabel('y')
Title('Campo de presiones')

5 Líneas de corriente

Para dibujar las líneas de corriente del campo [math]\overrightarrow{u}[/math], debemos tener en cuenta que estas son tangenes a [math]\overrightarrow{u}[/math] en cada apunto.

Procedemos a calcular [math]\overrightarrow{v}[/math] que es ortogonal a [math]\overrightarrow{u}[/math] ya que se comprueba que: [math]\overrightarrow{v}=\overrightarrow{e_{\theta }}\cdot \overrightarrow{u}[/math].

Como hemos hallado anteriormente, la divergencia de [math]\overrightarrow{u}[/math] es nula. Consequentemente, el rotacional del campo [math]\overrightarrow{v}[/math] es nulo también.

La función de corriente o potencial escalar viene definido como: [math]\psi [/math],([math]\bigtriangledown \psi =\overrightarrow{v}[/math])

[math]\overrightarrow{v}=\begin{vmatrix} \overrightarrow{e_{\rho }} & \overrightarrow{e_{\theta }}&\overrightarrow{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]:

[math]\overrightarrow{v}=\left ( \frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\right )\overrightarrow{e_{\rho }}[/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].

Obtenemos la relación: [math]\frac{d\psi }{d\rho }=\frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\ [/math].

Integrando el gradiente del potencial escalar: [math]\psi=\frac{p2-p1}{12\mu }\rho ^{3} +\frac{p1-p2}{\mu }\rho [/math].

El resultado de la función potencial quedará: [math] \psi= -\frac{1}{4} \rho ^{3} + 3\rho[/math].

Líneas de corriente

x=0:0.2:2; %Definimos 'ρ'
z=0:0.2:10; %Definimos 'z'
[X,Z]=meshgrid(x,z);
lineas=(-1./4).*X.^3+3.*X; %Definimos campo escalar  
contour(X,Z,lineas); 
axis([0,3,0,10]);
colorbar
title('Líneas de corriente');

6 Velocidad máxima del fluido y módulo de velocidad

Para obtener los máximos y mínimos de la velocidad vamos a optimizar la función de velocidad. Para ello, derivamos la misma y la igualamos a 0:

[math]\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4 } +3)\vec{e_z}[/math]

donde:

[math] \frac{\partial \vec{u}} {\partial \rho}=\frac{-3\rho}{2}\vec{e_z} [/math] [math] \frac{\partial \vec{u} }{\partial \theta}=0 [/math] [math] \frac{\partial \vec{u} }{\partial z}=0 [/math]

Consequentemente, el punto donde la velocidad es máxima es en [math] \rho =0 [/math]

Como podemos observar en la gráfica la velocidad de nuestro fluido no es constante a lo largo de la tubería. De hecho, disminuye a medida que nos acercamos a los bordes de la tubería. Esto se debe a las fuerzas viscosas potenciadas por la diferencia de presión que debe mantener el fluido para ser viscoso.

Velocidad de fluido

x=0:0.2:2;
z=0:0.2:10;
[X,Z]=meshgrid(x,z);
figure(1);
p=1-3.*Z;
surf(X,Z,p);
colorbar;
view(2);
axis([0,3,0,10]);
title('Velocidad de fluido');
xlabel('ρ');
ylabel('z');

7 Rotacional

Para hallar el rotacional, vamos a proceder con el cálculo:

[math]\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}[/math]

Sustituyendo en la fórmula:

[math]\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & (\frac{-3\rho ^{2}}{4 } +3)\end{vmatrix}[/math]

Obtenemos que [math]\nabla\times\vec{u}=(\frac{3\rho}{2})\vec{e_\theta}[/math],

Para visualizar esto de manera gráfica, nuevamente a Octave,

Rotacional del campo velocidad

x=0:0.1:2;
 z=0:0.1:10;
 [X,Z]=meshgrid(x,z);
 rot=abs(((9.*X/2+3./X);
 surf(X,Z,rot)
 colorbar
 view(2)
 axis([0,5,0,12])

8 Campo de temperaturas

En el enunciado nos viene dado la función Temperatura en coordenadas cilíndricas:

[math] T(ρ{,}θ{,}z)=log(1+ρ)e^{-(z-2)^2}[/math]

En la gráfica podemos observar como la temperatura es máxima para p=2, z=2; ya que los valores p y z son proporcionales a la temperatura.

Campo de temperaturas

p=0:0.01:2;
z=-2:0.05:10;
[P,Z]=meshgrid(p,z);
figure (1) 
a=log(1+P)*exp.^-(-Z-2).^2
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
10 colorbar 
11 figure(2)
12 contour(Y,Z,a,10,'k'); 
13 grid on
14 axis([0,8,-1,10]);

Curvas de nivel

Las curvas de nivel son líneas que concentran puntos con la misma temperatura. Se comprueba que tienen forma simétrica y logarítmica.

Curvas de nivel de campo de temperatura

y=0:0.05:2;
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

9 Gradiente de temperatura

Teniendo la función de temperatura, calculamos su gradiente:

[math]\nabla T(ρ,θ,z)=(\frac{e^{-(z-2)^2}}{1+ρ)}-2log(ρ+1)(z-2){e^{-(z-2)^2}})[/math]

La gráfica del gradiente de la temperatura, se ha representado utilizando el siguiente código:

rho=0:0.1:2;
z=0:0.1:10;
[RHO,Z]=meshgrid(rho,z);
figure(1)
T=log(1+RHO).*exp(-(Z-2).^2);
[TRHO,TZ]=gradient(T);
hold on
quiver(RHO,Z,TRHO,TZ)
xlabel('rho');
ylabel('z');
axis([0,2,0,10]);
title('Gradiente de Temperatura y Curvas de nivel')
shading flat
grid on
hold off

Gradiente ortogonal a las curvas de nivel

Código Octave utilizado:

Gradiente ortogonal a las curvas de nivel

Como podemos verificar, 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.

Por otro lado, podemos comprobar que si el gradiente es mayor, la temperatura tendrá una variación más rápida.

10 Caudal

El caudal es la cantidad de un fluido que atraviesa la tubería en un determinado tiempo. Para ello debemos tener en cuenta el flujo de volumen por unidad de tiempo. Vamos a representarlo de la siguiente forma:

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

donde [math]\vec{v}[/math] es el campo de velocidades.

El vector normal será perpendicular a la superficie por lo que hemos visto anteriormente de las curvas y gradientes.

El campo de velocidades finalmente será:

[math]\vec{u}(\rho,\theta,z)=(-\frac{3}{4}\rho ^{2} +3)\overrightarrow{e_{z}}[/math],

Finalmente resolviendo la integral doble obtendremos:

[math]Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left (\frac{-3\rho^{2}}{4}+{3\rho} \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}} dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{-3\rho^{2}}{4}+{3} \right )d_{\rho }d_{\theta }=8{\Pi}=25.133\left ( m/s \right )[/math].