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-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.1.1 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.1.2 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]

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

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

2.1.5 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.1.6 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.

1. Condiciones en la frontera Sabemos que la velocidad del fluido es: u(ρ,θ)=f(ρ) eθ,f(ρ)=aρ+bρ.\mathbf{u}(\rho,\theta) = f(\rho)\,\mathbf{e}_\theta,\quad f(\rho) = a\rho + \frac{b}{\rho}.u(ρ,θ)=f(ρ)eθ​,f(ρ)=aρ+ρb​. En las fronteras:

En el cilindro interior (ρ=1\rho = 1ρ=1):

uθ(1)=f(1)=a(1)+b1=a+b.u_\theta(1) = f(1) = a(1) + \frac{b}{1} = a + b.uθ​(1)=f(1)=a(1)+1b​=a+b. Esta debe coincidir con la velocidad tangencial del cilindro interior: uθ(1)=ωi⋅1=ωi.u_\theta(1) = \omega_i \cdot 1 = \omega_i.uθ​(1)=ωi​⋅1=ωi​.

En el cilindro exterior (ρ=2\rho = 2ρ=2):

uθ(2)=f(2)=a(2)+b2=2a+b2.u_\theta(2) = f(2) = a(2) + \frac{b}{2} = 2a + \frac{b}{2}.uθ​(2)=f(2)=a(2)+2b​=2a+2b​. Esta debe coincidir con la velocidad tangencial del cilindro exterior: uθ(2)=ωe⋅2=2ωe.u_\theta(2) = \omega_e \cdot 2 = 2\omega_e.uθ​(2)=ωe​⋅2=2ωe​.

2. Sistema de ecuaciones {a+b=ωi,2a+b2=2ωe.\begin{cases} a + b = \omega_i, \\ 2a + \dfrac{b}{2} = 2\omega_e. \end{cases}⎩⎨⎧​a+b=ωi​,2a+2b​=2ωe​.​ Multiplicamos la segunda ecuación por 2 para eliminar fracción: 4a+b=4ωe.4a + b = 4\omega_e.4a+b=4ωe​. Restamos la primera: (4a+b)−(a+b)=4ωe−ωi  ⟹  3a=4ωe−ωi.(4a + b) - (a + b) = 4\omega_e - \omega_i \implies 3a = 4\omega_e - \omega_i.(4a+b)−(a+b)=4ωe​−ωi​⟹3a=4ωe​−ωi​. Por tanto: a=4ωe−ωi3.a = \dfrac{4\omega_e - \omega_i}{3}.a=34ωe​−ωi​​. Luego: b=ωi−a=ωi−4ωe−ωi3=3ωi−(4ωe−ωi)3=4ωi−4ωe3.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}.b=ωi​−a=ωi​−34ωe​−ωi​​=33ωi​−(4ωe​−ωi​)​=34ωi​−4ωe​​.

Resultado a=4ωe−ωi3,b=4(ωi−ωe)3.\boxed{ a = \dfrac{4\omega_e - \omega_i}{3}, \quad b = \dfrac{4(\omega_i - \omega_e)}{3}. }a=34ωe​−ωi​​,b=34(ωi​−ωe​)​.​

3. Comprobación de incompresibilidad En coordenadas cilíndricas: ∇⋅u=1ρ∂(ρuρ)∂ρ+1ρ∂uθ∂θ+∂uz∂z.\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}.∇⋅u=ρ1​∂ρ∂(ρuρ​)​+ρ1​∂θ∂uθ​​+∂z∂uz​​. Como: uρ=0,  uz=0,  uθ=f(ρ),  independiente de θ,z,u_\rho = 0,\; u_z = 0,\; u_\theta = f(\rho),\; \text{independiente de }\theta,z,uρ​=0,uz​=0,uθ​=f(ρ),independiente de θ,z, entonces: ∇⋅u=1ρ∂(0)∂ρ+1ρ(0)+0=0.\nabla \cdot \mathbf{u} = \frac{1}{\rho}\frac{\partial (0)}{\partial \rho} + \frac{1}{\rho}(0) + 0 = 0.∇⋅u=ρ1​∂ρ∂(0)​+ρ1​(0)+0=0. Se cumple la condición de incompresibilidad.

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

%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.25,2.25]) 

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]