Diferencia entre revisiones de «Flujo de Poiseuille (Grupo 23)»

Revisión del 12:21 14 dic 2023

1 Introducción

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 [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 y como viene dado en coordenadas polares antes de representarlo debemos cambiarlo a coordenadas cartesianas.

centro

x=0:0.05:2;
y=0:0.05:10;
[xx,yy]=meshgrid(x,y);
ux=-3/4*xx.^2+3;
uy=-3/4*yy.^2+3;
hold on
quiver(xx,yy,ux,uy)
axis([0,3,0,11])
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

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,10])

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

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)}[/math]

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)}}{1+\rho} -2 \cdot log(\rho + 1) \cdot (z-2) \cdot{e^{-(z-2)}}[/math]

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

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]);

curvas de nivel

x=0:0.05:8;
y=0:0.05:1;
[X,Y]=meshgrid(x,y);
figure (1)
p=1+((Y.^2).*exp(-((X.^2)+(Y.^2)-1/2).^2));
[pX,pY]=gradient(p);
hold on
quiver(X,Y,pX,pY)
axis([0,8,-1,2]);
shading flat
grid on
hold off

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.

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

10 Caudal