From c776853a67d244af1669e8bbb7e96883d0851792 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 00:06:33 -0400 Subject: [PATCH 01/11] delete manifest.toml --- Manifest.toml | 351 -------------------------------------------------- 1 file changed, 351 deletions(-) delete mode 100644 Manifest.toml diff --git a/Manifest.toml b/Manifest.toml deleted file mode 100644 index 59f8db6..0000000 --- a/Manifest.toml +++ /dev/null @@ -1,351 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.9.3" -manifest_format = "2.0" -project_hash = "ccd43927e850a716519e856a10fcafa4145bbae1" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.1" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[deps.ChainRulesCore]] -deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "e0af648f0692ec1691b5d094b8724ba1346281cf" -uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.18.0" - - [deps.ChainRulesCore.extensions] - ChainRulesCoreSparseArraysExt = "SparseArrays" - - [deps.ChainRulesCore.weakdeps] - SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - -[[deps.CommonSolve]] -git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" -uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" -version = "0.2.4" - -[[deps.CommonSubexpressions]] -deps = ["MacroTools", "Test"] -git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" -uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" -version = "0.3.0" - -[[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "8a62af3e248a8c4bad6b32cbbe663ae02275e32c" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.10.0" -weakdeps = ["Dates", "LinearAlgebra"] - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" - -[[deps.ConstructionBase]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "c53fc348ca4d40d7b371e71fd52251839080cbc9" -uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.4" - - [deps.ConstructionBase.extensions] - ConstructionBaseIntervalSetsExt = "IntervalSets" - ConstructionBaseStaticArraysExt = "StaticArrays" - - [deps.ConstructionBase.weakdeps] - IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" - StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[deps.DiffResults]] -deps = ["StaticArraysCore"] -git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" -uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" -version = "1.1.0" - -[[deps.DiffRules]] -deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272" -uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.15.1" - -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" - -[[deps.ForwardDiff]] -deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] -git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad" -uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.36" - - [deps.ForwardDiff.extensions] - ForwardDiffStaticArraysExt = "StaticArrays" - - [deps.ForwardDiff.weakdeps] - StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" - -[[deps.Future]] -deps = ["Random"] -uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - -[[deps.JLLWrappers]] -deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.5.0" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.3" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "7.84.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.10.2+0" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.26" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[deps.MacroTools]] -deps = ["Markdown", "Random"] -git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48" -uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.11" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+0" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.10.11" - -[[deps.NaNMath]] -deps = ["OpenLibm_jll"] -git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" -uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.0.2" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.21+4" - -[[deps.OpenLibm_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+0" - -[[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" - -[[deps.PrecompileTools]] -deps = ["Preferences"] -git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f" -uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.2.0" - -[[deps.Preferences]] -deps = ["TOML"] -git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.1" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Random]] -deps = ["SHA", "Serialization"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[deps.Roots]] -deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"] -git-tree-sha1 = "06b5ac80ff1b88bd82df92c1c1875eea3954cd6e" -uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "2.0.20" - - [deps.Roots.extensions] - RootsForwardDiffExt = "ForwardDiff" - RootsIntervalRootFindingExt = "IntervalRootFinding" - RootsSymPyExt = "SymPy" - RootsSymPyPythonCallExt = "SymPyPythonCall" - - [deps.Roots.weakdeps] - ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" - IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" - SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" - SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.Setfield]] -deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] -git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" -uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" -version = "1.1.1" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[deps.SpecialFunctions]] -deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.1" -weakdeps = ["ChainRulesCore"] - - [deps.SpecialFunctions.extensions] - SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" - -[[deps.StaticArraysCore]] -git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" -uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.2" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.Unitful]] -deps = ["Dates", "LinearAlgebra", "Random"] -git-tree-sha1 = "a72d22c7e13fe2de562feda8645aa134712a87ee" -uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.17.0" - - [deps.Unitful.extensions] - ConstructionBaseUnitfulExt = "ConstructionBase" - InverseFunctionsUnitfulExt = "InverseFunctions" - - [deps.Unitful.weakdeps] - ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+0" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+0" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.48.0+0" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+0" From 3c16b0075790915ff628c7432bf3b94dae4be6f0 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 00:07:31 -0400 Subject: [PATCH 02/11] add Manifest.toml to .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index b0383ea..d5cc55a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.jl.*.cov *.jl.mem *.pdf +Manifest.toml \ No newline at end of file From 70347ec0cb91bb40e9f502463bc9551a4ee10919 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 00:08:29 -0400 Subject: [PATCH 03/11] add Th, Ts solvers for regions 3, 5 --- Project.toml | 6 +- README.md | 4 +- src/SteamTables.jl | 1277 ++++++++++++++++++++------------------------ 3 files changed, 578 insertions(+), 709 deletions(-) diff --git a/Project.toml b/Project.toml index f6e938f..07b2773 100644 --- a/Project.toml +++ b/Project.toml @@ -4,17 +4,13 @@ author = ["Braam van Dyk "] version = "1.4.2" [deps] -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] -ForwardDiff = "0.8, 0.10" -Roots = "1, 2" +PrecompileTools = "1" Unitful = "1" julia = "1" -PrecompileTools = "1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index 674ca26..2b3b6bc 100644 --- a/README.md +++ b/README.md @@ -34,11 +34,11 @@ DeltaHvap |J/kg |Latent heat of vaporisation #### P and h - SpecificG_Ph, SpecificF_Ph, SpecificV_Ph, SpecificU_Ph, SpecificS_Ph, SpecificH_Ph, SpecificCP_Ph, SpecificCV_Ph, SpeedOfSound_Ph, Quality_Ph + SpecificG_Ph, SpecificF_Ph, SpecificV_Ph, SpecificU_Ph, SpecificS_Ph, SpecificH_Ph, SpecificCP_Ph, SpecificCV_Ph, SpeedOfSound_Ph, Quality_Ph, Temperature_Ph #### P and s - SpecificG_Ps, SpecificF_Ps, SpecificV_Ps, SpecificU_Ps, SpecificS_Ps, SpecificH_Ps, SpecificCP_Ps, SpecificCV_Ps, SpeedOfSound_Ps, Quality_Ps + SpecificG_Ps, SpecificF_Ps, SpecificV_Ps, SpecificU_Ps, SpecificS_Ps, SpecificH_Ps, SpecificCP_Ps, SpecificCV_Ps, SpeedOfSound_Ps, Quality_Ps, Temperature_Ps #### T and h Quality_Th diff --git a/src/SteamTables.jl b/src/SteamTables.jl index 60bcc9f..8fcd252 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -23,12 +23,12 @@ P and T P and h SpecificG_Ph, SpecificF_Ph, SpecificV_Ph, SpecificU_Ph, SpecificS_Ph, SpecificH_Ph, SpecificCP_Ph, SpecificCV_Ph, SpeedOfSound_Ph - Quality_Ph + Quality_Ph, Temperature_Ph P and s SpecificG_Ps, SpecificF_Ps, SpecificV_Ps, SpecificU_Ps, SpecificS_Ps, SpecificH_Ps, SpecificCP_Ps, SpecificCV_Ps, SpeedOfSound_Ps - Quality_Ps + Quality_Ps, Temperature_Ps T and h Quality_Th @@ -87,14 +87,12 @@ export Psat, Tsat, SpecificG, SpecificF, SpecificV, SpecificU, SpecificS, SpecificH, SpecificCP, SpecificCV, SpeedOfSound, SpecificG_Ph, SpecificF_Ph, SpecificV_Ph, SpecificU_Ph, SpecificS_Ph, - SpecificCP_Ph, SpecificCV_Ph, SpeedOfSound_Ph, + SpecificCP_Ph, SpecificCV_Ph, SpeedOfSound_Ph, Temperature_Ph, SpecificG_Ps, SpecificF_Ps, SpecificV_Ps, SpecificU_Ps, SpecificH_Ps, - SpecificCP_Ps, SpecificCV_Ps, SpeedOfSound_Ps, + SpecificCP_Ps, SpecificCV_Ps, SpeedOfSound_Ps, Temperature_Ps, R, Tc, Pc, ρc, T3, P3, Mr -using Roots using Unitful -using ForwardDiff struct UnitsError <: Exception @@ -111,11 +109,11 @@ Returns the boundary between regions 2 and 3. InputType is either :T or :P to indicate that InputValue is temperature [K] or pressure [MPa]. The complimentary value is returned. """ function B23(InputType::Symbol, InputValue) - n = [0.348_051_856_289_69E3, + n = (0.348_051_856_289_69E3, -0.116_718_598_799_75E1, 0.101_929_700_393_26E-2, 0.572_544_598_627_46E3, - 0.139_188_397_788_70E2] + 0.139_188_397_788_70E2) if InputType == :T return n[1] + n[2] * InputValue + n[3] * InputValue^2 @@ -136,11 +134,11 @@ end Valid from saturation line at 554.485K & 6.54670MPa to 1019.32K & 100MPa """ function B2bc(InputType::Symbol, InputValue) - n = [0.905_842_785_147_23E3, + n = (0.905_842_785_147_23E3, -0.679_557_863_992_41, 0.128_090_027_301_36E-3, 0.265_265_719_084_28E4, - 0.452_575_789_059_48E1] + 0.452_575_789_059_48E1) if InputType == :h return n[1] + n[2] * InputValue + n[3] * InputValue^2 @@ -1150,6 +1148,102 @@ function Region2c_TPh(P, h) return sum([n[i]*ππ^I[i]*ηη^J[i] for i=1:23]) end +""" + Region3_TPh + + Returns T [K] from P[MPa] and h[kJ/kg] in Region 3. + Pressures in MPa and temperature in [K] +""" +function Region3_TPh(P, h) + Tlow = 623.15 + Thigh = B23(:P, P) + if 16.529164252604478 <= P <= 22.064 + #we are in the saturation boundary + Tsat = Tsat(P) + hl,hv = SatHL(Tsat),SatHV(Tsat) + if hl <= h <= hv + return Tsat + elseif h > hv #gas phase, interpolate with h(T_high) + hlow = hv + hhigh = Region3(:SpecificH, P, Thigh) + elseif h < hl #liquid phase, interpolate with h(T_low) + hlow = Region3(:SpecificH, P, Tlow) + hhigh = hl + end + else #22.064 < p <= 100 + hlow = Region3(:SpecificH, P, Tlow) + hhigh = Region3(:SpecificH, P, Thigh) + end + T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + ρ = Region3_ρPT(P, T) + Told = T + for i in 1:100 + h_i = Region3_ρ(:SpecificH, ρ, T) + if hlow < h_i + hlow = h_i + Tlow = T + else + hhigh = h_i + Thigh = T + end + Told = T + T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + if abs(T - Told) < 1e-12 + return T + end + ρ = Region3_ρPT(P, T) + end + throw(error("Region3_TPh: temperature iterations failed to converge.")) +end + +""" + Region3_TPs + + Returns T [K] from P[MPa] and s[kJ/kg/K] in Region 3. + Pressures in MPa and temperature in [K] +""" +function Region3_TPs(P, s) + Tlow = 623.15 + Thigh = B23(:P, P) + if 16.529164252604478 <= P <= 22.064 + #we are in the saturation boundary + Tsat = Tsat(P) + sl,sv = SatSL(Tsat),SatSV(Tsat) + if sl <= s <= sv + return Tsat + elseif s > sv #gas phase, interpolate with s(T_high) + slow = sv + shigh = Region3(:SpecificS, P, Thigh) + elseif s < sl #liquid phase, interpolate with s(T_low) + slow = Region3(:SpecificS, P, Tlow) + shigh = sl + end + else #22.064 < p <= 100 + slow = Region3(:SpecificS, P, Tlow) + shigh = Region3(:SpecificS, P, Thigh) + end + T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + ρ = Region3_ρPT(P, T) + Told = T + for i in 1:100 + s_i = Region3_ρ(:SpecificS, ρ, T) + if slow < s_i + slow = s_i + Tlow = T + else + shigh = s_i + Thigh = T + end + Told = T + T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + if abs(T - Told) < 1e-12 + return T + end + ρ = Region3_ρPT(P, T, ρ) + end + throw(error("Region3_TPs: temperature iterations failed to converge.")) +end + """ Region2a_TPs @@ -1589,7 +1683,7 @@ function Region3_ρ(Output::Symbol, ρ, T) ρstar = ρc #kg/m3 Tstar = Tc #K - n = [0.106_580_700_285_13E1, + n = (0.106_580_700_285_13E1, -0.157_328_452_902_39E2, 0.209_443_969_743_07E2, -0.768_677_078_787_16E1, @@ -1628,9 +1722,9 @@ function Region3_ρ(Output::Symbol, ρ, T) 0.323_089_047_037_11E-2, 0.809_648_029_962_15E-4, -0.165_576_797_950_37E-3, - -0.449_238_990_618_15E-4] + -0.449_238_990_618_15E-4) - I = [0, + I = (0, 0, 0, 0, @@ -1669,9 +1763,9 @@ function Region3_ρ(Output::Symbol, ρ, T) 9, 10, 10, - 11] + 11) - J = [0, + J = (0, 0, 1, 2, @@ -1710,7 +1804,7 @@ function Region3_ρ(Output::Symbol, ρ, T) 26, 0, 1, - 26] + 26) δ = ρ / ρstar τ = Tstar / T @@ -1745,23 +1839,158 @@ function Region3_ρ(Output::Symbol, ρ, T) end end +function Region3_p∂p∂ρ(ρ, T) + ρstar = ρc #kg/m3 + Tstar = Tc #K + + n = (0.106_580_700_285_13E1, + -0.157_328_452_902_39E2, + 0.209_443_969_743_07E2, + -0.768_677_078_787_16E1, + 0.261_859_477_879_54E1, + -0.280_807_811_486_20E1, + 0.120_533_696_965_17E1, + -0.845_668_128_125_02E-2, + -0.126_543_154_777_14E1, + -0.115_244_078_066_81E1, + 0.885_210_439_843_18, + -0.642_077_651_816_07, + 0.384_934_601_866_71, + -0.852_147_088_242_06, + 0.489_722_815_418_77E1, + -0.305_026_172_569_65E1, + 0.394_205_368_791_54E-1, + 0.125_584_084_243_08, + -0.279_993_296_987_10, + 0.138_997_995_694_60E1, + -0.201_899_150_235_70E1, + -0.821_476_371_739_63E-2, + -0.475_960_357_349_23, + 0.439_840_744_735_00E-1, + -0.444_764_354_287_39, + 0.905_720_707_197_33, + 0.705_224_500_879_67, + 0.107_705_126_263_32, + -0.329_136_232_589_54, + -0.508_710_620_411_58, + -0.221_754_008_730_96E-1, + 0.942_607_516_650_92E-1, + 0.164_362_784_479_61, + -0.135_033_722_413_48E-1, + -0.148_343_453_524_72E-1, + 0.579_229_536_280_84E-3, + 0.323_089_047_037_11E-2, + 0.809_648_029_962_15E-4, + -0.165_576_797_950_37E-3, + -0.449_238_990_618_15E-4) + + I = (0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1, + 1, + 1, + 1, + 2, + 2, + 2, + 2, + 2, + 2, + 3, + 3, + 3, + 3, + 3, + 4, + 4, + 4, + 4, + 5, + 5, + 5, + 6, + 6, + 6, + 7, + 8, + 9, + 9, + 10, + 10, + 11) + + J = (0, + 0, + 1, + 2, + 7, + 10, + 12, + 23, + 2, + 6, + 15, + 17, + 0, + 2, + 6, + 7, + 22, + 26, + 0, + 2, + 4, + 16, + 26, + 0, + 2, + 4, + 26, + 1, + 3, + 26, + 0, + 2, + 26, + 2, + 26, + 2, + 26, + 0, + 1, + 26) + + δ = ρ / ρstar + τ = Tstar / T + logδ = log(δ) + logτ = log(τ) + ϕ = n[1]*logδ+ sum([n[i]*exp(logδ*I[i]+logτ*J[i]) for i=2:40]) + ϕ_δ = n[1]/δ + sum([n[i]*I[i]*exp(logδ*(I[i]-1)+logτ*J[i]) for i=2:40]) + ϕ_δδ = -n[1]/(δ^2) + sum([n[i]*I[i]*(I[i]-1)*exp(logδ*(I[i]-2)+logτ*J[i]) for i=2:40]) + #ϕ_δ = n[1]/δ + sum([n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40]) + #ϕ_δδ = -n[1]/(δ^2) + sum([n[i]*I[i]*(I[i]-1)*(δ^(I[i]-2))*(τ^J[i]) for i=2:40]) + + k = R*T*ρstar/1000 + #p = ρ*R*T*δ*ϕ_δ/1000 + p = k*δ*δ*ϕ_δ + ∂p∂δ = k*(2*δ*ϕ_δ + δ*δ*ϕ_δδ) + ∂p∂ρ = ∂p∂δ/ρstar + return p, ∂p∂ρ +end """ Initialise from ideal gas law, then use root finder to calculate P by iterating on Region3ρ. - Pass through properties from Region3ρ + Pass through properties from Region3_ρ """ - -#TODO Remove Roots dependency here - function Region3(Output::Symbol, P, T) - # Start with a higher density to fix the problems in convergence. - # TODO Replace this with a density from Peng-Robinson or similar EoS # ρ0 = 1000*P/(R*T) #Starting value from ideal gas - ρ0 = 500.0 - - f(ρ) = Region3_ρ(:Pressure, ρ, T) - P - ρ = Roots.find_zero(f, ρ0) - + ρ = Region3_ρPT(P,T) if Output == :SpecificV return 1.0 / ρ else @@ -1769,6 +1998,44 @@ function Region3(Output::Symbol, P, T) end end +function Region3_ρ0(P, T) + if T <= 647.96 #we are inside saturation. use the saturation volume as initial guess + if 16.529164252604478 <= P <= Psat(T) #gas phase + #gas phase + ρ0 = 1000*P/(R*T) + else + #liquid phase + ρ0 = SatDensL(T) + end + else #over critical point,it does not matter we start, but a high density point is recomended + ρ0_liquid = 574.6703980963142 #volume at the intersection of regions 1-3-4 + _,∂p∂ρᵢ = Region3_p∂p∂ρ(ρ0_liquid, T) + if ∂p∂ρᵢ > 0 + return ρ0_liquid + else + return 1/Region2(:SpecificV,P,B23(:P,P)) + end + end +end + +Region3_ρPT(P,T) = Region3_ρPT(P,T,Region3_ρ0(P,T)) + +function Region3_ρPT(P,T,ρ0) + #this solver is used to solve the GERG2008 equation of state. + #is a modified newton raphson solver. + logρ = log(ρ0) + for i in 1:100 + ρᵢ = exp(logρ) + Pᵢ,∂p∂ρᵢ = Region3_p∂p∂ρ(ρᵢ, T) + ΔP = (P - Pᵢ) + Δ = ΔP/(ρᵢ*∂p∂ρᵢ) + logρ = logρ + Δ + abs(ΔP) < 3eps(Pᵢ) && return exp(logρ) + abs(Δ) < 1e-12 && return exp(logρ) + end + return zero(ρ0)/zero(ρ0) +end + """ Region4 @@ -1780,7 +2047,7 @@ end or pressure [MPa]. The complimentary value is returned. """ function Region4(InputType::Symbol, InputValue) - n = [0.116_705_214_527_67E4, + n = (0.116_705_214_527_67E4, -0.724_213_167_032_06E6, -0.170_738_469_400_92E2, 0.120_208_247_024_70E5, @@ -1789,7 +2056,7 @@ function Region4(InputType::Symbol, InputValue) -0.482_326_573_615_91E4, 0.405_113_405_420_57E6, -0.238_555_575_678_49, - 0.650_175_348_447_98E3] + 0.650_175_348_447_98E3) if InputType == :T Θ = InputValue + n[9]/(InputValue - n[10]) @@ -1813,7 +2080,6 @@ function Region4(InputType::Symbol, InputValue) end end - """ Region5 @@ -1911,6 +2177,54 @@ function Region5(Output::Symbol, P, T) end end +function Region5_TPh(P, h) + hhighT = Region5(:SpecificH, P, 1073.15) + hlowT = Region5(:SpecificH, P, 2273.15) + hhigh,hlow = minmax(hhighT,hlowT) + T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + Told = T + for i in 1:100 + h_i = Region5(:SpecificH, P, T) + if hlow < h_i + hlow = h_i + Tlow = T + else + hhigh = h_i + Thigh = T + end + Told = T + T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + if abs(T - Told) < 1e-12 + return T + end + end + throw(error("Region5_TPh: temperature iterations failed to converge.")) +end + +function Region5_TPs(P, s) + shighT = Region5(:SpecificH, P, 1073.15) + slowT = Region5(:SpecificH, P, 2273.15) + shigh,slow = minmax(shighT,slowT) + T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + Told = T + for i in 1:100 + s_i = Region5(:SpecificS, P, T) + if slow < s_i + slow = s_i + Tlow = T + else + shigh = s_i + Thigh = T + end + Told = T + T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + if abs(T - Told) < 1e-12 + return T + end + end + throw(error("Region5_TPs: temperature iterations failed to converge.")) +end + """ RegionID @@ -1996,6 +2310,9 @@ function RegionID_Ph(P, h)::Symbol 2a: P ≤ 4Mpa 2b: s ≥ 5.85 kJ/kgK (or use function B2bc if P,h specificed) 2c: s ≥ 5.85 kJ/kgK (or use function B2bc if P,h specificed) + + Region 3: + Psat(623.15) ≤ P ≤ 100MPa =# # Check overall region first @@ -2030,6 +2347,8 @@ function RegionID_Ph(P, h)::Symbol end elseif B23(:P, P) ≤ T ≤ 1073.15 return :Region2c #, Region2(:SpecificH, P, T) # Return forward h for consistency check + elseif h <= Region3_ρ(:SpecificH, P, B23(:P, P)) + return :Region3 else throw(DomainError((P, T), "Pressure/Temperature outside valid ranges.")) end @@ -2145,13 +2464,13 @@ end """ function Psat2(T) - a = [ -7.85951783000, + a = ( -7.85951783000, 1.84408259000, -11.78664970000, 22.68074110000, -15.96187190000, 1.80122502000 - ] + ) θ = T/Tc τ = 1 - θ @@ -2160,6 +2479,35 @@ function Psat2(T) + a[5]*τ^4 + a[6]*τ^7.5)) end +function ∂Psat∂T(T) + a = ( -7.85951783000, + 1.84408259000, + -11.78664970000, + 22.68074110000, + -15.96187190000, + 1.80122502000 + ) + θ = T/Tc + τ = 1 - θ + ∂τ = - 1/Tc + ∂T = one(T) + τhalf = sqrt(τ) + τ15 = τ*τhalf + τ2 = τ*τ + τ25 = τ2*τhalf + τ3 = τ2*τ + τ35=τ3*τhalf + τ4=τ2*τ2 + τ65 = τ35*τ3 + τ75 = τ35*τ4 + + x = Tc*(a[1]*τ + a[2]*τ15 + a[3]*τ3 + a[4]*τ35 + a[5]*τ4 + a[6]*τ75) + ∂x = Tc*∂τ*(a[1] + 1.5*a[2]*τhalf + 3.0*a[3]*τ2 + 3.5*a[4]*τ25 + 4.0*a[5]*τ3 + 7.5*a[6]*τ65) + xT = x/T + #res = Pc*exp(xT) + ∂res = Pc*exp(xT)*(∂x - xT*∂T)/T +end + #============================================================================ ============================================================================= @@ -2180,7 +2528,7 @@ function Psat(T) if T3 ≤ T ≤ Tc return Region4(:T, T) else - throw(DomainError(T, "Temperature not between tripple and critical points.")) + throw(DomainError(T, "Temperature not between triple and critical points.")) end end @@ -2195,7 +2543,7 @@ function Psat(T::Q) where Q <: Quantity if T3 ≤ T.val ≤ Tc return Region4(:T, T.val)*u"MPa" else - throw(DomainError(T, "Temperature not between tripple and critical points.")) + throw(DomainError(T, "Temperature not between triple and critical points.")) end end @@ -2211,7 +2559,7 @@ function Tsat(P) if 0.611_213E-3 ≤ P ≤ Pc # Note that botton limit not exactly P3 return Region4(:P, P) else - throw(DomainError(P, "Pressure not between tripple and critical points")) + throw(DomainError(P, "Pressure not between triple and critical points")) end end @@ -2226,10 +2574,23 @@ function Tsat(P::Q) where Q <: Quantity if 0.611_213E-3 ≤ P.val ≤ Pc return Region4(:P, P.val)*u"K" else - throw(DomainError(P, "Pressure not between tripple and critical points")) + throw(DomainError(P, "Pressure not between triple and critical points")) end end +function property_PT(property::Symbol,Region::Symbol,P,T) + if Region == :Region1 + return Region1(property, P, T) + elseif Region == :Region2a || Region == :Region2b || Region == :Region2c || Region == :Region2 + return Region2(property, P, T) + elseif Region == :Region3 + return Region3(property, P, T) + elseif Region == :Region5 + return Region5(property, P, T) + else + throw(DomainError(Region, "Invalid Region")) + end +end """ SpecificG @@ -2241,16 +2602,7 @@ end """ function SpecificG(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificG, P, T) - elseif Region == :Region2 - return Region2(:SpecificG, P, T) - elseif Region == :Region3 - return Region3(:SpecificG, P, T) - elseif Region == :Region5 - return Region5(:SpecificG, P, T) - end + return property_PT(:SpecificG,Region,P,T) end @@ -2259,51 +2611,93 @@ function SpecificG(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity P = 1.0*uconvert(u"MPa", P) T = 1.0*uconvert(u"K", T) Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificG, P.val, T.val)*u"kJ/kg" - elseif Region == :Region2 - return Region2(:SpecificG, P.val, T.val)*u"kJ/kg" - elseif Region == :Region3 - return Region3(:SpecificG, P.val, T.val)*u"kJ/kg" - elseif Region == :Region5 - return Region5(:SpecificG, P.val, T.val)*u"kJ/kg" - end + return property_PT(:SpecificG,Region,P.val,T.val)*u"kJ/kg" catch throw(UnitsError((P,T), "Invalid input units.")) end +end + +""" + Temperature_Ph(P, h) + Utility function that returns the Temperature [K] from P [MPa] and h [kJ/kg] + The explicit backwards equations are only available in regions 1 and 2. For Region 3 and 5 there is an implicit solver that may fail. Input outside these + will result in a DomainError exception. + If inputs have associated units, the value is returned with associated + units of K via Uniful.jl. +""" +function Temperature_Ph(P, h) + Region = RegionID_Ph(P, h) + return _Temperature_Ph(Region, P, h) end +function _Temperature_Ph(Region::Symbol, P, h) + if Region == :Region1 + return Region1_TPh(P, h) + elseif Region == :Region2a + return Region2a_TPh(P, h) + elseif Region == :Region2b + return Region2b_TPh(P, h) + elseif Region == :Region2c + return Region2c_TPh(P, h) + elseif Region == :Region3 + return Region3_TPh(P, h) + elseif Region == :Region5 + return Region5_TPh(P, h) + else + throw(DomainError(Region, "Invalid Region")) + end +end """ - SpecificG_Ph + Temperature_Ph(P, h) - Utility function that returns the Gibbs free energy [kJ/kgK] from P [MPa] and h [kJ/kg] - The explicit backwards equations are only available in regions 1 and 2. Input outside these + Utility function that returns the Temperature [K] from P [MPa] and s [kJ/kg/K] + The explicit backwards equations are only available in regions 1 and 2. For Region 3 and 5 there is an implicit solver that may fail. Input outside these will result in a DomainError exception. If inputs have associated units, the value is returned with associated - units of kJ/kg via Uniful.jl. + units of K via Uniful.jl. """ -function SpecificG_Ph(P, h) - Region = RegionID_Ph(P, h) +function Temperature_Ps(P, s) + Region = RegionID_Ps(P, s) + return _Temperature_Ps(Region, P, s) +end +function _Temperature_Ps(Region::Symbol, P, s) if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificG, P, T) + return Region1_TPs(P, s) elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificG, P, T) + return Region2a_TPs(P, s) elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificG, P, T) + return Region2b_TPs(P, s) elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificG, P, T) + return Region2c_TPs(P, s) + elseif Region == :Region3 + return Region3_TPs(P, s) + elseif Region == :Region5 + return Region5_TPs(P, s) + else + throw(DomainError(Region, "Invalid Region")) end end +""" + SpecificG_Ph + + Utility function that returns the Gibbs free energy [kJ/kgK] from P [MPa] and h [kJ/kg] + The explicit backwards equations are only available in regions 1 and 2. Input outside these + will result in a DomainError exception. + If inputs have associated units, the value is returned with associated + units of kJ/kg via Uniful.jl. +""" +function SpecificG_Ph(P, h) + Region = RegionID_Ph(P, h) + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificG,Region,P,T) +end + + function SpecificG_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity try P = 1.0*uconvert(u"MPa", P) @@ -2317,21 +2711,9 @@ function SpecificG_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificG, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificG, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificG, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificG, P.val, T)*u"kJ/kg" - end -end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificG,Region,P.val,T)*u"kJ/kg" +end """ @@ -2454,23 +2836,10 @@ end """ function SpecificF_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificF, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificF, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificF, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificF, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificF, Region, P, T) end - function SpecificF_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity try P = 1.0*uconvert(u"MPa", P) @@ -2484,20 +2853,8 @@ function SpecificF_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificF, Region, P.val, T)*u"kJ/kg" end @@ -2513,20 +2870,8 @@ end """ function SpecificF_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificF, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificF, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificF, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificF, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificF, Region, P, T) end @@ -2543,20 +2888,8 @@ function SpecificF_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificF, P.val, T)*u"kJ/kg" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificF, Region, P.val, T)*u"kJ/kg" end @@ -2569,16 +2902,7 @@ end """ function SpecificV(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificV, P, T) - elseif Region == :Region2 - return Region2(:SpecificV, P, T) - elseif Region == :Region3 - return Region3(:SpecificV, P, T) - elseif Region == :Region5 - return Region5(:SpecificV, P, T) - end + return property_PT(:SpecificV, Region, P, T) end @@ -2595,16 +2919,7 @@ function SpecificV(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificV, P.val, T.val)*u"m^3/kg" - elseif Region == :Region2 - return Region2(:SpecificV, P.val, T.val)*u"m^3/kg" - elseif Region == :Region3 - return Region3(:SpecificV, P.val, T.val)*u"m^3/kg" - elseif Region == :Region5 - return Region5(:SpecificV, P.val, T.val)*u"m^3/kg" - end + return property_PT(:SpecificV, Region, P.val, T.val)*u"m^3/kg" end @@ -2619,20 +2934,8 @@ end """ function SpecificV_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificV, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificV, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificV, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificV, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificV, Region, P, T) end @@ -2649,20 +2952,8 @@ function SpecificV_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificV, Region, P.val, T)*u"m^3/kg" end @@ -2678,20 +2969,8 @@ end """ function SpecificV_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificV, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificV, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificV, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificV, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificV, Region, P, T) end @@ -2708,20 +2987,8 @@ function SpecificV_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificV, P.val, T)*u"m^3/kg" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificV, Region, P.val, T)*u"m^3/kg" end @@ -2735,16 +3002,7 @@ end """ function SpecificU(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificU, P, T) - elseif Region == :Region2 - return Region2(:SpecificU, P, T) - elseif Region == :Region3 - return Region3(:SpecificU, P, T) - elseif Region == :Region5 - return Region5(:SpecificU, P, T) - end + return property_PT(:SpecificU, Region, P, T) end @@ -2761,16 +3019,7 @@ function SpecificU(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificU, P.val, T.val)*u"kJ/kg" - elseif Region == :Region2 - return Region2(:SpecificU, P.val, T.val)*u"kJ/kg" - elseif Region == :Region3 - return Region3(:SpecificU, P.val, T.val)*u"kJ/kg" - elseif Region == :Region5 - return Region5(:SpecificU, P.val, T.val)*u"kJ/kg" - end + return property_PT(:SpecificU, Region, P.val, T.val)*u"kJ/kg" end @@ -2786,20 +3035,8 @@ end """ function SpecificU_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificU, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificU, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificU, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificU, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificU, Region, P, T) end @@ -2816,20 +3053,8 @@ function SpecificU_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificU, Region, P.val, T)*u"kJ/kg" end @@ -2845,20 +3070,8 @@ end """ function SpecificU_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificU, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificU, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificU, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificU, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificU, Region, P, T) end @@ -2875,20 +3088,8 @@ function SpecificU_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificU, P.val, T)*u"kJ/kg" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificU, Region, P.val, T)*u"kJ/kg" end @@ -2902,16 +3103,7 @@ end """ function SpecificS(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificS, P, T) - elseif Region == :Region2 - return Region2(:SpecificS, P, T) - elseif Region == :Region3 - return Region3(:SpecificS, P, T) - elseif Region == :Region5 - return Region5(:SpecificS, P, T) - end + return property_PT(:SpecificS, Region, P, T) end @@ -2928,16 +3120,7 @@ function SpecificS(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificS, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region2 - return Region2(:SpecificS, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region3 - return Region3(:SpecificS, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region5 - return Region5(:SpecificS, P.val, T.val)*u"kJ/kg/K" - end + return property_PT(:SpecificS, Region, P.val, T.val)*u"kJ/kg/K" end @@ -2953,20 +3136,8 @@ end """ function SpecificS_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificS, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificS, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificS, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificS, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificS, Region, P, T) end @@ -2983,20 +3154,8 @@ function SpecificS_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificS, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificS, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificS, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificS, P.val, T)*u"kJ/kg/K" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificS, Region, P.val, T)*u"kJ/kg/K" end @@ -3015,16 +3174,7 @@ end """ function SpecificH(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificH, P, T) - elseif Region == :Region2 - return Region2(:SpecificH, P, T) - elseif Region == :Region3 - return Region3(:SpecificH, P, T) - elseif Region == :Region5 - return Region5(:SpecificH, P, T) - end + return property_PT(:SpecificH, Region, P, T) end @@ -3041,16 +3191,7 @@ function SpecificH(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificH, P.val, T.val)*u"kJ/kg" - elseif Region == :Region2 - return Region2(:SpecificH, P.val, T.val)*u"kJ/kg" - elseif Region == :Region3 - return Region3(:SpecificH, P.val, T.val)*u"kJ/kg" - elseif Region == :Region5 - return Region5(:SpecificH, P.val, T.val)*u"kJ/kg" - end + return property_PT(:SpecificH, Region, P.val, T.val)*u"kJ/kg" end @@ -3071,20 +3212,8 @@ end """ function SpecificH_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificH, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificH, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificH, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificH, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificH, Region, P, T) end @@ -3101,20 +3230,8 @@ function SpecificH_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificH, P.val, T)*u"kJ/kg" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificH, P.val, T)*u"kJ/kg" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificH, P.val, T)*u"kJ/kg" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificH, P.val, T)*u"kJ/kg" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificH, Region, P.val, T)*u"kJ/kg" end @@ -3128,16 +3245,7 @@ end """ function SpecificCP(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificCP, P, T) - elseif Region == :Region2 - return Region2(:SpecificCP, P, T) - elseif Region == :Region3 - return Region3(:SpecificCP, P, T) - elseif Region == :Region5 - return Region5(:SpecificCP, P, T) - end + return property_PT(:SpecificCP, Region, P, T) end @@ -3154,16 +3262,7 @@ function SpecificCP(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificCP, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region2 - return Region2(:SpecificCP, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region3 - return Region3(:SpecificCP, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region5 - return Region5(:SpecificCP, P.val, T.val)*u"kJ/kg/K" - end + return property_PT(:SpecificCP, Region, P.val, T.val)*u"kJ/kg/K" end @@ -3179,20 +3278,8 @@ end """ function SpecificCP_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificCP, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificCP, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificCP, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificCP, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificCP, Region, P, T) end @@ -3209,20 +3296,8 @@ function SpecificCP_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificCP, Region, P.val, T)*u"kJ/kg/K" end @@ -3238,20 +3313,8 @@ end """ function SpecificCP_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificCP, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificCP, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificCP, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificCP, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificCP, Region, P, T) end @@ -3268,20 +3331,8 @@ function SpecificCP_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificCP, P.val, T)*u"kJ/kg/K" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificCP, Region, P.val, T)*u"kJ/kg/K" end @@ -3295,16 +3346,7 @@ end """ function SpecificCV(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpecificCV, P, T) - elseif Region == :Region2 - return Region2(:SpecificCV, P, T) - elseif Region == :Region3 - return Region3(:SpecificCV, P, T) - elseif Region == :Region5 - return Region5(:SpecificCV, P, T) - end + return property_PT(:SpecificCV, Region, P, T) end @@ -3321,16 +3363,7 @@ function SpecificCV(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpecificCV, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region2 - return Region2(:SpecificCV, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region3 - return Region3(:SpecificCV, P.val, T.val)*u"kJ/kg/K" - elseif Region == :Region5 - return Region5(:SpecificCV, P.val, T.val)*u"kJ/kg/K" - end + return property_PT(:SpecificCV, Region, P.val, T.val)*u"kJ/kg/K" end @@ -3346,20 +3379,8 @@ end """ function SpecificCV_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpecificCV, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpecificCV, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpecificCV, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpecificCV, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpecificCV, Region, P, T) end @@ -3376,20 +3397,8 @@ function SpecificCV_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpecificCV, Region, P.val, T)*u"kJ/kg/K" end @@ -3405,20 +3414,8 @@ end """ function SpecificCV_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpecificCV, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpecificCV, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpecificCV, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpecificCV, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpecificCV, Region, P, T) end @@ -3435,20 +3432,8 @@ function SpecificCV_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpecificCV, P.val, T)*u"kJ/kg/K" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpecificCV, Region, P.val, T)*u"kJ/kg/K" end @@ -3461,16 +3446,7 @@ end """ function SpeedOfSound(P, T) Region = RegionID(P, T) - - if Region == :Region1 - return Region1(:SpeedOfSound, P, T) - elseif Region == :Region2 - return Region2(:SpeedOfSound, P, T) - elseif Region == :Region3 - return Region3(:SpeedOfSound, P, T) - elseif Region == :Region5 - return Region5(:SpeedOfSound, P, T) - end + return property_PT(:SpeedOfSound, Region, P, T) end @@ -3487,16 +3463,7 @@ function SpeedOfSound(P::Q1, T::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID(P.val, T.val) - - if Region == :Region1 - return Region1(:SpeedOfSound, P.val, T.val)*u"m/s" - elseif Region == :Region2 - return Region2(:SpeedOfSound, P.val, T.val)*u"m/s" - elseif Region == :Region3 - return Region3(:SpeedOfSound, P.val, T.val)*u"m/s" - elseif Region == :Region5 - return Region5(:SpeedOfSound, P.val, T.val)*u"m/s" - end + return property_PT(:SpeedOfSound, Region, P.val, T.val)*u"m/s" end @@ -3508,20 +3475,8 @@ end units of m/s via Uniful.jl.""" function SpeedOfSound_Ph(P, h) Region = RegionID_Ph(P, h) - - if Region == :Region1 - T = Region1_TPh(P, h) - return Region1(:SpeedOfSound, P, T) - elseif Region == :Region2a - T = Region2a_TPh(P, h) - return Region2(:SpeedOfSound, P, T) - elseif Region == :Region2b - T = Region2b_TPh(P, h) - return Region2(:SpeedOfSound, P, T) - elseif Region == :Region2c - T = Region2c_TPh(P, h) - return Region2(:SpeedOfSound, P, T) - end + T = _Temperature_Ph(Region, P, h) + return property_PT(:SpeedOfSound, Region, P, T) end @@ -3538,20 +3493,8 @@ function SpeedOfSound_Ph(P::Q1, h::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ph(P.val, h.val) - - if Region == :Region1 - T = Region1_TPh(P.val, h.val) - return Region1(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2a - T = Region2a_TPh(P.val, h.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2b - T = Region2b_TPh(P.val, h.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2c - T = Region2c_TPh(P.val, h.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - end + T = _Temperature_Ph(Region, P.val, h.val) + return property_PT(:SpeedOfSound, Region, P.val, T)*u"m/s" end @@ -3564,20 +3507,8 @@ end """ function SpeedOfSound_Ps(P, s) Region = RegionID_Ps(P, s) - - if Region == :Region1 - T = Region1_TPs(P, s) - return Region1(:SpeedOfSound, P, T) - elseif Region == :Region2a - T = Region2a_TPs(P, s) - return Region2(:SpeedOfSound, P, T) - elseif Region == :Region2b - T = Region2b_TPs(P, s) - return Region2(:SpeedOfSound, P, T) - elseif Region == :Region2c - T = Region2c_TPs(P, s) - return Region2(:SpeedOfSound, P, T) - end + T = _Temperature_Ps(Region, P, s) + return property_PT(:SpeedOfSound, Region, P, T) end @@ -3594,20 +3525,8 @@ function SpeedOfSound_Ps(P::Q1, s::Q2) where Q1 <: Quantity where Q2 <: Quantity end Region = RegionID_Ps(P.val, s.val) - - if Region == :Region1 - T = Region1_TPs(P.val, s.val) - return Region1(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2a - T = Region2a_TPs(P.val, s.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2b - T = Region2b_TPs(P.val, s.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - elseif Region == :Region2c - T = Region2c_TPs(P.val, s.val) - return Region2(:SpeedOfSound, P.val, T)*u"m/s" - end + T = _Temperature_Ps(Region, P.val, s.val) + return property_PT(:SpeedOfSound, Region, P.val, T)*u"m/s" end @@ -3633,7 +3552,7 @@ function SatDensL(T) return ρc*(1 + b[1]*τ^(1/3) + b[2]*τ^(2/3) + b[3]*τ^(5/3) + b[4]*τ^(16/3) + b[5]*τ^(43/3) + b[6]*τ^(110/3)) else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end #SatDensL @@ -3659,7 +3578,7 @@ function SatDensL(T::Q) where Q <: Quantity return ρc*(1 + b[1]*τ^(1/3) + b[2]*τ^(2/3) + b[3]*τ^(5/3) + b[4]*τ^(16/3) + b[5]*τ^(43/3) + b[6]*τ^(110/3))*u"kg/m^3" else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end #SatDensL @@ -3687,7 +3606,7 @@ function SatDensV(T) return ρc*exp(c[1]*τ^(2/6) + c[2]*τ^(4/6) + c[3]*τ^(8/6) + c[4]*τ^(18/6) + c[5]*τ^(37/6) + c[6]*τ^(71/6)) else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end #SatDensV @@ -3714,7 +3633,7 @@ function SatDensV(T::Q) where Q <: Quantity return ρc*exp(c[1]*τ^(2/6) + c[2]*τ^(4/6) + c[3]*τ^(8/6) + c[4]*τ^(18/6) + c[5]*τ^(37/6) + c[6]*τ^(71/6))*u"kg/m^3" else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end #SatDensV @@ -3739,9 +3658,9 @@ function SatHL(T) θ = T/Tc α = α0*(dα + d[1]*θ^(-19) + d[2]*θ + d[3]*θ^4.5 + d[4]*θ^5 + d[5]*θ^54.5) - return (α + T/SatDensL(T)*1e6*ForwardDiff.derivative(Psat2, T))/1000.0 + return (α + T/SatDensL(T)*1e6*∂Psat∂T(T))/1000.0 else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end @@ -3765,9 +3684,9 @@ function SatHL(T::Q) where Q <: Quantity θ = T.val/Tc α = α0*(dα + d[1]*θ^(-19) + d[2]*θ + d[3]*θ^4.5 + d[4]*θ^5 + d[5]*θ^54.5) - return (α + T.val/SatDensL(T.val)*1e6*ForwardDiff.derivative(Psat2, T.val))/1000.0*u"kJ/kg" + return (α + T.val/SatDensL(T.val)*1e6*∂Psat∂T(T.val))/1000.0*u"kJ/kg" else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end @@ -3792,9 +3711,9 @@ function SatHV(T) θ = T/Tc α = α0*(dα + d[1]*θ^(-19) + d[2]*θ + d[3]*θ^4.5 + d[4]*θ^5 + d[5]*θ^54.5) - return (α + T/SatDensV(T)*1e6*ForwardDiff.derivative(Psat2, T))/1000.0 + return (α + T/SatDensV(T)*1e6*∂Psat∂T(T))/1000.0 else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end @@ -3818,12 +3737,28 @@ function SatHV(T::Q) where Q <: Quantity θ = T.val/Tc α = α0*(dα + d[1]*θ^(-19) + d[2]*θ + d[3]*θ^4.5 + d[4]*θ^5 + d[5]*θ^54.5) - return (α + T.val/SatDensV(T.val)*1e6*ForwardDiff.derivative(Psat2, T.val))/1000.0*u"kJ/kg" + return (α + T.val/SatDensV(T.val)*1e6*∂Psat∂T(T.val))/1000.0*u"kJ/kg" else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end +function ϕS(T) + α0 = 1000 + ϕ0 = α0/Tc + + dϕ = 2319.5246 + d = ( -5.651_349_98E-8, + 2690.666_31, + 127.287_297, + -135.003_439, + 0.981_825_814 + ) + + θ = T/Tc + ϕ = ϕ0*(dϕ + 19/20*d[1]*θ^(-20) + d[2]*log(θ) + 9/7*d[3]*θ^3.5+ 5/4*d[4]*θ^4 + 109/107*d[5]*θ^53.5) +end + """ SatSL @@ -3834,22 +3769,10 @@ end """ function SatSL(T) if T3 ≤ T ≤ Tc - α0 = 1000 - ϕ0 = α0/Tc - - dϕ = 2319.5246 - d = [ -5.651_349_98E-8, - 2690.666_31, - 127.287_297, - -135.003_439, - 0.981_825_814 - ] - - θ = T/Tc - ϕ = ϕ0*(dϕ + 19/20*d[1]*θ^(-20) + d[2]*log(θ) + 9/7*d[3]*θ^3.5 + 5/4*d[4]*θ^4 + 109/107*d[5]*θ^53.5) - return (ϕ + 1/SatDensL(T)*1e6*ForwardDiff.derivative(Psat2, T))/1000.0 + ϕ = ϕS(T) + return (ϕ + 1/SatDensL(T)*1e6*∂Psat∂T(T))/1000.0 else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end @@ -3860,29 +3783,10 @@ function SatSL(T::Q) where Q <: Quantity catch throw(UnitsError(T, "Invalid input units.")) end - - if T3 ≤ T.val ≤ Tc - α0 = 1000 - ϕ0 = α0/Tc - - dϕ = 2319.5246 - d = [ -5.651_349_98E-8, - 2690.666_31, - 127.287_297, - -135.003_439, - 0.981_825_814 - ] - - θ = T.val/Tc - ϕ = ϕ0*(dϕ + 19/20*d[1]*θ^(-20) + d[2]*log(θ) + 9/7*d[3]*θ^3.5 - + 5/4*d[4]*θ^4 + 109/107*d[5]*θ^53.5) - return (ϕ + 1/SatDensL(T.val)*1e6*ForwardDiff.derivative(Psat2, T.val))/1000.0*u"kJ/kg/K" - else - throw(DomainError(T, "Temperature not between tripple and critical points")) - end + T3 ≤ T.val ≤ Tc || throw(DomainError(T, "Temperature not between triple and critical points")) + return SatSL(T.val)*u"kJ/kg/K" end - """ SatSV @@ -3892,23 +3796,10 @@ end """ function SatSV(T) if T3 ≤ T ≤ Tc - α0 = 1000 - ϕ0 = α0/Tc - - dϕ = 2319.5246 - d = [ -5.651_349_98E-8, - 2690.666_31, - 127.287_297, - -135.003_439, - 0.981_825_814 - ] - - θ = T/Tc - ϕ = ϕ0*(dϕ + 19/20*d[1]*θ^(-20) + d[2]*log(θ) + 9/7*d[3]*θ^3.5 - + 5/4*d[4]*θ^4 + 109/107*d[5]*θ^53.5) - return (ϕ + 1/SatDensV(T)*1e6*ForwardDiff.derivative(Psat2, T))/1000.0 + ϕ = ϕS(T) + return (ϕ + 1/SatDensV(T)*1e6*∂Psat∂T(T))/1000.0 else - throw(DomainError(T, "Temperature not between tripple and critical points")) + throw(DomainError(T, "Temperature not between triple and critical points")) end end @@ -3919,26 +3810,8 @@ function SatSV(T::Q) where Q <: Quantity catch throw(UnitsError(T, "Invalid input units.")) end - - if T3 ≤ T.val ≤ Tc - α0 = 1000 - ϕ0 = α0/Tc - - dϕ = 2319.5246 - d = [ -5.651_349_98E-8, - 2690.666_31, - 127.287_297, - -135.003_439, - 0.981_825_814 - ] - - θ = T.val/Tc - ϕ = ϕ0*(dϕ + 19/20*d[1]*θ^(-20) + d[2]*log(θ) + 9/7*d[3]*θ^3.5 - + 5/4*d[4]*θ^4 + 109/107*d[5]*θ^53.5) - return (ϕ + 1/SatDensV(T.val)*1e6*ForwardDiff.derivative(Psat2, T.val))/1000.0*u"kJ/kg/K" - else - throw(DomainError(T, "Temperature not between tripple and critical points")) - end + T3 ≤ T.val ≤ Tc || throw(DomainError(T, "Temperature not between triple and critical points")) + return SatSV(T.val)*u"kJ/kg/K" end From 374d5d1ed38cb78d81c5c95b14888af7ad4087b9 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 00:29:11 -0400 Subject: [PATCH 04/11] add region5 to RegionID_Ph and RegionID_Ps, add tests --- src/SteamTables.jl | 26 +++++++++++++++----------- test/testreg3.jl | 8 ++++++++ test/testreg5.jl | 8 ++++++++ 3 files changed, 31 insertions(+), 11 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index 8fcd252..44bfa73 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -1175,10 +1175,9 @@ function Region3_TPh(P, h) hhigh = Region3(:SpecificH, P, Thigh) end T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) - ρ = Region3_ρPT(P, T) Told = T for i in 1:100 - h_i = Region3_ρ(:SpecificH, ρ, T) + h_i = Region3(:SpecificH, P, T) if hlow < h_i hlow = h_i Tlow = T @@ -1191,7 +1190,6 @@ function Region3_TPh(P, h) if abs(T - Told) < 1e-12 return T end - ρ = Region3_ρPT(P, T) end throw(error("Region3_TPh: temperature iterations failed to converge.")) end @@ -1226,7 +1224,7 @@ function Region3_TPs(P, s) ρ = Region3_ρPT(P, T) Told = T for i in 1:100 - s_i = Region3_ρ(:SpecificS, ρ, T) + s_i = Region3(:SpecificS, P, T) if slow < s_i slow = s_i Tlow = T @@ -1239,7 +1237,6 @@ function Region3_TPs(P, s) if abs(T - Told) < 1e-12 return T end - ρ = Region3_ρPT(P, T, ρ) end throw(error("Region3_TPs: temperature iterations failed to converge.")) end @@ -2178,8 +2175,9 @@ function Region5(Output::Symbol, P, T) end function Region5_TPh(P, h) - hhighT = Region5(:SpecificH, P, 1073.15) - hlowT = Region5(:SpecificH, P, 2273.15) + Thigh,Tlow = (1073.15,2273.15) + hhighT = Region5(:SpecificH, P, Tlow) + hlowT = Region5(:SpecificH, P, Thigh) hhigh,hlow = minmax(hhighT,hlowT) T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) Told = T @@ -2202,9 +2200,9 @@ function Region5_TPh(P, h) end function Region5_TPs(P, s) - shighT = Region5(:SpecificH, P, 1073.15) - slowT = Region5(:SpecificH, P, 2273.15) - shigh,slow = minmax(shighT,slowT) + Thigh,Tlow = (1073.15,2273.15) + shigh = Region5(:SpecificS, P, Thigh) + slow = Region5(:SpecificS, P, Tlow) T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) Told = T for i in 1:100 @@ -2347,7 +2345,7 @@ function RegionID_Ph(P, h)::Symbol end elseif B23(:P, P) ≤ T ≤ 1073.15 return :Region2c #, Region2(:SpecificH, P, T) # Return forward h for consistency check - elseif h <= Region3_ρ(:SpecificH, P, B23(:P, P)) + elseif h <= Region3(:SpecificH, P, B23(:P, P)) return :Region3 else throw(DomainError((P, T), "Pressure/Temperature outside valid ranges.")) @@ -2361,6 +2359,8 @@ function RegionID_Ph(P, h)::Symbol end elseif B23(:P, P) ≤ T ≤ 1073.15 return :Region2b #, Region2(:SpecificH, P, T) # Return forward h for consistency check + elseif P < 50 && Region5(:SpecificH, P, 1073.15) <= h <= Region5(:SpecificH, P, 2273.15) + return :Region5 else throw(DomainError((P, T), "Pressure/Temperature outside valid ranges.")) end @@ -2429,6 +2429,8 @@ function RegionID_Ps(P, s)::Symbol end elseif B23(:P, P) ≤ T ≤ 1073.15 return :Region2c #, Region2(:SpecificS, P, T) # Return forward s for consistency check + elseif s <= Region3(:SpecificS, P, B23(:P, P)) + return :Region3 else throw(DomainError((P, s), "Pressure/entropy not in valid ranges.")) end @@ -2441,6 +2443,8 @@ function RegionID_Ps(P, s)::Symbol end elseif B23(:P, P) ≤ T ≤ 1073.15 return :Region2b #, Region2(:SpecificS, P, T) # Return forward s for consistency check + elseif P < 50 && Region5(:SpecificS, P, 1073.15) <= s <= Region5(:SpecificS, P, 2273.15) + return :Region5 else throw(DomainError((P, s), "Pressure/entropy not in valid ranges.")) diff --git a/test/testreg3.jl b/test/testreg3.jl index a3f581a..99a589b 100644 --- a/test/testreg3.jl +++ b/test/testreg3.jl @@ -51,3 +51,11 @@ @test SpeedOfSound(0.783_095_639E2, 750.0) ≈ 0.760_696_041E3 @test SpeedOfSound(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.760_696_041E3u"m/s" + +#issue #23 +let P = 30,T = 650 + H = SpecificH(P,T) + S = SpecificS(P,T) + @test S ≈ SpecificS_Ph(P,H) + @test H ≈ SpecificH_Ps(P,S) +end diff --git a/test/testreg5.jl b/test/testreg5.jl index c56c65c..8a64cc1 100644 --- a/test/testreg5.jl +++ b/test/testreg5.jl @@ -51,3 +51,11 @@ @test SpeedOfSound(30.0, 2000.0) ≈ 0.106_736_948E4 @test SpeedOfSound(30.0u"MPa", 2000.0u"K") ≈ 0.106_736_948E4u"m/s" + +#issue #23 +let P = 30,T = 2000 + H = SpecificH(P,T) + S = SpecificS(P,T) + @test S ≈ SpecificS_Ph(P,H) + @test H ≈ SpecificH_Ps(P,S) +end \ No newline at end of file From c5a7c4df13cc5e751d133502a21350d43408a67c Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 00:50:52 -0400 Subject: [PATCH 05/11] fix tolerance in Ps solvers, add tests in files --- src/SteamTables.jl | 10 ++++++++-- test/runtests.jl | 36 ++++++++++++++++++++++-------------- test/testreg3.jl | 8 -------- test/testreg3Px.jl | 6 ++++++ test/testreg5.jl | 8 -------- test/testreg5Px.jl | 6 ++++++ 6 files changed, 42 insertions(+), 32 deletions(-) create mode 100644 test/testreg3Px.jl create mode 100644 test/testreg5Px.jl diff --git a/src/SteamTables.jl b/src/SteamTables.jl index 44bfa73..3c23acf 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -1223,7 +1223,10 @@ function Region3_TPs(P, s) T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) ρ = Region3_ρPT(P, T) Told = T + s_i = (slow + shigh)/2 + sold = s_i for i in 1:100 + sold = s_i s_i = Region3(:SpecificS, P, T) if slow < s_i slow = s_i @@ -1234,7 +1237,7 @@ function Region3_TPs(P, s) end Told = T T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - if abs(T - Told) < 1e-12 + if abs(s_i - sold) < 1e-12 return T end end @@ -2205,7 +2208,10 @@ function Region5_TPs(P, s) slow = Region5(:SpecificS, P, Tlow) T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) Told = T + s_i = (shigh + slow)/2 + sold = s_i for i in 1:100 + sold = s_i s_i = Region5(:SpecificS, P, T) if slow < s_i slow = s_i @@ -2216,7 +2222,7 @@ function Region5_TPs(P, s) end Told = T T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - if abs(T - Told) < 1e-12 + if abs(s_i - sold) < 1e-12 return T end end diff --git a/test/runtests.jl b/test/runtests.jl index d63c3b2..576c0cb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,51 +12,51 @@ end const TestDigits = 6 # number of significant digits to test selected functions to -@testset "Region 1, forwards equations for P and T " begin +@testset "Region 1 - forwards equations for P and T " begin include("testreg1.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 1, backwards equations for P and h " begin +@testset "Region 1 - backwards equations for P and h " begin include("testreg1Ph.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 1, backwards equations for P and s " begin +@testset "Region 1 - backwards equations for P and s " begin include("testreg1Ps.jl") end -@testset "Region 2, forwards equations for P and T " begin +@testset "Region 2 - forwards equations for P and T " begin include("testreg2.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 2a, backwards equations for P and h " begin +@testset "Region 2a - backwards equations for P and h " begin include("testreg2aPh.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 2b, backwards equations for P and h " begin +@testset "Region 2b - backwards equations for P and h " begin include("testreg2bPh.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 2c, backwards equations for P and h " begin +@testset "Region 2c - backwards equations for P and h " begin include("testreg2cPh.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 2a, backwards equations for P and s " begin +@testset "Region 2a - backwards equations for P and s " begin include("testreg2aPs.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. #@test mysignif(SpecificG(8.0, 0.600_484_040E3), TestDigits) ≈ mysignif(SpecificG_Ps(8.0, 6.0), TestDigits) -@testset "Region 2b, backwards equations for P and s " begin +@testset "Region 2b - backwards equations for P and s " begin include("testreg2bPs.jl") end # The free energies have lower consistency between the forward and backwards equations. Check to fewer significant digits. -@testset "Region 2c, backwards equations for P and s " begin +@testset "Region 2c - backwards equations for P and s " begin include("testreg2cPs.jl") end @@ -64,18 +64,26 @@ end # for the pressure. This does introduce some small error, which means we should check to fewer significant digits. # For consistency the checks here are done to the same number of digits as for the free energies, although the # accuracy here will be better. -@testset "Region 3, forwards equations for P and T " begin +@testset "Region 3 - forwards equations for P and T " begin include("testreg3.jl") end -@testset "Region 4 - saturation pressure from temperature " begin +@testset "Region 3 - backwards equations for P and T " begin + include("testreg3Px.jl") +end + +@testset "Region 4 - saturation pressure from temperature " begin include("testreg4T.jl") end -@testset "Region 4 - saturation temperature from pressure " begin +@testset "Region 4 - saturation temperature from pressure " begin include("testreg4P.jl") end -@testset "Region 5 - forwards equations for P and T " begin +@testset "Region 5 - forwards equations for P and T " begin include("testreg5.jl") end + +@testset "Region 5 - backwards equations for P and T " begin + include("testreg5Px.jl") +end diff --git a/test/testreg3.jl b/test/testreg3.jl index 99a589b..a3f581a 100644 --- a/test/testreg3.jl +++ b/test/testreg3.jl @@ -51,11 +51,3 @@ @test SpeedOfSound(0.783_095_639E2, 750.0) ≈ 0.760_696_041E3 @test SpeedOfSound(0.783_095_639E2u"MPa", 750.0u"K") ≈ 0.760_696_041E3u"m/s" - -#issue #23 -let P = 30,T = 650 - H = SpecificH(P,T) - S = SpecificS(P,T) - @test S ≈ SpecificS_Ph(P,H) - @test H ≈ SpecificH_Ps(P,S) -end diff --git a/test/testreg3Px.jl b/test/testreg3Px.jl new file mode 100644 index 0000000..3cc1c78 --- /dev/null +++ b/test/testreg3Px.jl @@ -0,0 +1,6 @@ +P = 30 +T = 650 +H = SpecificH(P,T) +S = SpecificS(P,T) +@test S ≈ SpecificS_Ph(P,H) +@test H ≈ SpecificH_Ps(P,S) diff --git a/test/testreg5.jl b/test/testreg5.jl index 8a64cc1..c56c65c 100644 --- a/test/testreg5.jl +++ b/test/testreg5.jl @@ -51,11 +51,3 @@ @test SpeedOfSound(30.0, 2000.0) ≈ 0.106_736_948E4 @test SpeedOfSound(30.0u"MPa", 2000.0u"K") ≈ 0.106_736_948E4u"m/s" - -#issue #23 -let P = 30,T = 2000 - H = SpecificH(P,T) - S = SpecificS(P,T) - @test S ≈ SpecificS_Ph(P,H) - @test H ≈ SpecificH_Ps(P,S) -end \ No newline at end of file diff --git a/test/testreg5Px.jl b/test/testreg5Px.jl new file mode 100644 index 0000000..58ef52f --- /dev/null +++ b/test/testreg5Px.jl @@ -0,0 +1,6 @@ +let P = 30,T = 2000 + H = SpecificH(P,T) + S = SpecificS(P,T) + @test S ≈ SpecificS_Ph(P,H) + @test H ≈ SpecificH_Ps(P,S) +end From 3308ae4d426ccec2aa5eb6d7d13a355ed87770c8 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 12:47:17 -0400 Subject: [PATCH 06/11] use constant Tc,Pc instead of hardcoded value --- src/SteamTables.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index 3c23acf..fa19a68 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -1157,7 +1157,7 @@ end function Region3_TPh(P, h) Tlow = 623.15 Thigh = B23(:P, P) - if 16.529164252604478 <= P <= 22.064 + if 16.529164252604478 <= P <= Pc #we are in the saturation boundary Tsat = Tsat(P) hl,hv = SatHL(Tsat),SatHV(Tsat) @@ -1203,7 +1203,7 @@ end function Region3_TPs(P, s) Tlow = 623.15 Thigh = B23(:P, P) - if 16.529164252604478 <= P <= 22.064 + if 16.529164252604478 <= P <= Pc #we are in the saturation boundary Tsat = Tsat(P) sl,sv = SatSL(Tsat),SatSV(Tsat) @@ -1216,12 +1216,11 @@ function Region3_TPs(P, s) slow = Region3(:SpecificS, P, Tlow) shigh = sl end - else #22.064 < p <= 100 + else #Pc < p <= 100 slow = Region3(:SpecificS, P, Tlow) shigh = Region3(:SpecificS, P, Thigh) end T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - ρ = Region3_ρPT(P, T) Told = T s_i = (slow + shigh)/2 sold = s_i @@ -1999,7 +1998,7 @@ function Region3(Output::Symbol, P, T) end function Region3_ρ0(P, T) - if T <= 647.96 #we are inside saturation. use the saturation volume as initial guess + if T <= Tc #we are inside saturation. use the saturation volume as initial guess if 16.529164252604478 <= P <= Psat(T) #gas phase #gas phase ρ0 = 1000*P/(R*T) From 3458296d558bee30467700636e1608fa8c201d19 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 15:27:54 -0400 Subject: [PATCH 07/11] add bisection step if secant is out of bounds --- src/SteamTables.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index fa19a68..b17897f 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -1187,6 +1187,7 @@ function Region3_TPh(P, h) end Told = T T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) if abs(T - Told) < 1e-12 return T end @@ -1236,6 +1237,7 @@ function Region3_TPs(P, s) end Told = T T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) if abs(s_i - sold) < 1e-12 return T end @@ -2194,6 +2196,7 @@ function Region5_TPh(P, h) end Told = T T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) + Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) if abs(T - Told) < 1e-12 return T end @@ -2221,6 +2224,7 @@ function Region5_TPs(P, s) end Told = T T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) + Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) if abs(s_i - sold) < 1e-12 return T end From ee6d346d134907a751d352f1f4338430feb680a4 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 15:40:19 -0400 Subject: [PATCH 08/11] RegionID_Ph and RegionID_Ps: Region 1 is only in the liquid phase --- src/SteamTables.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index b17897f..a134ce3 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -2330,7 +2330,7 @@ function RegionID_Ph(P, h)::Symbol T = Region1_TPh(P, h) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) + if P ≥ Psat(T) && h <= SatHL(T) #Region1 is only the liquid phase return :Region1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check end end @@ -2415,7 +2415,7 @@ function RegionID_Ps(P, s)::Symbol T = Region1_TPs(P, s) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) + if P ≥ Psat(T) && s < SatSL(T) #Region 1 is only the liquid phase return :Region1 #, Region1(:SpecificS, P, T) # Return forward s for consistency check end end From 9814b6d3618237f5988435a20ffec8672ffee6b1 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Sun, 21 Apr 2024 16:49:35 -0400 Subject: [PATCH 09/11] RegionID_Ph/RegionID_Ps: use satXL(Region4(:P,P)) instead of satXL(Region1_T) --- src/SteamTables.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index a134ce3..7fbc8b3 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -2330,7 +2330,7 @@ function RegionID_Ph(P, h)::Symbol T = Region1_TPh(P, h) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) && h <= SatHL(T) #Region1 is only the liquid phase + if P ≥ Psat(T) && h <= SatHL(Region4(:P,P)) #Region1 is only the liquid phase return :Region1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check end end @@ -2415,7 +2415,7 @@ function RegionID_Ps(P, s)::Symbol T = Region1_TPs(P, s) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) && s < SatSL(T) #Region 1 is only the liquid phase + if P ≥ Psat(T) && s < SatSL(Region4(:P,P)) #Region 1 is only the liquid phase return :Region1 #, Region1(:SpecificS, P, T) # Return forward s for consistency check end end From d37e02c61f376bc0858dc111889ffe9263f89267 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Mon, 22 Apr 2024 17:37:00 -0400 Subject: [PATCH 10/11] use ITP solver, change crriteria for region1 again --- src/SteamTables.jl | 193 +++++++++++++++++++-------------------------- 1 file changed, 80 insertions(+), 113 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index 7fbc8b3..e29172d 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -78,7 +78,7 @@ const ρc = 322. #kg/m3 Critical density of water const T3 = 273.16 #K Triple point temperature of water const P3 = 611.657E-6 #MPa Triple point pressure of water const Mr = 18.01528 #kg/kmol Molecular weight of water - +const P_134 = 16.529164252604478 #intersection between region 1, 3 and 4. at T = 623.15 export Psat, Tsat, SatDensL, SatDensV, SatHL, SatHV, SatSL, SatSV, @@ -101,7 +101,46 @@ struct UnitsError <: Exception end Base.showerror(io::IO, e::UnitsError) = print(io, "UnitsError with ",e.var, ": ", e.message) +#ITP: interpolate-truncate-project. +#https://en.wikipedia.org/wiki/ITP_method +function itp_step(f::F,xa,xb,ya,yb,ϵ,κ₁,κ₂,j,nmid,nmax) where F + #calculate parameters + xmid = 0.5*(xa + xb) + r = ϵ*exp2(nmax-j) - 0.5*(xb - xa) + δ = κ₁*(xb - xa)^κ₂ + #interpolation: + xf = (yb*xa-ya*xb)/(yb-ya) + #truncate: + σ = sign(xmid-xf) + xt = δ <= abs(xmid-xf) ? xf + σ*δ : xmid + #project: + x_itp = abs(xt-xmid) <= r ? xt : xmid - σ*r + y_itp = f(x_itp) + #update interval + if y_itp > 0 + xb,yb = x_itp,y_itp + elseif y_itp < 0 + xa,ya = x_itp,y_itp + else + xa,xb = x_itp,x_itp + end + return xa,xb,ya,yb +end +function itp_find_zero(f::F,xa,xb;tol = 0.5e-12,max_iters = 100) where F + κ₁,κ₂ = 0.2,2.0 + n0 = 1 + ya,yb = f(xa),f(xb) + nmid = log2((xb-xa)/(2*tol)) + nmax = nmid + n0 + for i in 0:max_iters + xa,xb,ya,yb = itp_step(f::F,xa,xb,ya,yb,tol,κ₁,κ₂,i,nmid,nmax) + if xb - xa <= 2*tol + return 0.5*(xa+xb) + end + end + return zero(xb)/zero(xb) +end """ B23(InputType::Symbol, InputValue) @@ -1157,42 +1196,22 @@ end function Region3_TPh(P, h) Tlow = 623.15 Thigh = B23(:P, P) - if 16.529164252604478 <= P <= Pc + if P_134 <= P <= Pc #we are in the saturation boundary Tsat = Tsat(P) hl,hv = SatHL(Tsat),SatHV(Tsat) if hl <= h <= hv return Tsat elseif h > hv #gas phase, interpolate with h(T_high) - hlow = hv - hhigh = Region3(:SpecificH, P, Thigh) + Tlow = Tsat elseif h < hl #liquid phase, interpolate with h(T_low) - hlow = Region3(:SpecificH, P, Tlow) - hhigh = hl + Thigh = Tsat end - else #22.064 < p <= 100 - hlow = Region3(:SpecificH, P, Tlow) - hhigh = Region3(:SpecificH, P, Thigh) - end - T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) - Told = T - for i in 1:100 - h_i = Region3(:SpecificH, P, T) - if hlow < h_i - hlow = h_i - Tlow = T - else - hhigh = h_i - Thigh = T - end - Told = T - T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) - Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) - if abs(T - Told) < 1e-12 - return T - end - end - throw(error("Region3_TPh: temperature iterations failed to converge.")) + end #22.064 < p <= 100 + f(T) = Region3(:SpecificH, P, T) - h + T = itp_find_zero(f,Tlow, Thigh) + isnan(T) && throw(error("Region3_TPh: temperature iterations failed to converge.")) + return T end """ @@ -1204,45 +1223,22 @@ end function Region3_TPs(P, s) Tlow = 623.15 Thigh = B23(:P, P) - if 16.529164252604478 <= P <= Pc + if P_134 <= P <= Pc #we are in the saturation boundary Tsat = Tsat(P) sl,sv = SatSL(Tsat),SatSV(Tsat) if sl <= s <= sv return Tsat elseif s > sv #gas phase, interpolate with s(T_high) - slow = sv - shigh = Region3(:SpecificS, P, Thigh) + Tlow = Tsat elseif s < sl #liquid phase, interpolate with s(T_low) - slow = Region3(:SpecificS, P, Tlow) - shigh = sl - end - else #Pc < p <= 100 - slow = Region3(:SpecificS, P, Tlow) - shigh = Region3(:SpecificS, P, Thigh) - end - T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - Told = T - s_i = (slow + shigh)/2 - sold = s_i - for i in 1:100 - sold = s_i - s_i = Region3(:SpecificS, P, T) - if slow < s_i - slow = s_i - Tlow = T - else - shigh = s_i - Thigh = T - end - Told = T - T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) - if abs(s_i - sold) < 1e-12 - return T + Thigh = Tsat end end - throw(error("Region3_TPs: temperature iterations failed to converge.")) + f(T) = Region3(:SpecificS, P, T) - s + T = itp_find_zero(f,Tlow, Thigh) + isnan(T) && throw(error("Region3_TPs: temperature iterations failed to converge.")) + return T end """ @@ -2001,12 +1997,12 @@ end function Region3_ρ0(P, T) if T <= Tc #we are inside saturation. use the saturation volume as initial guess - if 16.529164252604478 <= P <= Psat(T) #gas phase + if P_134 <= P <= Psat(T) #gas phase #gas phase - ρ0 = 1000*P/(R*T) + return 1000*P/(R*T) else #liquid phase - ρ0 = SatDensL(T) + ρ0 = 574.6703980963142 end else #over critical point,it does not matter we start, but a high density point is recomended ρ0_liquid = 574.6703980963142 #volume at the intersection of regions 1-3-4 @@ -2179,57 +2175,19 @@ function Region5(Output::Symbol, P, T) end function Region5_TPh(P, h) - Thigh,Tlow = (1073.15,2273.15) - hhighT = Region5(:SpecificH, P, Tlow) - hlowT = Region5(:SpecificH, P, Thigh) - hhigh,hlow = minmax(hhighT,hlowT) - T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) - Told = T - for i in 1:100 - h_i = Region5(:SpecificH, P, T) - if hlow < h_i - hlow = h_i - Tlow = T - else - hhigh = h_i - Thigh = T - end - Told = T - T = Tlow + (Thigh - Tlow) / (hhigh - hlow) * (h - hlow) - Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) - if abs(T - Told) < 1e-12 - return T - end - end - throw(error("Region5_TPh: temperature iterations failed to converge.")) + Thigh,Tlow = 1073.15,2273.15 + f(T) = Region5(:SpecificH, P, T) - h + T = itp_find_zero(f, Tlow, Thigh) + isnan(T) && throw(error("Region5_TPh: temperature iterations failed to converge.")) + return T end function Region5_TPs(P, s) - Thigh,Tlow = (1073.15,2273.15) - shigh = Region5(:SpecificS, P, Thigh) - slow = Region5(:SpecificS, P, Tlow) - T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - Told = T - s_i = (shigh + slow)/2 - sold = s_i - for i in 1:100 - sold = s_i - s_i = Region5(:SpecificS, P, T) - if slow < s_i - slow = s_i - Tlow = T - else - shigh = s_i - Thigh = T - end - Told = T - T = Tlow + (Thigh - Tlow) / (shigh - slow) * (s - slow) - Tlow <= T <= Thigh || (T = 0.5*(Thigh + Tlow)) - if abs(s_i - sold) < 1e-12 - return T - end - end - throw(error("Region5_TPs: temperature iterations failed to converge.")) + Thigh,Tlow = 1073.15,2273.15 + f(T) = Region5(:SpecificS, P, T) - s + T = itp_find_zero(f, Tlow, Thigh) + isnan(T) && throw(error("Region5_TPs: temperature iterations failed to converge.")) + return T end @@ -2330,8 +2288,12 @@ function RegionID_Ph(P, h)::Symbol T = Region1_TPh(P, h) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) && h <= SatHL(Region4(:P,P)) #Region1 is only the liquid phase - return :Region1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check + if P ≥ Psat(T) + #Region1 is only the liquid phase + hlim = P <= P_134 ? SatHL(Region4(:P,P)) : Region1(:SpecificH, P,623.15) + if h <= hlim + return :Region1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check + end end end @@ -2415,7 +2377,12 @@ function RegionID_Ps(P, s)::Symbol T = Region1_TPs(P, s) if 273.15 ≤ T ≤ 623.15 # could be Region 1 - if P ≥ Psat(T) && s < SatSL(Region4(:P,P)) #Region 1 is only the liquid phase + if P ≥ Psat(T) + #Region1 is only the liquid phase + slim = P <= P_134 ? SatSL(Region4(:P,P)) : Region1(:SpecificS, P,623.15) + if s <= slim + return :Region1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check + end return :Region1 #, Region1(:SpecificS, P, T) # Return forward s for consistency check end end From 0a45cf5e39919d7cb165d780e405b3c20fa59daa Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Mon, 22 Apr 2024 17:46:34 -0400 Subject: [PATCH 11/11] Region5_TPs/Region5_TPh: fix inputs --- src/SteamTables.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/SteamTables.jl b/src/SteamTables.jl index e29172d..d0647b6 100644 --- a/src/SteamTables.jl +++ b/src/SteamTables.jl @@ -131,7 +131,8 @@ function itp_find_zero(f::F,xa,xb;tol = 0.5e-12,max_iters = 100) where F κ₁,κ₂ = 0.2,2.0 n0 = 1 ya,yb = f(xa),f(xb) - nmid = log2((xb-xa)/(2*tol)) + @assert ya < yb + nmid = ceil(log2(abs(xb-xa)/(2*tol))) nmax = nmid + n0 for i in 0:max_iters xa,xb,ya,yb = itp_step(f::F,xa,xb,ya,yb,tol,κ₁,κ₂,i,nmid,nmax) @@ -2175,7 +2176,7 @@ function Region5(Output::Symbol, P, T) end function Region5_TPh(P, h) - Thigh,Tlow = 1073.15,2273.15 + Tlow,Thigh = 1073.15,2273.15 f(T) = Region5(:SpecificH, P, T) - h T = itp_find_zero(f, Tlow, Thigh) isnan(T) && throw(error("Region5_TPh: temperature iterations failed to converge.")) @@ -2183,7 +2184,7 @@ function Region5_TPh(P, h) end function Region5_TPs(P, s) - Thigh,Tlow = 1073.15,2273.15 + Tlow,Thigh = 1073.15,2273.15 f(T) = Region5(:SpecificS, P, T) - s T = itp_find_zero(f, Tlow, Thigh) isnan(T) && throw(error("Region5_TPs: temperature iterations failed to converge."))