The VBA listing below calculates some of the steam table functions using the 1967 ASME formulation. Note that it calculates in only one direction>> Pressure and Temperature (psia, degF) to the steam properties. It will not calculate "in the other direction." You must iterate (or use Excel's Solver) to do that.
Hope it helps.
Dim qDataCache(100), Cnst(100)
'____________________________ Specific Volume _______________________________
Function fnVsh(P, F)
If Cnst(1) = 0 Then Initialize
B = P / 3208.234759 'Beta (ratio of P/Pcrit)
Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit)
Chi = 4.260321148 * Th / B ' Initialize Chi (ratio of v/vCrit)
Kount1 = 7
For M = 1 To 5
Call Sum2(Sum, Kount1, 0)
Chi = Chi - Sum * M * B ^ (M - 1)
Next
Kount1 = 36: Kount2 = 59
For M = 6 To 8
Call Sum2(Sum, Kount1, 0)
Top = Sum * (M - 2) * B ^ (1 - M)
Call Sum2(Sum, Kount2, 0)
Bot = (Sum + B ^ (2 - M)) ^ 2
Chi = Chi - Top / Bot
Next
Call Sum3(Sum, Bl, 0)
Chi = Chi + 11 * Sum * (B / Bl) ^ 10
fnVsh = Chi * 0.0507785289
End Function
'____________________________________________________________________________
'______________________________ Entropy _____________________________________
Function fnSsh(P, F)
' Private B, Sig, Sum, Kount1, Kount2, Kount3, Kount4, Top, Bot, Bl, M
If Cnst(1) = 0 Then Initialize
B = P / 3208.234759 'Beta (ratio of P/Pcrit)
Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit)
Sig = -4.260321148 * Log(B) + Cnst(1) * Log(Th) ' Initialize Sig
' which is Sigma
Call Sum1(Sum, 0) ' s/sCrit
Sig = Sig - Sum
Kount1 = 7
For M = 1 To 5
Call Sum2(Sum, Kount1, -1)
Sig = Sig - Sum * bLil * B ^ M
Next
Kount1 = 59: Kount2 = 36: Kount3 = 59: Kount4 = 36
For M = 6 To 8
Call Sum2(Top, Kount1, -1)
Call Sum2(Sum, Kount2, 0)
Top = Top + Sum
Call Sum2(Sum, Kount3, 0)
Bot = Sum + B ^ (2 - M)
Call Sum2(Sum, Kount4, -1)
Sig = Sig - bLil * (Sum - Top / Bot) / Bot
Next
Call Sum3(Sum, Bl, 1)
Sig = Sig + B * Sum * (B / Bl) ^ 10
fnSsh = Sig * 0.02587358228
End Function
'____________________________________________________________________________
'_________________________________ Enthalpy _________________________________
Function fnHsh(P, F)
'Private B, Eps, Sum, Kount1, Kount2, Kount3, Kount4, Top, Bot, Buv, Bl, M
If Cnst(1) = 0 Then Initialize
B = P / 3208.234759 'Beta (ratio of P/Pcrit)
Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit)
Eps = Cnst(1) * Th ' Initialize Eps which is Epselon (h/hCrit)
Call Sum1(Sum, -1)
Eps = Eps - Sum
Kount1 = 7: Kount2 = 7
For M = 1 To 5
Call Sum2(Sum, Kount1, 0)
Top = Sum
Call Sum2(Sum, Kount2, -1)
Eps = Eps - (Top + bLil * Th * Sum) * B ^ M
Next
Kount1 = 59: Kount2 = 59: Kount3 = 36: Kount4 = 36
For M = 6 To 8
Call Sum2(Top, Kount1, -1)
Call Sum2(Sum, Kount2, 0)
Bot = Sum + B ^ (2 - M)
Call Sum2(Buv, Kount3, 0)
Eps = Eps - (Buv + bLil * Th * Sum - Buv * bLil * Th * Top / Bot) / Bot
Next
Call Sum3(Sum, Bl, 2)
Eps = Eps + Sum * B * (B / Bl) ^ 10
fnHsh = Eps * 30.14634566
End Function
'____________________________________________________________________________
'______________________________ Sat Pressure ________________________________
Function fnPSat(F)
If F < 32 Then Beep: qPrint "Bad Data": Exit Function
B = P / 3208.234759 'Beta (ratio of P/Pcrit)
Th = (F + 459.67) / 1165.14 'Theta (ratio of degR/degRcrit)
temp = 1 - Th
Top = temp * (temp * (temp * (temp * (temp * -118.9646225 + 64.23285504) _
- 168.1706546) - 26.08023696) - 7.691234564)
Bot = 1 + temp * (temp * 20.9750676 + 4.16711732)
fnPSat = 3208.234759 * Exp(Top / (Th * Bot) - temp / (1000000000# * temp ^ 2 + 6))
End Function
'____________________________________________________________________________
'_______________________________ Sat Temperature ____________________________
Function fnTSat(P)
'Fit by Bize Not ASME Steam Tables
If P < 0.0886 Then Beep: qPrint "Bad Data : exit function"
temp = -268.633 - 0.0054 * P + 2.4417 * Sqr(P) + 32.558 * Sqr(Sqr(P)) _
+ 50.338 * Sqr(Sqr(Sqr(P))) + 285.04 * Sqr(Sqr(Sqr(Sqr(P))))
If P > 1400 Then
Del = P - 1400
temp = temp - Del * (0.000106 + Del * (0.000000343 + Del * 0.00000000013))
End If
fnTSat = temp
End Function
'____________________________________________________________________________
Sub Initialize() ' Reads ASME constants into Cnst() and bLil
'B sub 0,v ===> Cnst(1) thru (6)
qData 16.83599274, 28.56067796, -54.38923329
qData 0.4330662834, -0.6547711697, 0.08565182058
'B sub u,v and z sub u,v where u = 1 to 5
'========> Cnst(7) thru (35) (zero is flag to end loop)
qData 0.06670375918, 13, 1.388983801, 3, 0
qData 0.08390104328, 18, 0.02614670893, 2, -0.03373439453, 1, 0
qData 0.4520918904, 18, 0.1069036614, 10, 0
qData -0.5975336707, 25, -0.08847535804, 14, 0
qData 0.5958051609, 32, -0.5159303373, 28, 0.2075021122, 24, 0
'B sub u,v and z sub u,v where u = 6 to 8
'========> Cnst(36) thru (50) (zero is flag to end loop)
qData 0.1190610271, 12, -0.09867174132, 11, 0
qData 0.1683998803, 24, -0.05809438001, 18, 0
qData 0.006552390126, 24, 0.0005710218649, 14, 0
'B sub 9,v =======> Cnst(51) thru (57) bLil = Cnst(58)
qData 193.6587558, -1388.522425, 4126.607219, -6508.211677
qData 5745.984054, -2693.088365, 523.5718623, 0.7633333333 ' bLil
'b sub u,lambda and x sub u,lambda
'=====> Cnst(59) thru (69) (zero is flag to end loop)
qData 0.4006073948, 14, 0
qData 0.08636081627, 19, 0
qData -0.8532322921, 54, 0.3460208861, 27, 0
For I = 1 To 69
qRead Cnst(I)
Next I
bLil = Cnst(58)
End Sub
Sub NonDim(Psia, degF, Beta, Theta)
' Changes Psia and degF into Non-dimensional parameters
' Beta is ratio of pressure/pressure at critical point of steam
' Theta is ratio of degR/degR at critical point
Beta = Psia / 3208.234759
Theta = (degF + 459.67) / 1165.14
End Sub
Sub Sum1(Sum1, iFlag)
'Sum1 is sum of (N1*Beta(0,v)*Theta^N2) for v = 1 to 5
' where if
' iFlag = -1 then N1 = v-1
' N2 = v-2
' else N1 = v-2
' N2 = v-1
'Private LilV, N1, N2
Sum1 = 0
For LilV = 1 To 5
If iFlag Then
N1 = LilV - 2: N2 = LilV - 1
Else
N1 = LilV - 1: N2 = LilV - 2
End If
Sum1 = Sum1 + N1 * Cnst(LilV + 1) * Th ^ N2
Next
End Sub
Sub Sum2(Sum2, K, iFlag)
'Private Prod
Sum2 = 0
While Cnst(K) <> 0
Prod = Cnst(K) * Exp(Cnst(K + 1) * bLil * (1 - Th))
If iFlag Then Prod = Prod * Cnst(K + 1)
Sum2 = Sum2 + Prod
K = K + 2
Wend
K = K + 1
End Sub
Sub Sum3(Sum3, Bl, iFlag)
'Private LilV, Prod, temp, L1, L2, Blprime
L1 = -34.17061978: L2 = 19.31380707
Bl = 15.74373327 + Th * (L1 + Th * L2)
Blprime = L1 + 2 * L2 * Th
Sum3 = 0
For LilV = 0 To 6
Prod = Cnst(51 + LilV) * Exp(LilV * bLil * (1 - Th))
If iFlag <> 0 Then
temp = 10 * Blprime / Bl + LilV * bLil
If iFlag = 1 Then Prod = Prod * temp
If iFlag = 2 Then Prod = Prod * (1 + Th * temp)
End If
Sum3 = Sum3 + Prod
Next
End Sub