Visualización de campos escales y vectoriales en fluidos.6/19

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

Texto de encabezado	Texto de encabezado
Ejemplo	Ejemplo
Ejemplo	Ejemplo
Ejemplo	Ejemplo
Ejemplo	Ejemplo
Ejemplo	Ejemplo

1 Introducción

Antes de comenzar con el desarrollo del ejercicio, en necesario explicar qué es un fluido incompresible. Un fluido incompresible es todo aquel cuya masa y volumen no sufren variaciones, por lo que permanecen constantes en el tiempo, es decir, se opone a las compresiones ofrecidas bajo cualquier condición o circunstancia. En términos matemáticos, lo que esto significa es que dicho fluido tendrá densidad constante, lo que facilitará las operaciones por homogeneidad de sus propiedades en cualquier posición.

2 Región ocupada por un fluido

Para la representación del fluido utilizamos coordenadas polares con los límites establecidos en el enunciado.

                      x=u*cos(v)
                      y=u*sin(v)

El código para esta representación es:

3 Campo de velocidades de un fluido

Sabiendo que estamos ante un campo escalar en coordenadas cilíndricas, el cálculo del gradiente lo realizaremos a partir de:

                                          Y conociendo

                                         obtenemos:

Representamos la función potencial (de esta ahora vamos a presentar solo las curvas de nivel para poder observar en relación a las velocidades)

Y la función gradiente de ϕ que representa el campo de velocidades del fluido

Observamos en el gráfico que nos salen ortogonales a las curvas de nivel ϕ

4 Gradiente ortogonal a superficies de nivel

La ecuación nos muestra el producto escalar de dos vectores. Como este es igual a cero quiere decir que los vectores son perpendiculares. u es el vector gradiente y n es el vector perpendicular a las curvas de nivel como podemos deducir del enunciado. Lo que interpretamos de esta ecuación una vez conocidos los vectores, es que el gradiente es ortogonal a las superficies de nivel. Para demostrarlo, vamos a comprobar que el ∇f(ρ,θ,z) es ortogonal a cualquier curva contenida en la superficie. Consideramos un punto P de la superficie de nivel S: 'f(ρ,θ,z)= K y c(t) la parametrización de una curva. Se cumple que:

Puesto que hemos dicho que la curva está contenida en la superficie. Si derivamos, la expresión anterior respecto de t usando la regla de la cadena, tenemos que:

Esto significa que ∇f es ortogonal al vector tangente a la curva c, es decir, es ortogonal a la curva. Como hemos trabajado de forma genérica, queda demostrado que el ∇f es ortogonal a todas las curvas de la superficie y, por tanto, es ortogonal a S.

4.1 Caso particular

Nos dicen que ρ es muy grande, y que por lo tanto, 1/ρ≈0. En este caso, u tendrá el siguiente valor:

5 Demostración de rotacional y divergencia nulos

5.1 Rotacional nulo

Comprobamos analíticamente que el rotacional de u es nulo, el rotacional es el determinante de la siguiente matriz:

En nuestro caso al venir la función vectorial del gradiente de una función escalar, el rotacional sería:

De manera que hacer el determinante se anula de la siguiente manera

5.2 Demostración de divergencia nula

Sabiendo que

realizamos el cálculo:

Como nos dicen que esto es igual a cero:

Simplificamos y llegamos a una igualdad y por tanto se verifica lo propuesto.

6 Lineas de corriente del campo

Ahora vamos a dibujar las lineas de corriente v, que van ser tangentes de u. Para ello calculamos el campo que en cada punto es ortogonal a u:v = k x u

También necesitaremos la función de la corriente ψ que es la potencial de v:

Dibujaremos la función potencial inicial, la nueva función potencial que representa la corriente, las velocidades y sus ortogonales:

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

                               vmax=2
                               vmin=0.0248

Ahora dibujamos las líneas de nivel para observar los puntos de remanso, donde la velocidad es nula, con contour(x,y,t)

8 Presión del fluido. Ecuación de bernouilli

Suponiendo que la densidad de flujo es constante d=2 y verificando la ecuación:

Nos pide calcular la presión dando el valor 10 a la constante, dibujarla y en que puntos alcanza al mayor y menor valor.

Tenemos que el campo es:

despejando P en la ecuación de arriba nos queda: P=10-|u|^2, donde:

Por lo tanto.

Ahora lo dibujamos, el código utilizado es:

Puntos donde se alcanza la mayor y menor presión:

Comparándolo con los puntos donde se alcanza mayor y menor velocidad: La velocidad será mínima cuando la presión sea máxima y viceversa

Las gráficas conseguidas son:

9 Ecuación de Navier-Stokes

Haremos el gradiente de: 1/2*d|u| + P = cte, de esa forma demostraremos que cumple la ecuación.

Que el gradiente sea lineal, quiere decir:

Por lo que:

De esta manera queda demostrado que cumple la ecuación de Navier-Stokes.

10 Líneas de corriente

Las líneas de corriente son tangentes a la velocidad. Para poder dibujarlas primero calculamos el vector perpendicular a u:v=k x u

Por ser la divergencia nula como se ha dicho en el apartado anterior, el rotacional será 0.

Calculamos la función potencial φ:

Comprobación:

Ahora lo dibujamos, el código utilizado es:

También tenemos la presión:

Módulo de la velocidad:

¿Cómo cambia la velocidad y presión al rodear al obstáculo?

La presión y la velocidad bajan.

11 Teorema de Kutta-Joukowski.Paradoja D´Alembert

El Teorema Kutta-Joukowski es un teorema perteneciente a la aerodinámica. Se trata de un teorema que relaciona la fuerza que ejerce un fluido con su circulación a lo largo del cuerpo. La circulación es la integral de línea de la velacidad del fluido, en una curva cerrada que contiene al cuerpo. En las descripciones de dicho teorema, 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: (ρ=2,θ=t,z=0,t€(0,2π)

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.

12 Curvas de nivel de presión

Como se ha demostrado en apartados anteriores, al rodear al objeto las presiones se ven disminuidas. Teniendo esto en cuenta, y el intervalo en el que trabajamos, podemos deducir el siguiente código para escribir las curvas de nivel:

%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-(1-1/p.^2).^2*sin(V).^2+(1+1/p.^2).*(cos(V).^2)/(p.^2);
 %Dibujo de las curvas de nivel
 contour(X,Y,p,50,'linewidth',1.5)
 axis([-5,5-5,5])
 axis equal

13 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. 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-(1-1/p.^2).^2*sin(V).^2+(1+1/p.^2).*(cos(V).^2)/(p.^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

Con la aplicación de este código, podemos calcular el valor medio de presión, con lo que obtenemos un valor de