MateWiki - Contribuciones del usuario [es]

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-03T17:52:56Z

Ramón Castán Naval: /* Campo de fuerza en la superficie expuesta */

{{ TrabajoED |Torres de enfriamiento hiperbólicas (Grupo 06) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] |
*Álvaro Díaz Martín
*Ramón Castán Naval
*Manuel Álvarez-Campana Illescas
*Gonzalo Quemada Simón
*Fernando Barettino Moreno-Aurioles }}
[[Categoría:Teoría de Campos]][[Categoría:TC25/26]] 

== 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:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</math> unos valores a determinar
Sean en nuestro modelo <math>Rmáx=55m</math> ; <math>Rmín=30m</math> ; <math>H=150m</math>

== 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=150m</math>). Al ser la figura un solido de revolución, generado al girar una hipérbola alrededor del eje Z , tiene simetría radial.

=== 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>.

<center><math>z_0=\dfrac{\scriptsize 2}{\scriptsize 3}H=\dfrac{\scriptsize 2}{\scriptsize 3}150=100</math></center>
 
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>.

<center><math>a=30</math></center>
 
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>
 
<center><math>\dfrac{x^2+y^2}{30^2}-\dfrac{(z-100)^2}{c^2} = 1</math></center>
 
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.
 
<center><math>\dfrac{55^2+0^2}{30^2}-\dfrac{(0-100)^2}{c^2} = 1</math></center>
 
Simplificando la ecuación
 
<center><math>\dfrac{55^2}{30^2}-\dfrac{(-100)^2}{c^2} = 1</math></center>
 
Elevando al cuadrado:
<center><math>\dfrac{3025}{900}-\dfrac{10000}{c^2} = 1</math></center>
 
Despejando <math>c</math>:
 
<center><math>3,361-\dfrac{10000}{c^2} = 1</math> <math>\longrightarrow</math> <math>3,361-1 = \dfrac{10000}{c^2}</math></center>
 
<center><math>2,361= \dfrac{10000}{c^2}</math> <math>\longrightarrow</math> <math>c^2=\dfrac{10000}{2,361}</math></center>
 
<center><math>c=\sqrt{\dfrac{10000}{2,361}}</math> <math>\longrightarrow</math> <math>c=65,07914</math></center>
 
Los parámetros obtenidos son <math>a=30,z_0=100,c=65,07914</math>

===Ecuación de la superficie hiperbólica en coordenadas cilíndricas===
Para obtener la ecuación en coordenadas cilíndricas se parte de la expresión anterior en coordenadas cartesianas:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</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 de las coordenadas cartesianas <math>(x,y,z)</math> en coordenadas cilíndricas <math>(ρ,θ,z).</math>
 
<center><math>x=ρ\cdot\cos(θ)</math> 
<math>y=ρ\cdot\sin(θ)</math> 
<math>z=z</math></center>
 
Cambiando <math>x^2+y^2</math> a coordenadas cilíndricas:
 
<center><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></center>
 
Se saca factor común <math>ρ^2</math>:
 
<center><math>x^2+y^2=ρ^2\cdot(\cos^2(θ)+\sin^2(θ))=1</math></center>
 
Aplicando la siguiente identidad trigonométrica:
 
<center><math>\cos^2(θ)+\sin^2(θ)=1</math></center>
 
Sustituyendo en la ecuación
 
<center><math>x^2+y^2=ρ^2\cdot 1 \longrightarrow x^2+y^2=ρ^2</math></center>
 
Al ser los demás datos constantes, no hace falta cambiar nada más. La ecuación en coordenadas cilíndricas queda: 
<center><math>\dfrac{ρ^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
===Representación grafica de la torre hiperbólica ===
 
A partir de la ecuación anterior se modela la torre de enfriamiento estándar de altura <math>(H=150m)</math>, radio máximo <math>(Rmáx=55m)</math> y radio mínimo <math>(Rmín=30m)</math>

 
[[Archivo:Torre11.jpg|miniaturadeimagen|500px|derecha|Representación grafica de la torre de enfriamiento estándar ]]
<syntaxhighlight lang="matlab">
clc
clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
u=linspace(0,2*pi,80);
v=linspace(0,H,80);

[Mt,Mz]=meshgrid(u,v);
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

figure
Mx=Mr.*cos(Mt)
My=Mr.*sin(Mt);
Mz=Mz;

surf(Mx,My,Mz, 'FaceColor','white','EdgeColor','black')
axis equal
hold on
grid on
title('Torre de enfriamiento hiperbólica')
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
</syntaxhighlight>

== 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: 
 
<center><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></center>
 
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: 
 
<center><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></center>
 
<center><math>T(\rho, z) = 65 - 38 \left(\frac{z}{150}\right)^{1,8} - 8 + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
<center><math>T(\rho, z) = 57 - 38 \left(\frac{z}{150}\right)^{1,8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
===Representación del campo escalar 𝑇(𝜌,𝑧) como un mapa de colores en una sección vertical y varias secciones transversales===
Se procede a representar el campo escalar de temperaturas como un mapa de colores con una interfaz que simula la de una cámara térmica.
 
 
'''Interpretación'''
 
Como se aprecia, la temperatura más elevada se da en el origen de coordenadas cartesianas <math>(X,Y,Z)=(0,0,0)</math>.
 
* La temperatura en el interior de la torre disminuye según aumenta la coordenada z.
 
* A su vez, la temperatura disminuye según aumenta el alejamiento al eje z, es decir, según aumenta la coordenada ρ (en coordenadas cilíndricas).
 
Por tanto, se aprecia que la variación en las temperaturas es radial y en altura. Este resultado es coherente, teniendo en cuenta que la función escalar de temperatura depende únicamente de los parámetros en coordenadas cilíndricas: z y ρ (y no de θ).
 
 
[[Archivo:grupo_6_campo_temperaturas_vertical.jpg|500px|Representación grafica del campo escalar de temperatura (sección vertical)]]
[[Archivo:grupo_6_campo_temperaturas_transversal.jpg|500px|Representación grafica del campo escalar de temperatura (secciones transversales)]]
 
A continuación se muestra el código en matlab que muestra las dos ventanas gráficas:
<syntaxhighlight lang="matlab">
clc; clear
%se establecen los parámetros de la torre:

H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:

Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIÓN VERTICAL-
z = linspace(0,H,300); %Se establecen los puntos del plano de sección XOZ
x = linspace(-Rmax,Rmax,300);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=Mrho; %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se representa el campo escalar Temperatura interior del plano XOZ de la
%torre
figure('name','Campo de Temperatura en sección XOZ de la torre')
contourf(Mx,Mz,T,60,'LineColor','none');
colormap('hot');
colorbar;
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo escalar de temperaturas en la sección vertical XOZ de la torre:');
tit.Position=[0,156,0]; %se edita la posición del titulo
hold on;

%Se representa la silueta de la torre
contorno_izq = -Mrho; %Lado izquierdo de la silueta
contorno_dch = Mrho; %Lado derecho de la silueta
contorno_alt = z; %Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); %Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); %Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2);%Vista en 2D
xlim([-Rmax-15,Rmax+15]);%se ajusta la ventana gráfica
ylim([0,H]);
grid off;

%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='Temp. (ºC)';
hold off;

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIONES TRANSVERSALES-
figure('name','Campo de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
rango=linspace(-Rmax,Rmax,200);
[Mx, My] = meshgrid(rango,rango);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);
for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre.
Zi=(H/3)*i; %Cuatro cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));
% T(rho, z)=T=Tbase-dTz.*(Zi./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2)))

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>Rho_z)=NaN;

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1)
pcolor(Mx, My, Ti);
shading interp; %Se suavizar el color

colormap('hot');
title(['Temperatura (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2);%Vista en 2D
xlim([-70,70]);%se ajusta la ventana gráfica
ylim([-70,70]);

hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='T (ºC)';
colorbar;
end
hold off;
</syntaxhighlight>

=== Cálculo y representación del gradiente del campo escalar Temperatura ===

Con el fin de interpretar la variación de temperaturas dentro de la torre de enfriamiento hiperbólica en cada punto, se estudia el campo vectorial gradiente del campo escalar temperatura.
 
Este campo gradiente tiene la siguiente expresión:
<center><math>∇T(r,θ,z)=(\dfrac{dT}{dρ},\dfrac{1}{ρ}·\dfrac{dT}{dθ},\dfrac{dT}{dz})</math></center>
 
Como el campo escalar temperatura no depende de <math>θ</math>,el gradiente de temperatura tiene solo dos componentes:
 
<center>Componente radial: <math>\dfrac{dT}{dρ}</math></center>
 
<center>Componente vertical: <math>\dfrac{dT}{dz}</math></center>
 
De este modo, se calculan las derivadas parciales respecto a <math>ρ</math> y respecto a <math>z</math>:
 
<center><math>\dfrac{dT}{dρ}=-ΔT_r·\dfrac{-2ρ·Rmax^2}{(Rmax^2-ρ^2)^2}·e^{(\dfrac{-ρ^2}{Rmax^2-ρ^2})}</math></center>
 
<center><math>\dfrac{dT}{dz}=-ΔT_z·n·\dfrac{z^{n-1}}{H^n}</math></center>
 
 
A continuación se muestra la representación del campo vectorial gradiente del campo Temperaturas <math>T(ρ,z):</math>
 
[[Archivo:grupo_6_gradiente_vertical.jpg|500px|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_gradiente_transversal.jpg|500px|Representación grafica del campo escalar de temperatura ]]
 
 
'''Interpretación gráfica y analítica'''
 
* El vector gradiente en cada punto apunta hacia donde la temperatura aumenta más rápidamente (el vector gradiente en cada punto sigue la trayectoria más empinada de ascenso de temperatura). Es por eso que todas los vectores apuntan hacia abajo (coordenada z=0) y hacia el eje de revolución torre (coordenadas x,y=0,0), ya que, como se ha expuesto en el apartado anterior, la temperatura en la torre aumenta según se acerca a la base y al centro de la torre, y por eso todos los vectores gradiente apuntan hacia allí.
 
 
* Otra interpretación importante es que los vectores tienen mayor módulo en los puntos más alejados de la base de la torre. Esto es así debido según se aumenta la altura, hay mayor diferencia de temperatura para la misma distancia. En la cima de la torre el aire caliente entra en contacto con el aire atmosférico, por tanto, cuanto más cerca se esté de este punto, mayor será el enfriamiento. Análogamente, según se acerca a las paredes de la torre, el gradiente tiene mayor módulo, puesto que para un mismo recorrido radial, la variación de temperatura es mayor.
 
A continuación se muestra el código en 'matlab' que genera sendas ventanas gráficas:
<syntaxhighlight lang="matlab">
clc, clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se establecen los puntos del plano de sección XOZ
z = linspace(0,H,25);
x = linspace(-Rmax,Rmax,25);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=(Mrho+4.*ones(25)); %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se calcula el campo gradiente de T
[gradx,gradz]=gradient(T,x,z);

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN VERTICAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOZ')
quiver(x,z,gradx,gradz,'r','LineWidth',0.7)
hold on

%Se representa la silueta de la torre
contorno_izq = -Mrho; % Lado izquierdo de la silueta
contorno_dch = Mrho; % Lado derecho de la silueta
contorno_alt = z; % Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); % Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); % Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2); %Vista en 2D
xlim([-Rmax-15,Rmax+15]); %se ajusta la ventana gráfica
ylim([0,H]);
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo vectorial gradiente de temperatura (Sección Vertical):');
tit.Position=[0,153,0]; %se edita la posición del titulo
hold off

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN TRANSVERSAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
x=linspace(-Rmax-10,Rmax+10,20);
y=x;
[Mx,My] = meshgrid(x,y);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);

for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre
Zi=(H/3)*i; %4 cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>(Rho_z))=NaN;

%Se calcula el gradiente
[gradx,grady]=gradient(Ti,x,y);

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1);
quiver(Mx,My,gradx,grady,'r','LineWidth',0.7);

title(['Gradiente(T) (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2); %Vista en 2D
xlim([-70,70]); %se ajusta la ventana gráfica
ylim([-70,70]);
hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
end
hold off;
</syntaxhighlight>

=== Representación de superficies isotérmicas ===
Las curvas de nivel son superficies generadas a partir de la asignación de un valor constante a un campo escalar (en <math>ℝ^3</math>). En este caso, el campo escalar es la Temperatura, T(ρ,z), y sus curvas de nivel son las 'Superficies Isotermas', obtenidas al asignar a T un valor de temperatura constante.
 
A continuación se representan las superficies isotérmicas para los valores de 30,40,50 y 60 ºC.
 
[[Archivo:grupo_6_isotermas_vertical.jpg|500px|izquierda|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_isotermas_xyz.jpg|500px|derecho|Representación grafica del campo escalar de temperatura ]]
 
'''Interpretación de las superficies de nivel:'''
 
* Las superficies isotermas tienen una forma arqueada hacia abajo, es decir, que son convexas hacia el eje z. Esto se debe a que el campo escalar temperatura (<math>T(ρ,z)</math>), como se ha comentado anteriormente, varía radial y verticalmente, disminuyendo la temperatura según se alejan los puntos del eje z (de revolución) y de la base de la torre. De este modo, para una misma altura (fijando z), la temperatura será mayor en el centro de la sección transversal que en los extremos, como bien se puede apreciar en la segunda figura del apartado 3.1.
* Es importante destacar que, según se asciende en la torre, las curvas isotermas pasan de asemejarse a una esfera a asemejarse a un plano horizontal, gradualmente. La explicación es debida a que el campo escalar Temperatura depende de dos términos, que dependen respectivamente de z y ρ.
 
El término que depende de z es:
 
<center><math>T(z) = -\Delta T_z \left(\frac{z}{H}\right)^n</math></center>
El término que depende de ρ es:
<center><math>T(\rho) = - \Delta T_{\rho} \left(1 - e^{- \frac{\rho^2}{R_{\text{max}}^2 - \rho^2}}\right)</math></center>
 
Por tanto:
* En la base de la torre (z cercano a 0), la superficie isoterma se asemeja a una esfera. Esto se debe a que el término <math>T(z)</math>, es mínimo, por lo que <math>T(ρ)</math> domina el campo escalar. Así, este último término da a las superficies isotermas la forma esférica, ya que el radio (ρ) produce el efecto más importante sobre la temperatura.
* Análogamente, en la cima de la torre (z cercano a H), la superficie se puede asemejar a un plano horizontal, puesto que el término <math>T(z)</math> es máximo. Así, el efecto que produce el radio (ρ), es menor.
 
A continuación se muestra el programa de matlab que representa en una ventana gráfica las superficies isotermas.
<syntaxhighlight lang="matlab">
clc,clear;
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;
%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38; %variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8; %variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se crean los puntos x,y,z
x = linspace(-Rmax,Rmax,200);
y = linspace(-Rmax,Rmax,200);
z = linspace(0,H,200);

[Mx,My,Mz]=meshgrid(x,y,z); %se crea la malla de puntos
Rho=sqrt(Mx.^2+My.^2); %Coordenada rho (cilíndricas)

%se calcula el campo escalar en la extensión de Rho, Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se calcula el radio de la torre para cada z
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

%Se eliminan los puntos que no están DENTRO de la torre (rho<=R(z))
Mask_interior=(Rho<=Mr);

%Se aplica la máscara: asignamos NaN a las temperaturas fuera de la torre.
T(~Mask_interior) = NaN;

%Se representa
figure('Name', 'Superficies Isotérmicas dentro de la Torre Hiperbólica');
hold on;

%Se define el colormap para los colores
clim([Tcima-5,Tbase]);
cmap = hot(45);
colormap(cmap);
cbar.Direction='reverse';

%% BUCLE
Temperaturas=[30,40,50,60]; %Se definen las temperaturas de las superficies
for i=1:4;
T_nivel=Temperaturas(i);
view(3); % Vista 3D
colorbar
%Se representa la superficie isoterma
[caras, vertices] = isosurface(Mx, My, Mz, T, T_nivel);
curva_nivel=patch('faces', caras, 'Vertices', vertices, 'FaceColor', cmap(T_nivel-23,:), 'EdgeColor', 'none');

%Se representa el contorno de la torre (con transparencia)
contorno = Rho - Mr;
[caras_t, vertices_t] = isosurface(Mx, My, Mz, contorno, 0);
h_tower = patch('faces', caras_t, 'vertices', vertices_t, 'FaceColor', [0.5 0.5 0.5], ...
'EdgeColor', 'none', 'FaceAlpha', 0.1); % Gris muy transparente

%Se ajusta de la ventana gráfica
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
axis equal tight;
end
%vista final
hold off
title('Superficies Isotérmicas dentro de la Torre');
lighting gouraud
camlight;
</syntaxhighlight>

=== Cálculo de la masa del aire contenido en la torre ===
La densidad del aire 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:
 
<center><math>\rho(T) = \rho_{\text{aire}} \frac{T_0 + 273}{T + 273}</math></center>
 
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:
 
<center><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></center>
 
Si se desarrolla la expresión se llega a lo siguiente:
 
<center><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></center>
 
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:
<center><math> M_{\text{real}} \approx 781420 \text{ kg}</math></center>
 
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:
 
<center><math>V = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} \rho \, d\rho dz d\theta</math></center>
En la sección 2.2 se obtiene la expresión del hiperboloide de revolución en coordenadas cilíndricas:
 
<center><math>\frac{\rho^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1</math></center>
Despejando se obtiene una relación entre <math>\rho</math> y <math>z</math> 
 
<center><math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math></center>
Volviendo a la integral, integrando primero <math>\theta</math> y después <math>\rho</math> se obtiene lo siguiente:
 
<center><math>V = 2\pi \int_0^{150} \left[ \frac{\rho^2}{2} \right]_0^{\rho(z)} dz = \pi \int_0^{150} \rho^2(z) dz</math></center>
 
Sustituyendo <math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math> con <math>z_0=100</math>:
 
<center><math>V = \pi a^2 \int_0^{150} \left[ 1 + \frac{(z - 100)^2}{c^2} \right] dz</math>
 
<math>V = \pi a^2 \left[ \int_0^{150} dz + \frac{1}{c^2} \int_0^{150} (z - 100)^2 dz \right]</math></center>
 
Resolviendo ambas integrales:
 
<center><math>\int_0^{150} dz = [z]_0^{150} = 150</math></center>
 
Para la realización de <math>\int_0^{150} (z - 100)^2 dz</math> usamos el cambio de variable <math>u = z - 100, du=dz</math> con límites <math>u(0)=-100, u(150)=50</math>
 
<center><math>\int_{-100}^{50} u^2 du = \left[ \frac{u^3}{3} \right]_{-100}^{50} = \frac{1}{3} \left[ 50^3 - (-100)^3 \right] = \frac{1}{3} [125000 + 1000000] = \frac{1125000}{3} = 375000</math></center>
 
Sustituyendo en la integral, con <math>a=30</math> y <math>c^2 =\frac{10000}{2.361}</math>:
 
<center><math>V = 900\pi \left[ 150 + \frac{375000}{\frac{360000}{85}} \right] = 900\pi \left[ 150 + \frac{375000 \cdot 85}{360000} \right]</math></center>
 
<center><math>V = 900\pi \left[ 150 + \frac{31875000}{360000} \right] = 900\pi \left[ 150 + 88.54167 \right]</math></center>
 
<center><math>V \approx 900\pi (238.54167) \approx 674185.87 \text{ m}^3</math></center>
 
Por tanto, suponiendo que la masa de aire estuviera toda ella a 15ºC, esta tendría la siguiente masa:
 
<center><math>M_{\text{15ºC}} = V \cdot \rho_{\text{aire}} = 674185.87 \text{ m}^3 \cdot 1.225 \text{ kg/m}^3</math></center>
 
<center><math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math></center>
 
Comparando ambos valores: <math>M_{\text{real}} \approx 781420 \text{ kg}</math> y <math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math> se observa que el aire tiene una masa aproximadamente un '''5% menor''' que la masa que tendría si estuviera a temperatura constante <math>T=15ºC</math>.
 

===Cálculo de la temperatura media y la potencia térmica evacuada===
La temperatura media se define como:
 
<center><math>T_{\text{media}} = \frac{1}{V} \iiint_V T(\rho, z) dV = \frac{1}{V} \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math></center>
 
Se definirá <math>I_T=\iiint_V T(\rho, z) dV = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math>
Sustituyendo y definiendo los límites de integración se llega a la siguiente integral:
 
<center><math>I_T = 2\pi \int_0^{150} \int_0^{\rho(z)} \left[ 57 - 38 \left(\frac{z}{150}\right)^{1.8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}} \right] \rho \, d\rho dz</math></center>
 
Debido al término exponencial esta integral no es de fácil resolución, luego utilizando una herramienta de internet se llega al siguiente resultado:
 
<center><math>I_T \approx 19534420 \text{ m}^3 \cdot \text{°C}</math></center>
 
Luego la temperatura media será:
 
<center><math>T_{\text{media}} = \frac{I_T}{V} = \frac{19534420 \text{ m}^3 \cdot \text{°C}}{674185.87 \text{ m}^3}</math></center>
 
<center><math>T_{\text{media}} \approx 28.97 \text{ °C}</math></center>
 
La potencia térmica evacuada puede calcularse mediante la siguiente expresión:
 
<center><math>P = \rho_{\text{aire}} c_p Q (T_{\text{media}} - T_{\text{ambiente}})</math></center>
 
Siendo:
*<math>\rho_{\text{aire}} = 1.225 \text{ kg/m}^3</math> la densidad del aire a temperatura <math>T_{\text{ambiente}}=15ºC</math>;
*<math>c_p = 1005 \text{ J/(kg K)}</math> el calor específico a presión constante del aire y <math>T_{\text{ambiente}}=15ºC</math>;
*<math>Q = 1500 \text{ m}^3/\text{s}</math> el caudal de aire;
*<math>T_{\text{media}} = 28.97 \text{ °C}</math> la temperatura media calculada anteriormente;
*<math>T_{\text{ambiente}} = 15 \text{ °C}</math> la temperatura ambiente.
 
Sustituyendo:
 
<center><math>P = 1.225 \cdot 1005 \cdot 1500 \cdot (28.97 - 15)</math></center>
 
<center><math>P = 1845600 \cdot 13.97</math></center>
 
<center><math>P \approx 25781000 \text{ W}=25.781 \text{ MW}</math></center>
 

== Análisis de la 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:
 
<center><math>V(z)=V_0·(\dfrac{z}{Zref})^α</math></center>
 
 <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>.

===.-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:
 
<center><math>P(z)=\dfrac{1}{2}ρ·V(z)^2</math></center>
 
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:
 
<center><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></center>
 
Como vemos analíticamente y veremos a continuación gráficamente, la presión es máxima cuando la altura es máxima, es decir a 150 metros, en la parte superior de la torre de enfriamiento, siendo esta de:
 
<center><math>P(z=150)=\dfrac{1}{2}ρ·V(z=150)^2=\dfrac{1}{2}ρ·(V_0·(\dfrac{z}{Zref})^α)^2=\dfrac{1,225}{2}·(18·(\dfrac{z}{10})^{0,1429})^2≈102.79·150^{0,286}=430.837Pa</math></center>
 

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:

[[Archivo:Captura de pantalla 2025-12-01 125000.png|600px|miniaturadeimagen|derecha|Campo escalar de presiones sobre la mitad de la torre expuesta al viento]]

<syntaxhighlight lang="matlab">
%determinación de los parámetros

a = 30; z0 = 100; c = 65,07914;
V0 = 18;
zref = 10;
alpha = 0.1429;
rho_air = 1.225; %densidad en [kg/m3], (dato extraido de worldpress.com)

% Determinación de parámetros para hacer un mallado

v = linspace(0, 150, 100);
u = linspace(0, pi, 100);

%mallado de la superficie a través de los parámetros anteriores y cálculos

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

%representacíon de la superficie

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)');
</syntaxhighlight>

=== 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:
 
<center><math>\vec{F}(x, y, z)=-P(z)·\vec{n}</math></center>
 
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:

[[Archivo:Captura_de_pantalla_2025-12-01_124522.png|600px|miniaturadeimagen|derecha|Campo de fuerzas sobre la mitad de la torre de enfriamiento]]

<syntaxhighlight lang="matlab">
%determinacíon de constantes

a = 30; c = 65,07914; z0 = 100;
H = 150;zref=10;

%determinacíon de parametros

theta = linspace(0, pi, 18);
z = linspace(0, H, 18);

%mallado con los parámetros anteriores

[Theta, Z] = meshgrid(theta, z);
R = a * sqrt(1 + ((Z - z0).^2 / c^2));

%pasar de cordenadas cilindricas a cartesianas

X = R .* cos(Theta);
Y = R .* sin(Theta);
P_z = 102.79 * Z.^0.286

%derivadas parciales

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

% obtención del 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;
Norm = sqrt(NX.^2 + NY.^2 + NZ.^2);

% Paso a unitarios 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;
</syntaxhighlight>

=== Fuerza total y por unidad de superficie del viento sobre la mitad expuesta ===
Teniendo en cuenta la presión dinámica total calculada previamente P(z), se calcula la fuerza total sobre la mitad frontal partiendo de la expresión (integrando en θ la mitad frontal da un factor 2):
 
<center><math>F = 2 \int_0^H p(z)\, R(z)\, dz= 2 \int_0^H \frac{1}{2}\rho\, V(z)^2\, R(z)\, dz= \rho \int_0^H V(z)^2 R(z)\, dz .</math></center>
 
Sustituyendo el perfil de viento y el radio lineal:
 
<center><math>V(z)^2 = V_0^{\,2} \left( \frac{z}{Z_{\mathrm{ref}}} \right)^{2\alpha}, R(z) = R_{\min} + \frac{\Delta R}{H} z, </math></center>
 
Queda
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\int_0^H z^{2\alpha} \left( R_{\min} + \frac{\Delta R}{H}z \right)\, dz.</math></center>
 
Separando integrales:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min} \int_0^H z^{2\alpha}\, dz+\frac{\Delta R}{H} \int_0^H z^{2\alpha+1}\, dz\right].</math></center>
 
Las integrales elementales resueltas son:
 
<center><math>\int_0^H z^{2\alpha}\, dz= \frac{H^{2\alpha+1}}{2\alpha+1},\qquad\int_0^H z^{2\alpha+1}\, dz= \frac{H^{2\alpha+2}}{2\alpha+2}.</math></center>
 
Por tanto:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min}\,\frac{H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R}{H}\,\frac{H^{2\alpha+2}}{2\alpha+2}\right].</math></center>
 
Finalmente
<center><math>\boxed{F =\rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left(\frac{R_{\min}\, H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R\, H^{2\alpha+1}}{2\alpha+2}\right)}</math></center>
 
Los datos utilizados son:
 
*<math>H = 150\ \text{m}</math>,altura de la chimenea
*<math>R_{\min} = 30\ \text{m}, \qquad R_{\max} = 55\ \text{m}</math> radio minimo y máximo
*<math>R(z) = R_{\min} + \frac{\Delta R}{H} z</math> ley de variación lineal del radio
*<math>V_0 = 18\ \text{m/s}, \qquad Z_{\text{ref}} = 10\ \text{m}, \qquad\alpha = 0.1429</math> valores numéricos del perfil de viento
*<math>\rho = 1.225\ \text{kg/m}^3</math> densidad del aire
*<math>p(z) = \frac{1}{2}\,\rho\, V(z)^2</math> presión dinámica del viento
 
Calculando numéricamente la expresión se obtiene:
 
<center><math>F \approx 4.423835 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
<center><math>F = 4.423835 \times 10^{6}\ \text{N}.</math></center>
 
Es necesario calcular el área de la mitad de la torre de enfriamiento, la superficie lateral de revolución correspondiente a la mitad frontal es:
 
<center><math>A_{\text{half}}=\pi \int_0^{H} R(z)\, \sqrt{1 + [R'(z)]^{2}}\, dz.</math></center>
Dado que
<center><math>R'(z)=\frac{\Delta R}{H},</math></center>
es constante, el factor geométrico puede sacarse fuera de la integral:
 
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\int_0^{H} R(z)\, dz.</math></center>
 
La integral del radio es:
 
<center><math>\int_0^{H} R(z)\, dz=\int_0^{H} \left( R_{\min} + \frac{\Delta R}{H} z \right) dz=R_{\min} H + \frac{\Delta R}{2}H.</math></center>
 
Por tanto:
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\left[R_{\min} H+\frac{\Delta R}{2}H\right].</math></center>
 
<center><math>\boxed{A_{\text{half}}=\pi H \sqrt{1 + \left( \frac{\Delta R}{H} \right)^2}\left( R_{\min} + \frac{\Delta R}{2} \right)}.</math></center>
 
<center><math>A_{\text{half}} \approx 2.0303909708635903 \times 10^{4}\ \text{m}^2,\qquad(\approx 2.03\times 10^{4}\ \text{m}^2).</math></center>
 
La fuerza por unidad de superficie es:
 
<center><math>\boxed{p_{\text{avg}} = \frac{F}{A_{\text{half}}}}.</math></center>
 
Con los valores obtenidos:
 
<center><math>p_{\text{avg}} \approx 217.881\ \text{Pa}.</math></center>

===Comparación con una torre cilíndrica===
Torre cilíndrica de radio 𝑅 = 𝑅max = 55 m y misma altura H = 150m.
 
Realizando los mismos cálculos que en el apartado anterior se compara la forma hiperbólica de una torre de enfriamiento con una cilíndrica.
Pero hay que tener en cuenta que es un cilindro, por lo tanto:
 
<center><math>R(z) = R = \text{cte}</math></center> El radio es constante
 
Por lo que las fórmula que se deben utilizar son:
 
<center><math>\boxed{F_{\text{cil}} =\rho\,R\,V_0^{2}\,Z_{\text{ref}}^{-2\alpha}\,\frac{H^{\,2\alpha+1}}{2\alpha+1}}</math></center>
 
<center><math>\boxed{A_{\text{cil}} = \pi R H}</math></center>
 
<center><math>\boxed{p(z) = \frac{1}{2}\,\rho\,V_0^{2}\left( \frac{z}{Z_{\text{ref}}} \right)^{2\alpha}}</math></center>
 
Fuerza total sobre la mitad frontal
 
Tronco hiperbólico:
 
<center><math>F_{\text{hiperbólico}} \approx 4.4238 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
Cilindro:
 
<center><math>F_{\text{cil}} \approx 5.5219 \times 10^{6}\ \text{N}\quad (\approx 5.52\ \text{MN}).</math></center>
 
Relación:
 
<center><math>\frac{F_{\text{cil}}}{F_{\text{hiperbólico}}} \approx 1.248,</math></center>
 
es decir, el cilindro recibe aproximadamente un \(24.8\%\) más de fuerza total.
 
 
Área de la mitad expuesta
 
 
Tronco hiperbólico:
 
<center><math>A_{\text{half, hiperbólico}} \approx 2.0304 \times 10^{4}\ \text{m}^2.</math></center>
 
Cilindro:
 
<center><math>A_{\text{half, cil}} = \pi R H \approx 2.5918 \times 10^{4}\ \text{m}^2.</math></center>
 
Relación:
 
<center><math>\frac{A_{\text{half, cil}}}{A_{\text{half, hiperbólico}}}\approx 1.277,</math></center>
 
por lo que el cilindro tiene aproximadamente un \(27.7\%\) más de área expuesta, la misma relación indica que el cilindro requiere cerca de un \(27.7\%\) más de material superficial.
 
 
Fuerza por unidad de superficie (presión media)
 
 
Tronco hiperbólico:
 
<center><math>p_{\text{avg, hiperbólico}}= \frac{F_{\text{hiperbólico}}}{A_{\text{half, hiperbólico}}}\approx 217.88\ \text{Pa}.</math></center>
 
Cilindro:
 
<center><math>p_{\text{avg, cil}}= \frac{F_{\text{cil}}}{A_{\text{half, cil}}}\approx 213.05\ \text{Pa}.</math></center>
 
Aunque el cilindro soporta mayor fuerza total, su presión media es ligeramente menor, debido a que tiene una superficie expuesta significativamente mayor y la carga se reparte sobre un área más grande.
 
En conclusión, la torre de forma hiperbólica soporta menos fuerza lateral y usa menos material que si fuese de forma cilíndrica, aunque esta soporta menos presión media debido a sus mayor superficie. Esto muestra que el uso de una forma hiperbólica no es cuestión estética sino de mayor eficiencia en resistencia y en coste.

== Aplicaciones practicas del hiperboloide de revolución ==
[[Archivo:Cofrentes.png|250px||miniaturadeimagen|derecha|Torres de enfriamiento de la central nuclear de Cofrentes ]]
[[Archivo:Didcot-power-station-5196798.jpg|miniaturadeimagen|centro|Torres de enfriamiento hiperbólicas de la central térmica de Didcot ]]
[[Archivo:Imagentorres.jpg ‎|miniaturadeimagen|derecha|Torres de enfriamiento de la central termoeléctrica de Niederaussem (Alemania).]]
 
Tras la investigación hecha anteriormente sobre las ventajas que tiene este tipo de geometría, estructurales y económicas, se presentaran algunos ejemplos del uso de esta geometría en las torres de enfriamiento y de la aplicación en otras estructuras mas del ámbito arquitectónico.
 
===Ejemplo de torres de enfriamiento===

 
Las torres de enfriamiento de las centrales energéticas son la máxima expresión del hiperboloide de revolución en la construcción. Ejemplos muy representativos son la central nuclear de Cofrentes en España, Didcot en Reino Unido y Niederaussem en Alemania.
 
==== Central nuclear de Cofrentes ====
 
Esta central se sitúa en la provincia de Valencia. Cuenta con dos torres de enfriamiento de 129 metros de altura. Tiene un diámetro de base de 90 metros y gracias a su funcionamiento ayudan a reducir el consumo de agua del río Júcar.
 
==== Central térmica de Didcot ====
 
La central térmica de Didcot se encuentra en Reino Unido y esta compuesta de dos partes, Didoct A y Didosct B . La central es famosa por sus seis torres de enfriamiento hiperbólicas distribuidas en dos grupos de tres localizadas en la parte A. Esta parte ceso su funcionamiento en 2013 por no cumplir con los estatutos de la Unión Europea sobre grandes instalaciones de combustión. En 2014 empezó su demolición hasta 2019 que culmino con la destrucción de las tres torres del norte. Estas estructuras median 114 metros y eran parte de Oxfordshire.
 
====Central termoeléctrica de Niederaussem ====
 
Esta central se encuentra en Alemania y alberga una de las mayores torres de enfriamiento hiperbólico en el mundo. Tuvo el récord de la torre hiperbólica mas grande desde 2002 hasta 2012, cuando le fue arrebatado por la torre de la central térmica de Kalisindh, situada en Rajastán (India). La torre más importante de Niederaussem tiene una altura de <math>200 m</math> y cuenta con un diámetro en la base de <math>152 m</math>. A pesar de su altura el grosor del hormigón en algunas partes es de tan <math>22 cm</math>.

=== Ejemplos de estructuras hiperboloides en la arquitectura e ingeniera ===
[[Archivo:11Photo0jpg.jpg|miniaturadeimagen|derecha|Planetario James S. McDonnell (San Luis, EE. UU.)]]
[[Archivo:61jl35o0nyt61.jpg|miniaturadeimagen|derecha|Restaurante los Manantiales (Xochimilco, México)]]
 
Las torres de enfriamiento no es la única estructura que emplea geometría hiperbólicas. Esta geometría es empleada en todo tipo de construcciones como torres, puentes, tejados e incluso fachadas. El planetario James S. McDonnell (San Luis, EE. UU.) y el el restaurante los Manantiales (Xochimilco, México) son dos ejemplos del uso de la geometría hiperbólica.
 
==== El planetario James S. McDonnell (San Luis, EE. UU.) ====
 
Este edificio tiene estructura de hiperboloide de revolución, es la misma estructura que las torres de enfriamiento pero empleada para un edificio cultural. Fue construido en 1963 por el arquitecto Gyo Obata. Este eligió esta forma para buscar un edificio que conectara la tierra con el espacio, lo que simboliza el planetario. El planetario es un hiperboloide de una hoja, gracias a esta forma las paredes son delgadas, reduciendo el material empleado.
 
==== Restaurante Los Manantiales (Xochimilco, México) ====
==== Restaurante Los Manantiales (Xochimilco, México) ====

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-03T17:44:36Z

Ramón Castán Naval: /* .-Representación del campo escalar de presiones. */

{{ TrabajoED |Torres de enfriamiento hiperbólicas (Grupo 06) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] |
*Álvaro Díaz Martín
*Ramón Castán Naval
*Manuel Álvarez-Campana Illescas
*Gonzalo Quemada Simón
*Fernando Barettino Moreno-Aurioles }}
[[Categoría:Teoría de Campos]][[Categoría:TC25/26]] 

== 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:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</math> unos valores a determinar
Sean en nuestro modelo <math>Rmáx=55m</math> ; <math>Rmín=30m</math> ; <math>H=150m</math>

== 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=150m</math>). Al ser la figura un solido de revolución, generado al girar una hipérbola alrededor del eje Z , tiene simetría radial.

=== 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>.

<center><math>z_0=\dfrac{\scriptsize 2}{\scriptsize 3}H=\dfrac{\scriptsize 2}{\scriptsize 3}150=100</math></center>
 
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>.

<center><math>a=30</math></center>
 
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>
 
<center><math>\dfrac{x^2+y^2}{30^2}-\dfrac{(z-100)^2}{c^2} = 1</math></center>
 
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.
 
<center><math>\dfrac{55^2+0^2}{30^2}-\dfrac{(0-100)^2}{c^2} = 1</math></center>
 
Simplificando la ecuación
 
<center><math>\dfrac{55^2}{30^2}-\dfrac{(-100)^2}{c^2} = 1</math></center>
 
Elevando al cuadrado:
<center><math>\dfrac{3025}{900}-\dfrac{10000}{c^2} = 1</math></center>
 
Despejando <math>c</math>:
 
<center><math>3,361-\dfrac{10000}{c^2} = 1</math> <math>\longrightarrow</math> <math>3,361-1 = \dfrac{10000}{c^2}</math></center>
 
<center><math>2,361= \dfrac{10000}{c^2}</math> <math>\longrightarrow</math> <math>c^2=\dfrac{10000}{2,361}</math></center>
 
<center><math>c=\sqrt{\dfrac{10000}{2,361}}</math> <math>\longrightarrow</math> <math>c=65,07914</math></center>
 
Los parámetros obtenidos son <math>a=30,z_0=100,c=65,07914</math>

===Ecuación de la superficie hiperbólica en coordenadas cilíndricas===
Para obtener la ecuación en coordenadas cilíndricas se parte de la expresión anterior en coordenadas cartesianas:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</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 de las coordenadas cartesianas <math>(x,y,z)</math> en coordenadas cilíndricas <math>(ρ,θ,z).</math>
 
<center><math>x=ρ\cdot\cos(θ)</math> 
<math>y=ρ\cdot\sin(θ)</math> 
<math>z=z</math></center>
 
Cambiando <math>x^2+y^2</math> a coordenadas cilíndricas:
 
<center><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></center>
 
Se saca factor común <math>ρ^2</math>:
 
<center><math>x^2+y^2=ρ^2\cdot(\cos^2(θ)+\sin^2(θ))=1</math></center>
 
Aplicando la siguiente identidad trigonométrica:
 
<center><math>\cos^2(θ)+\sin^2(θ)=1</math></center>
 
Sustituyendo en la ecuación
 
<center><math>x^2+y^2=ρ^2\cdot 1 \longrightarrow x^2+y^2=ρ^2</math></center>
 
Al ser los demás datos constantes, no hace falta cambiar nada más. La ecuación en coordenadas cilíndricas queda: 
<center><math>\dfrac{ρ^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
===Representación grafica de la torre hiperbólica ===
 
A partir de la ecuación anterior se modela la torre de enfriamiento estándar de altura <math>(H=150m)</math>, radio máximo <math>(Rmáx=55m)</math> y radio mínimo <math>(Rmín=30m)</math>

 
[[Archivo:Torre11.jpg|miniaturadeimagen|500px|derecha|Representación grafica de la torre de enfriamiento estándar ]]
<syntaxhighlight lang="matlab">
clc
clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
u=linspace(0,2*pi,80);
v=linspace(0,H,80);

[Mt,Mz]=meshgrid(u,v);
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

figure
Mx=Mr.*cos(Mt)
My=Mr.*sin(Mt);
Mz=Mz;

surf(Mx,My,Mz, 'FaceColor','white','EdgeColor','black')
axis equal
hold on
grid on
title('Torre de enfriamiento hiperbólica')
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
</syntaxhighlight>

== 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: 
 
<center><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></center>
 
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: 
 
<center><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></center>
 
<center><math>T(\rho, z) = 65 - 38 \left(\frac{z}{150}\right)^{1,8} - 8 + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
<center><math>T(\rho, z) = 57 - 38 \left(\frac{z}{150}\right)^{1,8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
===Representación del campo escalar 𝑇(𝜌,𝑧) como un mapa de colores en una sección vertical y varias secciones transversales===
Se procede a representar el campo escalar de temperaturas como un mapa de colores con una interfaz que simula la de una cámara térmica.
 
 
'''Interpretación'''
 
Como se aprecia, la temperatura más elevada se da en el origen de coordenadas cartesianas <math>(X,Y,Z)=(0,0,0)</math>.
 
* La temperatura en el interior de la torre disminuye según aumenta la coordenada z.
 
* A su vez, la temperatura disminuye según aumenta el alejamiento al eje z, es decir, según aumenta la coordenada ρ (en coordenadas cilíndricas).
 
Por tanto, se aprecia que la variación en las temperaturas es radial y en altura. Este resultado es coherente, teniendo en cuenta que la función escalar de temperatura depende únicamente de los parámetros en coordenadas cilíndricas: z y ρ (y no de θ).
 
 
[[Archivo:grupo_6_campo_temperaturas_vertical.jpg|500px|Representación grafica del campo escalar de temperatura (sección vertical)]]
[[Archivo:grupo_6_campo_temperaturas_transversal.jpg|500px|Representación grafica del campo escalar de temperatura (secciones transversales)]]
 
A continuación se muestra el código en matlab que muestra las dos ventanas gráficas:
<syntaxhighlight lang="matlab">
clc; clear
%se establecen los parámetros de la torre:

H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:

Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIÓN VERTICAL-
z = linspace(0,H,300); %Se establecen los puntos del plano de sección XOZ
x = linspace(-Rmax,Rmax,300);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=Mrho; %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se representa el campo escalar Temperatura interior del plano XOZ de la
%torre
figure('name','Campo de Temperatura en sección XOZ de la torre')
contourf(Mx,Mz,T,60,'LineColor','none');
colormap('hot');
colorbar;
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo escalar de temperaturas en la sección vertical XOZ de la torre:');
tit.Position=[0,156,0]; %se edita la posición del titulo
hold on;

%Se representa la silueta de la torre
contorno_izq = -Mrho; %Lado izquierdo de la silueta
contorno_dch = Mrho; %Lado derecho de la silueta
contorno_alt = z; %Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); %Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); %Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2);%Vista en 2D
xlim([-Rmax-15,Rmax+15]);%se ajusta la ventana gráfica
ylim([0,H]);
grid off;

%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='Temp. (ºC)';
hold off;

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIONES TRANSVERSALES-
figure('name','Campo de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
rango=linspace(-Rmax,Rmax,200);
[Mx, My] = meshgrid(rango,rango);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);
for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre.
Zi=(H/3)*i; %Cuatro cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));
% T(rho, z)=T=Tbase-dTz.*(Zi./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2)))

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>Rho_z)=NaN;

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1)
pcolor(Mx, My, Ti);
shading interp; %Se suavizar el color

colormap('hot');
title(['Temperatura (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2);%Vista en 2D
xlim([-70,70]);%se ajusta la ventana gráfica
ylim([-70,70]);

hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='T (ºC)';
colorbar;
end
hold off;
</syntaxhighlight>

=== Cálculo y representación del gradiente del campo escalar Temperatura ===

Con el fin de interpretar la variación de temperaturas dentro de la torre de enfriamiento hiperbólica en cada punto, se estudia el campo vectorial gradiente del campo escalar temperatura.
 
Este campo gradiente tiene la siguiente expresión:
<center><math>∇T(r,θ,z)=(\dfrac{dT}{dρ},\dfrac{1}{ρ}·\dfrac{dT}{dθ},\dfrac{dT}{dz})</math></center>
 
Como el campo escalar temperatura no depende de <math>θ</math>,el gradiente de temperatura tiene solo dos componentes:
 
<center>Componente radial: <math>\dfrac{dT}{dρ}</math></center>
 
<center>Componente vertical: <math>\dfrac{dT}{dz}</math></center>
 
De este modo, se calculan las derivadas parciales respecto a <math>ρ</math> y respecto a <math>z</math>:
 
<center><math>\dfrac{dT}{dρ}=-ΔT_r·\dfrac{-2ρ·Rmax^2}{(Rmax^2-ρ^2)^2}·e^{(\dfrac{-ρ^2}{Rmax^2-ρ^2})}</math></center>
 
<center><math>\dfrac{dT}{dz}=-ΔT_z·n·\dfrac{z^{n-1}}{H^n}</math></center>
 
 
A continuación se muestra la representación del campo vectorial gradiente del campo Temperaturas <math>T(ρ,z):</math>
 
[[Archivo:grupo_6_gradiente_vertical.jpg|500px|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_gradiente_transversal.jpg|500px|Representación grafica del campo escalar de temperatura ]]
 
 
'''Interpretación gráfica y analítica'''
 
* El vector gradiente en cada punto apunta hacia donde la temperatura aumenta más rápidamente (el vector gradiente en cada punto sigue la trayectoria más empinada de ascenso de temperatura). Es por eso que todas los vectores apuntan hacia abajo (coordenada z=0) y hacia el eje de revolución torre (coordenadas x,y=0,0), ya que, como se ha expuesto en el apartado anterior, la temperatura en la torre aumenta según se acerca a la base y al centro de la torre, y por eso todos los vectores gradiente apuntan hacia allí.
 
 
* Otra interpretación importante es que los vectores tienen mayor módulo en los puntos más alejados de la base de la torre. Esto es así debido según se aumenta la altura, hay mayor diferencia de temperatura para la misma distancia. En la cima de la torre el aire caliente entra en contacto con el aire atmosférico, por tanto, cuanto más cerca se esté de este punto, mayor será el enfriamiento. Análogamente, según se acerca a las paredes de la torre, el gradiente tiene mayor módulo, puesto que para un mismo recorrido radial, la variación de temperatura es mayor.
 
A continuación se muestra el código en 'matlab' que genera sendas ventanas gráficas:
<syntaxhighlight lang="matlab">
clc, clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se establecen los puntos del plano de sección XOZ
z = linspace(0,H,25);
x = linspace(-Rmax,Rmax,25);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=(Mrho+4.*ones(25)); %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se calcula el campo gradiente de T
[gradx,gradz]=gradient(T,x,z);

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN VERTICAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOZ')
quiver(x,z,gradx,gradz,'r','LineWidth',0.7)
hold on

%Se representa la silueta de la torre
contorno_izq = -Mrho; % Lado izquierdo de la silueta
contorno_dch = Mrho; % Lado derecho de la silueta
contorno_alt = z; % Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); % Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); % Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2); %Vista en 2D
xlim([-Rmax-15,Rmax+15]); %se ajusta la ventana gráfica
ylim([0,H]);
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo vectorial gradiente de temperatura (Sección Vertical):');
tit.Position=[0,153,0]; %se edita la posición del titulo
hold off

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN TRANSVERSAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
x=linspace(-Rmax-10,Rmax+10,20);
y=x;
[Mx,My] = meshgrid(x,y);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);

for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre
Zi=(H/3)*i; %4 cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>(Rho_z))=NaN;

%Se calcula el gradiente
[gradx,grady]=gradient(Ti,x,y);

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1);
quiver(Mx,My,gradx,grady,'r','LineWidth',0.7);

title(['Gradiente(T) (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2); %Vista en 2D
xlim([-70,70]); %se ajusta la ventana gráfica
ylim([-70,70]);
hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
end
hold off;
</syntaxhighlight>

=== Representación de superficies isotérmicas ===
Las curvas de nivel son superficies generadas a partir de la asignación de un valor constante a un campo escalar (en <math>ℝ^3</math>). En este caso, el campo escalar es la Temperatura, T(ρ,z), y sus curvas de nivel son las 'Superficies Isotermas', obtenidas al asignar a T un valor de temperatura constante.
 
A continuación se representan las superficies isotérmicas para los valores de 30,40,50 y 60 ºC.
 
[[Archivo:grupo_6_isotermas_vertical.jpg|500px|izquierda|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_isotermas_xyz.jpg|500px|derecho|Representación grafica del campo escalar de temperatura ]]
 
'''Interpretación de las superficies de nivel:'''
 
* Las superficies isotermas tienen una forma arqueada hacia abajo, es decir, que son convexas hacia el eje z. Esto se debe a que el campo escalar temperatura (<math>T(ρ,z)</math>), como se ha comentado anteriormente, varía radial y verticalmente, disminuyendo la temperatura según se alejan los puntos del eje z (de revolución) y de la base de la torre. De este modo, para una misma altura (fijando z), la temperatura será mayor en el centro de la sección transversal que en los extremos, como bien se puede apreciar en la segunda figura del apartado 3.1.
* Es importante destacar que, según se asciende en la torre, las curvas isotermas pasan de asemejarse a una esfera a asemejarse a un plano horizontal, gradualmente. La explicación es debida a que el campo escalar Temperatura depende de dos términos, que dependen respectivamente de z y ρ.
 
El término que depende de z es:
 
<center><math>T(z) = -\Delta T_z \left(\frac{z}{H}\right)^n</math></center>
El término que depende de ρ es:
<center><math>T(\rho) = - \Delta T_{\rho} \left(1 - e^{- \frac{\rho^2}{R_{\text{max}}^2 - \rho^2}}\right)</math></center>
 
Por tanto:
* En la base de la torre (z cercano a 0), la superficie isoterma se asemeja a una esfera. Esto se debe a que el término <math>T(z)</math>, es mínimo, por lo que <math>T(ρ)</math> domina el campo escalar. Así, este último término da a las superficies isotermas la forma esférica, ya que el radio (ρ) produce el efecto más importante sobre la temperatura.
* Análogamente, en la cima de la torre (z cercano a H), la superficie se puede asemejar a un plano horizontal, puesto que el término <math>T(z)</math> es máximo. Así, el efecto que produce el radio (ρ), es menor.
 
A continuación se muestra el programa de matlab que representa en una ventana gráfica las superficies isotermas.
<syntaxhighlight lang="matlab">
clc,clear;
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;
%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38; %variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8; %variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se crean los puntos x,y,z
x = linspace(-Rmax,Rmax,200);
y = linspace(-Rmax,Rmax,200);
z = linspace(0,H,200);

[Mx,My,Mz]=meshgrid(x,y,z); %se crea la malla de puntos
Rho=sqrt(Mx.^2+My.^2); %Coordenada rho (cilíndricas)

%se calcula el campo escalar en la extensión de Rho, Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se calcula el radio de la torre para cada z
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

%Se eliminan los puntos que no están DENTRO de la torre (rho<=R(z))
Mask_interior=(Rho<=Mr);

%Se aplica la máscara: asignamos NaN a las temperaturas fuera de la torre.
T(~Mask_interior) = NaN;

%Se representa
figure('Name', 'Superficies Isotérmicas dentro de la Torre Hiperbólica');
hold on;

%Se define el colormap para los colores
clim([Tcima-5,Tbase]);
cmap = hot(45);
colormap(cmap);
cbar.Direction='reverse';

%% BUCLE
Temperaturas=[30,40,50,60]; %Se definen las temperaturas de las superficies
for i=1:4;
T_nivel=Temperaturas(i);
view(3); % Vista 3D
colorbar
%Se representa la superficie isoterma
[caras, vertices] = isosurface(Mx, My, Mz, T, T_nivel);
curva_nivel=patch('faces', caras, 'Vertices', vertices, 'FaceColor', cmap(T_nivel-23,:), 'EdgeColor', 'none');

%Se representa el contorno de la torre (con transparencia)
contorno = Rho - Mr;
[caras_t, vertices_t] = isosurface(Mx, My, Mz, contorno, 0);
h_tower = patch('faces', caras_t, 'vertices', vertices_t, 'FaceColor', [0.5 0.5 0.5], ...
'EdgeColor', 'none', 'FaceAlpha', 0.1); % Gris muy transparente

%Se ajusta de la ventana gráfica
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
axis equal tight;
end
%vista final
hold off
title('Superficies Isotérmicas dentro de la Torre');
lighting gouraud
camlight;
</syntaxhighlight>

=== Cálculo de la masa del aire contenido en la torre ===
La densidad del aire 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:
 
<center><math>\rho(T) = \rho_{\text{aire}} \frac{T_0 + 273}{T + 273}</math></center>
 
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:
 
<center><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></center>
 
Si se desarrolla la expresión se llega a lo siguiente:
 
<center><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></center>
 
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:
<center><math> M_{\text{real}} \approx 781420 \text{ kg}</math></center>
 
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:
 
<center><math>V = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} \rho \, d\rho dz d\theta</math></center>
En la sección 2.2 se obtiene la expresión del hiperboloide de revolución en coordenadas cilíndricas:
 
<center><math>\frac{\rho^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1</math></center>
Despejando se obtiene una relación entre <math>\rho</math> y <math>z</math> 
 
<center><math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math></center>
Volviendo a la integral, integrando primero <math>\theta</math> y después <math>\rho</math> se obtiene lo siguiente:
 
<center><math>V = 2\pi \int_0^{150} \left[ \frac{\rho^2}{2} \right]_0^{\rho(z)} dz = \pi \int_0^{150} \rho^2(z) dz</math></center>
 
Sustituyendo <math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math> con <math>z_0=100</math>:
 
<center><math>V = \pi a^2 \int_0^{150} \left[ 1 + \frac{(z - 100)^2}{c^2} \right] dz</math>
 
<math>V = \pi a^2 \left[ \int_0^{150} dz + \frac{1}{c^2} \int_0^{150} (z - 100)^2 dz \right]</math></center>
 
Resolviendo ambas integrales:
 
<center><math>\int_0^{150} dz = [z]_0^{150} = 150</math></center>
 
Para la realización de <math>\int_0^{150} (z - 100)^2 dz</math> usamos el cambio de variable <math>u = z - 100, du=dz</math> con límites <math>u(0)=-100, u(150)=50</math>
 
<center><math>\int_{-100}^{50} u^2 du = \left[ \frac{u^3}{3} \right]_{-100}^{50} = \frac{1}{3} \left[ 50^3 - (-100)^3 \right] = \frac{1}{3} [125000 + 1000000] = \frac{1125000}{3} = 375000</math></center>
 
Sustituyendo en la integral, con <math>a=30</math> y <math>c^2 =\frac{10000}{2.361}</math>:
 
<center><math>V = 900\pi \left[ 150 + \frac{375000}{\frac{360000}{85}} \right] = 900\pi \left[ 150 + \frac{375000 \cdot 85}{360000} \right]</math></center>
 
<center><math>V = 900\pi \left[ 150 + \frac{31875000}{360000} \right] = 900\pi \left[ 150 + 88.54167 \right]</math></center>
 
<center><math>V \approx 900\pi (238.54167) \approx 674185.87 \text{ m}^3</math></center>
 
Por tanto, suponiendo que la masa de aire estuviera toda ella a 15ºC, esta tendría la siguiente masa:
 
<center><math>M_{\text{15ºC}} = V \cdot \rho_{\text{aire}} = 674185.87 \text{ m}^3 \cdot 1.225 \text{ kg/m}^3</math></center>
 
<center><math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math></center>
 
Comparando ambos valores: <math>M_{\text{real}} \approx 781420 \text{ kg}</math> y <math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math> se observa que el aire tiene una masa aproximadamente un '''5% menor''' que la masa que tendría si estuviera a temperatura constante <math>T=15ºC</math>.
 

===Cálculo de la temperatura media y la potencia térmica evacuada===
La temperatura media se define como:
 
<center><math>T_{\text{media}} = \frac{1}{V} \iiint_V T(\rho, z) dV = \frac{1}{V} \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math></center>
 
Se definirá <math>I_T=\iiint_V T(\rho, z) dV = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math>
Sustituyendo y definiendo los límites de integración se llega a la siguiente integral:
 
<center><math>I_T = 2\pi \int_0^{150} \int_0^{\rho(z)} \left[ 57 - 38 \left(\frac{z}{150}\right)^{1.8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}} \right] \rho \, d\rho dz</math></center>
 
Debido al término exponencial esta integral no es de fácil resolución, luego utilizando una herramienta de internet se llega al siguiente resultado:
 
<center><math>I_T \approx 19534420 \text{ m}^3 \cdot \text{°C}</math></center>
 
Luego la temperatura media será:
 
<center><math>T_{\text{media}} = \frac{I_T}{V} = \frac{19534420 \text{ m}^3 \cdot \text{°C}}{674185.87 \text{ m}^3}</math></center>
 
<center><math>T_{\text{media}} \approx 28.97 \text{ °C}</math></center>
 
La potencia térmica evacuada puede calcularse mediante la siguiente expresión:
 
<center><math>P = \rho_{\text{aire}} c_p Q (T_{\text{media}} - T_{\text{ambiente}})</math></center>
 
Siendo:
*<math>\rho_{\text{aire}} = 1.225 \text{ kg/m}^3</math> la densidad del aire a temperatura <math>T_{\text{ambiente}}=15ºC</math>;
*<math>c_p = 1005 \text{ J/(kg K)}</math> el calor específico a presión constante del aire y <math>T_{\text{ambiente}}=15ºC</math>;
*<math>Q = 1500 \text{ m}^3/\text{s}</math> el caudal de aire;
*<math>T_{\text{media}} = 28.97 \text{ °C}</math> la temperatura media calculada anteriormente;
*<math>T_{\text{ambiente}} = 15 \text{ °C}</math> la temperatura ambiente.
 
Sustituyendo:
 
<center><math>P = 1.225 \cdot 1005 \cdot 1500 \cdot (28.97 - 15)</math></center>
 
<center><math>P = 1845600 \cdot 13.97</math></center>
 
<center><math>P \approx 25781000 \text{ W}=25.781 \text{ MW}</math></center>
 

== Análisis de la 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:
 
<center><math>V(z)=V_0·(\dfrac{z}{Zref})^α</math></center>
 
 <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>.

===.-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:
 
<center><math>P(z)=\dfrac{1}{2}ρ·V(z)^2</math></center>
 
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:
 
<center><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></center>
 
Como vemos analíticamente y veremos a continuación gráficamente, la presión es máxima cuando la altura es máxima, es decir a 150 metros, en la parte superior de la torre de enfriamiento, siendo esta de:
 
<center><math>P(z=150)=\dfrac{1}{2}ρ·V(z=150)^2=\dfrac{1}{2}ρ·(V_0·(\dfrac{z}{Zref})^α)^2=\dfrac{1,225}{2}·(18·(\dfrac{z}{10})^{0,1429})^2≈102.79·150^{0,286}=430.837Pa</math></center>
 

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:

[[Archivo:Captura de pantalla 2025-12-01 125000.png|600px|miniaturadeimagen|derecha|Campo escalar de presiones sobre la mitad de la torre expuesta al viento]]

<syntaxhighlight lang="matlab">
%determinación de los parámetros

a = 30; z0 = 100; c = 65,07914;
V0 = 18;
zref = 10;
alpha = 0.1429;
rho_air = 1.225; %densidad en [kg/m3], (dato extraido de worldpress.com)

% Determinación de parámetros para hacer un mallado

v = linspace(0, 150, 100);
u = linspace(0, pi, 100);

%mallado de la superficie a través de los parámetros anteriores y cálculos

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

%representacíon de la superficie

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)');
</syntaxhighlight>

=== 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:
 
<center><math>\vec{F}(x, y, z)=-P(z)·\vec{n}</math></center>
 
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:

[[Archivo:Captura_de_pantalla_2025-12-01_124522.png|600px|miniaturadeimagen|derecha|Campo de fuerzas sobre la mitad de la torre de enfriamiento]]

<syntaxhighlight lang="matlab">
a = 30; c = 65,07914; z0 = 100;
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;
</syntaxhighlight>

=== Fuerza total y por unidad de superficie del viento sobre la mitad expuesta ===
Teniendo en cuenta la presión dinámica total calculada previamente P(z), se calcula la fuerza total sobre la mitad frontal partiendo de la expresión (integrando en θ la mitad frontal da un factor 2):
 
<center><math>F = 2 \int_0^H p(z)\, R(z)\, dz= 2 \int_0^H \frac{1}{2}\rho\, V(z)^2\, R(z)\, dz= \rho \int_0^H V(z)^2 R(z)\, dz .</math></center>
 
Sustituyendo el perfil de viento y el radio lineal:
 
<center><math>V(z)^2 = V_0^{\,2} \left( \frac{z}{Z_{\mathrm{ref}}} \right)^{2\alpha}, R(z) = R_{\min} + \frac{\Delta R}{H} z, </math></center>
 
Queda
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\int_0^H z^{2\alpha} \left( R_{\min} + \frac{\Delta R}{H}z \right)\, dz.</math></center>
 
Separando integrales:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min} \int_0^H z^{2\alpha}\, dz+\frac{\Delta R}{H} \int_0^H z^{2\alpha+1}\, dz\right].</math></center>
 
Las integrales elementales resueltas son:
 
<center><math>\int_0^H z^{2\alpha}\, dz= \frac{H^{2\alpha+1}}{2\alpha+1},\qquad\int_0^H z^{2\alpha+1}\, dz= \frac{H^{2\alpha+2}}{2\alpha+2}.</math></center>
 
Por tanto:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min}\,\frac{H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R}{H}\,\frac{H^{2\alpha+2}}{2\alpha+2}\right].</math></center>
 
Finalmente
<center><math>\boxed{F =\rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left(\frac{R_{\min}\, H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R\, H^{2\alpha+1}}{2\alpha+2}\right)}</math></center>
 
Los datos utilizados son:
 
*<math>H = 150\ \text{m}</math>,altura de la chimenea
*<math>R_{\min} = 30\ \text{m}, \qquad R_{\max} = 55\ \text{m}</math> radio minimo y máximo
*<math>R(z) = R_{\min} + \frac{\Delta R}{H} z</math> ley de variación lineal del radio
*<math>V_0 = 18\ \text{m/s}, \qquad Z_{\text{ref}} = 10\ \text{m}, \qquad\alpha = 0.1429</math> valores numéricos del perfil de viento
*<math>\rho = 1.225\ \text{kg/m}^3</math> densidad del aire
*<math>p(z) = \frac{1}{2}\,\rho\, V(z)^2</math> presión dinámica del viento
 
Calculando numéricamente la expresión se obtiene:
 
<center><math>F \approx 4.423835 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
<center><math>F = 4.423835 \times 10^{6}\ \text{N}.</math></center>
 
Es necesario calcular el área de la mitad de la torre de enfriamiento, la superficie lateral de revolución correspondiente a la mitad frontal es:
 
<center><math>A_{\text{half}}=\pi \int_0^{H} R(z)\, \sqrt{1 + [R'(z)]^{2}}\, dz.</math></center>
Dado que
<center><math>R'(z)=\frac{\Delta R}{H},</math></center>
es constante, el factor geométrico puede sacarse fuera de la integral:
 
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\int_0^{H} R(z)\, dz.</math></center>
 
La integral del radio es:
 
<center><math>\int_0^{H} R(z)\, dz=\int_0^{H} \left( R_{\min} + \frac{\Delta R}{H} z \right) dz=R_{\min} H + \frac{\Delta R}{2}H.</math></center>
 
Por tanto:
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\left[R_{\min} H+\frac{\Delta R}{2}H\right].</math></center>
 
<center><math>\boxed{A_{\text{half}}=\pi H \sqrt{1 + \left( \frac{\Delta R}{H} \right)^2}\left( R_{\min} + \frac{\Delta R}{2} \right)}.</math></center>
 
<center><math>A_{\text{half}} \approx 2.0303909708635903 \times 10^{4}\ \text{m}^2,\qquad(\approx 2.03\times 10^{4}\ \text{m}^2).</math></center>
 
La fuerza por unidad de superficie es:
 
<center><math>\boxed{p_{\text{avg}} = \frac{F}{A_{\text{half}}}}.</math></center>
 
Con los valores obtenidos:
 
<center><math>p_{\text{avg}} \approx 217.881\ \text{Pa}.</math></center>

===Comparación con una torre cilíndrica===
Torre cilíndrica de radio 𝑅 = 𝑅max = 55 m y misma altura H = 150m.
 
Realizando los mismos cálculos que en el apartado anterior se compara la forma hiperbólica de una torre de enfriamiento con una cilíndrica.
Pero hay que tener en cuenta que es un cilindro, por lo tanto:
 
<center><math>R(z) = R = \text{cte}</math></center> El radio es constante
 
Por lo que las fórmula que se deben utilizar son:
 
<center><math>\boxed{F_{\text{cil}} =\rho\,R\,V_0^{2}\,Z_{\text{ref}}^{-2\alpha}\,\frac{H^{\,2\alpha+1}}{2\alpha+1}}</math></center>
 
<center><math>\boxed{A_{\text{cil}} = \pi R H}</math></center>
 
<center><math>\boxed{p(z) = \frac{1}{2}\,\rho\,V_0^{2}\left( \frac{z}{Z_{\text{ref}}} \right)^{2\alpha}}</math></center>
 
Fuerza total sobre la mitad frontal
 
Tronco hiperbólico:
 
<center><math>F_{\text{hiperbólico}} \approx 4.4238 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
Cilindro:
 
<center><math>F_{\text{cil}} \approx 5.5219 \times 10^{6}\ \text{N}\quad (\approx 5.52\ \text{MN}).</math></center>
 
Relación:
 
<center><math>\frac{F_{\text{cil}}}{F_{\text{hiperbólico}}} \approx 1.248,</math></center>
 
es decir, el cilindro recibe aproximadamente un \(24.8\%\) más de fuerza total.
 
 
Área de la mitad expuesta
 
 
Tronco hiperbólico:
 
<center><math>A_{\text{half, hiperbólico}} \approx 2.0304 \times 10^{4}\ \text{m}^2.</math></center>
 
Cilindro:
 
<center><math>A_{\text{half, cil}} = \pi R H \approx 2.5918 \times 10^{4}\ \text{m}^2.</math></center>
 
Relación:
 
<center><math>\frac{A_{\text{half, cil}}}{A_{\text{half, hiperbólico}}}\approx 1.277,</math></center>
 
por lo que el cilindro tiene aproximadamente un \(27.7\%\) más de área expuesta, la misma relación indica que el cilindro requiere cerca de un \(27.7\%\) más de material superficial.
 
 
Fuerza por unidad de superficie (presión media)
 
 
Tronco hiperbólico:
 
<center><math>p_{\text{avg, hiperbólico}}= \frac{F_{\text{hiperbólico}}}{A_{\text{half, hiperbólico}}}\approx 217.88\ \text{Pa}.</math></center>
 
Cilindro:
 
<center><math>p_{\text{avg, cil}}= \frac{F_{\text{cil}}}{A_{\text{half, cil}}}\approx 213.05\ \text{Pa}.</math></center>
 
Aunque el cilindro soporta mayor fuerza total, su presión media es ligeramente menor, debido a que tiene una superficie expuesta significativamente mayor y la carga se reparte sobre un área más grande.
 
En conclusión, la torre de forma hiperbólica soporta menos fuerza lateral y usa menos material que si fuese de forma cilíndrica, aunque esta soporta menos presión media debido a sus mayor superficie. Esto muestra que el uso de una forma hiperbólica no es cuestión estética sino de mayor eficiencia en resistencia y en coste.

== Aplicaciones practicas del hiperboloide de revolución ==
[[Archivo:Cofrentes.png|250px||miniaturadeimagen|derecha|Torres de enfriamiento de la central nuclear de Cofrentes ]]
[[Archivo:Didcot-power-station-5196798.jpg|miniaturadeimagen|centro|Torres de enfriamiento hiperbólicas de la central térmica de Didcot ]]
[[Archivo:Imagentorres.jpg ‎|miniaturadeimagen|derecha|Torres de enfriamiento de la central termoeléctrica de Niederaussem (Alemania).]]
 
Tras la investigación hecha anteriormente sobre las ventajas que tiene este tipo de geometría, estructurales y económicas, se presentaran algunos ejemplos del uso de esta geometría en las torres de enfriamiento y de la aplicación en otras estructuras mas del ámbito arquitectónico.
 
===Ejemplo de torres de enfriamiento===

 
Las torres de enfriamiento de las centrales energéticas son la máxima expresión del hiperboloide de revolución en la construcción. Ejemplos muy representativos son la central nuclear de Cofrentes en España, Didcot en Reino Unido y Niederaussem en Alemania.
 
==== Central nuclear de Cofrentes ====
 
Esta central se sitúa en la provincia de Valencia. Cuenta con dos torres de enfriamiento de 129 metros de altura. Tiene un diámetro de base de 90 metros y gracias a su funcionamiento ayudan a reducir el consumo de agua del río Júcar.
 
==== Central térmica de Didcot ====
 
La central térmica de Didcot se encuentra en Reino Unido y esta compuesta de dos partes, Didoct A y Didosct B . La central es famosa por sus seis torres de enfriamiento hiperbólicas distribuidas en dos grupos de tres localizadas en la parte A. Esta parte ceso su funcionamiento en 2013 por no cumplir con los estatutos de la Unión Europea sobre grandes instalaciones de combustión. En 2014 empezó su demolición hasta 2019 que culmino con la destrucción de las tres torres del norte. Estas estructuras median 114 metros y eran parte de Oxfordshire.
 
====Central termoeléctrica de Niederaussem ====
 
Esta central se encuentra en Alemania y alberga una de las mayores torres de enfriamiento hiperbólico en el mundo. Tuvo el récord de la torre hiperbólica mas grande desde 2002 hasta 2012, cuando le fue arrebatado por la torre de la central térmica de Kalisindh, situada en Rajastán (India). La torre más importante de Niederaussem tiene una altura de <math>200 m</math> y cuenta con un diámetro en la base de <math>152 m</math>. A pesar de su altura el grosor del hormigón en algunas partes es de tan <math>22 cm</math>.

=== Ejemplos de estructuras hiperboloides en la arquitectura e ingeniera ===
[[Archivo:11Photo0jpg.jpg|miniaturadeimagen|derecha|Planetario James S. McDonnell (San Luis, EE. UU.)]]
[[Archivo:61jl35o0nyt61.jpg|miniaturadeimagen|derecha|Restaurante los Manantiales (Xochimilco, México)]]
 
Las torres de enfriamiento no es la única estructura que emplea geometría hiperbólicas. Esta geometría es empleada en todo tipo de construcciones como torres, puentes, tejados e incluso fachadas. El planetario James S. McDonnell (San Luis, EE. UU.) y el el restaurante los Manantiales (Xochimilco, México) son dos ejemplos del uso de la geometría hiperbólica.
 
==== El planetario James S. McDonnell (San Luis, EE. UU.) ====
 
Este edificio tiene estructura de hiperboloide de revolución, es la misma estructura que las torres de enfriamiento pero empleada para un edificio cultural. Fue construido en 1963 por el arquitecto Gyo Obata. Este eligió esta forma para buscar un edificio que conectara la tierra con el espacio, lo que simboliza el planetario. El planetario es un hiperboloide de una hoja, gracias a esta forma las paredes son delgadas, reduciendo el material empleado.
 
==== Restaurante Los Manantiales (Xochimilco, México) ====
==== Restaurante Los Manantiales (Xochimilco, México) ====

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-03T17:40:52Z

Ramón Castán Naval: /* .-Representación del campo escalar de presiones. */

{{ TrabajoED |Torres de enfriamiento hiperbólicas (Grupo 06) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] |
*Álvaro Díaz Martín
*Ramón Castán Naval
*Manuel Álvarez-Campana Illescas
*Gonzalo Quemada Simón
*Fernando Barettino Moreno-Aurioles }}
[[Categoría:Teoría de Campos]][[Categoría:TC25/26]] 

== 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:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</math> unos valores a determinar
Sean en nuestro modelo <math>Rmáx=55m</math> ; <math>Rmín=30m</math> ; <math>H=150m</math>

== 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=150m</math>). Al ser la figura un solido de revolución, generado al girar una hipérbola alrededor del eje Z , tiene simetría radial.

=== 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>.

<center><math>z_0=\dfrac{\scriptsize 2}{\scriptsize 3}H=\dfrac{\scriptsize 2}{\scriptsize 3}150=100</math></center>
 
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>.

<center><math>a=30</math></center>
 
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>
 
<center><math>\dfrac{x^2+y^2}{30^2}-\dfrac{(z-100)^2}{c^2} = 1</math></center>
 
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.
 
<center><math>\dfrac{55^2+0^2}{30^2}-\dfrac{(0-100)^2}{c^2} = 1</math></center>
 
Simplificando la ecuación
 
<center><math>\dfrac{55^2}{30^2}-\dfrac{(-100)^2}{c^2} = 1</math></center>
 
Elevando al cuadrado:
<center><math>\dfrac{3025}{900}-\dfrac{10000}{c^2} = 1</math></center>
 
Despejando <math>c</math>:
 
<center><math>3,361-\dfrac{10000}{c^2} = 1</math> <math>\longrightarrow</math> <math>3,361-1 = \dfrac{10000}{c^2}</math></center>
 
<center><math>2,361= \dfrac{10000}{c^2}</math> <math>\longrightarrow</math> <math>c^2=\dfrac{10000}{2,361}</math></center>
 
<center><math>c=\sqrt{\dfrac{10000}{2,361}}</math> <math>\longrightarrow</math> <math>c=65,07914</math></center>
 
Los parámetros obtenidos son <math>a=30,z_0=100,c=65,07914</math>

===Ecuación de la superficie hiperbólica en coordenadas cilíndricas===
Para obtener la ecuación en coordenadas cilíndricas se parte de la expresión anterior en coordenadas cartesianas:
 
<center><math>\dfrac{x^2+y^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
Sean <math>a, c, z_0>0</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 de las coordenadas cartesianas <math>(x,y,z)</math> en coordenadas cilíndricas <math>(ρ,θ,z).</math>
 
<center><math>x=ρ\cdot\cos(θ)</math> 
<math>y=ρ\cdot\sin(θ)</math> 
<math>z=z</math></center>
 
Cambiando <math>x^2+y^2</math> a coordenadas cilíndricas:
 
<center><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></center>
 
Se saca factor común <math>ρ^2</math>:
 
<center><math>x^2+y^2=ρ^2\cdot(\cos^2(θ)+\sin^2(θ))=1</math></center>
 
Aplicando la siguiente identidad trigonométrica:
 
<center><math>\cos^2(θ)+\sin^2(θ)=1</math></center>
 
Sustituyendo en la ecuación
 
<center><math>x^2+y^2=ρ^2\cdot 1 \longrightarrow x^2+y^2=ρ^2</math></center>
 
Al ser los demás datos constantes, no hace falta cambiar nada más. La ecuación en coordenadas cilíndricas queda: 
<center><math>\dfrac{ρ^2}{a^2}-\dfrac{(z-z_0)^2}{c^2} = 1</math></center>
 
===Representación grafica de la torre hiperbólica ===
 
A partir de la ecuación anterior se modela la torre de enfriamiento estándar de altura <math>(H=150m)</math>, radio máximo <math>(Rmáx=55m)</math> y radio mínimo <math>(Rmín=30m)</math>

 
[[Archivo:Torre11.jpg|miniaturadeimagen|500px|derecha|Representación grafica de la torre de enfriamiento estándar ]]
<syntaxhighlight lang="matlab">
clc
clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
u=linspace(0,2*pi,80);
v=linspace(0,H,80);

[Mt,Mz]=meshgrid(u,v);
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

figure
Mx=Mr.*cos(Mt)
My=Mr.*sin(Mt);
Mz=Mz;

surf(Mx,My,Mz, 'FaceColor','white','EdgeColor','black')
axis equal
hold on
grid on
title('Torre de enfriamiento hiperbólica')
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
</syntaxhighlight>

== 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: 
 
<center><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></center>
 
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: 
 
<center><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></center>
 
<center><math>T(\rho, z) = 65 - 38 \left(\frac{z}{150}\right)^{1,8} - 8 + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
<center><math>T(\rho, z) = 57 - 38 \left(\frac{z}{150}\right)^{1,8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}}</math></center>
 
===Representación del campo escalar 𝑇(𝜌,𝑧) como un mapa de colores en una sección vertical y varias secciones transversales===
Se procede a representar el campo escalar de temperaturas como un mapa de colores con una interfaz que simula la de una cámara térmica.
 
 
'''Interpretación'''
 
Como se aprecia, la temperatura más elevada se da en el origen de coordenadas cartesianas <math>(X,Y,Z)=(0,0,0)</math>.
 
* La temperatura en el interior de la torre disminuye según aumenta la coordenada z.
 
* A su vez, la temperatura disminuye según aumenta el alejamiento al eje z, es decir, según aumenta la coordenada ρ (en coordenadas cilíndricas).
 
Por tanto, se aprecia que la variación en las temperaturas es radial y en altura. Este resultado es coherente, teniendo en cuenta que la función escalar de temperatura depende únicamente de los parámetros en coordenadas cilíndricas: z y ρ (y no de θ).
 
 
[[Archivo:grupo_6_campo_temperaturas_vertical.jpg|500px|Representación grafica del campo escalar de temperatura (sección vertical)]]
[[Archivo:grupo_6_campo_temperaturas_transversal.jpg|500px|Representación grafica del campo escalar de temperatura (secciones transversales)]]
 
A continuación se muestra el código en matlab que muestra las dos ventanas gráficas:
<syntaxhighlight lang="matlab">
clc; clear
%se establecen los parámetros de la torre:

H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:

Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIÓN VERTICAL-
z = linspace(0,H,300); %Se establecen los puntos del plano de sección XOZ
x = linspace(-Rmax,Rmax,300);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=Mrho; %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se representa el campo escalar Temperatura interior del plano XOZ de la
%torre
figure('name','Campo de Temperatura en sección XOZ de la torre')
contourf(Mx,Mz,T,60,'LineColor','none');
colormap('hot');
colorbar;
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo escalar de temperaturas en la sección vertical XOZ de la torre:');
tit.Position=[0,156,0]; %se edita la posición del titulo
hold on;

%Se representa la silueta de la torre
contorno_izq = -Mrho; %Lado izquierdo de la silueta
contorno_dch = Mrho; %Lado derecho de la silueta
contorno_alt = z; %Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); %Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); %Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2);%Vista en 2D
xlim([-Rmax-15,Rmax+15]);%se ajusta la ventana gráfica
ylim([0,H]);
grid off;

%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='Temp. (ºC)';
hold off;

%% -REPRESENTACIÓN DEL CAMPO ESCALAR TEMPERATURA EN SECCIONES TRANSVERSALES-
figure('name','Campo de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
rango=linspace(-Rmax,Rmax,200);
[Mx, My] = meshgrid(rango,rango);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);
for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre.
Zi=(H/3)*i; %Cuatro cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));
% T(rho, z)=T=Tbase-dTz.*(Zi./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2)))

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>Rho_z)=NaN;

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1)
pcolor(Mx, My, Ti);
shading interp; %Se suavizar el color

colormap('hot');
title(['Temperatura (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2);%Vista en 2D
xlim([-70,70]);%se ajusta la ventana gráfica
ylim([-70,70]);

hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
%se establecen los valores máximos y minimos de la escala de color
clim([Tcima-5,Tbase]);
%se invierte la dirección de la barra de colores
cbar=colorbar;
cbar.Direction='reverse';
cbar.Title.String='T (ºC)';
colorbar;
end
hold off;
</syntaxhighlight>

=== Cálculo y representación del gradiente del campo escalar Temperatura ===

Con el fin de interpretar la variación de temperaturas dentro de la torre de enfriamiento hiperbólica en cada punto, se estudia el campo vectorial gradiente del campo escalar temperatura.
 
Este campo gradiente tiene la siguiente expresión:
<center><math>∇T(r,θ,z)=(\dfrac{dT}{dρ},\dfrac{1}{ρ}·\dfrac{dT}{dθ},\dfrac{dT}{dz})</math></center>
 
Como el campo escalar temperatura no depende de <math>θ</math>,el gradiente de temperatura tiene solo dos componentes:
 
<center>Componente radial: <math>\dfrac{dT}{dρ}</math></center>
 
<center>Componente vertical: <math>\dfrac{dT}{dz}</math></center>
 
De este modo, se calculan las derivadas parciales respecto a <math>ρ</math> y respecto a <math>z</math>:
 
<center><math>\dfrac{dT}{dρ}=-ΔT_r·\dfrac{-2ρ·Rmax^2}{(Rmax^2-ρ^2)^2}·e^{(\dfrac{-ρ^2}{Rmax^2-ρ^2})}</math></center>
 
<center><math>\dfrac{dT}{dz}=-ΔT_z·n·\dfrac{z^{n-1}}{H^n}</math></center>
 
 
A continuación se muestra la representación del campo vectorial gradiente del campo Temperaturas <math>T(ρ,z):</math>
 
[[Archivo:grupo_6_gradiente_vertical.jpg|500px|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_gradiente_transversal.jpg|500px|Representación grafica del campo escalar de temperatura ]]
 
 
'''Interpretación gráfica y analítica'''
 
* El vector gradiente en cada punto apunta hacia donde la temperatura aumenta más rápidamente (el vector gradiente en cada punto sigue la trayectoria más empinada de ascenso de temperatura). Es por eso que todas los vectores apuntan hacia abajo (coordenada z=0) y hacia el eje de revolución torre (coordenadas x,y=0,0), ya que, como se ha expuesto en el apartado anterior, la temperatura en la torre aumenta según se acerca a la base y al centro de la torre, y por eso todos los vectores gradiente apuntan hacia allí.
 
 
* Otra interpretación importante es que los vectores tienen mayor módulo en los puntos más alejados de la base de la torre. Esto es así debido según se aumenta la altura, hay mayor diferencia de temperatura para la misma distancia. En la cima de la torre el aire caliente entra en contacto con el aire atmosférico, por tanto, cuanto más cerca se esté de este punto, mayor será el enfriamiento. Análogamente, según se acerca a las paredes de la torre, el gradiente tiene mayor módulo, puesto que para un mismo recorrido radial, la variación de temperatura es mayor.
 
A continuación se muestra el código en 'matlab' que genera sendas ventanas gráficas:
<syntaxhighlight lang="matlab">
clc, clear
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;

%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38;%variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8;%variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se establecen los puntos del plano de sección XOZ
z = linspace(0,H,25);
x = linspace(-Rmax,Rmax,25);
[Mx,Mz]=meshgrid(x,z); %se crea la malla de puntos
MR=abs(Mx); %Radio en valor absoluto

%se calcula el campo escalar en la extensión de X y Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-MR.^2./(Rmax^2-MR.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se eliminan los puntos sobrantes que no pertenecen al paraboloide
Mrho=a*sqrt(1+(((Mz-Z0).^2)/c.^2)); %Coordenadas cillindricas: rho
mask=abs(Mx)<=(Mrho+4.*ones(25)); %Recorte del campo escalar a la forma del hiperboloide
T(~mask)=NaN;

%Se calcula el campo gradiente de T
[gradx,gradz]=gradient(T,x,z);

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN VERTICAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOZ')
quiver(x,z,gradx,gradz,'r','LineWidth',0.7)
hold on

%Se representa la silueta de la torre
contorno_izq = -Mrho; % Lado izquierdo de la silueta
contorno_dch = Mrho; % Lado derecho de la silueta
contorno_alt = z; % Coordenadas verticales
plot(contorno_izq,contorno_alt,'k','LineWidth',3); % Silueta izquierda
plot(contorno_dch,contorno_alt,'k','LineWidth',3); % Silueta derecha

%Se establecen los parámetros de la ventana
xlabel('X(m)');
ylabel('Z(m)');
axis equal;
view(2); %Vista en 2D
xlim([-Rmax-15,Rmax+15]); %se ajusta la ventana gráfica
ylim([0,H]);
ax=gca; %se edita el tamaño del titulo.
ax.FontSize = 9.5;
tit=title('Campo vectorial gradiente de temperatura (Sección Vertical):');
tit.Position=[0,153,0]; %se edita la posición del titulo
hold off

%% -REPRESENTACIÓN DEL CAMPO GRADIENTE EN SECCIÓN TRANSVERSAL-
figure('name','Campo vectorial gradiente de Temperatura en sección XOY de la torre')
%Se crea una malla de coordenadas x,y (cartesianas)
x=linspace(-Rmax-10,Rmax+10,20);
y=x;
[Mx,My] = meshgrid(x,y);

%Se convierte a coordenadas cilíndricas para insertar los datos.
Rho=sqrt(Mx.^2+My.^2);

for i=[0,1,2,3];
%Se establece la altura a la que se hace el corte de la torre
Zi=(H/3)*i; %4 cortes equidistantes
Rho_z=a*sqrt(1+(((Zi-Z0).^2)/c.^2));

%Se calcula la temperatura:
Ti=Tbase-dTz*(Zi/H)^n-dTr*(1-exp(-(Rho.^2)./(Rmax^2-Rho.^2)));

%se eliminan los puntos que no pertenecen a la circunferencia.
Ti(Rho>(Rho_z))=NaN;

%Se calcula el gradiente
[gradx,grady]=gradient(Ti,x,y);

%Se representa el campo de temperaturas en la sección
subplot(2,2,i+1);
quiver(Mx,My,gradx,grady,'r','LineWidth',0.7);

title(['Gradiente(T) (z = ' num2str(Zi) ' m)']);
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
axis equal
view(2); %Vista en 2D
xlim([-70,70]); %se ajusta la ventana gráfica
ylim([-70,70]);
hold on

%Se representa el contorno circular
theta=linspace(0,2*pi,100);
x_cont=Rho_z.*cos(theta);
y_cont=Rho_z.*sin(theta);
plot(x_cont,y_cont,'k','LineWidth',2);
end
hold off;
</syntaxhighlight>

=== Representación de superficies isotérmicas ===
Las curvas de nivel son superficies generadas a partir de la asignación de un valor constante a un campo escalar (en <math>ℝ^3</math>). En este caso, el campo escalar es la Temperatura, T(ρ,z), y sus curvas de nivel son las 'Superficies Isotermas', obtenidas al asignar a T un valor de temperatura constante.
 
A continuación se representan las superficies isotérmicas para los valores de 30,40,50 y 60 ºC.
 
[[Archivo:grupo_6_isotermas_vertical.jpg|500px|izquierda|Representación grafica del campo escalar de temperatura ]]
[[Archivo:grupo_6_isotermas_xyz.jpg|500px|derecho|Representación grafica del campo escalar de temperatura ]]
 
'''Interpretación de las superficies de nivel:'''
 
* Las superficies isotermas tienen una forma arqueada hacia abajo, es decir, que son convexas hacia el eje z. Esto se debe a que el campo escalar temperatura (<math>T(ρ,z)</math>), como se ha comentado anteriormente, varía radial y verticalmente, disminuyendo la temperatura según se alejan los puntos del eje z (de revolución) y de la base de la torre. De este modo, para una misma altura (fijando z), la temperatura será mayor en el centro de la sección transversal que en los extremos, como bien se puede apreciar en la segunda figura del apartado 3.1.
* Es importante destacar que, según se asciende en la torre, las curvas isotermas pasan de asemejarse a una esfera a asemejarse a un plano horizontal, gradualmente. La explicación es debida a que el campo escalar Temperatura depende de dos términos, que dependen respectivamente de z y ρ.
 
El término que depende de z es:
 
<center><math>T(z) = -\Delta T_z \left(\frac{z}{H}\right)^n</math></center>
El término que depende de ρ es:
<center><math>T(\rho) = - \Delta T_{\rho} \left(1 - e^{- \frac{\rho^2}{R_{\text{max}}^2 - \rho^2}}\right)</math></center>
 
Por tanto:
* En la base de la torre (z cercano a 0), la superficie isoterma se asemeja a una esfera. Esto se debe a que el término <math>T(z)</math>, es mínimo, por lo que <math>T(ρ)</math> domina el campo escalar. Así, este último término da a las superficies isotermas la forma esférica, ya que el radio (ρ) produce el efecto más importante sobre la temperatura.
* Análogamente, en la cima de la torre (z cercano a H), la superficie se puede asemejar a un plano horizontal, puesto que el término <math>T(z)</math> es máximo. Así, el efecto que produce el radio (ρ), es menor.
 
A continuación se muestra el programa de matlab que representa en una ventana gráfica las superficies isotermas.
<syntaxhighlight lang="matlab">
clc,clear;
%se establecen los parámetros de la torre:
H=150; %altura de la torre de enfriamiento, en metros
Z0=100; %altura de mínimo radio del hiperboloide, en metros
a=30;
c=65.07913735;
Rmax=55;
Rmin=30;
%se establecen los parámetros de temperatura:
Tbase=65; %temperatura de la base en grados centígrados
dTz=38; %variación de temperatura a lo largo del eje z en grados centígrados
Tcima=Tbase-dTz; %temperatura en el punto más alto de la torre en grados centígrados
dTr=8; %variación de temperatura radial en grados centígrados
n=1.8; %exponente de convección

%Se crean los puntos x,y,z
x = linspace(-Rmax,Rmax,200);
y = linspace(-Rmax,Rmax,200);
z = linspace(0,H,200);

[Mx,My,Mz]=meshgrid(x,y,z); %se crea la malla de puntos
Rho=sqrt(Mx.^2+My.^2); %Coordenada rho (cilíndricas)

%se calcula el campo escalar en la extensión de Rho, Z
T=Tbase-dTz.*(Mz./H).^n-dTr.*(1-exp(-Rho.^2./(Rmax^2-Rho.^2))); %Cálculo del campo escalar Temperatura: T(rho,z)

%Se calcula el radio de la torre para cada z
Mr=a*sqrt(1+(((Mz-Z0).^2)/c.^2));

%Se eliminan los puntos que no están DENTRO de la torre (rho<=R(z))
Mask_interior=(Rho<=Mr);

%Se aplica la máscara: asignamos NaN a las temperaturas fuera de la torre.
T(~Mask_interior) = NaN;

%Se representa
figure('Name', 'Superficies Isotérmicas dentro de la Torre Hiperbólica');
hold on;

%Se define el colormap para los colores
clim([Tcima-5,Tbase]);
cmap = hot(45);
colormap(cmap);
cbar.Direction='reverse';

%% BUCLE
Temperaturas=[30,40,50,60]; %Se definen las temperaturas de las superficies
for i=1:4;
T_nivel=Temperaturas(i);
view(3); % Vista 3D
colorbar
%Se representa la superficie isoterma
[caras, vertices] = isosurface(Mx, My, Mz, T, T_nivel);
curva_nivel=patch('faces', caras, 'Vertices', vertices, 'FaceColor', cmap(T_nivel-23,:), 'EdgeColor', 'none');

%Se representa el contorno de la torre (con transparencia)
contorno = Rho - Mr;
[caras_t, vertices_t] = isosurface(Mx, My, Mz, contorno, 0);
h_tower = patch('faces', caras_t, 'vertices', vertices_t, 'FaceColor', [0.5 0.5 0.5], ...
'EdgeColor', 'none', 'FaceAlpha', 0.1); % Gris muy transparente

%Se ajusta de la ventana gráfica
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
axis equal tight;
end
%vista final
hold off
title('Superficies Isotérmicas dentro de la Torre');
lighting gouraud
camlight;
</syntaxhighlight>

=== Cálculo de la masa del aire contenido en la torre ===
La densidad del aire 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:
 
<center><math>\rho(T) = \rho_{\text{aire}} \frac{T_0 + 273}{T + 273}</math></center>
 
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:
 
<center><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></center>
 
Si se desarrolla la expresión se llega a lo siguiente:
 
<center><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></center>
 
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:
<center><math> M_{\text{real}} \approx 781420 \text{ kg}</math></center>
 
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:
 
<center><math>V = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} \rho \, d\rho dz d\theta</math></center>
En la sección 2.2 se obtiene la expresión del hiperboloide de revolución en coordenadas cilíndricas:
 
<center><math>\frac{\rho^2}{a^2} - \frac{(z - z_0)^2}{c^2} = 1</math></center>
Despejando se obtiene una relación entre <math>\rho</math> y <math>z</math> 
 
<center><math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math></center>
Volviendo a la integral, integrando primero <math>\theta</math> y después <math>\rho</math> se obtiene lo siguiente:
 
<center><math>V = 2\pi \int_0^{150} \left[ \frac{\rho^2}{2} \right]_0^{\rho(z)} dz = \pi \int_0^{150} \rho^2(z) dz</math></center>
 
Sustituyendo <math>\rho^2(z) = a^2 \left[ 1 + \frac{(z - z_0)^2}{c^2} \right]</math> con <math>z_0=100</math>:
 
<center><math>V = \pi a^2 \int_0^{150} \left[ 1 + \frac{(z - 100)^2}{c^2} \right] dz</math>
 
<math>V = \pi a^2 \left[ \int_0^{150} dz + \frac{1}{c^2} \int_0^{150} (z - 100)^2 dz \right]</math></center>
 
Resolviendo ambas integrales:
 
<center><math>\int_0^{150} dz = [z]_0^{150} = 150</math></center>
 
Para la realización de <math>\int_0^{150} (z - 100)^2 dz</math> usamos el cambio de variable <math>u = z - 100, du=dz</math> con límites <math>u(0)=-100, u(150)=50</math>
 
<center><math>\int_{-100}^{50} u^2 du = \left[ \frac{u^3}{3} \right]_{-100}^{50} = \frac{1}{3} \left[ 50^3 - (-100)^3 \right] = \frac{1}{3} [125000 + 1000000] = \frac{1125000}{3} = 375000</math></center>
 
Sustituyendo en la integral, con <math>a=30</math> y <math>c^2 =\frac{10000}{2.361}</math>:
 
<center><math>V = 900\pi \left[ 150 + \frac{375000}{\frac{360000}{85}} \right] = 900\pi \left[ 150 + \frac{375000 \cdot 85}{360000} \right]</math></center>
 
<center><math>V = 900\pi \left[ 150 + \frac{31875000}{360000} \right] = 900\pi \left[ 150 + 88.54167 \right]</math></center>
 
<center><math>V \approx 900\pi (238.54167) \approx 674185.87 \text{ m}^3</math></center>
 
Por tanto, suponiendo que la masa de aire estuviera toda ella a 15ºC, esta tendría la siguiente masa:
 
<center><math>M_{\text{15ºC}} = V \cdot \rho_{\text{aire}} = 674185.87 \text{ m}^3 \cdot 1.225 \text{ kg/m}^3</math></center>
 
<center><math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math></center>
 
Comparando ambos valores: <math>M_{\text{real}} \approx 781420 \text{ kg}</math> y <math>M_{\text{15ºC}} \approx 825377 \text{ kg}</math> se observa que el aire tiene una masa aproximadamente un '''5% menor''' que la masa que tendría si estuviera a temperatura constante <math>T=15ºC</math>.
 

===Cálculo de la temperatura media y la potencia térmica evacuada===
La temperatura media se define como:
 
<center><math>T_{\text{media}} = \frac{1}{V} \iiint_V T(\rho, z) dV = \frac{1}{V} \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math></center>
 
Se definirá <math>I_T=\iiint_V T(\rho, z) dV = \int_0^{2\pi} \int_0^{H} \int_0^{\rho(z)} T(\rho, z) \cdot \rho \, d\rho dz d\theta</math>
Sustituyendo y definiendo los límites de integración se llega a la siguiente integral:
 
<center><math>I_T = 2\pi \int_0^{150} \int_0^{\rho(z)} \left[ 57 - 38 \left(\frac{z}{150}\right)^{1.8} + 8 \cdot e^{-\frac{\rho^2}{55^2 - \rho^2}} \right] \rho \, d\rho dz</math></center>
 
Debido al término exponencial esta integral no es de fácil resolución, luego utilizando una herramienta de internet se llega al siguiente resultado:
 
<center><math>I_T \approx 19534420 \text{ m}^3 \cdot \text{°C}</math></center>
 
Luego la temperatura media será:
 
<center><math>T_{\text{media}} = \frac{I_T}{V} = \frac{19534420 \text{ m}^3 \cdot \text{°C}}{674185.87 \text{ m}^3}</math></center>
 
<center><math>T_{\text{media}} \approx 28.97 \text{ °C}</math></center>
 
La potencia térmica evacuada puede calcularse mediante la siguiente expresión:
 
<center><math>P = \rho_{\text{aire}} c_p Q (T_{\text{media}} - T_{\text{ambiente}})</math></center>
 
Siendo:
*<math>\rho_{\text{aire}} = 1.225 \text{ kg/m}^3</math> la densidad del aire a temperatura <math>T_{\text{ambiente}}=15ºC</math>;
*<math>c_p = 1005 \text{ J/(kg K)}</math> el calor específico a presión constante del aire y <math>T_{\text{ambiente}}=15ºC</math>;
*<math>Q = 1500 \text{ m}^3/\text{s}</math> el caudal de aire;
*<math>T_{\text{media}} = 28.97 \text{ °C}</math> la temperatura media calculada anteriormente;
*<math>T_{\text{ambiente}} = 15 \text{ °C}</math> la temperatura ambiente.
 
Sustituyendo:
 
<center><math>P = 1.225 \cdot 1005 \cdot 1500 \cdot (28.97 - 15)</math></center>
 
<center><math>P = 1845600 \cdot 13.97</math></center>
 
<center><math>P \approx 25781000 \text{ W}=25.781 \text{ MW}</math></center>
 

== Análisis de la 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:
 
<center><math>V(z)=V_0·(\dfrac{z}{Zref})^α</math></center>
 
 <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>.

===.-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:
 
<center><math>P(z)=\dfrac{1}{2}ρ·V(z)^2</math></center>
 
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:
 
<center><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></center>
 
Como vemos analíticamente y veremos a continuación gráficamente, la presión es máxima cuando la altura es máxima, es decir a 150 metros, en la parte superior de la torre de enfriamiento, siendo esta de:
 
<center><math>P(z=150)=\dfrac{1}{2}ρ·V(z=150)^2=\dfrac{1}{2}ρ·(V_0·(\dfrac{z}{Zref})^α)^2=\dfrac{1,225}{2}·(18·(\dfrac{z}{10})^{0,1429})^2≈102.79·150^{0,286}=430.837Pa</math></center>
 

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:

[[Archivo:Captura de pantalla 2025-12-01 125000.png|600px|miniaturadeimagen|derecha|Campo escalar de presiones sobre la mitad de la torre expuesta al viento]]

<syntaxhighlight lang="matlab">
%determinación de los parámetros
a = 30; z0 = 100; c = 65,07914;
V0 = 18;
zref = 10;
alpha = 0.1429;
rho_air = 1.225;
%densidad en [kg/m3], (dato extraido de worldpress.com)
v = linspace(0, 150, 100);
u = linspace(0, pi, 100);
%mallado y calculos
[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;
%representacíon de la superficie
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)');
</syntaxhighlight>

=== 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:
 
<center><math>\vec{F}(x, y, z)=-P(z)·\vec{n}</math></center>
 
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:

[[Archivo:Captura_de_pantalla_2025-12-01_124522.png|600px|miniaturadeimagen|derecha|Campo de fuerzas sobre la mitad de la torre de enfriamiento]]

<syntaxhighlight lang="matlab">
a = 30; c = 65,07914; z0 = 100;
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;
</syntaxhighlight>

=== Fuerza total y por unidad de superficie del viento sobre la mitad expuesta ===
Teniendo en cuenta la presión dinámica total calculada previamente P(z), se calcula la fuerza total sobre la mitad frontal partiendo de la expresión (integrando en θ la mitad frontal da un factor 2):
 
<center><math>F = 2 \int_0^H p(z)\, R(z)\, dz= 2 \int_0^H \frac{1}{2}\rho\, V(z)^2\, R(z)\, dz= \rho \int_0^H V(z)^2 R(z)\, dz .</math></center>
 
Sustituyendo el perfil de viento y el radio lineal:
 
<center><math>V(z)^2 = V_0^{\,2} \left( \frac{z}{Z_{\mathrm{ref}}} \right)^{2\alpha}, R(z) = R_{\min} + \frac{\Delta R}{H} z, </math></center>
 
Queda
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\int_0^H z^{2\alpha} \left( R_{\min} + \frac{\Delta R}{H}z \right)\, dz.</math></center>
 
Separando integrales:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min} \int_0^H z^{2\alpha}\, dz+\frac{\Delta R}{H} \int_0^H z^{2\alpha+1}\, dz\right].</math></center>
 
Las integrales elementales resueltas son:
 
<center><math>\int_0^H z^{2\alpha}\, dz= \frac{H^{2\alpha+1}}{2\alpha+1},\qquad\int_0^H z^{2\alpha+1}\, dz= \frac{H^{2\alpha+2}}{2\alpha+2}.</math></center>
 
Por tanto:
 
<center><math>F = \rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left[R_{\min}\,\frac{H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R}{H}\,\frac{H^{2\alpha+2}}{2\alpha+2}\right].</math></center>
 
Finalmente
<center><math>\boxed{F =\rho V_0^{\,2} Z_{\mathrm{ref}}^{-2\alpha}\left(\frac{R_{\min}\, H^{2\alpha+1}}{2\alpha+1}+\frac{\Delta R\, H^{2\alpha+1}}{2\alpha+2}\right)}</math></center>
 
Los datos utilizados son:
 
*<math>H = 150\ \text{m}</math>,altura de la chimenea
*<math>R_{\min} = 30\ \text{m}, \qquad R_{\max} = 55\ \text{m}</math> radio minimo y máximo
*<math>R(z) = R_{\min} + \frac{\Delta R}{H} z</math> ley de variación lineal del radio
*<math>V_0 = 18\ \text{m/s}, \qquad Z_{\text{ref}} = 10\ \text{m}, \qquad\alpha = 0.1429</math> valores numéricos del perfil de viento
*<math>\rho = 1.225\ \text{kg/m}^3</math> densidad del aire
*<math>p(z) = \frac{1}{2}\,\rho\, V(z)^2</math> presión dinámica del viento
 
Calculando numéricamente la expresión se obtiene:
 
<center><math>F \approx 4.423835 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
<center><math>F = 4.423835 \times 10^{6}\ \text{N}.</math></center>
 
Es necesario calcular el área de la mitad de la torre de enfriamiento, la superficie lateral de revolución correspondiente a la mitad frontal es:
 
<center><math>A_{\text{half}}=\pi \int_0^{H} R(z)\, \sqrt{1 + [R'(z)]^{2}}\, dz.</math></center>
Dado que
<center><math>R'(z)=\frac{\Delta R}{H},</math></center>
es constante, el factor geométrico puede sacarse fuera de la integral:
 
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\int_0^{H} R(z)\, dz.</math></center>
 
La integral del radio es:
 
<center><math>\int_0^{H} R(z)\, dz=\int_0^{H} \left( R_{\min} + \frac{\Delta R}{H} z \right) dz=R_{\min} H + \frac{\Delta R}{2}H.</math></center>
 
Por tanto:
<center><math>A_{\text{half}}=\pi \sqrt{1 + R'^2}\left[R_{\min} H+\frac{\Delta R}{2}H\right].</math></center>
 
<center><math>\boxed{A_{\text{half}}=\pi H \sqrt{1 + \left( \frac{\Delta R}{H} \right)^2}\left( R_{\min} + \frac{\Delta R}{2} \right)}.</math></center>
 
<center><math>A_{\text{half}} \approx 2.0303909708635903 \times 10^{4}\ \text{m}^2,\qquad(\approx 2.03\times 10^{4}\ \text{m}^2).</math></center>
 
La fuerza por unidad de superficie es:
 
<center><math>\boxed{p_{\text{avg}} = \frac{F}{A_{\text{half}}}}.</math></center>
 
Con los valores obtenidos:
 
<center><math>p_{\text{avg}} \approx 217.881\ \text{Pa}.</math></center>

===Comparación con una torre cilíndrica===
Torre cilíndrica de radio 𝑅 = 𝑅max = 55 m y misma altura H = 150m.
 
Realizando los mismos cálculos que en el apartado anterior se compara la forma hiperbólica de una torre de enfriamiento con una cilíndrica.
Pero hay que tener en cuenta que es un cilindro, por lo tanto:
 
<center><math>R(z) = R = \text{cte}</math></center> El radio es constante
 
Por lo que las fórmula que se deben utilizar son:
 
<center><math>\boxed{F_{\text{cil}} =\rho\,R\,V_0^{2}\,Z_{\text{ref}}^{-2\alpha}\,\frac{H^{\,2\alpha+1}}{2\alpha+1}}</math></center>
 
<center><math>\boxed{A_{\text{cil}} = \pi R H}</math></center>
 
<center><math>\boxed{p(z) = \frac{1}{2}\,\rho\,V_0^{2}\left( \frac{z}{Z_{\text{ref}}} \right)^{2\alpha}}</math></center>
 
Fuerza total sobre la mitad frontal
 
Tronco hiperbólico:
 
<center><math>F_{\text{hiperbólico}} \approx 4.4238 \times 10^{6}\ \text{N}\quad (\approx 4.42\ \text{MN}).</math></center>
 
Cilindro:
 
<center><math>F_{\text{cil}} \approx 5.5219 \times 10^{6}\ \text{N}\quad (\approx 5.52\ \text{MN}).</math></center>
 
Relación:
 
<center><math>\frac{F_{\text{cil}}}{F_{\text{hiperbólico}}} \approx 1.248,</math></center>
 
es decir, el cilindro recibe aproximadamente un \(24.8\%\) más de fuerza total.
 
 
Área de la mitad expuesta
 
 
Tronco hiperbólico:
 
<center><math>A_{\text{half, hiperbólico}} \approx 2.0304 \times 10^{4}\ \text{m}^2.</math></center>
 
Cilindro:
 
<center><math>A_{\text{half, cil}} = \pi R H \approx 2.5918 \times 10^{4}\ \text{m}^2.</math></center>
 
Relación:
 
<center><math>\frac{A_{\text{half, cil}}}{A_{\text{half, hiperbólico}}}\approx 1.277,</math></center>
 
por lo que el cilindro tiene aproximadamente un \(27.7\%\) más de área expuesta, la misma relación indica que el cilindro requiere cerca de un \(27.7\%\) más de material superficial.
 
 
Fuerza por unidad de superficie (presión media)
 
 
Tronco hiperbólico:
 
<center><math>p_{\text{avg, hiperbólico}}= \frac{F_{\text{hiperbólico}}}{A_{\text{half, hiperbólico}}}\approx 217.88\ \text{Pa}.</math></center>
 
Cilindro:
 
<center><math>p_{\text{avg, cil}}= \frac{F_{\text{cil}}}{A_{\text{half, cil}}}\approx 213.05\ \text{Pa}.</math></center>
 
Aunque el cilindro soporta mayor fuerza total, su presión media es ligeramente menor, debido a que tiene una superficie expuesta significativamente mayor y la carga se reparte sobre un área más grande.
 
En conclusión, la torre de forma hiperbólica soporta menos fuerza lateral y usa menos material que si fuese de forma cilíndrica, aunque esta soporta menos presión media debido a sus mayor superficie. Esto muestra que el uso de una forma hiperbólica no es cuestión estética sino de mayor eficiencia en resistencia y en coste.

== Aplicaciones practicas del hiperboloide de revolución ==
[[Archivo:Cofrentes.png|250px||miniaturadeimagen|derecha|Torres de enfriamiento de la central nuclear de Cofrentes ]]
[[Archivo:Didcot-power-station-5196798.jpg|miniaturadeimagen|centro|Torres de enfriamiento hiperbólicas de la central térmica de Didcot ]]
[[Archivo:Imagentorres.jpg ‎|miniaturadeimagen|derecha|Torres de enfriamiento de la central termoeléctrica de Niederaussem (Alemania).]]
 
Tras la investigación hecha anteriormente sobre las ventajas que tiene este tipo de geometría, estructurales y económicas, se presentaran algunos ejemplos del uso de esta geometría en las torres de enfriamiento y de la aplicación en otras estructuras mas del ámbito arquitectónico.
 
===Ejemplo de torres de enfriamiento===

 
Las torres de enfriamiento de las centrales energéticas son la máxima expresión del hiperboloide de revolución en la construcción. Ejemplos muy representativos son la central nuclear de Cofrentes en España, Didcot en Reino Unido y Niederaussem en Alemania.
 
==== Central nuclear de Cofrentes ====
 
Esta central se sitúa en la provincia de Valencia. Cuenta con dos torres de enfriamiento de 129 metros de altura. Tiene un diámetro de base de 90 metros y gracias a su funcionamiento ayudan a reducir el consumo de agua del río Júcar.
 
==== Central térmica de Didcot ====
 
La central térmica de Didcot se encuentra en Reino Unido y esta compuesta de dos partes, Didoct A y Didosct B . La central es famosa por sus seis torres de enfriamiento hiperbólicas distribuidas en dos grupos de tres localizadas en la parte A. Esta parte ceso su funcionamiento en 2013 por no cumplir con los estatutos de la Unión Europea sobre grandes instalaciones de combustión. En 2014 empezó su demolición hasta 2019 que culmino con la destrucción de las tres torres del norte. Estas estructuras median 114 metros y eran parte de Oxfordshire.
 
====Central termoeléctrica de Niederaussem ====
 
Esta central se encuentra en Alemania y alberga una de las mayores torres de enfriamiento hiperbólico en el mundo. Tuvo el récord de la torre hiperbólica mas grande desde 2002 hasta 2012, cuando le fue arrebatado por la torre de la central térmica de Kalisindh, situada en Rajastán (India). La torre más importante de Niederaussem tiene una altura de <math>200 m</math> y cuenta con un diámetro en la base de <math>152 m</math>. A pesar de su altura el grosor del hormigón en algunas partes es de tan <math>22 cm</math>.

=== Ejemplos de estructuras hiperboloides en la arquitectura e ingeniera ===
[[Archivo:11Photo0jpg.jpg|miniaturadeimagen|derecha|Planetario James S. McDonnell (San Luis, EE. UU.)]]
[[Archivo:61jl35o0nyt61.jpg|miniaturadeimagen|derecha|Restaurante los Manantiales (Xochimilco, México)]]
 
Las torres de enfriamiento no es la única estructura que emplea geometría hiperbólicas. Esta geometría es empleada en todo tipo de construcciones como torres, puentes, tejados e incluso fachadas. El planetario James S. McDonnell (San Luis, EE. UU.) y el el restaurante los Manantiales (Xochimilco, México) son dos ejemplos del uso de la geometría hiperbólica.
 
==== El planetario James S. McDonnell (San Luis, EE. UU.) ====
 
Este edificio tiene estructura de hiperboloide de revolución, es la misma estructura que las torres de enfriamiento pero empleada para un edificio cultural. Fue construido en 1963 por el arquitecto Gyo Obata. Este eligió esta forma para buscar un edificio que conectara la tierra con el espacio, lo que simboliza el planetario. El planetario es un hiperboloide de una hoja, gracias a esta forma las paredes son delgadas, reduciendo el material empleado.
 
==== Restaurante Los Manantiales (Xochimilco, México) ====
==== Restaurante Los Manantiales (Xochimilco, México) ====

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T12:08:01Z

Ramón Castán Naval: /* .-Representación del campo escalar de presiones. */

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T12:00:41Z

Ramón Castán Naval: /* Análisis de la presión del viento */

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T11:53:07Z

Ramón Castán Naval: /* Análisis de la presión del viento */

Archivo:Captura de pantalla 2025-12-01 125000.png

2025-12-01T11:49:59Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T11:48:45Z

Ramón Castán Naval: /* .-Representación del campo escalar de presiones. */

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T11:47:43Z

Ramón Castán Naval: /* Campo de fuerza en la superficie expuesta */

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-12-01T11:46:55Z

Ramón Castán Naval: /* Campo de fuerza en la superficie expuesta */

Archivo:Captura de pantalla 2025-12-01 124522.png

2025-12-01T11:46:12Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T16:05:17Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:57:27Z

Ramón Castán Naval:

Archivo:Captura de pantalla 2025-11-25 165346.png

2025-11-25T15:53:40Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:50:06Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:48:59Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:47:36Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:38:48Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:29:00Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:27:29Z

Ramón Castán Naval:

Torres de enfriamiento hiperbólicas (Grupo 06)

2025-11-25T15:27:10Z

Ramón Castán Naval: /* Campo de temperaturas y transferencia de calor */