Torres de enfriamiento hiperbólicas (grupo 33)

Trabajo realizado por estudiantes
Título	Torres de enfriamiento Hiperbólicas. Grupo 33
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Marcos Sanchez Martínez Guillermo Garrido Torres Carlos Aguado Esparrells Hector Perucho Conde
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Las torres de enfriamiento hiperbólicas son estructuras fundamentales en la industria energética, ampliamente utilizadas desde mediados del siglo XX debido a su alta eficiencia en la transferencia de calor. Estas torres, con su distintiva forma hiperbólica, combinan propiedades geométricas y mecánicas que las hacen tanto resistentes como funcionales para optimizar los procesos de enfriamiento en plantas termoeléctricas y nucleares.

En este trabajo, se analiza el diseño y comportamiento de una torre de enfriamiento hiperbólica típica, considerando su geometría, las fuerzas inducidas por el viento en su superficie, y el campo de temperatura dentro de la estructura.

derecha

Consideremos una torre de enfriamiento hiperbólica, caracterizada por su altura total H, su radio máximo en la base Rmax, y su radio mínimo Rmin alcanzado a 2/3 de la altura H de la torre.

La superficie de la torre sigue la forma de un hiperboloide hiperbólico con centro a la altura 2/3H, el cual, en coordenadas cartesianas, tiene la siguiente forma:

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

Se pueden suponer los siguientes parámetros:

[math]Rmax=50m,\qquad Rmin=50m,\qquad H=120m[/math]

Por otro lado, el viento ejerce una presión lateral que varía a lo largo de la superficie de la torre. Considerando que la velocidad escalar del viento aumenta con la altura, podemos representarla con la función:

[math]V(z)=V_o(\dfrac{z}{z_o})^α [/math]

Donde:

- [math]V_0[/math] es la velocidad de referencia del viento a una altura [math]z_0[/math]. Para una simulación de viento, podemos fijar [math]V_0 = 15 m/s[/math] como valor de referencia.

- α es un exponente que depende del terreno; para ´areas abiertas suele ser alrededor de 0.14.

Utilizando esta velocidad del viento, la presión del viento sobre la superficie de la torre se puede modelar como:

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

donde ρ es la densidad del aire estándar.

Para determinar la distribución de las fuerzas laterales, calculamos el campo vectorial de las fuerzas inducidas por el viento en la superficie de la torre:

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

donde [math]\vec{n}[/math] es el vector normal a la superficie.

Por ultimo, Suponemos que dentro de la torre de enfriamiento se da un campo de temperatura representado mediante la ecuación:

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

1 Encontrar los valores de a, c, [math]z_0[/math] para la ecuación de la torre

Con la ecuación: [math]\dfrac{x^2+y^2}{a^2} - \dfrac{(z - z_0)^2}{c^2} = 1[/math] en cartesianas, lo primero es pasar la fórmula del hiperboloide a cilíndricas que es [math]\dfrac{\rho^2}{a^2} - \dfrac{(z - z_0)^2}{c^2} = 1[/math]. Ya que en cilíndricas [math] \rho^2=x^2+y^2 [/math]

Ahora se tiene en cuenta los parámetros:

La altura (H) es 120 metros

En z=0 (la base) el radio máximo es 50 metros

En [math]z=\dfrac{2}{3} \cdot H=80[/math] el radio mínimo es 20 metros

Lo primero es saber que el hiperboloide está centrado a una altura [math]z_0[/math](dado por el enunciado)​, la cual se toma como el centro de simetría. Por la simetría del problema, sabemos que la curva alcanza su mínimo radio en [math]z=z_0[/math]​. Esto implica que:

[math]z_0=\dfrac{2}{3} \cdot H=80m[/math]

Así que la ecuación(1) situada en el centro de simetría con los parámetros:

[math] \rho=20, \qquad z=z_0=80[/math]

[math]\dfrac{20^2}{a^2} - \dfrac{(80-80)^2}{c^2} = 1 → \qquad \dfrac{20^2}{a^2}=1 → \qquad a=20 [/math]

Por otro lado tenemos la ecuación(2) situada en la base con parámetros:

[math] \rho=50, \qquad z=0, \qquad z_0=80[/math]

[math]\dfrac{50^2}{a^2} - \dfrac{(0-80)^2}{c^2} = 1[/math]

Por último se sustituye el valor de a en (2) y se consigue el valor de c

[math]c≈34,91[/math]

Así que la ecuación de la torre en cilíndricas con los valores ya sustituidos de a,c y [math]z_0[/math] es:

[math]\dfrac{\rho^2}{20^2} - \dfrac{(z-80)^2}{34,91^2} = 1 → \qquad \dfrac{\rho^2}{400} - \dfrac{(z-80)^2}{1219} = 1[/math]

2 Representación de la superficie parametrizada

Lo primero, es parametrizar la superficie conseguida en la tarea 1 qué es: [math]\dfrac{\rho^2}{20^2} - \dfrac{(z-80)^2}{34,91^2} = 1 → \qquad \dfrac{\rho^2}{400} - \dfrac{(z-80)^2}{1219} = 1[/math]

La parametrización que le damos nosotros para la representación es:

r(u,v)=(p(v)[math]\cdot[/math]cos(u) , p(v)[math]\cdot[/math]sen(u) , v)

Donde [math]p(v)=\sqrt{400 \cdot (1+ \dfrac{(v-80)^2}{1219})}[/math]

Y u pertenece a [0,2pi] y v pertenece a [0,120]

Por lo que es la representación en matlab es:

% Parámetros de la torre 
Rm=50% Parámetro a (radio max)
a = 20; % Parámetro a (radio mínimo)
z0 = 80; % Centro del hiperboloide
H = 120; % Altura total de la torre
c = sqrt(((6400*a^2)/(Rm^2-a^2))); % Parámetro c


% Dominio de los parámetros 
theta = linspace(0, 2*pi, 100); % Ángulo theta [0, 2pi] 
z = linspace(0, H, 100); % Altura z [0, H] 
% Crear mallas para theta y z 
[Theta, Z] = meshgrid(theta, z); 
% Radio \rho en función de z 
Rho = sqrt(a^2 * (1 + ((Z - z0).^2) / c^2)); 

% Coordenadas cartesianas 
X = Rho .* cos(Theta); % Coordenada x 
Y = Rho .* sin(Theta); % Coordenada y 

% Graficar la superficie 
figure; 
surf(X, Y, Z, 'EdgeColor', 'none'); % Representar superficie sin bordes 
colormap([0.5, 0.5, 0.5]); % Color único (gris) 
xlabel('X (m)'); 
ylabel('Y (m)'); 
zlabel('Z (m)'); 
title('Superficie de la Torre de Enfriamiento'); 
axis equal; % Ejes proporcionales 
view(3); % Vista 3D

En esta representación se puede observar la torre de enfriamiento hiperbólica cuya altura es 120 metros, cuyo Radio máximo es alcanzado en la base(z=0) como se puede observar y cuyo radio mínimo es alcanzado a 2/3H (z=80)

3 Ecuación de la torre como una superficie rigada

Una superficie rigada es una superficie generada al mover una recta a lo largo de una trayectoria. El hiperboloide de una hoja es una superficie reglada, y su representación puede expresarse como:

[math]r(u,v)=r_0(u)+v \cdot d(u)[/math]

donde [math]r_0(u)[/math] es una curva generadora, en cambio d(u) es una dirección de la recta generadora.

Para el hiperboloide:

1. Podemos parametrizar una curva generadora en [math]z_0[/math]=cte como: [math]r_0(u)= (a \cdot cos(u) , a \cdot sen(u) , z_0)[/math]

2. Una dirección de la recta generadora es: [math]d(u)= (sinh(v) \cdot cos(u) , sinh(v) \cdot sin(u) , cosh(v))[/math]

Por lo que la parametrización rigada es: [math]r(u,v)= (a \cdot cos(u)+v \cdot sinh(v) \cdot cos(u) , a \cdot sen(u)+v \cdot sinh(v) \cdot sin(u) , z_o+v \cdot cosh(v))[/math]

4 Representación del campo escalar de presión como un mapa de colores sobre la superficie parametrizada de la torre

Esta pregunta se centra en representar gráficamente el campo escalar de presión inducido por el viento en la superficie de la torre. Esto permite comprender cómo varía la presión a lo largo de la estructura en función de la altura, un factor crucial en el diseño y la estabilidad estructural.

% Crear la malla para parametrizar la superficie
angulo = linspace(0, 2*pi, 100); % Ángulo alrededor del eje z
altura = linspace(0, altura_torre, 200); % Altura
[angulo_malla, altura_malla] = meshgrid(angulo, altura); % Crear mallas 2D para ángulo y altura
% Parametrización del hiperboloide
coord_x = radio_min_torre * cos(angulo_malla) .* sqrt(1 + ((altura_malla - centro_hiperboloide) / parametro_c).^2);
coord_y = radio_min_torre * sin(angulo_malla) .* sqrt(1 + ((altura_malla - centro_hiperboloide) / parametro_c).^2);
coord_z = altura_malla;
% Perfil del viento
vel_viento = vel_viento_ref * (altura_malla / altura_ref_viento).^exponente_viento;
% Presión del viento
presion_viento = 0.5 * densidad_aire * vel_viento.^2;
% Máscara para dividir la torre: seleccionar puntos donde coord_x + coord_y < 0
mascara_division = (coord_x + coord_y < 0);
coord_x(~mascara_division) = NaN; % Asignar NaN para ocultar la otra mitad
coord_y(~mascara_division) = NaN;
coord_z(~mascara_division) = NaN;
presion_viento(~mascara_division) = NaN;
% Representar la torre y el campo de presión
figure;
surf(coord_x, coord_y, coord_z, presion_viento, 'EdgeColor', 'none'); % Campo de presión
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Presión del viento en la mitad expuesta de la torre');
axis equal;
grid on;
view(3);
colormap(jet); % Representación del campo de presiones
colorbar; % Leyenda de colores

(Descripción breve de lo visto en la imagen)

5 Representación del campo vectorial de la fuerza generada por la presión del viento en la superficie de la mitad de la torre expuesta.

En la quinta pregunta, se busca visualizar el campo vectorial de las fuerzas generadas por la presión del viento sobre la superficie de la torre. Este análisis es esencial para estudiar las cargas mecánicas que la estructura debe soportar bajo condiciones de viento.

% Parámetros de la torre
radio_max_base = 50; % Radio máximo de la base
radio_min_torre = 20; % Radio mínimo de la torre
centro_hiperboloide = 80; % Centro del hiperboloide
parametro_c = sqrt(((6400 * radio_min_torre^2) / (radio_max_base^2 - radio_min_torre^2))); % Parámetro c
altura_torre = 120; % Altura total de la torre
densidad_aire = 1.225; % Densidad del aire (kg/m^3)
vel_viento_ref = 15; % Velocidad del viento a altura de referencia (m/s)
altura_ref_viento = 10; % Altura de referencia para el viento (m)
exponente_viento = 0.14; % Exponente del perfil del viento
% Parámetros de la torre
radio_max_base = 50; % Radio máximo
radio_min_torre = 20; % Radio mínimo
centro_hiperboloide = 80; % Centro del hiperboloide
parametro_c = sqrt(((6400 * radio_min_torre^2) / (radio_max_base^2 - radio_min_torre^2))); % Parámetro c
altura_torre = 120; % Altura total
densidad_aire = 1.225; % Densidad del aire (kg/m^3)
vel_viento_ref = 15; % Velocidad del viento a altura de referencia (m/s)
altura_ref_viento = 10; % Altura de referencia para el viento (m)
exponente_viento = 0.14; % Exponente del perfil del viento
% Crear la malla para parametrizar la superficie
angulo = linspace(0, 2*pi, 100); % Ángulo alrededor del eje z
altura = linspace(0, altura_torre, 200); % Altura
[angulo_malla, altura_malla] = meshgrid(angulo, altura); % Mallas 2D para ángulo y altura
% Parametrización del hiperboloide
coord_x = radio_min_torre * cos(angulo_malla) .* sqrt(1 + ((altura_malla - centro_hiperboloide) / parametro_c).^2);
coord_y = radio_min_torre * sin(angulo_malla) .* sqrt(1 + ((altura_malla - centro_hiperboloide) / parametro_c).^2);
coord_z = altura_malla;
% Crear una máscara para cortar la torre según el eje -i-j (X + Y <= 0)
mascara_corte = (coord_x + coord_y) <= 0;
coord_x_corte = coord_x .* mascara_corte;
coord_y_corte = coord_y .* mascara_corte;
coord_z_corte = coord_z .* mascara_corte;
% Crear flechas apuntando al centro (eje Z)
fuerza_x_centro = -coord_x_corte; % Componente en X hacia el centro
fuerza_y_centro = -coord_y_corte; % Componente en Y hacia el centro
fuerza_z_centro = 0 * coord_z_corte; % Sin componente en Z (paralelas al eje Z)
% Escalar las flechas según la presión del viento
factor_escala = coord_z_corte / max(coord_z_corte(:));
factor_escala(isnan(factor_escala)) = 0; % Manejar valores NaN en la base
fuerza_x_centro = fuerza_x_centro .* factor_escala;
fuerza_y_centro = fuerza_y_centro .* factor_escala;
% Filtrar las flechas para evitar puntos nulos
puntos_validos = mascara_corte & (coord_x_corte ~= 0 | coord_y_corte ~= 0);
% --------- Gráfico de las flechas ---------
figure;
hold on;
% Dibujar las flechas
quiver3(coord_x_corte(puntos_validos), coord_y_corte(puntos_validos), coord_z_corte(puntos_validos), ...
      fuerza_x_centro(puntos_validos), fuerza_y_centro(puntos_validos), fuerza_z_centro(puntos_validos), ...
      2, 'Color', [0.5 0.7 1]); % Azul claro (RGB: 50% rojo, 70% verde, 100% azul)
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Campo vectorial de la presión del viento sobre la parte expuesta de la torre');
axis equal; % Mantener proporciones
grid on;
view(3); % Vista en 3D
hold off;

(Descripción breve de lo visto en la imagen)

6 Representación el campo de temperatura utilizando un mapa de colores en un plano vertical que corta la torre pasando por el eje de simetría.

En el apartado 6, se analiza y representa el campo de temperatura en la torre. Por un lado, se estudia cómo varía la temperatura en un plano vertical que corta la torre por el eje de simetría. Por otro lado, se examinan planos paralelos al suelo mediante una animación que permite visualizar la distribución de temperatura en diferentes niveles de la torre.

La animación en matlab es:

% Iniciar figura
figure;
for i = 1:length(altura_malla)
  % Altura de la capa actual
  altura_actual = altura_malla(i);
   % Radio en esta altura
  radio_capa = sqrt((1 + (altura_actual - centro_hiperboloide)^2 / parametro_c^2) * radio_min_torre^2);
   % Coordenadas de la capa
  coord_x_capa = radio_capa * cos(angulo_theta);
  coord_y_capa = radio_capa * sin(angulo_theta);
  coord_z_capa = altura_actual * ones(size(coord_x_capa));
   % Temperatura en esta capa
  temp_capa = temp_base - delta_temp_vertical * (altura_actual / altura_torre)^exponente_conveccion - ...
              delta_temp_radial * (1 - exp(-radio_capa^2 / (radio_max_base^2 - radio_capa^2)));
   % Dibujar la capa
  fill3(coord_x_capa, coord_y_capa, coord_z_capa, temp_capa, 'EdgeColor', 'none');
  hold on;
   % Dibujar las líneas discontinuas en las alturas predefinidas
  for altura_linea = alturas_lineas
      radio_linea = sqrt((1 + (altura_linea - centro_hiperboloide)^2 / parametro_c^2) * radio_min_torre^2); % Radio en esta altura
      coord_x_linea = radio_linea * cos(angulo_theta);
      coord_y_linea = radio_linea * sin(angulo_theta);
      plot3(coord_x_linea, coord_y_linea, altura_linea * ones(size(coord_x_linea)), 'm--', 'LineWidth', 1.5); % Línea discontinua magenta
  end
   % Mostrar el contador de la altura actual en el título
  title(sprintf('Campo de temperatura a z = %.2f m', altura_actual));
   % Configuración gráfica
  colormap('jet');
  colorbar;
  caxis([temp_tope, temp_base]); % Escala de colores
  xlabel('x (m)');
  ylabel('y (m)');
  zlabel('z (m)');
  xlim([-radio_max_base, radio_max_base]);
  ylim([-radio_max_base, radio_max_base]);
  zlim([0, altura_torre]);
  view(3);
   % Pausar para animación
  pause(0.1);
  hold off;
end

Por otro lado, el plano vertical en matlab es:

figure; % Ventana 2: Plano vertical
% Generar plano vertical en xz
x_plane = linspace(-Rm, Rm, 100); % Eje x del plano
z_plane = linspace(0, H, 100);        % Eje z del plano
[X_plane, Z_plane] = meshgrid(x_plane, z_plane);
% Cálculo de la temperatura en el plano
R_plane = abs(X_plane); % Distancia radial en el plano
T_plane = Tbase - DeltaTz * (Z_plane / H).^n - DeltaTr * (1 - exp(-R_plane.^2 ./ (Rmax^2 - R_plane.^2 + eps)));
% Representar el plano vertical
surf(X_plane, zeros(size(X_plane)), Z_plane, T_plane, 'EdgeColor', 'none', 'FaceAlpha', 0.8);
% Configuración gráfica
colormap('jet');
colorbar;
caxis([Ttop, Tbase]); % Escala de colores
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('Plano Vertical: Campo de Temperatura');
xlim([-Rm, Rm]);
ylim([-Rm, Rm]);
zlim([0, H]);
view(3);
hold off;

(Insertar imagenes)

(Breve descripción de las imagenes vistas)

7 Representación del campo del gradiente de temperatura en los puntos de un plano que corta la torre verticalmente pasando por el eje de simetría y los puntos de un plano que corta la torre horizontalmente en correspondencia con el radio mínimo

(Introducción 7)

7.1 Representación del campo del gradiente de temperatura en los puntos de un plano que corta la torre verticalmente pasando por el eje de simetría

El codigo en matlab es:

% Parámetros del modelo
radioBase = 20; % Semi-eje en las direcciones x e y (m)
radioAltura = 34.17; % Semi-eje en la dirección z (m)
alturaCentro = 80; % Altura del centro del hiperboloide
alturaTotal = 120; % Altura total de la torre (m)
radioMaximo = 50; % Radio máximo de la torre (m)
tempInicial = 70; % Temperatura en la base (°C)
tempFinal = 25; % Temperatura en la cima (°C)
variacionAlturaTemp = tempInicial - tempFinal;
variacionRadialTemp = 5; % Variación máxima radial (°C)
expConveccion = 1.5; % Exponente de convección
% Malla para el hiperboloide
anguloTheta = linspace(0, 2*pi, 50); % Ángulo de rotación
altura = linspace(0, alturaTotal, 50); % Coordenada vertical (z)
% Inicializar figura
figure;
% Ciclo de animación para construir el hiperboloide
for paso = 1:length(altura)
   % Limpiar figura de iteración anterior
   clf;
  
   % Dibujar cada capa hasta el paso actual
   hold on;
   for capa = 1:paso
       % Coordenada z de la capa actual
       alturaActual = altura(capa);
      
       % Calcular el radio en esta altura
       radioActual = sqrt((1 + ((alturaActual - alturaCentro)^2) / radioAltura^2) * radioBase^2);
      
       % Coordenadas de la capa en el plano x-y
       coordX = radioActual * cos(anguloTheta);
       coordY = radioActual * sin(anguloTheta);
       capaZ = alturaActual * ones(size(coordX));
      
       % Calcular la temperatura en la capa
       parteAltura = (alturaActual / alturaTotal)^expConveccion;
       variacionRadial = exp(-radioActual^2 / (radioMaximo^2 - radioActual^2));
       temperaturaCapa = tempInicial - variacionAlturaTemp * parteAltura - variacionRadialTemp * (1 - variacionRadial);
      
       % Dibujar la capa con temperatura
       fill3(coordX, coordY, capaZ, temperaturaCapa, 'EdgeColor', 'none');
   end
  
   % Configuración de la gráfica
   colormap('jet');
   colorbar;
   caxis([tempFinal, tempInicial]); % Escala de colores
   xlabel('x (m)');
   ylabel('y (m)');
   zlabel('z (m)');
   title('Construcción del hiperboloide con campo de temperatura');
   xlim([-radioMaximo, radioMaximo]);
   ylim([-radioMaximo, radioMaximo]);
   zlim([0, alturaTotal]);
   view(3);
end

(Insertar imagenes)

(Breve descripción de las imagenes vistas)

7.2 Puntos de un plano que corta la torre horizontalmente en correspondencia con el radio mínimo

(Apartado 2)

8 Animación que representa las superficies isotérmicas para varios valores de temperatura

Una superficie isoterma es aquella en la que todos sus puntos mantienen una misma temperatura. Para representarla, se han seleccionado varias temperaturas y se han visualizado estas superficies junto con sus gradientes correspondientes. Esto permitirá visualizar cómo cambian las superficies y el gradiente asociado de manera dinámica.

La animación en matlab es :

% Parámetros geométricos y térmicos de la torre
radioBase = 20; % Semi-eje horizontal en x e y (m)
alturaCentro = 80; % Altura del centro del hiperboloide (m)
alturaTotal = 120; % Altura completa de la estructura (m)
radioMaximo = 50; % Radio máximo permitido (m)
radioHiperboloide = 34.17; % Semi-eje vertical en z (m)
tempBase = 70; % Temperatura en la base (°C)
tempCima = 25; % Temperatura en la parte superior (°C)
variacionTempAltura = tempBase - tempCima; % Diferencia de temperatura por altura
variacionTempRadial = 5; % Variación máxima radial de temperatura (°C)
expConveccion = 1.5; % Exponente para cálculo térmico
% Niveles isotérmicos
nivelesTemp = linspace(tempCima, tempBase, 6); % Definimos 6 superficies isotérmicas
% Crear malla en coordenadas cilíndricas
angulos = linspace(0, 2*pi, 50); % Discretización angular
alturas = linspace(0, alturaTotal, 50); % Discretización vertical
[AnguloMalla, AlturaMalla] = meshgrid(angulos, alturas);
% Coordenadas del hiperboloide
radioMalla = sqrt((1 + ((AlturaMalla - alturaCentro).^2) / radioHiperboloide^2) * radioBase^2);
xMalla = radioMalla .* cos(AnguloMalla);
yMalla = radioMalla .* sin(AnguloMalla);
% Cálculo del campo de temperaturas en 3D
campoTemp = tempBase - variacionTempAltura * (AlturaMalla / alturaTotal).^expConveccion ...
           - variacionTempRadial * (1 - exp(-radioMalla.^2 ./ (radioMaximo^2 - radioMalla.^2)));
% Gradientes del campo de temperatura
[gradTempX, gradTempZ] = gradient(campoTemp, mean(diff(xMalla(1,:))), mean(diff(AlturaMalla(:,1))));
gradTempY = gradTempX; % Por simetría angular
% Iniciar figura para la animación
figure;
for nivel = 1:length(nivelesTemp)
   % Temperatura isotérmica actual
   tempNivel = nivelesTemp(nivel);
  
   % Máscara lógica para identificar puntos cercanos al nivel isotérmico
   mascaraIso = abs(campoTemp - tempNivel) < 0.5; % Margen de ±0.5 °C
  
   % Extraer puntos de la superficie isotérmica
   xIso = xMalla(mascaraIso);
   yIso = yMalla(mascaraIso);
   zIso = AlturaMalla(mascaraIso);
   coloresIso = campoTemp(mascaraIso); % Asociar colores a la temperatura
   % Graficar la superficie isotérmica
   scatter3(xIso, yIso, zIso, 20, coloresIso, 'filled');
   hold on;
  
   % Graficar el gradiente térmico en los puntos isotérmicos
   quiver3(xIso, yIso, zIso, gradTempX(mascaraIso), gradTempY(mascaraIso), gradTempZ(mascaraIso), ...
           'k', 'AutoScaleFactor', 0.5);
  
   % Configuración de visualización
   colormap('cool'); % Aquí cambiamos el mapa de colores a 'cool'
   colorbar;
   caxis([tempCima, tempBase]); % Escala de colores de temperatura
   xlabel('X (m)');
   ylabel('Y (m)');
   zlabel('Z (m)');
   title(['Superficie isotérmica: T = ', num2str(tempNivel), '°C']);
   xlim([-radioMaximo, radioMaximo]);
   ylim([-radioMaximo, radioMaximo]);
   zlim([0, alturaTotal]);
   view(3);
   grid on;
   % Pausa para animación
   pause(1);
   % Limpiar la figura para el siguiente cuadro
   hold off;
end

(Insertar imagenes)

(Breve descripción de las imagenes vistas)