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 Gaitán Tsukasov, Pablo Casado del Campo, Esteban Gabaldón Hermoso , Carlos Ascanio Martín, Aarón Pérez Luna
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Flujo de Couette entre dos tubos concéntricos

Introducción

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

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

1 Mallado de la sección transversal

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

Leyenda: Malla de elementos finitos para el flujo de Couette.

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

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

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

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

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

Siendo \mu el coeficiente de viscosidad del fluido.

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

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

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

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

2.2 Cálculo del gradiente de una divergencia

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

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

Si se particulariza al ejercicio se obtiene que:

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

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

2.3 Cálculo del rotacional del rotacional

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

[math]\nabla \times \vec{F}(\rho, \theta, z) = \frac{1}{\rho} \begin{vmatrix} \vec{e}_{\rho} & \rho \vec{e}_{\theta} & \vec{e}_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ F_{\rho} & \rho F_{\theta} & F_{z} \end{vmatrix}[/math]

Así, particularizando al ejercicio se obtiene:

[math]\nabla \times \vec{u} = \frac{1}{\rho} \begin{vmatrix} \vec{e}_{\rho} & \rho \vec{e}_{\theta} & \vec{e}_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ u_{\rho} & \rho u_{\theta} & u_{z} \end{vmatrix} = \left( \frac{1}{\rho} \frac{\partial}{\partial \rho}(\rho f(\rho)) \right) \vec{e}_{z}[/math]

Tras ello, se desarrolla el siguiente rotacional:

[math]\nabla \times (\nabla \times \vec{u}) = - \frac{1}{\rho} \det\begin{vmatrix}\vec{e}_{\rho} & \rho \vec{e}_{\theta} & \vec{e}_{z} \\\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\0 & \frac{1}{\rho} \frac{\partial}{\partial \rho} (\rho f(\rho)) & 0\end{vmatrix} = -\left[ \frac{\partial}{\partial \rho}\left(\frac{1}{\rho}\frac{\partial}{\partial \rho}(\rho f(\rho))\right) \right] \vec{e}_{\theta} = -\left[ \frac{\rho^2 \frac{\partial^2 f}{\partial \rho^2} + \rho \frac{\partial f}{\partial \rho} - f(\rho)}{\rho^2} \right] \vec{e}_{\theta}= -\left[ \frac{\partial^2 f}{\partial \rho^2} + \frac{1}{\rho} \frac{\partial f}{\partial \rho} - \frac{f(\rho)}{\rho^2} \right] \vec{e}_{\theta}[/math]

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

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

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

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

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

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

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

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

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

2.5 Comprobación de la ecuación obtenida

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

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

Para ello hay que derivar la expresión obtenida:

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

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

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

2.6 Calculo de los parámetros a y b

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

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

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

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

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

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

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

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

2.7 Condición de incompresibilidad

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

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

3 Representación del campo de velocidades

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

Leyenda: Representación gráfica campo de velocidades.

clear; clc;

Ri = 1;         % radio interior
Re = 2;         % radio exterior
omega_i = +1;   % interior: antihorario (convención +)
omega_e = -1;   % exterior: horario     (convención -)

B = (omega_i - omega_e) * (Ri^2 * Re^2) / (Re^2 - Ri^2);
A = omega_i - B / (Ri^2);

Ntheta = 20;     
Nr     = 15;     
thetas = linspace(0, 2*pi, Ntheta+1);
thetas(end) = [];           % evita duplicar 0 y 2*pi
rhos   = linspace(Ri, Re, Nr);

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 (tangenciales)
    U = -u_theta .* sin(th);
    V =  u_theta .* cos(th);

    % posiciones
    X = rhos .* cos(th);
    Y = rhos .* sin(th);

    % vectores del campo
    quiver(X, Y, U, V, scale_quiver, 'LineWidth', 1.1, ...
        'MaxHeadSize', 1.6, 'Color','r'); % azul verdoso
end

tc = linspace(0,2*pi,500);
plot(Ri*cos(tc), Ri*sin(tc), 'k', 'LineWidth', 2.0);
plot(Re*cos(tc), Re*sin(tc), 'k', 'LineWidth', 2.0);

axis equal;
xlim([-2.3, 2.3]); ylim([-2.3, 2.3]);

xlabel('X (m)'); 
ylabel('Y (m)');
title('Gráfica del campo de velocidades');

grid on;
holf off;

4 Líneas de Corriente

4.1 Campo v ortogonal a u

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

4.2 Irrotacionalidad de v

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

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

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

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

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

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

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

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

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

4.3 Función de corriente de u (𝜓)

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

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

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

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

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

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

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

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

El resultado sería:

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

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

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

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

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

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

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

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

Procedemos a sustituir las componentes por las aportadas anteriormente:

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

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

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

4.5 Líneas de corriente

Leyenda: Líneas de Corriente.

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

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

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

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

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

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

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

5 Puntos de velocidad máxima del fluido

Para enconntrar los puntos dentro de los cilindros donde la velocidad de nuestro fluido es máxima, debemos estudiar la primera derivada del campo de velocidades obtenido: [math]\vec{u}(\rho,\theta) = \left[ \left( -\frac{2}{3}\omega_e - \frac{1}{3}\omega_i \right)\rho + \frac{\frac{2}{3}\omega_e + \frac{4}{3}\omega_i}{\rho} \right]\vec{e}_\theta[/math]

6 Rotacional de u

6.1 Cálculo del Rotacional

Las componentes del campo [math]u[/math] son: [math]u_\rho = 0, u_\theta = f(\rho), u_z = 0[/math]

La fórmula del rotacional de un vector u es: [math]\nabla \times \mathbf{u} = \frac{1}{\rho} \left[ \frac{\partial u_z}{\partial \theta} - \frac{\partial (\rho u_\theta)}{\partial z} \right] \mathbf{e}_\rho + \left[ \frac{\partial u_\rho}{\partial z} - \frac{\partial u_z}{\partial \rho} \right] \mathbf{e}_\theta + \frac{1}{\rho} \left[ \frac{\partial (\rho u_\theta)}{\partial \rho} - \frac{\partial u_\rho}{\partial \theta} \right] \mathbf{e}_z[/math]

Al sustituir las componentes en la ecuación del rotacional de la componente [math]\vec{e}_{\rho}[/math] nos queda lo siguiente: [math]\frac{1}{\rho} \left[ \frac{\partial (0)}{\partial \theta} - \frac{\partial (\rho f(\rho))}{\partial z} \right][/math]

En este caso el resultado sería = 0 y con la componente [math]\vec{e}_{\theta}[/math] ocurre lo mismo.

Por último calculamos la componente [math]\vec{e}_{z}[/math]  : [math] (\nabla \times \mathbf{u})_{\mathbf{e}_z} = \frac{1}{\rho} \left[ \frac{\partial (\rho f(\rho))}{\partial \rho} - \frac{\partial (0)}{\partial \theta} \right] [/math]

El segundo valor es 0 así que solo debemos calcular el primero con [math] \rho f(\rho) = \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] = \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) \rho^2 + \left( \frac{2}{3} w_e + \frac{4}{3} w_i \right) [/math]

Calculamos la derivada de eso: [math] \frac{d}{d\rho}(\rho f(\rho)) = \frac{d}{d\rho} \left[ \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) \rho^2 + \text{Constante} \right] = \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) (2\rho) + 0 = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \rho [/math]

Ahora solo debemos sustituir en la componente [math]\vec{e}_{z}[/math]

[math] (\nabla \times \mathbf{u})_{\mathbf{e}_z} = \frac{1}{\rho} \left[ \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \rho \right] = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) [/math]

Por lo tanto, el rotacional de u es [math] \nabla \times \mathbf{u} = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \mathbf{e}_z [/math]

6.2 Representación del rotacional

Leyenda: Rotacional de u.

%% Parámetros y Cálculo de la Magnitud Escalar 
Ri = 1; Re = 2;
we = -1; wi = +1; 

% La magnitud del rotacional es un escalar constante:
Rot = abs((-4/3) * we - (2/3) * wi);

% Calcular el valor numérico
fprintf('La magnitud del rotacional es constante: %.4f\n', Rot);

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

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

% Debe estar en una matriz para la función pcolor
Rot = ones(size(R)) * MagRot_val;

% Enmascaramos para solo mostrar el anillo del rotacional
mask = (R < Ri | R > Re);
Rot(mask) = NaN; 

%% Representación Gráfica 
figure; 
hold on;
grid on;
% Usamos pcolor para mapear la magnitud del campo
pcolor(X, Y, Rot);
shading flat; % Suaviza la cuadrícula
colormap('hot'); % Utiliza un mapa de colores calido

% Dibujar los bordes
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); 

% Configuración del gráfico
axis equal;
axis([-L, L, -L, L]);
xlabel('x');
ylabel('y');
title('Rotacional de u');

7 Gradiente de la temperatura

7.1 =Representación grafica

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

Leyenda: Figura 8: Gradiente temperatura/curvas de nivel.

clear all % Limpiar la pantalla y espacio de trabajo
clc

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

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

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

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

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

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

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

7.2 Calculo del gradiente

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

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

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

8 Caudal de una sección longitudinal

Para comenzar calcularemos el caudal que atraviesa los dos cilindros por una sección exterior del cilindro.

Para ello tendremos en cuenta que los dos cilindros tiene 2 metros de profundidad y la velocidad del fluido esta condicionada por la ecuacion: [math]\vec{u}(\rho,\theta) = \left(\frac{4}{3}\rho - \frac{4}{3}\rho\right)\,\vec e_{\theta} [/math]

Donde depende de la variable se obtienen estos datos:

[math]\quad u_{\rho}=0,\quad u_{\theta}=\frac{4}{3}\rho - \frac{4}{3}\rho,\quad u_{z}=0 [/math]

Para calcular el caudal que recorre la superficie del cilindro(Imagen 9) se usa esta formula la cual consiste en integrar la derivada de la superficie por la velocidad:

[math] \int_{S} \vec{u}\cdot d\vec{S} \;=\; \iint_{D} \vec{u}\bigl(\Phi(u,v)\bigr)\,\cdot\, \left(\Phi_{u}\times\Phi_{v}\right)\,du\,dv [/math]

si elk caudal sale negativo significa que hemos tomado la dirreccion contraria al vector theta.

8.1 Parametrización de la superficie

Para poder calcular el caudal es necesario parametrizas la superficie por la cual circula el caudal (Imagen 9)

Esta es una de las posibles parametrizaciones de la superficie S:

[math] \vec{r}(\theta, z) = \begin{pmatrix} R \cos\theta \\ R \sin\theta \\ z \end{pmatrix}, \qquad 0 \le \theta \le 2\pi,\quad 0 \le z \le H. [/math]

calculando con esta parametrizacion nos da que el flujo es negativo por lo tanto buscamos otra que nos lo de correctamente:

8.2 Calculo del Caudal

Una vez obtenida la parametrizacion se puede obtener el vector u en funcion de esta:

Leyenda: Superficie cilindrica exterior recorrida por el caudal

R_exterior = 2.0;      % Radio exterior
altura     = 2.0;      

% Generar cilindro exterior
[Xe, Ye, Ze] = cylinder(R_exterior, 80);
Ze = Ze * altura;

% GRAFICA
figure;
hold on;
grid on;
axis equal;

color_superficie = [0.1 0.7 0.8];

% Superficie cilíndrica exterior
surf(Xe, Ye, Ze, ...
    'FaceColor', color_superficie, ...
    'FaceAlpha', 0.7, ...
    'EdgeColor', 'k');

xlabel('X');
ylabel('Y');
zlabel('Z');

title('Pared exterior del cilindro');

view(45, 25);

hold off;