1 #include "LogK_function_of_T.h"
2 #include "ThermoParameters.h"
3 #include "Reaction.h"
4
5 namespace ThermoFun {
6
thermoPropertiesReaction_LogK_fT(Reaktoro_::Temperature TK,Reaktoro_::Pressure Pbar,Reaction reaction,MethodCorrT_Thrift::type CE)7 auto thermoPropertiesReaction_LogK_fT(Reaktoro_::Temperature TK, Reaktoro_::Pressure Pbar, Reaction reaction, MethodCorrT_Thrift::type CE) -> ThermoPropertiesReaction
8 {
9 ThermoPropertiesReaction tpr;
10
11 auto Rln10 = R_CONSTANT * lg_to_ln;
12 auto R_T = TK * R_CONSTANT;
13 auto ref_tpr = reaction.thermoReferenceProperties();
14 auto A = reaction.thermoParameters().reaction_logK_fT_coeff;
15 auto CpCoeff = reaction.thermoParameters().reaction_Cp_fT_coeff;
16 auto dVr = ref_tpr.reaction_volume; //Gr = rc[q].Gs[0];
17 auto dHr = ref_tpr.reaction_enthalpy;
18 auto dSr = ref_tpr.reaction_entropy;
19 auto dCpr = ref_tpr.reaction_heat_capacity_cp;
20 auto dGr = ref_tpr.reaction_gibbs_energy;
21 auto lgK = ref_tpr.log_equilibrium_constant;
22 auto Tr = reaction.referenceT();
23
24 /// deal with Cp and A parameters conversion
25 switch (CE)
26 {
27 case MethodCorrT_Thrift::type::CTM_EK0: // 1-term lgK = const
28 // lgK = A[0];
29 dCpr = 0.0;
30 dHr = 0.0;
31 // dGr = -dSr * TK;
32 break;
33 case MethodCorrT_Thrift::type::CTM_EK1: // 1-term dGr = const
34 // lgK = A[2]/T;
35 dCpr = 0.0;
36 dSr = 0.0;
37 dHr = dGr;
38 // dGr_d = dHr;
39 lgK = - dHr / TK / Rln10;
40 break;
41 case MethodCorrT_Thrift::type::CTM_EK2: // 2-term or 1-term lgK=const at dHr=0
42 // lgK = A[0] + A[2]/T;
43 dCpr = 0.0;
44 // dGr_d = dHr - dSr * T;
45 lgK = (dSr - dHr/TK ) / Rln10;
46 break;
47 case MethodCorrT_Thrift::type::CTM_EK3: // 3-term
48 // lgK = A[0] + A[2]/T + A[3] * lnT;
49 lgK = ( dSr - dHr/TK - dCpr * ( 1 - Tr/TK - log( TK/Tr ))) / Rln10;
50 dSr += dCpr * log( TK / Tr );
51 dHr += dCpr * (TK - Tr );
52 // dGr_d = dHr - dSr * T;
53 break;
54 case MethodCorrT_Thrift::type::CTM_LGK: // full 7-term logK approx
55 case MethodCorrT_Thrift::type::CTM_LGX: // (derived from dCp=f(T))
56 if (A.size() == 0 && CpCoeff.size() > 0)
57 {
58 //from Cp to A
59 A.resize(7, 0.0);
60 CpCoeff.resize(5, 0.0);
61
62 // calculation of logK=f(T) coeffs (only first 5 Cp coefficients, conforming to Haas-Fisher function)
63 A[0] = (( dSr - CpCoeff[0] - CpCoeff[0]*log(TK) - CpCoeff[1]*TK + CpCoeff[2]/(2.0*TK*TK)
64 + 2.0*CpCoeff[3]/pow(TK,0.5) - CpCoeff[4]*TK*TK/2.0 ) / Rln10).val;
65 A[1] = CpCoeff[1]/(2.0*Rln10);
66 A[2] = (-( dHr - CpCoeff[0]*TK - CpCoeff[1]*TK*TK/2.0 + CpCoeff[2]/TK
67 - 2.0*CpCoeff[3]*pow(TK,0.5) - CpCoeff[4]*TK*TK*TK/3.0 ) / Rln10).val;
68 A[3] = CpCoeff[0]/Rln10;
69 A[4] = CpCoeff[2]/(2.0*Rln10);
70 A[5] = CpCoeff[4]/(6.0*Rln10);
71 A[6] = -4.0*CpCoeff[3]/Rln10;
72 }
73
74 if (CpCoeff.size() == 0 && A.size() > 0)
75 {
76 //from A to Cp
77 A.resize(7, 0.0);
78 CpCoeff.resize(5, 0.0);
79
80 CpCoeff[0] = Rln10 * A[3];
81 CpCoeff[1] = Rln10 * 2.0 * A[1];
82 CpCoeff[2] = Rln10 * 2.0 * A[4];
83 CpCoeff[3] = -Rln10 * 0.25 * A[6];
84 CpCoeff[4] = Rln10 * 6.0 * A[5];
85 }
86
87 lgK = A[0] + A[1] * TK + A[2]/TK + A[3] * log(TK) + A[4] / (TK*TK) +
88 A[5] * (TK*TK) + A[6] / pow(TK,0.5);
89 dHr = Rln10 *( A[1]*(TK*TK) - A[2] + A[3]*TK - 2.0*A[4]/TK
90 + 2.0*A[5]*(TK*TK*TK) - 0.5*A[6]*pow(TK,0.5) );
91 dSr = Rln10 * ( A[0] + 2.0*A[1]*TK + A[3]*(1.0+log(TK)) -
92 A[4]/(TK*TK) + 3.0*A[5]*(TK*TK) + 0.5*A[6]/pow(TK,0.5) );
93 // dGr_d = dHr - dSr * T;
94 // if( rc[q].DCp )
95 dCpr = CpCoeff[0] + CpCoeff[1]*TK + CpCoeff[2]/(TK*TK) + CpCoeff[4]*(TK*TK) + CpCoeff[3]/pow(TK,0.5);
96 // dHr = Rln10 * T_2 * ( A[1] - A[2]/T_2 + A[3]/T -
97 // 2.0*A[4]/T_3 + 2.0*A[5]*T - 0.5*A[6]/T_15 );
98 // dSr = Rln10 * ( A[0] + 2.0*A[1]*T + A[3]*(1.0+lnT) -
99 // A[4]/T_2 + 3.0*A[5]*T_2 + 0.5*A[6]/T_05 );
100 break;
101 default:
102 ; // error message ?
103 }
104 // Calculation of dGr
105 dGr = -R_T * lgK * lg_to_ln;
106 // if (dHr.val == 0) dHr = dGr + dSr * TK;
107 Reaktoro_::ThermoScalar dUr, dAr;
108 if (dHr.val != 0)
109 {
110 dUr = dHr - Pbar*dVr;
111 dAr = dUr - TK*dSr;
112 }
113
114 // Loading output data
115 tpr.log_equilibrium_constant = lgK;
116 tpr.ln_equilibrium_constant = lgK*lg_to_ln;
117 tpr.reaction_gibbs_energy = dGr;
118 tpr.reaction_enthalpy = dHr;
119 tpr.reaction_entropy = dSr;
120 tpr.reaction_heat_capacity_cp = dCpr;
121 tpr.reaction_volume = dVr;
122 tpr.reaction_helmholtz_energy = dAr;
123 tpr.reaction_internal_energy = dUr;
124
125 return tpr;
126 }
127
128 }
129