Prototipo 101
Contenido
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 grandes 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:
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 educació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:
Por lo tanto la ecuación en coordenadas cilíndricas es:
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]):
Segunda ecuación: [math]z=80[/math], en la que sabemos que el radio (ρ) es igual al radio mínimo ([math]Rmín[/math]=20m</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:
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:
2.2 Descripción como superficie 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]
La segunda directriz es un círculo superior en el plano [math]z=120[/math] con radio igual a [math]R≈30,414m[/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]
2.2.2 Generatriz
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:
% 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:
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:
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:
Aquí adjuntamos una visualización de como este viento lateral aplica presión en uno de los laterales de la superficie:
% Parámetros geométricos
a = 20; % Radio base menor
c = 34.915; % Altura característica
z_0 = 80; % Altura del radio mínimo
R_min = 20; % Radio mínimo
R_max = 50; % Radio máximo
z_max = 120; % Altura máxima
% Discretización
theta = linspace(0, pi, 100); % Ángulo para la mitad de la torre
z = linspace(0, z_max, 100); % Altura
% 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'); % Superficie con colores de presión
colormap('parula'); % Mapa de colores
colorbar; % Barra de colores
caxis([min(P(:)), max(P(:))]); % Ajustar rango de colores
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Torre de enfriamiento con presión del viento lateral');
view(3); % Vista 3D
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:
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:
% Parámetros del hiperboloide
a = 20; c = 34.91; z0 = 80; H = 120;zref=10;
% Parámetros del viento
V0 = 15; % Velocidad del viento (m/s)
alpha = 0.14; % Exponente del viento
rho = 1.225; % Densidad del aire (kg/m^3)
% Malla de puntos para la representación de la mitad de la torre
theta = linspace(0, pi, 30); % Solo la mitad de la torre
z = linspace(0, H, 30); % Altura de la torre
[Theta, Z] = meshgrid(theta, z);
% Coordenadas cilíndricas de la torre
R = a * sqrt(1 + ((Z - z0).^2 / c^2)); % Radio en función de Z
X = R .* cos(Theta);
Y = R .* sin(Theta);
%Funciones
V_z = V0 * (Z / zref).^alpha;
P_z = (1/2) * rho * V_z.^2;
% Cálculo de los vectores normales
% Derivadas parciales de la parametrización
dX_dtheta = -R .* sin(Theta);
dY_dtheta = R .* cos(Theta);
dZ_dtheta = zeros(size(Z)); % Z no depende de 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)); % 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'); % Escala ajustada para claridad
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:
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:
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:
% Parámetros del hiperboloide truncado
z_max = 120; % Altura truncada
r_base = 50; % Radio en la base
r_min = 20; % Radio mínimo en el cuello
r_top = 34; % Radio aproximado en la cima truncada
z_neck = 80; % Altura del cuello
% Crear coordenadas cilíndricas para el hiperboloide
z = linspace(0, z_max, 200); % Altura truncada
theta = linspace(0, 2*pi, 200); % Ángulo
[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; % Barra de colores para el campo de temperatura
colormap(jet); % Escala de colores
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); % Vista 3D
axis equal;
grid on;
hold off;
--
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:
El gradiente de temperatura tiene dos componentes:
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):
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].
4.3 Representación de superficies isotérmicas
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.
% Parámetros dados para el cilindro
Rmax = 50; % Radio constante del cilindro (m)
H = 120; % Altura total del cilindro (m)
% Parámetros para la velocidad del viento
V0 = 10; % Velocidad del viento en la altura de referencia (m/s)
zref = 10; % Altura de referencia (m)
alpha = 0.14; % Exponente (valor proporcionado)
% 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); % Puntos de altura (de 0 a H)
theta = linspace(0, pi, 100); % Puntos angulares (de 0 a 2*pi)
[Z, Theta] = meshgrid(z, theta); % Malla de puntos de altura y ángulo
% Cálculo del radio (en este caso es constante para todo z)
r = Rmax * ones(size(Z)); % Radio constante para el cilindro
% Calcular la velocidad del viento a cada altura z
V = V0 * (Z / zref).^alpha; % Velocidad del viento en función de z
% Calcular la presión del viento a cada altura z
P = 0.5 * rho * V.^2; % Presión del viento en función de la velocidad
% Cálculo del área diferencial de cada elemento del cilindro (dA = r * dz * dtheta)
dz = z(2) - z(1); % Diferencia de altura
dtheta = theta(2) - theta(1); % Diferencia angular
% Cálculo de la fuerza total sobre el cilindro (integración de la presión)
F_total = 0; % Inicializamos la fuerza total
% 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; % Área lateral del cilindro
% Cálculo de la presión total (fuerza total / área total)
P_total = F_total / A_total; % Fuerza total por área 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]
% Parámetros dados
Rmax = 50; % Radio máximo en la base (m)
Rmin = 20; % Radio mínimo a 2/3 de la altura (m)
H = 120; % Altura total de la torre (m)
a = 20; % Constante calculada
c = 34.91; % Constante calculada
% Parámetros para la velocidad del viento
V0 = 10; % Velocidad del viento en la altura de referencia (m/s)
zref = 10; % Altura de referencia (m)
alpha = 0.14; % Exponente (valor proporcionado)
% 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); % Puntos de altura (de 0 a H)
theta = linspace(0, pi, 100); % Puntos angulares (de 0 a 2*pi)
[Z, Theta] = meshgrid(z, theta); % Malla de puntos de altura y ángulo
% 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); % Radio en función de z
% Calcular la velocidad del viento a cada altura z
V = V0 * (Z / zref).^alpha; % Velocidad del viento en función de z
% Calcular la presión del viento a cada altura z
P = 0.5 * rho * V.^2; % Presión del viento en función de la velocidad
% 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); % Diferencia de altura
dtheta = theta(2) - theta(1); % Diferencia angular
% Cálculo de la fuerza total sobre la torre (integración de la presión)
F_total = 0; % Inicializamos la fuerza total
% 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); % Suma de las áreas diferenciales
% Cálculo de la presión total (fuerza total / área total)
P_total = F_total / A_total; % Fuerza total por área 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
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.
