Visualización de campos escalares y vectoriales en fluidos (Grupo A5)
El 24 de marzo de 2000 se anunció en el Collège de France de París el conjunto de los siete problemas que constituyen los Millennium Prize Problems, patrocinados por el Clay Mathematics Institute. En esta lista se recogieron siete de los problemas más importantes a los que se enfrentaría la ciencia en el nuevo siglo. Estos problemas recorren diversas áreas de las matemáticas, puras y aplicadas. La resolución de cada uno de ellos supondrá para el autor un premio de un millón de dólares. La cuestión número seis se enuncia como "Existencia y regularidad para las Ecuaciones de Navier-Stokes; Mecánica de Fluidos y EDPs[math]^{1}[/math]". Dicha ecuación es la que analizaremos a continuación.
[math]^{1}[/math] EDPs: ecuaciones en derivadas parciales.
Contenido
1 Introducción
Un fluidoes una sustancia que, debido a su poca cohesión intermolecular, carece de forma propia y adopta la forma del recipiente que lo contiene. Los fluidos están presentes en la naturaleza de manera sistemática. Su comportamiento influye en la manera en la que se producen distintos fenómenos; desde la manera en la que el agua se desplaza por el interior de una tubería a la formación de tornados o corrientes marinas.
El estudio y obtención de modelos matemáticos para fluidos surge del deseo humano del conocimiento y del control de los fenómenos que ocurren en la naturaleza.
Para establecer un modelo matemático adecuado es necesario conocer la densidad o viscosidad característica del fluido, además de los campos de presiones y velocidades a los que está sometido. Todo esto englobado en un determinado dominio ( de [math]\R²[/math] o [math]\R³[/math]).
Las ecuaciones fundamentales de los fluidos surgieron a partir de las leyes de la conservación de la masa y de la cantidad de movimiento de Newton. La identificación de las componentes de las fuerzas que actúan llevó un siglo y medio y en esta tarea participaron Johann, Bernoulli, Euler y Cauchy. Claude-Louis Navier y George Gabriel Stokes culminaron la modelización añadiendo la fórmula del efecto viscoso.
Los postulados de los fluidos perfectos tenían implícitas grandes simplificaciones. Esta situación fue solventada con el modelo de Navier-Stokes, que trataba fluidos más realistas al tener en cuenta el efecto de la viscosidad. Estas ecuaciones son utilizadas en numerosos campos, como la hidráulica, la meteorología o la aeronáutica, ya que describe el movimiento de las partículas del fluido.
Los matemáticos y los físicos creen que una explicación para la predicción tanto de la brisa como de las turbulencias que las corrientes de aire causan se puede encontrar a través de la comprensión de las soluciones a las ecuaciones de Navier-Stokes. Aunque estas ecuaciones fueron escritas en el siglo XIX, la comprensión de las mismas sigue siendo insuficiente. El reto es conseguir una teoría matemática que descubra los secretos escondidos en estas ecuaciones.
2.1 Campos de presiones y de velocidades
Vamos a considerar el flujo de un fluido incompresible a través de un canal con paredes rectas. Trabajaremos en el plano XY con coordenadas cartesianas.
Primero vamos a ver cómo dibujar un mallado que represente los puntos interiores del rectángulo [0,4]x[0,1] ocupado por un fluido, con ejes en [0,4]x[-1,2]. Sobre este rectángulo veremos los campos escalares y vectoriales que hemos estado viendo. Para ello utilizamos el siguiente código en MatLab:
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); % matrices de las coordenadas x e y
mesh(xx,yy,xx*0) % dibujamos el mallado
axis([0,4,-1,2]) % fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XY
La visualización de campos escalares y vectoriales es una herramienta muy utilizada en el estudio del comportamiento de los fluidos, ya que permiten apreciar de manera gráfica las acciones y efectos a los que está sometido un fluido. En esta imagen podemos ver la sección longitudinal del canal con paredes rectas que contiene al fluido incompresible sobre el que trata este trabajo.
A partir de ahora todo será referido al campo vectorial
[math] \vec u (x,y)= \frac{y(1-y)(p_{1}-p_{2})}{(2\mu)}\vec i[/math]
que representa la velocidad de las partículas de un fluido y al campo escalar [math] p(x,y)=p_{1}+(p_{2}-p_{1})(x-1) [/math] que representa su presión, 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.
Dibujamos ahora el campo de presiones y el de velocidades tomando [math] p_{1}=2[/math] [math], p_{2}=1[/math] y [math] \mu=1[/math], quedando los campos: [math] \vec u(x,y)=\frac{1}{2}(y-y^2) , p(x,y)=3-x [/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); % matrices de las coordenadas x e y
figure (1)
f=3-xx %campo de presiones
surf(xx,yy,f) %dibujamos la superficie
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XY
x=0:0.2:4; % creamos el vector x con valores entre 0 y 4 equidistantes 0.2 y 0.1 respectivamente
y=0:0.1:1; % creamos el vector y con valores entre 0 y 1 equidistantes 0.1
[xx,yy]=meshgrid(x,y); % matrices de las coordenadas x e y
figure (2)
fx=(1/2).*(yy-(yy.^2)); %campo de velocidades en x
fy=0.*xx; %campo de velocidades en y, en nuestro caso, nulo
quiver(xx,yy,fx,fy) %dibujamos el campo vectorial de velocidades
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XYLos colores cálidos representan los valores más altos y los fríos los más bajos, así vemos como la presión disminuye según el fluido avanza
Vemos que la velocidad es máxima en el centro y nula en las paredes
Vamos a comprobar que el campo vectorial verifica la ecuación de Navier-Stokes estacionaria: [math] \vec u \cdot \nabla \vec u + \nabla p = \mu Δ \vec u [/math] Calculamos ambos miembros de la ecuación: [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]:
[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]
Así pues, 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] En el segundo miembro calculamos el laplaciano de [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] El segundo miembro queda entonces [math] \mu Δ \vec u =(p_2-p_1) \vec i [/math] que vemos que es igual al primer miembro, luego se cumple la igualdad de Navier-Stokes.
2.3 Condición de incompresibilidad
La característica de incompresibilidad de un fluido implica que este tiene densidad constante en todo el conjunto a lo largo del tiempo. Por lo tanto, la masa de fluido que entra y la que sale a través de la sección de área de un canal con un volumen determinado será la misma. Esto equivale a decir que la divergencia del campo de velocidades de un fluido es nula y por consiguiente el fluido es incompresible: [math] \nabla \cdot \vec u = 0[/math]: [math]\nabla \cdot \vec u = \frac{ \partial \vec u_1 }{ \partial x_1} + \frac{ \partial \vec u_2 }{ \partial x_2} + \frac{ \partial \vec u_3 }{ \partial x_3} = \frac{ \partial}{\partial x_1}(\frac{y(1-y)(p_{1}-p_{2})}{2\mu}) = 0 [/math]
3 Velocidad del fluido
A continuación vamos a calcular los puntos en los que la velocidad es máxima sabiendo que la presión [math] p_1=2 [/math], la presión [math] p_2=1 [/math] y el coeficiente de viscosidad del fluido [math] \mu = 1 [/math].
Los puntos del fluido en que la velocidad es máxima serán aquellos en lo cuales el módulo del campo de velocidades ([math] \vec u [/math]) sea máximo:
[math] \vec u (x,y)= \frac{y(1-y)}{2}\vec i[/math]:
[math]| \vec u (x,y)|= \frac{y(1-y)}{2} [/math]
Para encontrar el máximo de la función del módulo de la velocidad, hacemos la derivada y vemos para qué valores se hace cero:
[math]| \vec u (x,y)|’= \frac{1}{2}-y [/math]:
[math]\frac{1}{2}-y=0 \implies y=\frac{1}{2}[/math]
Para ver si [math] y=\frac{1}{2}[/math] es un máximo o un mínimo, hacemos la segunda derivada:
[math]| \vec u (x,y)|’’= -1 [/math]
como la segunda derivada (que representa los cambios de la curvatura de la función) es negativa, sabemos que el valor es un máximo (el valor negativo indica una curvatura convexa, es decir, es un punto en el que los valores de su alrededor son menores que él).
Como podemos ver en la parábola, el punto [math] y=\frac{1}{2}[/math]( fluido de la parte central del canal) es el de velocidad máxima.
y=0:0.1:1; % creamos el vector y con valores entre 0 y 1 equidistantes 0.1
funcion = 1/2.* (y- y.^2) % creamos una función que represente la parábola
plot (y,funcion) % dibujamos la parábola
También comprobamos que se cumple la condición de que en las paredes la velocidad del fluido es nula, para nuestro caso las paredes son [math] y=0 , y=1 [/math]: [math]\vec u (x,0)= \frac{0(1-0)(p_{1}-p_{2})}{(2\mu)}\vec i = 0[/math]: [math]\vec u (x,1)= \frac{1(1-1)(p_{1}-p_{2})}{(2\mu)}\vec i = 0[/math] La viscosidad de un fluido hace que este se pegue a las paredes y por lo tanto que adquiera su velocidad. Como en nuestro caso las paredes no se mueven, la velocidad del fluido en contacto con ellas es cero, tal y como se demuestra en los cálculos anteriores.
4 Líneas de corriente del campo [math] \vec u [/math]
Vamos a dibujar ahora las líneas de corriente del campo [math]\vec{u}[/math], es decir, las líneas que son tangentes a [math]\vec{u}[/math] en cada punto. Para ello calculamos el campo [math]\vec{v}[/math] que en cada punto es ortogonal a [math]\vec{u}[/math]. Podemos comprobar que [math]\vec{v}[/math] es irrotacional (por ser [math]\vec{u}[/math] de divergencia nula) y tiene potencial escalar [math]\psi[/math], que se conoce como función de corriente de [math]\vec{u}[/math]. Vamos a calcular entonces [math]\psi[/math]:
[math]\vec{v}=\vec{k}\times \vec{u} =\frac{y(1-y)(p_1-p_2)}{2\mu}\vec{j}= \nabla \psi [/math]:
[math]\frac{ \partial \psi}{\partial y}=v_2=\frac{y(1-y)(p_1-p_2)}{2\mu} \implies\psi =\int\frac{y(1-y)(p_1-p_2)}{2\mu}dy =\frac{(p_1-p_2)}{2\mu}(\frac{y^2}{2}-\frac{y^3}{3})+ f(x) [/math]:
[math]\frac{ \partial \psi}{\partial x}= v_1=0\implies \frac{\partial}{\partial x}(\frac{(p_1-p_2)}{2\mu}(\frac{y^2}{2}-\frac{y^3}{3})) + f'(x)=0[/math]
[math]\psi =\frac{(p_1-p_2)}{2\mu}(\frac{y^2}{2}-\frac{y^3}{3}) + C[/math]:
Dibujaremos las líneas del potencial escalar [math]\psi =cte[/math] tomando [math] p_{1}=2[/math] [math] p_{2}=1[/math] y [math] \mu=1[/math] y considerando [math]C=0[/math]. Por tanto, [math]\psi =\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})[/math]
x=0:0.1:4; % creamos el vector x con valores entre 0 y 4 equidistantes 0.1
y=0:0.2:1; % creamos el vector y con valores entre 0 y 1 equidistantes 0.2
[xx,yy]=meshgrid(x,y); % matrices de las coordenadas x e y
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 XYDEMOSTRACIONES DE AFIRMACIONES DE ESTE APARTADO:
Si [math]\vec{u}[/math] es de divergencia nula, [math]\vec{v}[/math] es irrotacional:
[math]\nabla \cdot \vec u = 0\implies \frac{ \partial u_1 }{ \partial x} + \frac{\partial u_2 }{ \partial y} = 0 [/math]:
[math]\vec{v}=\vec{k}\times \vec{u}=det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ 0 & 0 & 1\\ u_1 & u_2 & u_3 \end{pmatrix} = -u_2\vec{i} +u_1\vec{j} [/math]:
[math]\nabla\times\vec{v}= det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ v_1 & v_2 & v_3 \end{pmatrix}= det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ -u_2 & u_1 & 0 \end{pmatrix}=-\frac{\partial}{\partial z}(u_1)\vec{i} - \frac{\partial}{\partial z}(-u_2)\vec{j} + (\frac{\partial}{\partial x}(u_1) -\frac{\partial}{\partial y}(-u_2))\vec{k} [/math]
[math]\frac{\partial}{\partial x}(-u_2)= -\frac{\partial}{\partial x}(u_2)[/math] por lo tanto: [math]\frac{\partial}{\partial x}(u_1) - \frac{\partial}{\partial y}(-u_2)= \frac{\partial}{\partial x}(u_1) + \frac{\partial}{\partial y}(u_2)=\nabla \cdot \vec u = 0 [/math] y así demostramos que si [math]\vec{u}[/math] es de divergencia nula, [math]\vec{v}[/math] es irrotacional.
La segunda parte del apartado nos dice que calcular el potencial escalar de [math]\vec{v}[/math] equivale a calcular la función de corriente de [math]\vec{u}[/math]. Vamos a ver por qué:
Sea
[math]\vec{u}(x,y) = u_1\vec{i} + u_2\vec{j}[/math] y [math]\vec{v}(x,y)=\vec{k}\times \vec{u}= v_1\vec{i} + v_2\vec{j} = -u_2\vec{i} +u_1\vec{j} [/math]
La función de corriente de [math]\vec{u}[/math] a la que llamaremos [math]\phi[/math] se calcula:
[math] u_1= \frac{ \partial \phi}{\partial y}[/math]: [math] u_2=-\frac{ \partial \phi}{\partial x}[/math]:
El potencial escalar de [math]\vec{v}[/math], anteriormente llamado [math]\psi[/math] lo calculamos:
[math] v_1= \frac{\partial \psi}{\partial x}[/math]: [math] v_2= \frac{\partial \psi}{\partial y}[/math]: Por la definición de [math]\vec{v}[/math] sabemos que [math]v_1=-u_2 [/math] y [math]v_2=u_1[/math]. Así pues:
[math]u_1=\frac{\partial \phi}{\partial y}= v_2=\frac{ \partial \psi}{\partial y}[/math]: [math]u_2=-\frac{\partial \phi}{\partial x}= -v_1=-\frac{ \partial \psi}{\partial x}[/math] Así pues, [math] \psi=\phi + C [/math] Tal y como queríamos demostrar.
5 Cálculo del rotacional [math] \vec u [/math]
Calcular el rotacional de [math]\vec{u}[/math]:
[math]\nabla\times\vec{u}= det\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}=det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \frac{y(1-y)(p_1-p_2)}{2\mu}& 0 & 0 \end{pmatrix}=\frac{\partial}{\partial y}(\frac{y(1-y)(p_1-p_2)}{2\mu})\vec{k} =\frac{ -(1-2y)(p_1-p_2)}{2\mu}\vec{k} [/math]
[math]|\nabla\times\vec{u}|=\frac{ (1-2y)(p_1-p_2)}{2\mu}[/math]
Para encontrar el máximo de la función del módulo del rotacional, hacemos la derivada y vemos para qué valores se hace cero: Como [math] p_1=2[/math], [math] p_2=1[/math] y [math] \mu=1 [/math]; así pues, caracterizamos el módulo para estos valores: [math]|\nabla\times\vec{u}|=\frac{ (1-2y)}{2} [/math]:
[math]y=0[/math]: [math]|\nabla\times\vec{u}|=|\frac{ (1-2y)}{2}|=|\frac{1}{2} |=\frac{1}{2}[/math]: [math]y=1[/math]: [math]|\nabla\times\vec{u}|=|\frac{ (1-2y)}{2}|=|- \frac{1}{2}| =\frac{1}{2}[/math] Por lo tanto, los puntos del canal con mayor rotacional serán los de las paredes. Esto es razonable porque las partículas del fluido del centro del canal tienen mayor velocidad.
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); % matrices de las coordenadas x e y
s=abs (1/2.*(1-2.*yy)) %definimos la función que representa el módulo del rotacional
surf(xx,yy,s) %dibujamos la superficie
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view (2) % para visualizar en el plano XYLa intuición nos dice que el fluido avanza en línea recta y esta se ve confirmada al ver las líneas de corriente de [math]\vec u [/math]. Sin embargo, si analizados de manera local el rotacional en las partículas cercanas a las paredes del canal parece que estas giran. Esto se debe a que las velocidad de las partículas de la parte central es mayor y la de las paredes es nula lo que "induciría" un giro.
6 Temperatura del fluido
Supongamos ahora que la temperatura del fluido viene dada por el campo escalar [math]T(x,y)=(x-1)^2– y^2[/math]. Vamos a verlo representado gráficamente:
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); % matrices de las coordenadas x e y
t=((xx-1).^2) -(yy.^2) %definimos en la variable t el campo de temperaturas
surf(xx,yy,t) %dibujamos la superficie
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
Se puede observar como la temperatura alcanza un mínimo hacia x=1 y luego crece conforme avanza hacia la derecha
Ahora vamos vamos a calcular el gradiente del campo de temperaturas para después representarlo. En él podremos ver las direcciones de máximo crecimiento:
[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] Para representarlo:
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); % matrices de las coordenadas x e y
fi=2.*(xx-1); %componente i del campo vectorial
fj=-2.*yy; %componente j del campo vectorial
quiver(xx,yy,fi,fj)%dibujamos el campo gradiente de temperaturas
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XY
Se ve como las temperaturas crecen más fuertemente hacia la parte final de avance del fluido
Por último, vamos a comprobar gráficamente que el gradiente de temperatura es ortogonal a las curvas de nivel de la temperatura
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); % matrices 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 poner encima el gradiente de temperaturas
fi=2.*(xx-1); %componente i del campo vectorial
fj=-2.*yy; %componente j del campo vectorial
quiver (xx,yy,fi,fj) %dibujamos el campo gradiente de temperaturas
axis([0,4,-1,2]) %fijamos los ejes en [0,4]x[-1,2]
view(2) %para visualizar el campo en el plano XYSe puede apreciar como el gradiente de temperaturas es ortogonal a las curvas de nivel, que se corresponden con valores del campo de temperaturas.
