Visualización de campos en un fluido incompresible. (Grupo 7-C)

De MateWiki
Saltar a: navegación, buscar


Trabajo realizado por estudiantes
Título Visualización de campos en un fluido incompresible. (Grupo 7-C)
Asignatura Teoría de Campos
Curso 2014-15
Autores Alejandro García Sainz,

Adela González Barbado,

Sergio Reig Vellón,

Diego Solano López,

Ramiro Torres Garófalo

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 .-Visualización de campos escalares y vectoriales en fluidos.

Croquizado de la situación a estudiar.

Vamos a estudiar la visualización de campos escalares y vectoriales en un fluido incompresible alrededor de un obstáculo con forma circular. Trabajaremos en el plano en coordenadas cilíndricas,es decir, polares, por trabajar en 2D.





2 .-Movimiento del fluido.-Representación de los puntos interiores.

Dibujamos en un software científico de visualización (MATLAB en este caso) la malla en la que vamos a trabajar. Hemos realizado el mallado del anillo comprendido entre los radios 2 y 6 con centro el origen. la visualización la hemos realizado con los ejes en el intervalo [-5,5]x[-5,5]. El dibujo se ha modificado añadiendo sobre la malla de MATLAB el obstáculo con un programa de edición de imagen.

Malla.


El código MATLAB utilizado es el siguiente: 


u=linspace(2,6,50);
v=linspace(0,2*pi,50);
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);
yy=uu.*sin(vv);
zz=zeros(50,50);
x=2*cos(v);
y=2*sin(v);
hold on
plot(x,y,'k','linewidth',2)
mesh(xx,yy,zz)
axis([-5,5,-5,5])
hold off


3 .-Función potencial y velocidad.

3.1 Curvas de nivel y Campo de velocidades.

La velocidad de las partículas viene dada por el gradiente de la función potencial φ=2(ρ+4/ρ)cosθ. Hemos representado, haciendo uso de MATLAB, las curvas de nivel de la función potencial φ y el campo de velocidades del fluido.

Los códigos MATLAB utilizados para ello son los que siguen:

Curvas de nivel de la función potencial φ.
  • Representación de las curvas de nivel.
u=linspace(2,6,50);
v=linspace(0,2*pi,50);
fi=inline('2*x+8*x./(x.^2+y.^2)','x','y');
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);
yy=uu.*sin(vv);
f=fi(xx,yy);
contour(xx,yy,f,50);
axis([-5,5,-5,5])


Representación del campo de velocidades.


  • Representación del campo de velocidades.
u=linspace(2,6,10);
v=linspace(0,2*pi,10);
fx=inline('2+(8*(y.^2)-8*(x.^2))./((y.^2+x.^2).^2)','x','y');
fy=inline('-16*(x.*y)./((y.^2+x.^2).^2)','x','y');
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);
yy=uu.*sin(vv);
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)
axis([-5,5,-5,5])


En las siguientes imágenes hemos superpuesto las representaciones de la función potencial y el campo de velocidades (\(\vec{u}\)=Δφ) para comprobar (gráficamente) que éstos son ortogonales. La imagen de la derecha es un zoom sobre la imagen de la izquierda para facilitar la visualización de la perpendicularidad entre las líneas de los dos campos.




Superposición de las representaciones de los dos campos (izquierda) y zoom para facilitar la visualización de la ortogonalidad entre φ y \(\vec{u}\)=Δφ (derecha).


Si \(\vec{n}\) es el vector normal a los puntos del obstáculo, [math]\nabla\vec{u}[/math]•\(\vec{n}\)=0, lo que significa que \(\vec{n}\) es ortogonal a [math]\nabla\vec{u}[/math].

3.2 Rotacional y divergencia.

Hemos definido la función potencial φ=2(ρ+4/ρ)cosθ. También hemos definido [math]\vec{u}=\nabla[/math]φ=(2cosθ-[math]\frac{8}{ρ^2}[/math]cosθ)\(\vec{g^ρ}\)+(-2ρsenθ-[math]\frac{8}{ρ}[/math]senθ)\(\vec{g^θ}\), campo de velocidades del fluido. Si estudiamos la situación alejándonos del obstáculo, ρ es muy grande, y por tanto [math]\frac{1}{ρ}[/math] (y todos sus múltiplos) es despreciable. En este caso:

[math]\vec{u} = (2-\frac{8}{ρ^2})cosθ\vec{g^ρ} - (2ρ + \frac{8}{ρ}) sen θ \vec{g^θ}[/math] → despreciamos los cocientes: [math]\frac{8}{ρ^2}[/math] [math]\frac{8}{ρ}[/math] → 2cosθ[math]\vec{g^ρ} - 2ρsen θ \vec{g^θ}[/math]

Para calcular el rotacional y la divergencia operamos en coordenadas polares:

3.2.1 Rotacional

[math]\vec{\nabla}\times\vec{u}[/math]=[math]\begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ { U }_{ \rho } & { U }_{ \theta } & 0 \end{vmatrix}[/math]= + ([math]\frac { \partial }{ \partial \rho }[/math]Uθ - [math]\frac { \partial }{ \partial \theta }[/math]Uρ)\(\vec{g_z}\)=((-2+[math]\frac{8}{ρ^2}[/math])senθ + (2-[math]\frac{8}{ρ^2}[/math])senθ)\(\vec{g_z}\)=\(\vec{0}\)


El rotacional da la dirección y la velocidad de giro. Un campo con rotacional nulo (o irrotacional), no tiene rotaciones. Extrapolando a nuestro caso, acabamos de comprobar que el rotacional es nulo, por lo que nuestro fluido no genera remolinos.

3.2.2 Divergencia

[math]\vec{\nabla}\vec{u}[/math]=[math]\frac{1}{\sqrt {g}}[/math]([math]\frac { \partial }{ \partial \rho }[/math]([math]\sqrt{g}[/math]Uρ))+[math]\frac{1}{\sqrt {g}}[/math]([math]\frac { \partial }{ \partial \theta }[/math]([math]\sqrt{g}[/math]Uθ))= →[math]\sqrt{g}[/math]=ρ→ =[math]\frac{1}{\rho}[/math]([math]\frac { \partial }{ \partial \rho }[/math](ρ(2-[math]\frac{8}{ρ^2}[/math])cosθ))+[math]\frac{1}{\rho}[/math]([math]\frac { \partial }{ \partial \theta }[/math](-ρ([math]\frac{2}{ρ}[/math]+[math]\frac{8}{ρ^3}[/math])senθ))=[math]\frac{1}{\rho}[/math]((2+[math]\frac{8}{ρ^2}[/math])cosθ+(-2-[math]\frac{8}{ρ^2}[/math])cosθ)=0

La divergencia es la tasa de flujo neto hacia el exterior por unidad de volumen. En términos más sencillos, es lo que se expande el fluido (lo que crece o decrece). Como la divergencia es 0, el fluido es incompresible, es decir, ni masa ni volumen cambian.

3.3 Líneas de corriente.

Las líneas de corriente del campo \(\vec{u}\) son las tangente a la velocidad, entendida como el gradiente de la función potencial φ. Para poder dibujar las líneas de corriente procedemos al cálculo de un vector perpendicular a \(\vec{u}\), el vector [math]\vec{v}=\vec{k}\times\vec{u}[/math]

[math]\vec{v}=\vec{k}\times\vec{u}= \sqrt { g }\begin{vmatrix} \bar { { g }^{ \rho } } & \bar { { g }^{ \theta } } & \bar { { g }^{ z } } \\ 0 & 0 & 1 \\ { u }^{ \rho } & { u }^{ \theta } & { u }^{ z } \end{vmatrix}= ρ (\vec{u^ρ} \vec{g^θ} - \vec{u^θ} \vec{g^ρ}) = \boxed {\vec{v} = 2ρ(1-\frac{4}{ρ^2})cosθ\vec{g^θ} + 2(1 + \frac{4}{ρ^2}) sen θ \vec{g^ρ}= 2(\frac{1}{ρ}-\frac{4}{ρ^3})cosθ \vec{g_θ} + 2(1 + \frac{4}{ρ^2}) sen θ \vec{g_ρ}} [/math]


Por ser la divergencia de \(\vec{u}\) nula, es decir, [math]\nabla·\vec{u}=0[/math] (como se ha calculado en el apartado anterior), el vector \(\vec{v}\) será irrotacional. Calculando [math]\nabla\times\vec{v}=0[/math] podremos comprobarlo.

[math]\nabla\times\vec{v}=\frac{1}{\sqrt { g }}\begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ { v }_{ \rho } & { v }_{ \theta } & { v }_{ z } \end{vmatrix}=\frac{2}{ρ}\begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ (1+\frac{4}{ρ^2}) senθ & ρ(1-\frac{4}{ρ^2}) cosθ & 0 \end{vmatrix}=\frac{1}{ρ}((2+\frac{8}{ρ^2})cosθ- (2+\frac{8}{ρ^2})cosθ)\vec{g_z}= \frac{2}{ρ}((1+\frac{4}{ρ^2}) - (1+\frac{4}{ρ^2}))cosθ\vec{g_z} = 0[/math] → [math]\boxed { \bar { rot } (\bar { v } )=\nabla\times\vec{v}=0 } [/math]


Además \(\vec{v}\) tendrá un potencial escalar ψ que se conoce como la función de corriente de \(\vec{u}\). Dicho potencial lo calcularemos gracias al campo \(\vec{v}\) pues es el gradiente de dicho potencial.

[math]\vec{v}=\vec{grad}ψ= \frac{∂ψ}{∂ρ}\vec{g^ρ} + \frac{∂ψ}{∂θ}\vec{g^θ} + \frac{∂ψ}{∂z}\vec{g^z} = 2(1 + \frac{4}{ρ^2})senθ \vec{g^ρ} + 2ρ(1 - \frac{4}{ρ^2})cosθ \vec{g^θ}[/math]


[math]\left\{\begin{matrix}\ \frac{∂ψ}{∂ρ}= 2(1 + \frac{4}{ρ^2})senθ → ψ=\int 2(1 + \frac{4}{ρ^2})senθ dρ= 2(ρ-\frac{4}{ρ})senθ + h(θ) \\ \frac{∂ψ}{∂θ}=2ρ(1-\frac{4}{ρ^2})cosθ=2(ρ-\frac{4}{ρ})cosθ+h'(θ) → h'(θ) = 0\ltbr /\gt → h(θ) = cte \end{matrix}\right\}[/math] → [math]\boxed {ψ=2(ρ-\frac{4}{ρ})senθ}[/math].


Se observa que h(θ) puede tomar el valor de cualquier constante por lo que podríamos encontrar infinitos potenciales escalares, pues no influye en su gradiente. Fácilmente comprobamos el cálculo de ψ hallando su gradiente [math]\vec{grad}ψ = \frac{∂ψ}{∂ρ}+\frac{∂ψ}{∂θ}= 2(1 + \frac{4}{ρ^2})senθ \vec{g^ρ} + 2(ρ-\frac{4}{ρ})cosθ \vec{g^θ}=2(1 + \frac{4}{ρ^2})senθ \vec{g^ρ} + 2ρ(1-\frac{4}{ρ^2})cosθ \vec{g^θ} [/math] que efectivamente corresponde al campo \(\vec{v}\).

Cabe destacar que al existir el potencial escalar ψ de \(\vec{v}\) podríamos haber anticipado que el rotacional de éste seria nulo y por tanto es que isorrotacional.

Finalmente se dibujan las líneas de corriente del campo de velocidades \(\vec{u}\) sobre las que añadiremos las curvas de nivel del potencial escalar ψ.


Curvas de nivel de las ctes del potencial escalar ψ sobre las líneas de corriente del campo de velocidades \(\vec{u}\).
%Generación de la líneas de corriente de la velocidad
u=linspace(2,6,50);
v=linspace(0,2*pi,50);
fi=inline('2*y-8*y./(x.^2+y.^2)','x','y');

%Creación de la malla y paso a polares
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);
yy=uu.*sin(vv);

%Definición y dibujo del campo de velocidades
f=fi(xx,yy);
hold on
contour(xx,yy,f,100)
axis([-5,5,-5,5])
axis equal

%Generación de las curvas de nivel
u=linspace(2,6,10);
v=linspace(0,2*pi,10);

%Definición y dibujo de función potencial 
fx=inline('2+(8*(y.^2)-8*(x.^2))./((y.^2+x.^2).^2)','x','y');
fy=inline('-16*(x.*y)./((y.^2+x.^2).^2)','x','y');
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);
yy=uu.*sin(vv);
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)
axis([-5,5,-5,5])



En la imagen podemos observar como las curvas de nivel son tangentes al campo de velocidades \(\vec{u}\) en cada uno de los puntos.

3.4 Puntos de frontera.

Gracias al módulo de \(\vec{u}\) podemos establecer una función dependiente de ρ y θ, u(ρ,θ), que nos permitirá calcular aquellos puntos con mayor y menor velocidad en la frontera S fijando el valor de ρ en 2 y 6.

Dado que el módulo no depende de la base escogida podremos hacerlo tanto en coordenadas contravariantes como covariantes, siendo éstas últimas las utilizadas para el cálculo. Además, se hallarán los puntos de remanso, donde la velocidad se anula, que veremos que además son los puntos de velocidad mínima. Tendremos que buscar los máximos y mínimos relativos y posteriormente encontrar entre ellos los absolutos. En los casos en los que la función a maximizar o minimizar ([math]|\vec u|[/math]) sea complicada debido a una raíz cuadrada, se estudiará el valor del interior de la raíz, cuyos puntos donde haya máximos y mínimos coincidirán con los de la velocidad.

[math] |\vec u| (ρ,θ)= 2 {\sqrt { (1-\frac{4}{ρ^2})^2 cos^2θ + (ρ+\frac{4}{ρ})^2 \frac{1}{ρ^2} sen^2θ} } [/math]

♦ ρ=2 → [math]|\vec u|(2,θ)= f(θ)= 4 senθ [/math]
[math] \ Puntos\ de\ remanso\ que\ coinciden\ con\ los\ de\ velocidad\ mínima\ → |\vec u|(2,θ)= 4 senθ=0 \\ θ_1=0 →\ |\vec u|(2,0)=0 \\ θ_2=π\ →\ |\vec u|(2,π)=0 \\ Máximos\ relativos\ → |\vec u|'(2,θ)= f'(θ)= 4 cosθ=0 \\ θ_3=\frac{π}{2} →\ |\vec u|(2,\frac{π}{2})=4 \\ θ_4=\frac{3π}{2} →\ |\vec u|(2,\frac{3π}{2})=4\[/math]


♦ ρ=6 → [math]|\vec u|(6,θ)= 2 {\sqrt { (\frac{64}{81}) cos^2θ + \frac{100}{81}sen^2θ} } [/math] → [math]f(θ)= { (\frac{64}{81}) cos^2θ + \frac{100}{81}sen^2θ} [/math] → [math] \left\{\begin{matrix}\ Mínimos\ relativos\ → f'(θ)= \frac{8}{9}\cosθsenθ=0 \\ θ_1=0 \ → |\vec u|(6,0)=1,78 \\ θ_2=π con\ |\vec u|(6,π)=1,78 \\ Máximos\ relativos\ → f'(θ)= \frac{8}{9}\cosθsenθ= 0 \\ θ_3=\frac{π}{2} \ →|\vec u|(6,\frac{π}{2})=2,22 \\ θ_4=\frac{3π}{2} \ → |\vec u|(6,\frac{3π}{2})=2,22 \\ Puntos\ de\ remanso\ →|\vec u|(6,θ) \ es\ siempre\ distinto\ de\ 0\ luego\ no\ existen .\end{matrix}\right.[/math]


Podemos ver que los puntos de velocidad máxima, mínima y de remanso son: → [math] \left\{\begin{matrix}\ Velocidad\ máxima\ → (ρ=2, θ=\frac{π}{2})\ y\ (ρ=2, θ=\frac{3π}{2})\ con\ |\vec u|=4 \\ Velocidad\ mínima\ y\ de\ remanso → (ρ=2, θ=0)\ y\ (ρ=2, θ=π)\ con\ |\vec u|=0 \end{matrix}\right.[/math]

4 .-Presión.

4.1 Ecuación de Bernouilli.

La ecuación de Bernouilli describe el movimiento del fluido con respecto a ciertas variables de éste que podemos describir como: 

 → p: Presión estática a la que se ve sometido el fluido p debido a las moléculas que lo rodean.
 → ρ: Densidad del fluido.
 → \(\vec{u}\): Velocidad de flujo del fluido.

Esta ecuación se aplica en la dinámica de fluidos, dicho fluido, que podemos considerar gas o líquido, adoptará la forma del recipiente que lo contiene pues sus partículas no están rígidamente unidas, como ocurre con los cuerpos sólidos

Al limitar el nivel de aplicabilidad con la suposiciones que citamos a continuación se puede llegar a la ecuación: [math] \frac{1}{2}ρ (|\vec{u}|)^2 + p = cte [/math]

♦ La velocidad de cada punto no será variante con el tiempo por lo que nuestro fluido se moverá en un régimen estacionario.

♦ Despreciaremos la viscosidad del fluido entendida como una fuerza interna de rozamiento.

♦ El efecto Bernouilli es consecuencia directa que surge a partir de su ecuación: Al aumentar la velocidad del punto disminuirá la presión de éste, es decir, habrá una relación inversa entre la presión y la velocidad del fluido. Dicho efecto nos facilitará el cálculo de los puntos de máxima y mínima presión.


Estableciendo una densidad constante para el fluido con ρ=2 y valor para la constante de la ecuación de 15 podremos calcular cuál es la presión de éste en cada uno de sus puntos.

[math] \frac{1}{2}ρ (|\vec{u}|)^2 + p = cte [/math] → ρ=2 y cte=15 → [math] (|\vec{u}|)^2 + p = 15 [/math] → \(\boxed { p= 15 - 4 [(1-\frac{4}{ρ^2})^2 cos^2θ + (ρ+\frac{4}{ρ})^2 \frac{1}{ρ^2} sen^2θ] } \)

Por último dibujaremos la presión estática P del fluido y cácularemos númericamente el valor máximo y mínimo ésta.

Dibujo en 2D de la superficie de la presión estática P del fluido.
u=linspace(2,6,50);
v=linspace(0,2*pi,50);

%Creación de la malla y paso a cartesianas
[U,V]=meshgrid(u,v);
x=U.*cos(V);
y=U.*sin(V);

%Definimos la presión estática del fluido
P=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2))));

%dibujo de la presión
figure(1)
surf(x,y,P)
axis([-5,5,-5,5])
view(2)

%dibujo en 3D 
figure(2)
surf(x,y,P)
axis([-5,5,-5,5])
view(3)

%Cálculo de la presión máxima y mínima 
Pmax=max(max(P))
Pmin=min(min(P))


Los valores que nos proporciona el programa para la presión máxima de 15 que podemos encontrar en los ángulos θ=0 y θ=π y un valor mínimo de la presión que podemos suponer en θ=π/2 y θ=3π/2 (por los colores de la gráfica, al igual que los máximos) con un valor de 0.

Matlab nos proporciona un valor negativo de -0'98 que físicamente no tendría sentido, debido al uso de un método numérico con un error que podemos reducir al aumentar tamaño del mallado, por ello vemos valores negativos en la gráfica 3D. Analíticamente si hemos obtenido un valor de presión mínima nulo para θ=π/2 y θ=3π/2.


Dibujo en 3D de la superficie de la presión estática P del fluido.

4.2 Curvas de nivel.

En este apartado únicamente se nos pide el dibujo de las curvas de nivel de la presión estática.

Dibujo de las curvas de nivel de la presión estática P del fluido.
%Dibujo de las curvas de nivel de la presión estática
u=linspace(2,6,100);
v=linspace(0,2*pi,100);

%Creación de la malla y paso a polares
[U,V]=meshgrid(u,v);
X=U.*cos(V);
Y=U.*sin(V);

%Dibujo y definición de P
P=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2))));

contour(X,Y,P,100)
axis([-5,5,-5,5])


4.3 Presión media.

Para el cálculo de la presión media hemos utilizado un método numérico en MATLAB basado en el método del trapecio. El resultado de la presión estática media del fluido es de 10'5619, valor que podríamos haber aproximado con la gráfica 3D de la función pues la presión se concentra entre los valores de 10 y 11. 

%Cálculo de la presión media 
i=1000; 
j=1000;                  
u=linspace(2,6,i+1);
v=linspace(0,2*pi,j+1);
h=(6-2)/i;
H=2*pi/j;

%Creación la malla.
[U,V]=meshgrid(u,v);

%Definición de la presión estática
p=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2))));
P=U.*p;

%Integración y cálculo de presión total
A=ones(i+1,1);                 
A(1)=1/2;
A(i+1)=1/2;

B=ones(j+1,1);                 
B(1)=1/2; 
w(j+1)=1/2;
PT=h*H*B'*P*A;

%Cálculo del área y de la presión media
ro=2;
roo=6;
Ar=pi*(roo^2-ro^2);
PM=PT/Ar


4.4 Variación de la velocidad y la presión.

Como se ha comentado en el primer apartado de esta sección existe una relación entre la velocidad de flujo y la presión estática del fluido determinada por el efecto Bernouilli que establece que en los puntos donde la velocidad aumenta disminuirá la presión y viceversa, por lo que podemos afirmar que aquellos donde la presión sea máxima la velocidad será mínima así como donde la primera sea mínima la velocidad será máxima.

La variación de la presión estática y de la velocidad de una partícula del fluido seguirá una relación inversa que podremos comprobar con el dibujo de las gráficas de la presión P y el módulo de la velocidad |\(\vec{u}\)|. Podremos ver en las gráfica como aquellos puntos que alcanzan el valor máximo corresponden al color rojo mientras que aquellos en los que es mínimo se muestran en azul.

Observando las gráficas de la superficie y curvas de nivel del módulo del campo de velocidades \(\vec{u}\) y de la presión, junto con las líneas de corriente (ψ), podemos ver cómo varían estos dos valores a medida que nos desplazamos por distintas líneas de corriente.

♦ Observando las gráficas de la superficie y curvas de nivel del módulo del campo de velocidades \(\vec{u}\), podemos ver como los puntos máximo se concentran puntos superiores e inferiores (θ=0 y θ=π y alrededores) en el obstáculo (ρ=2) y algunos ligeramente menores en los lados (θ=π/2 y θ=3π/2 y alrededores) en la frontera exterior de ρ=6. El resto de la gráfica la observamos en colores azulados que nos indican los puntos de menor velocidad.

♦ Sin embargo en el dibujo de las curvas de nivel de la presión, así como en su superficie, se observa como ocurre exactamente lo contrario, ésta tiende a colores rojizos en la mayor parte de la función salvo en aquellos puntos donde la velocidad era rojo que podemos ver en tonos azules.


Cuando nos aproximamos al obstáculo desde la línea de corriente central reducimos nuestra velocidad hasta llegar a un punto de remanso (θ=π y ρ=2). En el resto de líneas de corriente, rodeamos el obstáculo aumentando la velocidades a medida que rodeamos el obstáculo (disminuyendo la presión) hasta llegar a un punto (θ=π/2 o θ=3π/2), desde el que seguiremos rodeando el obstáculo dado que no hay velocidad normal, y desde el que empezará a disminuir la velocidad y a aumentar la presión. Al terminar la curva nuestra partícula saldrá con la misma velocidad y presión con la que entraba. 


ro=linspace(2,6,50);
th=linspace(0,2*pi,50);

%Creaación de la malla y paso a cartesianas
[U,V]=meshgrid(ro,th);
x=U.*cos(V);
y=U.*sin(V);

%Definición de las líneas de corriente
UU=2.*sqrt(((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2)))));

%Definimos la presión estática del fluido
P=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2))));

%dibujo de ambas superficies 
figure(1)
surf(x,y,UU)
axis([-5,5,-5,5])
view(2)

figure(2)
surf(x,y,P)
axis([-5,5,-5,5])
view(2)

figure(3)
contour(x,y,UU,100)
axis([-5,5,-5,5])
view(2)

figure(4)
contour(x,y,P,100)
axis([-5,5,-5,5])
view(2)


Superficie del módulo del campo de velocidades \(\vec{u}\).
Superficie de la presión estática P.


Curvas de nivel del módulo del campo de velocidades \(\vec{u}\).
Curvas de nivel de la presión estática P.




5 .- Teorema de Kutta-Joukowski y paradoja de D'Alembert.

Debemos comprobar que la circulación del campo \(\vec{u}\) a lo largo de la circunferencia de radio 2 centrada en el origen.

Antes de proceder al cálculo directo de la integral analizaremos nuestro caso detenidamente pues podemos anticiparnos al resultado. Nuestro campo \(\vec{u}\) posee un potencial escalar φ, pues su rotacional es cero lo cual nos lleva a dos posibles deducciones del valor de la circulación en la curva.

→ Por un lado el teorema de Stokes establece que la integral de un campo vectorial alrededor de una frontera de la superficie es igual a la del rotacional de éste mismo campo alrededor de toda la superficie. Al ser el rotacional de campo \(\vec{u}\) nulo podemos afirma que la circulación a lo largo del obstáculo entendido como una curva cerrada será nula. 

 [math] ∫_{əS} \vec {u} \vec {dr} = ∫_S (\nabla\times\vec{u}) \vec {dS}=0 [/math]

→ Por otro lado al poseer el campo un potencial escalar, la circulación no dependerá del camino recorrido pues únicamente influye en ella el origen y el extremo de éste pues todo campo que tenga asociado un potencial escalar será conservativo. 

[math] ∫_Ύ \vec {u} \vec {dr} = ∫_0^{2π} \vec u \vec {dr}= φ(2π)- φ(0)=0 [/math]

Además vemos que la curva sobre la que estamos trabajando tiene como origen y como extremo el mismo punto, es decir, se trata de una curva cerrada y por tanto la circulación será nula. 

[math] ∫_Ύ \vec {u} \vec {dr} = ∫_0^{2π} \vec u \vec {dr}=0 [/math]

El teorema de Kutta-Joukowski, que establece que la fuerza que ejerce el fluido sobre el obstáculo no es sólo proporcional a la circulación si no también a la densidad y la velocidad de fluido. La ecuación viene dado por: 

→l: La fuerza de sustentación.
→ρ: Densidad del fluido.
→v: Velocidad del fluido.
→Γ: La circulación.

\(\boxed { l=ρvΓ } \)


Al ser nula la circulación, el obstáculo no se ve sometido a ninguna fuerza por el fluido lo cual nos lleva a mencionar la siguiente paradoja pues el resultado obtenido se aleja por completo de aquel hubiéramos podido intuir con la descripción del movimiento del fluido, la paradoja de D'Alembert.

Si el resultado de las fuerzas ejercidas sobre el obstáculo no hubiesen sido nulas éstas serían ortogonales a las líneas de corriente.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_en_un_fluido_incompresible._(Grupo_7-C)&oldid=37488»