Visualización de campos escalares y vectoriales en fluidos (grupo A7)

De MateWiki
Saltar a: navegación, buscar

En el siguiente articulo vamos a estudiar el flujo de un fluido incompresible alrededor de un obstáculo circular. Trabajaremos con coordenadas cilindricas (polares ya que se reduce al plano) ya que debido a la forma del obstáculo los cálculos resultarán más sencillos. Las gráficas que se mostraran aqui han sido obtenidas con el programa MATLAB.

Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en fluidos. Grupo 7-A
Asignatura Teoría de Campos
Curso 2015-16
Autores Enrique Pellico Martín, Eduardo Moreda Meleiro, Rommel Beltran Carrero, Antonio Maya Hidalgo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En primer lugar debemos hacer referencia al concepto de fluido incompresible. Un fluido incompresible es aquel fluido cuya densidad permanece constante con el tiempo y se opone a la compresión bajo cualquier condición. Esto quiere decir que tanto su masa como su volumen permanecerán constantes. Su significado matemático es que dicho fluido tendrá una densidad constante, facilitando los cálculos necesarios.

2 Región ocupada por el fluido

El primer paso para la la visualización de campos escalares y vectoriales de un fluido es delimitar la región de estudio. La región de estudio es una corona circular de radio menor 2 y radio mayor 6.

La representamos mediante un mallado que represente los puntos interiores de la corona circular en el intervalo [−5, 5] × [−5, 5]

 1 %rho=u; phi=v;
 2 %Creamos el mallado.
 3 rho=linspace(2,6,50);
 4 th=linspace(0,2*pi,50);
 5 %generamos la reticula
 6 [U,V]=meshgrid(rho,th);
 7 %definimos la parametrización
 8 X=U.*cos(V);
 9 Y=U.*sin(V);
10 %Dibujamos la malla y definimos los ejes
11 mesh(X,Y,0*X)
12 axis([-5,5,-5,5]);
13 axis equal
Región que ocupa el fluido

3 Función potencial y campo de velocidades

La velocidad de un fluido viene dada por el campo de velocidades de este aplicado en cada punto.

El campo de velocidades deriva de un potencial escalar. Vamos a representar primero el campo escalar y después el campo de velocidades

3.1 Función Potencial

El campo escalar del que deriva es el siguiente:

Formula1A7.png

Se aprecia como la función potencial crece de izquierda a derecha según la dirección de movimiento del fluido

 1 %Dibujamos la malla 
 2 rho=linspace(2,6,40);
 3 th=linspace(0,2*pi,40);
 4 [U,V]=meshgrid(rho,th);
 5 %dibujamos la función potencial
 6 %primero la superficie
 7 X=U.*cos(V);
 8 Y=U.*sin(V);
 9 %Despues definimos la función
10 phi=inline('2*(rho+(4./rho)).*cos(th)','rho','th');
11 %aplicamos la función
12 Z=phi(U,V);
13 %Dibujamos
14 p=surf(X,Y,Z);
15 axis equal
16 view(2)
17 hold on
18 %por último dibujamos el obstáculo
19 plot(2*cos(th),2*sin(th),'k','linewidth',2)
Representación función potencial

3.2 Campo de velocidades

El campo de velocidades viene dado por el gradiente de la función potencial [math]\varphi[/math]

Formula2A7.png


Formula3A7.png


Su representación viene dada por:

 1 %Dibujamos la malla 
 2 rho=linspace(2,6,30);
 3 th=linspace(0,2*pi,30);
 4 [U,V]=meshgrid(rho,th);
 5 %dibujamos la función potencial
 6 %primero la superficie
 7 X=U.*cos(V);
 8 Y=U.*sin(V);
 9 %Despues definimos la función
10 phi=inline('2*(rho+(4./rho)).*cos(th)','rho','th');
11 %aplicamos la función
12 Z=phi(U,V);
13 %Dibujamos las curvas de nivel
14 contour(X,Y,Z,40)
15 %axis([-5,5,-5,5]);
16 axis equal
17 view(2)
18 hold on
19 %Dibujamos el campo de velocidades del fluido:
20 %aplicamos las funciones
21 ux=(2-8./(U.*U)).*cos(V).*cos(V)+U.*sin(V).*sin(V).*(2./U+8./(U.*U.*U));
22 uy=(2-8./(U.*U)).*cos(V).*sin(V)-U.*sin(V).*cos(V).*(2./U+8./(U.*U.*U));
23 %dibujamos el campo de velocidades
24 quiver(X,Y,ux,uy);
25 axis equal
26 view(2)
27 %por último dibujamos el obstáculo
28 plot(2*cos(th),2*sin(th),'k','linewidth',2)
29 hold off
Campo de velocidades [math]\vec{u}[/math]

Se ve perfectamente como el fluido evita el obstáculo rodeándolo. Disminuye de velocidad al chocar contra el obstáculo y aumenta al evitarlo.

También vamos a realizar un Zoom sobre la imagen para apreciar en detalle como el campo de velocidades [math]\vec{u}[/math] es ortogonal a las curvas de nivel de la función potencial [math]\varphi[/math] del que proviene.

Las curvas de nivel de la función potencial representan los puntos del plano en el que esta tiene el mismo valor. Según se ve en el apartado anterior, en la representación gráfica se observa que la función crece de izquierda a derecha evitando el obstáculo. El campo de velocidades viene dado por el gradiente de esa función. El gradiente representa el valor y la dirección de máximo crecimiento de de la función escalar. Como en las curvas de nivel la función es constante, es decir, en las curvas la variación es 0, la dirección de máximo crecimiento es ortogonal a las curvas de nivel y viene dada por el gradiente, es decir por el campo de velocidades del fluido.

Zoom

3.2.1 Campo de velocidades respecto al obstáculo

Si consideramos el obstaculo como la superficie circular de centro el origen y radio 2, llamamos [math]\vec{n}[/math] al vector normal a esta superficie, es decir a los puntos del obstáculo comprobamos que [math]\vec{n}\cdot\vec{u}=0[/math]

Como el vector normal a la superficie circular es:

Formula5A7.png

Y el campo de velocidades [math]\vec{u}[/math] es:

Formula4A7.png

Como [math]\vec{g_{\rho}}[/math] y [math]\vec{g_{\theta}}[/math] son ortogonales al vector [math]\vec{g_{z}}[/math] el producto escalar [math]\vec{n}\cdot\vec{u}=0[/math] es igual a 0.

Esto quiere decir que el fluido en todo momento las particulas del fluido al chocar con el obstáculo cambian de dirección de forma que rodean el obstáculo, y nunca lo pasan por encima

3.2.2 Campo de velocidades lejos del obstáculo

Si nos encontramos lejos del obstáculo el valor de [math]\rho[/math] es muy grande por lo que podemos considerar que [math]\frac{1}{ \rho }=0[/math] por lo tanto el campo de velocidades del fluido queda:


Formula4A7.png


Formula6A7.png


Por lo tanto como [math]|\vec{g_{\rho}}|=1[/math] entonces la velocidad del fluido en los puntos alejados del obstáculo es [math]|\vec{u}|= 2cos( \theta ) [/math]

4 Divergencia y rotacional

Con el estudio de la divergencia y de rotacional se demuestra la incompresibilidad del fluido. Al tratarse de una región simplemente conexa y con el campo de velocidades que tiene potencial escalar, podemos afirmar que el rotacional del campo será nulo: El rotacional en coordenadas cilíndricas:

Rotacional.PNG

Resolviendo:

Resolucion rotacional.PNG

La divergencia tendrá que ser nula demostrando así la imposibilidad de variación de volumen, como se interpreta de este operador.

Definimos la divergencia para coordenadas cilíndricas:

Divergencia.JPG

Resolviendo:

Divergencia resuelta.JPG

Observamos que igualmente se anula. Con la divergencia nula se demuestra la incompresibilidad del fluido como se había predicho.

5 Líneas de corriente del campo

Las líneas de corriente son líneas que representadas en el plano definen la trayectoria del campo, en nuestro caso, de las partículas de agua.

Son, a su vez, el potencial de un campo vectorial ortogonal al campo de velocidades. Luego para definir las líneas de corriente se calculará el campo ortogonal y su potencial.

5.1 Campo ortogonal al campo de velocidades

Este campo obtiene multiplicando escalarmente el vector de posicion [math]\vec{k}[/math] con el campo [math]\vec{u}[/math]. Para facilitar el proceso utilizamos el vector de posición de la coordenada curvilinea en z:

Campov.JPG

5.1.1 Comprobación gráfia

  • Campo v en coordenadas cartesianas
  • Código MATLAB
 1 %Reticula
 2 h=0.3;
 3 rho=2:h:6;
 4 th=0:h:2*pi+h;
 5 [U,V]=meshgrid(rho,th);
 6 
 7 %Cambio de coordenadas a cartesianas
 8 X=U.*cos(V);
 9 Y=U.*sin(V);
10 
11 %Componentes del campo de velocidades en cartesianas
12 ux=(2-8./(U.*U)).*cos(V).*cos(V)+U.*sin(V).*sin(V).*(2./U+8./(U.*U.*U));
13 uy=(2-8./(U.*U)).*cos(V).*sin(V)-U.*sin(V).*cos(V).*(2./U+8./(U.*U.*U));
14 
15 %Componentes del campo v en cartesianas
16 vx=(2+8./(U.*U)).*sin(V).*cos(V)-(2./U-8./(U.*U.*U)).*cos(V).*sin(V).*U;
17 vy=(2+8./(U.*U)).*sin(V).*sin(V)+(2./U-8./(U.*U.*U)).*cos(V).*cos(V).*U;
18 
19 %Dibujo campo de velocidades
20 quiver(X,Y,ux,uy)
21 axis([-5,5,-5,5])
22 axis equal
23 view(2)
24 hold on
25 
26 %Dibujo campo v
27 quiver(X,Y,vx,vy)
28 axis([-5,5,-5,5])
29 axis equal
30 view(2)
Campos u y v

5.1.2 Observación del rotacional

Al ser u de divergencia nula el rotacional de [math]\vec{v}[/math] será nulo:

Rotav.JPG

5.2 Potencial escalar de corriente

Una vez obtenido el campo [math]\vec{v}[/math] podremos interpretar las líneas de corriente como la función potencial de este campo, luego:

Potencialc.JPG

Calculamos el potencial escalar:

Cpotencial1.JPG

Donde C igualando:

Cpotencial3.JPG

5.3 Representación de las líneas de corriente

Una vez tenemos el potencial escalar de líneas de corriente se puede representar de la siguiente manera:

 1 %Reticula
 2 h=0.3;
 3 rho=2:h:6;
 4 th=0:h:2*pi+h;
 5 [U,V]=meshgrid(rho,th);
 6 
 7 %Cambio de coordenadas a cartesianas
 8 X=U.*cos(V);
 9 Y=U.*sin(V);
10 
11 %Potencial escalar. Lineas de corriente
12 chi=(2*U-8./U).*sin(V);
13 
14 %Componentes del campo de velocidades en cartesianas
15 ux=(2-8./(U.*U)).*cos(V).*cos(V)+U.*sin(V).*sin(V).*(2./U+8./(U.*U.*U));
16 uy=(2-8./(U.*U)).*cos(V).*sin(V)-U.*sin(V).*cos(V).*(2./U+8./(U.*U.*U));
17 
18 %Dibujo lineas de corriente
19 contour(X,Y,chi,50,'LineWidth',1.5)
20 axis([-5,5,-5,5])
21 axis equal
22 view(2)
23 hold on
24 
25 %Dibujo campo de velocidades
26 quiver(X,Y,ux,uy)
27 axis([-5,5,-5,5])
28 axis equal
29 view(2)
30 hold on
Líneas de corriente y campo de veocidades

5.3.1 Observaciones

  1. Podemos observar que la velocidad es tangente a las líneas de corriente realizando un zoom en la representación anterior
Velocidad y líneas de corriente. Zoom
  1. Observamos la ortogonalidad de las líneas de corriente con el campo equipotencial ejecutando el código MATLAB posterior.Ampliando observamos con detalle la ortogonalidad.
 1 %Reticula
 2 h=0.3;
 3 rho=2:h:6;
 4 th=0:h:2*pi+h;
 5 [U,V]=meshgrid(rho,th);
 6 
 7 %Cambio de coordenadas a cartesianas
 8 X=U.*cos(V);
 9 Y=U.*sin(V);
10 
11 %Potencial escalar. Lineas de corriente (chi).
12 chi=(2*U-8./U).*sin(V);
13 %Cuervas equipotenciales(po)
14 po=(2*U+8./U).*cos(V);
15 
16 %Dibujo lineas de corriente
17 contour(X,Y,chi,50,'LineWidth',1.5)
18 axis([-5,5,-5,5])
19 axis equal
20 view(2)
21 hold on
22 
23 %Dibujo equipotenciales
24 contour(X,Y,po,50,'LineWidth',1.5)
25 axis([-5,5,-5,5])
26 axis equal
27 view(2)
Equipotenciales y líneas de corriente. Zoom

6 Velocidad máxima y mínima del fluido

Los máximos y los mínimos, en el contorno del obstáculo, se pueden calcular maximizando o minimizando el módulo del campo vectorial [math]\vec{v}[/math], que es más sencillo, ya que sabemos que este campo tiene el mismo módulo que [math]\vec{u}[/math] (fácil de demostrar a partir de la definición del módulo del producto escalar).

Modulov.JPG

Donde el módulo de [math]\vec{v}[/math] sea máximo la velocidad será máxima y donde sea mínimo la velocidad será mínima:

Los puntos de velocidad máxima son: [0,2] y [0,-2].
Los puntos de velocidad mínima son: [-2,0] y [2,0].
  • Puntos de remanso

Son aquellos puntos donde la velocidad del fluido es nula. Observamos que coinciden con los puntos de velocidad mínima. En [-2,0] y [2,0] existiran dos puntos de remanso.

7 Presión del fluido. Ecuación de Bernouilli

7.1 Ecuación de Bernouilli

Como sabemos, un fluido se caracteriza por por no describir siempre una forma concreta, sino que su forma o dinámica depende principalmente de agentes externos (en el caso de la forma, depende de un recipiente). Aquí vamos a estudiar el comportamiento dinámico de un fluido, dependiendo de la presión a la que esté sometido y la densidad del mismo, ya que su movimiento no depende del tiempo, por lo que es estacionario. Así pues, se dedujo y describe la ecuación de Bernouilli:

Ecuacion bernouilli.jpg

Con nuestro campo vectorial u y teniendo como datos d=2 y cte=15, podremos obtener la presión del fluido:

Formula1 bernouilli.jpg
Formula2 bernouilli.jpg





De esta manera, operamos y nos queda la presión:

Formula3 bernouilli.jpg

Una vez hallada, procedemos a representarla gráficamente. Para ello, empleamos el siguiente comando en Matlab, dando:

 1 % vemos variables
 2 h=0.1;
 3 rho=2:h:6;
 4 th=0:h:2*pi+h
 5 
 6 %hacemos retícula
 7 [U,V]=meshgrid(rho,th);
 8 
 9 %describo coordenadas
10 X=U.*cos(V);
11 Y=U.*sin(V);
12 
13 %aplicamos funcion presion
14 f=15-((4*cos(V).^2)).*((1-(4./(U.^2))))-((4*((sin(V)).^2)).*(((1./U)+(4./(U.^3))).^2));
15 
16 %dibujo
17 surf(X,Y,f)
18 axis([-5,5,-5,5])
19 axis equal
20 view(2)
Campo de presiones


Con esto, nos damos cuenta de que las mayores presiones se concentran en los extremos superior e inferior del obstáculo, ya que es el momento posterior al choque del fluido contra él, y el desplazamiento a los lados. En cambio, las menores presiones se distinguen a los laterales izquierdo y derecho, siendo la dirección a la que se desplaza el fluido y donde no tiene contacto con el obstáculo, debiéndose así que sus presiones sean mucho menores.

7.2 Campo de presiones y lineas de corriente

La conclusión sacada anteriormente se puede corroborar observando las líneas de corriente del fluido, ya que a su vez representan la velocidad del mismo y los cambios que se producen cuando impacta con el obstáculo.

Lineas de corriente
Campo de presiones



Al comparar las dos gráficas, podemos ver que, efectivamente, en los extremos superior e inferior del obstáculo hay una mayor concentración de corrientes, y sus velocidades son mucho mayores, generando presiones significativamente altas. En cambio, en los lados izquierdo y derecho todo está `más tranquilo´ porque no hay nada que se interponga a la corriente normal del fluido. Por tanto, si fuéramos una partícula del fluido, inicialmente podríamos comenzar moviéndonos en cotas inferiores, frías, pero al llegar al obstáculo, se crea una barrera que nos desplazaría a los extremos, donde se acumularían muchas más partículas, generando más presión y haciendo que el rozamiento rápido entre éstas llevara a un mayor calor, impulsándonos a cotas más superiores en ese tramo. Posteriormente, una vez pasado el obstáculo, iríamos volviendo a nuestro estado inicial.

8 Ecuación de Navier-Stokes

Las ecuaciones de Navier-Stokes tratan de describir el movimiento de un fluido en la atmósfera, las corrientes oceánicas, y sobre cualquier objeto, como proyectiles. Está publicado como uno de los Siete Grandes Problemas del Milenio, ya que se busca su generalización para cualquier cosa donde estén involucrados los fluidos newtonianos. Actualmente se ha demostrado pero para casos muy particulares. De manera simplificada, la ecuación satisface, con los parámetros de campo u y presión ,que:

Ecuacion stokes.jpg


Nosotros vamos a comprobar que se va a cumplir en el caso de que la viscosidad del fluido sea nula, y donde además se verifica que es incompresible, debido a que siempre ocupará el mismo volumen. Así pues, quedaría demostrado que el gradiente por el campo u es nulo. Con todo, la ecuación de Navier-Stokes resultante es la siguiente:

Ecuacion1 stokes.jpg


Sabiendo que d es la densidad del fluido, y por tanto una constante, y que:
Ecuacion222 stokes.jpg

Desarrollamos más profundamente los cálculos:

Ecuacion33 stokes.jpg




Como podemos observar, nos queda una ecuación determinada por coordenadas cartesianas. Estudiándola detenidamente, y relacionando las leyes físicas que conllevaron a su elaboración, nos damos cuenta de que se asemeja a la ecuación de Bernouilli, ya que en las dos se menosprecia la viscosidad y tratan con los campos vectoriales del fluido y sus presiones. Así pues, al relacionarlas con las mismas condiciones matemáticas, podríamos establecer la igualdad de las mismas. Para ello, le aplicamos el gradiente a la ecuación de Bernouilli:

Ecuacion444 stokes.jpg
Ecuacion33 stokes.jpg






El resultado final equivale al de la ecuación de Navier-Stokes, por lo que si se cumple una, la otra también, y como la de Bernouilli vimos que sí, poodemos afirmar que Navier-Stokes se cumple cuando la viscosidad se considera nula.

9 Teorema de Kutta-Joukowski. Paradoja de D'Alembert

En primer lugar haremos una breve explicación de la Paradoja de D'Alembert. Dicha paradoja, es una contradicción a la que llegó D'Alembert después de estudiar matemáticamente el fenómeno de la resistencia producida sobre un cuerpo cuando una corriente de fluido circula sobre él.

D'Alembert aplicó la teoría de flujo potencial para modelar el fenómeno, y concluyó que la fuerza resultante sobre el cuerpo sobre el cual fluye el aire es cero, lo cual se contradice con la observación.

Por otra parte según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo es proporcional a la circulación.

A continuación comprobaremos la circulación a lo largo de la circunferencia de radio 2. En primer lugar parametrizamos dicho obstaculo:

GammaA7.png



Calculamos la circulacion en el intervalo (a,b)=(0,2π)

Circulacion1A7.png




Sustituyendo nuestros datos dentro de la integral obtenemos:

Circulacion4A7.png




Finalmente desarrolando dicha integral:

Circulacion3A7.png




Obtenemos el resultado de que la circulación es nula. Esto va en contra de la intuición ya que se supone que el fluido deberia ejercer una fuerza sobre el obstáculo, pero aplicando el Teorema de Kutta-Joukowski, que como hemos dicho antes dice que la fuerza es proporcional a la circulación, llegamos a la conclusión de que la fuerza deberia ser nula también. Encontrandonos con la conocida como Paradoja de D'Alembert.

10 Curvas de nivel de la presión

A continuación observamos el código MATLAB utilizado para obtener el gráfico de las curvas de presión del fluido.

 1 %Reticula
 2 ro=2:0.1:6;
 3 teta=0:0.1:2*pi;
 4 [U,V]=meshgrid(ro,teta);
 5 %Cambio de coordenadas
 6 X=U.*cos(V);
 7 Y=U.*sin(V);
 8 %Formula de la presión
 9 p=15-((4*((cos(V)).^2)).*((1-(4./(U.^2))).^2))-((4*((sin(V)).^2)).*(((1./U)+(4./(U.^3))).^2));
10 %Dibujo de las curvas de nivel
11 contour(X,Y,p,50,'linewidth',1.5)
12 axis([-5,5-5,5])
13 axis equal

Como anteriormente hemos observado los maximos valores para la presión los encontramos en la zona superior e inferior del óbstaculo correspondiendose con las zonas de mayor velocidad. Mientras que el los laterales las líneas son de un color azul que equivale a valores más bajos de la presión, justo donde las corrientes se alejan del óbstaculo.

Curvas de nivel de la presión del fluido

11 Presión media del fluido

Vamos a calcular la presión media del fluido. Para ello se hará un aaproximacion de la integral de la presion y lo dividiremos entre el área de la región que ocupa el fluido. En este caso es el área de una corona circularde radios 2 y 6.

AreaA7.png



Realizamos la aproximacion de la presión:

 1 %Reticula
 2 h=0.01;
 3 ro=2:h:6;
 4 teta=0:h:2*pi;
 5 [U,V]=meshgrid(ro,teta);
 6 %Presión del fluido
 7 p=15-((4*((cos(V)).^2)).*((1-(4./(U.^2))).^2))-((4*((sin(V)).^2)).*(((1./U)+(4./(U.^3))).^2));
 8 %Integral de la presión
 9 pres=U.*p;
10 %Integral
11 volumen=h^2*pres;
12 w=sum(sum(volumen));
13 %Área de la corona circular
14 area=(pi*6^2)-(pi*2^2);
15 %Presión media
16 presmedia=w/area

Realizando estos calculos obtenemos un resultado de 13.6441 que si contrastamos con una de las gráficas de la presión podemso ver que concuerda, ya que en la gráfica a simple vista podemos observar como los valores medios rondan este valor.

PresionA7.png