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;

2 Ecuación Navier-Stokes

2.1 Velocidad de las partículas del fluido y análisis de la ecuación

La velocidad de las partículas del fluido viene expresada según el campo vectorial [math]\vec{u}(ρ, θ) = f(ρ) \vec{e_θ} [/math], donde la presión ([math]p[/math]) es constante. Además, el campo de velocidades tiene que cumplir la ecuación de Navier-Stockes estacionaria definida por la siguiente expresión:

[math](\vec{u} · ∇)\vec{u} + ∇p = µ∆\vec{u} [/math]

Nota: El valor µ representa el coeficiente de viscosidad del fluido.

Analizando la expresión, aseguramos que al tratarse de una presión constante su gradiente es nulo ( [math]∇p = 0[/math] ) y además despreciando el término convectivo (primer término de la fórmula) finalmente obtenemos que:

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

Nota: [math]∆\vec{u} [/math] representa el laplaciano vectorial del campo de velocidades.

La simplificación obtenida de la ecuación de Navier-Stokes representa un flujo viscoso estacionario en el cual no hay gradientes de presión ni fuerzas externas con efectos inerciales significativos. Su resolución depende unicamente de las condiciones de frontera en el recinto que se indique.

2.2 Laplaciano del campo vectorial [math]\vec{u} [/math]

Por definición tenemos que el laplaciano de un campo vectorial viene dado por la siguiente ecuación:

[math]∆\vec{u} = ∆(u_1\vec{i} + u_2\vec{j}+ u_3\vec{k}) = ∆u_1\vec{i} + ∆u_2\vec{j}+ ∆u_3\vec{k}[/math]

Debido a que el campo de velocidades está proporcionado en coordenadas cilíndricas, resulta más sencillo calcularlo con dichas coordenadas y por lo tanto la ecuación pasaría a ser la siguiente:

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

2.2.1 Gradiente de la divergencia

Como podemos apreciar, debemos calcular el gradiente de la divergencia. De modo que primeramente realizaremos el cálculo de la divergencia.

Por otro lado, al tratarse de un fluido incomprensible, dicha divergencia debe ser nula.

[math] ∇ · \vec{u} = \frac{1}{ρ}[\frac{ \partial}{\partial ρ}(0) + \frac{ \partial}{\partial θ}(f(ρ)) + \frac{ \partial}{\partial z}(0) = 0 [/math]

Tras realizar los cálculos, observamos que como hemos deducido, la divergencia es nula y por lo tanto:

[math] ∇(∇ · \vec{u})= ∇(0)= \vec{0} [/math]

Finalmente tras los cálculos realizados, aseguramos que el primer término del laplaciano es nulo, simplificándose así la ecuación y quedando de la siguiente forma:

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

2.2.2 Rotacional del campo de velocidades

Realizamos los cálculos aplicando la fórmula para el rotacional definido en coordenadas cilíndricas.

[math](\nabla\times\vec{u}) = \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & ρf(ρ) & 0 \end{vmatrix} = \frac{1}{ρ}\frac{ \partial (ρf(ρ)) }{\partial ρ}\vec{e_z} = [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}]\vec{e_z}[/math]

2.2.3 Rotacional del rotacional del campo de velocidades

Una vez realizado el rotacional del campo de velocidades, proseguimos con el cálculo del rotacional sobre el vector hallado anteriormente:

[math]\nabla\times(\nabla\times\vec{u})= \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & 0 & [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}] \end{vmatrix} = -\frac{1}{ρ}\frac{ \partial }{\partial ρ}(\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}) ρ\vec{e_θ} = \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math]

2.2.4 Resultado del Laplaciano vectorial

Sustituimos en la ecuación obtenida del laplaciano vectorial.

[math]\Delta\vec{u} = - \nabla \times (\nabla \times \vec{u})[/math] [math]\Delta\vec{u} = - \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math] [math]\Delta\vec{u} = \frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math]

2.3 Resultado de la ecuación de Navier-Stokes

Una vez calculado el laplaciano, procedemos a sustituir el valor en la ecuación analizada previamente de Navier-Stokes obteniendo:

[math] µ∆\vec{u} = \vec{0}[/math] [math]∆\vec{u} = \vec{0}[/math] [math]\frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math] [math][-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math]

2.3.1 Ecuación diferencial

Analizando y simplificando la ecuación diferencial obtenida, llegamos a la siguiente expresión:

[math] [-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = [-\frac{f(ρ)}{ρ}+\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math]

Podemos ver que la simplificación está bien realizada debido a que al operar la derivada propuesta, obtendríamos la expresión de partida.

[math]\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = \frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})[/math]

Reordenamos la ecuación para la obtención de una ecuación diferencial de segundo orden:

[math]\frac{f(ρ)}{ρ}-\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = 0 [/math] [math]\frac{f(ρ)}{ρ} = \frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) [/math]

2.4 Análisis de la solución conocida

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

2.5 Obtención de [math]a[/math], [math]b[/math] y [math]\vec{u}[/math]

Para determinar [math]a[/math] y [math]b[/math], de manera que la velocidad lineal 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] , por lo que la velocidad lineal de los puntos del cilindro es [math]\rho.\omega_e=2\omega_e[/math]. El sentido de giro es contrario al creciente del parámetro [math]\theta[/math], 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] , por lo que la velocidad lineal de los puntos del cilindro es [math]\rho.\omega_i=1\omega_i[/math] . El sentido de giro es el creciente del parámetro [math]\theta[/math], 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 [math]a[/math], [math]b[/math], [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} = -2\omega_e[/math]

Una vez obtenidas ambas relaciones, se resuelve el sistema de manera trivial buscando dejar [math]a[/math] y [math]b[/math] 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} = -2\omega_e [/math]

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

[math]b = -\tfrac{2}{3}\,(-2\omega_e - 2\omega_i) \to b = \tfrac{4}{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{4}{3}\,\omega_e - \tfrac{1}{3}\,\omega_i [/math]

Tras hallar ambas constantes [math]a[/math] y [math]b[/math], la función [math]\vec{u}(\rho,\theta)[/math] 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{4}{3} w_e - \frac{1}{3} w_i \right)\rho + \left(\frac{4}{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].

2.6 Condición de incompresibilidad

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.

3 Campo de velocidades

clear; clc;
% --------------------------
% PARÁMETROS
% --------------------------
Ri = 1;       % radio interior
Re = 2;       % radio exterior
omega_i = +1;   % interior: antihorario
omega_e = -1;   % exterior: horario
% --------------------------
% COEFICIENTES A y B 
% --------------------------
B = (omega_i - omega_e) * (Ri^2 * Re^2) / (Re^2 - Ri^2);
A = omega_i - B / (Ri^2);
% --------------------------
% MALLA PARA FLECHAS 
% --------------------------
Ntheta = 20;     
Nr     = 15;     
thetas = linspace(0, 2*pi, Ntheta+1);
thetas(end) = [];    
rhos   = linspace(Ri, Re, Nr);
% --------------------------
% FIGURA
% --------------------------
figure('Color',[1 1 1]); hold on;
scale_quiver = 0.7;   
for th = thetas
    
    % velocidad tangencial u_theta en cada rho
    u_theta = A .* rhos + B ./ rhos;
    
    % componentes cartesianas
    U = -u_theta .* sin(th);
    V =  u_theta .* cos(th);
    % posiciones
    X = rhos .* cos(th);
    Y = rhos .* sin(th);
    %Vectores
    quiver(X, Y, U, V, scale_quiver, 'LineWidth', 1, ...
        'MaxHeadSize', 1.6, 'Color',[0.85 0.45 0]);
end
% --------------------------
% CILINDROS
% --------------------------
tc = linspace(0,2*pi,500);
plot(Ri*cos(tc), Ri*sin(tc), 'k', 'LineWidth', 2);
plot(Re*cos(tc), Re*sin(tc), 'k', 'LineWidth', 2);
axis equal;
xlim([-2.3 2.3]); ylim([-2.3 2.3]);
xlabel('x'); ylabel('y');
title('Campo de velocidades');
grid on;

4 Líneas de Corriente del campo

4.1 Campo [math]\vec{v}[/math] ostogonal a [math]\vec{u}[/math]

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{4}{3} \omega_e + \frac{1}{3} \omega_i \right) \rho - \left(\frac{4}{3} \omega_e + \frac{4}{3} \omega_i \right) \frac{1}{\rho} \right] \vec{e}_\rho[/math]

4.2 [math]\vec{v}[/math] irrotacional

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

4.3 Potencial escalar [math]\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{4}{3}w_{e}+\frac{1}{3}w_{i}\right)\rho - \left(\frac{4}{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{4}{3}w_{e}+\frac{1}{3}w_{i} \right)\rho - \left( \frac{4}{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{4}{3} w_e + \frac{1}{3} w_i \right) \rho \, d\rho - \int \left( \frac{4}{3} w_e + \frac{4}{3} w_i \right) \frac{1}{\rho} \, d\rho = \left( \frac{4}{3} w_e + \frac{1}{3} w_i \right) \int \rho \, d\rho - \left( \frac{4}{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{4}{3} w_e + \frac{1}{3} w_i \right) \frac{\rho^2}{2} - \left( \frac{4}{3} w_e + \frac{4}{3} w_i \right) \ln \rho + \alpha, \quad \alpha \in \mathbb{R} [/math]

4.4 Líneas de corriente del campo [math]\vec{u}[/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(-\frac{4}{3} w_e - \frac{1}{3} w_i \right)\rho + \left(\frac{4}{3} w_e + \frac{4}{3} w_i \right)\frac{1}{\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, se transforma la expresión vectorial [math]\vec{u}(\rho,\theta)[/math] en una función escalar: [math] u(\rho) = \left| \left[ \left(-\frac{4}{3} w_e - \frac{1}{3} w_i \right)\rho + \left(\frac{4}{3} w_e + \frac{4}{3} w_i \right)\frac{1}{\rho} \right] \vec{e}_\theta \right| [/math] de [math]\mathbb{R}[/math] en [math]\mathbb{R}[/math]. Si se continúa con el supuesto de que [math]|\omega_i|=|\omega_e|=1[/math], queda finalmente como:

[math]u(\rho)=-\tfrac{5}{3}\rho+\tfrac{8}{3\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)=-\tfrac{5}{3}\rho+\tfrac{8}{3\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) = -\frac{5}{3} - \frac{8}{3\rho^{2}}. [/math]

Para [math]\rho\gt0[/math], el término [math]\tfrac{8}{3\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 -\frac{5}{3}-\frac{8}{3\rho^{2}}=0 \Longrightarrow -5-\frac{8}{\rho^{2}}=0 \Longrightarrow \frac{8}{\rho^{2}}=-5 \Longrightarrow \rho^{2}=-\frac{8}{5}, [/math]

que no tiene solución real.

• Valores en los extremos del intervalo:

[math] u(1)=-\frac{5}{3}\cdot 1+\frac{8}{3\cdot 1}=\frac{-5+8}{3}=1, \qquad u(2)=-\frac{5}{3}\cdot 2+\frac{8}{3\cdot 2}=-\frac{10}{3}+\frac{4}{3}=-\frac{6}{3}=-2. [/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)=-2[/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 -\frac{5}{3}\rho+\frac{8}{3\rho}=0 \Longrightarrow -5\rho^{2}+8=0 \Longrightarrow 5\rho^{2}=8 \Longrightarrow \rho^{2}=\frac{8}{5} \Longrightarrow \rho=\sqrt{\tfrac{8}{5}}\approx 1{,}2649. [/math]

(Observación: [math]u(1)=1\gt0[/math] y [math]u(2)=-2\lt0[/math], por lo que efectivamente existe tal [math]\rho_0\in(1,2)[/math] y su valor aproximado es [math]\rho_0\approx1{,}2649.[/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 [math]\rho=2[/math] (si tomamos el valor absoluto), siendo ésta igual a dos.

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|-\tfrac{5}{3}\rho+\tfrac{8}{3\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(-(5/3)*rho + 8./(3*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, se utiliza 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 [math] \vec{u} [/math] :

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

Se puede determinar que depende únicamente de la componente [math] e_{\theta} [/math], 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, se obtiene:

[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] \vec e_{\theta} [/math] no depende de z, 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]

Se sustituye 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]

Se deriva 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

Véase una representación gráfica realizada con MATLAB: (Supóngase que tanto ωi como ωe son unitarios)

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

Rotacional = cte

%% Parámetros del problema (modulo = 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 = 2;   % radio exterior

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

%% 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;
colormap winter;

title('｜∇×u｜ con ｜ω_i｜=｜ω_e｜=1');
axis equal;

%% (Opcional) 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

% Mallado cilíndrico
rho   = 1:0.08:2;
theta = 0:0.08:2*pi;
[RHO,THETA] = meshgrid(rho,theta);

% Campo de temperatura T(rho,theta)
T = log(1 + 2*RHO) .* (cos(THETA).^2);

% Derivadas en coordenadas cilíndricas
Tr     = (2 ./ (1 + 2*RHO)) .* (cos(THETA).^2);   % dT/drho
Ttheta = -log(1 + 2*RHO) .* sin(2*THETA);         % dT/dtheta

% Componentes del gradiente (e_rho, e_theta)
Grad_rho   = Tr;
Grad_theta = Ttheta ./ RHO;

% Paso a cartesianas
X = RHO .* cos(THETA);
Y = RHO .* sin(THETA);

Ux = Grad_rho .* cos(THETA) - Grad_theta .* sin(THETA);
Uy = Grad_rho .* sin(THETA) + Grad_theta .* cos(THETA);

step  = 2;      % densidad de flechas
scale = 0.45;   % un poco más grandes

figure;

%% Subplot 1: curvas de nivel + gradiente
subplot(1,2,1);
contour(X, Y, T, 20, 'LineWidth', 1.2);
colormap('winter');
colorbar;
axis equal tight;
grid on;
xlabel('Eje X1');
ylabel('Eje X2');
title('Curvas de nivel de T y gradiente \nabla T');
hold on;
quiver(X(1:step:end,1:step:end), ...
       Y(1:step:end,1:step:end), ...
       Ux(1:step:end,1:step:end), ...
       Uy(1:step:end,1:step:end), ...
       scale, 'Color', [1 0.5 0], 'LineWidth', 1.1);
hold off;

%% Subplot 2: solo gradiente de temperatura
subplot(1,2,2);
quiver(X(1:step:end,1:step:end), ...
       Y(1:step:end,1:step:end), ...
       Ux(1:step:end,1:step:end), ...
       Uy(1:step:end,1:step:end), ...
       scale, 'Color', [1 0.5 0], 'LineWidth', 1.1);
axis equal tight;
grid on;
xlabel('Eje X1');
ylabel('Eje X2');
title('Campo del gradiente de temperatura \nabla T');

9 Caudal circulante en la sección longitudinal

Finalmente, el cálculo del caudal cuando [math]\theta = 0[/math] se abordará integrando el módulo de la velocidad antes calculado y analizado, de tal manera que, como la sección resultante es un cuadrado de área 1 m², y se supone que la velocidad viene expresada en m/s, dicho resultado es directamente el caudal [math]Q[/math], puesto [math]Q = v\cdot S[/math].

[math]\displaystyle \int_{1}^{2} u(\rho)\,d\rho = \int_{1}^{2}\left(-\tfrac{5}{3}\rho+\tfrac{8}{3\rho}\right)d\rho = \left[-\tfrac{5}{6}\rho^{2}+\tfrac{8}{3}\ln\rho\right]_{1}^{2} = \left(-\tfrac{5}{6}\cdot 2^{2}+\tfrac{8}{3}\ln 2\right) - \left(-\tfrac{5}{6}\cdot 1^{2}+\tfrac{8}{3}\ln 1\right) = -\tfrac{5}{2}+\tfrac{8}{3}\ln 2. [/math]

Cuya aproximación numérica es:

[math]\displaystyle \int_{1}^{2}\left(-\tfrac{5}{3}\rho+\tfrac{8}{3\rho}\right)d\rho = -\tfrac{5}{2}+\tfrac{8}{3}\ln 2 \approx -0.65176 \, m^{3}/s [/math]

Así pues, atendiendo a los signos, el caudal compensado, entre la velocidad que avanza según [math]\vec{e}_{\theta}[/math] y la que va en sentido contrario, queda en sentido -[math]\vec{e}_{\theta}[/math]. Dicho resultado es coherente con el análisis sobre el módulo de la velocidad, ya que se aprecia que el área bajo la curva del módulo de [math]u[/math] es mayor cuando [math]\rho \lt \sqrt{\tfrac{8}{5}}[/math]. Queda a continuación el dibujo que representa la sección y la velocidad del fluido a través de éste.

% Flujo de Couette entre cilindros - sección theta = 0
clear; close all; clc;
%% Parámetros geométricos y físicos
Ri = 1;        % radio interior
Re = 2;        % radio exterior
h  = 1;        % altura (profundidad)

% Velocidades angulares (módulo 1). 
% Ejemplo por defecto: interior -> horario (-1), exterior -> antihorario (+1)
wi = +1;   % omega interior (rad/s)
we = -1;   % omega exterior (rad/s)

%% Construcción de la solución analítica: v_theta(r) = A*r + B/r
% Condiciones: v_theta(Ri) = wi*Ri,  v_theta(Re) = we*Re

% Resolvemos para A y B
A = [Ri, 1/Ri; Re, 1/Re] \ [wi*Ri; we*Re];  % sistema 2x2
A = A(1); B = A; % overwrite fix below

% Above line was a mistake: recompute properly
M = [Ri, 1/Ri; Re, 1/Re];
coeff = M \ [wi*Ri; we*Re];
A = coeff(1);
B = coeff(2);

%% Malla en (r,z) para representar la sección theta = 0
Nr = 10; Nz = 10;
r = linspace(Ri,Re,Nr);
z = linspace(0,h,Nz);
[R_grid, Z_grid] = meshgrid(r,z);

% Campo tangencial v_theta
v_theta = A .* R_grid + B ./ R_grid;   % magnitud en dirección e_theta

% Convertir a componentes cartesianas (en theta = 0 -> e_theta = [0,1,0])

Ux = zeros(size(v_theta)); 
Uy = v_theta; 
Uz = zeros(size(v_theta));

%% Figura 3D
figure('Color','w','Position',[200 100 900 700]);
hold on; axis equal; 
view(35,18);

% Dibujar cilindros (superficies)
nCirc = 150;
[Xi,Yi,Zi] = cylinder(Ri,nCirc); Zi = Zi*h;
surf(Xi,Yi,Zi, 'FaceAlpha',0.25, 'EdgeColor','none', 'FaceColor',[0.85 0.8 0.6]);
[Xe,Ye,Ze] = cylinder(Re,nCirc); Ze = Ze*h;
surf(Xe,Ye,Ze, 'FaceAlpha',0.25, 'EdgeColor','none', 'FaceColor',[0.85 0.9 0.75]);

% Dibujar la sección theta = 0: es el plano y = 0 (mostramos solo la porción r∈[Ri,Re], z∈[0,h], x>=0)
% En coordenadas cartesianas, sobre theta=0 los puntos tienen (x = r, y = 0).
X_patch = [Ri Re Re Ri];  % x extents
Y_patch = [0 0 0 0];      % y = 0 (plano)
Z_patch = [0 0 h h];
patch(X_patch, Y_patch, Z_patch, 'r','FaceAlpha',0.12,'EdgeColor','r','LineWidth',1.8);

% Resaltar aristas de la sección en rojo
plot3([Ri Re],[0 0],[0 0],'r','LineWidth',2)   % borde inferior de la sección
plot3([Ri Re],[0 0],[h h],'r','LineWidth',2)   % borde superior
plot3([Ri Ri],[0 0],[0 h],'r','LineWidth',2)   % arista interior
plot3([Re Re],[0 0],[0 h],'r','LineWidth',2)   % arista exterior

% Dibujar campo de velocidades en el plano theta=0 (posiciones x=r, y=0)
[xq,zq] = meshgrid(r,z);
yq = zeros(size(xq));

% Aquí Ux = 0, Uy = v_theta
quiverScale = 2;
quiver3(xq, yq, zq, Ux, Uy, Uz, quiverScale, 'b', 'LineWidth', 1.5, 'MaxHeadSize',0.5);
title('Sección \theta = 0 del flujo de Couette entre cilindros');
xlim([-0.1, Re+0.5]);
ylim([-Re-0.5, Re+0.5]); 
zlim([0 h]);
axis off;

% Mejoras estéticas
camlight; lighting phong;
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