Flujo de Couette entre dos tubos concéntricos (Grupo 8-C)

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos. Grupo 8-C
Asignatura	Teoría de Campos
Curso	2022-23
Autores	Gabriel Moreno Pardo Daniel Alexandre Ferreira Patricio Francisco Javier Vela Cobos Juan Carlos Fernández Alonso
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 INTRODUCCIÓN

En este proyecto tenemos un fluido dado, el cual es incompresible, tiene la capacidad de oponerse a la compresión de sí mismo bajo cualquier condición. Se nos ha dado también un obstáculo, el cual son dos cilindros concéntricos, de manera que el exterior se mueve con velocidad angular constante en sentido antihorario mientras que el interior está fijo. Tenemos un radio distinto para cada cilindro (1 y 2). En la primera imagen adjuntada (Figura 1) podemos observar una representación del mallado del flujo de Couette entre los cilindros antes comentados, con su código de lenguaje m (MATLab).

clear, close all               
rho=1:0.1:2;              
theta=0:0.1:2*pi;           
[RHO,THETA]=meshgrid(rho,theta); 
figure(1)
x=RHO.*cos(THETA);        
y=RHO.*sin(THETA);
mesh(x,y,0*x)          
axis([-3,3,-3,3])      
view(2)

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

2 ECUACIÓN DE NAVIER-STOKES

El enunciado nos proporciona el campo vectorial de la velocidad de las partículas que recorren el canal [math]\vec u(\rho,\theta)=f(\rho)\vec e_\theta[/math] Sabemos que se cumple la ecuación de Navier-Stokes estacionaria: [math](\vec u \cdot \triangledown) \vec u + \triangledown p= \mu \triangle \vec u [/math]

Debido a que la presión es constante: [math]\triangledown p= \vec 0[/math]

Sabiendo esto, y que se desprecia la primera parte de la ecuación estacionaria, tenemos la siguiente igualdad:

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

Para calcular el laplaciano, y ahorrarnos el cambio de base a la base cartesiana, hemos decidido aplicar la siguiente fórmula: [math]\triangle\vec u = \triangledown(\triangledown\cdot\vec u)-\triangledown\times(\triangledown\times\vec u)[/math]

Aparte de la razón anterior, también hemos querido calcular el laplaciano de esta forma porque en esta fórmula se incluyen cálculos como la divergencia y el rotacional, los cuales son requeridos más adelante en el problema.

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]

Con esto también comprobamos la condición de incompresibilidad, al ser la divergencia igual a 0.

Como es conocido, el gradiente de 0 es igual a 0, así que:

[math]\triangledown(\triangledown\cdot\vec u)=\vec 0[/math]

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 } = -\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]

Con el rotacional del rotacional, tendríamos todos los cálculos necesarios para obtener el laplaciano, y por lo tanto la ecuación de Navier Stokes.

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

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

La ecuación sería la siguiente:

[math]\mu \triangle \vec{u}= 0 \rightarrow \frac{1}{\rho}\left [ -\frac{1}{\rho}f(\rho)+f'(\rho)+ \rho f''(\rho) \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]

Vamos a comprobar si [math]f(\rho)[/math] 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.

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

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

Valores de a y b:

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, e igualar la función a las velocidades que tienen los cilindros exterior(w) e interior(0):

[math]\ \rho=1\rightarrow \vec{u}=\vec{0} \rightarrow f(1)\vec{e_{\theta }}=\vec{0}\rightarrow (a\cdot 1+\frac{b}{1})\vec{e_{\theta }}=0\rightarrow a+b=0 [/math]

[math]\ \rho=2\rightarrow \vec{u}=w\vec{e_{\theta }} \rightarrow f(2)\vec{e_{\theta }}=w\vec{e_{\theta }} \rightarrow (a\cdot 2+\frac{b}{2})\vec{e_{\theta }}=w\vec{e_{\theta }}\rightarrow 2a+\frac{b}{2}=w [/math]

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

Con lo que tenemos finalmente: [math]\ a=\frac{2}{3}w [/math] y [math]\ b=-\frac{2}{3}w [/math]

3 CAMPO DE VELOCIDADES

Tenemos [math]\ \omega =1 [/math] y [math]\ \mu=1 [/math]

Campo de velocidades:[math]\ \vec{u}=f(\rho)\vec e_{\theta}= \frac{2}{3}w\rho -\frac{2w}{3\rho } \vec{e_{\theta }} = (\frac{2}{3}\rho-\frac{2}{3\rho})\vec{e_{\theta }} [/math]

Figura 2: Campo de velocidades

clear, close all                 
rho=1:0.1:2;              
theta=0:0.1:2*pi;           
[RHO,THETA]=meshgrid(rho,theta); 
figure(1)
x=RHO.*cos(THETA);        
y=RHO.*sin(THETA);
xx=sin(THETA).*(-2/3*(RHO-1./RHO));
yy=cos(THETA).*(2/3*(RHO-1./RHO));
quiver(x,y,xx,yy);

4 LÍNEAS DE CORRIENTE

Vamos a calcular las líneas de corriente del campo [math] \vec u [/math] ,es decir, las líneas que son tangentes a [math] \vec u [/math] en cada punto. Se empieza calculando 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{2}{3}\rho-\frac{2}{3\rho})\vec{e_{\theta }}=(\frac{2}{3\rho}-\frac{2\rho}{3})\vec e_\rho[/math]

A pesar de que se menciona en el enunciado, vamos a comprobar si el campo es irrotacional, por lo que tenemos que calcular el rotacional de [math]\vec v[/math]

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

Vamos a calcular [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{2}{3\rho}-\frac{2\rho}{3})\vec e_\rho[/math]

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

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

Una vez calculado el potencial escalar, vamos a representar gráficamente las líneas de corriente:

Figura 3: Líneas de corrientes

h=80;                       % Número de intervalos
rho=linspace(1,2,h);           % Parametro u [1,2]
theta=linspace(0,2*pi,h);        % Parametro  v [0,2*pi]
[RHO,THETA]=meshgrid(rho,theta);         % Matrices de coordenadas de U y V
f=((2*log(RHO)/3)-((2*RHO)/3));      % Campo escalar
X=RHO.*cos(THETA);                 % Parametrización
Y=RHO.*sin(THETA);
axis([-3,3,-3,3])            % Ejes
view(2)                     
contour(X,Y,f,20);colorbar            % Líneas de nivel
axis([-3,3,-3,3])            % Ejes del dibujo
colormap("cool")

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 que se demuestra que el máximo es alcanzado en [math]\ \rho = 2 [/math] y su valor es de 1

clear, close all

%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);

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

%Módulo de la velocidad
Modulo=2/3.*RHO-2/3./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');

Figura 4: Variación del módulo de la velocidad en función de [math]\ \rho [/math]

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

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

%Se definen las variables rho y theta
rho=1:0.33:2;
theta=0:0.33:2*pi;

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

%Paso a cartesianas
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
z=zeros(size(x));

% rotacional
xx=zeros(size(x));
yy=zeros(size(x));
zz=ones(size(y)).*4/3;

%Representación del módulo del rotacional
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');

Figura 5: Representación del rotacional

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 CAMPO DE TEMPERATURAS

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]

A continuación, se va a representar el campo escalar con la ayuda de matlab y se estudiará el gradiente de la temperatura.

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, close all
%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);
%Paso a cartesianas
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
%La temperatura del fluido viene dada por la expresión
z=1+RHO.^2.*sin(THETA).^2.*exp(-(RHO-3/2).^2);

%Representación del campo en 3D
figure(1)
mesh(x,y,z)
axis([-3,3,-3,3]);
surf(x,y,z);
colorbar;
title('Representación de la temperatura en 3D');
xlabel('Eje X1');
ylabel('Eje X2');
zlabel('Eje X3');


%Representación del campo en 2D
figure(2)
mesh(x,y,z)
axis([-3,3,-3,3]);
surf(x,y,z);
view(2);
colorbar;
title('Representación de la temperatura en 2D');
xlabel('Eje X1');
ylabel('Eje X2');

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

Figura 7: Representación de la temperatura en 2D.

Como se demuestra en las representaciones, los máximos de la temperatura se consiguen en [math]\ \rho = 2 [/math] y en [math]\ \theta = 0 [/math] y [math]\ \theta = \pi [/math].

Para dibujar las curvas de nivel se utiliza el comando contour con 15 líneas.

clear, close all

%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);

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

%La temperatura del fluido viene dada por la expresión
T=1+RHO.^2.*sin(THETA).^2.*exp(-(RHO-3/2).^2);

%Representación de las curvas de nivel
contour(x,y,T,15,'r');
title('Curvas de nivel');
xlabel('Eje X1');
ylabel('Eje X2');

Figura 8: Representación de las curvas de nivel.

7.2 Gradiente de la temperatura

Se define en coordenadas cilíndricas el campo vectorial gradiente de un campo escalar como:

[math]\ \bigtriangledown f(\rho ,\theta ,z) = \frac{\partial f}{\partial \rho }\vec e_\rho + \frac{1}{\rho}\frac{\partial f}{\partial \theta}\vec e_\theta + \frac{\partial f}{\partial z} \vec e_z [/math]

Si se aplica esta fórmula a nuestro campo de temperaturas se tiene [math]\ \bigtriangledown T [/math] :

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

clear, close all

%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);

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

%La temperatura viene dada por:
T=1+RHO.^2.*sin(THETA).^2.*exp(-(RHO-3/2).^2);

%Las componentes del gradiente son:
xx=sin(THETA).^2.*exp(-(RHO-3/2).^2).*(2.*RHO-2.*RHO.^2.*(RHO-3/2));
yy=2.*RHO.*sin(THETA).*cos(THETA).*exp(-(RHO-3/2).^2);

%Representación del gradiente
hold on
quiver(x,y,xx,yy)
contour(x,y,T,15,'r');
axis([-3,3,-3,3]);
title('Representación del gradiente de la temperatura');
xlabel('Eje X1');
ylabel('Eje X2');
hold off

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

8 CAUDAL

En el siguiente apartado se va a calcular el caudal que pasa por la sección longitudinal de los cilindros concéntricos (intersección del plano [math]\ x_2=0 [/math] con los cilindros).

Para ello, se supone que los cilindros tienen una profundidad de 1m y que la velocidad del fluido dada en [m/s] se corresponde con el campo:

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

De esta forma, el caudal que atraviesa la sección viene dado por la integral:

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

siendo [math]\ S [/math] la superficie que atraviesa el fluido (sección longitudinal en este caso) y [math]\ D [/math] el dominio de los parámetros [math]\ u [/math] y [math]\ v [/math] .

La superficie se va a suponer con orientación en la dirección de [math]\ \vec e_\theta [/math] , para obtener un caudal positivo.

En primer lugar, se va a parametrizar la superficie:

Figura 10: Representación de la sección longitudinal

Como se aprecia en la figura anterior, la sección resulta en dos superficies rectangulares, [math]\ S_1 [/math] y [math]\ S_2 [/math]. La superficie total será la suma de ambas, [math]\ S = S_1 + S_2[/math] La parametrización de [math]\ S_1 [/math] y [math]\ S_2 [/math] es la siguiente:

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

• vectores velocidad:

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

La orientación es "mala", es contraria a la supuesta en la superficie. Para corregir esto, después al calcular la integral de la superficie se le cambia el signo y se le pone negativo.

Estos vectores se han calculado 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]\ 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]

• vectores velocidad:

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

La orientación es "mala", es contraria a la supuesta en la superficie. Para corregir esto, después al calcular la integral de la superficie se le cambia el signo y se le pone negativo.

Una vez definimos las parametrizaciones se calcula la integral.

La integral sobre la superficie [math]\ S [/math] es la suma de las integrales sobre 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{2}{3} \int_{0}^{1}\int_{1}^{2}\left ( u-\frac{1}{u}\right ) \vec e_\theta(-\vec e_\theta)dudv - \frac{2}{3} \int_{0}^{1}\int_{1}^{2}\left ( u-\frac{1}{u}\right ) \vec e_\theta(-\vec e_\theta)dudv = [/math]

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

Como comentario, si en vez de haber cogido [math]\ S_1 [/math] y [math]\ S_2 [/math] como la superficie total, se hubiera tomado una de las dos el resultado del caudal sería [math]\ 1 - \frac{2}{3} \ln 2 \left [m^3/s \right ] [/math].