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

Debido a que el campo u solo depende de la coordenada θ, un campo ortoganl será el mismopero dependiendo de la coordenada ρ

[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 | R > Re);
Psi(mask) = NaN; 

%% Representación Gráfica
figure; 
hold on;

% Dibujar las líneas de corriente 
contour(X, Y, Psi, 20, 'r', 'LineWidth', 1.5); 

% Dibujar los cilindros (fronteras)
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('\Psi(\rho) = cte. Líneas de Corriente del campo U');
grid on;
hold off;

5 Puntos de velocidad máxima

5.1 Busqueda de la derivada

Para enconntrar los puntos dentro de los cilindros donde la velocidad de nuestro fluido es máxima, debemos estudiar el 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]

Realmente el campo solo depende de la variable ρ, podemos escribir: [math]\vec{u}(\rho) = \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]

Como lo que queremos analizar es el módulo, vamos a sacar nuestra función: [math]f(\rho) = \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}[/math]

Por último, antes de comenzar con el estudio, vamos a sustituir nuestros valores [math]|\omega_i|=|\omega_e|=1[/math] , para así obtener: [math]f(\rho) = -\rho + \frac{2}{\rho}[/math]

Ya con nuestra función terminada, podemos proceder a realizar su primera derivada para estudiarla: [math]f'(\rho)=\frac{df}{d\rho} = -1 - \frac{2}{\rho^2}[/math]

5.2 Estudio de la derivada

Primero, buscaremos los puntos criticos: [math]f'(\rho)=\frac{df}{d\rho} = -1 - \frac{2}{\rho^2}=0[/math]  ; [math]\rho^2 = -2[/math] . No existe solución real, por lo que podemos afirmar que la funcion no tendra puntos críticos en nuestro intervalo [1,2]. Además, si observamos nuestra función, nos damos cuenta que en nuestro intervalo [math]f'(\rho)\lt0[/math] siempre, por lo que el máximo estara en uno de los extremos del intervalo.

Sustituyendo en nuestra derivada obtenermos: [math]f'(1)=0[/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');

6.3 ¿Qué puntos tienen mayor rotacional?

[math] |\nabla \times \mathbf{u}| = \left| -\frac{4}{3} w_e - \frac{2}{3} w_i \right| [/math] y tomamos [math]w_e[/math]=-1 y [math]w_i[/math]= 1 como hemos hecho anteriormente.

El valor del rotacional es tan simple como calcular el módulo de este [math] |\nabla \times \mathbf{u}| = \left| -\frac{4}{3}(-1) - \frac{2}{3}(+1) \right| = \left| \frac{4}{3} - \frac{2}{3} \right| = \left| \frac{2}{3} \right| \approx 0.667 [/math] Esto significa que el resultado es constante, por lo que no hay puntos con mayor o menor rotacional sino que es el mismo en todos ellos.

6.4 ¿Es razonable?

Sí, es totalmente razonable. En este caso se trata de un flujo de couette cilindrico por lo que es entendible que no hayan puntos con mayor o menor rotacionalidad y que sea constante se debe a la uniformidad de la vorticidad. Esto indica que toda la masa encontrada entre ambas circunferencias sienten la misma intensidad de giro. También podemos dar por hecho que es constante ya que la gráfica tiene elmismo color en todas partes, y si fuera el caso contrario tendría distintas tonalidades.

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 por una sección

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:

Campo vectorial en función de [math]\rho[/math] y [math]\theta[/math]

Sea el campo dado por:

[math] \vec{u}(\rho,\theta) = \left(\frac{4}{3}\rho - \frac{4}{3}\rho\right)\,\vec{e}_{\theta}. [/math]

Simplificando:

[math] \vec{u}(\rho,\theta) = 0 \,\vec{e}_{\theta}. [/math]

Por componentes en coordenadas cilíndricas:

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

Así, el campo en función de [math]\rho[/math] y [math]\theta[/math] queda:

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

Calculamos el vector tangente:

[math] \vec{r}(\theta, z) = (R\cos\theta,\; R\sin\theta,\; z) [/math]

Calculamos los vectores tangentes derivando respecto a cada parámetro:

Derivada respecto a [math]\theta[/math] [math] \frac{\partial \vec{r}}{\partial \theta} = (-R\sin\theta,\; R\cos\theta,\; 0) [/math]

Derivada respecto a [math]z[/math] [math] \frac{\partial \vec{r}}{\partial z} = (0,\; 0,\; 1) [/math]

Por tanto, los vectores tangentes de la superficie son:

[math] \vec{r}_\theta = (-R\sin\theta,\; R\cos\theta,\; 0) [/math]

[math] \vec{r}_z = (0,\; 0,\; 1) [/math]

Al ser positivos es una parametrización correcta.

8.2 Calculo del Caudal

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

Se tiene la parametrización del cilindro de radio [math]R[/math] y altura [math]H[/math]:

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

y el campo vectorial dado:

[math] \vec{u}(\rho,\theta) = \left(\frac{4}{3}\rho - \frac{4}{3}\rho \right)\vec{e}_{\theta} = 0\,\vec{e}_{\theta}. [/math]

Vectores tangentes de la superficie

[math] \mathbf r_\theta = \frac{\partial \mathbf r}{\partial \theta} = (-R\sin\theta,\, R\cos\theta,\, 0) [/math]

[math] \mathbf r_z = \frac{\partial \mathbf r}{\partial z} = (0,\,0,\,1) [/math]

Vector normal a la superficie

[math] \mathbf n = \mathbf r_\theta \times \mathbf r_z = (R\cos\theta,\, R\sin\theta,\, 0) [/math]

[math] \hat{\mathbf n} = \frac{\mathbf n}{\|\mathbf n\|} = (\cos\theta,\,\sin\theta,\,0) [/math]

Cálculo del caudal

El caudal a través de la superficie lateral se calcula como:

[math] Q = \iint_S \vec u \cdot \hat{\mathbf n} \, dS [/math]

donde [math]dS = \|\mathbf r_\theta \times \mathbf r_z\|\,d\theta\,dz = R\,d\theta\,dz[/math].

Como [math]\vec u = 0[/math], se tiene:

[math] Q = 0 [/math]

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;