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	Carlos Azcona Ruiz Alonso Gómez Oliveira Marcos Gómez Rosillo María de las Mercedes Oriola Zamora Alba Zilu Martínez Toda
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 artículo vamos a considerar el flujo de un fluido incompresible a través de dos cilindros concéntricos, 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. 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

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

% 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

2 Ecuación Navier-Stokes y obtención de la ecuación diferencial a resolver

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:

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

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:

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

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:

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

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:

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

• 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:

[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, el gradiente de la propia divergencia saldrá 0, de igual forma:

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

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:

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

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

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

• 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:

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

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

• 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.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:

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

[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.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:

[math]f(ρ) = aρ +\frac{b}{ρ},\hspace{20pt}a,b \in \mathbb{R} [/math]

Para realizar la comprobación, es necesario derivar la expresión de forma que dichas derivadas aparezcan en la ecuación diferencial obtenida.

[math]\frac{\partial f(ρ)}{\partial ρ} = a -\frac{b}{ρ^2},\hspace{20pt}a,b \in \mathbb{R} [/math] [math]\frac{\partial }{\partial ρ}(ρ\frac{\partial f(ρ)}{\partial ρ}) = \frac{\partial }{\partial ρ}(ρ(a -\frac{b}{ρ^2})) = \frac{\partial }{\partial ρ}(aρ -\frac{b}{ρ})) = a +\frac{b}{ρ^2},\hspace{20pt}a,b \in \mathbb{R}[/math]

Una vez halladas las derivadas, sustituimos en la ecuación y comprobamos que:

[math] \frac{1}{ρ}(aρ +\frac{b}{ρ}) = a +\frac{b}{ρ^2} \Longrightarrow \frac{f(ρ)}{ρ} = \frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}),\hspace{20pt}a,b \in \mathbb{R} [/math]

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

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. Se detalla a continuación el cálculo de la incompresibilidad: observación: se tomarán 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: Observación: se considera que por el sentido, ω_i = 1 y ω_e = -1. [math] \boxed{ a = \dfrac{4\omega_e - \omega_i}{3}, \qquad b = \dfrac{4(\omega_i - \omega_e)}{3}. } [/math]

[math]a = -\frac{5}{3} [/math] y [math] b = \frac{8}{3}[/math]

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]

Campo de velocidades

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

• Obtención del campo de velocidades v

Sabemos que el campo de velocidades [math]\vec{v}[/math] ha de ser ortogonal a [math]\vec{u}[/math].

1-Tomamos el siguiente producto vectorial, para obtener el vector perpendicular [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.

2-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:

[math]\vec{v}(ρ,\theta) = f(ρ)\,(\vec{k}\times\vec{e_\theta}) = -\,f(ρ)\,\vec{e_ρ}[/math]

Cambiando a coordenadas cartesianas:

[math]\vec{v}(x,y)=\big(-f(ρ)\cos\theta,\ -f(ρ)\sin\theta\big)[/math]

con [math]ρ=\sqrt{x^2+y^2}[/math], y [math]\theta=\operatorname{arctan2}(y,x)[/math].

3-Se observa 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]:

[math] \vec{u}\cdot\vec{v} = f(ρ)\,\vec{e_\theta}\cdot\big(-f(ρ)\,\vec{e_ρ}\big) = -f(ρ)^2(\vec{e_\theta}\cdot\vec{e_ρ}) = 0. [/math]

4-En nuestro caso, la función seguiría la siguiente expresión [math]f(ρ)=a\,ρ+\frac{b}{ρ}[/math].

5- Se fijaron previamente las constantes [math]a,b[/math] para [math]\omega_i = +1[/math], [math]\omega_e = -1[/math]: [math]a = -\tfrac{5}{3}[/math], [math]b = \tfrac{8}{3}[/math], entonces: [math] f(ρ)=a\,ρ+\frac{b}{ρ} [/math] y el campo quedaría expresado en coordenadas cartesianas de la siguiente forma: [math] \vec{v}(x,y) = \Big(-(a\,ρ+\frac{b}{ρ})\,\frac{x}{ρ},\; -(a\,ρ+\frac{b}{ρ})\,\frac{y}{ρ}\Big) = \Big(-(a\,ρ+\frac{b}{ρ})\cos\theta,\; -(a\,ρ+\frac{b}{ρ})\sin\theta\Big) [/math]

• Campo v en coordenadas cilíndricas

En coordenadas cilíndricas [math](ρ,\theta,z)[/math], el campo [math]\vec{v} = \vec{k} \times \vec{u}[/math] se expresa de forma muy sencilla.

Sabemos que:

[math]\vec{u}(ρ,\theta) = f(ρ)\,\vec{e_\theta}[/math], con [math]f(ρ) = a\,ρ + \tfrac{b}{ρ}[/math].

El producto cruzado con [math]\vec{k} = \vec{e_z}[/math] da:

[math]\vec{k}\times\vec{e_\theta} = -\,\vec{e_ρ}[/math].

Por tanto:

[math]\vec{v}(ρ,\theta) = -\,f(ρ)\,\vec{e_ρ} = -\Big(a\,ρ + \tfrac{b}{ρ}\Big)\,\vec{e_ρ}[/math]

[math]\boxed{\;\vec{v}(ρ,\theta) = -\,f(ρ)\,\vec{e_ρ} = -\Big(a\,ρ + \tfrac{b}{ρ}\Big)\,\vec{e_ρ}\;}[/math]

---

• Componentes en la base cilíndrica

[math]v_ρ = -\Big(a\,ρ + \tfrac{b}{ρ}\Big), \quad v_\theta = 0, \quad v_z = 0[/math]

• Función de corriente [math]\psi(ρ)[/math]

Queremos calcular la función de corriente [math]\psi(ρ)[/math] tal que:

[math]\nabla \psi = \vec{v}[/math]

---

1. Observación previa

El campo [math]\vec{v} = \vec{k} \times \vec{u} = -f(ρ)\,\vec{e_ρ}[/math] es radial y apunta hacia el centro. Además:

- [math]\vec{v}[/math] es **irrotacional** (porque [math]\vec{u}[/math] tiene divergencia nula). - Existe un potencial escalar [math]\psi[/math] tal que [math]\nabla \psi = \vec{v}[/math]. Este potencial se conoce como **función de corriente** de [math]\vec{u}[/math].

---

2. Gradiente en coordenadas cilíndricas

En coordenadas cilíndricas [math](ρ,\theta,z)[/math]:

[math]\nabla \psi = \frac{\partial \psi}{\partial r}\,\vec{e_r} + \frac{1}{ρ}\frac{\partial \psi}{\partial \theta}\,\vec{e_\theta} + \frac{\partial \psi}{\partial z}\,\vec{e_z}[/math]

---

3. Igualación con el campo [math]\vec{v}[/math]

Sabemos que:

[math]\vec{v}(ρ,\theta) = -\Big(a\,ρ + \tfrac{b}{ρ}\Big)\,\vec{e_ρ}[/math]

Por tanto:

[math]\frac{\partial \psi}{\partial ρ} = -\Big(a\,ρ + \tfrac{b}{ρ}\Big), \quad \frac{\partial \psi}{\partial \theta} = 0, \quad \frac{\partial \psi}{\partial z} = 0[/math]

---

\subsection*{4. Integración respecto a [math]ρ\ltmath\gt} \ltmath\gt\psi(ρ) = -\int \Big(a\,ρ + \tfrac{b}{ρ}\Big)\,dρ[/math]

[math]\psi(ρ) = -\Big(\tfrac{a}{2}ρ^2 + b\ln ρ\Big) + C[/math]

---

\subsection*{5. Resultado final}

[math]\boxed{\psi(ρ) = -\tfrac{a}{2}ρ^2 - b\ln ρ + C}[/math]

---

\subsection*{6. Interpretación}

- Como [math]\psi[/math] depende únicamente de [math]ρ[/math], las curvas de nivel [math]\psi = \text{cte}[/math] son **circunferencias concéntricas** [math]ρ = \text{constante}[/math]. - Estas curvas son las **líneas de corriente del campo $\vec{u}$**, ya que se cumple:

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

• Irrotacionalidad de v

En coordenadas cilíndricas, para un campo puramente radial [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 ρ}(r 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.

• Potencial escalar

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]

• líneas 𝜓 = cte

Como $\psi$ 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].

Como $\psi$ depende únicamente de [math]ρ[/math], las líneas [math]\psi = \text{cte}[/math] son circunferencias concéntricas. Basta con elegir varios valores de [math]ρ[/math] en el intervalo [math][1,2][/math] y dibujar los círculos correspondientes.

center

Líneas de corriente

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

Esto se puede comprobar que siempre será negativo, pues ambos términos están precedidos por un signo negativo para [math]r \gt 0[/math]. Así que, es decreciente en este intervalo.

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 ρ:

Módulo de la velocidad

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 es puramente azimutal 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 vorticidad

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]

∇×u

%% 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

Campo de Temperaturas y curvas de nivel

%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 en [math]\theta = 0[/math] o [math]\theta = \pi[/math].

Por tanto, el máximo de [math]T[/math] en el anillo se da en: [math] \rho = 2,\qquad \theta = 0 \text{ o } \pi. [/math]

7.3 Máximo de la temperatura en el anillo [math]1 \le \rho \le 2[/math]

Trabajaremos en la sección transversal y dentro del fluido (anillo [math]1 \le \rho \le 2[/math]).

La temperatura está dada por: [math] T(\rho,\theta) = \log\!\big(1+\rho^2\big)\,\cos^2\theta. [/math]

7.4 Comportamiento de los factores

1) Para [math]\rho \ge 0[/math], 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, en [math]\theta = 0,\,\pi,\,2\pi,\,\dots[/math]

7.5 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 en [math]\theta = 0[/math] o [math]\theta = \pi[/math].

Por tanto, el 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

Para [math] T(r,\theta) = \log(1+r^2)\,\cos^2\theta, [/math] calculamos sus derivadas parciales:

[math] \frac{\partial T}{\partial r} = \frac{2r}{1+r^2}\,\cos^2\theta, \qquad \frac{\partial T}{\partial \theta} = \log(1+r^2)\,\frac{d}{d\theta}(\cos^2\theta) = -\,\log(1+r^2)\,\sin(2\theta). [/math]

8.1 En base cilíndrica

El gradiente se expresa como: [math] \nabla T = \frac{\partial T}{\partial r}\,\mathbf{e}_r + \frac{1}{r}\frac{\partial T}{\partial \theta}\,\mathbf{e}_\theta = T_r\,\mathbf{e}_r + T_\theta\,\mathbf{e}_\theta, [/math] donde: [math] T_r = \frac{2r}{1+r^2}\,\cos^2\theta, \qquad T_\theta = -\,\frac{1}{r}\,\log(1+r^2)\,\sin(2\theta). [/math]

8.2 Representación Gráfica

Gradiente de temperatura

Curvas de nivel y gradiente

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 Interpretación física: Cálculo del caudal

La sección [math]\theta = 0[/math] es una línea radial desde [math]r = 1[/math] hasta [math]r = 2[/math]. El fluido se mueve azimutalmente ([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_{r=1}^{r=2} u_\theta(r)\, h \, dr, [/math] porque el área diferencial es [math]dA = h \, dr[/math] (una tira de ancho [math]dr[/math] y profundidad [math]h[/math]).

9.2 Perfil de velocidad

[math] u_\theta(r) = f(r) = a r + \frac{b}{r}, \qquad a = -\frac{5}{3},\; b = \frac{8}{3}. [/math]

Entonces: [math] Q = \int_{1}^{2} \big(a r + \frac{b}{r}\big)\, dr = \Big[\frac{a}{2}r^2 + b \ln r\Big]_{1}^{2}. [/math]

9.3 Sustituimos valores

[math] \frac{a}{2}r^2 = \frac{-5/3}{2}r^2 = -\frac{5}{6}r^2,\qquad b \ln r = \frac{8}{3}\ln r. [/math]

Evaluamos: [math] Q = \Big(-\frac{5}{6}(4) + \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]

Como [math]\ln 2 \approx 0.693[/math]: [math] \frac{8}{3}\ln 2 \approx \frac{8}{3}\cdot 0.693 \approx 1.848. [/math]

Entonces: [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).

El módulo del caudal es: [math] |Q| \approx 0.65 \;\text|Q| \approx 0.65 \;\text{m}^3/\text{s}. ===Representación Gráfica===[/math]