Ver el código de «Fluido alrededor de un obstáculo circular (Grupo 20A)»


Un fluido es un medio continuo tal que en reposo no transmite tensión tangencial. No es capaz soportar su propio peso, y tiende a adoptar las características del recipiente en el que se encuentra.

Los fluidos pueden ser líquidos o gases, se van a estudiar los fluidos incompresibles, que son aquellos cuyo volumen no disminuye al ejercer una fuerza sobre él, este tipo de fluidos son los líquidos. 

En este artículo se estudiará el comportamiento de estos fluidos al fluir alrededor de un obstáculo con forma circular. 
Por conveniencia se trabajará con coordenadas cilíndricas.


{{ TrabajoED | Fluido alrededor de un obstáculo circular (Grupo 20A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Nicolás Martín de Andrade  Andrea Amo AllasㅤㅤㅤㅤJorge Mancheño EstévezㅤㅤPedro Pérez HernanzㅤㅤÁlvaro Román Aguilera }}

==Mallado==
Se empezará dibujando un mallado que represente los puntos interiores de la región ocupada por un fluido, que será el exterior del círculo unidad. Para llevar a cabo esta gráfica, se tomará un mallado del anillo comprendido entre los radios 1 y 5 y centro el origen.

La finalidad de realizar un mallado es dividir el dominio continuo (la región ocupada por el fluido) en subdominios más pequeños llamados elementos finitos o celdas. Estos elementos forman una malla que cubre toda la región y permite realizar cálculos y simulaciones numéricas en cada uno de ellos.
 
Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-5,5]×[-5,5]</math>.

[[Archivo:mallado56.jpg|thumb|right|Mallado del obstáculo circular]]
{{matlab|codigo=

rho=linspace(1,5,30);                        %Se definen rho y theta 
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                      %Se crea el mallado

hold on
X=U.*cos(V);                                 %Se parametriza la superficie
Y=U.*sin(V);
Z=0.*U;
mesh(X,Y,0*Z);                               %Dibujo de la placa 
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Se representa el obstáculo
axis([-5,5,-5,5]);                           %Se fijan los ejes 
view(2);
title ('Placa');
xlabel 'EJE X'
ylabel 'EJE Y'
axis equal
hold off
}}

==Función Potencial del Fluido==
Una vez definida la región que se va a estudiar, se analizará la velocidad de las partículas del fluido, que viene dada por su función potencial. 

La función potencial establecida es <math> \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta )- \theta  </math>. Se representará la función potencial, para luego examinar la velocidad de las partículas y el movimiento del fluido.
[[Archivo:funcionpotencial2.jpg|thumb|right|Mallado del obstáculo circular]]
{{matlab|codigo=
rho=linspace(1,5,30);                             %Se define rho y theta
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                           %Se crea el mallado

hold on
X=U.*cos(V);                                      %Se parametriza la superficie
Y=U.*sin(V);

f=@(rho,th)(rho+(1./rho)).*cos(th)-th;  %Se define la función potencial
Z=f(U,V);                                         %Se aplica la función potencial
surf(X,Y,Z);                                      %Se representa la función
plot(1*cos(th),1*sin(th),'k','lineWidth',3);      %Se representa el obstáculo

view(2);  
axis([-5,5,-5,5]);
colorbar;
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}

==Campo de Velocidades del Fluido==
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial. Esto determina el valor y la dirección de máximo crecimiento de la función potencial.

El gradiente de una función en coordenadas cilíndricas se obtiene a partir de la siguiente fórmula: 

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Siguiendo esta fórmula, el gradiente de la función potencial determinada es: 

<center><math> \vec u=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta) + \frac{1}{\rho}]\vec e_\theta  </math></center>


Como la representación se realiza en coordenadas cartesianas, se debe pasar el gradiente de la función potencial a estas. Para ello usamos la matriz de cambio de base de cilíndricas a cartesianas y la multiplicamos por el vector dado en coordenadas cilíndricas.


Matriz de cambio de base:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta  & 0\\
0 & 0 & 1
\end{pmatrix}


Tras el cambio obtenemos;

<center><math> \vec u=\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) +  (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)  + \frac{sin(\theta)}{\rho}\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - \frac{cos(\theta)}{\rho} \right) \vec j </math></center> 


A continuación, se representa el campo de velocidades.
[[Archivo:campodevelocidades2.jpg|thumb|right|Campo de Velocidades]]
{{matlab|codigo=
rho=linspace(1,5,30);                              %Definimos rho y theta
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                            %Creamos la malla

X=U.*cos(V);                                       %Parametrizamos la superficie
Y=U.*sin(V);

f=@(rho,th)(rho+(1./rho)).*cos(th)-th;  %Definimos la función potencial 

Z=f(U,V);                                          %Aplicamos la función

contour(X,Y,Z,30);                                 %Dibujamos las curvas de nivel
hold on
                                                   %Definimos las componentes X e Y del gradiente
Cx=(1-(1./U.^2)).*cos(V).^2+(1+(1./U.^2)).*sin(V).^2+(sin(V)./rho);           
Cy=(1-(1./U.^2)).*sin(V).*cos(V)-(1+(1./U.^2)).*sin(V).*cos(V)-(cos(V)./rho);
quiver(X,Y,Cx,Cy);                                 %Dibujamos el campo de velocidades 
plot(1*cos(th),1*sin(th),'k','lineWidth',1);       %Representación del obstáculo
axis([-5,5,-5,5]);
colorbar;                                          %Añadimos una barra de color
axis equal 
view(2);
hold off
}}
Ahora más de cerca comprobamos que \( u = \nabla\varphi \) es ortogonal a las curvas de nivel de \( \varphi \). Esto representa la dirección de máxima tasa de variación. 
[[Archivo:campodevelocidades2cerca.jpg|Zoom Campo de Velocidades]]




===Interpretación y Valor de  <math>\vec u </math>===
Utilizando la gráfica anterior, se puede visualizar como el campo de velocidades <math> \vec u=\nabla \varphi </math>, representado por las flechas, es ortogonal en todo punto a las curvas de nivel de <math> \varphi </math>. 

Para comprobar matemáticamente esta observación, se elige el vector <math> \vec n </math> normal a los puntos del obstáculo, y se verifica que cumple la fórmula <math>\vec{u}\cdot \vec{n}=0 </math>. Esto demuestra que, efectivamente los vectores son ortogonales entre sí.

Se evalúa el gradiente <math> \vec u </math> en el punto <math>\rho = 1</math>, y por ello se elimina la componente <math>\vec e_\rho </math> del gradiente. De esta forma, si multiplicamos el gradiente por <math> \vec n </math>, el vector normal, que en este caso es  <math>-\vec e_\rho </math> se tiene que:

<center><math>-\vec e_\rho \cdot \left ( -2\sin (\theta) -1 \right )\vec e_\theta = 0 </math></center>

Se evalúa <math> \vec u </math> en la frontera del obstáculo:

<center><math>\vec u(1,\theta,z)=\nabla \varphi(1,\theta,z)= \left (-2\sin(\theta)-1 \right)\vec e_\theta </math></center>

===Campo de velocidades lejos del obstáculo===
Para poder estudiar la función estando lejos del obstáculo, se hará la siguiente suposición: al estar lejos del obstáculo circular, <math>\rho </math> es muy grande, y se puede suponer que <math> \frac{1}{\rho } </math> es despreciable.

Trabajamos sobre la función potencial establecida: 

<center><math> \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) - \theta  </math></center>

Al despreciar <math> \frac{1}{\rho } </math> se obtiene:

<center><math> \varphi* = \rho \cos (\theta) - \theta </math></center>

Seguidamente calculamos el gradiente de la función potencial hallada:

<center><math>\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z= \nabla \varphi*=
</math>
<math>\cos(\theta) \vec{e}_{\rho} - (\sin(\theta) + \frac{1}{\rho}) \vec{e}_{\theta}</math></center>

===Rotacional Nulo y Divergencia Nula===
El rotacional y la divergencia del campo vectorial proporcionan información sobre las características físicas del fluido.

El operador rotacional indica la dirección y la velocidad del giro del campo en cada punto. Considera el movimiento del fluido y su tendencia a inducir rotación en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta) - \frac{1}{\rho}) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right)  \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) + 1 \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Obtenemos que el rotacional es 0, por lo que las partículas del fluido no están girando.



La divergencia en cambio, mide la tasa de expansión y contracción de un fluido en cada punto, calculando la diferencia entre el flujo entrante y saliente de la superficie.

A continuación se demuestra que la divergencia es nula:

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho^2})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) - \frac{1}{\rho} \right) + 0 \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

Obtenemos que al igual que el rotacional, la divergencia es también 0, lo que indica que el volumen del fluido es constante, es decir, es un fluido incompresible.

==Líneas de Corriente==
En primer lugar, se calcularán las lineas de corriente del campo <math>\vec{u}</math>, donde las líneas de corriente son aquellas que son tangentes al campo <math>\vec{u}</math> en cada punto. Dichas líneas son las cuales establecen la trayectoria que sigue las partículas del fluido. 

Para su realización, se calculará el campo que es ortogonal en cada punto a <math>\vec{u}</math>. Este vector se asignará como  <math>\vec{v}</math> y se calculará con la siguiente fórmula:  <math>\vec{v}=\vec{k}\times \vec{u}</math>

En el apartado 3.3, se puede apreciar que el rotacional es nulo y esto implica  <math>\vec{v}</math> sea irrotacional. A su vez, tiene un potencial escalar <math>\psi </math> el cual tiene  <math>\vec{v}</math> de gradiente. Este potencial es conocido como la función de corriente de <math>\vec{u}</math>. 

Para calcular dicho potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> para finalmente dibujar las líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta) - \frac{1}{\rho}) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta}+\frac{1}{\rho})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

Resolvemos el sistema del potencial sabiendo que:

<math>\frac{\partial \psi}{\partial \rho} = \vec{v}_\rho = ((1+\frac{1}{\rho^2})\sin{\theta}+\frac{1}{\rho}\vec{e}_\rho</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \vec{v}_\theta = (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<math>\frac{\partial \psi}{\partial Z} = \vec{v}_Z = 0</math>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta} + \frac{1}{\rho})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho + \int\frac{1}{\rho}\partial\rho= \sin{\theta}(\rho-\frac{1}{\rho})+ log({\rho}) + \delta(\theta,Z)</math>

<math>\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})+\frac{\partial\delta}{\partial\theta}=\rho\cos{\theta}(1-\frac{1}{\rho^2})</math>

<math>\frac{\partial \delta(\theta,Z)}{\partial \theta} = 0</math>

<math>\delta(\theta,Z) = cte</math>

<math>\psi(\rho,\theta,Z) = \log(\rho)+\sin{\theta}(\rho-\frac{1}{\rho})-cte</math>


Finalmente, obtenemos la ecuación <math>\psi: \log(\rho)+\sin{\theta}(\rho-\frac{1}{\rho})=cte</math>


[[Archivo:partepedro.jpg|miniaturadeimagen|derecha|Lineas de corriente]]
{{matlab|codigo= 
rho=linspace(1,5,30);
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                                %Creamos la malla

X=U.*cos(V);                                           %Parametrizamos la superficie
Y=U.*sin(V);

psi=@(rho,th)(log(U))+(sin(V)).*(U-(1./U));  %Definimos la función  ψ (psi)

Z=psi(U,V);                                            %Aplicamos la función

contour(X,Y,Z,30);                                     %Dibujamos las líneas de corriente 
hold on

%Dibujamos el campo de velocidades (gradiente de phi)

Cx=((cos(V).^2).*(1-(1./U.^2)))+(sin(V).^2).*(1+(1./U.^2)+(sin(V)./U));
Cy=sin(V).*cos(V).*(1-(1./U.^2))-sin(V).*cos(V).*(1+(1./U.^2))-(cos(V)./U);


quiver(X,Y,Cx,Cy); 

plot(1*cos(th),1*sin(th),'k','lineWidth',1);           %Representación del obstáculo
axis([-5,5,-5,5]);
colorbar;
xlabel 'Eje X';
ylabel 'Eje Y';
axis equal 
view(2);
hold off}}

Posteriormente comprobaremos que las líneas de corriente de <math>\vec{u}</math> son tangentes a la velocidad del fluido. 

[[Archivo:partepedro1.jpg|Líneas de corriente]]


Del  mismo modo comprobamos que las líneas de corriente <math>\psi</math> son ortogonales a las curvas equipotenciales.

[[Archivo:pedro2bien.jpg]]

==Puntos de la Frontera con Velocidades Máximas y Mínimas==
Seguidamente, donde el módulo de la velocidad es mayor y menor, se calcularán los puntos de la frontera del obstáculo S. Los puntos de remanso, serán en los que la velocidad es nula. En el borde del obstáculo de la gráfica se señalarán los puntos de remanso. 

Debido a que   <math> \vec{u} =\nabla\varphi </math> , se tiene que hallar el módulo de la velocidad,  <math> \left \| \vec u \right \|</math>. Este módulo es un escalar, y se analiza en la región S= <math> \left \{(\rho,θ): \rho=1 \right \}</math>. Al calcular el módulo en dicha región, se obtiene que: <math> \left \| \vec u \right \|= -2\sin(\theta)-1  </math>

Para encontrar los puntos de inflexión, es decir, los máximos y mínimos, se deriva la ecuación: 
<math> \frac {\partial}{\partial\theta}\left ( -2sin\theta - 1 \right )=0 \Rightarrow cos(\theta)=0 \Rightarrow \theta=3π/2 </math> y <math>\theta=π/2 </math>

Por tanto, obtenemos los siguientes puntos de inflexión; <math> \theta=3π/2 \text{ y } \theta=π/2</math>.

Ya que el módulo de la velocidad es dado en valor absoluto, su valor mínimo posible es <math> \left \| \vec u \right \|=0 \Rightarrow \left \| \vec u \right \|= -2\sin(\theta)-1=0 </math>. 

Con estos sacamos los mínimos: <math> \theta=7π/6 \text{ y } \theta=11π/6</math>

En conclusión, con ayuda del grafico: 
	el máximo absoluto se da en: <math> \theta=π/2 </math> y los mínimos, es decir los punto de remanso se dan en: <math> \theta=7π/6 \text{ y } \theta=11π/6</math>

[[Archivo:punto66.jpg|thumb|right|Como cambia el módulo de la velocidad en el borde del obstáculo al aumentar θ]]

{{matlab|codigo=

x=linspace(0,2*pi,1000);
y=-2.*sin(x)-1;
a=1;

%Se cambian de signo los valores negativos

for i=y
    if i<0
        i=-i;
    end
    y(a)=i;
    a=a+1;
end
plot(x,y)
xlabel('ángulo')
ylabel('módulo de la velocidad')
}}

==Ecuación de Bernoulli==
Para poder calcular los puntos de mayor y menor presión del fluido, se supondrá que la densidad del fluido es constante <math> d = 2 </math>, tal que se verifica la ecuación de Bernoulli: 
<math> 1/2 d \left \| \vec u \right \|{^2}=cte </math>

Se le asignará el valor de 10 a la constante anterior, para calcular la presión del fluido.

Aplicando los valores anteriores, la ecuación quedará de la forma: <math> p=10-\left \| \vec u \right \|{^2} </math>. 
Debido a que se está observando la región del obstáculo S, se tiene que <math> \left \| \vec u \right \|{^2} </math> queda en función del parámetro <math> \theta </math>.

<center><math> p=10- ( -2sinθ - 1 ){^2} = 9 - 4sin{^2}\theta -  4sin\theta</math></center>


Para hallar los máximos y mínimos, se deriva la función anterior, y se iguala a 0. 

<center><math>\frac {\partial}{\partial\theta}\left ( 9 - 4sin{^2}\theta -  4sin\theta\right )= -8sin\theta cos\theta - 4cos\theta = 0 \Rightarrow 2sin\theta cos\theta = -cos\theta \Rightarrow sin\theta = -1/2</math></center>


De esta forma, se tienen los puntos de inflexión: <math> 7π/6 </math> y <math> 11π/6 </math>. 
Puesto que <math> cos\theta = 0 </math>, puede haber un máximo o mínimo tanto en <math> 3π/2 </math> como en <math> π/2 </math>.


Los máximos se encuentran en <math> 7π/6 </math> y <math> 11π/6 </math>, mientras que el mínimo se da en <math> π/2 </math>.


A continuación, se proporciona una gráfica que demuestra la presión en el obstáculo circular.

[[Archivo:ecuacionbernuilli1.jpg|thumb|right|Representación de como cambia la presión en el borde del obstáculo al aumentar θ]]

{{matlab|codigo=
tt=linspace(0,2*pi,100)
p=10-(-2.*sin(tt)-1).^2
plot(tt,p)
}}

Los puntos del obstáculo donde se alcanzan máximos de presión son opuestos a los puntos donde el módulo de la velocidad es mínimo y viceversa. 
[[Archivo:ecuacionbernuilli2.jpg|thumb|right|Representación de como cambia la presión en el borde del obstáculo al aumentar θ]]

{{matlab|codigo=
f=@(x) -2*sin(x) - 1;
x=linspace(0,2*pi,1000);
y=f(x);
a=1;

%Se cambian de signo los valores negativos

for i=y
    if i<0
        i=-i;
    end
    y(a)=i;
    a=a+1;
end
hold on
plot(x,y);
tt=linspace(0,2*pi,100);
p=10-(2.*sin(tt+pi)-1).^2;
plot(tt,p)
xlabel('ángulo')
ylabel('módulo de la velocidad y presión')
hold off
}}

==Velocidades y presiones==
Tras hallar gráficamente el campo de velocidades y el campo de presiones, se puede observar que hay una relación entre presión y velocidad. 

Esto se puede deducir a partir del Trinomio de Bernoulli (que se cumple ya que es un fluido incompresible), en el que vemos que cuanto menor es su presión, mayor es el módulo de la velocidad, debido a que el Trinomio ha de ser constante. 


A continuación se expone el Trinomio de Bernoulli y se explican sus componentes;

<center><math> B=z+ \frac{p}{γ} + \frac{v{^2}}{2g} </math></center>

B = valor del trinomio de Bernoulli

p = presión 

z = altura

v = velocidad que lleva el fluido en un determinado punto

g = valor de la gravedad

γ = valor del peso específico

==Ecuación de Navier-Stokes==
La ecuación de Navier-Stokes es un conjunto de ecuaciones en derivadas parciales, que describen el movimiento de un fluido en función de su campo de velocidades y presiones.

Partiendo de la ecuación de Bernoulli, se busca comprobar que <math> (\vec{u},p) </math> satisface la ecuación de Navier-Stokes estacionaria, dada por la fórmula:

<center><math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math></center>


Debido a que los fluidos incompresibles ocupan siempre el mismo volumen, la viscosidad será nula, por lo que <math> \mu =0 </math>. Al igual que anteriormente, se considera que la densidad <math> d = 2 </math>.

<center><math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math></center>


Para resolver la ecuación establecida, se verifica la siguiente identidad vectorial:

<center><math> (\vec{u}·\nabla)\vec{u} = \frac{1}{2} \nabla \| \vec u  \|{^2} - \vec u × \nabla × \vec u </math></center>

Donde <math> \vec u × \nabla × \vec u = 0 </math> debido a que el rotacional es nulo, y al multiplicarlo vectorialmente por <math> \vec u </math> sigue siendo nulo. Por ello, se tiene que:

<center><math> (\vec{u}·\nabla)\vec{u} = \frac{1}{2} \nabla \| \vec u  \|{^2} </math></center>

<center><math> (\vec{u}·\nabla)\vec{u} = \frac{1}{2} \nabla (4 sin {^2} \theta + 4 sin \theta + 1) </math></center>

<center><math> (\vec{u}·\nabla)\vec{u} = \frac{1}{2 \rho}(8sin\theta cos\theta +4cos\theta)\vec{e}_{\theta} </math></center>

De la misma forma, se calcula el gradiente de la ecuación de Bernoulli establecida en el apartado anterior: <math> p = 9 - 4sin{^2}\theta -  4sin\theta </math>

<center><math> \nabla p= \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{1}{\rho} \frac{\partial p}{\partial \theta} \vec{e}_{\theta} </math></center>

<center><math> \nabla p = \frac{1}{\rho}(-8sin\theta cos\theta -4cos\theta)\vec{e}_{\theta} = -\frac{1}{\rho}(8sin\theta cos\theta +4cos\theta)\vec{e}_{\theta} </math></center>


Volviendo a la ecuación inicial, se verifica que <math> (\vec{u},p) </math> satisface la ecuación de Navier-Stokes:

<center><math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math></center>

<center><math> \frac{2}{2 \rho}(8sin\theta cos\theta +4cos\theta)\vec{e}_{\theta} -\frac{1}{\rho}(8sin\theta cos\theta +4cos\theta)\vec{e}_{\theta} = 0 </math></center>

==Paradoja de D'Alembert==

En primer lugar, <math> p · \vec {n} </math> es lo que se consideraría como la fuerza ejercida por el fluido en los diferentes puntos de la frontera, al ver esto, resulta que el sumatorio de todas las proyecciones de las fuerzas sobre la dirección <math> \vec {i} </math> da como respuesta que la resultante es nula. En particular, sobre la dirección horizontal del obstáculo el fluido no realiza ningún empuje.

Como ya se ha calculado en el apartado 6, el valor de <math> \ p = 10 - \left\|u\right \|^2 = 9 - 4sin{^2}\theta -  4sin\theta</math>


Para calcular el vector normal, <math> \vec {n} = -\vec e_\rho </math>, se utilizará la siguiente matriz de cambio de base de cilíndricas a cartesianas y en este caso será la matriz traspuesta: 


\begin{bmatrix}
\cos(\theta) & \sin(\theta) & 0 \\
-\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix} y con esto, se calculará las coordenadas del vector normal: 


\begin{bmatrix}
\cos(\theta) & \sin(\theta) & 0 \\
-\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix} \begin{bmatrix}
1\\
0\\
0
\end{bmatrix} = \begin{bmatrix}
\cos(\theta) \\
-\sin(\theta)\\
0 
\end{bmatrix}  
<math> = cos\theta \vec i - sin\theta \vec j </math>


Tenemos; <math> \vec {n} = cos\theta \vec i - sin\theta \vec j  </math> es decir, que el producto <math> \vec {n}·\vec {i} = -\cosθ </math>

Los límites de integración serán: <math> \left \{ (\rho,\theta)=(1,\theta) \right \} \theta\epsilon[0,2π) </math>

<math> \int_{0}^{2\pi} p·\vec {n}·\vec {i} = \int_{0}^{2\pi} p (-cosθ dθ) =\int_{0}^{2\pi} {(9 - 4sin{^2}\theta -  4sin\theta)(-cosθ)} dθ= \ [-9sen\theta]_{0}^{2\pi} + [\frac{4sen^3\theta}{3}]_{0}^{2\pi} +  [2sen^2\theta]_{0}^{2\pi} = 0 </math> 

Con ello sabemos que el fluido no realiza ningún empuje sobre el obstáculo en la dirección horizontal.

==Curvas de nivel de la presión==
Analizando las curvas de nivel de la presión, se puede deducir que la zona donde la presión es mínima es cuando <math> \theta = π/2 </math>, mientras que las zonas de presión máxima ocurren en <math> \theta = 7π/6 </math> y <math> \theta = 11π/6 </math>.

Se puede observar también, que la presión es menor en la zona superior que la que se da en la zona inferior de la gráfica.

[[Archivo:ULTIMOcorrecto.jpg|thumb|right|Curvas de nivel de la presión]]
{{matlab|codigo=

tt=0:pi/30:2*pi;
rho=1:0.1:2;
[R,T]=meshgrid(rho,tt);
X=R.*cos(T);
Y=R.*sin(T);
p=((1-(1./R.^2)).*(cos(T).^2)+(1+(1./R.^2)).*(sin(T).^2)+(sin(T)./R)).^2+((1-(1./R.^2)).*(cos(T).*sin(T))-(1+(1./R.^2)).*(sin(T).*cos(T))-(cos(T)./R)).^2;
contour(X,Y,p,500)
axis equal
xlabel('EjeX')
ylabel('EjeY')
}}

[[Categoría:Teoría de Campos]] [[Categoría:TC23/24|2023-24]]