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

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos tubos concéntricos (Grupo B9)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Ignacio Calmet, Roddyk Inocencio, Joel Huanca, Alejandro Molina, Daniel Zúñiga
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Flujo de Couette es un tipo de flujo que ocurre en mecánica de fluidos entre dos paredes paralelas, una en movimiento y otra con un movimiento estacionario. Este tipo de flujo se utiliza para estudiar la viscosidad y el comportamiento de corte de los fluidos. En nuestro caso 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 del cilindro interior es ω > 0

1 Representación de la sección transversal

Con la finalidad de ver los puntos ocupados por el fluido, vamos a realizar una sección transversal por el plano [math]x_3=0[/math]. Dando como resultado la siguiente sección:

centro

Código en Matlab:

h=0.08;                 % intervalo para dividir los parámetros 
u=1:h:2;                % intervalos de u [1,2]
v=0:h:2*pi;             % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v);  % creación de matrices de los parámetros
figure(1)
xx=uu.*cos(vv);         % parametrización
yy=uu.*sin(vv);
mesh(xx,yy,0*xx);       % dibujo de la figura con la coordenadas parametrizadas
axis([-3,3,-3,3])       % Ejes del dibujo
view(2)                 % Elección de perspectiva

2 Cálculo de las velocidades

2.1 Definición del campo de velocidades

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

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

2.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 como el campo de velocidades está dado en coordenadas cilíndricas utilizaremos:

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

2.2.1 Gradiente de la divergencia

En cuanto al 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]

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

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

2.3 Cálculos con de 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]

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

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

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

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

2.3.3 Búsqueda de las constantes "a" y "b"

Una vez confirmado la solución, para poder hallar las constantes "a" y "b" debemos imponer las condiciones dadas:

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]

Como resultado obtenemos un sistema de dos ecuaciones con dos incógnitas, es este caso lo resolveremos por el método de sustitución, despejando en la segunda ecuación tenemos que [math] (b=-4a) [/math]. Gracias a esto tenemos las dos constantes:?

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

Quedando la ecuación diferencial de siguiente forma:

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

3 Representación gráfica del campo de velocidades

Imponiendo las condiciones del enunciado ([math] \omega = 1 [/math] y [math] \mu = 1 [/math]) tenemos: [math] a = -\frac{1}{3}\omega [/math] [math] b = \frac{4}{3}\omega[/math]

Sustituyendo en la función obtenemos:

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

Por lo tanto tenemos el siguiente campo vectorial:

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

Su representación sería la siguiente:

Figura 2. Campo vectorial de las velocidades

Código en Matlab:

h=0.1;                 % intervalo para dividir los parámetros 
u=1:h:2;               % intervalos de u [1,2]
v=0:h:2*pi;            % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v); % creacion de matrices de los parámetros
xx=uu.*cos(vv);        % parametrización
yy=uu.*sin(vv);
dx=-sin(vv).*(1/3*(4./uu-uu)); % dirección de la velocidades 
dy=cos(vv).*(1/3*(4./uu-uu)); 
axis([-3,3,-3,3])      % Ejes del dibujo
view(2)                % Elección de perspectiva
quiver(xx,yy,dx,dy);   % dibujo del campo vectorial

4 Líneas de Corriente

4.1 Vector ortogonal a [math] \vec u [/math]

Queremos dibujar las líneas de corriente del campo vectorial [math] \vec u [/math], es decir las líneas tangentes en cada punto. Para ello calculamos el vector [math] \vec v [/math] que es el ortogonal a [math] \vec u [/math] en cada punto.

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

4.2 Demostración de que el campo [math] \vec v [/math] es irrotacional

Como [math] \vec u [/math] es de divergencia nula, el campo [math] \vec v [/math] es irrotacional.

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

4.3 Cálculo de las líneas de corriente

Por último, calculamos el potencial escalar [math]\psi [/math] de [math] \vec v [/math], siendo la función de corriente de [math] \vec u [/math].

[math]\triangledown\psi= \vec v [/math], [math]\psi= f(\rho, \theta)[/math], [math]\triangledown\psi= \frac{\partial \psi}{\partial \rho}\vec e_\rho = \frac{1}{3}\ ( \rho-\frac{4}{\rho} ) \vec{e}_\rho[/math]

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

4.3.1 Representación grafica de [math]\psi [/math]

Figura 3. Líneas de corrientes

Codigo de matlab:

h=0.1;                 % Intervalo para dividir los parámetros 
u=1:h:2;               % Intervalos de u [1,2]
v=0:h:2*pi+h;          % Intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v); % Creación de matrices de los parámetros
xx=uu.*cos(vv);        % Parametrización
yy=uu.*sin(vv);
f= 1/3*(-4log(uu)+(uu.^2)./2); % Función de las líneas de corriente
axis([-3,3,-3,3])      % Región de vista del dibujo
view(2)                % Elección de perspectiva
contour(xx,yy,f,20);   % Dibujo de la líneas de corriente
colorbar;

5 Velocidad máxima del fluido

5.1 Módulo del campo de velocidades

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

5.2 Puntos críticos

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

5.3 Valor máximo de[math]\left| \vec{u}(\rho) \right|[/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 [math] \rho\in(1,2) [/math]:

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

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

El módulo de la velocidad es máximo cuando [math]\rho=2[/math], es decir, en los puntos del tubo interior que gira con velocidad angular constante.

5.4 Representación gráfica

Representamos en una gráfica el comportamiento del módulo de la velocidad en función de [math]\rho[/math]:

Figura 4. Campo escalar: Módulo de la velocidad

Codigo de matlab:

h=0.1;                 % intervalo para dividir los parámetros
u=1:h:2;               % intervalos de u [1,2]
v=0:h:2*pi+h;          % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v); % Creación de matrices de los parámetros
xx=uu.*cos(vv);        % parametrization
yy=uu.*sin(vv);
f= 1/3*(4./uu-uu);     % módulo de campo de velociades
surf(xx,yy,f);         % Dibujo del campo escalar
axis([-3,3,-3,3])      % Ejes del dibujo
view(2)                % Elección de perpectiva

6 Rotacional

6.1 Cálculo del rotacional de [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}{ρ}\frac{ \partial (ρf(ρ)) }{\partial ρ}\vec{e_z} = [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}]\vec{e_z}[/math] De apartados anteriores sabemos que [math] f(\rho) = \frac{1}{3}\left ( \frac{4}{\rho}-\rho \right ) [/math]

Por lo que resulta

[math](\nabla\times\vec{u}) = -\frac{2}{3}\vec{e_z} [/math]

El rotacional es constante en todos los puntos, esto es razonable debido a que el rotacional mide el giro de un cuerpo ortogonal a [math]\vec{u}[/math]sobre sí mismo cuando flota sobre un fluido.

6.2 Representación gráfica

A continuación se representa el campo [math]\ \left | \bigtriangledown \times \vec u \right | = \frac{2}{3} [/math] :

Figura 5. Campo escalar: Modulo del rotacional

Código de Matlab:

h=0.1;                 % Intervalo para dividir los parámetros 
u=1:h:2;               % Intervalos de u [1,2]
v=0:h:2*pi;            % Intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v); % Creacion de matrices de los parámetros
xx=uu.*cos(vv);        % Parametrización
yy=uu.*sin(vv);
z=0.*uu;
fz=ones(size(vv));     
zz=fz*(2/3);           % valor del modulo del rotacional
surf(xx,yy,zz)
axis([-3,3,-3,3])      %  ejes del dibujo
view(2)                % Prespectiva del dibujo

7 Campo de temperaturas

La temperatura del fluido viene definida por el siguiente 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

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

Figura 6. Campo de temperaturas y curvas de nivel

Código de Matlab:

h=0.1;                    % intervalo para dividir los parámetros 
u=1:h:2;                  % intervalos de u [1,2]
v=0:h:2*pi+h;             % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v);    % creacion de matrices de los parámetros
T= 1+uu.^2.*(sin(vv.*exp((-(uu-3/2).^2)))).^2;
subplot(1,2,1);           % Ventanas 
surf(uu,vv,T);            % Dibujo de la superficie 
shading flat              % Difuminado
subplot(1,2,2)            % Ventana 2
pcolor(uu,vv,T);          % Proyección en planta
hold on                   % Mantener ventana
contour(uu,vv,T,10,'k');  % 10 curvas de nivel
hold off

7.2 Punto máximo de temperatura

En esto caso vamos a calcular el punto máximo gráficamente, por lo tanto, su representación gráfica es la siguiente:

Figura 7. Gráfico de la temperatura

Código en Matlab

h=0.1;                    % intervalo para dividir los parámetros 
u=1:h:2;                  % intervalos de u [1,2]
v=0:h:2*pi+h;             % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v);    % creacion de matrices de los parámetros
T= 1+uu.^2.*(sin(vv.*exp((-(uu-3/2).^2)))).^2;
plot3(uu,vv,T);           % dibujo de las lineas de campo

8 Gradiente de la temperatura

8.1 Cálculo del gradiente

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]

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

[math] T(\rho,\theta) = 1 + \rho^{2}sen^{2}\theta e^{ - (\rho - \frac{3}{2})^{2}}[/math]

Por tanto, el gradiente de la temperatura será:

[math]\triangledown T = [sen^{2}\theta[2\rho e^{ - (\rho - \frac{3}{2})^{2}} + \rho^{2}[ - 2(\rho - \frac{3}{2}) e^{ - (\rho - \frac{3}{2})^{2}}]]]\vec{e_\rho} + [sen(2\theta)\rho e^{ - (\rho - \frac{3}{2})^{2}}]\vec{e_\theta}[/math]

8.2 Representación gráfica del gradiente

Figura 8. Representación del gradiente de la temperatura

Codigo de matlab:

h=0.1;                     % intervalo para dividir los parámetros 
u=1:h:2;                   % intervalos de u [1,2]
v=0:h:2*pi+h;              % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v);     % creacion de matrices de los parámetros
T= 1+uu.^2.*(sin(vv.*exp((-(uu-3/2).^2)))).^2;
[dx,dy] = gradient(T,u,v); % Gradiente de la temperatura
figure;
quiver(uu,vv,dx,dy);       % Dibujo del campo vectorial

8.3 Comprobación de la ortogonalidad del gradiente

Esta comprobación se ha realizado gráficamente con la siguiente figura:

Figura 9. Ortogonalidad gradiente y curvas de nivel

h=0.1;                     % intervalo para dividir los parámetros 
u=1:h:2;                   % intervalos de u [1,2]
v=0:h:2*pi+h;              % intervalo de v [0,2*pi]
[uu,vv]=meshgrid(u,v);     % creacion de matrices de los parámetros
T= 1+uu.^2.*(sin(vv.*exp((-(uu-3/2).^2)))).^2;
[dx,dy] = gradient(T,u,v); % Gradiente de la temperatura
figure;
contour(uu,vv,T)           % curvas de nivel 
hold on
quiver(uu,vv,dx,dy);       % gradiente
hold off;

9 Caudal

Queremos conocer el caudal que atraviesa la sección longitudinal de los cilindros. Estos tienen una altura de 1 metro y la velocidad del fluido lo medimos en m/s. La velocidad de los puntos del fluido está determinada por el campo vectorial:

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

La sección a considerar se trata de la composición de dos superficies [math]S_1[/math] y [math]S_2[/math], tal y como representa la siguiente figura:

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

Para calcular el caudal, debemos calcular la siguiente integral:

[math]\int_{S}^{} \vec{u} \cdot d\vec{S} = \iint_{D}^{} \vec{u} ( \Phi (u,v))\cdot ( \Phi_u \times \Phi_v) du dv [/math]

[math]S[/math] es la composición de las superficies [math]S_1[/math] y [math]S_2[/math], [math]D[/math] es el dominio de los valores que adoptan [math]u[/math] y [math]v[/math], y escogeremos el vector [math]\vec{e_\theta}[/math] como sentido de orientación de las superficies.

9.1 parametrización de la superficie

El primer paso es determinar la parametrización de las superficies.

[math]\ S_1 [/math]: [math]\ \Phi (u,v) = (\rho (u,v),\theta (u,v),z(u,v)) = (u,\frac{\pi}{2},v) [/math] donde [math]\ \left\{\begin{matrix} u\in (1,2) \\ v\in (0,1) \end{matrix}\right. [/math]

[math]\ S_2 [/math]: [math]\ \Phi (u,v) = (\rho (u,v),\theta (u,v),z(u,v)) = (u,\frac{3\pi}{2},v) [/math] donde [math]\ \left\{\begin{matrix} u\in (1,2) \\ v\in (0,1) \end{matrix}\right. [/math]

Calculamos los vectores velocidad [math]\Phi_u[/math] y [math]\Phi_v[/math] para posteriormente obtener el vector normal a la superficie [math] \Phi_u \times \Phi_v [/math]. Primero con la superficie [math]\ S_1 [/math]:

[math]\Phi_u = \frac{ \partial}{\partial u} u \vec{e_\rho} + \rho\frac{ \partial}{\partial u}\frac{\pi}{2}\vec{e_\theta} + \frac{ \partial}{\partial u} v \vec{e_z} = \vec{e_\rho} [/math]

[math]\Phi_v = \frac{ \partial}{\partial v} u \vec{e_\rho} + \rho\frac{ \partial}{\partial v}\frac{\pi}{2}\vec{e_\theta} + \frac{ \partial}{\partial v} v \vec{e_z} = \vec{e_z} [/math]

[math] \Phi_u \times \Phi_v = \vec{e_\rho} \times \vec{e_z} = -\vec{e_\theta}[/math]

Como podemos observar en las operaciones, el valor de [math] \theta = \frac{\pi}{2} [/math] no influye en el cálculo del vector normal, ya que se trata de un valor fijo que, al derivarlo, se anula. Si realizáramos los mismos cálculos con la superficie [math]\ S_2 [/math], obtendríamos los mismos valores, puesto que el único valor en el que difieren en la parametrización es en [math] \theta [/math], que también es un valor fijo.

El vector normal para ambas superficies es [math] \Phi_u \times \Phi_v = -\vec{e_\theta} [/math], que mira hacia el lado opuesto de la orientación que hemos escogido ([math]\vec{e_\theta}[/math]). Por lo tanto, cuando calculemos el caudal, cambiaremos de signo al integral.

El vector [math]\vec{u} ( \Phi (u,v))[/math], tanto para [math]\ S_1 [/math] como a [math]\ S_2 [/math], será:

[math]\vec{u} ( \Phi (u,v)) = \vec{u} (u,\frac{\pi}{2},v) = \vec{u} (u,\frac{3\pi}{2},v) =\frac{1}{3}(\frac{4}{u}-u)\vec{e_\theta} [/math]

9.2 Cálculo del caudal

Finalmente, podemos calcular la integral para obtener el caudal:

[math]\ Caudal = \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,\frac{\pi}{2},v)\cdot (-\vec e_\theta)dudv - \int_{0}^{1}\int_{1}^{2} \vec u (u,\frac{3\pi}{2},v)\cdot (-\vec e_\theta)dudv = [/math] [math]\ = -2 \int_{0}^{1}\int_{1}^{2} (\frac{1}{3}(\frac{4}{u}-u)\vec{e_\theta}) \cdot (-\vec e_\theta)dudv = \frac{2}{3} \int_{0}^{1}\int_{1}^{2} (\frac{4}{u}-u)dudv = [/math] [math]\ = \frac{2}{3} \int_{1}^{2} (\frac{4}{u}-u)du = \frac{2}{3}[(4\ln{u} - \frac{u^{2}}{2}]^2 _1 = (\frac{8}{3} \ln{2} - 1) \frac{m^{3}}{s}[/math]