Flujo de Couette entre dos tubos concéntricos (Grupo 33)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Flujo de Couette entre dos tubos concéntricos. Grupo 33 |
| Asignatura | Teoría de Campos |
| Curso | 2025-26 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
El flujo de Couette es un fenómeno básico de la mecánica de fluidos que describe cómo se comporta un fluido situado entre dos superficies paralelas cuando una de ellas se desplaza respecto a la otra. Este flujo, sencillo y ampliamente aplicable, ha sido analizado en profundidad por su importancia en la modelización de sistemas mecánicos.
En este trabajo nos centramos en la forma más simple de este sistema: dos cilindros circulares concéntricos, uno interior y otro exterior, entre los cuales se encuentra un fluido viscoso newtoniano. Cuando los cilindros gira, el fluido responde desarrollando un campo de velocidad que depende del radio, generando un flujo que solo gira alrededor del eje, sin moverse hacia dentro, hacia fuera ni hacia arriba. Este fenómeno es un ejemplo de solución exacta de las ecuaciones de Navier–Stokes en geometría cilíndrica.
Lo interesante del flujo de Couette en cilindros concéntricos es que:
- Representa un flujo laminar estable para velocidades de rotación moderadas, lo que permite obtener soluciones analíticas explícitas.
- Sirve como modelo idealizado para numerosos sistemas reales, como rodamientos hidrodinámicos, reómetros rotacionales o dispositivos de mezcla.
- Ofrece una geometría donde las ecuaciones se simplifican, al no existir dependencia ni del ángulo θ ni de la coordenada axial z, lo que convierte al problema en esencialmente unidimensional.
El problema propuesto nos indica la existencia de dos cilindros concéntricos con un flujo de fluido incompresible a través de ellos, de manera que el exterior se mueve con velocidad angular constante 𝜔⃗𝑒 en sentido horario mientras que el interior se mueve con velocidad angular 𝜔⃗𝑖 en sentido contrario, es decir, antihorario. Si suponemos que ambos cilindros tienen su eje en 𝑂𝑋3 y pintamos la sección transversal (𝑥3 = 0), el cilindro exterior queda proyectado sobre la circunferencia 𝜌 = 2 y el interior sobre la circunferencia 𝜌 = 1
Contenido
- 1 Mallado de la sección transversal
- 2 Ecuación Navier-Stokes y obtención de la ecuación diferencial a resolver
- 3 Campo de velocidades
- 4 Líneas de corriente
- 5 Módulo de la velocidad del fluido
- 6 Rotacional de [math]\vec{u}[/math]
- 7 Dibujar el campo de temperaturas y las curvas de nivel.
- 8 Gradiente de la temperatura en coordenadas cilíndricas
- 9 Caudal que pasa por la sección [math]\theta = 0[/math]
- 10 Bibliografía
1 Mallado de la sección transversal
En la siguiente figura se representa la sección transversal de los dos tubos concéntricos con los ejes fijados en la región (𝜌, 𝜃) ∈ [0, 3] × [0, 2𝜋].
% MALLADO DE LA SECCIÓN TRANSVERSAL
rho=1:0.1:2; % INTERVALO DE RHO [1,2]
theta=0:0.05:2*pi+0.02; % INTERVALO DE THETA [0,2*PI]
[RHO,THETA]=meshgrid(rho,theta); % MATRIZ DE PARAMETROS
rho1=1;
rho2=2;
[RHO1,THETA1]=meshgrid(rho1,theta);
[RHO2,THETA2]=meshgrid(rho2,theta);
hold on
grid on % MALLADO EJES
x=RHO.*cos(THETA); % PARAM DE X
y=RHO.*sin(THETA); % PARAM DE Y
mesh(x,y,0*x,'EdgeColor','b') % DIBUJO DE LA FIGURA
t=rho1.*cos(theta); % CREACION CILINDRO INTERIOR
s=rho1.*sin(theta);
plot(t,s,'k','LineWidth',2)
t=rho2.*cos(theta); % CREACION CILINDRO EXTERIOR
s=rho2.*sin(theta);
plot(t,s,'k','LineWidth',2)
axis([-3,3,-2.5,2.5]) % EJES DEL DIBUJO
title('MALLADO DE LA SECCION TRANSVERSAL') % TITULO DE LA GÁFICA
xlabel('ρ');
ylabel('θ');
hold off
El principal objetivo de este apartado es la obtención de la ecuación diferencial a partir de la ecuación de Navier-Stokes. Una vez conocida esta expresión, será posible calcular las constantes a y b que permitan que la velocidad en la frontera del fluido coincida con la de los tubos interior y exterior.
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. Adicionalmente, el campo de velocidades tiene que cumplir la ecuación de Navier-Stokes estacionaria que sigue la siguiente expresión:
Observación: El valor µ se relaciona con el coeficiente de viscosidad del fluido.
Analizando la expresión, y teniendo en cuenta que la función tomada tiene de variable ρ, la cual es constante, el campo gradiente será nulo ( [math]∇p = 0[/math] ). Si además es despreciable el término convectivo (primer término de la fórmula) finalmente obtenemos que la expresión queda:
Observación: [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 únicamente de las condiciones de frontera en el recinto que se indique.
- 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:
Por la condición de que el campo de velocidades está proporcionado en coordenadas cilíndricas, es de mayor facilidad el cálculo con dichas coordenadas y por consecuente, la ecuación quedaría de la siguiente forma:
- 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. Por definición ésta se calcula en coordenadas cilíndricas mediante la siguiente expresión:
Tras realizar los cálculos, observamos que como hemos deducido, la divergencia es nula y por lo tanto, el gradiente de la propia divergencia saldrá 0, de igual forma:
Una vez deducido que el gradiente de la divergencia es 0, por lo cálculos previamente realizados, el primer término se elimina y la expresión queda simplificada de la siguiente forma:
- Rotacional del campo de velocidades
Realizamos los cálculos aplicando la fórmula para el rotacional definido en coordenadas cilíndricas. Siendo la expresión del campo de velocidades: [math]u(\rho,\theta)=f(\rho)\,e_{\theta}[/math]
- Rotacional del rotacional del campo de velocidades
Una vez realizado el rotacional del campo de velocidades, proseguimos con el cálculo del rotacional del rotacional. Para ello aplicaremos la fórmula en coordenadas cilíndricas sobre el vector hallado anteriormente:
- Resultado del Laplaciano vectorial
Sustituimos en la ecuación obtenida del laplaciano vectorial.
- 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:
2.2 Obtención de las constantes de la ecuación diferencial
2.2.1 Obtención ecuación diferencial segundo orden
Una vez conocida la solución de la ecuación de Navier-Stokes, procedemos a la obtención de la ecuación diferencial:
Reordenamos la ecuación para la obtención de una ecuación diferencial de segundo orden:
2.2.2 Comprobación ecuación diferencial
A continuación, vamos a comprobar si la ecuación diferencial que hemos obtenido satisface la solución conocida, que viene dada por:
Para realizar la comprobación, es necesario derivar la expresión de forma que dichas derivadas aparezcan en la ecuación diferencial obtenida.
Una vez halladas las derivadas, sustituimos en la ecuación y comprobamos que:
Queda por tanto comprobado que se cumple la ecuación diferencial.
2.2.3 Búsqueda de las constantes "a" y "b"
Finalmente, hemos solucionado la EDO y debemos conocer las soluciones específicas de ésta, a través de las constante a y b de tal forma que la velocidad en la frontera coincida con las de los tubo interior y exterior sen la misma.
Sabemos que la velocidad del fluido es: [math] \mathbf{u}(\rho,\theta) = f(\rho)\,\mathbf{e}_\theta, \qquad f(\rho) = a\rho + \frac{b}{\rho}. [/math]
-En el cilindro interior ( [math]\rho = 1[/math] ), la velocidad tangencial será: [math] u_\theta(1) = f(1) = a(1) + \frac{b}{1} = a + b, [/math] debe coincidir con la velocidad tangencial del cilindro interior: [math] u_\theta(1) = \omega_i \cdot 1 = \omega_i. [/math]
-En el cilindro exterior ( [math]\rho = 2[/math] ),la velocidad tangencial: [math] u_\theta(2) = f(2) = a(2) + \frac{b}{2} = 2a + \frac{b}{2}, [/math] debe coincidir con la velocidad tangencial del cilindro exterior: [math] u_\theta(2) = \omega_e \cdot 2 = 2\omega_e. [/math]
- Sistema de ecuaciones
De las condiciones de contorno: [math] \begin{cases} a + b = \omega_i, \\ 2a + \dfrac{b}{2} = 2\omega_e. \end{cases} [/math]
Multiplicamos la segunda ecuación por 2 para eliminar la fracción: [math] 4a + b = 4\omega_e. [/math]
Restamos la primera ecuación: [math] (4a + b) - (a + b) = 4\omega_e - \omega_i \;\;\Longrightarrow\;\; 3a = 4\omega_e - \omega_i. [/math]
Por tanto: [math] a = \dfrac{4\omega_e - \omega_i}{3}. [/math]
Luego: [math] b = \omega_i - a = \omega_i - \dfrac{4\omega_e - \omega_i}{3} = \dfrac{3\omega_i - (4\omega_e - \omega_i)}{3} = \dfrac{4\omega_i - 4\omega_e}{3} = \dfrac{4(\omega_i - \omega_e)}{3}. [/math]
Resultado: [math] \boxed{ a = \dfrac{4\omega_e - \omega_i}{3}, \qquad b = \dfrac{4(\omega_i - \omega_e)}{3}. } [/math]
2.3 Comprobación de incompresibilidad
Un fluido incompresible es aquel cuya densidad es constante, esto matemáticamente se analiza obteniendo la divergencia del campo de velocidades, esta debe de ser cero (∇ ⋅ v = 0).
Se ha comprobado con anterioridad, a la hora de la obtención de la ecuación diferencial, que la divergencia del campo de velocidades es igual a 0. Aun así, se detalla a continuación el cálculo de la incompresibilidad:
observación: se toman coordenadas cilíndricas.
[math]
\nabla \cdot \mathbf{u}
= \frac{1}{\rho}\frac{\partial(\rho u_\rho)}{\partial \rho}
+ \frac{1}{\rho}\frac{\partial u_\theta}{\partial \theta}
+ \frac{\partial u_z}{\partial z}.
[/math]
Como
[math]
u_\rho = 0,\quad u_z = 0,\quad
u_\theta = f(\rho)\ \text{(independiente de }\theta,z\text{)},
[/math]
entonces:
[math]
\nabla \cdot \mathbf{u}
= \frac{1}{\rho}\frac{\partial(0)}{\partial \rho}
+ \frac{1}{\rho}(0) + 0
= 0.
[/math]
3 Campo de velocidades
Con el objetivo de calcular el campo de velocidades establecemos que la velocidad angular del cilindro exterior e interior es 1 (ω = 1) y μ = 1 (coeficiente de viscosidad del fluido), y en consecuencia las constantes serán de la siguiente forma: [math] \boxed{ a = \dfrac{4\omega_e - \omega_i}{3}, \qquad b = \dfrac{4(\omega_i - \omega_e)}{3}. } [/math]
Observación: se considera que por el sentido, ωi = 1 y ωe = -1.
Conocida la siguiente expresión de la velocidad de las partículas en función de f(ρ):
[math]\vec{u}(ρ, θ) = f(ρ) \vec{e_θ} [/math].
Hemos calculado las constantes a y b, desarrollamos la ecuación. Quedando el campo de velocidades la las partículas como:
[math]\vec{u}(\rho, \theta) = \left(-\frac{5}{3} \rho + \frac{8}{3 \rho}\right) \hat{e}_\theta [/math]
%Apartado 3: Campo de velocidades
rho=linspace(1,2,30); %radio de 1 a 2
theta=linspace(0,2*pi,50); %theta de 0 a 2pi
a=-5/3;
b=8/3;
w1=1; %velocidad angular interior
w2=-1; %velocidad angular exterior
%Mallado coordenadas polares
[R,T]=meshgrid(rho,theta); % meshgrid para coordenadas polares
X=R.*cos(T); % Paso a cartesianas
Y=R.*sin(T);
%Campo de velocidades
Vtheta= a.*R+b./R; %Velocidad de theta
Vx= w2.*Vtheta .*sin(T);
Vy= w1.*Vtheta.*cos(T);
%Dibujo campo de velocidades
hold on
grid on
quiver(X,Y,Vx,Vy,0.5,'Color',[0 0 1],'LineWidth',1.2);
title('CAMPO DE VELOCIDADES')
axis([-3.5,3.5,-2.1,2.1])
hold off
4 Líneas de corriente
Para obtener las lineas de corriente del campo [math]\vec{u}[/math] conocemos que el campo de velocidades del fluido en la sección transversal es [math]\vec{u}(ρ,\theta)=f(ρ)\,\vec{e_\theta}[/math], donde [math]f(\rho) = a\rho + \frac{b}{\rho}[/math] y [math]\vec{e_\theta}=(-\sin\theta,\ \cos\theta)[/math]. Dichas lineas de corriente se obtienen de forma más sencilla recurriendo a un campo de velocidades ortogonal a [math]\vec{u}[/math], el cual, hemos llamado [math]\vec{v}[/math].
4.1 Obtención del campo de velocidades [math]\vec{v}[/math]
Sabemos que el campo de velocidades [math]\vec{v}[/math] ha de ser perpendicular a [math]\vec{u}[/math].
Para ello, tomamos el siguiente producto vectorial: [math]\vec{v}=\vec{k}\times\vec{u}[/math], siendo [math]\vec{k}=(0,0,1)[/math] el vector unitario correspondiente al eje Z en base cartesiana. Usando propiedades del producto vectorial: [math]\vec{k}\times\vec{e_\theta}=-\vec{e_ρ}[/math], donde [math]\vec{e_ρ}=(\cos\theta,\ \sin\theta)[/math].
Por tanto, en base cilíndrica, el campo de velocidades seguiría la siguiente expresión:
Cambiando a coordenadas cartesianas:
con [math]ρ=\sqrt{x^2+y^2}[/math], y [math]\theta=\operatorname{arctan}(y,x)[/math].
- Componentes en la base cilíndrica
[math]v_ρ = -\Big(a\,ρ + \tfrac{b}{ρ}\Big), \quad v_\theta = 0, \quad v_z = 0[/math]
4.2 Observación de [math]\vec{v}[/math] irrotacional
Se quiere comprobar la irrotacionalidad del campo v, ya que, al igual que el campo u, este no debe de rotar sobre si mismo. Conocemos entonces que [math]\vec{v}[/math] sigue la dirección radial, pues apunta hacia [math]-\vec{e_ρ}[/math]y es ortogonal a [math]\vec{u}[/math] porque [math]\vec{e_ρ}\cdot\vec{e_\theta}=0[/math]:
- Irrotacionalidad de v
Para calcular si el campo es rotacional o no, iremos a su definición. Un campo dado en coordenadas cilíndricas será irrotacional, siempre que el rotacional de 0:
[math] \nabla \times \vec{v}(\rho, \theta, z) = \frac{1}{\rho} \begin{vmatrix} \mathbf{e}_{\rho} & \rho \mathbf{e}_{\theta} & \mathbf{e}_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ u_{\rho} & \rho u_{\theta} & u_{z} \end{vmatrix} = 0 [/math]
En coordenadas cilíndricas, para un campo que tiene únicamente coordenada [math]\mathbf{e}_\rho[/math]
[math]v_ρ(ρ)[/math] con [math]v_\theta = v_z = 0[/math], el rotacional es:
[math]\nabla \times \vec{v} = \begin{cases} \text{componente $z$: } \tfrac{1}{ρ}\tfrac{\partial}{\partial ρ}(ρ v_\theta) - \tfrac{\partial v_ρ}{\partial \theta} = 0, \\ \text{otras componentes: también } 0 \end{cases}[/math]
Como [math]v_\theta = 0[/math] y [math]v_ρ[/math] no depende de [math]\theta[/math], se cumple: [math]\nabla \times \vec{v} = \vec{0}[/math]
Por tanto, [math]\vec{v}[/math] es irrotacional.
4.3 Función de corriente [math]\psi(ρ)[/math]
Existe un potencial escalar [math]\psi[/math] tal que:
[math]\nabla \psi = \vec{v}[/math]
En coordenadas cilíndricas:
[math]\nabla \psi = \frac{\partial \psi}{\partial ρ}\,\vec{e_ρ} + \frac{1}{ρ}\frac{\partial \psi}{\partial \theta}\,\vec{e_\theta} + \frac{\partial \psi}{\partial z}\,\vec{e_z}[/math]
Como [math]\vec{v} = v_ρ\,\vec{e_ρ}[/math] con [math]v_ρ = -f(ρ)[/math], tenemos:
[math]\frac{\partial \psi}{\partial ρ} = -f(ρ), \quad \frac{\partial \psi}{\partial \theta} = 0, \quad \frac{\partial \psi}{\partial z} = 0[/math]
Por tanto:
[math]\psi(ρ) = -\int f(ρ)\,dρ = -\int \left(a ρ + \frac{b}{ρ}\right)\,dρ = -\left(\tfrac{a}{2}ρ^2 + b\ln ρ\right) + C[/math]
Finalmente:
[math]\boxed{\psi(ρ) = -\tfrac{a}{2}ρ^2 - b\ln ρ + C}[/math]
4.4 Comprobación de que las líneas 𝜓 = cte
Como 𝜓 depende solo de [math]ρ[/math], las curvas [math]\psi = \text{cte}[/math] son circunferencias concéntricas: [math]ρ = \text{constante}[/math].
Estas son líneas de corriente de [math]\vec{u}[/math] porque: [math]\vec{u} \cdot \nabla \psi = f(ρ)\,\vec{e_\theta} \cdot \big(-f(ρ)\,\vec{e_ρ}\big) = 0[/math].
Además, se ha calculado: [math]\psi(ρ) = -\tfrac{a}{2}ρ^2 - b \ln ρ + C[/math], [math]f(ρ) = a ρ + \tfrac{b}{ρ}[/math].
Basta con elegir varios valores de [math]ρ[/math] en el intervalo [math][1,2][/math] y dibujar los círculos correspondientes.
%% Parámetros del problema
r_in = 1;
r_out = 2;
omega_i = +1;
omega_e = -1;
% Coeficientes del perfil
a = (4*omega_e - omega_i)/3; % = -5/3
b = 4*(omega_i - omega_e)/3; % = 8/3
%% Función de corriente psi(r) = -(a/2)*r^2 - b*ln(r)
psi = @(r) -(a/2)*r.^2 - b*log(r);
%% Elegimos varios niveles de psi (líneas de corriente)
n_lines = 25; % número de círculos
r_vals = linspace(r_in, r_out, n_lines);
%% Dibujar
figure('Color','w'); hold on; axis([-3 3 -2.5 2.5]); grid on;
title('Líneas de corriente (\psi = cte)');
xlabel('x'); ylabel('y');
% Dibujar círculos para cada r
theta = linspace(0, 2*pi, 400);
for r = r_vals
x = r*cos(theta);
y = r*sin(theta);
plot(x, y, 'Color', [0.2 0.6 0.8], 'LineWidth', 1.2);
end
5 Módulo de la velocidad del fluido
El módulo de la velocidad del fluido viene dado por: [math] |\mathbf{u}| = |u_\theta| = |f(ρ)| = \big| a\,ρ + \frac{b}{ρ} \big|. [/math]
Con [math] a = -\frac{5}{3}, \qquad b = \frac{8}{3}, [/math]
tenemos: [math] f(ρ) = -\frac{5}{3}ρ + \frac{8}{3ρ}. [/math]
- Buscar el máximo en [math]ρ \in [1,2][/math]
Derivamos la expresión para obtener los tramos de crecimiento y decrecimiento, y por tanto conocer los puntos máximos y/o mínimos que puedan existir. [math] \frac{d}{dρ} f(ρ) = -\frac{5}{3} - \frac{8}{3ρ^2}. [/math]
Por tanto, procedemos a calcular los máximos en los extremos del intervalo:
En [math]ρ = 1[/math]: [math] f(1) = -\frac{5}{3}(1)+\frac{8}{3(1)} = -\frac{5}{3}+\frac{8}{3} = 1. [/math]
En [math]ρ = 2[/math]: [math] f(2) = -\frac{10}{3}+\frac{8}{6} = -\frac{10}{3}+\frac{4}{3} = -2. [/math]
El módulo: [math] |f(1)| = 1,\quad |f(2)| = 2. [/math]
En conclusión, el máximo está en la frontera exterior [math]ρ = 2[/math], con: [math] |\mathbf{u}| = 2. [/math]
El módulo de la velocidad es máximo en |f(2)| = 2
- Representación gráfica del módulo de la velocidad
En la siguiente figura se representa el módulo de la velocidad con respecto en función del parámetro ρ:
clear; clc;
% Parámetros
r1 = 1;
r2 = 2;
w1= +1;
w2 = -1;
% Coeficientes del perfil
a = (4*w2- w1)/3; % = -5/3
b = 4*(w1 - w2)/3; % = 8/3
% Vector radial y cálculo del módulo
r = linspace(r1, r2, 200);
u = abs(a*r + b./r); %
6 Rotacional de [math]\vec{u}[/math]
6.1 Expresión del rotacional y cálculo del mismo
Por definición el rotacional de un campo vectorial en coordenadas cilíndricas se calcula como: [math] \nabla \times \vec{u} = \frac{1}{\rho} \begin{vmatrix} \mathbf{e}_{\rho} & \rho \mathbf{e}_{\theta} & \mathbf{e}_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ u_{\rho} & \rho u_{\theta} & u_{z} \end{vmatrix} [/math] De datos de partida conocemos que: [math] u_ρ = 0,\; u_\theta = f(ρ),\; u_z = 0,\; \text{y } f(ρ) = a ρ + \frac{b}{ρ}, [/math], y sustituimos en la expresión anterior: [math] \nabla \times \vec{u} = \frac{1}{\rho} \begin{vmatrix} \mathbf{e}_{\rho} & \rho \mathbf{e}_{\theta} & \mathbf{e}_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ 0 & \rho\left(a\rho + \frac{b}{\rho}\right) & 0 \end{vmatrix} [/math] Desarollando el determinante, obtenemos que el rotacional sigue la siguiente expresión: [math] (\nabla \times \vec{u}) = \frac{1}{\rho} \frac{\partial (\rho f(\rho))}{\partial \rho} \mathbf{e}_z [/math]
- Calculamos [math]ρ f(ρ)[/math]
[math] ρf(ρ) = ρ\big(a ρ + \frac{b}{ρ}\big) = a ρ^2 + b. [/math]
Derivamos respecto de ρ: [math] \frac{\partial}{\partial ρ}(a ρ^2 + b) = 2 a ρ. [/math]
Por tanto, obtenemos que: [math] (\nabla \times \vec{u}) = \frac{1}{\rho} (2 a \rho) \mathbf{e}_z = 2 a \mathbf{e}_z [/math]
- Resultado del rotacional
El rotacional del propio campo de velocidades quedaría expresado de la siguiente manera: [math] \boxed{\nabla \times \mathbf{u} = (0,0,2a)}. [/math]
Observamos que es constante y únicamente depende de [math]a[/math].
Con [math]a = -\frac{5}{3}[/math], obtenemos que el rotacional es: [math] \nabla \times \mathbf{u} = (0,0,-\frac{10}{3}) \approx (0,0,-3.33). [/math]
6.2 ¿Qué puntos tienen mayor rotacional?
Todos los puntos tienen el mismo rotacional, porque el flujo depende únicamente de θ en la base de coordenadas cilíndricas y la derivada radial genera un valor constante.
6.3 ¿Es razonable que salgan esos puntos?
Sí, porque este flujo es un flujo de Couette cilíndrico con velocidad angular variable, pero su rotacional depende únicamente del gradiente de la velocidad angular. el rotacional en el anillo es uniforme, representando que el fluido gira en el interior de los tubos concéntricos con un gradiente constante.
6.4 Representación gráfica del rotacional del campo de velocidades de las partículas.
En el gráfico posterior se expresa el rotacional del campo de velocidades de partículas [math]\vec{u}[/math]
%% Cálculo del rotacional (constante)
A = 10/3;
curl_value = abs(2*A); % ∇×u
%% Mallado para representación
Nradios = 250;
Ntheta = 350;
rho = linspace(1, 2, Nradios);
theta = linspace(0, 2*pi, Ntheta);
[Rho, Theta] = meshgrid(rho, theta);
% El rotacional es constante → matriz uniforme
CURL = curl_value * ones(size(Rho));
%% Conversión a cartesianas
X = Rho .* cos(Theta);
Y = Rho .* sin(Theta);
%% Gráfica
figure;
pcolor(X, Y, CURL);
shading interp;
colormap parula;
title('rotacional de u');
axis([-3 3 -2.5 2.5]);
7 Dibujar el campo de temperaturas y las curvas de nivel.
Trabajaremos en la sección transversal y dentro del fluido (anillo [math]1 \le \rho \le 2[/math]).
La temperatura está dada por la siguiente fórmula: [math] T(\rho,\theta) = \log\!\big(1+\rho^2\big)\,\cos^2\theta. [/math]
7.1 Representación del campo de temperaturas y obtención gráfica del punto máximo
%Apartado 7: Campo de temperaturas
rho=linspace(1,2,200); %radio de 1 a 2
theta=linspace(0,2*pi,200); %theta de 0 a 2pi
[R,Theta]=meshgrid(rho,theta);
% función temperatura
T=log(1+R.^2).*cos(Theta).^2;
X=R.*cos(Theta); % Paso a cartesianas
Y=R.*sin(Theta);
%Dibujar campos de temperatura
hold on
surf(X, Y, T, 'EdgeColor', 'none');
colorbar;
xlabel('eje x');
ylabel('eje y');
zlabel('Temperatura');
title('Campo de Temperatura');
view(3);
axis([-3,3,-2.1,2.1]);
%Dibujar las curvas de nivel
contourf(X, Y, T, 20); % 20 niveles de contorno
title('Curvas de Nivel del Campo de Temperaturas');
xlabel('Eje X');
ylabel('Eje Y');
colorbar;
%Encontrar temperatura máxima
T_max=max(T(:));
[f,c] =find(T == max_T, 1);
% Coordenadas cartesianas Temp máx
x_max= X(f,c);
y_max= Y(f,c);
plot(x_max, y_max, 'ro', 'MarkerSize', 10, 'LineWidth', 2);
text(x_max,y_max,'Temperatura máxima');Se ha obtenido el MÁXIMO DE TEMPERATURA de forma gráfica, para su mejor comprensión y para su comprobación se realizará igualmente de froma analítica.
7.2 Obtención del punto máximo de temperatura analíticamente
- Comportamiento de los factores
Debido a la complejidad de derivación de la expresión resulta de mayor utilidad analizar los componentes que la conforman.
1) Para [math]\rho \ge 0[/math], (sabemos que ρ está entre 1 y 2) la función: [math] \log(1+\rho^2) [/math] es estrictamente creciente en [math]\rho[/math] (pues [math]\frac{d}{d\rho}\log(1+\rho^2) = \frac{2\rho}{1+\rho^2} \ge 0[/math] y es >0 para [math]\rho\gt0[/math]).
2) El factor angular: [math] \cos^2\theta [/math] alcanza su máximo igual a [math]1[/math] cuando [math]\cos\theta = \pm 1[/math], es decir, para aquellos valores: [math]\theta = 0,\,\pi,\,2\pi,\,\dots[/math]
- Máximo en el dominio [math]1 \le \rho \le 2[/math]
Para maximizar [math]T(\rho,\theta) = \log(1+\rho^2)\cos^2\theta[/math]:
- Tomamos el mayor valor radial permitido: [math]\rho = 2[/math] (porque [math]\log(1+\rho^2)[/math] crece con [math]\rho[/math]).
- Tomamos el mayor valor angular: [math]\cos^2\theta = 1[/math], que ocurre únicamente en los valores [math]\theta = 0[/math] o [math]\theta = \pi[/math].
Por tanto, el punto máximo de [math]T[/math] en el anillo se da en: [math] \rho = 2,\qquad \theta = 0 \text{ o } \pi. [/math]
8 Gradiente de la temperatura en coordenadas cilíndricas
8.1 Cálculo del gradiente de T
En dicho apartado se procederá al cálculo del gradiente de el campo de temperaturas dado por la siguiente expresión:
[math]
T(ρ,\theta) = \log(1+ρ^2)\,\cos^2\theta,
[/math]
El gradiente de un campo vectorial, como es el caso, en coordenadas cilíndricas, viene expresado por la posterior definición:
[math]
\nabla T =
\frac{\partial T}{\partial \rho} \mathbf{e}_{\rho}
+ \frac{1}{\rho} \frac{\partial T}{\partial \theta} \mathbf{e}_{\theta}
+ \frac{\partial T}{\partial z} \mathbf{e}_{z}
[/math]
Por tanto, es necesario el cálculo de las derivadas parciales del campo vectorial:
[math] \frac{\partial T}{\partial ρ} = \frac{2ρ}{1+ρ^2}\,\cos^2\theta, \qquad \frac{\partial T}{\partial \theta} = \log(1+ρ^2)\,\frac{d}{d\theta}(\cos^2\theta) = -\,\log(1+ρ^2)\,\sin(2\theta). [/math]
En base cilíndrica, el gradiente, al no depender el campo de temperaturas de [math]\mathbf{e}_z[/math] se expresa como:
[math]
\nabla T = \frac{\partial T}{\partial ρ}\,\mathbf{e}_ρ + \frac{1}{ρ}\frac{\partial T}{\partial \theta}\,\mathbf{e}_\theta
= T_ρ\,\mathbf{e}_ρ + T_\theta\,\mathbf{e}_\theta,
[/math]
Con los cálculos obtenidos previamente de las derivadas parciales, sustituimos en la expresión:
[math] \nabla T = T_{\rho} \mathbf{e}_{\rho} + T_{\theta} \mathbf{e}_{\theta} = \left( 2 \rho (1 + \rho^2 \cos 2\theta) \right) \mathbf{e}_{\rho} + \left( -\frac{1}{\rho} \log(1 + \rho^2) \sin(2\theta) \right) \mathbf{e}_{\theta} [/math]
8.2 Representación gráfica del campo gradiente y comprobación de ortogonalidad
En las siguientes imágenes se representa en primer lugar el campo gradiente; y posteriormente, la comprobación de perpendicularidad del campo gradiente con respecto a las líneas de nivel de la temperatura.
r_in = 1; % Geometría del dominio (anillo)
r_out = 2;
xL = -2.2; xR = 2.2; % Malla cartesiana para calcular/visualizar
N = 401;
yB = -2.2; yT = 2.2;
[x, y] = meshgrid(linspace(xL,xR,N), linspace(yB,yT,N));
r = hypot(x, y); % Coordenadas polares
theta = atan2(y, x);
T = log(1 + r.^2) .* (cos(theta).^2); % Temperatura
Tr = (2.*r ./ (1 + r.^2)) .* (cos(theta).^2); % dT/dr
Ttheta = - log(1 + r.^2) .* sin(2*theta); % dT/dtheta
% Componentes de e_r y e_theta
er_x = cos(theta); er_y = sin(theta);
et_x = -sin(theta); et_y = cos(theta);
gx = Tr .* er_x + (Ttheta ./ r) .* et_x;
gy = Tr .* er_y + (Ttheta ./ r) .* et_y;
gx(r==0) = 0; gy(r==0) = 0; %(fuera del anillo no dibujaremos)
mask = (r >= r_in) & (r <= r_out);
% Figura 1:
figure('Color','w'); hold on; axis equal; box on; grid on;
levels = 10; % nº de niveles de contorno
% Curvas de nivel de T en el anillo
T_plot = T; T_plot(~mask) = NaN;
contour(x, y, T_plot, levels, 'LineWidth', 1.0);
% Campo gradiente (submuestreo para legibilidad)
step = 16;
qx = x(1:step:end, 1:step:end);
qy = y(1:step:end, 1:step:end);
qgx = gx(1:step:end, 1:step:end);
qgy = gy(1:step:end, 1:step:end);
qmask = (hypot(qx, qy) >= r_in) & (hypot(qx, qy) <= r_out);
qgx(~qmask) = 0; qgy(~qmask) = 0;
quiver(qx, qy, qgx, qgy, 0.8, 'Color', [1 0 0], 'LineWidth', 0.8); % gradiente
title('Curvas de nivel de T y gradiente');
xlabel('x'); ylabel('y');
xlim([xL xR]); ylim([yB yT]);
%figura 2
figure('Color','w'); hold on; axis equal; box on; grid on;
step = 16;
qx = x(1:step:end, 1:step:end);
qy = y(1:step:end, 1:step:end);
qgx = gx(1:step:end, 1:step:end);
qgy = gy(1:step:end, 1:step:end);
qmask = (hypot(qx, qy) >= r_in) & (hypot(qx, qy) <= r_out);
qgx(~qmask) = 0; qgy(~qmask) = 0;
quiver(qx, qy, qgx, qgy, 0.8, 'Color', [1 0 0], 'LineWidth', 0.8); % gradiente
title('Gradiente');
xlabel('x'); ylabel('y');
xlim([xL xR]); ylim([yB yT]);
9 Caudal que pasa por la sección [math]\theta = 0[/math]
9.1 Cálculo del caudal
La sección [math]\theta = 0[/math] es una línea radial desde [math]ρ = 1[/math] hasta [math]ρ = 2[/math].
El fluido se mueve dependiente de [math]u_\theta[/math], así que el flujo atraviesa esa línea perpendicularmente.
El caudal [math]Q[/math] a través de esa sección (profundidad [math]h = 1 \,\text{m}[/math]) se calcula como:
[math]
Q = \int_{ρ=1}^{ρ=2} u_\theta(ρ)\, h \, dρ,
[/math]
porque el área diferencial es [math]dA = h \, dρ[/math] (una tira de ancho [math]dρ[/math] y profundidad [math]h[/math]).
[math]
u_\theta(ρ) = f(ρ) = a ρ + \frac{b}{ρ}, \qquad a = -\frac{5}{3},\; b = \frac{8}{3}.
[/math]
[math]
Q = \int_{1}^{2} \big(a ρ + \frac{b}{ρ}\big)\, dρ
= \Big[\frac{a}{2}ρ^2 + b \ln ρ\Big]_{1}^{2}.
[/math]
Evaluamos en los valores de a y b obtenidos en el apartado 2:
[math]
\left[ \frac{a}{2} \rho^2 + b \ln \rho \right]_{2}^{1} =
\left[ -\frac{5}{6} \rho^2 + \frac{8}{3} \ln \rho \right]_{2}^{1}
[/math]
[math]
Q = \Big(-\frac{5}{6}(2^2) + \frac{8}{3}\ln 2\Big) - \Big(-\frac{5}{6}(1) + \frac{8}{3}\ln 1\Big)
= -\frac{20}{6} + \frac{8}{3}\ln 2 + \frac{5}{6}
= -\frac{15}{6} + \frac{8}{3}\ln 2
= -2.5 + \frac{8}{3}\ln 2.
[/math]
El resultado final del caudal quedaría de la siguiente forma:
[math] Q \approx -2.5 + 1.848 = -0.652 \;\text{m}^3/\text{s}. [/math]
El signo negativo indica que el flujo atraviesa la sección en sentido horario (hacia abajo respecto a la orientación positiva).
9.2 Representación Gráfica
En las posteriores figuras se representa el caudal del campo visto desde planta (con objetivo de conocer el sentido del movimiento) y vista de perfil (para poder observar el movimiento real del caudal dentro de los dos tubos concéntricos).
%% Parámetros geométricos
r_in = 1; % radio interior
r_out = 2; % radio exterior
h = 1; % profundidad (m)
%% Velocidades angulares (interior antihorario, exterior horario)
omega_i = +1;
omega_e = -1;
%% Perfil azimutal: u_theta(r) = a*r + b/r
a = (4*omega_e - omega_i)/3; % = -5/3
b = 4*(omega_i - omega_e)/3; % = 8/3
f = @(r) a*r + b./r;
%% Caudal a través de la sección theta=0: Q = ∫_{r=1}^{2} u_theta(r) * h dr
Q = (a/2)*(r_out^2 - r_in^2) + b*log(r_out/r_in); % integral analítica
%% Geometría para dibujar
th = linspace(0, 2*pi, 500);
figure('Color','w'); hold on; axis equal; grid on; box on;
xlim([-2.4 2.4]); ylim([-2.0 2.0]);
xlabel('x'); ylabel('y');
title(sprintf('Sección \\theta=0 y flujo azimutal atravesándola (Q = %.3f m^3/s)', Q));
% Dibujar el anillo del fluido
plot(r_in*cos(th), r_in*sin(th), 'k', 'LineWidth', 1.5);
plot(r_out*cos(th), r_out*sin(th), 'k', 'LineWidth', 1.5);
%% Sección theta=0: línea radial sobre el eje x desde r=1 a r=2
plot([r_in, r_out], [0, 0], 'r', 'LineWidth', 2); % sección en magenta
%% Flechas del flujo atravesando la sección (u_theta es perpendicular a la sección)
n_arrows = 10; % nº de flechas sobre la sección
r_samples = linspace(r_in, r_out, n_arrows);
x_s = r_samples; y_s = zeros(size(r_samples));
u_x = zeros(size(r_samples));
u_y = f(r_samples); % positivo: flecha hacia +y; negativo: hacia -y
% Escalado de flechas para visibilidad
scale = 0.9;
quiver(x_s, y_s, u_x, u_y, scale, 'Color', [0 0 1], 'LineWidth', 1.4);
%% Referencias y anotaciones
% Marcar valores en los bordes
plot(r_in, 0, 'ko', 'MarkerFaceColor','k'); text(r_in, -0.12, 'r=1', 'HorizontalAlignment','center');
plot(r_out,0, 'ko', 'MarkerFaceColor','k'); text(r_out, -0.12, 'r=2', 'HorizontalAlignment','center');
% Leyenda conceptual
%legend({'rlegend({'r=1','r=2','sección \theta=0','flujo azimutal'}, 'Location','bestoutside');
10 Bibliografía
Microsoft. (2025). Microsoft Copilot. Recuperado de https://copilot.microsoft.com
Wikimedia Foundation. (2025). Wikipedia. Recuperado de https://www.wikipedia.org
Apuntes Tema 2. Campos tensoriales y operadores diferenciales. Proporcionados por los profesores de la asignatura
