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

Contenido

1 .-Superficie Mallada
2 .-Función Potencial y Campo de Velocidades del Fluido.
- 2.1 .-Representación de la Función Potencial
- 2.2 .-Representación del Campo de Velocidades
3 .-Rotacional y Divergencia
- 3.1 .-Rotacional nulo
- 3.2 .-Divergencia nula
4 .-Líneas de corriente del campo [math]\vec{u}[/math]
5 .-Puntos de Frontera S y Remanso
- 5.1 .-Puntos de la frontera S
- 5.2 .-Puntos de remanso
6 .-Ecuación de Bernoulli
7 .-Trayectoria de la Partícula
8 .-Paradoja de D´Alembert
9 .-Reanalisis de los apartados 2,3 y 4

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 ofrecen información fundamental sobre las propiedades físicas del fluido que dicho campo describe.

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]

FALTA EDITAR 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 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

El módulo de la velocidad sobre la frontera [math]\rho = 1[/math] es:

[math]|u(1,\theta)| = 2 |\sin\theta|[/math]

Para hallar los máximos y mínimos, derivamos respecto a [math]\theta[/math]:

[math]\frac{d}{d\theta}|u| = 2 \cos\theta[/math]
Igualando a cero:

[math]\cos\theta = 0 \Rightarrow \theta = \pi/2, 3\pi/2[/math]

Por tanto:

Velocidad máxima: [math]|u| = 2[/math] en [math]\theta = \pi/2, 3\pi/2[/math]

 Coordenadas cartesianas: [math](0,1)[/math] y [math](0,-1)[/math]

Velocidad mínima (excepto puntos de remanso): [math]|u| = 0[/math] en [math]\theta = 0, \pi[/math]

 (estos ya se identificarán como puntos de remanso en la sección 5.2).

Para representar cómo varía el módulo de la velocidad con el ángulo [math]\theta[/math]:

theta = linspace(0,2*pi,1000);
u_mod = 2*abs(sin(theta));
plot(theta,u_mod,'LineWidth',2);
xlabel('\theta [rad]');
ylabel('

5.2 .-Puntos de remanso

En el borde del obstáculo [math]\rho = 1[/math], las componentes de la velocidad en coordenadas polares son:

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

Por tanto, el vector de velocidad es tangencial al borde y su módulo es:

[math]|u(1,\theta)| = 2 |\sin\theta|[/math]

Los puntos de remanso se encuentran donde [math]|u| = 0 \Rightarrow \sin\theta = 0[/math], es decir:

[math]\theta = 0 \quad y \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]

Se puede observar estos puntos visualmente codificando el flujo en MATLAB:

clear; clc;

% Mallado
r = linspace(1,5,50);   
theta_vec = linspace(0,2*pi,50);
[rho,theta] = meshgrid(r,theta_vec);

% Coordenadas cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

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

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

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

% Curvas de nivel y campo de velocidades
figure;
contour(X,Y,phi,50); hold on;
quiver(X,Y,DX,DY,'k');

% Obstáculo circular
plot(cos(theta_vec), sin(theta_vec), 'r', 'LineWidth',1.3);

% Puntos de remanso
plot([1,-1],[0,0],'b*','LineWidth',2,'MarkerSize',10);

axis equal; axis([-5 5 -5 5]);
colorbar;
title('Puntos de remanso y campo de velocidades');
xlabel('X'); ylabel('Y');
hold off;

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

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:

clear; clc; I = linspace(1,5,50); J = linspace(0,2*pi,80); [rho,theta] = meshgrid(I,J);

X = rho.*cos(theta); Y = rho.*sin(theta);

% Potencial f = @(rho,theta)((rho + 1./rho).*cos(theta)); contour(X,Y,f(rho,theta),50) hold on

% Componentes de la velocidad DX = (1 - 1./(rho.^2)).*cos(theta).*cos(theta) ...

  - (-(1 + 1./(rho.^2)).*sin(theta)).*sin(theta);

DY = (1 - 1./(rho.^2)).*cos(theta).*sin(theta) ...

  + (-(1 + 1./(rho.^2)).*sin(theta)).*cos(theta);

quiver(X,Y,DX,DY);

plot(cos(J),sin(J),'k','lineWidth',2);

axis([-5,5,-5,-5]); colorbar; title('Variacion de la Velocidad alrededor del Obstaculo'); xlabel('X'); ylabel('Y'); axis equal; view(2); hold off

8 .-Paradoja de D´Alembert

Para demostrar que el fluido no ejerce ningún empuje sobre el obstáculo en dirección horizontal, se resolverá la siguiente integral:

[math]\int p\vec n \vec i\,d\theta=0[/math]

El sumatorio de esta proyección calculado en dirección i es nulo, por lo que, el fluido no realiza ningún empuje sobre el obstáculo en la dirección horizontal. Sean:

p: presión particularizada en [math]\rho=1[/math]

[math]\vec n[/math]: vector perpendicular a la curva

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

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:

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