Visualización con Matlab de campos escalares y vectoriales en fluidos (Grupo 4C)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización con Matlab de campos escalares y vectoriales en fluidos Grupo 4C |
| Asignatura | Teoría de Campos |
| Curso | 2015-16 |
| Autores |
Juan Carlos Durán |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Vamos a considerar el flujo de un fluido incompresible alrededor de un obstáculo con forma circular con forma de anillo comprendido entre las circunferencias centradas en el origen de radios 2 y 6. Trabajaremos en el plano con coordenadas polares [math](\rho,\theta)[/math]. Para ilustrar que el fluido ocupa el exterior de un círculo consideraremos [math](\rho,\theta) \in [-5,5]*[-5,5][/math]
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial [math]\varphi=2(\rho+\frac{4}{\rho})cos(\theta)[/math]
Contenido
- 1 º Mallado del sólido
- 2 º Función potencial [math]\varphi[/math] y velocidad de las partículas del fluido
- 3 º Ortogonalidad [math]\nabla\vec u[/math] y [math]\vec n[/math]
- 4 º Rotacional y Divergencia
- 5 º Líneas de corriente
- 6 º Puntos frontera
- 7 º Presión del fluido
- 8 º Ecuación de Navier-Stokes
- 9 º Trayectoria de la línea de corriente
- 10 º Paradoja D'Alembert
- 11 º Curvas de nivel de la presión
- 12 º Presión media
1 º Mallado del sólido
Nuestro mallado se ha centrado en la representación de los puntos interiores del anillo circular por el que circulará el fluido. Se ha elegido un paso de muestreo similar a otros trabajos, de [math]h=\frac{1}{10}[/math] para las variables [math]\rho[/math] y [math]\theta[/math].
Código Matlab:
h=0.1; % Paso de muestreo.
u=2:h:6; % Intervalo [2,6].
v=0:h:2*pi+h; % Intervalo [0,2*pi].
[uu,vv]=meshgrid(u,v); % Matrices.
xx=uu.*cos(vv); % Parametrización X.
yy=uu.*sin(vv); % Parametrización Y.
zz=zeros(size(xx)); % Parametrización Z.
subplot(1,1,1); % Muestra imágenes.
mesh(xx,yy,zz) % Mallado.
axis([-5,5,-5,5]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
2 º Función potencial [math]\varphi[/math] y velocidad de las partículas del fluido
2.1 º Función potencial [math]\varphi(\rho,\theta,{z})[/math]
[math]\varphi=2(\rho+\frac{4}{\rho})cos(\theta)[/math]
Código Matlab de la función potencial [math]\varphi[/math] mediante 30 líneas de nivel:
f=2.*(uu+4./uu).*cos(vv); % Campo Escalar de la Temperatura.
subplot(1,1,1); % Muestra una imagen.
contour(xx,yy,f,30); colorbar % Define 30 líneas de nivel.
2.2 º Campo de velocidades del fluido [math]\nabla\varphi[/math]
La velocidad de las partículas viene dada por el gradiente de la función [math]\varphi[/math] potencial. Es decir, [math]\vec u=\nabla\varphi[/math].:
[math]\vec{u}=\nabla\varphi=\frac{\partial\varphi}{\partial {x}^{i}}{\vec {g}}^{i}=\frac{\partial\varphi}{\partial\rho}{\vec{g}}^{\rho}+\frac{\partial\varphi}{\partial\rho}{\vec{g}}^{\theta}+\frac{\partial\varphi}{\partial{z}}{\vec{g}}^{z}=\frac{\partial\varphi}{\partial\rho}{\vec{g}}_{\rho}+\frac{1}{{\rho}^{2}}\frac{\partial\varphi}{\partial\theta}{\vec{g}}_{\theta}+\frac{\partial\varphi}{\partial{z}}{\vec{g}}_{z}[/math]:
[math]\vec{g}_{\rho}=\vec{e}_{\rho}\qquad\vec{g}_{\theta}=\rho\vec{e}_{\theta}[/math]
En este caso::
[math]\nabla\varphi=2cos(\theta)[1-\frac{4}{\rho^2}]{\vec{g}}^{\rho}-2sen(\theta)[\rho+\frac{4}{\rho}]{\vec{g}}^{\theta}=2cos(\theta)[1-\frac{4}{{\rho}^2}]{\vec{g}}_{\rho}-2sen(\theta)[\frac{1}{\rho}+\frac{4}{{\rho}^3}]{\vec{g}}_{\theta}=2cos(\theta)[1-\frac{4}{{\rho}^2}]{\vec{e}}_{\rho}-2sen(\theta)[{1}+\frac{4}{{\rho}^2}]{\vec{e}}_{\theta}[/math]:
[math]\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})[/math]:
Código Matlab:
syms x y % Variables definidas para la función "diff".
F=2.*(sqrt(x.^2+y.^2)+4./sqrt(x.^2+y.^2)).*(x./(sqrt(x.^2+y.^2))); % Campo Potencial en coordenadas cartesianas.
FX=diff(F,x) % Derivada del campo respecto de x (Simbólico).
FY=diff(F,y) % Derivada del campo respecto de y (Simbólico).Mediante este código obtenemos las derivadas parciales que componen la función Gradiente del campo potencial [math]\varphi[/math]::
[math]FX=\frac{\partial\varphi(x,y)}{\partial x}=2+\frac{8y^2-8x^2}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{16xy}{(x^2+y^2)^2}[/math]
Código Matlab:
fx=inline('2+(8*(y.^2)-8*(x.^2))./((y.^2+x.^2).^2)','x','y'); % Función derivada del campo respecto de x.
fy=inline('-16*(x.*y)./((y.^2+x.^2).^2)','x','y'); % Función derivada del campo respecto de y.
d1=fx(xx,yy); % Matrices de la función derivada respecto de x.
d2=fy(xx,yy); % Matrices de la función derivada respecto de y.
hold on % Inicio superposición de gráficos.
contour(xx,yy,f,30); colorbar % Define 30 líneas de nivel.
quiver(xx,yy,d1,d2) % Representación de los vectores gradientes del campo escalar.
axis([-5,5,-5,5]) % Selecciona la región a dibujar.
hold off % Fin superposición de gráficos.
Se comprueba en la imagen que la función potencial [math]\varphi[/math] y el campo de velocidades [math]\nabla\varphi[/math] son ortogonales.:
3 º Ortogonalidad [math]\nabla\vec u[/math] y [math]\vec n[/math]
Se denomina a [math]\vec n[/math] como vector normal a los puntos del obstáculo, es decir [math]\vec n= \vec k[/math]: [math]\nabla\varphi=\begin{vmatrix} \frac{\partial{u}_{1}}{\partial{\rho}} & \frac{\partial{u}_{1}}{\partial{\theta}} & \frac{\partial{u}_{1}}{\partial{z}} \\ \frac{\partial{u}_{2}}{\partial{\rho}} & \frac{\partial{u}_{2}}{\partial{\theta}} & \frac{\partial{u}_{2}}{\partial{z}} \\ \frac{\partial{u}_{3}}{\partial{\rho}} & \frac{\partial{u}_{3}}{\partial{\theta}} & \frac{\partial{u}_{3}}{\partial{z}} \end{vmatrix}=\begin{vmatrix} \frac{16cos({\theta})}{\rho^3} & -2sen(\theta)(1-\frac{4}{\rho^2}) & {0}\\ 2sen(\theta)(\frac{1}{\rho^2}+\frac{12}{\rho^4}) & -2cos(\theta)(\frac{1}{\rho}+\frac{4}{\rho^3})&{0}\\ {0}&{0}&{0}\end{vmatrix}[/math]:
[math]\nabla\varphi\cdot\vec n=\begin{vmatrix} \frac{\partial{u}_{1}}{\partial{\rho}} & \frac{\partial{u}_{1}}{\partial{\theta}} & \frac{\partial{u}_{1}}{\partial{z}} \\ \frac{\partial{u}_{2}}{\partial{\rho}} & \frac{\partial{u}_{2}}{\partial{\theta}} & \frac{\partial{u}_{2}}{\partial{z}} \\ \frac{\partial{u}_{3}}{\partial{\rho}} & \frac{\partial{u}_{3}}{\partial{\theta}} & \frac{\partial{u}_{3}}{\partial{z}} \end{vmatrix}=\begin{vmatrix} \frac{16cos({\theta})}{\rho^3} & -2sen(\theta)(1-\frac{4}{\rho^2}) & {0}\\ 2sen(\theta)(\frac{1}{\rho^2}+\frac{12}{\rho^4}) & -2cos(\theta)(\frac{1}{\rho}+\frac{4}{\rho^3})&{0}\\ {0}&{0}&{0}\end{vmatrix}\begin{vmatrix}{0}\\{0}\\{1}\end{vmatrix}=\begin{vmatrix}{0}\\{0}\\{0}\end{vmatrix}[/math]
Como interpretación podríamos decir que el vector velocidad [math]\vec v[/math], además de ser perpendicular a las líneas de nivel de la función potencial, está únicamente definido en el plano 2D y no en un espacio 3D. Como resultado evidente, el producto escalar del gradiente de la velocidad y el vector normal siempre será cero.
Supongamos ahora que nos encontramos lejos del obstáculo siendo así [math]\rho[/math] muy grande, por tanto suponemos que [math]\frac{1}{\rho}≃{0}[/math] en esos puntos::
[math]\vec{u}=2cos(\theta){\vec{e}_{\rho}}-2sen(\theta){\vec{e}_{\theta}}[/math]
Por tanto::
[math]|\vec{u}| =2\sqrt{cos(\theta)^2+sen(\theta)^2}=2[/math]
4 º Rotacional y Divergencia
[math]\nabla\times\vec u=\begin{vmatrix} \vec {g}_{\rho}&\vec {g}_{\theta}&\vec {g}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ 2cos(\theta)[1-\frac{4}{\rho^2}] & -2sen(\theta)[1-\frac{4}{\rho^2}] & {0} \end{vmatrix}=\frac{\partial}{\partial{\rho}}[-2sen(\theta)(1-\frac{4}{\rho^2})]\vec {g}_{z}-\frac{\partial}{\partial{\theta}}[-2cos(\theta)(1-\frac{4}{\rho^2})]\vec {g}_{z}=[-2sen(\theta)(1-\frac{4}{\rho^2})]\vec {g}_{z}+[2sen(\theta)(1-\frac{4}{\rho^2})]\vec {g}_{z}=\vec {0}[/math]
[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}(2{\rho}cos(\theta)(\rho -\frac{4}{\rho^2}))+\frac{\partial}{\partial{\theta}}(-2sen(\theta)({1}+\frac{4}{\rho^2})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[2cos(\theta)(1+\frac{4}{\rho^2})-2cos(\theta)(1+\frac{4}{\rho^2})]=0 [/math]
El rotacional obtenido es nulo, quedando demostrado que las partículas del fluido no giran. La divergencia obtenida también es nula, indicando que localmente el fluido no varía su volumen, es decir, ni se expande ni se contrae. Teniendo en cuenta este resultado, podemos confirmar que el fluido es incompresible.
5 º Líneas de corriente
Vamos a dibujar ahora las líneas de corriente del campo [math]\vec u[/math], es decir las líneas que son tangentes a [math]\vec u[/math] en cada punto. Para ello calculamos el campo [math]\vec v[/math] que en cada punto es ortogonal a [math]\vec u[/math] (tomar [math]\vec v=\vec k\times\vec u[/math]).:
[math]\vec v=\frac{1}{\rho}=\begin{vmatrix} \vec {g}_{\rho}&\vec {g}_{\theta}&\vec {g}_{z} \\ {0} & {0} & {1} \\ 2cos(\theta)[1-\frac{4}{\rho^2}] & -2sen(\theta)[\rho+\frac{4}{\rho}] & {0} \end{vmatrix}=2sen(\theta)[1+\frac{4}{\rho^2}]\vec {g}_{\rho}+2cos(\theta)[\frac{1}{\rho}-\frac{4}{\rho^3}]\vec {g}_{\theta}[/math]: [math]\vec v=\vec {grad}\psi[/math]:
[math]\vec v=\vec {grad}\psi=\frac{\partial\psi}{\partial\rho}\frac{{\vec {g}}_{\rho}}{|{\vec {g}}_{\rho}|^2}+\frac{\partial\psi}{\partial\theta}\frac{{\vec {g}}_{\theta}}{|{\vec {g}}_{\theta}|^2}+\frac{\partial\psi}{\partial{z}}\frac{{\vec {g}}_{z}}{|{\vec {g}}_{z}|^2}=\frac{\partial\psi}{\partial\rho}\frac{{\vec {g}}_{\rho}}{1}+\frac{\partial\psi}{\partial\theta}\frac{{\vec {g}}_{\theta}}{\rho^2}+\frac{\partial\psi}{\partial{z}}\frac{{\vec {g}}_{z}}{1}[/math]:
[math]\psi=\int 2sen(\theta)[1+\frac{4}{\rho^2}] \, d\rho \,\!=2sen(\theta)[\rho+\frac{4}{\rho}]+f(\theta) [/math]: [math]\frac{\partial\psi}{\partial\theta}=2cos(\theta)[\rho-\frac{4}{\rho}+f'(\theta)\qquad f'(\theta)=0\qquad f(\theta)=cte[/math]: [math]\psi=2sen(\theta)[\rho-\frac{4}{\rho}][/math]:
[math]\nabla\times\vec v=\begin{vmatrix} \vec {g}_{\rho}&\vec {g}_{\theta}&\vec {g}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ 2sen(\theta)[1+\frac{4}{\rho^2}] & 2cos(\theta)[\rho-\frac{4}{\rho}] & {0} \end{vmatrix}=\frac{\partial}{\partial{\rho}}[2cos(\theta)[\rho-\frac{4}{\rho}]]\vec {g}_{z}-\frac{\partial}{\partial{\theta}}[2sen(\theta)[1+\frac{4}{\rho^2}]]\vec {g}_{z}=[2cos(\theta)(1+\frac{4}{\rho^2})]\vec {g}_{z}-[2cos(\theta)(1+\frac{4}{\rho^2})]\vec {g}_{z}=\vec {0}[/math]
Como [math]\vec u[/math] es de divergencia nula, queda comprobado que [math]\vec v[/math] es irrotacional.
Código Matlab:
fi=inline('2*y-8*y./(x.^2+y.^2)','x','y'); % Función línea de corriente
f=fi(xx,yy); % Línea de corriente
hold on % Inicio superposición de gráficos.
contour(xx,yy,f,100) % Define 100 líneas de corriente.
fx=inline('2+(8*(y.^2)-8*(x.^2))./((y.^2+x.^2).^2)','x','y'); % Función componente horizontal del campo de velocidades.
fy=inline('-16*(x.*y)./((y.^2+x.^2).^2)','x','y'); % Función componente vertical del campo de velocidades.
u1=fx(xx,yy); % Componente horizontal del campo de velocidades.
u2=fy(xx,yy); % Componente horizontal del campo de velocidades.
quiver(xx,yy,u1,u2) % Representación de los vectores del campo de velocidad.
axis([-5,5,-5,5]) % Selecciona la región a dibujar.
hold off % Fin superposición de gráficos.
6 º Puntos frontera
Definido el campo de velocidades anterior, se calcula su módulo en función de [math]\theta[/math].
[math]\vec u=2sen(\theta)[1+\frac{4}{\rho^2}]\vec {g}_{z}-2cos(\theta)[1+\frac{4}{\rho^2}]\vec {g}_{z}=2cos(\theta)[1-\frac{4}{{\rho}^2}]{\vec{e}}_{\rho}-2sen(\theta)[{1}+\frac{4}{{\rho}^2}]{\vec{e}}_{\theta}[/math]
Para [math]\rho=2\qquad|\vec u|=\sqrt{4cos(\theta)^2[1-\frac{4}{{2}^2}]^2+4sen(\theta)^2[{1}+\frac{4}{{2}^2}]^2}=\sqrt{16sen(\theta)^2}=4|sen(\theta)|[/math]
Los puntos de velocidad máxima son [math]\theta=\frac{\pi}{2}[/math] y [math]\theta=\frac{3\pi}{2}[/math]
Los puntos de velocidad mínima o puntos de remanso son [math]\theta=0[/math] y [math]\theta=\pi[/math].
7 º Presión del fluido
La ecuación de Bernoulli describe como se comporta un fluido, relacionando las siguientes variables:
[math]p =[/math]presión estática a la que está sometido el fluido debida a las moléculas que lo rodean.
[math]d =[/math]densidad del fluido.
[math]\vec { u } =[/math]velocidad de flujo del fluido. :
[math]\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte[/math]
Para que esta expresión pueda aplicarse, tienen que cumplirse las siguientes condiciones:
- - El fluido se mueve en un régimen de tipo estacionario, es decir, la velocidad en un punto determinado no varía con el tiempo.
- - Se desprecia la viscosidad del fluido (que es una fuerza de rozamiento interna).
En nuestro caso, se cumplen la primera condición, y aplicamos la segunda; por lo que se puede aplicar esta ecuación.
En nuestro flujo, se cumple que:
- - La densidad del fluido es [math]d = 2[/math].
- - La constante de la ecuación toma el valor [math]15[/math].
Sustituyendo en la ecuación anterior, obtenemos que: :
[math] |\vec{u}|^2 + p = 15 [/math] :
[math] p = 15- |\vec{u}|^2 [/math]
Por otro lado, hallamos [math]|\vec{u}|^2 [/math]: :
[math] |\vec{u}|^2 = 4(1-\frac{4}{ρ^2})^2 cos^2(θ) + 4(ρ+\frac{4}{ρ})^2 sen^2(θ) [/math]
De esta forma, la presión ejercida por el fluido vendrá determinada por la siguiente ecuación: :
[math] p= 15 - [ 4(1-\frac{4}{ρ^2})^2 cos^2(θ) + 4(ρ+\frac{4}{ρ})^2 sen^2(θ) ] [/math]
Código Matlab:
u=linspace(2,6,50); % Intervalo [2,6].
v=linspace(0,2*pi,50); % Intervalo [0,2*pi].
[U,V]=meshgrid(u,v); % Matrices.
x=U.*cos(V); % Parametrización X.
y=U.*sin(V); % Parametrización Y.
P=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2)))); % Presión estática del fluido.
surf(x,y,P);colorbar; % Visualización de superficie en 3D más leyenda en color.
axis([-5,5,-5,5]) % Selecciona la región a dibujar.
view(2) % Vista desde arriba.
Pmax=max(max(P)) % Presión Máxima.
Pmin=min(min(P)) % Presión Mínima.
El programa nos devuelve que los puntos de presión máxima se encuentran en [math] \theta = 0 [/math] y [math] \theta = \pi [/math]; mientras que los puntos de presión mínima los sitúa en [math] \theta =\frac{\pi}{2} [/math] y [math] \theta =\frac{3\pi}{2} [/math]. Observando los resultados del apartado anterior comprobamos que los puntos de presión máxima coinciden con los puntos de velocidad mínima y viceversa, lo cual es un resultado coherente que verifica la ecuación de Bernoulli.
La ecuación de Navier-Stokes describe el movimiento de los fluidos reales a partir de su campo de velocidades y su campo de presiones. :
[math] d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} [/math]
En nuestro caso, tenemos un fluido incompresible (es decir, que siempre ocupa el mismo volumen, como por ejemplo el agua), cuya viscosidad es [math] \mu =0 [/math] y cuya densidad es [math]d = 2[/math]. Sustituyendo: :
[math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math]
Considerando los campos [math] (\vec{u},p) [/math] hallados anteriormente: :
[math] (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = [/math] :
[math] = 2(1-\frac{4}{\rho^2})cos(\theta)[2(\frac{8}{\rho^3})cos(\theta) \vec{e}_{\rho} - 2(1+\frac{4}{\rho^2})sin(\theta) \vec{e}_{\theta}] - 2(1+\frac{4}{\rho^2})sin(\theta)[2(\frac{8}{\rho^3})sin(\theta) \vec{e}_{\rho} - 2(1+\frac{4}{\rho^2})cos(\theta) \vec{e}_{\theta}] = [/math] :
[math] (\vec{u}·\nabla)\vec{u} = [\frac{32}{\rho^3}(1-\frac{4}{\rho^2})cos^2(\theta) - \frac{32}{\rho^3}(1+\frac{4}{\rho^2})sin^2(\theta)] \vec{e}_{\rho} + [-4(1-\frac{4}{\rho^2})^2 sin(\theta)cos(\theta) + 4(1+\frac{4}{\rho^2})^2 sin(\theta)cos(\theta)] \vec{e}_{\theta} [/math]
De la misma forma: :
[math] \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = [/math] :
[math] \nabla p = [- \frac{64}{\rho^3}(1-\frac{4}{\rho^2})cos^2(\theta) + \frac{64}{\rho^3}(1+\frac{4}{\rho^2})sin^2(\theta)] \vec{e}_{\rho} + [8(1-\frac{4}{\rho^2})^2 sin(\theta)cos(\theta) - 8(1+\frac{4}{\rho^2})^2 sin(\theta)cos(\theta)] \vec{e}_{\theta} [/math]
Si sustituimos los valores obtenidos en la ecuación de Navier-Stokes, se obtiene de forma inmediata que se cumple [math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math], por lo que queda demostrado que se satisface dicha ecuación.
9 º Trayectoria de la línea de corriente
Si fuésemos una partícula del fluido, seguiríamos la trayectoria de una línea de corriente, esto es, según avanza por esta y a medida que se acerca al obstáculo sigue el contorno del mismo para que una vez superado, siga la dirección original de su línea de corriente de nuevo. En cuanto una partícula se va aproximando al obstáculo, la velocidad que ésta tiene es mínima, a diferencia de su presión que es máxima.: Una vez que la partícula ha llegado al obstáculo en donde su velocidad es mínima (nula) y a medida que va avanzando, su presión disminuye a favor de aumentar su velocidad hasta alcanzar su valor máximo. Posteriormente la partícula experimenta una subida de presión a la vez que su velocidad disminuye de nuevo hasta ser nula.
10 º Paradoja D'Alembert
[math]\vec u=2cos(\theta)[1-\frac{4}{\rho^2}]\vec {g}_{\rho}-2sen(\theta)[\frac{1}{\rho}+\frac{4}{\rho^3}]\vec {g}_{\theta}[/math]: Para [math]\rho=2[/math]: [math]\vec u=-2sen(\theta)[\frac{1}{\rho}+\frac{4}{\rho^3}]\vec {g}_{\theta}[/math]: [math]\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} -2sen(\theta)\vec{g}_{\theta}\vec{g}_{\theta}\,d\theta =\int_{0}^{2\pi} -8sen(\theta)\,d\theta=8[cos(2\pi)-cos(0)]=0 [/math]:
El teorema de Kutta-Joukowski establece que la fuerza que ejerce el fluido sobre el obstáculo, en este caso, el círculo de radio 2, es proporcional a esta circulación. Al ser nula, el fluido no ejerce ninguna fuerza sobre el obstáculo, en contra de la intuición. Esto es lo que se conoce como la paradoja de D’Alembert. Para ilustrar esta fuerza es necesario tomar soluciones de modelos de fluidos más sofisticados (con viscosidad µ > 0). Se demuestra así que la circulación del campo [math]\vec u[/math] a lo largo de la circunferencia [math]\rho=2[/math] se anula.
11 º Curvas de nivel de la presión
Por medio de las curvas de nivel podemos comprender de una manera mejor la conclusión a la que habíamos llegado en la cuestión 8 del trabajo. Recordamos que la presión máxima se alcanzaba en [math]\theta=0[/math] y [math]\theta=\pi[/math] y la presión mínima en [math]\theta=\frac{\pi}{2}[/math] y [math]\theta=\frac{3\pi}{2}[/math]. También podemos saber la cantidad de variación de presión, puesto que según estén más juntas las líneas de nivel, la velocidad a la que varíe la presión será mayor.
Código Matlab:
contour(x,y,P,100) % Curvas de nivel
axis([-5,5,-5,5]) % Selecciona la región a dibujar
12 º Presión media
Para calcular la presión media en los puntos del fluido de manera numérica, se aproxima la integral de la presión en todo el fluido y se divide por el área total en el anillo comprendido entre [math]2\lt\rho\lt6[/math]
Código Matlab:
i=1000; % Muestreo
j=1000; % Muestreo
u=linspace(2,6,i+1); % Intervalo [2,6]
v=linspace(0,2*pi,j+1); % Intervalo [0,2*pi]
h=(6-2)/i; H=2*pi/j; % Pasos
[U,V]=meshgrid(u,v); % Mallado
p=15-4.*((((1-(4./(U.^2))).^2).*((cos(V)).^2))+(((U+(4./U)).^2).*((sin(V)).^2).*(1./(U.^2)))); % Función de presiones
P=U.*p; % Presión estática
A=ones(i+1,1); % Matriz auxiliar para integración
A(1)=1/2; A(i+1)=1/2; % Ajuste en los extremos
B=ones(j+1,1); % Matriz auxiliar para integración
B(1)=1/2; B(j+1)=1/2; % Ajuste en los extremos
PT=h*H*B'*P*A; % Presión Total
ro=2; % Radio Interior
roo=6; % Radio Exterior
Ar=pi*(roo^2-ro^2); % Área
PresionMedia=PT/Ar % Presión Media
Una vez ejecutado este código, observamos que la presión media en los puntos del fluido es 10,556.
