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

Trabajo realizado por estudiantes
Título	Flujo alrededor de un obstáculo circular (Grupo 29)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Rodrigo de Pelayo García García Sergio Resino Velayo Cayetano Gilabert Castejón Hugo Moreno Peregrina
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El estudio del flujo alrededor de un cilindro circular es uno de los problemas más clásicos y analizados en la mecánica de fluidos teórica y aplicada. Su importancia en ingeniería civil es notable, pues muchas estructuras presentan formas aproximadamente cilíndricas: pilas de puentes, torres, tirantes, pilotes o tuberías expuestas a corrientes o viento.

Comprender el comportamiento del flujo permite predecir fenómenos como la distribución de presiones alrededor del obstáculo, las fuerzas de arrastre y sustentación, el desprendimiento periódico de vórtices, los fenómenos de resonancia estructural y la erosión local en estructuras hidráulicas. El flujo alrededor de un cilindro combina varios conceptos clave de Teoría de Campos: los campos vectoriales, el gradiente, el rotacional, la divergencia y las líneas de corriente. A lo largo de este trabajo se desarrolla progresivamente la formulación matemática completa.

1 Región del fluido y mallado en coordenadas polares

Para comenzar, se genera una malla que representa los puntos pertenecientes a la región donde se encuentra el fluido, situada en el exterior del cilindro de radio unidad. Para describir esta zona se emplean coordenadas polares dentro del anillo comprendido entre los radios 1 y 5, con centro en el origen. La visualización final se mostrará en coordenadas cartesianas dentro del dominio [-4,4] x [-4,4], de modo que se aprecie correctamente el área de trabajo.

Mediante el siguiente código elaborado en MATLAB es posible representar dicha malla:

Mallado que representa los puntos interiores de la región que ocupa un fluido

% Parámetros del anillo
rho_inner = 1;      % radio interior (círculo unidad)
rho_outer = 5;      % radio exterior
n_rho = 40;         % número de curvas rho = cte (radiales concéntricas)
n_theta = 40;       % número de curvas theta = cte (radiales)

% Generar malla en coordenadas polares
rho = linspace(rho_inner, rho_outer, n_rho);
theta = linspace(0, 2*pi, n_theta);
[U,V] = meshgrid(rho,theta);

%Prarametrizar la superficie
hold on
X=U.*cos(V);
Y=U.*sin(V);
Z=0.*U;

%Representación del mallado y dibujo del círculo
mesh(X,Y,Z);
plot(1*cos(theta), 1*sin(theta), 'LineWidth',2);
%Delimitar los ejes
axis equal;
axis([-4,4,-4,4]);

% Marcas de eje 
box on;
grid off;

% Etiquetas y título
xlabel('x');
ylabel('y');
title('Mallado del anillo');

hold off

2 Función potencial y campo de velocidades

Una vez delimitada la región en la que vamos a trabajar, el siguiente paso consiste en analizar el movimiento del fluido. Para ello estudiaremos la velocidad asociada a la función potencial, recordando que la velocidad en este tipo de modelos se obtiene aplicando el gradiente a dicha función.

Antes de continuar, resulta útil recordar brevemente qué entendemos por campo escalar y campo vectorial. Un campo escalar definido sobre un dominio [math] D \subset \mathbb{R}^{3} [/math] es una función que asigna a cada punto de D un valor numérico. En cambio, un campo vectorial en un dominio [math] D \subset \mathbb{R}^{3} [/math] asocia a cada punto de D un veto, describiendo así tanto magnitud como direccional en cada posición del espacio.

2.1 Función potencial

La función potencial a representar es la siguiente:

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

Funcion Potencial

I = linspace(1,5,200);      % mallado en rho
J = linspace(0,2*pi,300);   % mallado en theta

[rho,theta] = meshgrid(I,J);

hold on

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

%  Funcion potencial correcta 
% phi(rho,theta) = (rho + 1/rho)*cos(theta)
phi = (rho + 1./rho) .* cos(theta);

% Contornos de la función de nivel 
contour(X, Y, phi, 50, 'LineWidth', 1)

% Círculo de radio 1
plot(cos(J), sin(J), 'LineWidth', 2);

view(2);
axis([-5 5 -5 5]);
colorbar;

title('Función Potencial');
xlabel('X');
ylabel('Y');

axis equal
hold off

2.2 Campo de velocidades

La velocidad de un campo escalar se representa a través del gradiente de esa función. Por eso, el campo de velocidades se expresa como [math]\vec u=\nabla\varphi[/math]. El gradiente escrito en coordenadas polares, queda:

[math]\vec u=\nabla\varphi=\frac{∂φ}{∂\rho}\vec{e}_\rho + \frac{1}{\rho}\frac{∂φ}{∂\theta}\vec{e}_\theta [/math]

Por lo tanto, se tiene que el campo de velocidades de la función potencial dada será:

[math]\vec u=((1-\frac{1}{\rho^2})cos(\theta)\vec{e}_\rho-(1+\frac{1}{\rho})sin(\theta)\vec{e}_\theta[/math]

Campo de velocidades

Zoom campo de velocidades

rho = linspace(1, 5, 40);          % Rango radial
th  = linspace(0, 2*pi, 80);       % Rango angular más fino para mejor detalle

[U, V] = meshgrid(rho, th);

% Conversión a coordenadas cartesianas
X = U .* cos(V);
Y = U .* sin(V);

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

%  Derivadas parciales en polares
dphi_dr  = (1 - 1./U.^2) .* cos(V);
dphi_dth = -(U + 1./U) .* sin(V);

% Gradiente transformado a coordenadas cartesianas 
Cx = cos(V).*dphi_dr - (sin(V)./U).*dphi_dth;
Cy = sin(V).*dphi_dr + (cos(V)./U).*dphi_dth;


% FIGURA COMPLETA

figure;
subplot(1,2,1)
contour(X, Y, phi, 30)
hold on
quiver(X, Y, Cx, Cy, 'r')
plot(cos(th), sin(th), 'k', 'LineWidth', 1)
axis equal axis([-5 5 -5 5])
title('Función potencial y campo de velocidades')
xlabel('X')
ylabel('Y')
colorbar
hold off

% ZOOM ORTOGONALIDAD

subplot(1,2,2)
contour(X, Y, phi, 30)
hold on
quiver(X, Y, Cx, Cy, 'r')
plot(cos(th), sin(th), 'k', 'LineWidth', 1)

% Ajuste del zoom 
axis equal
zoom_x = [-2 0];
zoom_y = [0 2];
axis([zoom_x zoom_y])
title('Zoom: ∇\phi ortogonal a curvas de nivel')
xlabel('X')
ylabel('Y')
hold off

3 Divergencia y Rotacional

3.1 Divergencia del campo de velocidades

Sea [math]\vec u=u_\rho\vec{e}_\rho+u_\theta\vec{e}_\theta+u_z\vec{e}_z[/math]

La divergencia de un campo viene dada por la siguiente fórmula:

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

La divergencia de un campo vectorial representa una magnitud escalar en la que se compara el flujo entrante y saliente en una superficie. Ahora calculamos la divergencia para nuestra función:

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

En conclusión, hemos obtenido que [math] {\nabla \cdot \vec u = 0} [/math] , por lo que sabemos que el fluido ni se expande ni se contrae, así que podemos garantizar que el fluido es incompresible.

3.2 Rotacional del campo de velocidades

Sea [math]\vec u=u_\rho\vec{e}_\rho+u_\theta\vec{e}_\theta+u_z\vec{e}_z[/math]

El rotacional de un campo se define como:

[math]\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\cdot\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ u_\rho & \rho\cdot\ u_\theta & {u_z} \end{vmatrix}[/math]

El rotacional indica tanto la dirección como la velocidad del giro del campo estudiado en cada punto. Muestra la tendencia de un campo a inducir rotación alrededor de un punto. Aplicando esto a nuestra función, obtenemos que:

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

Tras los cálculos, obtenemos que [math]{\nabla\times\bar{u}=0}[/math], por tanto las partículas del flujo no giran.

4 Líneas de Corriente

Las líneas de corriente de [math]\vec{u}[/math] determinan la trayectoria que siguen las partículas del fluido. Son curvas tangentes a los vectores velocidad de cada partícula. Para calcularlas, simplemente basta con hallar el campo ortogonal a [math]\vec{u}[/math], o sea, [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\\(1-\frac{1}{\rho ^2})cos(\theta)&-(1+\frac{1}{\rho^2})sin(\theta)&0\end{vmatrix}=[(1+\frac{1}{\rho^2})sin(\theta)]e_\rho + [(1-\frac{1}{\rho ^2})cos(\theta)]e_\theta[/math]

Sabiendo que la divergencia de [math]\vec{u}[/math] es nula, podemos decir que [math]\vec{v}[/math] es irrotacional. Además existe un potencial escalar [math]\psi [/math], cuyo gradiente es [math]\vec{v}[/math], que se conoce como función de corriente de [math]\vec{u}[/math]. Para calcularla, se resuelve el siguiente sistema de ecuaciones diferenciales ([math]V_z[/math] se puede despreciar ya que se está trabajando en el plano):

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

Tras despejar [math]\frac{1}{\rho }[/math] e integrar ambas ecuaciones, obtenemos la función [math]\psi=(\rho-\frac{1}{\rho})sin(\theta)[/math]

A continuación creamos una representación gráfica en la que comprobamos que los vectores que componen [math]\vec{v}[/math] son ortogonales a las líneas de corriente de [math]\psi[/math]. Este ha sido el código utilizado:

5 Puntos de Frontera

Analizamos la velocidad en la frontera es decir sobre la superficie del cilindro ([math]\rho = 1[/math]).

• La componente radial es nula: [math]u_\rho(1, \theta) = 0[/math].

• La componente tangencial es: [math]u_\theta(1, \theta) = -2\sin\theta[/math].

• El módulo de la velocidad es: [math]|\vec{u}| = 2|\sin\theta|[/math].

5.1 Puntos de Remanso (Velocidad mínima = 0):

Se dan cuando [math]\sin\theta = 0[/math].

• [math]\theta = 0 \rightarrow (x, y) = (1, 0)[/math]

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

En estos puntos, el fluido choca contra el obstáculo y se detiene momentáneamente.

5.2 Puntos de Velocidad Máxima (Velocidad = 2):

Se dan cuando [math]|\sin\theta| = 1[/math].

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

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

El fluido se acelera en los "hombros" del cilindro debido al estrechamiento de las líneas de corriente.

6 Presión según Bernouilli

Presiones y puntos de remanso y velocidad máxima

Partimos de la ecuación de Bernoulli para un fluido con densidad [math]\rho = 2[/math]:

[math]\frac{1}{2}\rho|\vec{u}|^2 + p = \text{cte}[/math]

Sustituyendo [math]\rho = 2[/math], la ecuación se simplifica:

[math]|\vec{u}|^2 + p = C \implies p = C - |\vec{u}|^2[/math]

Para el cálculo y la representación gráfica, elegiremos un valor arbitrario para la constante, por ejemplo C = 5. Así, la ecuación del campo de presiones es:

[math]p(r, \theta) = 5 - |\vec{u}(r, \theta)|^2[/math]

Análisis de Extremos:

Observando la ecuación, la presión y la velocidad tienen una relación inversa (si una crece, la otra disminuye):

• Presión Máxima: Se alcanza donde la velocidad es mínima ([math]|\vec{u}|=0[/math]).

  Coincide con los Puntos de Remanso en [math](1, 0)[/math] y [math](-1, 0)[/math].

  Valor máximo: [math]p_{max} = 5 - 0 = 5[/math].

• Presión Mínima: Se alcanza donde la velocidad es máxima ([math]|\vec{u}|=2[/math]).

  Ocurre en la parte superior [math](0, 1)[/math] e inferior [math](0, -1)[/math] del cilindro.

  Valor mínimo: [math]p_{min} = 5 - 2^2 = 1[/math].

Esto explica el fenómeno físico: el fluido pierde presión al acelerarse para rodear el obstáculo (efecto Venturi) y recupera presión al frenarse en la zona de estancamiento.

9 Flujo con circulación impuesta

I = linspace(1,5,200);          % mallado en rho
J = linspace(0,2*pi,300);       % mallado en theta

[rho,theta] = meshgrid(I,J);

hold on

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

% --- Función potencial corregida ---
% phi(rho,theta) = (rho + 1/rho)*cos(theta) + theta/(4*pi)
phi = (rho + 1./rho).*cos(theta) + theta/(4*pi);

% Contornos de la función de nivel
contour(X, Y, phi, 50, 'LineWidth', 1)

% Círculo de radio 1
plot(cos(J), sin(J), 'LineWidth', 2);

view(2);
axis([-5 5 -5 5]);
colorbar;

title('Función Potencial');
xlabel('X');
ylabel('Y');

axis equal
hold off