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

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

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

25.01.2023 10:22

hydroargerum 4 766 Источник

При составлении библиотеки использованы документы Международной ассоциации по свойствам воды и водяного пара (МАСПВ, анг.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

Комментарии (4)

ViktorAPT
25.01.2023 16:59
#25153990
+1
Я тоже делаю этот расчет. :) Вот вам ещё литература. Ещё стесняюсь спросить, для каких расчетов будете использовать? ;)
1. hydroargerum Автор
  25.01.2023 18:46
  #25154396
  Спасибо, знаком с этим справочником, начинал с него. В нём Александров и Григорьев в предисловии ссылаются на первоисточник, это те же документы МАСПВ (IAPWS), что и послужили основой для написания данной библиотеки функций. Только в первоисточнике хоть и на английском, но проще: там коэффициенты даны уже готовые для всех частных производных энергий Гиббса и Гельмгольца, а по справочнику нужно получать дифференцируя. Я пробовал это: долго, нудно и нужна внимательность, а в итоге получается то же самое.
  По вопросу "для чего", то использую для:
  
  проверки показаний приборов учёта тепловой энергии;
  
  гидравлических расчётов тепловых сетей и сетей горячего водоснабжения. Кстати, формулы плотности и энтальпии полностью идентичны приложению А ГОСТ Р ЕН 1434-1-2011 "Теплосчетчики. Часть 1. Общие требования"
  1. ViktorAPT
    25.01.2023 22:37
    #25155034
    А что вы думаете пор сооружения на сетях и прочую тепломеханику (мои схемы)? Я вообще сейчас интересуюсь паровыми преобразователями тепла. Если всякое такое интересно пишите в личку. (у меня низкий рейтинг могу комментировать 1раз в час) В ответ пришлю мои расчеты, не закончил ещё проверку. Как у вас с проверкой? Как находите удельный объём/плотность в 3 диапазоне?
    
    hydroargerum Автор
    26.01.2023 18:03
    #25158000
    В третьем диапазоне есть функция Плотность3(температура гр.Цельсия, давление МПа, точность знаков после запятой). А в ней банально методом половинного деления (колесо изобретать не стал).