Flujo de Couette entre dos tubos concéntricos GRUPO 50

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos (Grupo 50)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Brisa Mora, Alejandro Morales, Nicolás López, Adrián Muñoz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Póster del artículo

Mediante el siguiente enlace podrá acceder al póster del siguiente artículo: Descargar póster

2 Introducción

En este proyecto 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 [math]\vec{\mathbf{\omega_e}}[/math] en sentido horario mientras que el interior se mueve con velocidad angular [math]\vec{\mathbf{\omega_i}}[/math] en sentido contrario. Si suponemos que ambos cilindros tienen su eje en [math]𝑂𝑋_3[/math] y pintamos la sección transversal [math](𝑥_{3} = 0)[/math], el cilindro exterior queda proyectado sobre la la circunferencia 𝜌 = 2 y el interior sobre la circunferencia 𝜌 =1. Trabajaremos en coordenadas cilíndricas

3 Representación de la sección trasversal

Como punto de comienzo, vamos a realizar una sección trasversal de los tubos que represente los tubos que representen los puntos ocupado por el fluido. Para ello hemos dibujado un mallado en Matlab que se ve de la siguiente forma:

Mallado Flujo de Couette entre dos tubos concéntricos

Código Matlab:

rho=1:0.1:2;
theta=linspace(0,2*pi,100);
[RHO,THETA]=meshgrid(rho,theta);
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
hold on;
grid on
axis equal;          
axis([-3 3 -3 3]);
plot(x,y,'b')
plot(x',y','b')
rin=1;       
rex=2;     
th=linspace(0, 2*pi, 400);
plot(rin*cos(th),rin*sin(th),'k','LineWidth',2);
plot(rex*cos(th),rex*sin(th),'k','LineWidth',2);
title('Mallado del flujo de Couette entre dos tubos concéntricos')
hold off;

4 Calculo de las velocidades

4.1 Definición del campo de velocidades

Se sabe que la velocidad de las partículas viene dada por [math]\vec{u}(𝜌, θ) = f(𝜌) \vec{e_θ} [/math] y que su presión 𝜌 es constante. Además el campo de velocidades tiene que cumplir la ecuación de Navier-Stokes estacionaria:

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

en la que µ es el coeficiente de viscosidad del fluido, y donde vamos a despreciar el primer término (parte convectiva). Obteniendo así:

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

4.2 Cálculo del Laplaciano del campo de velocidades

Para el cálculo del laplaciano vectorial en coordenadas cartesianas tenemos la siguiente formula:

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

Pero en este ejercicio el campo de velocidades está dado en la base cilíndrica, así que utilizaremos:

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

4.2.1 Gradiente de la divergencia

Si lo separamos por pasos, en el primer sumando tenemos el gradiente de la divergencia, pero para ello necesitamos calcular primero la divergencia. En nuestro caso, tenemos un fluido incompresible y dado que la divergencia mide el cambio en la densidad de un fluido moviéndose de acuerdo con un campo vectorial, será nulo.

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

Como la divergencia es nula el gradiente también lo es.

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

4.2.2 Rotacional del campo de velocidades

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

Tras haber calculado el rotacional del campo de velocidades, continuamos con el procedimiento calculando el rotacional del rotacional del campo de velocidades.

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

4.2.3 Calculo final

Tras haber calculado todas las partes de la ecuación, sustituimos en la definición de Laplaciano vectorial

[math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math] [math]\Delta\vec{u} = \vec{0} - \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]

4.3 Cálculos con la ecuación de Navier-Stocks

Dado que ya conocemos todas los partes de la ecuación, podemos resolverla, para hallar así el campo de velocidades:

[math](\vec{u} · ∇)\vec{u} + ∇p = µ∆\vec{u} [/math] [math] \vec{0} + {0}-µ∆\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]

4.3.1 Comprobación de que f(ρ) satisface una ecuación diferencial

La ecuación diferencial dada es:

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

si desarrollamos la primera parte de la ecuación tenemos:

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

El resultado obtenido coincide con el segundo y tercer sumando de la ecuación de Navier-Stocks, por lo tanto sustituyendo tenemos:

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

Y reordenando la ecuación, comprobamos lo que se nos pedía:

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

4.3.2 Comprobación de una solución conocida

Dada una solución posible solución:

[math]f(ρ) = aρ +\frac{b}{ρ},\hspace{20pt}a,b \in \mathbb{R} [/math]

Para poder comprobar que esta solución es válida, 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]

Tras hallar dichas derivadas, las introducimos en la ecuación obtenida anteriormente y verificamos 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]

4.3.3 Valores de "a" y "b"

Tenemos que buscar los valores "a" y "b" tal que [math]\vec{\mathbf{u}}[/math] coincide en las fronteras. Para ello se imponen las siguientes condiciones:

[math]\begin{cases}\vec{\mathbf{u}}(ρ=1)=\vec{\mathbf{\omega_i}} \times ρ\vec{\mathbf{e_\rho}} = \omega \rho \vec{\mathbf{e_\theta}} = \omega \vec{\mathbf{e_\theta}} \\ \vec{\mathbf{u}}(ρ=2)=\vec{\mathbf{\omega_e}} \times ρ\vec{\mathbf{e_\rho}} = -n \omega \rho \vec{\mathbf{e_\theta}} = - 2n \omega \vec{\mathbf{e_\theta}}\end{cases} \Longrightarrow \begin{cases}\vec{\mathbf{u}}(ρ=1)= (a\rho + \frac{b}{\rho}) \vec{\mathbf{e_\theta}} = (a+b) \vec{\mathbf{e_\theta}} \\ \vec{\mathbf{u}}(ρ=2)= (a\rho + \frac{b}{\rho}) \vec{\mathbf{e_\theta}}= (2a+\frac{b}{2})\vec{\mathbf{e_\theta}} \end{cases}[/math]

Que haciendo un resolución de un sistema de 2 ecuaciones con 2 incognitas al tener dos funciones ([math]\vec{\mathbf{u}}[/math] y [math]\vec{\mathbf{\omega}} \times \rho\vec{\mathbf{e_\rho}}[/math] ) nos tienes que salir que:

[math]\begin{cases} \omega = (a+b) \Longrightarrow a= \omega -b \\ -2n \omega = 2a + \frac{b}{2} \end{cases} \Longrightarrow [/math] [math]\Longrightarrow 2n\omega = -2\omega-2b+ \frac{b}{2} \Longrightarrow -(2n+2)\omega= \frac{-3b}{2} \Longrightarrow [/math] [math]\Longrightarrow b=\frac{4}{3}(n+1)\omega [/math]

Entonces sustituyo b en a y sale:

[math]\Longrightarrow a=\omega-\frac{4}{3}(n+1)\omega \Longrightarrow a=(\frac{4}{3}(n+1))\omega = (-\frac{4}{3}n -\frac{4}{3}+1)\omega =(-\frac{4}{3}n -\frac{1}{3})\omega \Longrightarrow a= -\frac{1}{3}(4n+1)\omega [/math]

quedando la ecuación diferencial como:

[math] \vec{\mathbf{u_\rho}}= -[ \frac{1}{3}(4n+1)\omega\rho-\frac{4}{3}(n+1)\omega\frac{1}{\rho} ]\vec{\mathbf{e_\theta}} \text{ tal que} \begin{cases} |\vec{\mathbf{\omega_i}}| = \omega \text{ (antihorario)}\\ |\vec{\mathbf{\omega_e}}| = n\omega \text{ (horario) }\end{cases}[/math]

5 Representación del campo de velocidades

suponiendo las condiciones del enunciado [math]|\vec{\mathbf{\omega_i}}|=|\vec{\mathbf{\omega_e}}|=1 \text{ y } \mu=1 [/math] y además tenemos que:

[math]b=\frac{4}{3}(n+1)\omega[/math] [math]a= -\frac{1}{3}(4n+1)\omega [/math]

Sustituimos en la función y obtenemos:

[math] \vec{\mathbf{u_\rho}}= -[ \frac{5}{3}\rho-\frac{8}{3}\cdot\frac{1}{\rho} ]\vec{\mathbf{e_\theta}} \text{ tal que } \rho \in [1,2] \begin{cases} -\frac{5}{3} \cdot 1 + \frac{8}{3}\cdot \frac{1}{1}= 1 =1\cdot1=\vec{\mathbf{\omega_i}}\times \rho\vec{\mathbf{e_\rho}} \\ -\frac{5}{3} \cdot 2 + \frac{8}{3}\cdot \frac{1}{2}= -2 =-2\cdot1=\vec{\mathbf{\omega_e}}\times \rho\vec{\mathbf{e_\rho}} \end{cases}[/math]

y su representación sería la siguiente:

Campo de velocidades

Código de Matlab:

rho=linspace(1,2,15);
theta=linspace(0,2*pi,20);
[RHO,THETA]=meshgrid(rho,theta);
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
f=(1/3)*((8./RHO)-(5.*RHO));                 
Vrho=zeros(size(RHO));    
Vx=-f.*sin(THETA);
Vy=f.*cos(THETA);
hold on;
quiver(x,y,Vx,Vy);
axis equal;
xlabel('x'); ylabel('y');
title('Campo de velocidades');
axis([-3 3 -3 3]);
hold off;

6 Representación de las líneas de corriente del campo

Para ello, necesitamos calcular el campo [math]\vec{\mathbf{v}}[/math] que en cada punto es ortogonal a [math]\vec{\mathbf{u}}[/math]

[math]\vec{\mathbf{v}}=\vec{\mathbf{k}}\times\vec{\mathbf{u}} \rightarrow \vec{\mathbf{v}}(\rho)= \vec{\mathbf{k}}\times -[\frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega\frac{1}{\rho}]\vec{\mathbf{e_\theta}}=[/math]

[math]=-[\frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega\frac{1}{\rho}]\vec{\mathbf{(-e_\rho)}}=[/math]

[math]=-[\frac{4}{3}(n+1)\omega\frac{1}{\rho}-\frac{1}{3}(4n+1)\omega\rho ]\vec{\mathbf{e_\rho}}[/math]

6.1 Demostración de que el campo [math]\vec{\mathbf{v}}[/math] es irrotacional

Comprobamos que, como dice el enunciado, el campo [math]\vec{\mathbf{v}}[/math] es irrotacional debido a que la divergencia de [math]\vec{\mathbf{u}}[/math] es nula, cosa que ya hemos calculado con anterioridad.

[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}(n+1)\omega\frac{1}{\rho}-\frac{1}{3}(4n+1)\omega\rho ]\vec{\mathbf{e_\rho}} & 0 & 0\end{vmatrix}=\vec 0[/math]

Se comprueba la irrotacionalidad de [math]\vec{\mathbf{v}}[/math] .

6.2 Cálculo de las líneas de corriente

Conocemos que:

[math] \vec{v} = \nabla\psi \Longrightarrow -[\frac{4}{3}(n+1)\omega\frac{1}{\rho}-\frac{1}{3}(4n+1)\omega\rho ]\vec{\mathbf{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]

Despejamos [math]\psi[/math]:

[math] \frac{\partial\psi}{\partial\rho} = -[\frac{4}{3}(n+1)\omega\frac{1}{\rho}-\frac{1}{3}(4n+1)\omega\rho ]\vec{\mathbf{e_\rho}} \Longrightarrow \psi = \int -[\frac{4}{3}(n+1)\omega\frac{1}{\rho}-\frac{1}{3}(4n+1)\omega\rho ]\vec{\mathbf{e_\rho}} \text{ } \partial \rho [/math]

Calculamos la integral y nos sale que [math]\psi[/math] es: [math] \psi= -\frac{4}{3}(n+1)\omega \cdot ln(\rho) + \frac{1}{3}(4n+1)\omega \cdot\frac{\rho^2}{2} + C, C\in R[/math]

6.3 Representación de las líneas de corriente de [math]\vec{\mathbf{u}}[/math]

Campo de velocidades

Código Matlab:

rho=linspace(1,2,200); 
theta=linspace(0, 2*pi,200);
[RHO,THETA]=meshgrid(rho,theta);
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
t=(5/6).*RHO.^2-(8/3).*log(RHO);
hold on;
axis equal;
grid on;
title('Líneas de corriente de u');
contour(X,Y,t,20,'LineWidth',1.5);
xlabel('x');
ylabel('y');
hold off;

7 Velocidad máxima del fluido

En este apartado, identificaremos los puntos donde la velocidad del fluido es máxima. Para ello, primero calcularemos el módulo del campo vectorial que corresponde a la velocidad: [math]\ \left | \vec u (\rho ) \right | = \left | -\left [\frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega \cdot \frac{1}{\rho} \right ] \vec e_\theta \right | = \frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega \cdot \frac{1}{\rho} [/math]

A continuación buscamos los valores máximos que dicho módulo puede tomar, es decir, estamos buscando puntos críticos. Por lo tanto, lo derivamos e igualamos a 0.

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

Como la ecuación obtenida en el punto anterior no puede cumplirse con valores reales, concluimos que la función del módulo no tiene máximos ni mínimos; es estrictamente creciente o estrictamente decreciente. Para determinar el máximo valor que puede tomar el módulo, evaluamos la función en los extremos del intervalo ρ∈(1,2):

[math] \left| \vec{u}(\rho=1) \right| = \frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega \cdot \frac{1}{1} = \omega [/math]

[math] \left| \vec{u}(\rho=2) \right| = \frac{1}{3}(4n+1)\omega\rho - \frac{2}{3}(n+1)\omega = 2n\omega [/math]

Entnces podemos concluir que si: [math] \left| \vec{u}(\rho=1) \right| \gt \left| \vec{u}(\rho=2) \right| \text{si } n\lt \frac{1}{2} \Rightarrow \vec{\mathbf{\omega_e}} \lt \frac{1}{2}\vec{\mathbf{\omega_i}} \Rightarrow \text{ Máx en } \left| \vec{u}(\rho=1) \right| [/math]

[math] \left| \vec{u}(\rho=1) \right| \lt \left| \vec{u}(\rho=2) \right| \text{si } n \gt \frac{1}{2} \Rightarrow \vec{\mathbf{\omega_e}} \gt \frac{1}{2}\vec{\mathbf{\omega_i}} \Rightarrow \text{ Máx en } \left| \vec{u}(\rho=2) \right| [/math]

[math] \left| \vec{u}(\rho=1) \right| = \left| \vec{u}(\rho=2) \right| \text{si } n = \frac{1}{2} \Rightarrow \vec{\mathbf{\omega_e}} = \frac{1}{2}\vec{\mathbf{\omega_i}} \Rightarrow \text{ Máx en } \left| \vec{u}(\rho=1) \right| \text{y } \left| \vec{u}(\rho=2) \right| [/math]

7.1 Representación grafica

Campo de velocidades

Código Matlab:

rho=linspace(1,2,400);          
theta=linspace(0,2*pi,200);
f=abs((1/3)*(-5.*rho+8./rho));
rhomin=sqrt(8/5);
fmin=abs((1/3)*(-5*rhomin+8/rhomin));
rhomax=2;
fmax=abs((1/3)*(-5*rhomax+8/rhomax));
[RHO,THETA]=meshgrid(rho,theta);
F=abs((1/3)*(-5.*RHO + 8./RHO));
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
Z=F;
figure;
subplot(1,2,1)
cmap=turbo(length(rho));
idx=round(rescale(f,1,length(rho)));
hold on
for i=1:length(rho)-1
   plot(rho(i:i+1),f(i:i+1),'Color',cmap(idx(i),:),'LineWidth',2);
end
h1=plot(rho_min,fmin,'ko','MarkerSize',8,'LineWidth',2);
h2=plot(rho_max,fmax,'ko','MarkerSize',8,'LineWidth',2);
text(rhomin,fmin,sprintf('   Mínimo = %.f',fmin),'VerticalAlignment','bottom','Color','k','FontSize',10,'FontWeight','bold');
text(rhomax,fmax,sprintf('  Máximo = %.f',fmax),'VerticalAlignment','top','Color','k','FontSize',10,'FontWeight','bold');
xlabel('\rho');
ylabel('

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

8.1 Cálculo del rotacional de [math] \vec{u} [/math]

Para calcular el rotacional de un campo vectorial en coordenadas cilíndricas, se utiliza la siguiente expresión:

[math] \nabla \times \vec{u} = -\frac{1}{\rho} \begin{vmatrix} \vec e_{\rho} & (\rho \vec e_{\theta}) & \vec e_{z} \\ \dfrac{\partial}{\partial \rho} & \dfrac{\partial}{\partial \theta} & \dfrac{\partial}{\partial z} \\ u_{\rho} & {\rho}u_{\theta} & u_{z} \end{vmatrix} [/math]

En este caso, conociendo [math] \vec{u} [/math] : Siendo [math]\vec{u}(ρ, θ) = f(ρ) \vec{e_θ} [/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}{ρ} \left [\vec{-e_\rho} \left (\frac{\partial \left [\frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega \cdot \frac{1}{\rho} \right ] \vec e_\theta}{\partial z} \right) + \vec{e_z} \left (\frac{\partial \left [\frac{1}{3}(4n+1)\omega\rho - \frac{4}{3}(n+1)\omega \cdot \frac{1}{\rho} \right ] \vec e_\theta}{\partial \rho} \right) \right ]=[/math] [math] = -\frac{1}{ρ} \left [\vec{e_z} \left (\frac{2}{3} (4n+1)\omega\rho - 0 \right ) \right ] = \left [-\frac{2}{3} (4n+1)\omega \right] \vec{e_z} [/math]

Por tanto, el rotacional es constante en todo el dominio.

8.2 Representación grafica

Módulo del rotacional

Código Matlab:

a=5/3;   
mr=2*a;
rho=linspace(1,2,300);        
theta=linspace(0,2*pi,300);
[RHO,THETA]=meshgrid(rho,theta);
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
C=mr*ones(size(RHO));
figure;
pcolor(X,Y,C);
shading interp;
colorbar;
axis equal;
axis([-3 3 -3 3]);
title('Módulo del rotacional de U');
xlabel('x');
ylabel('y');

9 Campo de Temperaturas

La temperatura del fluido viene definida por el siguiente campo:

[math] T(\rho,\theta) = log(1+ \rho^2)\cos^2 \theta [/math]

Y la temperatura máxima que puede alcanzar se puede ver en la grafica como la zona mas amarilla.

9.1 Representación gráfica

Dado el campo de Temperaturas, la representación del campo y de las curvas de nivel es el siguiente:

Campo de temperatura

Código Matlab:

rho=linspace(1,2,400);
theta=linspace(0,2*pi,400);
[RHO,THETA]=meshgrid(rho,theta);
T=log10(1+RHO.^2).*(cos(THETA)).^2; 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
rhomax=2;
thetamax=[0, pi];
xmax=rhomax*cos(thetamax);
ymax=rhomax*sin(thetamax);
Tmax=log10(1+rhomax^2);
fprintf('Temperatura máxima = %.6f\n',Tmax);
for i=1:2
   fprintf('Máx %d: rho = %.6f, theta = %.6f rad, (x,y) = (%.6f, %.6f)\n',i,rhomax,thetamax(i),xmax(i),ymax(i));
end
figure;
contour(X,Y,T,40); 
colorbar;
axis equal tight;
title('Curvas de nivel de T(\rho,\theta)');
xlabel('x');
ylabel('y');
hold on;
clim([min(T(:)),max(T(:))]); 
colormap(turbo); 
for i=1:2
   plot(xmax(i),ymax(i),'ko','MarkerSize',10,'LineWidth',2);
  
   if xmax(i)>0
       text(xmax(i)-0.8,ymax(i),sprintf(' Máx %.3f',Tmax),'Color','k','FontSize',11,'FontWeight','bold');
   else
       text(xmax(i)+0.08,ymax(i),sprintf(' Máx %.3f',Tmax),'Color','k','FontSize',11,'FontWeight','bold');
   end
end
grid on
figure;
surf(X,Y,T,'EdgeColor','none'); 
colormap(turbo); 
colorbar;
axis equal;
xlabel('x');
ylabel('y');
zlabel('T');
title('Campo escalar T(\rho,\theta)');
view(25,15);
clim([min(T(:)), max(T(:))]);
hold on;
for i=1:2
   plot3(xmax(i),ymax(i),Tmax,'go','MarkerSize',8,'LineWidth',2);
   text(xmax(i),ymax(i),Tmax,sprintf('  Máx %.3f',Tmax),'Color','g','FontSize',11,'FontWeight','bold');
end

10 Gradiente de la temperatura

Primero definimos en coordenadas cilíndricas el campo vectorial gradiente de un campo escalar de la siguiente manera:

[math]\triangledown T = \frac{ \partial T}{\partial \rho}\vec{e_\rho} + \frac{1}{\rho}\frac{ \partial T}{\partial \theta}\vec{e_\theta} + \frac{ \partial T}{\partial z}\vec{e_z}[/math]

Pero nuestra función de temperatura T solo depende de [math] \rho [/math] y [math] \theta [/math]

[math]\triangledown T(\rho,\theta)= \cos^{2}\theta \cdot \frac{2\rho}{(1 + \rho^{2}) ln10} \cdot \vec{e}_\rho - \frac{\log(1 + \rho^{2})\cdot2\sin\theta\cdot\cos\theta}{\rho} \cdot \vec{e}_\theta [/math]

10.1 Demostración que la gráfica del gradiente de temperatura sea ortogonal a las curvas de nivel de la temperatura

En la primera imagen se visualiza el campo gradiente, en la segunda las curvas de nivel de la temperatura y finalmente, la comprobación de que el campo gradiente es ortogonal a las líneas de nivel de la temperatura.

• 

  Gradiente

• 

  Líneas de nivel

• 

  Gradiente y líneas de nivel

Codigo Matlab:

[x,y]=meshgrid(linspace(-2.2,2.2,401),linspace(-2.2,2.2,401));
rho=hypot(x,y);     
theta=atan2(y,x);
T=log10(1+rho.^2).*(cos(theta).^2);
Tr=(2.*rho.*(cos(theta).^2))./((1+rho.^2)*log(10));
Ttheta=-2*cos(theta).*sin(theta).*log10(1+rho.^2);
gx=Tr.*cos(theta)-(Ttheta./rho).*sin(theta); 
gy=Tr.*sin(theta)+(Ttheta./rho).*cos(theta);
mask=(rho>=1)&(rho<=2);
T(~mask)=NaN;
gx(~mask)=0; 
gy(~mask)=0;  
figure('Color','w');

hold on;
axis equal;
box on;
grid on;
contour(x,y,T,20);
colormap(turbo);
step=16;
quiver(x(1:step:end,1:step:end),y(1:step:end,1:step:end),gx(1:step:end,1:step:end),gy(1:step:end,1:step:end),0.8,'Color',[1 0 0],'LineWidth',0.8);
title('Curvas de nivel y gradiente');
xlabel('x');
ylabel('y');
xlim([-2.2 2.2]);
ylim([-2.2 2.2]);
figure('Color','w');

hold on;
axis equal;
box on;
grid on;
quiver(x(1:step:end,1:step:end),y(1:step:end,1:step:end),gx(1:step:end,1:step:end),gy(1:step:end,1:step:end),0.8,'Color',[1 0 0],'LineWidth',0.8);
title('Gradiente');
xlabel('x');
ylabel('y');
xlim([-2.2 2.2]);
ylim([-2.2 2.2]);

11 Caudal que pasa por la sección θ=0

11.1 Cálculo del caudal

La sección [math]\theta = 0[/math] es una línea radial desde [math]ρ = 1[/math] hasta [math]ρ = 2[/math]. El fluido se mueve dependiente de [math]u_\theta[/math], así que el flujo atraviesa esa línea perpendicularmente.
El caudal [math]Q[/math] a través de esa sección (profundidad [math]h = 1 \,\text{m}[/math]) se calcula como: [math]Q = \int_{A} u_\theta(ρ)\, dA = \int_{z=0}^{z=1}\int_{ρ=1}^{ρ=2} u_\theta(ρ)\, dρ dz, [/math] y haciendo la integral donde Z=altura y [math]\rho =base[/math] nos sale que: [math]Q=-\frac{1}{2}(4n+1)\omega+\frac{4}{3}(n+1)\omega\cdot ln(2) [/math]