/[ascend]/trunk/models/johnpye/fprops/sat.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/sat.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2301 - (show annotations) (download) (as text)
Sat Aug 21 13:28:35 2010 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 12783 byte(s)
Regen toluene model working, next water.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2009 Carnegie Mellon University
3
4 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 */
27
28 #include "sat.h"
29 #include "helmholtz_impl.h"
30
31 //#define SAT_DEBUG
32
33 #ifdef SAT_DEBUG
34 # include <assert.h>
35 #else
36 # define assert(ARGS...)
37 #endif
38
39 #include <math.h>
40 #include <stdio.h>
41
42 #define SQ(X) ((X)*(X))
43
44 //#define THROW_FPE
45
46 #ifdef THROW_FPE
47 #define _GNU_SOURCE
48 #include <fenv.h>
49 int feenableexcept (int excepts);
50 int fedisableexcept (int excepts);
51 int fegetexcept (void);
52 #endif
53
54 #ifdef SAT_DEBUG
55 # define MSG(STR,...) fprintf(stderr,"%s:%d: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
56 #else
57 # define MSG(ARGS...)
58 #endif
59 #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
60
61 /**
62 Estimate of saturation pressure using H W Xiang ''The new simple extended
63 corresponding-states principle: vapor pressure and second virial
64 coefficient'', Chemical Engineering Science,
65 57 (2002) pp 1439-1449.
66 */
67 double fprops_psat_T_xiang(double T, const HelmholtzData *d){
68 double p_c = fprops_pc(d);
69 double Zc = p_c / (8314. * d->rho_c * d->T_c);
70
71 #ifdef TEST
72 fprintf(stderr,"Zc = %f\n",Zc);
73 #endif
74
75 double theta = SQ(Zc - 0.29);
76
77 #ifdef TEST
78 fprintf(stderr,"theta = %f\n",theta);
79 #endif
80
81 double aa[] = {5.790206, 4.888195, 33.91196
82 , 6.251894, 15.08591, -315.0248
83 , 11.65859, 46.78273, -1672.179
84 };
85
86 double a0 = aa[0] + aa[1]*d->omega + aa[2]*theta;
87 double a1 = aa[3] + aa[4]*d->omega + aa[5]*theta;
88 double a2 = aa[6] + aa[7]*d->omega + aa[8]*theta;
89
90 double Tr = T / d->T_c;
91 double tau = 1 - Tr;
92
93 double taupow = pow(tau, 1.89);
94 double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow;
95
96 double logpr = temp * log(Tr);
97 double p_r = exp(logpr);
98
99 #ifdef TEST
100 fprintf(stderr,"a0 = %f\n", a0);
101 fprintf(stderr,"a1 = %f\n", a1);
102 fprintf(stderr,"a2 = %f\n", a2);
103 fprintf(stderr,"temp = %f\n", temp);
104 fprintf(stderr,"T_r = %f\n", Tr);
105 fprintf(stderr,"p_r = %f\n", p_r);
106 #endif
107
108 return p_r * p_c;
109 }
110
111 /**
112 Estimate saturation pressure using acentric factor. This algorithm
113 is used for first estimates for later refinement in the program REFPROP.
114 */
115 double fprops_psat_T_acentric(double T, const HelmholtzData *d){
116 /* first guess using acentric factor */
117 double p_c = fprops_pc(d);
118 double p = p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
119 return p;
120 }
121
122
123 /**
124 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
125 see http://dx.doi.org/10.1002/aic.690250412
126 */
127 double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
128 double p_c = fprops_pc(D);
129 double Zc = D->rho_c * D->R * D->T_c / p_c;
130 double Tau = 1. - T/D->T_c;
131 double vf = (D->R * D->T_c / p_c) * pow(Zc, -1 - pow(Tau, 2./7));
132
133 return 1./vf;
134 }
135
136 /**
137 Inverse of fprops_rhof_T_rackett. FIXME this need checking.
138 */
139 double fprops_T_rhof_rackett(double rhof, const HelmholtzData *D){
140 double p_c = fprops_pc(D);
141 double Zc = D->rho_c * D->R * D->T_c / p_c;
142 double f1 = p_c / D->R / D->T_c / rhof;
143 double f2 = -log(f1)/log(Zc);
144 return pow(f2 -1, 3./2);
145 }
146
147 /**
148 Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
149 see http://dx.doi.org/10.1016/j.tca.2004.05.017
150 */
151 double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
152 double Tau = 1. - T/D->T_c;
153 #if 0
154 double p_c = fprops_pc(D);
155 double Zc = D->rho_c * D->R * D->T_c / p_c;
156 # define N1 -0.1497547
157 # define N2 0.6006565
158 # define P1 -19.348354
159 # define P2 -41.060325
160 # define P3 1.1878726
161 double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */
162 //MMM = 2.4686277;
163 double PPP = Zc / (P1 + P2*Zc*log(Zc) + P3/Zc);
164 fprintf(stderr,"PPP = %f\n",PPP);
165 //PPP = -0.6240188;
166 double NNN = PPP + 1./(N1*D->omega + N2);
167 #else
168 /* exact values from Chouaieb for CO2 */
169 # define MMM 2.4686277
170 # define NNN 1.1345838
171 # define PPP -0.6240188
172 #endif
173
174 double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM));
175 return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
176 }
177
178 /**
179 Solve saturation condition for a specified temperature.
180 @param T temperature [K]
181 @param psat_out output, saturation pressure [Pa]
182 @param rhof_out output, saturated liquid density [kg/m^3]
183 @param rhog_out output, saturated vapour density [kg/m^3]
184 @param d helmholtz data object for the fluid in question.
185 @return 0 on success, non-zero on error (eg algorithm failed to converge, T out of range, etc.)
186 */
187 int fprops_sat_T(double T, double *psat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
188 double tau = d->T_c / T;
189 double delf = 1.1 * fprops_rhof_T_rackett(T,d) / d->rho_c;
190 double delg = 0.9 * fprops_rhog_T_chouaieb(T,d) / d->rho_c;
191 #ifdef THROW_FPE
192 feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
193 #endif
194 #ifdef SAT_DEBUG
195 fprintf(stderr,"%s: calculating for %s, T = %.12e\n",__func__,d->name,T);
196 #endif
197
198 if(T < d->T_t - 1e-8){
199 ERRMSG("Input temperature is below triple-point temperature");
200 return 1;
201 }
202
203 if(fabs(T - d->T_c) < 1e-9){
204 *psat_out = fprops_pc(d);
205 *rhof_out = d->rho_c;
206 *rhog_out = d->rho_c;
207 return 0;
208 }
209
210 int i = 0;
211 while(i++ < 20){
212 assert(!__isnan(delg));
213 #ifdef SAT_DEBUG
214 fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delf*d->rho_c, delg*d->rho_c);
215 #endif
216 double phirf = helm_resid(tau,delf,d);
217 double phirf_d = helm_resid_del(tau,delf,d);
218 double phirf_dd = helm_resid_deldel(tau,delf,d);
219 double phirg = helm_resid(tau,delg,d);
220 double phirg_d = helm_resid_del(tau,delg,d);
221 double phirg_dd = helm_resid_deldel(tau,delg,d);
222 assert(!__isnan(phirf));
223 assert(!__isnan(phirf_d));
224 assert(!__isnan(phirf_dd));
225 assert(!__isnan(phirg));
226 assert(!__isnan(phirg_d));
227 assert(!__isnan(phirg_dd));
228
229 #define J(FG) (del##FG * (1. + del##FG * phir##FG##_d))
230 #define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG))
231 #define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd)
232 #define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG)
233 double Jf = J(f);
234 double Jg = J(g);
235 double Kf = K(f);
236 double Kg = K(g);
237 double Jf_del = J_del(f);
238 double Jg_del = J_del(g);
239 double Kf_del = K_del(f);
240 double Kg_del = K_del(g);
241
242 double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
243 assert(!__isnan(DELTA));
244
245 #define gamma 1.0
246 delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
247 delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);
248
249 assert(!__isnan(delg));
250 assert(!__isnan(delf));
251
252 if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){
253 //fprintf(stderr,"%s: CONVERGED\n",__func__);
254 *rhof_out = delf * d->rho_c;
255 *rhog_out = delg * d->rho_c;
256 if(__isnan(*rhog_out)){
257 fprintf(stderr,"%s: T = %.12e\n",__func__,T);
258 }
259 *psat_out = helmholtz_p_raw(T, *rhog_out, d);
260 return 0;
261 }
262 if(delg < 0)delg = -0.5*delg;
263 if(delf < 0)delf = -0.5*delf;
264 }
265 *rhof_out = delf * d->rho_c;
266 *rhog_out = delg * d->rho_c;
267 *psat_out = helmholtz_p_raw(T, *rhog_out, d);
268 ERRMSG("Not converged: '%s' with T = %e (rhof=%f, rhog=%f).",d->name,T,*rhof_out,*rhog_out);
269 return 1;
270
271 }
272
273 /**
274 Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.
275 */
276 double fprops_pc(const HelmholtzData *d){
277 static const HelmholtzData *d_last = NULL;
278 static double p_c = 0;
279 if(d == d_last){
280 return p_c;
281 }
282 p_c = helmholtz_p_raw(d->T_c, d->rho_c,d);
283 d_last = d;
284 return p_c;
285 }
286
287 /**
288 Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.
289 */
290 int fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const HelmholtzData *d){
291 static const HelmholtzData *d_last = NULL;
292 static double p_t, rhof_t, rhog_t;
293 if(d == d_last){
294 *p_t_out = p_t;
295 *rhof_t_out = rhof_t;
296 *rhog_t_out = rhog_t;
297 return 0;
298 }
299 MSG("Calculating saturation for T = %f",d->T_t);
300 int res = fprops_sat_T(d->T_t, &p_t, &rhof_t, &rhog_t,d);
301 if(res)return res;
302 else{
303 d_last = d;
304 *p_t_out = p_t;
305 *rhof_t_out = rhof_t;
306 *rhog_t_out = rhog_t;
307 return 0;
308 }
309 }
310
311
312 /**
313 Solve saturation properties in terms of pressure.
314 This function makes calls to fprops_sat_T, and solves for temperature using
315 a Newton solver algorith. Derivatives dp/dT are calculated using the
316 Clapeyron equation.
317 @return 0 on success.
318 */
319 int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
320 MSG("Calculating for %s at p = %.12e Pa",d->name,p);
321 double T1;
322 double p_c = fprops_pc(d);
323 if(fabs(p - p_c)/p_c < 1e-6){
324 MSG("Very close to critical pressure: using critical temperature without iteration.");
325 T1 = d->T_c;
326 double p1, rhof, rhog;
327 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
328 *Tsat_out = T1;
329 *rhof_out = rhof;
330 *rhog_out = rhog;
331 return res;
332 }else{
333 /*
334 Estimate of saturation temperature using definition of acentric factor and
335 the assumed p(T) relationship:
336 log10(p)=A + B/T
337 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
338 */
339 T1 = d->T_c / (1. - 3./7. / (1.+d->omega) * log10(p / p_c));
340 MSG("Estimated using acentric factor: T = %f",T1);
341 if(T1 < d->T_t){
342 T1 = d->T_t;
343 MSG("Estimate moved up to T_t = %f",T1);
344 }
345 }
346 double p1, rhof, rhog;
347 int i = 0;
348 while(i++ < 50){
349 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
350 if(res){
351 MSG("Got error %d from fprops_sat_T at T = %.12e", res,T1);
352 return 1;
353 }
354 MSG("T1 = %f ������> p = %f bar\trhof = %f\trhog = %f",T1, p1/1e5, rhof, rhog);
355 if(fabs(p1 - p) < 1e-5){
356 *Tsat_out = T1;
357 *rhof_out = rhof;
358 *rhog_out = rhog;
359 return 0;
360 }
361 double hf = helmholtz_h_raw(T1, rhof, d);
362 double hg = helmholtz_h_raw(T1, rhog, d);
363 double dpdT_sat = (hg - hf) / T1 / (1./rhog - 1./rhof);
364 //fprintf(stderr,"\t\tdpdT_sat = %f bar/K\n",dpdT_sat/1e5);
365 double delta_T = -(p1 - p)/dpdT_sat;
366 if(T1 + delta_T < d->T_t - 1e-2){
367 MSG("Correcting sub-triple-point temperature guess");
368 T1 = 0.5 * (d->T_t + T1);
369 }
370 else T1 += delta_T;
371 }
372 MSG("Exceeded iteration limit, returning last guess with error code");
373 *Tsat_out = T1;
374 *rhof_out = rhof;
375 *rhog_out = rhog;
376 return 1;
377 }
378
379
380
381 /**
382 Calculate Tsat based on a value of hf. This value is useful in setting
383 first guess Temperatures when solving for the coordinates (p,h).
384 This function uses the secant method for the iterative solution.
385 */
386 int fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
387 double T1 = 0.4 * d->T_t + 0.6 * d->T_c;
388 double T2 = d->T_t;
389 double h1, h2, p, rhof, rhog;
390 int res = fprops_sat_T(T2, &p, &rhof, &rhog, d);
391 if(res){
392 ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,d->name);
393 return 1;
394 }
395 double tol = 1e-6;
396 h2 = helmholtz_h(T2,rhof,d);
397 if(hf < h2){
398 ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
399 return 2;
400 }
401
402 int i = 0;
403 while(i++ < 60){
404 assert(T1 >= d->T_t - 1e-4);
405 assert(T1 <= d->T_c);
406 MSG("T1 = %f\n",T1);
407 res = fprops_sat_T(T1, &p, &rhof, &rhog, d);
408 if(res){
409 ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,d->name);
410 return 1;
411 }
412 h1 = helmholtz_h(T1,rhof, d);
413 if(fabs(h1 - hf) < tol){
414 *Tsat_out = T1;
415 *psat_out = p;
416 *rhof_out = rhof;
417 *rhog_out = rhog;
418 return 0;
419 }
420 if(h1 == h2){
421 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",d->name,h1,hf);
422 return 2;
423 }
424
425 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
426 T2 = T1;
427 h2 = h1;
428 while(T1 + delta_T > d->T_c)delta_T *= 0.5;
429 T1 += delta_T;
430 if(T1 < d->T_t)T1 = d->T_t;
431 if(i==20 || i==30)tol*=100;
432 }
433 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
434 *Tsat_out = T1;
435 *psat_out = p;
436 *rhof_out = rhof;
437 *rhog_out = rhog;
438 return 1;
439 }
440
441
442

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22