Torres de Enfriamiento Hiperbólicas (Grupo 15, Retiro)

Trabajo realizado por estudiantes
Título	Torres de Enfriamiento Hiperbolicas Grupo 15
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Miguel Cervigón, Ernesto Irazabal, Alejandro Barranquero, Guillermo Monreal
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción a las torres de enfriamiento hiperbólicas

Central Nuclear de Temelín, en la República Checa

Las torres de enfriamiento constituyen un elemento distintivo en el ámbito energético global, teniendo como principal objetivo la disipación de calor producido por centrales termoeléctricas y nucleares. La forma hiperbólica hace que la refrigeración se produzca mediante un tiro natural, aspirando el aire frio y expulsando el calor, sin necesidad de ventiladores, su geometría ofrece ademas una gran resistencia al viento.
El diseño hiperbólico para torres de enfriamiento se introdujo en la década de 1950, cuando la industria energética buscaba soluciones más eficientes para la refrigeración en centrales térmicas y nucleares. Antes de esto, se utilizaban estructuras cilíndricas, que requerían grandes cantidades de material y ofrecían menor estabilidad frente al viento. Desde entonces, este diseño se convirtió en el estándar mundial, consolidándose como símbolo de la ingeniería moderna en el sector energético.

2 Modelo geométrico del hiperboloide

Una torre de enfriamiento hiperbólica, se caracteriza por su altura máxima [math]H[/math], su radio minimo [math]Rmin[/math], y su radio máximo [math]Rmax[/math] alcanzado en la base de la torre. La superficie de la torre sigue la forma de un hiperboloide hiperbólico con centro a la altura [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math] , y se modela con la ecuación de un hiperboloide de revolución de una hoja:

[math]\dfrac{x^2 + y^2}{a^2} - \dfrac{(z - z_0)^2}{c^2} = 1[/math]

Donde

[math]Rmax=55m,\qquad Rmin=30m,\qquad H=150m[/math]

y

[math] a=Rmin,\qquad c= Parámetro \ de \ curvatura,\qquad z_0= Centro \ del \ hiperboloide[/math]

2.1 Ecuación del hiperboloide y parámetros

Para poder definir el hiperboloide es imprescindible hallar los parámetros a, c y z₀, para ello nos apoyaremos en el cambio de coordenadas cilíndricas:

[math]\begin{cases} x = \rho \cos\theta \\ y = \rho \sin\theta \\ z = z \end{cases}[/math]

La ecuación de un hiperboloide en coordenadas cartesianas es:

[math]\frac{x^2 + y^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1[/math]

Sabiendo que en coordenadas cilíndricas x² + y² = ρ², la ecuación se reescribira de la siguiente manera:

[math]\frac{\rho^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1[/math]

• Se conocen los siguientes datos:

Centro en z₀ = 100 m, altura total H = 150 m, radio mínimo ρ_min = 30 m, radio máximo ρ_max = 55 m.

• Empezaremos calculando el paremetro a:

Sabemos que en el plano z = z₀, el segundo término de la ecuación:

[math]\frac{\rho^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1[/math]

se anula porque (z - z₀) = 0, lo que implica que estamos en la sección más estrecha del hiperboloide. Por tanto, la ecuación se simplifica a:

[math]\frac{\rho^2}{a^2} = 1[/math]

Sustituimos ρ = 30 m:

[math]\frac{30^2}{a^2} = 1 \implies a = 30 \text{ m}[/math]

Con el parámetro a ya calculado, el cálculo del parámetro c se facilita enormemente:
Con a = 30 m y ρ = 55 m:

[math]\frac{55^2}{30^2} - \frac{(z - z_0)^2}{c^2} = 1[/math]

Despejamos:

[math]\frac{(z - z_0)^2}{c^2} = \frac{55^2}{30^2} - 1[/math]

Calculamos con z - z₀ = 100 m:

[math]\frac{100^2}{c^2} = \frac{3025}{900} - 1 = 2.361[/math] [math]\frac{10000}{c^2} = 2.361 \implies c^2 = \frac{10000}{2.361} \approx 4236.7 \implies c \approx 65.079 \text{ m}[/math]

• Los parametros obtenidos tras los calculos son los siguientes:

[math]a = 30 \text{ m}, \quad c \approx 65.079 \text{ m}, \quad z_0 = 100 \text{ m}[/math]

2.2 Representación grafica de la superficie parametrizada

Primero se establecen los parámetros con sus distintos rangos (se define [math] v=\dfrac{z-z_0}{c} [/math]) , después se definen las ecuaciones paramétricas (despejando el radio en la expresión en coordenadas cilíndricas) y por último se grafica la función.

% 1
a = 30;       % radio mínimo
c = 65.1;     % parámetro vertical
z0 = 100;     % centro del hiperboloide
H = 150;      % altura total

% 2
u = linspace(0, 2*pi, 100);   % ángulo
v_min = (0 - z0)/c;           % para z = 0
v_max = (H - z0)/c;           % para z = H
v = linspace(v_min, v_max, 100);

[U,V] = meshgrid(u,v);

% 3
X = a * sqrt(1 + V.^2) .* cos(U);
Y = a * sqrt(1 + V.^2) .* sin(U);
Z = z0 + c * V;

% 4
figure;
surf(X,Y,Z);
axis equal;
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Torre de Enfriamiento Hiperbólica');
shading faceted;

3 Torres de enfriamiento hiperbólicas famosas en el mundo

A lo largo del mundo existen una serie de famosas centrales nucleares con sus propias torres de enfriamiento hiperbólicas, las más destacadas son las siguientes:

Torre de Cantón

• Central nuclear de Cofrentes (España)

Se sitúa en la confluencia de los ríos Juncar y Cabriel. Sus dos torres de enfriamiento de tiro natural son vitales porque la central opera en un ciclo cerrado para minimizar el impacto térmico en el río y permitir la vida en él su geometría está dictada por la necesidad de enfriamiento frente a las altas temperaturas de la zona.

• Central nuclear de Kashiwazaki-Kariwa:

La central nuclear más grande del mundo en capacidad instalada (unos 8000MW). La estructura hiperbólica de sus torres de enfriamiento han permitido que resista terremotos de magnitud 6.8 en la escala de Richter.

• Central térmica de Niederaussem (Alemania):

La principal característica de dicha central es que alberga la torre de enfriamiento más alta del mundo (200 metros). La forma hiperbólica es la única viable para una estructura de tal magnitud, ya que una cilíndrica no aguantaría los momentos de vuelco. Esta torre también tiene funcionalidad como “chimenea”. Los gases de combustión desulfurados se inyectan en la atmósfera, de forma que la estructura se usa para dos cosas a la vez, enfriar y expulsar gases tratados.

• Las torres de Didcot (Reino Unido):

La demolición de esta estructura mostró con claridad la superioridad de esta forma. Al detonar los explosivos en la base, las torres se derrumbaron sobre su propia planta, mostrando que no depende de los materiales sino de su forma.
[math]\textbf{Otras estructuras hiperbólicas}[/math]

1. Torre Shabolovka (Moscú)
Pese a la escasez de acero en la Rusia post revolucionaria, la gran eficiencia de las torres hiperbólicas a nivel de materiales permitió su construcción. La torre de comunicaciones de Shújov (coloquialmente conocida como torre Shabolovka) fue la primera torre hiperbólica en construirse.
2. Torre de Cantón (Guangzhou, China):
Con 600 metros de altura. Es un hiperboloide generado por columnas de acero rectas pero inclinadas, creando una cintura estrecha. Esta forma se utiliza para resistencia sísmica y de tifones.

4 Superioridad de la forma hiperbólica ante la cilíndrica

Para abordar esto con la profundidad requerida, se contrastará desde los siguientes dos aspectos físicos: estabilidad estructural, resistencia al viento y termodinámica. La carga crítica para una torre de enfriamiento no es su propio peso sino la presión dinámica del viento, cuya velocidad viene definida por la ley de potencia:

[math]V(z) = V_0 \left( \dfrac{z}{z_{ref}} \right)^{\alpha}[/math]

Y la presión dinámica del viento:

[math]P(z) = \dfrac{1}{2} \varrho_{aire} V(z)^2[/math]

[math]\textbf{Problemas del cilindro}[/math]

La "doble curvatura" del hiperboloidele otorga mayor resistencia a esfuerzos. Nótese que el punto de mayor resistencia está en el estrechamiento, donde la curvatura es más pronunciada.

El cilindro tiene varias desventajas frente a cargas laterales. Al ser una superficie con curvatura nula, no aporta rigidez geométrica, por lo que depende del grosor del material para resistir esfuerzos. Cuando actúa el viento, la sección circular puede deformarse en una elipse, lo que obliga a usar paredes muy gruesas o refuerzos adicionales. Además, su diámetro constante limita la estabilidad frente al vuelco, ya que la base no ofrece una superficie amplia.

El hiperboloide, en cambio, presenta curvatura negativa, lo que le da ventajas importantes. Su doble curvatura aporta rigidez, reduciendo esfuerzos de flexión sin aumentar el espesor. También desvía las fuerzas hacia las generatrices y la base, mejorando la resistencia al pandeo. La diferencia entre el radio mayor y el menor crea una base más amplia, aumentando la estabilidad y reduciendo el consumo de material.

Desde el punto de vista aerodinámico, el hiperboloide reduce la presión del viento; los cálculos muestran que soporta hasta un 40 % menos carga lateral que un cilindro de igual altura y radio máximo. Además, su forma favorece corrientes ascendentes de aire, útiles para enfriamiento mediante el efecto Venturi y el principio de Bernoulli.

5 Análisis en profundidad de la presión del viento

5.1 Representación del campo escalar P(z)

El viento sopla en dirección del vector:

[math]-\dfrac{1}{\sqrt{2}} (\vec{i} + \vec{j})[/math]

Mapa de presión sobre la torre

   Este cuadro de texto está debajo de la imagen alineada a la derecha. Aquí puedes añadir una descripción o explicación adicional.

%% Mapa de presión del viento sobre la torre hiperboloidal
clear; clc; close all;

% Parámetros del viento
V0 = 18;          % velocidad a zref [m/s]
zref = 10;        % altura de referencia [m]
alpha = 1/7;      % exponente terreno abierto
rho_air = 1.2;    % densidad aire [kg/m^3]

% Geometría de la torre
a = 30;           % radio mínimo [m]
c = 65.1;         % parámetro vertical [m]
z0 = 100;         % centro del hiperboloide [m]
H = 150;          % altura total [m]

% Mallado paramétrico (media torre)
u = linspace(0, pi, 100);    % mitad expuesta
v_min = (0 - z0)/c;
v_max = (H - z0)/c;
v = linspace(v_min, v_max, 200);
[U,V] = meshgrid(u,v);

% Coordenadas
X = a * sqrt(1 + V.^2) .* cos(U);
Y = a * sqrt(1 + V.^2) .* sin(U);
Z = z0 + c * V;

% Calcular velocidad y presión en cada altura
Vwind = V0 * (Z / zref).^alpha;           % ley de potencia
Pwind = 0.5 * rho_air .* Vwind.^2;        % presión dinámica

% Encontrar punto de máxima presión
[maxP, idx] = max(Pwind(:));
[maxRow, maxCol] = ind2sub(size(Pwind), idx);
Xmax = X(maxRow, maxCol);
Ymax = Y(maxRow, maxCol);
Zmax = Z(maxRow, maxCol);

% Graficar superficie con mapa de presión
figure('Color','w');
surf(X,Y,Z,Pwind,'EdgeColor','none');
axis equal;
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Mapa de presión del viento sobre la torre hiperboloidal');
colormap(jet);
colorbar;

% Marcar punto de máxima presión
hold on;
plot3(Xmax,Ymax,Zmax,'ko','MarkerFaceColor','b','MarkerSize',10);
text(Xmax,Ymax,Zmax+5,sprintf('Max P = %.1f Pa',maxP),'Color','k','FontSize',10);

5.2 Representación del campo vectorial de fuerza sobre la superficie expuesta

La presión de viento genera cierta fuerza sobre la superficie expuesta su representación gráfica es la siguiente:

Mapa de presión sobre la torre

   Este cuadro de texto está debajo de la imagen alineada a la derecha. Aquí puedes añadir una descripción o explicación adicional.

% --- Parámetros de la torre (ya calculados)
H = 150;
Rmax = 55; Rmin = 30;
h = 2/3*H; a = Rmin; z0 = h;
c = h / sqrt((Rmax^2 / a^2) - 1);   % parámetro vertical del hiperboloide

% --- Definir V(z) y P(z)
V   = @(z) 18*(z./10).^(1/7);   % z > 0
rho = 1.204;
Pfun = @(z) 0.5*rho.*(V(z)).^2;

% --- Malla paramétrica (theta, z) COMPLETA (0..2pi)
nt = 120; nz = 200;
theta = linspace(0,2*pi,nt);
z     = linspace(0, H, nz);
[TH, Z] = meshgrid(theta, z);    % size nz x nt

% --- radio r(z)
r = a * sqrt(1 + (Z - z0).^2 ./ c^2);   % nz x nt

% --- coordenadas cartesianas de la superficie (transpuestas para nt x nz)
X  = (r .* cos(TH))';   % nt x nz
Y  = (r .* sin(TH))';
Zt = Z';                % nt x nz

% --- valores P(z) (igual para todas las theta a una z dada)
Pz = Pfun(Zt);          % nt x nz

% --- Gradiente de G: G = x^2/a^2 + y^2/a^2 - (z-z0)^2/c^2 - 1
%     ∇G = [2x/a^2, 2y/a^2, -2(z-z0)/c^2]
Gx = 2*X./a^2;
Gy = 2*Y./a^2;
Gz = -2*(Zt - z0)./c^2;

% --- Normal unitario (hacia exterior)
Nmag = sqrt(Gx.^2 + Gy.^2 + Gz.^2);
Nx = Gx ./ Nmag;
Ny = Gy ./ Nmag;
Nz = Gz ./ Nmag;

% --- Fuerza: F = -P(z) * n
Fx = -Pz .* Nx;
Fy = -Pz .* Ny;
Fz = -Pz .* Nz;

% === Dirección del viento (horizontal) ===
% azimut en grados medido desde +X hacia +Y (sentido antihorario).
az_deg = 45;  % ejemplo: viento del Suroeste hacia el Noreste (apunta a NE)
az = deg2rad(az_deg);
w_hat = [cos(az); sin(az); 0];  % vector unitario

% Producto escalar n·w_hat (usar la normal exterior)
dotNW = Nx*w_hat(1) + Ny*w_hat(2) + Nz*w_hat(3);  % nt x nz

% --- Máscara de la mitad expuesta a barlovento:
% puntos donde la presión empuja hacia dentro: F = -P n (⇒ n·w_hat < 0)
mask = (dotNW < 0);

% Aplicar máscara (ocultamos la otra mitad con NaN)
X(~mask)   = NaN;  Y(~mask)   = NaN;  Zt(~mask)  = NaN;
Fx(~mask)  = NaN;  Fy(~mask)  = NaN;  Fz(~mask)  = NaN;

% Representación
figure('Color','w')
Fmag = sqrt(Fx.^2 + Fy.^2 + Fz.^2);
surf(X, Y, Zt, Fmag, 'EdgeColor','none', 'FaceAlpha',0.9)
axis equal
xlabel('x (m)'), ylabel('y (m)'), zlabel('z (m)')
title(sprintf('Campo F = -P(z)·n (mitad expuesta, azimut viento = %d°)', az_deg))
colormap(parula)
cb = colorbar;                      % evitar conflicto con 'c' geométrico
cb.Label.String = '

5.3 Cálculo de la fuerza total y la fuerza por unidad de superficie

Queremos la fuerza total debida al viento sobre la mitad de la superficie lateral expuesta de la torre.

• El viento genera una presión que solo depende de la altura:

[math] V(z) = V_0 \left(\frac{z}{z_{\text{ref}}}\right)^\alpha, \quad V_0 = 18\,\text{m/s},\ z_{\text{ref}}=10\,\text{m},\ \alpha=\tfrac17. [/math]

• La presión dinámica del viento es:

[math] P(z) = \tfrac12\,\rho_{\text{aire}}\,V(z)^2, [/math]

Donde se toma como densidad del aire estandar, [math]\rho_{\text{aire}} \approx 1{,}204\ \text{kg/m}^3[/math]

Como se pide “integrar la presión sobre la mitad de la superficie expuesta”, matemáticamente es una integral de superficie de un campo escalar:

[math] F_{\text{total}} = \iint_{S_{\text{media}}} P(z)\,dS. [/math]

A partir de la geometría y parametrización de la torre, un hiperboloide de revolución de una hoja, con ecuación:

[math] \frac{x^2 + y^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1, [/math]

y con:

• Altura total: [math]H = 150\,\text{m}[/math].
• Radio máximo en la base: [math]R_{\max} = 55\,\text{m}[/math].
• Radio mínimo (estrangulamiento): [math]R_{\min} = 30\,\text{m}[/math], a la altura
• Elegimos [math]z_0 = h = 100\,\text{m}[/math] y [math]a = R_{\min} = 30\,\text{m}[/math].

Obtendremos el radio en función de la altura (coordenadas cilíndricas, [math]\rho[/math]):

[math] \rho(z) = \sqrt{x^2 + y^2} = a\,\sqrt{1 + \frac{(z - z_0)^2}{c^2}}. [/math]

Imponiendo que en la base [math]z=0[/math] el radio sea [math]R_{\max} = 55\,\text{m}[/math]:

[math] R_{\max}^2 = a^2\left(1 + \frac{z_0^2}{c^2}\right) \;\Longrightarrow\; c^2 = \frac{z_0^2}{\dfrac{R_{\max}^2}{a^2}-1}. [/math]

Con los valores numéricos: [math] H = 150,\quad z_0 = 100,\quad a = 30,\quad R_{\max} = 55, [/math] se obtiene aproximadamente [math]c \approx 65{,}08\,\text{m}[/math].

Ahora parametrizamos la superficie de revolución usando [math](\theta,z)[/math] como parámetros:

[math] \mathbf r(\theta,z) = \bigl(\rho(z)\cos\theta,\ \rho(z)\sin\theta,\ z\bigr). [/math]

• [math]z[/math] recorre la altura: [math]0 \le z \le H[/math].
• [math]\theta[/math] recorre el ángulo alrededor del eje. Para toda la torre [math]0 \le \theta < 2\pi[/math].

Para la mitad expuesta al viento tomamos un intervalo de longitud de media circunferencia por ejemplo: [math] -\frac{\pi}{2} \le \theta \le \frac{\pi}{2}, [/math] que corresponde a la “media circunferencia” que mira hacia el viento.

El elemento de área es:

[math] dS = \bigl|\mathbf r_\theta\times \mathbf r_z\bigr|\,d\theta\,dz. [/math]

Calculamos las derivadas parciales:

• Derivada respecto a [math]\theta[/math]:

[math] \mathbf r_\theta = \bigl(-\rho(z)\sin\theta,\ \rho(z)\cos\theta,\ 0\bigr). [/math]

• Derivada respecto a [math]z[/math]:

[math] \mathbf r_z = \bigl(\rho'(z)\cos\theta,\ \rho'(z)\sin\theta,\ 1\bigr), [/math]

El producto vectorial (simplificado) es:

[math] \mathbf r_\theta \times \mathbf r_z = \bigl(\rho(z)\cos\theta,\ \rho(z)\sin\theta,\ -\rho(z)\rho'(z)\bigr). [/math]

Su módulo:

[math] \bigl|\mathbf r_\theta \times \mathbf r_z\bigr| = \sqrt{\rho^2\cos^2\theta + \rho^2\sin^2\theta + \rho^2(\rho')^2} = \rho(z)\sqrt{1 + \rho'(z)^2}. [/math]

Por lo tanto:

[math] dS = \rho(z)\,\sqrt{1 + \rho'(z)^2}\,d\theta\,dz. [/math]

Esta es la fórmula general del área de una superficie de revolución.

La derivada de [math]\rho(z)[/math] a partir de la fórmula del hiperboloide es:

[math] \rho(z) = a\sqrt{1 + \frac{(z-z_0)^2}{c^2}} \;\Longrightarrow\; \rho'(z) = \frac{a(z-z_0)}{c^2\sqrt{1 + \dfrac{(z-z_0)^2}{c^2}}}. [/math]

Finalmente laIntegral de la fuerza total se define como:

[math] F_{\text{total}} = \iint_{S_{\text{media}}} P(z)\,dS = \int_{z=0}^{H} \int_{\theta=-\pi/2}^{\pi/2} P(z)\,\rho(z)\sqrt{1+\rho'(z)^2}\,d\theta\,dz. [/math]

Como [math]P(z)[/math] no depende de [math]\theta[/math], la integral en [math]\theta[/math] se puede sacar:

[math] \int_{-\pi/2}^{\pi/2} d\theta = \pi. [/math]

Por tanto:

[math] F_{\text{total}} = \pi \int_0^{H} P(z)\,\rho(z)\,\sqrt{1+\rho'(z)^2}\,dz, [/math]

[math] \rho(z) = a\sqrt{1 + \frac{(z-z_0)^2}{c^2}},\quad \rho'(z) = \frac{a(z-z_0)}{c^2\sqrt{1 + \dfrac{(z-z_0)^2}{c^2}}}. [/math]

Área de la mitad de la torre y fuerza por unidad de superficie La superficie de la mitad de la torre expuesta al viento es la integral de [math]1[/math] sobre esa superficie:

[math] A_{\text{media}} = \iint_{S_{\text{media}}} dS = \int_0^H \int_{-\pi/2}^{\pi/2} \rho(z)\sqrt{1+\rho'(z)^2}\,d\theta\,dz = \pi\int_0^H \rho(z)\sqrt{1+\rho'(z)^2}\,dz. [/math]

La fuerza media por unidad de superficie (presión media) es:

[math] P_{\text{media}} = \frac{F_{\text{total}}}{A_{\text{media}}} = \frac{\displaystyle \int_0^H P(z)\,\rho(z)\sqrt{1+\rho'(z)^2}\,dz} {\displaystyle \int_0^H \rho(z)\sqrt{1+\rho'(z)^2}\,dz}. [/math]

Obsérvese que el factor [math]\pi[/math] se simplifica.

Valor numérico aproximado si se evalúan las integrales numéricamente:

• [math]H = 150\,\text{m}[/math],
• [math]a = 30\,\text{m}[/math], [math]z_0 = 100\,\text{m}[/math],
• [math]c \approx 65{,}08\,\text{m}[/math],
• [math]\rho_{\text{aire}} = 1{,}204\,\text{kg/m}^3[/math],
• [math]V(z) = 18\,(z/10)^{1/7}[/math],

Se obtiene aproximadamente:

• Fuerza total sobre la mitad de la torre:

[math] F_{\text{total}} \approx 5{,}7\times 10^6\ \text{N}. [/math]

• Área de la mitad de la torre:

[math] A_{\text{media}} \approx 1{,}8\times 10^4\ \text{m}^2. [/math]

• Fuerza por unidad de superficie (presión media):

[math] P_{\text{media}} \approx 3{,}1\times 10^2\ \text{Pa}. [/math]

6 Analisis del flujo de aire en el interior de la torre

6.1 Representación del campo vectorial de velocidad [math]\vec{v}[/math]

Mapa de presión sobre la torre

   Este cuadro de texto está debajo de la imagen alineada a la derecha. Aquí puedes añadir una descripción o explicación adicional.

%% Sección vertical de la torre hiperboloidal y campo de velocidad v_z(x,z)
clear; close all; clc;

%% Parámetros geométricos de la torre
H    = 150;   % altura total (m)
Rmax = 55;    % radio en la base (m)
Rmin = 30;    % radio mínimo (estrangulamiento) (m)
z0   = 100;   % altura del estrangulamiento (m)

% Parámetro c del hiperboloide calculado a partir de la condición de la base:
c = sqrt( (z0^2) / ( (Rmax^2 / Rmin^2) - 1 ) );

% Función radio de la torre a la altura z (sección vertical)
rho_torre = @(z) Rmin .* sqrt( 1 + ((z - z0).^2) ./ c.^2 );  % a = Rmin

%% Parámetros del campo de velocidad
vmax = 4;     % m/s, velocidad máxima en el eje central en el tope

% v_z(x,z) = vmax * (1 - rho^2 / rho_torre(z)^2) * (z/H)^0.6, con rho = abs(x)
vz_fun = @(x,z) vmax .* ( 1 - (abs(x).^2) ./ (rho_torre(z).^2) ) .* (z./H).^(0.6);

%% Mallado del plano x-z (sección que contiene el eje)
Nx = 161; Nz = 161;
x  = linspace(-Rmax, Rmax, Nx);
z  = linspace(0, H, Nz);
[X,Z] = meshgrid(x, z);

% Máscara: puntos dentro de la torre (abs(x) <= rho_torre(z))
RHOz   = rho_torre(Z);
inside = abs(X) <= RHOz;

%% Campo de velocidad vertical
Vz = vz_fun(X, Z);
Vz(~inside) = NaN;  % fuera de la torre, NaN para transparencia

% Componente horizontal (cero) para quiver
Ux = zeros(size(Vz));
Ux(~inside) = NaN;

%% Figura: mapa de color + contorno + flechas con fondo blanco
figure('Color','w');
h = imagesc(x, z, Vz);
set(h, 'AlphaData', ~isnan(Vz));  % Transparencia fuera de la torre
set(gca,'YDir','normal','Color','w');  % Fondo blanco en el área fuera de datos
axis equal tight;
colormap(turbo);
cb = colorbar;
cb.Label.String = 'v_z (m/s)';

hold on;
% Contorno de la sección de la torre
plot( rho_torre(z), z, 'k-', 'LineWidth', 2);
plot(-rho_torre(z), z, 'k-', 'LineWidth', 2);

% Flechas del campo (muestreo para no saturar)
stepX = 6; stepZ = 6;
Xq = X(1:stepZ:end, 1:stepX:end);
Zq = Z(1:stepZ:end, 1:stepX:end);
Uq = Ux(1:stepZ:end, 1:stepX:end);
Vq = Vz(1:stepZ:end, 1:stepX:end);

quiver(Xq, Zq, Uq, Vq, 'k', 'LineWidth', 1, 'AutoScale', 'on', 'AutoScaleFactor', 0.8);

xlabel('x (m)');
ylabel('z (m)');
title('Sección vertical (x-z) y campo de velocidad v_z dentro de la torre');
grid on;

Animación mostrando el mapa de velocidades en distintas secciones horizontales a distintas alturas

   Este cuadro de texto está debajo de la imagen alineada a la derecha. Aquí puedes añadir una descripción o explicación adicional.

%% Animación 3D del campo de velocidad en planos horizontales con la torre
clear; close all; clc;

%% Parámetros geométricos
H    = 150; Rmax = 55; Rmin = 30; z0 = 100;
c    = sqrt( (z0^2) / ( (Rmax^2 / Rmin^2) - 1 ) );
rho_torre = @(z) Rmin .* sqrt( 1 + ((z - z0).^2) ./ c.^2 );

%% Campo de velocidad
vmax = 4;
vz_fun = @(rho,z) vmax .* (1 - (rho.^2) ./ (rho_torre(z).^2)) .* (z./H).^(0.6);

%% Mallado XY para planos
Nx = 101; Ny = 101;
x = linspace(-Rmax,Rmax,Nx);
y = linspace(-Rmax,Rmax,Ny);
[X,Y] = meshgrid(x,y);
RHO = sqrt(X.^2 + Y.^2);

%% Mallado para la torre (hiperboloide)
theta = linspace(0,2*pi,100);
z_vals_torre = linspace(0,H,100);
[Theta,Zt] = meshgrid(theta,z_vals_torre);
Rt = rho_torre(Zt);
Xt = Rt .* cos(Theta);
Yt = Rt .* sin(Theta);

clf
%% Configuración de animación
% figure('Color','w')
% Establecer el color de fondo de la figura a blanco
set(gca, 'Color', 'w');
axis equal tight;
colormap(turbo);
axis([-Rmax Rmax -Rmax Rmax 0 H]);
xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
title('Animación 3D del campo de velocidad en planos horizontales');
view(45,25);  % Vista con acimut 45° y cenit 25°
hold on

% Dibujar la torre como superficie transparente
surf(Xt,Yt,Zt,'FaceAlpha',0.1,'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);

frames = 50;
z_vals = linspace(0,H,frames);
F(frames) = struct('cdata',[],'colormap',[]);

t=0;
while t<= frames
    for k = 1:frames
        t=k;
        zc = z_vals(k);
        Rz = rho_torre(zc);
        
        % Campo de velocidad en este plano
        Vz_plane = vz_fun(RHO,zc);
        Vz_plane(RHO > Rz) = NaN;  % fuera de la torre
        
        % Dibujar plano horizontal en z = zc
        h = surf(X,Y,zc*ones(size(X)),Vz_plane,'EdgeColor','none');
        
        % Ajustes visuales
        axis([-Rmax Rmax -Rmax Rmax 0 H]);
        view(45,30);
        colorbar; caxis([0 vmax]);
        
        title(sprintf('Plano horizontal a z = %.1f m',zc));
        drawnow;
        
        % Capturar frame
        F(k) = getframe(gcf);
        
        % Eliminar el plano anterior para el siguiente frame
        delete(h);
    end

    break
end

6.2 Divergencia del campo de velocidades

Campo de velocidad:

[math] v_z(\rho, z) = v_{max} \left( 1 - \frac{\rho^2}{R(z)^2} \right) \left( \frac{z}{H} \right)^{0,6} [/math]

Divergencia en cilíndricas:

[math] \nabla \cdot \vec{v} = \frac{1}{\rho} \left( \frac{\partial}{\partial \rho}(\rho v_\rho) + \frac{\partial}{\partial \theta}(v_\theta) + \frac{\partial}{\partial z}(\rho v_z) \right) [/math]

Al solo tener componente vertical ([math]v_z[/math]), la expresión se simplifica:

[math] \nabla \cdot \vec{v} = \frac{1}{\rho} \left( \frac{\partial}{\partial z}(\rho v_z) \right) [/math]

Sustituimos la expresión de la velocidad:

[math] \nabla \cdot \vec{v} = \frac{1}{\rho} \cdot \frac{\partial}{\partial z} \left[ \rho \cdot v_{max} \left( 1 - \frac{\rho^2}{R(z)^2} \right) \left( \frac{z}{H} \right)^{0,6} \right] [/math]

Al ser derivada respecto de [math]z[/math], se puede sacar el [math]\rho[/math] de dentro de la derivada y simplificarla con el [math]\frac{1}{\rho}[/math] de fuera. Aplicando la regla del producto para derivar respecto a [math]z[/math]:

[math] \nabla \cdot \vec{v} = v_{max} \left[ \frac{2\rho^2 R'(z)}{R(z)^3} \cdot \left( \frac{z}{H} \right)^{0,6} + \left( 1 - \frac{\rho^2}{R(z)^2} \right) \cdot \frac{0,6}{z} \left( \frac{z}{H} \right)^{0,6} \right] [/math]

• Zona 1: Del suelo hasta la garganta ([math]z=0 \to z=100[/math])
  • Al subir, la torre se estrecha, el radio disminuye (Negativo).
  • Al subir, la velocidad aumenta por el tiro térmico (Positivo).

• Zona 2: De la garganta a la parte más alta ([math]z=100 \to z=150[/math])
  • Al subir, la torre se ensancha, el radio aumenta (Positivo).
  • Al subir, la velocidad aumenta por el tiro térmico (Positivo).

[math] \Rightarrow \text{Es positiva} [/math]

Físicamente, una divergencia positiva indica que el fluido se expande a medida que asciende. Esto se debe a:

1. Expansión térmica: El aire absorbe el calor, aumentando su temperatura y disminuyendo su densidad. Para conservar la masa, el gas debe aumentar su volumen.
2. Tiro vertical: Se acelera más rápido que lo que se ensancha (o estrecha), por tanto también se expande verticalmente ("estiramiento" del flujo).

6.3 Caudal volumétrico de aire en una sección a altura Z

Nos piden calcular:

[math] Q(z) = \iint_S v_z(\rho, z) dS [/math]

En coordenadas cilíndricas [math] dS = \rho d\rho d\theta [/math]:

• [math] 0 \le \theta \le 2\pi [/math]
• [math] 0 \le \rho \le R(z) [/math]

Planteamos la integral:

[math] Q(z) = \int_{0}^{2\pi} \int_{0}^{R(z)} \underbrace{v_{max} \left( \frac{z}{H} \right)^{0,6}}_{\text{constante}} \cdot \left( 1 - \frac{\rho^2}{R(z)^2} \right) \underbrace{\rho d\rho d\theta}_{\text{Jacobiano}} [/math]

Separamos los términos constantes respecto a la integral radial y angular:

[math] Q(z) = v_{max} \left( \frac{z}{H} \right)^{0,6} \cdot \int_{0}^{2\pi} d\theta \cdot \int_{0}^{R(z)} \left( 1 - \frac{\rho^2}{R(z)^2} \right) \rho d\rho [/math]

Resolvemos la integral radial:

[math] \int_{0}^{R(z)} \left( \rho - \frac{\rho^3}{R(z)^2} \right) d\rho = \left[ \frac{\rho^2}{2} - \frac{\rho^4}{4 R(z)^2} \right]_{0}^{R(z)} [/math]

Evaluando en los límites:

[math] = \left[ \frac{R(z)^2}{2} - \frac{R(z)^4}{4 R(z)^2} \right] = \frac{R(z)^2}{2} - \frac{R(z)^2}{4} = \frac{R(z)^2}{4} [/math]

Sustituimos el resultado en la expresión original del caudal:

[math] Q(z) = v_{max} \left( \frac{z}{H} \right)^{0,6} \cdot 2\pi \cdot \frac{R(z)^2}{4} [/math]

Resultado final:

[math] Q(z) = \frac{\pi}{2} v_{max} R(z)^2 \left( \frac{z}{H} \right)^{0,6} [/math]

1. ¿Es el caudal constante o varía con la altura?

   Sí varía con la altura. En la expresión del caudal aparecen [math]z[/math] y [math]R(z)[/math], ambas variables dependientes de la altura.
   

2. ¿Qué implica esto físicamente?

   Implica que la densidad del aire está disminuyendo. Ya que si el caudal aumenta por la altura, por la ley de la conservación de la masa, la densidad disminuye.
   

3. ¿Dónde es máximo y por qué?

   La expresión del caudal se maximiza cuando el radio es máximo y la altura es máxima, esto se da en la cima de la torre.

6.3.1 Representación grafica de Q(z) en función de z

% Parámetros
Vmax = 10;       % Velocidad máxima (m/s)
H = 150;         % Altura total (m)
R0 = 30;         % Radio mínimo en la base (m)
Rmax = 55;       % Radio máximo en la parte superior (m)

% Vector de alturas
z = linspace(0, H, 200);

% Función del radio en función de z (lineal)
Rz = R0 + (Rmax - R0) * (z / H);

% Cálculo del caudal Q(z)
Q = (pi/2) * Vmax .* (Rz.^2) .* ((z / H).^0.6);

% Gráfica
figure;
plot(z, Q, 'LineWidth', 2);
xlabel('Altura z (m)');
ylabel('Caudal Q(z) (m^3/s)');
title('Caudal volumétrico en función de la altura');
grid on;

% Mostrar el valor máximo y su posición
[maxQ, idx] = max(Q);
fprintf('Caudal máximo: %.2f m^3/s en z = %.2f m\n', maxQ, z(idx));

6.4 Cálculo de la potencia disipada

Calculamos la potencia térmica disipada:

[math] \dot{Q} = \rho_{aire} c_p Q(H) (T_{base} - T_{top}) [/math]

Datos:

• [math] c_p = 1005 \, J/(kg \cdot K) [/math]
• [math] \Delta T_z = 38^\circ C [/math]
• [math] \rho_{aire} = 1,225 \, kg/m^3 [/math] (Buscado)
• [math] v_{max} = 4 \, m/s [/math]

Calculamos Q(H) usando la fórmula hallada antes:

[math] Q(H) = \frac{\pi}{2} v_{max} R(H)^2 \cdot \underbrace{\left( \frac{H}{H} \right)^{0,6}}_{1} \Rightarrow Q(H) = 2\pi R(H)^2 [/math]

Necesitamos el radio en la cima ([math]R(z)[/math]):

[math] R(z) = a \sqrt{1 + \frac{(z-z_0)^2}{c^2}} \rightarrow R(150) = 30 \sqrt{1 + \frac{(150-100)^2}{4235,4}} = 37,82 \, m [/math]

Sustituimos en el caudal:

[math] Q(H) = \frac{\pi}{2} (4) (37,82)^2 \approx 8987 \, m^3/s [/math]

Calculamos la potencia:

[math] \dot{Q} = 1,225 \cdot 1005 \cdot 8987 \cdot 38 = 420.450.000 \, W = 420,45 \, MW [/math]

Análisis con Radio mínimo = 20m

Habrá que calcular una nueva [math]c[/math] ([math]c'[/math]) con [math]a=20[/math]:

[math] 55 = 20 \sqrt{1 + \frac{100^2}{c'^2}} \Rightarrow c' = 1523,8 [/math]

Nuevo radio en la cima:

[math] R'(150) = 20 \sqrt{1 + \frac{(150-100)^2}{1523,8}} = 32,5 \, m [/math]

Nuevo caudal:

[math] Q'(H) = \frac{\pi}{2} (4) (32,5)^2 = 6636 \, m^3/s [/math]

Es menor debido a que se cumple [math]Q = V \cdot A[/math]. Si se mantiene la velocidad y el área disminuye, el caudal también lo hace.

Conclusiones sobre el estrangulamiento

• ¿Cómo afecta un estrangulamiento mayor al flujo del aire?

  Como se ha explicado antes con la relación [math]Q = V \cdot A[/math], al reducir el área, se reduce el flujo.

• ¿Existe un compromiso entre resistencia estructural (menor Rmin) y eficiencia de enfriamiento?

  Sí. Estrechar el cuello tiene de positivo una mayor resistencia estructural y rigidez con una capa de poco espesor.

  Tiene de negativo, como se ha demostrado antes, que disminuye el caudal, y por tanto, disminuye su capacidad para evacuar calor.

6.5 Análisis del número de Reynolds en el flujo

El número de Reynolds caracteriza el flujo:

[math] Re = \frac{\rho v_{tipica} L}{\mu} [/math]

Datos:

• [math] \mu = 1,8 \cdot 10^{-5} \, Pa \cdot s [/math]
• [math] \rho_{aire} = 1,225 \, kg/m^3 [/math]
• Longitud característica: [math] L = 2 \cdot R_{min} = 2 \cdot 30 = 60 \, m [/math]

Cálculo de velocidad típica: Evaluamos en el centro ([math]\rho=0[/math]) y a la altura de la garganta ([math]z=2/3 H[/math]):

[math] v_z(\rho, z) = v_{max} \left( 1 - \frac{\rho^2}{R(z)^2} \right) \left( \frac{z}{H} \right)^{0,6} [/math]

Al ser [math]\rho=0[/math] se simplifican los cálculos:

[math] v_{tipica} = 4 \cdot 1 \cdot (0,667)^{0,6} = 4 \cdot 0,784 = 3,136 \, m/s [/math]

Sustituimos en el número de Reynolds:

[math] Re = \frac{1,225 \cdot 3,136 \cdot 60}{1,8 \cdot 10^{-5}} = 12.805.555 = 1,28 \cdot 10^7 [/math]

Régimen del flujo

Por convención:

• Si [math] Re \lt 2000 \to [/math] Laminar
• Si [math] Re \gt 2000 \to [/math] Turbulento

El régimen es extremadamente turbulento.

¿Por qué es tan alto? En general esto es debido a la escala "L" (diámetro). Aunque la viscosidad del aire intente ordenar el flujo, la inercia de mover una masa de 60 metros de ancho a 4 m/s es millones de veces superior. Lejos de ser un problema, esto es beneficioso, ya que el caos en el aire genera vórtices y mezclas constantes, las cuales ayudan a la disipación del calor.