Flujo de Couette entre dos tubos cocéntricos. Grupo 14

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos. Grupo 14
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Oscar García Caballero Daniel García-Alcaide Presa Eduardo Juarranz del Valle Juan Holgado Sánchez Jose Antonio Calvo de las Heras
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 fundamental en la mecánica de fluidos que describe el comportamiento de un fluido entre dos superficies paralelas, una de las cuales se mueve en relación con la otra mientras el fluido permanece confinado en el espacio intermedio. Este flujo, caracterizado por su sencillez y su amplia aplicabilidad, ha sido estudiado extensamente debido a su relevancia en la modelización de sistemas mecánicos.

En este artículo se va a proceder al análisis del flujo de un fluido incomprensible a través de dos cilíndros concéntricos de manera que el cilindro exterior se mueve con velocidad angular constante en sentido horario, mientras que el cilindro interior está fijo.

1 Mallado de la sección transversal

Primeramente vamos a representar una sección transversal del fluido respecto a los parámetros [math] \rho=2 [/math] para el tubo exterior, [math] \rho=1 [/math] para el tubo interior y [math] \theta ∈ [0,2\pi][/math] . Seccionamos ambos tubos según el plano [math]x_3=0[/math], representamos los ejes mediante los [math] (x,y) ∈ [-4,4] × [-4,4][/math].

Para poder representar la sección seguimos el siguiente código en MatLab:

% MALLADO DE LA SECCIÓN TRANSVERSAL
rho=1:0.08:2;                                % INTERVALO DE RHO [1,2]
theta=0:0.08:2*pi;                           % 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([-4,4,-4,4])                            % EJES DEL DIBUJO  
title('MALLADO DE LA SECCION TRANSVERSAL')   % TITULO DE LA GÁFICA
hold off

2 Ecuación Navier-Stokes

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. Además, el campo de velocidades tiene que cumplir la ecuación de Navier-Stockes estacionaria definida por la siguiente expresión:

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

Nota: El valor µ representa el coeficiente de viscosidad del fluido.

Analizando la expresión, aseguramos que al tratarse de una presión constante su gradiente es nulo ( [math]∇p = 0[/math] ) y además despreciando el término convectivo (primer término de la fórmula) finalmente obtenemos que:

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

Nota: [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 unicamente de las condiciones de frontera en el recinto que se indique.

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

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.1.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.1.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.1.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.1.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.2 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.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.2.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.2.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

Para el cálculo del campo de velocidades suponemos que ω = 1 y μ = 1, y por lo tanto las constantes serán:

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

De tal forma 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]

A continuación, representamos este campo mediante una gráfica de mathlab:

% 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,-3,3]) 
quiver(x,y,X,Y);                          % DIBUJO CAMPO VECTORIAL
hold off

En la imagen del campo de velocidades podemos apreciar lo siguiente: las flechas indican que las partículas del fluido se mueven en trayectorias cerradas; el tamaño de las flechas puede ser proporcional a la magnitud de la velocidad en cada punto, las flechas son más grandes en regiones exteriores lo que indica que la velocidad aumenta al alejarse del centro y el centro muestra una zona con flechas pequeñas o ausentes que indica una baja velocidad, al alejarnos del eje la velocidad aumenta.

4 Líneas de Corriente del campo

Las líneas de corriente son tangentes al campo [math] \vec u [/math] en cada punto, de modo que vamos a buscar el campo [math] \vec v [/math] que es perpendicular al campo [math] \vec u [/math].

[math] \vec{v} = {e_z}\times\vec{u} \Longrightarrow \vec{v} = \frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) (\vec{e_z}\times\vec{e_\theta}) \Longrightarrow \vec{v} = \frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) (-\vec{e_\rho}) \Longrightarrow \vec{v} = -\frac{4}{3}\left ( \rho - \frac{1}{\rho} \right )\cdot\vec{e_\rho} [/math]

4.1 Irrotacionalidad de [math] \vec v [/math]

Vamos a comprobar si el campo [math] \vec v [/math] es irrotacional, el campo seguramente lo sea ya que la divergencia de [math] \vec u [/math] es nula. Esto se comprueba al calcular el rotacional del campo [math] \vec v [/math].

[math] \nabla\times\vec{v} = \begin{vmatrix} \vec{e_\rho} & \rho\vec{e_\theta} & \vec{e_z} \\ \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\ -\frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) & 0 & 0 \end{vmatrix} \Longrightarrow \nabla\times\vec{v} = \vec{0} [/math]

Efectivamente el campo [math] \vec v [/math] es irrotacional.

4.2 Cálculo de la función de corriente de [math] \vec u [/math] y de [math] \psi [/math]

Sabemos que el vector [math] \vec v [/math] es:

[math] \vec{v} = \nabla\psi \Longrightarrow -\frac{4}{3}\left ( \rho - \frac{1}{\rho} \right )\cdot\vec{e}_\rho = \frac{\partial\psi}{\partial\rho}\vec{e}_\rho + \frac{1}{\rho}\frac{\partial\psi}{\partial\theta}\vec{e}_\theta + \frac{\partial\psi}{\partial z}\vec{e}_z [/math]

De esta igualdad despejamos [math] \psi [/math]:

[math] \frac{\partial\psi}{\partial\rho} = -\frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) \Longrightarrow \psi = \int -\frac{4}{3}\left ( \rho - \frac{1}{\rho} \right ) \partial\rho =-\frac{4}{3}\int-\frac{1}\rho {\partial\rho} - \frac{4}{3}\int\rho {\partial\rho}=\frac{4}{3}[ln(\rho)-\frac{\rho^2}{2}] [/math]

4.3 Representación de las líneas [math] \psi=cte [/math]

Las líneas del potencial se van a representar con el siguiente código:

% LINEAS DE CORRIENTE
h=90;                             % NÚMERO DE PARTES DEL INTERVALO
rho=linspace(1,2,h);              % PARAMERTO RHO [1,2]
theta=linspace(0,2*pi,h);         % PARAMETRO THETA [0,2*pi]
[RHO,THETA]=meshgrid(rho,theta);  % MATRIZ DE PARAMETROS
fv=((4*log(RHO)/3)-((4*RHO)/3));  % CAMPO ESCALAR
X=RHO.*cos(THETA);                % PARAMETRIZACIÓN
Y=RHO.*sin(THETA);
hold on
grid on
title('LINEAS DE CORRIENTE')
axis([-3,3,-3,3])         % EJES
view(2)                     
contour(X,Y,fv,20);colorbar       % LÍNEAS DE NIVEL
colormap("default")
hold off

En cuanto a la gráfica, podemos interpretar que el flujo es estacionario, dado que las líneas de corriente son cerradas y regulares; el color de las líneas corresponde a un campo escalar donde los colores más cálidos indican mayor magnitud relativa, y los fríos menor magnitud; la región interior podría representar una zona de baja velocidad y el gradiente de colores sugiere que la magnitud de la velocidad disminuye hacia el centro.

5 Velocidad máxima del fluido

Se obtiene la velocidad máxima del fluido con el módulo del campo vectorial:

[math]\ \left | \vec u (\rho ) \right | = \left | \left ( \frac{4}{3}\rho w - \frac{4w}{3\rho } \right ) \vec e_\theta \right | = \frac{4}{3} w \left ( \rho -\frac{1}{\rho } \right ) [/math]

A continuación se deriva e iguala a 0, teniendo en cuenta que la velocidad angular es [math]\ w=1 [/math]

[math]\ \frac{\partial \left | \vec u (\rho ) \right | }{\partial \rho } = \frac{4}{3} \left ( 1 + \frac{1}{\rho^2} \right ) = 0 [/math]

Esta solución sería: [math]\ \sqrt{-1} [/math]; por lo que no es Real

Por ello, hay que analizarlo en los extremos del intervalo, nuestros extremos están en [math]\ \rho \in (1,2) [/math]

El resultado del campo en estos dos extremos es:

[math]\ \left\{\begin{matrix} \left | \vec u (1) \right | = \left | \left ( \frac{4}{3} \rho w - \frac{4w}{3\rho } \right ) \vec e_\theta \right | = \frac{4}{3} \left ( 1 -\frac{1}{1 } \right ) = 0 \\ \left | \vec u (2) \right | = \left | \left ( \frac{4}{3} \rho w - \frac{4w}{3\rho } \right ) \vec e_\theta \right | = \frac{4}{3} \left ( 2 -\frac{1}{2 } \right ) = 1 \end{matrix}\right. [/math]

Por lo tanto se comprueba [math]\ \rho = 2 [/math] y su valor es de 1

derecha

clear, close all
rho=1:0.1:2;
theta=0:0.1:2*pi; 
[RHO,THETA]=meshgrid(rho,theta); 
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
Modulo=2/3.*RHO-2/3./RHO; 
mesh(x,y,Modulo)
axis([-3,3,-3,3]);
surf(x,y,Modulo);
view(2);
colorbar;
title('Variación del módulo de la velocidad en función de RHO');
xlabel('Eje X1');
ylabel('Eje X2');

En la anterior figura 4 se demuestra el comportamiento del módulo. Se puede observar, como ya se ha explicado, que el máximo obtenido se encuentra en [math]\ \rho = 2 [/math] y su valor es de 1. Además el comportamiento del fluido muestra un decrecimiento de velocidad hasta llegar a [math]\ \rho = 1 [/math] dónde la velocidad es nula.

6 Rotacional de [math] \vec{u} [/math]

Se calcula el rotacional de [math] \vec u [/math] en coordenadas cilíndricas:

[math]\ \vec u(\rho ,\theta ) = f(\rho )\cdot \vec e_\theta [/math] [math]\ \bigtriangledown \times \vec u = \frac{1}{\rho} \cdot \left|\begin{matrix} \vec e_\rho & \rho\cdot \vec e_\theta & \vec e_z \\ \frac{\partial}{\partial \rho } & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ 0 & \rho\cdot f(\rho ) & 0 \end{matrix}\right| = \frac{1}{\rho} \cdot \left (\left ( f(\rho) + \rho\cdot f'(\rho ) \right ) \vec e_z \right ) = \left (\frac{f(\rho )}{\rho} + f'(\rho ) \right )\vec e_z [/math]

Como se demuestra en apartados anteriores [math]\ f(\rho ) [/math] es [math]\ f(\rho ) = \frac{4}{3} w\rho - \frac{4w}{3\rho } [/math] con la velocidad angular [math]\ w=1 [/math]

Por lo tanto, el campo vectorial es [math]\ \vec u(\rho ,\theta ) = f(\rho )\cdot \vec e_\theta = \frac{4}{3}\left ( \rho - \frac{1}{\rho }\right ) \vec e_\theta[/math]

Quedando el rotacional:

[math]\ \bigtriangledown \times \vec u = \left (\frac{f(\rho )}{\rho} + f'(\rho ) \right )\vec e_z = \left ( \frac{ \frac{4}{3} \left ( \rho - \frac{1}{\rho } \right ) }{\rho} + \frac{4}{3} \left (1 + \frac{1}{\rho ^2} \right ) \right ) \vec e_z = \frac{4}{3} \left ( 1- \frac{1}{\rho ^2} + 1 + \frac{1}{\rho ^2} \right )\vec e_z = \frac{8}{3} \vec e_z [/math]

Siendo la distancia del vector [math]\ \left | \bigtriangledown \times \vec u \right | = \sqrt{\left ( \frac{8}{3} \right )^2} = \frac{8}{3} [/math]

Este es el código con el que se representa el campo [math]\ \left | \bigtriangledown \times \vec u \right | [/math] :

derecha

%MÓDULO ROTACIONAL

 rho=1:0.08:2;
 theta=0:0.08:2*pi;
  
  %se genera una retícula rectangular donde se representan rho y theta
  [RHO,THETA]=meshgrid(rho,theta);
    
  %Pasamos a cartesianas
  x=RHO.*cos(THETA);
  y=RHO.*sin(THETA);
  z=ones(size(y)).*(8/3);
 
 %Representación del módulo del rotacional
 subplot(1,2,1)
 surf(x,y,z);
 axis([-3,3,-3,3]);
 title('Módulo del Rotacional');
 xlabel('Eje X1');
ylabel('Eje X2');
zlabel('Eje X3');



%Representación del rotacional
 rho=1:0.5:2;
 theta=0:0.5:2*pi;
 [RHO,THETA]=meshgrid(rho,theta);
 x=RHO.*cos(THETA);
 y=RHO.*sin(THETA);
 z=zeros(size(x));
 xx=zeros(size(x));
 yy=zeros(size(x));
zz=ones(size(y)).*4/3;

subplot(1,2,2)
quiver3(x,y,z,xx,yy,zz);
axis([-3,3,-3,3]);
title('Representación del rotacional');
xlabel('Eje X1');
ylabel('Eje X2');
zlabel('Eje X3');

El módulo del rotacional es constante e igual a [math] \frac{8}{3} [/math] para cualquier [math] \rho [/math] y [math] \theta [/math]. Además en la gráfica se puede ver que el rotacional es igual en todos los puntos en dirección [math] \vec e_z [/math].

7 Representación del Campo de Temperaturas

Sabiendo que la temperatura del fluido viene dada por la expresión [math] T(\rho,\theta) = log(1+ \rho^2)\cos^2 \theta\ e^{−\left ( \rho − \frac{5}{2} \right ) ^2} [/math] se dibuja el campo de temperaturas y las curvas de nivel:

7.1 Representación del campo y curvas de nivel

Mediante el siguiente código de matlab mostramos la interpretación gráfica del campo escalar que representa la temperatura del fluido. Se presentan dos vistas, una en tres dimensiones y en dos dimensiones. Las gráficas tienen una escala de colores que representa el valor de la temperatura en cada punto de la región estudiada.

%MALLADO DE LA SECCIÓN TRANSVERSAL
rho=1:0.08:2;
theta=0:0.08:2*pi;
[RHO,THETA]=meshgrid(rho,theta);
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
z=log(1+RHO.^2).*cos(THETA).*exp(-(RHO-(5/2)).^2);
%Representación del campo en 3D
subplot(1,3,1)
surf(x,y,z);
axis([-3,3,-3,3]);
title('Representación de la temperatura en 3D')
xlabel('Eje X1');
ylabel('Eje X2');
zlabel('Eje X3');
%REPRESENTACIÓN DEL CAMPO EN 2D
colorbar;
subplot(1,3,2)
pcolor(x,y,z);
axis([-3,3,-3,3]);
title('Representación de la temperatura en 2D')
xlabel('Eje X1');
ylabel('Eje X2');
colorbar;

%REPRESENTACIÓN DE LAS LÍNEAS DE NIVEL
subplot(1,3,3);
pcolor(x,y,z);
%axis([-3,3,-3,3]);
contour(x,y,z,"LineWidth",2)
axis([-3,3,-3,3]);
colorbar;
title('Representación de las lineas de nivel de la temperatura')
xlabel('Eje X1');
ylabel('Eje X2');

8 Gradiente de la temperatura

8.1 Cálculo del gradiente de T

El gradiente de la temperatura se define mediante la siguiente fórmula:

[math] \nabla T(\rho,\theta) = \frac{\partial T}{\partial\rho}\vec{e}_\rho + \frac{1}{ρ}\frac{\partial T}{\partial\theta}\vec{e}_\theta + \frac{\partial T}{\partial z}\vec{e}_z [/math]

De manera que:

[math]\ \bigtriangledown T(\rho ,\theta ,z) = \cos^{2}\theta \ e^{-\left(\rho-\frac{5}{2}\right)^{2}} \left (\frac{( 2\rho +(2\rho ^3 -5\rho^2+2\rho-5)\left (\log(1+\rho^2) \right) }{(1+\rho^2)} \right )\vec e_\rho - \frac{1}{\rho}\log(1+\rho^2) 2\sin \theta \cos \theta \ e^{-\left(\rho-\frac{5}{2}\right)^{2}} \vec e_\theta [/math]

clear
clc
rho=1:0.13:2;
theta=0:0.13:2*pi;
[RHO,THETA]=meshgrid(rho,theta);
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
z=log(1+RHO.^2).*cos(THETA).^2.*exp(-(RHO-(5/2)).^2);
[DX,DY]=gradient(z);
hold on;
quiver(x,y,DX,DY);
axis([-3,3,-3,3]);
title('Gradiente de la temperatura');
xlabel('Eje X1');
ylabel('Eje X2');
hold off

8.2 Demostración gráfica de la ortogonalidad del gradiente de la temperatura y las curvas de nivel

Para ello, se ha utilizado el siguiente código de Matlab, donde se pueden ver las líneas de nivel en colores y el campo gradiente.

clear
clc
rho=1:0.13:2;
theta=0:0.13:2*pi;
[RHO,THETA]=meshgrid(rho,theta);
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
z=log(1+RHO.^2).*cos(THETA).^2.*exp(-(RHO-(5/2)).^2);
[DX,DY]=gradient(z);
figure;
quiver(x,y,DX,DY);
axis([-3,3,-3,3]);
title('Comprobación ortogonalidad del gradiente');
xlabel('Eje X1');
ylabel('Eje X2');
%curvas de nivel
hold on
contour(x,y,z,"LineWidth",2)
hold off

El gradiente, por definición, es la dirección de máximo crecimiento de la función. Las líneas de nivel de la temperatura representan donde la función es constante. Por lo tanto, si recorres una línea de nivel el valor de la función no cambia, por lo que el movimiento está siendo ortogonal a la dirección en la que más cambia ,que es el gradiente.

9 Caudal circulante en la sección longitudinal

Partiendo del campo de velocidades hallado anteriormente,

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

Introducimos la siguiente fórmula para el cálculo del caudal del campo a través de una superficie [math]\ S [/math] (la sección longitudinal):

[math]\ \int_{S}^{} \vec u \cdot d\vec S = \iint_{D}^{} \vec u \left ( \Phi (u,v) \right )\cdot \left ( \Phi _u\times \Phi _v \right ) du dv [/math]

Como se puede ver en la figura, la superficie [math]\ S [/math] está compuesta por la suma de dos superficies rectangulares [math]\ S_1 [/math] y [math]\ S_2 [/math]. En primer lugar, debemos orientar dicha superficie y parametrizarla consecuentemente. La orientación decidida es aquella que coincide con el sentido de crecimiento del campo de velocidades.

La parametrización de [math]\ S_1 [/math] sería:

[math]\ S_1 [/math]: [math]\ \Phi (u,v) = (\rho (u,v),\theta (u,v),z(u,v)) = (v,\pi/2,u) [/math] con [math]\ \left\{\begin{matrix} u\in (0,1) \\ v\in (1,2) \end{matrix}\right. [/math]

Procedemos al cálculo del vector normal a la superficie:

[math]\ \left.\begin{matrix} \Phi _u = \vec e_\rho \\ \Phi _v = \vec e_z \end{matrix}\right\} [/math] [math]\ \rightarrow [/math] [math]\ \Phi _u\times \Phi _v = \vec e_\theta [/math]

Esta orientación no sirve ya que va al contrario a la que se tomó, por lo que en el cálculo de la integral habrá que poner un menos para cambiar el sentido.

Y la parametrización de [math]\ S_2 [/math] sería:

[math]\ S_2 [/math]: [math]\ \Phi (u,v) = (\rho (u,v),\theta (u,v),z(u,v)) = (v,3\pi/2,u) [/math] con [math]\ \left\{\begin{matrix} u\in (0,1) \\ v\in (1,2) \end{matrix}\right. [/math]

Procedemos al cálculo del vector normal a la superficie:

[math]\ \left.\begin{matrix} \Phi _u = \vec e_\rho \\ \Phi _v = \vec e_z \end{matrix}\right\} [/math] [math]\ \rightarrow [/math] [math]\ \Phi _u\times \Phi _v = \vec e_\theta [/math]

Finalmente, el caudal se calcula como la integral de la sección longitudinal, que se puede poner como la suma de las integrales de las superficies [math]\ S_1 [/math] y [math]\ S_2 [/math].

[math]\ \int_{S}^{} \vec u \cdot d\vec S = \int_{S_1}^{} \vec u \cdot d\vec S_1 + \int_{S_2}^{} \vec u \cdot d\vec S_2 = [/math]

[math]\ = \int_{1}^{2}\int_{0}^{1} \vec u (\rho,\theta,z)\cdot (\vec e_\theta)dudv + \int_{1}^{2}\int_{0}^{1} \vec u (\rho,\theta,z)\cdot (\vec e_\theta)dudv = [/math]

[math]\ = \frac{4}{3} \int_{1}^{2}\int_{0}^{1}\left ( v-\frac{1}{v} \right )dudv + \frac{4}{3} \int_{1}^{2}\int_{0}^{1}\left ( v-\frac{1}{v} \right )dudv = [/math]

[math]\ = \frac{8}{3} \int_{1}^{2}\int_{0}^{1}\left ( v-\frac{1}{v} \right )dudv = [/math]

[math]\ = \frac{8}{3} \int_{1}^{2} \left [vu-\frac{1}{v}u\right]^1_0 dv = \frac{8}{3}\int_{1}^{2}\left (v-\frac{1}{v}\right)dv = \frac{8}{3}\left [(\frac{v^2}{2})-lnv\right]^2_1 = 4-\frac{8}{3}ln2 = 2.15 [/math]

Este resultado es coherente, ya que el valor obtenido es mayor que cero, lo cual concuerda con la hipótesis inicial de orientar la superficie en el sentido del crecimiento de la velocidad.

10 Ejemplos prácticos

A continuación podemos observar algunos ejemplos de este flujo de couette entre tubos cocéntricos:

10.1 Medición de viscosidad en reómetros

Un reómetro con geometría de cilindros concéntricos se utiliza para medir la viscosidad de fluidos como el asfalto modificado, utilizado en carreteras.

10.2 Sistemas de lubricación

El lubricante entre cojinetes cilíndricos concéntricos experimenta un flujo de Couette. Este diseño reduce el desgaste por fricción y mejora la eficiencia mecánica del sistema.

10.3 Diseño de sistemas de drenaje

En túneles largos, el agua se transporta a través de sistemas de drenaje diseñados con tubos concéntricos. El flujo de Couette modela el comportamiento del agua en espacios anulares cuando hay un gradiente de velocidad entre el flujo en la pared del tubo interno y el externo. Este análisis asegura un drenaje eficiente y evita acumulación de agua.

10.4 Revestimiento protectores en tuberías

Cuando se aplican recubrimientos líquidos o plásticos en la superficie interna de tuberías, el flujo generado en el espacio entre el cilindro aplicador y la tubería sigue patrones similares al flujo de Couette. Esto asegura un espesor uniforme del revestimiento.

11 Bibliografía

Wikipedia. (s.f.). Flujo de Couette: https://es.wikipedia.org/wiki/Flujo_de_Couette#:~:text=En%20Din%C3%A1mica%20de%20fluidos%2C%20Flujo,relativo%20con%20respecto%20al%20otro.

StudySmarter. (s.f.). Flujo de Couette: https://www.studysmarter.es/resumenes/ingenieria/mecanica-de-fluidos-en-ingenieria/flujo-de-couette/

Wikipedia. (s.f.). Flujo Taylor-Couette: https://es.wikipedia.org/wiki/Flujo_Taylor-Couette

Apuntes proporcionados en ETSI Caminos, Canales y Puertos, de la Universidad Politécnica de Madrid