MateWiki - Contribuciones del usuario [es]

Flujo alrededor de un obstáculo circular (Grupo 26)

2025-11-28T11:16:48Z

Aaron Garcia: /* . Representación de la Función potencial */

{{ TrabajoED |Flujo alrededor de un obstáculo circular (Grupo 26) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] |
*Gonzalo Gallego Fulgencio
*Andrea García Carrasco
*Aarón García Martín
*Miryam Sánchez-Ferragut Samalea
*Guillermo Rodríguez Navadijos }}
Vamos a estudiar el flujo de un fluido incompresible alrededor de un obstáculo circular, trabajando en el plano y utilizando coordenadas cilíndricas (polares) para describir el campo de velocidades y las condiciones en la superficie del cilindro. Este enfoque permite formular de manera directa las ecuaciones del flujo potencial y analizar cómo la presencia del obstáculo modifica la distribución de velocidades y presiones. A partir de este planteamiento se desarrollarán las cuestiones que se piden a continuación.

==. Representación del mallado==
En este primer apartado representaremos la región ocupada por el fluido, que corresponde al exterior del círculo unidad. Para ello construiremos un mallado en coordenadas polares que cubra el anillo comprendido entre los radios 1 y 5, con centro en el origen. Este mallado permitirá visualizar los puntos interiores de la zona de estudio y establecer la geometría sobre la que se formulará posteriormente el problema del flujo. Para completar la representación, dibujaremos también los ejes cartesianos en el dominio
[
−
4
,
4
]
×
[
−
4
,
4
]
[−4,4]×[−4,4], lo que facilitará interpretar la posición del obstáculo circular y la extensión del fluido respecto al sistema de referencia.

[[Archivo:apartado1G26.png|400px|miniaturadeimagen| Mallado que representa los puntos de la región ocupada por un
fluido]]

{{matlab|codigo=
% Trabajo P - Apartado (1)
% Mallado del anillo 1 <= r <= 5 en coordenadas polares

clear; clc; close all;

R1 = 1; % radio interior (obstáculo)
R2 = 5; % radio exterior del fluido
Nr = 25; % número de divisiones radiales
Nth = 80; % número de divisiones angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);

[RHO, TH] = meshgrid(rho, theta);

X = RHO .* cos(TH);
Y = RHO .* sin(TH);
Z = 0.*RHO;
figure; hold on;

% Líneas radiales (theta = constante)
for i = 1:Nth
plot(X(i,:), Y(i,:), 'g');
end

% Circunferencias (r = constante)
for j = 1:Nr
plot(X(:,j), Y(:,j), 'g');
end

% Obstáculo circular (r = 1) representado solo con contorno
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % obstáculo circular

axis equal;
xlim([-4 4]);
ylim([-4 4]);
xlabel('x');
ylabel('y');
title('Mallado en el anillo 1 \leq r \leq 5 (flujo alrededor de un cilindro)');
grid off;
hold off;
}}

==. Función potencial y campo de velocidades del fluido==
En este apartado analizaremos la velocidad de las partículas dada por el gradiente de la siguiente función potencial:

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

===. Representación de la Función potencial===

Primero representaremos la función potencial que describe el flujo asociado al movimiento de un fluido incompresible alrededor de un obstáculo circular. Representaremos gráficamente la función potencial en el dominio exterior al círculo unidad para visualizar cómo varía en el plano y cómo organiza la estructura del flujo alrededor del cilindro.
[[Archivo:Curvasnivel26.png|400px|miniaturadeimagen| Función Potencial]]
{{matlab|codigo=
% Trabajo P - Apartado (2)
% Función potencial y campo de velocidades para
% phi(rho,theta) = (rho + 1/rho) * cos(theta)

clear; clc; close all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 40; % puntos radiales
Nth = 120; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Función potencial phi(rho,theta) = (rho + 1/rho) cos(theta)

phi = (RHO + 1./RHO) .* cos(TH);

% Campo de velocidades u = grad(phi)
% En polares: u_rho = dphi/drho, u_th = (1/rho) dphi/dth

u_rho = (1 - 1./RHO.^2) .* cos(TH);
u_th = -(1 + 1./RHO.^2) .* sin(TH);

% Pasamos a componentes cartesianas:
u_x = u_rho .* cos(TH) - u_th .* sin(TH);
u_y = u_rho .* sin(TH) + u_th .* cos(TH);

% Puntos del contorno del obstáculo (rho = 1)
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);

%Dibujamos las curvas de nivel del potencial
figure;
contour(X, Y, phi, 30); % 30 niveles de phi
hold on;
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % cilindro
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Curvas de nivel de la función potencial \phi');
hold off;
}}

===. Representación del campo de velocidades===
A partir de la función potencial, la velocidad del fluido se determina mediante su gradiente, <math>\vec{u}</math>
=∇φ.

<center><math>
\vec{u}=
\left(1-\frac{1}{\rho^2}\right)\cos\theta\,\vec{e}_\rho
-
\left(1+\frac{1}{\rho^2}\right)\sin\theta\,\vec{e}_\theta
</math></center>

Aquí representaremos el campo de velocidades resultante y analizaremos la dirección y magnitud del movimiento de las partículas del fluido, donde podremos observar que la velocidad es ortogonal a las curvas de nivel de
φ.

El campo <math> \vec u </math> lo hemos pasado manualmente a coordenadas cartesianas con la matriz de cambio de base para añadirlo directamente a nuestro código de Matlab. Dándonos como resultado:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\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)\right) \vec j </math></center>

[[Archivo:Campovelocidades26.png|400px|miniaturadeimagen| Campo de velocidades]]
[[Archivo:Campovelocidadesampl26.png |400px|miniaturadeimagen| Campo de velocidades ampliado]]
{{matlab|codigo=
clear; clc;clear all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 10; % puntos radiales
Nth = 70; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

%Definimos función potencial y la aplicamos a Z
f=@(rho,theta)(rho+(1./rho)).*cos(theta);
Z=f(RHO,TH);

%Dibujamos las curvas de nivel
contour(X,Y,Z,15);
hold on

%Definimos las componentes X e Y del gradiente
Gx=(1-(1./RHO.^2)).*cos(TH).^2+(1+(1./RHO.^2)).*sin(TH).^2;
Gy=(1-(1./RHO.^2)).*sin(TH).*cos(TH)-(1+(1./RHO.^2)).*sin(TH).*cos(TH);

%Dibujamos el campo de velocidades
quiver(X,Y,Gx,Gy);

% Representamos nuestro obstáculo
plot(1*cos(theta),1*sin(theta),'k','lineWidth',1);
axis([-4,4,-4,4]);
colorbar;
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}

==. Comprobación rotacional y divergencia nulos==
A partir del campo de velocidades calculado en el apartado anterior, calculamos su rotacional y su divergencia para conocer las características del fluido.
<center><math>
\vec{u}=
\left(1-\frac{1}{\rho^2}\right)\cos\theta\,\vec{e}_\rho
-
\left(1+\frac{1}{\rho^2}\right)\sin\theta\,\vec{e}_\theta
</math></center>
===. Comprobación del rotacional nulo===

Conociendo la fórmula del rotacional calculamos:
<center><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 & {0} \end{vmatrix}</math></center>

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

</math></center>

Obtenemos un rotacional nulo, es decir, se trata de un fluido irrotacional, por lo tanto, podemos deducir que las partículas de fluido no giran.

===. Comprobación de la divergencia nula===

Conociendo la fórmula de la divergencia calculamos:
<center><math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}(\rho(u_\rho))+\frac{\partial}{\partial{\theta}}(u_\theta)+\frac{\partial}{\partial{z}}(\rho(u_z))]</math></center>

<center><math>
\nabla\cdot\vec{u}
=
\frac{1}{\rho}
\left[
\frac{\partial}{\partial\rho}
\Bigl( \left(1-\frac{1}{\rho^{2}}\right)\cos\theta \; \rho\,\vec{e}_{\rho} \Bigr)
\;-\;
\frac{\partial}{\partial\theta}
\Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr)
\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></center>

Obtenemos una divergencia nula, es decir, significa que el fluido mantiene su volumen constante (ni se expande ni se contrae).

==. Líneas de corriente==

Primero calcularemos el campo <math>\vec{v}</math>, que en cada punto es ortogonal a <math>\vec{u}</math>, (<math>\vec{v}</math> = <math>\vec{k}\times\vec{u}</math>, donde <math>\vec{k}</math>=<math>\vec {e}_{z}</math>).

<center><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+\frac{1}{\rho^2})sin(\theta)\vec {e}_{\rho} + [(1-\frac{1}{\rho^2})cos(\theta)]\vec {e}_{\theta} =\vec v</math></center>

Comprobamos que <math>\vec{v}</math> es irrotacional:
<center><math>\nabla\times\vec v= \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}} \\ v_\rho & \rho v_\theta & {0} \end{vmatrix}=\frac{1}{\rho}[[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}-[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}]=\vec {0}</math></center>

A continuación calculamos <math>\psi</math>, para ello resolveremos el sistema de ecuaciones <math>\nabla\cdot\psi=\vec v</math>

<center><math>\frac{\partial\psi}{\partial\rho}=v_\rho=\int (1+\frac{1}{\rho^2})sen(\theta)\,d\rho=sen(\theta) (\rho-\frac{1}{\rho})+f(\theta)</math></center>

<center><math>\frac{\partial\psi}{\partial\theta}=\rho v_\theta=\int (\rho-\frac{1}{\rho})cos(\theta),d\theta=sen(\theta) (\rho-\frac{1}{\rho})</math></center>

De ello se obtendrá la siguiente igualdad, que representa el potencial escalar, y se conoce como función de corriente de <math>\vec{u}</math>

<center><math>\psi = \sin(\theta)\left(\rho - \frac{1}{\rho}\right)</math></center>

A continuación, se representa el campo y el potencial escalar.

[[Archivo:apartado44G26.png|400px|thumb|left|Lineas de corriente ]]

{{matlab|codigo=

u = linspace(1,5,250);
v = linspace(0,2*pi,250);
[rho,th] = meshgrid(u,v);

Mx = rho.*cos(th);
My = rho.*sin(th);

% CirculaciÃ³n
Gamma = 1/2;

psi = (rho - 1./rho).*sin(th) - (Gamma/(2*pi))*log(rho);

% Velocidades en polares
u_r = (1 - 1./rho.^2).*cos(th);
u_th = -(1 + 1./rho.^2).*sin(th) + Gamma./(2*pi*rho);

% Velocidades en cartesianas
Ux = u_r.*cos(th) - u_th.*sin(th);
Uy = u_r.*sin(th) + u_th.*cos(th);

% quitar flechas
step = 12;

Mx_q = Mx(1:step:end, 1:step:end);
My_q = My(1:step:end, 1:step:end);
Ux_q = Ux(1:step:end, 1:step:end);
Uy_q = Uy(1:step:end, 1:step:end);

figure;
hold on;

contour(Mx, My, psi, 80);
quiver(Mx_q, My_q, Ux_q, Uy_q, 'k');
axis equal;
xlabel

}}

Gracias a esta representación observamos que el potencial escalar y campo de velocidades son paralelos, concluyendo que son efectivamente líneas de corriente de <math>\vec{u}</math>.

==. Puntos de la frontera S==
En la frontera del cilindro se tiene <math>\rho = 1</math>.

Las componentes del campo de velocidades son:

<center><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></center>

Sustituyendo <math>\rho = 1</math>:

<center><math>
u_\rho(1,\theta) = 0, \qquad
u_\theta(1,\theta) = -2\sin\theta.
</math></center>

La rapidez en la frontera es:

<center><math>
\left\lvert \vec u(1,\theta)\right\rvert
= \sqrt{u_\rho^2 + u_\theta^2}
= 2\lvert \sin\theta\rvert.
</math></center>

===. Velocidad máxima ===

La velocidad es máxima cuando:

<center><math>
\lvert \sin\theta\rvert = 1 \quad\Longrightarrow\quad
\theta = \frac{\pi}{2},\ \frac{3\pi}{2}.
</math></center>

Coordenadas sobre el cilindro:<math>
(0,1), \qquad (0,-1).
</math>

===. Velocidad mínima ===

La rapidez es mínima cuando:

<center><math>
\lvert \sin\theta\rvert = 0 \quad\Longrightarrow\quad
\theta = 0,\ \pi.
</math></center>

Coordenadas:<math>
(1,0) \quad\text{y}\quad (-1,0).
</math>

===. Puntos de remanso ===

Los puntos de remanso son aquellos donde la velocidad es nula:

<center><math>
u_\rho = 0, \qquad u_\theta = 0.
</math></center>

Esto ocurre cuando:

<center><math>
\sin\theta = 0 \quad\Longrightarrow\quad
\theta = 0,\ \pi.
</math></center>

Por tanto, los puntos de remanso son:<math>
(1,0) \quad\text{y}\quad (-1,0).
</math>

==Presión del fluido==

En este apartado se calcula el campo de presiones del flujo usando la ecuación de Bernoulli.

Las componentes del campo de velocidades son:

<center><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></center>

La rapidez viene dada por

<center><math>
\lvert\vec u\rvert^2 = u_\rho^2 + u_\theta^2
</math></center>

Sustituyendo:

<center><math>
\begin{aligned}
\lvert\vec u\rvert^2
&= \left(1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)\cos^2\theta
+ \left(1 + \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)\sin^2\theta \\
&= 1 + \frac{1}{\rho^4} - \frac{2}{\rho^2}\cos 2\theta
\end{aligned}
</math></center>

La ecuación de Bernoulli para un flujo incompresible e irrotacional es

<center><math>
p + \frac{1}{2}\rho_f \lvert\vec u\rvert^2
= p_\infty + \frac{1}{2}\rho_f U_\infty^2,
</math></center>

donde <math>U_\infty = 1</math>. Por tanto,

<center><math>
p(\rho,\theta)
= p_\infty + \frac{1}{2}\rho_f\left(1 - \lvert\vec u\rvert^2\right).
</math></center>

Sustituyendo la expresión de la velocidad:

<center><math>
\begin{aligned}
p(\rho,\theta)
&= p_\infty + \frac{1}{2}\rho_f\left[
1 - \left(1 + \frac{1}{\rho^4} - \frac{2}{\rho^2}\cos 2\theta\right)
\right] \\
&= p_\infty + \rho_f\left(
\frac{\cos 2\theta}{\rho^2} - \frac{1}{2\rho^4}
\right).
\end{aligned}
</math></center>

Por tanto, el campo de presiones es:

<center><math>
{
p(\rho,\theta)
= p_\infty
+ \rho_f\left(\frac{\cos 2\theta}{\rho^2}
.⁠ ⁠\frac{1}{2\rho^4}\right)
}
</math></center>

=== Presión sobre la superficie del cilindro (ρ = 1) ===

En <math>\rho = 1</math> la rapidez es

<center><math>
\lvert\vec u(1,\theta)\rvert^2 = 4\sin^2\theta,
</math></center>

y entonces Bernoulli da

<center><math>
p(1,\theta)
= p_\infty + \frac{1}{2}\rho_f(1 - 4\sin^2\theta)
= p_\infty - \frac{1}{2}\rho_f + \rho_f \cos 2\theta.
</math></center>

=== Puntos de máxima y mínima presión ===

En los puntos de remanso (<math>\theta = 0,\pi</math>), donde la velocidad es nula:

<center><math>
p = p_\infty + \frac{1}{2}\rho_f.
</math></center>

En los puntos de máxima velocidad (<math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>):

<center><math>
p = p_\infty - \frac{3}{2}\rho_f.
</math></center>

La presión disminuye donde aumenta la velocidad, en concordancia con la ecuación de Bernoulli.

===Código de la presión de fluido===

[[Archivo:apartado6G26.png|400px|thumb|left|Presión de fluido]]

{{matlab|codigo=
clear; clc;clear all;

function presion_apartado6()

p_inf = 0;
rho_f = 1;

R1 = 1;
R2 = 5;
Nr = 180;
Nth = 360;

r = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RR, TT] = meshgrid(r, theta);

P = p_inf + rho_f*( cos(2*TT)./RR.^2 - 1./(2*RR.^4) );

% Pasar a coordenadas cartesianas para dibujar
X = RR .* cos(TT);
Y = RR .* sin(TT);

th = linspace(0, 2*pi, 400);
xc = cos(th);
yc = sin(th);

figure;
contourf(X, Y, P, 40, "LineColor", "none");
hold on;
plot(xc, yc, 'k', 'LineWidth', 2);
colorbar;
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Campo de presiones p(r,\theta) del Apartado 6');
hold off;
end
}}

En la imagen vemos el campo de presiones alrededor del cilindro de radio 1 cuando el flujo pasa a velocidad uniforme u=1
Las zonas amarillas representan las presiones mas altas que se consideran los puntos theta=0 y theta=pi. Los consideramos puntos de remanso, lugar donde la velocidad cae a cero y la presión sube al máximo .
Las zonas azul y verde representan la zona de menor presión que son theta=pi/2 y theta=3pi/2 .El fluido en este caso acelera para bordear el cilindro luego llega a velocidad máxima y presión mínima.

==Partícula del fluido==

En este apartado analizamos la trayectoria que seguiría una partícula del fluido y cómo cambian la
velocidad y la presión al rodear el obstáculo.

=== Trayectorias y líneas de corriente ===

En un flujo estacionario e incompresible, las trayectorias de las partículas coinciden con las líneas de corriente.
Por tanto, una partícula del fluido seguirá exactamente las curvas que verifican:

<center><math>
\psi(\rho,\theta) = \text{cte},
</math></center>

donde la función corriente es

<center><math>
\psi(\rho,\theta)
= \left(\rho - \frac{1}{\rho}\right)\sin\theta.
</math></center>

Estas líneas describen las trayectorias del fluido alrededor del cilindro y muestran cómo la partícula se desvía
en torno al obstáculo siguiendo la geometría del flujo potencial.

=== Variación de la velocidad al rodear el cilindro ===

La rapidez del fluido viene dada por

<center><math>
\lvert\vec u\rvert
= \sqrt{
\left(1 - \frac{1}{\rho^2}\right)^2 \cos^2\theta
+
\left(1 + \frac{1}{\rho^2}\right)^2 \sin^2\theta
}.
</math></center>

En la superficie del cilindro (<math>\rho = 1</math>) se simplifica a

<center><math>
\lvert\vec u(1,\theta)\rvert = 2\lvert\sin\theta\rvert.
</math></center>

Por tanto:

* La velocidad es máxima en <math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>.
* La velocidad se anula en los puntos de remanso: <math>\theta = 0, \pi</math>.

La partícula acelera al desplazarse hacia la parte superior e inferior del cilindro y se frena al
pasar por los puntos frontales y traseros del obstáculo.

=== Variación de la presión ===

Según Bernoulli,

<center><math>
p + \frac{1}{2}\rho_f\lvert\vec u\rvert^2
= p_\infty + \frac{1}{2}\rho_f.
</math></center>

La presión disminuye cuando la velocidad aumenta. Aplicando esto sobre la superficie del cilindro:

* En los puntos de remanso (<math>\theta = 0, \pi</math>) la presión es máxima.
* En los puntos de mayor velocidad (<math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>) la presión es mínima.

==Circulación del campo==

En este apartado comprobamos que la circulación del campo de velocidades alrededor de una circunferencia
de radio 1 es nula en el caso sin circulación añadida. Asimismo, se explica la relación entre este hecho,
la fuerza ejercida por el fluido y la paradoja de D’Alembert.

=== Circulación del campo de velocidades ===

La circulación alrededor de una curva cerrada <math>C</math> se define como:

<center><math>
\Gamma = \oint_C \vec u \cdot d\vec s.
</math></center>

Tomamos como curva de referencia la circunferencia de radio 1: <math>\rho = 1</math>.
El campo de velocidades sobre el cilindro es:

<center><math>
u_\rho(1,\theta) = 0, \qquad
u_\theta(1,\theta) = -2\sin\theta.
</math></center>

Sobre la circunferencia, el elemento de arco es

<center><math>
d\vec s = \hat{e}\theta \, \rho \, d\theta = \hat{e}\theta \, d\theta.
</math></center>

Por tanto, la circulación es:

<center><math>
\Gamma = \int_0^{2\pi} u_\theta(1,\theta)\, d\theta
= \int_0^{2\pi} -2\sin\theta \, d\theta.
</math></center>

Como la integral de <math>\sin\theta</math> en un período completo es cero, obtenemos:

<center><math>
\Gamma = 0.
</math></center>

=== Relación con la fuerza sobre el cilindro ===

Lo relacionamos directamente con la circulación mediante el teorema de Kutta–Joukowski:

<center><math>
F = -\rho_f U_\infty \Gamma.
</math></center>

Como en este caso <math>\Gamma = 0</math>, la fuerza resultante es:

<center><math>
F = 0.
</math></center>

Es decir, a pesar de que el fluido se desvía alrededor del cilindro no aparece fuerza neta sobre él.

Este resultado es una manifestación de la paradoja de D’Alembert en un flujo potencial ideal y sin viscosidad la fuerza sobre un obstáculo fijo es exactamente cero, lo cual contradice la experiencia real.

==. Nueva función potencial==
===. Nueva representación del Potencial y del campo de velocidades===
Ahora supondremos que la velocidad de las partículas del fluido viene dada por el gradiente de la siguiente función potencial:

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

A continuación repetiremos el mismo proceso anterior. Primero representaremos la función potencial.

[[Archivo:Curvasnivel926.png|400px|miniaturadeimagen| Función Potencial 2]]
{{matlab|codigo=
% Función potencial y campo de velocidades para
% phi(rho,theta) = (rho + 1/rho) * cos(theta) + theta/4*pi

clear; clc; close all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 40; % puntos radiales
Nth = 120; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Función potencial phi(r,theta) = (r + 1/r) cos(theta)+ theta/4*pi

phi = (RHO + 1./RHO) .* cos(TH) + TH./(4*pi);

% Campo de velocidades u = grad(phi)
% En polares: u_rho = dphi/drho, u_th = (1/rho) dphi/dth

u_rho = (1 - 1./RHO.^2) .* cos(TH);
% dphi/dtheta
dphi_dtheta = -(RHO + 1./RHO) .* sin(TH) + 1/(4*pi);

% u_theta = (1/r)*dphi/dtheta
u_th = (1./RHO) .* dphi_dtheta;

% Pasamos a componentes cartesianas:
u_x = u_rho .* cos(TH) - u_th .* sin(TH);
u_y = u_rho .* sin(TH) + u_th .* cos(TH);

% Puntos del contorno del obstáculo (r = 1)
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);

%Dibujamos las curvas de nivel del potencial
figure;
contour(X, Y, phi, 30); % 30 niveles de phi
hold on;
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % cilindro
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Curvas de nivel de la función potencial \phi');
hold off;
}}

A partir de la función potencial, la velocidad del fluido se determina mediante su gradiente, <math>\vec{u}</math>
=∇φ.

<center><math>
\vec{u}(\rho,\theta,z)=\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\mathbf e_{\rho}+\frac{1}{\rho}\left[-\left(\rho+\frac{1}{\rho}\right)\sin\theta+\frac{1}{4\pi}\right]\mathbf e_{\theta}
</math></center>

Aquí representaremos el campo de velocidades resultante y analizaremos la dirección y magnitud del movimiento de las partículas del fluido.

El campo <math> \vec u </math> lo hemos pasado manualmente a coordenadas cartesianas con la matriz de cambio de base para añadirlo directamente a nuestro código de Matlab. Dándonos como resultado:
<center><math> \vec u (\vec i,\vec j,\vec k)=\left[\cos^{2}\theta\left(1-\frac{1}{\rho^{2}}\right)-\sin\theta\cos\theta+\frac{\sin\theta}{4\pi\rho}\right]\mathbf i+\left[\sin\theta\cos\theta\left(1-\frac{1}{\rho^{2}}\right)+\cos^{2}\theta+\frac{\cos\theta}{4\pi\rho}\right]\mathbf j </math></center>

[[Archivo:Campovelocidades926.png|400px|miniaturadeimagen| Campo de Velocidades resultante]]
[[Archivo:Campovelocidadesampl926.png|400px|miniaturadeimagen| Campo de Velocidades resultante]]
{{matlab|codigo=
clear; clc;clear all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 10; % puntos radiales
Nth = 70; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

%Definimos función potencial y la aplicamos a Z
f=@(rho,theta)(rho+(1./rho)).*cos(theta) + theta./(4.*pi);
Z=f(RHO,TH);

% Dibujamos las curvas de nivel
contour(X,Y,Z,15);
hold on

%Definimos las componentes X e Y del gradiente en coordenadas cartesianas
Gx=(cos(TH).^2).*(1 - 1./RHO.^2) - sin(TH).*cos(TH) + sin(TH)./(4*pi.*RHO);
Gy=sin(TH).*cos(TH).*(1 - 1./RHO.^2) + cos(TH).^2 + cos(TH)./(4*pi.*RHO);

%Dibujamos el campo de velocidades
quiver(X,Y,Gx,Gy);

% Representamos nuestro obstáculo
plot(1*cos(theta),1*sin(theta),'k','lineWidth',1);
axis([-4,4,-4,4]);
colorbar;
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}

Observamos como nuestro campo es ortogonal a las curvas de nivel.

===. Comprobación rotacional y divergencia nulos===

A continuación, comprobaremos que el rotacional y la divergencia sean nulos:

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

Conociendo el campo de velocidades calculado anteriormente:

<center><math>
\vec{u}
= u_{\rho}\,\vec{e}_{\rho}+u_{\theta}\,\vec{e}_{\theta}
=\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\vec{e}_{\rho}
-\left(1+\frac{1}{\rho^{2}}\right)\sin\theta\,\vec{e}_{\theta}
+\frac{1}{4\pi\rho}\,\vec{e}_{\theta}.
</math></center>

Comprobamos el rotacional nulo:

<center><math>
\nabla\times\vec{u}
=\frac{1}{\rho}
\begin{vmatrix}
\vec{e}_{\rho} & \rho\vec{e}_{\theta} & \vec{e}_{z} \\
\dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial z} \\
u_{\rho} & \rho u_{\theta} & 0
\end{vmatrix},
</math></center>

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

</math></center>

Obtenemos un rotacional nulo, por lo que el flujo sigue siendo irrotacional y las partículas de fluido no giran localmente.

Comprobamos la divergencia nula:

<center><math>
\nabla\cdot\vec{u}
=\frac{1}{\rho}\frac{\partial}{\partial\rho}\bigl(\rho u_{\rho}\bigr)
+\frac{1}{\rho}\frac{\partial u_{\theta}}{\partial\theta}
+\frac{\partial u_{z}}{\partial z},
\qquad u_{z}=0.
</math></center>

<center><math>
\nabla\cdot\vec{u}
=\frac{1}{\rho}\left[
\frac{\partial}{\partial\rho}\left(\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\rho\right)
+\frac{\partial}{\partial\theta}\left(-\left(1+\frac{1}{\rho^{2}}\right)\sin\theta+\frac{1}{4\pi\rho}\right)
+\frac{\partial}{\partial z}(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></center>

Obtenemos una divergencia nula, es decir, significa que el fluido mantiene su volumen constante (ni se expande ni se contrae), de modo que se trata de un flujo incompresible.

===Lineas de corriente===

==Póster==

[[Categoría:Matemáticas I]]
[[Categoría:MatI/19]]

Flujo alrededor de un obstáculo circular (Grupo 26)

2025-11-28T10:09:22Z

Aaron Garcia: /* Partícula del fluido */

{{ TrabajoED | Flujo alrededor de un obstáculo circular (Grupo 26) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] |
*Gonzalo Gallego Fulgencio
*Andrea García Carrasco
*Aarón García Martín
*Miryam Sánchez-Ferragut Samalea
*Guillermo Rodríguez Navadijos }}
Vamos a estudiar el flujo de un fluido incompresible alrededor de un obstáculo circular, trabajando en el plano y utilizando coordenadas cilíndricas (polares) para describir el campo de velocidades y las condiciones en la superficie del cilindro. Este enfoque permite formular de manera directa las ecuaciones del flujo potencial y analizar cómo la presencia del obstáculo modifica la distribución de velocidades y presiones. A partir de este planteamiento se desarrollarán las cuestiones que se piden a continuación.

==. Representación del mallado==
En este primer apartado representaremos la región ocupada por el fluido, que corresponde al exterior del círculo unidad. Para ello construiremos un mallado en coordenadas polares que cubra el anillo comprendido entre los radios 1 y 5, con centro en el origen. Este mallado permitirá visualizar los puntos interiores de la zona de estudio y establecer la geometría sobre la que se formulará posteriormente el problema del flujo. Para completar la representación, dibujaremos también los ejes cartesianos en el dominio
[
−
4
,
4
]
×
[
−
4
,
4
]
[−4,4]×[−4,4], lo que facilitará interpretar la posición del obstáculo circular y la extensión del fluido respecto al sistema de referencia.

[[Archivo:apartado1G26.png|400px|miniaturadeimagen| Mallado que representa los puntos de la región ocupada por un
fluido]]

{{matlab|codigo=
% Trabajo P - Apartado (1)
% Mallado del anillo 1 <= r <= 5 en coordenadas polares

clear; clc; close all;

R1 = 1; % radio interior (obstáculo)
R2 = 5; % radio exterior del fluido
Nr = 25; % número de divisiones radiales
Nth = 80; % número de divisiones angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);

[RHO, TH] = meshgrid(rho, theta);

X = RHO .* cos(TH);
Y = RHO .* sin(TH);
Z = 0.*RHO;
figure; hold on;

% Líneas radiales (theta = constante)
for i = 1:Nth
plot(X(i,:), Y(i,:), 'g');
end

% Circunferencias (r = constante)
for j = 1:Nr
plot(X(:,j), Y(:,j), 'g');
end

% Obstáculo circular (r = 1) representado solo con contorno
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % obstáculo circular

axis equal;
xlim([-4 4]);
ylim([-4 4]);
xlabel('x');
ylabel('y');
title('Mallado en el anillo 1 \leq r \leq 5 (flujo alrededor de un cilindro)');
grid off;
hold off;
}}

==. Función potencial y campo de velocidades del fluido==
En este apartado analizaremos la velocidad de las partículas dada por el gradiente de la siguiente función potencial:

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

===. Representación de la Función potencial===

Primero representaremos la función potencial que describe el flujo asociado al movimiento de un fluido incompresible alrededor de un obstáculo circular. Representaremos gráficamente la función potencial en el dominio exterior al círculo unidad para visualizar cómo varía en el plano y cómo organiza la estructura del flujo alrededor del cilindro.
[[Archivo:Curvasnivel26.png|400px|miniaturadeimagen| Función Potencial]]
{{matlab|codigo=
% Trabajo P - Apartado (2)
% Función potencial y campo de velocidades para
% phi(r,theta) = (r + 1/r) * cos(theta)

clear; clc; close all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 40; % puntos radiales
Nth = 120; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Función potencial phi(r,theta) = (r + 1/r) cos(theta)

phi = (RHO + 1./RHO) .* cos(TH);

% Campo de velocidades u = grad(phi)
% En polares: u_rho = dphi/drho, u_th = (1/rho) dphi/dth

u_rho = (1 - 1./RHO.^2) .* cos(TH);
u_th = -(1 + 1./RHO.^2) .* sin(TH);

% Pasamos a componentes cartesianas:
u_x = u_rho .* cos(TH) - u_th .* sin(TH);
u_y = u_rho .* sin(TH) + u_th .* cos(TH);

% Puntos del contorno del obstáculo (r = 1)
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);

%Dibujamos las curvas de nivel del potencial
figure;
contour(X, Y, phi, 30); % 30 niveles de phi
hold on;
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % cilindro
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Curvas de nivel de la función potencial \phi');
hold off;
}}

===. Representación del campo de velocidades===
A partir de la función potencial, la velocidad del fluido se determina mediante su gradiente, <math>\vec{u}</math>
=∇φ.

<center><math>
\vec{u}=
\left(1-\frac{1}{\rho^2}\right)\cos\theta\,\vec{e}_\rho
-
\left(1+\frac{1}{\rho^2}\right)\sin\theta\,\vec{e}_\theta
</math></center>

Aquí representaremos el campo de velocidades resultante y analizaremos la dirección y magnitud del movimiento de las partículas del fluido, donde podremos observar que la velocidad es ortogonal a las curvas de nivel de
φ.

El campo <math> \vec u </math> lo hemos pasado manualmente a coordenadas cartesianas con la matriz de cambio de base para añadirlo directamente a nuestro código de Matlab. Dándonos como resultado:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\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)\right) \vec j </math></center>

[[Archivo:Campovelocidades26.png|400px|miniaturadeimagen| Campo de velocidades]]
[[Archivo:Campovelocidadesampl26.png |400px|miniaturadeimagen| Campo de velocidades ampliado]]
{{matlab|codigo=
clear; clc;clear all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 10; % puntos radiales
Nth = 70; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

%Definimos función potencial y la aplicamos a Z
f=@(rho,theta)(rho+(1./rho)).*cos(theta);
Z=f(RHO,TH);

%Dibujamos las curvas de nivel
contour(X,Y,Z,15);
hold on

%Definimos las componentes X e Y del gradiente
Gx=(1-(1./RHO.^2)).*cos(TH).^2+(1+(1./RHO.^2)).*sin(TH).^2;
Gy=(1-(1./RHO.^2)).*sin(TH).*cos(TH)-(1+(1./RHO.^2)).*sin(TH).*cos(TH);

%Dibujamos el campo de velocidades
quiver(X,Y,Gx,Gy);

% Representamos nuestro obstáculo
plot(1*cos(theta),1*sin(theta),'k','lineWidth',1);
axis([-4,4,-4,4]);
colorbar;
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}

==. Comprobación rotacional y divergencia nulos==
A partir del campo de velocidades calculado en el apartado anterior, calculamos su rotacional y su divergencia para conocer las características del fluido.
<center><math>
\vec{u}=
\left(1-\frac{1}{\rho^2}\right)\cos\theta\,\vec{e}_\rho
-
\left(1+\frac{1}{\rho^2}\right)\sin\theta\,\vec{e}_\theta
</math></center>
===. Comprobación del rotacional nulo===

Conociendo la fórmula del rotacional calculamos:
<center><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 & {0} \end{vmatrix}</math></center>

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

</math></center>

Obtenemos un rotacional nulo, es decir, se trata de un fluido irrotacional, por lo tanto, podemos deducir que las partículas de fluido no giran.

===. Comprobación de la divergencia nula===

Conociendo la fórmula de la divergencia calculamos:
<center><math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}(\rho(u_\rho))+\frac{\partial}{\partial{\theta}}(u_\theta)+\frac{\partial}{\partial{z}}(\rho(u_z))]</math></center>

<center><math>
\nabla\cdot\vec{u}
=
\frac{1}{\rho}
\left[
\frac{\partial}{\partial\rho}
\Bigl( \left(1-\frac{1}{\rho^{2}}\right)\cos\theta \; \rho\,\vec{e}_{\rho} \Bigr)
\;-\;
\frac{\partial}{\partial\theta}
\Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr)
\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></center>

Obtenemos una divergencia nula, es decir, significa que el fluido mantiene su volumen constante (ni se expande ni se contrae).

==. Líneas de corriente==

Primero calcularemos el campo <math>\vec{v}</math>, que en cada punto es ortogonal a <math>\vec{u}</math>, (<math>\vec{v}</math> = <math>\vec{k}\times\vec{u}</math>, donde <math>\vec{k}</math>=<math>\vec {e}_{z}</math>).

<center><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+\frac{1}{\rho^2})sin(\theta)\vec {e}_{\rho} + [(1-\frac{1}{\rho^2})cos(\theta)]\vec {e}_{\theta} =\vec v</math></center>

Comprobamos que <math>\vec{v}</math> es irrotacional:
<center><math>\nabla\times\vec v= \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}} \\ v_\rho & \rho v_\theta & {0} \end{vmatrix}=\frac{1}{\rho}[[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}-[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}]=\vec {0}</math></center>

A continuación calculamos <math>\psi</math>, para ello resolveremos el sistema de ecuaciones <math>\nabla\cdot\psi=\vec v</math>

<center><math>\frac{\partial\psi}{\partial\rho}=v_\rho=\int (1+\frac{1}{\rho^2})sen(\theta)\,d\rho=sen(\theta) (\rho-\frac{1}{\rho})+f(\theta)</math></center>

<center><math>\frac{\partial\psi}{\partial\theta}=\rho v_\theta=\int (\rho-\frac{1}{\rho})cos(\theta),d\theta=sen(\theta) (\rho-\frac{1}{\rho})</math></center>

<center><math>\psi = \sin(\theta)\left(\rho - \frac{1}{\rho}\right)</math></center>

==. Puntos de la frontera S==
En la frontera del cilindro se tiene <math>\rho = 1</math>.

Las componentes del campo de velocidades son:

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

Sustituyendo <math>\rho = 1</math>:

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

La rapidez en la frontera es:

<math>
\left\lvert \vec u(1,\theta)\right\rvert
= \sqrt{u_\rho^2 + u_\theta^2}
= 2\lvert \sin\theta\rvert.
</math>

===. Velocidad máxima ===

La velocidad es máxima cuando:

<math>
\lvert \sin\theta\rvert = 1 \quad\Longrightarrow\quad
\theta = \frac{\pi}{2},\ \frac{3\pi}{2}.
</math>

Coordenadas sobre el cilindro:

<math>
(0,1), \qquad (0,-1).
</math>

===. Velocidad mínima ===

La rapidez es mínima cuando:

<math>
\lvert \sin\theta\rvert = 0 \quad\Longrightarrow\quad
\theta = 0,\ \pi.
</math>

Coordenadas:

<math>
(1,0), \qquad (-1,0).
</math>

===. Puntos de remanso ===

Los puntos de remanso son aquellos donde la velocidad es nula:

<math>
u_\rho = 0, \qquad u_\theta = 0.
</math>

Esto ocurre cuando:

<math>
\sin\theta = 0 \quad\Longrightarrow\quad
\theta = 0,\ \pi.
</math>

Por tanto, los puntos de remanso son:

<math>
(1,0) \quad\text{y}\quad (-1,0).
</math>

==Presión del fluido==

En este apartado se calcula el campo de presiones del flujo usando la ecuación de Bernoulli.

Las componentes del campo de velocidades son:

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

La rapidez viene dada por

<math>
\lvert\vec u\rvert^2 = u_\rho^2 + u_\theta^2.
</math>

Sustituyendo:

<math>
\begin{aligned}
\lvert\vec u\rvert^2
&= \left(1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)\cos^2\theta
+ \left(1 + \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)\sin^2\theta \\
&= 1 + \frac{1}{\rho^4} - \frac{2}{\rho^2}\cos 2\theta.
\end{aligned}
</math>

La ecuación de Bernoulli para un flujo incompresible e irrotacional es

<math>
p + \frac{1}{2}\rho_f \lvert\vec u\rvert^2
= p_\infty + \frac{1}{2}\rho_f U_\infty^2,
</math>

donde <math>U_\infty = 1</math>. Por tanto,

<math>
p(\rho,\theta)
= p_\infty + \frac{1}{2}\rho_f\left(1 - \lvert\vec u\rvert^2\right).
</math>

Sustituyendo la expresión de la velocidad:

<math>
\begin{aligned}
p(\rho,\theta)
&= p_\infty + \frac{1}{2}\rho_f\left[
1 - \left(1 + \frac{1}{\rho^4} - \frac{2}{\rho^2}\cos 2\theta\right)
\right] \\
&= p_\infty + \rho_f\left(
\frac{\cos 2\theta}{\rho^2} - \frac{1}{2\rho^4}
\right).
\end{aligned}
</math>

Por tanto, el campo de presiones es:

<math>
\boxed{
p(\rho,\theta)
= p_\infty
+ \rho_f\left(\frac{\cos 2\theta}{\rho^2}
•⁠ ⁠\frac{1}{2\rho^4}\right)
}
</math>

=== Presión sobre la superficie del cilindro (ρ = 1) ===

En <math>\rho = 1</math> la rapidez es

<math>
\lvert\vec u(1,\theta)\rvert^2 = 4\sin^2\theta,
</math>

y entonces Bernoulli da

<math>
p(1,\theta)
= p_\infty + \frac{1}{2}\rho_f(1 - 4\sin^2\theta)
= p_\infty - \frac{1}{2}\rho_f + \rho_f \cos 2\theta.
</math>

=== Puntos de máxima y mínima presión ===

En los puntos de remanso (<math>\theta = 0,\pi</math>), donde la velocidad es nula:

<math>
p = p_\infty + \frac{1}{2}\rho_f.
</math>

En los puntos de máxima velocidad (<math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>):

<math>
p = p_\infty - \frac{3}{2}\rho_f.
</math>

La presión disminuye donde aumenta la velocidad, en concordancia con la ecuación de Bernoulli.

===Código de la presión de fluido===

[[Archivo:apartado6G26.png|400px|thumb|left|texto alternativo]]

{{matlab|codigo=
clear; clc;clear all;

function presion_apartado6()

p_inf = 0;
rho_f = 1;

R1 = 1;
R2 = 5;
Nr = 180;
Nth = 360;

r = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RR, TT] = meshgrid(r, theta);

P = p_inf + rho_f*( cos(2*TT)./RR.^2 - 1./(2*RR.^4) );

% Pasar a coordenadas cartesianas para dibujar
X = RR .* cos(TT);
Y = RR .* sin(TT);

th = linspace(0, 2*pi, 400);
xc = cos(th);
yc = sin(th);

figure;
contourf(X, Y, P, 40, "LineColor", "none");
hold on;
plot(xc, yc, 'k', 'LineWidth', 2);
colorbar;
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Campo de presiones p(r,\theta) del Apartado 6');
hold off;
end
}}

En la imagen vemos el campo de presiones alrededor del cilindro de radio 1 cuando el flujo pasa a velocidad uniforme u=1
Las zonas amarillas representan las presiones mas altas que se consideran los puntos theta=0 y theta=pi. Los consideramos puntos de remanso, lugar donde la velocidad cae a cero y la presión sube al máximo .
Las zonas azul y verde representan la zona de menor presión que son theta=pi/2 y theta=3pi/2 .El fluido en este caso acelera para bordear el cilindro luego llega a velocidad máxima y presión mínima.

==Partícula del fluido==

En este apartado analizamos la trayectoria que seguiría una partícula del fluido y cómo cambian la
velocidad y la presión al rodear el obstáculo.

=== Trayectorias y líneas de corriente ===

En un flujo estacionario e incompresible, las trayectorias de las partículas coinciden con las líneas de corriente.
Por tanto, una partícula del fluido seguirá exactamente las curvas que verifican:

<math>
\psi(\rho,\theta) = \text{cte},
</math>

donde la función corriente es

<math>
\psi(\rho,\theta)
= \left(\rho - \frac{1}{\rho}\right)\sin\theta.
</math>

Estas líneas describen las trayectorias del fluido alrededor del cilindro y muestran cómo la partícula se desvía
en torno al obstáculo siguiendo la geometría del flujo potencial.

=== Variación de la velocidad al rodear el cilindro ===

La rapidez del fluido viene dada por

<math>
\lvert\vec u\rvert
= \sqrt{
\left(1 - \frac{1}{\rho^2}\right)^2 \cos^2\theta
+
\left(1 + \frac{1}{\rho^2}\right)^2 \sin^2\theta
}.
</math>

En la superficie del cilindro (<math>\rho = 1</math>) se simplifica a

<math>
\lvert\vec u(1,\theta)\rvert = 2\lvert\sin\theta\rvert.
</math>

Por tanto:

* La velocidad **es máxima** en <math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>.
* La velocidad **se anula** en los puntos de remanso: <math>\theta = 0, \pi</math>.

La partícula acelera al desplazarse hacia la parte superior e inferior del cilindro y se frena al
pasar por los puntos frontales y traseros del obstáculo.

=== Variación de la presión ===

Según Bernoulli,

<math>
p + \frac{1}{2}\rho_f\lvert\vec u\rvert^2
= p_\infty + \frac{1}{2}\rho_f.
</math>

La presión disminuye cuando la velocidad aumenta. Aplicando esto sobre la superficie del cilindro:

* En los puntos de remanso (<math>\theta = 0, \pi</math>) la presión es **máxima**.
* En los puntos de mayor velocidad (<math>\theta = \frac{\pi}{2}, \frac{3\pi}{2}</math>) la presión es **mínima**.

==Circulación del campo==

En este apartado comprobamos que la circulación del campo de velocidades alrededor de una circunferencia
de radio 1 es nula en el caso sin circulación añadida. Asimismo, se explica la relación entre este hecho,
la fuerza ejercida por el fluido y la paradoja de D’Alembert.

=== Circulación del campo de velocidades ===

La circulación alrededor de una curva cerrada <math>C</math> se define como:

<math>
\Gamma = \oint_C \vec u \cdot d\vec s.
</math>

Tomamos como curva de referencia la circunferencia de radio 1: <math>\rho = 1</math>.
El campo de velocidades sobre el cilindro es:

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

Sobre la circunferencia, el elemento de arco es

<math>
d\vec s = \hat{e}\theta \, \rho \, d\theta = \hat{e}\theta \, d\theta.
</math>

Por tanto, la circulación es:

<math>
\Gamma = \int_0^{2\pi} u_\theta(1,\theta)\, d\theta
= \int_0^{2\pi} -2\sin\theta \, d\theta.
</math>

Como la integral de <math>\sin\theta</math> en un período completo es cero, obtenemos:

<math>
\Gamma = 0.
</math>

=== Relación con la fuerza sobre el cilindro ===

En aerodinámica potencial, el esfuerzo o fuerza lateral que ejerce el fluido sobre un cuerpo se relaciona
directamente con la circulación mediante el teorema de Kutta–Joukowski:

<math>
F = -\rho_f U_\infty \Gamma.
</math>

Como en este caso <math>\Gamma = 0</math>, la fuerza resultante es:

<math>
F = 0.
</math>

Es decir, a pesar de que el fluido se desvía alrededor del cilindro, *no aparece fuerza neta sobre él*.

Este resultado es una manifestación de la *paradoja de D’Alembert*:
en un flujo potencial ideal, sin viscosidad y sin separación de la capa límite, la fuerza sobre un
obstáculo fijo es exactamente cero, lo cual contradice la experiencia real.

=== Interpretación física ===

•⁠ ⁠En el modelo idealizado de fluido potencial no viscoso, el flujo se ajusta perfectamente a la superficie del cilindro.
•⁠ ⁠No hay formación de estela ni separación de flujo.
•⁠ ⁠Las presiones alrededor del cilindro se distribuyen de manera simétrica.
•⁠ ⁠Esta simetría implica que las fuerzas horizontales se cancelan exactamente.

Sin embargo, en la realidad:

•⁠ ⁠La viscosidad provoca la separación de la capa límite.
•⁠ ⁠Aparece una estela turbulenta detrás del cilindro.
•⁠ ⁠La presión en la parte posterior es menor que en la delantera.
•⁠ ⁠Esto genera una fuerza de arrastre real, que el modelo potencial no es capaz de predecir.

Por ello, este fenómeno se denomina la *paradoja de D’Alembert*, ya que muestra las limitaciones del modelo de
flujo potencial para describir fuerzas sobre cuerpos sumergidos.

==apartado 9==
===.Nueva representación del Potencial y del campo de velocidades===
Repetimos el apartado 2 con el nuevo potencial:
<center><math> \varphi (\rho ,\theta, z)=(\rho +\frac{1}{\rho})\cos (\theta )+\frac{\theta}{4\pi} </math></center>
[[Archivo:Curvasnivel926.png|400px|miniaturadeimagen| Función Potencial 2]]
{{matlab|codigo=
% Función potencial y campo de velocidades para
% phi(rho,theta) = (rho + 1/rho) * cos(theta) + theta/4*pi

clear; clc; close all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 40; % puntos radiales
Nth = 120; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Campo de velocidades u = grad(phi)
% En polares: u_rho = dphi/drho, u_th = (1/rho) dphi/dth

u_rho = (1 - 1./RHO.^2) .* cos(TH);
% dphi/dtheta
dphi_dtheta = -(RHO + 1./RHO) .* sin(TH) + 1/(4*pi);

% u_theta = (1/r)*dphi/dtheta
u_th = (1./RHO) .* dphi_dtheta;

% Pasamos a componentes cartesianas:
u_x = u_rho .* cos(TH) - u_th .* sin(TH);
u_y = u_rho .* sin(TH) + u_th .* cos(TH);

% Puntos del contorno del obstáculo (r = 1)
th_circ = linspace(0, 2*pi, 400);
x_circ = R1 * cos(th_circ);
y_circ = R1 * sin(th_circ);

%Dibujamos las curvas de nivel del potencial
figure;
contour(X, Y, phi, 30); % 30 niveles de phi
hold on;
plot(x_circ, y_circ, 'k', 'LineWidth', 2); % cilindro
axis equal;
xlim([-4 4]); ylim([-4 4]);
xlabel('x'); ylabel('y');
title('Curvas de nivel de la función potencial \phi');
hold off;
}}
A partir de la función potencial, la velocidad del fluido se determina mediante su gradiente, <math>\vec{u}</math>
=∇φ.

<center><math>
\vec{u}(\rho,\theta,z)=\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\mathbf e_{\rho}+\frac{1}{\rho}\left[-\left(\rho+\frac{1}{\rho}\right)\sin\theta+\frac{1}{4\pi}\right]\mathbf e_{\theta}
</math></center>

Aquí representaremos el campo de velocidades resultante y analizaremos la dirección y magnitud del movimiento de las partículas del fluido.

El campo <math> \vec u </math> lo hemos pasado manualmente a coordenadas cartesianas con la matriz de cambio de base para añadirlo directamente a nuestro código de Matlab. Dándonos como resultado:
<center><math> \vec u (\vec i,\vec j,\vec k)=\left[\cos^{2}\theta\left(1-\frac{1}{\rho^{2}}\right)-\sin\theta\cos\theta+\frac{\sin\theta}{4\pi\rho}\right]\mathbf i+\left[\sin\theta\cos\theta\left(1-\frac{1}{\rho^{2}}\right)+\cos^{2}\theta+\frac{\cos\theta}{4\pi\rho}\right]\mathbf j </math></center>

[[Archivo:Campovelocidades926.png|400px|miniaturadeimagen| Campo de Velocidades resultante]]
[[Archivo:Campovelocidadesampl926.png|400px|miniaturadeimagen| Campo de Velocidades resultante]]
{{matlab|codigo=
clear; clc;clear all;

% Parámetros del dominio
R1 = 1; % radio del cilindro
R2 = 5; % radio exterior
Nr = 10; % puntos radiales
Nth = 70; % puntos angulares

rho = linspace(R1, R2, Nr);
theta = linspace(0, 2*pi, Nth);
[RHO, TH] = meshgrid(rho, theta);

% Coordenadas cartesianas
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Dibujamos las curvas de nivel
contour(X,Y,Z,15);
hold on

%Definimos las componentes X e Y del gradiente en coordenadas cartesianas
Gx=(cos(TH).^2).*(1 - 1./RHO.^2) - sin(TH).*cos(TH) + sin(TH)./(4*pi.*RHO);
Gy=sin(TH).*cos(TH).*(1 - 1./RHO.^2) + cos(TH).^2 + cos(TH)./(4*pi.*RHO);

%Dibujamos el campo de velocidades
quiver(X,Y,Gx,Gy);

% Representamos nuestro obstáculo
plot(1*cos(theta),1*sin(theta),'k','lineWidth',1);
axis([-4,4,-4,4]);
colorbar;
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}

===. Comprobación rotacional y divergencia nulos===

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

Las componentes de velocidad:

<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}
= \frac{1}{\rho}\left(-\left(\rho+\frac{1}{\rho}\right)\sin\theta+\frac{1}{4\pi}\right)
= -\left(1+\frac{1}{\rho^{2}}\right)\sin\theta + \frac{1}{4\pi\rho}.
</math>

El campo de velocidades es:

<math>
\vec{u}
= u_{\rho}\,\vec{e}_{\rho}+u_{\theta}\,\vec{e}_{\theta}
=\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\vec{e}_{\rho}
-\left(1+\frac{1}{\rho^{2}}\right)\sin\theta\,\vec{e}_{\theta}
+\frac{1}{4\pi\rho}\,\vec{e}_{\theta}.
</math>

Comprobación Rotacional nulo:

<math>
\nabla\times\vec{u}
=\frac{1}{\rho}
\begin{vmatrix}
\vec{e}_{\rho} & \rho\vec{e}_{\theta} & \vec{e}_{z} \\
\dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial z} \\
u_{\rho} & \rho u_{\theta} & 0
\end{vmatrix},
</math>

<math>
\nabla\times\vec{u}
=\frac{1}{\rho}
\begin{vmatrix}
\vec{e}_{\rho} & \rho\vec{e}_{\theta} & \vec{e}_{z} \\
\dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial z} \\
\left(1-\dfrac{1}{\rho^{2}}\right)\cos\theta
&
-\left(\rho+\dfrac{1}{\rho}\right)\sin\theta+\dfrac{1}{4\pi}
&
0
\end{vmatrix}
=
-\left(1-\dfrac{1}{\rho^{2}}\right)\sin\theta\,\vec{e}_{z}
+
\left(\dfrac{1}{\rho^{2}}\right)\sin\theta\,\vec{e}_{\theta}

</math>

Esta vez, nuestro rotacional no es nulo, por lo que el flujo no es irrotacional y las partículas de fluido giran localmente.

Comprobamos la divergencia nula:

<math>
\nabla\cdot\vec{u}
=\frac{1}{\rho}\frac{\partial}{\partial\rho}\bigl(\rho u_{\rho}\bigr)
+\frac{1}{\rho}\frac{\partial u_{\theta}}{\partial\theta}
+\frac{\partial u_{z}}{\partial z},
\qquad u_{z}=0.
</math>

<math>
\nabla\cdot\vec{u}
=\frac{1}{\rho}\left[
\frac{\partial}{\partial\rho}\left(\left(1-\frac{1}{\rho^{2}}\right)\cos\theta\,\rho\right)
+\frac{\partial}{\partial\theta}\left(-\left(1+\frac{1}{\rho^{2}}\right)\sin\theta+\frac{1}{4\pi\rho}\right)
+\frac{\partial}{\partial z}(0)
\right]
=\frac{1}{\rho}\left[
\left(1+\frac{1}{\rho^{2}}\right)\cos\theta]
</math>

Al igual, obtenemos una divergencia no nula, es decir, significa que el volumen del fluido varía

[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]