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

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 Gaitan, Pablo Casado, Esteban Gabaldon, Carlos Ascanio, Aarón Perez
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].

% MALLADO DE LA SECCIÓN TRANSVERSAL
figure; 
rho=1:0.1:2;                                  % SE ESTABLECEN PARAMETROS ALGULARES Y CARTESIANOS                    
theta_malla = linspace(0, 2*pi, 50);             
[RHO,THETA]=meshgrid(rho,theta_malla);        % SE OBTIENE 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');                % SE REPRESENTA LA MALLA
theta_suave = theta_malla; 
x_ext = rho_ext.*cos(theta_suave);            % SE OBTIENEN LOS CILINDROS INTERIOR Y EXTERIOR EN PLANTA
y_ext = rho_ext.*sin(theta_suave);
x_int = rho_int.*cos(theta_suave); 
y_int = rho_int.*sin(theta_suave);
axis([-2, 2, -2, 2]);                         % SE LIMITAL LOS EJES                  
axis equal;                                      
plot(x_ext, y_ext, 'b', 'LineWidth', 3);      % SE REPRESENTAN AMBOS CILINDROS  
plot(x_int, y_int, 'b', 'LineWidth', 3);         
view(2)                                       % MALLADO DE LA SECCIÓN TRANSVERSAL - MALLA CIRCULAR Y BAJA DENSIDAD
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], llegamos a nuestro campo de velocidades:

[math]\vec u(\rho, \theta)=f(\rho)\vec e_\theta=\left ( -\rho -\frac{2}{\rho}\right )\vec e_\theta[/math]

4 Líneas de Corriente

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

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

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

[math]\psi(\rho, \theta) = \int v_\rho \, d\rho[/math]

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