[CHNOSZ-commits] r34 - in pkg/CHNOSZ: . R inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 27 13:23:52 CET 2012
Author: jedick
Date: 2012-12-27 13:23:52 +0100 (Thu, 27 Dec 2012)
New Revision: 34
Added:
pkg/CHNOSZ/inst/tests/test-water.IAPWS95.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/water.Rd
Log:
add some tests of water.IAPWS95()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-12-18 13:25:21 UTC (rev 33)
+++ pkg/CHNOSZ/DESCRIPTION 2012-12-27 12:23:52 UTC (rev 34)
@@ -1,6 +1,6 @@
-Date: 2012-12-18
+Date: 2012-12-27
Package: CHNOSZ
-Version: 0.9.8-8
+Version: 0.9.8-9
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
@@ -17,4 +17,4 @@
License: GPL (>= 2)
BuildResaveData: no
ZipData: no
-URL: http://www.chnosz.net/
+URL: http://www.chnosz.net/, http://chnosz.r-forge.r-project.org/
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2012-12-18 13:25:21 UTC (rev 33)
+++ pkg/CHNOSZ/R/water.R 2012-12-27 12:23:52 UTC (rev 34)
@@ -51,174 +51,174 @@
return(t)
}
-water.IAPWS95 <- function(property,T=298.15,rho=1000) {
- # the IAPWS-95 formulation for pure H2O from
- # Wagner and Pruss, 2002
- property <- tolower(property)
- # triple point
- T.triple <- 273.16 # K
- P.triple <- 611.657 # Pa
- rho.triple.liquid <- 999.793
- rho.triple.vapor <- 0.00485458
- # normal boiling point
- T.boiling <- 373.124
- P.boiling <- 0.101325
- rho.boiling.liquid <- 958.367
- rho.boiling.vapor <- 0.597657
- # critical point constants
- T.critical <- 647.096 # K
- rho.critical <- 322 # kg m-3
- # specific and molar gas constants
- R <- 0.46151805 # kJ kg-1 K-1
- # R.M <- 8.314472 # J mol-1 K-1
- # molar mass
- M <- 18.015268 # g mol-1
-
- ### idealgas
- idealgas <- function(p) {
- # Table 6.1
- n <- c( -8.32044648201, 6.6832105268, 3.00632, 0.012436,
- 0.97315, 1.27950, 0.96956, 0.24873 )
- gamma <- c( NA, NA, NA, 1.28728967,
- 3.53734222, 7.74073708, 9.24437796, 27.5075105 )
- # Equation 6.5
- phi <- function() log(delta) + n[1] + n[2]*tau + n[3]*log(tau) +
- sum( n[4:8] * log(1-exp(-gamma[4:8]*tau)) )
- phi.delta <- function() 1/delta+0+0+0+0
- phi.delta.delta <- function() -1/delta^2+0+0+0+0
- phi.tau <- function() 0+0+n[2]+n[3]/tau+sum(n[4:8]*gamma[4:8]*((1-exp(-gamma[4:8]*tau))^-1-1))
- phi.tau.tau <- function() 0+0+0-n[3]/tau^2-sum(n[4:8]*gamma[4:8]^2 *
- exp(-gamma[4:8]*tau)*(1-exp(-gamma[4:8]*tau))^-2)
- phi.delta.tau <- function() 0+0+0+0+0
- return(get(p)())
- }
-
- #### residual
- residual <- function(p) {
- # Table 6.2
- c <- c(rep(NA,7),rep(1,15),rep(2,20),rep(3,4),4,rep(6,4),rep(NA,5))
- d <- c(1,1,1,2,2,3,4,1,1,1,2,2,3,4,
- 4,5,7,9,10,11,13,15,1,2,2,2,3,4,
- 4,4,5,6,6,7,9,9,9,9,9,10,10,12,
- 3,4,4,5,14,3,6,6,6,3,3,3,NA,NA)
- t <- c(-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4,2,
- 13,9,3,4,11,4,13,1,7,1,9,10,10,3,
- 7,10,10,6,10,10,1,2,3,4,8,6,9,8,
- 16,22,23,23,10,50,44,46,50,0,1,4,NA,NA)
- n <- c( 0.12533547935523E-1, 0.78957634722828E1 ,-0.87803203303561E1 ,
- 0.31802509345418 ,-0.26145533859358 ,-0.78199751687981E-2,
- 0.88089493102134E-2,-0.66856572307965 , 0.20433810950965 ,
- -0.66212605039687E-4,-0.19232721156002 ,-0.25709043003438 ,
- 0.16074868486251 ,-0.40092828925807E-1, 0.39343422603254E-6,
- -0.75941377088144E-5, 0.56250979351888E-3,-0.15608652257135E-4,
- 0.11537996422951E-8, 0.36582165144204E-6,-0.13251180074668E-11,
- -0.62639586912454E-9,-0.10793600908932 , 0.17611491008752E-1,
- 0.22132295167546 ,-0.40247669763528 , 0.58083399985759 ,
- 0.49969146990806E-2,-0.31358700712549E-1,-0.74315929710341 ,
- 0.47807329915480 , 0.20527940895948E-1,-0.13636435110343 ,
- 0.14180634400617E-1, 0.83326504880713E-2,-0.29052336009585E-1,
- 0.38615085574206E-1,-0.20393486513704E-1,-0.16554050063734E-2,
- 0.19955571979541E-2, 0.15870308324157E-3,-0.16388568342530E-4,
- 0.43613615723811E-1, 0.34994005463765E-1,-0.76788197844621E-1,
- 0.22446277332006E-1,-0.62689710414685E-4,-0.55711118565645E-9,
- -0.19905718354408 , 0.31777497330738 ,-0.11841182425981 ,
- -0.31306260323435E2 , 0.31546140237781E2 ,-0.25213154341695E4 ,
- -0.14874640856724 , 0.31806110878444)
- alpha <- c(rep(NA,51),20,20,20,NA,NA)
- beta <- c(rep(NA,51),150,150,250,0.3,0.3)
- gamma <- c(rep(NA,51),1.21,1.21,1.25,NA,NA)
- epsilon <- c(rep(NA,51),1,1,1,NA,NA)
- a <- c(rep(NA,54),3.5,3.5)
- b <- c(rep(NA,54),0.85,0.95)
- B <- c(rep(NA,54),0.2,0.2)
- C <- c(rep(NA,54),28,32)
- D <- c(rep(NA,54),700,800)
- A <- c(rep(NA,54),0.32,0.32)
+idealgas.IAPWS95 <- function(p, delta, tau) {
+ ## the ideal gas part in the IAPWS-95 formulation
+ # from Table 6.1 of Wagner and Pruss, 2002
+ n <- c( -8.32044648201, 6.6832105268, 3.00632, 0.012436,
+ 0.97315, 1.27950, 0.96956, 0.24873 )
+ gamma <- c( NA, NA, NA, 1.28728967,
+ 3.53734222, 7.74073708, 9.24437796, 27.5075105 )
+ # Equation 6.5
+ phi <- function() log(delta) + n[1] + n[2]*tau + n[3]*log(tau) +
+ sum( n[4:8] * log(1-exp(-gamma[4:8]*tau)) )
+ # derivatives from Table 6.4
+ phi.delta <- function() 1/delta+0+0+0+0
+ phi.delta.delta <- function() -1/delta^2+0+0+0+0
+ phi.tau <- function() 0+0+n[2]+n[3]/tau+sum(n[4:8]*gamma[4:8]*((1-exp(-gamma[4:8]*tau))^-1-1))
+ phi.tau.tau <- function() 0+0+0-n[3]/tau^2-sum(n[4:8]*gamma[4:8]^2 *
+ exp(-gamma[4:8]*tau)*(1-exp(-gamma[4:8]*tau))^-2)
+ phi.delta.tau <- function() 0+0+0+0+0
+ return(get(p)())
+}
-
- # Table 6.5
- i1 <- 1:7
- i2 <- 8:51
- i3 <- 52:54
- i4 <- 55:56
- # deriviatives of distance function
- Delta <- function(i) { Theta(i)^2 + B[i] * ((delta-1)^2)^a[i] }
- Theta <- function(i) { (1-tau) + A[i] * ((delta-1)^2)^(1/(2*beta[i])) }
- Psi <- function(i) { exp ( -C[i]*(delta-1)^2 - D[i]*(tau-1)^2 ) }
- dDelta.bi.ddelta <- function(i) { b[i]*Delta(i)^(b[i]-1)*dDelta.ddelta(i) }
- d2Delta.bi.ddelta2 <- function(i) { b[i]*( Delta(i)^(b[i]-1) * d2Delta.ddelta2(i) +
- (b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i)^2 ) }
- dDelta.bi.dtau <- function(i) { -2*Theta(i)*b[i]*Delta(i)^(b[i]-1) }
- d2Delta.bi.dtau2 <- function(i) { 2*b[i]*Delta(i)^(b[i]-1) + 4*Theta(i)^2*b[i]*(b[i]-1)*Delta(i)^(b[i]-2) }
- d2Delta.bi.ddelta.dtau <- function(i) { -A[i]*b[i]*2/beta[i]*Delta(i)^(b[i]-1)*(delta-1) *
- ((delta-1)^2)^(1/(2*beta[i])-1) - 2*Theta(i)*b[i]*(b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i) }
- dDelta.ddelta <- function(i) { (delta-1) * ( A[i]*Theta(i)*2/beta[i]*((delta-1)^2)^(1/(2*beta[i])-1) +
- 2*B[i]*a[i]*((delta-1)^2)^(a[i]-1) ) }
- d2Delta.ddelta2 <- function(i) { 1/(delta-1)*dDelta.ddelta(i) + (delta-1)^2 * (
- 4*B[i]*a[i]*(a[i]-1)*((delta-1)^2)^(a[i]-2) + 2*A[i]^2*(1/beta[i])^2 *
- (((delta-1)^2)^(1/(2*B[i])-1))^2 + A[i]*Theta(i)*4/beta[i]*(1/(2*B[i])-1) *
- ((delta-1)^2)^(1/(2*beta[i])-2) ) }
- # derivatives of exponential function
- dPsi.ddelta <- function(i) { -2*C[i]*(delta-1)*Psi(i) }
- d2Psi.ddelta2 <- function(i) { ( 2*C[i]*(delta-1)^2 - 1 ) * 2*C[i]*Psi(i) }
- dPsi.dtau <- function(i) { -2*D[i]*(tau-1)*Psi(i) }
- d2Psi.dtau2 <- function(i) { (2*D[i]*(tau-1)^2 - 1) * 2*D[i]*Psi(i) }
- d2Psi.ddelta.dtau <- function(i) { 4*C[i]*D[i]*(delta-1)*(tau-1)*Psi(i) }
- # dimensionless Helmholtz free energy and derivatives
- phi <- function() {
- sum(n[i1]*delta^d[i1]*tau^t[i1]) +
- sum(n[i2]*delta^d[i2]*tau^t[i2]*exp(-delta^c[i2])) +
- sum(n[i3]*delta^d[i3]*tau^t[i3] *
- exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) ) +
- sum(n[i4]*Delta(i4)^b[i4]*delta*Psi(i4))
- }
- phi.delta <- function() {
- sum(n[i1]*d[i1]*delta^(d[i1]-1)*tau^t[i1]) +
- sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-1)*tau^t[i2]*(d[i2]-c[i2]*delta^c[i2]))) +
- sum(n[i3]*delta^d[i3]*tau^t[i3] *
- exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) *
- (d[i3]/delta - 2 * alpha[i3]*(delta-epsilon[i3])) ) +
- sum(n[i4] * ( Delta(i4)^b[i4] * (Psi(i4)+delta*dPsi.ddelta(i4)) + dDelta.bi.ddelta(i4)*delta*Psi(i4) ) )
- }
- phi.delta.delta <- function() {
- sum(n[i1]*d[i1]*(d[i1]-1)*delta^(d[i1]-2)*tau^t[i1]) +
- sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-2)*tau^t[i2]*((d[i2]-c[i2]*delta^c[i2]) *
- (d[i2]-1-c[i2]*delta^c[i2])-c[i2]^2*delta^c[i2]))) +
- sum(n[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2) * (
- -2*alpha[i3]*delta^d[i3]+4*alpha[i3]^2*delta^d[i3]*(delta-epsilon[i3])^2 -
- 4*d[i3]*alpha[i3]*delta^(d[i3]-1)*(delta-epsilon[i3])+d[i3]*(d[i3]-1)*delta^(d[i3]-2) ) ) +
- sum(n[i4]*( Delta(i4)^b[i4]*(2*dPsi.ddelta(i4)+delta*d2Psi.ddelta2(i4)) +
- 2*dDelta.bi.ddelta(i4)*(Psi(i4)+delta*dPsi.ddelta(i4)) + d2Delta.bi.ddelta2(i4)*delta*Psi(i4) ) )
- }
- phi.tau <- function() {
- sum(n[i1]*t[i1]*delta^d[i1]*tau^(t[i1]-1)) +
- sum(n[i2]*t[i2]*delta^d[i2]*tau^(t[i2]-1)*exp(-delta^c[i2])) +
- sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
- (t[i3]/tau-2*beta[i3]*(tau-gamma[i3]))) +
- sum(n[i4]*delta*(dDelta.bi.dtau(i4)*Psi(i4)+Delta(i4)^b[i4]*dPsi.dtau(i4)))
- }
- phi.tau.tau <- function() {
- sum(n[i1]*t[i1]*(t[i1]-1)*delta^d[i1]*tau^(t[i1]-2)) +
- sum(n[i2]*t[i2]*(t[i2]-1)*delta^d[i2]*tau^(t[i2]-2)*exp(-delta^c[i2])) +
- sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
- (((t[i3]/tau)-2*beta[i3]*(tau-gamma[i3]))^2-t[i3]/tau^2-2*beta[i3])) +
- sum(n[i4]*delta*(d2Delta.bi.dtau2(i4)*Psi(i4)+2*dDelta.bi.dtau(i4)*dPsi.dtau(i4) +
- Delta(i4)^b[i4]*d2Psi.dtau2(i4)))
- }
- phi.delta.tau <- function() {
- sum(n[i1]*d[i1]*t[i1]*delta^(d[i1]-1)*tau^(t[i1]-1)) +
- sum(n[i2]*t[i2]*delta^(d[i2]-1)*tau^(t[i2]-1)*(d[i2]-c[i2]*delta^c[i2])*exp(-delta^c[i2])) +
- sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
- ((d[i3]/delta)-2*alpha[i3]*(delta-epsilon[i3]))*(t[i3]/tau-2*beta[i3]*(tau-gamma[i3])) ) +
- sum(n[i4]*(Delta(i4)^b[i4]*(dPsi.dtau(i4)+delta*d2Psi.ddelta.dtau(i4)) +
- delta*dDelta.bi.ddelta(i4)*dPsi.dtau(i4)+dDelta.bi.dtau(i4) * (Psi(i4)+delta*dPsi.ddelta(i4)) +
- d2Delta.bi.ddelta.dtau(i4)*delta*Psi(i4) ))
- }
-
- return(get(p)())
+residual.IAPWS95 <- function(p, delta, tau) {
+ ## the residual part in the IAPWS-95 formulation
+ # from Table 6.2 of Wagner and Pruss, 2002
+ c <- c(rep(NA,7),rep(1,15),rep(2,20),rep(3,4),4,rep(6,4),rep(NA,5))
+ d <- c(1,1,1,2,2,3,4,1,1,1,2,2,3,4,
+ 4,5,7,9,10,11,13,15,1,2,2,2,3,4,
+ 4,4,5,6,6,7,9,9,9,9,9,10,10,12,
+ 3,4,4,5,14,3,6,6,6,3,3,3,NA,NA)
+ t <- c(-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4,2,
+ 13,9,3,4,11,4,13,1,7,1,9,10,10,3,
+ 7,10,10,6,10,10,1,2,3,4,8,6,9,8,
+ 16,22,23,23,10,50,44,46,50,0,1,4,NA,NA)
+ n <- c( 0.12533547935523E-1, 0.78957634722828E1 ,-0.87803203303561E1 ,
+ 0.31802509345418 ,-0.26145533859358 ,-0.78199751687981E-2,
+ 0.88089493102134E-2,-0.66856572307965 , 0.20433810950965 ,
+ -0.66212605039687E-4,-0.19232721156002 ,-0.25709043003438 ,
+ 0.16074868486251 ,-0.40092828925807E-1, 0.39343422603254E-6,
+ -0.75941377088144E-5, 0.56250979351888E-3,-0.15608652257135E-4,
+ 0.11537996422951E-8, 0.36582165144204E-6,-0.13251180074668E-11,
+ -0.62639586912454E-9,-0.10793600908932 , 0.17611491008752E-1,
+ 0.22132295167546 ,-0.40247669763528 , 0.58083399985759 ,
+ 0.49969146990806E-2,-0.31358700712549E-1,-0.74315929710341 ,
+ 0.47807329915480 , 0.20527940895948E-1,-0.13636435110343 ,
+ 0.14180634400617E-1, 0.83326504880713E-2,-0.29052336009585E-1,
+ 0.38615085574206E-1,-0.20393486513704E-1,-0.16554050063734E-2,
+ 0.19955571979541E-2, 0.15870308324157E-3,-0.16388568342530E-4,
+ 0.43613615723811E-1, 0.34994005463765E-1,-0.76788197844621E-1,
+ 0.22446277332006E-1,-0.62689710414685E-4,-0.55711118565645E-9,
+ -0.19905718354408 , 0.31777497330738 ,-0.11841182425981 ,
+ -0.31306260323435E2 , 0.31546140237781E2 ,-0.25213154341695E4 ,
+ -0.14874640856724 , 0.31806110878444)
+ alpha <- c(rep(NA,51),20,20,20,NA,NA)
+ beta <- c(rep(NA,51),150,150,250,0.3,0.3)
+ gamma <- c(rep(NA,51),1.21,1.21,1.25,NA,NA)
+ epsilon <- c(rep(NA,51),1,1,1,NA,NA)
+ a <- c(rep(NA,54),3.5,3.5)
+ b <- c(rep(NA,54),0.85,0.95)
+ B <- c(rep(NA,54),0.2,0.2)
+ C <- c(rep(NA,54),28,32)
+ D <- c(rep(NA,54),700,800)
+ A <- c(rep(NA,54),0.32,0.32)
+ # from Table 6.5
+ i1 <- 1:7
+ i2 <- 8:51
+ i3 <- 52:54
+ i4 <- 55:56
+ # deriviatives of distance function
+ Delta <- function(i) { Theta(i)^2 + B[i] * ((delta-1)^2)^a[i] }
+ Theta <- function(i) { (1-tau) + A[i] * ((delta-1)^2)^(1/(2*beta[i])) }
+ Psi <- function(i) { exp ( -C[i]*(delta-1)^2 - D[i]*(tau-1)^2 ) }
+ dDelta.bi.ddelta <- function(i) { b[i]*Delta(i)^(b[i]-1)*dDelta.ddelta(i) }
+ d2Delta.bi.ddelta2 <- function(i) { b[i]*( Delta(i)^(b[i]-1) * d2Delta.ddelta2(i) +
+ (b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i)^2 ) }
+ dDelta.bi.dtau <- function(i) { -2*Theta(i)*b[i]*Delta(i)^(b[i]-1) }
+ d2Delta.bi.dtau2 <- function(i) { 2*b[i]*Delta(i)^(b[i]-1) + 4*Theta(i)^2*b[i]*(b[i]-1)*Delta(i)^(b[i]-2) }
+ d2Delta.bi.ddelta.dtau <- function(i) { -A[i]*b[i]*2/beta[i]*Delta(i)^(b[i]-1)*(delta-1) *
+ ((delta-1)^2)^(1/(2*beta[i])-1) - 2*Theta(i)*b[i]*(b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i) }
+ dDelta.ddelta <- function(i) { (delta-1) * ( A[i]*Theta(i)*2/beta[i]*((delta-1)^2)^(1/(2*beta[i])-1) +
+ 2*B[i]*a[i]*((delta-1)^2)^(a[i]-1) ) }
+ d2Delta.ddelta2 <- function(i) { 1/(delta-1)*dDelta.ddelta(i) + (delta-1)^2 * (
+ 4*B[i]*a[i]*(a[i]-1)*((delta-1)^2)^(a[i]-2) + 2*A[i]^2*(1/beta[i])^2 *
+ (((delta-1)^2)^(1/(2*B[i])-1))^2 + A[i]*Theta(i)*4/beta[i]*(1/(2*B[i])-1) *
+ ((delta-1)^2)^(1/(2*beta[i])-2) ) }
+ # derivatives of exponential function
+ dPsi.ddelta <- function(i) { -2*C[i]*(delta-1)*Psi(i) }
+ d2Psi.ddelta2 <- function(i) { ( 2*C[i]*(delta-1)^2 - 1 ) * 2*C[i]*Psi(i) }
+ dPsi.dtau <- function(i) { -2*D[i]*(tau-1)*Psi(i) }
+ d2Psi.dtau2 <- function(i) { (2*D[i]*(tau-1)^2 - 1) * 2*D[i]*Psi(i) }
+ d2Psi.ddelta.dtau <- function(i) { 4*C[i]*D[i]*(delta-1)*(tau-1)*Psi(i) }
+ # dimensionless Helmholtz free energy and derivatives
+ phi <- function() {
+ sum(n[i1]*delta^d[i1]*tau^t[i1]) +
+ sum(n[i2]*delta^d[i2]*tau^t[i2]*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3] *
+ exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) ) +
+ sum(n[i4]*Delta(i4)^b[i4]*delta*Psi(i4))
}
-
- # relation of thermodynamic properties to Helmholtz free energy
+ phi.delta <- function() {
+ sum(n[i1]*d[i1]*delta^(d[i1]-1)*tau^t[i1]) +
+ sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-1)*tau^t[i2]*(d[i2]-c[i2]*delta^c[i2]))) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3] *
+ exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) *
+ (d[i3]/delta - 2 * alpha[i3]*(delta-epsilon[i3])) ) +
+ sum(n[i4] * ( Delta(i4)^b[i4] * (Psi(i4)+delta*dPsi.ddelta(i4)) + dDelta.bi.ddelta(i4)*delta*Psi(i4) ) )
+ }
+ phi.delta.delta <- function() {
+ sum(n[i1]*d[i1]*(d[i1]-1)*delta^(d[i1]-2)*tau^t[i1]) +
+ sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-2)*tau^t[i2]*((d[i2]-c[i2]*delta^c[i2]) *
+ (d[i2]-1-c[i2]*delta^c[i2])-c[i2]^2*delta^c[i2]))) +
+ sum(n[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2) * (
+ -2*alpha[i3]*delta^d[i3]+4*alpha[i3]^2*delta^d[i3]*(delta-epsilon[i3])^2 -
+ 4*d[i3]*alpha[i3]*delta^(d[i3]-1)*(delta-epsilon[i3])+d[i3]*(d[i3]-1)*delta^(d[i3]-2) ) ) +
+ sum(n[i4]*( Delta(i4)^b[i4]*(2*dPsi.ddelta(i4)+delta*d2Psi.ddelta2(i4)) +
+ 2*dDelta.bi.ddelta(i4)*(Psi(i4)+delta*dPsi.ddelta(i4)) + d2Delta.bi.ddelta2(i4)*delta*Psi(i4) ) )
+ }
+ phi.tau <- function() {
+ sum(n[i1]*t[i1]*delta^d[i1]*tau^(t[i1]-1)) +
+ sum(n[i2]*t[i2]*delta^d[i2]*tau^(t[i2]-1)*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ (t[i3]/tau-2*beta[i3]*(tau-gamma[i3]))) +
+ sum(n[i4]*delta*(dDelta.bi.dtau(i4)*Psi(i4)+Delta(i4)^b[i4]*dPsi.dtau(i4)))
+ }
+ phi.tau.tau <- function() {
+ sum(n[i1]*t[i1]*(t[i1]-1)*delta^d[i1]*tau^(t[i1]-2)) +
+ sum(n[i2]*t[i2]*(t[i2]-1)*delta^d[i2]*tau^(t[i2]-2)*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ (((t[i3]/tau)-2*beta[i3]*(tau-gamma[i3]))^2-t[i3]/tau^2-2*beta[i3])) +
+ sum(n[i4]*delta*(d2Delta.bi.dtau2(i4)*Psi(i4)+2*dDelta.bi.dtau(i4)*dPsi.dtau(i4) +
+ Delta(i4)^b[i4]*d2Psi.dtau2(i4)))
+ }
+ phi.delta.tau <- function() {
+ sum(n[i1]*d[i1]*t[i1]*delta^(d[i1]-1)*tau^(t[i1]-1)) +
+ sum(n[i2]*t[i2]*delta^(d[i2]-1)*tau^(t[i2]-1)*(d[i2]-c[i2]*delta^c[i2])*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ ((d[i3]/delta)-2*alpha[i3]*(delta-epsilon[i3]))*(t[i3]/tau-2*beta[i3]*(tau-gamma[i3])) ) +
+ sum(n[i4]*(Delta(i4)^b[i4]*(dPsi.dtau(i4)+delta*d2Psi.ddelta.dtau(i4)) +
+ delta*dDelta.bi.ddelta(i4)*dPsi.dtau(i4)+dDelta.bi.dtau(i4) * (Psi(i4)+delta*dPsi.ddelta(i4)) +
+ d2Delta.bi.ddelta.dtau(i4)*delta*Psi(i4) ))
+ }
+ return(get(p)())
+}
+
+water.IAPWS95 <- function(property,T=298.15,rho=1000) {
+ ## the IAPWS-95 formulation for pure H2O from
+ ## Wagner and Pruss, 2002
+ property <- tolower(property)
+ # triple point
+ T.triple <- 273.16 # K
+ P.triple <- 611.657 # Pa
+ rho.triple.liquid <- 999.793
+ rho.triple.vapor <- 0.00485458
+ # normal boiling point
+ T.boiling <- 373.124
+ P.boiling <- 0.101325
+ rho.boiling.liquid <- 958.367
+ rho.boiling.vapor <- 0.597657
+ # critical point constants
+ T.critical <- 647.096 # K
+ rho.critical <- 322 # kg m-3
+ # specific and molar gas constants
+ R <- 0.46151805 # kJ kg-1 K-1
+ # R.M <- 8.314472 # J mol-1 K-1
+ # molar mass
+ M <- 18.015268 # g mol-1
+ ## define functions idealgas and residual, supplying arguments delta and tau
+ idealgas <- function(p) idealgas.IAPWS95(p, delta, tau)
+ residual <- function(p) residual.IAPWS95(p, delta, tau)
+ ## relation of thermodynamic properties to Helmholtz free energy
a <- function() {
x <- idealgas('phi')+residual('phi')
return(x*R*T)
@@ -275,69 +275,22 @@
(idealgas('phi.tau.tau')+residual('phi.tau.tau'))*(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) )
return(x/(R*rho))
}
-
- if(property[1]=='test') {
- # Table 6.6
- pr <- c('phi','phi.delta','phi.delta.delta','phi.tau','phi.tau.tau','phi.delta.tau')
- i1 <- r1 <- i2 <- r2 <- numeric()
- print('expt 1: T=500,rho=838.025')
- T <- 500; rho <- 838.025
- #T <- 354.467; rho <- 970.942
- delta <- rho / rho.critical
- tau <- T.critical / T
- for(i in 1:length(pr)) i1 <- c(i1,idealgas(pr[i]))
- for(i in 1:length(pr)) r1 <- c(r1,residual(pr[i]))
- print('expt 2: T=647,rho=358')
- T <- 647; rho <- 358
- #T <- 354.467; rho <- 0.30864
- delta <- rho / rho.critical
- tau <- T.critical / T
- for(i in 1:length(pr)) i2 <- c(i2,idealgas(pr[i]))
- for(i in 1:length(pr)) r2 <- c(r2,residual(pr[i]))
- t <- data.frame(idealgas.1=i1,residual.1=r1,idealgas.2=i2,residual.2=r2)
- rownames(t) <- pr
- print(t)
- # properties at triple, boiling, and critical points
- tl <- tv <- bl <- bv <- cr <- numeric()
- pr <- c('p','s','u','h','g','cv','cp','mu')
- T <- T.triple; rho <- rho.triple.liquid
- delta <- rho / rho.critical; tau <- T.critical / T
- for(i in 1:length(pr)) tl <- c(tl,get(pr[i])())
- T <- T.triple; rho <- rho.triple.vapor
- delta <- rho / rho.critical; tau <- T.critical / T
- for(i in 1:length(pr)) tv <- c(tv,get(pr[i])())
- T <- T.boiling; rho <- rho.boiling.liquid
- delta <- rho / rho.critical; tau <- T.critical / T
- for(i in 1:length(pr)) bl <- c(bl,get(pr[i])())
- T <- T.boiling; rho <- rho.boiling.vapor
- delta <- rho / rho.critical; tau <- T.critical / T
- for(i in 1:length(pr)) bv <- c(bv,get(pr[i])())
- T <- T.critical; rho <- rho.critical
- delta <- rho / rho.critical; tau <- T.critical / T
- for(i in 1:length(pr)) cr <- c(cr,get(pr[i])())
- T <- c(T.triple,T.triple,T.boiling,T.boiling,T.critical)
- rho <- c(rho.triple.liquid,rho.triple.vapor,rho.boiling.liquid,rho.boiling.vapor,rho.critical)
- t <- as.data.frame(matrix(c(tl,tv,bl,bv,cr),byrow=TRUE,nrow=5))
- t <- data.frame(T=T,rho=rho,t)
- rownames(t) <- c('triple.liquid','triple.vapor','boiling.liquid','boiling.vapor','critical')
- colnames(t) <- c('T','rho',pr)
- print(t)
- return()
- } else {
- ww <- NULL
- my.T <- T; my.rho <- rho
- for(j in 1:length(property)) {
- t <- numeric()
- for(i in 1:length(my.T)) {
- T <- my.T[i]; rho <- my.rho[i]
- # Equation 6.4
- delta <- rho / rho.critical
- tau <- T.critical / T
- t <- c(t,get(property[j])())
- }
- t <- data.frame(t)
- if(j==1) ww <- t else ww <- cbind(ww,t)
+ ## run the calculations
+ ww <- NULL
+ my.T <- T
+ my.rho <- rho
+ for(j in 1:length(property)) {
+ t <- numeric()
+ for(i in 1:length(my.T)) {
+ T <- my.T[i]
+ rho <- my.rho[i]
+ # Equation 6.4
+ delta <- rho / rho.critical
+ tau <- T.critical / T
+ t <- c(t,get(property[j])())
}
+ t <- data.frame(t)
+ if(j==1) ww <- t else ww <- cbind(ww,t)
}
colnames(ww) <- property
return(ww)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2012-12-18 13:25:21 UTC (rev 33)
+++ pkg/CHNOSZ/inst/NEWS 2012-12-27 12:23:52 UTC (rev 34)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9.8-8 (2012-12-18)
+CHANGES IN CHNOSZ 0.9.8-9 (2012-12-27)
--------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
@@ -181,6 +181,9 @@
their own definitions, with an attribute indicating whether the
function is minimized or maximized.
+- Separate idealgas.IAPWS95() and residual.IAPWS95() from
+ water.IAPWS95() (in order to write some test_that tests.)
+
DATA UPDATES:
- Add 148 liquid and 148 crystalline acyclic isoprenoids, polycyclic
Added: pkg/CHNOSZ/inst/tests/test-water.IAPWS95.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-water.IAPWS95.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-water.IAPWS95.R 2012-12-27 12:23:52 UTC (rev 34)
@@ -0,0 +1,83 @@
+context("water.IAPWS95")
+
+test_that("calculations of Helmholtz free energy and its derivatives are consistent with reference cases", {
+ ## reference values of these terms are listed Table 6.6 of Wagner and Pruss, 2002
+ p <- c("phi", "phi.delta", "phi.delta.delta", "phi.tau", "phi.tau.tau", "phi.delta.tau")
+ # reference values for case 1: T=500 K, rho=838.025 kg m-3
+ idealgas.ref.1 <- c(0.204797734e1, 0.384236747, -0.147637878, 0.904611106e1, -0.193249185e1, 0)
+ residual.ref.1 <- c(-0.342693206e1, -0.364366650, 0.856063701, -0.581403435e1, -0.223440737e1, -0.112176915e1)
+ # reference values for case 1: T=647 K, rho=358 kg m-3
+ idealgas.ref.2 <- c(-0.156319605e1, 0.899441341, -0.808994726, 0.980343918e1, -0.343316334e1, 0)
+ residual.ref.2 <- c(-0.121202657e1, -0.714012024, 0.475730696, -0.321722501e1, -0.996029507e1, -0.133214720e1)
+ ## set up the problem
+ # critical point constants
+ T.critical <- 647.096 # K
+ rho.critical <- 322 # kg m-3
+ # T and rho for cases 1 and 2
+ T <- c(500, 647)
+ rho <- c(838.025, 358)
+ # delta and tau for cases 1 and 2
+ delta <- rho / rho.critical
+ tau <- T.critical / T
+ # calculated ideal gas and residual parts for case 1
+ idealgas.calc.1 <- sapply(p, idealgas.IAPWS95, delta[1], tau[1])
+ residual.calc.1 <- sapply(p, residual.IAPWS95, delta[1], tau[1])
+ # calculated ideal gas and residual parts for case 2
+ idealgas.calc.2 <- sapply(p, idealgas.IAPWS95, delta[2], tau[2])
+ residual.calc.2 <- sapply(p, residual.IAPWS95, delta[2], tau[2])
+ ## perform the tests
+ # we almost get away without increasing the tolerance in any test ...
+ expect_equal(idealgas.calc.1, idealgas.ref.1, check.attributes=FALSE)
+ expect_equal(residual.calc.1, residual.ref.1, check.attributes=FALSE)
+ expect_equal(idealgas.calc.2, idealgas.ref.2, check.attributes=FALSE)
+ # ... however an offset is apparent in the value of the residual phi.delta.delta for case 2
+ expect_equal(residual.calc.2, residual.ref.2, check.attributes=FALSE, tol=1e-5)
+})
+
+test_that("calculations of thermodynamic properties are consistent with reference values", {
+ ## these are the properties we test - from Table 13.1 of Wagner and Pruss, 2002
+ ## (speed of sound omitted as it's not currently implemented)
+ p <- c("P", "H", "S", "Cv", "Cp")
+ ## a selection of T and rho at vapor-liquid boundary
+ ## (NOTE: excluding triple and critical points; we have some unresolved issues there)
+ T <- c(274, 320, 368, 416, 464, 512, 560, 608, 647)
+ rho.liquid <- c(999.843, 989.387, 961.984, 923.577, 875.125, 814.982, 737.831, 626.74, 357.34)
+ rho.vapor <- c(0.00514, 0.07166, 0.50231, 2.1203, 6.5107, 16.409, 37.147, 84.173, 286.51)
+ ## reference values
+ P.ref <- c(0.000650, 0.010546, 0.084142, 0.39166, 1.2788, 3.2798, 7.1062, 13.681, 22.038)
+ H.liquid.ref <- c(3.544, 196.170, 397.457, 601.396, 811.225, 1032.06, 1273.11, 1558.42, 2029.44)
+ H.vapor.ref <- c(2502.46, 2585.71, 2667.37, 2737.09, 2785.91, 2803.05, 2771.24, 2646.01, 2148.56)
+ S.liquid.ref <- c(0.0130, 0.6629, 1.2487, 1.7687, 2.2436, 2.6914, 3.1319, 3.6034, 4.3224)
+ S.vapor.ref <- c(9.1331, 8.1302, 7.4169, 6.9025, 6.4994, 6.1504, 5.8071, 5.3922, 4.5065)
+ cv.liquid.ref <- c(4.2156, 4.0416, 3.7950, 3.5560, 3.3523, 3.1884, 3.0722, 3.0617, 6.2344)
+ cv.vapor.ref <- c(1.4191, 1.4627, 1.5428, 1.7138, 2.0015, 2.3697, 2.8225, 3.4538, 6.2740)
+ cp.liquid.ref <- c(4.2171, 4.1807, 4.2101, 4.2892, 4.4513, 4.7616, 5.4239, 7.6169, 3905.2)
+ cp.vapor.ref <- c(1.8852, 1.9417, 2.0601, 2.3334, 2.8561, 3.7263, 5.4099, 10.805, 5334.1)
+ ## calculated values
+ liquid.calc <- water.IAPWS95(p, T, rho.liquid)
+ vapor.calc <- water.IAPWS95(p, T, rho.vapor)
+ ## the tests
+ # take P to 5 significant digits but not more than 6 decimals
+ expect_equal(round(signif(liquid.calc$p, 5), 6), P.ref, tol=1e-3)
+ expect_equal(round(signif(vapor.calc$p, 5), 6), P.ref, tol=1e-4)
+ # take H to 6 significant digits but not more than 3 decimals
+ expect_equal(round(signif(liquid.calc$h, 6), 3), H.liquid.ref, tol=1e-5)
+ expect_that(round(signif(vapor.calc$h, 6), 3), equals(H.vapor.ref)) # spot on!
+ # round S to 4 decimals
+ expect_that(round(liquid.calc$s, 4), equals(S.liquid.ref)) # spot on!
+ expect_equal(round(vapor.calc$s, 4), S.vapor.ref, tol=1e-4)
+ # round cv to 4 decimals
+ expect_equal(round(liquid.calc$cv, 4), cv.liquid.ref, tol=1e-4)
+ expect_equal(round(vapor.calc$cv, 4), cv.vapor.ref, tol=1e-4)
+ # take cp to 5 significant digits but not more than 4 decimals
+ # note high tolerance setting: the highest temperature is the challenge
+ expect_equal(round(signif(liquid.calc$cp, 5), 4), cp.liquid.ref, tol=1e0)
+ expect_equal(round(signif(vapor.calc$cp, 5), 4), cp.vapor.ref, tol=1e-1)
+})
+
+# reference
+
+# Wagner, W. and Pruss, A. (2002)
+# The IAPWS formulation 1995 for the thermodynamic properties of
+# ordinary water substance for general and scientific use
+# J. Phys. Chem. Ref. Data 31, 387--535. http://dx.doi.org/10.1063/1.1461829
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2012-12-18 13:25:21 UTC (rev 33)
+++ pkg/CHNOSZ/man/water.Rd 2012-12-27 12:23:52 UTC (rev 34)
@@ -3,6 +3,8 @@
\alias{water}
\alias{water.AW90}
\alias{water.IAPWS95}
+\alias{idealgas.IAPWS95}
+\alias{residual.IAPWS95}
\alias{water.WP02}
\alias{water.SUPCRT92}
\title{Properties of Water}
@@ -14,16 +16,21 @@
water(property = NULL, T = thermo$opt$Tr, P = "Psat")
water.SUPCRT92(property, T = 298.15, P = 1, isat = 0)
water.IAPWS95(property, T = 298.15, rho = 1000)
+ idealgas.IAPWS95(p, delta, tau)
+ residual.IAPWS95(p, delta, tau)
water.WP02(property, T = 298.15)
water.AW90(T = 298.15, rho = 1000, P = 0.1)
}
\arguments{
- \item{property}{character, name(s) of property(s) to calculate.}
- \item{T}{numeric, temperature (K).}
- \item{P}{character or numeric, \samp{Psat}; or pressure (bar for \code{water}, \code{water.SUPCRT92}; MPa for \code{water.AW90}, \code{water.IAPWS95}, \code{water.WP02}).}
+ \item{property}{character, name(s) of property(s) to calculate}
+ \item{T}{numeric, temperature (K)}
+ \item{P}{character or numeric, \samp{Psat}; or pressure (bar for \code{water}, \code{water.SUPCRT92}; MPa for \code{water.AW90}, \code{water.IAPWS95}, \code{water.WP02})}
\item{isat}{numeric, if \eqn{1}, calculate values of \samp{Psat}}
- \item{rho}{numeric, density (kg m\eqn{^{-3}}{^-3}).}
+ \item{rho}{numeric, density (kg m\eqn{^{-3}}{^-3})}
+ \item{p}{character, name of property (Helmholtz free energy or its derivatives)}
+ \item{delta}{numeric, density divided by critical density}
+ \item{tau}{numeric, critical temperature divided by temperature}
}
\details{
@@ -76,10 +83,13 @@
\code{water.SUPCRT92} interfaces to the FORTRAN subroutine taken from the \acronym{SUPCRT92} package (H2O92D.F) for calculating properties of water. These calculations are based on data and equations of Levelt-Sengers et al., 1983, Haar et al., 1984, and Johnson and Norton, 1991, among others (see Johnson et al., 1992). If \code{isat} is equal to \eqn{1} (or \code{TRUE}), the values of \code{P} are ignored and values of \samp{Psat} are returned. \samp{Psat} refers to one bar below 100 \eqn{^{\circ}}{°}C, otherwise to the vapor-liquid saturation pressure at temperatures below the critical point (\samp{Psat} is not available at temperatures above the critical point). \code{water.SUPCRT92} function provides a limited interface to the FORTRAN subroutine; some functions provided there are not made available here (e.g., using variable density instead of pressure, or calculating the properties of steam). The properties of steam in CHNOSZ, as in \acronym{SUPCRT92}, are calculated using general equations for crystalline, gaseous and liquid species (\code{\link{cgl}}). The IAPWS-95 formulation also has provisions for computing the properties of steam, but these are currently not used by CHNOSZ.
- \code{water.IAPWS95} provides an implementation of the IAPWS-95 formulation for properties (including pressure) calculated as a function of temperature and density. To compute the thermodynamic and electrostatic properties of water as a function of temperature and pressure using \code{water.IAPWS95}, \code{water} applies a root-finding function (\code{\link{uniroot}}) to determine the corresponding values of density. Electrostatic properties in this case are derived from values of the static dielectric constant (\code{epsilon}) calculated using equations given by Archer and Wang, 1990 and coded in \code{water.AW90}. Note that the \code{water.AW90} computes the static dielectric constant at given temperatures and pressures, so \code{water} contains routines for calculating its derivatives with respect to temperature and pressure. A keyword, \samp{test}, may be given as \code{property} to \code{water.IAPWS95}, which causes the printing of two tables, one representing the ideal-gas and residual contributions to the Helmholtz free energy (Table 6.6 of Wagner and Pruss, 2002), and a second with a selection of calculated properties for the liquid and vapor at the triple and boiling points.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 34
More information about the CHNOSZ-commits
mailing list