Flujo de Couette. Grupo 8. Planos paralelos horizontales

Trabajo realizado por estudiantes
Título	Flujo de Couette entre dos planos paralelos horizontales (Grupo 8)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Carlos Reillo González, Carlos Sayalero Gómez, Daniel Portincasa Navarro, Diego Vásquez Peñaloza, Irene Blanes Grau
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El objetivo principal de este trabajo es analizar el Flujo de Couette, el cual podemos definir como un tipo de flujo viscoso confinado entre dos planos paralelos horizontales, en el que uno de los planos está en movimiento constante mientras que el otro permanece estacionario. El fluido se encuentra entre estos dos planos y se mueve en respuesta al movimiento relativo de los planos. En nuestro caso vamos a considerar un fluido incompresible (aquel fluido cuya densidad no varía significativamente bajo la influencia de cambios de presión), de manera que el plano inferior es el que se mueve con velocidad constante v[math]\vec{j}[/math]; v>0. Trabajaremos en el plano con coordenadas cartesianas (y,z). Nos hemos apoyado en el programa informático MATLAB para poder comprender y visualizar el comportamiento del fluido.

2 Mallado del fluido

En la siguiente figura representamos la superficie que ocupa el fluido en los puntos interiores del rectángulo de dimensiones [0,8] x [0,1], y con ejes (y,z) definidos en los intervalos [0,8] x [-1,2]:

Mallado

Para la obtención de esta gráfica hemos usado MATLAB, con un paso de muestreo de 0.1. Adjunto el código utilizado:

y=0:0.1:8;       
z=0:0.1:1;            
[yy,zz]=meshgrid(y,z); 
figure(1)
mesh(yy,zz,0*yy)       
axis([0,8,-1,2])     
view(2)

3 Ecuación de Navier-Stokes estacionaria

La ecuación de Navier-Stokes estacionaria es fundamental en la mecánica de fluidos, describe el comportamiento de fluidos estables o en movimiento constante. Se elimina la dependencia temporal y nos resulta útil en los sistemas en equilibrio. Para nuestro caso en particular la vamos a aplicar para modelar el flujo de fluidos viscosos entre dos placas sólidas, proporcionando detalles sobre la distribución de velocidad y presión en el sistema.

Sabiendo que la ecuación estacionaria de Navier-Stokes es la siguiente:

[math](\vec{u} \cdot \nabla) \vec{u} + \nabla p=μ \Delta\vec{u} [/math]

Y siendo la representación de un fluido viscoso en sus componentes. Su descripción se hará de la siguiente forma:

[math](\vec{u} \cdot \nabla) \vec{u}=u_j\frac{\partial u_i}{\partial x_j}\vec{e_i} [/math] ; [math]\Delta\vec {u}=u_i\vec{e_i}[/math] ; [math]\vec{u}=u_i\vec{e_i}[/math]

Sabiendo que [math](\vec{u},p)[/math], satisface la ecuación de Navier-Stokes y que en este caso en concreto, trabajaremos con: [math]\vec{u}(y,z)=f(z)\vec{j} [/math], que es el campo vectorial de la velocidad de las partículas del fluido ; [math]\ p(x,y,z)=p_1+(p_2-p_1)(y-1)[/math] es el campo de presiones del fluido y [math]μ[/math] es la viscosidad del fluido. Calcularemos la ecuación diferencial que satisfaga [math] f(z)[/math].

1.) Vamos a resolver cada término de la ecuación por separado:

[math](\vec{u} \cdot \nabla) \vec{u}=\begin{pmatrix}{\frac{\partial u_1}{\partial y}}&{\frac{\partial u_1}{\partial z}}\\{\frac{\partial u_2}{\partial y}}&{\frac{\partial u_2}{\partial z}}\end{pmatrix}\cdot \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} [/math]

y sustituyendo con nuestros datos:

[math](\vec{u} \cdot \nabla) \vec{u}=\begin{pmatrix} 0 & f'(z) \\ 0 & 0 \end{pmatrix} \cdot \begin{pmatrix} f(z)\\ 0 \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \end{pmatrix} [/math]

2.) Gradiente del campo de presiones:

[math]\nabla p=\begin{pmatrix} \frac{\partial p}{\partial y}\\ \frac{\partial p}{\partial z} \end{pmatrix}=\begin{pmatrix} p_2-p_1 \\ 0 \end{pmatrix}= (p2-p1)\vec{j}[/math]

3.) Laplaciano [math]\Delta\vec{u} [/math] , del campo de velocidades. El cual calculamos con la siguiente fórmula:

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

Calculamos la divergergencia:

[math]\nabla \cdot \vec{u}=\frac{\partial u_1}{\partial y}+\frac{\partial u_2}{\partial z} =0[/math]

Así [math]\nabla(\nabla \cdot \vec{u})=0[/math]

La divergencia de un campo vectorial [math]\vec{u}[/math] describe el comportamiento del fluido en diferentes regiones. Divergencia positiva [math](\nabla \cdot \vec{u})\gt0[/math] implica expansión local, mientras que divergencia negativa [math](\nabla \cdot \vec{u})\lt0[/math] indica contracción. Con divergencia cero [math](\nabla \cdot \vec{u})=0[/math], el fluido es incompresible y conserva la masa. Estos conceptos son esenciales en la dinámica de fluidos y se aplican en el modelado fluidodinámico, donde la divergencia se utiliza para analizar flujos y evaluar la variación local de la densidad del fluido.

• Calculamos el rotacional:

[math]\nabla\times\vec u= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ 0 & f(z) & 0 \end{vmatrix}=-f'(z)\vec{i}[/math]

• Y finalmente calculamos [math]\nabla\times(\nabla\times\vec u)[/math]:

[math]\nabla\times(\nabla\times\vec u)= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ -f'(z) & 0 & 0 \end{vmatrix}= -f''(z)\vec{j}[/math]

Sustituyendo en [math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math], obtenemos:

[math]\Delta\vec{u}=f''(z)\vec{j}[/math]

Sustituyendo en la ecuación de Navier-Stokes obtenemos:

[math]\begin{pmatrix} 0 \\ 0 \end{pmatrix}+\begin{pmatrix} p_2-p_1 \\ 0 \end{pmatrix}=μ\begin{pmatrix} f''(z)\\ 0 \end{pmatrix}[/math]

[math] p2-p1=μ\cdot f''(z)[/math]

[math] f''(z)= \frac{(p_2-p_1)}{μ}\vec{j} [/math]

Teniendo ya calculada la segunda derivada del campo vectorial de velocidades, podemos calcular [math]f'(z)[/math] y [math]f(z)[/math] facilmente integrando. Obteniendo:

[math] f'(z)= \frac{(p_2-p_1)\cdot z}{μ} [/math]

[math] f(z)= \frac{(p_2-p_1)\cdot z^2}{2\cdotμ}[/math]

Suponemos que la velocidad del fluido en [math]y=0,1[/math] coincide con los planos inferior y superior respectivamente; y aplicamos las condiciones de contorno f(z=0)=v y f(z=1)=0 lo que finalmente hace que lleguemos al campo vectorial de velocidades:

[math]f(z)=(1-z)(v+ \frac{(p_2-p_1)z}{2μ})[/math]

4 Campo de presiones y campo de velocidades

Suponiendo que [math]\ p_1=1 ; p_2=2 ; v=1 ; μ=1[/math], dibujaremos los respectivos campos de presiones y de velocidades.

CAMPO DE VELOCIDADES:

Para la expresión del campo de velocidades, hemos obtenido:

[math](1-\frac{z}{2}-\frac{z^2}{2})[/math]

De la fórmula anterior, podemos interpretar que la velocidad va a cambiar en función de z ,es decir, que variará en función de la altura. Los vectores del campo de velocidades son paralelos a las paredes del canal, dado que el campo cambia únicamente en la dirección j.

Mediante MATLAB, hemos logrado representarlo usando el siguiente código:

[Y,Z] = meshgrid([0:0.25:8],[0:0.1:1]); 
 uy=inline('((1-z./2-z.^2./2))','y','z');
 uz=inline('0.*y','y','z'); 
 U=uy(Y,Z);
 V=uz(Y,Z);
 quiver(Y,Z,U,V)
 axis([0,8,-1,2])

Campo de velocidad del fluido

Interpretación de la figura

Analizando la figura, podemos observar que la velocidad es nula en la pared del canal en z=1. También, observamos que la velocidad máxima se alcanza z=0 y como habíamos deducido anteriormente, la velocidad es paralela a las paredes del canal.

CAMPO DE PRESIONES:

Utilizando la fórmula del campo de presiones:

[math] p(x,y,z)=p_1+(p_2-p_1)(y-1)[/math]

hemos obtenido:

[math] p(x,y)=y[/math]

Logrando representarlo mediante el siguiente código:

y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z); 
figure (1);
p=Y;
surf(Y,Z,p);
view(2);
axis([0,8,-1,2]);

Campo de presiones del fluido

Interpretación de la figura

Como podemos observar en la figura a medida que el fluido avanza por el canal la presión va aumentando. Las presiones altas están representadas en amarillo y las más bajas en azul. El aumento de presión en una tubería o canal es el aumento de energía dinámica del fluido debido al uso de dispositivos específicos diseñados para este propósito (como por ejemplo: bombas o compresores).

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

Ahora vamos a dibujar las líneas de corriente del campo [math]\vec{u}[/math] (las líneas que son tangentes a [math]\vec{u}[/math] en cada punto). Para esto calculamos el campo [math]\vec{v}[/math] que en cada punto es ortogonal a [math]\vec{u}[/math], utilizando:

[math]\vec{v}=\vec{i}\times\vec{u}[/math]

Que al sustituir en nuestro campo de velocidades, nos queda:

[math]\vec{i}\times\vec u=\vec{i}\times\ f(z)\vec{j}=\ f(z)\vec{k}[/math]

Ahora vamos a demostrar que [math]\vec{v}[/math] es irrotacional:

[math]\nabla\times\vec v= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ 0 & 0 & \ f(z)\end{vmatrix}[/math]

Operando obtenemos, que: [math]\nabla\times\vec{v}=0[/math], por lo que podemos confirmar que el campo es irrotacional.

Ya que sabemos que es irrotacional, vamos a calcular [math]\psi[/math], el cual es un potencial escalar que se conoce como función de corriente de [math]\vec{u}[/math]

• Aplicamos la definición:

[math] \nabla \psi (y,z)= \vec{v} \longleftrightarrow \frac{\partial \psi }{\partial y}=\ v_1 \ {,} \ \frac{\partial \psi}{\partial z}=\ v_2 [/math]

• Operamos:

[math]\psi (y,z)=\int 0 \ dy=0+f(z)[/math]

• Derivamos con respecto a [math]z[/math], y teniendo en cuenta [math]p1=1,p2=2,μ=1,v=1[/math]

[math] \frac{\partial \psi}{\partial z}=(1-z)(1+\frac{z}{2})=(1-\frac{z}{2}-\frac{z^2}{2}) \rightarrow\frac{\partial }{\partial z}(0+f(z))\rightarrow f'(z)=1-\frac{z}{2}-\frac{z^2}{2}[/math]

• Integrando:

[math]f(z)=\int (1-\frac{z}{2}-\frac{z^2}{2})dz = z -\frac{z^2}{4}-\frac{z^3}{6}+C \ {,} \ C\in\mathbb{R}[/math]

• Simplificando:

[math] \psi (y,z)=z -\frac{z^2}{4}-\frac{z^3}{6}[/math]

Representando las líneas de corriente a través del siguiente código de MATLAB:

y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z); % Mallado
figure (1);
lineas=Z-(Z.^2)/4-(Z.^3)/6;
contour (Y,Z,lineas);
axis([0,8,-1,2]);
view (2);

Líneas de corriente de [math]\vec{u}[/math]

Interpretación de la figura:

Podemos ver que las líneas de corriente de [math]\vec{u}[/math] son tangentes a la velocidad del fluido. En régimen estacionario las líneas de corriente están fijas y nos indican la trayectoria seguida por los elementos del fluido. El hecho de que las lineas de corriente no se corten, nos indica que nos encontramos ante un flujo que no tiene puntos singulares, es decir, no existen fuentes o sumideros.

6 Velocidad máxima del fluido

Ahora vamos a obtener los puntos de máxima velocidad, mediante el cálculo de las derivadas parciales del campo:

[math]\vec{u}(y,z) = (1-\frac{z}{2}-\frac{z^2}{2})\vec{j}[/math] [math]\frac{\partial \vec{u}} {\partial y}=0 \ {,} \ \frac{\partial \vec{u} }{\partial z}=-\frac{1}{2}-z=0 [/math]

Procedemos a igualar dichas derivadas parciales a 0, obteniendo así que los puntos de la recta donde la velocidad es máxima es [math]z=-\frac{1}{2}[/math]

Con los resultados obtenidos nos hemos percatado que la velocidad máxima que el fluido es capaz de alcanzar es en z=-0.5, debido a este resultado es normal que aparezcan nuevas interrogantes al respecto: ¿Es posible que la velocidad máxima esté fuera del recinto? Esto cobra sentido ya que si el recinto fuese más amplio, sería z=-0.5 la velocidad máxima, pero como el recinto está acotado en z=0, damos por sentado que la velocidad máxima dentro del espacio confinado entre las dos placas se encuentra en z=0.

7 Rotacional de [math]\vec{u}[/math]

En el contexto de un fluido viscoso entre dos placas contenidas en el plano YZ, el rotacional del campo de velocidad \(\nabla \times \mathbf{\vec{u}}\) resalta la rotación local del fluido. Para un movimiento bidimensional dentro del plano \(\mathbf{\vec{u}} = \langle 0, u_y, u_z \rangle\) y el rotacional se simplifica a \(\langle \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z}, 0, 0 \rangle\). La componente X describe la rotación alrededor del eje de propagación, indicando la generación de esfuerzos cortantes y la deformación del fluido entre las placas. Este concepto es crucial en la mecánica de fluidos para analizar fluidos viscosos en el plano YZ.

De acuerdo con lo anterior nos preguntamos ¿Cuáles son los puntos que poseen mayor rotacional?

Recordando que:

[math]f(z)\vec{j} = (1-z)(v+ \frac{(p_2-p_1)z}{2μ}) \vec{j} [/math]

Pasamos al calculo del rotacional utilizando la siguiente expresión:

[math]\nabla\times\vec u= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ 0 & f(z) & 0 \end{vmatrix} = - \frac{\partial }{\partial z}(f(z))\vec{i} [/math]

Operando, obtenemos:

[math]\nabla\times\vec{u}= v-\frac{(p_2-p_1)(1/2-z)}{μ}\vec{i}[/math]

Y sustituyendo los valores [math] p_1=1, p_2=2, μ=1, v=1 [/math]. Tenemos que [math]\nabla\times\vec u= z +\frac{1}{2}\ \vec{i}[/math], siendo su módulo [math]\left | \nabla\times\vec u \right |=z+\frac{1}{2}[/math]. Como podemos obervar solo depende del parámetro [math]z[/math].

Con el siguiente código de MATLAB,

y=0:0.05:8;
 z=0:0.05:1;
 [Y,Z]=meshgrid(y,z); 
 rota=abs(1/2+Z); 
 surf(Y,Z,rota)
 axis([0,8,-1,2])
 view(2)

Rotacional del campo velocidad

Como podemos ver en la gráfica el máximo del rotacional se alcanza en el punto [math]z=0[/math], alcanzando el mínimo en el punto [math]z=1[/math]. Se entiende por rotacional al operador vectorial que muestra la tendencia de un campo a inducir rotación alrededor de un punto, en nuestro caso alrededor del mínimo [math]z=0[/math], que a su vez en ese punto, ya se ha demostrado que la velocidad es máxima. El punto con menor tendencia a la rotación es [math]z=0[/math].Para el valor z=−1/2 (valor fuera del recinto), el rotacional es nulo y este punto coincide con la velocidad máxima. En el gráfico los puntos con mayor tendencia a la rotación están representados con colores cálidos (amarillo) y los puntos con menor tendencia a la rotación aparecen representados con colores fríos (azules).

8 Temperatura máxima

La temperatura a la que está sometido el fluido viene determinada por el campo escalar:

[math] T(ρ{,}θ)=log(1+ρ)*cos^2θ[/math]

La expresión del campo está en coordenadas cilíndricas. Trabajaremos a partir de ahora en coordenadas cilíndricas:

CAMPO DE TEMPERATURA

El campo de temperatura describe la distribución espacial de las temperaturas en un sistema físico y éste depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión:[math]T(ρ{,}θ)=log(1+ρ)*cos^2θ[/math], por lo tanto no depende del tiempo (es estacionario) y depende exclusivamente de las componentes espaciales [math](ρ{,}z)[/math]

Logrando representar el campo, mediante el siguiente código de MATLAB:

x=0:0.25:8;         %Tomamos la x como rho
y=0:0.05:2*pi;         %Tomamos la y como theta
[Mx,My]=meshgrid(x,y);         %Matriz del campo
Mz= log(1+Mx).*(cos(My)).^2;                                     
surf(Mx,My,Mz);                                    
shading flat

Campo de temperaturas del fluido

Interpretación de la figura:

Observamos que la temperatura es mayor cuando ρ aumenta (la temperatura mayor se representa con colores cálidos y la menor temperatura con colores fríos).

CURVAS DE NIVEL

Mediante un cambio en el programa, representamos las curvas de nivel:

x=0:0.25:8;                                         
y=0:0.05:2*pi;                                      
[Mx,My]=meshgrid(x,y);                              
%Matriz del campo
Mz= log(1+Mx).*(cos(My)).^2;                                                                        
pcolor(Mx,My,Mz);                                  
hold on                                           
contour(Mx,My,Mz,8,'w')                           
hold off

Líneas de nivel del campo T

Interpretación de la figura:

En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica.

9 Gradiente de la temperatura

Ahora procedemos a demostrar que el gradiente es ortogonal a las curvas de nivel de la temperatura:

El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente [math]\nabla T[/math].

GRÁFICA DEL GRADIENTE:

Mediante el siguiente código, hemos logrado representar el gráfico del gradiente de la temperatura:

x=0:0.1:8;                                         
y=0:0.1:2*pi;                                      
%Seguimos usando los parametros anteriores
[Mx,My]=meshgrid(x,y);                             
Mz= log(1+Mx).*(cos(My)).^2;
%Calculamos el gradiente
[Dx,Dy]=gradient(Mz,x,y);                        
figure;
hold on
axis equal
quiver(Mx,My,Dx,Dy);
hold off

Campo vectorial gradiente

En el gráfico se puede observar que hemos representado el gradiente como un campo vectorial

REPRESENTACIÓN DEL GRADIENTE Y LAS CURVAS DE NIVEL

Código MATLAB utilizado:

x=0:0.1:8;                                         
y=0:0.1:2*pi;                                      
[Mx,My]=meshgrid(x,y);                             
Mz= log(1+Mx).*(cos(My)).^2;      
[Dx,Dy] = gradient(Mz,x,y);                        
figure;
%Combinamos las líneas de nivel con el gradiente
contour(Mx,My,Mz)                                  
hold on 
quiver(Mx,My,Dx,Dy);                               
hold off

Gradiente ortogonal a las curvas de nivel

Interpretación de la figura:

El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente que indica la dirección de crecimiento de la función temperatura en cada punto. En el gráfico aparece el gradiente como un campo vectorial, como ya se ha mencionado, junto con las curvas de nivel. Se observa que en torno a los puntos [math]y=1[/math] e [math]y=-1[/math] el módulo del gradiente aumenta por lo que la temperatura varía más rapidamente, en torno a esos puntos. Por otra parte, se observa que el gradiente es ortogonal a las curvas de nivel, la razón de esta ortogonalidad se duduce de la definición de gradiente y curvas de nivel, el gradiente indica la dirección de crecimiento y sentido de los valores del campo en cada punto, y las curvas de nivel son el lugar geométrico de los puntos equipotenciales. Si el gradiente tuviera componentes paralelas a las curvas las definiciones anteriores serían contradictorias.

10 Caudal del canal

Suponiendo que la velocidad está dada en m/s y que el canal tiene una profundidad de 1 metro, ¿qué caudal lleva el canal?

El caudal es el volumen de fluido que pasa a través del canal por unidad de tiempo y se calcula mediante la integral de superficie siguiente:

[math]\int_{S}\vec{u} \cdot d\vec{S}[/math]

donde [math]\vec{v}[/math] es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.

En nuestro caso, el campo de velocidades es:

[math](1-\frac{z}{2}-\frac{z^2}{2})\vec{j}[/math]

Esa expresión resulta de sustituir los valores [math]p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1[/math] en la expresión incial del campo:

[math]f(z)=(1-z)(v+ \frac{(p_2-p_1)z}{2μ})[/math]

CÁLCULO DEL CAUDAL

[math]\int_{S}\vec{u} \cdot d\vec{S}=\int_{0}^{1} \int_{0}^{1} (1-\frac{z^2}{2}-\frac{z^2}{2})dzdy=0.583\frac{m^3}{s}[/math]