Next: subcrt
Up: CHNOSZ examples
Previous: hkf
water> ## Don't show:
water> data(thermo)
thermo: loaded 1997 aqueous, 3089 total species to thermo$obigt
thermo: loaded 5264 proteins to thermo$ECO
thermo: loaded 6717 proteins to thermo$SGD
thermo: loaded 4155 localizations and 3570 abundances to thermo$yeastgfp
water> ## End Don't show
water> ## set temperature, density
water> T <- 500; rho <- 838.0235
water> # calculate pressure
water> P <- as.numeric(water.IAPWS95('P',T=T,rho=rho))
water> # output table of test values
water> water.IAPWS95('test')
[1] "expt 1: T=500,rho=838.025"
[1] "expt 2: T=647,rho=358"
idealgas.1 residual.1 idealgas.2 residual.2
phi 2.0479773 -3.4269321 -1.5631960 -1.2120266
phi.delta 0.3842367 -0.3643667 0.8994413 -0.7140120
phi.delta.delta -0.1476379 0.8560637 -0.8089947 0.4757806
phi.tau 9.0461111 -5.8140344 9.8034392 -3.2172250
phi.tau.tau -1.9324919 -2.2344074 -3.4331633 -9.9602951
phi.delta.tau 0.0000000 -1.1217691 0.0000000 -1.3321472
T rho p s u h
triple.liquid 273.160 9.99793e+02 0.0015548952 6.331605e-08 1.727836e-05 1.572495e-03
triple.vapor 273.160 4.85458e-03 0.0006116553 9.155493e+00 2.374920e+03 2.500915e+03
boiling.liquid 373.124 9.58367e+02 0.0998144792 1.306919e+00 4.189512e+02 4.190554e+02
boiling.vapor 373.124 5.97657e-01 0.1013249528 7.354426e+00 2.505992e+03 2.675529e+03
critical 647.096 3.22000e+02 NaN NaN NaN NaN
g cv cp mu
triple.liquid 1.555200e-03 4.217427 4.219907 -0.0002415825
triple.vapor 7.226521e-04 1.418384 1.884352 0.7045275101
boiling.liquid -6.858736e+01 3.768298 4.215647 -0.0002040517
boiling.vapor -6.858408e+01 1.555764 2.079938 0.0981709913
critical NaN NaN NaN NaN
NULL
water> # calculate dielectric constant
water> water.AW90(T=T,rho=rho,P=P)
[1] 30.81207
water> # find water density for this T, P
water> water('rho',T=T,P=convert(P,'bar'))
rho
1 838.1686
water> ## density along saturation curve
water> T <- seq(273.15,623.15,25)
water> water.WP02(T=T) # liquid from WP02
[1] 999.7885 996.9994 987.9908 974.8091 958.3468 939.0253 917.0095 892.2849 864.6701 833.7879
[11] 798.9850 759.1563 712.3294 654.4707 574.4712
water> water.WP02('rho.vapor',T) # steam from WP02
[1] 4.850914e-03 2.307363e-02 8.315054e-02 2.421845e-01 5.980992e-01 1.298500e+00 2.547898e+00
[8] 4.617113e+00 7.860972e+00 1.275428e+01 1.996592e+01 3.051875e+01 4.617059e+01 7.050918e+01
[15] 1.135604e+02
water> water('rho',T=T,P='Psat') # liquid from SUPCRT92
rho
1 999.8289
2 997.0614
3 988.0295
4 974.8643
5 958.3926
6 939.0726
7 917.0577
8 892.3427
9 864.7434
10 833.8733
11 799.0719
12 759.2362
13 712.4075
14 654.5772
15 574.6875
water> # values of the density, Psat, Gibbs energy
water> water(c('rho','psat','G'),T=T,P='Psat')
rho Psat G
1 999.8289 1.000000 -56289.50
2 997.0614 1.000000 -56687.71
3 988.0295 1.000000 -57123.89
4 974.8643 1.000000 -57594.93
5 958.3926 1.013220 -58098.40
6 939.0726 2.320144 -58631.71
7 917.0577 4.757169 -59193.26
8 892.3427 8.918049 -59781.38
9 864.7434 15.536499 -60394.50
10 833.8733 25.478603 -61031.25
11 799.0719 39.736493 -61690.35
12 759.2362 59.431251 -62370.65
13 712.4075 85.837843 -63071.13
14 654.5772 120.457572 -63790.84
15 574.6875 165.211289 -64528.89
water> # derivatives of the dielectric constant (Born functions)
water> water(c('Q','Y','X','U'),T=T)
Q Y X U
1 6.356099e-07 -4.983317e-05 -3.702958e-07 -67885.49
2 6.634152e-07 -5.798650e-05 -3.055586e-07 -67434.49
3 7.601678e-07 -6.564514e-05 -3.122511e-07 -66984.28
4 9.141878e-07 -7.376119e-05 -3.393917e-07 -66533.88
5 1.136993e-06 -8.272774e-05 -3.804005e-07 -66081.54
6 1.453281e-06 -9.290949e-05 -4.390711e-07 -65625.71
7 1.906248e-06 -1.048443e-04 -5.249732e-07 -65165.09
8 2.568927e-06 -1.193620e-04 -6.547650e-07 -64698.24
9 3.570729e-06 -1.378067e-04 -8.588469e-07 -64223.03
10 5.156325e-06 -1.624786e-04 -1.198071e-06 -63736.36
11 7.831114e-06 -1.976426e-04 -1.809707e-06 -63233.72
12 1.277509e-05 -2.521721e-04 -3.057376e-06 -62708.54
13 2.328664e-05 -3.479667e-04 -6.160726e-06 -62150.13
14 5.175063e-05 -5.565210e-04 -1.719160e-05 -61537.93
15 1.836321e-04 -1.292075e-03 -1.075086e-04 -60816.67
water> # now at constant pressure
water> water(c('Q','Y','X','U'),T=T,P=2000)
Q Y X U
1 3.008917e-07 -4.758098e-05 -2.234814e-07 -67896.54
2 3.635263e-07 -5.305377e-05 -2.146022e-07 -67491.35
3 4.298437e-07 -5.840948e-05 -2.156297e-07 -67080.52
4 5.099209e-07 -6.387878e-05 -2.221501e-07 -66669.54
5 6.091712e-07 -6.951973e-05 -2.291804e-07 -66258.66
6 7.327542e-07 -7.534558e-05 -2.371544e-07 -65847.43
7 8.867312e-07 -8.139605e-05 -2.473538e-07 -65436.01
8 1.078433e-06 -8.774143e-05 -2.609294e-07 -65024.81
9 1.316934e-06 -9.447920e-05 -2.789023e-07 -64614.02
10 1.613599e-06 -1.017308e-04 -3.021737e-07 -64203.50
11 1.982701e-06 -1.096385e-04 -3.315172e-07 -63792.86
12 2.442138e-06 -1.183624e-04 -3.675474e-07 -63381.60
13 3.014248e-06 -1.280750e-04 -4.106622e-07 -62969.21
14 3.726675e-06 -1.389555e-04 -4.609576e-07 -62555.26
15 4.613237e-06 -1.511802e-04 -5.181166e-07 -62139.48
water> ## NaCl dissocation logK f(T,P)
water> # after Shock et al., 1992, Fig. 1
water> # make note of the warning in the subcrt help page
water> species <- c('NaCl','Na+','Cl-')
water> coeffs <- c(-1,1,1)
water> # start a new plot with the experimental data
water> thermo.plot.new(xlim=c(0,1000),ylim=c(-5.5,1),
water+ xlab=axis.label("T"),ylab=axis.label("logK"))
water> expt <- thermo$expt$SOJSH92
water> points(expt$T,expt$logK,pch=expt$pch)
water> T <- list(seq(0,370,25),seq(265,465,25),
water+ seq(285,760,25),seq(395,920,25))
water> for(i in 5:9) T[[i]] <- seq(400,1000,25)
water> P <- list("Psat",500,1000,1500,2000,2500,3000,3500,4000)
water> for(i in 1:length(T)) {
water+ s <- subcrt(species,coeffs,T=T[[i]],P=P[[i]])
water+ lines(s$out$T,s$out$logK)
water+ }
subcrt: 3 species at 15 values of T and P (wet)
subcrt: 3 species at 9 values of T and P (wet)
subcrt: 3 species at 20 values of T and P (wet)
subcrt: 3 species at 22 values of T and P (wet)
subcrt: 3 species at 25 values of T and P (wet)
subcrt: 3 species at 25 values of T and P (wet)
subcrt: 3 species at 25 values of T and P (wet)
subcrt: 3 species at 25 values of T and P (wet)
subcrt: 3 species at 25 values of T and P (wet)
water> legend("bottomleft",pch=unique(expt$pch),
water+ legend=unique(expt$source))
water> title(main=paste('NaCl(aq) = Na+ + Cl-\n',
water+ 'Psat and 500-4000 bar, after Shock et al., 1992'))
water> ## comparing the computational options
water> prop <- c('A','G','S','U','H','Cv','Cp','w','epsilon',
water+ 'Y','Q','X','rho','Psat')
water> thermo$opt$water <- 'SUPCRT'
water> print(water(prop,T=convert(c(25,100,200,300),'K')))
A G S U H Cv Cp w epsilon
1 -55814.06 -56687.71 16.71228 -67434.49 -68316.76 17.81996 18.01160 149727.79 78.24514
2 -57224.77 -58098.40 20.75956 -66081.54 -66963.78 16.23669 18.15793 154171.14 55.49238
3 -59528.16 -60394.50 25.16818 -64223.03 -65097.99 14.28333 19.32884 132927.66 34.90937
4 -62248.93 -63071.13 29.14072 -62150.13 -62980.94 13.18363 24.73943 90574.94 20.40832
Y Q X rho Psat
1 -5.798650e-05 6.634152e-07 -3.055586e-07 997.0614 1.00000
2 -8.272774e-05 1.136993e-06 -3.804005e-07 958.3926 1.01322
3 -1.378067e-04 3.570729e-06 -8.588469e-07 864.7434 15.53650
4 -3.479667e-04 2.328664e-05 -6.160726e-06 712.4075 85.83784
water> thermo$opt$water <- 'IAPWS'
water> print(water(c(prop,'N','UBorn'),T=convert(c(25,100,200,300),'K')))
A G S U H Cv Cp w epsilon
1 -55814.06 -56687.71 16.71230 -67434.50 -68316.76 17.81535 18.00372 4797.468 78.38084
2 -57224.74 -58098.37 20.75976 -66081.45 -66963.69 16.22477 18.15165 6318.160 55.51178
3 -59527.94 -60394.28 25.16593 -64223.90 -65098.85 14.28620 19.35798 9042.855 34.76655
4 -62248.96 -63071.12 29.14717 -62146.48 -62977.24 13.14546 24.75981 10434.110 20.07577
Y Q X rho Psat N UBorn
1 -5.848294e-05 5.934840e-07 -2.773788e-07 997.0470 1.00000 -1.290808e-10 4.071396e-09
2 -8.298249e-05 1.159833e-06 -3.876482e-07 958.3491 1.01418 -4.967554e-10 1.196612e-08
3 -1.408906e-04 3.830012e-06 -9.164440e-07 864.6581 15.54939 -3.593861e-09 5.491354e-08
4 -3.623493e-04 2.480856e-05 -6.426547e-06 712.1356 85.87867 -9.333187e-08 7.777806e-07
water> # fixme: things seem to be working except speed of
water> # sound in our IAPWS calculations
water>
water> # calculating Q Born function
water> # after Table 22 of Johnson and Norton, 1991
water> thermo$opt$water <- 'SUPCRT'
water> T <- rep(c(375,400,425,450,475),each=5)
water> P <- rep(c(250,300,350,400,450),5)
water> t <- water('Q',T=convert(T,'K'),P=P)
water> # the rest is to make a readable table
water> t <- as.data.frame(matrix(t[[1]],nrow=5))
water> colnames(t) <- T[1:5*5]
water> rownames(t) <- P[1:5]
water> print(t)
375 400 425 450 475
250 0.0004550761 0.0038535648 0.0024521064 0.0019850838 0.001722244
300 0.0001576272 0.0024803336 0.0026610046 0.0019996332 0.001690020
350 0.0000942922 0.0004264376 0.0020475555 0.0019040904 0.001611918
400 0.0000668617 0.0001945159 0.0007988531 0.0014967134 0.001446774
450 0.0000515948 0.0001206553 0.0003588294 0.0008965782 0.001165101
Next: subcrt
Up: CHNOSZ examples
Previous: hkf
