/[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 2225 - (show annotations) (download) (as text)
Mon Jul 26 03:40:18 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 15503 byte(s)
Added python bindings for new function. 
We can see that it's not converging for about 10% of cases.
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
30 #include "helmholtz_impl.h"
31 # include <assert.h>
32
33 #include <math.h>
34 #ifdef TEST
35 # include <stdio.h>
36 # include <stdlib.h>
37 #endif
38
39 #include <gsl/gsl_multiroots.h>
40
41 #define SQ(X) ((X)*(X))
42
43 /**
44 Estimate of saturation pressure using H W Xiang ''The new simple extended
45 corresponding-states principle: vapor pressure and second virial
46 coefficient'', Chemical Engineering Science,
47 57 (2002) pp 1439-1449.
48 */
49 double fprops_psat_T_xiang(double T, const HelmholtzData *d){
50
51 double Zc = d->p_c / (8314. * d->rho_c * d->T_c);
52
53 #ifdef TEST
54 fprintf(stderr,"Zc = %f\n",Zc);
55 #endif
56
57 double theta = SQ(Zc - 0.29);
58
59 #ifdef TEST
60 fprintf(stderr,"theta = %f\n",theta);
61 #endif
62
63 double aa[] = {5.790206, 4.888195, 33.91196
64 , 6.251894, 15.08591, -315.0248
65 , 11.65859, 46.78273, -1672.179
66 };
67
68 double a0 = aa[0] + aa[1]*d->omega + aa[2]*theta;
69 double a1 = aa[3] + aa[4]*d->omega + aa[5]*theta;
70 double a2 = aa[6] + aa[7]*d->omega + aa[8]*theta;
71
72 double Tr = T / d->T_c;
73 double tau = 1 - Tr;
74
75 double taupow = pow(tau, 1.89);
76 double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow;
77
78 double logpr = temp * log(Tr);
79 double p_r = exp(logpr);
80
81 #ifdef TEST
82 fprintf(stderr,"a0 = %f\n", a0);
83 fprintf(stderr,"a1 = %f\n", a1);
84 fprintf(stderr,"a2 = %f\n", a2);
85 fprintf(stderr,"temp = %f\n", temp);
86 fprintf(stderr,"T_r = %f\n", Tr);
87 fprintf(stderr,"p_r = %f\n", p_r);
88 #endif
89
90 return p_r * d->p_c;
91 }
92
93 /**
94 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
95 see http://dx.doi.org/10.1002/aic.690250412
96 */
97 double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
98
99 double Zc = D->rho_c * D->R * D->T_c / D->p_c;
100 double Tau = 1. - T/D->T_c;
101 double vf = (D->R * D->T_c / D->p_c) * pow(Zc, -1 - pow(Tau, 2./7));
102
103 return 1./vf;
104 }
105
106
107 /**
108 Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
109 see http://dx.doi.org/10.1016/j.tca.2004.05.017
110 */
111 double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
112 double Zc = D->rho_c * D->R * D->T_c / D->p_c;
113 double Tau = 1. - T/D->T_c;
114 #if 0
115 # define N1 -0.1497547
116 # define N2 0.6006565
117 # define P1 -19.348354
118 # define P2 -41.060325
119 # define P3 1.1878726
120 double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */
121 //MMM = 2.4686277;
122 double PPP = Zc / (P1 + P2*Zc*log(Zc) + P3/Zc);
123 fprintf(stderr,"PPP = %f\n",PPP);
124 //PPP = -0.6240188;
125 double NNN = PPP + 1./(N1*D->omega + N2);
126 #else
127 /* exact values from Chouaieb for CO2 */
128 # define MMM 2.4686277
129 # define NNN 1.1345838
130 # define PPP -0.6240188
131 #endif
132
133 double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM));
134 return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
135 }
136
137
138 /**
139 Maxwell phase criterion, first principles
140 */
141 void phase_criterion(double T, double rho_f, double rho_g, double p_sat, double *eq1, double *eq2, double *eq3, const HelmholtzData *D){
142 #ifdef TEST
143 fprintf(stderr,"PHASE CRITERION: T = %f, rho_f = %f, rho_g = %f, p_sat = %f\n", T, rho_f, rho_g, p_sat);
144 #endif
145 double delta_f, delta_g, tau;
146 tau = D->T_c / T;
147
148 delta_f = rho_f / D->rho_c;
149 delta_g = rho_g/ D->rho_c;
150
151 #ifdef TEST
152 assert(!isnan(delta_f));
153 assert(!isnan(delta_g));
154 assert(!isnan(p_sat));
155 assert(!isinf(delta_f));
156 assert(!isinf(delta_g));
157 assert(!isinf(p_sat));
158 #endif
159
160 *eq1 = (p_sat - helmholtz_p(T, rho_f, D));
161 *eq2 = (p_sat - helmholtz_p(T, rho_g, D));
162 *eq3 = helmholtz_g(T, rho_f, D) - helmholtz_g(T,rho_g, D);
163
164 #ifdef TEST
165 fprintf(stderr,"eq1 = %e\t\teq2 = %e\t\teq3 = %e\n", *eq1, *eq2, *eq3);
166 #endif
167 }
168
169
170 typedef struct PhaseSolve_struct{
171 double tau;
172 const HelmholtzData *D;
173 const double *scaling;
174 } PhaseSolve;
175
176 /**
177 Maxwell phase criterion, from equations from IAPWS-95.
178
179 We define here pi = p_sat / (R T rho_c)
180 */
181 static int phase_resid(const gsl_vector *x, void *params, gsl_vector *f){
182 const double tau = ((PhaseSolve *)params)->tau;
183 const double pi = gsl_vector_get(x,0);
184 double delta_f = gsl_vector_get(x,1);
185 double delta_g = exp(gsl_vector_get(x,2));
186 const HelmholtzData *D = ((PhaseSolve *)params)->D;
187 const double *scaling = ((PhaseSolve *)params)->scaling;
188
189 assert(!isinf(pi));
190 assert(!isinf(delta_f));
191 assert(!isinf(delta_g));
192
193 gsl_vector_set(f, 0, scaling[0] * (delta_f + delta_f * delta_f * helm_resid_del(tau, delta_f, D) - pi));
194 gsl_vector_set(f, 1, scaling[1] * (1 + delta_g * helm_resid_del(tau, delta_g, D) - pi/delta_g));
195 gsl_vector_set(f, 2, scaling[2] * (pi * (delta_f - delta_g) - delta_f * delta_g * (log(delta_f / delta_g) + helm_resid(delta_g,tau,D) - helm_resid(delta_f,tau,D))));
196
197
198 if(isnan(gsl_vector_get(f,2)) || isnan(gsl_vector_get(f,2)) || isnan(gsl_vector_get(f,2))){
199 fprintf(stderr,"NaN encountered... ");
200 }
201
202 //gsl_vector_set(f, 2, pi * (delta_f - delta_g) * delta_f*delta_g* (log(delta_f/delta_g) + helm_resid(delta_g,tau,D) - helm_resid(delta_f,tau,D)));
203
204 return GSL_SUCCESS;
205 }
206
207 static int print_state(size_t iter, gsl_multiroot_fsolver * s){
208 fprintf(stderr,"iter = %3u: pi = %.3f delf = %.3f delg = %.3f "
209 "E = % .3e % .3e %.3e\n",
210 iter,
211 gsl_vector_get(s->x, 0),
212 gsl_vector_get(s->x, 1),
213 gsl_vector_get(s->x, 2),
214 gsl_vector_get (s->f, 0),
215 gsl_vector_get (s->f, 1),
216 gsl_vector_get (s->f, 2));
217 }
218
219 #if 0
220 int phase_solve(double T, double *p_sat, double *rho_f, double *rho_g, const HelmholtzData *D){
221 gsl_multiroot_fsolver *s;
222 int status;
223 const size_t n = 3;
224 double tau = D->T_c / T;
225 size_t i, iter = 0;
226 const gsl_multiroot_fsolver_type *ftype;
227
228 const double scaling[3] = {1., 1., 0.1};
229
230 PhaseSolve p = {
231 tau
232 , D
233 , scaling
234 };
235
236 gsl_multiroot_function f = {
237 &phase_resid,
238 n, &p
239 };
240
241 /* TODO use the first guess 'provided' if psat, rho_f, rho_g less than zero? */
242
243 double x_init[3] = {
244 /* first guesses, such as they are... */
245 fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c
246 ,fprops_rhof_T_rackett(T, D) / D->rho_c
247 ,log(fprops_rhog_T_chouaieb(T, D) / D->rho_c)
248 };
249
250 gsl_vector *x = gsl_vector_alloc (n);
251 for(i=0; i<3; ++i)gsl_vector_set(x, i, x_init[i]);
252 ftype = gsl_multiroot_fsolver_broyden;
253 s = gsl_multiroot_fsolver_alloc(ftype, n);
254 gsl_multiroot_fsolver_set(s, &f, x);
255
256 //print_state(iter, s);
257
258 double pi_xiang = fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c;
259 double delg_chouaieb = fprops_rhog_T_chouaieb(T, D) / D->rho_c;
260 double delf_rackett = fprops_rhof_T_rackett(T, D) / D->rho_c;
261
262 assert(!isnan(pi_xiang));
263 assert(!isnan(delg_chouaieb));
264 assert(!isnan(delf_rackett));
265
266 do{
267 iter++;
268
269 status = gsl_multiroot_fsolver_iterate(s);
270
271 print_state(iter, s);
272 if(status)break;
273
274 status = gsl_multiroot_test_residual(s->f, 1e-7);
275
276 double pi = gsl_vector_get(s->x,0);
277 double delf = gsl_vector_get(s->x,1);
278 double delg = exp(gsl_vector_get(s->x,2));
279
280
281 #define VAR_PI 0
282 #define VAR_DELF 1
283 #define VAR_DELG 2
284 #define CHECK_RESET(COND, VAR, NEWVAL)\
285 if(COND){\
286 gsl_vector_set(s->x,VAR,NEWVAL);\
287 fprintf(stderr,"RESET %s to %s = %f (FAILED %s)\n", #VAR, #NEWVAL, NEWVAL, #COND);\
288 }
289
290 #define CHECK_RESET_CONTINUE(COND, VAR, NEWVAL)\
291 if(COND){\
292 gsl_vector_set(s->x,VAR,NEWVAL);\
293 fprintf(stderr,"RESET %s to %s = %f (FAILED %s)\n", #VAR, #NEWVAL, NEWVAL, #COND);\
294 continue;\
295 }
296
297 //CHECK_RESET_CONTINUE(pi < 0, VAR_PI, pi_xiang);
298 //CHECK_RESET_CONTINUE(delg < 0, VAR_DELG, delg_chouaieb);
299 //CHECK_RESET_CONTINUE(delf < 0, VAR_DELF, delf_rackett);
300
301 #if 0
302 CHECK_RESET_CONTINUE(delf < 1, VAR_DELF, delf_rackett);
303 CHECK_RESET_CONTINUE(delg > 1, VAR_DELG, delg_chouaieb);
304
305 CHECK_RESET_CONTINUE(pi > 2. * pi_xiang, VAR_PI, pi_xiang)
306 else CHECK_RESET_CONTINUE(pi < 0.5 * pi_xiang, VAR_PI, pi_xiang)
307
308 CHECK_RESET_CONTINUE(delg < 0.8 * delg_chouaieb, VAR_DELG, delg_chouaieb)
309 else CHECK_RESET_CONTINUE(delg > 1.5 * delg_chouaieb, VAR_DELG, delg_chouaieb);
310
311 //if(gsl_vector_get(s->x,1) < D->rho_c)gsl_vector_set(s->x,1,2 * D->rho_c);
312 //if(gsl_vector_get(s->x,2) > D->rho_c)gsl_vector_set(s->x,2,0.5 * D->rho_c);
313
314 #endif
315 }while(status == GSL_CONTINUE && iter < 40);
316
317 if(status!=GSL_SUCCESS)fprintf(stderr,"%s:%d: warning: (%d) status = %s\n", __FILE__,__LINE__, status, gsl_strerror (status));
318
319 //fprintf(stderr,"SOLUTION: pi = %f\n", gsl_vector_get(s->x, 0));
320
321
322 //fprintf(stderr," p = %f\n", p_sat);
323
324 *p_sat = gsl_vector_get(s->x, 0) * T * D->R * D->rho_c;
325 *rho_f = gsl_vector_get(s->x, 1) * D->rho_c;
326 *rho_g = gsl_vector_get(s->x, 2) * D->rho_c;
327
328 gsl_multiroot_fsolver_free(s);
329 gsl_vector_free(x);
330
331 return status;
332 }
333 #endif
334
335 #if 0
336
337
338 #ifdef TEST
339 fprintf(stderr,"PHASE CRITERION: T = %f, rho_f = %f, rho_g = %f, p_sat = %f\n", T, rho_f, rho_g, p_sat);
340 #endif
341 double delta_f, delta_g, tau;
342 tau = D->T_c / T;
343
344 delta_f = rho_f / D->rho_c;
345 delta_g = rho_g/ D->rho_c;
346
347 #ifdef TEST
348 assert(!isnan(delta_f));
349 assert(!isnan(delta_g));
350 assert(!isnan(p_sat));
351 assert(!isinf(delta_f));
352 assert(!isinf(delta_g));
353 assert(!isinf(p_sat));
354 #endif
355
356 *eq1 = (p_sat - helmholtz_p(T, rho_f, D));
357 *eq2 = (p_sat - helmholtz_p(T, rho_g, D));
358 *eq3 = helmholtz_g(T, rho_f, D) - helmholtz_g(T,rho_g, D);
359
360 #ifdef TEST
361 fprintf(stderr,"eq1 = %e\t\teq2 = %e\t\teq3 = %e\n", *eq1, *eq2, *eq3);
362 #endif
363 }
364
365
366
367 #if 0
368 static int phase_deriv(const gsl_vector * x, void *params, gsl_matrix * J){
369 const double tau = ((PhaseSolve *)params)->tau;
370 const double pi = gsl_vector_get(x,0);
371 const double delta_f = gsl_vector_get(x,1);
372 const double delta_g = gsl_vector_get(x,2);
373 const HelmholtzData *D = ((PhaseSolve *)params)->D;
374
375 double dE1ddelf = helm_resid_del(tau, delta_f, D) + delta_f * helm_resid_deldel(tau, delta_f, D) - pi ;
376 double dE1ddelg = 0;
377 double dE1dpi = - 1./ delta_f;
378
379 double dE2ddelf = 0;
380 double dE2ddelg = helm_resid_del(tau, delta_g, D) + delta_g * helm_resid_deldel(tau, delta_g, D) - pi ;
381 double dE2dpi = - 1./ delta_g;
382
383 double dE3ddelf = pi - delta_g * (log(delta_f/delta_g) + helm_resid(delta_g, tau, D) - helm_resid(delta_f, tau, D) + 1 - delta_f * helm_resid_del(delta_f, tau, D));
384 double dE3ddelg = -pi - delta_f * (log(delta_f/delta_g) + helm_resid(delta_g, tau, D) - helm_resid(delta_f, tau, D) - 1 + delta_g * helm_resid_del(delta_g, tau, D));
385 double dE3dpi = delta_f - delta_g;
386
387 gsl_matrix_set (J, 0, 0, dE1ddelf);
388 gsl_matrix_set (J, 0, 1, dE1ddelg);
389 gsl_matrix_set (J, 0, 2, dE1dpi);
390 gsl_matrix_set (J, 1, 0, dE2ddelf);
391 gsl_matrix_set (J, 1, 1, dE2ddelg);
392 gsl_matrix_set (J, 1, 2, dE2dpi);
393 gsl_matrix_set (J, 2, 0, dE3ddelf);
394 gsl_matrix_set (J, 2, 1, dE3ddelg);
395 gsl_matrix_set (J, 2, 2, dE3dpi);
396
397 return GSL_SUCCESS;
398 }
399
400 static int phase_residderiv(
401 const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J
402 ){
403 phase_resid(x, params, f);
404 phase_deriv(x, params, J);
405 return GSL_SUCCESS;
406 }
407
408 static int print_state_fdf(size_t iter, gsl_multiroot_fdfsolver * s){
409 fprintf(stderr,"iter = %3u: delf = %.3f delg = %.3f pi = %.3f "
410 "E = % .3e % .3e %.3e\n",
411 iter,
412 gsl_vector_get (s->x, 0),
413 gsl_vector_get (s->x, 1),
414 gsl_vector_get (s->x, 2),
415 gsl_vector_get (s->f, 0),
416 gsl_vector_get (s->f, 1),
417 gsl_vector_get (s->f, 2));
418 }
419
420
421 double phase_solve_fdf(double T, const HelmholtzData *D){
422 gsl_multiroot_fdfsolver *s;
423 int status;
424 const size_t n = 3;
425 double tau = D->T_c / T;
426 size_t i, iter = 0;
427 const gsl_multiroot_fdfsolver_type *fdftype;
428
429 PhaseSolve p = {tau, D};
430 gsl_multiroot_function_fdf f = {
431 &phase_resid,
432 &phase_deriv,
433 &phase_residderiv,
434 n, &p
435 };
436
437 double x_init[3] = {
438 /* first guesses, such as they are... */
439 fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c
440 ,fprops_rhof_T_rackett(T, D) / D->rho_c
441 ,fprops_rhog_T_chouaieb(T, D) / D->rho_c
442 };
443
444 gsl_vector *x = gsl_vector_alloc (n);
445 for(i=0; i<3; ++i)gsl_vector_set(x, i, x_init[i]);
446 fdftype = gsl_multiroot_fdfsolver_hybridsj;
447 s = gsl_multiroot_fdfsolver_alloc (fdftype, n);
448 gsl_multiroot_fdfsolver_set (s, &f, x);
449
450 print_state_fdf(iter, s);
451
452 do{
453 iter++;
454 status = gsl_multiroot_fdfsolver_iterate (s);
455 print_state_fdf(iter, s);
456 if(status)break;
457 if(gsl_vector_get(s->x,2) > D->rho_c)gsl_vector_set(s->x,2,D->rho_c);
458 status = gsl_multiroot_test_residual(s->f, 1e-7);
459 }while(status == GSL_CONTINUE && iter < 1000);
460
461 fprintf(stderr,"status = %s\n", gsl_strerror (status));
462
463 fprintf(stderr,"SOLUTION: pi = %f\n", gsl_vector_get(s->f, 0));
464
465 double p_sat = gsl_vector_get(s->f, 0) * T * D->R * D->rho_c;
466
467 fprintf(stderr," p = %f\n", p_sat);
468 gsl_multiroot_fdfsolver_free (s);
469 gsl_vector_free (x);
470
471 return 0;
472 }
473 #endif
474
475
476 void solve_saturation(double T, const HelmholtzData *D){
477 double rho_f, rho_g, p_sat;
478
479 /* first guesses, such as they are... */
480 p_sat = fprops_psat_T_xiang(T, D);
481 rho_f = fprops_rhof_T_rackett(T, D);
482 rho_g = fprops_rhog_T_chouaieb(T, D);
483
484
485
486 #ifdef TEST
487 fprintf(stderr,"Rackett liquid density: %f\n", rho_f);
488 fprintf(stderr,"Chouaieb vapour density: %f\n", rho_g);
489 #endif
490
491 double delta_f, delta_g, eq1, eq2, eq3, tau;
492
493 int i = 40;
494 while(--i > 0){
495 delta_f = rho_f / D->rho_c;
496 delta_g = rho_g/ D->rho_c;
497
498 #ifdef TEST
499 assert(!isnan(delta_f));
500 assert(!isnan(delta_g));
501 assert(!isnan(p_sat));
502 assert(!isinf(delta_f));
503 assert(!isinf(delta_g));
504 assert(!isinf(p_sat));
505 #endif
506
507 phase_criterion(T, rho_f, rho_g, p_sat, &eq1, &eq2, &eq3, D);
508
509 #ifdef TEST
510 fprintf(stderr,"p_sat = %f\trho_f = %f\trho_g = %f\teq1 = %e\t\teq2 = %e\t\teq3 = %e\n",p_sat, rho_f, rho_g, eq1, eq2, eq3);
511 assert(!isnan(eq1));
512 assert(!isnan(eq2));
513 assert(!isnan(eq3));
514 assert(!isinf(eq1));
515 assert(!isinf(eq2));
516 assert(!isinf(eq3));
517 #endif
518
519 double p1 = D->R * T * D->rho_c * delta_f * delta_g / (delta_f - delta_g) * (helm_resid(delta_f,tau,D) - helm_resid(delta_g,tau,D) + log(delta_f/delta_g));
520 if(p1/p_sat > 1.5){
521 p1 = p_sat * 1.5;
522 }else if(p1/p_sat < 0.7){
523 p1 = p_sat * 0.7;
524 }
525 p_sat = p1;
526 if(p_sat < D->p_t) p_sat = D->p_t;
527 if(p_sat > D->p_c) p_sat = D->p_c;
528
529 double rho_f1 = rho_f - (helmholtz_p(T, rho_f, D) - p_sat) / helmholtz_dpdrho_T(T, rho_f, D);
530 double rho_g1 = rho_g - (helmholtz_p(T, rho_g, D) - p_sat) / helmholtz_dpdrho_T(T, rho_g, D);
531
532 if(rho_g1 > 0){
533 rho_g = rho_g1;
534 }
535 if(rho_f1 > 0){
536 rho_f = rho_f1;
537 }
538 if(fabs(rho_f - rho_g) < 1e-5){
539 #ifdef TEST
540 fprintf(stderr,"%s:%d: rho_f = rho_g\n", __FILE__, __LINE__);
541 exit(1);
542 #else
543 break;
544 #endif
545 }
546
547
548
549 }
550
551 //return eq3;
552 }
553
554
555 #endif
556
557
558
559

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