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

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

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

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

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

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]

2.2.1 Gradiente de la divergencia

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

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

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

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

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

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

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

2.2.2 Rotacional del campo de velocidades

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

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

2.2.3 Rotacional del rotacional del campo de velocidades

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

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

2.2.4 Resultado del Laplaciano vectorial

Sustituimos en la ecuación obtenida del laplaciano vectorial.

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

2.3 Resultado de la ecuación de Navier-Stokes

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

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

2.3.1 Ecuación diferencial

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

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

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

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

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

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

2.3.2 Análisis de la solución conocida

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]

2.3.3 Búsqueda de las constantes "a" y "b"

Por último, para hallar el campo de velocidades concreto, tenemos que imponer las condiciones que conocemos, ya que al resolver la ecuación diferencial de segundo orden han aparecido dos constantes desconocidas.

Condición 1: Si [math]ρ=1 \Longrightarrow \vec{u}=\vec{0} \Longrightarrow f(1)\vec{e_θ}=\vec{0} \Longrightarrow (a+b)\vec{e_θ}=\vec{0}\Longrightarrow a+b=0 [/math]

Condición 2: Si [math]ρ=2 \Longrightarrow \vec{u}= \omega \vec{e_z} \Longrightarrow f(2)\vec{e_z}=\omega \vec{e_z} \Longrightarrow (2a+\frac{b}{2})\vec{e_θ}=2\omega\vec{e_θ} \Longrightarrow 2a+\frac{b}{2}=2\omega [/math]

Una vez impuestas estas dos condiciones, resulta un sistema de dos ecuaciones con dos incógnitas donde

[math] a = \frac{4}{3}\omega [/math] [math] b = -\frac{4}{3}\omega[/math]

concluyendo por lo tanto que:

[math]f(ρ) = \frac{4}{3}\omega ρ -\frac{4}{3ρ}\omega[/math]

3 Campo de velocidades

Con el objetivo de calcular el campo de velocidades establecemos que ω = 1 y μ = 1, y en consecuencia las constantes serán de la siguiente forma:

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

Sustituimos en la función:

[math] f(\rho) = a\rho + \frac{b}{\rho} \Longrightarrow f(\rho) = \frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) [/math]

Es por tanto que nuestro campo de velocidades quedaría expresado por:

[math] \vec{u} = f(\rho)\vec{e}_\theta \Longrightarrow \vec{u} = \frac{4}{3}\left ( \rho - \frac{1}{\rho} \right )\cdot\vec{e}_\theta [/math]

En la siguiente figura, se representa el campo mediante una gráfica de Matlab:

Campo de velocidades

% CAMPO DE VELOCIDADES
u=1:0.15:2;                               % INTERVALO DE U
v=0:0.3:2*pi;                            % INTERVALO DE V
[U,V]=meshgrid(u,v);                      % MATRIZ DE PARAMETROS
x=U.*cos(V);                              % PARAMETRIZACIÓN
y=U.*sin(V); 

X=sin(V).*((4/3)*(U-1./U));               % FUNCION
Y=cos(V).*(-(4/3)*(U-1./U));              % FUNCION

hold on
grid on 
title('CAMPO DE VELOCIDADES')
axis([-3,3,-2.5,2.5]) 
quiver(x,y,X,Y); 
                         % DIBUJO CAMPO VECTORIAL
hold off

4 Líneas del campo de corriente

Antes de representar las líneas de corriente, hemos de verificar primero la irrotacionalidad del campo, antes de ello lo que haremos será multiplicar vectorialmente por el vector “k” de la base ortonormal cartesiana (0,0,1).

[math]\vec{v}=\vec{k}\times\vec{u}=\begin{vmatrix} \vec{e}_{\rho }&\rho \cdot \vec{e} _{\theta } & \vec{e}_{z}\\ 0 & \frac{4}{3}\left (\rho-\frac{4}{3\rho}\right ) & 0 \\ 0 & 0 & 1\end{vmatrix}=\frac{4}{3}\left (\rho-\frac{4}{3\rho}\vec{e}_\rho\right )[/math]

4.1 Comprobación de irrotacionalidad

Procedemos al cálculo la irrotacionalidad para posteriormente obtener las líneas de corriente (Función Potencial).

[math]\triangledown\times\vec v= \frac{1}{\rho}\begin{vmatrix} \vec{e}_{\rho }&\rho \cdot \vec{e} _{\theta } & \vec{e}_{z}\\ \frac{\partial }{\partial \rho } & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z} \\ \frac{4}{3}(\rho-\frac{4}{3\rho}) & 0 & 0\end{vmatrix}=\vec 0[/math]

Por lo tanto, es posible afirmar que el campo es irrotacional.

4.2 Calculo del potencial escalar

Ahora con la teoría del potencial vectorial procedemos a calcular la función potencial, que es la línea de corriente: (Hacemos uso de la expresión del gradiente en coordenadas cilíndricas)

[math]\vec ∇v=\frac{\partial \psi}{\partial \rho}\vec e_\rho+\frac{1}{\rho}\frac{\partial \psi}{\partial θ}\vec e_θ+\frac{\partial\psi}{\partial z} = \frac{4}{3}\rho-\frac{4}{3\rho}\vec e_\rho[/math]

Igualamos términos y resolvemos la función potencial.

[math]\frac{4}{3}\rho-\frac{4}{3\rho} = \frac{\partial\psi}{\partial\rho} \rightarrow \partial\psi = \left (\frac{4}{3}\rho-\frac{4}{3\rho} \right )\partial\rho \rightarrow \psi = \int\left (\frac{4}{3}\rho-\frac{4}{3\rho} \right )\partial\rho \rightarrow \psi = \frac{4}{6}\rho^2-\frac{4}{3}ln(\rho)+C [/math]

Tenemos la línea de corriente.

4.3 Representación gráfica

líneas de corriente

El siguiente gráfico representa las líneas de corriente de la función potencial.

clear all 
clc
puntos = 100; % Definir el número de puntos en la malla
 
den_radial = linspace(1,2,puntos); % Generacion de valores de densidad radial 'ro'
Theta = linspace(0, 2*pi, puntos); % Generacion de valores para theta
[RO, THETA] = meshgrid(den_radial, Theta); % Genera una malla en 2D
funcion_psi = (4/3)*RO^2-(4/3)*log(RO); % Línea de corriente
x = RO.*cos(THETA);                 % Convertir coordenadas polares a coordenadas cartesianas
y = RO.*sin(THETA);
 
view(2); % Vista en 2D
 
contour(x,y,funcion_psi, 20); % Gráfico para la funcion en coordenadas cartesianas
colorbar; % Mostrar la barra de colores
 
axis([-3, 3, -2.5, 2.5]);
colormap("cool");

5 Velocidad máxima del fluido

Para poder calcular los puntos del fluido donde su módulo de velocidad es máxima tenemos que calcular el módulo del campo de velocidades. Sabiendo que [math]\vec u(\rho,\theta)=f(\rho)\vec e_\theta[/math] (los campos de velocidades se calcularon en el apartado anterior a este) Entonces:

[math]\vec u(\rho)=(\frac{4}{3}\rho\omega-\frac{4}{3\rho}\omega)\vec e_\theta[/math]

Tras esto calculamos su módulo

[math]|\vec u(\rho)=(\frac{4}{3}\rho\omega-\frac{4}{3\rho}\omega)\vec e_\theta|=\frac{ω}{3}(4\rho-\frac{4}{\rho})[/math]

Ahora para calcular los puntos máximos derivamos el módulo respecto a \rho y se iguala a 0

[math]\frac {(∂|u⃗ (ρ)|)}{∂ρ}=\frac{ω}{3}(4ρ-\frac{4}{ρ})=0[/math]

Si despejamos \rho nos da un número imaginario.

[math]\rho=\sqrt{-1}[/math]