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^2})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 = -(1 + 1./U^2) .* 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

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:

Líneas de corriente, velocidad y vector u

Líneas de corriente, velocidad y vector u con zoom

% Definición del mallado
nr = 40; nt = 80;
r = linspace(1, 4, nr);
theta = linspace(0, 2*pi, nt);
[R, T] = meshgrid(r, theta);

X = R .* cos(T);
Y = R .* sin(T);

% Cálculo función de corriente(Psi)
Psi = (R - 1./R) .* sin(T);

% Cálculo vector velocidad (u)
% Componentes en Polares
Ur = (1 - 1./R.^2) .* cos(T);
Ut = -(1 + 1./R.^2) .* sin(T);

% Convertimos a Cartesianas 
Ux = Ur .* cos(T) - Ut .* sin(T);
Uy = Ur .* sin(T) + Ut .* cos(T);
Vx = -Uy; 
Vy = Ux;

% Representación gráfica
figure(4); clf;
set(gcf, 'Color', 'w');

[C, h] = contour(X, Y, Psi, 50, 'LineWidth', 1); 
hold on;

% Superponemos el campo de velocidades (u)
quiver(X , Y , Ux, Uy, 'LineWidth', 0.45 );

quiver(X, Y, Vx, Vy,'color','r', 'LineWidth', 0.45)

% Dibujamos el obstáculo
plot(cos(theta), sin(theta), 'LineWidth', 2);

axis equal; axis([-3.5 3.5 -3 3]);
%si le ponemos zoom: axis([-1.5 -0.5 -0 1])

title('Apartado 4: Líneas de Corriente (\psi = cte) y Velocidad (\vec{u})');
xlabel('x'); ylabel('y');

fprintf('Observa cómo las flechas negras (\vec{u}) son siempre TANGENTES a las líneas azules (\psi).\n');

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

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(\rho, \theta) = 5 - |\vec{u}(\rho, \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.

Presiones y puntos de remanso y velocidad máxima

%% Parámetros
U = 1;              % Velocidad del flujo uniforme
a = 1;              % Radio del cilindro
N = 400;            % Resolución del mallado

rho = 1;            % Densidad del fluido (puedes cambiarla)

%% Malla
x = linspace(-3,3,N);
y = linspace(-3,3,N);
[X,Y] = meshgrid(x,y);

R = sqrt(X.^2 + Y.^2);
Theta = atan2(Y,X);

%% Campo de velocidades potencial
ur = U .* cos(Theta) .* (1 - (a^2 ./ R.^2));
ut = -U .* sin(Theta) .* (1 + (a^2 ./ R.^2));

Ux = ur .* cos(Theta) - ut .* sin(Theta);
Uy = ur .* sin(Theta) + ut .* cos(Theta);

% Magnitud
V = sqrt(Ux.^2 + Uy.^2);

% En el interior del cilindro no hay flujo
mask = (R < a);
V(mask) = NaN;

%% Cálculo del coeficiente de presión
% Bernoulli ideal:
% C_p = 1 - (V/U)^2
Cp = 1 - (V./U).^2;

%% Presión absoluta relativa (solo escala)
p = rho * U^2 * Cp / 2;
p(mask) = NaN;

%% Gráfica
figure("Color","w");
contourf(X, Y, Cp, 60, "LineStyle", "none");
colorbar
hold on

% Contorno del cilindro
th = linspace(0,2*pi,300);
plot(a*cos(th), a*sin(th), 'k', 'LineWidth', 2);

axis equal
xlabel("x");
ylabel("y");
title("Campo de presiones (Coeficiente C_p) alrededor de un cilindro");

7 Trayectorias de Partículas

Observando las representaciones de las líneas de corriente, la distribución de presiones y el módulo de la velocidad (gráficas correspondientes a los apartados 2.1, 6 y 4), se aprecia que conforme el flujo se acerca al obstáculo circular, tanto la presión como la velocidad experimentan cambios significativos.

• 
• 
• 

Al observar la gráfica de las líneas de corriente, se puede notar un estrechamiento al rodear el cuerpo. Debido a que el área de paso se reduce, el fluido debe acelerar para mantener el caudal constante. Mientras que en la zonas donde las líneas de corriente están más espaciadas, el área de paso es mayor y por ello la velocidad disminuye.

Utilizando el principio de Bernoulli ([math]\frac{1}{2}\rho|\vec{u}|^2 + p = \text{cte}[/math] ) explicado en el apartado anterior, un aumento de velocidad supone una disminución de presión. Por lo tanto, al rodear un cuerpo circular, la velocidad aumenta y la presión disminuye. Por el contrario, en el punto de impacto frontal se produce una velocidad nula, provocando una presión máxima. Los puntos de velocidad máxima y los puntos de velocidad mínima (puntos de remanso) se ven representados en la grafica de presiones del apartado 6.

8 Circulación y Paradoja de D' Alembert

A continuación, vamos a comprobar que la circulación del campo [math]\vec{u}[/math] alrededor del obstáculo circular se anula. Aplicando el teorema de Kutta-Joukowski, el cual nos dice que la fuerza que ejerce el fluido sobre el obstáculo es proporcional a esta circulación, y sabiendo que es nula, podemos concluir que el fluido no ejerce ninguna fuerza sobre el obstáculo. Esto es lo que se conoce como la Paradoja de D'Alembert. Mediante la resolución de la siguiente integral, demostramos este suceso.

Sean: p = presión particularizada en [math]\rho=1[/math], [math]\vec n[/math] = vector perpendicular a la curva

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

Puesto que [math]p\cdot\vec{n}[/math] es la fuerza que ejerce el fluido en cada punto de la frontera, esto significa que, al sumar la proyección de todas estas fuerzas sobre la dirección [math]\vec{i}[/math] la resultante es nula.

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

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

[math]\vec {i} = 1\cdot\vec {i}[/math] ; para hallar el valor de este vector en cilíndricas se usará la matriz de cambio de base previamente utilizada pero traspuesta: [math] \begin{pmatrix} cos(\theta) & sin(\theta) & 0 \\ -sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = cos(\theta)\vec {e}_\rho-sin(\theta)\vec {e}_\theta[/math]

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

[math]\int_{0}^{2\pi}p·\vec{n}·\vec{i}d\theta=\int_{0}^{2\pi}(5 - 4sin^2(\theta))\cdot(-cos(\theta))d\theta=-[\frac{1}{3}(-4sin^3(\theta)+15sin(\theta))]_{0}^{2\pi}=0[/math]

Finalmente, queda demostrado que la resultante de las fuerzas que ejerce el flujo en dirección horizontal es nula, cumpliendo así la Paradoja de D'Alembert.

9 Flujo cambiando la función

9.1 Función potencial y campo de velocidades

A continuación, observáremos como la función

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

Influye en:

9.1.1 El flujo potencial

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

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

9.1.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{\partial \varphi}{\partial \rho}\vec{e}_\rho + \frac{1}{\rho}\frac{\partial \varphi}{\partial \theta}\vec{e}_\theta[/math]

Para este caso, la función potencial incluye un término de circulación, definida como:

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

Calculando las derivadas parciales correspondientes:

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

Por lo tanto, sustituyendo en la definición del gradiente (recordando dividir el término angular por [math]\rho[/math]), se tiene que el campo de velocidades será:

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

Funcion Potencial

Funcion Potencial

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) + V/(4*pi); 

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

% 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), 'LineWidth', 2)
axis equal; axis([-5 5 -5 5])
title('Función potencial y campo de velocidades')
xlabel('X')
ylabel('Y')
colorbar
hold off

% Zoom

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

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

9.2 Divergencia y Rotacional

9.2.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 en coordenadas cilíndricas:

[math]\nabla \cdot \vec{u} = \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho}(\rho \cdot u_\rho) + \frac{\partial}{\partial \theta} u_\theta + \frac{\partial}{\partial z}(\rho \cdot u_z) \right][/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 con circulación:

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

Observamos que la derivada respecto a [math]\theta[/math] del término de circulación ([math]\frac{1}{4\pi\rho}[/math]) es cero.

[math]\nabla \cdot \vec{u} = \frac{1}{\rho}\left[ \left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right] + \frac{1}{\rho}\left[ -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right] = 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.

9.2.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, donde [math]\rho u_\theta = -(\rho + 1/\rho)\sin(\theta) + 1/4\pi[/math], 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) + \frac{1}{4\pi} & 0 \end{vmatrix}[/math]

Calculando la componente en [math]\vec{e}_z[/math] (única no nula en 2D):

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

Tras los cálculos, obtenemos que [math]\nabla \times \vec{u} = 0[/math], por tanto las partículas del fluido no giran (es irrotacional), a pesar de la existencia de circulación global alrededor del cilindro.

9.3 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} = \vec{k} \times \vec{u}[/math]:

[math]\vec{v} = \vec{k} \times \vec{u} = \begin{vmatrix} \vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\ 0 & 0 & 1 \\ (1 - \frac{1}{\rho^2})\cos(\theta) & -\left[(1 + \frac{1}{\rho^2})\sin(\theta) - \frac{1}{4\pi\rho}\right] & 0 \end{vmatrix}[/math] [math]= \left[ \left(1 + \frac{1}{\rho^2}\right)\sin(\theta) - \frac{1}{4\pi\rho} \right]\vec{e}_\rho + \left[ \left(1 - \frac{1}{\rho^2}\right)\cos(\theta) \right]\vec{e}_\theta[/math]

Sabiendo que la divergencia de [math]\vec{u}[/math] es nula, podemos decir que [math]\vec{v}[/math] es irrotacional (y viceversa). 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}, \quad V_\theta = \frac{1}{\rho}\frac{\partial \psi}{\partial \theta}[/math]

Al integrar [math]V_\rho[/math] respecto a [math]\rho[/math], el término [math]1/\rho[/math] da lugar a un logaritmo natural. Tras despejar e integrar ambas ecuaciones, obtenemos la función [math]\psi[/math]:

[math]\psi = \left(\rho - \frac{1}{\rho}\right)\sin(\theta) - \frac{1}{4\pi}\ln(\rho)[/math]

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

% Configuración del Modelo
nr = 60; nt = 100;              % Aumentamos la resolución para curvas suaves
r = linspace(1, 4.5, nr);       % Radio desde 1 hasta 4.5
theta = linspace(0, 2*pi, nt);  % 0 a 2pi
[R, T] = meshgrid(r, theta);

% Coordenadas Cartesianas
X = R .* cos(T);
Y = R .* sin(T);

% Cálculos Físicos
% 1. Función de Corriente (Psi)
Psi = (R - 1./R) .* sin(T) - (1/(4*pi)) * log(R);

% 2. Campo de Velocidades (u)
Ur = (1 - 1./R.^2) .* cos(T);
Ut = -(1 + 1./R.^2) .* sin(T) + 1./(4*pi*R);

% Conversión de vectores a Cartesianas
Ux = Ur .* cos(T) - Ut .* sin(T);
Uy = Ur .* sin(T) + Ut .* cos(T);

% GENERACIÓN DE GRÁFICOS
figure('Name', 'Líneas de Corriente y Velocidad', 'Color', 'w');

% GRÁFICA 1: VISTA NORMAL
subplot(1, 2, 1)

% Dibujamos líneas de corriente
contour(X, Y, Psi, 40, 'LineWidth', 1.0); 
hold on

% Dibujamos vectores de velocidad (u)
q1 = quiver(X, Y, Ux, Uy, 'Color', [0.2 0.2 0.2], 'LineWidth', 0.8);

% Dibujamos el cilindro
plot(cos(theta), sin(theta), 'k', 'LineWidth', 2.5);
axis equal
axis([-3.5 3.5 -3 3])
title('Vista General: Líneas de Corriente \psi')
xlabel('x'); ylabel('y');
hold off

% GRÁFICA 2: ZOOM
subplot(1, 2, 2)

% Dibujamos líneas de corriente
contour(X, Y, Psi, 60, 'LineWidth', 1.2);
hold on

% Dibujamos vectores de velocidad
q2 = quiver(X, Y, Ux, Uy, 'Color', 'r', 'LineWidth', 1.2, 'AutoScaleFactor', 0.6);

% Dibujamos el cilindro
plot(cos(theta), sin(theta), 'k', 'LineWidth', 3);
axis equal
axis([-2 0.5 0 2]) 
title('Zoom: Detalle cerca de la superficie')
xlabel('x'); ylabel('y');
hold off