二氧化碳物性计算源代码-VBA版-span算法

Option Explicit
    '计算Φr的相关数组
    Public n(1 To 42) As Double
    Public d(1 To 39) As Double
    Public T(1 To 39) As Double
    Public c(8 To 34) As Double
    Public aerf(35 To 39) As Double
    Public beita(35 To 42) As Double
    Public gama(35 To 39) As Double
    Public ebsl(35 To 39) As Double
    Public a(40 To 42) As Double
    Public b(40 To 42) As Double
    Public aa(40 To 42) As Double
    Public bb(40 To 42) As Double
    Public cc(40 To 42) As Double
    Public dd(40 To 42) As Double
    Public i As Integer
    '常数的值
    Const M = 44.0098: Const Rm = 8.31451
    Const R = 0.1889241
    Const Tc = 304.1282: Const pc = 7.3773
    Const rouc = 467.6
    Const Tt = 216.592: Const pt = 0.51795
    Const T0 = 298.15: Const p0 = 0.101325
    'Const ho0 = 0: so0 = 0
    
    
    
'Φr初始化系数
Public Sub init()

    n(1) = 0.38856823203161: n(2) = 2.938547594274: n(3) = -5.5867188534934
    n(4) = -0.76753199592477: n(5) = 0.31729005580416: n(6) = 0.54803315897767
    n(7) = 0.12279411220335: n(8) = 2.165896154322: n(9) = 1.5841735109724
    n(10) = -0.23132705405503: n(11) = 0.058116916431436: n(12) = -0.55369137205382
    n(13) = 0.48946615909422: n(14) = -0.024275739843501: n(15) = 0.062494790501678
    n(16) = -0.12175860225246: n(17) = -0.37055685270086: n(18) = -0.016775879700426
    n(19) = -0.11960736637987: n(20) = -0.045619362508778: n(21) = 0.035612789270346
    n(22) = -7.4427727132052E-03: n(23) = -1.7395704902432E-03: n(24) = -0.021810121289527
    n(25) = 0.024332166559236: n(26) = -0.037440133423463: n(27) = 0.14338715756878
    n(28) = -0.13491969083286: n(29) = -0.02315122505348: n(30) = 0.012363125492901
    n(31) = 0.002105832197294: n(32) = -3.3958519026368E-04: n(33) = 5.5993651771592E-03
    n(34) = -3.0335118055646E-04: n(35) = -213.6548868832: n(36) = 26641.569149272
    n(37) = -24027.212204557: n(38) = -283.41603423999: n(39) = 212.47284400179
    n(40) = -0.66642276540751: n(41) = 0.72608632349897: n(42) = 0.055068668612842
    

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

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

    c(8) = 1: c(9) = 1: c(10) = 1
    c(11) = 1: c(12) = 1: c(13) = 1: c(14) = 1: c(15) = 1
    c(16) = 1: c(17) = 2: c(18) = 2: c(19) = 2: c(20) = 2
    c(21) = 2: c(22) = 2: c(23) = 2: c(24) = 3: c(25) = 3
    c(26) = 3: c(27) = 4: c(28) = 4: c(29) = 4: c(30) = 4
    c(31) = 4: c(32) = 4: c(33) = 5: c(34) = 6
    

    aerf(35) = 25: aerf(36) = 25: aerf(37) = 25
    aerf(38) = 15: aerf(39) = 20
    

    beita(35) = 325: beita(36) = 300: beita(37) = 300: beita(38) = 275: beita(39) = 275
    beita(40) = 0.3: beita(41) = 0.3: beita(42) = 0.3
    

    gama(35) = 1.16: gama(36) = 1.19: gama(37) = 1.19
    gama(38) = 1.25: gama(39) = 1.22
    

    ebsl(35) = 1: ebsl(36) = 1: ebsl(37) = 1
    ebsl(38) = 1: ebsl(39) = 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 fr(ByVal derta As Double, ByVal tao As Double) As Double
  Call init
  
   Dim fr1, fr2 As Double
   fr1 = 0: fr2 = 0
   
   For i = 1 To 7
        fr1 = fr1 + n(i) * (derta ^ d(i)) * tao ^ T(i)
    Next
    For i = 8 To 34
        fr2 = fr2 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * (2.718282 ^ (-derta * c(i)))
    Next
    fr = fr1 + fr2
End Function

'Φrδ(δ,τ) Eq.(6.5) a Table 31 P34
Public Function frd(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    
    Dim frd1, frd2 As Double
    'frd1 = 0: frd2 = 0
    For i = 1 To 7
        frd1 = frd1 + n(i) * d(i) * (derta ^ (d(i) - 1)) * tao ^ T(i)
    Next
    
    For i = 8 To 34
        frd2 = frd2 + n(i) * (2.718282 ^ (-derta ^ c(i))) * ((derta ^ (d(i) - 1)) * (tao ^ T(i)) * (d(i) - c(i) * derta ^ c(i)))
    Next
    frd = frd1 + frd2
End Function
'Φrδδ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function frdd(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim frdd1, frdd2 As Double
    For i = 1 To 7
        frdd1 = frdd1 + n(i) * d(i) * (d(i) - 1) * (derta ^ (d(i) - 2)) * (tao ^ T(i))
    Next
    For i = 8 To 34
        frdd2 = frdd2 + 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
    frdd = frdd1 + frdd2
End Function
'Φrτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function frt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim frt1, frt2 As Double
    For i = 1 To 7
        frt1 = frt1 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1))
    Next
    For i = 8 To 34
        frt2 = frt2 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1)) * 2.718282 ^ (-derta ^ c(i))
    Next
    frt = frt1 + frt2
    
End Function
'Φrττ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function frtt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim frtt1, frtt2 As Double
    For i = 1 To 7
        frtt1 = frtt1 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2))
    Next
    For i = 8 To 34
        frtt2 = frtt2 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2)) * 2.718282 ^ (-derta ^ c(i))
    Next
    frtt = frtt1 + frtt2
    
End Function
'Φrδτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function frdt(ByVal derta As Double, ByVal tao As Double) As Double
    Call init
    Dim frdt1, frdt2 As Double
    For i = 1 To 7
        frdt1 = frdt1 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 1))
    Next
    For i = 8 To 34
        frdt2 = frdt2 + 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
    frdt = frdt1 + frdt2
    
End Function
'the expression for Φ0(δ,τ) Eq.(6.3)
Function f0(derta, tao) As Double

    Dim a1, a2, a3, a4, a5 As Double
    Dim a6, a7, a8 As Double

        a1 = 8.37304456
        a2 = -3.70454304
        a3 = 2.5
        a4 = 1.99427042
        a5 = 0.62105248
        a6 = 0.41195293
        a7 = 1.04028922
        a8 = 0.08327678
        
    Dim a() As Double
    ReDim a(1 To 8)
    a(1) = a1: a(2) = a2: a(3) = a3: a(4) = a4
    a(5) = a5: a(6) = a6: a(7) = a7: a(8) = a8
    
    Dim st() As Double
    ReDim st(4 To 8)
    st(4) = 3.15163
    st(5) = 6.1119
    st(6) = 6.77708
    st(7) = 11.32384
    st(8) = 27.08792
    
    
 
    
    Dim Tc, rouc, R As Double
        Tc = 304.1282
        rouc = 467.6
        R = 0.1889241
        
    Dim i As Integer
    Dim z As Double
    
    
        z = 0
        z = z + Log(derta)
        z = z + a(1)
        z = z + a(2) * tao + a(3) * Log(tao)
    
    For i = 4 To 8
       
        z = z + a(i) * Log(1 - Exp(-tao * st(i)))
        
        
    Next
    
  
    f0 = z

End Function

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

End Function

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

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

End Function

'Φ0τ Eq.(6.3) Table 28 P34
Function f0t(tao) As Double
   
    Dim a1, a2, a3, a4, a5 As Double
    Dim a6, a7, a8 As Double

        a1 = 8.37304456
        a2 = -3.70454304
        a3 = 2.5
        a4 = 1.99427042
        a5 = 0.62105248
        a6 = 0.41195293
        a7 = 1.04028922
        a8 = 0.08327678
        
    Dim a() As Double
    ReDim a(1 To 8)
    a(1) = a1: a(2) = a2: a(3) = a3: a(4) = a4
    a(5) = a5: a(6) = a6: a(7) = a7: a(8) = a8
    
    Dim st() As Double
    ReDim st(4 To 8)
    st(4) = 3.15163
    st(5) = 6.1119
    st(6) = 6.77708
    st(7) = 11.32384
    st(8) = 27.08792

    Dim i As Integer
    Dim z As Double
    
    z = 0
    z = z + a(2)
    z = z + a(3) / tao
    For i = 4 To 8
        z = z + a(i) * st(i) * (1 / (1 - 2.718282 ^ (-st(i) * tao)) - 1)
    Next
    
    f0t = z
    
End Function
'Φ0ττ(τ) Eq.(6.3) Table 28 P34
Function f0tt(tao) As Double
   
    Dim a1, a2, a3, a4, a5 As Double
    Dim a6, a7, a8 As Double

        a1 = 8.37304456
        a2 = -3.70454304
        a3 = 2.5
        a4 = 1.99427042
        a5 = 0.62105248
        a6 = 0.41195293
        a7 = 1.04028922
        a8 = 0.08327678
        
    Dim a() As Double
    ReDim a(1 To 8)
    a(1) = a1: a(2) = a2: a(3) = a3: a(4) = a4
    a(5) = a5: a(6) = a6: a(7) = a7: a(8) = a8
    
    Dim st() As Double
    ReDim st(4 To 8)
    st(4) = 3.15163
    st(5) = 6.1119
    st(6) = 6.77708
    st(7) = 11.32384
    st(8) = 27.08792

    Dim i As Integer
    Dim z As Double
    
    z = 0
  
    z = z - a(3) / (tao * tao)
    
    For i = 4 To 8
        z = z - a(i) * st(i) * st(i) * (2.718282 ^ (-st(i) * tao)) * (1 - 2.718282 ^ (-st(i) * tao)) ^ (-2)
    Next
    
    f0tt = z
    
End Function




'Pressure根据密度，温度计算压力
Public Function rouT_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 frd(derta, tao)
    Dim z As Double
    z = frd(derta, tao)
    rouT_P = rou * R * T * (1 + derta * z) / 1000


End Function
'Entropy 根据密度、温度计算熵
Public Function rouT_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 * (f0t(tao) + frt(derta, tao)) - f0(derta, tao) - fr(derta, tao)
    
    rouT_S = z * R


End Function
'Internal energy 根据密度、温度计算内能
Public Function rouT_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 * (f0t(tao) + frt(derta, tao))
    
    rouT_U = z * R * T


End Function
'Enthalpy 根据密度、温度计算焓
Public Function rouT_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 * (f0t(tao) + frt(derta, tao)) + derta * frd(derta, tao)
    
    rouT_H = z * R * T


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

Sub test()

Debug.Print "rouT_U="; rouT_U(5.4054, 203.314)

End Sub