diff --git a/chemicals/dippr.py b/chemicals/dippr.py index 8a0e455..e647ae3 100644 --- a/chemicals/dippr.py +++ b/chemicals/dippr.py @@ -50,7 +50,9 @@ ----------------------- .. autofunction:: chemicals.dippr.EQ101_fitting_jacobian .. autofunction:: chemicals.dippr.EQ102_fitting_jacobian - +.. autofunction:: chemicals.dippr.EQ105_fitting_jacobian +.. autofunction:: chemicals.dippr.EQ106_fitting_jacobian +.. autofunction:: chemicals.dippr.EQ107_fitting_jacobian """ @@ -58,7 +60,9 @@ __all__ = ['EQ100', 'EQ101', 'EQ102', 'EQ104', 'EQ105', 'EQ106', 'EQ107', 'EQ114', 'EQ115', 'EQ116', 'EQ127', - 'EQ101_fitting_jacobian', 'EQ102_fitting_jacobian'] + 'EQ101_fitting_jacobian', 'EQ102_fitting_jacobian', + 'EQ106_fitting_jacobian', 'EQ105_fitting_jacobian', + 'EQ107_fitting_jacobian'] from chemicals.utils import log, exp, sinh, cosh, atan, atanh, sqrt, tanh from cmath import log as clog @@ -404,6 +408,125 @@ def EQ102_fitting_jacobian(Ts, A, B, C, D): out[i][3] = -x2*x5 return out +def EQ105_fitting_jacobian(Ts, A, B, C, D): + r'''Compute and return the Jacobian of the property predicted by + DIPPR Equation # 105 with respect to all the coefficients. This is used in + fitting parameters for chemicals. + + Parameters + ---------- + Ts : list[float] + Temperatures of the experimental data points, [K] + A-D : float + Parameter for the equation; chemical and property specific [-] + + Returns + ------- + jac : list[list[float, 4], len(Ts)] + Matrix of derivatives of the equation with respect to the fitting + parameters, [various] + + ''' + N = len(Ts) +# out = np.zeros((N, 4)) # numba: uncomment + out = [[0.0]*4 for _ in range(N)] # numba: delete + for i in range(N): + r = out[i] + x0 = 1.0 - Ts[i]/C + + if D < 1.0 and x0 < 0.0: + r[0] = 1.0/B + r[1] = -A/(B*B) + else: + x1 = x0**D + x2 = x1 + 1.0 + x3 = A*B**(-x1 - 1.0) + x4 = x1*x3*log(B) + r[0] = B**(-x2) + r[1] = -x2*x3/B + r[2] = -D*Ts[i]*x4/(C*C*x0) + r[3] = -x4*log(x0) + return out + +def EQ106_fitting_jacobian(Ts, Tc, A, B, C, D, E): + r'''Compute and return the Jacobian of the property predicted by + DIPPR Equation # 106 with respect to all the coefficients. This is used in + fitting parameters for chemicals. + + Parameters + ---------- + Ts : list[float] + Temperatures of the experimental data points, [K] + Tc : float + Critical temperature, [K] + A-E : float + Parameter for the equation; chemical and property specific [-] + + Returns + ------- + jac : list[list[float, 5], len(Ts)] + Matrix of derivatives of the equation with respect to the fitting + parameters, [various] + + ''' + N = len(Ts) +# out = np.zeros((N, 5)) # numba: uncomment + out = [[0.0]*5 for _ in range(N)] # numba: delete + for i in range(N): + x0 = Ts[i]/Tc + if x0 != 1.0: + x1 = 1.0 - x0 + x2 = x1**(B + x0*(C + x0*(D + E*x0))) + x3 = A*x2*log(x1) + r = out[i] + r[0] = x2 + r[1] = x3 + r[2] = x0*x3 + r[3] = x0*x0*x3 + r[4] = x0*x0*x0*x3 + return out + +def EQ107_fitting_jacobian(Ts, A, B, C, D, E): + r'''Compute and return the Jacobian of the property predicted by + DIPPR Equation # 107 with respect to all the coefficients. This is used in + fitting parameters for chemicals. + + Parameters + ---------- + Ts : list[float] + Temperatures of the experimental data points, [K] + A-E : float + Parameter for the equation; chemical and property specific [-] + + Returns + ------- + jac : list[list[float, 5], len(Ts)] + Matrix of derivatives of the equation with respect to the fitting + parameters, [various] + + ''' + N = len(Ts) +# out = np.zeros((N, 5)) # numba: uncomment + out = [[0.0]*5 for _ in range(N)] # numba: delete + for i in range(N): + r = out[i] + x1 = 1.0/Ts[i] + x0 = x1*x1 + x2 = C*x1 + x3 = sinh(x2) + x3_inv = 1.0/x3 + x4 = x0*x3_inv*x3_inv + x5 = E*x1 + x6 = cosh(x5) + x6_inv = 1.0/x6 + x7 = x0*x6_inv*x6_inv + r[0] = 1.0 + r[1] = C*C*x4 + r[2] = 2.0*B*C*x4*(-x2*cosh(x2)*x3_inv + 1.0) + r[3] = E*E*x7 + r[4] = 2.0*D*E*x7*(-x5*sinh(x5)*x6_inv + 1.0) + return out + def EQ104(T, A, B, C=0.0, D=0.0, E=0.0, order=0): r'''DIPPR Equation #104. Often used in calculating second virial coefficients of gases. All 5 parameters are required. @@ -788,7 +911,11 @@ def EQ107(T, A=0, B=0, C=0, D=0, E=0, order=0): doi:10.1016/0378-3812(81)85002-9. ''' if order == 0: - return A + B*((C/T)/sinh(C/T))**2 + D*((E/T)/cosh(E/T))**2 + C_T = C/T + t0 = 2.0*C_T/(trunc_exp(C_T) - trunc_exp(-C_T)) + E_T = E/T + t1 = 2.0*E_T/(trunc_exp(-E_T) + trunc_exp(E_T)) + return A + B*t0*t0 + D*t1*t1 elif order == 1: return (2*B*C**3*cosh(C/T)/(T**4*sinh(C/T)**3) - 2*B*C**2/(T**3*sinh(C/T)**2) diff --git a/chemicals/numba.py b/chemicals/numba.py index 592240c..87e6024 100644 --- a/chemicals/numba.py +++ b/chemicals/numba.py @@ -75,9 +75,13 @@ def transform_complete_chemicals(replaced, __funcs, __all__, normal, vec): 'vapor_pressure.Wagner_fitting_jacobian', 'vapor_pressure.Wagner_original_fitting_jacobian', + 'vapor_pressure.Antoine_fitting_jacobian', 'dippr.EQ102_fitting_jacobian', 'dippr.EQ101_fitting_jacobian', + 'dippr.EQ106_fitting_jacobian', + 'dippr.EQ105_fitting_jacobian', + 'dippr.EQ107_fitting_jacobian', 'vapor_pressure.Yaws_Psat_fitting_jacobian', 'viscosity.mu_Yaws_fitting_jacobian', diff --git a/chemicals/vapor_pressure.py b/chemicals/vapor_pressure.py index 3c50f24..8235366 100644 --- a/chemicals/vapor_pressure.py +++ b/chemicals/vapor_pressure.py @@ -57,6 +57,7 @@ ----------------------- .. autofunction:: chemicals.vapor_pressure.Wagner_fitting_jacobian .. autofunction:: chemicals.vapor_pressure.Wagner_original_fitting_jacobian +.. autofunction:: chemicals.vapor_pressure.Antoine_fitting_jacobian .. autofunction:: chemicals.vapor_pressure.Yaws_Psat_fitting_jacobian Vapor Pressure Estimation Correlations @@ -165,6 +166,7 @@ 'Antoine_coeffs_from_point', 'Antoine_AB_coeffs_from_point', 'DIPPR101_ABC_coeffs_from_point', 'Wagner_original_fitting_jacobian', 'Wagner_fitting_jacobian', 'Yaws_Psat_fitting_jacobian', + 'Antoine_fitting_jacobian', 'TDE_PVExpansion'] import os @@ -1327,8 +1329,8 @@ def Wagner_original_fitting_jacobian(Ts, Tc, Pc, a, b, c, d): Parameters ---------- - T : float - Temperature of fluid, [K] + Ts : list[float] + Temperatures of fluid data points, [K] Tc : float Critical temperature, [K] Pc : float @@ -1369,8 +1371,8 @@ def Wagner_fitting_jacobian(Ts, Tc, Pc, a, b, c, d): Parameters ---------- - T : float - Temperature of fluid, [K] + Ts : list[float] + Temperatures of fluid data points, [K] Tc : float Critical temperature, [K] Pc : float @@ -1402,6 +1404,50 @@ def Wagner_fitting_jacobian(Ts, Tc, Pc, a, b, c, d): row[3] = x4*x5 return out +def Antoine_fitting_jacobian(Ts, A, B, C, base=10.0): + r'''Calculates the jacobian of the Antoine vapor pressure equation + for use in fitting these parameters when experimental values are known. + + Requires three coefficients specific to each chemical. + + Parameters + ---------- + Ts : list[float] + Temperatures of fluid data points, [K] + A : float + Antoine `A` parameter, [-] + B : float + Antoine `B` parameter, [K] + C : float + Antoine `C` parameter, [K] + base : float, optional + Optional base of logarithm; 10 by default, [-] + + Returns + ------- + jac : list[list[float, 3], len(Ts)] + Matrix of derivatives of the equation with respect to the fitting + parameters, [various] + ''' + N = len(Ts) +# out = np.zeros((N, 3)) # numba: uncomment + out = [[0.0]*3 for _ in range(N)] # numba: delete + ln_base = log(base) + for i in range(N): + row = out[i] + x0 = C + Ts[i] + if x0 <= 0.0: + row[0] = 0.0 + row[1] = 0.0 + row[2] = 0.0 + else: + x1 = 1.0/x0 + x2 = base**(A - B*x1)*ln_base + row[0] = x2 + row[1] = -x1*x2 + row[2] = B*x2*x1*x1 + return out + def Wagner(T, Tc, Pc, a, b, c, d): r'''Calculates vapor pressure using the Wagner equation (2.5, 5 form). diff --git a/tests/test_dippr.py b/tests/test_dippr.py index 7d6fd92..cd56e5d 100644 --- a/tests/test_dippr.py +++ b/tests/test_dippr.py @@ -227,6 +227,10 @@ def test_EQ107_more(): with pytest.raises(Exception): EQ107(20., *coeffs, order=1E100) + + # Case that requires overflow handling + EQ107(**{'T': 377.77777777777777, 'A': 1539249.2020718465, 'B': -46807441.804555826, 'C': -409401.9169728528, 'D': -2164118.45731599, 'E': 339.5030595758336, 'order': 0}) + def test_EQ114_more(): # T derivative @@ -367,3 +371,56 @@ def test_EQ102_fitting(): der_analytical = EQ102_fitting_jacobian([T], A, B, C, D) assert_close1d(der_analytical, [der_num]) assert_close1d(der_analytical, der_expect, rtol=1e-13) + +def test_EQ105_fitting(): + T, A, B, C, D = 300., 0.70824, 0.26411, 507.6, 0.27537 + der_num = [derivative(lambda A: EQ105(T, A, B, C, D), A, dx=A*1e-5), + derivative(lambda B: EQ105(T, A, B, C, D), B, dx=B*1e-5), + derivative(lambda C: EQ105(T, A, B, C, D), C, dx=C*1e-5), + derivative(lambda D: EQ105(T, A, B, C, D), D, dx=D*1e-5)] + der_expect = [[10.721182221195127, -51.225752844842916, 0.006195731841673147, -7.0661094492142285]] + der_analytical = EQ105_fitting_jacobian([T], A, B, C, D) + assert_close1d(der_analytical, [der_num]) + assert_close1d(der_analytical, der_expect, rtol=1e-13) + + + T, A, B, C, D = 304, 3733.087734888731, 0.30552803535622014, 301.6993863907116, 0.3415888512743092 + der_num = [derivative(lambda A: EQ105(T, A, B, C, D), A, dx=A*1e-5), + derivative(lambda B: EQ105(T, A, B, C, D), B, dx=B*1e-5), + derivative(lambda C: EQ105(T, A, B, C, D), C, dx=C*1e-5), + derivative(lambda D: EQ105(T, A, B, C, D), D, dx=D*1e-5)] + der_analytical = EQ105_fitting_jacobian([T], A, B, C, D) + assert_close1d(der_analytical, [der_num]) + + +def test_EQ106_fitting(): + T, Tc, A, B, C, D, E = 300, 647.096, 0.17766, 2.567, -3.3377, 1.9699, 0.25 + der_num = [derivative(lambda A: EQ106(T, Tc, A, B, C, D, E), A, dx=A*1e-5), + derivative(lambda B: EQ106(T, Tc, A, B, C, D, E), B, dx=B*1e-5), + derivative(lambda C: EQ106(T, Tc, A, B, C, D, E), C, dx=C*1e-5), + derivative(lambda D: EQ106(T, Tc, A, B, C, D, E), D, dx=D*1e-5), + derivative(lambda E: EQ106(T, Tc, A, B, C, D, E), E, dx=E*1e-5), + ] + + der_expect = [[0.4007741423076755, -0.04435095583995359, -0.020561534535812425, -0.009532527415937865, -0.004419372434354963]] + der_analytical = EQ106_fitting_jacobian([T], Tc, A, B, C, D, E) + assert_close1d(der_analytical, [der_num]) + assert_close1d(der_analytical, der_expect, rtol=1e-13) + + # Test case with errors + EQ106_fitting_jacobian(Ts=[466.],Tc = 466.0, A = 47700.0, B = 0.37, C = 0.0, D = 0.0, E = 0.0) + +def test_EQ107_fitting(): + T, A, B, C, D, E = 250.0, 33363., 26790., 2610.5, 8896., 1169 + der_num = [derivative(lambda A: EQ107(T, A, B, C, D, E), A, dx=A*1e-5), + derivative(lambda B: EQ107(T, A, B, C, D, E), B, dx=B*3e-3, order=3), + derivative(lambda C: EQ107(T, A, B, C, D, E), C, dx=C*1e-3, order=9), + derivative(lambda D: EQ107(T, A, B, C, D, E), D, dx=D*1e-5, order=3), + derivative(lambda E: EQ107(T, A, B, C, D, E), E, dx=E*1e-5), + ] + + der_expect = [[1.0, 3.7138247806865474e-07, -7.197214038362036e-05, 0.00758947296962729, -0.42452325330497365]] + der_analytical = EQ107_fitting_jacobian([T], A, B, C, D, E) + assert_close1d(der_analytical, [der_num]) + assert_close1d(der_analytical, der_expect, rtol=1e-13) + diff --git a/tests/test_vapor_pressure.py b/tests/test_vapor_pressure.py index cb5b254..5c0e97c 100644 --- a/tests/test_vapor_pressure.py +++ b/tests/test_vapor_pressure.py @@ -378,6 +378,27 @@ def test_Wagner_fitting_jacobian(): assert_close1d(der_expect, der_analytical, rtol=1e-13) assert_close1d(der_analytical, [der_num], rtol=1e-7) +def test_Antoine_fitting_jacobian(): + + T, A, B, C = 100.0, 8.7687, 395.744, -6.469 + der_num = [derivative(lambda A: Antoine(T, A, B, C), A, dx=A*1e-5), + derivative(lambda B: Antoine(T, A, B, C), B, dx=B*1e-5), + derivative(lambda C: Antoine(T, A, B, C), C, dx=C*1e-5)] + der_expect = [[79389.37469005348, -848.802800034785, 3591.414774748115]] + der_analytical = Antoine_fitting_jacobian([T], A, B, C) + assert_close1d(der_analytical, [der_num], rtol=1e-7) + assert_close1d(der_expect, der_analytical, rtol=1e-13) + + # Zero point + T, A, B, C = 30, 3.45604+5, 1044.038, -53.893 + der_num = [derivative(lambda A: Antoine(T, A, B, C), A, dx=A*1e-5), + derivative(lambda B: Antoine(T, A, B, C), B, dx=B*1e-5), + derivative(lambda C: Antoine(T, A, B, C), C, dx=C*1e-5)] + der_expect = [[0.0, 0.0, 0.0]] + der_analytical = Antoine_fitting_jacobian([T], A, B, C) + assert_close1d(der_analytical, [der_num], rtol=1e-7) + assert_close1d(der_expect, der_analytical, rtol=1e-13) + def test_Yaws_Psat(): assert_close(Yaws_Psat(T=400.0, A=28.588+ log10(101325/760), B=-2469, C=-7.351, D=2.8025E-10, E=2.7361E-6), 708657.0891069275, rtol=1e-14)