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

Representación de la sección transversal 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]

Debido a que el campo de velocidades está proporcionado en coordenadas cilíndricas, resulta más sencillo calcularlo con dichas coordenadas y por lo tanto la ecuación pasaría a ser la siguiente:

[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

Campo de velocidades

=
clear all
clc
hold on

radiodemalla = 0.1:0.1:2;           % radios desde 0.1 hasta 2
angulodemalla = 0:0.1:2*pi;
[R, A] = meshgrid(radiodemalla, angulodemalla);

%paso a cartesianas
Xpos = R .* cos(A);
Ypos = R .* sin(A);


U = zeros(size(R));
V = zeros(size(R));

%región antihorario
region1 = (R < 1);

% Componente tangencial (dirección antihoraria: +θ̂)
% θ̂ = (-sinθ, cosθ)
U(region1) = -sin(A(region1));  
V(region1) =  cos(A(region1));

% Magnitud opcional según (1 - r):
magnitud1 = (1 - R(region1));
U(region1) = U(region1) .* magnitud1;
V(region1) = V(region1) .* magnitud1;

%=========================================================
%   CAMPO CLOCKWISE (HORARIO) PARA 1 < r < 2
%=========================================================
region2 = (R > 1 & R <= 2);

% Dirección horaria = -θ̂ = (sinθ, -cosθ)
U(region2) =  sin(A(region2));
V(region2) = -cos(A(region2));

% Magnitud opcional según (r - 1):
magnitud2 = (R(region2) - 1);
U(region2) = U(region2) .* magnitud2;
V(region2) = V(region2) .* magnitud2;

%=========================================================
%   Campo nulo exactamente en r = 1
%=========================================================
region3 = (abs(R - 1) < 1e-6);
U(region3) = 0;
V(region3) = 0;

%--- Gráfico ---
quiver(Xpos, Ypos, U, V, 'AutoScale','on');

axis equal
axis([-2.5 2.5 -2.5 2.5]);
title('Campo antihorario para r<1 y horario para 1<r<2');
hold off