Flujo de Couette entre dos tubos concéntricos (Grupo 25)

De MateWiki
Revisión del 19:54 6 dic 2025 de Andres.gaitan (Discusión | contribuciones) (→‎Líneas de corriente)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Flujo de Couette entre dos tubos concéntricos . Grupo 25
Asignatura Teoría de Campos
Curso 2025-26
Autores Andrés Gaitán Tsukasov, Pablo Casado del Campo, Esteban Gabaldón Hermoso , Carlos Ascanio Martín, Aarón Pérez Luna
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Flujo de Couette entre dos tubos concéntricos

Introducción

El flujo de Couette entre dos tubos concentricos, designa el movimiento de un fluido viscoso entre dos cilindros coaxiales de diferente radio, siendo dominado así el movimiento por la rotación de los cilindros alrederor de una eje común. En ingenieria, el flujo de Couette es sumamente importante ya que es uno de los pocos fluidos viscosos de los que se conoce la solución exacta de la ecuacion de Navier-Stokes, volviendolo fundamental en el estudio de propiedades de los fluidos y estabilidad de los flujos.

En este trabajo partiremos de un flujo generado por la rotación en sentido antihorario del cilindro interior con velocidad angular constante [math]\vec{w}_i[/math] y la rotación en sentido horario del cilindro exterior con velocidad angular [math]\vec{w}_e[/math]. Ambos cilindros tienen su eje [math]OX_3[/math] y sus radios [math]\rho[/math]=2 para el cilindro exterior y [math]\rho[/math]=1 para el interior. A partir de ello, con el campo de velocidades entre ambos cilindros, estudiaremos el flujo de temperaturas resultante, el gradiente y el caudal, entre otras.

1 Mallado de la sección transversal

En primer lugar, se toma la sección de ambos cilindros en el plano [math]x_3[/math]=0 y se representa en matlab de acuerdo al parámetro angular [math]\theta \in [0,2\pi][/math] y a los parámetros cartesianos [math](x,y) \in [-2,2] \times [-2,2][/math].

Leyenda: Malla de elementos finitos para el flujo de Couette.
% MALLADO DE LA SECCIÓN TRANSVERSAL
figure; 
rho=1:0.1:2;                                %SE ESTABLECEN PARÁMETROS                     
theta_malla = linspace(0, 2*pi, 50);             
[RHO,THETA]=meshgrid(rho,theta_malla);      %OBTENCION DE LA MALLA       
rho_int=1;
rho_ext=2;
hold on
grid on
x=RHO.*cos(THETA);                               
y=RHO.*sin(THETA);                               
mesh(x,y,0*y,'EdgeColor','g');              
x_ext = rho_ext.*cos(theta_malla);          %OBTENCIÓN DE COORDENADAS
y_ext = rho_ext.*sin(theta_malla);
x_int = rho_int.*cos(theta_malla); 
y_int = rho_int.*sin(theta_malla);
axis([-2, 2, -2, 2]);                       %AJUSTE DE LOS EJES           
axis equal;                                      
plot(x_ext, y_ext, 'b', 'LineWidth', 3);    %SE REPRESENTA EL CILINDRO EXTERIOR
plot(x_int, y_int, 'b', 'LineWidth', 3);    %SE REPRESENTA EL CILINDRO EXTERIOR      
view(2)                                     %VISTA EN 2D 
hold off


2 Calculo de velocidades. Ecuación de Navier-Stokes

2.1 Velocidad de las partículas del fluido y particularización de la ecuación

La velocidad de las partículas del fluido viene dada por la expresión [math]\vec{u}(\rho, \theta) = f(\rho)\vec{e}_{\theta}[/math], donde la presión p es constante. Sabiendo que (\vec{u},p) cumple la ecuacion de Navier-Stokes estacionaria se obtiene:

[math](\vec{u} \cdot \nabla)\vec{u} + \nabla p = \mu \Delta\vec{u}[/math]

Siendo \mu el coeficiente de viscosidad del fluido.

Además, como la presión p es constante, su gradiente sera nulo: [math]\nabla p = 0[/math], y si aparte despreciamos el primer término de la fórmula anterior obtenemos que:

[math]µ∆\vec{u} =\vec{0}[/math]

Así, se llega al Laplaciano de un campo vectorial que se desarrollara en la base cilíndrica de acuerdo a la siguiente expresión:

[math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math]

2.2 Cálculo del gradiente de una divergencia

Se desarrolla el primer termino, que es el gradiente de una divergencia en coordenadas cilindricas:

[math]\nabla \cdot \vec{u} = \frac{1}{\rho} \frac{\partial}{\partial \rho}(\rho u_{\rho}) + \frac{1}{\rho} \frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial u_{z}}{\partial z}[/math]

Si se particulariza al ejercicio se obtiene que:

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

Llegandose asi a la conclusión de que el primer término es igual al vector nulo: [math]\nabla(\nabla \cdot \vec{u}) = \vec{0}[/math]

2.3 Cálculo del rotacional del rotacional

Ahora se desarrolla el segundo termino, el rotacional del rotacional de [math]\vec{u}[/math], para ello se empieza en primer lugar con el primer rotacional:

[math]\nabla \times \vec{F}(\rho, \theta, z) = \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_{\rho} & \rho F_{\theta} & F_{z} \end{vmatrix}[/math]

Así, particularizando al ejercicio se obtiene:

[math]\nabla \times \vec{u} = \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} = \left( \frac{1}{\rho} \frac{\partial}{\partial \rho}(\rho f(\rho)) \right) \vec{e}_{z}[/math]

Tras ello, se desarrolla el siguiente rotacional:

[math]\nabla \times (\nabla \times \vec{u}) = - \frac{1}{\rho} \det\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 & \frac{1}{\rho} \frac{\partial}{\partial \rho} (\rho f(\rho)) & 0\end{vmatrix} = -\left[ \frac{\partial}{\partial \rho}\left(\frac{1}{\rho}\frac{\partial}{\partial \rho}(\rho f(\rho))\right) \right] \vec{e}_{\theta} = -\left[ \frac{\rho^2 \frac{\partial^2 f}{\partial \rho^2} + \rho \frac{\partial f}{\partial \rho} - f(\rho)}{\rho^2} \right] \vec{e}_{\theta}= -\left[ \frac{\partial^2 f}{\partial \rho^2} + \frac{1}{\rho} \frac{\partial f}{\partial \rho} - \frac{f(\rho)}{\rho^2} \right] \vec{e}_{\theta}[/math]

Por tanto, tras desarrollar todos los terminos se llega a que:

[math]\Delta\vec{u} = \left[ \frac{\partial^2 f}{\partial \rho^2} + \frac{1}{\rho} \frac{\partial f}{\partial \rho} - \frac{f(\rho)}{\rho^2} \right] \vec{e}_{\theta}[/math]

2.4 Solución de la ecuacion de Navier-Stokes de forma diferencial

Ya, una vez conocidos todos los términos de la ecuación de Navier-Stokes, se puede resolver:

[math](\vec{u} \cdot \nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} \Rightarrow 0 + \vec{0} - \mu \Delta \vec{u} = \vec{0} \Rightarrow \Delta\vec{u} = \vec{0}[/math]

De lo que se deduce la siguiente ecuación diferencial:

[math]\left[ \frac{\partial^2 f}{\partial \rho^2} + \frac{1}{\rho} \frac{\partial f}{\partial \rho} - \frac{f(\rho)}{\rho^2} \right] \vec{e}_{\theta} = 0[/math]

Primero se despeja [math]f(\rho)[/math] y luego se simplifica llegando a la ecuación final:

[math]\frac{\partial^2 f}{\partial \rho^2} + \frac{1}{\rho} \frac{\partial f}{\partial \rho} - \frac{f(\rho)}{\rho^2} = 0 \quad \Rightarrow \quad \frac{1}{\rho} \frac{\partial}{\partial \rho} \left[ \rho \frac{\partial f}{\partial \rho} \right] = \frac{f(\rho)}{\rho^2} [/math]
[math]\frac{f(\rho)}{\rho} = \frac{\partial}{\partial \rho}\left(\rho \frac{\partial (f(\rho))}{\partial \rho}\right)[/math]

2.5 Comprobación de la ecuación obtenida

A continuación se comprueba si la ecuación obtenida cumple la solución obtenida:

[math]f(\rho) = a\rho + \frac{b}{\rho}, \quad a, b \in \mathbb{R}[/math]

Para ello hay que derivar la expresión obtenida:

[math]\frac{\partial}{\partial \rho}\left(\rho \frac{\partial f(\rho)}{\partial \rho}\right) = \frac{\partial}{\partial \rho}\left(\rho \left(a - \frac{b}{\rho^2}\right)\right) = \frac{\partial}{\partial \rho}\left(a\rho - \frac{b}{\rho}\right) = a + \frac{b}{\rho^2}, \quad a, b \in \mathbb{R}[/math]

Tras haber realizado la derivada y sustituir en la ecuación, se comprueba que la ecuación cumple la solución obtenida:

[math]\frac{f(\rho)}{\rho} = \frac{\partial}{\partial \rho}\left(\rho \frac{\partial (f(\rho))}{\partial \rho}\right) \implies \frac{1}{\rho}\left(a\rho + \frac{b}{\rho}\right) = a + \frac{b}{\rho^2}, \quad a, b \in \mathbb{R}[/math]

2.6 Calculo de los parámetros a y b

Ahora queda determinar a y b para que las velocidades en la frontera del fluido coincidan con las de los cilindros. Para ello, se aplican las condiciones obtenidas sobre los parámetros obtenidos. Primero se particulariza en el cilindro interior, donde [math]\rho[/math]=1:

[math]\vec{u} = \omega_i \vec{e}_\theta = f(1)\vec{e}_\theta = \left(a \cdot 1 + \frac{b}{1}\right) \vec{e}_\theta \implies a + b = \omega_i[/math]

Segundo se particulariza en el cilindro exterior, donde [math]\rho[/math]=2:

[math]\vec{u} = -\omega_e \vec{e}_\theta = -f(2)\vec{e}_\theta = -\left(2a + \frac{b}{2}\right)\vec{e}_\theta \implies 2a + \frac{b}{2} = -\omega_e[/math]

Solucionando el sistema en función de las velocidades angulares se obtiene que:

[math]b = \frac{2}{3} \omega_e + \frac{4}{3} \omega_i[/math]
[math]a = -\frac{2}{3} \omega_e - \frac{1}{3} \omega_i[/math]

Tras averiguar los valores de a y b, [math]\vec{u}[/math] queda así:

[math]\vec{u}(\rho, \theta) = \left[\left(\left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right)\rho\right) + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho}\right] \vec{e}_{\theta}[/math]

2.7 Condición de incompresibilidad

Para comprobar la condición de incompresibilidad, se hace la divergencia del [math]\vec{u}[/math] obtenido en coordenadas cilindricas y se comprueba que la divergencia del campo de velocidades es nulo:

[math]\vec{u}(\rho, \theta) = \left[\left(\left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right)\rho\right) + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho}\right] \vec{e}_{\theta} \\\nabla \cdot \vec{u} = \frac{1}{\rho}\left(\frac{\partial}{\partial \rho}(\rho u_{\rho}) + \frac{\partial}{\partial \theta}(u_{\theta}) + \frac{\partial}{\partial z}(\rho u_z)\right) = \frac{1}{\rho}\left(0 + \frac{\partial}{\partial \theta}(u_{\theta}) + 0\right) = 0[/math]


3 Representación del campo de velocidades

Partiendo de las siguientes expresiones; [math]\vec{u}(\rho,\theta) = f(\rho)\,\vec{e}_\theta,\quad f(\rho) = a\rho + \frac{b}{\rho}[/math]; y los valores anteriormente calculados; [math]b = \frac{2}{3} \omega_e + \frac{4}{3} \omega_i[/math], [math]a = - \frac{2}{3} \omega_e - \frac{1}{3} \omega_i[/math], si suponemos que [math]|\vec{\omega}_e| = |\vec{\omega}_i| = 1 \text{ y } \mu = 1[/math], podemos representar graficamente nuestro campo de velocidades:

Leyenda: Representación gráfica campo de velocidades.
% Parametros
Ri = 1;            % Radio interno
Ro = 2;            % Radio externo
Nr = 28;           % puntos radiales
Ntheta = 40;       % puntos angulares
omega1 = 1;        % velocidad angular interior
omega2 = -1;       % velocidad angular exterior

% Mallado polar
r = linspace(Ri, Ro, Nr);
theta = linspace(0, 2*pi, Ntheta);
[R, Theta] = meshgrid(r, theta);
X = -R .* cos(Theta);    % signo negativo para mantener estilo del ejemplo
Y =  R .* sin(Theta);

% Campo de velocidades (Couette)
A = (omega2*Ro - omega1*Ri) / (Ro - Ri);
B = (omega1*Ro - omega2*Ri) * Ri * Ro / (Ro - Ri);
Vtheta = A .* R + B ./ R;
Vx = -Vtheta .* sin(Theta);
Vy =  Vtheta .* cos(Theta);

% Dibujo del campo (color cambiado a azul intenso)
quiver(X, Y, Vx, Vy, 0.6, 'Color', [0 0.4 0.9], 'LineWidth', 1.3);
hold on;

% Lineas negras que simulan los cilindros (contornos en Ri y Ro)
t = linspace(0, 2*pi, 400);
plot(Ri*cos(t), Ri*sin(t), 'k-', 'LineWidth', 1.2);  % cilindro interior
plot(Ro*cos(t), Ro*sin(t), 'k-', 'LineWidth', 1.2);  % cilindro exterior

% Ajuste de ejes
padding = 0.25 * (Ro - Ri);
xlim([-(Ro+padding), Ro+padding]);
ylim([-(Ro+padding), Ro+padding]);
axis equal;

% Texto (cambiado respecto al original)
title('Distribución del campo de velocidades');
xlabel('Coordenada X');
ylabel('Coordenada Y');

hold off;


4 Líneas de Corriente

4.1 Campo v ortogonal a u

[math]\vec{v}(\rho, \theta) = \left[ \left( \left(-\frac{2}{3}w_{e}-\frac{1}{3}w_{i}\right)\rho \right) + \left( \frac{2}{3}w_{e}+\frac{4}{3}w_{i} \right) \frac{1}{\rho} \right] \vec{e}_{\rho}[/math]


4.2 Irrotacionalidad de v

[math]\text{Campo vectorial: } \vec{v}(\rho, \theta) = f(\rho)\vec{e}_{\rho} \implies v_{\rho} = f(\rho), \ v_{\theta} = 0, \ v_{z} = 0[/math]


[math]\text{Fórmula del Rotacional en Coordenadas Cilíndricas:}[/math]

[math]\nabla \times \vec{v} = \frac{1}{\rho}\left[\frac{\partial v_{z}}{\partial \theta} - \frac{\partial (\rho v_{\theta})}{\partial z}\right]\vec{e}_{\rho} + \left[\frac{\partial v_{\rho}}{\partial z} - \frac{\partial v_{z}}{\partial \rho}\right]\vec{e}_{\theta} + \frac{1}{\rho}\left[\frac{\partial (\rho v_{\theta})}{\partial \rho} - \frac{\partial v_{\rho}}{\partial \theta}\right]\vec{e}_{z}[/math]


\text{Sustitución de Componentes y Derivadas:}

[math]* \text{Componente } \vec{e}_{\rho}: \frac{1}{\rho}\left[\frac{\partial (0)}{\partial \theta} - \frac{\partial (\rho \cdot 0)}{\partial z}\right] = 0[/math]


[math]* \text{Componente } \vec{e}_{\theta}: \left[\frac{\partial f(\rho)}{\partial z} - \frac{\partial (0)}{\partial \rho}\right] = 0 - 0 = 0 \quad (\text{porque } f(\rho) \text{ no depende de } z)[/math]


[math]* \text{Componente } \vec{e}_{z}: \ltmath\gt\frac{1}{\rho}\left[\frac{\partial (\rho \cdot 0)}{\partial \rho} - \frac{\partial f(\rho)}{\partial \theta}\right] = \frac{1}{\rho}[0 - 0] = 0 \quad (\text{porque } f(\rho) \text{ no depende de } \theta)[/math]


\text{Conclusión: El campo es irrotacional.}

[math]\nabla \times \vec{v} = \vec{0}[/math]

4.3 Función de corriente de u (𝜓)

La función de corriente de 𝑢⃗ es el potencia escalar de dicho campo, lo que significa: ∇𝜓 = v


[math]\nabla \psi = \frac{\partial \psi}{\partial \rho} \vec{e}_{\rho} + \frac{1}{\rho} \frac{\partial \psi}{\partial \theta} \vec{e}_{\theta} = \vec{v}[/math]


[math]\frac{\partial \psi}{\partial \rho} = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho}[/math]


[math]\frac{1}{\rho} \frac{\partial \psi}{\partial \theta} = 0 \implies \frac{\partial \psi}{\partial \theta} = 0[/math]

Podemos ver que ψ depende únicamente de ρ, por lo que integramos la ecuación respecto a ρ


[math]\psi(\rho) = \int \left[ \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho} \right] d\rho[/math]


[math]\psi(\rho) = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \int \rho \, d\rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \int \frac{1}{\rho} \, d\rho[/math]


[math]\psi(\rho) = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \frac{\rho^2}{2} + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \ln(\rho) + C[/math]


El resultado sería:

[math]\psi(\rho, \theta) = \left(-\frac{1}{3}w_e - \frac{1}{6}w_i\right) \rho^2 + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \ln(\rho) + C[/math]


4.4 Comprobación de que son líneas de corriente

Como ya mencionamos anteriormente, para que v represente las líneas de correinte de u estos dos campos deben ser ortogonales entre sí.

[math]\vec{u} \cdot \nabla \psi = 0 \implies \vec{u} \cdot \vec{v} = 0[/math]


Componentes del campo tangencial [math]\vec{u}:[/math] [math]\quad u_{\rho}=0, \quad u_{\theta}=f(\rho), \quad u_{z}=0[/math]


Componentes del radial ortogonal [math]\vec{v}:[/math] [math]\quad v_{\rho}=f(\rho), \quad v_{\theta}=0, \quad v_{z}=0[/math]


La multiplicación escalar sin sustituir los valores quedaría tal que asi:


[math]\vec{u} \cdot \vec{v} = u_{\rho}v_{\rho} + u_{\theta}v_{\theta} + u_{z}v_{z}[/math]


Procedemos a sustituir las componentes por las aportadas anteriormente:


[math]\vec{u} \cdot \vec{v} = (0) \cdot (f(\rho)) + (f(\rho)) \cdot (0) + (0) \cdot (0)[/math]

[math]\vec{u} \cdot \vec{v} = 0 + 0 + 0 = 0[/math]


Como podemos ver se cumple la ortogonalidad ya que el resultado es = 0

4.5 Líneas de corriente

Archivo:Figura apartado 4.png 

%% Parámetros
Ri = 1; Re = 2;
we = -1; wi = +1; 

% Coeficientes de la función f(rho) = A_f*rho + B_f/rho
A_f = (-2/3)*we - (1/3)*wi; 
B_f = (2/3)*we + (4/3)*wi;

A_Psi = A_f / 2; % Coeficiente de R^2
B_Psi = B_f;     % Coeficiente de log(R)

%% Mallado y Cálculo de Psi
L = 3; N = 300;
rho = linspace(0.001, L, N); 
theta = linspace(0, 2*pi, N);
[R,T] = meshgrid(rho, theta);

X = R.*cos(T);
Y = R.*sin(T);

% Cálculo del potencial escalar
Psi = A_Psi*(R.^2) + B_Psi*log(R);

% Enmascaramiento de zonas fuera del anillo
mask = (R < Ri


5 Rotacional de u

6 Gradiente de la temperatura

6.1 =Representación grafica

En la figura 8 se puede apreciar como las curvas de nivel que representan la temperatura son ortogonales al gradiente.

Leyenda: Figura 8: Gradiente temperatura/curvas de nivel.
clear all % Limpiar la pantalla y espacio de trabajo
clc

r=1:0.1:2; % parámetros
th=0:0.1:2*pi;
[R,theta]=meshgrid(r,th);

temp=log10(1+R.^2).*(((cos(theta)).^2).*exp(-(R-(5/2)).^2)); % Expresión de la temperatura

[Dx,Dy]=gradient(temp,r,th); % Gradiente de la temperatura

figure; % Crea una figura y dibuja las curvas de nivel

quiver(R,theta,Dx,Dy,'color','r'); % Dibuja un campo de vectores sobre las curvas de nivel

xlabel('r'), ylabel('\theta'), colorbar;
title('curvas de nivel');

hold on
contour(R,theta,temp,'b');
legend('vectores', 'curvas de nivel'); % Leyenda para indicar el campod e vectores
hold off

6.2 Calculo del gradiente

El gradiente de la temperatura en coordenadas cilíndricas queda definido por la ecuacuacion:

[math]\nabla T(\rho,\theta) = \frac{\partial T}{\partial \rho}\,\vec e_{\rho} + \left(\frac{1}{\rho}\frac{\partial T}{\partial \theta}\right)\vec e_{\theta} + \frac{\partial T}{\partial z}\,\vec e_{z}[/math]

Cada coordenada tendrá su derivada parcial respecto a ella misma,las cuales sumadas son el gradiente:

[math]\quad (\nabla T)_{\rho}=\frac{\partial T}{\partial \rho},\quad (\nabla T)_{\theta}= \frac{1}{\rho}\frac{\partial T}{\partial \theta},\quad (\nabla T)_{z}= \frac{\partial T}{\partial z} [/math]

7 Caudal de una sección longitudinal

Para finalizar calcularemos el caudal que atraviesa los dos cilindros por una sección longitudinal producida por el plano x2=0 .

Para ello tendremos en cuenta que los dos cilindros tiene 1 metro de profundidad y la velocidad del fluido (dada en m/s) viene dada por la forma:

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Flujo_de_Couette_entre_dos_tubos_concéntricos_(Grupo_25)&oldid=102364»