Visualización de campos escalares y vectoriales en fluidos 9-A
Este artículo pretende mostrar la visualización de Campos Escalares y Vectoriales de un fluido incompresible en contacto con un obstáculo circular.
Se ha decidido trabajar en coordenadas cilíndricas (polares) siempre que ésto sea posible.
Se pide realizar el siguiente estudio analítico:
Contenido
- 1 Mallado
- 2 Función Potencial y Campo de Velocidades
- 3 Interpretación de nulidad del producto:
- 4 Valor del gradiente en función de distancia al obstáculo
- 5 Rotacional y Divergencia
- 6 Lineas de Corriente
- 7 Puntos de forntera de la superficie
- 8 Ecuacion de Bernouilli
- 9 Partícula de fluido
- 10 Teorema de Kutta-Joukowski
- 11 Curvas de nivel de la Presion
- 12 Presión media en los puntos del fluido
1 Mallado
El primer paso a realizar es un mallado que represente los puntos interiores de la región ocupada por un fluido. Dicho fluido sera el exterior del círculo unidad (obstáculo). Mallamos un anillo de radio 2 a 6 centrado en el origen. Quedará patente que el fluido ocupa el exterior de un círculo dibujando los ejes en el intervalo [-5,5] x [-5x5]. El siguiente comando de Matlab muestra el mallado generado:
1.1 Comando Matlab de Mallado
- Recogemos aquí los comandos empleados en el programa Matlab para realizar el mallado circular:
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
figure(1)
hold on
plot(cos(t),sin(t),'k','linewidth',2)
mesh(xx,yy,0.*xx)
axis([-5,5,-5,5])
view(2)
hold off
Aunque se nos ha requerido que las operaciones sean realizadas en polares, para poder trabajar en Octave hemos recurrido a la parametrización definida en el comando.
1.2 Gráfica obtenida
- Adjuntamos la gráfica donde queda recogida la región antes citada:
Los puntos interiores del mallado representan puntos del fluido, mientras que el círculo centrado de radio unidad representa el obstáculo.
2 Función Potencial y Campo de Velocidades
2.1 Funcion potencial
Se ha establecido que la velocidad de las partículas del fluido a tratar vengan dadas por el gradiente de la función potencial:
ϕ = 2(ρ + 4/ρ) cos θ
2.2 Campo de velocidades del fluido
Dibujaremos las función gradiente (campo de velocidades) de ϕ, para ello realizaremos los siguientes cálculos, expresando en este apartado el gradiente en covas. Si en algún apartado necesitáramos trabajar en contras, simplemente realizaríamos un cambio de base con la matriz de Gram.
zoom=5;
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
Ur=(2-8./rr.^2).*cos(tt);
Ut=-(2+8./rr.^2).*sin(tt);
Ux=Ur.*cos(tt)-Ut.*sin(tt);
Uy=Ur.*sin(tt)+Ut.*cos(tt);
ff=2.*(rr+4./rr).*cos(tt);
figure(2)
hold on
axis([-zoom,zoom,-zoom,zoom])
plot(cos(t),sin(t),'k','linewidth',2)
contour(xx,yy,ff,50);
quiver(xx,yy,Ux,Uy)
hold off
2.2.1 Gráfica del gradiente del campo
- En la siguiente gráfica se recoge el campo de velocidades de la función potencial(gradiente):
- Como se puede observar en la siguiente imagen ampliada, el gradiente es ortogonal a las curvas de nivel:
3 Interpretación de nulidad del producto: 
Tras el estudio del producto escalar, queda demostrada su nulidad.
La interpretación física de lo anterior quiere decir que a medida que las partículas del fluido avanzan y chocan con el obstáculo fijo, el vector de movimiento de cada una de ellas queda en perpendicular al vector normal del punto del sólido obstáculo con el que chocan produciendo la nulidad en el producto escalar. En consecuencia, las partículas de fluido que entran en contacto con el obstáculo se detienen; adquiriendo momentáneamente velocidad cero. De forma que todas las partículas de fluido van rodeando al obstáculo para seguir avanzando.
4 Valor del gradiente en función de distancia al obstáculo
- A continuación,se pide que nos situemos a una distancia lejana del obstáculo, por lo que podemos suponer despreciable el valor de 1/ro.
En este caso concreto el valor de la función potencial varía, al igual que el valor del gradiente de la función citada.
Podemos comprobar que el valor del gradiente ha aumentado, por lo tanto la velocidad aumentará a medida que nos alejemos del obstáculo.
5 Rotacional y Divergencia
En este apartado, nos disponemos a comprobar analíticamente, la nulidad de la divergencia:
Que la divergencia sea nula,esta asociado a que localmente el fluido mantiene su volumen, siendo este un fluido incompresible.
Comprobamos que el campo también es irrotacional:
Por lo tanto, como el rotacional indica la cantidad de giro neto del fluido al rededor de nuestra superficie, al ser este nulo, no hay giro y entonces la circulación es nula.
6 Lineas de Corriente
Ahora vamos a dibujar las líneas de corriente des campo u, es decir, las líneas que son tangentes a este en cada punto. Para esto calculamos un campo v, que es ortogonal a u.
Una vez que tenemos calculado el campo v, comprobamos analíticamente que es irrotacional (por ser u de divergencia nula).
Además vamos a calcular el potencial escalar ψ, que se conoce como función de corriente de u.
Las líneas ψ=cte, son las líneas de corriente de u que el enunciado nos pide que dibujemos.
El siguiente código de Matlab explica el proceso resolutorio
zoom=5;
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
gg=2.*(rr-4./rr).*sin(tt);
Ur=(2-8./rr.^2).*cos(tt);
Ut=-(2+8./rr.^2).*sin(tt);
Ux=Ur.*cos(tt)-Ut.*sin(tt);
Uy=Ur.*sin(tt)+Ut.*cos(tt);
figure(300)
axis([-zoom,zoom,-zoom,zoom])
plot(cos(t),sin(t),'k','linewidth',2)
hold on
contour(xx,yy,gg,50);
quiver(xx,yy,Ux,Uy)
hold off
En la gráfica, que nos representa las líneas ψ=cte, podemos comprobar, que en efecto son tangentes al campo u, por lo tanto son sus líneas de corriente.
Para verlo mejor hacemos un zoom de la imagen anterior:
Ahora para comprobar que las líneas de corriente son ortogonales a las curvas equipotenciales, realizamos un comando Matlab que nos dibuje en la misma grafica las dos líneas. Así podremos comprobar gráficamente su ortogonalidad.
zoom=5;
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
gg=2.*(rr-4./rr).*sin(tt);
ff=2.*(rr+4./rr).*cos(tt);
Ur=(2-8./rr.^2).*cos(tt);
Ut=-(2+8./rr.^2).*sin(tt);
Ux=Ur.*cos(tt)-Ut.*sin(tt);
Uy=Ur.*sin(tt)+Ut.*cos(tt);
figure(301)
axis([-zoom,zoom,-zoom,zoom])
plot(2.*cos(t),2.*sin(t),'k','linewidth',2)
hold on
contour(xx,yy,gg,50);
contour(xx,yy,ff,50);
hold off
Para apreciar mejor el dibujo hacemos un zoom, en el que se ve mejor que las líneas son ortogonales entre sí.
7 Puntos de forntera de la superficie
Las velocidades del fluido en la frontera S con el obstáculo circular, van variando. Para ver las velocidades máximas y mínimas en esta frontera, entramos en la función del gradiente y damos valor ρ=2.
Viendo la formula se ve claramente que el gradiente (velocidad del fluido), toma valores mínimos para ϴ = 0,Π (en valor absoluto). Los valores máximos los toma para ϴ =Π/2, 3Π/2.
Además podemos comprobarlo mirando el dibujo del gradiente del apartado 2.
8 Ecuacion de Bernouilli
Vamos a suponer que la densidad del fluido es ρ=2, y que se verifica la ecuación de Bernouilli:
Primero sacamos el módulo de u, que al estar en coordenadas polares no se hace directamente.
Además vamos a dar valor 15 a la constante, y a partir de ahí sacamos la presión del fluido:
Con el siguiente comando Matlab, dibujamos dicha presión:
zoom=5;
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
Ur=(2-8./rr.^2).*cos(tt);
Ut=-(2+8./rr.^2).*sin(tt);
Ux=Ur.*cos(tt)-Ut.*sin(tt);
Uy=Ur.*sin(tt)+Ut.*cos(tt);
pp=15-[Ux.^2+Uy.^2]
figure(6)
plot(2.*cos(t),2.*sin(t),'k','linewidth',2)
axis([-zoom,zoom,-zoom,zoom])
hold on
contour(xx,yy,pp,50);
hold off
zoom=5;
dr=0.1;
dt=2*pi/40;
r=2:dr:6;
t=0:dt:2*pi+dt;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
pp=15-[(sin(tt).^2).*((2+8./rr.^2).^2)+((cos(tt).^2).*((2-8./rr.^2).^2))];
figure(63)
surf(xx,yy,pp)
Las dos gráficas nos representan la presión del fluido. Una esta en 3D, que es mas visual a la hora de apreciar las máximos y mínimos de las presiones, y la otra es en 2D y nos representa sus curvas de nivel (que se piden en el apartado 10)
Usando estas graficas, se ve claramente que los puntos de mayor y menor presión se alcanzan en la frontera S.
Los de menor presión son para valores ϴ=Π/2, 3Π/2, que a su vez son los puntos de velocidad máxima.
Los de mayor presión son para valores de ϴ=0, Π, que a su vez son los puntos de velocidad mínima.
Esta conclusión tiene sentido, ya que la propia formula de Bernoulli, nos dice que cuanta mas presión menos velocidad, y viceversa.
9 Partícula de fluido
Si fuéramos una partícula de fluido, seguiríamos una línea de corriente. Fijándonos en las graficas que hemos dibujado anteriormente, nuestra presión y velocidad iría variando según recorremos el obstáculo.
A medida que nos fuéramos acercando al obstáculo, nuestra velocidad iría disminuyendo, hasta llegar a este, en donde se anularía. Este punto a su vez es de presión máxima, por lo que la partícula de fluido tendería a ir a zonas de menor presión, lo que provoca que esta rodee el obstáculo.
Por lo tanto nos trasladaríamos de izquierda a derecha rodeando el obstáculo.
10 Teorema de Kutta-Joukowski
El teorema de Green establece, que la integral de un campo vectorial alrededor de la frontera de la superficie es igual a la del rotacional de este campo alrededor de toda la superficie, por lo tanto la circulación nos la define el rotacional. Este ha sido calculado anteriormente y es nulo, con lo que queda demostrado que la circulación es nula.
Además el teorema de Kutta-Joukowski, establece que la fuerza que el fluido ejerce sobre el obstáculo, es propocional a esta circulación. Si analizamos los resultados que hemos obtenido anteriormente, podemos comprobar que este teorema se cumple, ya que el fluido rodea al obstáculo y no genera ninguna fuerza sobre él.
11 Curvas de nivel de la Presion
Las curvas de nivel de la presión ya han sido representadas en apartados anteriores. De todas formas, las volvemos a dibujar a continuación. (El comando Matlab queda representado en el apartado 8)
12 Presión media en los puntos del fluido
Para calcular la presión media de los puntos del fluido, generamos un comando Matlab, y así calcularlo numéricamente.
Para aproximar la integral de presióne en todo el fluido, utilizamos el método del trapecio. Esta expresión hay que dividirla por el área total del anillo (2<ρ<6)
N1=1000; N2=500; %Number of points
a=2; b=6; c=0; d=2*pi; %Extremes of the interval
h1=(b-a)/N1; h2=(d-c)/N2;
u=a:h1:b; v=c:h2:d; %coordinates of the partition
[uu,vv]=meshgrid(u,v); %coordinates of the rectangle
f=uu.*(15-[(sin(vv).^2).*((2+8./uu.^2).^2)+((cos(vv).^2).*((2-8./uu.^2).^2))]);
w1=ones(N1+1,1); %weights vector
w1(1)=1/2; w1(N1+1)=1/2;
w2=ones(N2+1,1); %weights vector
w2(1)=1/2; w2(N2+1)=1/2;
result=(h1*h2*w2'*f*w1)/(32*pi) %dividimos por el área total del anillo
El resultado que nos devuelve el programa es 10,556, que se corresponde con la presión media de los puntos del fluido.
