氦气物性计算-VBA代码

'氦气
'常数
Public Const R = 2077.1 'J/kg.K，气体常数
Public Const M = 0.004003 'kg/mol，氦气的分子量
Public Const Pcr = 0.229 'MPa，临界压力
Public Const roucr = 69.3 'kg/m3，临界密度
Public Const tb = -268.94 '℃，沸点
Public Const tm = -267.96 '℃，熔点
Public Const cp0 = 5198 'J/kg.K
Public Const cv0 = 3121 'J/kg.K
Public Const h0 = 5557 'J/kg 温度315.56℃，压力0.1013MPa时氦气的焓
Public Const s0 = 28016 'J/kg.K 温度315.56℃，压力0.1013MPa时氦气的熵
Public Const T0 = 273.15 'K
Public Const p0 = 100000 'Pa
Public Const rou0 = 0.1762 'kg/m3


'维里系数B(T)
Function b(T)
    Dim c(5)
    c(1) = 0.0009489433
    c(2) = 0.0009528079
    c(3) = 0.0342068
    c(4) = 0.00273947
    c(5) = 0.000940912
    b = c(1) + c(2) / (1 - c(3) * T) + c(4) / (1 + c(5) * T)
End Function
'比容，m3/kg
Function v(p, T)
    v = R * T / p + b(T)
End Function
'密度，kg/m3
Function rou(p, T)
    rou = p / (R * T + b(T) * p)
End Function
'维里系数的一阶导数
Function BP(T)
    Dim c(5)
    c(1) = 0.0009489433
    c(2) = 0.0009528079
    c(3) = 0.0342068
    c(4) = 0.00273947
    c(5) = 0.000940912
    BP = c(2) * c(3) / (1 - c(3) * T) ^ 2 - c(4) * c(5) / (1 + c(5) * T) ^ 2
End Function
'焓值计算
Function h(p, T)
    h = cp0 * T + h0 + (b(T) - BP(T) * T) * p
End Function
'熵计算
Function s(p, T)
    s = cp0 * Log(T / T0) - R * Log(p / p0) - BP(T) * p + s0
    
End Function
'粘度
'压力在100000-140*100000Pa范围内，粘度只是温度的函数
Function yita(T)
    yita = 0.00001855 * (T / T0) ^ 0.68
End Function
'导热系数
Function lamda(p, T)
    Dim a, n, lamda0
    a = 0.0001838
    n = 1.17
    lamda0 = 0.1448 * (T / T0) ^ 0.68
    lamda = lamda0 * (1 + 0.0001665 * ((p / p0) ^ 1.17) / (T / T0) ^ 1.85)
    
    
End Function
'普朗特数
Function Pr(p, T)
    Pr = 0.666 * (1 + 0.0001665 * ((p / p0) ^ 1.17) / (T / T0) ^ 1.85) ^ (-1)
End Function
'压缩因子
Function Z(p, T)
    Z = p * v(p, T) / (R * T)
End Function
'声速 m/s
Function Vs(p, T)
    Vs = 1.2903 * Z(p, T) * (R * T) ^ 0.5
End Function


Sub bb()
    Debug.Print Vs(0.1, 173.15)
    
    
    
End Sub