Visualización de campos escalares y vectoriales en fluidos. (Grupo 24-C)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en fluidos.
Asignatura	Teoría de Campos
Curso	2022-23
Autores	Erick Morales Pruna Hugo Sacristán de Agustín Yago Amador Page Ismael James Sala Dillon Ignacio del Río Cárdenas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción.

Un fluido es una sustancia que presenta fuerzas de atracción débiles entre partículas y solo es capaz de soportar tensiones normales de compresión. Un fluido ideal es incompresible por lo que no cambia su densidad y por lo tanto tampoco su viscosidad. En el siguiente trabajo presentaremos sus propiedades a partir de un caso en el que suponemos el flujo de un fluido ideal alrededor de un obstáculo circular de radio uno. Haremos uso de campos escalares y vectoriales para su visualización.

2 Representación del mallado de la región ocupada por un fluido incompresible

Dibujaremos un mallado que represente los puntos interiores de la región ocupada por un fluido, que será el exterior del círculo unidad. Para ello tomaremos un mallado del anillo comprendido entre los radios 1 y 5 y centro el origen. Para ilustrar que el fluido ocupa el exterior de un círculo dibujaremos los ejes en el intervalo [−5, 5] × [−5, 5].

%Numero de particiones
n = 40;
%Variables polares
ro = linspace(1,5,n);
tt = linspace(0,2*pi,n);
%Mallado en coordenadas polares
[Mro, Mtt] = meshgrid(ro,tt); 
%Pasamos el mallado a coordenadas cartesianas
x = Mro.*cos(Mtt);
y = Mro.*sin(Mtt);
%Representamos
mesh(x,y,0*x);
axis([-5,5,-5,5]);
view(2);

Figura 1: Mallado de la región ocupada por un fluido.

3 Velocidad del fluido y función potencial

3.1 Representación de la velocidad

Figura 2: Gradiante de la velocidad de las partículas de agua.

La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial

[math] φ(ρ,θ,z) = (ρ + \frac{1}{ρ}) \cos(θ)+\sqrt{2}θ[/math]

A continuación, calcularemos y representaremos la función potencial [math]φ[/math] y el campo de velocidades del fluido. Observar gráficamente (dibujando un zoom de la gráfica) que [math]\vec u = ∇φ[/math] es ortogonal a las curvas de nivel de [math]φ[/math].

La fórmula del gradiente de cualquier campo escalar expresado en coordenadas de la base física cilíndrica es:

[math]∇f (ρ,θ,z) = (\frac{∂f}{∂ρ})\vec e_ρ + \frac{1}{ρ}(\frac{∂f}{∂θ})\vec e_θ + (\frac{∂f}{∂z})\vec e_z [/math]

Sustitiyendo el campo escalar de velocidades [math]φ[/math] en la fórmula del gradiente obtenemos lo siguiente:

[math] ∇φ(ρ,θ,z) = \frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ + \frac{\sqrt{2}ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ [/math]

Esto es el campo de velocidades del fluido sobre y alrededor del obstáculo circular representadas en el gráfico a la derecha generada con el código de MATLAB:

%Definimos nuestra función y sus gradientes
f = @(ro,tt) (cos(tt).*(ro + 1./ro)+sqrt(2).*tt);
Gro = @(ro,tt) (cos(tt)-(cos(tt)./ro.^2));
Gtt = @(ro,tt) (-sin(tt)).*(ro.^2 + 1)./ro.^2 + sqrt(2)./ro;
%Definimos variables en polares y su mallado
n =20; m=60;
ro = linspace(1,50,n);
tt = linspace(0,2*pi,n);
[Mro,Mtt] = meshgrid(ro,tt); 
%Pasamos a coordenadas cartesianas
x = Mro.*cos(Mtt);
y = Mro.*sin(Mtt);
z = f(Mro,Mtt);
fro = Gro(Mro,Mtt);
ftt = Gtt(Mro,Mtt);
fx = (fro.*x - ftt.*y);
fy = (fro.*y + ftt.*x);
%Representamos
 quiver(x,y,fx,fy);
 hold on
 contour(x,y,z,m)

3.2 Gradiante ortogonal a curvas de nivel

Figura 3: Acercamiento Figura 2.

Podemos observar en el acercamiento de la gráfica, que los vectores velocidad del gradiente del campo escalar [math]φ[/math] son ortogonales a las curvas de nivel. Esto se debe a que el flujo del fluido de una curva de nivel a otra sigue la trayectoria de la línea de máxima pendiente. Recorre así la menor distancia posible para cambiar de curva de nivel.

3.3 Comprobación de perpendicularidad

Considerando [math]∇φ(ρ,θ,z) = \vec u[/math] y si denotamos [math]\vec n[/math] al vector normal a los puntos del obstáculo, comprobaremos que [math]\vec u*\vec n=0[/math]. En la Figura 3 se puede apreciar la perpendicularidad de los vectores gradiente a las curvas de nivel de la función potencial. La trayectoria del gradiente viene definida por los puntos de la gráfica los cuales tienen mayor incremento. Gráficamente son las rectas de máxima pendiente entre curvas de nivel. Por lo tanto, el producto de [math]\vec u*\vec n[/math] es igual a cero.

[math] \vec n = \vec e_z [/math]
[math] \vec u = ∇φ(ρ,θ,z) = \frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ + \frac{\sqrt{2}ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ + 0\vec e_z [/math]

[math] \vec n·\vec u = (0\vec e_ρ+0\vec e_θ+1\vec e_z)·(\frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ + \frac{\sqrt{2}ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ + 0\vec e_z) = \vec 0[/math]

4 Aproximación de [math]∇φ[/math]

Si estamos lejos del obstáculo, [math]ρ[/math] es muy grande y podemos suponer [math]\frac{1}{ρ}[/math] despreciable. ¿Cuánto vale aproximadamente [math]\vec u[/math] (escrito en coordenadas cartesianas) en esos puntos?
Suponiendo que [math]ρ[/math] es muy grande:

[math]\frac{1}{ρ} = 0[/math]

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

Expresamos el vector [math]\vec u[/math] en coordenadas cartesianas:

[math]\begin{pmatrix} v_1 \\ v_2 \\ v_3\end{pmatrix} = \begin{pmatrix}cos(θ) & -sen(θ) & 0 \\ sen(θ) & cos(θ) & 0 \\ 0 & 0 & 1\end{pmatrix} \begin{pmatrix}cos(θ) \\ -sen(θ) \\ 0\end{pmatrix} = \begin{pmatrix}cos^2(θ)+sen^2(θ) \\ cos(θ)sen(θ)-cos(θ)sen(θ) \\ 0\end{pmatrix} = \vec u(x_1,x_2,x_3)=\vec i[/math]

Esto significa que el vector gradiente, y por lo tanto los vectores velocidad del fluido irían en la dirección generalizada de [math]\vec i[/math].

5 Rotacional y Divergencia de [math]∇φ[/math]

A continuación, comprobaremos analíticamente que [math]\vec u[/math] tiene rotacional y divergencia nulos comprobando así la incompresibilidad del fluido a nivel local.

Rotacional en coordenadas cilíndricas:[math]∇×\vec u = \frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ u_ρ & ρu_θ & u_z \end{pmatrix}[/math]

5.1 Rotacional de [math]∇φ[/math]

[math]∇×\vec u(ρ,θ,z) = \frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ \frac{(ρ^2-1)cos(θ)}{ρ^2} & ρ\frac{\sqrt{2}ρ-(ρ^2+1)sen(θ)}{ρ^2} & 0\end{pmatrix} = 0\vec e_ρ-0\vec e_θ+[\frac{∂}{∂ρ}(\frac{\sqrt2ρ}{ρ}-\frac{(ρ^2+1)sen(θ)}{ρ})-\frac{∂}{∂θ}(\frac{(ρ^2-1)cos(θ)}{ρ^2})]\vec e_z[/math]

[math] = [0-sen(θ)\frac{∂}{∂ρ}(\frac{ρ^2+1}{ρ})+\frac{(ρ^2-1)sen(θ)}{ρ^2}]\vec e_z = [-sen(θ)\frac{(ρ^2-1)}{ρ^2}+sen(θ)\frac{(ρ^2-1)}{ρ^2}]\vec e_z = 0\vec e_z[/math]

5.2 Divergencia de [math]∇φ[/math]

[math]∇·\vec u(ρ,θ,z) = \frac{1}{ρ}[\frac{∂}{∂ρ}(ρ\vec u_ρ)+\frac{∂}{∂θ}(\vec u_θ)+\frac{∂}{∂z}(ρ\vec u_z)][/math]

[math]\vec u = \frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ+\frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ [/math]

[math]∇·\vec u(ρ,θ,z) = \frac{1}{ρ}[cos(θ)\frac{∂}{∂ρ}(\frac{ρ^2-1}{ρ})+(\frac{-(ρ^2+1)}{ρ^2})\frac{∂}{∂θ}sen(θ)] = \frac{1}{ρ}[cos(θ)\frac{ρ^2+1}{ρ}-cos(θ)\frac{ρ^2+1}{ρ}] = 0[/math]

6 Lineas de corriente

Definimos como líneas de corriente como aquellas que son tangentes a cada punto del campo [math]\vec u = ∇φ (ρ,θ,z)[/math], es decir, ortogonales a un vector [math]\vec v[/math] que lo calculamos de forma [math]\vec v = \vec k × \vec u[/math].
Siendo: [math]\vec k = \vec e_z[/math]

[math]\vec u = ∇φ (ρ,θ,z) = \frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ+\frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ+0\vec e_z[/math]

Realizamos el producto vectorial:

[math]\vec v = \vec k×\vec u = \begin{pmatrix} \vec{e_ρ} & \vec{e_θ} & \vec{e_z} \\ 0 & 0 & 1 \\ \frac{(ρ^2-1)cos(θ)}{ρ^2} & \frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2} & 0 \end{pmatrix} = ... = [/math] [math]\vec v = -\frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_ρ+\frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_θ [/math]

Sabemos que el vector [math]\vec u[/math] tiene divergencia nula ([math]∇∙\vec u=0[/math]), por lo que el vector [math]\vec v[/math] sera irrotacional ([math]∇×\vec u=0[/math]):
[math]∇×\vec v=0[/math], como ya hemos visto en el apartado anterior.


La función potencial, definimos: [math]∇ψ = \vec v[/math]

[math] ∇ψ = \frac{∂ψ}{∂ρ}+\frac{1}{ρ}\frac{∂ψ}{∂θ}+\frac{∂ψ}{∂z}[/math]

• [math] \frac{∂ψ}{∂ρ} = -\frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2} [/math] → [math] ψ=-∫\frac{\sqrt2ρ-(ρ^2+1)sen(θ)}{ρ^2}dρ = ... = \frac{-sen(θ)(\sqrt2ρ*ln(ρ)-ρ^2+1)}{ρ}+f(θ)[/math]
• [math] \frac{1}{ρ}\frac{∂ψ}{∂θ} : ψ = ∫(\frac{-sen(θ)(\sqrt2ρ*ln(ρ)-ρ^2+1)}{ρ}+f(θ))dθ = \frac{cos(θ)(\sqrt2ρ*ln(ρ)-ρ^2+1}{ρ}+f'(θ) [/math] →
  

→ [math] \frac{1}{ρ}ψ = \frac{cos(θ)(\sqrt2ρ*ln(ρ)-ρ^2+1}{ρ^2}+\frac{f'(θ)}{ρ} [/math]

• [math] \frac{∂ψ}{∂z} = 0 [/math]

Nos damos cuenta que: [math]\frac{f'(θ)}{ρ}=0[/math] → [math]f'(θ)=0[/math] → [math]f(θ)=cte[/math]

[math] ψ=cte-\frac{-sen(θ)(\sqrt2ρ*ln(ρ)-ρ^2+1)}{ρ} [/math]

Con todo esto, nos damos cuenta que las líneas de corriente son tangentes al campo.

Figura 4: Lineas de corriente respecto al campo.

Código de MATLAB creado para observar el resultado de nuestro campo:

%Definimos nuestra función y sus gradientes
 f = @(ro,tt) (cos(tt).*(ro + 1./ro)+sqrt(2).*tt);
 Gro = @(ro,tt) (cos(tt)-(cos(tt)./ro.^2));
 Gtt = @(ro,tt) (-sin(tt)).*(ro.^2 + 1)./ro.^2 + sqrt(2)./ro;
 fc = @(ro,tt) ((sin(tt)./ro).*(((sqrt(2).*ro).*log(ro))-ro.^2+1));
 %Definimos variables en polares y su mallado
 n =30;
 ro = linspace(1,5,n);
 tt = linspace(0,2*pi,n);
 [Mro,Mtt] = meshgrid(ro,tt);
 %Pasamos a Coordenadas cartesianas
 x = Mro.*cos(Mtt);
 y = Mro.*sin(Mtt);
 z = f(Mro,Mtt);
 z2 = fc(Mro,Mtt);
 fro = Gro(Mro,Mtt);
 ftt = Gtt(Mro,Mtt);
 fx = (fro.*x - ftt.*y);
 fy = (fro.*y + ftt.*x);
 %Representamos
 quiver(x,y,fx,fy);
 hold on
 contour(x,y,z2,n)

7 Velocidades máximas y mínimas en frontera

Consideramos como "S" a la frontera entre el obstáculo y el fluido. Vamos a calcular los puntos de "S" donde el módulo de la velocidad del fluido es máxima, mínima y los puntos de remanso (donde el módulo de la velocidad es nulo).

Consideramos:

[math]ρ = 1[/math]
[math]\vec u = \frac{(ρ^2-1)cos(θ)}{ρ^2}\vec e_ρ + \frac{\sqrt{2}ρ-(ρ^2+1)sen(θ)}{ρ^2}\vec e_θ + 0\vec e_z[/math]

Figura 5: Gráfico puntos máximo, minimo y remanso

Entonces: [math]\vec u = (\sqrt2-2sen(θ))\vec e_θ[/math]

Podemos tomar en cuenta el intervalo [math]-1 ≤ cos(θ) ≤ 1[/math] ya que el valor del coseno tiene como máximo 1 y mínimo -1. Los valores extremos se encuentran en [math]θ_1 = 0[/math] , [math]θ_2 = π[/math]. Para encontrar los puntos de remanso consideramos [math]cos(θ) = 0[/math] ➔ [math]θ_3 = \frac{π}{2}[/math] , [math]θ_4 = \frac{3π}{2}[/math]. Todo esto se puede apreciar en la gráfica adyacente.

8 Ecuación de Bernoulli

8.1 Presion

Supondremos que la densidad del fluido es constante [math]d = 2[/math], y que se verifica la ecuación de Bernoulli [math]\frac{1}{2}d|\vec u|^2+p=cte[/math] Calcularemos y representaremos la presión del fluido dando el valor 10 a la constante anterior.

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

Por lo tanto

[math]p=-(1+\frac{4cos^2(θ)}{ρ^2}-\frac{2\sqrt{2}sen(θ)(ρ^2-1)}{ρ^3}+\frac{1}{ρ^4})+10[/math]

Con el siguiente código de MATLAB podemos graficar las presiones del fluido alrededor del obstáculo. La presión máxima es 10 y la mínima es 4.

Figura 6: Presiones alrededor del obstáculo en 3D

Figura 7: Vista aérea presiones alrededor del obstáculo

%Definimos nuestra función y sus gradientes
 f = @(ro,tt) (cos(tt).*(ro + 1./ro)+sqrt(2).*tt);
 Gro = @(ro,tt) (cos(tt)-(cos(tt)./ro.^2));
 Gtt = @(ro,tt) (-sin(tt)).*(ro.^2 + 1)./ro.^2 + sqrt(2)./ro;
 fc = @(ro,tt) ((sin(tt)./ro).*(((sqrt(2).*ro).*log(ro))-ro.^2+1));
 %Definimos variables en polares y su mallado
 n =30;
 ro = linspace(1,5,n);
 tt = linspace(0,2*pi,n);
 [Mro,Mtt] = meshgrid(ro,tt);
 %Pasamos a Coordenadas cartesianas
 x = Mro.*cos(Mtt);
 y = Mro.*sin(Mtt);
 
 fro = Gro(Mro,Mtt);
 ftt = Gtt(Mro,Mtt);
 fx = (fro.*x - ftt.*y);
 fy = (fro.*y + ftt.*x);
 % Definimos la función de la presión
  p =-(1.+(4.*cos(Mtt).^2./Mro.^2)-(2.*sqrt(2).*sin(Mtt).*((Mro.^2)-1)./Mro.^3)+(1./Mro.^4))+10;
% Representamos la presión del fluido
 figure(1)
 surf(x,y,p)
 axis([-5,5;-5,5])
 axis equal
 view(2)
 figure(2)
 surf(x,y,p)
 axis([-5,5;-5,5])
 axis equal
 view(3)

9 Navier-Stokes estacionaria

Partiendo de la ecuación de Bernouilli comprobar que [math](\vec u, p)[/math] satisface la ecuación de Navier-Stokes estacionaria

[math]d(\vec u·∇)\vec u+∇p = µ∆\vec u[/math] para [math]µ=0[/math] ([math]µ[/math] esta asociado a la viscocidad del fluido), donde [math](\vec u·∇)\vec u = u_j\frac{∂u_i}{∂x_j}\vec e_i[/math]
[math]∆\vec u = ∆u_i\vec e_i[/math]
[math]\vec u = u_i\vec e_i[/math].

Para ello calcularemos el gradiente de la ecuación de Bernouilli y verificaremos la siguiente identidad vectorial usando coordenadas cartesianas que además se verifica la condición de incompresibilidad ([math]∇·\vec u = 0[/math])

[math](\vec u·∇)\vec u = \frac{1}{2}∇(|\vec u|^2)-\vec u×∇×\vec u[/math]

Las ecuaciones de Navier-Stokes son ecuaciones diferenciales parciales que describen el movimiento de cuerpos fluidos viscosos. Consisten en una ecuación de energía, otra de conservación del momento, de la masa y la continuidad. Sin embargo al estar trabajando con un fluido ideal, despreciamos las ecuaciones que tienen que ver con la conservación de la masa y la cantidad de movimiento para describir la mecánica de el fluido en cuestión.

[math]\frac{1}{2}d|\vec u|^2+p=cte[/math] ➔ [math]\frac{1}{2}d∇(|\vec u|^2)+∇p=∇(cte)=0[/math] ➔
➔ [math]\frac{1}{2}d∇(\vec u·\vec u)+∇p=0[/math] ➔ [math]\frac{1}{2}2d(∇·\vec u)\vec u+∇p=0[/math] ➔ [math]d(∇·\vec u)\vec u+∇p=0[/math]

10 Cambios cinemáticos y de presión sobre el obstáculo

Las partículas de fluido siguen la trayectoria de las líneas de corriente. No es hasta que se presenta un obstáculo que cambia su dirección, velocidad y presión. Vamos a Deducir esto a partir de las gráficas de las líneas de corriente, presión y módulo de la velocidad [math]|\vec u|[/math].

Partimos entonces de la expresión de la presión:

[math]p=-(1+\frac{4cos^2(θ)}{ρ^2}-\frac{2\sqrt{2}sen(θ)(ρ^2-1)}{ρ^3}+\frac{1}{ρ^4})+10[/math]

[math]∇p=[\frac{8cos^2(θ)}{ρ^3}+\frac{4\sqrt{2}ρ^2sen(θ)-6\sqrt{2}sen(θ)(ρ^2-1)}{ρ^4}+\frac{4}{ρ^5}]\vec e_ρ + [\frac{4sen(2θ)}{ρ^3}+\frac{2\sqrt{2}cos(θ)(ρ^2-1)}{ρ^4}]\vec e_θ[/math]

Cuando la presión aumenta, la velocidad disminuye con:

[math]\frac{1}{2}d|\vec u|^2+p = cte ➔ |\vec u|=\sqrt{\frac{2(cte-p)}{d}} [/math]

Visto lo anterior podemos definir que el aumento de la presión del fluido es inverso al aumento de la velocidad del mismo. La máxima presión se encuentra en [math]P_1(ρ,θ)=(1,\frac{π}{2})[/math] y [math]P_2(ρ,θ)=(1,\frac{3π}{2})[/math], por lo tanto la velocidad mínima se encuentra en los mismos puntos.

11 Paradoja de D'Alembert

La integral del campo p [math]\vec n \vec i[/math] a lo largo de la circunferencia p = 1 se anula ya que al ser p [math]\vec n[/math] la fuerza que ejerce el fluido en cada punto de la frontera, este vector nos da como resultado un vector compuesto por tres componentes: a[math]\vec i[/math] + b [math]\vec j[/math] + c [math]\vec k[/math]; centrándonos en la a es fácil comprobar que esta va a dar 0, ya que contiene dos partes iguales de distinto signo. En cuanto a las otras también se nos van al multiplicar dicho vector por [math]\vec i[/math] ya que solo queremos saber la fuerza que ejerce en esta dirección. Como resultado nos da la paradoja de D’Alembert.

[math]∫_S p\vec n\vec i dS = ∫_S (p\vec n)\vec i dS = ∫_S \vec 0 dS = 0[/math]

Otra manera de llegar a esta paradoja sería aplicando el teorema de Stokes lo que con un rotacional nulo nos daría el mismo resultado.

[math] ∫_S (∇ × \vec u) \vec n · dS = ∫_S (0)\vec n · dS = 0 [/math]

La paradoja de D’Alembert es como bien indica su nombre una paradoja creada por el matemático francés Jean Le Rond D'Alambert la cual llega a la contradicción escrita antes. Es imposible que a un cuerpo que se mueve a velocidad constante en un fluido no le afecte ninguna fuerza impuesta por dicho fluido, es decir que no actuarán fuerzas en la dirección del movimiento sobre el móvil; como se puede comprobar a simple vista esto no tiene ningún tipo de sentido, es como decir que el aire no ejerce ninguna fuerza sobre un avión en movimiento por ejemplo. Sin embargo, el consejo científico de la época dictó que este resultado se debió a la falta de viscosidad del fluido en los cálculos. Para obtener esta fuerza es necesario tomar soluciones de la ecuación de Navier Stokes con viscosidad μ > 0.

[math]d(\vec u·∇)\vec u+∇p = µ∆\vec u ≠ 0[/math]

12 Representación de presiones

El dibujo de las curvas de nivel de presión se pondrá a continuación junto al código en lenguaje M que lo proporciona.

Figura 8: Representación presiones por curvas de nivel

%Definimos la función y su gradiante
f = @(ro,tt) (cos(tt).*(ro + 1./ro)+sqrt(2).*tt);
Gro = @(ro,tt) (cos(tt)-(cos(tt)./ro.^2));
Gtt = @(ro,tt) (-sin(tt)).*(ro.^2 + 1)./ro.^2 + sqrt(2)./ro;
%Discretizamos las variables
o =100;
ro = linspace(1,10,o);
tt = linspace(0,2*pi,o);
[Mro,Mtt] = meshgrid(ro,tt); 
%Pasamos los valores y su gradiante a cartesianas
x = Mro.*cos(Mtt);
y = Mro.*sin(Mtt);
fro = Gro(Mro,Mtt);
ftt = Gtt(Mro,Mtt);
%Definimos y representamos
fx = (fro.*x - ftt.*y)./sqrt(x.^2 + y.^2);
fy = (fro.*y + ftt.*x)./sqrt(x.^2 + y.^2);
p = -(1+(4.*cos(Mtt).^2./Mro.^2)-(2.*sqrt(2).*sin(Mtt).*((Mro.^2)-1)./Mro.^3)+(1./Mro.^4));
z = 10-v.^2;
contour(x,y,z,o)

Gracias a este sencillo gráfico es muy fácil comprobar que donde la velocidad es menor se aprecia un notable aumento de la presión. Como podemos comprobar en la gráfica, las líneas de flujo van reduciendo su velocidad según se acercan al cilindro de radio 1, aumentando la presión; pero en el momento en el que llegan a este radio empiezan a incrementar su velocidad mientras lo rodean alcanzando puntos máximos en [math]ρ=1[/math] y [math]θ=\frac{π}{2}[/math] para posteriormente volver a reducir su velocidad terminando de rodear este objeto.