Library of functions for calculation of water and water vapor properties (repost)

1

Vote

Original post

https://habr.com/ru/articles/712656/ (Russian)

Description

These documents of the International Association for the Properties of Water and Water Vapor (IASPW, ang. IASPW) were used in compiling the library:

• Revised Release on the IAPWS Industrial Formulation 1997 For the Thermodynamic Properties of Water and Steam. August 2007. (Further on in the text IF97)
• Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance. September 2008.

Regions and equations of IAPWS-IF97:


1. Main functions

• Gibbs_energy(t As Double, p As Double, Optional reg) As Double, [J/kg]
• Specific_volume(t As Double, p As Double, Optional reg) As Double, [m³/kg]
• Density(t As Double, p As Double, Optional reg) As Double, [kg/m³]
• Specific_energy(t As Double, p As Double, Optional reg) As Double, [J/kg]
• Specific_entropy(t As Double, p As Double, Optional reg) As Double, [J/kg]
• Specific_enthalpy(t As Double, p As Double, Optional reg) As Double, [J/kg]
• Heat_capacity_isobaric(t As Double, p As Double, Optional reg) As Double, [J/(kg∙°C)]
• Heat_capacity_isochoric(t As Double, p As Double, Optional reg) As Double, [J/(kg∙°C)]
• Speed_sound(t As Double, p As Double, Optional reg) As Double, [m/s]

All functions have temperature in degrees Celsius and pressure in MPa as a parameter. The functions independently define the region (using the Region function) and the formula to be applied in this region.

If necessary, the functions allow manual selection of the region. This may be necessary, for example, to calculate parameters in a metastable zone. The values of the reg parameter are numbers with the region numbers of the corresponding regions in the diagram given in the IF97 formula. In addition, for the metastable vapor the parameter reg=21 (since this region overlaps with region 1 - water).

In region 3 in IF97, density is used as an argument instead of pressure, so the functions determine the density by means of the function Density3 with an accuracy of 3 decimal places. As a consequence, the accuracy may decrease slightly, and in the area of the critical point the error may reach 3-4%. Therefore, it is more appropriate to use special functions of the 3rd region.

Note: function Gibbs_energy in the 3rd region returns Helmholtz energy as a result instead of Gibbs energy itself.

2. Functions of the 3rd region

• Helmholtz_energy(t As Double, ro As Double) As Double, [J/kg]
• Pressure3(t As Double, ro As Double) As Double, [Pa]
• Specific_energy3(t As Double, ro As Double) As Double, [J/kg]
• Specific_entropy3(t As Double, ro As Double) As Double, [J/kg]
• Specific_enthalpy3(t As Double, ro As Double) As Double, [J/kg]
• Heat_capacity_isobaric3(t As Double, ro As Double) As Double, [J/(kg∙°C)]
• Heat_capacity_isochoric3(t As Double, ro As Double) As Double, [J/(kg∙°C)]
• Speed_Sound3(t As Double, ro As Double) As Double, [m/s]

The functions of the 3rd region have temperature in degrees Celsius and density in kg/m3 as a parameter. These functions are completely identical to the IF97 formulas and give the most accurate results in this region.

3. Auxiliary functions

• Density3(t As Double, p As Double, Optional accuracy) As Double

The function for the 3rd region searches for the density corresponding to the obtained parameters: temperature in degrees Celsius and pressure in MPa. By default, the accuracy of the result is 3, which corresponds to 3 decimal places ±0.001 kg/m3. Near the critical point (P=22.064MPa and t=373.946°C) the function cannot accurately determine the value due to the fact that the graph of dependence ρ=f(p) becomes almost horizontal and the error increases to 3-4% and more. To increase the accuracy, it is recommended to enter a higher value of accuracy. The maximum allowable value of the parameter is limited to accuracy=13 !!! This is caused by the accuracy with which VBA operates with values of double type.
The function works in the range of densities from 113.6 to 762.4 kg/m3.

4. Auxiliary functions

• Saturation_temperature(p As Double) As Double, [°C]
• Saturation_pressure(t As Double) As Double, [MPa]

The functions describe a plot of the saturation line curve and return a second parameter corresponding to their temperature or pressure. The arguments and results of the functions are expressed: pressure in MPa, temperature in degrees Celsius.

• Temperature_boundary(p As Double) As Double, [°C]
• Pressure_boundary(t As Double) As Double, [MPa]

The functions describe the graph of the curve of the line limiting region 2 (vapor) from regions 1 and 3. The parameters of the functions are identical to the saturation line functions.

Note: at temperatures up to 350°C, the functions coincide with the saturation line functions.

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

The function returns the number of the region corresponding to the entered pressure and temperature. The region numbers correspond to the diagram given in IF97.

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

The function returns for all regions the value of the given characteristic energies of the region reg. For region 3 it is the Helmholtz free energy function f(ρ,t)/RT, and for the other regions the Gibbs energy g(p,t)/RT.
The function has 4 mandatory parameters:

1. Temperature in degrees Celsius
2. Pressure in MPa
3. The region has acceptable values - 1,2,3,4,5,21
4. Trigger has the following possible options:
   • a. The function itself is 0
   • b. The partial derivative of temperature - 1
   • c. Second partial derivative of temperature - 2
   • d. Mixed partial derivative of temperature and pressure (or density for region 3) - 3
   • e. The partial derivative of pressure (or density for the 3rd region) - 4
   • f. Second partial derivative of pressure (or density for the 3rd region) - 5

• Viscosity(t As Double, p As Double, Optional reg) As Double, [Pa∙s]

The function returns the value of dynamic viscosity expressed in Pa∙s calculated by the formula given in “Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance” September 2008. But since the parameter μ2 is taken equal to one, the function is not applicable in the critical region (P=22,064MPa and t=373,946°C). In other regions the function gives values with acceptable accuracy up to 1-2%. Otherwise, the definition area of the function completely coincides with the definition area of the density function for the IF97 document.

5. Functions of MI 2412-97
Density and enthalpy functions calculated by formulas P.1 and P.2 given in MI 2412-97 “Water heating systems. Equations of measurement of heat energy and quantity of heat carrier”. (https://meganorm.ru/Data2/1/4294845/4294845246.htm - Recommendation from the Russian State System for Ensuring Uniformity of Measurements Water heating systems equations for measuring heat energy and quantity of heat transfer fluid ).

• Density_MI(t As Double, p As Double) As Double, [kg/m3]
• Specific_enthalpy_MI(t As Double, p As Double) As Double, [J/kg]

VBA code

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 ' Triggers of the Gibbs energy function
Const R = 461.526, Default_accuracy = 3
Sub Auto_Open()
    MsgBox "The functions are written using :" & vbCrLf & "1. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam (2007 edition)." & vbCrLf & "2. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance." & vbCrLf & "3. MI 2412-97 Water heat supply systems equations of heat energy and heat carrier quantity measurements. (Recommendation from the Russian State System for Ensuring Uniformity of Measurements Water heating systems equations for measuring heat energy and quantity of heat transfer fluid )" 
    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 = Pressure_boundary(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 Density3(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 Saturation_pressure(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 Density3 = 1 / 0 Else Density3 = ro_mid
End Function
Function Helmholtz_energy(t As Double, ro As Double) As Double
    Helmholtz_energy = JF(t, ro, non, 3) * (t + 273.15) * R
End Function
Function Pressure3(t As Double, ro As Double) As Double
    Pressure3 = ro ^ 2 * R * (t + 273.15) * JF(t, ro, dp, 3) / 322
End Function
Function Specific_energy3(t As Double, ro As Double) As Double
    Specific_energy3 = R * 647.096 * JF(t, ro, dt, 3)
End Function
Function Specific_entropy3(t As Double, ro As Double) As Double
    Specific_entropy3 = R * (647.096 / (t + 273.15) * JF(t, ro, dt, 3) - JF(t, ro, non, 3))
End Function
Function Specific_enthalpy3(t As Double, ro As Double) As Double
    Specific_enthalpy3 = R * (647.096 * JF(t, ro, dt, 3) + (t + 273.15) * JF(t, ro, dp, 3) * ro / 322)
End Function
Function Heat_capacity_isobaric3(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)
    Heat_capacity_isobaric3 = (-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 Heat_capacity_isochoric3(t As Double, ro As Double) As Double
    Heat_capacity_isochoric3 = -(647.096 / (t + 273.15)) ^ 2 * JF(t, ro, dtt, 3) * R
End Function
Function Speed_Sound3(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)
    Speed_Sound3 = (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 Gibbs_energy(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
        Gibbs_energy = JF(t, p, non, Trigger) * (t + 273.15) * R
    Case Is = 3
        ro = Density3(t, p)
        Gibbs_energy = JF(t, ro, non, 3) * (t + 273.15) * R
    End Select
End Function
Function Specific_volume(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
        Specific_volume = (t + 273.15) * R * JF(t, p, dp, Trigger) / 16530000
    Case Is = 2, 4, 5, 21
        Specific_volume = (t + 273.15) * R * JF(t, p, dp, Trigger) / 1000000
    Case Is = 3
        Specific_volume = 1 / Density3(t, p)
    End Select
End Function
Function Density(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
        Density = 1 / ((t + 273.15) * R * JF(t, p, dp, Trigger) / 16530000)
    Case Is = 2, 4, 5, 21
        Density = 1 / ((t + 273.15) * R * JF(t, p, dp, Trigger) / 1000000)
    Case Is = 3
        Density = Density3(t, p)
    End Select
End Function
Function Specific_energy(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
        Specific_energy = R * (1386 * JF(t, p, dt, Trigger) - p / 16.53 * (t + 273.15) * JF(t, p, dp, Trigger))
    Case Is = 2, 4, 21
        Specific_energy = R * (540 * JF(t, p, dt, Trigger) - p * (t + 273.15) * JF(t, p, dp, Trigger))
    Case Is = 5
        Specific_energy = R * (1000 * JF(t, p, dt, Trigger) - p * (t + 273.15) * JF(t, p, dp, Trigger))
    Case Is = 3
        ro = Density3(t, p)
        Specific_energy = R * 647.096 * JF(t, ro, dt, Trigger)
    End Select
End Function
Function Specific_entropy(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
        Specific_entropy = R * (1386 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))
    Case Is = 2, 4, 21
        Specific_entropy = R * (540 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))
    Case Is = 5
        Specific_entropy = R * (1000 / (t + 273.15) * JF(t, p, dt, Trigger) - JF(t, p, non, Trigger))
    Case Is = 3
        ro = Density3(t, p)
        Specific_entropy = R * (647.096 / (t + 273.15) * JF(t, ro, dt, Trigger) - JF(t, ro, non, Trigger))
    End Select
End Function
Function Specific_enthalpy(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
        Specific_enthalpy = R * 1386 * JF(t, p, dt, Trigger)
    Case Is = 2, 4, 21
        Specific_enthalpy = R * 540 * JF(t, p, dt, Trigger)
    Case Is = 5
        Specific_enthalpy = R * 1000 * JF(t, p, dt, Trigger)
    Case Is = 3
        ro = Density3(t, p)
        Specific_enthalpy = R * (647.096 * JF(t, ro, dt, Trigger) + (t + 273.15) * JF(t, ro, dp, Trigger) * ro / 322)
    End Select
End Function
Function Heat_capacity_isobaric(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
        Heat_capacity_isobaric = -R * (1386 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)
    Case Is = 2, 4, 21
        Heat_capacity_isobaric = -R * (540 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)
    Case Is = 5
        Heat_capacity_isobaric = -R * (1000 / (t + 273.15)) ^ 2 * JF(t, p, dtt, Trigger)
    Case Is = 3
        ro = Density3(t, p)
        tf = 647.096 / (t + 273.15)
        rof = ro / 322
        fp = JF(t, ro, dp, Trigger)
        Heat_capacity_isobaric = (-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 Heat_capacity_isochoric(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)
        Heat_capacity_isochoric = 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)
        Heat_capacity_isochoric = 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)
        Heat_capacity_isochoric = 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 = Density3(t, p)
        Heat_capacity_isochoric = -(647.096 / (t + 273.15)) ^ 2 * JF(t, ro, dtt, 3) * R
    End Select
End Function
Function Speed_sound(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)
        Speed_sound = 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)
        Speed_sound = 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)
        Speed_sound = 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 = Density3(t, p)
        tf = 647.096 / (t + 273.15)
        rof = ro / 322
        fp = JF(t, ro, dp, Trigger)
        Speed_sound = (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 Saturation_temperature(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)
    Saturation_temperature = (n10 + D - ((n10 + D) ^ 2 - 4 * (n9 + n10 * D)) ^ 0.5) / 2 - 273.15
End Function
Function Saturation_pressure(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
    Saturation_pressure = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4
End Function
Function Temperature_boundary(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)
        Temperature_boundary = (n10 + D - ((n10 + D) ^ 2 - 4 * (n9 + n10 * D)) ^ 0.5) / 2 - 273.15
    Else
        Temperature_boundary = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5 - 273.15
    End If
End Function
Function Pressure_boundary(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
        Pressure_boundary = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4
    Else
        Pressure_boundary = 348.05185628969 - 1.1671859879975 * (t + 273.15) + 1.0192970039326E-03 * (t + 273.15) ^ 2
    End If
End Function
Function Viscosity(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 = Density(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
    Viscosity = 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 Density_MI(t As Double, p As Double) As Double
    tau = (t + 273.15) / 647.14
    Pi = p / 22.064
    Density_MI = 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 Specific_enthalpy_MI(t As Double, p As Double) As Double
    tau = (t + 273.15) / 647.14
    Pi = p / 22.064
    Specific_enthalpy_MI = (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

Sort by date Sort by score

Enable comment auto-refresher

Add your comment

petrofaq welcomes all comments. If you do not want to be anonymous, register or log in. It is free.