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

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos cocé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

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

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

3 Campo de velocidades

Para este apartado vamos a suponer que ω = 1 y μ = 1:

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

Sustituimos en la función:

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

De tal forma que nuestro campo quedaría tal que:

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

Ahora representamos este campo en una gráfica de mathlab:

% CAMPO DE VELOCIDADES

u=1:0.08:2;                               % INTERVALO DE U
v=0:0.08: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).*((2/3)*(U-1./U));               % FUNCION
Y=cos(V).*(-(2/3)*(U-1./U));              % FUNCION
hold on
grid on 
title('CAMPO DE VELOCIDADES')
axis([-2.5,2.5,-2.5,2.5]) 
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. 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{2}{3}\left ( \rho - \frac{1}{\rho} \right ) (\vec{e_z}\times\vec{e_\theta}) \Longrightarrow \vec{v} = \frac{2}{3}\left ( \rho - \frac{1}{\rho} \right ) (-\vec{e_\rho}) \Longrightarrow \vec{v} = -\frac{2}{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{2}{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{2}{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{2}{3}\left ( \rho - \frac{1}{\rho} \right ) \partial\rho =-\frac{2}{3}\int-\frac{1}\rho {\partial\rho} - \frac{2}{3}\int\rho {\partial\rho}=\frac{2}{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=((2*log(RHO)/3)-((2*RHO)/3));  % CAMPO ESCALAR
X=RHO.*cos(THETA);                % PARAMETRIZACIÓN
Y=RHO.*sin(THETA);
hold on
grid on
title('LINEAS DE CORRIENTE')
axis([-2.5,2.5,-2.5,2.5])         % 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 dado que las líneas de corriente son cerradas y regulares, el flujo es estacionario; 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 modulo del campo vectorial:

[math]\ \left | \vec u (\rho ) \right | = \left | \left ( \frac{2}{3}\rho w - \frac{2w}{3\rho } \right ) \vec e_\theta \right | = \frac{2}{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{2}{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{2}{3} \rho w - \frac{2w}{3\rho } \right ) \vec e_\theta \right | = \frac{2}{3} \left ( 1 -\frac{1}{1 } \right ) = 0 \\ \left | \vec u (2) \right | = \left | \left ( \frac{2}{3} \rho w - \frac{2w}{3\rho } \right ) \vec e_\theta \right | = \frac{2}{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

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{2}{3} w\rho - \frac{2w}{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{2}{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{2}{3} \left ( \rho - \frac{1}{\rho } \right ) }{\rho} + \frac{2}{3} \left (1 + \frac{1}{\rho ^2} \right ) \right ) \vec e_z = \frac{2}{3} \left ( 1- \frac{1}{\rho ^2} + 1 + \frac{1}{\rho ^2} \right )\vec e_z = \frac{4}{3} \vec e_z [/math]

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

Con el siguiente código matlab se muestra el campo [math]\ \left \bigtriangledown \times \vec u \right [/math] :

clear, close all


rho=1:0.33:2;
theta=0:0.33: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;

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

El módulo del rotacional es constante e igual a [math] \frac{4}{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.

clear
clc
%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
figure (1)
surf(x,y,z);
shading flat
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;
figure (2);
pcolor(x,y,z);
shading flat 
axis([-3,3,-3,3]);
title('Representación de la temperatura en 2D');
xlabel('Eje X1');
ylabel('Eje X2');
colorbar;

Para las curvas de nivel se utiliza la misma herramienta

clear
clc
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);
contour(x,y,z,15);
axis([-3,3,-3,3]);
title('Representación líneas de nivel');
xlabel('Eje X1');
ylabel('Eje X2');
colorbar;

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 - ln(1+\rho^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).*exp(-(RHO-(5/2)).^2);
[DX,DY]=gradient(z);
figure;
quiver(x,y,DX,DY);
axis([-3,3,-3,3]);
title('Gradiente de la temperatura');
xlabel('Eje X1');
ylabel('Eje X2');

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).*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

9 Caudal circulante en la sección longitudinal

Partiendo del campo de velocidades hallado anteriormente,

[math]\ \vec u(\rho ,\theta ) = \frac{2}{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,0,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 tomo, por lo que en el calculo 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,2\pi,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,\tetha,z)\cdot (\vec e_\theta)dudv = [/math]

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

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

[math]\ = \frac{4}{3} \int_{0}^{1} \left [lnu-\frac{u^2}{2}\right]^2_1 dv = \frac{4}{3}\int_{0}^{1}\left (ln2-ln1-(\frac{4}{2}-\frac{1}{2})\right)dv= \frac{4}{3}\int_{0}^{1}\left (ln2-\frac{3}{2}\right)dv = \frac{4}{3}\left [ln2v-(\frac{3v}{2})\right]^1_0 = \frac{4}{3} (ln2-\frac{3}{2}) = \frac{4}{3}ln2-2 [/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.