Torres de enfriamientos hiperbólicas (Grupo 46)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | TORRES DE ENFRIAMIENTO HIPERBOLICAS |
| Asignatura | Teoría de Campos |
| Curso | 2025-26 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Modelo geométrico del hiperboloide
- 3 Determinación de los parámetros del modelo
- 4 Análisis de presión del viento
- 4.1 Opción 1: Análisis de presión del viento
- 4.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?
- 4.3 Calcula y representa el campo vectorial de fuerza \(\vec{F}\) sobre la superficie expuesta.
- 4.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.
- 4.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)?
- 5 Campo de temperatura y transferencia de calor
- 5.1 Representación del campo escalar de temperatura en MATLAB
- 5.2 Calculo del gradiente de temperatura [math]\nabla T ^\circ \text{C}[/math]
- 5.3 Representación del campo vectorial gradiente sobre el hiperboloide
- 5.4 Representación de superficies isotérmicas
- 5.5 Simulación de un caso en el que entra aire caliente a la base y cálculo numérico de la masa total de aire dentro de la torre
- 5.6 Cálculo de la temperatura media del aire dentro de la torre
- 6 Estructuras hiperbólicas famosas y ventajas frente a una superficie cilíndrica
- 7 Link al Poster
- 8 Bibliografía
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
El parámetro a controla la escala radial del hiperboloide en el plano horizontal xy. Específicamente, determina el radio de las secciones circulares que se obtienen al cortar la superficie con planos paralelos al plano xy. Cuanto mayor sea el valor de a, más ancho será el hiperboloide en cada nivel de altura z, lo que implica una expansión más pronunciada en las direcciones x y y. En términos geométricos, a define el tamaño mínimo de las circunferencias que forman la “cintura” del hiperboloide.
El parámetro c regula la apertura del hiperboloide en la dirección vertical z. Actúa como un factor de escala que controla cuán rápido se ensancha la superficie conforme nos alejamos del plano central z = z₀. Un valor grande de c produce una apertura más suave y extendida, mientras que un valor pequeño genera una forma más estrecha y cerrada. En esencia, c define la curvatura vertical del hiperboloide y afecta la forma de las hipérbolas que aparecen en las secciones longitudinales.
El parámetro z₀ representa el desplazamiento vertical del centro del hiperboloide. Es el valor de z en el cual se encuentra la sección más estrecha de la superficie, también conocida como el cuello o el eje central. Al modificar z₀, se traslada el hiperboloide hacia arriba o hacia abajo en el espacio tridimensional sin alterar su forma. Este parámetro es clave para posicionar la superficie en el contexto de un sistema de coordenadas, especialmente cuando se modelan estructuras físicas como torres de refrigeración o cubiertas arquitectónicas.
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. [math] \begin{cases} x = \rho \cos \theta \\ y = \rho \sin \theta \\ z = z \end{cases}[/math]
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]:
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:
% 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 Análisis de presión del viento
4.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).
4.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);4.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);4.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]
4.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.
5 Campo de temperatura y transferencia de calor
5.1 Representación del campo escalar de temperatura en MATLAB
(2.4) 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 en Matlab necesitamos hallar la temperatura en el centro de la altura máxima de la torre, es decir: [math]\rho=0 \;y\; z=150 [/math]
[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]
[math]T_{\text{top}} = T(0, 150)= 65 - 38 \left( \frac{150}{150} \right)^{1,8} - 8 \left( 1 - e^{-\frac{0^2}{55^2 - 0^2}} \right) = 65 - 38 = 27^ºC[/math]
5.1.1 Localización del punto más 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 (m)
a = 30;
c = 65.079;
z0 = 100;
H = 150;
Rmax = 55;
%Parámetros de temperatura (ºC)
Tbase = 65; % Temperatura en la base
Ttop = 27; % Temperatura en la parte superior
DeltaTz = 38; % Caída de temperatura desde la base hasta el tope
DeltaTr = 8; % Variación radial de temperatura
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;5.2 Calculo del gradiente de temperatura [math]\nabla T ^\circ \text{C}[/math]
(2.5) Calcula el gradiente de temperatura [math]\nabla T[/math]. Representa este campo vectorial sobre la misma sección vertical y sobre las secciones transversales seleccionadas en el apartado anterior. ¿En qué dirección se orienta de manera predominante?
Siendo el campo escalar: [math] T(\rho, z) = T_{\text{base}} - \Delta T_z \left( \frac{z}{H} \right)^n - \Delta T_{\rho} \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: [math] \nabla T = \frac{\partial T}{\partial \rho} \, \vec{e}_\rho + \frac{1}{\rho} \frac{\partial T}{\partial \theta} \, \vec{e}_\theta + \frac{\partial T}{\partial z} \, \vec{e}_z [/math]
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: [math]\nabla T = \frac{\partial T}{\partial \rho} \, \vec{e}_\rho + \frac{\partial T}{\partial z} \, \vec{e}_z[/math]
Derivada respecto a [math]\rho[/math]: [math]\frac{\partial T}{\partial \rho} = -\Delta T_{\rho} \cdot \frac{d}{d\rho} \left( 1 - e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}} \right) = \Delta T_{\rho} \cdot \frac{d}{d\rho} \left( e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}} \right) [/math]
Aplicamos regla de la cadena:[math]\frac{\partial T}{\partial \rho} = \Delta T_{\rho} \cdot e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}} \cdot \frac{d}{d\rho} \left( -\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}\right) = \Delta T_{\rho} \cdot e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}}\cdot \left(- \frac{ (2\rho)(R_{\text{max}}^2 - \rho^2) - (\rho^2)(-2\rho) }{ (R_{\text{max}}^2 - \rho^2)^2 }\right) = \Delta T_{\rho} \cdot e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}}\cdot \left( - \frac{ 2\rho R_{\text{max}}^2 }{ (R_{\text{max}}^2 - \rho^2)^2 } \right)
[/math]
Derivada respecto a z:
[math]
\frac{\partial T}{\partial z} = -\Delta T_z \cdot n \left( \frac{z}{H} \right)^{n-1} \cdot \frac{1}{H}
[/math]
Gradiente final: [math] \nabla T = \left[ \Delta T_r \cdot e^{-\frac{\rho^2}{R_{\text{max}}^2 - \rho^2}} \cdot \left( - \frac{ 2\rho R_{\text{max}}^2 }{ (R_{\text{max}}^2 - \rho^2)^2 } \right) \right] \vec{e}_\rho + \left[ -\Delta T_z \cdot n \left( \frac{z}{H} \right)^{n-1} \cdot \frac{1}{H} \right] \vec{e}_z [/math]
Gradiente habiendo sustituido los datos:
[math] \boxed{\nabla T = \Bigg[ 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}} \cdot \left( - \frac{ 2\rho \cdot 55^2 }{ (55^2 - \rho^2)^2 } \right) \Bigg] \vec{e}_\rho + \Bigg[ -38 \cdot 1,8 \cdot \left( \frac{z}{150} \right)^{0,8} \cdot \frac{1}{150} \Bigg] \vec{e}_z} [/math]
5.3 Representación del campo vectorial gradiente sobre el hiperboloide
En el siguiente recuadro se muestra el gradiente del campo de temperaturas en una sección vertical que pasa por el eje de simetría de la torre, más su respectivo código de Matlab para llevarlo a cabo.
% Parámetros del hiperboloide y el campo de temperatura
a = 30;
c = 65.079;
z0 = 100;
H = 150;
Rmax = 55;
%Parámetros de temperatura
Tbase = 65;
Ttop = 27;
DeltaTz = 38;
DeltaTr = 8;
n = 1.8;
% Definición de la malla
rho = linspace(-Rmax, Rmax, 100);
z = linspace(0, H, 100);
[R, Z] = meshgrid(rho, 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 = Tbase - DeltaTz * (Z / H).^n - DeltaTr*(1 - exp(-R.^2 ./ (Rmax^2 - R.^2)));
% Gradiente completo del campo de temperatura
[Tr, Tz] = gradient(T, rho, 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 rho y z)
figure;
quiver(R_sub, Z_sub, Tr_sub, Tz_sub, 'b-', '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('\nablaT respecto al plano vertical que pasa por el eje de simetría');
axis equal;
xlim([-Rmax ,Rmax]);
ylim([0 H]);
grid on;
A continuación se muestra el código de Matlab utilizado para representar, mediante una animación, cómo varía el gradiente del campo de temperaturas a través de secciones transversales.
% Parámetros geométricos (m)
a = 30;
c = 65.079;
z0 = 100;
H = 150;
Rmax = 55;
% Parámetros de temperatura (ºC)
Tbase = 65;
Ttop = 27;
DeltaTz = 38;
DeltaTr = 8;
n = 1.8;
%% Animación de planos transversales con gradiente 3D y borde circular
figure('Name','Animación transversal del gradiente con borde circular');
z_levels = linspace(0, H, 40); % cortes horizontales
% --- Superficie completa del hiperboloide ---
theta = linspace(0,2*pi,80);
z_surf = linspace(0,H,80);
[Theta,Zsurf] = meshgrid(theta,z_surf);
Rsurf = a*sqrt(1+(Zsurf-z0).^2/c^2);
Xsurf = Rsurf.*cos(Theta);
Ysurf = Rsurf.*sin(Theta);
surf(Xsurf,Ysurf,Zsurf,'FaceAlpha',0.1,'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);
hold on;
% Configuración fija de la vista
xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
view(45,30); axis equal; grid on;
xlim([-Rmax,Rmax]); ylim([-Rmax,Rmax]); zlim([0,H]);
axis vis3d;
for k = 1:length(z_levels)
% Borrar objetos anteriores del plano
delete(findobj(gca,'Tag','slice'));
delete(findobj(gca,'Tag','quivers'));
delete(findobj(gca,'Tag','border'));
zi = z_levels(k); % altura actual
% Radio máximo en esa altura según hiperboloide
rmax_slice = sqrt((1 + (zi - z0)^2 / c^2) * a^2);
[X,Y] = meshgrid(linspace(-rmax_slice,rmax_slice,40), ...
linspace(-rmax_slice,rmax_slice,40));
Rxy = sqrt(X.^2 + Y.^2);
% Máscara para limitar al círculo de radio rmax_slice
mask = Rxy <= rmax_slice;
% Campo de temperatura en el plano
T = Tbase - DeltaTz*(zi/H)^n ...
- DeltaTr*(1 - exp(-(Rxy.^2)./(Rmax^2 - Rxy.^2)));
T(~mask) = NaN;
% Gradiente completo en 3D
dx = X(1,2)-X(1,1);
dy = Y(2,1)-Y(1,1);
[Tx,Ty] = gradient(T,dx,dy);
dz = H/200;
T_up = Tbase - DeltaTz*((zi+dz)/H)^n ...
- DeltaTr*(1 - exp(-(Rxy.^2)./(Rmax^2 - Rxy.^2)));
T_down = Tbase - DeltaTz*((zi-dz)/H)^n ...
- DeltaTr*(1 - exp(-(Rxy.^2)./(Rmax^2 - Rxy.^2)));
Tz = (T_up - T_down)/(2*dz);
Tz(~mask) = NaN;
% Escalar para visibilidad
scale = 2 / max(abs([Tx(:); Ty(:); Tz(:)]));
Tx = Tx*scale; Ty = Ty*scale; Tz = Tz*scale;
% Dibujar plano coloreado
surf(X,Y,zi*ones(size(X)),T,'EdgeColor','none','FaceAlpha',0.3,'Tag','slice');
% Dibujar vectores del gradiente en 3D
skip = 4;
quiver3(X(1:skip:end,1:skip:end), ...
Y(1:skip:end,1:skip:end), ...
zi*ones(size(X(1:skip:end,1:skip:end))), ...
Tx(1:skip:end,1:skip:end), ...
Ty(1:skip:end,1:skip:end), ...
Tz(1:skip:end,1:skip:end), ...
'b','LineWidth',1.2,'Tag','quivers');
% --- Borde circular negro ---
theta_border = linspace(0,2*pi,200);
xb = rmax_slice*cos(theta_border);
yb = rmax_slice*sin(theta_border);
zb = zi*ones(size(xb));
plot3(xb,yb,zb,'k','LineWidth',2,'Tag','border');
title(sprintf('Plano transversal en z = %.2f m',zi));
drawnow;
pause(0.15);
endComo se puede observar en la imagen y cumpliendo con su propiedad, el campo gradiente apunta al punto de máximo crecimiento del campo escalar [math]T(\rho,z)[/math], que como se comprobó, en nuestro caso es [math]z=0[/math] y [math]\rho=0[/math].
5.4 Representación de superficies isotérmicas
(2.6) Representa varias superficies isotérmicas para distintos valores de temperatura (por ejemplo, [math]T = 30,40,50,60,\ldots \; ^\circ C[/math]). ¿Qué forma presentan? ¿Son aproximadamente horizontales o muestran una curvatura significativa?
Las superficies isotermas se describen a partir de la siguiente expresión del campo de temperatura:
Cada valor de temperatura define una superficie dependiente de la coordenada vertical y del radio. El resultado de esta ecuación genera formas que, en gran medida, se asemejan a una semiesfera para cada nivel térmico considerado, es decir, presentan una curvatura significativa.
Con el fin de facilitar la interpretación, se ha elaborado un código de Matlab que muestra una vista lateral de la figura con 6 superficies isotérmicas representadas como líneas.
% Parámetros del hiperboloide y el campo de temperatura
a = 30;
c = 65.079;
z0 = 100;
H = 150;
Rmax = 55;
%Parámetros de temperatura
Tbase = 65;
Ttop = 27;
DeltaTz = 38;
DeltaTr = 8;
n = 1.8;
% Definición de la malla
rho = linspace(-Rmax, Rmax, 100);
z = linspace(0, H, 100);
[R, Z] = meshgrid(rho, z);
% Campo de temperatura T(rho, z)
T = Tbase - DeltaTz * (Z / H).^n - DeltaTr*(1 - exp(-R.^2 ./ (Rmax^2 - R.^2)));
% Definir 5 valores isotérmicos (valores de temperatura)
T_values = linspace(min(T(:)), max(T(:)), 8);
% 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('\rho (m)');
ylabel('z (m)');
title('Superficies isotérmicas dentro del hiperboloide');
axis equal;
xlim([-Rmax Rmax]);
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]);5.5 Simulación de un caso en el que entra aire caliente a la base y cálculo numérico de la masa total de aire dentro de la torre
(2.7) En la base de la torre entra aire caliente. Considera que la densidad del aire depende de la temperatura según: [math] \delta(T) = \delta_{\text{aire}} \frac{T_0 + 273}{T + 273} [/math]
donde [math]\delta_{\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].
Calcula (analíticamente o numéricamente) la masa total de aire dentro de la torre: [math] M = \iiint_V \delta(T(\rho, z)) \, dV [/math]
y compárala con la masa si el aire estuviera uniformemente a temperatura [math]T_{\text{ambiente}} = 15^\circ \text{C}[/math].
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 \, d\rho \, d\theta \, dz [/math]
Quedando la siguiente integral:
Donde:
[math] \delta(T) = 1,225 \, \frac{15 + 273}{T(r,z) + 273} [/math]
[math] T(\rho,z) = 65 - 38 \left( \frac{z}{150} \right)^{1.8} - 8 \left( 1 - e^{-\frac{\rho^2}{55^2 - \rho^2}} \right) [/math]
Para la resolución de la integral se realizó un programa de Matlab, ya que la integral no tiene primitiva directa.
% Parámetros
z_max = 150;
% Definición de rho(z)
rho_lim = @(z) 30 * sqrt(1 + ((z - 100).^2) / (65.079^2));
% Función de temperatura T(rho, z)
T = @(rho, z) 65 - 38 * (z / 150).^1.8 ...
- 8 * (1 - exp(-rho.^2 ./ (552 - rho.^2)));
% Función de densidad delta(T)
delta = @(rho, z) 1.225 * 288 ./ (T(rho, z) + 273);
% Integrando: delta(T) * rho
integrando = @(rho, z) delta(rho, z) .* rho;
% Integral doble con límite superior dependiente de z
M = 2 * pi * integral(@(z) arrayfun(@(zz) ...
integral(@(rho) integrando(rho, zz), 0, rho_lim(zz)), z), ...
0, z_max);
% Mostrar resultado
fprintf('La masa total del aire dentro de la torre es %.6f kg\n', M);5.5.1 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] \delta(T) = 1,225 \, \frac{15 + 273}{15 + 273} = 1,225 \, \text{kg/m}^3 [/math]
Volumen del hiperboloide truncado (de 0 a 150 m):
Obteniendo como resultado:
[math]
V \approx 674461,7295\,\text{m}^3
[/math]
Entonces:
Debido a que [math] M_{\text{aire caliente}} \lt M_{15} [/math] , podemos decir que en el caso de temperatura uniforme a 15ºC se sobreestima la masa y en el otro caso se captura la variación espacial de densidad y ofrece un valor más realista.
5.6 Cálculo de la temperatura media del aire dentro de la torre
Para ello tomamos:
donde [math]V[/math] es el volumen interior.
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 \, d\rho \, d\theta \, dz [/math]
Quedando la siguiente integral:
Donde:
[math] T(\rho,z) = 65 - 38 \left( \frac{z}{150} \right)^{1.8} - 8 \left( 1 - e^{-\frac{\rho^2}{55^2 - \rho^2}} \right) [/math]
Para la resolución de la integral se realizó un programa de Matlab, ya que la integral no tiene primitiva directa.
% Límites
z_min = 0;
z_max = 150;
rho_cap = 55; % singularidad a rho = 55
eps_r = 1e-6; % margen de seguridad
% Radio límite dependiente de z
rho_lim_nom = @(z) 30 .* sqrt(1 + ((z - 100).^2) ./ (65.079^2));
rho_lim = @(z) min(rho_lim_nom(z), rho_cap - eps_r); % evitar tocar 55
% Campo de Temperatura
T = @(rho, z) 65 - 38 * (z / 150).^1.8 ...
- 8 * (1 - exp(-rho.^2 ./ (3025 - rho.^2)));
% Volumen V = pi * ∫ rho(z)^2 dz
rho2 = @(z) (rho_lim(z)).^2;
V = pi * integral(@(z) rho2(z), z_min, z_max, ...
'RelTol',1e-8,'AbsTol',1e-10);
% Integral triple de T sobre el volumen (en cilíndricas)
integrando = @(rho, z) T(rho, z) .* rho;
IT = 2 * pi * integral(@(z) arrayfun(@(zz) ...
integral(@(rho) integrando(rho, zz), 0, rho_lim(zz), ...
'RelTol',1e-8,'AbsTol',1e-10), z), ...
z_min, z_max, 'RelTol',1e-8,'AbsTol',1e-10);
Tmedia = IT / V;
fprintf('V = %.9f\n', V);
fprintf('Integral de T sobre V (IT) = %.9f\n', IT);
fprintf('Tmedia = %.9f ºC\n', Tmedia);Como resultado final: [math]\boxed{T_{\text{media}} \approx 51,7845^ºC}[/math]
5.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].
Una vez calculada la temperatura media se calcula la potencia, ajustamos las unidades:
La variacion de temperatura es igual en kelvin que en centígrados por lo que no nos hace falta hacer el cambio.
El resto de unidades concuerdan ya que según la formula:
[math]\frac{\text{kg}}{\text{m}^3} \times \frac{\text{J}}{\text{kg}\cdot\text{K}} \times \frac{\text{m}^3}{\text{s}} \times \text{K} = \frac{\text{J}}{\text{s}} = \text{W}.[/math]
6 Estructuras hiperbólicas famosas y ventajas frente a una superficie cilíndrica
6.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.
6.2 Torres de enfriamiento hiperbólicas famosas
Central Nuclear de Cofrentes (España): Sus torres hiperbólicas son un referente en el paisaje industrial español, diseñadas para enfriar grandes volúmenes de agua mediante tiro natural.
Didcot Power Station (Reino Unido): Las icónicas torres, demolidas en 2014, fueron un símbolo de la ingeniería británica y ejemplo del uso estandarizado de hiperboloides.
Niederaussem Power Station (Alemania): Alberga algunas de las torres de enfriamiento más grandes del mundo, con alturas superiores a 200 metros, que aprovechan la geometría hiperbólica para maximizar la eficiencia.
6.3 Otras estructuras hiperboloides en arquitectura e ingeniería
• 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.
7 Link al Poster
8 Bibliografía
Apuntes: TeoríadeCampos,cursoacadémico2025–2026,ETSICaminos,UniversidadPolitécnicadeMadrid. ©2025JuanAntonioBarceló,CarlosCastro,DavidGonzález,NicoMicheleSchiavone
