Прежде чем перейти к статье, хочу вам представить, экономическую онлайн игру Brave Knights, в которой вы можете играть и зарабатывать. Регистируйтесь, играйте и зарабатывайте!
При составлении библиотеки использованы документы Международной ассоциации по свойствам воды и водяного пара (МАСПВ, анг.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.
![](https://habrastorage.org/r/w780q1/getpro/habr/upload_files/b71/af3/fda/b71af3fda61c88ee3a6f3d3e6149554b.jpg)
1. Основные функции
Энергия_Гиббса(t As Double, p As Double, Optional reg) As Double, [Дж/кг]
Удельный_объем(t As Double, p As Double, Optional reg) As Double, [м3/кг]
Плотность(t As Double, p As Double, Optional reg) As Double, [кг/ м3]
Удельная_энергия(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-й области в качестве параметра имеют температуру в градусах Цельсия и плотность в кг/ м3. Данные функции полностью идентичны формулам IF97 и дают наиболее точные результаты в этой области.
3. Вспомогательные функции
Плотность3(t As Double, p As Double, Optional accuracy) As Double
Функция для 3-ей области осуществляет поиск плотности соответствующей полученным параметрам : температуре в градусах Цельсия и давлению в МПа. По умолчанию точность результата 3, что соответствует 3 знакам после запятой ±0,001 кг/ м3. В районе критической точки (P=22,064МПа и t=373,946°С) функция не может точно определить значение из-за того, что график зависимости ρ=f(p) становится практически горизонтальным и погрешность возрастает до 3-4% и более. Для повышения точности рекомендуется вводить более высокое значение точности. Максимально допустимое значение параметра ограничено accuracy=13 !!! Это вызвано точностью с которой оперирует VBA значениями типа double. При вводе точности 14 и более функция зависает.
Функция работает в диапазоне плотностей от 113,6 до 762,4 кг/ м3. В случае если результат выходит за указанные рамки функция возвращает ошибку неправильное значение #ЗНАЧ!
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. Но так как параметр μ2 принят равным единице, то в районе критической области (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, [кг/ м3]
Удельная_энтальпия_МИ(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