Diferencia entre revisiones de «Flujo de Couette grupo 21»

Revisión actual del 11:54 15 dic 2023

1 INTRODUCCIÓN

Este proyecto tiene como objetivo analizar y representar visualmente tanto campos escalares como vectoriales, que fueron abordados durante el curso de grado, en el contexto de fluidos. En este estudio, nos enfocaremos en un fluido incompresible con fluctuaciones de caudal a lo largo de un canal con paredes horizontales rectas, sometido simultáneamente a tres campos: dos escalares (presión y temperatura) y uno vectorial (velocidad).

La cinemática de los fluidos se centra en el movimiento de estos sin considerar sus causas subyacentes, destacando aspectos como trayectorias, velocidades y aceleraciones. Además, reconocemos que un fluido incompresible mantiene una densidad constante a lo largo del tiempo y posee la capacidad de resistir la compresión en diversas condiciones.

Para llevar a cabo este trabajo, empleamos el programa informático Matlab, el cual nos facilitó la visualización del comportamiento del fluido a través de gráficos.

2 SUPERFICIE DE TRABAJO

La superficie en la que nos vamos a basar va a ser [0,8]x[0,1] en [math]\{\vec{j} {,} \vec{k} \}[/math], por lo que trabajamos en el plano x=0, entonces definimos la y en [0,8] y la z en [0,1], aunque el eje z lo definimos como [-1,2]. Para obsérvalo vamos a recurrir a octave:

y=0:0.1:8; %vector y
z=0:0.1:1; %vector z
[yy,zz]=meshgrid(y,z); %mallado ZY
hold on
mesh(yy,zz,0.*yy); %representar el canal
axis([0,8,-1,2]);
axis equal
view(2);
title ('Mallado del canal');
hold off

La superficie de trabajo es la que se muestra en la siguiente imagen:

Mallado

3 ECUACIÓN DE NAVIER-STOKES ESTACIONARIA

En el ámbito de la dinámica de fluidos, las ecuaciones de Navier-Stokes son expresiones matemáticas que describen el movimiento tridimensional de sustancias fluidas viscosas. Estas ecuaciones encuentran aplicaciones diversas, siendo utilizadas para prever fenómenos como el clima, las corrientes oceánicas, el flujo de agua en sistemas de tuberías o reactores, así como en el análisis del flujo sanguíneo, entre otros usos, como el diseño de submarinos. En esta sección específica, nuestro propósito es demostrar que el campo de velocidad y el campo de presión al que está sujeto el fluido cumplen con la ecuación estacionaria de Navier-Stokes. Esto implicaría que el fluido es incompresible, ya que estas ecuaciones modelan el comportamiento de los fluidos newtonianos, es decir, aquellos en los que la resistencia a deformaciones puede considerarse constante a lo largo del tiempo.

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

Trabajando en componentes tenemos que:

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

En nuestro caso, [math]\vec{u}(y,z)=f(z)\vec{j} [/math] es el campo vectorial velocidad , [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.

Para la demostración, se trabajará en componentes respecto de la base cartesiana [math]\{\vec{j} {,} \vec{k} \}[/math].

El primer término de la ecuación será:

Sustituyendo, nuestro campo y operando tenemos:

Cálculamos ahora, el gradiente del campo de presiones, [math]\nabla p[/math], siendo [math]\ p(x,y,z)=p_1+(p_2-p_1)(y-1)[/math]

Por último, calculamos el Laplaciano,[math]\Delta\vec{u} [/math] , del campo de velocidades, [math]\vec{u}(y,z)=f(z)\vec{j} [/math]. Para calcularlo, utilizaremos la siguiente fórmula:

1. Cálculo de la divergencia:

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

Obtenemos que el fluido es incompresible, pues la condición de incompresibilidad es [math]\nabla \cdot \vec{u}=0[/math]. Teniendo en cuenta que la divergencia mide el cambio de volumen del fluido inducido por el campo, ésta al ser nula, conlleva que el movimiento de las partículas no afecta al volumen provocando que la densidad permanezca constante en el tiempo, coincidiendo con la definición dada de fluido incompresible ya mencionada anteriormente.

2. Cálculo del rotacional:

3. Cálculo de [math]\nabla\times(\nabla\times\vec u)[/math] (doble rotacional)

Sustituimos los resultados anteriores en [math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math], obtenemos que:

Sustituimos por último, todos los términos en la expresión de la ecuación de Navier-Stokes estacionaria [math](\vec{u} \cdot \nabla) \vec{u} + \nabla p=μ \Delta\vec{u} [/math]; resultando:

Por lo que:

Para averiguar f(z) tenemos que integrar 2 veces sobre [math]\frac{p_1-p_2}{-μ}[/math]. Al hacerlo y aplicando las condiciones de contorno proporcionadas por el enunciado donde f(z=0)=v y f(z=1)=0:

4 CAMPO DE VELOCIDADES

Para la representación del campo se han tomado los siguientes valores:

Sustituyendo en la expresión del campo obtenemos:

Para su representación se ha implementado en Matlab, el siguiente código:

% Definir el dominio de y y z
[y, z] = meshgrid(0:0.25:8, 0:0.1:1);

% Definir las funciones vectoriales uy(y, z) y uz(y, z)
uy = (4-3*z-z.^2);
uz = zeros(size(y));

% Visualizar los resultados con quiver
quiver(y, z, uy, uz);
xlabel('y');
ylabel('z');
title('campo velocidades');
axis([0, 8, -1, 2]);
grid on;

Resultando el siguiente gráfico:

Lo primero que observamos, es que la velocidad es nula la pared superior del canal pues en [math]z=1[/math] no existen líneas de campo, como ya habíamos demostrado analíticamente, mientras que en y [math]z=0[/math] si. Además, en las deducciones previas habíamos asegurado que la velocidad sería paralela a las paredes del canal, cosa que también observamos en la gráfica, pues las líneas de campo son paralelas a las rectas [math]z=1[/math] y [math]z=0[/math] que coinciden con las paredes del canal a la vez que se forma la mencionada ya curva parabólica que sigue nuestra función solución de la ecuación diferencial al sustituir los valores mencionados.

Por otro lado, se intuye que la velocidad máxima se alcanza en [math]z=-1/2[/math] como a continuación se demuestra analíticamente.

4.1 CÁLCULO DE LA VELOCIDAD MÁXIMA DEL FLUIDO

Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:

Igualando a 0, las derivadas anteriores, obtenemos que la velocidad máxima se obtienen en los puntos de la recta [math]z=-1/2[/math], interpretando así que la velocidad será máxima fuera del recinto propuesto coincidiendo así con el máximo de dicha función parabólica.

5 CAMPO DE PRESIONES

El campo de presiones al que está sometido nuestro fluido es un campo escalar, con las siguiente expresión:

Para representarlo gráficamente se consideran nuevamente los valores: [math]p_1=1 \ {,} \ p_2=2 \ {,} \ μ=1\ {,} \ v=1 [/math] .Obteniendo, por tanto:

Implementado el código siguiente en Matlab:

y=0:0.1:8; %vector y
z=0:0.1:1; %vector z
[Y,Z]=meshgrid(y,z); % Mallado
p=Y;
surf(Y,Z,p)
title('Campo de presiones');
view(2)
axis([0,8,-1,2]);

Obtenemos el siguiente gráfico:

Campo de presiones

Observamos en el gráfico que a medida que el fluido avanza por el canal la presión va disminuyendo, las presiones altas están representadas con colores azules y derivados de éste, y las más bajas con tonos anaranjados. Esto es así, porque para mantener el caudal de un fluido viscoso estable debe mantenerse una diferencia de presiones entre las paredes del canal. Esta diferencia de presión es necesaria debida a la fuerza de arrastre o frenada que ejerce el canal sobre la capa de fluido en contacto con él y la que ejerce cada capa de fluido sobre la adyacente que se está moviendo con distinta velocidad. A estas fuerzas las denominamos fuerzas viscosas. El resultado de su presencia, hace que la velocidad del fluido no sea constante a lo largo del canal siendo mayor cerca del centro y menor cerca de las paredes tal como se explicó y demostró en el apartado anterior.

6 ROTACIONAL

Para calcular el rotacional del campo, utilizamos la siguiente expresión:

Operando obtenemos, que: [math]\nabla\times\vec{u}=-f'(z)\vec{i}[/math], sustituyendo los valores:

Tenemos que: [math]\nabla\times\vec u=-z+C_1\vec{i}[/math], donde la constante 1 ya ha sido calculada anteriormente aplicando las condiciones de contorno y siendo igual a -1/2

Finalmente obtenemos el módulo del rotacional, [math]\left | \nabla\times\vec u \right |=z+1/2[/math], como se observa en la expresión, depende del valor de [math]z[/math]. Tomando los valores máximos en los extremos del intervalo en el que está definido [math]z[/math], [math] [0,1][/math], y el valor mínimo en [math]z=-1/2[/math]

Para visualizar esto de manera gráfica, nuevamente recurrimos a Octave,

y=0:0.01:8;
z=0:0.01:1;
[Y,Z]=meshgrid(y,z); %Mallado
rota=abs(Z); %Rotacional
surf(Y,Z,rota)
axis([0,8,-1,2])
view(2)

Resultando el siguiente gráfico, donde se representa el rotacional:

Rotacional del campo velocidad

El rotacional muestra la tendencia del campo a inducir rotación alrededor de un punto. Por otra parte, el rotacional caracteriza la rotación de un fluido, por lo que en los extremos del canal en efecto de giro será máximo, particularizando en [math]z=0[/math] tendrá el valor máximo en el sentido positivo de [math]\vec{i}[/math] , lo que implica rotaciones en sentido antihorario, y en [math]z=1[/math] no alcanzará ningun valor debido a las hipótesis tomadas en las condiciones de contorno atribuyendo a este 0 mencionadas por el enunciado ,lo que provocará rotaciones en sentido horario. Para el valor [math]z=-1/2[/math], que es un valor fuera de nuestro recinto, el rotacional es nulo y este punto coincide con la velocidad máxima. Con esto quedan demostradas, las hipótesis anteriores, pues en el gráfico los puntos con mayor tendencia a la rotación aparecen representados con colores cálidos, amarillos hacia verde, y los puntos con menor tendencia a la rotación aparecen representados con colores fríos, azules.

7 CAMPO DE TEMPERATURAS

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

Se dibuja el campo de temperaturas y las curvas de nivel:

En Matlab hemos empleado el siguiente programa:

x=1:0.05:2;                                         %Definición de la x (rho)
y=0:0.05:2*pi;                                      %Definicón de la y (teta)
[Mx,My]=meshgrid(x,y);                              %Matriz de los parámetros
Mz= log(1+Mx).*(cos(My)).^2;                        %Matriz de la z
subplot(1,2,1);                                     %Ventanas 
surf(Mx,My,Mz);                                     %Dibujo de la superficie 
shading flat                                        %Difuminado
subplot(1,2,2)                                      %Ventana 2
pcolor(Mx,My,Mz);                                   %Proyección en planta
hold on                                                  %Mantener ventana
contour(Mx,My,Mz,7,'k')                             %7 curvas de nivel
hold off

A continuación, representamos una gráfica donde observamos en que punto la temperatura es máxima.

En Matlab hemos utilizado el siguiente código:

x=1:0.05:2;                                     %Definición de la x (rho)
y=0:0.05:2*pi;                                  %Definicón de la y (teta)
[Mx,My]=meshgrid(x,y);                          %Matriz de los parámetros
Mz= log(1+Mx).*(cos(My)).^2;             %Matriz de la z
plot3(Mx,My,Mz);                                %Dibujo en líneas del campo

8 GRADIENTE DE LAS TEMPERATURAS

Tomando los operadores "log" del campo de temperaturas como "ln" obtenemos un gradiente:

[math]\nabla T=(Cos^2(θ)*(\frac{1}{1+ρ}))\vec{e_ρ}+((\frac{1}{ρ})*ln(1+ρ)*2*Cos(θ)*(-1Sin(θ)))\vec{e_θ}[/math]

Finalmente obtenemos la siguiente gráfica:

Mediante la siguiente linea de código:

% Definir la función T y sus derivadas parciales
T = @(rho, theta) log(1 + rho).*cos(theta).^2;
dT_dRho = @(rho,theta) cos(theta).^2./(1+rho);
dT_dTheta = @(rho,theta) (2.*cos(theta).*log(1+rho)./rho).*(rho>0); % Evitar división por cero
% Definir el rango de valores para rho y theta
rho=1:0.1:2;
theta=0:0.1:(2*pi);
% Crear una malla de valores para rho y theta
[Rho, Theta]=meshgrid(rho,theta);
% Calcular las componentes del campo vectorial en la malla
dT_dRho_values=dT_dRho(Rho,Theta);
dT_dTheta_values=dT_dTheta(Rho,Theta);
% Graficar el campo de vectores usando quiver
figure;
quiver(Rho,Theta,dT_dRho_values,dT_dTheta_values);
xlabel('\rho');
ylabel('\theta');
title('Campo de Vectores \nabla T');
axis tight;
grid on;

Comprobamos así que el gradiente de temperaturas es ortogonal a las curvas de nivel del campo de temperaturas.

9 LÍNEAS DE CORRIENTE

Las líneas de corriente, son rectas tangentes al campo de velocidad en cada punto. Para calcular estas rectas, calculamos el campo [math]\vec{v}=\vec{i}\times\vec u[/math]. Dicho campo es:

Las líneas de corriente corresponden con [math]\psi=cte[/math], siendo [math]\psi[/math] la función de corriente de [math]\vec{u}[/math]. Esta función, es la función potencial o el potencial escalar del que deriva el campo [math]\vec{v}[/math], para calcularlo, lo primero que debemos comprobar es que [math]\vec{v}[/math] es irrotacional, pues en caso contrario, dicho campo no admitiría función potencial.

Demostración [math]\nabla\times\vec v=0[/math]

Operando obtenemos, que: [math]\nabla\times\vec{v}=0[/math], por lo que el campo es irrotacional y ahora ya podemos calcular su función potencial.

Cálculo de la función potencial [math]\psi (x,y)[/math]: Busco un campo escalar u tal que [math]\nabla\times\vec{u}=\vec{v}[/math]

1. Aplicamos la definición:

2. Operamos:

Teniendo en cuenta que las lineas de corriente dependerán de un parámetro, asociaremos dicho parámetro a la Cte obtenida mediante integración

Representación gráfica

Representamos las líneas [math]\psi=cte[/math]

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

Como se observa en el gráfico, estas líneas son tangentes al campo vectorial [math]\vec{u}[/math] en cada punto quedando así demostrado que son las líneas de corriente de [math]\vec{u}[/math]

Código Octave utilizado para la representación

y=0:0.01:8;
z=0:0.01:1;
[Y,Z]=meshgrid(y,z); %Mallado
lineas=(-((1-Z).^2)./2).*v+((p1-p2)./(2*mu)).*((Z.^2./2)-(Z.^3./3)); %ecuación 
contour(Y,Z,lineas)
axis([0,8,-1,2])
view(2)

10 CAUDAL

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:

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:

Esa expresión resulta de sustituir los valores siguientes en la expresión inicial del campo:

La profundidad del caudal es de 1 metro, por lo que los intervalos de integración serán 0 y 1 en ambos casos ya que es de sección cuadrada. Finalmente para la parametrización de esta tomamos como valores:

CÁLCULO DEL CAUDAL