Introducción

Este documento presenta el desarrollo, simulación y análisis de un modelo térmico dinámico aplicado a un receptor solar de concentración.
El modelo se inspira en la estructura metodológica descrita en el artículo:

Artículo original:
Thermal performance modelling of solar receivers: A comparative review,
Renewable and Sustainable Energy Reviews, 2026.
Disponible en: https://www.sciencedirect.com/science/article/pii/S136403212600568X

El artículo revisa distintos enfoques de modelización térmica y concluye que los modelos dinámicos de primer orden, basados en balances de energía, son adecuados para estudiar:

la evolución temporal de la temperatura del receptor,
el equilibrio térmico,
la influencia de parámetros físicos y ambientales,
y la comparación entre escenarios de operación.

El modelo implementado en este documento corresponde exactamente a esa categoría.

La aplicación Shiny asociada permite explorar el comportamiento térmico de forma interactiva:

URL de la app Shiny:
https://r-teresa-fraile.shinyapps.io/shiny/

Función del modelo en el artículo original

El artículo clasifica los modelos térmicos en varias categorías: estacionarios, dinámicos, CFD y modelos híbridos.
El modelo utilizado aquí pertenece a los modelos dinámicos simplificados, cuya función principal es:

predecir la evolución temporal de la temperatura del receptor,
comparar escenarios de operación y diseño,
analizar la sensibilidad a parámetros como masa, emisividad, coeficiente convectivo,
optimizar el rendimiento térmico del sistema.

Este tipo de modelo es especialmente útil cuando se busca un equilibrio entre precisión y coste computacional, tal como se discute en el artículo.

Modelo térmico

El modelo parte de la ecuación diferencial:

\[ m \cdot C_p \cdot \frac{dT}{dt} = Q_{\text{in}} - Q_{\text{loss}} \]

donde:

\(m\): masa del receptor
\(C_p\): calor específico
\(Q_{\text{in}}\): potencia térmica absorbida
\(Q_{\text{loss}}\): pérdidas térmicas por convección y radiación

La potencia absorbida se calcula como:

\[ Q_{\text{in}} = \eta \cdot DNI \cdot A \]

Las pérdidas se modelan como:

\[ Q_{\text{loss}} = h \cdot A \cdot (T - T_{\text{amb}}) + \epsilon \sigma A (T^4 - T_{\text{amb}}^4) \]

Código del modelo

# Parámetros
m <- 50          # kg
Cp <- 500        # J/kgK
A <- 2           # m2
eta <- 0.85
h <- 15          # W/m2K
epsilon <- 0.9
sigma <- 5.67e-8
Tamb <- 300      # K
DNI <- 800       # W/m2

# Ecuación diferencial
modelo <- function(t, T, parms){
  Qin <- eta * DNI * A
  Qloss <- h * A * (T - Tamb) + epsilon * sigma * A * (T^4 - Tamb^4)
  dTdt <- (Qin - Qloss) / (m * Cp)
  list(dTdt)
}

library(deSolve)
tiempo <- seq(0, 3600, by = 10)
sol <- ode(y = 300, times = tiempo, func = modelo, parms = NULL)

Resultados: gráfico principal

plot(sol[,1], sol[,2],
     type = "l",
     col = "red",
     lwd = 2,
     xlab = "Tiempo (s)",
     ylab = "Temperatura (K)",
     main = "Evolución temporal de la temperatura del receptor")

# Gráfico de sensibilidad

DNI_vals <- c(600, 800, 1000)
colores <- c("blue", "red", "darkgreen")

plot(NULL, xlim = c(0, 3600), ylim = c(300, 900),
     xlab = "Tiempo (s)", ylab = "Temperatura (K)",
     main = "Sensibilidad a la irradiancia")

for(i in seq_along(DNI_vals)){
  DNI <- DNI_vals[i]
  sol_i <- ode(y = 300, times = tiempo, func = modelo, parms = NULL)
  lines(sol_i[,1], sol_i[,2], col = colores[i], lwd = 2)
}

legend("bottomright", legend = DNI_vals, col = colores, lwd = 2)

library(ggplot2)

sigma <- 5.67e-8

# Resistencias
R_conv <- function(h) 1 / h
R_wall <- function(e, k) e / k
R_rad <- function(Twall, Tsky, emiss, F = 1) 1 / (F * emiss * sigma * (Twall^4 - Tsky^4))

# Pérdidas
q_loss_conv <- function(Twall, Tair, h_ext) h_ext * (Twall - Tair)
q_loss_rad  <- function(Twall, Tsky, emiss, F = 1) emiss * sigma * F * (Twall^4 - Tsky^4)

# Eficiencia
eta_receptor <- function(alpha, q_sun, q_loss_total) alpha * (1 - q_loss_total / q_sun)

# Modelo indirecto
modelo_indirecto <- function(q_sun, alpha = 0.9, emiss = 0.85,
                             Tfluid = 700 + 273, Tair = 20 + 273,
                             Tsky = -20 + 273, h_int = 1200,
                             h_ext = 20, e = 0.003, k = 30) {
  
  Twall_int <- Tfluid + q_sun * R_conv(h_int)
  Twall_ext <- Twall_int + q_sun * R_wall(e, k)
  
  q_conv <- q_loss_conv(Twall_ext, Tair, h_ext)
  q_rad  <- q_loss_rad(Twall_ext, Tsky, emiss)
  q_loss_total <- q_conv + q_rad
  
  eficiencia <- eta_receptor(alpha, q_sun, q_loss_total)
  
  return(data.frame(
    q_sun, Twall_int, Twall_ext, q_conv, q_rad, q_loss_total, eficiencia,
    tipo = "Indirecto"
  ))
}

# Modelo directo
modelo_directo <- function(q_sun, alpha = 0.95, emiss = 0.85,
                           Tpart = 800 + 273, Tair = 20 + 273,
                           Tsky = -20 + 273, h_adv = 150) {
  
  Twall <- Tpart
  
  q_conv <- q_loss_conv(Twall, Tair, h_adv)
  q_rad  <- q_loss_rad(Twall, Tsky, emiss)
  q_loss_total <- q_conv + q_rad
  
  eficiencia <- eta_receptor(alpha, q_sun, q_loss_total)
  
  return(data.frame(
    q_sun, Twall_int = Twall, Twall_ext = Twall,
    q_conv, q_rad, q_loss_total, eficiencia,
    tipo = "Directo"
  ))
}

# Rango de flujos solares
flujos <- seq(200, 2000, length.out = 200)

# Simulaciones
datos_ind <- do.call(rbind, lapply(flujos, modelo_indirecto))
datos_dir <- do.call(rbind, lapply(flujos, modelo_directo))
datos <- rbind(datos_ind, datos_dir)

ggplot(datos, aes(q_sun, eficiencia, color = tipo)) +
  geom_line(size = 1.4) +
  scale_color_manual(
    values = c("Directo" = "#0072B2", "Indirecto" = "#D55E00"),
    name = "Tipo de receptor",
    labels = c("Directo (absorción en pared)", "Indirecto (fluido interno)")
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(size = 14, face = "bold"),
    legend.text  = element_text(size = 12)
  ) +
  labs(
    title = "Eficiencia del receptor vs Flujo solar",
    x = "Flujo solar (W/m²)",
    y = "Eficiencia"
  )

ggplot(datos, aes(q_sun, q_conv, color = tipo)) +
  geom_line(size = 1.4) +
  scale_color_manual(
    values = c("Directo" = "#0072B2", "Indirecto" = "#D55E00"),
    name = "Tipo de receptor",
    labels = c("Directo (absorción en pared)", "Indirecto (fluido interno)")
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(size = 14, face = "bold"),
    legend.text  = element_text(size = 12)
  ) +
  labs(
    title = "Pérdidas por convección/advección",
    x = "Flujo solar (W/m²)",
    y = "q_conv (W/m²)"
  )

ggplot(datos, aes(q_sun, q_rad, color = tipo)) +
  geom_line(size = 1.4) +
  scale_color_manual(
    values = c("Directo" = "#0072B2", "Indirecto" = "#D55E00"),
    name = "Tipo de receptor",
    labels = c("Directo (absorción en pared)", "Indirecto (fluido interno)")
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(size = 14, face = "bold"),
    legend.text  = element_text(size = 12)
  ) +
  labs(
    title = "Pérdidas por radiación",
    x = "Flujo solar (W/m²)",
    y = "q_rad (W/m²)"
  )

ggplot(datos, aes(q_sun, Twall_ext, color = tipo)) +
  geom_line(size = 1.4) +
  scale_color_manual(
    values = c("Directo" = "#0072B2", "Indirecto" = "#D55E00"),
    name = "Tipo de receptor",
    labels = c("Directo (absorción en pared)", "Indirecto (fluido interno)")
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(size = 14, face = "bold"),
    legend.text  = element_text(size = 12)
  ) +
  labs(
    title = "Temperatura de pared vs Flujo solar",
    x = "Flujo solar (W/m²)",
    y = "Temperatura pared (K)"
  )

Modelo térmico de receptor solar: simulación y análisis

Teresa Fraile

2026-07-29

Introducción

Función del modelo en el artículo original

Modelo térmico

Código del modelo

Resultados: gráfico principal