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

Trabajo realizado por estudiantes
Título	Flujo alrededor de un obstáculo circular (Grupo 26)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Gonzalo Gallego Fulgencio Andrea García Carrasco Aarón García Martín Miryam Sánchez-Ferragut Samalea Guillermo Rodríguez Navadijos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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.

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

Mallado que representa los puntos de la región ocupada por un fluido

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

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

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

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

Mallado que representa los puntos de la región ocupada por un fluido

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

2.2 . 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] =∇φ.

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

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:

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

Mallado que representa los puntos de la región ocupada por un fluido

Mallado que representa los puntos de la región ocupada por un fluido

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

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

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

3.1 . Comprobación del rotacional nulo

Conociendo la fórmula del rotacional calculamos:

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

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

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.

3.2 . Comprobación de la divergencia nula

Conociendo la fórmula de la divergencia calculamos:

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

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

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

4 . 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]).

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

Comprobamos que [math]\vec{v}[/math] es irrotacional:

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

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

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

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

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

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

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

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

9 apartado 9

Repetimos el apartado 2 con el nuevo potencial:

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

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

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(1-\dfrac{1}{\rho^{2}}\right)\sin\theta\,\vec{e}_{z} =0 [/math]

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:

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

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.