Fluido Alrededor de un Obstáculo Circular (Grupo 45)

Trabajo realizado por estudiantes
Título	Fluido Alrededor de un Obstáculo Circular (Grupo 45)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Jose Antonio Martín-Caro Xavier Grimalt Roig Uriel Hidalgo Marcos Emilio Tavío
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El objetivo de este artículo se estudiará el comportamiento de un fluido alrededor de un sólido circular.

Un fluido es una sustancia que puede deformarse continuamente bajo la aplicación de una fuerza de cizallamiento (es decir, una fuerza que actúa paralela a una superficie) sin mostrar resistencia permanente. A nivel físico, los fluidos pueden ser líquidos y gases, ya que ninguno de los dos puede conservar una forma estable. La diferencia entre ellos es que los primeros toman la forma del recipiente donde están, mientras que los segundos tienen tan poca unión entre sus partículas que pueden comprimirse y no tienen ni forma ni volumen propios.

1 .-Superficie Mallada

Se comienza realizando un mallado que describe los puntos interiores de la región ocupada por el fluido. Para llevar a cabo la representación de esta región se emplean coordenadas cilíndricas, definidas en el intervalo radial 1 ≤ r ≤ 5. que posteriormente se transforman a coordenadas cartesianas. Tras esta transformación, el dominio queda incluido en: (x,y) ∈ [−4,4] × [−4,4].

Mediante el siguiente código elaborado en Matlab, se podrá visualizar la superficie de trabajo:

Representación de Superficie Mallada

clear; clc;

%Mallado
r = linspace(1,5,60); %Radios entre 1 y 5
theta = linspace(0,2*pi,80); %Ángulos

% Mallado en coordenadas polares
[R,TH] = meshgrid(r,theta);

%Conversión a coordenadas cartesianas
X = R .* cos(TH);
Y = R .* sin(TH);
Z = zeros(size(R));
% Dibujar el mallado
figure; hold on;
mesh(X,Y,Z);

%Dibujar el círculo unidad
plot(1*cos(theta), 1*sin(theta), 'k', 'LineWidth', 2);

%Vista
axis equal;
axis([-4 4 -4 4]);
view(2); % Vista en planta
title('Mallado del Fluido (Región Exterior al Círculo Unidad)');
xlabel('X');
ylabel('Y');
hold off;

2 .-Función Potencial y Campo de Velocidades del Fluido.

Una vez delimitada la zona que se va a estudiar, se procederá a examinar cómo se comporta el fluido. Para ello, será necesario determinar el campo de velocidades. El dominio del campo escalar estará definido por la función potencial siendo este [math] D\subset \mathbb{R}^{3} [/math] . La función potencial con la que se trabajará para conocer su comportamiento será: [math] \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) [/math]

2.1 .-Representación de la Función Potencial

En primer lugar, se presentará la representación de la función potencial. Todo esto se llevará a cabo utilizando funciones implementadas con la ayuda del programa Matlab:

Lineas de Nivel de la Función Potencial

clear;clc; 
I=linspace(1,5,30);   
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);  
hold on
X=rho.*cos(theta);      
Y=rho.*sin(theta);
%Se introduce la función potencial
f=@(x,y)(x+(1./x)).*cos(y)-y;
%Se generan las líneas de nivel de la función
contour(X,Y,f(rho,theta),50)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
view(2);
axis([-5,5,-5,5]);
colorbar;
title('Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
hold off

2.2 .-Representación del Campo de Velocidades

Sabiendo que la velocidad de las partículas de este fluido es el gradiente de la función potencial tal que: [math] \nabla \varphi=\vec u= \frac{\partial \vec{u}}{\partial \rho }\vec e_\rho + \frac{1}{\rho }\frac{\partial \vec{u}}{\partial \theta }\vec e_\theta+\frac{\partial \vec{u}}{\partial z}\vec e_z = \left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta)\vec e_\rho - [\left ( 1+\frac{1}{\rho^2} \right )\cdot\sin(\theta)-\frac{1}{\rho}]\vec e_\theta [/math]

A continuación, se mostrará una representación de los vectores del campo de velocidades, los cuales son ortogonales a las curvas de nivel de la función potencial. Todo ello se realizará mediante funciones programadas en Matlab. Para poder trabajar en este programa, el gradiente se expresa en coordenadas cartesianas, lo cual se obtiene aplicando la siguiente matriz de cambio de base: [math] \begin{pmatrix} cos(\theta) & - sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math]

Velocidad de las Partículas

clear;clc;
I=linspace(1,5,30);   
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);  
X=rho.*cos(theta);      
Y=rho.*sin(theta);
f=@(x,y)((x+(1./x)).*cos(y))-y;
contour(X,Y,f(rho,theta),50);
hold on
%Gradiente de la función potencial
DX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
quiver(X,Y,DX,DY);
plot(1*cos(J),1*sin(J),'k','lineWidth',1);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de la Función Potencial');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Si se ampliara la imagen en cualquier punto del campo de velocidades, se vería que los vectores que representan la velocidad de las partículas del fluido son siempre perpendiculares a las curvas de nivel de la función potencial durante todo su recorrido alrededor del obstáculo circular. Esto puede apreciarse en la siguiente imagen:

centro

3 .-Rotacional y Divergencia

Es posible verificar que el campo de velocidades [math] \nabla \varphi=\vec {u} [/math] es ortogonal a las curvas de nivel examinando el producto escalar mediante el vector normal [math] \vec n [/math]. Si [math] \vec n [/math] es un vector normal a los puntos del obstáculo entonces [math] \vec u \cdot \vec n = 0 [/math] , esto implica que el campo de velocidades es ortogonal a dichas curvas. Al expresar el problema en coordenadas cilíndricas, el vector normal se obtiene a partir del gradiente [math]\nabla \phi = \frac{\partial \phi}{\partial \rho}\,\vec{e}_{\rho}[/math]
Suponiendo el gradiente en ρ=1 y multiplicándolo por el vector ortogonal anteriormente nombrado, se puede comprobar mediante un breve cálculo la condición de ortogonalidad:
[math](\vec u(\rho,\theta,z)=\left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta)\vec e_\rho - [\left ( 1+\frac{1}{\rho^2} \right )\cdot\sin(\theta)-\frac{1}{\rho}]\vec e_\theta) \cdot (-\vec e_\rho) = 0 [/math]

Si calculamos el valor de [math] \vec u [/math] en la frontera del obstáculo se obtendría: [math]\vec u(1,\theta,z)=\nabla \varphi(1,\theta,z)= (- 2\cdot\sin(\theta) - 1 )\vec e_\theta [/math]

4 .-Líneas de corriente del campo [math]\vec{u}[/math]

-Las líneas de corriente [math]\vec{u}[/math] son las curvas tangentes a los vectores de velocidad de las partículas del fluido. Para hallarlas simplemente hay que calcular el campo ortogonal a [math]\vec{u}[/math] resolviendo el siguiente producto vectorial: [math]\vec{v}[/math] = [math]\vec{k}[/math] x [math]\vec{u}[/math]:

[math]\vec{v}[/math] = [math]\vec{k}[/math] x [math]\vec{u}[/math] = [math]\begin{vmatrix} e_\rho & e_\theta & e_z \\ 0 & 0 & 1\\ cos(\theta)\cdot (1-\frac{1}{\rho ^2}) & sin(\theta) \cdot (1-\frac{1}{\rho^2}) -\frac{1}{\rho} & 0 \end{vmatrix}=[/math] [math] [(1+\frac{1}{\rho ^2}) \cdot sin(\theta) + \frac{1}{\rho}]e_\rho + [(1-\frac{1}{\rho ^2}) \cdot cos(\theta)]e_\theta[/math].


Como se ha comprobado en el apartado 3.1, el rotacional de [math]\vec{u}[/math] es nulo. Esto quiere decir que [math]\vec{v}[/math] es irrotacional. Asimismo, [math]\vec{v}[/math] tiene un potencial llamado [math]\psi [/math], gradiente del propio [math]\vec{v}[/math], y que será descifrado resolviendo el siguiente sistema de ecuaciones diferenciales (se desprecia [math]V_z [/math] ya que se está trabajando en el plano):
[math]\frac{\partial \psi }{\partial \rho }=V_\rho [/math]

[math]\frac{1}{\rho }\frac{\partial \psi }{\partial \theta }=V_\theta [/math]
Tras despejar [math]\frac{1}{\rho }[/math] e integrar ambas ecuaciones obtenemos la función [math]\psi[/math] = [math] (\rho - \frac{1}{\rho}) \cdot sin(\theta) + Ln(\rho) [/math].
Una vez se obtiene tanto [math]\psi[/math] como [math]\vec{u}[/math], se introducen en MATLAB, así comprobando que efectivamente los vectores que componen [math]\vec{v}[/math] son ortogonales a las línes de corriente de [math]\psi[/math]. Para ello se ha empleado el siguiente código:

Lineas de Corriente

{{matlab|codigo=

5 .-Puntos de Frontera y Remanso

El rotacional y la divergencia de un campo vectorial proporcionan información sobre las características físicas del fluido al que estudia y define.

5.1 .-Frontera S

Por su parte, el rotacional indica la dirección y la velocidad del giro del campo en cada punto. En él se considera el movimiento del fluido y su tendencia a inducir rotación en cada punto. Se calcula el rotacional, tal que:

[math]\nabla\times \vec{u}=\frac{1}{\rho }\begin{vmatrix} \vec e_\rho &\rho\vec e_\theta &\vec e_z \\ \frac{\partial }{\partial \rho }& \frac{\partial }{\partial \theta }& \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}=\frac{1}{\rho }\begin{vmatrix} \vec e_\rho &\rho\vec e_\theta &\vec e_z \\ \frac{\partial }{\partial \rho }& \frac{\partial }{\partial \theta }& \frac{\partial }{\partial z}\\ cos\theta (1-\frac{1}{\rho^2}) & -(1+\frac{1}{\rho^2}) \rho sin\theta + 1 & 0 \end{vmatrix} = \frac{1}{\rho }[ (\vec e_\rho \cdot 0) + (\vec e_\theta \cdot 0) + \vec e_z \cdot (-sin\theta(1-\frac{1}{\rho^2}) + sin\theta(1-\frac{1}{\rho^2}))] = 0 [/math]

Por tanto, se puede concluir que es un campo irrotacional , dado que [math]\boxed { \nabla \times \bar { u } =0 } [/math] ,y, por ende, que las partículas del fluido no giran.

5.2 .-Velocidad Nula

La divergencia de un campo vectorial es una magnitud escalar que compara el flujo entrante y el flujo saliente en una superficie.
[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\}= \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\}= \frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0 [/math]

Como se puede observar también se demuestra que la divergencia es nula, dado que [math] \boxed {\nabla \cdot \vec u = 0} [/math] , lo que indica que el volumen del fluido es constante, es decir el fluido no se expande ni se contrae.

6 .-Ecuación de Bernoulli

Las líneas de corriente de [math]\vec{u}[/math] son las curvas tangentes a los vectores de velocidad de las párticulas del fluido. Para hallarlas simplemente hay que calcular el campo ortogonal a [math]\vec{u}[/math] resolviendo el siguiente producto vectorial: [math]\vec{v}[/math] = [math]\vec{k}[/math] x [math]\vec{u}[/math]:

[math]\vec{v}[/math] = [math]\vec{k}[/math] x [math]\vec{u}[/math] = [math]\begin{vmatrix} e_\rho & e_\theta & e_z \\ 0 & 0 & 1\\ cos(\theta)\cdot (1-\frac{1}{\rho ^2}) & sin(\theta) \cdot (1-\frac{1}{\rho^2}) -\frac{1}{\rho} & 0 \end{vmatrix}=[/math] [math] [(1+\frac{1}{\rho ^2}) \cdot sin(\theta) + \frac{1}{\rho}]e_\rho + [(1-\frac{1}{\rho ^2}) \cdot cos(\theta)]e_\theta[/math].


Como se ha comprobado en el apartado 5, el rotacional de [math]\vec{u}[/math] es nulo. Esto quiere decir que [math]\vec{v}[/math] es irrotacional. Asimismo, [math]\vec{v}[/math] tiene un potencial llamado [math]\psi [/math], gradiente del propio [math]\vec{v}[/math], y que será descifrado resolviendo el siguiente sistema de ecuaciones diferenciales (se desprecia [math]V_z [/math] ya que se está trabajando en el plano):
[math]\frac{\partial \psi }{\partial \rho }=V_\rho [/math]

[math]\frac{1}{\rho }\frac{\partial \psi }{\partial \theta }=V_\theta [/math]
Tras despejar [math]\frac{1}{\rho }[/math] e integrar ambas ecuaciones obtenemos la función [math]\psi[/math] = [math] (\rho - \frac{1}{\rho}) \cdot sin(\theta) + Ln(\rho) [/math].
Una vez se obtiene tanto [math]\psi[/math] como [math]\vec{u}[/math], se introducen en MATLAB, así comprobando que efectivamente los vectores que componen [math]\vec{v}[/math] son ortogonales a las línes de corriente de [math]\psi[/math]. Para ello se ha empleado el siguiente código:

Lineas de Corriente

clear;clc;
I=linspace(1,5,30);   
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);
X=rho.*cos(theta);      
Y=rho.*sin(theta);
%Función PSI
C=@(rho,theta)((sin(theta).*(rho-1./rho))+log(rho));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Grandiente de PSI
DX=((1+(1./(rho.^2))).*cos(theta).*sin(theta))+(cos(theta)./rho)+((-1+(1./(rho.^2))).*cos(theta).*sin(theta));
DY=((sin(theta).^2).*(1+(1./(rho.^2))))+(sin(theta)./rho)+((cos(theta).^2).*(1+(1./(rho.^2))));
quiver(X,Y,DX,DY);
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica Líneas de Corriente');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

A la vez que los vectores de [math]\vec{v}[/math] son ortogonales a las líneas de corriente, los de [math]\vec{u}[/math] deberían ser tangentes a estas, y a su vez perpendiculares a los ya mencionados vectores de [math]\vec{v}[/math]. Para demostrar esta afirmación gráficamente, se ha diseñado un nuevo código que permite observar los ángulos rectos que se forman:

Comparación entre [math] \vec v [/math] y [math] \vec u [/math]

clear;clc;
I=linspace(1,5,30);   
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);
X=rho.*cos(theta);      
Y=rho.*sin(theta);
C=@(x,y)((sin(theta).*(rho-1./rho))+log(rho));
Z=C(I,J);
contour(X,Y,Z,30);
hold on
X=rho.*cos(theta);      
Y=rho.*sin(theta);
%Grandiente de PSI
DX=((1+(1./(rho.^2))).*cos(theta).*sin(theta))+(cos(theta)./rho)+((-1+(1./(rho.^2))).*cos(theta).*sin(theta));
DY=((sin(theta).^2).*(1+(1./(rho.^2))))+(sin(theta)./rho)+((cos(theta).^2).*(1+(1./(rho.^2))));
quiver(X,Y,DX,DY);
%Gradiente de la función potencial
DXX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DYY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
quiver(X,Y,DXX,DYY);
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
view(2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de la Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
hold off

Para una mayor apreciación, de las tangencias que forma [math] \vec u [/math] a las líneas de corriente y de la ortogonalidad formada por [math] \vec v [/math] también con las las líneas de corriente, se ha ampliado la fotografía de la gráfica anterior en un punto cualquiera, dado que se cumple a lo largo de todo el campo. Como se comprueba en la siguiente fotografía:

centro

7 .-Trayectoria de la Partícula

En el borde del obstáculo [math]\rho = 1[/math] lo que quiere decir que [math]\vec{u}(1,\theta) = (-2 \cdot sin(\theta) - 1)[/math].
El módulo de [math]\vec{u}(1,\theta)[/math] es |[math]\vec{u}[/math]| = [math]\sqrt(4 \cdot sin(\theta)^2 + 4 \cdot sin(\theta) + 1)[/math]
Para que |[math]\vec{u}[/math]| alcance su valor máximo, [math]sin(\theta)[/math] tiene que ser igual a 1. Esto se alcanza en el ángulo [math]\frac{\pi}{2}[/math].
Para calcular el punto de remanso se ha de igualar |[math]\vec{u}[/math]| a 0 y despejar el valor de [math]\theta[/math]:

[math]\sqrt(4 \cdot sin(\theta)^2 + 4 \cdot sin(\theta) + 1) = 0[/math]

Finalmente [math]\theta = \frac{11 \pi}{6}[/math] y [math]\frac{7 \pi}{6}[/math] ; existen 2 puntos de remanso.

Para observar estos cálculos visualmente, una vez más se a codificado todo en MATLAB para así observar como en los puntos de remanso efectivamente no se aprecia ningún vector:

Punto de Remanso

clear;clc;
I=linspace(1,5,50);   
J=linspace(0,2*pi,50);
[rho,theta]=meshgrid(I,J);  
X=rho.*cos(theta);      
Y=rho.*sin(theta);
f=@(x,y)(x+(1./x)).*cos(y)-y;
contour(X,Y,f(rho,theta),50);
hold on
DX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
quiver(X,Y,DX,DY);
plot(1*cos(J),1*sin(J),'k','lineWidth',1);
%Punto de remanso previamente calculado
plot(0.866025,-0.5,'r*','LineWidth',3);
plot(-0.866025,-0.5,'r*','LineWidth',3)
axis([-5,5,-5,5]);
colorbar;
title('Punto de Remanso');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Para una mayor claridad se ha ampliado la gráfica anterior, para así conseguir una mayor nitidez de uno de los puntos de remanso:

centro

8 .-Paradoja de D´Alembert

Con el fin de calcular los puntos donde se alcanza mayor y menor presión. Para ello, se supondrá que la densidad es igual a 2 (d=2). Además, se ha de satisfacer la ecuación de Bernouilli, que estipula : [math] \frac{1}{2} d |\vec u|^2 + p = cte [/math].
También, se le asignará a la cte anterior el valor 10 para permitir el cálculo de las presión del fluido.
Una vez aplicados los cambios nombrados anteriormente, queda la siguiente expresión: [math] p=10-|\vec u|{^2} [/math]

El módulo del gradiente es: |[math]\vec{u}[/math]| = [math]\sqrt((\left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta))^2 + (\left ( 1+\frac{1}{\rho^2} \right )\cdot\sin(\theta)-\frac{1}{\rho})^2)[/math]

El módulo del gradiente al cuadrado, es decir, [math] |\vec u|^2 = cos^{2} \theta \cdot (1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}) + (sin^2 \theta + \frac{2sin^2 \theta}{\rho^2} + \frac{2sin^2 \theta}{\rho^4}) + \frac{2}{\rho} \cdot (sin \theta + \frac{sin \theta}{\rho^2}) + \frac{1}{\rho^2} [/math]

Por lo tanto, la presión es: [math] p = 10 - cos^{2} \theta \cdot (1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}) + (sin^2 \theta + \frac{2sin^2 \theta}{\rho^2} + \frac{2sin^2 \theta}{\rho^4}) + \frac{2}{\rho} \cdot (sin \theta + \frac{sin \theta}{\rho^2}) + \frac{1}{\rho^2}. [/math]

Este campo de presiones queda representado en la siguiente gráfica que se obtuvo con la ayuda del siguiente código desarrollado en MATLAB:

Mallado de Presiones

clear;clc;
I=linspace(1,5,50);   
J=linspace(0,2*pi,50);
[rho,theta]=meshgrid(I,J);  
X=rho.*cos(theta);      
Y=rho.*sin(theta);
Z=zeros(size(X));
%Ecuación de la presión
p=10-(((cos(theta).^2).*(1-(2./rho.^2)+(1./rho.^4)))+((sin(theta).^2)+((2.*(sin(theta).^2))./rho.^2)+((sin(theta).^2)./rho.^4))+((2./rho).*(sin(theta)+(sin(theta)./(rho.^2))))+(1./(rho.^2)));
surf(X,Y,p);
hold on 
axis([-5,5,-5,5]);
colorbar;
view(2); 
title('Mallado de Presiones');
xlabel 'X'
ylabel 'Y'
axis equal
hold off

Para finalizar se empleará de nuevo un código de Matlab, en el que se obtiene el punto con mayor presión que es: [math] p_{max}= 18.9938 [/math] uds y el punto con menor presión siendo este: [math] p_{min}=9.0890[/math] uds.
A continuación se muestra el código empleado:

p=10-((cos(theta).^2).*(1-(2./rho.^2)+(1./rho.^4)))+((sin(theta).^2)+((2.*(sin(theta).^2))./rho.^2)+((sin(theta).^2)./rho.^4))+((2./rho).*(sin(theta)+(sin(theta)./(rho.^2))))+(1./(rho.^2));
max(max(p))
min(min(p))

9 .-Reanalisis de los apartados 2,3 y 4

Las llamadas ecuaciones de Navier-Stokes describen matemáticamente el movimiento de un fluido. Estas ecuaciones deberían determinar el movimiento futuro del fluido a partir de su estado inicial. En este apartado,se pretende comprobar que partiendo de la ecuación de Bernouilli, que [math] (\vec{u},p) [/math] satisfacen la ecuación de Navier-Stokes estacionaria, que viene dada por la siguiente expresión:

[math] d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} [/math]

Para este cálculo, se supondrá que µ = 0, es decir, viscosidad nula; y que d(densidad) [math] = [/math] 2:

[math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math]

Aplicando las propiedades teóricas algebraicas se produce la siguiente igualdad:

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

En consecuencia, a que el rotacional es nulo, al multiplicarlo por [math] \vec u [/math] sigue siendo nulo y por lo tanto se comprueba que [math] \vec u × \nabla × \vec u = 0 [/math] . Por tanto sustituyendo [math] \vec u × \nabla × \vec u = 0 [/math] en el paso anterior obtenemos:

[math] (\vec{u}·\nabla)\vec{u} = \frac{1}{2} \nabla |\vec u|{^2} = (\vec{u}·\nabla)\vec{u} = \frac{1}{2} \nabla (4 sin {^2} \theta + 4 sin \theta + 1) = (\vec{u}·\nabla)\vec{u} = \frac{1}{2 \rho}(8sin\theta cos\theta +4cos\theta)\vec{e}_{\theta} [/math]

A continuación se calcula el gradiente de la ecuación de Bernouilli:

[math] p = 9 - 4sin{^2}\theta - 4sin\theta [/math]

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

Retrocediendo hasta el inicio de este apartado, e introduciendo en [math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math] las variables calculadas, se concluye finalmente con que [math] (\vec{u},p) [/math] satisface la ecuación de Navier-Stokes:

[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]

9.1 .-Función Potencial y Campo de Velocidades

Función Potencial

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
hold on
X = rho.*cos(theta);
Y = rho.*sin(theta);
f = @(rho,theta)((rho + 1./rho).*cos(theta) + theta./(4*pi));
contour(X,Y,f(rho,theta),50)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
view(2);
axis([-5,5,-5,5]);
colorbar;
title('Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
hold off

Función Velocidad

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
X = rho.*cos(theta);
Y = rho.*sin(theta);
f = @(rho,theta)((rho + 1./rho).*cos(theta) + theta./(4*pi));
contour(X,Y,f(rho,theta),50)
hold on
dphi_r = (1 - 1./(rho.^2)).*cos(theta);
dphi_t = -(1 + 1./(rho.^2)).*sin(theta) + 1./(4*pi*rho);
DX = dphi_r.*cos(theta) - dphi_t.*sin(theta);
DY = dphi_r.*sin(theta) + dphi_t.*cos(theta);
quiver(X,Y,DX,DY)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de la Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
view(2);
hold off