next up previous
Next: [2] wjd Up: CHNOSZ examples Previous: [1] anim

[2] EOSregress

## Don't show: data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species
## End(Don't show) ## 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/Cp.CH4.HW97.csv", package="CHNOSZ") d <- read.csv(f) # have to 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) # the result is within 10% of the accepted # values of c1, c2 and omega for CH4(aq) CH4coeffs <- EOScoeffs("CH4", "Cp") dcoeffs <- EOSlm$coefficients - CH4coeffs stopifnot(all(abs(dcoeffs/CH4coeffs) < 0.1)) ## make plots comparing the regressions ## here with the accepted EOS parameters of CH4 par(mfrow=c(2,2)) EOSplot(d, T.max=600)
EOSplot: plotting line for P=279.96 bar
title("Cp of CH4(aq), fit to 600 K") legend("bottomleft", pch=1, legend="Hnedkovsky and Wood, 1997") EOSplot(d, coefficients=CH4coeffs)
EOSplot: plotting line for P=279.96 bar
title("Cp from EOS parameters in database") EOSplot(d, T.max=600, T.plot=600)
EOSplot: plotting line for P=279.96 bar
title("Cp fit to 600 K, plot to 600 K") EOSplot(d, coefficients=CH4coeffs, T.plot=600)
EOSplot: plotting line for P=279.96 bar
title("Cp from EOS parameters in database")

Image EOSregress1

 

# continuing from above, with user-defined variables invTTTheta3 <- function(T, P) (2*T)/(T-T*thermo$opt$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)
[1] -1.850630e-12 -1.225139e-12 -8.502064e-13 -6.243590e-13 -5.151951e-13 -4.401916e-13 -4.062540e-13 [8] -3.902133e-13 -3.855398e-13 -3.818611e-13 -3.795887e-13 -3.763294e-13 -3.706987e-13 -3.605415e-13 [15] -3.452356e-13
# use invTTTheta and invTX in a regression var <- c("invTTTheta3", "invTX") EOSregress(d, var)
Call: lm(formula = fmla) Coefficients: (Intercept) invTTTheta3 invTX -1.224e+03 -9.378e+14 -2.379e+09
# 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/V.CH4.HWM96.csv", package="CHNOSZ") d <- read.csv(f) d$P <- convert(d$P, "bar") # the HKF equation varHKF <- c("invTTheta", "QBorn") # alpha is the expansivity coefficient of water varal <- c("invTTheta", "alpha") 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)
EOSplot: plotting line for P=280.010526315789 bar
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)
EOSplot: plotting line for P=280.010526315789 bar
title("V of CH4(aq), expansivity equation") EOSplot(d, varHKF, T.max=663)
EOSplot: plotting line for P=280.010526315789 bar
title("V of CH4(aq), HKF equation") EOSplot(d, varal, T.max=663)
EOSplot: plotting line for P=280.010526315789 bar
title("V of CH4(aq), expansivity equation") # 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
EOSrgr>

Image EOSregress2

 


next up previous
Next: [2] wjd Up: CHNOSZ examples Previous: [1] anim