Ajuste de la capacidad calorífica

Capacidad calorifica específica

La capacidad calorífica es una propiedad relevante desde el punto de vista teórico, pero también para el diseño de procesos en los que este involucrados equipos intercambio de calor y reactores químicos.

La capacidad calorífica específica se define como la cantidad de calor que se debe suministar a una unidad de masa para aumentar su temperatura en \(\Delta T\) °C. Esto se escribe matemáticamente así:

\[ Q = m c_p \Delta T,\]

donde \(Q\) es el calor suministrado, \(m\) es la masa y \(c_p\) es la capacidad calorífica específica. Los experimentos para determinar \(c_p\) controlan la masa, el calor suministrado, la presión y la temperatura. Entonces, \(c_p\) se halla despejando de la expresión anterior,

\[c_p = \frac{Q}{m\Delta T}.\]

Normalmente hemos visto como los libros de termodinámica nos presentan tablas de \(c_p\), en las cuales se reportan los coeficientes de un polinomio dado un rango de temperatura, aqui les va un ejemplo extraido de colby, Tabla 7.2.3:

En este pequeño notebook vamos a realizar un ejercicio guiado para ajustar la capacidad calorífica específica del gas metano (\(CH_4\)).

Los datos

Vamos a importar un conjunto de datos experimentales del paquete CHNOSZ.

library(CHNOSZ)

file <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
Cpdat <- read.csv(file)

# Agradecimientos para el alumno Ignacio Barrio por esa actualización ----------
# Filtramos el dataframe para que se quede solo CH4
Cpdat <- Cpdat[Cpdat$species == 'CH4',]

# Eliminamos la columna 'species', que solo aporta el nombre de la especie CH4
Cpdat <- subset (Cpdat, select = -c(species))
#-------------------------------------------------------------------------------
head(Cpdat)

Este paquete tiene una función para cambiar unidades. Por ejemplo, podemos cambiar \(c_p\) de Julios a calorías, y la presión de mega pascales a bares:

Cpdat$Cp <- convert(Cpdat$Cp, "cal")
Cpdat$P <- convert(Cpdat$P, "bar")
head(Cpdat)

Visualicemos los datos:

library(ggplot2)

Use suppressPackageStartupMessages() to eliminate package startup messages

ggplot(mapping = aes(x=T, y=Cp), data = Cpdat) + geom_point()

Podemos observar que hay un cambio d comportamiento alrededor de una temperatura de 669 K. Por lo tanto, requerimos de dos ecuaciones. Dividamos el conjunto de datos para luego ajustar las dos ecuaciones:

Cpdat1 <- Cpdat[Cpdat$T < 660,]
head(Cpdat1)

Cpdat2 <- Cpdat[Cpdat$T >= 660,]
head(Cpdat2)

Regresión clásica

Pregunta 1

Grafique Cpdat1 y Cpdat2 por aparte, y utilice estas gráficas como punto de partida para decidir el grado de su polinomio.

# Grafique en esta celda

---

Pregunta 2

Ahora, debes ajustar una ecuación polinomial para la capacidad calorífica del metano según su temperatura.

# Ajuste en esta celda un modelo con lm() para Cpdat1

# Ajuste en esta celda un modelo con lm() para Cpdat2

---

Herrramientas novedosas

El propio paquete CHNOSZ tiene una función llamada EOSregress(). Esta función considera tras bambalinas dos aspectos super interesantes:

1. Utiliza la ecuación de Helgeson para \(c_p\)

\[c_p = a + \frac{b}{(T-T_0)^2} + \omega TX,\]

donde \(a\) y \(b\) son parámetros relacionados con fenómenos de no solvatación, y \(\omega\) ers un parámetro de solvatación. \(T_0\) es una temperatura de referencia (228 K) y \(X\) es la función de Born que relaciona los procesos de solvatación con las propiedades dieléctricas del agua.

¡Vaya! En este ejercicio no requerimos ninguna de estas sofisticaciones porque el metano es apolar y es gaseoso, pero lo interesante es que podríamos usar este paquete para abordar casos mucho más complejos.

Pregunta 3

Bueno, aqui es importante hacer una pregunta: ¿es la ecuación de Hegelson lineal? Responda, y justifique su respuesta.

Respuesta (escribe a continuación)

---

2. Su algoritmo de ajuste utiliza la misma función lm(). Además, usa internamente la función water() para encontrar aquellas propiedades electrostáticas del agua relacionadas con la solvatación que se mencionó anteriormente.

Por ejemplo, se puede ajustar una ecuación para los datos completos así:

var <- c("invTTheta2", "TXBorn")
Cplm_high <- EOSregress(Cpdat, var)
Cplm_high$coefficients

  (Intercept)    invTTheta2        TXBorn 
     329.3733 -1247884.4469  2052194.0779 

y el paquete trae la siguiente función para graficar. Ejecute la siguiente celda y observe la gráfica.

EOSplot(Cpdat, coeficients = round(Cplm_high$coefficients, 1))

También se le puede especificar el punto de cambio, observe el siguiente código:

Cplm_low <- EOSregress(Cpdat, var, T.max = 600)
Cplm_low$coefficients

  (Intercept)    invTTheta2        TXBorn 
-1.660861e+00  4.459696e+04 -1.295538e+05 

Pregunta 4

Ahora, grafica el modelo Cplm_low por ti mismo usando la función EOSplot()y observa la diferencia con el modelo Cplm_high

# Grafica en esta celda
EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))

summary(Cplm_low)

Call:
lm(formula = fmla)

Residuals:
      1       2       3       4       5 
-0.3177  1.7028 -0.5365 -2.2518  1.4032 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -1.661e+00  5.242e+00  -0.317   0.7814  
invTTheta2   4.460e+04  2.342e+04   1.905   0.1971  
TXBorn      -1.296e+05  3.982e+04  -3.254   0.0829 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.272 on 2 degrees of freedom
Multiple R-squared:  0.8456,    Adjusted R-squared:  0.6911 
F-statistic: 5.476 on 2 and 2 DF,  p-value: 0.1544

summary(Cplm_high)

Call:
lm(formula = fmla)

Residuals:
     Min       1Q   Median       3Q      Max 
-1137.04  -123.48    31.21   196.93  1420.14 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)      329.4      718.1   0.459    0.655
invTTheta2  -1247884.4  5464255.3  -0.228    0.823
TXBorn       2052194.1  2980374.9   0.689    0.504

Residual standard error: 672.7 on 12 degrees of freedom
Multiple R-squared:  0.04475,   Adjusted R-squared:  -0.1145 
F-statistic: 0.2811 on 2 and 12 DF,  p-value: 0.7598

---

Pregunta 5

Como mencionamos antes, EOSregress() usa a lm(). Por lo tanto, puedes usar la función summary() para indagar la calidad del ajuste. Compara los modelos Cplm_high y Cplm_low con los modelos que has ajustado antes. ¿Cuál obtiene mejor calidad de ajuste?

Respuesta (escribe a continuación)