From: Linus Torvalds <[email protected]> Date: Wed, 9 Mar 2016 16:41:45 -0800 Subject: [PATCH 1/2] 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 <[email protected]> --- So doing the least-square fitting to *also* look for a middle ground between 273K and 323K may be questionable, but it works, and has the advantage of directly using the experimental data that we have rather than doing any other intermediate munging of the experimental data. The He compressibility curve itself doesn't change noticeably. Yes, you can see it is two different curves, but just barely. subsurface-core/compressibility.r | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r index 15aca6ea27c9..66310f3aa8d4 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) -- 2.8.0.rc2 _______________________________________________ subsurface mailing list [email protected] http://lists.subsurface-divelog.org/cgi-bin/mailman/listinfo/subsurface
