Diferencia entre revisiones de «Torres de enfriamiento hiperbólicas (Grupo 3)»

Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/mat/public_html/w/includes/diff/DairikiDiff.php on line 434

Revisión actual del 13:15 26 nov 2025

Trabajo realizado por estudiantes
Título	Torres de enfriamiento hiperbólicas Grupo 3
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Asier Fernández Torrijos Jorge Mayoral Bote Samuel Portela Trujillo Alejandro Santos Álvarez
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 estructuras muy importantes en las instalaciones industriales y su característico diseño optimiza la transferencia de calor. Su sección hiperbólica es ideal para reducir el uso de materiales de construcción y garantizar su resistencia frente a fuertes vientos.

La geometría hiperbólica se ha estandarizado como la mejor solución para cumplir con su cometido. Su base ancha proporciona un área extensa para favorecer la evaporación del agua, mientras que el estrechamiento central incrementa la velocidad de disipación del aire. Además, la forma hiperbólica contribuye significativamente a la estabilidad estructural.

En este artículo se analizará un modelo estandarizado de torre de enfriamiento, definido por una altura total ([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 [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math]. La superficie de la torre se describe matemáticamente mediante la siguiente ecuación en coordenadas cartesianas:

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

Sean [math]a, c, z_0\gt0[/math] unos valores a determinar

Sean en nuestro modelo [math]Rmáx=50m[/math] ; [math]Rmín=20m[/math] ; [math]H=120m[/math]

2 Torre de enfriamiento y su geometría hiperbólica

Como se ha mencionado, la torre de enfriamiento sigue una geometría hiperbólica. Es importante aclarar que la torre es simétrica respecto del plano de corte horizontal en [math]z_0[/math], y el dominio abarca desde [math]z=0[/math] hasta [math]z=120[/math]. El [math]z_0[/math] marca la altura a la que se encuentra el radio mínimo ([math]Rmín[/math]) y se encuentra a [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math], es decir, a [math]80m[/math] del suelo. Es por ello que el radio de la base no coincide con el radio de la cúspide de la torre.

2.1 Parámetros y ecuación del hiperboloide

2.1.1 Significado de los valores

El valor de [math]a[/math] está asociado con las dimensiones de la torre en el plano [math]z[/math], que es el plano horizontal. Dicho valor controla el ancho horizontal de la figura y fija el radio mínimo ([math]Rmín[/math]).

El valor [math]c[/math] determina cómo cambia la forma del hiperboloide a lo largo del eje [math]z[/math]. Regula la curvatura de la hipérbole en la dirección vertical.

El valor [math]z_0[/math] representa la altura a la que se encuentra el radio mínimo ([math]Rmín[/math]) del hiperboloide respecto de la base. En nuestra torre de enfriamiento, esta altura se encuentra a [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math].

2.1.2 Paso de coordenadas cartesianas a cilíndricas

Para encontrar los valores de [math]a, c, z_0\gt0[/math] primero pasaremos nuestra ecuación del hiperboloide de coordenadas cartesianas [math](x, y, z)[/math] a coordenadas cilíndricas [math](r, θ, z)[/math]. Para ello, se tomará la ecuación en cartesianas y se realizarán las siguientes relaciones para transformarla en coordenadas cilíndricas:

[math]x=ρ·cos(θ) ; y=ρ·sen(θ) ; z=z[/math]

[math]sen(θ)^2+cos(θ)^2=1[/math]

[math]x^2+y^2=ρ^2·cos(θ)^2+ρ^2·sen(θ)^2=ρ^2·(sen(θ)^2+cos(θ)^2)=ρ^2[/math]

Por lo tanto la ecuación en coordenadas cilíndricas es:

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

2.1.3 Cálculo de valores

Para el cálculo de valores se deben resolver las siguientes dos ecuaciones con sus dos incógnitas [math]a[/math] y [math]c[/math] :

Primera ecuación: [math]z=0[/math], en la que sabemos que el radio (ρ) es igual al radio máximo ([math]Rmáx=50m[/math]):

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

Segunda ecuación: [math]z=80[/math], en la que sabemos que el radio (ρ) es igual al radio mínimo ([math]Rmín=20m[/math]):

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

De la segunda ecuación se obtiene directamente el valor de [math]a[/math]. Después se obtiene el valor de [math]c[/math] sustituyendo [math]a[/math] en la primera ecuación:

[math]1=\dfrac{20^2}{a^2} ; a=20[/math]

[math]1=\dfrac{50^2}{a^2}-\dfrac{(-80)^2}{c^2}=\dfrac{50^2}{20^2}-\dfrac{(-80)^2}{c^2} ; c≈34,915[/math]

2.1.4 Conclusión

Los valores finales que describen nuestra torre de enfriamiento modelo son: [math]a=20, c=34,915[/math] y [math]z_0=80[/math]. La ecuación final en coordenadas cilíndricas que modeliza nuestra torre es:

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

2.2 Descripción como superficie reglada

Estructura representativa de superficie hiperbólica reglada

Una superficie reglada es una superficie generada por la estela de una recta, la generatriz, que se desplaza por una o varias curvas, las directrices.

2.2.1 Directrices

La primera directriz es un círculo inferior en el plano [math]z=0[/math] con radio igual a [math]Rmáx=50m[/math]

[math]r_1(t)=(Rmáx·cos(t), Rmáx·sen(t), z), tε[0,2π][/math]

La segunda directriz es un círculo superior en el plano [math]z=120[/math] con radio igual a [math]R≈30,414m[/math]

[math]r_2(t)=(R·cos(t), R·sen(t), z), tε[0,2π][/math]

La tercera directriz es un círculo intermedio en el plano [math]z=80[/math] con radio igual a [math]Rmín=20m[/math]

[math]r_3(t)=(Rmín·cos(t), Rmín·sen(t), z), tε[0,2π][/math]

2.2.2 Generatriz

Esta es la generatriz:

[math]r(t,s)=(1-s)·r_1(t)+s·r_3(t)[/math]

El parámetro de [math]s[/math] varía entre [math]0[/math] y [math]1[/math]

Para [math]s=0[/math] la generatriz está en el punto [math]r_1(t)[/math] que corresponde a [math]z=0[/math]

Para [math]s=1[/math] la generatriz está en el punto [math]r_3(t)[/math] que corresponde a [math]z=120[/math]

2.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=120[/math] para una mejor visualización:

Visualización torre de enfriamiento

% Parámetros de la torre:
a = 20;
c = 34.915;
z_0 = 80;
R_max = 50;
R_min = 20;
z_min = 0;
z_max = 120;

% 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;

3 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]15m/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].

3.1 Mapa de presión por viento

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}·(15·(\dfrac{z}{10})^{0,14})^2≈70,849·z^{0,28}[/math]

Aquí adjuntamos una visualización de como este viento lateral aplica presión en uno de los laterales de la superficie:

Torre de enfriamiento con presión del viento lateral

% Parámetros geométricos
a = 20;
c = 34.915;
z_0 = 80;
R_min = 20;
R_max = 50;
z_max = 120;

% Discretización
theta = linspace(0, pi, 100);
z = linspace(0, z_max, 100);

% Generación de la superficie del hiperboloide
[Theta, Z] = meshgrid(theta, z);
Ro = sqrt(a^2 * (1 + ((Z - z_0).^2 / c^2)));

% Restricción a los radios mínimos y máximos
Ro(Z < z_0) = min(R_max, Ro(Z < z_0));
Ro(Z > z_0) = max(R_min, Ro(Z > z_0));

% Conversión a coordenadas cartesianas
X = Ro .* cos(Theta);
Y = Ro .* sin(Theta);

% Cálculo de la presión del viento
P = 70.849 * Z.^0.28;

% Gráfico de la torre
figure;
surf(X, Y, Z, P, 'EdgeColor', 'none');
colormap('parula');
colorbar;
caxis([min(P(:)), max(P(:))]);
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Torre de enfriamiento con presión del viento lateral');
view(3);
axis equal;

3.2 Campo de fuerza en la superficie expuesta

La fuerza que ejerce el viento sobre la torre actúan en la dirección normal a la superficie, es decir, que el campo vectorial se representa de la siguiente forma:

[math]\vec{F}(x, y, z)=-P(z)·\vec{n}[/math]

A continuación se adjunta una imagen con la visualización de la fuerza del viento sobre la superficie de la torre, y su correspondiente código en MATLAB:

Campo vectorial de fuerza sobre la mitad expuesta de la torre

% Parámetros del hiperboloide
a = 20; c = 34.91; z0 = 80; H = 120;zref=10;

% Parámetros del viento
V0 = 15;
alpha = 0.14;
rho = 1.225;

% Malla de puntos para la representación  de la mitad de la torre
theta = linspace(0, pi, 30);
z = linspace(0, H, 30);
[Theta, Z] = meshgrid(theta, z);

% Coordenadas cilíndricas de la torre
R = a * sqrt(1 + ((Z - z0).^2 / c^2));
X = R .* cos(Theta);
Y = R .* sin(Theta);

%Funciones
V_z = V0 * (Z / zref).^alpha;
P_z = (1/2) * rho * V_z.^2;

% Derivadas parciales de la parametrización
dX_dtheta = -R .* sin(Theta);
dY_dtheta = R .* cos(Theta);
dZ_dtheta = zeros(size(Z));
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)); % Derivada de Z respecto a sí misma

% Producto vectorial para obtener el vector normal
NX = dY_dtheta .* dZ_dz - dZ_dtheta .* dY_dz;
NY = dZ_dtheta .* dX_dz - dX_dtheta .* dZ_dz;
NZ = dX_dtheta .* dY_dz - dY_dtheta .* dX_dz;

% Normalización de los vectores normales
Norm = sqrt(NX.^2 + NY.^2 + NZ.^2);
NX = NX ./ Norm;  NY = NY ./ Norm;  NZ = NZ ./ Norm;

% Fuerza generada por la presión componenete por componente
FX = -P_z .* NX;  FY = -P_z .* NY;  FZ = -P_z .* NZ;

% Graficar el campo vectorial de la fuerza
quiver3(X, Y, Z, FX, FY, FZ, 1.5, 'k');
colorbar;
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Campo vectorial de fuerza sobre la mitad expuesta de la torre');
axis equal; view(3);

4 Campo de temperatura

Dentro de la torre de enfriamiento tenemos un campo de temperaturas modelizado mediante la siguiente expresión:

[math]T(r,z)=T_{base}-ΔT_z·(\dfrac{z}{H})^n-ΔT_r·(1-e^{\dfrac{-(r^2)}{Rmáx^2-r^2}})[/math]

Sea [math]T_{base}[/math] la temperatura del punto central de la base, es decir, cuando [math]r=z=0[/math].

Sea [math]ΔT_z=T_{base}-T_{top}[/math], siendo [math]T_{top}[/math] la temperatura en la parte superior.

Sea [math]ΔT_r=T_{base}-T_{Rmáx}[/math].

Sea [math]n\gt1[/math] un exponente que modela la convección.

En nuestro modelo utilizaremos los siguientes parámetros: [math]T_{base}=70ºC[/math]; [math]T_{top}=25ºC[/math]; [math]n=1,5[/math]; [math]ΔT_r=5ºC[/math].

4.1 Representación del campo de temperatura a través de diferentes planos de corte

Primeramente se trabaja la ecuación inicial introduciendo los valores que podemos calcular:

[math]ΔT_z=T_{base}-T_{top}=70-25=45ºC[/math]

[math]H=120m[/math] y [math]Rmáx=50m[/math]

[math]T(r,z)=70-45·(\dfrac{z}{120})^{1,5}-5·(1-e^{\dfrac{-(r^2)}{50^2-r^2}})[/math]

A continuación se muestran dos representaciones del campo de temperatura en nuestra torre de enfriamiento: La primera representación que muestra el campo de temperatura respecto de un corte de un plano vertical que pasa por el eje; y la segunda representación que consisten en una animación de varios cortes horizontales que ascienden desde la base hacia la cima de nuestra torre:

Campo de Temperaturas respecto de un plano vertical

% Parámetros del hiperboloide truncado
z_max = 120;
r_base = 50;
r_min = 20;
r_top = 34;
z_neck = 80;

% Crear coordenadas cilíndricas para el hiperboloide
z = linspace(0, z_max, 200);
theta = linspace(0, 2*pi, 200);
[Theta, Z_cyl] = meshgrid(theta, z);

% Calcular el radio en cada altura
R_hyper = sqrt(r_min^2 + ((r_base^2 - r_min^2) / z_neck^2) * (z - z_neck).^2);

% Expandir a coordenadas cartesianas
R_hyper_grid = repmat(R_hyper', 1, length(theta));
X_cyl = R_hyper_grid .* cos(Theta);
Y_cyl = R_hyper_grid .* sin(Theta);

% Crear el plano de corte
[X_plane, Z_plane] = meshgrid(linspace(-r_base, r_base, 200), z);
R_plane = abs(X_plane); % Radio en el plano
T_plane = 70 - 45 * (Z_plane / z_max).^1.5 - 5 * (1 - exp(-R_plane.^2 ./ (50^2 - R_plane.^2)));

% Reemplazar valores no definidos en el plano
T_plane(isnan(T_plane) | isinf(T_plane)) = NaN;

% Gráfica 3D
figure;
hold on;

% Dibujar el hiperboloide truncado (torre de enfriamiento)
surf(X_cyl, Y_cyl, Z_cyl, 'FaceColor', [0.8, 0.8, 0.8], 'EdgeColor', 'none', 'FaceAlpha', 0.3);

% Dibujar el plano con el mapa de colores del campo de temperaturas
surf(X_plane, zeros(size(X_plane)), Z_plane, T_plane, 'EdgeColor', 'none');

% Configuración de la gráfica
colorbar;
colormap(jet);
caxis([min(T_plane(:)), max(T_plane(:))]); % Ajustar escala de colores al rango de T
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Altura z (m)');
title('Campo de Temperaturas respecto de un plano vertical');
view(3);
axis equal;
grid on;
hold off;

Animación de planos de corte paralelos

Figura completa con los planos solapados

% Parámetros
a = 20;
c = 34.17;
z0 = 80;
H = 120;
Rmax = 50;
Rmin = 10;
Tbase = 70;
Ttop = 25;
DeltaTz = Tbase - Ttop;
DeltaTr = 5;
n = 1.5;

% Malla para el hiperboloide
theta = linspace(0, 2*pi, 50);
z = linspace(0, H, 50);

% Iniciar figura
figure;
for i = 1:length(z)
% Altura de la capa actual
zi = z(i);

% Radios en esta altura
r_ext = sqrt((1 + (zi - z0)^2 / c^2) * a^2);
r_int = max(Rmin, r_ext - 5);

% Coordenadas de la superficie externa
x_ext = r_ext * cos(theta);
y_ext = r_ext * sin(theta);
z_ext = zi * ones(size(x_ext));

% Coordenadas de la superficie interna
x_int = r_int * cos(theta);
y_int = r_int * sin(theta);
z_int = zi * ones(size(x_int));

% Temperatura en esta capa
T_layer = Tbase - DeltaTz * (zi / H)^n - DeltaTr * (1 - exp(-r_ext^2 / (Rmax^2 - r_ext^2)));

% Representar la capa externa
fill3(x_ext, y_ext, z_ext, T_layer, 'EdgeColor', 'none');
hold on;

% Representar la capa interna
fill3(x_int, y_int, z_int, T_layer, 'EdgeColor', 'none');

% Conexión entre superficies (laterales)
for j = 1:length(theta) - 1
patch([x_ext(j) x_ext(j+1) x_int(j+1) x_int(j)], ...
[y_ext(j) y_ext(j+1) y_int(j+1) y_int(j)], ...
[z_ext(j) z_ext(j+1) z_int(j+1) z_int(j)], ...
T_layer, 'EdgeColor', 'none');
end

% Configuración gráfica
colormap('jet');
colorbar;
caxis([Ttop, Tbase]);
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('Construcción el hiperboloide hueco con el campo de temperatura');
xlim([-Rmax, Rmax]);
ylim([-Rmax, Rmax]);
zlim([0, H]);
view(3);
pause(0.1);
axis equal;
end
hold off;

% Torre de enfriamiento en 3D
surf(X_Tower, Y_Tower, Z_Tower, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', [0.5, 0.5, 0.5]);
hold on;

% Campo vectorial en el plano (r, z)
[X_Plane, Y_Plane] = meshgrid(r, zeros(size(r)));% Plano contenido en el eje (x=radio)
quiver3(X_Plane, Y_Plane, Z, Grad_R, zeros(size(Grad_R)), Grad_Z, 'r', 'AutoScale', 'on', 'LineWidth', 0.5);

% Configuración de la gráfica
xlabel('X (radio)');
ylabel('Y (simetría)');
zlabel('Z (altura)');
title('Campo vectorial y geometría de la torre de enfriamiento');
grid on;
view(3);

4.2 Representación del gradiente a través de diferentes planos de corte

El gradiente de temperatura en un campo [math]T(r,z)[/math] que viene dado por el vector gradiente de la función de temperatura:

[math]GradT(r,z)=(\dfrac{dT}{dr},\dfrac{dT}{dz})[/math]

El gradiente de temperatura tiene dos componentes:

Componente radial:[math]\dfrac{dT}{dr}[/math]

Componente vertical:[math]\dfrac{dT}{dz}[/math]

Respecto la ecuación del campo de temperatura anterior, Calculamos las derivadas respecto a la [math]r[/math] (gradiente radial) y respecto a la [math]z[/math] (gradiente vertical):

[math]\dfrac{dT}{dr}=-ΔT_r·\dfrac{-2r·Rmáx^2}{(Rmáx^2-r^2)^2}·e^{\dfrac{-r^2}{Rmáx^2-r^2}}[/math]

[math]\dfrac{dT}{dz}=-ΔT_z·n·\dfrac{z^{n-1}}{H^n}[/math]

A la hora de representar los planos los dividimos en un plano vertical (que corta la torre por el eje de simetría) y en un plano horizontal (a la altura del radio mínimo).

El plano vertical es el que pasa por [math]r=0[/math], lo que significa que el gradiente radial será cero en todo momento (ya que no hay cambio de temperatura respecto al radio en el centro de la torre). Por tanto, el gradiente será únicamente en la dirección vertical.

El plano horizontal estará ubicado a una altura [math]z[/math], donde se tomará el radio mínimo (Rmín).

El gradiente de temperatura en este caso estará dado por la derivada radial, y la derivada vertical dependerá de la altura [math]z[/math].

A continuación se muestra una representación del gradiente respecto de un plano vertical que corta por el eje:

Plano vertical que corta por el eje

% Parámetros del hiperboloide y el campo de temperatura
a = 20;
c = 34.915;
z0 = 80;
H = 120;
r_max = 50;

% Definición de la malla (ejes r y z)
r = linspace(-r_max, r_max, 100);
z = linspace(0, H, 100);
[R, Z] = meshgrid(r, z);

% Ecuación del hiperboloide (en coordenadas 2D r-z)
hyperboloid_mask = (R.^2 / a^2 - (Z - z0).^2 / c^2) <= 1;

% Campo de temperatura T(r, z)
T = 70 - 45 * (Z / H).^1.5 - 5 * (1 - exp(-R.^2 ./ (50^2 - R.^2)));

% Gradiente completo del campo de temperatura
[Tr, Tz] = gradient(T, r, z); % Gradiente en r y z

% Enmascarar valores fuera del hiperboloide
Tr(~hyperboloid_mask) = NaN;
Tz(~hyperboloid_mask) = NaN;

% Escalar el gradiente para aumentar la visibilidad de las flechas (ajuste de tamaño)
Tz_scaled = Tz / max(abs(Tz(:))) * 40;
Tr_scaled = Tr / max(abs(Tr(:))) * 40;

% Tomar una submalla de las coordenadas para reducir el número de flechas
skip = 3;
R_sub = R(1:skip:end, 1:skip:end);
Z_sub = Z(1:skip:end, 1:skip:end);
Tr_sub = Tr(1:skip:end, 1:skip:end);
Tz_sub = Tz(1:skip:end, 1:skip:end);

% Representación del gradiente completo (componentes en r y z)
figure;
quiver(R_sub, Z_sub, Tr_sub, Tz_sub, 'r-', 'LineWidth', 1.5, 'MaxHeadSize', 5);
hold on;

% Graficar el contorno del hiperboloide a ambos lados
theta = linspace(0, 2*pi, 100);
z_contour = linspace(0, H, 100);
r_contour = a * sqrt((1 + (z_contour - z0).^2 / c^2));

% Graficar el contorno positivo (lado derecho)
plot(r_contour, z_contour, 'k', 'LineWidth', 2);

% Graficar el contorno negativo (lado izquierdo)
plot(-r_contour, z_contour, 'k', 'LineWidth', 2);

% Configuración de la gráfica
xlabel('x (m)');
ylabel('z (m)');
title('Gradiente de temperatura respecto del plano vertical que pasa por el eje');
axis equal;
xlim([-r_max r_max]);
ylim([0 H]);
grid on;

A continuación se muestra una representación del gradiente respecto de un plano horizontal que corta en [math]z_0[/math] a [math]80m[/math] de altura, donde se encuentra el [math]Rmín[/math]:

Plano horizontal que corta en [math]z_0[/math]

% Parámetros de la torre
Rmin = 20;
Rmax = 50;
H = 120;
z_horizontal = (2/3) * H;

% Parámetros para la malla dentro de un cuadrado [-60, 60] en X y Y, con mayor resolución
res_x = 50;
res_y = 50;
x = linspace(-60, 60, res_x); % Rango de X de -60 a 60
y = linspace(-60, 60, res_y); % Rango de Y de -60 a 60

% Generación de malla
[X, Y] = meshgrid(x, y);

% Aplicar filtro para obtener puntos dentro del círculo de radio 20 metros
in_circle = X.^2 + Y.^2 <= Rmin^2;

% Filtrar las coordenadas X y Y
X_filtered = X(in_circle);
Y_filtered = Y(in_circle);

% Direcciones del gradiente (flechas hacia abajo)
Tz_grad = -0.8 * ones(size(X_filtered)); % Componente Z de las flechas (fijas, hacia abajo)

% Gráfica
figure;
hold on;

% Dibujar flechas del gradiente en el plano horizontal (z = z_horizontal)
quiver3(X_filtered, Y_filtered, z_horizontal * ones(size(X_filtered)), ...
   zeros(size(X_filtered)), zeros(size(Y_filtered)), Tz_grad, 0.5, ...
   'Color', [1, 0, 0], 'LineWidth', 1, 'MaxHeadSize', 2);

% Contorno del círculo con radio 20 metros
theta = linspace(0, 2 * pi, 100);
x_circle = Rmin * cos(theta);
y_circle = Rmin * sin(theta);
z_circle = z_horizontal * ones(size(theta));

% Dibujar el contorno en el plano horizontal
plot3(x_circle, y_circle, z_circle, 'k-', 'LineWidth', 2);

% Configuración de la gráfica
title('Gradiente de temperatura en el plano horizonatal donde el radio es mínimo');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Altura Z (m)');
axis([-60 60 -60 60 0 150]);
grid on;

% Ajustar la escala de los ejes para exagerar el eje Z
daspect([1 1 0.2]);% Relación entre los ejes X:Y:Z

% Vista 3D inclinada para mejor percepción
view(75, 30);

% Animación para exagerar el eje Z dinámicamente
for scale = 1:0.1:2
   daspect([1 1 scale]);
   pause(0.1);
end

4.3 Representación de superficies isotérmicas

Para representar las superficies isotermas hacemos uso de la siguiente ecuación del campo de temperatura:

[math]T(r,z)=70-45·(\dfrac{z}{120})^{1,5}-5·(1-e^{\dfrac{-(r^2)}{50^2-r^2}})[/math]

Para cada valor de temperatura obtenemos una superficie que depende de la altura y del radio. Con esta ecuación las superficies que salen son relativamente parecidas a media esfera para cada valor de temperatura.

Para una mayor facilidad de entendimiento, hemos creado un código que representa una vista lateral de la figura con la representación de 5 superficies isotermas que se verán como líneas y otro código que representa una animación en 3D de las mismas 5 superficies isotermas.

Superficies isotermas dentro del hiperboloide

% Parámetros del hiperboloide y el campo de temperatura
a = 20;
c = 34.915;
z0 = 80;
H = 120;
r_max = 50;

% Definición de la malla (ejes r y z)
r = linspace(-r_max, r_max, 100);
z = linspace(0, H, 100);
[R, Z] = meshgrid(r, z);

% Campo de temperatura T(r, z)
T = 70 - 45 * (Z / H).^1.5 - 5 * (1 - exp(-R.^2 ./ (50^2 - R.^2)));

% Definir 5 valores isotérmicos (valores de temperatura)
T_values = linspace(min(T(:)), max(T(:)), 7);
 
% Crear la figura
figure;
hold on;

% Representar las superficies isotermas
for i = 1:length(T_values)-1
   contour(R, Z, T, [T_values(i) T_values(i)], 'LineWidth', 2); 
end

% Ecuación del hiperboloide (en coordenadas 2D r-z)
hyperboloid_mask = (R.^2 / a^2 - (Z - z0).^2 / c^2) <= 1;

% Graficar el contorno del hiperboloide a ambos lados
theta = linspace(0, 2*pi, 100);
z_contour = linspace(0, H, 100);  
r_contour = a * sqrt((1 + (z_contour - z0).^2 / c^2));
  
% Graficar el contorno positivo (lado derecho)
plot(r_contour, z_contour, 'k-', 'LineWidth', 2);
  
% Graficar el contorno negativo (lado izquierdo)
plot(-r_contour, z_contour, 'k-', 'LineWidth', 2); 
 
% Configuración de la gráfica
xlabel('r (Radio)');
ylabel('z (Altura)');
title('Superficies isotermas dentro del hiperboloide');
axis equal;
xlim([-r_max r_max]);   
ylim([0 H]);            
grid on;

% Ajustar la apariencia del fondo y las líneas de la cuadrícula para mejor visualización
set(gca, 'Color', [1 1 1]);

Animación de superficies isotermas

% Parámetros del hiperboloide y el campo de temperatura
a = 20;
c = 34.915;
z0 = 80;
H = 120;
r_max = 50;

% Definición de la malla cilíndrica
theta = linspace(0, 2*pi, 100);
z = linspace(0, H, 100);
r = linspace(0, r_max, 50);
[Theta, Z, R] = meshgrid(theta, z, r);

% Conversiones a coordenadas cartesianas para el campo
X = R .* cos(Theta);
Y = R .* sin(Theta);

% Campo de temperatura T(r, z)
T = 70 - 45 * (Z / H).^1.5 - 5 * (1 - exp(-R.^2 ./ (50^2 - R.^2)));

% Defino valores isotérmicos (valores de temperatura)
T_values = linspace(min(T(:)), max(T(:)), 7);
T_values = [max(T(:)) + 10, T_values];

% Ecuación del hiperboloide
[Theta_h, Z_h] = meshgrid(theta, z);
R_h = a * sqrt(1 + ((Z_h - z0).^2) / c^2);
X_h = R_h .* cos(Theta_h);
Y_h = R_h .* sin(Theta_h);

% Creao figura
figure;
hold on;

% Graficar el hiperboloide en 3D
surf(X_h, Y_h, Z_h,'FaceAlpha', 0.4, 'EdgeColor', 'none', 'FaceColor', [0.5, 0.5, 0.5]);

% Configuración de la gráfica
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Animación de superficies isotermas');
grid on;
axis equal;
view(3);

% Mejorar visualización
light;
lighting phong;
camlight('headlight');

% Crear animación
cmap = jet(length(T_values));
for i = 1:length(T_values)
   T_iso = T_values(i);
   iso_surface = isosurface(X, Y, Z, T, T_iso);
   p = patch(iso_surface);
   set(p, 'FaceColor', cmap(i, :), 'EdgeColor', 'none', 'FaceAlpha', 0.7);
   drawnow;
   pause(0.5);
end

5 Comparativa de eficiencia en otros diseños

Cuando implementamos la condición de que el [math]Rmín=Rmáx=50[/math], calculamos las nuevas constantes [math]a[/math] y [math]c[/math], que calculándolas de la misma forma que anteriormente, el valor de [math]a[/math] no cambia pues no depende del [math]Rmin[/math] y el valor de [math]c[/math] sin embargo tiende a infinito, lo que nos lleva a la conclusión de que se trata de un cilindro.

Se ha creado en MATLAB dos programa que nos arroja la fuerza por unidad de superficie que tendría que soportar cada uno de los modelos. El primer código será el encargado de calcular la fuerza por unidad de superficie que tendría que soportar el cilindro. El segundo código es el encargado de calcular la fuerza por unidad de superficie que soportaría una torre hiperbólico.

Geometría cilíndrica: Silo de grano

% Parámetros dados para el cilindro
Rmax = 50;
H = 120;

% Parámetros para la velocidad del viento
V0 = 10;
zref = 10;
alpha = 0.14;

% Densidad del aire (kg/m^3) a nivel del mar
rho = 1.225;

% Generar puntos para la superficie del cilindro
z = linspace(0, H, 100);
theta = linspace(0, pi, 100);
[Z, Theta] = meshgrid(z, theta);

% Cálculo del radio (en este caso es constante para todo z)
r = Rmax * ones(size(Z));
% Calcular la velocidad del viento a cada altura z
V = V0 * (Z / zref).^alpha;

% Calcular la presión del viento a cada altura z
P = 0.5 * rho * V.^2;

% Cálculo del área diferencial de cada elemento del cilindro (dA = r * dz * dtheta)
dz = z(2) - z(1);
dtheta = theta(2) - theta(1);

% Cálculo de la fuerza total sobre el cilindro (integración de la presión)
F_total = 0;

% Realizamos la integración para calcular la fuerza total
for i = 1:length(z)
   for j = 1:length(theta)
       % Calcular la fuerza diferencial en cada punto
       F_total = F_total + P(i,j) * r(i,j) * dz * dtheta;
   end
end

% Cálculo del área total del cilindro
A_total = pi * Rmax * H;

% Cálculo de la presión total (fuerza total / área total)
P_total = F_total / A_total;

% Mostrar el resultado de la presión total
disp(['Presión total ejercida por el viento sobre el cilindro: ', num2str(P_total), ' Pa (N/m^2)']);

La torre cilíndrica soportaría una fuerza por unidad de superficie igual a [math]97,4457Pa[/math]

Geometría hiperbólica: Torre de enfriamiento

% Parámetros dados
Rmax = 50;
Rmin = 20;
H = 120;
a = 20;
c = 34.91;

% Parámetros para la velocidad del viento
V0 = 10;
zref = 10;
alpha = 0.14;

% Densidad del aire (kg/m^3) a nivel del mar
rho = 1.225;

% Generar puntos para la superficie de la torre
z = linspace(0, H, 100);
theta = linspace(0, pi, 100);
[Z, Theta] = meshgrid(z, theta);

% Cálculo del radio en función de z (es función de la altura z)
r = a * sqrt(1 + (Z - 80).^2 / c^2);

% Calcular la velocidad del viento a cada altura z
V = V0 * (Z / zref).^alpha;

% Calcular la presión del viento a cada altura z
P = 0.5 * rho * V.^2;

% Coordenadas cartesianas de la superficie de la torre (X, Y)
X = r .* cos(Theta);
Y = r .* sin(Theta);

% Cálculo del área diferencial de cada elemento de la torre (dA = r * dz * dtheta)
dz = z(2) - z(1);
dtheta = theta(2) - theta(1);

% Cálculo de la fuerza total sobre la torre (integración de la presión)
F_total = 0;

% Realizamos la integración para calcular la fuerza total
for i = 1:length(z)
   for j = 1:length(theta)
       % Calcular la fuerza diferencial en cada punto
       F_total = F_total + P(i,j) * r(i,j) * dz * dtheta;
   end
end

% Cálculo del área total de la torre
A_total = sum(r(:) * dz * dtheta);

% Cálculo de la presión total (fuerza total / área total)
P_total = F_total / A_total;

% Mostrar el resultado de la presión total
disp(['Presión total ejercida por el viento sobre la torre: ', num2str(P_total), ' Pa (N/m^2)']);

La torre hiperbólica tendría que soportar una fuerza por unidad de superficie igual a [math]89,4393Pa[/math]

Una vez observados los resultados de los cálculos, podemos afirmar que el modelo de torre más eficaz en este estudio sería una torre hiperbólica. En condiciones iguales de viento lateral, la geometría hiperbólica consigue obtener valores de presión menores en su superficie que una geometría cilíndrica.

6 Aplicaciones de estructuras hiperboloides en ingeniería

Torre de Kobe

Las estructuras hiperboloides se emplean en ingeniería debido principalmente a sus características de resistencia y ligereza, lo que las hace muy útiles en aplicaciones arquitectónicas y en el ámbito de la ingeniería civil.

Entre sus características más importantes destacamos que estas estructuras son especialmente valiosas por ser geométricamente eficientes, lo que las hace ideales para soportar grandes cargas con un material mínimo (como acero o concreto) consiguiendo una relación óptima entre peso y resistencia, y consiguiendo así el mayor rendimiento posible.

Otra característica a destacar es la estabilidad estructural ya que debido a sus formas distribuyen las cargas de manera uniforme, aumentando la durabilidad y la seguridad sin olvidar que su diseño aerodinámico también minimiza la resistencia al viento y mejora la estabilidad en zonas sísmicas lo que convierte a estas estructuras en una obra maestra en la ingeniería.

Estas características las hacen ideales para aplicaciones como torres, cubiertas de edificios, puentes y otras construcciones donde se requiere una estructura ligera pero resistente y, a día de hoy, continúan siendo un área de innovación en el diseño de construcciones más sostenibles y económicas a nivel mundial.

Puente con geometría hiperbólica

7 Bibliografía

Torre de enfriamiento

[1]

Hiperboloide

[2]

Superficie reglada

[3]