氮气物性计算-源代码-VBA-加载宏-可以在Excel单元格里直接调用：在Excel后台里插入模块，然后复制粘贴输入下面代码，就可以在前面Excel单元格里调用计算氩气物性的函数进行物性计算了

'氮气

'常数

Const R = 296.78 'J/kg.K，气体常数

Const M = 0.028016 'kg/mol，氮气的分子量

Const Pcr = 3.396 'MPa，临界压力

Const roucr = 304 'kg/m3，临界密度

Const tb = -195.8 '℃，沸点

Const tm = -210#  '℃，熔点

'Const cp0 = 5198 'J/kg.K

'Const cv0 = 3121 'J/kg.K

'Const h0 = 5557 'J/kg 温度315.56℃，压力0.1013MPa时氮气的焓

'Const s0 = 28016 'J/kg.K 温度315.56℃，压力0.1013MPa时氮气的熵

Const T0 = 273.15 'K

Const p0 = 100000 'Pa

'Const rou0 = 0.1762 'kg/m3

'1.比容计算，m3/kg

Function f(V, P, T) 'p-MPa,T-K,v-m3/kg

    Dim A0, a, B0, b, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        b = -0.000246645

        c = 1499.14

        'p = p * 1000000

       

    f = ((R * T * (1 - c / (V * T ^ 3))) / (V ^ 2)) * (V + B0 * (1 - b / V))

    f = f - (1 - a / V) * A0 / V ^ 2

    f = f - P

    'f = f

    On Error Resume Next

      

End Function

'比容，m3/kg

Function V_N2(P, T)

    Dim vi(10000)

    vi(1) = 0.01

    vi(2) = 4

    WuCha = 0.00001

   

    If f(vi(1), P, T) <> 0 And f(vi(2), P, T) Then

    vi(3) = vi(2) - f(vi(2), P, T) * (vi(2) - vi(1)) / (f(vi(2), P, T) - f(vi(1), P, T))

    'Debug.Print "vi(3)="; vi(3)

    'v = vi(3)

    End If

   

   

    If f(vi(2), P, T) <> 0 And f(vi(3), P, T) Then

    vi(4) = vi(3) - f(vi(3), P, T) * (vi(3) - vi(1)) / (f(vi(3), P, T) - f(vi(1), P, T))

    'Debug.Print "vi(4)="; vi(4)

    'v = vi(4)

    End If

    Dim i As Integer

    i = 4

   

    'Debug.Print "f(vi(1),p,T)="; f(vi(1), p, T)

    'Debug.Print "f(vi(2),p,T)="; f(vi(2), p, T)

    Do Until Abs(f(vi(i), P, T)) < 0.00001

   

        vi(i + 1) = vi(i) - f(vi(i), P, T) * (vi(i) - vi(1)) / (f(vi(i), P, T) - f(vi(1), P, T))

       

        i = i + 1

        'Debug.Print "vi("; i; ")="; vi(i)

        'Debug.Print "f("; i; ","; p; ","; T; ")="; f(vi(i), p, T)

        If i >= 10000 Then

        Exit Function

        End If

       

    Loop

   

    V_N2 = vi(i)

   

   

    On Error Resume Next

   

   

End Function

'密度，kg/m3

Function rou_N2(P, T)

    rou_N2 = 1 / V_N2(P, T)

   

End Function

'2.氮气的比热

Function cp0_N2(T)

    Dim b()

    ReDim Preserve b(0 To 4)

    b(0) = 3.771078

    b(1) = -0.001939065

    b(2) = 0.000004287664

    b(3) = -0.000000002810956

    b(4) = 6.260205E-13

    cp0_N2 = 0

    For i = 0 To 4

        cp0_N2 = cp0_N2 + b(i) * T ^ (i)

    Next

    cp0_N2 = R * cp0_N2

   

End Function

Function cv0_N2(T)

    cv0_N2 = cp0_N2(T) - R

End Function

Function cv_N2(P, T) 'J/kg.K

    Dim A0, a, B0, b, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        b = -0.000246645

        c = 1499.14

        Dim V

        V = V_N2(P, T)

       

        cv_N2 = cv0_N2(T) + (6 * R * c / ((T ^ 3) * V)) * (1 + B0 / V * (0.5 - b / (3 * V)))

End Function

Function cp_N2(P, T)

    Dim A0, a, B0, b, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        b = -0.000246645

        c = 1499.14

        Dim V

        V = V_N2(P, T)

    Dim z1, z2, z3, m1, m2, m3, m4, m5

        z1 = V + B0 * (1 - b / V)

        z2 = 1 + 2 * c / (T * T * T * V)

        z3 = R * (z1 * z2) ^ 2

        m1 = (2 * P * V ^ 3) / (R * T)

        m2 = -(c * B0 / T ^ 3) * (1 - 2 * b / V)

        m3 = -V ^ 2 - B0 * b

        m4 = A0 * a / (R * T)

        m5 = m1 + m2 + m3 + m4

    cp_N2 = cv_N2(P, T) + z3 / m5

End Function

'3 焓值计算

Function H_N2(P, T) 'k/kg

      Dim A0, a, B0, bb, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        bb = -0.000246645

        c = 1499.14

        Dim V

        V = V_N2(P, T)

        'Debug.Print "v="; v

        Dim b()

        ReDim Preserve b(0 To 4)

        b(0) = 3.771078

        b(1) = -0.001939065

        b(2) = 0.000004287664

        b(3) = -0.000000002810956

        b(4) = 6.260205E-13

       

    Dim f1, f2, f3, f4, f5, f6

        f1 = 0

        For N = 0 To 4

            f1 = f1 + (b(N) * T ^ N) / (N + 1)

        Next

       

        f1 = R * T * (f1 - 1)

       

        f2 = (A0 / V) * (a / (2 * V) - 1)

        f3 = -3 * R * c / (T * T * c)

        f4 = 1 + (B0 / V) * (0.5 - bb / (3 * V))

        f5 = P * V - 8082.61

        f6 = f1 + f2 + f3 * f4 + f5

        H_N2 = f6

End Function

'熵计算

Function S_N2(P, T)

    's = cp0 * Log(T / T0) - R * Log(p / p0) - BP(T) * p + s0

   

     Dim A0, a, B0, bb, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        bb = -0.000246645

        c = 1499.14

        Dim V

        V = V_N2(P, T)

        'Debug.Print "v="; v

        Dim b()

        ReDim Preserve b(0 To 4)

        b(0) = 3.771078

        b(1) = -0.001939065

        b(2) = 0.000004287664

        b(3) = -0.000000002810956

        b(4) = 6.260205E-13

       

    Dim f1, f2, f3, f4, f5, f6, f7, f8

        f1 = 0

        For N = 1 To 4

            f1 = f1 + (b(N) * T ^ N) / N

        Next

       

        f1 = R * ((b(0) - 1) * Log(T) + f1)

       

        f2 = Log(Abs(V))

        f3 = 1 - bb / (2 * V)

        f4 = 2 / B0 + 1 / V - 2 * bb / (3 * V * V)

        f5 = c * f4 / T ^ 3

        f6 = -(f5 + f3) * B0 / V

        f7 = R * (f2 + f6)

        f8 = f1 + f7 + 2311.59

        S_N2 = f8

   

   

End Function

'粘度

'压力在100000-140*100000Pa范围内，粘度只是温度的函数'kg/m.s

Function yita_N2(P, T)

    Dim D(5)

    D(1) = 0.0000013974

    D(2) = 107

    D(3) = 0.00115

    D(4) = 1

    D(5) = 1.7

    Dim yita1

    yita1 = (D(1) * T ^ 0.5) / (1 + D(2) / T)

    Dim P1, T1

    P1 = 101325

    T1 = 273.15

   

    yita_N2 = yita1 * (1 + D(3) * ((P / P1 - 1) ^ D(4)) * (T1 / T) ^ D(5))

   

   

End Function

'导热系数'W/m.K

Function lamda_N2(P, T)

     Dim a()

     ReDim Preserve a(0 To 4)

        a(0) = 0.001310641

        a(1) = 0.000092366747

        a(2) = -0.000000034716902

        a(3) = -5.089499E-12

        a(4) = 1.0818806E-14

    Dim lamda0

    lamda0 = 0

    Dim N As Integer

    For N = 0 To 4

        lamda0 = lamda0 + a(N) * T ^ N

    Next

    Dim rr

    rr = rou_N2(P, T)

    lamda_N2 = lamda0 + 0.00001617573 * rr ^ 1.23

End Function

'普朗特数

Function Pr_N2(P, T)

    Dim cp, yita, lamda

    cp = cp_N2(P, T)

    yita = yita_N2(P, T)

    lamda = lamda_N2(P, T)

    Pr_N2 = cp * yita / lamda

   

End Function

'声速 m/s

Function Vs_N2(P, T)

    '等熵指数K计算

   

    Dim k, cp, cv, V

    Dim A0, a, B0, b, c

        A0 = 173.571

        a = 0.000934109

        B0 = 0.00177511

        b = -0.000246645

        c = 1499.14

        cp = cp_N2(P, T)

        cv = cv_N2(P, T)

        V = V_N2(P, T)

       

    Dim f1, f2, f3, f4, f5, f6, f7

        f1 = cp / cv

        f2 = 296.78 * T / (P * V)

        f3 = c * B0 / (T * T * T * V * V)

        f4 = 1 - 2 * b / V

        f5 = 1 + B0 * b / (V * V)

        f6 = A0 * a / (V * V * 296.78 * T)

        f7 = f1 * (2 - f2 * (f3 * f4 + f5 - f6))

        k = f7

   

    Vs_N2 = (k * P * V) ^ 0.5

End Function

'压缩因子

Function Z_N2(P, T)

    Z_N2 = P * V_N2(P, T) / (296.78 * T)

End Function

'反算，根据压力，熵反算焓

Function f_pTs(P, T, S)

    f_pTs = S_N2(P, T) - S

End Function

Function H_PS_N2(pp, ss) '根据压力和熵求焓

   Dim Ti(10000), f1, f2, f3

    Ti(1) = 100

    Ti(2) = 1300

    Do Until Abs(Ti(2) - Ti(1)) < 0.00001

            Ti(3) = (Ti(1) + Ti(2)) / 2

            f1 = f_pTs(pp, Ti(1), ss)

            f2 = f_pTs(pp, Ti(2), ss)

            f3 = f_pTs(pp, Ti(3), ss)

        If f1 * f3 < 0 Then

            Ti(2) = Ti(3)

        Else

            Ti(1) = Ti(3)

        End If

    Loop

        H_PS_N2 = H_N2(pp, Ti(3))

       

End Function

Function f_pTh(P, T, H)

    f_pTh = H_N2(P, T) - H

End Function

Function S_PH_N2(pp, hh) '根据压力和焓求熵

   Dim Ti(10000), f1, f2, f3

    Ti(1) = 100

    Ti(2) = 1300

    Do Until Abs(Ti(2) - Ti(1)) < 0.00001

            Ti(3) = (Ti(1) + Ti(2)) / 2

            f1 = f_pTh(pp, Ti(1), hh)

            f2 = f_pTh(pp, Ti(2), hh)

            f3 = f_pTh(pp, Ti(3), hh)

        If f1 * f3 < 0 Then

            Ti(2) = Ti(3)

        Else

            Ti(1) = Ti(3)

        End If

    Loop

        S_PH_N2 = S_N2(pp, Ti(3))

End Function

Function T_PH_N2(pp, hh)

    Dim Ti(10000), f1, f2, f3

    Ti(1) = 100

    Ti(2) = 1300

    Do Until Abs(Ti(2) - Ti(1)) < 0.00001

            Ti(3) = (Ti(1) + Ti(2)) / 2

            f1 = f_pTh(pp, Ti(1), hh)

            f2 = f_pTh(pp, Ti(2), hh)

            f3 = f_pTh(pp, Ti(3), hh)

        If f1 * f3 < 0 Then

            Ti(2) = Ti(3)

        Else

            Ti(1) = Ti(3)

        End If

    Loop

        T_PH_N2 = Ti(3)

End Function

Function T_PS_N2(pp, ss)

    Dim Ti(10000), f1, f2, f3

    Ti(1) = 99

    Ti(2) = 1300

    Do Until Abs(Ti(2) - Ti(1)) < 0.00001

            Ti(3) = (Ti(1) + Ti(2)) / 2

            f1 = f_pTs(pp, Ti(1), ss)

            f2 = f_pTs(pp, Ti(2), ss)

            f3 = f_pTs(pp, Ti(3), ss)

            'Debug.Print "f1="; f1

            'Debug.Print "f2="; f2

            'Debug.Print "f3="; f3

        If f1 * f3 < 0 Then

            Ti(2) = Ti(3)

        Else

            Ti(1) = Ti(3)

        End If

    Loop

        T_PS_N2 = Ti(3)

End Function

Function P_TH_N2(T, H)

    Dim pi1, pi2, pi3, f1, f2, f3

    pi1 = 99999

    pi2 = 1100000

    Do Until Abs(pi2 - pi1) < 0.00001

            pi3 = (pi1 + pi2) / 2

            f1 = f_pTh(pi1, T, H)

            f2 = f_pTh(pi2, T, H)

            f3 = f_pTh(pi3, T, H)

        If f1 * f3 < 0 Then

            pi2 = pi3

        Else

            pi1 = pi3

        End If

    Loop

        P_TH_N2 = pi3

End Function

Function S_TH_N2(T, H)

    S_TH_N2 = S_N2(P_TH_N2(T, H), T)

End Function

Function P_TS_N2(T, S)

    Dim pi1, pi2, pi3, f1, f2, f3

    pi1 = 99999

    pi2 = 1100000

    Do Until Abs(pi2 - pi1) < 0.00001

            pi3 = (pi1 + pi2) / 2

            f1 = f_pTs(pi1, T, S)

            f2 = f_pTs(pi2, T, S)

            f3 = f_pTs(pi3, T, S)

        If f1 * f3 < 0 Then

            pi2 = pi3

        Else

            pi1 = pi3

        End If

    Loop

        P_TS_N2 = pi3

End Function

Function H_TS_N2(T, S)

    H_TS_N2 = H_N2(P_TS_N2(T, S), T)

End Function

Function V_PS_N2(P, S) '反算，根据压力和熵计算比容 V代表比容m3/kg

    V_PS_N2 = V_N2(P, T_PS_N2(P, S))

End Function

Function V_PH_N2(P, H)

    V_PH_N2 = V_N2(P, T_PH_N2(P, H))

End Function

Function V_TS_N2(T, S)

    V_TS_N2 = V_N2(P_TS_N2(T, S), T)

End Function

Function V_TH_N2(T, H) '误差有点大

    V_TH_N2 = V_N2(P_TH_N2(T, H), T)

End Function

Function f_pTv(P, T, V)

 f_pTv = V_N2(P, T) - V

End Function

Function P_VT_N2(V, T)

    Dim pi1, pi2, pi3, f1, f2, f3

    pi1 = 99999

    pi2 = 1100000

    Do Until Abs(pi2 - pi1) < 0.00001

            pi3 = (pi1 + pi2) / 2

            f1 = f_pTv(pi1, T, V)

            f2 = f_pTv(pi2, T, V)

            f3 = f_pTv(pi3, T, V)

        If f1 * f3 < 0 Then

            pi2 = pi3

        Else

            pi1 = pi3

        End If

    Loop

        P_VT_N2 = pi3

End Function

Function H_VT_N2(V, T)

    H_VT_N2 = H_N2(P_VT_N2(V, T), T)

End Function

Function S_VT_N2(V, T)

    S_VT_N2 = S_N2(P_VT_N2(V, T), T)

End Function

Function T_VP_N2(V, P)

    Dim Ti(10000), f1, f2, f3

    Ti(1) = 99

    Ti(2) = 1300

    Do Until Abs(Ti(2) - Ti(1)) < 0.00001

            Ti(3) = (Ti(1) + Ti(2)) / 2

            f1 = f_pTv(P, Ti(1), V)

            f2 = f_pTv(P, Ti(2), V)

            f3 = f_pTv(P, Ti(3), V)

        If f1 * f3 < 0 Then

            Ti(2) = Ti(3)

        Else

            Ti(1) = Ti(3)

        End If

    Loop

        T_VP_N2 = Ti(3)

End Function

Function H_VP_N2(V, P)

    H_VP_N2 = H_N2(P, T_VP_N2(V, P))

End Function

Function S_VP_N2(V, P)

    S_VP_N2 = S_N2(P, T_VP_N2(V, P))

End Function

Private Sub CommandButton1_Click()

    Debug.Print "="; S_VP_N2(0.51198, 100000)

End Sub