Torres de enfriamiento hiperbólicas (Grupo 06)

Trabajo realizado por estudiantes
Título	Torres de enfriamiento hiperbólicas (Grupo 06)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Álvaro Díaz Martín Ramón Castán Naval Manuel Álvarez-Campana Illescas Gonzalo Quemada Simón Fernando Barettino Moreno-Aurioles
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Las torres de enfriamiento hiperbólicas, son estructuras emblemáticas del paisaje energético. Estas se pueden encontrar en instalaciones como la Central Nuclear de Cofrentes (España) o en Niederaussem (Alemania). Son un elemento clave para la disipación térmica en el ciclo de producción de energía. Su diseño de hiperboloide de revolución de una hoja acelera el flujo de aire ascendente y enfría el agua de manera natural sin necesidad de usar ventiladores. La doble curvatura proporcionada por esta estructura hiperbólica otorga a la torre una enorme resistencia frente al viento y su propio peso.

En el desarrollo de este articulo vamos a analizar una torre de enfriamiento estándar de altura máxima ([math]H[/math]) con radio máximo ([math]Rmáx[/math]) y radio mínimo ([math]Rmín[/math]) el cual se encuentra a [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math].La superficie de la torre se describe matemáticamente mediante la siguiente ecuación en coordenadas cartesianas:

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

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

Sean en nuestro modelo [math]Rmáx=55m[/math] ; [math]Rmín=30m[/math] ; [math]H=150m[/math]

2 La geometría hiperbólica en las torres de enfriamiento

Como se menciona anteriormente, estas estructuras siguen una morfología hiperbólica. La coordenada [math]z[/math] esta comprendida desde [math]z=0[/math] y [math]z=150[/math] al ser la altura máxima ([math]H=150[/math]). Al ser la figura un solido de revolución, generado al girar una hipérbola alrededor del eje Z , tiene simetría radial.

2.1 Obtención de los parámetros

Como se ha visto anteriormente el parámetro [math]z_0[/math] es la posición vertical donde ocurre el estrechamiento, el cual está localizado a una altura de [math]\dfrac{\scriptsize 2}{\scriptsize 3}H[/math].

[math]z_0=\dfrac{\scriptsize 2}{\scriptsize 3}H=\dfrac{\scriptsize 2}{\scriptsize 3}150=100[/math]

El parámetro [math]a[/math] representa el radio de la estrechez del hiperboloide, el cual es su radio mas pequeño [math]Rmín=30m[/math].

[math]a=30[/math]

El parámetro [math]c[/math] es el factor de estiramiento vertical. Esta variable controla la curvatura de las paredes. Determina que tan rápido se ensancha la torre a medida que se alejas de la garganta.

Para calcular el parámetro [math]c[/math] partimos de la siguiente ecuación [math](a=30,z_0=100)[/math]

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

Se necesita otro punto conocida para la obtención de [math]c[/math]. Se conoce que la base se encuentra a una altura [math]z=0[/math] y que tiene un radio [math]R=55[/math].

Cualquier punto del borde de la base cumple que [math]x^2+y^2=55^2[/math]. Para facilitar los cálculos se toma el siguiente punto de la base [math](x,y)=(55,0)[/math]. Se sustituye con este punto en la ecuación inicial.

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

Simplificando la ecuación

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

Elevando al cuadrado:

[math]\dfrac{3025}{900}-\dfrac{10000}{c^2} = 1[/math]

Despejando [math]c[/math]:

[math]3,361-\dfrac{10000}{c^2} = 1[/math] [math]\longrightarrow[/math] [math]3,361-1 = \dfrac{10000}{c^2}[/math]

[math]2,361= \dfrac{10000}{c^2}[/math] [math]\longrightarrow[/math] [math]c^2=\dfrac{10000}{2,361}[/math]

[math]c=\sqrt{\dfrac{10000}{2,361}}[/math] [math]\longrightarrow[/math] [math]c=65,07914[/math]

Los parámetros obtenidos son [math]a=30,z_0=100,c=65,07914[/math]

2.2 Ecuación de la superficie hiperbólica en coordenadas cilíndricas

Para obtener la ecuación en coordenadas cilíndricas se parte de la anterior expresión en coordenadas cartesianas:

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

Sean [math]a, c, z_0\gt0[/math] los valores determinados en el apartado anterior.

Con el objetivo de pasar la ecuación a coordenadas cilíndricas se definen los parámetros [math](x,y,z)[/math] en [math](ρ,θ,z).[/math]

[math]x=ρ\cdot\cos(θ)[/math]

[math]y=ρ\cdot\sin(θ)[/math]

[math]z=z[/math]

Cambiando [math]x^2+y^2[/math] a coordenadas cilíndricas:

[math]x^2+y^2=(ρ\cdot\cos(θ))^2+(ρ\cdot\sin(θ))^2[/math]
[math]x^2+y^2=ρ^2\cdot\cos^2(θ)+ρ^2\cdot\sin^2(θ)[/math]

Se saca factor común [math]ρ^2[/math]:

[math]x^2+y^2=ρ^2\cdot(\cos^2(θ)+\sin^2(θ))=1[/math]

Aplicando la siguiente identidad trigonométrica:

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

Sustituyendo en la ecuación

[math]x^2+y^2=ρ^2\cdot 1 \longrightarrow x^2+y^2=ρ^2[/math]

3 Campo de temperaturas y transferencia de calor

El campo de temperaturas del aire en el interior de la torre se modeliza mediante la siguiente expresión:

[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]

Siendo:

• [math]T_{\text{base}}[/math] la temperatura en el centro de la base, en este caso [math]T_{\text{base}}=65ºC[/math];
• [math]\Delta T_z[/math] la caída de temperatura de base a tope, en este caso [math]\Delta T_z=38ºC[/math];
• [math]\Delta T_{\rho}[/math] la variación radial de temperatura, en este caso [math]\Delta T_{\rho}=8ºC[/math];
• [math]n[/math] el exponente de convección, en este caso [math]n=1,8[/math].

Sustituyendo estos valores en el campo y simplificando se llega a la siguiente expresión:

[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]

[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]

3.1 Cálculo de la masa del aire contenido en la torre

La densidad del aire cambia varía con la temperatura, por ello y dadas las dimensiones de la torre de enfriamiento, hay que tener en cuenta esta variación para el correcto cálculo de la masa. La expresión que relaciona la densidad con la temperatura es la siguiente:

[math]\rho(T) = \rho_{\text{aire}} \frac{T_0 + 273}{T + 273}[/math]

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

De este modo, se calcula la masa del aire dentro de la torre mediante la siguiente expresión:

[math]M_{\text{real}} = \iiint_{V} \rho(T(\rho, z)) \, dV=\int_{0}^{2\pi} \int_{0}^{150} \int_{0}^{\rho(z)} \rho(T(\rho, z)) \cdot \rho \, d\rho \, dz \, d\theta[/math]

Si se desarrolla la siguiente expresión se llega a lo siguiente:

[math]M_{\text{real}} = 2\pi \cdot 352.8 \int_{0}^{150} \int_{0}^{\rho(z)} \frac{\rho}{330 - 38 \left(\frac{z}{150}\right)^{1.8} + 8 \cdot e^{- \frac{\rho^2}{55^2 - \rho^2}}} \, d\rho \, dz[/math]

Esta integral no es de fácil resolución, luego utilizando una herramienta de internet para realizar el cálculo numérico se llega al siguiente resultado:

[math] M_{\text{real}} \approx 781,420kg[/math]

Por otro lado, si se considera que todo el aire está a una temperatura [math]T=15ºC[/math] su densidad es constante y de valor [math]\rho_{\text{aire}}=1,225kg/m^3[/math]. Por lo tanto, la masa se calcula multiplicando dicha densidad por el volumen de aire contenido en la torre. Por lo que el problema pasa a ser calcular dicho volumen.

El volumen de la torre se calcula integrando el jacobiano en coordenadas cilíndricas en la región requerida:

[math]V = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} \rho \, d\rho dz d\theta[/math]

Si se escribe la ecuación de un hiperboloide de revolución en coordenadas cilíndricas se obtiene esta expresión:

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

Despejando se obtiene una relación entre [math]\rho[/math] y [math]z[/math]

[math]\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right][/math]

4 Presión del viento

El viento es un factor que ejerce presión lateral sobre la torre de enfriamiento que varían a lo largo de la superficie en función de la altura. La velocidad escalar se puede representar de la siguiente forma:

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

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

[math]Zref[/math] es la altura de referencia sobre el suelo, que es [math]10m[/math].

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

4.1 .-Representación del campo escalar de presiones.

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

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

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

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

[math]P(z)=\dfrac{1}{2}ρ·V(z)^2=\dfrac{1}{2}ρ·(V_0·(\dfrac{z}{Zref})^α)^2=\dfrac{1,225}{2}·(18·(\dfrac{z}{10})^{0,1429})^2≈102.79·z^{0,286}[/math]

Haremos el supuesto de que la torre es golpeada por un viento paralelo al vector [math]-\frac{1}{\sqrt{2}}(\vec{i} + \vec{j})[/math] para la mitad expuesta de la torre, se adjunta a continuación el código usado en Matlab y la representación obtenida: IMAGENPRESIONES

a = 20; z0 = 80; c = 34.91;
V0 = 18;
zref = 10;
alpha = 0.1429;
rho_air = 1.225; %[kg/m3], (dato extraido de worldpress.com)
v = linspace(0, 150, 100);
u = linspace(0, pi, 100); 
[U, V] = meshgrid(u, v);
Rho = a * sqrt(1 + ((V - z0) / c).^2);
X = Rho .* cos(U);
Y = Rho .* sin(U);
P_z = 102.79 * V.^0.286;
figure;
surf(X, Y, V, P_z, 'EdgeColor', 'none'); 
colormap(jet); 
colorbar;axis equal; 
title('Mapa de Presión del Viento en la Mitad de la Torre Expuesta');
xlabel('X (m)');ylabel('Y (m)');zlabel('Z (m)');

4.2 Campo de fuerza en la superficie expuesta

La fuerza ejercida por el viento sobre la superficie actúa en la dirección normal a ella y el campo vectorial se representa de la siguiente forma:

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

A continuación se adjunta la imagen con la representación que la fuerza creada por la presión del viento sobre la superficie de la torre, y el código de MATLAB con el que se he representado:

Torre de enfriamiento con presión del viento lateral

a = 20; c = 34.91; z0 = 80; 
H = 150;zref=10;
theta = linspace(0, pi, 18);
z = linspace(0, H, 18);
[Theta, Z] = meshgrid(theta, z);
R = a * sqrt(1 + ((Z - z0).^2 / c^2));
X = R .* cos(Theta);
Y = R .* sin(Theta);
P_z = 102.79 * Z.^0.286
dX_dtheta = -R .* sin(Theta);% Derivadas parciales de la parametrización
dY_dtheta = R .* cos(Theta);
dZ_dtheta = zeros(size(Z));
dX_dz = a * (Z - z0) ./ (c^2 * sqrt(1 + ((Z - z0).^2 / c^2))) .* cos(Theta);
dY_dz = a * (Z - z0) ./ (c^2 * sqrt(1 + ((Z - z0).^2 / c^2))) .* sin(Theta);
dZ_dz = ones(size(Z)); % Derivada de Z respecto de Z
NX = dY_dtheta .* dZ_dz - dZ_dtheta .* dY_dz;% Producto vectorial para obtener el vector normal
NY = dZ_dtheta .* dX_dz - dX_dtheta .* dZ_dz;
NZ = dX_dtheta .* dY_dz - dY_dtheta .* dX_dz;
Norm = sqrt(NX.^2 + NY.^2 + NZ.^2);% Normalización de los vectores normales
NX = NX ./ Norm;  NY = NY ./ Norm;  NZ = NZ ./ Norm;
% Fuerza descompuesta generada por la presión del aire
FX = -P_z .* NX;  FY = -P_z .* NY;  FZ = -P_z .* NZ;
%representación del campo vectorial de fuerzas sobre dicha superfie:
quiver3(X, Y, Z, FX, FY, FZ, 1.5, 'r');
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Campo vectorial de fuerza sobre media torre');
axis equal;