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.

5 Puntos de Remanso

Es importante marcar que la divergencia del campo vectorial [math]\vec{v}[/math] es en todo caso nula. La demostración es trivial, pues el campo depende exclusivamente de la variable [math]z[/math] mientras que todos los vectores resultantes son paralelos al plano marcado por [math]\vec{i}[/math] y [math]\vec{j}[/math]. En otras palabras, esto resulta que en que la derivada [math]\frac{\partial}{\partial z}\vec{v}_z = 0[/math], pues [math]\vec{v}_z[/math] es en sí ya nulo.

En el contexto en el que nos encontramos, no está de más preguntarse la razón de esta nulidad de la divergencia. Si entendemos la divergencia como un fenómeno físico, es correspondiente al cambio de volumen inducido por un campo. Siendo el agua el fluido incompresible que es, la velocidad de sus moléculas adquirida gracias al efecto del viento es totalmente independiente del volumen. Para todo campo velocidad de un fluido incompresible, sería por tanto de esperar que [math]\nabla · \vec{v} = 0 [/math]. Se trata de una consecuencia directa del Teorema de la divergencia de Gauss.

6 Presión según Bernouilli

Animación de la espiral de Ekman (curva roja) vista desde el cenit y definida a partir del vector velocidad [math]\vec{v}[/math] (flecha de color cambiante) con su rotacional [math]\nabla\times\vec{v}[/math] (flecha rosa).

En línea con el análisis de curvas, el rotacional de la espiral de Ekman, sirve para comprender la condición de espiral de la que goza. A partir del campo de velocidades [math]\vec{v}= u \vec{i} + v \vec{j}[/math], se da que:

[math] \nabla\times\vec{v}(z) = {\begin{vmatrix} \vec i & \vec j & \vec k \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ \vec v_x & \vec v_y & \vec v_z \end{vmatrix}} [/math]

Siendo [math]\vec v_x[/math], [math]\vec v_y[/math] y [math]\vec v_z[/math] [math]u[/math], [math]v[/math] y 0 respectivamente. Resolviendo el determinante, el proceso queda simplificado a:

[math] \nabla\times\vec{v}(z) = -\frac{\partial v}{\partial z}\vec{i} + \frac{\partial u}{\partial z}\vec{j} [/math]

Por simplicidad, llamamos [math]\varphi[/math] a la fase angular dentro de las funciones trigonométricas ([math]z/d_E + \vartheta[/math]). Sustituyendo:

[math] \nabla\times\vec{v}(z) = \frac{V_0}{d_E}·e^{\frac{z}{d_E}}·\left[-\left(\sin\varphi +\cos\varphi\right)\vec{i} + \left(\cos\varphi - \sin\varphi\right)\vec{j}\right] [/math]

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.