氮气物性计算源代码-python版：文章正文全为python代码，复制粘贴即可，最好在pycharm里应用

本文链接：https://blog.csdn.net/weixin_43441028/article/details/136577781

import math

# 常数

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

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

Pcr = 3.396 # MPa，临界压力

roucr = 304 # kg/m3，临界密度

tb = -195.8 # ℃，沸点

tm = -210 # ℃，熔点

T0 = 273.15 # K

p0 = 100000 # Pa

# 比容计算，m3/kg

def f(V, P, T):

A0 = 173.571

a = 0.000934109

B0 = 0.00177511

b = -0.000246645

c = 1499.14

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

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

f -= P

return f

# 比容，m3/kg

def V_N2(P, T):

vi = [0.01, 4]

WuCha = 0.00001

if f(vi[0], P, T) != 0 and f(vi[1], P, T) != 0:

vi[2] = vi[1] - f(vi[1], P, T) * (vi[1] - vi[0]) / (f(vi[1], P, T) - f(vi[0], P, T))

if f(vi[1], P, T) != 0 and f(vi[2], P, T) != 0:

vi[3] = vi[2] - f(vi[2], P, T) * (vi[2] - vi[0]) / (f(vi[2], P, T) - f(vi[0], P, T))

i = 3

while abs(f(vi[i], P, T)) >= 0.00001:

vi.append(vi[i] - f(vi[i], P, T) * (vi[i] - vi[0]) / (f(vi[i], P, T) - f(vi[0], P, T)))

i += 1

if i >= 10000:

break

return vi[i]

# 密度，kg/m3

def rou_N2(P, T):

return 1 / V_N2(P, T)

# 氮气的比热

def cp0_N2(T):

b = [3.771078, -0.001939065, 0.000004287664, -0.000000002810956, 6.260205E-13]

cp0 = sum([b[i] * T ** i for i in range(5)])

cp0 *= R

return cp0

def cv0_N2(T):

return cp0_N2(T) - R

def cv_N2(P, T):

A0 = 173.571

a = 0.000934109

B0 = 0.00177511

b = -0.000246645

c = 1499.14

V = V_N2(P, T)

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

return cv_N2

def cp_N2(P, T):

A0 = 173.571

a = 0.000934109

B0 = 0.00177511

b = -0.000246645

c = 1499.14

V = V_N2(P, T)

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

return cp_N2

# 焓值计算

def H_N2(P, T):

A0 = 173.571

a = 0.000934109

B0 = 0.00177511

bb = -0.000246645

c = 1499.14

V = V_N2(P, T)

b = [3.771078, -0.001939065, 0.000004287664, -0.000000002810956, 6.260205E-13]

f1 = 0

for N in range(5):

f1 += (b[N] * T ** N) / (N + 1)

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

return f6

# 熵计算

def S_N2(P, T):

A0 = 173.571

a = 0.000934109

B0 = 0.00177511

bb = -0.000246645

c = 1499.14

V = V_N2(P, T)

b = [3.771078, -0.001939065, 0.000004287664, -0.000000002810956, 6.260205E-13]

f1 = 0

for N in range(1, 5):

f1 += (b[N] * T ** N) / N

f1 = R * ((b[0] - 1) * math.log(T) + f1)

f2 = math.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

return f8

# 粘度

def yita_N2(P, T):

D = [0.0000013974, 107, 0.00115, 1, 1.7]

yita1 = (D[0] * T ** 0.5) / (1 + D[1] / T)

P1 = 101325

T1 = 273.15

yita_N2 = yita1 * (1 + D[2] * ((P / P1 - 1) ** D[3]) * (T1 / T) ** D[4])

return yita_N2

# 导热系数

def lamda_N2(P, T):

a = [0.001310641, 0.000092366747, -0.000000034716902, -5.089499E-12, 1.0818806E-14]

lamda0 = 0

for N in range(5):

lamda0 += a[N] * T ** N

rr = rou_N2(P, T)

lamda_N2 = lamda0 + 0.00001617573 * rr ** 1.23

return lamda_N2

# 普朗特数

def Pr_N2(P, T):

cp = cp_N2(P, T)

yita = yita_N2(P, T)

lamda = lamda_N2(P, T)

Pr_N2 = cp * yita / lamda

return Pr_N2

# 声速 m/s

def Vs_N2(P, T):

k = cp_N2(P, T) / cv_N2(P, T)

V = V_N2(P, T)

f1 = k * P * V

Vs_N2 = f1 ** 0.5

return Vs_N2

# 压缩因子

def Z_N2(P, T):

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

return Z_N2

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

def f_pTs(P, T, S):

f_pTs = S_N2(P, T) - S

return f_pTs

def H_PS_N2(pp, ss):

Ti = [100, 1300]

while 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:

Ti[2] = Ti[3]

else:

Ti[1] = Ti[3]

H_PS_N2 = H_N2(pp, Ti[3])

return H_PS_N2

def f_pTh(P, T, H):

f_pTh = H_N2(P, T) - H

return f_pTh

def S_PH_N2(pp, hh):

Ti = [100, 1300]

while 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:

Ti[2] = Ti[3]

else:

Ti[1] = Ti[3]

S_PH_N2 = S_N2(pp, Ti[3])

return S_PH_N2

def T_PH_N2(pp, hh):

Ti = [100, 1300]

while 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:

Ti[2] = Ti[3]

else:

Ti[1] = Ti[3]

T_PH_N2 = Ti[3]

return T_PH_N2

def T_PS_N2(pp, ss):

Ti = [99, 1300]

while 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:

Ti[2] = Ti[3]

else:

Ti[1] = Ti[3]

T_PS_N2 = Ti[3]

return T_PS_N2

def P_TH_N2(T, H):

pi1 = 99999

pi2 = 1100000

while 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:

pi2 = pi3

else:

pi1 = pi3

P_TH_N2 = pi3

return P_TH_N2

def S_TH_N2(T, H):

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

return S_TH_N2

def P_TS_N2(T, S):

pi1 = 99999

pi2 = 1100000

while 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:

pi2 = pi3

else:

pi1 = pi3

P_TS_N2 = pi3

return P_TS_N2

def H_TS_N2(T, S):

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

return H_TS_N2

def V_PS_N2(P, S):

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

return V_PS_N2

def V_PH_N2(P, H):

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

return V_PH_N2

def V_TS_N2(T, S):

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

return V_TS_N2

def V_TH_N2(T, H):

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

return V_TH_N2

def f_pTv(P, T, V):

f_pTv = V_N2(P, T) - V

return f_pTv

def P_VT_N2(V, T):

pi1 = 99999

pi2 = 1100000

while 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:

pi2 = pi3

else:

pi1 = pi3

P_VT_N2 = pi3

return P_VT_N2

def H_VT_N2(V, T):

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

return H_VT_N2

def S_VT_N2(V, T):

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

return S_VT_N2

def T_VP_N2(V, P):

Ti = [99, 1300]

while 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:

Ti[2] = Ti[3]

else:

Ti[1] = Ti[3]

T_VP_N2 = Ti[3]

return T_VP_N2

def H_VP_N2(V, P):

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

return H_VP_N2

def S_VP_N2(V, P):

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

return S_VP_N2

# 调用示例

print("S_VP_N2:", S_VP_N2(0.51198, 100000))