Visualización de campos escalares y vectoriales en fluidos (Grupo A-17)
Contenido
1 Introducción
Se denomina fluido a un tipo de medio continuo formado por alguna sustancia entre cuyas moléculas hay una fuerza de atracción débil. La característica fundamental que define a los fluidos es su incapacidad para resistir esfuerzos cortantes (lo que provoca que carezcan de forma definida). El término, por tanto, engloba a los líquidos y a los gases.
1.1 Hipótesis básica
El estudio de la la mecánica de fluidos y las interacciones con el medio que lo limita se basa en la hipótesis del medio continuo: en esta hipótesis se considera que el fluido es continuo a lo largo del espacio que ocupa, ignorando por tanto su estructura molecular y las discontinuidades asociadas a esta. Con esta hipótesis se puede considerar que las propiedades del fluido (densidad, temperatura, etc.) son funciones continuas.
1.2 Aproximaciones al movimiento de un fluido
A la hora de describir el movimiento de un fluido existen dos puntos de vista:
Una primera forma de hacerlo es seguir a cada partícula fluida en su movimiento, de manera que buscaremos unas funciones que nos den la posición, así como las propiedades de la partícula fluida en cada instante. Esta es la descripción Lagrangiana del movimiento.
Una segunda forma es asignar a cada punto del espacio y en cada instante, un valor para las propiedades o magnitudes sin importar qué partícula fluida ocupa ese volumen diferencial en ese instante. Esta es la descripción Euleriana del movimiento, que resulta más útil por lo que será la usada en este estudio.
Las tres ecuaciones fundamentales que rigen los fluidos son:
- la ecuación de continuidad
- la ecuación de la cantidad de movimiento, y
- la ecuación de conservación de la energía.
Estas ecuaciones pueden darse en forma integral o en forma diferencial. Para llegar a su formulación diferencial, Navier y Stokes las unificaron aplicando ciertas consideraciones (principalmente la Ley de viscosidad de Newton, que establece una relación lineal entre los esfuerzos tangenciales y el gradiente de la velocidad).
De esta manera se obtiene la formulación diferencial, conocida como Ecuaciones de Navier-Stokes, muy útiles para la resolución de problemas de mecánica de fluidos:
[math] \rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} \right) = -\nabla p + \nabla \cdot\boldsymbol{\mathsf{T}} + \mathbf{f} [/math]
En el caso de fluidos ([math]\mathsf{T}=\mu \nabla \vec v[/math]) incompresibles ([math]\rho[/math] constante), circulando en régimen estacionario ([math]\frac{\partial \mathbf{v}}{\partial t}=0[/math]), y despreciando los efectos de la gravedad ([math]\rho=1 , \mathbf{f}=0[/math]), la ecuación queda simplificada de la siguiente forma:
[math] \vec u \cdot \nabla \vec u + \nabla p = \mu \Delta \vec u, [/math] donde [math]\vec u[/math] representa el campo de velocidades del fluido, [math]p[/math] el campo de presiones y [math]\mu[/math] la viscosidad del fluido.
Dichas ecuaciones son un conjunto de ecuaciones en derivadas parciales no lineales. No se dispone de una solución general para este conjunto de ecuaciones, y salvo ciertos tipos de flujo y situaciones muy concretas, no es posible hallar una solución analítica; por ello en muchas ocasiones es preciso recurrir al análisis numérico computacional para determinar una solución aproximada.
2 Estudio de un caso particular
En este artículo estudiaremos el caso particular de un fluido incompresible a través de un canal con paredes rectas.
La velocidad del fluido viene dada por la expresión: [math] \vec u (x,y) = y (1-y)(p_{1}-p_{2})/(2\mu) \vec i[/math] y la presión por: [math] p(x,y) = p_{1}+(p_{2}-p_{1})(x-1) [/math] siendo [math] p_{1}[/math]la presión en los puntos [math]x=1[/math], [math] p_{2}[/math]la presión de los puntos en [math]x=2[/math] y [math] \mu[/math] el coeficiente de viscosidad del fluido.
Es fácil comprobar que la velocidad en las paredes [math]\vec u (y=0)[/math] y [math]\vec u (y=1)[/math] es nula, simplemente sustituyendo los valores de [math]y[/math] en la ecuación del campo de velocidades.
2.1 Estacionalidad
Para saber si este campo de velocidades corresponde a un régimen estacionario, comprobaremos que verifica la ecuación de Navier-Stokes estacionaria: [math] \vec u \cdot \nabla \vec u + \nabla p = \mu Δ \vec u[/math] Para ello desarrollaremos por separado ambos miembros de la igualdad, verificando que llegamos al mismo resultado. Así: [math] \nabla \vec u = \frac{ \partial u_{i}}{ \partial x^{j} } = \begin{bmatrix} \frac{\partial u_1}{\partial x} & \frac{\partial u_1}{\partial y} \\\frac{\partial u_2}{\partial x} & \frac{\partial u_2}{\partial y} \end{bmatrix} = \frac{(1-2y)(p_{1}-p_{2})}{(2\mu)}\vec j [/math]
resultando [math] \vec u \cdot \nabla \vec u = \frac{y(1-y)(p_{1}-p_{2})}{(2\mu)}\vec i \cdot \frac{(1-2y)(p_{1}-p_{2})}{(2\mu)}\vec j = 0 [/math] por ser el producto escalar de dos vectores ortogonales. El primer miembro queda reducido entonces a: [math] \nabla p = \frac{ \partial p }{ \partial x } \vec i + \frac{ \partial p }{ \partial y} \vec j = (p_{2}-p_{1}) \vec i [/math] Por otro lado, en el segundo miembro calculamos [math]Δ \vec u[/math]: [math] Δ \vec u = \nabla \nabla \vec u = \frac{ \partial^2 \vec u }{ \partial x^2 } + \frac{ \partial^2 \vec u }{ \partial y^2 }= \frac{ \partial^2 u_1 }{ \partial x^2 } \vec i + \frac{ \partial^2 u_2 }{ \partial x^2 } \vec j + \frac{ \partial^2 u_1 }{ \partial y^2 } \vec i + \frac{ \partial^2 u_2 }{ \partial y^2 } \vec j = \frac { p_2-p_1}{\mu} \vec i [/math] Multiplicando por la viscosidad, nos queda [math] \mu Δ \vec u =(p_2-p_1) \vec i [/math] , verificándose la ecuación estacionaria de Navier-Stokes.
2.2 Condición de incompresibilidad
Todos los fluidos, al igual que los sólidos, son -desde el punto de vista estrictamente físico- compresibles en cierto grado. No obstante, los líquidos son prácticamente incompresibles, a diferencia de los gases, que son altamente compresibles.
Para este análisis matemático, tomaremos el fluido como incompresible, es decir, su densidad se mantiene constante a lo largo del tiempo y del espacio. Para ello, la masa de fluido que entra al volumen de control debe ser igual a la que sale. Esto equivale a decir que la divergencia del campo de velocidades de un fluido incompresible es nula:
[math] \nabla \cdot \vec u = 0 [/math]
La demostración es trivial:
[math] \nabla \cdot \vec u = \frac{ \partial \vec u_1 }{ \partial x_1} + \frac{ \partial \vec u_2 }{ \partial x_2} = \frac{ \partial (y(1-y)(p_{1}-p_{2})/(2\mu))}{ \partial x} = 0 [/math]
3 Análisis del campo. Representación gráfica
A partir de ahora, supondremos los valores de [math]p_{1}=2[/math], [math]p_{2}=1[/math] y [math]\mu=1[/math]. Por lo tanto, los campos de velocidades y presiones nos quedan: [math] \vec u(x,y)= \frac{y(1-y)}{2}\vec{i} [/math]: [math] \vec p(x,y)= 3-x[/math]
3.1 Mallado
Para todo el análisis siguiente usaremos una malla cartesiana que nos represente los puntos interiores de un plano horizontal del canal. Asociada a esta malla, con su misma distancia entre puntos (0.1), se creará la matriz de los puntos sobre los que aplicaremos el modelo de cálculo numérico.
En la figura de la derecha podemos ver una representación de dicha malla, y debajo el código en OCTAVE para crearla:
x=0:0.1:4; % intervalo de X de 0 a 4
y=0:0.1:1; % intervalo de Y de 0 a 1
[xx,yy]=meshgrid(x,y); % creamos la malla
figure(1)
mesh(xx,yy,0*xx) % dibujamos la malla
axis([0,4,-1,2]) % elegimos la región a representar
view(2) % visualizar la imagen en el plano XY
3.2 Campo de velocidades y campo de presiones
Para el análisis de campos, la representación gráfica es una herramienta muy intuitiva que nos permite hacernos una idea rápida y visual de cómo es el campo. Existen dos tipos de representaciones:
- la representación de campos escalares: En estos campos, a cada punto del plano ($\mathbb{R^2}$) o del espacio ($\mathbb{R^3}$) le corresponde un valor escalar; al representarlo gráficamente, este valor se representa por variables escalares como la altura (resultando mapas 3D) o el color (asignando mediante una función un rango de colores para un determinado rango de valores).
- la representación de campos vectoriales: En estos campos, a cada punto del plano o del espacio le corresponde un valor vectorial. Por tanto, la representación gráfica debe ser orientada, por lo que se usan flechas de distinto tamaño y orientación para representar la magnitud vectorial.
Vamos a representar gráficamente a continuación los campos de velocidades (vectorial) y presiones (escalar)
La representación del campo de velocidades nos indica gráficamente que el fluido va ganando velocidad a medida que se acerca al plano central del canal, por lo que en dicho plano la velocidad será máxima. Esto se debe a las propiedades de viscosidad explicadas en el primer apartado del artículo: el fluido, al desplazarse en un espacio confinado genera un rozamiento con las paredes hasta el punto que llegará a generarse una finísima capa de fluido que tendrá la misma velocidad que el recipiente, en este caso nula.
En el caso de la presión, los colores más cálidos representan valores mayores de presión, y los más fríos valores menores.
Observamos cómo la presión disminuye en el sentido de avance del fluido. Al contrario de lo que ocurría con la velocidad, cuyo valor variaba con la Y pero se mantenía constante según X, la presión varía según nos movemos a lo largo de X, pero permanece constante al movernos a lo largo de Y
A continuación se adjuntan los códigos de OCTAVE correspondientes a las representaciones de los campos de velocidades y de presiones:
x=linspace(0,4,8); %División del dominio de definición de X (perfil longitudinal)
y=linspace(0,1,10); %División del dominio de definición de Y (perfil transversal)
[X,Y]=meshgrid(x,y); %Matriz asociada a la división del dominio
U=(Y.*(1-Y))/2; %Componente horizontal (i) de los vectores del campo
V=0*U; %Componente vertical (j) de los vectores del campo
quiver(X,Y,U,V) %Representación del campo de velocidades del fluido
axis([0,4,-1,2]) %Extensión del intervalo mostrado
print('velocidades.png') %Imagen del gráfico mostrado en formato .png
x=linspace(0,4,400); %División del dominio de definición de X (perfil longitudinal)
y=linspace(0,1,100); %División del dominio de definición de Y (perfil transversal)
[X,Y]=meshgrid(x,y); %Matriz asociada a la división del dominio
Z=3-X; %Valor del campo de presiones
h=pcolor(X,Y,Z); %Representación del campo de presiones mediante una gama de colores
axis([0,4,-1,2]) %Extensión del intervalo mostrado
print('presiones.png') %Imagen del gráfico mostrado en formato .png
3.3 Lineas de corriente del campo
Para hallar las lineas de corriente del campo [math]\vec{u}[/math] tenemos que tener en cuenta que son tangentes en todo momento al propio campo [math]\vec{u}[/math]. Nos vamos a ayudar para su cálculo de un campo [math]\vec{v}[/math] que será ortogonal al campo y por lo tanto también a las líneas de corriente. En nuestro caso llevará dirección [math]\vec{j}[/math]
Sabemos que nuestro campo [math]\vec{u}[/math] es adivergente (ver apartado 2.1), de lo que deducimos que el campo [math]\vec{v}[/math], que es ortogonal a él, es irrotacional, por lo que [math]\vec{u}[/math] tiene un potencial escalar [math]\psi[/math]. Para calcularlo, usaremos la metodología explicada en clase.
[math]\frac{ \partial \psi}{\partial y}=\frac{y(1-y)}{2} [/math] ; [math]\psi =\int\frac{y(1-y)}{2}dy =\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})+ f(x) [/math]
[math]\frac{ \partial \psi}{\partial x}=0 [/math] ; [math]\frac{\partial}{\partial x}(\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})) + f'(x)=0 \implies f'(x)=0 \implies f(x)=cte[/math]:
[math] \psi =\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3}) + C [/math]
A continuación mostramos el código de OCTAVE para representar las lineas de corriente, que serán las del potencial escalar para [math]\psi = cte[/math]
x=0:0.1:4; % creamos el vector x con valores entre 0 y 4 equidistantes 0.1
y=0:0.1:1; % creamos el vector y con valores entre 0 y 1 equidistantes 0.1
[xx,yy]=meshgrid(x,y); % malla para las coordenadas
potencial= (1/2).* ((yy.^2/2)-(yy.^3/3) ); % función del potencial escalar
contour(xx,yy,potencial) %Dibujamos las líneas de nivel en el plano XY
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XY
3.4 Rotacional
El rotacional de [math]\vec{u}[/math] se calcula como:
[math]\nabla\times\vec{u}= \begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ u_1 & u_2 & u_3 \end{pmatrix}=\frac{\partial}{\partial y}(\frac{y(1-y)}{2})\vec{k} =\frac{1-2y}{2}\vec{k} [/math]
Su módulo es: [math]|\nabla\times\vec{u}|=\frac{1-2y}{2}[/math]
Observamos que en [math]y=0[/math] el módulo del rotacional vale [math]1/2[/math], mientras que en [math]y=1[/math] vale [math]-1/2[/math]. Esto se refleja en el gráfico de la derecha, donde el color rojo representa el valor máximo y el azul el valor mínimo. En el centro, el rotacional vale 0.
Aunque todas las partículas del fluido tengan su velocidad orientada en la misma dirección -lo cual podría hacernos pensar que el rotacional vale 0 en todos los puntos-, la diferencia (gradiente) de magnitudes de esas velocidades hace que el rotacional no se anule, salvo en los puntos del centro del canal, donde por simetría la velocidad vale lo mismo a un lado y a otro. En cambio, junto a las paredes del canal, el gradiente de velocidades es máximo, por lo que el valor del rotacional también lo es.
A continuación se muestra el código de OCTAVE mediante el cual se ha representado gráficamente el rotacional de [math]\vec{u}[/math].
x=0:0.1:4; % definimos el intervalo en el eje X
y=0:0.1:1; % definimos el intervalo en el eje Y
[xx,yy]=meshgrid(x,y); % creamos la malla
figure(1)
f=1/2.*(1-2.*yy); % definimos el valor del rotacional
surf(xx,yy,f) % dibujamos el campo
axis([0,4,-1,2]) % elegimos región en la que dibujar
view(2) % visualizamos la representación
4 Temperatura del fluido
Por último, vamos a estudiar el campo de temperaturas asociado al fluido, que suponemos dado por:
[math]
T(x,y)=(x-1)^2-y^2
[/math]
Al igual que hemos hecho con el campo de presiones, la temperatura, al ser otra variable escalar, puede representarse también con una representación de campo escalar, en este caso asignando colores a valores de temperatura, como se muestra en la gráfica de la derecha.
Debajo se muestra el código de OCTAVE usado para obtener dicha representación.
x=0:0.05:4; %Vector en X de la malla
y=0:0.05:1; %Vector en Y de la malla
[xx,yy]=meshgrid(x,y); %Malla
figure(1)
ff=(xx-1).^2-yy.^2; %Campo de temperaturas
surf(xx,yy,ff) %Dibujo del campo
axis([0,4,-1,2]) %Ejes
view(2) %Visualizar el campo
print('Temperaturas.jpeg')%Imprime la imagen en formato .jpeg
4.1 Gradiente de temperaturas
El gradiente de temperaturas es un campo vectorial asociado al campo escalar de temperaturas. Este campo vectorial muestra la dirección, sentido y magnitud de la máxima variación del campo escalar para cada punto; y se obtiene derivando la expresión de este campo escalar respecto de cada una de las variables: :[math]\nabla \ T = \frac{ \partial T}{ \partial x}\vec i + \frac{ \partial T}{ \partial y}\vec j= (2x-2) \vec i - (2y) \vec j[/math]
A continuación se muestra el código de OCTAVE para obtener la representación del gradiente de la temperatura:
x=0:0.2:4; %Vector en X de la malla
y=0:0.2:1; %Vector en Y de la malla
[xx,yy]=meshgrid(x,y); %Malla
figure(1)
fx=2.*xx-2; %Gradiente de temperaturas en i
fy=-2.*yy; %Gradiente de temperaturas en j
quiver(xx,yy,fx,fy) %Dibujo del gradiente
axis([0,4,-1,2]) %Ejes
view(2)
print('gradienTemp.jpg') %Imprime la imagen en formato .jpg
Por último, vamos a comprobar gráficamente que el gradiente de temperaturas es ortogonal a las curvas de nivel. Para ello, en la imagen de la derecha, hemos representado simultáneamente gradiente y curvas.
Debajo tenemos el código de OCTAVE que nos ha permitido obtener dicha representación.
x=0:0.1:4; %definimos los valores de X
y=0:0.1:1; %definimos los valores de Y
[xx,yy]=meshgrid(x,y); %matriz de las coordenadas X e Y
t=((xx-1).^2) -(yy.^2) %definimos en la variable t el campo de temperaturas
contour (xx,yy,t) %dibujamos las lineas de nivel
hold on %mantenemos el dibujo de las curvas para añadir el del gradiente
fi=2.*xx-2; %componente i del campo vectorial gradiente
fj=-2.*yy; %componente j del campo vectorial gradiente
quiver (xx,yy,fi,fj) %dibujamos el campo gradiente de temperaturas
axis([0,4,-1,2]) %fijamos el rango de representación
view(2) %mostramos ambas gráficas a la vez
