Diferencia entre revisiones de «Torres de enfriamientos hiperbólicas (Grupo 46)»

Revisión del 17:55 2 dic 2025

1 Introducción

Las torres de enfriamiento hiperbólicas son estructuras icónicas de plantas de energía (nucleares y térmicas) diseñadas para disipar calor mediante convección natural. Su forma hiperbólica no es solo estética: proporciona resistencia estructural óptima frente al viento y maximiza el flujo de aire ascendente. Estas estructuras se popularizaron desde mediados del siglo XX y son un ejemplo perfecto de diseño ingenieril donde forma y función se unen.

2 Modelo geométrico del hiperboloide

Consideramos una torre de enfriamiento con altura total ([math]H[/math]), radio máximo en la base ([math]R_{\text{max}}[/math]), y radio mínimo ([math]R_{\text{min}}[/math]) (estrangulamiento) alcanzado a una altura ([math]h = \dfrac{2}{3}H[/math]) desde la base. La superficie de la torre es un hiperboloide de revolución de una hoja, descrito en coordenadas cartesianas como:

donde [math]a, c, z_0 \gt 0[/math] son parámetros a determinar en función de la geometría de la torre. Para nuestro modelo, usamos los siguientes parámetros de una torre típica:

2.1 Presion del viento

El viento ejerce una presión lateral sobre la torre. La velocidad del viento aumenta con la altura según la ley de potencia:

donde [math]V_0 = 18 \,\text{m/s}[/math] es la velocidad a altura de referencia [math]z_{\text{ref}} = 10 \,\text{m}[/math], y [math]\alpha = \dfrac{1}{7}[/math] es el exponente para terreno abierto.

La presión dinámica del viento es:

donde [math]\rho_{\text{aire}}[/math] es la densidad del aire estándar.

El campo vectorial de fuerza por unidad de superficie sobre la torre es:

donde [math]\vec{n}[/math] es el vector normal a la superficie (apuntando hacia el exterior).

2.2 Campo de temperaturas

Dentro de la torre, el aire caliente asciende por convección natural. Modelamos el campo de temperatura (en °C) como:

donde:

[math]T_{\text{base}} = 65^\circ \text{C}[/math]: temperatura en el centro de la base;

[math]\Delta T_z = 38^\circ \text{C}[/math]: caída de temperatura desde base hasta tope;

[math]\Delta T_\rho = 8^\circ \text{C}[/math]: variación radial de temperatura;

[math]n = 1.8[/math]: exponente de convección.

3 Determinación de los parámetros del modelo

3.1 Significado de cada parámetro

3.2 calculo de los parámetros

Para la determinación de los parámetros del modelo se utilizan lo datos iniciales y se realiza el cambio de la ecuación del hiperboloide a coordenadas cilíndricas.

Obteniendo la ecuación de la superficie:

El primer parámetro a determinar va a ser el de [math]z_0[/math]:

Conociendo que [math]z_0[/math] se alcanza a la altura [math]h = \dfrac{2}{3}H[/math] y que el valor de [math]H=150[/math] despejamos y obtenemos que:

Ahora se calcula el parámetro [math]a[/math] sabiendo que a la altura [math]z=100[/math] [math]ρ=30[/math] se sustituye en la ecuación obteniendo:

Por último se calcula el parámetro [math]c[/math], sustituyendo todos los valores previamente obtenidos y sabiendo que en la altura [math]z=0[/math] el radio es [math]ρ=55[/math]:

De la cual se despeja el parametro [math]c[/math]

3.3 Representación en MATLAB

En la siguiente imagen se muestra la visualización de la superficie de nuestra torre de enfriamiento en color gris a través de MATLAB. A la visualización se le han acompañado dos círculos azules en los planos [math]z=0[/math] y [math]z=150[/math] para una mejor visualización:

Visualización de la superficie parametrizada

% Parámetros
a  = 30;
c  = 65.079;
z0 = 100;
zmin = 0;
zmax = 150;

% Mallado de u y v
u = linspace(0, 2*pi, 100);
v  = linspace(zmin, zmax, 100);
[Mu, Mv] = meshgrid(u,v);

% Parametrización
rho = a * sqrt(1 + ((Mv - z0).^2) / (c^2));
theta = Mu;
z = Mv;

% Conversión a coordenadas cartesianas
X = rho .* cos(theta);
Y = rho .* sin(theta);
Z = z; 

% Gráfica sin mallado y en gris
surf(X, Y, Z, 'FaceColor', [0.7 0.7 0.7], 'EdgeColor', 'none');
hold on;

% Dibujar borde inferior (z = zmin)
Rmin = a * sqrt(1 + ((zmin - z0)^2) / (c^2));
x_min = Rmin * cos(u);
y_min = Rmin * sin(u);
plot3(x_min, y_min, zmin*ones(size(u)), 'b', 'LineWidth', 2);

% Dibujar borde superior (z = zmax)
Rmax = a * sqrt(1 + ((zmax - z0)^2) / (c^2));
x_max = Rmax * cos(u);
y_max = Rmax * sin(u);
plot3(x_max, y_max, zmax*ones(size(u)), 'b', 'LineWidth', 2);

axis equal;
xlim([-60, 60]);
ylim([-60, 60]);
xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
title('Hiperboloide de revolución de una hoja');
view(3)
hold off;

4 Estructuras hiperbólicas famosas y ventajas frente a una superficie cilíndrica

4.1 Razón por la cual se usan hiperboloides

La forma hiperbólica en torres de enfriamiento ofrece ventajas estructurales y funcionales frente a la superficie cilíndrica. En cuanto a la eficiencia aerodinámica, la curvatura hiperbólica favorece el tiro natural, generando una corriente ascendente de aire que mejora la refrigeración sin necesidad de ventiladores. Por otra parte, la geometría hiperbólica distribuye mejor las cargas de viento, reduce los esfuerzos y aumenta la estabilidad frente a condiciones climáticas extremas. Además, su forma permite construir estructuras altas y resistentes con menor cantidad de hormigón y acero, lo que reduce costes. Por último, presenta una mayor superficie de intercambio térmico, es decir, el ensanchamiento superior incrementa el área de contacto entre aire y agua, mejorando la transferencia de calor.

4.2 Torres de enfriamiento hiperbólicas famosas

4.3 Otras estructuras hiperboloides en arquitectura e ingeniería

Torre de Polibino

• Torres de Shújov (Rusia): Vladimir Shújov patentó en 1896 las primeras torres hiperboloides de celosía en acero, como la torre de Polibino y la famosa Torre Shabolovka en Moscú (1920–1922).

Aplicaciones arquitectónicas modernas:

• Antonio Gaudí utilizó geometrías hiperbólicas en la Sagrada Familia.

• Eduardo Torroja diseñó cubiertas hiperbólicas como la del Hipódromo de la Zarzuela en Madrid.

• Oscar Niemeyer y Le Corbusier también exploraron estas formas en proyectos emblemáticos.

Ejemplos contemporáneos: La Torre de Kōbe (Japón) y diversas cubiertas culturales emplean hiperboloides por su estética y eficiencia estructural.

La forma hiperbólica en torres de enfriamiento no es un mero capricho estético, sino una solución técnica que combina eficiencia térmica, resistencia estructural y economía de materiales. Su éxito ha trascendido el ámbito energético, inspirando a arquitectos e ingenieros en la creación de torres, cubiertas y estructuras icónicas que aprovechan las propiedades geométricas del hiperboloide para unir funcionalidad y belleza.

5 Análisis de presión del viento

5.1 Opción 1: Análisis de presión del viento

Supón que el viento sopla en dirección del vector \(-\frac{1}{\sqrt{2}}\left(\vec{i}+\vec{j}\right)\) (45° desde el noreste).

5.2 Calcula y representa el campo escalar \(P(z)\) como un mapa de colores sobre la mitad de la torre expuesta al viento. Añade colorbar. ¿Dónde es máxima la presión?

Para calcular el campo escalar \(P(z)\) sustituimos el campo escalar de la velocidad del viento [math]V(z) = V_0 \left( \frac{z}{z_{\text{ref}}} \right)^{\alpha}[/math] en la presión dinámica del viento [math]P(z) = \frac{1}{2} \rho_{\text{aire}} V(z)^2[/math], obteniendo la expresión: [math]P(z) = \frac{1}{2} \rho_{\text{aire}} \Bigl( V_0 \left( \frac{z}{z_{\text{ref}}} \right)^{\alpha}\Bigr)^2[/math]

Sustituyendo todos los parámetros tenemos: [math]P(z) = \tfrac{1}{2} \cdot 1{,}225 \cdot \Bigl( 18 \cdot \left(\tfrac{z}{10}\right)^{1/7} \Bigr)^2 \;\;\;\Rightarrow\;\;\; \boxed{P(z)\; \approx \; 102{,}787 \cdot z^{2/7}}[/math]

Como el exponente del campo escalar \(P(z)\) es [math]\tfrac{2}{7}\gt0 [/math] entonces la presión es directamente proporcional a la altura. Por tanto, la presión va a ser máxima en [math]z=150[/math]

¿Qué parte de la superficie está expuesta al viento?

Tenemos que el viento tiene como vector dirección [math]\vec{v}=-\frac{1}{\sqrt{2}}\left(\vec{i}+\vec{j}\right) [/math] y el vector normal \(\vec{n}\) a la superficie apunta hacia el exterior.

Si tomamos el plano [math]z=0[/math] y consideremos la circunferencia de centro en el origen: [math]x^2 + y^2 = r^2[/math]

En cualquier punto de esa circunferencia, el vector que va desde el centro hasta el punto es: [math]\vec{r} = (x, y)[/math]

Este vector radial siempre es perpendicular (normal) a la circunferencia, porque la circunferencia es un conjunto de puntos a una distancia fija del centro.

Un punto genérico de la circunferencia se puede escribir como: [math](x, y) = (r\cos(\theta),\ r\sin(\theta))[/math]

Entonces el vector radial(normal sin normalizar) es: [math]\vec{r} = (r\cos(\theta),\ r\sin(\theta))[/math]

Al normalizar el vector radial tenemos [math] \vec{r} = \frac{\vec{r}}{|\vec{r}|} \;\;\;\Rightarrow\;\;\; \vec{r} = \frac{r \cdot (\cos\theta, \sin\theta)}{r \cdot \sqrt{\cos^2\theta + \sin^2\theta}} \;\;\;\Rightarrow\;\;\; \vec{r} = (\cos\theta, \sin\theta) [/math]

Extrapolando del plano al espacio, el vector normal en [math]z=0[/math] es [math]\vec{n}=(cos(\theta),sin(\theta),0)[/math]

Según las propiedades del producto escalar, aplicamos la condición necesaria para identificar que parte de la superficie se ve afectada por el viento:

[math]\vec{v} \cdot \vec{n} \lt 0 \Rightarrow -\frac{1}{\sqrt{2}} (1,\,1,\,0) \cdot (cos(\theta),sin(\theta),0) \lt 0 \Rightarrow -\frac{1}{\sqrt{2}} \cdot (cos(\theta) + sin(\theta)) \lt 0 \Rightarrow cos(\theta) + sin(\theta) \gt 0 \Rightarrow \sqrt{2} \cdot \sin\left(\theta + \tfrac{\pi}{4}\right) \gt 0 \Rightarrow \sin\left(\theta + \tfrac{\pi}{4}\right) \gt 0 \Rightarrow 0 \lt \theta + \tfrac{\pi}{4} \lt \pi [/math]

[math]\Rightarrow \boxed{-\tfrac{\pi}{4} \lt \theta \lt \tfrac{3\pi}{4}} [/math]

En el siguiente recuadro se muestra como el campo escalar \(P(z)\) afecta sobre la mitad de la superficie expuesta al viento. También se incluye el código de Matlab utilizado para la representación.

% Parámetros
a  = 30;
c  = 65.079;
z0 = 100;
zmin = 0;
zmax = 150;

% Parámetros del viento
rho_aire = 1.225;         % densidad del aire [kg/m^3]
V0  = 18;            % velocidad referencia [m/s]
zref = 10;           % altura referencia [m]
alpha = 1/7;

% Mallado de u y v
u = linspace(0, 2*pi, 100);
v  = linspace(zmin, zmax, 100);
[Mu, Mv] = meshgrid(u,v);

% Parametrización
rho = a * sqrt(1 + ((Mv - z0).^2) / (c^2));
theta = Mu;
z = Mv;

% Conversión a coordenadas cartesianas
X = rho .* cos(theta);
Y = rho .* sin(theta);
Z = z; 

% Dirección del viento
vx = -1/sqrt(2);
vy = -1/sqrt(2);

% Vector normal
Nx = cos(theta);  
Ny = sin(theta);  

% Mascara: zona expuesta al viento
mask = (Nx*vx + Ny*vy) < 0;

% PRESIÓN DINÁMICA P(z)
Vz = V0 * (Z/zref).^alpha;   % velocidad del viento
P  = 0.5 * rho_aire .* Vz.^2;     % presión dinámica

% Aplicar máscara (ocultar lado no expuesto)
P(~mask) = NaN;

% Representación
figure(1)
surf(X, Y, Z, P, 'EdgeColor','none');
colormap(jet);
colorbar;
xlabel('x(m)'); ylabel('y(m)'); zlabel('z(m)');
title('Presión del viento sobre la mitad expuesta de la torre');
axis equal;
view(20,20);

figure(2)
surf(X, Y, Z, P, 'EdgeColor','none');
colormap(jet);
colorbar;
xlabel('x(m)'); ylabel('y(m)'); zlabel('z(m)');
title('Vista en planta');
axis equal;
view(2);

5.3 Calcula y representa el campo vectorial de fuerza \(\vec{F}\) sobre la superficie expuesta.

Así queda representado el campo vectorial de fuerzas sobre la superficie expuesta:

% Parámetros
a  = 30;
c  = 65.079;
z0 = 100;
zmin = 0;
zmax = 150;

% Parámetros del viento
rho_aire = 1.225;         % densidad del aire [kg/m^3]
V0  = 18;            % velocidad referencia [m/s]
zref = 10;           % altura referencia [m]
alpha = 1/7;

% Mallado de u y v
u = linspace(0, 2*pi, 40);
v  = linspace(zmin, zmax, 40);
[Mu, Mv] = meshgrid(u,v);

% Parametrización
rho = a * sqrt(1 + ((Mv - z0).^2) / (c^2));
theta = Mu;
z = Mv;

% Conversión a coordenadas cartesianas
X = rho .* cos(theta);
Y = rho .* sin(theta);
Z = z; 

% Dirección del viento
vx = -1/sqrt(2);
vy = -1/sqrt(2);

% Derivadas parciales
dX_dtt = -rho .* sin(theta);
dY_dtt = rho .* cos(theta);
dZ_dtt = zeros(size(theta));

dX_dz = a * (Z - z0) ./ (c^2 * sqrt(1 + ((Z - z0).^2 / c^2))) .* cos(theta);
dY_dz = a * (Z - z0) ./ (c^2 * sqrt(1 + ((Z - z0).^2 / c^2))) .* sin(theta);
dZ_dz = ones(size(Z)); 

% Normal exterior
Nx = dY_dtt .* dZ_dz - dZ_dtt .* dY_dz;
Ny = dZ_dtt .* dX_dz - dX_dtt .* dZ_dz;
Nz = dX_dtt .* dY_dz - dY_dtt .* dX_dz;

% Mascara: zona expuesta al viento
mask = (Nx*vx + Ny*vy) < 0;

% PRESIÓN DINÁMICA P(z)
Vz = V0 * (Z/zref).^alpha;   % velocidad del viento
P  = 0.5 * rho_aire .* Vz.^2;     % presión dinámica

% Fuerza F = -P * n
Fx = -P .* Nx;
Fy = -P .* Ny;
Fz = -P .* Nz;

% Aplicar máscara
Fx(~mask) = NaN;
Fy(~mask) = NaN;
Fz(~mask) = NaN;

% Representación
figure(1)
quiver3(X, Y, Z, Fx, Fy, Fz, 2, 'b');  
axis equal;
xlabel('x(m)'); ylabel('y(m)'); zlabel('z(m)');
title('Campo vectorial de fuerza F = -P(z) n sobre la superficie expuesta');
view(0,30);

figure(2)
quiver3(X, Y, Z, Fx, Fy, Fz, 2, 'b');   
axis equal;
xlabel('x(m)'); ylabel('y(m)'); zlabel('z(m)');
title('Vista en planta');
view(2);

5.4 Calcula la fuerza total y la fuerza por unidad de superficie que el viento ejerce sobre la torre integrando la presión sobre la mitad de superficie expuesta. Puedes hacerlo numéricamente o analíticamente. Puedes hacer esto numéricamente discretizando la superficie.

5.5 Compara con una torre cilíndrica de radio \(R = R_{\text{max}} = 55\ \text{m}\) y misma altura. Calcula la fuerza por unidad de superficie del viento sobre la mitad del cilindro expuesta. ¿Cuál estructura soporta menos fuerza lateral? ¿Cuál usa menos material (compara áreas superficiales)?

6 Campo de temperatura y transferencia de calor

6.1 Representación del campo escalar de temperatura en MATLAB

% apartado_2_4_torre_hiperbolica_animacion.m
% Torre hiperboloidal real con animación de corte horizontal

clear; close all; clc;

%% ------------------ PARÁMETROS FÍSICOS ------------------
Tbase = 57;         
Ttope = 22.43;
dTz   = Tbase - Ttope;   % = 34.57 °C
dTr   = 8;               % variación radial
n     = 1.8;

%% ------------------ GEOMETRÍA REAL ------------------
H = 150;           % altura de la torre
z0 = 100;          % centro del hiperboloide

a = 30;            % semieje radial mínimo
c = 68.079;        % semieje vertical
Rmax_base = 55;    % radio máximo en la base

%% ------------------ Superficie hiperboloide ------------------
NzS = 180;
Ntheta = 160;

zS = linspace(0, H, NzS);
theta = linspace(0, 2*pi, Ntheta);

[TH, Zsurf] = meshgrid(theta, zS);

% radio hiperboloidal exacto
Rsurf = a .* sqrt(1 + ((Zsurf - z0).^2)/c^2);

% coordenadas cartesianas
Xsurf = Rsurf .* cos(TH);
Ysurf = Rsurf .* sin(TH);

%% ------------------ TEMPERATURA EN LA SUPERFICIE ------------------
Tsurf = Tbase ...
      - dTz * (Zsurf/H).^n ...
      - dTr .* (1 - exp(-(Rsurf.^2)/(Rmax_base^2)));

Tmin = min(Tsurf(:));
Tmax = max(Tsurf(:));

%% ------------------ FIGURA ------------------
figure('Units','normalized','Position',[0.03 0.1 0.9 0.8]);
colormap(parula);

%% ----------- SUBPLOT IZQUIERDO: SECCIÓN VERTICAL ------------
subplot(1,2,1);

rho_vec = linspace(0, Rmax_base, 300);
z_vec   = linspace(0, H, 300);
[RR, ZZ] = meshgrid(rho_vec, z_vec);

Tmap = Tbase ...
     - dTz*(ZZ/H).^n ...
     - dTr*(1 - exp(-(RR.^2)/(Rmax_base^2)));


hImg = imagesc(rho_vec, z_vec, Tmap);
set(gca,'YDir','normal');
xlabel('\rho (m)'); ylabel('z (m)');
title('Sección vertical \rho-z');
clim([Tmin Tmax]);

hold on;
hLine = plot([0 Rmax_base],[0 0],'w--','LineWidth',2);

% Punto más caliente (rho=0, z=0)
hHot2D = plot(0, 0, 'ro', 'MarkerFaceColor','r', 'MarkerSize',8);
text(2, 3, 'Punto más caliente', 'Color','white', 'FontSize',10);

%% ----------- SUBPLOT DERECHO: HIPERBOLOIDE REAL 3D ---------------
subplot(1,2,2);

hSurf = surf(Xsurf, Ysurf, Zsurf, Tsurf);


shading interp;
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
title('Torre hiperboloidal en 3D con corte horizontal');
clim([Tmin Tmax]);
view(35,25);
hold on;

% Plano horizontal móvil
[XP, YP] = meshgrid(linspace(-Rmax_base, Rmax_base, 70));
ZP = zeros(size(XP));
hPlane = surf(XP, YP, ZP, 'FaceAlpha',0.25, 'EdgeColor','none');



colorbar;

%% ------------------ ANIMACIÓN ------------------
for k = 1:length(z_vec)
    zcur = z_vec(k);

    % actualizar línea en sección vertical
    set(hLine, 'YData', [zcur zcur]);

    % actualizar plano en 3D
    set(hPlane, 'ZData', zcur * ones(size(ZP)));

    drawnow;
    pause(0.03);
end

6.2 Calculo del gradiente de temperatura [math]\Delta T ^\circ \text{C}[/math]

texto alternativo

texto alternativo

6.3 Representación del campo vectorial gradiente sobre el hiperboloide

%% APARTADO 2.5 — Sección vertical con campo completo -∇T en el interior
clear; close all; clc;

%% ---- Parámetros de la torre ----
H  = 150;
z0 = 100;
a  = 30;
c  = 98.079;

%% ---- Mallado vertical ----
Nr = 160;   % resolución radial
Nz = 300;   % resolución vertical

rho = linspace(0, 60, Nr);
z   = linspace(0, H, Nz);

[RR, ZZ] = meshgrid(rho, z);

%% ---- Radio de la pared ----
r_wall = a .* sqrt(1 + ((z - z0).^2)/c^2);  % r(z)
r_wall_2D = repmat(r_wall(:), 1, Nr);       % expandido a matriz

%% ---- Máscara interior ----
inside = RR <= r_wall_2D;

%% ---- Gradiente EN CILÍNDRICAS (TU FÓRMULA EXACTA) ----
Grad_r = -(912 * (RR.^2) .* exp(-(RR.^2)./(3249 - RR.^2))) ./ (3249 - RR.^2).^2;
Grad_z = -(0.00829 * ZZ.^0.8);

% eliminar valores fuera de la torre
Grad_r(~inside) = NaN;
Grad_z(~inside) = NaN;

%% ---- FIGURA ----
figure('Color','w','Position',[100 100 850 700]); hold on;

% Contorno de la figura (lado derecho e izquierdo)
plot(r_wall, z, 'k', 'LineWidth', 2);
plot(-r_wall, z, 'k', 'LineWidth', 2);

% Eje central
plot(0*z, z, 'k--', 'LineWidth', 1);

%% ---- Campo completo de flechas ----
skip = 6;  % densidad del campo

quiver( RR(1:skip:end,1:skip:end), ...
        ZZ(1:skip:end,1:skip:end), ...
        Grad_r(1:skip:end,1:skip:end), ...
        Grad_z(1:skip:end,1:skip:end), ...
        'k');

%% ---- Estética ----
xlabel('\rho (m)');
ylabel('z (m)');
title('Campo completo -\nablaT en sección vertical (delimitado por la torre)');
xlim([-max(r_wall) max(r_wall)]);
ylim([0 H]);
axis equal;
grid on;
set(gca,'FontSize',12);


%% ---- Estética ----
xlabel('\rho (m)');
ylabel('z (m)');
title('Campo completo -\nablaT en sección vertical (delimitado por la torre)');
xlim([-max(r_wall) max(r_wall)]);
ylim([0 H]);
axis equal;
grid on;
set(gca,'FontSize',12);

6.4 2.6

texto alternativo