reset()
## Fit experimental heat capacities of CH4
## using revised Helgeson-Kirkham-Flowers equations
# Read the data from Hnedkovsky and Wood, 1997
f <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
d <- read.csv(f)
# Use data for CH4
d <- d[d$species == "CH4", ]
d <- d[, -1]
# Convert J to cal and MPa to bar
d$Cp <- convert(d$Cp, "cal")
d$P <- convert(d$P, "bar")
# Specify the terms in the HKF equations
var <- c("invTTheta2", "TXBorn")
# Perform regression, with a temperature limit
EOSlm <- EOSregress(d, var, T.max = 600)
# Calculate the Cp at some temperature and pressure
EOScalc(EOSlm$coefficients, 298.15, 1)
# Get the database values of c1, c2 and omega for CH4(aq)
CH4coeffs <- EOScoeffs("CH4", "Cp")
## Make plots comparing the regressions
## with the accepted EOS parameters of CH4
opar <- par(mfrow = c(2,2))
EOSplot(d, T.max = 600)
title("Cp of CH4(aq), fit to 600 K")
legend("bottomleft", pch = 1, legend = "Hnedkovsky and Wood, 1997")
EOSplot(d, coefficients = CH4coeffs)
title("Cp from EOS parameters in database")
EOSplot(d, T.max = 600, T.plot = 600)
title("Cp fit to 600 K, plot to 600 K")
EOSplot(d, coefficients = CH4coeffs, T.plot = 600)
title("Cp from EOS parameters in database")
par(opar)
# Continuing from above, with user-defined variables
Theta <- 228 # K
invTTTheta3 <- function(T, P) (2*T) / (T-T*Theta) ^ 3
invTX <- function(T, P) 1 / T * water("XBorn", T = T, P = P)[,1]
# Print the calculated values of invTTTheta3
EOSvar("invTTTheta3", d$T, d$P)
# Use invTTTheta and invTX in a regression
var <- c("invTTTheta3", "invTX")
EOSregress(d, var)
# Give them a "label" attribute for use in the legend
attr(invTTTheta3, "label") <-
quote(phantom()%*%2 * italic(T) / (italic(T) - italic(T) * Theta) ^ 3)
attr(invTX, "label") <- quote(phantom() / italic(T * X))
# Uncomment the following to make the plot
#EOSplot(d, var)
## Model experimental volumes of CH4
## using HKF equation and an exploratory one
f <- system.file("extdata/cpetc/HWM96_V.csv", package = "CHNOSZ")
d <- read.csv(f)
# Use data for CH4 near 280 bar
d <- d[d$species == "CH4", ]
d <- d[, -1]
d <- d[abs(d$P - 28) < 0.1, ]
d$P <- convert(d$P, "bar")
# The HKF equation
varHKF <- c("invTTheta", "QBorn")
# alpha is the expansivity coefficient of water
varal <- c("invTTheta", "alpha")
opar <- par(mfrow = c(2, 2))
# For both HKF and the expansivity equation,
# we'll fit up to a temperature limit
EOSplot(d, varHKF, T.max = 663, T.plot = 625)
legend("bottomright", pch = 1, legend = "Hnedkovsky et al., 1996")
title("V of CH4(aq), HKF equation")
EOSplot(d, varal, T.max = 663, T.plot = 625)
title("V of CH4(aq), expansivity equation")
EOSplot(d, varHKF, T.max = 663)
title("V of CH4(aq), HKF equation")
EOSplot(d, varal, T.max = 663)
title("V of CH4(aq), expansivity equation")
par(opar)
# Note that the volume regression using the HKF gives
# a result for omega (coefficient on Q) that is
# not consistent with the high-T heat capacities
Run the code above in your browser using DataLab