Diferencia entre revisiones de «Flujo de Poiseuille Grupo 30»

Revisión actual del 10:42 9 dic 2024

1 Introducción

La ley de Poiseuille, también conocida como ley de Hagen-Poiseuille, describe el flujo laminar estacionario de un líquido incompresible. En este caso, analizaremos el flujo de un líquido incompresible a través de una tubería cilíndrica con un radio de 3, lo que implica una sección transversal circular constante. Este flujo depende del gradiente de presión y del radio de la tubería.

Para desarrollar este análisis, hemos utilizado el software Matlab, que nos ha permitido representar gráficamente los resultados, como secciones transversales y gradientes, de manera visual. Esto facilita al lector una mejor comprensión de la Ley de Poiseuille, ayudándole a interpretar y entender sus implicaciones de forma clara y didáctica.

2 Dibujar un mallado 2D de la sección longitudinal

Líneas de corriente

Se trata de graficar la sección longitudinal de la tubería en coordenadas [math] x_{1} = 0 [/math], [math] \left ( \rho,z \right )\epsilon \left [ 0,4 \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,4,0,10]);  %Rango de los ejes
xlabel('ρ') ;
ylabel('z') ;
view(2);
title ('Malla de la Sección Longitudinal');

3 Resolver la ecuación diferencial para f(ρ)

3.1 Ecuación de Navier-Stokes

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)/4 [/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=5 [/math].

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 = 3 [/math] la velocidad es cero, por tanto, como la velocidad es [math] f\left ( \rho \right )\vec{e_{z}}, [/math] entonces [math] f\left ( 3 \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 )}{4\mu} \cdot \left (\frac{\rho^{2}}{1} - 9 \right ) [/math]

3.2 Verificación de la condición de incompresibilidad

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

[math]\triangledown\cdot\vec{u} \left(\rho,\theta,z\right)=\frac{1}{\rho}\left ( \frac{\partial }{\partial z}\left(\rho \cdot f\left(\rho\right)\right)\right)[/math]

[math]\triangledown\cdot\vec{u} \left(\rho,\theta,z\right)= \frac{1}{\rho} \frac{\partial}{\partial \rho} \left[ \frac{(p_2 - p_1)}{4\mu} \rho \left( \rho^2 - 9 \right) \right] = 0.[/math] (para cualquier valor de [math]\rho[/math]).

4 Campo de velocidades y campo de presiones

Vamos a dar como dato: p1=1, p2=6 y μ=1

4.1 Campo de velocidades

Anteriormente hemos obtenido la función de velocidad.

Vamos a sustituir los valores dados:

Como podemos observar, a medida que nos acercamos al borde de la tubería, entra menos velocidad mientras que en el centro, hay un mayor flujo de concentración de velocidades.

Líneas de corriente

% Crear la malla 2D
[X, Z] = meshgrid(x, z);
% Calcular el campo de velocidades
p1 = 1;        % Presión inicial
p2 = 6;        % Presión final
mu = 1;        % Viscosidad dinámica
ux = ((p2 - p1) / (4 * mu)) .* (X.^2 - 9);  % Componente en ρ (x)
uz = 0 .* Z;                                % Componente en z (constante)
% Corregir los valores de ux donde X^2 > 9 para evitar valores inválidos
ux(X.^2 > 9) = 0;
% Representación del campo de velocidades
figure;
hold on;
quiver(X, Z, ux, uz, 'b');  % Campo de velocidades con flechas
axis([0, 4, 0, 10]);        % Ajustar los límites del gráfico
xlabel('\rho');             % Etiqueta del eje ρ
ylabel('z');                % Etiqueta del eje z
title('Campo de velocidades');
grid on;
hold off;
view(2);

4.2 Campo de presiones

La expresión del campo de presiones viene dada como:

Campo de presiones

Como podemos comprobar, la presión no aumenta con el radio, sino con el parámetro "z".

% Parámetros del problema
p1 = 2;        % Presión inicial
p2 = 6;        % Presión final
% Intervalo de altura 'z'
z = 0:0.1:5;   % Altura en el eje z
% Cálculo del campo de presiones
f = p1 + ((p2 - p1) / 4) .* (z - 1);  % Campo de presiones
% Graficar el campo de presiones
figure;
plot(z, f, 'r', 'LineWidth', 1.5);  % Gráfica con línea roja
grid on;                            % Activar cuadrícula
xlabel('Altura (z)');               % Etiqueta del eje z
ylabel('Presión (p)');              % Etiqueta del eje p
title('Campo de Presiones');        % Título de la gráfica
xlim([0, 5]);                       % Limitar el eje z para claridad

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 }}\times \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])

Líneas de corriente

% Parámetros
rho = linspace(0, 5, 100); % Valores de rho (radio)
theta = linspace(0, 2*pi, 100); % Valores de theta (ángulo)
[R, T] = meshgrid(rho, theta); % Crear la malla en coordenadas polares
% Función de corriente
psi = (5/12)*R.^3 - (45/4)*R; 
% Conversión a coordenadas cartesianas para graficar
X = R .* cos(T);
Y = R .* sin(T);
% Gráfico de líneas de corriente
contour(X, Y, psi, 20, 'LineWidth', 1.5); % Dibujar las líneas de nivel (líneas de corriente)
colorbar; % Barra de color para las líneas de corriente
title('Líneas de corriente');
xlabel('x (m)');
ylabel('y (m)');
axis equal; % Mantener proporción
grid on;

6 velocidad máxima del fluido y su gráfica de comportamiento .

Para encontrar los puntos donde la velocidad del fluido es maxima, derivamos respecto de ρ:

Si [math] \frac{\partial \vec{u} }{\partial ρ}=0 [/math] entonces ρ=0, que corresponde al eje del tubo.Es decir, la velocidad maxima en el eje del tubo.

La siguiente gráfica es la comportamiento de la velocidad,la cual nos muestra el comportamiento en ρ=0, en la que se obseva que cuanto mas nos acerquemos alos bordes de la tuberia la velocidad disminuye, entonces deducimos que la velocidad máxima es en el eje

Comportamiento módulo de máxima velocidad

x=0:0.1:3;
f=((5/4)*(x.^2))-9; 
plot(x,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,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. Para llevar a cabo el rotacional del campo [math]\overrightarrow{u}[/math] utilizaremos la fórmula:

Con el campo vectorial que sacamos del apartado 3 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 & \left(\frac{5}{4}\rho^{2}-{9}\right)\vec{e_{z}} \end{vmatrix}=\left(\frac{-5}{2}\rho\right)\vec{e_{\theta}} [/math]

El rotacional nos queda: [math]\triangledown\times\vec{u}\left(\rho,\theta,z\right)=\frac{-5\rho}{2}[/math]

Para poder ver la gráfica vamos a utilizar Matlab para visualizarlo gráficamente:[math]|\triangledown\times\vec{u}|\left(\rho,\theta,z\right)=\frac{5\rho}{2}[/math]

Rotacional del Campo

clear;
clc;
x=0:0.1:2;
y=0:0.1:10;
[x,y]=meshgrid(x,y);
rot=abs((5./2).*x);
surf(x,y,rot)
colorbar
view(2)
axis([0,3,0,10])
title('Rotacional del campo');
hold off

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.

8 Temperatura del Fluido

La Temperatura del fluido viene dada por [math] T(ρ,θ,z)=e^{(1+ρ)}2-(z-2)^2. [/math]

El campo de temperaturas las curvas de nivel los podemos observar en la siguiente gráfica realizada en Matlab.

En la que observamos gracias a los tonos más cálidos que la temperatura es máxima en el ρ=2 y z=2.

rho=0:0.01:2;
z=0:0.05:10;
[X,Y]=meshgrid(rho,z);
figure (1)
p=(exp(1+X).*2) - (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]);

izquierda

derecha

9 Gradiente de la temperatura

El gradiente de la temperatura lo calcularíamos mediante la siguiente fórmula:

Pero ya que nos lo piden gráficamente lo obtenemos mediante Matlab, y observamos que el gradiente es ortogonal a las curvas de nivel, gracias al siguiente código.

derecha

rho=0:0.1:2; 
z=0:0.1:10;
[X,Y]=meshgrid(rho,z);
figure(1)
T=(exp(1+X).*2) - (Y-2).^2;
[TRHO,TZ]=gradient(T); 
hold on
quiver(X,Y,TRHO,TZ)
contour(X,Y,T,'k') 
axis([0,2,0,10]);
title('Gradiente de Temperatura y Curvas de nivel')
shading flat
grid on
hold off

10 Caudal

El caudal es la cantidad de un fluido que atraviesa la tubería, rio, canal ,…, en un determinado tiempo. Para ello debemos tener en cuenta el flujo de volumen por unidad de tiempo.

Para poder calcular el caudal que llevara la tubería, se estudia el volumen de fluido que pasa a través de la tubería. La integral de la superficie es la siguiente:

donde [math]\vec{e_{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:

Sustituimos los siguientes valores:

[math]Q\left(m/s\right)=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left (\frac{p2-p1}{4\mu}\right) \left(\frac{\rho^2}{1} -9\right)\cdot \rho \overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}} dS=\int_{0}^{2\Pi }\int_{0}^{3}\left (\frac{p2-p1}{4\mu}\right) \left(\frac{\rho^2}{1} -9\right)\cdot \rho d_{\rho }d_{\theta }=\int_{0}^{2\Pi }\int_{0}^{3} (\rho^{2}-9 ) \cdot \rho =(-81/2)\cdot \prod_{}^{} =-127,23\left ( m/s \right )[/math].

El caudal nos a salido negativo, lo cual significa que interpolando el caudal va en sentido contrario, esto depende de nuestra normal y de nuestro signo de la resultante.

Según los datos del problema, p1=2 en z=1 y p2=6 en z=5, esto significa que la presión es mayor en z=5 y menor en z=1. El fluido siempre se mueve de mayor presión a menor presión. En este caso, el flujo se dirige desde z=5 hacia z=1 es decir, en el sentido contrario al eje z. Por lo tanto, la componente de velocidad será negativa por lo cual va en sentido contrario .