Visualización de campos escalares y vectoriales en un fluido. Grupo A13

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en un fluido
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Alberto Rodríguez Soto, Ignacio Pérez Abellán
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Nuestro trabajo consiste en el estudio de un fluido incompresible, esto es, que su densidad es constante en el tiempo y que se opone a cualquier fuerza externa de compresión. Para realizar el trabajo hemos sometido al fluido a diversos campos escalares y vectoriales según lo aprendido en la asignatura.

2 Mallado de la 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 1 y radio mayor 5.

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

%rho=u; phi=v;
%Creamos el mallado.
rho=linspace(1,5,50);
th=linspace(0,2*pi,50);
%generamos la reticula
[U,V]=meshgrid(rho,th);
%definimos la parametrización
X=U.*cos(V);
Y=U.*sin(V);
%Dibujamos la malla y definimos los ejes
mesh(X,Y,0*X)
axis([-5,5,-5,5]);
axis equal

Mallado

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: Se aprecia como la función potencial crece de abajo a arriba según el movimiento del fluido

%Dibujamos la malla 
 rho=linspace(1,5,25);
 th=linspace(0,2*pi,25);
 [U,V]=meshgrid(rho,th);
 %dibujamos la función potencial
 %primero la superficie
 X=U.*cos(V);
 Y=U.*sin(V);
 %Despues definimos la función
 phi=inline('(rho+(1./rho)).*sin(th)','rho','th');
 %aplicamos la función
 Z=phi(U,V);
 %Dibujamos
 p=surf(X,Y,Z);
 axis equal
 view(2)
 hold on
 %por último dibujamos el obstáculo
 plot(cos(th),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

Su representación viene dada por:

%Dibujamos la malla 
rho=linspace(1,5,30);
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th);
%dibujamos la función potencial
%primero la superficie
X=U.*cos(V);
Y=U.*sin(V);
%Despues definimos la función
phi=inline('(rho+(1./rho)).*sin(th)','rho','th');
%aplicamos la función
Z=phi(U,V);
%Dibujamos las curvas de nivel
contour(X,Y,Z,40)
%axis([-5,5,-5,5]);
axis equal
view(2)
hold on
%Dibujamos el campo de velocidades del fluido:
ux=sin(V).*cos(V)-(sin(V).*cos(V))./(U.^2)-(sin(V).*cos(V)+(sin(V).*cos(V))./(U.^2));
uy=(sin(V).*sin(V)-(sin(V).*sin(V))./(U.^2))+cos(V).*cos(V)+(cos(V).*cos(V))./(U.^2);
%dibujamos el campo de velocidades
quiver(X,Y,ux,uy);
axis equal
view(2)
hold on
%por último dibujamos el obstáculo
plot(cos(th),sin(th),'k','linewidth',2)

Representación campo de velocidades

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 se aprecia la ortogonalidad a las curvas de nivel de la función potencial φ=(ρ+1/ρ)·sen(θ) 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 abajo a arriba 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.

4 Vector normal a los puntos del obstáculo

Si consideramos nuestro obstáculo como la superficie circular de radio 1 y centro origen, cuyo vector normal a su superficie es y dado nuestro campo de velocidades Observamos claramente cómo al hacer el producto escalar entre ambos vectores, el resultado es nulo. La interpretación física que podemos dar a esto, es que cuando las partículas del fluido chocan contra el obstáculo, éstas cambian de dirección de modo que siempre van a rodear el obstáculo y nunca pasarlo por encima

5 Valor de la velocidad al alejarnos del obstáculo

Debido a que estamos muy lejos del obstáculo, podemos despreciar el valor de 1/ρ, nos quedará el campo escalar

El gradiente del campo escalar cuando nos encontramos muy lejos del centro será

6 Divergencia y rotacional

Con el estudio de la divergencia y del rotacional se demuestra la incompresibilidad del fluido. Se trata 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: Resolviendolo, podemos comprobar que es irrotacional.

Así mismo, la divergencia se define como: y tambien comprobamos que es 0

7 Lineas de corriente del campo de velocidades

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 del fluido. 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. Multiplicando vectorialmente el vector k por la velocidad, obtenemos:

7.1 Comprobación gráfica

• Campo v en coordenadas cartesianas
• Código MATLAB

%Reticula
rho=linspace(1,5,30);
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th);
%Cambio de coordenadas a cartesianas
X=U.*cos(V);
Y=U.*sin(V);
%Componentes del campo de velocidades en cartesianas
ux=sin(V).*cos(V)-(sin(V).*cos(V))./(U.^2)-(sin(V).*cos(V)+(sin(V).*cos(V))./(U.^2));
uy=(sin(V).*sin(V)-(sin(V).*sin(V))./(U.^2))+cos(V).*cos(V)+(cos(V).*cos(V))./(U.^2); 
%Componentes del campo v en cartesianas
vx=(-cos(V).*cos(V))./U-U.*sin(V).*sin(V);
vy=(-cos(V).*sin(V))./U+U.*cos(V).*sin(V);
%Dibujo campo de velocidades
quiver(X,Y,ux,uy)
axis([-5,5,-5,5])
axis equal
view(2)
hold on
%Dibujo campo v
quiver(X,Y,vx,vy)
axis([-5,5,-5,5])
axis equal
view(2)

7.2 Cálculo de la función de corriente

La función de corriente, será la función potencial del campo v obtenido previamente.

La representaremos en MATlab mediante el siguiente código:

%Reticula
%Dibujamos la malla 
rho=linspace(1,5,25);
th=linspace(0,2*pi,25);
[U,V]=meshgrid(rho,th);

%Cambio de coordenadas a cartesianas
X=U.*cos(V);
Y=U.*sin(V);

%Potencial escalar. Lineas de corriente
chi=log(U).*sin(V)-U.*U.*cos(V);

%Componentes del campo de velocidades en cartesianas
ux=sin(V).*cos(V)-(sin(V).*cos(V))./(U.^2)-(sin(V).*cos(V)+(sin(V).*cos(V))./(U.^2));
uy=(sin(V).*sin(V)-(sin(V).*sin(V))./(U.^2))+cos(V).*cos(V)+(cos(V).*cos(V))./(U.^2); 
 
%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

En la figura observamos claramente que las de corriente son tangentes a la velocidad y ortogonales al campo equipotencial como bien dice su definición.

8 Velocidad máxima y mínima en la frontera del obstáculo

En nuestro caso, vamos a calcular cuál es la velocidad máxima y mínima en el entorno del obstáculo. Para ello, dado que estamos en la frontera, suponemos que ρ=1, y el campo u queda

|u|=2·cos(θ)

Dado que el valor del campo vectorial nos queda en función del ángulo θ, podemos afirmar lo siguiente:

• La velocidad será máxima para los valores de θ que sean iguales a 0 y -Π, alcanzando el valor de 2 unidades.
• La velocidad será mínima para los valores de θ que sean iguales a Π/2 y -Π/2, alcanzando el valor de 0 unidades. Por tanto, estos puntos además de ser mínimos serán puntos de remanso.

9 Presión del fluido para una densidad constante

Como hemos dicho previamente, estamos tratando con fluidos incompresibles, y por tanto, van a tener una densidad constante, que en este caso es 2. Sustituyendo en la fórmula de Bernoulli e igualándolo a la constante de 10, nos queda El código empleado en MATlab para su visualización es:

%Dibujamos la malla 
rho=linspace(1,5,25);
th=linspace(0,2*pi,25);

%hacemos retícula
[U,V]=meshgrid(rho,th);

%describo coordenadas
X=U.*cos(V);
Y=U.*sin(V);
 
%aplicamos funcion presion
p=10-(sin(V)-(1./(U.^2)).*sin(V)).^2-(U.^2).*(((1./U)+(1./U.^3)).^2).*cos(V).*cos(V);

%dibujo
surf(X,Y,p)
axis([-5,5,-5,5])
axis equal
view(2)

10 Ecuación de Navier-Stokes

Las ecuaciones de Navier-Stokes reciben su nombre de Claude-Louis Navier y George Gabriel Stokes. Se trata de un conjunto de ecuaciones en derivadas parciales no lineales que describen el movimiento de un fluido. Estas ecuaciones gobiernan la atmósfera terrestre, las corrientes oceánicas y el flujo alrededor de vehículos o proyectiles y, en general, cualquier fenómeno en el que se involucren fluidos newtonianos.

Estas ecuaciones se obtienen aplicando los principios de conservación de la mecánica y la termodinámica a un volumen fluido. Haciendo esto se obtiene la llamada formulación integral de las ecuaciones. Para llegar a su formulación diferencial se manipulan aplicando ciertas consideraciones, principalmente aquella en la que los esfuerzos tangenciales guardan una relación lineal con el gradiente de velocidad (ley de viscosidad de Newton), obteniendo de esta manera la formulación diferencial que generalmente es más útil para la resolución de los problemas que se plantean en la mecánica de fluidos.

Desarrollando con los cálculos pertinentes obtenemos que: 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.

11 Interpretación según la velocidad y la presión

Analizando los gráficos que hemos obtenido en apartados previos, vamos a calcular ahora las velocidades que tendrían las partículas del fluido en la frontera del obstáculo. Las velocidades menores, las obtenemos en los extremos inferior y superior, y debido a esto, es aquí donde se acumulará un gran número de partículas que es lo que nos supone que la presión sea tan elevada. Por contra, en los extremos izquierdo y derecho, las partículas tendrán una velocidad máxima y por lo tanto la presión aquí será mínima.

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

La paradoja de D'Alembert, es una contradicción a la que llegó D'Alembert luego de 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 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 cómo la circulación a lo largo de la circunferencia de radio 2 se anula.

Dada la curva: γ=[(ρ,θ)/ρ=2; θ=t; t€(0,2Π) y parametrizando el campo vectorial y desarrollando, llegamos a la integral:

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.

13 Curvas de presión

En la figura observamos que la presión es mucho más alta en los extremos superior e inferior, siendo esto coherente con lo desarrollado en apartados anteriores.

%Dibujamos la malla 
rho=linspace(1,5,25);
th=linspace(0,2*pi,25);

%hacemos retícula
[U,V]=meshgrid(rho,th);

%describo coordenadas
X=U.*cos(V);
Y=U.*sin(V);
 
%aplicamos funcion presion
p=10-(sin(V)-(1./(U.^2)).*sin(V)).^2-(U.^2).*(((1./U)+(1./U.^3)).^2).*cos(V).*cos(V);

%dibujo
contour(X,Y,p,50,'linewidth',1.5)
axis([-5,5,-5,5])
axis equal

Curvas de nivel de la presión del fluido

14 Presión media del fluido

Vamos a calcular la presión media del fluido. Para ello se hará un aproximación de la integral de la presión 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 5 y 1. El área es igual a π·(25-1)=24π

Con esto, ya podemos calcular la presión media del fluido con el siguiente código

%Dibujamos la malla 
h=0.01;
rho=1:h:5;
th=0:h:2*pi;
[U,V]=meshgrid(rho,th);
%Presión del fluido
p=10-(sin(V)-(1./(U.^2)).*sin(V)).^2+(U.^2).*(((1./U)+(1./U.^3)).^2).*cos(V).*cos(V);
%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/area

El resultado de la presión media ha resultado ser P=10.3065