Visualización de campos escalares y vectoriales en fluidos (Grupo 12A)

En el artículo desarrollado a continuación vamos a estudiar el flujo de un fluido incompresible alrededor de un obstáculo con forma circular. Vamos a trabajar con coordenadas cilíndricas (polares ya que se reduce al plano) ya que con la forma del obstáculo, los cálculos nos resultarán más sencillos. Las gráficas que mostremos aquí han sido obtenidas con el programa MATLAB / OCTAVE.

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en fluidos. Grupo 12A
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Jorge Barrena Blázquez, Rommel Beltrán Carrero, Pablo Domingo Escudero Bravo, Carlos Maroto Escribano
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Antes de comenzar, cabe explicar el concepto de fluido incompresible. Un fluido incompresible es aquel cuya masa y volumen no variarán, permaneciendo constantes en el tiempo; y se opone a las compresiones ofrecidas bajo cualquier condición o circunstancia. Su significado matemático, por tanto, es que dicho fluido tendrá una densidad constante, facilitando los cálculos necesarios por homogeneidad de sus propiedades en cualquier posición.

2 Región ocupada por un fluido

De forma previa a la realización de cálculo alguno respecto al comportamiento del fluido habrá que delimitar la región del estudio. Se nos da una corona circular de radios 1 y 5, por lo que representamos su mallado en el intervalo [-5,5] x [-5,5]:

izquierda

izquierda

3 Función potencial y su campo de velocidades

La velocidad de las partículas del fluido viene dada por el gradiente de :

La representación gráfica de la función potencial y su gradiente (campo de velocidades) es la siguiente:

izquierda

derecha

derecha

izquierda

Ampliando la imagen del campo de velocidades puede apreciarse que es ortogonal a la función potencial. Además, para mayor claridad se han superpuesto los gráficos de la función potencial y del campo de velocidades.

izquierda

derecha

3.1 Velocidades paralelas a la superficie

Vamos a experimentar con las direcciones de la velocidad al entrar en contacto con el obstáculo. Probando a hacer el producto escalar entre la velocidad del fluido y el vector normal del osbtáculo obtenemos:

Con lo que puede deducirse que la velocidad del fluido siempre será perpendicular al vector normal de la superficie, es decir, que las partículas del fluido únicamente rodearán el espacio, no lo pasarán por encima, pues las componentes verticales de la velocidad se anulan.

3.2 Campo de velocidades y función potencial lejos del obstáculo

Vamos a ver qué sucede cuando nos alejamos mucho de nuestro obstáculo y entonces \(\rho\) es muy grande. Para ello vamos a hacer el límite de cuando \(\rho\) tiende a infinito de nuestra función potencial \(\varphi= \big( \rho + \frac{1}{ \rho } \big) sen \big( \theta \big)\) veremos como \(\frac{1}{ \rho }\) será igual a cero.

\(\lim_{ \rho \rightarrow \infty }\big( \rho + \frac{1}{ \rho } \big) sen \big( \theta \big)=\rho sen \big( \theta \big)\)

Y nuestro gradiente

\(\nabla\varphi = \frac{\partial\varphi}{\partial \rho } g_{p}+ \frac{1}{ \rho^{2} } \frac{\partial \varphi }{\partial \theta } g_{ \theta } + \frac{\partial \varphi }{\partial z} g_{z} =sen \big( \theta \big) g_{ \rho } + \frac{1}{ \rho } cos \big( \theta \big) g_{ \theta }=sen \big( \theta \big) g_{ \rho } \)

Y por ello, este campo de velocidades se dará en los puntos donde \(\rho\) es muy grande.

4 Rotacional y divergencia del campo

En este punto estudiaremos analiticamente que \( \overline{u}\) su rotacional como su divergencia es nula y trataremos el tema de la incompresibilidad del fluido.

Nuestro rotacional en coordenadas polares:

\(\nabla \times \overline{u} = \frac{1}{ \rho } \begin{bmatrix} g_{ \rho } & g_{ \theta }& g_{ z }\\ \frac{\partial}{\partial\rho } & \frac{\partial}{\partial \theta } & \frac{\partial}{\partial z} \\ u g_{ \rho } & ug_{ \theta } & ug_{ z} \end{bmatrix}= \frac{1}{ \rho } \begin{bmatrix} g_{ \rho } & g_{ \theta }& g_{ z }\\ \frac{\partial}{\partial\rho } & \frac{\partial}{\partial \theta } & \frac{\partial}{\partial z} \\ sen \big( \theta \big) \big(1- \frac{1}{ \rho ^{3} } \big) & cos \big( \theta \big) \big( \rho + \frac{1}{ \rho ^{2} } \big) & 0 \end{bmatrix}= \frac{1}{ \rho } \big\{ \big(cos \big( \theta \big)- \frac{1}{ \rho ^{2} } cos \big( \theta \big)\big) - \big(cos \big( \theta \big)-cos \big( \theta \big) \frac{1}{ \rho ^{2} } \big) \big\} g_{z} =0\)

Nuestra divergencia en coordenadas polares:

\(\nabla . \overline{u} = \frac{1}{ \sqrt{g} } \frac{\partial}{\partial x_{i} } \big( \sqrt{g} u_{i} \big) = \frac{1}{ \rho } \big\{ \frac{\partial}{\partial \rho } \big(\rho sen \big( \theta \big) -sen \big( \theta \big) \frac{1}{\rho} \big) +\frac{\partial}{\partial \theta } \big(cos \big( \theta \big)+cos \big( \theta \big) \frac{1}{ \rho ^{2} } \big) \big\}= \)

\(\frac{1}{ \rho } \big\{ \frac{\partial}{\partial \rho } \big(\rho sen \big( \theta \big) -sen \big( \theta \big) \frac{1}{\rho} \big) +\frac{\partial}{\partial \theta } \big(cos \big( \theta \big)+cos \big( \theta \big) \frac{1}{ \rho ^{2} } \big) \big\}=\)

\(\frac{1}{ \rho } \big\{sen \big( \theta \big)+sen \big( \theta \big)\frac{1}{ \rho ^{2} }-sen \big( \theta \big) -sen\big( \theta \big) \frac{1}{ \rho ^{2} } \big\} =0\)

Vemos que ambos quedan reducidos a 0.

Demostrando que la divergencia es nula, vemos que nuestro fluido es un líquido incompresible.

5 Líneas de corriente del campo

En este punto vamos a dibujar las líneas de corriente que nos definirán la trayectoria del plano y son tangentes a \( \overline{u}\). Para ello necesitaremos calcular el campo ortogonal y su potencial.

5.1 Calculo analítico

Este campo lo obtendremos del resultado de multiplicar el vector \( \overline{k}\) con el campo \( \overline{u}\). Que al estar en coordenadas polares nuestro \( \overline{k}\) será \(\overline{ g_{z} }\).

\(\overline{v} = \overline{k} \times \overline{u} =\overline{g_{z}} \times \overline{ u }=\begin{bmatrix} g_{ \rho } & g_{ \theta }& g_{ z }\\ 0 & 0 & 1 \\ sen \theta \big(1- \frac{1}{ \rho ^{2} } \big) & \rho ^{2}cos \theta \big( \frac{1}{ \rho } + \frac{1}{ \rho ^{3} } \big) & 0 \end{bmatrix}=\frac{1}{ \rho }\big\{-\rho^{2} cos \theta \big( \frac{1}{ \rho } + \frac{1}{ \rho ^{3} } \big) \overline{g _{ \rho } } +sen \theta\big( 1-\frac{1}{ \rho ^{2} } \big) \overline{ g_{ \theta } } \big\} =-cos \theta \big( 1-\frac{1}{ \rho ^{2} } \big) \overline{ g_{ \rho } }+ sen \theta \big( \frac{1}{ \rho } + \frac{1}{ \rho ^{3} } \big) \overline{ g_{ \theta } }\)

5.2 Comprobación por grafica de matlab

clear all; clf;
%datos
p=linspace(1,5,25);
theta=linspace(0,2*pi,25);
%creamos la malla
[U,V]=meshgrid(p,theta);
%Le damos la forma que queremos
X=U.*cos(V);
Y=U.*sin(V);
%doy la función
phi=inline('(p+(1./p)).*sin(theta)','p','theta');
%consehuimos las alturas
Z=phi(U,V);
%el grafico
%fibujamoslas lineas de nivel
% contour(X,Y,Z,75);
% axis equal
 
hold on
%las derivadas parciales
ux=(1-1./(U.*U)).*sin(V).*cos(V)-U.*sin(V).*cos(V).*((1./U)+(1./(U.*U.*U)));
uy=(1-1./(U.*U)).*sin(V).*sin(V)+U.*cos(V).*cos(V).*((1./U)+(1./(U.*U.*U)));
%dibujamos el gradiente
quiver(X,Y,ux,uy);
%el ortogonal
vx=-(1+1./(U.*U)).*cos(V).*cos(V)-U.*sin(V).*sin(V).*((1./U)-(1./(U.*U.*U)));
vy=-(1+1./(U.*U)).*cos(V).*sin(V)+U.*cos(V).*sin(V).*((1./U)-(1./(U.*U.*U)));
%dibujar el otrtogonal
quiver(X,Y,vx,vy);

centro

5.3 Demuestro que este campo es irrotacional

Al ser de divergencia nula el rotacional también lo será:

\(\nabla \times \overline{v} = \frac{1}{ \rho } \begin{bmatrix} g_{ \rho } & g_{ \theta }& g_{ z }\\ \frac{\partial}{\partial\rho } & \frac{\partial}{\partial \theta } & \frac{\partial}{\partial z} \\ v g_{ \rho } & vg_{ \theta } & vg_{ z} \end{bmatrix}= \frac{1}{ \rho } \big\{sen \theta \big(1- \frac{1}{ \rho ^{2} } \big)- \big(sen \theta \big(1- \frac{1}{ \rho ^{2} } \big) \big) \big\} =0 \)

Así queda demostrado que es irrotacional.

5.4 Potencial escalar de corriente

Interpretaremos las líneas de corriente como la función potencial de este campo \(\overline{v} = \nabla \psi\)

Calcularemos el potencial escalar:

\(\frac{\partial \psi }{\partial \rho } =-cos \theta \big(1- \frac{1}{ \rho ^{2}} \big) \rightarrow \psi =- \int cos \theta d \rho + \int \frac{cos \theta }{ \rho ^{2} }d \rho \rightarrow \psi = - \rho cos \theta - \frac{1}{ \rho }cos \theta+C \big( \theta \big)\)

Donde C igualando \(\rightarrow cos \big( \theta \big)=cte \rightarrow 0\)

5.5 Representamos las líneas de corriente con Matlab

El código respectivo:

clear all; clf;
%datos
p=linspace(1,5,50);
theta=linspace(0,2*pi,50);
%creamos la malla
[U,V]=meshgrid(p,theta);
%Le damos la forma que queremos
X=U.*cos(V);
Y=U.*sin(V);
%potencial escalar. Lineas de corriente
chi=(1./U-U).*cos(V)
 
%las derivadas parciales
ux=(1-1./(U.*U)).*sin(V).*cos(V)-U.*sin(V).*cos(V).*((1./U)+(1./(U.*U.*U)));
uy=(1-1./(U.*U)).*sin(V).*sin(V)+U.*cos(V).*cos(V).*((1./U)+(1./(U.*U.*U)));
%Dibujo lineas de corriente
contour(X,Y,chi,50,'LineWidth',1.5)
 axis([-5,5,-5,5])
 axis equal
 view(2)
 hold on
 
%Dibujo campo de velocidades
 quiver(X,Y,ux,uy)
axis([-5,5,-5,5])
 axis equal
view(2)
hold on

centro

En ella se ve como la velocidad es tangente a las líneas de corriente.

5.6 Las lineas de corriente son ortogonales a las curvas equipotenciales

Se presenta gráficamente que son ortogonales:

clear all; clf;
%datos
p=linspace(1,5,50);
theta=linspace(0,2*pi,50);
%creamos la malla
[U,V]=meshgrid(p,theta);
%Le damos la forma que queremos
X=U.*cos(V);
Y=U.*sin(V);
%potencial escalar. Lineas de corriente
chi=(1./U-U).*cos(V)
%curvas equipotenciales
po=(1./U+U).*sin(V)
 %Dibujo lineas de corriente
 contour(X,Y,chi,50,'LineWidth',1.5)
 axis([-5,5,-5,5])
 axis equal
 view(2)
 hold on
 
%Dibujo equipotenciales
 contour(X,Y,po,50,'LineWidth',1.5)
 axis([-5,5,-5,5])
 axis equal
 view(2)

centro

6 Velocidades máxima y mínima en frontera del obstáculo

Como hemos demostrado anteriormente, el fluido va a rodear al obstáculo S, sin llegar a pasar por encima ni atravesarlo. Ahora bien, sí que va a rodearlo, encontrando así diferentes velocidades alrededor de la frontera que lo conforma. Así pues, como sabemos el campo de velocidades del fluido y las coordenadas de la frontera (Radio=1), procedemos a calcular cómo será el campo en todo el contorno de S:

Ya tenemos la dirección y valores del campo en todos esos puntos. Ahora bien, para hallar los puntos máximos y mínimos, necesitamos llevar a esta función a sus extremos. Para ello, vemos que, usando el módulo de esta, tenemos::

-Adquirirá sus valores máximos cuando:          
   
   

-Y sus valores mínimos cuando:
   

Por tanto, tenemos que: (en cartesianas)

 - Las velocidades máximas se darán en los puntos: [1,0] , [-1,0]
 - Las velocidades mínimas se darán en los puntos: [0,1] , [0,-1]

Los puntos de remanso, que son donde se adquiere velocidad nula, corresponden con los de las velocidades mínimas, como podemos comprobar.

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

7.1 Ecuación de Bernouilli

Sabemos que un fluido se caracteriza, entre otras cosas, por tener carácter compresible o incompresible, determinando así su dinámica. Además de ello, también describirán su comportamiento agentes externos como la temperatura o la presión. En este caso, la ecuación de Bernouilli nos permite relacionar la velocidad de un fluido en función de la presión a la que esté sometido, considerando también una densidad:

Consideramos el fluido incompresible (d=2), el campo de velocidades u y tomando la cte=10. Sustituimos y obtenemos el siguiente campo de presiones:

centro

Así pues, introducimos la ecuación de presiones resultante en Matlab y dibujamos:

% vemos variables
rho=1:0.1:6;
th=0:0.1:2*pi+0.1;
%hacemos retícula
[U,V]=meshgrid(rho,th);
%describo coordenadas
X=U.*cos(V);
Y=U.*sin(V);
%aplicamos funcion presion
f=10-((sin(V).^2)).*((1-(1./(U.^2))))-((((cos(V)).^2)).*(((1./U)+(1./(U.^3))).^2)).*(U.^2);
%dibujo
surf(X,Y,f)
axis([-5,5,-5,5])
axis equal
view(2)

centro

Así pues, viendo el gráfico nos damos cuenta de que las mayores presiones se darán cuando el fluido choque con el obstáculo verticalmente. Y las menores presiones se encontrarán, debido al desplazamiento producido por el choque, en los laterales horizontales del obstáculo. En tanto, lejos del obstáculo se mantendrán más o menos presiones constantes.

7.2 Campo de presiones y velocidad del fluido

Una vez hallado el campo de presiones de nuestro fluido, procedemos a encontrar la relación que, según Bernouilli, lo vincula al campo de velocidades el mismo. De este modo, comparamos los gráficos y vemos:

izquierda

derecha

El fluido, al desplazarse de sur a norte, cuando se encuentra con el fluido choca, viendo su velocidad reducida al mínimo, coincidiendo con la zona donde la presión es máxima. Esto hace que las partículas se desplacen a ambos lados del obstáculo, aumentando de nuevo su velocidad y reduciendo de nuevo la presión a mínimos, ya que no están sometidos a otro choque. Aún así, al llegar al extremo superior del obstáculo, de nuevo las partículas del fluido se chocan y se frenan entre ellas, volviendo a aumentar bastante la presión en esa zona. Por tanto, encontramos una clara relación entre velocidades y presiones.

8 Incompresibilidad del fluido. Ecuación de Navier-Stokes

Uno de los Siete Grandes Problemas del Milenio son las ecuaciones de Navier-Stokes, que tratan de describir el movimiento de un fluido en la atmósfera, las corrientes oceánicas, y sobre cualquier objeto, tales como las aeronaves. De esta forma, se intenta encontrar una relación general entre dichos movimientos y en todos donde estén involucrados fluidos newtonianos, aplicando las leyes elementales de la Mecánica y Termodinámica física. Hasta ahora solo se ha podido demostrar para casos muy particulares. En esta sección vamos a tratar una de las ecuaciones de manera simplificada, donde podamos tomar los parámetros de campo u y presión ,quedando:

Vamos a considerar una viscosidad nula, por lo que lo idealizamos como fluido perfecto. Además, sabemos que es incompresible por tener divergencia nula, lo que hace que tomemos la densidad d=cte. Así, queda: Además, nos dan la siguiente condición: Con todo ello, desarrollamos la ecuación y nos queda (teniendo e cuenta que estamos trabajando en 2D):

centro

Ahora bien, vamos a comprobar ahora que dicha ecuación guarda relación con la ecuación de Bernouilli, ya que para elaborarlas se han tenido en consideración leyes físicas semejantes, y se ha operado con el campo de presiones y de velocidades del fluido, el cual es incompresible y perfecto (sin viscosidad). Para ello, debemos aplicar el gradiente a la misma, y desarrollar en cartesianas partiendo de ahí:

centro

Como podemos comprobar, las soluciones de Bernouilli y Navier-Stokes coinciden, haciendo evidente su conexión tanto matemática como física, aunque solo cuando la viscosidad del fluido es nula.

9 Teorema de Kutta-Joukowski. Paradoja D'Alembert

El Teorema Kutta-Joukowski es un teorema fundamental de la aerodinámica. El teorema relaciona la fuerza que ejerce un fluido con su circulación a lo largo de un cuerpo. La circulación es la integral de línea de la velocidad del fluido, en una curva cerrada que contiene al cuerpo. En las descripciones del teorema Kutta-Joukowski, tomó como sólido un cilindro circular.

En resumen podemos decir que este teorema concluye que la fuerza que ejerce el fluido sobre el obstáculo es proporcional a la circulación. Por otro lado la Paradoja de D´Alembert es una contradicción a la que se llegó tras estudiar matemáticamente el fenómeno de la resistencia producida sobre un cuerpo cuando una corriente de fluido (líquido o gas) 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 es cero, lo cual se contradice con la observación.

En nuestro caso, empezaremos parametrizando el obstáculo:

centro

Calculamos la circulación en el intervalo (a,b)=(0,2π):

centro

centro

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 debería ejercer una fuerza sobre el obstáculo, pero aplicando el Teorema de Kutta-Joukowski, llegamos a la conclusión de que la fuerza es nula también. Encontrándonos con la conocida Paradoja D'Alembert.

10 Curvas de nivel de la presión

Como anteriormente hemos observado, los máximos valores para la presión los encontramos en la zona superior e inferior del óbstaculo, correspondiéndose con las zonas de menor velocidad. Mientras que en 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 obstáculo.

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

centro

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 circular de radios 1 y 5. (Agregar imagen 4) Realizamos la aproximación de la presión mediante Matlab:

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

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

centro