[CHNOSZ-commits] r52 - in pkg/CHNOSZ: . R inst inst/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 6 16:36:01 CEST 2013
Author: jedick
Date: 2013-05-06 16:35:56 +0200 (Mon, 06 May 2013)
New Revision: 52
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-IAPWS95.R
pkg/CHNOSZ/inst/tests/test-subcrt.R
Log:
IAPWS95: handle lower temperatures, and improve operation with subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/DESCRIPTION 2013-05-06 14:35:56 UTC (rev 52)
@@ -1,6 +1,6 @@
-Date: 2013-03-28
+Date: 2013-05-06
Package: CHNOSZ
-Version: 1.0.0
+Version: 1.0.0-1
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/R/subcrt.R 2013-05-06 14:35:56 UTC (rev 52)
@@ -363,6 +363,7 @@
# logK
if('logk' %in% prop) {
for(i in 1:length(out)) {
+ # NOTE: the following depends on the water function renaming g to G
out[[i]] <- cbind(out[[i]],data.frame(logK=convert(out[[i]]$G,'logK',T=T)))
colnames(out[[i]][ncol(out[[i]])]) <- 'logK'
}
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/R/water.R 2013-05-06 14:35:56 UTC (rev 52)
@@ -22,9 +22,9 @@
rho.0 <- 1000 # kg m-3
# Equation 1
epsilon.0 <- 8.8541878E-12 # permittivity of vacuum, C^2 J-1 m-1
- epsfun.lhs <- function(e) (e-1)*(2*e+1)/(9*e)
+ #epsfun.lhs <- function(e) (e-1)*(2*e+1)/(9*e)
epsfun.rhs <- function(T,V.m) N.A*(alpha+mufun()/(3*epsilon.0*k*T))/(3*V.m)
- epsfun <- function(e,T,V.m) epsfun.lhs(e) - epsfun.rhs(T,V.m)
+ #epsfun <- function(e,T,V.m) epsfun.lhs(e) - epsfun.rhs(T,V.m)
mufun <- function() gfun()*mu^2
gfun <- function() rhofun()*rho/rho.0 + 1
# Equation 3
@@ -32,7 +32,9 @@
b[4]*(T-215)^-0.5 + b[5]*(T-215)^-0.25 +
exp(b[6]*T^-1 + b[7]*T^-2 + b[8]*P*T^-1 + b[9]*P*T^-2)
epsilon <- function(T,rho) {
- tu <- try(uniroot(epsfun,c(1E-1,1E3),T=T,V.m=M/rho)$root,TRUE)
+ #tu <- try(uniroot(epsfun,c(1E-1,1E3),T=T,V.m=M/rho)$root,TRUE)
+ epspoly <- function() epsfun.rhs(T,V.m=M/rho)
+ tu <- (9*epspoly() + 1 + ((9*epspoly()+1)*(9*epspoly()+1) + 8)^0.5) / 4 #Marc Neveu added 4/24/2013
if(!is.numeric(tu)) {
warning('water.AW90: no root for density at ',T,' K and ',rho,' kg m-3.',call.=FALSE,immediate.=TRUE)
tu <- NA
@@ -71,7 +73,9 @@
} else {
# here we get properties using IAPWS-95
w.out <- water.IAPWS95(property, T, P)
- colnames(w.out) <- property
+ # normalize the names to use upper case (expected by subcrt())
+ iprop <- match(tolower(property), tolower(water.props("IAPWS95")))
+ colnames(w.out) <- water.props("IAPWS95")[iprop]
return(w.out)
}
}
@@ -87,6 +91,7 @@
"V", "rho", "Psat", "E", "kT")
else if(formulation=="IAPWS95")
props <- c("A", "G", "S", "U", "H", "Cv", "Cp",
+ "Speed", "diel",
"YBorn", "QBorn", "XBorn", "NBorn", "UBorn",
"V", "rho", "Psat", "de.dT", "de.dP", "P")
return(props)
@@ -198,11 +203,17 @@
rho <- numeric()
for(i in 1:length(T)) {
if(T[i] < 647.096) {
- rho.lower <- WP02.auxiliary("rho.liquid", T=T[i]) - 2
+ rho.lower <- WP02.auxiliary('rho.liquid',T=T[i])-2
rho.upper <- rho.lower + 400
if(P[i] < 5000) rho.upper <- rho.lower + 300
if(P[i] < 1000) rho.upper <- rho.lower + 200
- if(P[i] < 300) rho.upper <- rho.lower + 30
+ if(P[i] < 300) {
+ rho.upper <- rho.lower + 30
+ if(T[i] < 250) { #Marc Neveu added 4/23/2013
+ rho.lower <- rho.lower - 10
+ rho.upper <- rho.lower + 40
+ }
+ }
} else {
rho.lower <- 0.01
rho.upper <- 1200
@@ -259,10 +270,10 @@
return(convert(IAPWS95('cv',T=T,rho=my.rho)$cv*M,'cal'))
cp <- function()
return(convert(IAPWS95('cp',T=T,rho=my.rho)$cp*M,'cal'))
- w <- function()
+ speed <- function()
return(IAPWS95('w',T=T,rho=my.rho)$w*100) # to cm/s
## electrostatic properties
- epsilon <- function() return(water.AW90(T=T,rho=my.rho,P=convert(P,'MPa')))
+ diel <- function() return(water.AW90(T=T,rho=my.rho,P=convert(P,'MPa')))
de.dt <- function() {
p <- numeric()
for(i in 1:length(T)) {
@@ -378,6 +389,5 @@
if(i > 1) w.out <- cbind(w.out, wnew) else w.out <- wnew
}
msgout("\n")
- names(w.out) <- property
return(w.out)
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/inst/NEWS 2013-05-06 14:35:56 UTC (rev 52)
@@ -1,3 +1,16 @@
+CHANGES IN CHNOSZ 1.0.0-1 (2013-05-06)
+--------------------------------------
+
+- Fix IAPWS-95 calculations for subcrt(): in water.IAPWS95(), rename
+ internal function from 'epsilon' to 'diel'; in water(), return
+ upper-case versions of names of properties.
+
+- rho.IAPWS95() has modified search interval to give results at
+ T < 250 K and P < 300 bar. Thanks to Marc Neveu.
+
+- In water.AW90(), analytical solution to polynomial equation is
+ now used to calculate dielectric constant. Thanks to Marc Neveu.
+
CHANGES IN CHNOSZ 1.0.0 (2013-03-28)
------------------------------------
Modified: pkg/CHNOSZ/inst/tests/test-IAPWS95.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-IAPWS95.R 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/inst/tests/test-IAPWS95.R 2013-05-06 14:35:56 UTC (rev 52)
@@ -75,6 +75,13 @@
expect_equal(round(signif(vapor.calc$cp, 5), 4), cp.vapor.ref, tol=1e-1)
})
+test_that("calculations are possible at low temperatures", {
+ # the sequences start at the lowest whole-number temperature (K)
+ # where the function returns a value of density at the given pressure
+ expect_that(any(is.na(water.IAPWS95("rho", T=seq(234, 274, 3), P=rep(1, 14)))), is_false())
+ expect_that(any(is.na(water.IAPWS95("rho", T=seq(227, 272, 3), P=rep(1000, 16)))), is_false())
+})
+
# reference
# Wagner, W. and Pruss, A. (2002)
Modified: pkg/CHNOSZ/inst/tests/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-subcrt.R 2013-03-28 15:57:07 UTC (rev 51)
+++ pkg/CHNOSZ/inst/tests/test-subcrt.R 2013-05-06 14:35:56 UTC (rev 52)
@@ -87,6 +87,15 @@
expect_equal(s.H2O$T[22], 0.01)
})
+test_that("calculations using IAPWS-95 are possible", {
+ thermo$opt$water <<- "IAPWS95"
+ sb <- subcrt(c("H2O", "Na+"), T=c(-30, -20, 0, 10), P=1)$out
+ # the test is not a demanding numerical comparison, more that we got numbers and no error
+ expect_that(all(sb$`Na+`$G < sb$water$G), is_true())
+ # clean up
+ suppressMessages(data(thermo))
+})
+
# references
# Amend, J. P. and Shock, E. L. (2001)
More information about the CHNOSZ-commits
mailing list