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.

7 Trayectorias de Particulas

Gran parte de la problemática y motivación de Ekman era establecer el flujo de agua resultante de todo el fenómeno a partir de una columna de agua. Para calcularlo, conviene imaginarse un plano genérico [math]S[/math] definido por el vector [math]\vec{n}_S=\cos\alpha\vec{i}+\sin\alpha\vec{j}[/math], para un ángulo fijo [math]\alpha \in [0, 2\pi)[/math]. Este plano sobre el que mediremos el flujo tendría una profundidad infinita ([math]z \rightarrow \infty [/math]), aunque solo tendríamos que considerarlo hasta la profundidad [math]d_E[/math], pues es hasta donde el efecto tiene sentido estudiarlo.

Asimismo, diremos por conveniencia que el plano tiene una anchura [math]L[/math], de dimensión muy inferior a la profundidad total del océano.

Para hallar el flujo se hace uso del cálculo integral ateniendo a las definiciones pertinentes:

[math] \Phi = \int_S \vec{v}·d\vec{S} = \int_S (\vec{v}·\vec{n}_S)\ dS = \iint_D (\vec{v}·\vec{n})\,dl\,dz [/math]

Definimos el dominio [math]D[/math] de integración a partir de las dimensiones en las que existe el plano (lado [math]l[/math] y profundidad [math]z[/math]). La profundidad vamos a considerarla entre [math]-\infty[/math] y 0 por fines educativos, mientras que el lado simplemente entre 0 y [math]L[/math].

Asimismo, el resultado del producto vectorial [math]\vec{v}·\vec{n}_S [/math] es:

[math]\vec{v}·\vec{n}_S= \cos \alpha \cos\left(\frac{z}{d_E}+\vartheta\right)+\sin \alpha \sin\left(\frac{z}{d_E}+\vartheta\right) [/math]

Convenientemente, a partir de la fórmula de la suma de ángulos para el coseno, intuimos que esto es igual a:

[math]\vec{v}·\vec{n}_S=\cos\left(\frac{z}{d_E}+\vartheta - \alpha\right)[/math]

Integramos cual integral doble:

[math] \Phi = \int_{0}^{L} V_0 ·\left(\int_{-\infty}^{0} e^{\frac{z}{d_E}} ·\cos\left(\frac{z}{d_E}+\vartheta - \alpha\right) dz \right) dl [/math]

La integral interior se resuelve por partes, resultando:

[math] \Phi = \int_{0}^{L}\frac{V_0 \cdot d_E \cdot e^{\frac{z}{d_E}} \left(\cos\left(\frac{z}{d_E} + \vartheta - \alpha\right) + \sin\left(\frac{z}{d_E} + \vartheta - \alpha\right)\right)}{2} \, dl = [/math]

[math] = \int_{0}^{L} \frac{V_0}{2}· d_E \cdot \left(\sin\left(\vartheta - \alpha\right) + \cos\left(\vartheta - \alpha\right)\right) \, dl = [/math]

[math] = L \cdot \frac{V_0}{2} \cdot d_E \cdot \left(\sin\left(\vartheta - \alpha\right) + \cos\left(\vartheta - \alpha\right)\right) [/math]

Para hallar el flujo resultante, es cuestión de maximizar la expresión arriba descrita según el ángulo [math]\alpha[/math]. Atendiendo a procesos de optimización, para [math]\vartheta = 3\pi/4[/math], resulta que el flujo es máximo para un plano con el ángulo [math]\alpha = 3\pi/4[/math]. Esto significa, tal y como hemos definido el plano, que el flujo resultante es hacia el oeste. Generalizado, se demuestra así que, con el fenómeno de Ekman, el flujo resultante es siempre perpendicular a la dirección del viento, siendo su sentido dependiente precisamente el parámetro de Coriolis y el hemisferio en el que nos encontremos.

8 Circulaón y Paradoja D' Alembert

Espiral de Ekman graficada en coordenadas cilíndricas para valores comprendidos entre la superficie y [math]d_E[/math]. Naturalmente, la representación es idéntica a cualquier otro sistema de coordenadas.

La espiral de Ekman, tal y como convenientemente su nombre indica, se trata de una espiral y por tanto tiene cierto carácter curvo. Para este tipo de construcciones, siempre es interesante conocer cómo se pueden observar en otro tipo de coordenadas. Es de particular interés para esta curva el planteamiento en coordenadas cilíndricas.

Tal y como hemos descrito la curva anteriormente, la espiral de Ekman definida por el campo de velocidades respecto a la profundidad queda como:

[math] \vec{v}(t) = (u(z), v(z), z) = V_0 · e^{\frac{z}{d_E}} \ \left(sgn(f)·\cos\left( \frac{z}{d_E}+\vartheta\right), \ \sin\left( \frac{z}{d_E}+\vartheta \right), z\right) [/math]

Para poder efectuar el cambio, partimos de las definiciones de las coordenadas cilíndricas:

[math] \begin{cases} \rho(z)=\sqrt{u(z)^2 +v(z)^2} \\ \theta(z)=\arctan\left(\frac{v(z)}{u(z)}\right)\\ z(z)=z \end{cases}\ [/math]

Una vez aplicado todo, obtenemos:

\begin{cases} \rho(z)=V_0·e^{\frac{z}{d_E}} \\ \theta(z)=\frac{z}{d_E}+\vartheta\\ z(z)=z \end{cases}

Destacar que, puesto que [math]\rho[/math] queda definido para valores positivos de la recta real, la función [math]sgn(f)[/math] es irrelevante. Ídem para [math]\theta[/math].

Para expresar todo en la base física cilíndrica podemos usar la matriz de cambio de coordenadas cartesianas a cilíndricas [math]M_{C\ \rightarrow \ \mathcal{C}}[/math]. Sencillamente:

[math] [\vec{v}]_\mathcal{C} = M_{C\ \rightarrow \ \mathcal{C}} · [\vec{v}]_C = \rho · \begin{pmatrix} \cos \theta & \sin \theta & 0\\ -\sin\theta & \cos\theta & 0\\ 0 & 0 & 1 \end{pmatrix} · \begin{pmatrix} sgn(f)·\cos\theta\\ \sin\theta\\ z \end{pmatrix} \Rightarrow [/math]

[math] \Rightarrow \vec{v}(\rho, \theta, z) = \rho\left( \cdot \cos^2 \theta + \sin^2 \theta \right)\vec{e_\rho} + \rho\left(-\sin \theta \cdot \cos \theta + \cos \theta \cdot \sin \theta\right)\vec{e_\theta} + z\vec{e_z} = \rho\vec{e_\rho} + z\vec{e_z}[/math],

Lo cual era de esperar atendiendo a la definición del vector posición en coordenadas cilíndricas [math]\vec{r}=\rho\vec{e_\rho} + z\vec{e_z} \quad [/math].

En MATLAB, se trabaja exclusivamente con coordenadas cartesianas, de manera que para poder representar la curva en coordenadas cilíndricas, habría que deshacer todos los cambios hechos. Se trata de un paso trivial que ya queda definido en el siguiente código:

• Código para la representación tridimensional de la espiral de Ekman a partir de su definición en coordenadas cilíndricas:

% Parámetros iniciales
dE = sqrt(2 * 0.1 / 10^-4); % Cálculo de dE
z = linspace(0, dE, 500); % Valores de z
v=0.2;

% Definición de la curva
%Valores de rho, theta y z
rho = v .* exp(z ./ dE);
theta = (z ./ dE) + 3*pi/4;
Z=z-dE;
%Paso a cilíndricas para graficar en MATLAB
X=rho.*cos(theta);
Y=rho.*sin(theta);
Z=z-dE;

figure;
hold on;
%axis([-0.4,0.4,-0.4,0.4,-150,0]);
axis tight;
grid on;
zlabel('Profundidad');  %Z

view(3);
plot3(X, Y, Z, 'r', 'LineWidth', 1);

9 Flujo con circulación impuesta

Curvatura [math]k(t)[/math] y torsión [math]\tau (t)[/math] graficadas para profundidades entre 0 y [math]d_E[/math].

En este tipo de curvas, conocer cuál es su torsión y curvatura resulta más que relevante, puesto que nos proporcionan importante información acerca de esta, su forma y su comportamiento. Para hallar cada una:

• Curvatura [math]k(t)[/math]: expresa cuánto se parece la curva en un punto a su

circunferencia osculatriz de radio [math]1/k(t)[/math], es decir, cuánto comparte la curva con la circunferencia. Para una curva [math]\gamma (t)[/math] se define de la siguiente manera:

[math] k(t) = \frac{|\dot{\gamma}(t) \times \ddot{\gamma}(t)|}{|\dot{\gamma}(t)|^3} [/math]

• Torsión [math]\tau (t)[/math]: expresa cuánto se separa la curva del plano osculador, es decir, de aquel dado por la circunferencia osculatriz o, lo que es lo mismo, dado por el vector tangente y normal en un punto. Igualmente, se define a partir de:
  
  

[math] \tau (t) = \frac{[\dot{\gamma}(t), \ddot{\gamma}(t), \dddot{\gamma}(t)]}{|\dot{\gamma}(t) \times \ddot{\gamma}(t)|^2} [/math]

Concretamente, graficarlas resulta muy útil.

9.1 Análisis mediante el triedro de Frenet

derecha

Una curiosa forma en la que podemos observar la verdadera naturaleza curva de la espiral de Ekman es mediante la visualización animada del triedro de Frenet. Gastar líneas en calcular el triedro es irrelevante, pero el código empleado para mostrar la animación mostrada a partir de la función frenet() es el siguiente:

• La función frenet()

function [T,N,B,k,t] = frenet(x,y,z),

if nargin == 2,
    z = zeros(size(x));
end

% CONVERT TO COLUMN VECTOR
x = x(:);
y = y(:);
z = z(:);

% SPEED OF CURVE
dx = gradient(x);
dy = gradient(y);
dz = gradient(z);
dr = [dx dy dz];

ddx = gradient(dx);
ddy = gradient(dy);
ddz = gradient(dz);
ddr = [ddx ddy ddz];

% TANGENT
T = dr./mag(dr,3);


% DERIVIATIVE OF TANGENT
dTx =  gradient(T(:,1));
dTy =  gradient(T(:,2));
dTz =  gradient(T(:,3));

dT = [dTx dTy dTz];

% NORMAL
N = dT./mag(dT,3);
% BINORMAL
B = cross(T,N);
% CURVATURE
% k = mag(dT,1);
k = mag(cross(dr,ddr),1)./((mag(dr,1)).^3);
% TORSION
t = dot(-B,N,2);

function N = mag(T,n),
% MAGNATUDE OF A VECTOR (Nx3)
%  M = mag(U)
N = sum(abs(T).^2,2).^(1/2);
d = find(N==0); 
N(d) = eps*ones(size(d));
N = N(:,ones(n,1));

• El programa de animación:

% Parámetros iniciales
dE = sqrt(2 * 0.1 / 10^-4); % Cálculo de dE
z = linspace(-dE, 0, 100); % Valores de z de 0 a 250 m

% Definición de la curva
x = 0.2 * exp(z / dE) .* cos(z / dE + 3/4 * pi)*100;
y = 0.2 * exp(z / dE) .* sin(z / dE + 3/4 * pi)*100;


% Cálculo del triedro de Frenet usando tu función
[T, N, B, k, t] = frenet(x, y, z); % Asume que la función frenet.m está en el mismo directorio
T = 10*T; N=10*N; B=10*B;
% Crear la animación
figure;
hold on;
zlabel('Profundidad');
axis equal;
axis([-0.3*100 0.1*100 -0.1*100 0.25*100 -50 0]);
view(3);

% Plot de la espiral
plot3(x, y, z, 'k', 'LineWidth', 1); % La curva 
quiver3(10, 0, -1, -12, 0, 0, 'k', 'LineWidth', 2, 'MaxHeadSize',0.751) % Vector del viento
hT = quiver3(0, 0, 0, 0, 0, 0, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 0.5); % Vector T
hN = quiver3(0, 0, 0, 0, 0, 0, 'g', 'LineWidth', 1.5, 'MaxHeadSize', 0.5); % Vector N
hB = quiver3(0, 0, 0, 0, 0, 0, 'b', 'LineWidth', 1.5, 'MaxHeadSize', 0.5); % Vector B

% Animación
for i =1:length(z)
    % Actualizar los vectores del triedro
    set(hT, 'XData', x(length(z)+1-i), 'YData', y(length(z)+1-i), 'ZData', z(length(z)+1-i), ...
        'UData', T(length(z)+1-i, 1), 'VData', T(length(z)+1-i, 2), 'WData', T(length(z)+1-i, 3));
    set(hN, 'XData', x(length(z)+1-i), 'YData', y(length(z)+1-i), 'ZData', z(length(z)+1-i), ...
        'UData', N(length(z)+1-i, 1), 'VData', N(length(z)+1-i, 2), 'WData', N(length(z)+1-i, 3));
    set(hB, 'XData', x(length(z)+1-i), 'YData', y(length(z)+1-i), 'ZData', z(length(z)+1-i), ...
        'UData', B(length(z)+1-i, 1), 'VData', B(length(z)+1-i, 2), 'WData', B(length(z)+1-i, 3));
     exportgraphics(gca,"Ekman11Frenet.gif", "Append",true)
    % Pausa para animación
    pause(0.02);
end

• 

  Escalera en el Vaticano siguiendo la forma de una espiral logarítmica.

• 

  Antena de doble polarización en forma de espiral logarítmica.

• 

  Las borrascas son fenómenos naturales en forma de espiral.

• 

  Las las galaxias también siguen la forma de una espiral logarítmica por naturaleza.