1 |
jpye |
2122 |
/* ASCEND modelling environment |
2 |
|
|
Copyright (C) 2008-2009 Carnegie Mellon University |
3 |
jpye |
2114 |
|
4 |
jpye |
2122 |
This program is free software; you can redistribute it and/or modify |
5 |
|
|
it under the terms of the GNU General Public License as published by |
6 |
|
|
the Free Software Foundation; either version 2, or (at your option) |
7 |
|
|
any later version. |
8 |
|
|
|
9 |
|
|
This program is distributed in the hope that it will be useful, |
10 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
11 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
12 |
|
|
GNU General Public License for more details. |
13 |
|
|
|
14 |
|
|
You should have received a copy of the GNU General Public License |
15 |
|
|
along with this program; if not, write to the Free Software |
16 |
|
|
Foundation, Inc., 59 Temple Place - Suite 330, |
17 |
|
|
Boston, MA 02111-1307, USA. |
18 |
|
|
*//** @file |
19 |
|
|
Routines to calculate saturation properties using Helmholtz correlation |
20 |
|
|
data. We first include some 'generic' saturation equations that make use |
21 |
|
|
of the acentric factor and critical point properties to predict saturation |
22 |
|
|
properties (pressure, vapour density, liquid density). These correlations |
23 |
|
|
seem to be only very rough in some cases, but it is hoped that they will |
24 |
|
|
be suitable as first-guess values that can then be fed into an iterative |
25 |
|
|
solver to converge to an accurate result. |
26 |
jpye |
2114 |
*/ |
27 |
|
|
|
28 |
|
|
#include "sat.h" |
29 |
jpye |
2262 |
#include "helmholtz_impl.h" |
30 |
jpye |
2114 |
|
31 |
jpye |
2124 |
# include <assert.h> |
32 |
jpye |
2114 |
#include <math.h> |
33 |
jpye |
2262 |
#include <stdio.h> |
34 |
jpye |
2114 |
|
35 |
|
|
#define SQ(X) ((X)*(X)) |
36 |
|
|
|
37 |
jpye |
2262 |
#if 0 |
38 |
|
|
#define _GNU_SOURCE |
39 |
|
|
#include <fenv.h> |
40 |
|
|
#endif |
41 |
|
|
|
42 |
jpye |
2122 |
/** |
43 |
|
|
Estimate of saturation pressure using H W Xiang ''The new simple extended |
44 |
|
|
corresponding-states principle: vapor pressure and second virial |
45 |
|
|
coefficient'', Chemical Engineering Science, |
46 |
|
|
57 (2002) pp 1439-1449. |
47 |
|
|
*/ |
48 |
jpye |
2114 |
double fprops_psat_T_xiang(double T, const HelmholtzData *d){ |
49 |
|
|
|
50 |
|
|
double Zc = d->p_c / (8314. * d->rho_c * d->T_c); |
51 |
|
|
|
52 |
jpye |
2117 |
#ifdef TEST |
53 |
jpye |
2114 |
fprintf(stderr,"Zc = %f\n",Zc); |
54 |
jpye |
2117 |
#endif |
55 |
jpye |
2114 |
|
56 |
|
|
double theta = SQ(Zc - 0.29); |
57 |
|
|
|
58 |
jpye |
2117 |
#ifdef TEST |
59 |
jpye |
2114 |
fprintf(stderr,"theta = %f\n",theta); |
60 |
jpye |
2117 |
#endif |
61 |
jpye |
2114 |
|
62 |
|
|
double aa[] = {5.790206, 4.888195, 33.91196 |
63 |
|
|
, 6.251894, 15.08591, -315.0248 |
64 |
|
|
, 11.65859, 46.78273, -1672.179 |
65 |
|
|
}; |
66 |
|
|
|
67 |
|
|
double a0 = aa[0] + aa[1]*d->omega + aa[2]*theta; |
68 |
|
|
double a1 = aa[3] + aa[4]*d->omega + aa[5]*theta; |
69 |
|
|
double a2 = aa[6] + aa[7]*d->omega + aa[8]*theta; |
70 |
|
|
|
71 |
|
|
double Tr = T / d->T_c; |
72 |
|
|
double tau = 1 - Tr; |
73 |
|
|
|
74 |
|
|
double taupow = pow(tau, 1.89); |
75 |
|
|
double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow; |
76 |
|
|
|
77 |
|
|
double logpr = temp * log(Tr); |
78 |
|
|
double p_r = exp(logpr); |
79 |
|
|
|
80 |
|
|
#ifdef TEST |
81 |
|
|
fprintf(stderr,"a0 = %f\n", a0); |
82 |
|
|
fprintf(stderr,"a1 = %f\n", a1); |
83 |
|
|
fprintf(stderr,"a2 = %f\n", a2); |
84 |
|
|
fprintf(stderr,"temp = %f\n", temp); |
85 |
|
|
fprintf(stderr,"T_r = %f\n", Tr); |
86 |
|
|
fprintf(stderr,"p_r = %f\n", p_r); |
87 |
|
|
#endif |
88 |
|
|
|
89 |
|
|
return p_r * d->p_c; |
90 |
|
|
} |
91 |
|
|
|
92 |
jpye |
2115 |
/** |
93 |
jpye |
2230 |
Estimate saturation pressure using acentric factor. This algorithm |
94 |
|
|
is used for first estimates for later refinement in the program REFPROP. |
95 |
|
|
*/ |
96 |
|
|
double fprops_psat_T_acentric(double T, const HelmholtzData *d){ |
97 |
|
|
/* first guess using acentric factor */ |
98 |
|
|
double p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.)); |
99 |
|
|
return p; |
100 |
|
|
} |
101 |
|
|
|
102 |
|
|
|
103 |
|
|
/** |
104 |
jpye |
2120 |
Saturated liquid density correlation of Rackett, Spencer & Danner (1972) |
105 |
|
|
see http://dx.doi.org/10.1002/aic.690250412 |
106 |
|
|
*/ |
107 |
jpye |
2121 |
double fprops_rhof_T_rackett(double T, const HelmholtzData *D){ |
108 |
jpye |
2120 |
|
109 |
|
|
double Zc = D->rho_c * D->R * D->T_c / D->p_c; |
110 |
|
|
double Tau = 1. - T/D->T_c; |
111 |
|
|
double vf = (D->R * D->T_c / D->p_c) * pow(Zc, -1 - pow(Tau, 2./7)); |
112 |
|
|
|
113 |
|
|
return 1./vf; |
114 |
|
|
} |
115 |
|
|
|
116 |
|
|
|
117 |
|
|
/** |
118 |
|
|
Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi |
119 |
|
|
see http://dx.doi.org/10.1016/j.tca.2004.05.017 |
120 |
|
|
*/ |
121 |
jpye |
2121 |
double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){ |
122 |
jpye |
2120 |
double Zc = D->rho_c * D->R * D->T_c / D->p_c; |
123 |
|
|
double Tau = 1. - T/D->T_c; |
124 |
|
|
#if 0 |
125 |
|
|
# define N1 -0.1497547 |
126 |
|
|
# define N2 0.6006565 |
127 |
|
|
# define P1 -19.348354 |
128 |
|
|
# define P2 -41.060325 |
129 |
|
|
# define P3 1.1878726 |
130 |
|
|
double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */ |
131 |
jpye |
2122 |
//MMM = 2.4686277; |
132 |
|
|
double PPP = Zc / (P1 + P2*Zc*log(Zc) + P3/Zc); |
133 |
|
|
fprintf(stderr,"PPP = %f\n",PPP); |
134 |
|
|
//PPP = -0.6240188; |
135 |
jpye |
2120 |
double NNN = PPP + 1./(N1*D->omega + N2); |
136 |
|
|
#else |
137 |
jpye |
2122 |
/* exact values from Chouaieb for CO2 */ |
138 |
jpye |
2120 |
# define MMM 2.4686277 |
139 |
|
|
# define NNN 1.1345838 |
140 |
|
|
# define PPP -0.6240188 |
141 |
|
|
#endif |
142 |
|
|
|
143 |
|
|
double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM)); |
144 |
jpye |
2121 |
return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha))); |
145 |
jpye |
2120 |
} |
146 |
|
|
|
147 |
jpye |
2262 |
int fprops_sat_T(double T, double *psat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){ |
148 |
|
|
double tau = d->T_c / T; |
149 |
|
|
double delf = 1.1 * fprops_rhof_T_rackett(T,d) / d->rho_c; |
150 |
|
|
double delg = 0.9 * fprops_rhog_T_chouaieb(T,d) / d->rho_c; |
151 |
|
|
//feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); |
152 |
jpye |
2120 |
|
153 |
jpye |
2262 |
int i = 0; |
154 |
|
|
while(i++ < 50){ |
155 |
|
|
//fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delf*d->rho_c, delg*d->rho_c); |
156 |
|
|
double phirf = helm_resid(tau,delf,d); |
157 |
|
|
double phirf_d = helm_resid_del(tau,delf,d); |
158 |
|
|
double phirf_dd = helm_resid_deldel(tau,delf,d); |
159 |
|
|
double phirg = helm_resid(tau,delg,d); |
160 |
|
|
double phirg_d = helm_resid_del(tau,delg,d); |
161 |
|
|
double phirg_dd = helm_resid_deldel(tau,delg,d); |
162 |
jpye |
2119 |
|
163 |
jpye |
2262 |
#define J(FG) (del##FG * (1. + del##FG * phir##FG##_d)) |
164 |
|
|
#define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG)) |
165 |
|
|
#define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd) |
166 |
|
|
#define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG) |
167 |
|
|
double Jf = J(f); |
168 |
|
|
double Jg = J(g); |
169 |
|
|
double Kf = K(f); |
170 |
|
|
double Kg = K(g); |
171 |
|
|
double Jf_del = J_del(f); |
172 |
|
|
double Jg_del = J_del(g); |
173 |
|
|
double Kf_del = K_del(f); |
174 |
|
|
double Kg_del = K_del(g); |
175 |
jpye |
2119 |
|
176 |
jpye |
2262 |
double DELTA = Jg_del * Kf_del - Jf_del * Kg_del; |
177 |
jpye |
2122 |
|
178 |
jpye |
2262 |
#define gamma 1 |
179 |
|
|
delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del); |
180 |
|
|
delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del); |
181 |
jpye |
2119 |
|
182 |
jpye |
2262 |
if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){ |
183 |
|
|
//fprintf(stderr,"%s: CONVERGED\n",__func__); |
184 |
|
|
*rhof_out = delf * d->rho_c; |
185 |
|
|
*rhog_out = delg * d->rho_c; |
186 |
|
|
*psat_out = helmholtz_p_raw(T, *rhog_out, d); |
187 |
|
|
return 0; |
188 |
|
|
} |
189 |
jpye |
2124 |
} |
190 |
jpye |
2262 |
*rhof_out = delf * d->rho_c; |
191 |
|
|
*rhog_out = delg * d->rho_c; |
192 |
|
|
*psat_out = helmholtz_p_raw(T, *rhog_out, d); |
193 |
|
|
fprintf(stderr,"%s: NOT CONVERGED for '%s' with T = %e (rhof=%f, rhog=%f)\n",__func__,d->name,T,*rhof_out,*rhog_out); |
194 |
|
|
return 1; |
195 |
jpye |
2124 |
|
196 |
jpye |
2122 |
} |
197 |
|
|
|
198 |
jpye |
2262 |
int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){ |
199 |
|
|
fprintf(stderr,"%s: NOT YET IMPLEMENTED\n",__func__); |
200 |
|
|
return 1; |
201 |
jpye |
2123 |
} |
202 |
|
|
|
203 |
|
|
|
204 |
jpye |
2124 |
|
205 |
|
|
|
206 |
jpye |
2123 |
|
207 |
jpye |
2124 |
|