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]);                           % Centrar y ampliar

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

3 Campo de velocidades

% Parámetros
Ri = 1;               % Radio interior
Ro = 2;               % Radio exterior
Nr = 20;              % nodos radiales
Ntheta = 40;          % nodos angulares
omega1 = 1;           % velocidad angular cilindro 
                        interior
omega2 = -1;          % velocidad angular cilindro 
                        exterior
% Mallado polar
r = linspace(Ri, Ro, Nr);
theta = linspace(0, 2*pi, Ntheta);
[R, Theta] = meshgrid(r, theta);
X = R .* cos(Theta);
Y = R .* sin(Theta);

% Campo de velocidades (Couette)
A = (omega2*Ro - omega1*Ri) / (Ro - Ri);
B = (omega1*Ri*Ro - omega2*Ro^2) * Ri / (Ro - Ri);
Vtheta = A * R + B ./ R;
Vx = -Vtheta .* sin(Theta);
Vy = Vtheta .* cos(Theta);
quiver(X, Y, Vx, Vy, 0.6, 'Color', [1 0.5 0], 'LineWidth', 1.3); % naranja

% Encadrar y ampliar espacio en los ejes
padding = 0.25 * (Ro - Ri);
xlim([-Ro-padding, Ro+padding]);
ylim([-Ro-padding, Ro+padding]);
axis equal;

title('CAMPO DE VELOCIDADES');
xlabel('x');
ylabel('y');
hold off;

Para comprobar que la función [math] f(\rho) = a \rho + \frac{b}{\rho}, \quad a,b \in \mathbb{R} [/math], es solución de la ecuación diferencial [math] \frac{\partial}{\partial \rho} \Big( \rho \frac{\partial f(\rho)}{\partial \rho} \Big) = \frac{f(\rho)}{\rho} [/math], se calculan las expresiones de ambos miembros con el fin de llegar a una expresión común:

El miembro de la izquierda es:

[math] \frac{f(\rho)}{\rho} = a + \frac{b}{\rho^2} [/math]

Por otro lado, el miembro de la derecha es el siguiente:

[math] \frac{\partial}{\partial \rho} \Big( \rho(a - \frac{b}{\rho^2}) \Big) = \frac{\partial}{\partial \rho} (a \rho - \frac{b}{\rho}) = a + \frac{b}{\rho^2} [/math]

Como puede comprobarse, en ambos miembros se obtiene la misma expresión, por lo que se puede afirmar que [math] f(\rho) = a \rho + \frac{b}{\rho}, \quad a, b \in \mathbb{R} [/math] es solución de la ecuación diferencial [math] \frac{\partial}{\partial \rho} \Big( \rho \frac{\partial f(\rho)}{\partial \rho} \Big) = \frac{f(\rho)}{\rho}. [/math]

Para determinar a y b, de manera que la velocidad del fluido en la frontera coincida con la de los cilindros interior y exterior es necesario introducir una serie de condiciones sobre los parámetros:

El cilindro exterior queda proyectado en [math]z=0[/math] sobre [math]\rho=2[/math] y gira en sentido horario con una velocidad angular constante [math]\omega_e[/math] , el sentido de giro es contrario al creciente del parámetro θ, por lo que en la velocidad se considerará [math]-\vec{e}_\theta[/math].

En cambio, el cilindro interior queda proyectado en [math]z=0[/math] sobre [math]\rho=1[/math] y gira en sentido antihorario con una velocidad angular constante [math]\omega_i[/math] , el sentido de giro es el creciente del parámetro θ, por lo que en la velocidad se considerará [math]\vec{e}_\theta[/math].

Una vez impuestas las condiciones sobre los parámetros, se obtiene, para cada valor de [math]\rho[/math] una relación entre a, b, [math]\omega_i[/math] y [math]\omega_e[/math]:

[math]\rho = 1 \to \vec u = \omega_i \,\vec e_\theta = f(\rho)\,\vec e_\theta = f(1)\,\vec e_\theta = (a\cdot 1 + \tfrac{b}{1})\,\vec e_\theta = \omega_i \,\vec e_\theta \;\;\Rightarrow\;\; a + b = \omega_i[/math] [math]\rho = 2 \to \vec u = -\omega_e \,\vec e_\theta = -f(\rho)\,\vec e_\theta = -f(2)\,\vec e_\theta = -(2a + \tfrac{b}{2})\,\vec e_\theta = -\omega_e\,\vec e_\theta \;\;\Rightarrow\;\; 2a + \frac{b}{2} = -\omega_e[/math]

Una vez obtenidas ambas relaciones, se resuelve el sistema de manera trivial buscando dejar a y b en función de [math]\omega_i[/math] y [math]\omega_e[/math]. De la primera relación de despeja [math]a[/math] de la siguiente manera:

[math] a = \omega_i - b [/math]

A continuación, se introduce en la segunda relación y se despeja [math]b[/math]:

[math]2(\omega_i - b) + \frac{b}{2} = -\omega_e [/math]

[math]b(-2 + \tfrac{1}{2}) = -\omega_e - 2\omega_i [/math]

[math]b = -\tfrac{2}{3}\,(-\omega_e - 2\omega_i) \to b = \tfrac{2}{3}\,\omega_e + \tfrac{4}{3}\,\omega_i [/math]

Una vez obtenida [math]b[/math], se sustituye en la primera relación para obtener [math]a[/math]:

[math]a = \omega_i - \left( \tfrac{2}{3}\,\omega_e + \tfrac{4}{3}\,\omega_i \right) \to a = -\tfrac{2}{3}\,\omega_e - \tfrac{1}{3}\,\omega_i [/math]

Tras hallar ambas constantes [math]a[/math] y [math]b[/math], la función \vec{u} queda de la siguiente manera:

[math]\vec{u}(\rho,\theta) = f(\rho)\, \vec{e}_\theta = \left(a \rho + \frac{b}{\rho}\right) \vec{e}_\theta[/math] [math]\vec{u}(\rho,\theta) = \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] \vec{e}_\theta[/math]

Es importante comentar que el campo solo depende de [math]\rho[/math] y las velocidades [math]\omega_i[/math] y [math]\omega_e[/math]; y no de [math]\theta[/math] o [math]z[/math]. Este suceso era de esperar ya que el flujo es constante en cualquier altura de los cilindros y en cualquier ángulo [math]\theta[/math].

También es interesante comprobar la condición de incompresibilidad o divergente nulo; para ello es necesario introducir la forma en la que se calcula la divergencia de un campo vectorial en la base cilíndrica: [math]\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)[/math].

A continuación, se calcula la divergencia del campo de velocidades del fluido entre los cilindros coaxiales:

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

La divergencia es nula porque la componente [math]u_\theta[/math], que en este caso es el campo en su totalidad, no depende de [math]\theta[/math] y por lo tanto, la derivada parcial / total es nula.

4 Líneas de Corriente del campo

Para trazar las líneas de corriente del campo [math]\vec{u}[/math] es necesario calcular y dibujar las líneas [math]\psi = \text{cte}[/math] del potencial escalar [math]\psi[/math] de un campo ortogonal, [math]\vec{v}[/math], a [math]\vec{u}[/math]; que se conocen como función de corriente de [math]\vec{u}[/math].

Para calcular un campo [math]\vec{v}[/math] ortogonal a [math]\vec{u}[/math] es necesario conocer cómo es el campo [math]\vec{u}[/math]; éste es un campo que, en la base cilíndrica, no tiene ninguna componente radial, solo tangencial a la velocidad (en cada punto). Por lo que cualquier campo que solo tenga componentes radiales / perpendiculares a la velocidad en cada punto ya es ortogonal a [math]\vec{u}[/math]. En este caso, para simplificar los cálculos, como el campo de velocidades (en la base cartesiana con el eje z coincidente con el eje de los cilindros) solo tiene componentes [math]\vec{i}[/math] y [math]\vec{j}[/math], el campo [math]\vec{v}[/math] puede considerarse como: [math]\vec{v} = \vec{k} \times \vec{u}[/math].

Por lo tanto, el campo [math]\vec{v}[/math] se obtiene de la siguiente manera:

[math]\vec{v} = \begin{vmatrix} \vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\ 0 & 0 & 1 \\ 0 & u_\theta & 0 \end{vmatrix} = -u_\theta \, \vec{e}_\rho \to \vec{v} = \left[ \left(\frac{2}{3} \omega_e + \frac{1}{3} \omega_i \right) \rho - \left(\frac{2}{3} \omega_e + \frac{4}{3} \omega_i \right) \frac{1}{\rho} \right] \vec{e}_\rho[/math]

Como ya se ha demostrado, [math]\vec{u}[/math] tiene divergencia nula; por lo que el rotacional de [math]\vec{v}[/math] también será nulo. Para comprobarlo es necesario exponer previamente la fórmula del rotacional de un campo vectorial en la base cilíndrica:

[math] \nabla \times \vec{v} = \frac{1}{\rho} \begin{vmatrix} \vec e_{\rho} & (\rho \vec e_{\theta}) & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial \theta} & \dfrac{\partial}{\partial z} \\ u_{\rho} & \rho u_{\theta} & u_{z} \end{vmatrix} [/math]

En este caso el rotacional de campo [math]\vec{v}[/math] es el siguiente:

[math] \nabla \times \vec{v} = \frac{1}{\rho} \begin{vmatrix} \vec e_{\rho} & (\rho \vec e_{\theta}) & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial \theta} & \dfrac{\partial}{\partial z} \\ u_{\rho} & 0 & 0 \end{vmatrix} = 0 [/math]

Se observa que [math]u_{\rho}[/math] solo depende de [math]\rho[/math], como cabría esperar ya que el flujo es constante [math]\forall \theta \in [0,2\pi),\ \forall z \gt 0[/math], por lo que todas las derivadas cruzadas son nulas. Esto permite afirmar que el campo [math]\vec{v}[/math] es conservativo, por lo que tiene un potencial escalar [math]\psi[/math] tal que [math]\vec{v} = \nabla \psi[/math].

El cálculo del potencial escalar [math]\psi[/math] de [math]\vec{v}[/math] debe abordarse como la igualación de componentes de dos campos. El primero de ellos es [math]\nabla \psi[/math] que se expresa, en componentes, de la siguiente manera: [math]\psi'_{\rho}\,\vec e_{\rho} + \rho\,\psi'_{\theta}\,\vec e_{\theta} + \psi'_{z}\,\vec e_{z}[/math]; mientras que el campo [math]\vec{v}[/math] se escribe: [math]v_{\rho}\,\vec e_{\rho} + v_{\theta}\,\vec e_{\theta} + v_{z}\,\vec e_{z}[/math]. De esta forma, y sabiendo que están en la misma base, se igualan las componentes:

[math] \psi'_{\rho}\,\vec e_{\rho} + \rho\,\psi'_{\theta}\,\vec e_{\theta} + \psi'_{z}\,\vec e_{z} = v_{\rho}\,\vec e_{\rho} + v_{\theta}\,\vec e_{\theta} + v_{z}\,\vec e_{z} [/math]

En este caso, el vector [math]\vec{v}[/math] solo tiene componente [math]\vec e_{\rho}[/math] por lo que la igualdad queda de la siguiente manera:

[math] \psi'_{\rho}\,\vec e_{\rho} = \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]\vec e_{\rho} [/math]

Igualando las componentes:

[math] \psi'_{\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]

Como se observa, y se ha comentado anteriormente, las componentes [math]v_{\theta}\,\vec e_{\theta} [/math] y [math]v_{z}\,\vec e_{z}[/math] son ambas nulas, por lo que el potencial escalar es [math]\psi = \psi(\rho)[/math]. Entonces se debe afirmar la siguiente igualdad:

[math] \psi = \int v_{\rho} \, d\rho [/math]

Donde la constante de integración no será [math]f(\theta, z)[/math] sino una constante escalar arbitraria [math]\alpha[/math].

De esta manera, [math]\psi[/math] se obtiene de la siguiente forma:

[math] \psi = \int \left( \frac{2}{3} w_e + \frac{1}{3} w_i \right) \rho \, d\rho - \int \left( \frac{2}{3} w_e + \frac{4}{3} w_i \right) \frac{1}{\rho} \, d\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]

Y el potencial escalar [math]\psi[/math] queda:

[math] \psi = \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 + \alpha, \quad \alpha \in \mathbb{R} [/math]

Finalmente, al representar las líneas [math]\psi=cte[/math] ; como [math]\psi = \psi(\rho)[/math], en realidad se estarían representando las líneas [math]\rho=cte[/math]. Y éstas son directamente las líneas de corriente del campo [math]\vec{u}[/math].

Para la visualización gráfica es necesario apoyarse en el siguiente código de MATLAB:

% Parámetros geométricos
Ri = 1;
Re = 2;
% Valores de rotación, ambos módulos de we y wi = 1
% Sentido horario , signo negativo
% Sentido antihorario , signo positivo
we = -1;   % cilindro exterior
wi = +1;   % cilindro interior
% Coeficientes de la función de corriente (según tu fórmula)
A = (2/3)*we + (1/3)*wi;
B = (2/3)*we + (4/3)*wi;
% Mallado polar
N = 300;
rho = linspace(0,3,N);
theta = linspace(0,2*pi,N);
[R,T] = meshgrid(rho,theta);
% Coordenadas cartesianas
X = R.*cos(T);
Y = R.*sin(T);
% Función de corriente
Psi = A*(R.^2)/2 - B*log(R);
% Enmascaramos zona que NO es fluido
mask = (R < Ri / R > Re);
Psi(mask) = NaN; 
% Dibujar líneas psi = cte
figure; hold on;
contour(X, Y, Psi, 20, 'LineWidth', 1.3);
% Dibujar los cilindros
th = linspace(0,2*pi,400);
plot(Ri*cos(th), Ri*sin(th), 'k', 'LineWidth', 2);
plot(Re*cos(th), Re*sin(th), 'k', 'LineWidth', 2);
axis equal;
xlabel('x');
ylabel('y');
title('\psi(\rho) = cte . Líneas de corriente del campo U');
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]

El código empleado para el dibujo de la gráfica se presenta a continuación.

% Intervalo de trabajo
rho = linspace(1, 2, 500);

% Definimos la función en valor absoluto
u = abs(-rho + 2./rho);

% Dibujar la función
figure;
plot(rho, u, 'LineWidth', 2);
grid on;

% Etiquetas y título
xlabel('\rho');
ylabel('｜u(\rho)｜');
title('Módulo de la velocidad');

6 Rotacional de [math] \vec{u} [/math]

6.1 Cálculo del rotacional

Para calcular el rotacional de un campo vectorial en coordenadas cilíndricas, podemos utilizar la siguiente expresión:

[math] \nabla \times \vec{u} = \frac{1}{\rho} \begin{vmatrix} \vec e_{\rho} & (\rho \vec e_{\theta}) & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial \theta} & \dfrac{\partial}{\partial z} \\ u_{\rho} & {\rho}u_{\theta} & u_{z} \end{vmatrix} [/math]

En este caso, conociendo ‘u’ :

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

Podemos determinar que depende únicamente de la componente ēθ , siendo las otras dos componentes nulas.

[math] u_{\rho}=0,\qquad u_{z}=0,\qquad u_{\theta}=\left(-\tfrac{1}{3}\omega_i-\tfrac{2}{3}\omega_e\right)\rho +\frac{\tfrac{4}{3}\omega_i+\tfrac{2}{3}\omega_e}{\rho} [/math]

Aplicando en el determinante alguna de las fórmulas para su resolución se obtendrá el resultado. Se trabajará en este caso con menores, para facilitar los cálculos sabiendo que dos de las tres componentes son nulas. Se calculará desarrollando por la fila de abajo, pasando ⍴ al otro lado, multiplicando, quedará:

[math] \rho(\nabla\times\vec{u}) = -\,\rho u_{\theta} \begin{vmatrix} \vec e_{\rho} & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial z} \end{vmatrix} [/math]

Aplicando el determinante y desarrollando por la fila inferior, obtenemos:

[math] \rho(\nabla\times\vec{u}) = -\,\rho u_{\theta} \begin{vmatrix} \vec e_{\rho} & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial z} \end{vmatrix} [/math]

El determinante 2×2 anterior es:

[math] \begin{vmatrix} \vec e_{\rho} & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial z} \end{vmatrix} = \vec e_{\rho}\,\dfrac{\partial}{\partial z} - \vec e_{z}\,\dfrac{\partial}{\partial \rho} [/math]

Por lo que, como [math]u_{\theta}[/math] no depende de z, nos queda únicamente

[math] \rho(\nabla\times\vec u)=\vec e_{z}\,\dfrac{\partial}{\partial\rho}(\rho u_{\theta}) [/math]

Volviendo a dividir entre 𝜌:

[math] \nabla\times\vec u = \vec e_{z}\, \dfrac{1}{\rho}\, \dfrac{\partial}{\partial\rho}(\rho u_{\theta}) [/math]

Sustituimos la expresión explícita de [math]u_{\theta}[/math]:

[math] \rho u_{\theta} = \left(-\tfrac{1}{3}\omega_i-\tfrac{2}{3}\omega_e\right)\rho^{2} +\left(\tfrac{4}{3}\omega_i+\tfrac{2}{3}\omega_e\right) [/math]

Derivamos con respecto a 𝜌:

[math] \dfrac{\partial}{\partial\rho}(\rho u_{\theta}) = 2\left(-\tfrac{1}{3}\omega_i-\tfrac{2}{3}\omega_e\right)\rho [/math]

Finalmente, sustituyendo:

[math] \nabla\times\vec u = \left(-\tfrac{2}{3}\omega_i-\tfrac{4}{3}\omega_e\right)\vec e_{z} [/math]

Por tanto, el rotacional es constante en todo el dominio. No existen puntos con mayor rotacional: el flujo de Couette presenta una vorticidad uniforme.

6.2 Representación gráfica del rotacional

Veamos una representación gráfica realizada con MATLAB: (Supondremos que tanto ωi como ωe son unitarios)

Donde podemos apreciar que es constante. A continuación el código empleado para su representación:

Rotacional = cte

%% Parámetros del problema (módulo = 1)
omega_i = 1;      % +1 o -1 según el sentido del cilindro interior
omega_e = -1;     % +1 o -1 según el sentido del cilindro exterior

a = 1;            % radio interior
b = 3;            % radio exterior

%% Cálculo del rotacional (constante)
A = -1/3*omega_i - 2/3*omega_e;
curl_value = abs(2*A);    

%% Mallado para representación
Nr = 200;
Ntheta = 300;

rho = linspace(a, b, Nr);
theta = linspace(0, 2*pi, Ntheta);
[RR, TT] = meshgrid(rho, theta);

% El rotacional es constante → matriz uniforme
CURL = curl_value * ones(size(RR));

%% Conversión a cartesianas
X = RR .* cos(TT);
Y = RR .* sin(TT);

%% Gráfica
figure;
pcolor(X, Y, CURL);
shading interp;
colorbar;
colormap turbo;

xlabel('x');
ylabel('y');
axis equal;

%% Dibujar los cilindros
hold on;
th = linspace(0, 2*pi, 400);
plot(a*cos(th), a*sin(th), 'k', 'LineWidth', 2);
plot(b*cos(th), b*sin(th), 'k', 'LineWidth', 2);
hold off;

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;

9 Caudal circulante en la sección longitudinal

Suponiendo que las velocidades calculadas están expresadas en m/s, y que los cilindros tienen una profundidad de 1 metro; el caudal por la sección [math]\theta=0[/math] se calculará de la siguiente mnera:

El caudal se interpreta físicamente como el flujo del campo [math]\vec{u}[/math] a través de la sección de estudio [math]\theta=0[/math]; por lo que queda, en forma genérica, de la siguiente manera: [math]Q = \int_{S} \vec{u}\cdot \vec{N}\, dS[/math]. Donde [math]\vec{u}[/math] es el campo de velocidades del fluido en cada punto, [math]\vec{N}[/math] es el vector normal a la sección (que en este caso es [math]\vec e_{\theta}[/math]) y [math]S[/math] es la sección de estudio.

Particularizando en los vectores y la sección de estudio:

[math]Q = \int_{S} \vec{u}\cdot \vec{N}\, dS= \int_{S} \vec{u}\cdot \vec{e}_{\theta}\, dS = \int_{0}^{z} \int_{\rho_{1}}^{\rho_{2}} \vec{u}\cdot \vec{e}_{\theta}\, d\rho\, dz[/math]

Sustituyendo finalmente los valores, la expresión queda de la siguiente manera:

[math]Q=\int_{S}\left[\left(-\tfrac13 w_i-\tfrac23 w_e\right)\rho+ \left(\tfrac43 w_i+\tfrac23 w_e\right)\frac{1}{\rho}\right] \, \vec e_\theta\cdot\vec e_\theta \, ds =\int_{S}\left[\left(-\tfrac13 w_i-\tfrac23 w_e\right)\rho+ \left(\tfrac43 w_i+\tfrac23 w_e\right)\frac{1}{\rho}\right] ds.[/math]

[math]Q=\int_{0}^{1}\int_{1}^{2} \left[\left(-\tfrac13 w_i-\tfrac23 w_e\right)\rho+ \left(\tfrac43 w_i+\tfrac23 w_e\right)\frac{1}{\rho}\right] \, d\rho\, dz =\int_{0}^{1}\left[\left(-\tfrac13 w_i-\tfrac23 w_e\right) \int_{1}^{2}\rho\, d\rho+ \left(\tfrac43 w_i+\tfrac23 w_e\right) \int_{1}^{2}\frac{1}{\rho}\, d\rho\right] dz.[/math]

[math]Q=\int_0^1 \left[ \tfrac32\left(-\tfrac13 w_i-\tfrac23 w_e\right) +\ln 2\left(\tfrac43 w_i+\tfrac23 w_e\right) \right] dz = \int_0^1 \left[ w_i\left(-\tfrac12+\tfrac43\ln 2\right) +w_e\left(-1+\tfrac23\ln 2\right) \right] dz.[/math]

[math]Q=\left(-\tfrac12+\tfrac43\ln 2\right) w_i +\left(-1+\tfrac23\ln 2\right) w_e.[/math]

Es importante comentar que el flujo es constante para cualquier sección definida de la manera [math]\theta=\alpha, \quad \alpha \in \mathbb{R}[/math], ya que la velocidad del fluido es constante para [math]\theta=cte[/math] y [math]\rho=cte[/math], y el vector normal a las secciones siempre va a ser [math]\vec e_{\theta}[/math]. Ergo, el flujo calculado para la sección [math]\theta=0[/math] es, en realidad, el flujo del fluido a través de cualquier sección vertical (que comprenda todo el espacio entre los cilindros y [math]\theta=cte[/math]).

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