Diferencia entre revisiones de «Flujo de Poiseuille (Grupo 36)»
(→Representación del campo de presiones.) |
(→Comportamiento del módulo de la velocidad.) |
||
| Línea 158: | Línea 158: | ||
{{matlab|codigo= | {{matlab|codigo= | ||
rho=0:0.1:2; | rho=0:0.1:2; | ||
| − | f= | + | f=((-3/4)*(rho.^2))+3; |
plot(rho,f); | plot(rho,f); | ||
title('Comportamiento del módulo del campo de velocidades'); | title('Comportamiento del módulo del campo de velocidades'); | ||
Revisión del 22:56 14 dic 2023
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Flujo de Poiseuille (Grupo 36) |
| Asignatura | Teoría de Campos |
| Curso | 2023-24 |
| Autores | Lucía Domínguez Álvarez; Eduardo Juarranz del Valle; Adrián Díaz Gadea; Pablo Amado Silva; Carmen Pardos Martínez |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Sección Longitudinal de la Tubería
- 3 Ecuación de Navier-Stokes y divergencia nula
- 4 Campo de presiones y campo de velocidades.
- 5 Líneas de corriente del campo.
- 6 Puntos con módulo de velocidad máxima.
- 7 Rotacional
- 8 Campo de temperaturas y curvas de nivel
- 9 Gradiente de temperatura
- 10 Caudal
1 Introducción
La ley de Poiseuille o también llamada ley de Hagen-Poiseuille es una ley que permite determinar el flujo laminar estacionario de un líquido incomprensible, en nuestro caso consideraremos el flujo de un líquido incomprensible en una tubería cilíndrica de radio 2, teniendo así una sección circular constante. Por lo tanto, dicho flujo vendrá determinado por el gradiente de presión y el radio.
Para la realización y creación de este artículo, hemos hecho uso del programa informático Matlab para poder representar gráficamente todos los resultados, como por ejemplo gráficos de secciones, gradientes... de una manera visual para así poder ayudar al usuario a comprender la Ley de Poiseuille ayudándole a comprender las interpretaciones de dicha ley.
2 Sección Longitudinal de la Tubería
Mallado de la represenación de la sección longitudinal de la tubería [math] x_{1} = 0 [/math], [math] \left ( \rho,z \right )\epsilon \left [ 0,3 \right ]\times \left [ 0,10 \right ]. [/math]
x=0:0.05:2; %Creamos Vectores
y=0:0.2:10;
[XX,YY]=meshgrid(x,y); %Creamos Malla
mesh(XX,YY,0*XX); %Representamos la sección
axis([0,3,0,10]); %Rango de los ejes
xlabel('ρ') ;
ylabel('z') ;
view(2);
title ('Malla de la Sección Longitudinal');
La velocidad de las partículas de nuestro fluido viene dada por el campo [math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}[/math], y la presión 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].
Ambas magnitudes, [math] \left ( \vec{u},\rho \right ) [/math], cumplen la ecuación estacionaria de Navier-Stokes, independiente del tiempo, donde [math] \mu [/math] es el coeficiente de viscosidad de fluido:- 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 } = \frac{\left ( p_{2}-p_{1} \right )}{\mu} \cdot \frac{\rho^{2}}{2} [/math]
- [math]\frac{\partial f\left ( \rho \right )}{\partial \rho } =\frac{p_{2}-p_{1}}{2\mu}\cdot \rho[/math]
- 3)Integramos por segunda vez
- [math] \int \left ( \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 )}{2\mu} \cdot \left (\frac{\rho^{2}}{2} + c \right ) [/math]
Para darle valor a la constante, usamos el dato de que en [math] \rho = 2 [/math] la velocidad es cero, por tanto, como la velocidad es [math] f\left ( \rho \right )\vec{e_{z}}, [/math] entonces [math] f\left ( 2 \right ) [/math] debe ser cero. La [math] f\left ( \rho \right ) [/math] nos queda: [math] f\left(\rho \right )=\frac{\left ( p_{2}-p_{1} \right )}{2\mu} \cdot \left (\frac{\rho^{2}}{2} - 2 \right ) [/math]
3.2 Divergencia Nula
Como el agua es un líquido incompresible, debe ocupar siempre el mismo volumen, es decir ser su densidad constante, para ello se comprueba que la divergencia del campo de velocidades es nula, ya que la divergencia de un campo de velocidades de un fluido en un punto se puede interpretar como el cambio en la densidad del fluido en ese punto.[math]\triangledown\cdot\vec{u} \left(\rho,\theta,z\right)=\frac{1}{\rho} \left ( \frac{\partial }{\partial \rho}\left ( \rho \cdot u_{\rho } \right )+\frac{\partial }{\partial \theta}\left ( u_{\theta } \right ) + \frac{\partial }{\partial z}\left (\rho \cdot u_{z}\right)\right) =\frac{1}{\rho}\left ( \frac{\partial }{\partial z}\left(\rho \cdot f\left(\rho\right)\right)\right)=0 [/math] (para cualquier valor de [math]\rho[/math]).
4 Campo de presiones y campo de velocidades.
A continuación calcularemos el campo de presiones (campo escalar) y el campo de velocidades (campo vectorial). Suponiendo que [math] p_{1}=4, p_{2}=1 [/math] y [math]\mu=1. [/math] Donde [math] p_{1} [/math] es la presión en los puntos [math] z=4 [/math], [math] p_{2} [/math] la presión en los puntos [math] z=1 [/math] y [math] \mu [/math] el coeficiente de viscosidad del fluido.
4.1 Campo de presiones.
Para calcular el Campo de presiones hacemos uso de la ecuación de presión; [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}=\frac{-3z}{2}+\frac{11}{2} [/math]
4.1.1 Representación del campo de presiones.
Para poder representar el campo de presiones debemos estudiar su comportamiento de la presión frente a la altura. Como podemos ver hay una relación lineal entre ambas. Cuanto más aumenta la profundidad, más aumenta la presión y de igual forma, si disminuye una, también lo hace la otra.
clc;
clear all;
z=0:0.1:10;
f=(-3)*z/2+(11/2);
plot(z,f)
xlabel('Variación de altura');
ylabel('Variación de presión');
title(' Gráfica del campo de presiones');
4.2 Campo de velocidades.
Para calcular el Campo de velocidades haremos uso de la ecuación que representa la velocidad de las partículas del fluido. Es decir por:
[math]\vec{u}(\rho,\theta,z)= f\left(\rho\right)\vec{e_{z}}= [/math] [math]\left(\frac{\left(p_{2}-p_{1}\right)}{4\mu}\rho^{2}+(\frac{\left(p_{2}-p_{1}\right)}{2\mu}.({-2}))\right)\vec{e_{z}} =\left(\frac{-3}{4}\rho^{2}+{3}\right)\vec{e_{z}}[/math]
4.2.1 Representación del campo de velocidades.
Para proceder con la representación del campo de velocidades debemos tener en cuenta que este es de carácter vectorial. Por lo que para representarlo usando Matlab o Octave deberemos pasar primero de coordenadas cilíndricas (ya que se trata de una tubería cilíndrica) a coordenadas cartesianas.
clear all;
clc;
x=0:0.1:3;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=-3/4*xx^2 + 3;%Cambio a Coord. Cartesianas
uy=-3/4*yy^2 + 3;%Cambio a Coord. Cartesianas
hold on % Debo rellenar la parte de arriba del código
quiver(xx,yy,ux,uy)
axis([0,5,0,12])
hold off
colorbar;
view(2)
title('GRÁFICA DEL CAMPO DE VELOCIDADES')
5 Líneas de corriente del campo.
Para dibujar las líneas de corriente del campo [math]\overrightarrow{u}[/math], es decir, las líneas tangenes a [math]\overrightarrow{u}[/math] en cada punto, es necesario calcular el campo [math]\overrightarrow{v}[/math] que es ortogonal a [math]\overrightarrow{u}[/math] en cada punto. Sabiendo que [math]\overrightarrow{v}=\overrightarrow{e_{\theta }}\times \overrightarrow{u}[/math], siendo el campo [math]\overrightarrow{v}[/math] irrotacional por ser nula la divergencia de [math]\overrightarrow{u}[/math] y ademas teniendo un potencial escalar [math]\psi [/math],([math]\bigtriangledown \psi =\overrightarrow{v}[/math]), siendo esta la función corriente de [math]\overrightarrow{u}[/math].
Primero calcularemos [math]\psi[/math] y dibujaremos las líneas de [math]\psi=cte[/math] y comprobaremos que son corrientes de [math]\vec{u}[/math] viendo 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)\overrightarrow{e_{\rho}}\rightarrow \overrightarrow{v}=f\left (\rho\right)\overrightarrow{e_{\rho}}[/math].
Sustituyendo [math]f\left(\rho\right)[/math] se obtiene:
[math]\overrightarrow{v}=\left(\frac{p2-p1}{4\mu}\rho^{2} + \frac{p1-p2}{\mu}\right)\overrightarrow{e_{\rho}}[/math].
Para calcular [math]\psi[/math], obtenemos su gradiente e igualamos al campo [math]\overrightarrow{v}[/math], [math]\bigtriangledown\psi=\overrightarrow{v}[/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].
Por lo tanto, [math]\frac{d\psi}{d\rho}=\frac{p2-p1}{4\mu}\rho ^{2} +\frac{p1-p2}{\mu}\ [/math].
Tras integrar nos queda finalmente que [math]\psi=\frac{p2-p1}{12\mu}\rho^{2} + \frac{p1-p2}{\mu}\rho[/math].
Sustituyendo en la anterior función los valores asignados, el resultado de la función potencial será: [math]\psi=-\frac{1}{4}\rho^{3} + 3\rho[/math]
5.1 Representación de las líneas de corriente del campo.
Mediante el programa informático Matlab se representan las líneas [math]\psi=cte [/math], para posteriormente comprobar que son líneas de corriente de [math]\overrightarrow{u}[/math]. Es decir, que son tangentes a dicho campo en cada punto. Siguen la dirección de [math]\overrightarrow{e_{z}}[/math] por lo que sigue el sentido de la tubería cilíndrica.
x=0:0.2:2; %generamos las lineas de corriente en x
y=0:0.2:10; %generamos las lineas de corriente en y
[X,Y]=meshgrid(x,y); %dibujamos el mallado
lineas=(-1./4).*X.^3+3.*X; %incluimos el resultado de la funcion potencial
contour(X,Y,lineas); %creamos un diagrama de contorno
axis([0,3,0,10]); %especificamos los limites de los ejes
colorbar %mostramos una barra de colores vertical a la derecha de los ejes
title('Líneas de corriente'); %escribimos el titulo
6 Puntos con módulo de velocidad máxima.
La ley de Poiseuille establece que el caudal de un fluido incompresible en un tubo cilíndrico de radio constante es proporcional a la diferencia de presión entre los extremos del tubo, al radio del tubo y a la inversa de la viscosidad del fluido. La velocidad del fluido en un punto del tubo se puede calcular a partir del caudal y del radio del tubo. El módulo de la velocidad del fluido es:
Para encontrar los puntos donde el módulo de la velocidad del fluido es máximo, debemos derivar la expresión anterior con respecto a ρ. Obtenemos:
La derivada es igual a cero cuando ρ = 0. Por lo tanto, el módulo de la velocidad del fluido es máximo en el punto ρ = 0, que corresponde al eje del tubo. Es decir, el módulo de la velocidad del fluido es máximo en el eje del tubo y su valor es proporcional al radio del tubo y a la inversa de la viscosidad del fluido.
6.1 Comportamiento del módulo de la velocidad.
La siguiente gráfica muestra el comportamiento del módulo de la velocidad del fluido en función de ρ. Se muestra la curva de la variación de la velocidad del fluido en función del radio de la sección de la tubería cilíndrica. Podemos ver como a medida que aumenta el radio, disminuye la velocidad. Por lo que concluimos que la velocidad máxima se situa en el eje de la tubería.
rho=0:0.1:2;
f=((-3/4)*(rho.^2))+3;
plot(rho,f);
title('Comportamiento del módulo del campo de velocidades');
xlabel('Radio de la sección');
ylabel('Variación de la velocidad');
axis([0,2,0,3.5])
7 Rotacional
El rotacional de un campo vectorial nos enseña la tendencia que tiene este campo vectorial a girar en torno a un punto. Al tener el fluido un coeficiente de viscosidad el fluido no es perfecto lo que conlleva que al rotar no van a ser todas la velocidades iguales y el rotacional no va ha ser nulo. Para llevar a cabo el rotacional del campo [math]\overrightarrow{u}[/math] utilizaremos la fórmula:
Con el campo vectorial que sacamos del apartado 5 el rotacional nos quedaría: [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}=-f'\left ( \rho \right )\overrightarrow{e_{\theta }}[/math]
La derivada de la función respecto de [math] \rho [/math] nos queda: [math]f'\left ( \rho \right )= \frac{\left ( p_{2}-p_{1} \right )\rho}{2\mu} [/math]
Sustituyendo los valores:El rotacional nos queda: [math] \triangledown \times \vec{u}\left ( \rho, \theta,z\right )=\frac{3\rho}{2} [/math]
Como se puede observar el rotacional del campo nos queda respecto de [math] \rho [/math], que a su vez este se encuentra en el intervalo cerrado [math] [0,2][/math].
Para visualizarlo mejor vamos a utilizar Matlab para visualizarlo gráficamente:
clear;
clc;
x=0:0.1:2;
y=0:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs((3./2).*xx);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])
title('Rotacional del campo');
hold off
Como el rotacional se quedaba en función de [math] \rho [/math] los valores máximos del rotacional se encuentran en [math] \rho=2 [/math]. En el gráfico los puntos con menor tendencia a rotar se van a encontrar con unos colores fríos (azul) y los puntos con mayor rotación con colores más cálidos (naranja/amarillo). Como se puede observar en el gráfico los colores cálidos se encuentran en las paredes de [math] \rho=2 [/math].
8 Campo de temperaturas y curvas de nivel
El campo de temperaturas es un campo escalar ya que cada punto tiene una temperatura diferente al resto. Suele depender de la posición del punto y del tiempo para determinar la temperatura de ese punto. La temperatura del fluido nos viene determinada por el campo escalar:[math] T(ρ{,}θ{,}z)=log(1+ρ)e^{-(z-2)^2}[/math]
Para determinar donde la temperatura es máxima gráficamente vale con ver el gráfico de matlab, también se puede calcular con el gradiente igualando al 0. En este apartado no es necesario pasar el campo a coordenadas cartesianas, con lo cual trabajaremos con las coordenadas polares.
8.1 Campo de temperaturas
x=0:0.1:2;
y=0:0.1:10;
[X,Y]=meshgrid(x,y);
figure(1);
hold on
T=log(X+1).*(exp(-(Y-2).^2));
pcolor(X,Y,T);
shading flat
grid on
axis([0,3,0,10]);
colorbar;
title('Campo de temperaturas')
hold off
Como se puede observar la temperatura sube cuando se encuentra entre [math] z\lt3 [/math] y [math] z\gt2 [/math] pero cuando mas avanza a la derecha de [math] \rho [/math] mayor es, por lo tanto el máximo se encuentra en las inmediaciones de [math] z=2 [/math] y [math] \rho=2 [/math].
8.1.1 Curvas de nivel de campo
Las curvas de nivel de un campo de temperaturas unen los puntos que se encuentran a la misma temperatura dentro del campo estas curvas son conocidas como isotermas. Estas las vamos a representar por colores las temperaturas al igual que el campo de temperaturas siendo los colores azules los más fríos y los amarillos/rojos los más cálidos.
x=0:0.2:2;
y=0:0.2:10;
[X,Y]=meshgrid(x,y);
hold on
T=log(X+1).*(exp(-(Y-2).^2));
contour(X,Y,T,10);
axis([0,3,0,10]);
view(2);
title('Curvas de nivel de temperatura')
colorbar
hold off9 Gradiente de temperatura
El gradiente de la temperatura será:
[math] \bigtriangledown T\left ( \rho, \theta, z \right )= e^{-\left (z-1 \right )^{2}} \left(\left(2 ρ - 1\right)\vec{e_ρ}-\left(2z-2\right) \left(ρ^{2}-\frac{1}{4}\right)\vec{e_z}\right) [/math]
Para ver el gradiente gráficamente usaremos matlab:
x=0:0.2:2;
y=0:0.2:10;
[X,Y]=meshgrid(x,y);
figure(1);
rho=sqrt(X.^2+Y.^2);
T=1+((rho-(1/2)).^2).*(exp(-(Y-1).^2));
[TX,TY]=gradient(T);
hold on
quiver(X,Y,TX,TY)
axis([0,3,0,10]);
title('Gradiente de temperatura');
shading flat
grid on
hold off
A continuación se expone gráficamente la ortogonalidad del gradiente de temperatura con las curvas de nivel:
10 Caudal
Llamaremos Q al caudal que atraviesa una sección de la tubería y lo podremos calcular gracias a la siguiente integral:
[math]Q\left ( m^{3}/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left ( \frac{1}{4}\frac{p2-p1}{\mu }\rho ^{2}+\frac{p1-p2}{\mu }\rho \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}} dS=...=\frac{\Pi }{3}\left ( \frac{p2-p1}{\mu } \right )+4\Pi \frac{p1-p2}{\mu }[/math].
Que tomando los valores mencionados anteriormente:
[math]p1=4, p2=1 ,\mu=1[/math],
[math]Q=11\pi\left ( m^{3}/s \right )[/math].
