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 Pedro Comas Payeras
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El objetivo de este artículo es 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:

Figura 1 — 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 definida la región donde se estudiará el comportamiento del fluido, se procede a analizar su dinámica a partir del campo de velocidades. Para ello, se utiliza una función potencial escalar: [math]\phi = \phi(\rho,\theta)[/math], definida en un dominio [math]D \subset \mathbb{R}^2[/math] expresado en coordenadas polares.

En este caso, la función potencial viene dada por: [math]\phi(\rho,\theta) = \left( \rho + \frac{1}{\rho} \right)\cos\theta[/math].

A partir de esta función se obtiene el campo de velocidades aplicando: [math]\vec{u} = \nabla\phi[/math], lo que permite describir la dirección y magnitud del movimiento del fluido, así como estudiar la relación entre las curvas de nivel del potencial y las líneas de corriente.

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

Para estudiar con mayor claridad la naturaleza del flujo, es útil examinar la forma que adopta la función potencial [math]\phi(\rho,\theta)[/math] dentro del dominio considerado. La representación gráfica de esta función permite identificar zonas donde el potencial crece o disminuye con mayor rapidez, así como patrones característicos que influyen en el comportamiento del fluido.

La construcción de estas gráficas se realiza mediante herramientas de visualización numérica, en este caso, MATLAB, que posibilitan generar superficies del potencial. Estas representaciones facilitan la interpretación del campo y sirven como apoyo previo al análisis del gradiente y de las velocidades resultantes.

Figura 2.1 — Curvas de nivel de la función potencial 𝜙 ( 𝜌 , 𝜃 )

clear; clc;

%Mallado en coordenadas polares
R = linspace(1,5,80); %Rho
T = linspace(0,2*pi,180); %Theta
[rho,theta] = meshgrid(R,T);

%Transformación a coordenadas cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Definición de la función potencial
phi = (rho+1./rho).*cos(theta);

%Representación de la función potencial (curvas de nivel)
figure;
contour(X, Y, phi, 60, 'LineWidth', 1.5);
hold on;
plot(1*cos(T), 1*sin(T), 'k', 'LineWidth', 2); %Círculo interior
axis equal;
colorbar;
title('Función Potencial');
xlabel('X');
ylabel('Y');
hold off;

2.2 .-Representación del Campo de Velocidades

Con la función potencial es posible determinar el campo de velocidades que describe el movimiento del fluido. Recordemos que la velocidad se calcula a partir del gradiente de la función potencial: [math]\vec{u} = \nabla\phi[/math].

Para ello se expresan las derivadas en coordenadas polares, donde el operador gradiente viene dado por: [math]\nabla\phi = \frac{\partial\phi}{\partial\rho}\,\vec{e}_\rho + \frac{1}{\rho}\frac{\partial\phi}{\partial\theta}\,\vec{e}_\theta[/math].

• La función potencial es: [math]\phi(\rho,\theta) = \left(\rho + \frac{1}{\rho}\right)\cos\theta[/math].

• Sus derivadas parciales son: [math]\frac{\partial\phi}{\partial\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]        [math]\frac{\partial\phi}{\partial\theta} = -\left(\rho + \frac{1}{\rho}\right)\sin\theta[/math].

Por tanto, las componentes del campo de velocidades en polares son: [math]u_{\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]         [math]u_{\theta} = -\left(1 + \frac{1}{\rho^{2}}\right)\sin\theta[/math].

Para representar el campo en MATLAB es más conveniente trabajar en coordenadas cartesianas. Usamos la transformación:

- [math]u_x = u_\rho\cos\theta - u_\theta\sin\theta[/math],

- [math]u_y = u_\rho\sin\theta + u_\theta\cos\theta[/math].

Tras realizar los cálculos:

- [math]u_x = (1 - \frac{1}{\rho^2})\cos^2\theta + (1 + \frac{1}{\rho^2})\sin^2\theta[/math],

- [math]u_y = -\frac{2}{\rho^2}\sin\theta\cos\theta[/math].

Figura 2.2 - Campo de velocidades alrededor del obstáculo circular

clear; clc;

%Mallado en polares
R = linspace(1,5,40);
T = linspace(0,2*pi,40);
[rho,theta] = meshgrid(R,T);

%Transformación a cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Función potencial
phi = (rho + 1./rho).*cos(theta);

%Componentes de la velocidad en polares
ur = (1 - 1./(rho.^2)).*cos(theta);
ut = -(1 + 1./(rho.^2)).*sin(theta);

%Conversión a componentes cartesianas
DX = ur.*cos(theta) - ut.*sin(theta);
DY = ur.*sin(theta) + ut.*cos(theta);

%Gráfica
figure;
contour(X, Y, phi, 40); hold on;
quiver(X, Y, DX, DY, 'k');
plot(cos(T), sin(T), 'r', 'LineWidth', 1.3); %Obstáculo circular
axis equal; axis([-5 5 -5 5]);
title('Campo de Velocidades');
xlabel('X'); ylabel('Y'); colorbar;
hold off;

Una característica interesante de este flujo es que las líneas de corriente coinciden con las trayectorias que seguirían partículas sin inercia, y estas son siempre perpendiculares a las curvas de nivel de la función potencial. Esto es una propiedad general de los flujos potenciales: el gradiente siempre apunta en la dirección de máxima variación del potencial, mientras que las curvas de nivel representan zonas donde el potencial es constante. Si se hace un zoom en cualquier área del diagrama se aprecia claramente que los vectores del campo de velocidades mantienen esta ortogonalidad en todo el dominio, especialmente alrededor del obstáculo circular donde el cambio de dirección es más brusco.

3 .-Rotacional y Divergencia

El rotacional y la divergencia de un campo vectorial proporcionan información esencial sobre el comportamiento físico del fluido que representan. La divergencia permite identificar si el fluido se comprime o se expande localmente, mientras que el rotacional muestra si las partículas experimentan algún tipo de giro o movimiento

3.1 .-Rotacional nulo

El rotacional muestra la dirección y la intensidad del giro del fluido en cada punto. Para analizar si el flujo induce rotación, se calcula el rotacional del campo de velocidades del siguiente modo:

• El campo de velocidades es: [math]\vec{u}=\nabla\phi=\frac{\partial\phi}{\partial\rho}\vec{e}_\rho+\frac{1}{\rho}\frac{\partial\phi}{\partial\theta}\vec{e}_\theta[/math] y sus componentes son: [math]u_\rho=\left(1-\frac{1}{\rho^2}\right)\cos\theta[/math], [math]u_\theta=-\left(1+\frac{1}{\rho^2}\right)\sin\theta[/math] y [math]u_z=0[/math].

• El rotacional en coordenadas cilíndricas es:[math]\nabla\times\vec{u}=\frac{1}{\rho}\begin{vmatrix}\vec e_\rho & \rho\vec e_\theta & \vec e_z \\[6pt] \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\[6pt] \left(1-\frac{1}{\rho^2}\right)\cos\theta & -\rho\left(1+\frac{1}{\rho^2}\right)\sin\theta & 0 \end{vmatrix}=\frac{1}{\rho}\left[\vec e_\rho\cdot 0+\vec e_\theta\cdot 0+\vec e_z\left(\frac{\partial}{\partial\rho}\left(-\rho\left(1+\frac{1}{\rho^2}\right)\sin\theta\right)-\frac{\partial}{\partial\theta}\left(\left(1-\frac{1}{\rho^2}\right)\cos\theta\right)\right)\right]=0[/math]

Por lo tanto [math]\boxed { \nabla \times \bar { u } =0 } [/math] y esto significa que el campo es irrotacional por lo que las particulas del fluido no rotan sobre sí mismas al moverse por el flujo.

3.2 .-Divergencia nula

La divergencia de un campo vectorial es la cantidad que mide la diferencia entre el flujo que entra y el flujo que sale del volumen.

• [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(\left(\rho-\frac{1}{\rho}\right)\cos\theta\right)+\frac{\partial}{\partial\theta}\left(-\left(1+\frac{1}{\rho^2}\right)\sin\theta\right)+0\right\}=\frac{1}{\rho}\left\{\left(1+\frac{1}{\rho^2}\right)\cos\theta-\left(1+\frac{1}{\rho^2}\right)\cos\theta\right\}=0[/math]

Como se observa, la divergencia resulta ser nula, es decir, [math] \boxed {\nabla \cdot \vec u = 0} [/math]. Esto implica que el volumen del fluido permanece constante: el flujo no se expande ni se contrae, por lo que el movimiento es incompresible.

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

Las líneas de corriente de [math]\vec{u}[/math] son las curvas tangentes a los vectores velocidad de las partículas del fluido. Por lo que se comenzará calculando el campo [math]\vec{v}[/math], que en cada punto es ortogonal a [math]\boxed{\vec{u} : \vec{v} = \vec{k} \times \vec{u}}[/math] :

[math]\vec{v} =\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 & 0 \end{vmatrix} =(1+\tfrac{1}{\rho^2})\sin\theta\,\vec e_\rho +\left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta [/math]

En el apartado 3.1 se ha comprobado que , el rotacional de [math]\vec{u}[/math] es nulo. Este resultado implica que [math]\vec{v}[/math] es irrotacional. Asimismo, [math]\vec{v}[/math] posee un potencial [math]\psi[/math], cuyo gradiente coincide con el propio campo.

Para determinar dicho potencial debe resolverse el sistema (se desprecia [math]V_z[/math]):

• [math]\frac{\partial \psi }{\partial \rho }=V_\rho [/math]

• [math]\frac{1}{\rho }\frac{\partial \psi}{\partial \theta}=V_\theta [/math]

Al haber despejado e integrado ambas ecuaciones [math]\frac{1}{\rho }[/math] se obtiene la función:[math] \boxed{ \psi=\sin\theta\left(\rho-\frac{1}{\rho}\right) } [/math]

Una vez calculado [math]\psi[/math] del campo [math]\vec{u}[/math], con la ayuda de MATLAB verificamos que los vectores de [math]\vec{v}[/math] son ortogonales a las líneas de corriente de [math]\psi[/math]. Se puede ver gráficamente, gracias a el siguiente codigo de matlab:

Líneas 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)((rho-1./rho).*sin(theta));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Gradiente de PSI (campo v)
DX= ((1+1./(rho.^2)).*sin(theta));
DY= ((1-1./(rho.^2)).*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 Líneas de Corriente');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Aun cuando los vectores de [math]\vec{v}[/math] se presentan como ortogonales respecto a las líneas de corriente, los vectores de [math]\vec{u}[/math] deben ser tangentes a dichas curvas, manteniendo, a su vez, una perpendicularidad rigurosa con los vectores anteriormente citados de [math]\vec{v}[/math].Para poder demostrar esta afirmación graficamente, nos apoyamos en un codgio MATLAB.

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);
%Función PSI
C=@(rho,theta)((rho-1./rho).*sin(theta));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Gradiente de PSI (campo v)
DX= ((1+1./(rho.^2)).*sin(theta));
DY= ((1-1./(rho.^2)).*cos(theta));
quiver(X,Y,DX,DY);
%Gradiente de la función potencial phi (campo u)
DXX= (1-1./(rho.^2)).*cos(theta);
DYY= -(1+1./(rho.^2)).*sin(theta);
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

centro

Con el objetivo de una apreciación más detallada de la tangencialidad que forma [math]\vec{u}[/math] con las líneas de corriente y de ortogonalidad formada por [math]\vec{v}[/math] también con las las líneas de corriente, se ha ampliado la representación gráfica en un punto particular. Cabe destacar que esta relación se mantiene de manera uniforme en todo el campo, como se comprueba en la siguiente gráfica:

5 .-Puntos de Frontera S y Remanso

En este apartado se analizará la frontera [math]S[/math] del obstáculo circular de radio unidad.Se determinarán los puntos donde el módulo de la velocidad es máximo y mínimo.
Asimismo, se identificarán los puntos en los que la velocidad es nula, conocidos como puntos de remanso. Finalmente, se representarán gráficamente los puntos de remanso sobre el borde del obstáculo para comprobar su ubicación en el flujo.

5.1 .-Puntos de la frontera S

Sobre la frontera se cumple: [math]\rho = 1 . [/math]

Las componentes del campo de velocidades (ya calculadas previamente) eran:

[math] u_\rho = \left( 1 - \frac{1}{\rho^2} \right)\cos\theta, \qquad u_\theta = -\left( 1 + \frac{1}{\rho^2} \right)\sin\theta. [/math]

En la frontera [math] \rho = 1 . [/math]

[math] u_{\rho}\big|_{\rho=1} = 0, \qquad u_{\theta}\big|_{\rho=1} = -2\sin\theta. [/math]

La rapidez sobre 𝑆 es:

[math] |\vec{u}| = 2\,|\sin\theta|. [/math]

• Máxima velocidad en 𝑆:

Ocurre cuando:

[math] |\sin\theta| = 1 \quad\Rightarrow\quad \theta = \frac{\pi}{2},\ \frac{3\pi}{2}. [/math]       En cartesianas: [math] (\,x,y\,)_{\text{max}} = (0,\,1),(0,\,-1). [/math]

• Mínima velocidad en 𝑆:

Ocurre cuando:

[math] |\sin\theta| = 0 \quad\Rightarrow\quad \theta = 0,\ \pi. [/math]           En cartesianas: [math] (\,x,y\,)_{\text{min}} = (1,\,0),(-1,\,0). [/math]

También se podría haber derivado la expresión de la rapidez e igualado a cero para obtener los extremos, pero ambos procedimientos —derivar o razonarlo directamente a partir de 2∣sin𝜃∣— conducen exactamente al mismo resultado.

5.2 .-Puntos de remanso

Los puntos de remanso sobre 𝑆 se obtienen imponiendo: [math] u_{\rho}\big|_{\rho=1}=0,\quad u_{\theta}\big|_{\rho=1}=0 \;\Longrightarrow\; -2\sin\theta=0 \;\Longrightarrow\; \sin\theta=0. [/math]
De aquí: [math] \theta = 0,\quad \theta=\pi. [/math]

En coordenadas cartesianas:

• [math]\theta = 0 \Rightarrow (x,y) = (1,0)[/math]
  
• [math]\theta = \pi \Rightarrow (x,y) = (-1,0)[/math]

Observaciones físicas:

En estos puntos la velocidad del fluido es nula respecto al sólido (puntos de estancamiento en la superficie del cilindro).
Según la ecuación de Bernoulli (en flujo potencial incompresible), un punto de remanso corresponde localmente a un máximo de presión estática.

Se pueden observar mejor estos puntos calculados visualmente a través de MATLAB:

Figura 5 - Visualización de los puntos de remanso y de velocidades máxima y mínima

clear; clc;

%Mallado en polares
R = linspace(1,5,40);
T = linspace(0,2*pi,40);
[rho,theta] = meshgrid(R,T);

%Transformación a cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Función potencial
phi = (rho + 1./rho).*cos(theta);

%Componentes de la velocidad en polares
ur = (1 - 1./(rho.^2)).*cos(theta);
ut = -(1 + 1./(rho.^2)).*sin(theta);

%Conversión a componentes cartesianas
DX = ur.*cos(theta) - ut.*sin(theta);
DY = ur.*sin(theta) + ut.*cos(theta);

%%PUNTOS DE LA FRONTERA S
%Velocidad máxima: (0,1) y (0,-1)
p_max = [0 1;
0 -1];

%Velocidad mínima: (1,0) y (-1,0)
p_min = [1 0;
-1 0];

%Gráfica
figure;
contour(X, Y, phi, 40); hold on;
quiver(X, Y, DX, DY, 'k');
plot(cos(T), sin(T), 'r', 'LineWidth', 1.3); %Obstáculo circular

%Dibujar puntos
plot(p_max(:,1), p_max(:,2), 'bo', 'MarkerFaceColor','b', 'MarkerSize',8);
plot(p_min(:,1), p_min(:,2), 'mo', 'MarkerFaceColor','m', 'MarkerSize',8);

legend('Líneas de potencial','Campo de velocidades',...
'Obstáculo','Vel. máxima','Vel. mínima');

axis equal; axis([-5 5 -5 5]);
title('Campo de Velocidades con Puntos Característicos');
xlabel('X'); ylabel('Y'); colorbar;
hold off;

6 .-Ecuación de Bernoulli

Con el fin de determinar los puntos donde el fluido alcanza mayor y menor presión, se considera una densidad constante igual a 2 (d=2). Además, debe cumplirse la ecuación de Bernoulli, de modo que [math]\tfrac{1}{2}d|\vec{u}|^2 + p = \text{cte}[/math] y para simplificar los cálculos, se asigna a dicha constante el valor 10.

Sustituyendo y despejando la presión, obtenemos finalmente: [math]p = 10 - |\vec{u}|^2[/math].

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

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

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

Este campo de presiones se representa en la siguiente gráfica, obtenida mediante el código desarrollado en MATLAB.

Figura 7 - 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

Como los vectores [math]\vec{v}[/math] son perpendiculares a las líneas de corriente, los vectores [math]\vec{u}[/math] deben ser tangentes a ellas y, por tanto, ortogonales a [math]\vec{v}[/math]. Para visualizar esta propiedad se ha generado un código que muestra claramente dichos ángulos rectos.

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

Si fuéramos una partícula del fluido, seguiríamos la trayectoria de una línea de corriente. Para analizar cómo cambiarían nuestra velocidad y presión al rodear el obstáculo, partimos del potencial:

[math] \varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta. [/math]

Las componentes del campo de velocidades se obtienen derivando:

[math] u_\rho=\frac{\partial\varphi}{\partial\rho} = \left(1-\frac{1}{\rho^{2}}\right)\cos\theta, [/math] [math] u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta} = -\left(1+\frac{1}{\rho^{2}}\right)\sin\theta. [/math]

Al considerarlo sobre el borde del obstáculo [math]\rho=1[/math], queda:

[math] u_\rho(1,\theta)=0,\qquad u_\theta(1,\theta)=-2\sin\theta. [/math]

Por lo tanto, el módulo de la velocidad en la superficie del cuerpo viene dado por:

[math] |\vec{u}(1,\theta)| = | -2\sin\theta | = 2|\sin\theta|. [/math]

Esto significa que la velocidad es máxima cuando [math]|\sin\theta|=1[/math], es decir, en las posiciones laterales del obstáculo:

[math] \theta=\frac{\pi}{2},\qquad \theta=\frac{3\pi}{2}. [/math]

En estas zonas el fluido se acelera debido a la geometría del obstáculo. Según el principio de Bernoulli, donde la velocidad aumenta, la presión disminuye, por lo que la presión es mínima en los lados del cilindro.

Por el contrario, cuando [math]\sin\theta=0[/math], la velocidad es mínima:

[math] |\vec{u}(1,\theta)| = 0, [/math]

lo cual ocurre en:

[math] \theta=0,\qquad \theta=\pi. [/math]

En estas zonas el fluido prácticamente se detiene al encontrarse de frente con el obstáculo, por lo que la presión es máxima.

Para representar visualmente estas variaciones de velocidad alrededor del obstáculo, se utiliza el siguiente código en MATLAB:

8 .-Paradoja de D´Alembert

Como p [math]\vec{n}[/math] es la fuerza que ejerce el fluido en cada punto de la frontera, al sumar la proyección de todas estas fuerzas sobre la dirección [math]\vec{i}[/math] la resultante es nula. Por lo que, el fluido no realizara ninguna fuerza sobre el obstáculo en la dirección horizontal.

Para poder demostrarlo, hay que resolver la siguiente integral:

[math] \int_{0}^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta [/math]

[math] p = 1 - 4 \sin^2\theta [/math]

[math] \vec{n} = - \vec{e}_\rho [/math]

[math]\vec {i} = 1\cdot\vec {i}[/math] ; para hallar el valor de este vector en cilíndricas se usará la matriz de cambio de base previamente utilizada pero traspuesta:

y el vector unitario horizontal en cilíndricas se obtiene mediante la matriz de cambio de base transpuesta:

[math] \begin{pmatrix} cos(\theta) & sin(\theta) & 0 \\ -sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = cos(\theta)\vec {e}_\rho-sin(\theta)\vec {e}_\theta[/math]

Finalmente se llega a [math]\vec {n}\cdot\vec {i}=-cos(\theta)[/math] pudiéndose así realizar la integral:

[math] \int_0^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta = \int_0^{2\pi} (1 - 4 \sin^2\theta) \cdot (-\cos\theta) , d\theta = \int_0^{2\pi} (-\cos\theta + 4 \cos\theta \sin^2\theta) , d\theta [/math]

Resolviendo término a término:

[math] \int_0^{2\pi} -\cos\theta , d\theta = [-\sin\theta]_0^{2\pi} = 0 [/math]
[math] \int_0^{2\pi} 4 \cos\theta \sin^2\theta , d\theta = [\frac{4}{3} \sin^3\theta]_0^{2\pi} = 0 [/math]

Por lo tanto, la integral total es nula:

[math] \int_{0}^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta = 0 [/math]

9 .-Reanálisis 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

9.2 .-Rotacional y Divergencia

-Rotacional- 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 & \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} [/math]

Sustituyendo las componentes del campo obtenidas a partir del nuevo potencial:

[math] u_\rho=(1-\frac{1}{\rho^2})\cos\theta,\qquad u_\theta=-(1+\frac{1}{\rho^2})\sin\theta+\frac{1}{4\pi\rho},\qquad u_z=0, [/math]

tenemos:

[math] \nabla\times \vec{u} = \frac{1}{\rho} \left[ \vec{e}_\rho(0) + \vec{e}_\theta(0) + \vec{e}_z \left( \frac{\partial}{\partial\rho}(\rho u_\theta) - \frac{\partial u_\rho}{\partial\theta} \right) \right]. [/math]

Calculamos los términos:

[math] \frac{\partial}{\partial\rho}(\rho u_\theta) = \frac{\partial}{\partial\rho} \left[ \rho\left( -(1+\tfrac{1}{\rho^2})\sin\theta+\frac{1}{4\pi\rho} \right) \right] = -(1-\tfrac{1}{\rho^2})\sin\theta-\frac{1}{4\pi\rho^2}, [/math] [math] \frac{\partial u_\rho}{\partial\theta} = \frac{\partial}{\partial\theta} \left[ (1-\tfrac{1}{\rho^2})\cos\theta \right] = -(1-\tfrac{1}{\rho^2})\sin\theta. [/math]

Entonces:

[math] \frac{\partial}{\partial\rho}(\rho u_\theta) - \frac{\partial u_\rho}{\partial\theta} = -(1-\tfrac{1}{\rho^2})\sin\theta-\frac{1}{4\pi\rho^2} + (1-\tfrac{1}{\rho^2})\sin\theta = -\frac{1}{4\pi\rho^2}. [/math]

Dado que este término está multiplicado por [math]\vec{e}_z[/math] y además dividido por [math]\rho[/math], queda:

[math] \nabla\times\vec{u} = \frac{1}{\rho}\left(0+0-\frac{1}{4\pi\rho^2}\vec{e}_z\right) = -\frac{1}{4\pi\rho^3}\vec{e}_z. [/math]

Por tanto, el campo no tiene rotación salvo por una pequeña contribución del término angular, y se hace despreciable lejos del cilindro. En zonas de estudio donde [math]\rho \gg 1[/math], se considera aproximadamente irrotacional.

Divergencia La divergencia de un campo vectorial es una magnitud escalar que compara el flujo entrante y el flujo saliente en una superficie.

En coordenadas cilíndricas:

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

Como trabajamos en 2D, [math]u_z=0[/math] y su derivada también.

Sustituimos:

[math] \rho u_\rho = \rho(1-\tfrac{1}{\rho^2})\cos\theta = (\rho-\tfrac{1}{\rho})\cos\theta, [/math] [math] \frac{\partial}{\partial\rho}(\rho u_\rho) = (1+\tfrac{1}{\rho^2})\cos\theta, [/math] [math] \frac{\partial u_\theta}{\partial\theta} = -\,(1+\tfrac{1}{\rho^2})\cos\theta. [/math]

Por tanto:

[math] \nabla\cdot\vec{u} = \frac{1}{\rho} \left[ (1+\tfrac{1}{\rho^2})\cos\theta - (1+\tfrac{1}{\rho^2})\cos\theta \right] = 0. [/math]

Como se puede observar, también se demuestra que la divergencia es nula, dado que [math]\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: se cumple la condición de incompresibilidad.

9.3 .-Lineas de Corriente de campo u

Las líneas de corriente de [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}=\vec{k}\times \vec{u}[/math]:

[math] \vec{v}= \begin{vmatrix} \vec{e}_\rho & \vec{e}_\theta & \vec{e}_z\\ 0 & 0 & 1\\ \frac{\partial \varphi}{\partial\rho} & \frac{1}{\rho}\frac{\partial\varphi}{\partial\theta} & 0 \end{vmatrix} = -\left(\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}\right)\vec{e}_\rho + \left(\frac{\partial\varphi}{\partial\rho}\right)\vec{e}_\theta. [/math]

La función potencial dada es:

[math] \varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}. [/math]

Calculamos sus derivadas parciales:

[math] \frac{\partial\varphi}{\partial\rho} = \left(1-\frac{1}{\rho^{2}}\right)\cos\theta, [/math] [math] \frac{1}{\rho}\frac{\partial\varphi}{\partial\theta} = -(1+\frac{1}{\rho^{2}})\sin\theta + \frac{1}{4\pi\rho}. [/math]

Sustituyendo en la expresión de [math]\vec{v}[/math] obtenemos:

[math] V_\rho = (1+\frac{1}{\rho^{2}})\sin\theta - \frac{1}{4\pi\rho}, [/math] [math] V_\theta = (1-\frac{1}{\rho^{2}})\cos\theta. [/math]

Como se ha comprobado en el apartado anterior, el rotacional de [math]\vec{v}[/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(\rho,\theta)= (\rho-\frac{1}{\rho})\sin\theta - \frac{1}{4\pi}\ln(\rho). [/math]

Una vez se obtiene tanto [math]\psi[/math] como [math]\vec{v}[/math], se introducen en MATLAB, así comprobando que efectivamente los vectores que componen [math]\vec{v}[/math] son ortogonales a las líneas 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);
C = @(rho,theta)( sin(theta).*(rho - 1./rho) + (theta/(4*pi)).*log(rho) );
Z = C(rho,theta);
contour(X,Y,Z,30)
hold on
dpsi_r = sin(theta).*(1 + 1./(rho.^2)) + (theta./(4*pi)).*(1./rho);
dpsi_t = cos(theta).*(rho - 1./rho) + (1/(4*pi))*log(rho);
DX = dpsi_r.*cos(theta) - (dpsi_t./rho).*sin(theta);
DY = dpsi_r.*sin(theta) + (dpsi_t./rho).*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 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:

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
C = @(rho,theta)( sin(theta).*(rho - 1./rho) + (theta/(4*pi)).*log(rho) );
Z = C(rho,theta);
X = rho.*cos(theta);
Y = rho.*sin(theta);
contour(X,Y,Z,30)
hold on
dpsi_r = sin(theta).*(1 + 1./(rho.^2)) + (theta./(4*pi))./(rho);
dpsi_t = cos(theta).*(rho - 1./rho) + (1/(4*pi))*log(rho);
DX = dpsi_r.*cos(theta) - (dpsi_t./rho).*sin(theta);
DY = dpsi_r.*sin(theta) + (dpsi_t./rho).*cos(theta);
quiver(X,Y,DX,DY)
dphi_r = (1 - 1./(rho.^2)).*cos(theta);
dphi_t = -(1 + 1./(rho.^2)).*sin(theta) + 1./(4*pi*rho);
DXX = dphi_r.*cos(theta) - dphi_t.*sin(theta);
DYY = dphi_r.*sin(theta) + dphi_t.*cos(theta);
quiver(X,Y,DXX,DYY)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
title('Comparación entre v y u');
xlabel('X');
ylabel('Y');
axis([-5,5,-5,5]);
colorbar;
axis equal;
view(2);
hold off