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:
y su presión [math]p\left(x,y\right)[/math] que viene dada por:
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:
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.
Cabe comentar que hemos utilizado Matlab para este trabajo a la hora de obtener los mallados y gráficas pedidas.
Contenido
- 1 Mallado. Sección transversal de la tubería.
- 2 Velocidad y presión
- 3 Temperatura
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]
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
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):
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
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} + c [/math]
Sabemos que en [math] \rho=2 [/math] la velocidad es nula. Como la velocidad es [math] f\left ( \rho \right )\vec{e_{z}}, [/math] entonces [math] f\left ( 2 \right ) [/math] debe ser nulo. Matemáticamente:
Sustituimos con los datos [math] p_{1}=2, p_{2}=1 [/math] y [math]\mu=1: [/math]
Por lo tanto,
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:
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}}=\left ( \frac{ p_{2}-p_{1}}{4\mu}\rho^{2} + 1 \right )\vec{e_{z}} =\left ( -\frac{1}{4}\rho ^{2}+1\right )\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 Matlab.
clear all;
clc;
x=0:0.1:3;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=(-(1/4)*xx.^2)+1; %debido al cambio de coordenadas cartesianas
uy=(-(1/4)*yy.^2)+1; %í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 }}=\left ( \frac{p_{2}-p_{1}}{4\mu}\rho ^{2}+1\right )\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}}=\left ( \frac{p_{2}-p_{1}}{4\mu}\rho ^{2}+1\right )\vec{e_{\rho}} [/math]
[math] \psi=\frac{p_{2}-p_{1}}{4\mu}\int\rho ^{2}d\rho + \int d\rho =\frac{p_{2}-p_{1}}{12\mu}\rho ^{3} + \rho [/math]
Sustituimos [math] p_{1}=2, p_{2}=1 [/math] y [math]\mu=1 [/math] y nos queda:
[math] \psi=-\frac{1}{12}\rho^{3} + \rho [/math]
2.3.1 Representación de las líneas de corriente del campo.
clear;
clc;
%Líneas de corriente
rho=0:0.1:3,
z=0:0.1:10;
[R,Z]=meshgrid(rho,z);
F=(-1*(1/12).*(R.^3))+R;
contour(R,Z,F,10);
Ahora comprobamos que las líneas de nivel son tangentes al campo de velocidades.
2.4 Puntos con módulo de 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))+1;
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');
axis([0,3,0,3.5])
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{u}\left ( \rho, \theta,z\right )=\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} = -\frac{\partial f\left ( \rho \right )}{\partial \rho } \vec{e_{\theta}}[/math]
Sustituímos con los valores [math] p_{1}=2, p_{2}=1 [/math] y [math]\mu=1 [/math]
[math] \triangledown \times \vec{u}\left ( \rho, \theta,z\right )=\frac{\left ( p_{2}-p_{1} \right )\rho}{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.
Para ello, primero parametrizaremos la superficie:
Ahora, calculamos el vector normal a la superficie [math] \left.\begin{matrix} \vec{r_u}=\vec{e_\rho}\\ \vec{r_v}=u\vec{e_\theta} \end{matrix}\right\} \vec{r_u}\times\vec{r_v}=u\vec{e_z} [/math], y procedemos a calcular el caudal:
[math]Q=\iint_{S}^{}\vec{u}\left ( \vec{r}\left ( u, v \right )\cdot \left ( \vec{r_u}\times\vec{r_v} \right ) \right )dudv=\iint_{S}^{}\left ( -\frac{1}{4}u^{2}+1 \right )\vec{e_z}\cdot u \vec{e_z}dudv= \int_{0}^{2\pi}\int_{0}^{3}\left ( -\frac{1}{4}u^{3}+u \right )dudv= -\frac{9}{8}\pi=-3,5343\left [ m/s \right ][/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, el cual viene dado por:
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.
Representamos el campo de temperaturas con la escala cromática predeterminada y podemos observar como en un punto en las inmediaciones del punto (8,1) nuestro fluido alcanza la máxima temperatura.
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 de la temperatura.
El gradiente de la función de temperatura será:
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 Comparación entre el gradiente y las curvas de nivel de temperatura.
Por útimo, procedemos a representar las curvas de nivel de temperatura junto al gradiente, para posteriormente (y haciendo zoom en la gráfica generada por Matlab) comprobar que, efectivamente, en cada punto el gradiente es ortogonal a la isoterma que pasa por este.
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
