Flujo de Poiseuille Grupo 18-B

Trabajo realizado por estudiantes
Título	Flujo de Poiseuille Grupo 18-B
Asignatura	Teoría de Campos
Curso	2022-23
Autores	Alejandro Ramos García Alberto Monge Barrera Alejandro López Pérez Louciana M. Contreras Rivera
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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, de la cual conocemos que la velocidad [math]\vec{u}(\rho,\theta,z)[/math] de las partículas viene dada por:

[math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}},[/math]

y su 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).[/math]

Posteriormente, representaremos el campo de velocidades y presiones gráficamente. Además, analizaremos los vectores ortogonales a [math]\vec{u}[/math], los puntos donde este es máximo, su rotacional y calcularemos el caudal que pasa por una de sus secciones transversales.

También, analizaremos una segunda cantidad física, la temperatura; la cual viene dada por el campo:

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

de la cual, dibujaremos el campo de temperaturas y las curvas de nivel del campo. Luego, calcularemos su gradiente y analizaremos sus vectores gráficamente.

1 Mallado. Sección transversal de la tubería.

Dibujaremos un mallado en dimensión 2, de la sección longitudinal de la tubería [math] x_{1} = 0 [/math]. Fijaremos los ejes en la región [math] \left ( \rho,z \right )\epsilon \left [ 0,3 \right ]\times \left [ 0,10 \right ]. [/math] Cabe comentar que hemos usado Matlab para este trabajo a la hora de obtener los mallados y gráficas pedidas.

Sección transversal

clc;
clear all;
x=0:0.3:3;
y=0:0.3:10;
[xx,yy]=meshgrid(x,y)
figure(1)
mesh(xx,yy,0*xx)
axis([-2,5,0,10])
view(2);

2 Velocidad y presión

2.1 Ecuación de Navier-Stokes. Divergencia nula.

2.1.1 Ecuación de Navier-Stokes.

Como ya hemos mencionado, la velocidad de las partículas del fluido viene dada por [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}[/math] y su presión viene dada por [math]p\left(x,y\right)=p_{1}+\left (p_{2}-p_{1}\right)(z-1) [/math], donde [math] p_{1} [/math] es la presión en los puntos [math] z=1 [/math], [math] p_{2} [/math] la presión en los puntos [math] z=2 [/math] y [math] \mu [/math] el coeficiente de viscosidad del fluido. Sabiendo que [math] \left ( \vec{u},\rho \right ) [/math] satisface la ecuación de Navier-Stokes estacionaria (independiente del tiempo):

[math] \left ( \vec{u}\cdot \triangledown \right )\vec{u}+\triangledown p=\mu\Delta\vec{u}, [/math]

en la que vamos a despreciar el primer término (parte convectiva), comprobar que [math]f\left ( \rho \right ) [/math] satisface la ecuación 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]

La resolveremos multiplicando por [math] \rho [/math] e integrando dos veces. Supondremos que la velocidad del fluido en [math] \rho=2 [/math] es nula y que no se hace infinito en [math] \rho=0 [/math].

Procedemos a operar:

1) Multiplicamos por [math] \rho [/math]

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

2) Integramos

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

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

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

3)Integramos por segunda vez

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

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

2.1.2 Divergencia nula.

En este apartado, comprobaremos si se verifica la condición de incompresibilidad, es decir, que el agua ocupa el mismo volumen:

[math] \triangledown \cdot \vec{u}=0 [/math]

La divergencia de una campo vectorial en un punto, más concretamente, de un campo de velocidades de partículas; se puede interpretar como el cambio en la densidad en ese punto.

Al ser el agua un fluido incompresible, la divergencia debería ser nula en cualquier punto, pues su densidad permanece constante como ya se ha comentado. Matemáticamente:

[math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}[/math]

[math]\triangledown\cdot\vec{u} \left(\rho,\theta,z\right)=\frac{1}{\rho} \left ( \frac{\partial }{\partial \rho}\left ( \rho u_{\rho } \right )+\frac{\partial }{\partial \theta}\left ( u_{\theta } \right ) + \frac{\partial }{\partial z}\left (\rho u_{z}\right)\right) =\frac{1}{\rho}\left ( \frac{\partial }{\partial z}\left(\rho f\left(\rho\right)\right)\right)=0 [/math] (para cualquier punto).

2.2 Campo de presiones y campo de velocidades.

Suponiendo que [math] p_{1}=2, p_{2}=1 [/math] y [math]\mu=1, [/math] calcularemos el campo de presiones y el campo de velocidades, es decir, un campo escalar y uno vectorial, respectivamente.

El campo de presiones nos queda:

[math]p\left ( x,y \right )=p_{1}+\left ( p_{2}-p_{1} \right )\left ( z-1 \right )=2+\left ( 1-2 \right )\left ( z-1 \right )=3-z.[/math]

El campo de velocidades nos queda:

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

2.2.1 Representación del campo de presiones.

Como ya se comentó, procedemos a representar la variación de la presión frente a la altura. Como puede verse, se trata de una relación lineal: la presión aumenta con la profundidad. O lo que es lo mismo, que está disminuye en estas condiciones solo si también disminuye la cantidad de fluido "por encima" del punto en cuestión.

clc;
clear all;
%La función de presiones gráficamente es una recta.
z=0:0.1:10;
f=3-z;
plot(z,f)
xlabel('Variación de altura');
ylabel('Variación de presión');
title('Campo de presiones del fluido');

2.2.2 Representación del campo de velocidades.

Nuevamente procedemos a representar un campo, esta vez de carácter vectorial: el campo de velocidades de las partículas del fluido. Eso sí, un breve inciso, al venir dado en coordenadas cilíndricas (como es normal al tratarse de una tubería cilíndrica),tendremos que pasarlo a cartesianas para poder representarlo en M.

clear;
clc;
x=0:0.1:3;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=-(1/4)*xx.^2;           %debido al cambio a coordenadas cartesianas
uy=-(1/4)*yy.^2;           %ídem
hold on
quiver(xx,yy,ux,uy)
axis([0,5,0,12])
hold off
view(2)
title('CAMPO DE VELOCIDADES DEL FLUIDO')

2.3 Líneas de corriente del campo.

Vamos a dibujar ahora las líneas de corriente del campo [math] \vec{u} [/math], es decir, las líneas que son tangentes a [math] \vec{u} [/math] en cada punto.

Para ello, calcularemos el campo [math] \vec{v} [/math], que en cada punto es ortogonal a [math]\vec{u}[/math]. Observaremos que [math] \vec{v} [/math] es irrotacional (por ser [math]\vec{u}[/math] de divergencia nula) y que tiene potencial escalar [math]\psi[/math], que se conoce como función corriente de [math]\vec{u}[/math].

Posteriormente, calcularemos [math]\psi [/math] y dibujaremos las líneas [math]\psi=cte[/math]. Comprobaremos que son efectivamente corrientes de [math]\vec{u}[/math], observando que son tangentes al campo 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\\ 1 & 0 & f\left ( \rho \right ) \end{vmatrix} = f\left ( \rho \right )\vec{e_{\rho }}=\frac{p_{2}-p_{1}}{4\mu}\rho ^{2}\vec{e_{\rho}} [/math]

El vector [math] \vec{v} [/math] es irrotacional al tener [math] \vec{u} [/math] divergencia nula. Matemáticamente se puede comprobar:

[math] \triangledown \times \vec{v}=\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} = \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}\\ f\left(\rho\right) & 0 & 0 \end{vmatrix} =\frac{1}{\rho}\left ( \frac{\partial f\left ( \rho \right ) }{\partial z} \rho \vec{e_{\theta}}- \frac{\partial f\left ( \rho \right ) }{\partial \theta}\vec{e_{z}}\right )=0 [/math]

Ahora calculamos [math]\psi[/math], función corriente de [math] \vec{u}. \left (\bigtriangledown\psi =\vec{v}\right) [/math]

[math] \triangledown \psi=\frac{\partial \psi}{\partial \rho} \vec{e_{\rho}} + \frac{\partial \psi}{\partial \theta} \frac{1}{\rho}\vec{e_{\theta}}+\frac{\partial \psi}{\partial z} \vec{e_{z}}=\frac{p_{2}-p_{1}}{4\mu}\rho ^{2}\vec{e_{\rho}} [/math]

[math] \psi=\frac{p_{2}-p_{1}}{4\mu}\int\rho ^{2}d\rho=\frac{p_{2}-p_{1}}{12\mu}\rho ^{3}=-\frac{1}{12}\rho^{3} [/math]

2.4 Puntos con velocidad máxima.

Representamos la variación de la velocidad de las partículas en función del radio, y la gráfica que nos sale tiene una naturaleza exponencial, es decir, la velocidad crece exponencialmente entre puntos situados a diferente distancia respecto al eje de la tubería cilíndrica.

%Módulo del campo de velocidades
clc
clear all
rho=0:0.1:3;
f=(1/4)*(rho.^2);
plot(rho,f);
%El módulo es máximo cuando el rho es máximo, es decir, en rho=3
xlabel('aumento del radio')
ylabel('aumento de la velocidad')
title('Módulo del campo de velocidades');

2.5 Rotacional.

Podemos considerar que la interpretación del rotacional es el giro alrededor de un punto del fluido. Si el fluido fuera perfecto, el rotacional sería nulo, ya que todas las velocidades serían iguales. En este caso, estamos trabajando con un fluido no perfecto (o no newtoniano), pues tiene coeficiente de viscosidad [math] \mu [/math], entonces el rotacional será distinto de cero.

Procedemos a calcularlo:

[math]\triangledown \times \vec{v}=\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} = \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 & f\left(\rho\right)\end{vmatrix} = -{f}'\left ( \rho \right )\vec{e_{\theta}} [/math]

Sustituímos con los valores de [math] p_{1}, [/math] [math] p_{2} [/math] y [math] \mu.[/math]

[math] \frac{p_{2}-p_{1}}{2\mu} = -\frac{\rho}{2} [/math]

Con esto podemos concluir que los puntos con mayor rotacional serán los puntos en los que la función tenga mayor pendiente.

2.5.1 Representación del campo del rotacional

clear;
clc;
x=0:0.1:3;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs(-(1/4)*xx.^2-(1/4)*yy.^2);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])

2.6 Caudal por una sección transversal

Podemos calcular el caudal de una sección transversal integrando el campo mediante una integral (doble) de superficie.

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

3 Temperatura

3.1 Campo de temperaturas y curvas de nivel.

Ahora procederemos de nuevo a representar y analizar un campo escalar: el campo de temperaturas de un fluido. Eso sí, nuevamente se insiste en la necesidad de pasar a cartesianas la función temperatura al venir esta dada en cilíndricas.

3.1.1 Representación del campo de temperaturas.

y=0:0.05:10;
z=-2:0.05:11;
[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 
title('CAMPO DE TEMPERATURAS')

3.1.2 Representación de las curvas de nivel del campo de temperaturas.

Las curvas de nivel de un campo de temperatura unen los puntos de igual temperatura, es decir, estamos ante las famosas líneas isotermas.

y=0:0.05:8;
z=-2:0.05:10;
[Y,Z]=meshgrid(y,z);
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 
contour(Y,Z,p,10,'k');
title('CURVAS DE NIVEL')

3.2 Gradiente.

El gradiente de la función de temperatura será:

[math] \bigtriangledown T\left ( \rho, \theta, z \right )=\frac{1}{4}e^{-\left ( z-1 \right )^{2}}\vec{e_\rho}+\left ( \rho-\frac{1}{2} \right )^{2}\cdot e^{-\left ( z-1 \right )^{2}}\cdot -2\left ( z-1 \right )\vec{e_z} [/math]

3.2.1 Representación del gradiente.

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
title('GRADIENTE')
hold off

3.2.2 Representación curvas de nivel del gradiente.

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
title('CURVAS DE NIVEL')
hold off