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 e realizar una sección transversal por el plano [math]x_3=0[/math]. Dando como resultado la siguiente sección.

Sección trasversal de los dos cilindros

Para realizar esta sección se ha usado el siguiente programa 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])       % región de vista del dibujo
view(2)                 % Elección de perspectiva

2 Cálculo de las velocidades

2.1 Definición del campo de velocidades

Por el enunciado del problema, sabemos que la velocidad de las partículas viene dada por [math]\vec{u}(ρ, θ) = f(ρ) \vec{e_θ} [/math], además de que su presión (p) es constante. Por otro lado, 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 cilindricas, utilizaremos la siguiente fórmula:

[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. ¿A modo de hipótesis, la divergencia ha de ser nula, ya que el fluido es incompresible?.

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

A la vista de los resultados podemos afirmar que que se trata de un fluido incompresible. Además de esto, 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:

Campo vectorial de las velocidades

La siguiente representación se ha realizado con el siguiente programa 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);
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])      
view(2)                % región de vista del dibujo
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

Además, comprobamos que que el campo [math] \vec v [/math] es irrotacional, por lo que calculamos su 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{1}{3}\ ( \rho-\frac{4}{\rho} ) & 0 & 0\end{vmatrix}=\vec 0[/math]

4.3 Cálculo de las líneas de corriente

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

Luego de calcular [math]\psi [/math], representamos graficamente las líneas de corriente de [math]\vec u [/math]

Lineas 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); % creacion de matrices de los parámetros
xx=uu.*cos(vv);        % parametrization
yy=uu.*sin(vv);
f= 1/3*(-log(uu)+(uu.^2)./2); % función de las lineas de corriente
axis([-3,3,-3,3])      % región de vista del dibujo
view(2)
contour(xx,yy,f,20);   % dibujo de la lineas de corriente

5 Velocidad máxima del fluido

En este apartado, identificaremos los puntos donde la velocidad del fluido es máxima. Para ello, primero calcularemos el módulo del campo vectorial que corresponde a la velocidad:

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

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]

Como dicha ecuación 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.

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

7 Cálculo del 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:

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.

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]

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]

8 Cálculo del gradiente de la temperatura

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

[math]\triangledown\T = \ltmath\gt\frac{ \partial\T}{\partial \rho}\left\vec{e_\rho} + \frac{1}{\rho}frac{ \partial\T}{\partial \theta}\left\vec{e_\theta} + \frac{ \partial\T}{\partial \z}\left\vec{e_\z}\lt\math\gt Como nuestra función de temperatura T solo depende de \ltmath\gt \rho \lt\math\gt y \ltmath\gt \theta \lt\math\gt \ltcenter\gt\ltmath\gt T(\rho,\theta) = 1 + \rho^{2}sen^{2}\theta e^{ - (\rho - frac{3}{2})^{2}}\lt\math\gt Por tanto, el gradiente de la temperatura será: \ltcenter\gt\ltmath\gt\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}}]]]\left\vec{e_\rho} + [sen2\theta\rho e^{ - (\rho - frac{3}{2})^{2}}]\left\vec{e_\theta} [[Categoría:Teoría de Campos]] [[Categoría:TC23/24]][/math]