Diferencia entre revisiones de «Flujo de Poiseuille (GRUPO 10A)»

Revisión del 19:21 13 dic 2023

1 Introducción

La ley de Poiseuille se utiliza para describir el flujo estacionario y laminar de un líquido incomprensible.

En el estudio de esta ley nos enfocamos en el flujo de un líquido incomprensible 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, 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. 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. 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 el valor es nulo. Por ello, podemos despreciar el primer término de la ecuación Navier-Stokes obteniendo:

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:

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]C=0[/math]

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

Para comprobar la condición de incomprensibilidad (el agua siempre ocupa el mismo volumen), vamos a hallar la divergencia del campo:

En conclusión, como la divergencia es nula, podemos afirmar que no hay variaciones de volumen.

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:´

Para su representación se ha introducido en Matlab, el siguiente código:

Campo de velocidades.

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

[math]p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1[/math] .Obteniendo, por tanto:

Implementado el código siguiente en Octave:

Campo de presiones

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

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

Sustituyendo [math]f\left ( \rho \right )[/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].

Con ayuda del programa Matlab se muestran las líneas de corriente de la siguiente manera:

Líneas de corriente

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 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:

donde:

[math] \frac{-\rho}{2}=0[/math] 

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

Módulo de la velocidad

Campo de presiones.

x=0:0.2:2;
z=0:0.2:10;
[X,Z]=meshgrid(x,z);
figure(1);
p=4-3/2*(Z-1);
surf(X,Z,p);
colorbar;
view(2);
axis([0,3,0,10]);
title('CAMPO DE PRESIONES');
xlabel('ρ');
ylabel('z');

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.

7 ROTACIONAL

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

Sustituyendo en la fórmula:

Obtenemos que [math]\nabla\times\vec{u}=(\frac{3\rho}{2}+\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].

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

Rotacional del campo velocidad

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

8 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:

Campo de temperaturas del fluido

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

Curvas de nivel de campo

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.

8.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:

Campo vectorial 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
hold off

Código Octave utilizado:

Gradiente ortogonal a las curvas de nivel

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.

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