Flujo de Couette entre dos tubos concéntricos (Grupo 35)

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos. Grupo 35
Asignatura	Teoría de Campos
Curso	2023-24
Autores	David Astorga Trejos Abigail Perez Segura Carmen Efa Nchama Julia Gómez-Nieves Carlos Carramiñana
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos el flujo de un fluido incompresible a través de dos cilindros concéntricos de manera que el interior se mueve con una velocidad angular constante en sentido antihorario mientras que el exterior se encuentra fijo. Suponiendo que ambos cilindros tienen su eje en [math]OX_3[/math] y pintamos la sección transversal [math](x_3 = 0)[/math] el cilindro exterior queda proyectado sobre la la circunferencia ρ = 2 y el interior sobre la circunferencia ρ = 1. La velocidad angular cilindro interior es ω > 0.

1 Mallado de la sección transversal

Primero debemos representar el fluido definido en el enunciado acorde a los parámetros especificados en el enunciado. Las directrices requeridas son las siguientes: Tomar los ejes (comando axis) en el rectángulo [math](x, y) ∈ [-3;3] × [-3;3][/math] y para ver desde arriba view(2).

Figura 1: Mallado del flujo de Couette entre dos cilindros concéntricos.

Para representar este campo se sigue el código dado a continuación:

rho=1:0.1:2;                       % Intervalos de rho [1,2]
theta=0:0.1:2*pi;                  % Intervalo de theta [0,2*pi]
[RHO,THETA]=meshgrid(rho,theta);   % Creación de matrices de los parámetros
figure(1)
x=RHO.*cos(THETA);                 % Parametrización de x
y=RHO.*sin(THETA);                 % Parametrización de y
mesh(x,y,0*x)                      % Dibujo de la figura con la coordenadas ya parametrizadas
axis([-3,3,-3,3])                  % Región de vista del dibujo
view(2)

2 Ecuación de Navier-Stokes

2.1 Definición del campo de velocidades

La velocidad de las partículas del fluido viene dada por [math]\vec u(\rho,\theta)=f(\rho)\vec e_\theta[/math] y su presión p (constante). Sabemos que [math](\vec u, p)[/math] satisface la ecuación de Navier-Stokes estacionaria:

[math](\vec u \cdot \triangledown) \vec u + \triangledown p= \mu \triangle \vec u [/math] donde µ es el coeficiente de viscosidad del fluido. Debido a que la presión es constante: [math]\triangledown p= \vec 0[/math]

y que se desprecia el primer termino, parte de la ecuación estacionaria, tenemos la siguiente igualdad: [math]\mu\triangle\vec u=\vec 0[/math] .

2.2 Cálculo del Laplaciano

En esta igualdad vamos a usar el Laplaciano [math]\triangle \vec u[/math] en cilíndricas ya que es así como nos lo da el enunciado, por lo tanto la formula se nos queda en: [math]\triangle\vec u = \triangledown(\triangledown\cdot\vec u)-\triangledown\times(\triangledown\times\vec u)[/math]

Usamos esta alternativa ya que no hemos visto el laplaciano de un campo de vectores cuando ese campo esta definido en la base cartesiana, tambien porque con esta ecuacion podemos usar la divergencia y el rotacional, los cuales son requeridos a lo largo del problema.

2.2.1 Cálculo de la divergencia

Calculamos la divergencia de [math]\vec u[/math]: [math]\triangledown\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial\rho}(0)+{\frac{\partial}{\partial\theta}(f(\rho))}+\frac{\partial}{\partial z}(0)=0][/math] Esto nos indica la condición de incompresibilidad, que sucede al ser la divergencia igual a 0.

2.2.2 Cálculo del rotacional

Ahora procedemos a calcular el rotacional de [math]\vec u[/math], y posteriormente el rotacional del rotacional:

[math]\ \triangledown \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| = \left [ \frac{f\left ( \rho \right )}{\rho } + \frac{\partial \left ( f\left ( \rho \right ) \right )}{\partial \rho }\right ] \bar{e_{z}}= \left ( \frac{f\left ( \rho \right )}{\rho } + {f}'\left ( \rho \right )\right ) \bar{e_{z}} [/math] ;

[math]\ \triangledown \times \left ( \triangledown \times \vec{u} \right )=\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} \\ 0& 0 &\frac{f\left ( \rho \right )}{\rho }+{f}'\left ( \rho \right ) \end{vmatrix}= -\frac{1}{\rho }\left [ -\frac{1}{\rho ^{2}}f\left ( \rho \right ) +\frac{{f}'\left ( \rho \right )}{\rho }+{f}''\left ( \rho \right )\right ]\left ( \rho \right ) \vec{e}_{\theta } = \left [\frac{1}{\rho ^{2}}f\left ( \rho \right )-{\frac{1}{\rho}f'\left ( \rho \right )} - f''\left ( \rho \right )\right ] \vec{e}_{\theta } [/math]

2.2.3 Cálculo final

Con el rotacional y el rotacional del rotacional, tendríamos todos los cálculos necesarios para obtener el Laplaciano [math]\triangle \vec u[/math] , y por lo tanto la ecuación de Navier Stokes.

[math]\ \triangle \vec{u}=-\triangledown \times \triangledown \times \vec{u} =0-\left [\frac{1}{\rho ^{2}}f\left ( \rho \right )-{\frac{1}{\rho}f'\left ( \rho \right )} - f''\left ( \rho \right )\right ] \vec{e}_{\theta } = -\frac{1}{\rho }\left [ -\frac{1}{\rho ^{2}}f\left ( \rho \right ) +\frac{{f}'\left ( \rho \right )}{\rho }+{f}''\left ( \rho \right )\right ]\left ( \rho \right ) \vec{e}_{\theta } = -\frac{1}{\rho }\left [-\frac{1}{\rho}f\left ( \rho \right )+{f'\left ( \rho \right )} +\rho f''\left ( \rho \right )\right ] \vec{e}_{\theta } [/math]

[math] → \ \triangle \vec{u}= -\frac{1}{\rho }\left [ -\frac{1}{\rho ^{2}}f\left ( \rho \right ) +\frac{{f}'\left ( \rho \right )}{\rho }+{f}''\left ( \rho \right )\right ]\left ( \rho \right ) \vec{e}_{\theta } = \left [\frac{1}{\rho ^{2}}f\left ( \rho \right )-{\frac{1}{\rho}f'\left ( \rho \right )} - f''\left ( \rho \right )\right ] \vec{e}_{\theta } [/math]

2.3 Ecuación de Navier-Stokes

La ecuación sería la siguiente:

[math]\mu \triangle \vec{u}= 0 \rightarrow \left [\frac{1}{\rho ^{2}}f\left ( \rho \right )-{\frac{1}{\rho}f'\left ( \rho \right )} - f''\left ( \rho \right )\right ]\vec{e}_{\theta }=\vec{0} [/math]

Después de calcular la Ec. de Navier-Stokes:

Tenemos la siguiente ecuación diferencial: [math]\ \frac{\partial }{\partial \rho }\left ( \rho \frac{\partial f\left ( \rho \right )}{\partial \rho } \right )= \frac{f(\rho)}{\rho} [/math]

2.3.1 Comprobante de que f(ρ) cumple la ecuación diferencial

[math]\ \frac{\partial }{\partial \rho }\left ( \rho \frac{\partial f\left ( \rho \right )}{\partial \rho } \right )= \frac{\partial }{\partial \rho } (\rho\cdot f'(\rho))=f'(\rho)+\rho(f''(\rho)) \rightarrow f'(\rho)+\rho(f''(\rho))=\frac{f(\rho)}{\rho} [/math]

Si pasamos el cociente [math]\frac{f(\rho)}{\rho}[/math] al otro lado de la igualdad, vemos que coincide con la ecuación calculada anteriormente.

2.3.2 Comprobante que la función es solución

Ahora vamos a comprobar si la función [math]f(\rho)=a\rho+\frac{b}{\rho}[/math] es solución de la ecuación diferencial mencionada anteriormente, en el caso de ser así, la velocidad de la frontera del fluido coincida con la frontera de los cilindros, tanto interior como exterior.

[math]\ \frac{\partial }{\partial \rho }\left ( \rho \frac{\partial f\left ( \rho \right )}{\partial \rho } \right )= \frac{\partial }{\partial \rho } (\rho\cdot (a-\frac{b}{\rho^2}))= \frac{\partial }{\partial \rho }(\rho\cdot a-\frac{b}{\rho }) = a + \frac{b}{\rho ^2} [/math]

[math]f'(\rho)= a-\frac{b}{\rho^2}[/math]

[math]f''(\rho)=\frac{2b}{\rho^3}[/math]

[math]\ \frac{1}{\rho}(a\cdot\rho +\frac{b}{\rho })=a-\frac{b}{\rho ^2}+\rho (\frac{2b}{\rho^3}) = a+\frac{b}{\rho^2} [/math]

Se comprueba que la función mencionada si que es solución, ya que coincide con la ecuación diferencial calculada:

[math]\ \frac{1}{\rho}(a\cdot\rho +\frac{b}{\rho })=a+\frac{b}{\rho ^2}\rightarrow \frac{f(\rho)}{\rho}=f'(\rho)+\rho f''(\rho) [/math]

2.3.3 Determinar los valores "a" y "b"

Recordamos que la formula de la velocidad viene dada por [math]\vec u(\rho,\theta)=f(\rho)\vec e_\theta[/math] y [math]f(\rho)=a\rho+\frac{b}{\rho}[/math] .

Para hallar estos valores, tenemos que adaptar la función a las características del problema, es decir, darle los valores de [math]\rho[/math] de las paredes de la tubería en las dos condiciones, los cuales son para [math]\rho =1[/math] y para [math]\rho =2[/math], e igualar la función a las velocidades que tienen los cilindros interior(w) e exterior(0):

Condición 1: Si [math]ρ=1 \Longrightarrow \vec{u}=\omega \vec{e_θ} \Longrightarrow f(1)\vec{e_θ}=\omega \vec{e_θ} \Longrightarrow (a+b)\vec{e_θ}=\omega \vec{e_θ}\Longrightarrow a+b=\omega[/math]

Condición 2: Si [math]ρ=2 \Longrightarrow \vec{u}= \vec{0} \Longrightarrow f(2)\vec{e_θ}=\vec{0} \Longrightarrow (2a+\frac{b}{2})\vec{e_θ}=\vec{0}\Longrightarrow 2a+\frac{b}{2}= 0[/math]

Obtenemos como solucion de cada ecuación, un sistema de dos ecuaciones, iniciando con que [math]\ (a = -\frac{b}{4}\omega)[/math]

Con lo que tenemos finalmente obtenemos que: [math]\ a = -\frac{1}{3}\omega [/math] y [math]\ b = \frac{4}{3}\omega [/math]

Quedando la ecuación diferencial de Navier-Stokes de siguiente forma:

[math]f(ρ) = -\frac{1}{3}\omega ρ +\frac{4}{3ρ}\omega[/math]

3 Campo de velocidades

Cuando [math]\ \omega =1 [/math] y [math]\ \mu=1 [/math] el campo de velocidades resulta ser el siguiente:

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

Figura 2. Representación gráfica del campo de velocidades

Para representar este campo se sigue el código dado a continuación:

r=1:0.1:2;                     %Rangos en que se mueven las variables.
t=0:0.1:2*pi;         
[R,T]=meshgrid(r,t);           %Matrices de coordenadas.
x=R.*cos(T);                   %Parametrización en coordenadas cartesianas.
y=R.*sin(T);
X=sin(T).*(-(1/3)*(-R+4./R));  %Evaluación del campo en la parametrización.
Y=cos(T).*((1/3)*(-R+4./R));
axis([-3,3,-3,3])      
view(2)
quiver(x,y,X,Y);               %Representación del campo.

4 Líneas de corriente

Estas líneas serán tangentes a las del campo de velocidades en cada punto por lo que buscamos el campo vectorial [math] \vec v [/math],perpendicular a [math]\vec u[/math]:

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

En el enunciado se indica que este campo irrotacional, lo cual se puede ver también al calcular el rotacional:

[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{\rho}{3} & 0 & 0\end{vmatrix}=\vec 0[/math]

La función de corriente del campo de velocidades es una función escalar cuyo gradiente es [math]\vec v[/math]:

[math]\psi: \vec v=\triangledown\psi=\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{\rho}{3})\vec e_\rho[/math]

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

[math]\psi=-\frac{4ln(\rho)}{3}+\frac{\rho^2}{6}[/math]

Figura 3. Líneas tangentes a las del campo de velocidades.

Una vez calculado el potencial podemos representarlo con el siguiente código:

int=100;                       % Número de intervalos
r=linspace(1,2,int);           % Parametro r [1,2]
t=linspace(0,2*pi,int);        % Parametro t [0,2*pi]
[R,T]=meshgrid(r,t);         % Matrices de coordenadas 
potencial=((-4*log(R)./3)+((R.*R)./6));      % Campo escalar 
x=R.*cos(T);                 % Parametrización
y=R.*sin(T);
axis([-3,3,-3,3])            % Ejes
view(2)                     
contour(x,y,potencial,20);colorbar   % Líneas de nivel
axis([-3,3,-3,3])            % Ejes del dibujo
colormap("cool")

5 Velocidad máxima del fluido

Primero tomamos el modulo de la velocidad: [math]\ \left | \vec u (\rho ) \right | = \left | \left ( \frac{1}{3}(-\rho + \frac{4}{\rho}) \right ) \vec e_\theta \right | = \frac{1}{3} \left ( -\rho +\frac{4}{\rho } \right ) [/math] Después derivamos e igualamos a 0 para hallar la velocidad máxima: [math]\ \frac{\partial \left | \vec u (\rho ) \right | }{\partial \rho } = \frac{1}{3} \left ( -1 - \frac{4}{\rho^2} \right ) = 0 [/math] La solución sale: [math]\ \sqrt{-4} [/math] , que no es una solución real, así que analizamos los extremos del intervalo: [math]\ \left\{\begin{matrix} \left | \vec u (1) \right | = \left | \left ( \frac{1}{3} (-\rho + \frac{4}{\rho}) \right ) \vec e_\theta \right | = \frac{1}{3} \left | ( 1 -\frac{4}{1 } \right ) | = 1 \\ \left | \vec u (2) \right | = \left | \left ( \frac{1}{3} (-\rho + \frac{4}{\rho}) \right ) \vec e_\theta \right | = \frac{1}{3} \left | ( 2 -\frac{4}{2 } \right ) | = 0 \end{matrix}\right. [/math] Con estos resultados vemos que la velocidad máxima se alcanza en [math]\ \rho = 1 [/math] con valor máximo 1

Figura 4: Variación del módulo de la velocidad en función de ρ

clear, close all
%definimos rho y theta
rho=1:0.1:2;
theta=0:0.1:2*pi;

%Se genera una retícula rectangular a partir de las divisiones rho y theta
[RHO,THETA]=meshgrid(rho,theta);

%pasamos a cartesianas
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);

%Módulo de la velocidad
Modulo = abs((1/3)*(-RHO + 4./RHO));

%Representación del módulo de la velocidad
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 vemos el comportamiento del módulo velocidad en función de rho y se ve lo que se demostró antes, que la velocidad máxima se obtiene en [math]\ \rho = 1 [/math] y su valor es 1. También vemos que en [math]\ \rho = 2 [/math] la velocidad es nula.

6 ROTACIONAL

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

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

Por lo tanto, el rotacional queda: [math]\ \bigtriangledown \times \vec u = \left (\frac{f(\rho )}{\rho} + f'(\rho ) \right )\vec e_z = \left ( \frac{ \frac{1}{3} \left ( -\rho + \frac{4}{\rho } \right ) }{\rho} + \frac{1}{3} \left (-1 - \frac{4}{\rho ^2} \right ) \right ) \vec e_z = \frac{1}{3} \left ( 1 - \frac{4}{\rho ^2} - 1 - \frac{4}{\rho ^2} \right )\vec e_z = -\frac{2}{3} \vec e_z [/math]

Así, el módulo del vector queda: [math]\ \left | \bigtriangledown \times \vec u \right | = \sqrt{\left ( -\frac{2}{3} \right )^2} = \frac{2}{3} [/math]

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

Figura 5: Representación del rotacional

rho=1:0.33:2;
theta=0:0.33:2*pi;

%se genera una reticula rectancular donde se representan rho y theta
[RHO,THETA]=meshgrid(rho,theta);
  
%Pasoamos a cartesianas
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
z=ones(size(y)).*(2/3);

%Representación del módulo del rotacional
surf(x,y,z);
axis([-3,3,-3,3]);
view(2)
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{2}{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 Campo de Termperaturas

Suponemos que la temperatura del fluido viene dada por el campo escalar: [math]\ T(\rho,\theta) = 1 + \rho^{2} \sin^{2} \theta e^{-\left(\rho-\frac{3}{2}\right)^{2}} [/math]

Usando el código de matlab mostraremos la interpretación gráfica del campo escalar que representa la temperatura del fluido. Para representar el campo escalar se hará con la ayuda de Matlab y así estudiará el gradiente de la temperatura. Hemos hecho dos graficas para poder visualizar correctamente la grafica. La primera la veremos en 3D, observando donde se puede situar el máximo de la temperatura. En la segunda imagen podemos observar la misma figura pero representada en el eje x,y (usando pcolor), teniendo una vista desde arriba donde contemplamos los tonos tanto de frio como de calor. Para poner la barra de temperatura a la derecha ponemos (colorbar). Y terminamos insertando los ejes (_label) y el titulo (title)

Representamos primero el campo T aunque lo he denominado Temp en el código:

%Se definen las variables rho y theta
rho=1:0.1:2;
theta=0:0.1:2*pi;
%Se genera una retícula rectangular a partir de las divisiones rho y theta
[RHO,THETA]=meshgrid(rho,theta);
%La temperatura del fluido viene dada por la expresión
Temp=1+RHO.^2.*sin(THETA).^2.*exp(-(RHO-3/2).^2);

Figura 6: Representación de la temperatura en 3D.

7.1 Representación del campo 3D

%Representación del campo en 3D
figure(1)
mesh(RHO,THETA,Temp)
surf(RHO,THETA,Temp); % Dibujo de la superficie 
shading flat
colorbar;             %Barra de temperatura a la derecha
title('Representación de la temperatura en 3D');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');

Figura 7: Representación del campo y de las curvas de nivel en 2D.

7.2 Representación del campo y las curvas de nivel en 2D

%Representación del campo y de las curvas de nivel en 2D
figure(2)
mesh(RHO,THETA,Temp)
pcolor(uu,vv,T);       % Proyección en planta
hold on                % Mantener ventana
contour(uu,vv,T,'k');
hold off
colorbar;              %Barra de temperatura a la derecha
title('Representación de la temperatura en 2D');
xlabel('Eje X');
ylabel('Eje Y');

7.3 Punto máximo de la Temperatura

Figura 8: Representación del campo y de las curvas de nivel en 2D.

Considerando el rango de temperaturas, la descripción visual del campo y de las líneas que conectan puntos de igual temperatura, se puede ver la temperatura máxima en la siguiente imagen, representado con matlab :

rho = 1:0.1:2; % la variable rho definida en el intervalo [1,2]
theta =0:0.1:2*pi;
%la variable teta descrita en el intervalo [0,2*pi]
[RHO,THETA]=meshgrid(rho,theta);%Se genera una retícula a partir de las divisiones rho y theta
Temp=1+RHO.^2.*sin(THETA).^2.*exp(-(RHO-3/2).^2);
plot3(RHO,THETA, Temp);
title('gráfica del maximo punto de temperatura');

8 Gradiente de la Temperatura

En primer lugar, establecemos la representación en coordenadas cilíndricas del campo vectorial gradiente de un campo escalar de la siguiente forma:

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

Podemos observar en la formula de la temperatura, que esta solo depende de ρ y θ: [math] T(\rho,\theta) = 1 + \rho^{2}sen^{2}\theta e^{ - (\rho - \frac{3}{2})^{2}}[/math]

Por lo tanto, como no dependemos de z, sabemos que [math]\frac{ \partial T}{\partial z}\vec{e_z}=0[/math]

Haciendo todos los calculos, el gradiente nos queda como: [math]\triangledown T = [sen^{2}\theta[(- 2\rho^{3} + 3\rho^{2}) e^{ - (\rho - \frac{3}{2})^{2}}+ 2\rho e^{ - (\rho - \frac{3}{2})^{2}} ]]\vec{e_\rho} + [sen(2\theta)\rho^{2} e^{ - (\rho - \frac{3}{2})^{2}}]\vec{e_\theta}[/math]

Posteriormente seguimos con la representación gráfica del gradiente en Matlab:

Figura 9: Representación del gradiente de la temperatura.

% intervalo para dividir los parámetros 
r=1:0.1:2;                   % intervalos de u [1,2]
t=0:0.1:2*pi;              % intervalo de v [0,2*pi]
[rr,tt]=meshgrid(r,t);     % creacion de matrices de los parámetros
T= 1+rr.^2.*(sin(tt.*exp((-(rr-3/2).^2)))).^2;
[dx,dy] = gradient(T,r,t); % Gradiente de la temperatura

subplot(1,2,1);             % Ventana 1
quiver3(rr,tt,dx,dy);       % Dibujo del campo vectorial
title('gradiente de la temperatura en 3D') %Titulo de la grafica
xlabel('Eje X');    %Eje X de la gráfica
ylabel('Eje Y');    %Eje Y de la gráfica
zlabel('Eje Z');    %Eje Z de la gráfica

subplot(1,2,2);            % Ventana 2
quiver(rr,tt,dx,dy);       % Dibujo del campo vectorial
title('gradiente de la temperatura en 2D') %Titulo de la grafica
xlabel('Eje X');    %Eje X de la gráfica
ylabel('Eje Y');    %Eje Y de la gráfica

8.1 Ortogonalidad del gradiente a las curvas de nivel de la temperatura

Usando el codigo de matlab queda demostrado gráficamente q el gradiente es perpendicular a las curvas de nivel de la temperatura dada en el apartado anterior:

Figura 10:Ortogonalidad del gradiente a las curvas de nivel de la temperatura .

clear all                                                
  clc

  r = 1:0.1:2; % Genera una cuadrícula en coordenadas polares a partir de las variables                                          
  th = 0:0.1:2*pi;
  [R, TH] = meshgrid(r, th);

  % Calcula la función de temperatura en coordenadas polares
  Z = 1 + R.^2 .* (sin(TH .* exp((-(R - 3/2).^2)))).^2;  

 [Gx, Gy] = gradient(Z, r, th);  % Calcula el gradiente de Z con respecto a r y th

                                                                        
 figure;  % Crea una figura y dibuja las curvas de nivel de Z
 contour(R, TH, Z);
 title('Curvas de Nivel y campo de vectores');
 xlabel('r');
 ylabel('\theta');
 hold on;
 quiver(R, TH, Gx, Gy, 'Color', 'r');  % Dibuja un campo de vectores sobre las curvas de nivel

9 Caudal

Por último se calculará el caudal que atraviesa los dos cilindros por una sección longitudinal, en este caso la intersección del plano [math]\ x_2=0 [/math].

Suponiendo que los dos cilindros tiene 1 metro de profundidad y la velocidad del fluido, con el campo calculado en el apartado 3, en m/s vine dado por la forma: [math]\ \vec u(\rho ,\theta ) = (-\frac{1}{3}\rho+\frac{4}{3\rho})\vec{e_{\theta }} [/math]

Para el caudal del campo a través de una superficie [math]\ S [/math] (la sección longitudinal) se integrara según al formula:

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

Para conseguir un caudal positivo se tomara como orientación de la superficie la dirección de [math]\ \vec e_\theta [/math].

Para comenzar hay que visualizar la superficie a parametrizar.

Figura 11: Representación de la superficie

%cilindro interior
zz=0:0.1:1;              
theta=0:0.1:2*pi;
[ZZ,THETA]=meshgrid(zz,theta);
a=cos(THETA);        
b=sin(THETA);
c=ZZ;
surf(a,b,c,'FaceAlpha',0.3)   
hold on

%cilindro exterior 
d=2*cos(THETA);        
e=2*sin(THETA);
f=ZZ;
surf(d,e,f,'FaceAlpha',0.2)

%Superficie 1
rr=1:0.1:2;
hh=0:0.1:1;
[RR,HH]=meshgrid(rr,hh);
h=RR;
i=0*RR;
j=HH;
surf(h,i,j)

%superficie 2
k=-RR;
l=0*RR;
m=HH;
surf(k,l,m)
hold off

%ejes
xlabel('Eje X1');
ylabel('Eje X2');
zlabel('Eje X3');

Como se puede ver en la figura, la superficie [math]\ S [/math] está compuesta por la suma de dos rectangulos [math]\ S_1 [/math] y [math]\ S_2 [/math]. Comenzamos parametrizando cada una de llas.

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)) = (u,0,v) [/math] con [math]\ \left\{\begin{matrix} u\in (1,2) \\ v\in (0,1) \end{matrix}\right. [/math]

Los vectores velocidad se calculan con la fórmula:

[math]\ \Phi _t=\rho '(t)\cdot \vec e_\rho + \rho (t)\cdot \vec e_\theta + z (t)\cdot \vec e_z [/math]

[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)) = (u,\pi,v) [/math] con [math]\ \left\{\begin{matrix} u\in (1,2) \\ v\in (0,1) \end{matrix}\right. [/math]

Para los vectores velocidad se vuelve a usar la misma fórmula:

[math]\ \Phi _t=\rho '(t)\cdot \vec e_\rho + \rho (t)\cdot \vec e_\theta + z (t)\cdot \vec e_z [/math]

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

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

Después 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_{0}^{1}\int_{1}^{2} \vec u (u,0,v)\cdot (-\vec e_\theta)dudv - \int_{0}^{1}\int_{1}^{2} \vec u (u,\pi,v)\cdot (-\vec e_\theta)dudv = [/math]

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

[math]\ = -\frac{2}{3} \int_{0}^{1}dv \int_{1}^{2}\left ( u-\frac{4}{u}\right )du = -\frac{2}{3} \left [ \frac{u^2}{2} - 4\ln \left | u \right | \right ]^2_1 = -\frac{2}{3} \left ( \frac{4}{2} - 4\ln 2 - \frac{1}{2} \right ) = -1 + \frac{8}{3} \ln 2 \left [m^3/s \right ][/math]