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

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos. Grupo 67
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Carazo Ureña, Laura García de veas González, Marcos Molina Amigo, Pablo Ronchas Martin, Luis Alfonso
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El flujo de Couette describe el movimiento de un fluido viscoso confinado entre dos superficies que se desplazan tangencialmente una respecto de la otra. En geometría cilíndrica, este fenómeno aparece cuando el fluido ocupa el espacio entre dos cilindros concéntricos de gran longitud y el movimiento está dominado por la rotación alrededor del eje común. Esta configuración se utiliza con frecuencia como modelo sencillo para estudiar cómo se transmiten los esfuerzos cortantes en lubricantes, rodamientos y dispositivos de medida de propiedades reológicas.​

En el problema que se aborda en este trabajo, el fluido se encuentra entre dos tubos concéntricos que giran con velocidades angulares constantes y de sentido opuesto. Bajo la hipótesis de fluido incompresible, régimen laminar y simetría axial, el campo de velocidades resultante es puramente azimutal y depende únicamente del radio, determinado por la condición de no deslizamiento en ambas paredes cilíndricas. Estas condiciones de contorno contrapuestas generan un gradiente de velocidad más intenso en el espesor del anillo, lo que hace especialmente interesante analizar la distribución de temperatura, las curvas de nivel del campo escalar y el comportamiento del gradiente, así como discutir el significado físico de las magnitudes obtenidas y su posible relación con regímenes más complejos de tipo Taylor–Couette a mayores velocidades de rotación.

1 Mallado de la sección transversal

En primer lugar se representa la sección transversal del fluido correspondiente a los parámetros [math]\rho = 2[/math] para el tubo exterior, [math]\rho = 1[/math] para el tubo interior y [math]\theta \in [0,2\pi][/math]. Seccionamos ambos cilindros con el plano [math]x_3 = 0[/math] y describimos la región resultante en coordenadas cartesianas como [math](x,y) \in [-2,2] \times [-2,2][/math].

La sección se representa mediante el siguiente código de matlab.

% Parámetros
Ri = 1; % Radio interior
Ro = 2; % Radio exterior
Nr = 50; % Número de nodos en dirección radial
Ntheta = 100; % Número de nodos en dirección angular

% Vectores de coordenadas en polares
r = linspace(Ri, Ro, Nr);
theta = linspace(0, 2*pi, Ntheta);
% Construir el mallado con coordenadas polares
[R, Theta] = meshgrid(r, theta);

% Convertir a coordenadas cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);
% Graficar la malla interna en azul
figure;                                      
plot(X, Y, 'b');
hold on;
plot(X', Y', 'b');

plot(Ro * cos(theta), Ro * sin(theta), 'k', 'LineWidth', 3);   % Exterior
plot(Ri * cos(theta), Ri * sin(theta), 'k', 'LineWidth', 3);   % Interior
                                                               % Encuadrar y ampliar espacio en los ejes
padding = 0.25 * (Ro - Ri);                                    % Margen extra
xlim([-Ro - padding, Ro + padding]);
ylim([-Ro - padding, Ro + padding]);
axis equal;
set(gca, 'Position', [0.13 0.13 0.75 0.75]); % Centra la figura y la hace amplia

title('Mallado del dominio entre tubos concéntricos');
xlabel('x');
ylabel('y');
hold off;

5 Velocidad máxima del fluido

El estudio de los puntos donde el módulo del campo de velocidades [math]\vec{u}(\rho,\theta)=\left[\left(-\tfrac{1}{3}\omega_i-\tfrac{2}{3}\omega_e\right)\rho+\tfrac{\tfrac{4}{3}\omega_i+\tfrac{2}{3}\omega_e}{\rho}\right]\vec{e}_{\theta}[/math] , es máximo, se abordará mediante dos condiciones simplificadoras, que agilizan el problema:

Primera, puesto que en la expresión de [math]\vec{u}(\rho,\theta)[/math] solo aparece la variable [math]\rho[/math], puede tratarse el campo de velocidades como un campo vectorial dependiente de una única variable, así: [math]\vec{u}(\rho,\theta)=\vec{u}(\rho)[/math].

Segunda, como la función vectorial devuelve vectores dirigidos únicamente según la dirección [math]\vec{e}_{\theta}[/math], el resultado numérico del campo resulta ser el módulo de la velocidad.

Así pues, transformamos la expresión vectorial [math]\vec{u}(\rho,\theta)[/math] en una función escalar: [math] u(\rho) = \left| \left[\left(-\tfrac{1}{3}\omega_i-\tfrac{2}{3}\omega_e\right)\rho + \tfrac{\tfrac{4}{3}\omega_i+\tfrac{2}{3}\omega_e}{\rho} \right]\vec{e}_{\theta} \right| [/math] de [math]\mathbb{R}[/math] en [math]\mathbb{R}[/math]. Si además continuamos con el supuesto de que [math]|\omega_i|=|\omega_e|=1[/math], queda finalmente como:

[math] u(\rho)=-\,\rho+\frac{2}{\rho}, [/math]

y representa el módulo del campo de velocidades, incluyendo en su signo información sobre el sentido de la velocidad.

5.1 Estudio de la función [math] u(\rho)=-\,\rho+\frac{2}{\rho} [/math]

Para encontrar los puntos de máxima velocidad, puede estudiarse analíticamente hallando los valores de [math]\rho[/math] tales que su derivada se anule, y analizar ahí, junto con los valores en los extremos del dominio, cuál de ellos proporciona los valores máximos. Procedamos a su estudio en el intervalo [math]\rho\in[1,2][/math].

• Cálculo de la derivada:

[math] \frac{\partial u(\rho)}{\partial\rho} = u'(\rho) = -1 - \frac{2}{\rho^{2}}. [/math]

Para [math]\rho\gt0[/math], el término [math]\tfrac{2}{\rho^{2}}\gt0[/math], luego [math]u'(\rho)[/math] es la suma de dos términos negativos; por tanto:

[math] u'(\rho)\lt0,\qquad \forall\rho\in[1,2]. [/math]

La función es estrictamente decreciente en ese intervalo, característica que se confirma ante la inexistencia de puntos críticos, como demuestra:

[math] u'(\rho)=0 \Longrightarrow -1-\frac{2}{\rho^{2}}=0 \Longrightarrow \frac{2}{\rho^{2}}=-1 \Longrightarrow \rho^{2}=-2, [/math]

que no tiene solución real.

• Valores en los extremos del intervalo:

[math] u(1)=-1+\frac{2}{1}=1, \qquad u(2)=-2+\frac{2}{2}=-1. [/math]

Se concluye así que, como [math]u'(\rho)\lt0[/math] en todo [math]\rho\in[1,2][/math], el módulo de la velocidad es estrictamente decreciente en el espacio entre sendos tubos; así pues, el máximo se alcanza en [math]\rho=1[/math] con [math]u(1)=1[/math], y el mínimo en [math]\rho=2[/math] con [math]u(2)=-1[/math].

Por otro lado, si la función cambia su signo entre un extremo y otro, y es continua en [math][1,2][/math], el Teorema de Bolzano garantiza la existencia de un [math]\rho_0[/math] tal que: [math] u(\rho_0)=0. [/math]

Resolviendo:

[math] u(\rho)=0 \Longrightarrow -\rho+\frac{2}{\rho}=0 \Longrightarrow \rho=\frac{2}{\rho} \Longrightarrow \rho^{2}=2 \Longrightarrow \rho=\sqrt{2}\approx 1,41. [/math]

5.2 Interpretación física del resultado

No obstante, no hay que olvidar que el resultado obtenido es puramente matemático, y es necesario verlo desde un punto de vista físico para no interpretarlo erróneamente. La función estudiada analíticamente representa el módulo del campo de velocidades del fluido que se encuentra entre los dos tubos, de radios uno y dos, y proviene de calcular el módulo de su expresión vectorial. Así pues, el signo de la función indica si el fluido se mueve en dirección del vector [math]\vec{e}_{\theta}[/math] si es positivo, o, por el contrario, en sentido opuesto si es negativo. Así, no hay que entender el signo como una longitud literal (puesto que no existen longitudes negativas), sino como información relativa al sentido de circulación del fluido respecto a la dirección [math]\vec{e}_{\theta}[/math].

Con ello, el módulo, entendido como la longitud del vector, es máximo en dos ocasiones: en [math]\rho=1[/math] y en [math]\rho=2[/math] (si tomamos el valor absoluto), siendo ésta igual a uno.

5.3 Gráfica del comportamiento del módulo en función de [math]\rho[/math]

A continuación, se presenta el dibujo de las longitudes del vector, tomando solamente los valores absolutos de la función [math]u(\rho)[/math], que queda entonces como:

[math] |u(\rho)|=\left|-\rho+\frac{2}{\rho}\right|. [/math]

7 Representación del Campo de Temperaturas

Sea la temperatura del fluido dada por la expresión [math] T(\rho,\theta) = log(1+ \rho^2)\cos^2 \theta [/math] se dibuja el campo de temperaturas y las curvas de nivel.

7.1 Representación del campo y curvas de nivel

A partir del siguiente código de MATLAB se obtiene una representación gráfica del campo escalar que describe la temperatura del fluido. Se muestran dos vistas, una tridimensional y otra bidimensional, ambas con una escala de colores que indica el valor de la temperatura en cada punto de la región considerada.

Las figuras obtenidas muestran el campo de temperatura del fluido en la sección transversal comprendida entre los dos tubos concéntricos de radio máximo \(\rho = 2\), representado en coordenadas cartesianas mediante un mapa de colores. Las zonas más claras dentro de la escala empleada corresponden a temperaturas más altas, mientras que las regiones más oscuras indican valores más bajos de la función \(T(\rho,\theta) = \log(1+\rho^{2})\cos^{2}\theta\). Las curvas de nivel recogen líneas de igual temperatura y permiten identificar con claridad el punto de máximo absoluto, marcado en rojo, situado en \(\rho = 2\), \(\theta = 0\), es decir, sobre el cilindro exterior en la dirección horizontal positiva, donde se combinan el mayor valor radial de \(\log(1+\rho^{2})\) con el máximo de \(\cos^{2}\theta\).

% Mallado
rho   = 1:0.08:2;
theta = 0:0.08:2*pi;
[RHO,THETA] = meshgrid(rho,theta);
x = RHO .* cos(THETA);
y = RHO .* sin(THETA);

% Tu campo de temperatura (sustituye por el correcto si es otro)
T = log(1 + 2*RHO) .* (cos(THETA).^2);

figure;

%% 1) Temperatura en 3D
subplot(1,3,1)
surf(x,y,T);
shading faceted           % intermedio
colormap(gca,'winter');
grid on
xlabel('Eje X1'); ylabel('Eje X2'); zlabel('Temperatura; Eje X3');
title('Representación de la temperatura en 3D');
axis tight
colorbar;

%% 2) Temperatura en 2D
subplot(1,3,2)
surf(x,y,T);
shading faceted
view(2)
colormap(gca,'winter');
grid on
xlabel('Eje X1'); ylabel('Eje X2');
title('Representación de la temperatura en 2D');
axis equal tight
colorbar;

%% 3) Líneas de nivel de la temperatura
subplot(1,3,3)
contour(x,y,T,20,'LineWidth',1.0);   % solo curvas de nivel
colormap(gca,'winter');
grid on
xlabel('Eje X1'); ylabel('Eje X2');
title('Representación de las líneas de nivel de la temperatura');
axis equal tight
colorbar;

% + PTO. MÁX
hold on;
plot(x_max, y_max, 'ro', 'MarkerFaceColor','r', 'MarkerSize',8);
hold off;

8 Gradiente de la temperatura

8.1 Cálculo del gradiente de T

El gradiente de la temperatura se define mediante la siguiente fórmula:

[math] \nabla T(\rho,\theta) = \frac{\partial T}{\partial\rho}\vec{e}_\rho + \frac{1}{ρ}\frac{\partial T}{\partial\theta}\vec{e}_\theta [/math]

De manera que:

[math] \frac{\partial T}{\partial \rho} = \frac{2\rho}{1 + \rho^2} \cos^2 \theta [/math]

[math] \frac{1}{\rho} \frac{\partial T}{\partial \theta} = -\frac{2 \sin \theta \cos \theta}{\rho} \log(1 + \rho^2) [/math]

[math] \nabla T(\rho, \theta) = \cos^{2}\theta \cdot \frac{2\rho}{1 + \rho^{2}} \vec{e}_\rho - \frac{2\sin\theta\cos\theta}{\rho} \log(1 + \rho^{2}) \vec{e}_\theta [/math]

8.2 Demostración gráfica de la ortogonalidad del gradiente de la temperatura y las curvas de nivel

El gradiente de una función escalar indica, en cada punto, la dirección en la que dicha función crece más rápidamente y su módulo mide la rapidez de ese incremento. Las líneas de nivel de la temperatura son curvas sobre las que el valor de la función permanece constante, de modo que al desplazarse a lo largo de ellas no hay variación de temperatura. Esto implica que el movimiento tangencial a una línea de nivel es siempre perpendicular a la dirección de máximo aumento, es decir, el vector gradiente es ortogonal a las curvas de nivel en todos los puntos del dominio.

Esto se visualiza gráficamente mediante el siguiente código de MATLAB.

clear
clc
% Definir el rango de rho y theta
rho = linspace(0, 2, 300);
theta = linspace(0, 2*pi, 300);

% Crear malla para rho y theta
[R, T] = meshgrid(rho, theta);

% Calcular temperatura T
Temperature = log(1 + R.^2) .* cos(T).^2;

% Calcular el gradiente analítico en coordenadas polares
dT_drho = (2 .* R) ./ (1 + R.^2) .* cos(T).^2;
dT_dtheta = -2 .* log(1 + R.^2) .* cos(T) .* sin(T);

% Evitar división por cero en rho=0
epsilon = 1e-10;
R_mod = R + epsilon;

% Convertir gradiente a coordenadas cartesianas
grad_x = cos(T) .* dT_drho - (1 ./ R_mod) .* sin(T) .* dT_dtheta;
grad_y = sin(T) .* dT_drho + (1 ./ R_mod) .* cos(T) .* dT_dtheta;

% Convertir coordenadas polares a cartesianas para graficar
X = R .* cos(T);
Y = R .* sin(T);

% Crear figura con dos subplots
figure;

% Primer subplot: curvas de nivel + gradiente
subplot(1,2,1);
contour(X, Y, Temperature, 30, 'LineWidth', 1.5);
colormap('winter');
colorbar;
hold on;
quiver(X(1:10:end, 1:10:end), Y(1:10:end, 1:10:end), grad_x(1:10:end, 1:10:end), grad_y(1:10:end, 1:10:end), 'Color', [1 0.5 0], 'LineWidth', 1);
title('Curvas de nivel y gradiente');
axis equal;
grid on;
hold off;

% Segundo subplot: solo gradiente
subplot(1,2,2);
quiver(X(1:10:end, 1:10:end), Y(1:10:end, 1:10:end), grad_x(1:10:end, 1:10:end), grad_y(1:10:end, 1:10:end), 'Color', [1 0.5 0], 'LineWidth', 1);
title('Gradiente de temperatura');
axis equal;
grid on;
axis equal;
hold off;

10 Bibliografía

A low-cost, open-source cylindrical couette rheometer. (2024). En Scientific Reports (N.o 30187). https://www.nature.com/articles/s41598-024-76494-8

Couette flow. (2018). En Science Direct. https://www.sciencedirect.com/topics/physics-and-astronomy/couette-flow

Wikipedia contributors. (2025, 15 noviembre). Taylor–Couette flow. Wikipedia. https://en.wikipedia.org/wiki/Taylor%E2%80%93Couette_flow

Apuntes proporcionados en ETSI Caminos, Canales y Puertos, de la Universidad Politécnica de Madrid