氩气物性计算源代码-VBA版

Option Explicit
    '计算Φr的相关数组
    Dim n(1 To 41) As Double
    Dim d(1 To 41) As Double
    Dim T(1 To 41) As Double
    Dim c(13 To 37) As Double
    'Dim aerf(35 To 39) As Double
    Dim yita(38 To 41) As Double
    Dim beita(38 To 41) As Double
    Dim gama(38 To 41) As Double
    Dim ebsl(38 To 41) As Double
    'Dim a(40 To 42) As Double
    'Dim b(40 To 42) As Double
    'Dim aa(40 To 42) As Double
    'Dim bb(40 To 42) As Double
    'Dim cc(40 To 42) As Double
    'Dim dd(40 To 42) As Double
    Dim i As Integer
    '常数的值
    Const M = 39.948: Const Rm = 8.31451
    Const R = 0.2081333
    Const Tc = 150.687: Const pc = 4.863
    Const rouc = 535.6
    Const Tt = 83.8058: Const pt = 68.891 'kPa
    Const T0 = 298.15: Const p0 = 0.101325 ' MPa
    'Const ho0 = 0: so0 = 0
    
    
    
'Φr初始化系数
Public Sub init()

    n(1) = 0.088722304990011: n(2) = 0.70514805167298: n(3) = -1.682011565409
    n(4) = -0.14909014431486: n(5) = -0.1202480460094: n(6) = -0.12164978798599
    
    n(7) = 0.40035933626752: n(8) = -0.27136062699129: n(9) = 0.24211924579645
    n(10) = 0.005788958318557: n(11) = -0.041097335615341: n(12) = 0.024710761541614
    
    n(13) = -0.32181391750702: n(14) = 0.33230017695794: n(15) = 0.031019986287345
    n(16) = -0.0307770860024: n(17) = 0.093891137419581: n(18) = -0.090643210682031
    
    n(19) = -4.5778349276654E-04: n(20) = -8.2659729025197E-05: n(21) = 1.3013415603147E-04
    n(22) = -0.011397840001996: n(23) = -0.024455169960535: n(24) = -0.064324067175955
    
    n(25) = 0.058889471093674: n(26) = -6.4933552112965E-04: n(27) = -0.013889862158435
    n(28) = 0.4048983929691: n(29) = -0.38612519594749: n(30) = -0.18817142332233
    
    n(31) = 0.15977647596482: n(32) = 0.053985518513856: n(33) = -0.028953417958014
    n(34) = -0.013025413381384: n(35) = 2.8948696775778E-03: n(36) = -2.2647134304796E-03
    
    n(37) = 1.7616456196368E-03: n(38) = 5.8552454482774E-03: n(39) = -0.69251908270028
    n(40) = 1.5315490030516: n(41) = -2.7380447449783E-03
    

    d(1) = 1: d(2) = 1: d(3) = 1: d(4) = 1: d(5) = 1
    d(6) = 2: d(7) = 2: d(8) = 2: d(9) = 2: d(10) = 3
    d(11) = 3: d(12) = 4: d(13) = 1: d(14) = 1: d(15) = 3
    d(16) = 4: d(17) = 4: d(18) = 5: d(19) = 7: d(20) = 10
    d(21) = 10: d(22) = 2: d(23) = 2: d(24) = 4: d(25) = 4
    d(26) = 8: d(27) = 3: d(28) = 5: d(29) = 5: d(30) = 6
    d(31) = 6: d(32) = 7: d(33) = 7: d(34) = 8: d(35) = 9
    d(36) = 5: d(37) = 6: d(38) = 2: d(39) = 1: d(40) = 2
    d(41) = 3
    

    T(1) = 0: T(2) = 0.25: T(3) = 1: T(4) = 2.75: T(5) = 4
    T(6) = 0: T(7) = 0.25: T(8) = 0.75: T(9) = 2.75: T(10) = 0
    
    T(11) = 2: T(12) = 0.75: T(13) = 3: T(14) = 3.5: T(15) = 1
    T(16) = 2: T(17) = 4: T(18) = 3: T(19) = 0: T(20) = 0.5
    
    T(21) = 1: T(22) = 1: T(23) = 7: T(24) = 5: T(25) = 6
    T(26) = 6: T(27) = 10: T(28) = 13: T(29) = 14: T(30) = 11
    T(31) = 14: T(32) = 8: T(33) = 14: T(34) = 6: T(35) = 7
    T(36) = 24: T(37) = 22: T(38) = 3: T(39) = 3: T(40) = 0
    T(41) = 0
    

    
    c(13) = 1: c(14) = 1: c(15) = 1
    c(16) = 1: c(17) = 1: c(18) = 1: c(19) = 1: c(20) = 1
    c(21) = 1: c(22) = 2: c(23) = 2: c(24) = 3: c(25) = 2
    c(26) = 2: c(27) = 3: c(28) = 3: c(29) = 3: c(30) = 3
    c(31) = 3: c(32) = 3: c(33) = 3: c(34) = 3: c(35) = 3
    c(36) = 4: c(37) = 4
    
    yita(38) = 20: yita(39) = 20: yita(40) = 20: yita(41) = 20
    

    'aerf(35) = 25: aerf(36) = 25: aerf(37) = 25
    'aerf(38) = 15: aerf(39) = 20
      
    beita(38) = 250: beita(39) = 375: beita(40) = 300: beita(41) = 225
    
    gama(38) = 1.11: gama(39) = 1.14: gama(40) = 1.17: gama(41) = 1.11
    

    
    ebsl(38) = 1: ebsl(39) = 1: ebsl(40) = 1: ebsl(41) = 1
    

    'a(40) = 3.5: a(41) = 3.5: a(41) = 3.5
    

    'b(40) = 0.875: b(41) = 0.925: b(42) = 0.875
    

    'aa(40) = 0.7: aa(41) = 0.7: aa(42) = 0.7
    

    'bb(40) = 0.3: bb(41) = 0.3: bb(42) = 1
    

    'cc(40) = 10: cc(41) = 10: cc(42) = 12.5
    
 
    'dd(40) = 275: dd(41) = 275: dd(42) = 275
End Sub

'Φr(δ,τ) Eq.(6.5) a Table 31 P34
Public Function afr(ByVal derta As Double, ByVal tao As Double) As Double
  Call init
  
   Dim afr1, afr2, afr3 As Double
   afr1 = 0: afr2 = 0: afr3 = 0
   
   For i = 1 To 12
        afr1 = afr1 + n(i) * (derta ^ d(i)) * tao ^ T(i)
    Next
    For i = 13 To 37
        afr2 = afr2 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-derta ^ c(i))
    Next
    For i = 38 To 41
        afr3 = afr3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)
    Next
    
    afr = afr1 + afr2 + afr3
End Function

'Φrδ(δ,τ) Eq.(6.5) a Table 31 P34
Public Function afrd(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    
    Dim afrd1, afrd2, afrd3 As Double
    afrd1 = 0: afrd2 = 0: afrd3 = 0
    
    For i = 1 To 12
        afrd1 = afrd1 + n(i) * d(i) * (derta ^ (d(i) - 1)) * tao ^ T(i)
    Next
    
    For i = 13 To 37
        afrd2 = afrd2 + n(i) * (2.718282 ^ (-derta ^ c(i))) * ((derta ^ (d(i) - 1)) * (tao ^ T(i)) * (d(i) - c(i) * derta ^ c(i)))
    Next
    For i = 38 To 41
        afrd3 = afrd3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * (2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (d(i) / derta - 2 * yita(i) * (derta - ebsl(i)))
    Next
    
    afrd = afrd1 + afrd2 + afrd3
    
End Function
'Φrδδ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function afrdd(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim afrdd1, afrdd2, afrdd3 As Double
    afrdd1 = 0: afrdd2 = 0: afrdd3 = 0
    For i = 1 To 12
        afrdd1 = afrdd1 + n(i) * d(i) * (d(i) - 1) * (derta ^ (d(i) - 2)) * (tao ^ T(i))
    Next
    For i = 13 To 37
        afrdd2 = afrdd2 + n(i) * (2.718282 ^ (-derta ^ c(i))) * ((derta ^ (d(i) - 2)) * (tao ^ T(i)) * ((d(i) - c(i) * derta ^ c(i)) * (d(i) - 1 - c(i) * derta ^ c(i)) - c(i) * c(i) * derta ^ c(i)))
    Next
    'Dim i As Integer
    For i = 38 To 41
        afrdd3 = afrdd3 + (n(i) * (tao ^ T(i)) * 2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (-2 * yita(i) * derta ^ d(i) + 4 * yita(i) * yita(i) * (derta ^ d(i)) * (derta - ebsl(i)) * (derta - ebsl(i)) ^ 2 - 4 * d(i) * yita(i) * (derta ^ (d(i) - 1)) * (derta - ebsl(i)) + d(i) * (d(i) - 1) * derta ^ (d(i) - 2))
    Next
    
    afrdd = afrdd1 + afrdd2 + afrdd3
End Function
'Φrτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function afrt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim afrt1, afrt2, afrt3 As Double
    afrt1 = 0: afrt2 = 0: afrt3 = 0
    For i = 1 To 12
        afrt1 = afrt1 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1))
    Next
    For i = 13 To 37
        afrt2 = afrt2 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1)) * 2.718282 ^ (-derta ^ c(i))
    Next
    For i = 38 To 41
        afrt3 = afrt3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ ((-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (T(i) / tao - 2 * beita(i) * (tao - gama(i)))
    
    afrt = afrt1 + afrt2 + afrt3
    Next
    
End Function
'Φrττ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function afrtt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim afrtt1, afrtt2, afrtt3 As Double
    afrtt1 = 0: afrtt2 = 0: afrtt3 = 0
    For i = 1 To 12
        afrtt1 = afrtt1 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2))
    Next
    For i = 13 To 37
        afrtt2 = afrtt2 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2)) * 2.718282 ^ (-derta ^ c(i))
    Next
    For i = 38 To 41
        afrtt3 = afrtt3 + (n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * ((T(i) / tao - 2 * beita(i) * (tao - gama(i))) ^ 2 - T(i) / (tao * tao) - 2 * beita(i))
    Next
    afrtt = afrtt1 + afrtt2 + afrtt3
    
End Function
'Φrδτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function afrdt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim afrdt1, afrdt2, afrdt3 As Double
    afrdt1 = 0: afrdt2 = 0: afrdt3 = 0
    For i = 1 To 12
        afrdt1 = afrdt1 + n(i) * d(i) * T(i) * (derta ^ (T(i) - 1)) * (tao ^ (T(i) - 1))
    Next
    For i = 13 To 37
        afrdt2 = afrdt2 + n(i) * (2.718282 ^ (-derta ^ c(i))) * (derta ^ (d(i) - 1)) * T(i) * (tao ^ (T(i) - 1)) * (d(i) - c(i) * derta ^ c(i))
    Next
    For i = 38 To 41
        afrdt3 = afrdt3 + (n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (d(i) / derta - 2 * yita(i) * (derta - ebsl(i))) * (T(i) / tao - 2 * beita(i) * (tao - gama(i)))
    Next
    
    afrdt = afrdt1 + afrdt2 + afrdt3
    
End Function
'Φrδττ(δ,τ) Eq.(6.5) a Table 31 P32
Public Function afrdtt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim afrdtt1, afrdtt2, afrdtt3 As Double
    afrdtt1 = 0: afrdtt2 = 0: afrdtt3 = 0
    For i = 1 To 12
        afrdtt1 = afrdt1 + n(i) * d(i) * T(i) * (T(i) - 1) * (derta ^ (T(i) - 1)) * (tao ^ (T(i) - 2))
    Next
    For i = 13 To 37
        afrdtt2 = afrdt2 + n(i) * (2.718282 ^ (-derta ^ c(i))) * T(i) * (T(i) - 1) * (tao ^ (T(i) - 2)) * (derta ^ (d(i) - 1)) * (d(i) - c(i) * derta ^ c(i))
    Next
    For i = 38 To 41
        afrdtt3 = afrdt3 + (n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-yita(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * ((T(i) / tao - 2 * beita(i) * (tao - gama(i))) ^ 2 - T(i) / (tao * tao) - 2 * beita(i)) * (d(i) / derta - 2 * yita(i) * (derta - ebsl(i)))
            
    afrdtt = afrdtt1 + afrdtt2 + afrdtt3
        
End Function

'the expression for Φ0(δ,τ) Eq.(6.3)
Function af0(derta, tao) As Double

    Dim a1, a2 As Double
        a1 = 8.31666243
        a2 = -4.94651164
        
        
   af0 = Log(derta) + a1 + a2 * tao + 1.5 * Log(tao)
   

End Function

'Φ0δ Eq.(6.3) Table 28 P34
Function af0d(derta) As Double
    Dim derta As Double
        af0d = 1 / derta

End Function

'Φ0δδ Eq.(6.3) Table 28 P34
Function af0dd(derta) As Double
    Dim derta As Double
        af0dd = -1 / (derta * derta)

End Function
'Φ0δτ(δ,τ) Eq.(6.3) Table 28 P34
Function af0dt(derta, tao) As Double
 
        af0dt = 0

End Function

'Φ0τ Eq.(6.3) Table 28 P34
Function af0t(tao) As Double
   
    
    af0t = -4.94651164 + 1.5 / tao
    
End Function
'Φ0ττ(τ) Eq.(6.3) Table 28 P34
Function af0tt(tao) As Double
   
    
    af0tt = -1.5 / (tao * tao)
    
End Function

'Pressure根据密度，温度计算压力
Public Function arouT_P(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    Call afrd(derta, tao)
    Dim z As Double
    z = afrd(derta, tao)
    arouT_P = rou * R * T * (1 + derta * z) / 1000

End Function
'Entropy 根据密度、温度计算熵
Public Function arouT_S(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    
    Dim z As Double
    z = tao * (af0t(tao) + afrt(derta, tao)) - af0(derta, tao) - afr(derta, tao)
    
    arouT_S = z * R

End Function
'Internal energy 根据密度、温度计算内能
Public Function arouT_U(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    
    Dim z As Double
    z = tao * (af0t(tao) + afrt(derta, tao))
    
    arouT_U = z * R * T

End Function
'Internal energy 根据密度、温度计算定容比热
Public Function arouT_Cv(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    
    Dim z As Double
    z = -tao * tao * (af0tt(tao) + afrtt(derta, tao))
    
    arouT_Cv = z * R

End Function

'Enthalpy 根据密度、温度计算焓
Public Function arouT_H(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    
    Dim z As Double
    z = 1 + tao * (af0t(tao) + afrt(derta, tao)) + derta * afrd(derta, tao)
    
    arouT_H = z * R * T

End Function
'Isobaric heat capacity 根据密度、温度计算定压热容
Public Function arouT_Cp(ByVal rou As Double, ByVal T As Double) As Double
    Dim derta As Double
    Dim tao As Double
    derta = rou / rouc: tao = Tc / T
    
    Dim z As Double
    z = -tao * tao * (af0tt(tao) + afrtt(derta, tao)) + ((1 + derta * afrd(derta, tao) - derta * tao * afrdt(derta, tao)) ^ 2) / (1 + 2 * derta * afrd(derta, tao) + derta * derta * afrdd(derta, tao))
    
    arouT_Cp = z * R

End Function

'反算，根据压力和焓计算熵
'Public Function erPH_S(ByVal P As Double, ByVal H As Double) As Double
'    Dim derta As Double
'    Dim tao As Double
 '   derta = rou / rouc: tao = Tc / T
    
'    Dim z As Double
'    z = 1 + tao * (f0t(tao) + frt(derta, tao)) + derta * frd(derta, tao)
    
 '   rouT_H = z * R * T

'End Function

'melting-pressure
Public Function aT_Pm(T) As Double

    Dim a1, a2, z As Double
        a1 = -7476.2665
        a2 = 9959.0613
        
        
   z = 1 + a1 * ((T / Tt) ^ 1.05 - 1) + a2 * ((T / Tt) ^ 1.275 - 1)
   aT_Pm = z * pt

End Function
' sublimation-pressure
Public Function aT_Psub(T) As Double

    Dim a1, a2, z As Double
        a1 = -11.391604
        a2 = -0.39513431
        
        
   z = (Tt / T) * (a1 * (1 - T / Tt) + a2 * (1 - T / Tt) ^ 2.7)
   aT_Psub = Exp(z) * pt

End Function

' The simple vapor-pressure
Public Function aT_Ps(T) As Double

    Dim a1, a2, a3, a4, z, S As Double
        a1 = -5.9409785
        a2 = 1.3553888
        a3 = -0.46497607
        a4 = -1.5399043
        S = 1 - T / Tc
        
   z = (Tc / T) * (a1 * S + a2 * S ^ 1.5 + a3 * S * S + a4 * S ^ 4.5)
   
   aT_Ps = Exp(z) * pc

End Function

' the saturated liquid
Public Function aT_roup(T) As Double

    Dim a1, a2, a3, a4, z, S As Double
        a1 = 1.5004262
        a2 = -0.3138129
        a3 = 0.086461622
        a4 = -0.041477525
        S = 1 - T / Tc
        
   z = a1 * S ^ 0.334 + a2 * S ^ (2 / 3) + a3 * S ^ (7 / 3) + a4 * S ^ 4
   
   aT_roup = Exp(z) * rouc

End Function

' the saturated vapor density
Public Function aT_roupp(T) As Double

    Dim a1, a2, a3, a4, z, S As Double
        a1 = -1.70695656
        a2 = -4.02739448
        a3 = 1.55177558
        a4 = -2.30683228
        S = 1 - T / Tc
        
    z = (Tc / T) * (a1 * S ^ 0.345 + a2 * S ^ (5 / 6) + a3 * S ^ 1 + a4 * S ^ (13 / 3))
   
    aT_roupp = Exp(z) * rouc

End Function

Private Sub CommandButton1_Click()
    Dim T, rou As Double
    T = CDbl(TextBox1.Value) + 273.15
    rou = aT_roup(T)
    TextBox2.Value = rou
    

End Sub

Private Sub CommandButton2_Click()
   Dim T, rou As Double
    T = CDbl(TextBox3.Value) + 273.15
    rou = aT_roupp(T)
    TextBox4.Value = rou
End Sub

Private Sub CommandButton3_Click()
    
    Dim T, rou As Double
    T = CDbl(TextBox5.Value) + 273.15
    rou = CDbl(TextBox6.Value)
    Dim p, S, U, Cv, H, Cp As Double
    
    p = arouT_P(rou, T)
    S = arouT_S(rou, T)
    U = arouT_U(rou, T)
    Cv = arouT_Cv(rou, T)
    H = arouT_H(rou, T)
    Cp = arouT_Cp(rou, T)
    
    TextBox7.Value = p
    TextBox8.Value = S
    TextBox9.Value = U
    TextBox10.Value = Cv
    TextBox11.Value = H
    TextBox12.Value = Cp
    

End Sub