Расчёт свойств воды и водяного пара

[VADR]

Недавно довелось столкнуться с такой задачей - расчёт тепловой энергии пара и воды в разных диапазонах давлений и температур. Помня такую книгу, как "Термодинамические свойства воды и водяного пара", Вукалович М.П., предполагал, что тут всё непросто и аппроксимировать кривые можно до потери пульса. Но откуда-то из закоулков памяти вылезло и другое: существуют приближённые (тем не менее, весьма точные) формулы расчёта. По информации, которую дал гугль, разработал эти формулы "The International Association for the Properties of Water and Steam", текст существует в открытом доступе и лежит здесь: http://www.iapws.org/relguide/IF97-Rev.pdf. Документ, по сравнению с Вукаловичем, "всего" на 49 страниц :).
Формулы на первый взгляд кажутся многоэтажными, но это только на первый взгляд. Важно точно вбить таблицы (и для этого язык, на котором вы будете это реализовывать, должен поддерживать массивы), а сами формулы... впрочем, ниже я приведу то, что у меня получилось - видно, что ничего принципиально сложного нет.
В моей ситуации я был ограничен в выборе языка программирования: целевая система поддерживает расчёты:
1. На своём внутреннем скриптовом языке (не поддерживает массивов и циклов).
2. В виде формул в Excel (назначаются ячейки для входящих данных и результата, можно прописать любые расчёты, но тупит это просто неимоверно - раз в секунду вызывать не получится никак).
3. COM-объект. Даже не стал влезать в это дело: мало того, что саму программу надо написать, её ещё надо специальным образом зарегистрировать, о ней и её регистрации надо не забывать при резервном копировании... Я на такие жертвы пойду только в том случае, когда других вариантов не останется совсем.
4. EXE-файл. Что-то типа предыдущего пункта, так что тоже - самый крайний случай.
5. VBScript. Массивы и циклы поддерживаются, текст функций хранится непосредственно в базе данных целевой системы и бэкапится вместе с ней без каких-либо дополнительных действий, текст функций всегда виден во встроенном редакторе.
В общем, хоть я и недолюбливаю Visual Basic вообще и VBScript в частности, в данном случае пришлось остановиться именно на нём. Это я тут налепил лирическое отступление, чтобы потом не отвечать на вопросы, почему именно VBS.
В общем, переходя от лирики к технике:
1. Из всех рассчитываемых параметров мне нужны были только удельный объём и энтальпия для воды и пара. Если посмотреть на функции в первоисточнике, нет ничего сложного на базе этих функций сделать аналогичные для остальных параметров (энтропия, внутренняя энергия, изохорная и изобарная теплоёмкости, скорость звука)
2. Для линии насыщения - нужны зависимость температуры от давления и наоборот.
3. Если есть возможность, следует использовать переменные двойной точности (double). По крайней мере, в тестовых примерах так рекомендуется.
4. На странице 4 документа есть график, где отмечены зоны, в которых применима та или иная формула. В моём случае в технологии нет пара и воды с давлением выше ~16МПа либо температурой выше 1073,15K, так что меня во всех этой красоте интересуют только регионы 1 и 2 и линия насыщения 4. Принцип расчёта в других регионах ничем не отличается, так что желающие могут реализовать эти формулы самостоятельно.
Итак, что у меня получилось:
1. Удельный объём воды (там, где температура ниже линии насыщения), по графику - регион 1:

Код: Выделить всё

Function wspV1(ByVal P, ByVal T)
	Dim Pi, Tau, RES, n, I, J, k, R
	R = 0.461526
	n = Array (0, 0.14632971213167, -0.84548187169114, -0.37563603672040E1, 0.33855169168385E1, -0.95791963387872, 0.15772038513228, -0.16616417199501E-1, 0.81214629983568E-3, 0.28319080123804E-3, -0.60706301565874E-3, -0.18990068218419E-1, -0.32529748770505E-1, -0.21841717175414E-1, -0.52838357969930E-4, -0.47184321073267E-3, -0.30001780793026E-3, 0.47661393906987E-4, -0.44141845330846E-5, -0.72694996297594E-15, -0.31679644845054E-4, -0.28270797985312E-5, -0.85205128120103E-9, -0.22425281908000E-5, -0.65171222895601E-6, -0.14341729937924E-12, -0.40516996860117E-6, -0.12734301741641E-8, -0.17424871230634E-9, -0.68762131295531E-18, 0.14478307828521E-19, 0.26335781662795E-22, -0.11947622640071E-22, 0.18228094581404E-23, -0.93537087292458E-25)
	I = Array (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32)
	J = Array (0, -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41)
 	Pi = P / 16.53 
	Tau = 1386 / T
	RES = 0
	for k = 1 to 34
		RES = RES + (-1 * n(k) * I(k) * ((7.1 - Pi) ^ (I(k)-1)) * ((Tau - 1.222)^J(k)))
	Next
	wspV1 = RES * Pi * R * T / (P * 1000)
end Function

2. Удельная энтальпия воды:

Код: Выделить всё

Function wspH1(ByVal P, ByVal T)
	Dim Pi, Tau, RES, n, I, J, k, R
	R = 0.461526
	n = Array (0, 0.14632971213167, -0.84548187169114, -0.37563603672040E1, 0.33855169168385E1, -0.95791963387872, 0.15772038513228, -0.16616417199501E-1, 0.81214629983568E-3, 0.28319080123804E-3, -0.60706301565874E-3, -0.18990068218419E-1, -0.32529748770505E-1, -0.21841717175414E-1, -0.52838357969930E-4, -0.47184321073267E-3, -0.30001780793026E-3, 0.47661393906987E-4, -0.44141845330846E-5, -0.72694996297594E-15, -0.31679644845054E-4, -0.28270797985312E-5, -0.85205128120103E-9, -0.22425281908000E-5, -0.65171222895601E-6, -0.14341729937924E-12, -0.40516996860117E-6, -0.12734301741641E-8, -0.17424871230634E-9, -0.68762131295531E-18, 0.14478307828521E-19, 0.26335781662795E-22, -0.11947622640071E-22, 0.18228094581404E-23, -0.93537087292458E-25)
	I = Array (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32)
	J = Array (0, -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41)
 	Pi = P / 16.53 
	Tau = 1386 / T
	RES = 0
	for k = 1 to 34
		RES = RES + ( n(k) * ((7.1 - Pi) ^ I(k)) * J(k) * ((Tau - 1.222)^(J(k)-1)))
	Next
	wspH1 = RES * Tau * R * T 
end Function

3. Удельный объём пара (температура выше линии насыщения), по графику - регион 2:

Код: Выделить всё

Function wspV2(ByVal P, ByVal T)
	Dim n0, J0, RES, R, Pi, Tau, k, n, I, J
	R = 0.461526
	n0 = Array (0, -0.96927686500217E1, 0.10086655968018E2, -0.56087911283020E-2, 0.71452738081455E-1, -0.40710498223928, 0.14240819171444E1, -0.43839511319450E1, -0.28408632460772, 0.21268463753307E-1)
    J0 = Array (0, 0, 1, -5, -4, -3, -2, -1, 2, 3)
	I = Array (0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24)
	J = Array (0, 0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58)
	n = Array (0, -0.17731742473213E-2, -0.17834862292358E-1, -0.45996013696365E-1, -0.57581259083432E-1, -0.50325278727930E-1, -0.33032641670203E-4, -0.18948987516315E-3, -0.39392777243355E-2, -0.43797295650573E-1, -0.26674547914087E-4, 0.20481737692309E-7, 0.43870667284435E-6, -0.32277677238570E-4, -0.15033924542148E-2, -0.40668253562649E-1, -0.78847309559367E-9, 0.12790717852285E-7, 0.48225372718507E-6, 0.22922076337661E-5, -0.16714766451061E-10, -0.21171472321355E-2, -0.23895741934104E2, -0.59059564324270E-17, -0.12621808899101E-5, -0.38946842435739E-1, 0.11256211360459E-10, -0.82311340897998E1, 0.19809712802088E-7, 0.10406965210174E-18, -0.10234747095929E-12, -0.10018179379511E-8, -0.80882908646985E-10, 0.10693031879409, -0.33662250574171, 0.89185845355421E-24, 0.30629316876232E-12, -0.42002467698208E-5, -0.59056029685639E-25, 0.37826947613457E-5, -0.12768608934681E-14, 0.73087610595061E-28, 0.55414715350778E-16, -0.94369707241210E-6)
	Pi = P
	Tau = 540 / T
	RES = 1 / Pi
	For k = 1 to 43
		RES = RES + (n(k) * I(k) * (Pi ^ (I(k) - 1)) * ((Tau - 0.5) ^ J(k)))
	Next
	wspV2 = RES * Pi * R * T / (P * 1000)
end Function

4. Удельная энтальпия пара:

Код: Выделить всё

Function wspH2(ByVal P, ByVal T)
	Dim n0, J0, RES, R, Pi, Tau, k, n, I, J
	R = 0.461526
	n0 = Array (0, -0.96927686500217E1, 0.10086655968018E2, -0.56087911283020E-2, 0.71452738081455E-1, -0.40710498223928, 0.14240819171444E1, -0.43839511319450E1, -0.28408632460772, 0.21268463753307E-1)
    J0 = Array (0, 0, 1, -5, -4, -3, -2, -1, 2, 3)
	I = Array (0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24)
	J = Array (0, 0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58)
	n = Array (0, -0.17731742473213E-2, -0.17834862292358E-1, -0.45996013696365E-1, -0.57581259083432E-1, -0.50325278727930E-1, -0.33032641670203E-4, -0.18948987516315E-3, -0.39392777243355E-2, -0.43797295650573E-1, -0.26674547914087E-4, 0.20481737692309E-7, 0.43870667284435E-6, -0.32277677238570E-4, -0.15033924542148E-2, -0.40668253562649E-1, -0.78847309559367E-9, 0.12790717852285E-7, 0.48225372718507E-6, 0.22922076337661E-5, -0.16714766451061E-10, -0.21171472321355E-2, -0.23895741934104E2, -0.59059564324270E-17, -0.12621808899101E-5, -0.38946842435739E-1, 0.11256211360459E-10, -0.82311340897998E1, 0.19809712802088E-7, 0.10406965210174E-18, -0.10234747095929E-12, -0.10018179379511E-8, -0.80882908646985E-10, 0.10693031879409, -0.33662250574171, 0.89185845355421E-24, 0.30629316876232E-12, -0.42002467698208E-5, -0.59056029685639E-25, 0.37826947613457E-5, -0.12768608934681E-14, 0.73087610595061E-28, 0.55414715350778E-16, -0.94369707241210E-6)
	Pi = P
	Tau = 540 / T
	RES = 0
	For k = 1 to 9
		RES = RES + n0(k) * J0(k) * (Tau ^ (J0(k) - 1))
	Next
	For k = 1 to 43
		RES = RES + (n(k) * (Pi ^ I(k)) * J(k) * ((Tau - 0.5) ^ (J(k) - 1)))
	Next
	wspH2 = RES * Tau * R * T
end Function

5. Температура на линии насыщения в зависимости от давления и наоборот:

Код: Выделить всё

Function wspTsat(ByVal P)
	Dim Beta, D, E, F, G, n
	n = Array (0, 0.11670521452767E4, -0.72421316703206E6, -0.17073846940092E2, 0.12020824702470E5, -0.32325550322333E7, 0.14915108613530E2, -0.48232657361591E4, 0.40511340542057E6, -0.23855557567849, 0.65017534844798E3)
	Beta = Sqr( Sqr (P))
	E = Beta^2 + (n(3)*Beta) + n(6)
	F = (n(1)*Beta^2) + (n(4)*Beta) + n(7)
	G = (n(2)*Beta^2) + (n(5)*Beta) + n(8)
	D = 2*G / (-1*F - Sqr(F^2 - 4*E*G))
    wspTsat = (n(10) + D - Sqr((n(10)+D)^2 - 4*(n(9)+n(10)*D)))/2
end Function

Код: Выделить всё

Function wspPsat(ByVal T)
	Dim Theta, A, B, C, n
	n = Array (0, 0.11670521452767E4, -0.72421316703206E6, -0.17073846940092E2, 0.12020824702470E5, -0.32325550322333E7, 0.14915108613530E2, -0.48232657361591E4, 0.40511340542057E6, -0.23855557567849, 0.65017534844798E3)
	Theta = T + (n(9) / (T - n(10)))
	A = Theta^2 + (n(1) * Theta) + n(2)
	B = (n(3) * Theta^2) + (n(4) * Theta) + n(5)
	C = (n(6) * Theta^2) + (n(7) * Theta) + n(8)
	wspPsat = (2*C/(-1*B + Sqr(B^2 - 4*A*C))) ^ 4
end Function

Единицы измерения - можно в первоисточнике посмотреть. Там же есть тестовые данные для расчёта. Кстати, в формулах для удельного объёма в предпоследней строке я дополнительно разделил результат на 1000, хотя в описании формул такого не было. В общем, есть там какая-то нестыковка с единицами измерения, в итоге, чтобы совпали результаты по тестовым выражениям, пришлось их поделить на 1000 :).

Вот как-то так, собственно. Кому вдруг понадобится - пользуйтесь на здоровье.

[hell_boy]

Кроме IAPWS IF-97 существуют также ГСССД 187-89 и МИ 2451-98. Макросы VB лежали на сайте IAPWS

steam2010.xls

[VADR]

существуют также ГСССД 187-89 и МИ 2451-98

Блин, не знал.

Макросы VB лежали на сайте IAPWS

А этого дела что-то так и не нашёл. Или плохо искал. Знаю, что есть софтина под названием WaterSteamPro (и у них даже имена функций начинаются так же, как у меня - с wsp, однако под последней "p" я имел в виду "properties", а не "pro") и она сделана на основе того же самого документа, но в моём случае применить её вряд ли было возможно.