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

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos. Grupo 22
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Miguel Ángel Abad Robles Rafael Moreno Orellana Juan Perestrello Gilabert Lucía Rodríguez Montes Miguel Rubio Arraztio
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Flujo de Couette entre dos tubos concéntricos. Vamos a considerar el flujo de un fluido incompresible a través de dos cilindros concéntricos de manera que el interior se mueve con velocidad angular constante en sentido antihorario mientras que el exterior está fijo. Si suponemos que ambos cilindros tienen su eje en OX3 y pintamos la sección transversal (x3 = 0) el cilindro exterior queda proyectado sobre la la circunferencia ρ = 2 y el interior sobre la circunferencia ρ = 1. La velocidad angular cilindro interior es ω > 0. En el archivo adjuntado podemos observar una representación gráfica de los puntos ocupados por el fluido (Figura 1).

Para su obtención, hemos utilizado el siguiente código de Matlab.

clear all   % Limpiar la pantalla, espacio de trabajo y habilitar el modo de retención de gráficos
clc
hold on
radio_malla = 1:0.1:2;                          % Definir los parámetros para la malla de círculos
angulo_malla = 0:0.1:2*pi;           
[R_mesh, A_mesh] = meshgrid(radio_malla, angulo_malla); 
x_mesh = R_mesh .* cos(A_mesh);                    % Generar coordenadas para la malla de círculos
y_mesh = R_mesh .* sin(A_mesh);
mesh(x_mesh, y_mesh, 0 * x_mesh)                     % Dibujar la malla de círculos en el plano xy
axis([-3, 3, -3, 3])      
radio_externo = 2;                             % Parámetros para los cilindros internos y externos
radio_interno = 1;
X_externo = radio_externo .* cos(angulo_malla);             % Coordenadas para el cilindro externo  
Y_externo = radio_externo .* sin(angulo_malla);
X_interno = radio_interno .* cos(angulo_malla);             % Coordenadas para el cilindro interno 
Y_interno = radio_interno .* sin(angulo_malla);
plot(X_externo, Y_externo, 'k', 'LineWidth', 3);                     % Dibujar el cilindro externo
plot(X_interno, Y_interno, 'k', 'LineWidth', 3);                     % Dibujar el cilindro interno
view(2)                                                                % Establecer la vista en 2D
hold off                                             % Desactivar el modo de retención de gráficos

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

2 Calculo de velocidades. Ecuación de Navier-Stokes

Por el enunciado del problema, sabemos que la velocidad de las partículas del fluido viene dada por [math]\vec u(\rho,\theta)=f(\rho)\vec e_\theta[/math]. También sabemos que [math](\vec u, \rho)[/math] tiene que cumplir la ecuación de Navier-Stokes estacionaria:

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

Despreciamos el primer término de la ecuación y como la presión es constante su gradiente es nulo ([math]\triangledown p= \vec 0[/math]) por lo tanto obtenemos:

[math]\mu\triangle\vec u=\vec 0[/math]

Siendo [math]\mu[/math] el coeficiente de viscosidad del fluido y [math]\triangle\vec u[/math] el laplaciano.

Para nuestro ejercicio, utilizaremos la fórmula [math]\triangle\vec u = \triangledown(\triangledown\cdot\vec u)-\triangledown\times(\triangledown\times\vec u)[/math] ya que estamos trabajando en coordenadas cilíndricas.

2.1 Calculo del laplaciano vectorial

Primero calculamos la divergencia de [math]\vec u[/math] utilizando: [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]

Por lo tanto el gradiente de la divergencia también es nulo. ([math]\triangledown(0)=\vec 0[/math]).

Comprobado que el primer término de nuestra fórmula es nulo, pasamos a calcular el rotacional de [math]\vec u[/math]:[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}}= \frac{1}{\rho} \cdot \left [ f(\rho ) + \rho \cdot {f}'\left ( \rho \right )\right ] \bar{e_{z}} [/math]

Ahora calculamos el rotacional del rotacional de [math]\vec u[/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{{f}'\left ( \rho \right )}{\rho }+{f}''\left ( \rho \right )\right ]\vec{e}_{\theta }[/math]

Finalmente, podemos calcular el laplaciano: [math]\ \triangle \vec{u}=-\triangledown \times \triangledown \times \vec{u} =0-\frac{1}{\rho}\left [ \frac{1}{\rho} f\left ( \rho \right )-f'\left ( \rho \right )- \rho f''(\rho)\right ]\vec{e}_{\theta } [/math]

Una vez obtenido todos los términos, sustituimos en la ecuación de Navier-Stokes, obteniendo: [math]\vec{0} + 0 - µ∆\vec{u} = \vec{0} [/math] [math]∆\vec{u} = \vec{0} [/math] [math]\ \triangle \vec{u}= \frac{1}{\rho }\left [ -\frac{1}{\rho }f(\rho)+f'(\rho)+ \rho f''(\rho) \right ] \vec{e}_{\theta } [/math]

2.2 Ecuación diferencial propuesta

A continuación, el problema nos propone comprobar que [math]f(\rho)[/math] satisface la ecuación diferencial: [math]\frac{\partial}{\partial \rho}(\rho\frac{\partial(f(ρ)}{\partial ρ})=\frac{f(\rho)}{\rho} [/math]

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

Podemos afirmar que se cumple. Posteriormente, el enunciado nos pide comprobar si la función [math]f(\rho)=a\rho+\frac{b}{\rho}[/math] es solución. Para ello sustituimos la solución dentro de nuestra ecuación diferencial.

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

Después calculamos [math]f'(\rho)[/math] y [math]f''(\rho)[/math]

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

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

Sustituimos estos datos en [math]f'(\rho)+\rho (f''(\rho))= \frac{f(\rho)}{\rho} [/math] y obtenemos: [math]a-\frac{b}{\rho^2} + \rho \frac{2b}{\rho^3}= \frac{a\rho+\frac{b}{\rho}}{\rho}[/math]. De esta manera, demostramos que es solución.

Para finalizar, el necesitamos hallar los valores de a y b concretos de nuestro ejercicio. Para ello utilizaremos las condiciones de nuestro enunciado, dando los valores de la frontera interior y exterior a [math]\rho[/math] llegando a un sistema de dos ecuaciones.

• En el interior: [math]\ \rho=1\rightarrow \vec{u}=w\vec{e_{\theta }} \rightarrow f(1)\vec{e_{\theta }}=w\vec{e_{\theta }}\rightarrow (a\cdot 1+\frac{b}{1})\vec{e_{\theta }}=w\vec{e_{\theta }}\rightarrow a+b=w [/math]
• En el exterior: [math]\ \rho=2\rightarrow \vec{u}=\vec{0} \rightarrow f(2)\vec{e_{\theta }}=\vec{0} \rightarrow (a\cdot 2+\frac{b}{2})\vec{e_{\theta }}=0\rightarrow 2a+\frac{b}{2}=0 [/math]

[math]\ a=w-b \rightarrow 2(w-b)+\frac{b}{2}=0\rightarrow \frac{3}{2}b=2w [/math]

• [math]\ a=\frac{-1}{3}w [/math]
• [math]\ b=\frac{4}{3}w [/math]

2.3 Condición de incompresibilidad

Sabemos que en los fluidos incompresibles, el gradiente del vector velocidad tiene que ser igual a 0. Como hemos demostrado al inicio del apartado 2.1, nuestro fluido cumple la condición de incompresibilidad dado que [math]∆\vec{u} =\vec{0}[/math] .

3 Representación del campo de velocidades

Suponiendo que, [math] \omega = 1 [/math] y [math] \mu = 1 [/math] obtenemos:

• [math]\ a=\frac{-1}{3} [/math]
• [math]\ b=\frac{4}{3} [/math]

Sustituimos a y b en la ecuación [math]f(\rho)=a\rho+\frac{b}{\rho}[/math] dando lugar a [math]f(\rho)=\frac{-1}{3}\rho+\frac{4}{3\rho}\rightarrow f(\rho)=\frac{-1}{3}(\rho-\frac{4}{\rho})[/math]

Por lo tanto, nuestro campo es:

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

3.1 Representación gráfica

Para representarlo, hemos utilizado el siguiente código de Matlab, dando lugar a lo mostrado en la Figura 2.

clear all % Limpiar la pantalla, espacio de trabajo y habilitar el modo de retención de gráficos
clc
hold on

radio_malla = 1:0.1:2;                        % Definir los parámetros para la malla de círculos        
angulo_malla = 0:0.1:2*pi;           
[R_malla, A_malla] = meshgrid(radio_malla, angulo_malla); 

x = R_malla .* cos(A_malla);                       % Convertir coordenadas polares a cartesianas
y = R_malla .* sin(A_malla);

X = sin(A_malla) .* (4/3 * (R_malla - 1./R_malla));   % Calcular componentes del campo vectorial
Y = cos(A_malla) .* (-1/3 * (R_malla - 1./R_malla));

quiver(x, y, X, Y);                             % Dibujar el campo vectorial usando quiver en 2D

Figura 2: Campo de velocidades.

4 Representación de las líneas de corriente

Para dibujar las líneas de corriente del campo [math]\vec u [/math] tenemos que calcular el campo [math]\vec v [/math] que en cada punto es ortogonal a [math]\vec u [/math] (tomar [math]\vec v =\vec k \times \vec u[/math])

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

4.1 Comprobación de irrotacionalidad

Comprobamos que el campo [math]\vec v[/math] es irrotacional debido a que la divergencia de [math]\vec u[/math] es nula. Aplicando la fórmula del rotacional obtenemos:

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

De esta manera hemos comprobado la irrotacionalidad.

4.2 Calculo de [math]\psi[/math]

A continuación vamos a calcular [math]\psi[/math]:

[math]\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{1}{3}(\rho-\frac{4}{\rho})\vec e_\rho[/math]

Despejamos [math]\psi[/math] y resolvemos la integral:

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

4.3 Representación gráfica

Para realizar la representación de las líneas de corriente del campo (Figura 3), hemos utilizado el siguiente código de Matlab.

clear all                                                                             % Limpiar la pantalla
clc
ptos = 100;                                                       % Definir el número de puntos en la malla

densidad_radial = linspace(1, 2, ptos);  % Generar valores para la densidad radial (rho) en el rango [1, 2]

theta = linspace(0, 2*pi, ptos);               % Generar valores para el ángulo theta en el rango [0, 2*pi]

[RHO, THETA] = meshgrid(densidad_radial, theta);   % Crear una malla 2D de coordenadas polares (rho, theta)

funcion_f = ((RHO.^2/6) - ((4*log(RHO))/3));                 % Definir una función f en coordenadas polares

X = RHO .* cos(THETA);                            % Convertir coordenadas polares a coordenadas cartesianas
Y = RHO .* sin(THETA);

view(2);                                                                        % Establecer la vista en 2D

contour(X, Y, funcion_f, 20);     % Crear un gráfico de contorno de la función f en coordenadas cartesianas
colorbar;                                                                     % Mostrar la barra de colores

axis([-3, 3, -3, 3]);                                                      % Establecer los límites del eje
colormap("cool");                                                          % Establecer la colormap a 'cool'

Figura 3: Líneas de corriente del campo [math] \vec{u} [/math]

5 Velocidad máxima del fluido

Queremos hallar los puntos del fluido donde el módulo de la velocidad sea máxima. Para ello primero calculamos el módulo del campo de velocidades:

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

A continuación derivamos [math]|\vec{u}(\rho)|[/math] respecto [math]\rho[/math] e igualamos a 0:

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

La solución obtenida no es real de manera que al ser continua debemos analizar los extremos de nuestro intervalo ([1,2])

• [math]|\vec{u}(1)|= 1[/math]
• [math]|\vec{u}(2)|= 0[/math]

Por lo tanto nuesto máximo se encuentra cuando [math]\rho = 1[/math], cosa que tiene sentido debido a que en el enunciado nos decían que el exterior se encontraba fijo mientras que el interior se movía a velocidad constante.

5.1 Representación gráfica

Para representar el campo escalar utilizaremos el siguiente código de Matlab, siendo la representación la Figura 4

clear all                           % Limpiar la pantalla, espacio de trabajo
clc
hold on

w = 1;                                                 % Parámetros iniciales
a = -1/3 * w;
b = 4/3 * w;

n = 35;                                                    % Número de puntos

r = linspace(1, 2, n);                                         % Crear mallas
theta = linspace(0, 2*pi, n);
[R, TH] = meshgrid(r, theta);

x = R .* cos(TH);                  % Coordenadas en el espacio tridimensional
y = R .* sin(TH);

f = a * R + b ./ R;   % Función que define la superficie en función de R y TH

colorbar;              % Añade una barra de color para representar la magnitud de los vectores
colormap("cool");

figure;               % Crear y visualizar la superficie tridimensional en 2D
surf(x, y, f);
axis([-3, 3, -3, 3]);
view(2);

Figura 4: Campo escalar del módulo de la velocidad.

6 Representación del rotacional de [math]\vec u[/math]

Como hemos calculado previamente, el rotacional de [math]\vec u[/math] es: [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}}= \frac{1}{\rho} \cdot \left [ f(\rho ) + \rho \cdot {f}'\left ( \rho \right )\right ] \bar{e_{z}} [/math]

También sabemos de apartados anteriores que:

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

Sustituyendo en la fórmula del rotacional: [math]\frac{1}{\rho} \cdot \left [ a\rho+\frac{b}{\rho} + \rho \cdot ( a-\frac{b}{\rho^2}) \right] \bar{e_{z}} = \left [ a+\frac{b}{\rho^2} + a-\frac{b}{\rho^2} \right] \bar{e_{z}} = 2a \cdot \bar{e_{z}} = \frac {-2}{3} \cdot \bar{e_{z}} [/math] [math]\ |\triangledown \times \vec{u}| = \frac {2}{3}[/math]

6.1 Representación gráfica

A continuación mostramos la representación gráfica del campo [math]\ \left | \bigtriangledown \times \vec u \right | [/math] (Figura 5), con el correspondiente código de Matlab utilizado para su obtención.

clear all                                          % Limpiar la pantalla, espacio de trabajo
clc

r_val = 1:0.33:2;                                % Definición de valores para radio y ángulo
th_val = 0:0.33:2*pi;

[R, TH] = meshgrid(r_val, th_val);   % Creación de una malla de puntos en coordenadas polares

x = R .* cos(TH);               % Conversión de coordenadas polares a coordenadas cartesianas
y = R .* sin(TH);
z = zeros(size(x));

r_x = zeros(size(x));                                      % Definición del campo de rotación
r_y = zeros(size(x));
r_z = ones(size(y)) .* (-2/3);

quiver3(x, y, z, r_x, r_y, r_z);         % Representación gráfica del campo de rotación en 3D

axis([-3, 3, -3, 3]);                                  % Configuración de límites en los ejes

xlabel('EJE X');                                                    % Etiquetas para los ejes
ylabel('EJE Y');
zlabel('EJE Z');

Figura 5: Campo vectorial [math]\ \left | \bigtriangledown \times \vec u \right | [/math]

7 Campo de temperaturas

La temperatura del fluido viene dada por el campo [math] T(\rho,\theta) = 1 + \rho^2 sin^2 \theta\cdot e^{−\left ( \rho − \frac{3}{2} \right ) ^2} [/math]

7.1 Representación gráfica

Dibujamos el campo de temperaturas y las curvas de nivel utilizando Matlab.

clear all                                                     % Limpiar la pantalla y el espacio de trabajo 
clc

r = 1:0.1:2;                                                           % Definir las variables alpha y beta
ang = 0:0.1:2*pi;

[R, ANG] = meshgrid(r, ang);     % Generar una retícula rectangular a partir de las divisiones alpha y beta

x_coord = R .* cos(ANG);                                                   % Paso a coordenadas cartesianas
y_coord = R .* sin(ANG);

temp = 1 + R.^2 .* sin(ANG).^2 .* exp(-(R - 3/2).^2);                         % Expresión de la temperatura

figure(1)                                                                     % Representación del campo 3D
mesh(x_coord, y_coord, temp)
axis([-3, 3, -3, 3]);
colormap("cool");
surf(x_coord, y_coord, temp);
colorbar;
xlabel('Coordenada X');
ylabel('Coordenada Y');
zlabel('Temperatura');

figure(2)                                           % Representación de las lineas de nivel del campo en 2D
contour(x_coord, y_coord, temp, 30, 'b');
axis([-3, 3, -3, 3]);
colorbar;
colormap("cool");
xlabel('Coordenada X');
ylabel('Coordenada Y');

Figura 6: Campo de temperaturas 3D

Figura 7: Campo de temperaturas 2D

Figura 8: Curvas de nivel 2D

8 Gradiente de la temperatura

8.1 Cálculo del gradiente

El gradiente de la temperatura en coordenadas cilíndricas queda definido por : [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] Por el apartado 7 sabemos que [math] T(\rho,\theta)[/math] es: [math] T(\rho,\theta) = 1 + \rho^2 sin^2 \theta\cdot e^{−\left ( \rho − \frac{3}{2} \right ) ^2} [/math] Finalmente operando obtenemos: [math] \nabla T(\rho,\theta) = \left [ \rho sen^2 \theta\cdot e^{-\left ( \rho - \frac{3}{2}\right ) ^2} \cdot (-2\rho^2 + 3\rho + 2) \right ] \vec{e}_\rho + 2\rho \sin \theta \cos \theta \cdot e^{-\left(\rho-\frac{3}{2}\right)^{2}} \vec e_\theta = \left [ \rho sen^2 \theta\cdot e^{-\left ( \rho - \frac{3}{2}\right ) ^2} \cdot (-2\rho^2 + 3\rho + 2) \right ] \vec{e}_\rho + \left [ \rho sen(2\theta)\cdot e^{-\left ( \rho - \frac{3}{2}\right ) ^2} \right ] \vec{e}_\theta [/math]

8.2 Representación gráfica

Utilizamos el siguiente código de Matlab para obtener la representación grafica del gradiente de la temperatura, dando lugar a la Figura 9

clear all                                               % Limpiar la pantalla y el espacio de trabajo
clc

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

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

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

figure;                                                                             % Crea una figura

quiver(R, TH, Dx, Dy, 'Color', 'r');  % Dibuja un campo de vectores (quiver plot) con colores cálidos

title('Campo de vectores en coordenadas polares');                      % Añade un título a la figura

xlabel('r');                                                             % Añade etiquetas a los ejes
ylabel('\theta');

colorbar;                     % Añade una barra de color para representar la magnitud de los vectores

Figura 9: Representación del gradiente

Posteriormente, comprobamos que las curvas de nivel de la temperatura y el gradiente son ortogonales. Queda representado en la figura 10.

clear all                                                % Limpiar la pantalla y el espacio de trabajo
clc

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

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

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

figure;                                                                              % Crea una figura

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

figure;                                            % Crea una figura y dibuja las curvas de nivel de Z
contour(R, TH, Z);
title('Curvas de Nivel en Coordenadas Polares');
xlabel('r');
ylabel('\theta');
hold on;

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

quiver(R, TH, Gx, Gy, 'Color', 'r');           % Dibuja un campo de vectores sobre las curvas de nivel

legend('Curvas de Nivel', 'Campo de Vectores');            % Leyenda para indicar el campo de vectores

hold off;

Figura 10: Ortogonalidad del gradiente y las curvas de nivel.

9 Caudal que pasa por una sección longitudinal

El caudal que circula por una superficie queda definido por la siguiente expresión: [math] C = \int\int_S \vec{u}\cdot\vec{dS} [/math]

Nosotros queremos calcular el caudal que circula por la sección longitudinal de los cilindros concéntricos (de altura 1). Para ello, conocemos la velocidad [math] \vec u [/math] dada en [m/s]: [math] \vec u=\frac{-1}{3}(\rho-\frac{4}{\rho})\vec e_\theta[/math]

También conocemos que [math] \vec{dS} = \frac{\partial \vec{r}}{\partial u} \times \frac{\partial \vec{r}}{\partial v}[/math] Siendo [math]\vec{r}(u,v)[/math] una parametrización de la superficie S.

9.1 Parametrización de la superficie

Primero se procede a parametrizar nuestra superficie (Figura 11):

Figura 11: Representación de la sección.

Podemos dividir nuestra superficie S en dos subsuperficies, de manera que [math]S = S_1 + S_2[/math]. De esta manera, la suma de las integrales de cada uno de las subsuperficies será igual a la integral de la superficie S. Hacemos la parametrización:

[math]\left\{\begin{matrix} S_1:\hspace{10pt} \vec{r}_1 = u \vec{e_\rho} + \frac{\pi}{2} \vec{e_\theta} + v\vec{e_v} \hspace{10pt} u \in[1,2]\hspace{5pt} v \in [0,1]\\ S_2:\hspace{10pt} \vec{r}_2 = u \vec{e_\rho} + \frac{3\pi}{2} \vec{e_\theta} + v\vec{e_v}\hspace{10pt} u \in[-2,-1] \hspace{5pt}v \in [0,1] \end{matrix}\right. [/math]

9.2 Cálculo del vector [math] \vec{dS} [/math]

A continuación calculamos el vector [math] \vec{dS} [/math]:

[math]\frac{\partial \vec{r_1}}{\partial u} =\frac{\partial \vec{r_2}}{\partial u} = \vec{e_\rho} [/math] [math] \frac{\partial \vec{r_1}}{\partial v} =\frac{\partial \vec{r_2}}{\partial v} = \vec{e_z} [/math]

[math] \vec{dS_1} = \vec{dS_2} = \frac{\partial \vec{r}}{\partial u} \times \frac{\partial \vec{r}}{\partial v} = \vec{e_\rho} \times \vec{e_z} = -\vec{e_\theta }[/math]

9.3 Cálculo del caudal

Para calcular el caudal primero debemos escribir el vector [math]\vec u[/math] en función de nuestra parametrización:

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

A continuación, resolvemos la integral:

[math] C = \int\int_S \vec{u}\cdot\vec{dS} [/math]

[math] C = \int\int_{S_1} \vec{u}\cdot\vec{dS_1} + \int\int_{S_2} \vec{u}\cdot\vec{dS_2}[/math]

• [math] S_1=\int_{1}^{2}\int_{0}^{1} \frac{-1}{3}(u-\frac{4}{u})\vec e_\theta \cdot (-\vec{e_\theta })du \hspace{1pt} dv = \frac{1}{3}[\frac {u^2}{2}-4log u]^1 _0 = \frac{1}{3}[\frac {1}{2}-4log 1][/math]
• [math] S_2=\int_{-1}^{-2}\int_{0}^{1} \frac{-1}{3}(u-\frac{4}{u})\vec e_\theta \cdot (-\vec{e_\theta })du \hspace{1pt} dv = \frac{1}{3}[\frac {u^2}{2}-4log u]^1 _0 = \frac{1}{3}[\frac {1}{2}-4log 1][/math]

Por lo tanto el caudal final es:

[math]S= S_1+S_2 \rightarrow S= \frac{1}{3}[\frac {1}{2}-4log 1] + \frac{1}{3}[\frac {1}{2}-4log 1] = \frac{2}{3}[\frac {1}{2}-4log 1][/math]