Библиотека функций для расчета свойств воды и водяного пара

При составлении библиотеки использованы документы Международной ассоциации по свойствам воды и водяного пара (МАСПВ, анг.IASPW):

Revised Release on the IAPWS Industrial Formulation 1997 For the Thermodynamic Properties of Water and Steam. August 2007. (Далее в тексте IF97)
Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance. September 2008.

1. Основные функции

Энергия_Гиббса(t As Double, p As Double, Optional reg) As Double, [Дж/кг]

Удельный_объем(t As Double, p As Double, Optional reg) As Double, [м³/кг]

Плотность(t As Double, p As Double, Optional reg) As Double, [кг/ м³]

Удельная_энергия(t As Double, p As Double, Optional reg) As Double, [Дж/кг]

Удельная_энтропия(t As Double, p As Double, Optional reg) As Double, [Дж/кг]

Удельная_энтальпия(t As Double, p As Double, Optional reg) As Double, [Дж/кг]

Теплоемкость_изобарная(t As Double, p As Double, Optional reg) As Double, [Дж/(кг∙°С)]

Теплоемкость_изохорная(t As Double, p As Double, Optional reg) As Double, [Дж/(кг∙°С)]

Скорость_звука(t As Double, p As Double, Optional reg) As Double, [м/с]

Все функции в качестве параметра имеют температуру в градусах Цельсия и давление в МПа. Функции самостоятельно определяют область (с помощью функции Region) и формулу, которую необходимо применить в этом регионе.

Если имеется необходимость, то функции допускают ручное выбор региона. Это может потребоваться, например, для расчета параметров в метастабильном регионе. Значения параметра reg это числа с номерами регионов соответствующих областей на диаграмме приведенной в формуляции IF97. Кроме этого для области метастабильных паров параметр reg=21 (так как эта область перекрывается областью 1 – вода).

В области 3 в IF97 в в качестве аргумента вместо давления используется плотность поэтому функции определяют плотность посредством функции Плотность3 с точностью до 3 знаков после запятой. Вследствие этого точность может несколько снижаться, а в районе критической точки погрешность может достигать 3-4%. Поэтому целесообразней использование специальных функций 3-й области.

Замечание

функция Энергия_Гиббса в 3-ей области вместо собственно энергии Гиббса возвращает в качестве результата энергию Гельмгольца.

2. Функции 3-й области

Энергия_Гельмгольца(t As Double, ro As Double) As Double, [Дж/кг]

Давление3(t As Double, ro As Double) As Double, [Па]

Удельная_энергия3(t As Double, ro As Double) As Double, [Дж/кг]

Удельная_энтропия3(t As Double, ro As Double) As Double, [Дж/кг]

Удельная_энтальпия3(t As Double, ro As Double) As Double, [Дж/кг]

Теплоемкость_изобарная3(t As Double, ro As Double) As Double, [Дж/(кг∙°С)]

Теплоемкость_изохорная3(t As Double, ro As Double) As Double, [Дж/(кг∙°С)]

Скорость_звука3(t As Double, ro As Double) As Double, [м/с]

Функции 3-й области в качестве параметра имеют температуру в градусах Цельсия и плотность в кг/ м³. Данные функции полностью идентичны формулам IF97 и дают наиболее точные результаты в этой области.

3. Вспомогательные функции

Плотность3(t As Double, p As Double, Optional accuracy) As Double

Функция для 3-ей области осуществляет поиск плотности соответствующей полученным параметрам : температуре в градусах Цельсия и давлению в МПа. По умолчанию точность результата 3, что соответствует 3 знакам после запятой ±0,001 кг/ м³. В районе критической точки (P=22,064МПа и t=373,946°С) функция не может точно определить значение из-за того, что график зависимости ρ=f(p) становится практически горизонтальным и погрешность возрастает до 3-4% и более. Для повышения точности рекомендуется вводить более высокое значение точности. Максимально допустимое значение параметра ограничено accuracy=13 !!! Это вызвано точностью с которой оперирует VBA значениями типа double. При вводе точности 14 и более функция зависает.

Функция работает в диапазоне плотностей от 113,6 до 762,4 кг/ м³. В случае если результат выходит за указанные рамки функция возвращает ошибку неправильное значение #ЗНАЧ!

4. Вспомогательные функции

Температура_насыщения(p As Double) As Double, [°С]

Давление_насыщения(t As Double) As Double, [МПа]

Функции в описывают график кривой линии насыщения и возвращают соответствующее их температуре или давлению второй параметр. Аргументы и результаты функций выражены: давление в МПа, температура в градусах Цельсия.

Температура_границы(p As Double) As Double, [°С]

Давление_границы(t As Double) As Double, [МПа]

Функции в описывают график кривой линии ограничивающей область 2 (пар) от областей 1 и 3. Параметры функций идентичны функциям линии насыщения.

Замечание

При температурах до 350°С функции совпадают с функциями линии насыщения.

Region(t As Double, p As Double) As Byte

Функция возвращает номер области соответствующей введенным давлению и температуре. Номера областей соответствуют диаграмме приведенной в IF97

JF(t As Double, p As Double, Trigger As Byte, reg As Byte) As Double

Функция возвращает для всех областей значение приведенных характерных энергий региона reg. Для 3-й области это функция свободной энергии Гельмгольца , а для остальных областей энергия Гиббса

Функция имеет 4 обязательных параметра:

1. Температура в градусах Цельсия

2. Давление в МПа

3. Trigger имеет следующие возможные варианты:

a. Собственно сама функция – 0

b. Частная производная по температуре – 1

c. Вторая частная производная по температуре – 2

d. Смешанная частная производная по температуре и давлению (или плотности для 3-го региона) — 3

e. Частная производная давлению (или плотности для 3-го региона) — 4

f. Вторая частная производная давлению (или плотности для 3-го региона) — 5

4. Регион имеет допустимые значения – 1,2,3,4,5,21

Вязкость(t As Double, p As Double, Optional reg) As Double, [Па∙с]

Функция возвращает значение динамической вязкости выраженную в Па∙с посчитанную по формуле приведенной в «Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance» Сентябрь 2008. Но так как параметр μ₂принят равным единице, то в районе критической области (P=22,064МПа и t=373,946°С) функция не применима. В остальных областях функция дает значения с приемлемой точностью до 1-2%. В остальном область определения функции полностью совпадает с областью определения функции плотности для документа IF97.

5. Функции МИ 2412-97

В версии 1.1. добавлены функции плотности и энтальпии вычисляемые по формулам П.1 и П.2 приведенным в МИ 2412-97 «Водяные системы теплоснабжения. Уравнения измерений тепловой энергии и количества теплоносителя».

Плотность_МИ(t As Double, p As Double) As Double, [кг/ м³]

Удельная_энтальпия_МИ(t As Double, p As Double) As Double, [Дж/кг]

Attribute VB_Name = "IF97_1_1" Option Base 1 Public N1, I1, J1 'I phase Public N02, J02, Nr2, Ir2, Jr2 'II phase Public N02m, Nr2m, Ir2m, Jr2m 'metastable-vapor region Public N3, I3, J3 'III phase Public N05, J05, Nr5, Ir5, Jr5 'V phase Public Vi0, VHi, Vi, Vj, VHij ' Viscosity Const non = 0, dt = 1, dtt = 2, dtp = 3, dp = 4, dpp = 5 ' Триггеры функции энергии Гиббса Const R = 461.526, Default_accuracy = 3 Sub Auto_Open()     MsgBox "При написании функций использованы :" & vbCrLf & "1. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam (Редакция 2007 года)." & vbCrLf & "2. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance." & vbCrLf & "3. МИ 2412-97 Водяные системы теплоснабжения уравнения измерений тепловой энергии и количества теплоносителя", , "ООО «Рога и Копыта»"     N1 = Array(0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26)     I1 = Array(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)     J1 = Array(-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)     N02 = Array(-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)     N02m = Array(-9.6937268393049, 10.087275970006, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307)     J02 = Array(0, 1, -5, -4, -3, -2, -1, 2, 3)     Nr2 = Array(-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -0.082311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07)     Ir2 = Array(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)     Jr2 = Array(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)     Nr2m = Array(-7.3362260186506E-03, -0.088223831943146, -0.072334555213245, -4.0813178534455E-03, 2.0097803380207E-03, -0.053045921898642, -0.007619040908697, -6.3498037657313E-03, -0.086043093028588, 0.007532158152277, -7.9238375446139E-03, -2.2888160778447E-04, -0.002645650148281)     Ir2m = Array(1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5)     Jr2m = Array(0, 2, 5, 11, 1, 7, 16, 4, 16, 7, 10, 9, 10)     N3 = Array(-15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05)     I3 = Array(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)     J3 = Array(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)     N05 = Array(-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917)     J05 = Array(0, 1, -3, -2, -1, 2)     Nr5 = Array(1.5736404855259E-03, 9.0153761673944E-04, -5.0270077677648E-03, 2.2440037409485E-06, -4.1163275453471E-06, 3.7919454822955E-08)     Ir5 = Array(1, 1, 1, 2, 2, 3)     Jr5 = Array(1, 2, 3, 3, 9, 7)     VHi = Array(1.67752, 2.20462, 0.6366564, -0.241605)     Vi = Array(0, 1, 2, 3, 0, 1, 2, 3, 5, 0, 1, 2, 3, 4, 0, 1, 0, 3, 4, 3, 5)     Vj = Array(0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 6, 6)     VHij = Array(0.520094, 0.0850895, -1.08374, -0.289555, 0.222531, 0.999115, 1.88797, 1.26613, 0.120573, -0.281378, -0.906851, -0.772479, -0.489837, -0.25704, 0.161913, 0.257399, -0.0325372, 0.0698452, 0.00872102, -0.00435673, -0.000593264) End Sub Function Region(t As Double, p As Double) As Byte     Dim pf, tf As Double     If ((t <= 590) And (p <= 100)) Then         pf = Давление_границы(t)     End If     Region = 0     If ((590 < t) And (t <= 800) And (0 <= p) And (p <= 100)) Then         Region = 2     ElseIf ((800 < t) And (t <= 2000) And (0 <= p) And (p <= 50)) Then         Region = 5     ElseIf ((0 <= t) And (t <= 350) And (pf < p) And (p <= 100)) Then         Region = 1     ElseIf ((350 < t) And (t <= 590) And (pf < p) And (p <= 100)) Then         Region = 3     ElseIf ((0 <= t) And (t <= 590) And (0 <= p) And (p < pf)) Then         Region = 2     ElseIf ((0 <= t) And (t <= 590) And (p = pf)) Then         Region = 4     End If End Function Function Плотность3(t As Double, p As Double, Optional accuracy) As Double Dim ro_max As Double, ro_mid As Double, ro_min As Double, p_mid As Double, faccuracy As Double     If t < 373.946 Then         If Давление_насыщения(t) < p Then             ro_min = 322             ro_max = 762.4         Else             ro_min = 113.6             ro_max = 322         End If     Else             ro_min = 113.6             ro_max = 762.4     End If     ro_mid = (ro_max + ro_min) / 2     p_mid = ro_mid ^ 2 * R * (t + 273.15) * JF(t, ro_mid, dp, 3) / 322000000 - p     If IsMissing(accuracy) Then faccuracy = 10 ^ (-Default_accuracy) Else faccuracy = 10 ^ (-accuracy)     Do While (Abs(p_mid) > faccuracy) And (Abs(ro_max - ro_mid) > faccuracy)         If p_mid < 0 Then             ro_min = ro_mid         Else             ro_max = ro_mid         End If         ro_mid = (ro_max + ro_min) / 2         p_mid = ro_mid ^ 2 * R * (t + 273.15) * JF(t, ro_mid, dp, 3) / 322000000 - p     Loop     If Abs(p_mid) > faccuracy Then Плотность3 = 1 / 0 Else Плотность3 = ro_mid End Function Function Энергия_Гельмгольца(t As Double, ro As Double) As Double     Энергия_Гельмгольца = JF(t, ro, non, 3) * (t + 273.15) * R End Function Function Давление3(t As Double, ro As Double) As Double     Давление3 = ro ^ 2 * R * (t + 273.15) * JF(t, ro, dp, 3) / 322 End Function Function Удельная_энергия3(t As Double, ro As Double) As Double     Удельная_энергия3 = R * 647.096 * JF(t, ro, dt, 3) End Function Function Удельная_энтропия3(t As Double, ro As Double) As Double     Удельная_энтропия3 = R * (647.096 / (t + 273.15) * JF(t, ro, dt, 3) - JF(t, ro, non, 3)) End Function Function Удельная_энтальпия3(t As Double, ro As Double) As Double     Удельная_энтальпия3 = R * (647.096 * JF(t, ro, dt, 3) + (t + 273.15) * JF(t, ro, dp, 3) * ro / 322) End Function Function Теплоемкость_изобарная3(t As Double, ro As Double) As Double     Dim tf As Double, rof As Double, fp As Double     tf = 647.096 / (t + 273.15)     rof = ro / 322     fp = JF(t, ro, dp, 3)     Теплоемкость_изобарная3 = (-tf ^ 2 * JF(t, ro, dtt, 3) + (fp - tf * JF(t, ro, dtp, 3)) ^ 2 / (2 * fp / rof + JF(t, ro, dpp, 3))) * R End Function Function Теплоемкость_изохорная3(t As Double, ro As Double) As Double     Теплоемкость_изохорная3 = -(647.096 / (t + 273.15)) ^ 2 * JF(t, ro, dtt, 3) * R End Function Function Скорость_звука3(t As Double, ro As Double) As Double     Dim tf As Double, rof As Double, fp As Double     tf = 647.096 / (t + 273.15)     rof = ro / 322     fp = JF(t, ro, dp, 3)     Скорость_звука3 = (R * (t + 273.15) * rof ^ 2 * (2 * fp / rof + JF(t, ro, dpp, 3) - (fp - tf * JF(t, ro, dtp, 3)) ^ 2 / (tf ^ 2 * JF(t, ro, dtt, 3)))) ^ 0.5 End Function Function Энергия_Гиббса(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1, 2, 4, 5, 21         Энергия_Гиббса = JF(t, p, non, Trigger) * (t + 273.15) * R     Case Is = 3         ro = Плотность3(t, p)         Энергия_Гиббса = JF(t, ro, non, 3) * (t + 273.15) * R     End Select End Function Function Удельный_объем(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Удельный_объем = (t + 273.15) * R * JF(t, p, dp, Trigger) / 16530000     Case Is = 2, 4, 5, 21         Удельный_объем = (t + 273.15) * R * JF(t, p, dp, Trigger) / 1000000     Case Is = 3         Удельный_объем = 1 / Плотность3(t, p)     End Select End Function Function Плотность(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Плотность = 1 / ((t + 273.15) * R * JF(t, p, dp, Trigger) / 16530000)     Case Is = 2, 4, 5, 21         Плотность = 1 / ((t + 273.15) * R * JF(t, p, dp, Trigger) / 1000000)     Case Is = 3         Плотность = Плотность3(t, p)     End Select End Function Function Удельная_энергия(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Удельная_энергия = R * (1386 * JF(t, p, dt, Trigger) - p / 16.53 * (t + 273.15) * JF(t, p, dp, Trigger))     Case Is = 2, 4, 21         Удельная_энергия = R * (540 * JF(t, p, dt, Trigger) - p * (t + 273.15) * JF(t, p, dp, Trigger))     Case Is = 5         Удельная_энергия = R * (1000 * JF(t, p, dt, Trigger) - p * (t + 273.15) * JF(t, p, dp, Trigger))     Case Is = 3         ro = Плотность3(t, p)         Удельная_энергия = R * 647.096 * JF(t, ro, dt, Trigger)     End Select End Function Function Удельная_энтропия(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Удельная_энтропия = R * (1386 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))     Case Is = 2, 4, 21         Удельная_энтропия = R * (540 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))     Case Is = 5         Удельная_энтропия = R * (1000 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))     Case Is = 3         ro = Плотность3(t, p)         Удельная_энтропия = R * (647.096 / (t + 273.15) * JF(t, ro, dt, Trigger) - JF(t, ro, non, Trigger))     End Select End Function Function Удельная_энтальпия(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Удельная_энтальпия = R * 1386 * JF(t, p, dt, Trigger)     Case Is = 2, 4, 21         Удельная_энтальпия = R * 540 * JF(t, p, dt, Trigger)     Case Is = 5         Удельная_энтальпия = R * 1000 * JF(t, p, dt, Trigger)     Case Is = 3         ro = Плотность3(t, p)         Удельная_энтальпия = R * (647.096 * JF(t, ro, dt, Trigger) + (t + 273.15) * JF(t, ro, dp, Trigger) * ro / 322)     End Select End Function Function Теплоемкость_изобарная(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double, tf As Double, rof As Double, fp As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         Теплоемкость_изобарная = -R * (1386 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)     Case Is = 2, 4, 21         Теплоемкость_изобарная = -R * (540 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)     Case Is = 5         Теплоемкость_изобарная = -R * (1000 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)     Case Is = 3         ro = Плотность3(t, p)         tf = 647.096 / (t + 273.15)         rof = ro / 322         fp = JF(t, ro, dp, Trigger)         Теплоемкость_изобарная = (-tf ^ 2 * JF(t, ro, dtt, Trigger) + (fp - tf * JF(t, ro, dtp, Trigger)) ^ 2 / (2 * fp / rof + JF(t, ro, dpp, Trigger))) * R     End Select End Function Function Теплоемкость_изохорная(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double, tf As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         tf = 1386 / (t + 273.15)         Теплоемкость_изохорная = R * (-tf ^ 2 * JF(t, p, dtt, Trigger) + (JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / JF(t, p, dpp, Trigger))     Case Is = 2, 4, 21         tf = 540 / (t + 273.15)         Теплоемкость_изохорная = R * (-tf ^ 2 * JF(t, p, dtt, Trigger) + (JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / JF(t, p, dpp, Trigger))     Case Is = 5         tf = 1000 / (t + 273.15)         Теплоемкость_изохорная = R * (-tf ^ 2 * JF(t, p, dtt, Trigger) + (JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / JF(t, p, dpp, Trigger))     Case Is = 3         ro = Плотность3(t, p)         Теплоемкость_изохорная = -(647.096 / (t + 273.15)) ^ 2 * JF(t, ro, dtt, 3) * R     End Select End Function Function Скорость_звука(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double, tf As Double, rof As Double, fp As Double     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     Select Case Trigger     Case Is = 1         tf = 1386 / (t + 273.15)         Скорость_звука = Sqr(R * (t + 273.15) * JF(t, p, dp, Trigger) ^ 2 / ((JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / (tf ^ 2 * JF(t, p, dtt, Trigger)) - JF(t, p, dpp, Trigger)))     Case Is = 2, 4, 21         tf = 540 / (t + 273.15)         Скорость_звука = Sqr(R * (t + 273.15) * JF(t, p, dp, Trigger) ^ 2 / ((JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / (tf ^ 2 * JF(t, p, dtt, Trigger)) - JF(t, p, dpp, Trigger)))     Case Is = 5         tf = 1000 / (t + 273.15)         Скорость_звука = Sqr(R * (t + 273.15) * JF(t, p, dp, Trigger) ^ 2 / ((JF(t, p, dp, Trigger) - tf * JF(t, p, dtp, Trigger)) ^ 2 / (tf ^ 2 * JF(t, p, dtt, Trigger)) - JF(t, p, dpp, Trigger)))     Case Is = 3         ro = Плотность3(t, p)         tf = 647.096 / (t + 273.15)         rof = ro / 322         fp = JF(t, ro, dp, Trigger)         Скорость_звука = (R * (t + 273.15) * rof ^ 2 * (2 * fp / rof + JF(t, ro, dpp, Trigger) - (fp - tf * JF(t, ro, dtp, Trigger)) ^ 2 / (tf ^ 2 * JF(t, ro, dtt, Trigger)))) ^ 0.5     End Select End Function Function Температура_насыщения(p As Double) As Double     Dim pf As Double, E As Double, F As Double, G As Double, D As Double     Const N1 = 1167.0521452767, n2 = -724213.16703206, N3 = -17.073846940092, n4 = 12020.82470247, n5 = -3232555.0322333, n6 = 14.91510861353, n7 = -4823.2657361591, n8 = 405113.40542057, n9 = -0.23855557567849, n10 = 650.17534844798     pf = p ^ 0.25     E = pf ^ 2 + N3 * pf + n6     F = N1 * pf ^ 2 + n4 * pf + n7     G = n2 * pf ^ 2 + n5 * pf + n8     D = 2 * G / (-F - (F ^ 2 - 4 * E * G) ^ 0.5)     Температура_насыщения = (n10 + D - ((n10 + D) ^ 2 - 4 * (n9 + n10 * D)) ^ 0.5) / 2 - 273.15 End Function Function Давление_насыщения(t As Double) As Double     Dim tf As Double, a As Double, B As Double, C As Double     Const N1 = 1167.0521452767, n2 = -724213.16703206, N3 = -17.073846940092, n4 = 12020.82470247, n5 = -3232555.0322333, n6 = 14.91510861353, n7 = -4823.2657361591, n8 = 405113.40542057, n9 = -0.23855557567849, n10 = 650.17534844798     tf = (t + 273.15) + n9 / (t + 273.15 - n10)     a = tf ^ 2 + N1 * tf + n2     B = N3 * tf ^ 2 + n4 * tf + n5     C = n6 * tf ^ 2 + n7 * tf + n8     Давление_насыщения = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4 End Function Function Температура_границы(p As Double) As Double     If p < 16.5292 Then         Dim pf As Double, E As Double, F As Double, G As Double, D As Double         Const N1 = 1167.0521452767, n2 = -724213.16703206, N3 = -17.073846940092, n4 = 12020.82470247, n5 = -3232555.0322333, n6 = 14.91510861353, n7 = -4823.2657361591, n8 = 405113.40542057, n9 = -0.23855557567849, n10 = 650.17534844798         pf = p ^ 0.25         E = pf ^ 2 + N3 * pf + n6         F = N1 * pf ^ 2 + n4 * pf + n7         G = n2 * pf ^ 2 + n5 * pf + n8         D = 2 * G / (-F - (F ^ 2 - 4 * E * G) ^ 0.5)         Температура_границы = (n10 + D - ((n10 + D) ^ 2 - 4 * (n9 + n10 * D)) ^ 0.5) / 2 - 273.15     Else         Температура_границы = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5 - 273.15     End If End Function Function Давление_границы(t As Double) As Double     If t < 350 Then         Dim tf As Double, a As Double, B As Double, C As Double         Const N1 = 1167.0521452767, n2 = -724213.16703206, N3 = -17.073846940092, n4 = 12020.82470247, n5 = -3232555.0322333, n6 = 14.91510861353, n7 = -4823.2657361591, n8 = 405113.40542057, n9 = -0.23855557567849, n10 = 650.17534844798         tf = (t + 273.15) + n9 / (t + 273.15 - n10)         a = tf ^ 2 + N1 * tf + n2         B = N3 * tf ^ 2 + n4 * tf + n5         C = n6 * tf ^ 2 + n7 * tf + n8         Давление_границы = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4     Else         Давление_границы = 348.05185628969 - 1.1671859879975 * (t + 273.15) + 1.0192970039326E-03 * (t + 273.15) ^ 2     End If End Function Function Вязкость(t As Double, p As Double, Optional reg) As Double     Dim Trigger As Byte, ro As Double, tf As Double, rof As Double, fp As Double     Dim i As Integer     If IsMissing(reg) Then Trigger = Region(t, p) Else Trigger = reg     tf = (t + 273.15) / 647.096     rof = Плотность(t, p, Trigger) / 322     mu0 = 0     For i = 1 To 4         mu0 = mu0 + VHi(i) / tf ^ (i - 1)     Next     mu0 = 100 * tf ^ 0.5 / mu0     mu1 = 0     For i = 1 To 21         mu1 = mu1 + VHij(i) * (1 / tf - 1) ^ Vi(i) * (rof - 1) ^ Vj(i)     Next     mu1 = Exp(rof * mu1)     mu2 = 1     Вязкость = mu0 * mu1 * mu2 * 0.000001 End Function Function JF(t As Double, p As Double, Trigger As Byte, reg As Byte) As Double Dim i As Integer, tf As Double, pf As Double, rof As Double Select Case reg Case 1     pf = p / 16.53     tf = 1386 / (t + 273.15)     JF = 0     Select Case Trigger     Case non         For i = 1 To 34             JF = JF + N1(i) * (7.1 - pf) ^ I1(i) * (tf - 1.222) ^ J1(i)         Next     Case dt         For i = 1 To 34             JF = JF + N1(i) * J1(i) * (7.1 - pf) ^ I1(i) * (tf - 1.222) ^ (J1(i) - 1)         Next     Case dtt         For i = 1 To 34             JF = JF + N1(i) * J1(i) * (J1(i) - 1) * (7.1 - pf) ^ I1(i) * (tf - 1.222) ^ (J1(i) - 2)         Next     Case dtp         For i = 1 To 34             JF = JF - N1(i) * I1(i) * J1(i) * (7.1 - pf) ^ (I1(i) - 1) * (tf - 1.222) ^ (J1(i) - 1)         Next     Case dp         For i = 1 To 34             JF = JF - N1(i) * I1(i) * (7.1 - pf) ^ (I1(i) - 1) * (tf - 1.222) ^ J1(i)         Next     Case dpp         For i = 1 To 34             JF = JF + N1(i) * I1(i) * (I1(i) - 1) * (7.1 - pf) ^ (I1(i) - 2) * (tf - 1.222) ^ J1(i)         Next     End Select Case 2, 4     pf = p     tf = 540 / (t + 273.15)     Select Case Trigger     Case non         JF = Log(p)         For i = 1 To 9             JF = JF + N02(i) * tf ^ J02(i)         Next         For i = 1 To 43             JF = JF + Nr2(i) * pf ^ Ir2(i) * (tf - 0.5) ^ Jr2(i)         Next     Case dt         JF = 0         For i = 1 To 9             JF = JF + N02(i) * J02(i) * tf ^ (J02(i) - 1)         Next         For i = 1 To 43             JF = JF + Nr2(i) * pf ^ Ir2(i) * Jr2(i) * (tf - 0.5) ^ (Jr2(i) - 1)         Next     Case dtt         JF = 0         For i = 1 To 9             JF = JF + N02(i) * J02(i) * (J02(i) - 1) * tf ^ (J02(i) - 2)         Next         For i = 1 To 43             JF = JF + Nr2(i) * pf ^ Ir2(i) * Jr2(i) * (Jr2(i) - 1) * (tf - 0.5) ^ (Jr2(i) - 2)         Next     Case dtp         JF = 0         For i = 1 To 43             JF = JF + Nr2(i) * Ir2(i) * pf ^ (Ir2(i) - 1) * Jr2(i) * (tf - 0.5) ^ (Jr2(i) - 1)         Next     Case dp         JF = 1 / pf         For i = 1 To 43             JF = JF + Nr2(i) * Ir2(i) * pf ^ (Ir2(i) - 1) * (tf - 0.5) ^ Jr2(i)         Next     Case dpp         JF = -1 / pf ^ 2         For i = 1 To 43             JF = JF + Nr2(i) * Ir2(i) * (Ir2(i) - 1) * pf ^ (Ir2(i) - 2) * (tf - 0.5) ^ Jr2(i)         Next     End Select Case 21     pf = p     tf = 540 / (t + 273.15)     Select Case Trigger     Case non         JF = Log(p)         For i = 1 To 9             JF = JF + N02m(i) * tf ^ J02(i)         Next         For i = 1 To 13             JF = JF + Nr2m(i) * pf ^ Ir2m(i) * (tf - 0.5) ^ Jr2m(i)         Next     Case dt         JF = 0         For i = 1 To 9             JF = JF + N02m(i) * J02(i) * tf ^ (J02(i) - 1)         Next         For i = 1 To 13             JF = JF + Nr2m(i) * pf ^ Ir2m(i) * Jr2m(i) * (tf - 0.5) ^ (Jr2m(i) - 1)         Next     Case dtt         JF = 0         For i = 1 To 9             JF = JF + N02m(i) * J02(i) * (J02(i) - 1) * tf ^ (J02(i) - 2)         Next         For i = 1 To 13             JF = JF + Nr2m(i) * pf ^ Ir2m(i) * Jr2m(i) * (Jr2m(i) - 1) * (tf - 0.5) ^ (Jr2m(i) - 2)         Next     Case dtp         JF = 0         For i = 1 To 13             JF = JF + Nr2m(i) * Ir2m(i) * pf ^ (Ir2m(i) - 1) * Jr2m(i) * (tf - 0.5) ^ (Jr2m(i) - 1)         Next     Case dp         JF = 1 / pf         For i = 1 To 13             JF = JF + Nr2m(i) * Ir2m(i) * pf ^ (Ir2m(i) - 1) * (tf - 0.5) ^ Jr2m(i)         Next     Case dpp         JF = -1 / pf ^ 2         For i = 1 To 13             JF = JF + Nr2m(i) * Ir2m(i) * (Ir2m(i) - 1) * pf ^ (Ir2m(i) - 2) * (tf - 0.5) ^ Jr2m(i)         Next     End Select Case 3     rof = p / 322     tf = 647.096 / (t + 273.15)     Select Case Trigger     Case non         JF = 1.0658070028513 * Log(rof)         For i = 1 To 39             JF = JF + N3(i) * rof ^ I3(i) * tf ^ J3(i)         Next     Case dt         JF = 0         For i = 1 To 39             JF = JF + N3(i) * rof ^ I3(i) * J3(i) * tf ^ (J3(i) - 1)         Next     Case dtt         JF = 0         For i = 1 To 39             JF = JF + N3(i) * rof ^ I3(i) * J3(i) * (J3(i) - 1) * tf ^ (J3(i) - 2)         Next     Case dtp         JF = 0         For i = 1 To 39             JF = JF + N3(i) * I3(i) * rof ^ (I3(i) - 1) * J3(i) * tf ^ (J3(i) - 1)         Next     Case dp         JF = 1.0658070028513 / rof         For i = 1 To 39             JF = JF + N3(i) * I3(i) * rof ^ (I3(i) - 1) * tf ^ J3(i)         Next     Case dpp         JF = -1.0658070028513 / rof ^ 2         For i = 1 To 39             JF = JF + N3(i) * I3(i) * (I3(i) - 1) * rof ^ (I3(i) - 2) * tf ^ J3(i)         Next     End Select Case 5     pf = p     tf = 1000 / (t + 273.15)     Select Case Trigger     Case non         JF = Log(p)         For i = 1 To 6             JF = JF + N05(i) * tf ^ J05(i)         Next         For i = 1 To 6             JF = JF + Nr5(i) * pf ^ Ir5(i) * tf ^ Jr5(i)         Next     Case dt         JF = 0         For i = 1 To 6             JF = JF + N05(i) * J05(i) * tf ^ (J05(i) - 1)         Next         For i = 1 To 6             JF = JF + Nr5(i) * pf ^ Ir5(i) * Jr5(i) * tf ^ (Jr5(i) - 1)         Next     Case dtt         JF = 0         For i = 1 To 6             JF = JF + N05(i) * J05(i) * (J05(i) - 1) * tf ^ (J05(i) - 2)         Next         For i = 1 To 6             JF = JF + Nr5(i) * pf ^ Ir5(i) * Jr5(i) * (Jr5(i) - 1) * tf ^ (Jr5(i) - 2)         Next     Case dtp         JF = 0         For i = 1 To 6             JF = JF + Nr5(i) * Ir5(i) * pf ^ (Ir5(i) - 1) * Jr5(i) * tf ^ (Jr5(i) - 1)         Next     Case dp         JF = 1 / pf         For i = 1 To 6             JF = JF + Nr5(i) * Ir5(i) * pf ^ (Ir5(i) - 1) * tf ^ Jr5(i)         Next     Case dpp         JF = -1 / pf ^ 2         For i = 1 To 6             JF = JF + Nr5(i) * Ir5(i) * (Ir5(i) - 1) * pf ^ (Ir5(i) - 2) * tf ^ Jr5(i)         Next     End Select End Select End Function Function Плотность_МИ(t As Double, p As Double) As Double     tau = (t + 273.15) / 647.14     Pi = p / 22.064     Плотность_МИ = 1000 / (114.332 * tau - 431.6382 + 706.5474 / tau - 641.9127 / tau ^ 2 + 349.4417 / tau ^ 3 - 113.8191 / tau ^ 4 + 20.5199 / tau ^ 5 - 1.578507 / tau ^ 6 + Pi * (-3.117072 + 6.589303 / tau - 5.210142 / tau ^ 2 + 1.819096 / tau ^ 3 - 0.2365448 / tau ^ 4) + Pi ^ 2 * (-6.417443 * tau + 19.84842 - 24.00174 / tau + 14.21655 / tau ^ 2 - 4.13194 / tau ^ 3 + 0.4721637 / tau ^ 4)) End Function Function Удельная_энтальпия_МИ(t As Double, p As Double) As Double     tau = (t + 273.15) / 647.14     Pi = p / 22.064     Удельная_энтальпия_МИ = (7809.096 * tau - 13868.72 + 12725.22 / tau - 6370.893 / tau ^ 2 + 1595.86 / tau ^ 3 - 159.9064 / tau ^ 4 + Pi * (9.488789 / tau + 1) + Pi ^ 2 * (-148.1135 * tau + 224.3027 - 111.4602 / tau + 18.15823 / tau ^ 2)) * 1000 End Function

ссылка на оригинал статьи https://habr.com/ru/post/712656/