Diferencia entre revisiones de «Espiral de Ekman (Grupo 55, Retiro)»

Revisión del 21:48 1 dic 2025

1 Resolucion de [math]\vartheta[/math]

El valor de [math]\vartheta[/math] es importante para la definición de la espiral de Ekman, ya que determina la fase inicial fijada por la dirección del viento respecto a la fuerza de Coriolis.

Considerando una localidad de las Islas Canarias, con coordenadas aproximadas 30º10′24.2″N, 15º30′26.5″W, asumiendo una viscosidad turbulenta de 0.05 m²/s, una velocidad del viento de 12 m/s que sopla de Norte a Sur, una velocidad superficial inducida de aproximadamente 0.15 m/s y un desvío aproximado de 45º (hacia la derecha respecto a la dirección del viento) del flujo supercial. Usamos los siguientes valores:

• [math]z = 0[/math].
• [math]\frac{v(z)}{u(z)} = 1 = \tan(45º) [/math].
• [math]f \gt 0[/math] ya que estamos en el norte, y por lo tanto [math]sgn(f) = 1[/math].

Resolviendo nos queda:

[math] \frac{v(z)}{u(z)} = \frac{V_0 e^{\frac{0}{d_E}} \sin(\vartheta)} {V_0 e^{\frac{0}{d_E}} \cos(\vartheta)} = \tan(\vartheta) = 1 \Rightarrow [/math]

[math] \Rightarrow \vartheta = \left\{ \begin{aligned} \frac{\pi}{4}\\ \frac{3\pi}{4} \end{aligned} \right. [/math]

Como tanto [math]u(z)[/math] como [math]v(z)[/math] son negativos, la solución correcta para [math]\vartheta[/math] es:

2 Campo vectorial v

Este apartado, analiza la representación de un campo vectorial V evaluado en puntos de coordenadas cartesianas (0,0,𝑧)(a lo largo del eje vertical z), desde la superficie del mar hasta la profundidad de Ekman. La finalidad es visualizar la formación de la espiral de Ekman, que describe el patrón de circulación de las corrientes oceánicas debido a la acción del viento. A medida que se desciende en la columna de agua, el campo de velocidad cambia de dirección y magnitud, creando un perfil espiral en el cual los vectores de velocidad se representan en distintas capas de agua a lo largo del eje 𝑧. En este contexto, cada vector mostrado representa la dirección y la intensidad del flujo en una capa particular, permitiendo observar la transición desde las corrientes superficiales hasta las más profundas que componen la espiral de Ekman.

%% Apartado 4

clc;close all;clear;
omega = 7.2921*(10)^-5;  % Velocidad angular de la Tierra (rad/s)
phi = 45;  % Latitud en grados
V0 = 0.15;  % Velocidad superficial inducida por el viento (m/s)
nu = 0.05;  % Viscosidad turbulenta (m^2/s)
theta0 = 3*pi/4;  % Ángulo inicial en radianes

% Calcular el parámetro de Coriolis f
f = 2 * omega * sind(phi);

% Calcular la profundidad de Ekman
delta = sqrt(2 * nu / abs(f));

% Parámetros de la simulación
z_max = 5 * delta;  % Profundidad máxima para la animación (2 veces la profundidad de Ekman)
n_frames = 100;  % Número de valores de profundidad
z_vals = linspace(0, z_max, n_frames);  % Valores de profundidad
figure(1);
hold on;
view(3)
xlim([-0.3, 0.3]);
ylim([-0.3, 0.3]);
zlim([-z_max 0]);
xlabel('Este - Oeste (m)');
ylabel('Norte - Sur (m)');
u_m = (f/abs(f))*V0 * exp(-z_vals / delta)  .* cos(-z_vals / delta+theta0);  % Componente u(z)
v_m = V0 * exp(-z_vals / delta)  .* sin(-z_vals / delta+theta0);  % Componente v(z)
plot3(u_m,v_m,-z_vals,'Color','r')
hold on
% Inicializar la animación
for k = 1:n_frames
    z = z_vals(k);  % Profundidad actual
    % Calcular las componentes de la velocidad
    u = u_m(k); % Componente u(z)
    v = v_m(k);% Componente v(z)
    quiver3(u, v,-z, u, v,0, 'AutoScale', 'off', 'LineWidth', 0.5, 'Color', "b",'HandleVisibility', 'off');
    title(sprintf('Campo Vectorial de Ekman'))
    grid on
    view([45 45])
end

3 Visualización tridimensional de la espiral de Ekman

Representación tridimensional de los vectores [math]\vec V(z)[/math] evaluados en [math](0,0,z)[/math] para [math]z \in [0,3d_E][/math]. La curva roja une las puntas de los vectores y forma la espiral de Ekman.

5.1. Vectores de velocidad a lo largo del eje vertical

El campo de velocidad horizontal del agua viene dado por

[math]\vec V(z) = u(z)\,\vec i + v(z)\,\vec j,[/math]

donde las soluciones de las ecuaciones de Ekman (obtenidas en el apartado 3) son

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

siendo:

• [math]V_0 = 0{,}15\ \text{m/s}[/math] la velocidad superficial inducida por el viento (dato del modelo);
• [math]d_E[/math] la profundidad de Ekman calculada en el apartado (1);
• [math]\vartheta[/math] la fase inicial encontrada en el apartado (2);
• [math]\operatorname{sgn}(f)=1[/math] porque [math]f\gt0[/math] en el hemisferio norte.

Para visualizar la espiral tridimensional seguimos exactamente lo que indica el enunciado. Consideramos puntos del eje vertical [math](0,0,z)[/math] con

[math]z \in [0,-3d_E].[/math]

Tomamos entre 30 y 40 puntos igualmente espaciados (por ejemplo, [math]N=40[/math]). En cada profundidad [math]z[/math] calculamos [math]\vec V(z)[/math] y dibujamos un vector con origen en [math](0,0,z)[/math] y punta en [math](u(z),v(z),z)[/math]. Las puntas de todos esos vectores describen una curva en el espacio tridimensional: la espiral de Ekman.

5.2. Código de MATLAB para la visualización 3D

El siguiente script de MATLAB implementa la construcción descrita anteriormente, utilizando los valores de [math]f[/math], [math]d_E[/math] y [math]\vartheta[/math] calculados en los apartados (1) y (2). Se generan entre 30 y 40 vectores de velocidad a lo largo de la columna de agua y se representa la espiral de Ekman en tres dimensiones.

%% Datos conocidos de los apartados (1) y (2)
V0     = 0.15;          % Velocidad superficial inducida [m/s]
f      = 7.3e-5;        % Parámetro de Coriolis [1/s] (valor hallado en (1))
nu_e   = 0.05;          % Viscosidad turbulenta [m^2/s]
dE     = sqrt(2*nu_e/abs(f));   % Profundidad de Ekman (calculada en (1))
theta0 = 3*pi/4;        % Fase inicial ϑ (hallada en (2))

%% Discretización de la profundidad física: prof ∈ [0, 3 dE]
N    = 40;                          % 30–40 puntos
prof = linspace(0, 3*dE, N);        % prof ≥ 0, 0 = superficie

%% z del MODELO (negativo), tal como se definió en el enunciado
z_model = -prof;                    % z ≤ 0

%% Cálculo de u(z) y v(z) usando z_model
sgnf = sign(f);                     % = 1 en nuestro caso (hemisferio norte)

u = sgnf * V0 .* exp(z_model/dE) .* cos(z_model/dE + theta0);
v =        V0 .* exp(z_model/dE) .* sin(z_model/dE + theta0);

%% Representación 3D de los vectores V(z)
figure;

% Vectores con origen en (0,0,prof) y punta en (u,v,prof)
h = quiver3( zeros(size(prof)), zeros(size(prof)), prof, ...
             u,                v,                zeros(size(prof)), ...
             'LineWidth', 1.5 );
set(h,'AutoScale','on','AutoScaleFactor',5);   % Agranda las flechas para verlas bien

hold on;

% Curva que une las puntas de los vectores: la espiral de Ekman
plot3(u, v, prof, '--k', 'LineWidth', 1);

xlabel('Componente este u(z) [m/s]');
ylabel('Componente norte v(z) [m/s]');
zlabel('Profundidad z [m]');
title('Espiral de Ekman en 3D (vectores V evaluados en (0,0,z))');

grid on;
set(gca, 'ZDir','reverse');   % Profundidad positiva hacia abajo
view(45, 25);                 % Vista tridimensional agradable
% No usamos axis equal para que las flechas no queden aplastadas

hold off;

5.3. Interpretación del resultado numérico

La figura obtenida con el código anterior muestra la espiral de Ekman correspondiente a los parámetros del modelo. En la superficie ([math]z = 0[/math]) el vector de velocidad tiene módulo [math]V_0 = 0{,}15\,\text{m/s}[/math] y dirección fijada por la fase [math]\vartheta[/math], desviada aproximadamente [math]45^\circ[/math] a la derecha del viento, tal como se dedujo en el apartado (2).

A medida que aumenta la profundidad (valores crecientes de la variable [math]\text{prof}[/math]):

• la **dirección** de los vectores va rotando debido al efecto de Coriolis, describiendo una espiral;
• la **magnitud** de la velocidad decrece exponencialmente con el término [math]e^{z/d_E}[/math], por lo que los vectores se hacen cada vez más cortos.

La curva que une las puntas de los vectores [math](u(z),v(z),z)[/math] representa así la estructura tridimensional de la espiral de Ekman, y muestra cómo la influencia directa del viento se concentra en la capa superior de espesor del orden de la profundidad de Ekman [math]d_E[/math], desapareciendo progresivamente a mayor profundidad.

4 Divergencia de [math]\vec{v}[/math]

Es importante señalar que la divergencia del campo vectorial [math]\vec{v}[/math] es en todo caso nula. La justificación matemática es sencilla, pues el campo depende exclusivamente de la variable [math]z[/math], y todos sus vectores se orientan de forma paralela al plano generado por [math]\vec{i}[/math] y [math]\vec{j}[/math]. En otras palabras, esto implica que la derivada [math]\frac{\partial}{\partial z}\vec{v}_z = 0[/math], ya que la componente [math]\vec{v}_z[/math] es, de por sí, nula.

En el contexto en el que nos encontramos, no está de más preguntarse la causa de esta divergencia cero. 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.

5 Flujo Neto

El objetivo de este apartado es determinar el flujo neto de agua que atraviesa una pared vertical y demostrar que dicho flujo se orienta hacia el oeste, perpendicular al viento predominante. Para ello, consideramos una pared de agua cuya normal está dada por el vector genérico [math] \vec{n} = \cos(\alpha)\,\vec{i} + \sin(\alpha)\,\vec{j}[/math],donde [math] \alpha \in [0,2\pi) [/math] es un ángulo fijo. La pared tendrá una anchura [math]L = 10\ \text{m}[/math] y la profundidad será infinita [math] z \in [0,-\infty)[/math]. Para hallar el flujo usaremos la fórmula del flujo conociendo el campo vectorial [math]\vec{v} = u(z)\,\vec{i} + v(z)\,\vec{j}[/math],y un vector perpendicular a la superficie [math]\vec{n} = \cos(\alpha)\,\vec{i} + \sin(\alpha)\,\vec{j}[/math]: [math]\Phi = \int_0^{L} \int_{-\infty}^{0} (\vec{v} \cdot \vec{n}) \, dz \, dL[/math]

El resultado del producto vectorial[math](\vec{v} \cdot \vec{n})[/math]es:

[math]\vec{v}\cdot\vec{n} = u(z)\cos(\alpha) + v(z)\sin(\alpha)[/math]

[math]\vec{v}\cdot\vec{n} = V_0 e^{z/d_E} \left[ \cos(\alpha)\cos\left( \frac{z}{d_E} + \vartheta \right) + \sin(\alpha)\sin\left( \frac{z}{d_E} + \vartheta \right) \right][/math]

[math]\vec{v}\cdot\vec{n} = V_0 e^{z/d_E} \cos\left( \frac{z}{d_E} + \vartheta - \alpha \right).[/math]

Ahora la integral doble sería:

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

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

Para resolver lo que queda de integral usaremos la integración por partes:

[math] \int u \, dv = uv - \int v \, du [/math]

[math] u = \cos\left(\frac{z}{d_E} + ϑ - \alpha\right), \quad \implies \quad du = -\frac{1}{d_E} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) dz [/math]

[math]dv = e^{\frac{z}{d_E}} dz \implies \quad v = d_E e^{\frac{z}{d_E}} [/math]

Aplicamos la integración por partes: [math] \Phi = L V_0 \left[ d_E e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) - \int d_E e^{\frac{z}{d_E}} \left(-\frac{1}{d_E} \sin\left(\frac{z}{d_E} + ϑ -\alpha\right)\right) dz \right] [/math]

Si simplificamos: [math] \Phi = L V_0 \left[ d_E e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) + \int e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) dz \right] [/math]

Ahora nos queda otra integral:

[math] I_2 = \int e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) dz [/math]

Utilizando integral por partes otra vez [math] I_2 = d_E e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) - \int d_E e^{\frac{z}{d_E}} \frac{1}{d_E} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) dz[/math]

[math] I_2 = d_E e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) - \int e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) dz [/math]

Como la ultima integral es igual a la primera:

[math] I_2 = d_E e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right)-I [/math]

Esta la sustituimos en [math]\Phi [/math]

[math] \Phi = L V_0 \left[ d_E e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) + d_E e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) - I \right] [/math]

Como [math] \Phi=I = \int_{-\infty}^0 e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) dz [/math]

Entonces [math] 2I =V_{0} d_E e^{\frac{z}{d_E}} \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) + d_E e^{\frac{z}{d_E}} \sin\left(\frac{z}{d_E} + ϑ - \alpha\right)[/math]

Solo nos queda evaluar la integral entre [math] z=0 [/math] y [math] z=-\infty[/math]

[math] z =0:[/math]

[math]e^{\frac{z}{d_E}} = 1, \quad \cos\left(\frac{z}{d_E} + ϑ - \alpha\right) = \cos(ϑ - \alpha), \quad \sin\left(\frac{z}{d_E} + ϑ - \alpha\right) = \sin(ϑ - \alpha) [/math]

Y para [math] z =-\infty [/math] todos los [math] e^{\frac{z}{d_E}} \to 0 [/math] entonces se anulan

Entonces el flujo nos quedaría:

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

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

6 Curvatura y torsión para la espiral de Ekman

La curvatura y la torsión son elementos de una curva que son muy útiles para entender como se comporta. La curvatura expresa cuanto cambia la dirección de la trayectoria en cada punto, mientras que la torsión mide cuanto cambia el plano de curvatura a lo largo de la curva. Matemáticamente se expresan de la siguiente manera:

Representación de la curvatura y la torsión

Torsión:

[math]\tau(z) = \frac{[\vec{v}(z),\, \vec{a}(z),\, \vec{a}'(z)]}{|\vec{v}(z) \times \vec{a}(z)|^{2}}[/math]

Curvatura:

[math]\kappa(z) = \frac{|\vec{v}(z) \times \vec{a}(z)|}{|\vec{v}(z)|^{3}}[/math]

Consideramos la curva de Ekman parametrizada por la profundidad:

[math]\vec r(z)=\big(u(z),\,v(z),\,z\big),[/math]

donde

[math]u(z)=\operatorname{sgn}(f)\,V_0 e^{z/d_E}\cos\!\left(\frac{z}{d_E}+\theta\right),v(z)=V_0 e^{z/d_E}\sin\!\left(\frac{z}{d_E}+\theta\right),[/math]

[math]d_E=\sqrt{\frac{2\nu_e}{|f|}},A(z)=V_0 e^{z/d_E}.[/math]

Derivadas de la curva:

Definimos [math]\varphi(z)=\frac{z}{d_E}+\theta.[/math]

Entonces:

[math]\vec{r}'(z) = \left( \frac{\text{sgn}(f) \, A}{d_E} (\cos \varphi - \sin \varphi), \quad \frac{A}{d_E} (\sin \varphi + \cos \varphi),\quad 1 \right),[/math]

[math]\vec{r}''(z) = \left( \frac{2 \, \text{sgn}(f) \, A}{d_E^2} \sin \varphi, \quad \frac{2A}{d_E^2} \cos \varphi,\quad 0 \right),[/math]

[math]\vec{r}'''(z) = \left( \frac{2 \, \text{sgn}(f) \, A}{d_E^3} (\sin \varphi + \cos \varphi), \quad \frac{2A}{d_E^3} (\cos \varphi - \sin \varphi), \quad 0 \right)[/math]

Conclusión:

• El signo de la torsión depende del hemisferio (\(\operatorname{sgn}(f)\)).
• Hacia gran profundidad \(A(z)\to 0\), por lo que:

[math]\kappa(z)\to 0,\qquad\tau(z)\to \operatorname{sgn}(f)\frac{1}{d_E}.[/math]

• Cerca de la superficie las funciones crecen porque la velocidad inducida por el viento es mayor.

7 Triedro de Frenet

En el enunciado (10) se parametrizó la espiral de Ekman en 3D mediante

[math] \gamma(z) = \bigl(u(z), v(z), z\bigr), \qquad z \le 0, [/math]

y en el enunciado (11) se calcularon la curvatura [math]\kappa(z)[/math] y la torsión [math]\tau(z)[/math] de esta curva a partir de las expresiones generales

[math] \kappa(z) = \dfrac{\lVert \gamma'(z) \times \gamma''(z) \rVert} {\lVert \gamma'(z) \rVert^3}, \qquad \tau(z) = \dfrac{\bigl[\gamma'(z),\gamma''(z),\gamma'''(z)\bigr]} {\lVert \gamma'(z) \times \gamma''(z) \rVert^2}. [/math]

A partir de estos invariantes definimos el triedro de Frenet [math]\{\mathbf T(z),\mathbf N(z),\mathbf B(z)\}[/math] asociado a la espiral:

• El vector tangente unitario

[math] \mathbf T(z) = \dfrac{\gamma'(z)}{\lVert \gamma'(z)\rVert}, [/math]

• el vector normal principal

[math] \mathbf N(z) = \dfrac{\mathbf T'(z)}{\lVert \mathbf T'(z)\rVert}, [/math]

• y el vector binormal

[math] \mathbf B(z) = \mathbf T(z) \times \mathbf N(z). [/math]

En la práctica, las derivadas se calculan de forma numérica a partir de las coordenadas discretas de la espiral. A continuación se detalla una implementación en MATLAB que permite obtener [math]\mathbf T[/math], [math]\mathbf N[/math], [math]\mathbf B[/math], así como la curvatura [math]\kappa[/math] y la torsión [math]\tau[/math] de la curva.

7.1 12.1. Función frenet (cálculo numérico del triedro)

La función siguiente recibe las coordenadas [math](x_i,y_i,z_i)[/math] de una curva en forma de vectores y devuelve los vectores del triedro de Frenet en cada punto:

function [T,N,B,k,t] = frenet(x,y,z)
% FRENET  Triedro de Frenet de una curva dada por (x,y,z)

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

% Convertir a vectores columna
x = x(:);
y = y(:);
z = z(:);

% Primera derivada (velocidad de la curva)
dx = gradient(x);
dy = gradient(y);
dz = gradient(z);
dr = [dx dy dz];

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

% Vector tangente unitario
T = dr ./ mag(dr,3);

% Derivada del tangente
dTx = gradient(T(:,1));
dTy = gradient(T(:,2));
dTz = gradient(T(:,3));
dT  = [dTx dTy dTz];

% Vector normal unitario
N = dT ./ mag(dT,3);

% Binormal
B = cross(T,N);

% Curvatura (fórmula general)
k = mag(cross(dr,ddr),1) ./ (mag(dr,1).^3);

% Torsión (aproximada)
t = dot(-B,N,2);

end

function N = mag(T,n)
% MAG  Magnitud de un vector (Nx3)
N = sum(abs(T).^2,2).^(1/2);
d = find(N==0);
N(d) = eps*ones(size(d));
N = N(:,ones(n,1));
end

7.2 12.2. Representación del triedro en 5–6 puntos de la espiral

Para representar gráficamente el triedro de Frenet en la espiral de Ekman, utilizamos los mismos parámetros físicos fijados en los enunciados (1) y (2):

• velocidad superficial inducida [math]V_0 = 0{,}15\ \text{m/s}[/math],
• viscosidad turbulenta [math]\nu_e = 0{,}05\ \text{m}^2/\text{s}[/math],
• parámetro de Coriolis local [math]f[/math], de donde se obtiene la profundidad de Ekman [math]d_E = \sqrt{2\nu_e/|f|}[/math],
• fase inicial [math]\vartheta[/math] (ángulo de desviación superficial), calculada previamente.

Para mantener la coherencia con el enunciado, trabajamos con una variable de profundidad negativa [math]z \in [0,-3d_E][/math], donde [math]z=0[/math] corresponde a la superficie y [math]z\lt0[/math] a niveles más profundos.

La espiral de Ekman se describe entonces mediante

[math] \gamma(z) = \bigl(u(z), v(z), z\bigr), [/math]

con

[math] u(z) = \operatorname{sgn}(f)\,V_0 e^{z/d_E} \cos\!\left(\frac{z}{d_E}+\vartheta\right),\quad v(z) = V_0 e^{z/d_E} \sin\!\left(\frac{z}{d_E}+\vartheta\right),\qquad z \le 0. [/math]

El siguiente código de MATLAB construye la espiral, calcula el triedro de Frenet y lo representa en 5–6 puntos. En la gráfica el eje vertical muestra explícitamente profundidades negativas, desde [math]z=0[/math] (superficie) hasta [math]z=-3d_E[/math]:

%% Parámetros físicos (de los enunciados (1) y (2))
Omega  = 7.2921e-5;                         % Velocidad angular de la Tierra [rad/s]
phi    = deg2rad(30 + 10/60 + 24.2/3600);   % Latitud de la zona de estudio
f      = 2*Omega*sin(phi);                  % Parámetro de Coriolis [1/s]

nu_e   = 0.05;                              % Viscosidad turbulenta [m^2/s]
V0     = 0.15;                              % Velocidad superficial [m/s]
dE     = sqrt(2*nu_e/abs(f));               % Profundidad de Ekman
theta0 = 3*pi/4;                            % Fase inicial ϑ

%% Curva de la espiral en coordenadas (x,y,z)
Npts = 200;                                 % puntos para describir bien la curva
z    = linspace(0, -3*dE, Npts);            % profundidad negativa (0 = superficie)

u =  V0 .* exp(z/dE) .* cos(z/dE + theta0);
v =  V0 .* exp(z/dE) .* sin(z/dE + theta0);

% Escala horizontal para mejorar la visualización
escala = 80;
x = escala * u;
y = escala * v;

%% Triedro de Frenet a lo largo de la curva
[T,N,B,kappa,tau] = frenet(x,y,z);

%% Selección de 5–6 puntos
idx = round(linspace(10, Npts-10, 6));

%% Representación 3D
figure;
hold on;

% Curva de la espiral
plot3(x, y, z, 'k', 'LineWidth', 1.5);

% Longitud gráfica de los vectores del triedro
L = dE/5;

for j = 1:length(idx)
    i  = idx(j);
    xb = x(i); yb = y(i); zb = z(i);

    TT = L * T(i,:);
    NN = L * N(i,:);
    BB = L * B(i,:);

    quiver3(xb, yb, zb, TT(1), TT(2), TT(3), 'r', 'LineWidth', 1.5);
    quiver3(xb, yb, zb, NN(1), NN(2), NN(3), 'g', 'LineWidth', 1.5);
    quiver3(xb, yb, zb, BB(1), BB(2), BB(3), 'b', 'LineWidth', 1.5);
end

xlabel('Componente este u(z) (escalada)');
ylabel('Componente norte v(z) (escalada)');
zlabel('Profundidad z [m]');
title('Triedro de Frenet en varios puntos de la espiral de Ekman');

grid on;
view(35,25);                  % Vista tridimensional
% Aquí NO se invierte el eje Z: se muestran las profundidades negativas tal cual

hold off;

En la figura resultante se observan 5–6 sistemas de vectores ortogonales [math]\{\mathbf T,\mathbf N,\mathbf B\}[/math] distribuidos a lo largo de la espiral, con el eje vertical indicando explícitamente los valores negativos de profundidad.

7.3 12.3. Animación del triedro de Frenet

Podemos completar el estudio construyendo una animación en la que el triedro recorra la espiral de Ekman. La idea es la misma que en el apartado anterior, pero actualizando en cada fotograma los vectores [math]\mathbf T[/math], [math]\mathbf N[/math] y [math]\mathbf B[/math] en un único punto que se desplaza a lo largo de la curva.

El siguiente script de MATLAB genera la animación.

%% Animación del triedro de Frenet en la espiral de Ekman

% Parámetros físicos (los mismos del trabajo)
Omega  = 7.2921e-5;                         % [rad/s]
phi    = deg2rad(30 + 10/60 + 24.2/3600);   % Latitud Canarias
f      = 2*Omega*sin(phi);                  % Parámetro de Coriolis [1/s]

nu_e   = 0.05;                              % Viscosidad turbulenta [m^2/s]
V0     = 0.15;                              % Velocidad superficial [m/s]
dE     = sqrt(2*nu_e/abs(f));               % Profundidad de Ekman
theta0 = 3*pi/4;                            % Fase inicial ϑ

%% Curva de la espiral (x,y,z) con z negativa
Npts = 200;                                 % puntos de la espiral
z    = linspace(0, -3*dE, Npts);            % profundidad negativa

u =  V0 .* exp(z/dE) .* cos(z/dE + theta0);
v =  V0 .* exp(z/dE) .* sin(z/dE + theta0);

% Escala horizontal para mejorar la visualización
escala = 80;
x = escala * u;
y = escala * v;

%% Triedro de Frenet a lo largo de toda la curva
[T,N,B,kappa,tau] = frenet(x,y,z);

%% Preparar la figura para la animación
figure;
hold on;

% Curva de la espiral
plot3(x, y, z, 'k', 'LineWidth', 1.5);

% Vector del viento (opcional, solo decorativo)
quiver3(escala*0.15, 0, 0, -escala*0.3, 0, 0, ...
        'k', 'LineWidth', 2, 'MaxHeadSize', 0.7);

% Inicializar los manejadores de T, N, B
hT = quiver3(0,0,0, 0,0,0, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 0.5);
hN = quiver3(0,0,0, 0,0,0, 'g', 'LineWidth', 1.5, 'MaxHeadSize', 0.5);
hB = quiver3(0,0,0, 0,0,0, 'b', 'LineWidth', 1.5, 'MaxHeadSize', 0.5);

xlabel('Componente este u(z) (escalada)');
ylabel('Componente norte v(z) (escalada)');
zlabel('Profundidad z [m]');
title('Animación del triedro de Frenet en la espiral de Ekman');

grid on;
axis vis3d;
view(35,25);                  % z se muestra negativo, coherente con el modelo

hold off;

%% Parámetros de la animación
L = dE/5;                                   % longitud de los vectores del triedro
nombre_gif = 'Ekman_triedro_Frenet.gif';    % nombre del GIF (opcional)

% Si se quiere empezar un GIF nuevo, se borra el anterior
if exist(nombre_gif,'file')
    delete(nombre_gif);
end

%% Bucle de animación
for i = 1:Npts
    % Punto actual en la espiral
    xb = x(i); yb = y(i); zb = z(i);

    % Vectores del triedro escalados
    TT = L * T(i,:);
    NN = L * N(i,:);
    BB = L * B(i,:);

    % Actualizar T, N, B en el punto actual
    set(hT, 'XData', xb, 'YData', yb, 'ZData', zb, ...
            'UData', TT(1), 'VData', TT(2), 'WData', TT(3));
    set(hN, 'XData', xb, 'YData', yb, 'ZData', zb, ...
            'UData', NN(1), 'VData', NN(2), 'WData', NN(3));
    set(hB, 'XData', xb, 'YData', yb, 'ZData', zb, ...
            'UData', BB(1), 'VData', BB(2), 'WData', BB(3));

    drawnow;
    pause(0.02);                            % velocidad de la animación

    % Guardar cada frame en un GIF (opcional)
    frame = getframe(gcf);
    im    = frame2im(frame);
    [A,map] = rgb2ind(im,256);
    if i == 1
        imwrite(A,map,nombre_gif,'gif','LoopCount',Inf,'DelayTime',0.02);
    else
        imwrite(A,map,nombre_gif,'gif','WriteMode','append','DelayTime',0.02);
    end
end

En la animación se observa cómo el triedro de Frenet se desplaza a lo largo de la espiral de Ekman.