Torres de enfriamiento hiperbólicas (grupo 35)

Trabajo realizado por estudiantes
Título	Torres de enfriamiento hiperbólicas Grupo 35
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Miguel Álvarez Penabad Alejandro Jiménez García Pedro Miguel Jaume Méndez Rodrigo Martínez Villén Noah González Becerra
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción:

Las torres de enfriamiento hiperbólicas son elementos característicos debido a su forma y tamaño, tanto de las centrales nucleares, como de las termoeléctricas. Su geometría, tiene una explicación estructural y física: la curvatura hiperbólica les brinda gran estabilidad frente a las presiones causadas por el viento y también favorece el ascenso del aire caliente, mediante convección natural. Desde la segunda mitad del siglo XX, este tipo de torres se ha vuelto muy común y constituye un claro ejemplo de cómo la ingeniería combina de manera eficiente la forma con la función.

2 Geometría hiperbólica:

En este artículo, se analizará un modelo estándar de torre de enfriamiento, definido por unas magnitudes dadas, una altura máxima ([math]H[/math]), un radio máximo en la base ([math]Rmáx[/math]), y un radio mínimo ([math]Rmín[/math]), que se alcanza a una altura dada [math]Zo[/math], que es igual a, [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math]. La superficie de la torre, es la de un hiperboliode de revolución de una hoja, que viene descrita por la siguiente ecuación en coordenadas cartesianas:

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

Siendo [math]a, c, z_0\gt0[/math] unos valores a calcular

Siendo en este caso [math]Rmáx=55m[/math] ; [math]Rmín=30m[/math] ; [math]H=150m[/math]

Para completar la ecuacion de la torre hiperbólica, se debe obtener los valores de: [math]Zo[/math], [math]c[/math] y [math]a[/math].

El [math]Zo[/math] marca, como ya se ha comentado la altura de [math]Rmín[/math], que con la operación matemática:

[math]Zo[/math] = [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math]

Se obtiene el valor de [math]Zo(H=150)=100[/math].

A continuación, para hallar los valores de '[math]c[/math]' y '[math]a[/math]' se debe cambiar el sistema de referencia, de sistema cartesiano a sistema cilíndrico, para ello, se ultilizarán las siguientes fórmulas:

[math]x=\rho cos(\theta)[/math]

[math]y=\rho sin(\theta)[/math]

[math]z=z[/math]

Dando como resultado la siguiente ecuación:

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

Para hallar el valor de '[math]a[/math]', se sistituye en la ecuación, con:[math]\rho = Rmín[/math], y por lo tanto [math]Z = Zo[/math]. Que sustituyendo queda:

[math]\dfrac{30^2}{a^2}-\dfrac{(100-100)^2}{c^2}=1[/math]

Y de donde se saca el valor de la incógnita '[math]a=30[/math]'

Con este valor, y sustituyendo en la misma ecuación se obtiene el valor de '[math]c[/math]', pero esta vez, con [math]\rho = Rmáx[/math], y por lo tanto [math]Z = 0[/math]:

[math]\dfrac{55^2}{a^2}-\dfrac{(0-100)^2}{c^2}=1[/math]

De donde se saca el valor de la incógnita '[math]c = 65,0791 [/math]'

Siendo la ecuación final, con los valores '[math]a[/math]' y '[math]c[/math]' ya sustituidos la siguiente:

[math]\dfrac{ρ^2}{30^2}-\dfrac{(z-z_0)^2}{65,0791^2}=1[/math]

3 Visualizació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 en los planos [math]z=0[/math] y [math]z=150[/math] para una mejor visualización:

Superficie parametrizada de la torre de enfriamiento

% Parámetros de la torre:
a = 30;
c = 65.079;
z_0 = 100;
R_max = 55;
R_min = 30;
z_min = 0;
z_max = 150;

% Coordenadas cilíndricas:
theta = linspace(0, 2*pi, 100);
z = linspace(z_min, z_max, 100);
[THETA, Z] = meshgrid(theta, z);

% Definición de ro:
ro = sqrt(a^2 * (1 + ((Z - z_0).^2 / c^2)));

% Asegurar que los valores están en el dominio:
ro = min(max(ro, R_min), R_max);

% Conversión a coordenadas cartesianas:
X = ro .* cos(THETA);
Y = ro .* sin(THETA);

% Crear la figura:
figure;
surf(X, Y, Z, 'FaceColor', [0.8, 0.8, 0.8], 'EdgeColor', 'none');
hold on;

% Dibujar los círculos en los planos z=0 y z=120:
circle_theta = linspace(0, 2*pi, 100);
circle_x_z0 = R_max * cos(circle_theta);
circle_y_z0 = R_max * sin(circle_theta);
circle_x_z120 = sqrt(a^2 * (1 + ((z_max - z_0)^2 / c^2))) * cos(circle_theta);
circle_y_z120 = sqrt(a^2 * (1 + ((z_max - z_0)^2 / c^2))) * sin(circle_theta);

% Mostrar los círculos:
plot3(circle_x_z0, circle_y_z0, z_min * ones(size(circle_theta)), 'k', 'LineWidth', 1);
plot3(circle_x_z120, circle_y_z120, z_max * ones(size(circle_theta)), 'k', 'LineWidth', 1);

% Configurar límites y etiquetas:
axis equal;
xlim([-60, 60]);
ylim([-60, 60]);
zlim([z_min, z_max]);
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Superficie Hiperbólica en Coordenadas Cilíndricas');

% Mostrar la gráfica
view(3);
grid on;
hold off;

4 Análisis de presión del viento:

El viento es un agente externo que ejerce presiones laterales que varían a lo largo de la superficie en función de la altura. La velocidad escalar se puede representar de la siguiente forma:

[math]V(z)=V_0(\dfrac{z}{Zref})^α[/math]

Sea [math]V_0[/math] la velocidad de referencia a la altura [math]Zref[/math] que en nuestro modelo podemos fijar en [math]18m/s[/math].

Sea [math]Zref[/math] la altura de referencia sobre el suelo, que generalmente equivale a [math]10m[/math].

Sea [math]α[/math] un exponente que depende del terreno, en nuestro modelo tendrá un valor de [math]0,14[/math].

Utilizando la velocidad del viento previamente mencionada, la presión del viento sobre la superficie de la torre se puede representar como:

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

Sea ρ la densidad del aire estándar, que tiene un valor de [math]1,2Kg/m^3[/math].

Sustituyendo en la fórmula obtenemos que la presión en función de la altura es:

[math]P(z)=\dfrac{1}{2}ρV(z)^2=\dfrac{1}{2}ρ(V_0(\dfrac{z}{Zref})^α)^2=\dfrac{1,2}{2}(18(\dfrac{z}{10})^{0,14})^2≈102,0226z^{0,28}[/math]

Suponiendo que la torre está sometida a una presión de viento paralelo al vector [math]-\frac{1}{\sqrt{2}}(\overrightarrow{i}+\overrightarrow{j})[/math] desde el noreste hacia suroeste.

Mapa de presiones del viento

%% Parámetros del problema
H = 150; 
Rmax = 55;
Rmin = 30;
rho_air = 1.225;
V0 = 18; 
zref = 10;
alpha = 1/7;

% Altura del estrangulamiento
h = 2*H/3;

% Parámetros del hiperboloide
a = Rmin;          % radio mínimo
z0 = h;            % lugar del estrangulamiento
c = sqrt( (Rmax^2/a^2 - 1) * (0 - z0)^2 );  % derivado de ecuaciones geométricas

%% Mallado de la superficie
nz = 200; ntheta = 300;
z = linspace(0, H, nz);
theta = linspace(0, 2*pi, ntheta);
[TH, ZZ] = meshgrid(theta, z);

% Radio en cada altura
RR = a * sqrt(1 + ((ZZ - z0).^2)/c^2);

% Coordenadas cartesianas
XX = RR .* cos(TH);
YY = RR .* sin(TH);

%% Selección de la mitad expuesta al viento
w = [-1/sqrt(2), -1/sqrt(2), 0];  % dirección del viento
normal_component = cos(TH) + sin(TH); 
mask = normal_component > 0;

%% Presión del viento
VZ = V0 * (ZZ / zref).^alpha;
PZ = 0.5 * rho_air .* VZ.^2;

% Aplicar máscara: solo la mitad expuesta
PX = PZ;
PX(~mask) = NaN;

%% Gráfica
figure
surf(XX, YY, ZZ, PX, 'EdgeColor', 'none')
colormap jet
colorbar
xlabel('x (m)')
ylabel('y (m)')
zlabel('z (m)')
title('Mapa de presiones del viento sobre la torre (mitad expuesta)')
view(60, 30)

4.1 Campo vectorial de fuerza:

El campo vectorial de fuerza por unidad de superficie sobre la torre se puede expresar como [math]\overrightarrow{F}(x,y,z)=-P(z) \overrightarrow{n}[/math]

El campo vectorial de fuerza resultante es: [math]\overrightarrow{F}(x,y,z)=(-P(z) r(z) cos(\theta),-P(z) r(z) sin(\theta),P(z) r(z) r'(z))[/math]
Siendo [math]r(z)[/math] la superficie parametrizada.

A continuación se muestra la representación en formato MATLAB:

Campo vectorial de fuerza sobre la superficie de la torre

% REPRESENTACIÓN DEL CAMPO DE FUERZA F SOBRE LA SUPERFICIE EXPUESTA
 R_max = 55; % Radio de la base de la torre (m)
 R_min = 30; % Radio de estrangulamiento (m)
 z_max = 150; % Altura de la torre (m)
 z_0 = 100; % Centro hiperboloide (m)
 c = 65.079; % Parámetro c
 v_0 = 18; % Velocidad del viento a la altura de referencia (m/s)
 z_ref = 10; % Altura de referencia (m)
 alfa = 1/7; % Exponente para terreno abierto
 densidad_aire = 1.225; % Densidad del aire (Kg/m^3) a nivel del mar y a 15ºC
 % Creación de la malla:
 angulo_theta = linspace(0, 2*pi, 100); % Ángulo alrededor del eje Z
 altura = linspace(0, z_max, 200);
 [X, Z] = meshgrid(angulo_theta, altura); % Malla
 % Hiperbole parametrizado:
 coord_x = R_min * cos (X).* sqrt(1 + ((Z - z_0) / c).^2);
 cood_y = R_min * sin(X).* sqrt(1 + ((Z - z_0) / c).^2);
 coord_z = Z;
 % Corte para el eje de la torre:
 corte = (coord_x + cood_y) <= 0;
 coord_x_corte = coord_x .* corte;
 coord_y_corte = cood_y .* corte;
 coord_z_corte = coord_z .* corte;
 % Flechas (apuntando al centro del eje Z) en color azul:
 F_x_centro = -coord_x_corte; % Representación del campo en componente x hacia el centro
 F_y_centro = -coord_y_corte; % Representación del campo en componente y hacia el centro
 F_z_centro = 0 * coord_z_corte; % No hay componente en z por ser paralelas al eje Z
 % Flechas según la presión del viento:
 factor_escala = coord_z_corte / max(coord_z_corte(:));
 factor_escala(isnan(factor_escala)) = 0; % Valores NaN en la base
 F_x_centro = F_x_centro .* factor_escala;
 F_y_centro = F_y_centro .* factor_escala;
 % Evitar puntos nulos:
 puntos_validos = corte & (coord_x_corte ~= 0 | coord_y_corte ~= 0);
 % Dibujo de las flechas:
 figure;
 hold on;
 % Gráfico:
 quiver3(coord_x_corte(puntos_validos), coord_y_corte(puntos_validos), coord_z_corte(puntos_validos), F_x_centro(puntos_validos), F_y_centro(puntos_validos), F_z_centro(puntos_validos), 2, 'Color', [0.5 0.7 1]); % Azul claro
 xlabel('x (m)'); 
 ylabel('y (m)');
 zlabel('z (m)');
 title('Campo vectorial de la presión del viento sobre la mitad expuesta');
 axis equal; % Para mantener las proporciones iguales
 grid on;
view(3); % Para la vista en 3D
hold off;

4.2 Fuerza total y fuerza por unidad de superficie que se ejerce sobre la torre:

Para calcular la fuerza total sobre la mitad de la torre se parte de la expresión de la presión [math]P(z)[/math], el perfil de viento [math]V(z)[/math] y el radio lineal [math]R(z)[/math]:

[math]P(z)=\frac{1}{2} \rho V(z)^{2}[/math] [math]V(z)=V_{0} (\frac{z}{z_{ref}})^{\alpha}[/math] [math]R(z) = R_{\min} + \frac{\Delta R}{H} z[/math]

La integración respecto de [math]\theta[/math] en el dominio de media torre da un valor de 2. Sustituyendo en la integral de [math]F[/math]:

[math]F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\int_0^H z^{2\alpha} \left( R_{\min} + \frac{\Delta R}{H}z \right)\, dz[/math]

Integrando quedaría la siguiente función:

[math]F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min}\,\frac{H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R}{H}\,\frac{H^{2\alpha+2}}{2\alpha+2}\right][/math]

Calculando con los datos anteriormente obtenidos, la fuerza total sería:

[math]F \approx 4.423835 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN})[/math]

[math]F = 4.423835 \times 10^{6}\ \text{N}[/math]

Es necesario calcular el área de la mitad de la torre de enfriamiento, la superficie lateral de revolución correspondiente a la mitad frontal es:

[math]A_{\text{mitad}}=\pi \int_0^{H} R(z)\, \sqrt{1 + [R'(z)]^{2}}\, dz[/math]

Dado que

[math]R'(z)=\frac{\Delta R}{H},[/math]

es constante, el factor geométrico puede sacarse fuera de la integral:

[math]A_{\text{mitad}}=\pi \sqrt{1 + R'^2}\int_0^{H} R(z)\, dz[/math]

La integral del radio es:

[math]\int_0^{H} R(z)\, dz=\int_0^{H} \left( R_{\min} + \frac{\Delta R}{H} z \right) dz=R_{\min} H + \frac{\Delta R}{2}H[/math]

Finalmente:

[math]{A_{\text{mitad}}=\pi H \sqrt{1 + \left( \frac{\Delta R}{H} \right)^2}\left( R_{\min} + \frac{\Delta R}{2} \right)}[/math]

[math]A_{\text{mitad}} \approx 2.0303909708635903 \times 10^{4}\ \text{m}^2,\qquad(\approx 2.03\times 10^{4}\ \text{m}^2)[/math]

La fuerza por unidad de superficie es:

[math]P = \frac{F}{A_{\text{mitad}}}[/math]

Con los valores obtenidos:

[math]P \approx 217.881\ \text{Pa} [/math]

4.3 Comparación torre cilíndrica con hiperbólica:

Para calcular la fuerza por unidad de superficie sobre la mitad del cilindro cuya area es:

[math]\frac{1}{2}A_{cilindro} = πr(r+h)[/math] [math]\frac{1}{2}A_{cilindro} = 35421'45713 m^2[/math]

La fuerza que recibe el cilindro es:

[math]F = \rho \, V_{0}^{2} \, Z_{\text{ref}}^{-2\alpha} \, (\frac{R \times \, H^{2\alpha + 1}}{2\alpha + 1})[/math]

[math]F = 5349039'287 N [/math]

[math]P = \frac{F}{S}[/math] [math]P = 151'0112717 Pa[/math]

Observamos que [math] P_{hiperboloide}\gt P_{cilindro}[/math] lo que nos indica que el hiperboloide soporta más fuerza lateral.

En cuanto a material si comparamos las dos áreas superficiales, observamos que el cilindro consume más material [math]A_{cilindro}\gt A_{hiperboloide}[/math] y es por ello que podemos concluir que a nivel coste-utilidad el hiperboloide es favorable

5 Campo de temperatura y transferencia de calor

En este apartado vamos a estudiar la variación de la temperatura dentro de la torre, donde el aire caliente asciende por convección natural y se puede describir el campo escalar de temperatura (en °C) como: [math]T(\rho,z)=T_{base}-\Delta T_{z}*(\frac{z}{H})^{n}-\Delta T_{\rho}*(1-e^{-\frac{\rho^{2}}{R_{max}^{2}-\rho^{2}}})[/math] siendo:

• [math]\Delta T_{z}=38(°C)[/math] , correspondiente a la temperatura en el centro de la base;

• [math]T_{base}=65 (ºC)[/math] , es la caida de temperatura desde el tope hasta la base;

• [math]\Delta T_{\rho}=8(ºC)[/math] , variación radial de temperatura;

• [math]n=1,8[/math], es el exponente de convección.

Además, la temperatura del tope de la torre se puede calcular como [math]T_{tope}= T_{base}-\Delta T_{z}[/math]. En nuestro caso, [math]T_{tope}= 27 (ºC).[/math]

5.1 Representación del campo de temperatura respecto a un plano vertical:

Sección vertical del campo de temperatura

% REPRESENTACIÓN DEL CAMPO ESCALAR DE TEMPERATURA POR UN CORTE VERTICAL POR EL EJE DE SIMETRÍA
 % Parámetros:
 R_max = 55; % Base (m)
 R_min = 30; % Estrangulamiento (m)
 z_max = 150; % Altura de la torre (m)
 c = 65.079; % Parámetro calculado
 centro_hiperboloide = 100; % Centro del hiperboloide (m)
 temp_base = 65; % Temperatura de la base (°C)
 delta_temp_vertical = 38; % Caida de temperatura desde el tope hasta la base (°C)
 temp_top = temp_base - delta_temp_vertical; % Temperatura del tope de la torre (°C)
 delta_temp_radial = 8; % Variación de temperatura radial (°C)
 exponente_conveccion = 1.8; 
 figure( 'Name', 'Sección vertical del campo de temperatura');
 z = linspace(0, z_max,300); 
 x = linspace(-R_max - 10, R_max + 10, 300);
 [X, Z] = meshgrid(x, z); % Malla
 R = sqrt((1 + ((Z - centro_hiperboloide) / c).^2) * R_min^2); % Radio del hiperboloide
 r = abs(X); % Distancia radial
 T = temp_base - delta_temp_vertical * (Z / z_max).^exponente_conveccion - delta_temp_radial * (1 - exp(-r.^2 ./ (R_max^2 - r.^2))); % Campo escalar de temperatura
 mask = abs(X) <= R; % Limitar el hiperboloide
 T(~mask) = NaN; % No mostrar fuera de la silueta
 % Dibujo del campo de temperatura:
 contourf(X, Z, T, 55, 'LineColor', 'none'); 
 colormap('jet');
 colorbar;
 caxis([temp_top, temp_base]); % Escala de colores
 hold on;
 % Dibujo de la silueta de la torre:
 x_sil_left = -R; % Lado izquierdo de la silueta
 x_sil_right = R; % Lado derecho de la silueta
 z_sil = z; % Coordenadas verticales
 plot(x_sil_left, z_sil, 'k-', 'LineWidth', 2); % Silueta izquierda
 plot(x_sil_right, z_sil, 'k-', 'LineWidth', 2); % Silueta derecha
 % Configuración de la gráfica:
 xlabel('x (m)');
 ylabel('z (m)');
 title('Campo escalar de temperatura (°C) corte vertical por su eje');
 xlim([-R_max - 10, R_max + 10]);
 ylim([0, z_max]);
 grid on;
 view(2); % Vista en 2D
 hold off;

Como se puede observar, la zona donde la masa de aire se encuentra más caliente corresponde con la base de la torre.

5.2 Gradiente y Representación del campo de temperatura en secciones transversales

[math]T(\rho,z) = 65-38(\frac{z}{150})^{1,8}-8(1-e^{\frac{-p^2}{55^2-p^2}}))[/math]

[math]\bigtriangledown T(\rho,z) = \frac{\partial T}{\partial \rho}e_{\rho} + \frac{1}{\rho}\frac{\partial T}{\partial \theta}e_{\theta} + \frac{\partial T}{\partial z}e_{z}[/math]

[math] \nabla T(\rho, z) = -\frac{16\rho}{3025} e^{-\frac{\rho^2}{55^2}} \mathbf{e}_{\rho} - 0.456 \left(\frac{z}{150}\right)^{0.8} \mathbf{e}_{z}[/math]

Campo de temperatura en secciones transversales

       
% SECCIONES TRANSVERSALES DE LA TORRE      
% Parámetros:
R_max = 55; % Radio de la base de la torre
R_min = 30; % Radio de estrangulamiento
c = 65.079; % Parámetro calculado en el primer apartado
centro_hiperboloide = 100; % 2/3 de la altura total de la torre
z_max = 150; % Altura máxima de la torre
temp_base = 65; % Temperatura en la base de la torre
delta_temp_vertical = 38; % Diferencia de temperatura vertical
temp_top = temp_base - delta_temp_vertical; % Temperatura del tope de la torre
delta_temp_radial = 8; 
expo_con = 1.8; % Exponente de convección

% Creación de la malla para la figura
angulo_theta = linspace(0, 2*pi, 100); % Ángulo alrededor del eje z
altura_malla = linspace(0, z_max, 50); % Coordenada vertical
alturas_lineas = 0:20:z_max; % Altura de las lineas discontinuas

% Nombre del GIF
filename = 'animacion.gif';

% Figura
figure;

for i = 1:length(altura_malla)
    % Altura de la capa actual:
    altura_actual = altura_malla(i);

    % Radio en la altura actual
    radio_capa = sqrt((1 + (altura_actual - centro_hiperboloide)^2 / c^2) * R_min^2);

    % Coordenadas de la capa
    coord_x_capa = radio_capa * cos(angulo_theta);
    coord_y_capa = radio_capa * sin(angulo_theta);
    coord_z_capa = altura_actual * ones(size(coord_x_capa));

    % Temperatura de la capa actual:
    temp_capa = temp_base - delta_temp_vertical * (altura_actual / z_max)^expo_con - delta_temp_radial * (1 - exp(-radio_capa^2 / (R_max^2 - radio_capa^2)));

    % Dibujo de la capa
    fill3(coord_x_capa, coord_y_capa, coord_z_capa, temp_capa, 'EdgeColor', 'none');
    hold on;

    % Lineas discontinuas para las diferentes alturas
    for altura_linea = alturas_lineas
        radio_linea = sqrt((1+ (altura_linea - centro_hiperboloide)^2 / c^2) *R_min^2); % Radio en cada altura
        coord_x_linea = radio_linea * cos(angulo_theta);
        coord_y_linea = radio_linea * sin(angulo_theta);
        plot3(coord_x_linea, coord_y_linea, altura_linea * ones(size(coord_x_linea)), 'm--', 'LineWidth',1.5);
    end

    % Mostrar el contador en el titulo
    title(sprintf('Campo de temperatura a z = %.2f m', altura_actual));

    % Configuración gráfica
    colormap('jet');
    colorbar;
    caxis([temp_top, temp_base]); % Escala de colores
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('z (m)');
    xlim([-R_max, R_max]);
    ylim([-R_max, R_max]);
    zlim([0, z_max]);
    view(3); % Vista 3D

    % Captura del fotograma para el GIF
    drawnow; % Asegura que la figura se actualice
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);

    % Guardar en GIF
    if i == 1
        imwrite(imind, cm, filename, 'gif', 'Loopcount', inf, 'DelayTime', 0.15);
    else
        imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.15);
    end

    hold off;
end

5.3 Gradiente de la temperatura:

[math]\nabla T=(-8+8e^{\frac{-\rho^{2}}{55^{2}-\rho^{2}}})e_{\rho}-0.456(\frac{Z}{150})^{0.8}e_{z}[/math]

Al hacer el gradiente [math]\nabla T[/math] del campo escalar de temperatura obtenemos un campo vectorial que se puede visualizar como:

Sección vertical del campo de temperatura

 % REPRESENTACIÓN DEL GRADIENTE EN LA SECCIÓN VERTICAL
 % Parámetros:
 R_max = 55; % Radio de la base de la torre
 R_min = 30; % Radio de estrangulamiento
 z_max = 150; % Este es el mismo valor que H, altura de la torre
 c = 65.079; 
 z_0 = 100; % 2/3 de la altura total de la torre
 temp_base = 65; % Temperatura de la base
 delta_temp_vertical = 38; % Diferencia de temperatura de la base y el tope de la torre
 temp_top = temp_base - delta_temp_vertical; % Temperatura del tope
 delta_temp_radial = 8; 
 exponente_conveccion = 1.8; 
 % Creación de la malla:
 z = linspace(0, z_max,50);
 x = linspace(-R_max - 100, R_max + 100, 50);
 [X, Z] = meshgrid(x,z); 
 R = sqrt((1 + ((Z - z_0) / c).^2) * R_min^2); % Radio del hiperboloide
 r = abs(X); % Distancia radial
 % Campo de temperatura:
 T = temp_base - delta_temp_vertical * (Z / z_max).^exponente_conveccion - delta_temp_radial * (1 - exp(-r.^2 ./ (R_max^2 - r.^2)));
 % Limitación del hiperboloide:
 mask = abs(X) <= R; 
 T(~mask) = NaN; % No mostrar lo que hay fuera de la torre
 % Calculo del gradiente:
 [Tx,Tz] = gradient(T, x(2) - x(1), z(2) - z(1)); % Derivada del campo escalar 
 % Figura:
 figure('Name', 'Gradiente de temperatura en la sección vertical');
 % Silueta de la torre:
 angulo_theta = linspace(0, 2*pi, 100);
 x_sil_left = -R; % Lado izquierdo
 x_sil_right = R; % Lado derecho
 z_sil = z; % Vertical
 plot(x_sil_left, z_sil, 'g-', 'LineWidth',2); % Lado izquierdo de la torre
 hold on;
 plot(x_sil_right, z_sil, 'g-', 'LineWidth',2); % Lado derecho de la torre
 % Representación de las flechas en color azul:
 quiver(X(mask), Z(mask), Tx(mask), Tz(mask), 'b', 'LineWidth', 1, 'MaxHeadSize',3);
 % Configuración de la gráfica en 2D:
 xlabel('x (m)');
 ylabel('z (m)');
 title('Gradiente de temperatura en sección vertical');
 xlim([-R_max - 10, R_max + 10]);
 ylim([0, z_max]);
  grid on;
 axis equal; % La proporción de los objetos igual
 view(2); % Vista en 2D
 hold off;

5.4 Representar superficies isotermicas:

Superficies isotermicas

% SUPERFICIE ISOTÉRMICAS
 % Parámetros:
 radioBase = 30; % Radio de estrangulamiento
 alturaCentro = 100; % 2/3 de la altura total de la torre
 alturaTotal = 150; % Altura de la torre
 radioMaximo = 55; % Radio de la base
 radioHiperboloide = 65.079; % Parámetro c
 tempBase = 65; % Temperatura de la base de la torre
 tempCima = 27; % Temperatura del tope de la torre
 variacionTempAltura = tempBase - tempCima; % Delta_temp_vertical, diferencia de temperatura entre la base y el tope de la torre
 variacionTempRadial = 8; % Delta_temp_radial
 expConveccion = 1.8; % Exponente de convección
  
 % Discretización optimizada:
   theta = linspace(0, 2*pi, 50); % n menor para una mejor visualización
   z = linspace(0, alturaTotal, 50); % n igual a la de arriba por las mismas razones
  [Theta, Z] = meshgrid(theta, z); % Creación de la malla
 
 % Radio rho en función de z:
  Rho = sqrt(radioBase^2 * (1 + ((Z - alturaCentro).^2) / radioHiperboloide^2));
  X = Rho .* cos(Theta);
  Y = Rho .* sin(Theta);
 
 % Campo de temperaturas en 3D:
  campoTemp = tempBase - variacionTempAltura * (Z / alturaTotal).^expConveccion ...
          - variacionTempRadial * (1 - exp(-Rho.^2 ./ (radioMaximo^2 - Rho.^2))); % Campo escalar de temperatura
  
 % Niveles isotérmicos:
  nivelesTemp = linspace(tempCima, tempBase, 20); % Menos niveles isotérmicos
  
%% Calculo del gradiente térmico:
  [gradTempX, gradTempZ] = gradient(campoTemp, mean(diff(X(1,:))), mean(diff(Z(:,1))));
  gradTempY=gradTempX;
  figure;
  filename = 'animacion.gif'; % Para guardar la animación como gif
  for nivel = 1:length(nivelesTemp)

    tempNivel = nivelesTemp(nivel);
    mascaraIso = abs(campoTemp - tempNivel) < 0.5; % Margen ±0.5 °C

    % Gráfica:
    contour3(X, Y, Z, 6, 'm'); % m color en magenta
    xlabel('X (m)');
    ylabel('Y (m)');
    zlabel('Z (m)');
    hold on;

    % Extracción de  puntos isotérmicos:
    xIso = X(mascaraIso);
    yIso = Y(mascaraIso);
    zIso = Z(mascaraIso);
    coloresIso = campoTemp(mascaraIso);

    % Superficie isotérmica:
    scatter3(xIso, yIso, zIso, 10, coloresIso, 'filled');
    colormap('jet');
    colorbar;
    clim([tempCima, tempBase]);

    % Gradientes:
    quiver3(xIso, yIso, zIso, ...
        gradTempX(mascaraIso), gradTempY(mascaraIso), gradTempZ(mascaraIso), ...
        'm', 'AutoScaleFactor', 0.5);

    drawnow;

    % Guardar la animación como gif:
    frame = getframe(gcf); % Captura la animación
    img = frame2im(frame); % Conversión de la imagen
    [A, map] = rgb2ind(img, 256); % Convierte a índice GIF

    if nivel == 1
        imwrite(A, map, filename, "gif", ...
                "LoopCount", Inf, "DelayTime", 0.05);
    else
        imwrite(A, map, filename, "gif", ...
                "WriteMode", "append", "DelayTime", 0.05);
    end
    
    % Limpieza:
    delete(findall(gca, 'Type', 'scatter'));
    delete(findall(gca, 'Type', 'quiver'));

end

5.5 Masa total aire dentro de la torre:

Teniendo en cuenta que la densidad del aire varía según la función: [math]\varrho(T) = \varrho(aire)\frac{Tº+273}{T+273}[/math] (temperatura ambiente [math]Tº=15\mathcal{ºC}[/math]/densidad del aire es [math]\varrho(aire) = 1225\frac{kg}{m^3}[/math])

Nos calculamos la temperatura dentro de la torre analíticamente resolviendo esta integral triple mediante este código: [math]Masa = \oint_{V}^{}\varrho(T(\rho,z))dV[/math]

[math]Masa = \int_{0}^{2π}\int_{30}^{55}\int_{0}^{150}1225\frac{15+273}{338-38(\frac{z}{150})^{1,8}-8(1-e^{-\frac{\rho^2}{55^2-\rho^2}})}\rho dzd\rho d\theta [/math]

[math] Masa = 971206636.3383 kg [/math]

Cabe destacar que si mantenemos el aire a temperatura constante cómo por ejemplo [math]T=15\mathcal{ºC}[/math], la masa se vería enormemente afectada [math]\varrho(15) = \varrho(aire)\frac{Tº+273}{Tº+273} = 1225 [/math]

Entonces al calcular la masa nos quedaría:

[math]Masa(15) = \int_{0}^{2π}\int_{30}^{55}\int_{0}^{150}1225\rho dzd\rho d\theta [/math]

[math]Masa(15) = 1117590170.588 kg [/math]

Que al compararlo con la masa anterior observamos que hay una diferencia entre ambos de [math] 146383534.2498 kg [/math]

% Masa de aire dentro de la torre 


% Función con densidad variable
f3 = @(z,p,t) 1225 * ((15+273)*p ./ ...
   (338 - 38*(z/150).^1.8 - 8*(1 - exp((p.^2)./-(55^2 - p.^2)))));


% Función con aire uniforme
f2 = @(z,p,t) 1225 * p;


% Límites de integración
t_min = 0;  t_max = 2*pi;
p_min = 33; p_max = 55;
z_min = 0;  z_max = 150;


% Integral triple (orden: t, p, z)
I = integral3(f3, t_min, t_max, p_min, p_max, z_min, z_max);
W = integral3(f2, t_min, t_max, p_min, p_max, z_min, z_max);


% Diferencia entre ambas
X = I - W;




% Resultados 
disp(['La masa total de aire dentro de la torre es ', num2str(I)])
disp(['La masa del aire con temperatura uniforme a 15º es ', num2str(W)])
disp(['La resta de ambos da como resultado ', num2str(X)])

5.6 Potencia térmica evacuada:

Para calcular la potencia térmica evacuada primero necesitamos saber la temperatura media, para eso resolvemos la siguiente integral: [math]T_{\text{media}} = \frac{1}{V} \iiint_V T(\rho,z) \, dV [/math]

[math] Tmedia = \frac{\displaystyle \int_{0}^{2\pi} \int_{33}^{55} \int_{0}^{150} \left[ 65 - 38 \left(\frac{z}{150}\right)^{1.8} - 8 \left(1 - e^{-\frac{\rho^2}{55^2 -\rho^2}}\right)\right] \rho \, dz \, d\rho \, d\theta}{\displaystyle \int_{0}^{2\pi} \int_{33}^{55} \int_{0}^{150} \rho \, dz \, d\rho \, d\theta} [/math]

[math] Tmedia = 58.4149ºC [/math]

Una vez obtenida, calculamos la potencia que se evacua mediante la resolución de la siguiente fórmula: [math] P = \rho_{\text{aire}} \, c_p \, Q \, (T_{\text{media}} - T_{\text{ambiente}}) [/math]


[math] P = 1225\times1005\times1500\times(58'4149 - 15) [/math]

[math]P = 8,017375314 \times (10^{10}) \frac{J}{s} [/math]

% Función temperatura (dentro de la integral) 
f3 = @(z,p,t) p .* (65 - 38*(z/150).^1.8 - 8*(1 - exp((p.^2)./-(55^2 - p.^2))));

% Función con aire uniforme
f2 = @(z,p,t) p;


% Límites de integración
t_min = 0;  t_max = 2*pi;
p_min = 33; p_max = 55;
z_min = 0;  z_max = 150;


% Integral triple (orden: t, p, z)
I = integral3(f3, t_min, t_max, p_min, p_max, z_min, z_max);
W = integral3(f2, t_min, t_max, p_min, p_max, z_min, z_max);


% Diferencia entre ambas
X = 1/W * I ;




% Resultados
disp(['Volumen es ', num2str(W)])
disp(['Resultado de la integral temperatura ', num2str(I)])
disp(['La temperatura media dentro de la torre ', num2str(X)])

6 Ventajas frente a una forma cilíndrica

Central nuclear de Cofrentes (España)

A nivel estructural, la forma hiperboloide es una superficie doblemente reglada, es decir, cada punto puede obtenerse como la intersección de dos líneas rectas. Esto permite que la estructura sea mucho más resistente (especialmente frente a cargas de viento y esfuerzos de compresión ) y al mismo tiempo más rígida y barata que un cilindro. Su geometría hace posible construirla utilizando únicamente vigas rectas, lo que simplifica el encofrado, reduce tiempos de ejecución y disminuye la cantidad de material necesario. Además, esta forma logra una mayor resistencia con un espesor menor.

Desde el punto de vista aerodinámico, el hiperboloide también presenta ventajas frente a una torre cilíndrica. La “cintura” característica de su forma genera una constricción que acelera el paso del flujo de aire por efecto Venturi: al estrecharse el paso, aumenta la velocidad del aire y disminuye la presión. A esto se suma el efecto chimenea, que mejora la circulación natural del aire en el interior de la torre.

Como resultado de ambos fenómenos, la forma hiperboloide proporciona un tiro de aire mayor que el de una torre cilíndrica, incrementando así la eficiencia general de la planta.

7 Otras estructuras hiperbólicas

Las principales ventajas de usar una estructura hiperboloide en estructuras y edificación son su gran resistencia estructural, la facilidad para construirla con vigas rectas (al tratarse de una superficie doblemente reglada) y su singularidad estética. Esta geometría permite crear cubiertas amplias con pocos puntos de apoyo, optimizando el espacio y generando un impacto visual dinámico y moderno.

Varios ejemplos de estructuras hiperbólicas son los siguientes:

1. Torre del Puerto de Kobe (Japón): Se trata de una torre construida con un enrejado de tubería que forma una superficie hiperbólica.

Torre del Puerto de Kobe (Japón)

2. Cubierta del Hipódromo de la Zarzuela (Madrid, España): Esta cubierta se compone de láminas de hormigón armado con forma de sectores de hiperbóloides que consiguen una estructura ligera y resistente.

Hipódromo de la Zarzuela (Madrid, España)

.

3. Torre de Control Aeropuerto El Prat (Barcelona, España): La torre se encuentra recubierta por una estructura metálica en forma de superficie hiperbólica.

Torre de Control El Prat (Barcelona, España)

8 Póster

Enlace al pdf del póster: Medio:Torres de enfriamiento hiperbólicas (G35).pdf