这段VBA代码实现了超临界二氧化碳的物性计算，包括压力、熵、内能、定容比容和焓等参数。通过定义常数和系数，然后使用多个函数分别计算不同物理量，如`CO2_rouT_P`计算压力，`CO2_rouT_S`计算熵等。此外，还提供了饱和液体和饱和蒸汽的相关计算。
Option Explicit
'计算Φr的相关数组
Dim n(1 To 42) As Double
Dim d(1 To 39) As Double
Dim T(1 To 39) As Double
Dim c(8 To 34) As Double
Dim aerf(35 To 39) As Double
Dim beita(35 To 42) As Double
Dim gama(35 To 39) As Double
Dim ebsl(35 To 39) 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 = 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 'MPa
Const T0 = 298.15: Const p0 = 0.101325 ' MPa
'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
'1.1 the expression for Φ0(δ,τ) Eq.(6.3)
Public Function CO2_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 z As Double
z = Log(derta) + a(1) + a(2) * tao + a(3) * Log(tao)
For i = 4 To 8
z = z + a(i) * Log(1 - Exp(-tao * st(i)))
Next
CO2_f0 = z
End Function
'1.2 Φ0δ Eq.(6.3) Table 28 P34
Function CO2_f0d(derta) As Double
CO2_f0d = 1 / derta
End Function
'1.3 Φ0δδ Eq.(6.3) Table 28 P34
Function CO2_f0dd(derta) As Double
CO2_f0dd = -1 / (derta * derta)
End Function
'1.4 Φ0δτ(δ,τ) Eq.(6.3) Table 28 P34
Function CO2_f0dt(derta, tao) As Double
CO2_f0dt = 0
End Function
'1.5 Φ0τ Eq.(6.3) Table 28 P34
Function CO2_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 z As Double
z = a(2) + a(3) / tao
For i = 4 To 8
z = z + a(i) * st(i) * (1 / (1 - 2.718282 ^ (-st(i) * tao)) - 1)
Next
CO2_f0t = z
End Function
'1.6 Φ0ττ(τ) Eq.(6.3) Table 28 P34
Function CO2_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 z As Double
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
CO2_f0tt = z
End Function
'2.1Φr(δ,τ) Eq.(6.5) a Table 31 P34
Public Function CO2_fr(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 7
afr1 = afr1 + n(i) * (derta ^ d(i)) * tao ^ T(i)
Next
For i = 8 To 34
afr2 = afr2 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-derta ^ c(i))
Next
For i = 35 To 39
afr3 = afr3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-aerf(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)
Next
CO2_fr = afr1 + afr2 + afr3
End Function
'2.2Φrδ(δ,τ) Eq.(6.5) a Table 31 P34
Public Function CO2_frd(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 7
afrd1 = afrd1 + n(i) * d(i) * (derta ^ (d(i) - 1)) * tao ^ T(i)
Next
For i = 8 To 34
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 = 35 To 39
afrd3 = afrd3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * (2.718282 ^ (-aerf(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (d(i) / derta - 2 * aerf(i) * (derta - ebsl(i)))
Next
CO2_frd = afrd1 + afrd2 + afrd3
End Function
'2.3 Φrδδ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function CO2_frdd(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 7
afrdd1 = afrdd1 + n(i) * d(i) * (d(i) - 1) * (derta ^ (d(i) - 2)) * (tao ^ T(i))
Next
For i = 8 To 34
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
For i = 35 To 39
afrdd3 = afrdd3 + (n(i) * (tao ^ T(i)) * 2.718282 ^ (-aerf(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (-2 * aerf(i) * derta ^ d(i) + 4 * aerf(i) * aerf(i) * (derta ^ d(i)) * (derta - ebsl(i)) * (derta - ebsl(i)) - 4 * d(i) * aerf(i) * (derta ^ (d(i) - 1)) * (derta - ebsl(i)) + d(i) * (d(i) - 1) * derta ^ (d(i) - 2))
Next
CO2_frdd = afrdd1 + afrdd2 + afrdd3
End Function
'2.4 Φrτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function CO2_frt(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 7
afrt1 = afrt1 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1))
Next
For i = 8 To 34
afrt2 = afrt2 + n(i) * T(i) * (derta ^ d(i)) * (tao ^ (T(i) - 1)) * 2.718282 ^ (-derta ^ c(i))
Next
For i = 35 To 39
afrt3 = afrt3 + n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ ((-aerf(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (T(i) / tao - 2 * beita(i) * (tao - gama(i)))
Next
CO2_frt = afrt1 + afrt2 + afrt3
End Function
'2.5 Φrττ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function CO2_frtt(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 7
afrtt1 = afrtt1 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2))
Next
For i = 8 To 34
afrtt2 = afrtt2 + n(i) * T(i) * (T(i) - 1) * (derta ^ d(i)) * (tao ^ (T(i) - 2)) * 2.718282 ^ (-derta ^ c(i))
Next
For i = 35 To 39
afrtt3 = afrtt3 + (n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-aerf(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
CO2_frtt = afrtt1 + afrtt2 + afrtt3
End Function
'2.6 Φrδτ(δ,τ) Eq.(6.5) a Table 31 P38
Public Function CO2_frdt(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 7
afrdt1 = afrdt1 + n(i) * d(i) * T(i) * (derta ^ (T(i) - 1)) * (tao ^ (T(i) - 1))
Next
For i = 8 To 34
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 = 35 To 39
afrdt3 = afrdt3 + (n(i) * (derta ^ d(i)) * (tao ^ T(i)) * 2.718282 ^ (-aerf(i) * (derta - ebsl(i)) ^ 2 - beita(i) * (tao - gama(i)) ^ 2)) * (d(i) / derta - 2 * aerf(i) * (derta - ebsl(i))) * (T(i) / tao - 2 * beita(i) * (tao - gama(i)))
Next
CO2_frdt = afrdt1 + afrdt2 + afrdt3
End Function
'1.Pressure根据密度，温度计算压力
Public Function CO2_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
Dim z As Double
z = CO2_frd(derta, tao)
CO2_rouT_P = rou * R * T * (1 + derta * z) / 1000
End Function
'2.Entropy 根据密度、温度计算熵
Public Function CO2_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 * (CO2_f0t(tao) + CO2_frt(derta, tao)) - CO2_f0(derta, tao) - CO2_fr(derta, tao)
CO2_rouT_S = z * R
End Function
'3.Internal energy 根据密度、温度计算内能
Public Function CO2_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 * (CO2_f0t(tao) + CO2_frt(derta, tao))
CO2_rouT_U = z * R * T
End Function
'4.Cv 根据密度、温度计算定容比容
Public Function CO2_rouT_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 * (CO2_f0tt(tao) + CO2_frtt(derta, tao))
CO2_rouT_Cv = z * R
End Function
'5.Enthalpy 根据密度、温度计算焓
Public Function CO2_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 * (CO2_f0t(tao) + CO2_frt(derta, tao)) + derta * CO2_frd(derta, tao)
CO2_rouT_H = z * R * T
End Function
'6.Isobaric heat capacity 根据密度、温度计算定压热容
'Public Function CO2_rouT_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 * (CO2_f0tt(tao) + CO2_frtt(derta, tao)) + ((1 + derta * CO2_frd(derta, tao) - derta * tao * CO2_frdt(derta, tao)) ^ 2) / (1 + 2 * derta * CO2_frd(derta, tao) + derta * derta * CO2_frdd(derta, tao))
'CO2_rouT_Cp = z * R
'End Function
'饱和液体和饱和蒸汽计算
'3.1 melting-pressure
Public Function CO2_T_Pm(T) As Double
Dim a1, a2, z As Double
a1 = 1955.539
a2 = 2055.4593
z = 1 + a1 * (T / tt - 1) + a2 * (T / tt - 1) ^ 2
CO2_T_Pm = z * pt
End Function
'3.2 sublimation-pressure
Public Function CO2_T_Psub(T) As Double
Dim a1, a2, a3, z As Double
a1 = -14.740846
a2 = 2.4327015
a3 = -5.3061778
z = (tt / T) * (a1 * (1 - T / tt) + a2 * (1 - T / tt) ^ 1.9 + a3 * (1 - T / tt) ^ 2.9)
CO2_T_Psub = Exp(z) * pt
End Function
'3.3 The simple vapor-pressure
Public Function CO2_T_Ps(T) As Double
Dim aaa(1 To 4), tt(1 To 4), z As Double
aaa(1) = -7.0602087: aaa(2) = 1.9391218: aaa(3) = -1.6463597: aaa(4) = -3.2995634
tt(1) = 1: tt(2) = 1.5: tt(3) = 2: tt(4) = 4
z = 0
For i = 1 To 4
z = z + aa(i) * (1 - T / Tc) ^ tt(i)
Next
z = (Tc / T) * z
CO2_T_Ps = Exp(z) * pc
End Function
'3.4 the saturated liquid 饱和液体密度
Public Function CO2_T_roup(T) As Double
Dim aaa(1 To 4), tt(1 To 4), z As Double
aaa(1) = 1.9245108: aaa(2) = -0.62385555: aaa(3) = -0.32731127: aaa(4) = 0.39245142
tt(1) = 0.34: tt(2) = 1 / 2: tt(3) = 10 / 6: tt(4) = 11 / 6
z = 0
For i = 1 To 4
z = z + aaa(i) * (1 - T / Tc) ^ tt(i)
Next
CO2_T_roup = Exp(z) * rouc
End Function
'3.5 the saturated vapor density 饱和蒸汽密度
Public Function CO2_T_roupp(T) As Double
Dim aaa(1 To 4), tt(1 To 4), z As Double
aaa(1) = -1.7074879: aaa(2) = -0.8227467: aaa(3) = -4.6008549
aaa(4) = -10.111178: aaa(5) = -29.742252
tt(1) = 0.34: tt(2) = 1 / 2: tt(3) = 1: tt(4) = 7 / 3: tt(5) = 14 / 3
z = 0
For i = 1 To 5
z = z + aaa(i) * (1 - T / Tc) ^ tt(i)
CO2_T_roupp = Exp(z) * rouc
End Function
'3.6 动力粘度计算
Function yita0(T) As Double
Dim ke, Tx, Ge, a0, a1, a2, a3, a4, e As Double
a0 = 0.235156
a1 = -0.491266
a2 = 0.05211155
a3 = 0.05347906
a4 = -0.01537102
ke = 1 / 251.196
Tx = ke * T
e = a0 + a1 * Log(Tx) + a2 * Log(Tx) ^ 2 + a3 * Log(Tx) ^ 3 + a4 * Log(Tx) ^ 4
Ge = 2.718282 ^ e
yita0 = 1.00697 * T ^ (1 / 2) / Ge
End Function
Function dyita(rou, T) As Double
Dim ke, Tx, d11, d21, d64, d81, d82 As Double
ke = 1 / 251.196
Tx = ke * T
d11 = 0.004071119
d21 = 0.00007198037
d64 = 2.411697E-17
d81 = 2.971072E-23
d82 = -1.627888E-23
dyita = d11 * rou + d21 * rou ^ 2 + d64 * (rou ^ 6) / (Tx * Tx * Tx) + d81 * (rou ^ 8) + d82 * (rou ^ 8) / Tx
End Function
Function CO2_rouT_Yita(rou, T) As Double
CO2_rouT_Yita = yita0(T) + dyita(rou, T)
End Function
Sub test()
Dim rou, T As Double
rou = 407.828
T = 800
Debug.Print "CO2_rouT_P="; CO2_rouT_P(rou, T)
Debug.Print "CO2_rouT_S="; CO2_rouT_S(rou, T)
Debug.Print "CO2_rouT_U="; CO2_rouT_U(rou, T)
Debug.Print "CO2_rouT_Cv="; CO2_rouT_Cv(rou, T)
Debug.Print "CO2_rouT_H="; CO2_rouT_H(rou, T)
'Debug.Print "CO2_rouT_Cp="; CO2_rouT_Cp(rou, T)
Debug.Print "CO2_rouT_Yita="; CO2_rouT_Yita(rou, T)
End Sub
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
