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

Revisión del 18:20 5 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 Cálculo 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

En nuestro modelo 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.

Para calcular el campo vectorial [math]\vec{F}(x,y,z) = -P(z)\,\vec{n}[/math] necesitamos calcular primero el vector normal.

Parametrizamos la superficie: \begin{cases} \rho(u,v) = 30 \cdot \sqrt{\,1 + \frac{(v - 100)^2}{65{,}079^2}\,} \\ \theta(u,v) = u \\ z(u,v) = v \\ \end{cases} [math](u,v) \in D = \left[-\tfrac{\pi}{4}, \tfrac{3\pi}{4}\right] \times [0,150][/math]

Tenemos que:

[math] \vec{n}(u, v) = \frac{\vec{r}_u(u, v) \times \vec{r}_v(u, v)}{|\vec{r}_u(u, v) \times \vec{r}_v(u, v)|} [/math]

Ahora calculamos las derivadas parciales necesarias:

\begin{aligned} \frac{\partial \vec{r}}{\partial u} = \frac{\partial \rho}{\partial u} \vec{e_\rho} + \rho \cdot \frac{\partial \theta}{\partial u} \vec{e_\theta} + \frac{\partial z}{\partial u} \vec{e_z} = \rho \cdot \vec{e_\theta} \end{aligned}

\begin{aligned} \frac{\partial \vec{r}}{\partial v} = \frac{\partial \rho}{\partial v} \vec{e_\rho} + \rho \cdot \frac{\partial \theta}{\partial v} \vec{e_\theta} + \frac{\partial z}{\partial v} \vec{e_z} = \rho' \cdot \vec{e_\rho} + \vec{e_z} \end{aligned}

\begin{aligned} \vec{n} = \frac{\partial \vec{r}}{\partial u} \times \frac{\partial \vec{r}}{\partial v} = \begin{vmatrix} \vec{e_\rho} & \vec{e_\theta} & \vec{e_z} \\[6pt] 0 & \rho & 0 \\[6pt] \rho' & 0 & 1 \end{vmatrix} = \rho \cdot \vec{e_\rho} - \rho \rho' \cdot \vec{e_z} \end{aligned} (este es el vector normal sin normalizar)

Normalizando el vector nos queda:

[math] \vec{n} = \frac{\vec{n}}{|\vec{n}|} = \frac{\rho \cdot (1,\,0,\,\rho')}{\rho \cdot \sqrt{1 + (\rho')^2}} = \frac{(1,\,0,\,\rho')}{\sqrt{1 + (\rho')^2}}[/math]

Entonces el campo vectorial de fuerzas, en la nueva parametrización queda:

[math] \boxed{\vec{F}(u,v) = - P(v) \cdot \vec{n}(u,v) = - 102,787 \cdot v^\frac{2}{7} \cdot \frac{(1,\,0,\,\rho'(u,v))}{\sqrt{1 + \rho'(u,v)^2}}}[/math]

La derivada: [math]\rho'(u,v) = \rho'(v) [/math] porque solo depende de v. Esta se calcula paso a paso usando la regla de la cadena.

[math] \begin{aligned} \rho(v) &= 30\sqrt{1+\frac{(v-100)^2}{65,079^2}}\rightarrow A(v)=1+\frac{(v-100)^2}{65,079^2}\rightarrow \rho(v)=30\,A(v)^{1/2} \\[10pt] \frac{d\rho}{dv} &= 30 \cdot \frac{1}{2} A^{-1/2}\cdot A'(v)\rightarrow A'(v)=\frac{2(v-100)}{65,079^2} \\[10pt] \frac{d\rho}{dv} &= 30\cdot \frac{1}{2} A^{-1/2}\cdot \frac{2(v-100)}{65,079^2} \\[10pt] &\boxed{\rho'(v) = \frac{30(v-100)}{65,079^2\sqrt{1+\frac{(v-100)^2}{65,079^2}}}} \end{aligned} [/math]

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.

La fuerza total se calcula integrando [math] \vec{F}(u,v) [/math] sobre la mitad de la superficie del hiperboloide expuesta al viento, es decir, calcular el flujo.

Otra cosa a tener en cuenta es la orientación del vector normal para realizar la integral, en este caso es negativa

[math]F_{TH} = \int_{S} \vec{F} \cdot (-\vec{n}) \,dS = \iint_{D} (-P(v) \cdot \vec{n}) \cdot (-\vec{n}) \cdot |\vec{r_u} \times \vec{r_v}|\,du\,dv = \iint_{D} P(v) \cdot |\vec{r_u} \times \vec{r_v}|\,du\,dv = \iint_{D} P(v) \cdot \rho(u,v) \cdot \sqrt{1 + \rho'(u,v)^2}\,du\,dv[/math]

Como [math]\rho(u,v) [/math] solo depende de v entonces:

[math]F_{TH} = \int_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \,du \int_{0}^{150}P(v) \cdot \rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv = [u]_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \cdot \int_{0}^{150}P(v) \cdot \rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv = \frac{3\pi}{4} -(-\frac{\pi}{4}) \cdot \int_{0}^{150}P(v) \cdot \rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv [/math] [math]\boxed{F_{TH} = \pi \cdot \int_{0}^{150}P(v) \cdot \rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv} [/math]

Para calcular la fuerza por unidad de superficie ([math] F_{usH} = \frac{F_{TH}}{A_{H}}[/math]) se necesita calcular primero el área de la superficie hiperbólica:

[math] A_{H} = \int_{S} \,dS = \iint_{D} |\vec{r_u} \times \vec{r_v}|\,du\,dv = \int_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \,du \int_{0}^{150} \rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv = [u]_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \cdot \int_{0}^{150}\rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv = \frac{3\pi}{4} -(-\frac{\pi}{4}) \cdot \int_{0}^{150}\rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv [/math] [math] \boxed{A_{H} = \pi \cdot \int_{0}^{150}\rho(v) \cdot \sqrt{1 + \rho'(v)^2}\,dv}[/math]

Debido a la complejidad de estas integrales, para resolverlas utilizamos un programa de Matlab que se adjunta a continuación:

% Integral doble en dominio D = [-pi/4, 3pi/4] x [0,150]
% Parámetros
a = 30;
c = 65.079;
z0 = 100;

% Definición de rho(v) y su derivada
rho  = @(v) a .* sqrt(1 + ((v - z0).^2) ./ (c^2));
rho_p = @(v) a*(v - z0) ./ (c.^2 .* sqrt(1 + ((v - z0).^2) ./ c.^2));

% Integrando original
integrando_fuerza = @(v) -102.787 .* v.^(2/7) .* (- rho(v) .* sqrt(1 + (rho_p(v)).^2));
integrando_area = @(v) rho(v) .* sqrt(1 + (rho_p(v)).^2);

% Integral definida en v
IF_v = integral(integrando_fuerza, 0, 150);
IA_v = integral(integrando_area, 0, 150);

% Integral doble (multiplicando por el ancho en u)
IF_total = (3*pi/4 - (-pi/4)) * IF_v;
IA_total = (3*pi/4 - (-pi/4)) * IA_v;

%Fuerza por unidad de superficie
F_usup = IF_total / IA_total;

% Mostrar resultado
fprintf('La fuerza total en el hiperboloide es: %.10f N \n', IF_total);
fprintf('El área del hiperbolide es: %.10f m^2 \n', IA_total);
fprintf('La fuerza por unidad de superfice es: %.10f N/m^2 o Pa \n', F_usup);

Los resultados son:

[math] \boxed{F_{TH} \approx 5787054,768 N \\ A_{H} \approx 18099,151 m^2 \\ F_{usH} = \frac{F_{TH}}{A_{H}} \approx 319,741 Pa} \\ [/math]

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)?

Para calcular la fuerza por unidad de superficie del viento sobre la mitad del cilindro expuesta y el área de este mismo, realizamos una parametrización para simplificar los cálculos:

Parametrizamos la superficie: \begin{cases} \rho(u,v) = 55 \\ \theta(u,v) = u \\ z(u,v) = v \\ \end{cases} [math](u,v) \in D = \left[-\tfrac{\pi}{4}, \tfrac{3\pi}{4}\right] \times [0,150][/math]

Ahora calculamos las derivadas parciales necesarias:

\begin{aligned} \frac{\partial \vec{r}}{\partial u} = \frac{\partial \rho}{\partial u} \vec{e_\rho} + \rho \cdot \frac{\partial \theta}{\partial u} \vec{e_\theta} + \frac{\partial z}{\partial u} \vec{e_z} = \rho \cdot \vec{e_\theta} \end{aligned}

\begin{aligned} \frac{\partial \vec{r}}{\partial v} = \frac{\partial \rho}{\partial v} \vec{e_\rho} + \rho \cdot \frac{\partial \theta}{\partial v} \vec{e_\theta} + \frac{\partial z}{\partial v} \vec{e_z} = \vec{e_z} \end{aligned}

\begin{aligned} \vec{n} = \frac{\partial \vec{r}}{\partial u} \times \frac{\partial \vec{r}}{\partial v} = \begin{vmatrix} \vec{e_\rho} & \vec{e_\theta} & \vec{e_z} \\[6pt] 0 & \rho & 0 \\[6pt] 0 & 0 & 1 \end{vmatrix} = \rho \cdot \vec{e_\rho} \end{aligned} (este es el vector normal sin normalizar)

En este caso necesitamos [math]|\vec{r_u} \times \vec{r_v}| = \rho [/math]

Ahora calculamos la fuerza total:

[math]F_{TC} = \int_{S} \vec{F} \cdot (-\vec{n}) \,dS = \iint_{D} (-P(v) \cdot \vec{n}) \cdot (-\vec{n}) \cdot |\vec{r_u} \times \vec{r_v}|\,du\,dv = \iint_{D} P(v) \cdot |\vec{r_u} \times \vec{r_v}|\,du\,dv = 102,787 \cdot \iint_{D} v^{2/7} \cdot \rho(u,v)\,du\,dv[/math]

Como [math]\rho(u,v)=55 [/math] entonces:

[math]F_{TC} = 102,787 \cdot 55\cdot \int_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \,du \int_{0}^{150}v^{2/7}\,dv = 5653,285 \cdot [u]_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \cdot \left[\frac{7 \cdot \sqrt[7]{v^9}}{9}\right]_{0}^{150} =5653,285 \cdot \left(\frac{3\pi}{4} -(-\frac{\pi}{4})\right) \cdot \frac{7 \cdot \sqrt[7]{150^9}}{9} = 5653,285 \pi \cdot 488,301 [/math]

[math]\boxed{F_{TC} \approx 8672381,345 N} [/math]

Para calcular el área lateral del cilindro afectada por el viento hacemos la integral de superficie:

[math] A_{C} = \int_{S} \,dS = \iint_{D} |\vec{r_u} \times \vec{r_v}|\,du\,dv = \int_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \,du \int_{0}^{150} \rho(u,v)\,dv = [u]_{-\frac{\pi}{4}}^{\frac{3\pi}{4}} \cdot 55 \cdot [v]_{0}^{150} = \frac{3\pi}{4} -(-\frac{\pi}{4}) \cdot 55 \cdot 150 = 55\pi \cdot 150 [/math]

[math] \boxed{A_{C} \approx 25918,139 m^2}[/math]

La fuerza por unidad de superficie en la mitad de la superficie cilíndrica es:

[math] \boxed{F_{usC} = \frac{F_{TC}}{A_{C}} \approx 334,607 Pa } [/math]

Tomando los resultados del caso de la superficie hiperbólica y la cilíndrica tenemos que:

[math]F_{usH} \lt F_{usC} [/math] , lo que indica que la superficie del cilindro soporta menos fuerza lateral que la superficie del hiperboloide.

[math]A_{H} \lt A_{C} [/math] , esto implica que la superficie hiperbólica utiliza menos material que la superficie cilíndrica.

Por tanto, la superficie hiperbólica es estructuralmente más eficiente y estable frente a cargas laterales.

6 Campo de temperatura y transferencia de calor

6.1 Representación del campo escalar de temperatura en MATLAB

Representa el campo escalar 𝑇(𝜌,𝑧) como un mapa de colores en una sección vertical que corte la torre por el eje de simetría y en varias secciones transversales. Añade colorbar¿Dónde está el aire más caliente?

Para la representación del campo escalar sustituimos :

El radio [math]\rho[/math] en la base superior, el cual se calcula sustituyendo en nuestra ecuación del hiperboloide: [math]S(z = 150)[/math] nos da: [math]\rho = 37,222[/math]

La temperatura en la base ([math]z = 0[/math]) y en la tapa superior ([math]z = 150[/math]) con sus respectivos radios en cada altura, que sustituyendo en nuestra ecuación de la temperatura nos da:

6.1.1 Localización del punto mas caliente de la torre

El punto más caliente se encuentra en la base [math]z = 0[/math] centrado en el eje [math]\rho = 0[/math], lo cual se comprueba observando la ecuación de la temperatura [math]T(\rho,z)[/math], siendo esta temperatura:

Esto se debe a que la torre no ha empezado a cumplir con su enfriamiento en ese punto, ya que esta enfría mientras más alto te encuentres y mas alejado del eje.

A continuación se muestra el código de la representación del campo de temperatura a través de secciones transversales y de una sección vertical q pasa por eje de simetría:

% Parámetros geométricos
a = 30;      
c = 65.079;  
z0 = 100;   
H = 150;     
Rmax = 55; 

%Parámetros de temperatura
Tbase = 65;  % Temperatura en la base (°C)
Ttop = 27;   % Temperatura en la parte superior (°C)
DeltaTz = Tbase - Ttop; % Diferencia de temperatura vertical
DeltaTr = 8; % Variación máxima radial (°C)
n = 1.8;     % Exponente de convección

%% 1. Animación de la torre
figure('Name', 'Animación de la Torre');
z = linspace(0, H, 50); % Coordenada vertical

for i = 1:length(z)
    clf; % Limpiar figura
    hold on;
    
    for j = 1:i
        zi = z(j); % Altura actual
        r = sqrt((1 + (zi - z0)^2 / c^2) * a^2); % Radio en esta altura
        
        % Malla radial en el plano
        theta = linspace(0, 2*pi, 100);
        rvec  = linspace(0, r, 50);
        [RR,TT] = meshgrid(rvec, theta);
        X = RR .* cos(TT);
        Y = RR .* sin(TT);
        Z = zi * ones(size(X));
        
        % Temperatura en cada punto de la malla
        Rxy = sqrt(X.^2 + Y.^2);
        T = Tbase - DeltaTz * (zi/H)^n ...
            - DeltaTr * (1 - exp(-(Rxy.^2)./(Rmax^2 - Rxy.^2)));
        
        % Dibujar capa con gradiente de color
        surf(X, Y, Z, T, 'EdgeColor','none');
    end
    
    % Configuración gráfica
    colormap('jet');
    colorbar;
    cb = colorbar; % Crear colorbar
    ylabel(cb, 'Temperatura (°C)', 'FontWeight', 'bold'); 
    clim([Ttop, Tbase]); % Escala de colores
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('z (m)');
    title('Construcción del hiperboloide con T(\rho,z)');
    xlim([-Rmax, Rmax]);
    ylim([-Rmax, Rmax]);
    zlim([0, H]);
    view(3); % Vista en 3D
    axis equal;
    pause(0.1); % Controlar velocidad de animación
end

%% 2. Campo de temperatura en el plano vertical
figure('Name', 'Plano Vertical del Campo de Temperatura');
z = linspace(0, H, 300); % Coordenada vertical
x = linspace(-Rmax - 10, Rmax + 10, 300); % Coordenada horizontal
[X, Z] = meshgrid(x, z); % Crear malla
R = sqrt((1 + ((Z - z0) / c).^2) * a^2); % Radio según el hiperboloide
r = abs(X); % Distancia radial
T = Tbase - DeltaTz * (Z / H).^n - DeltaTr * (1 - exp(-r.^2 ./ (Rmax^2 - r.^2))); % Campo de temperatura
mask = abs(X) <= R; % Máscara para limitar al hiperboloide
T(~mask) = NaN; % Fuera de la silueta, no mostrar
% Dibujar el campo de temperatura
contourf(X, Z, T, 50, 'LineColor', 'none');
colormap('jet');
colorbar;
cb = colorbar; % Crear colorbar
ylabel(cb, 'Temperatura (°C)', 'FontWeight', 'bold'); 
clim([Ttop, Tbase]); % Escala de colores
hold on;
% Dibujar 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 gráfica
xlabel('x (m)');
ylabel('z (m)');
title('Sección vertical de T(\rho,z) en el eje central de la torre');
xlim([-Rmax - 10, Rmax + 10]);
ylim([0, H]);
grid on;
view(2); % Vista 2D
hold off;

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

Siendo el campo escalar: [math] T(\rho, z) = T_{\text{base}} - \Delta T_z \left( \frac{z}{H} \right)^n - \Delta T_r \left( 1 - e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}} \right) [/math]

Calculamos el gradiente utilizando la respectiva formula del gradiente para las coordenadas cilíndricas:

Como el campo [math]T(\rho,z)[/math] no depende de la coordenada [math]\theta[/math], entonces el término asociado a [math]\vec{e}_\theta[/math] desaparece del gradiente. Obteniendo:

Realizando todo:

Gradiente habiendo sustituido los datos:

6.3 Representación del campo vectorial gradiente sobre el hiperboloide

texto alternativo

texto alternativo

%% 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 Representación de superficies isotermas

texto alternativo

%% Apartado 2.6 – Superficies isotermas en la sección r–z de la torre

clear; close all; clc;

%% ---------------- PARÁMETROS DEL PROBLEMA ----------------
Tbase = 65;         
Ttope = 27;
dTz   = 38;  
dTr   = 8;               % variación radial
n     = 1.8;

H = 150;                 % altura total
r_max = 55;              % radio máximo
a = 30;                  % hiperboloide (no afecta T)
c = 65.079;
z0 = 100;

%% ---------------- MALLA EN r–z ----------------
r = linspace(-r_max, r_max, 200);
z = linspace(0, H, 200);
[R, Z] = meshgrid(r, z);

%% ---------------- CAMPO DE TEMPERATURA ----------------
T = Tbase ...
    - dTz * (Z/H).^n ...
    - dTr * (1 - exp(-(R.^2)/(r_max^2)));

%% ---------------- VALORES DE ISOTERMAS ----------------
T_values = linspace(min(T(:)), max(T(:)), 10);

%% ---------------- FIGURA ----------------
figure; hold on;

% Dibujar isotermas
contour(R, Z, T, T_values, 'LineWidth', 1.4);

%% ---------------- CONTORNO DEL HIPERBOLOIDE ----------------
z_contour = linspace(0, H, 300);
r_contour = a * sqrt(1 + (z_contour - z0).^2 / c^2);

plot( r_contour, z_contour, 'k', 'LineWidth', 2);
plot(-r_contour, z_contour, 'k', 'LineWidth', 2);

%% ---------------- AJUSTES DE GRÁFICA ----------------
xlabel('r (m)');
ylabel('z (m)');
title('Isotermas T(r,z) dentro del hiperboloide – Apartado 2.6');
ylim([0 H]);
xlim([-r_max r_max]);
axis equal;
grid on;
set(gca,'Color',[1 1 1]);

6.5 Simulación de un caso en el que entra aire caliente a la base

Consideramos que la densidad del aire depende de la temperatura según:

donde [math]\rho_{\text{aire}} = 1.225 \,\text{kg/m}^3[/math] es la densidad del aire estándar a temperatura [math]T_0 = 15^\circ \text{C}[/math].

6.5.1 Cálculo numérico de la masa total de aire dentro de la torre

Para realizar la integral primero despejamos de la ecuación del hiperboloide la componente [math]\rho[/math] y creamos una función [math]\rho(z)[/math] la cual se usara como limite de integración:

Luego por el jacobiano cambiamos [math]dV[/math] [math] = dx\,dy\,dz = r \, dr \, d\theta \, dz [/math]

Quedando la siguiente integral:

Donde:

[math] \rho(T) = 1.225 \, \frac{15 + 273}{T(r,z) + 273} [/math]

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

Para la resolución de la integral se realizaron aproximaciones por el método de Simpson, ya que la integral no tiene primitiva directa. dando como resultado numérico:

6.5.2 comparación con la masa si el aire estuviera uniformemente a temperatura [math]T_{\text{ambiente}} = 15^\circ \text{C}[/math]

Primero, la densidad del aire a 15°C:

[math] \rho_{15} = 1.225 \,\text{kg/m}^3 [/math]

Volumen del elipsoide truncado (de 0 a 150 m):

Ese volumen para el mismo dominio es:

pero podemos usar directamente la fórmula ejecutada:

Entonces:

6.6 Cálculo de la temperatura media del aire dentro de la torre

Para ello tomamos:

donde [math]V[/math] es el volumen interior.

6.6.1 ¿Cuánta potencia térmica se evacua?

Si el caudal de aire es [math]Q = 1500 \,\text{m}^3/\text{s}[/math]

con [math]c_p = 1005 \,\text{J}/(\text{kg}\cdot\text{K})[/math] el calor específico a presión constante del aire y [math]T_{\text{ambiente}} = 15^\circ \text{C}[/math].