diff options
| author |  Linus Torvalds <torvalds@linux-foundation.org> | 2016-03-09 16:41:45 -0800 |
|---|---|---|
| committer |  Dirk Hohndel <dirk@hohndel.org> | 2016-03-13 12:43:30 -0700 |
| commit | 6b5e22a68e8c86be9eb7c297d353191b655dbc25 (patch) | |
| tree | 5c38b6b015be87ccef4d63fa34ca84fbd9be5add /subsurface-core | |
| parent | 4f935595d94c7207fa1892cc44a92e48f081b13a (diff) | |
| download | subsurface-6b5e22a68e8c86be9eb7c297d353191b655dbc25.tar.gz | |
gas model: add proper He compressibility data and do a least-squares fit
Lubomir pointed to exactly where he got his data from, so I added that
raw Helium data to the R script, and let the least-squares fit just take
care of the interpolation between 273K and 323K.
Signed-off-by: Linus Torvalds <torvalds@linux-foundation.org>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>
Diffstat (limited to 'subsurface-core')
| -rw-r--r-- | subsurface-core/compressibility.r | 39 |
1 files changed, 38 insertions, 1 deletions
diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r index 15aca6ea2..66310f3aa 100644 --- a/subsurface-core/compressibility.r +++ b/subsurface-core/compressibility.r @@ -39,7 +39,6 @@ o2 = c(0.9994, 0.9968, 0.9941, 0.9884, 0.9771, 0.9676, 0.9597, 0.9542, 0.9560, 0 n2 = c(0.9998, 0.9990, 0.9983, 0.9971, 0.9964, 0.9973, 1.0000, 1.0052, 1.0559, 1.1422, 1.2480, 1.3629) he = c(1.0005, 1.0024, 1.0048, 1.0096, 1.0191, 1.0286, 1.0381, 1.0476, 1.0943, 1.1402, 1.1854, 1.2297) - options(digits=15) # @@ -76,3 +75,41 @@ summary(hefit) new = data.frame(x = seq(min(x),max(x),len=200)) lines(new$x,predict(hefit,newdata=new)) + +# +# Raw data from VOLUMETRIC BEHAVIOR OF HELIUM-ARGON MIXTURES [..] +# T=323.15K (50 C) +p323atm = c(674.837, 393.223, 237.310, 146.294, 91.4027, 57.5799, 36.4620, 23.1654, 14.7478, 9.4017, 5.9987, 3.8300, + 540.204, 319.943, 195.008, 120.951, 75.8599, 47.9005, 30.3791, 19.3193, 12.3080, 7.8495, 5.0100, 3.1992) + +Hez323 = c(1.28067, 1.16782, 1.10289, 1.06407, 1.04028, 1.02548, 1.01617, 1.01029, 1.00656, 1.00418, 1.00267, 1.00171, + 1.22738, 1.13754, 1.08493, 1.05312, 1.03349, 1.02122, 1.01349, 1.00859, 1.00548, 1.00349, 1.00223, 1.00143) + + +# T=273.15 (0 C) +p273atm = c(683.599, 391.213, 233.607, 143.091, 89.0521, 55.9640, 35.3851, 22.4593, 14.2908, 9.1072, 5.8095, 3.7083, + 534.047, 312.144, 188.741, 116.508, 72.8529, 45.9194, 29.0883, 18.4851, 11.7702, 7.5040, 4.7881, 3.0570) + +Hez273 = c(1.33969, 1.19985, 1.12121, 1.07494, 1.04689, 1.02957, 1.01874, 1.01191, 1.00758, 1.00484, 1.00309, 1.00197, + 1.26914, 1.16070, 1.09837, 1.06118, 1.03843, 1.02429, 1.01541, 1.00980, 1.00625, 1.00398, 1.00254, 1.00162) + +p323 = p323atm * 1.01325 +p273 = p273atm * 1.01325 + +x2=append(p323,p273) +he2=append(Hez323,Hez273) + +plot(x2,he2) + +hefit2 = nls(he2 ~ 1.0 + p1*x2 + p2*x2^2 + p3*x2^3, + start=list(p1=0,p2=0,p3=0)) +summary(hefit2) + +he3 = function(bar) +{ + 1.0 +0.00047961098687979363 * bar -0.00000004077670019935 * bar^2 +0.00000000000077707035 * bar^3 +} + +new = data.frame(x2 = seq(min(x2),max(x2),len=200)) +lines(new$x2,predict(hefit2,newdata=new)) +curve(he3, min(x2),max(x2),add=TRUE) |