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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2260 - (show annotations) (download) (as text)
Thu Aug 5 06:28:01 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 12729 byte(s)
Error messages to detect temperatures below critical point.
1 #include "sat2.h"
2 #include "helmholtz_impl.h"
3
4 #include <assert.h>
5 #include <math.h>
6 #include <stdio.h>
7
8 #define PHASE_DEBUG
9 #define PHASE_ERRORS
10
11 #ifndef PHASE_DEBUG
12 # define MSG(...)
13 #else
14 # define MSG(ARGS...) fprintf(stderr,"%s:%d: ",__func__,__LINE__);fprintf(stderr,ARGS)
15 #endif
16
17 #ifndef PHASE_DEBUG
18 # define ERR(...)
19 #else
20 # define ERR(ARGS...) fprintf(stderr,"%s: ",__func__);fprintf(stderr,ARGS)
21 #endif
22
23
24 /**
25 solve Maxwell phase criterion using successive substitution
26
27 @return 0 on success, else error code.
28 @param p output, pressure
29 @param rhof output, liquid density
30 @param rhof output, vapour density
31 */
32 int fprops_sat_T(double T, double *p_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
33 double p, p_new, delta_p, delta_p_old;
34
35 if(T < d->T_t){
36 ERR("Temperature T = %e K is below triple point, unable to calculate saturation properties.\n",T);
37 return 5;
38 }
39
40 /*
41 Estimate of saturation temperature using definition of acentric factor and
42 the assumed p(T) relationship:
43 log10(p)=A + B/T
44 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
45 */
46 p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
47
48 MSG("Initial estimate: psat(%f) = %f bar\n",T,p/1e5);
49
50 int i,j;
51 char use_guess = 0;
52
53 #define FPROPS_MAX_SUCCSUBS 12
54
55 double rhof = -1, rhog = -1;
56
57 /* this approach follows approximately the approach of Soave, page 204:
58 http://dx.doi.org/10.1016/0378-3812(86)90013-0
59 */
60 for(i=0;i<FPROPS_MAX_SUCCSUBS;++i){
61 if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof)){
62 MSG(" increasing p, error with rho_f\n");
63 p *= 1.005;
64 continue;
65 }
66
67 if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
68 MSG(" decreasing p, error with rho_g\n");
69 p *= 0.95;
70 continue;
71 }
72
73 MSG("Iter %d: p = %f bar, rhof = %f, rhog = %f\n",i, p/1e5, rhof, rhog);
74
75 use_guess = 1; /* after first run, start re-using current guess */
76
77 if(fabs(rhof - rhog) < 1e-8){
78 ERR("densities converged to same value\n");
79 *p_out = p;
80 *rhof_out = rhof;
81 *rhog_out = rhog;
82 return 1;
83 }
84
85 double delta_a = helmholtz_a(T, rhof, d) - helmholtz_a(T, rhog, d);
86 MSG(" delta_a = %f\n",delta_a);
87
88 p_new = delta_a / (1./rhog - 1./rhof);
89 delta_p = p_new - p;
90
91 MSG(" delta_p = %f bar\n",delta_p/1e5);
92
93 /* convergence test */
94 if(fabs(delta_p/p) < 1e-6){
95 MSG("CONVERGED to p = %f bar...\n",p/1e5);
96 /* note possible need for delp non-change test for certain fluids */
97
98 p = p_new;
99
100 /* find vapour density, using guess */
101 if(fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
102 ERR("failed to estimate vapour density\n");
103 *rhof_out = rhof;
104 *rhog_out = rhog;
105 *p_out = p;
106 return 2;
107 }
108 /* find liquid phase, using guess */
109 fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof);
110
111 /* re-evaluate exact p for the value of rhof calculated */
112 p = helmholtz_p_raw(T,rhof,d);
113
114 if(p > d->p_c || rhof < d->rho_c){
115 ERR("invalid converged value of p, beyond p_crit\n");
116 *p_out = p;
117 *rhof_out = rhof;
118 *rhog_out = rhog;
119 return 3;
120
121 }
122
123 MSG(" recalc p = %f bar\n",p/1e5);
124 *p_out = p;
125 *rhof_out = rhof;
126 *rhog_out = rhog;
127 return 0;
128 }
129
130 delta_p_old = delta_p;
131
132 if(fabs(delta_p) > 0.4 * p){
133 /* reduce the size of delta_p steps */
134 for(j=0;j<10;++j){
135 delta_p *= 0.5;
136 if(fabs(delta_p) < 0.4 * p)break;
137 }
138 MSG("FPROPS: delta_p was too large, has been shortened\n");
139 }
140
141 p = p + delta_p;
142 }
143
144 if(i==FPROPS_MAX_SUCCSUBS){
145 ERR("too many iterations (%d) in %s(T = %f), or need to try alternative method\n",i,__func__,T);
146 *p_out = p;
147 *rhof_out = rhof;
148 *rhog_out = rhog;
149 return 4;
150 }
151
152 ERR("unknown error with T = %f in %s\n",T,__func__);
153 *p_out = p;
154 *rhof_out = rhof;
155 *rhog_out = rhog;
156 return 9;
157 }
158
159 /**
160 Routine to calculate saturation condition at points close to saturation.
161 The above successive substitution method fails in those cases.
162 */
163 int fprops_sat_T_crit(double T, double *p_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
164 /* nothing yet */
165 return 99;
166 }
167
168 /**
169 Routine to calculate saturation properties for a set value of pressure.
170 @return 0 on success.
171 */
172 int fprops_sat_p(double p, double *T_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
173 double T, T_new, delta_T, delta_T_old, rhof = -1, rhog = -1;
174 int res, i, j;
175 char use_guess = 0;
176
177 MSG("p_c = %f bar\n",d->p_c/1e5);
178
179 /*
180 Estimate of saturation temperature using definition of acentric factor and
181 the assumed p(T) relationship:
182 log10(p)=A + B/T
183 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
184 */
185 T = d->T_c / (1. - 3./7. / (1.+d->omega) * log10(p / d->p_c));
186
187 /* following code is based on successive substitution */
188 MSG("Initial estimate: Tsat(%f bar) = %f K\n",p/1e5,T);
189
190 #if 0
191 int res = fprops_sat_T(T, &T_new, &rhof, &rhog, d);
192 MSG("...gives rhof = %f, rhog = %f and p = %f bar\n",rhof, rhof, p/1e5);
193 #endif
194
195
196 #define FPROPS_MAX_P_SUCCSUBS 40
197
198 /* this approach follows approximately the approach of Soave, page 204:
199 http://dx.doi.org/10.1016/0378-3812(86)90013-0
200 */
201 double Tfactor = 1.002;
202 for(i=0;i<FPROPS_MAX_P_SUCCSUBS;++i){
203 MSG("calculating rho_f estimate\n");
204
205 if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof) || rhof < d->rho_c){
206 T=T/Tfactor;
207 Tfactor = 1 + (Tfactor-1)/1.2;
208 MSG(" adjusting T, error with rho_f. new T = %f\n",T);
209 continue;
210 }
211 MSG("rho_f = %f, T = %f\n",rhof,T);
212
213 MSG("calculating rho_g estimate\n");
214 if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog) || rhog > d->rho_c){
215 T=T*1.001;
216 if(T > d->T_c)T = T/1.001 * 1.00005;
217 if(rhog > d->rho_c)rhog=d->rho_c / 2.;
218 MSG(" adjusting T, error with rho_g = %f, new T = %f\n",rhog,T);
219 continue;
220 }
221 MSG("rho_g = %f\n",rhog);
222
223 MSG("Iter %d: T = %f K, rhof = %f, rhog = %f\n",i, T, rhof, rhog);
224
225 use_guess = 1; /* after first run, start re-using current guess */
226
227 double delta_g = helmholtz_g(T, rhof, d) - helmholtz_g(T, rhog, d);
228 double delta_s = helmholtz_s(T, rhof, d) - helmholtz_s(T, rhog, d);
229 MSG(" delta_g = %f\n",delta_g);
230
231 /* equation from REFPROP */
232 delta_T = delta_g/delta_s;
233
234 MSG(" delta_T = %f K\n",delta_T);
235
236 assert(!isnan(T_new));
237
238 /* convergence test */
239 if(fabs(delta_T/T) < 1e-6){
240 MSG("CONVERGED...\n");
241 /* note possible need for delp non-change test for certain fluids */
242
243 T = T + delta_T;
244
245 /* find liquid phase, using guess */
246 if(
247 fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof)
248 || T > d->T_c || rhof < d->rho_c
249 ){
250 ERR("failed to converge liquid density\n");
251 if(p > 0.9999 * d->p_c){
252 /* no converged, but let's return critical point data */
253 *T_out = d->T_c;
254 *rhof_out = d->rho_c;
255 *rhog_out = d->rho_c;
256 }else{
257 *T_out = T;
258 *rhof_out = rhof;
259 *rhog_out = rhog;
260 }
261 return 3;
262
263 }
264
265 /* find vapour density, using guess */
266 if(
267 fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)
268 || T > d->T_c || rhog > d->rho_c
269 ){
270 ERR("failed to converge vapour density\n");
271 if(p > 0.9999 * d->p_c){
272 /* no converged, but let's return critical point data */
273 *T_out = d->T_c;
274 *rhof_out = d->rho_c;
275 *rhog_out = d->rho_c;
276 }else{
277 *T_out = T;
278 *rhof_out = rhof;
279 *rhog_out = rhog;
280 }
281 return 2;
282 }
283
284 MSG(" recalc T = %f K\n",T);
285 *T_out = T;
286 *rhof_out = rhof;
287 *rhog_out = rhog;
288 return 0;
289 }
290
291 delta_T_old = delta_T;
292
293 if(fabs(delta_T) > 0.5 * T){
294 /* reduce the size of delta_p steps */
295 for(j=0;j<10;++j){
296 delta_T *= 0.25;
297 if(fabs(delta_T) < 0.5 * p)break;
298 }
299 MSG(" delta_T was too large, has been shortened\n");
300 }
301 assert(T > 0);
302 if(T + delta_T < 0)delta_T = -0.2 * T;
303 while(T + delta_T > 10000)delta_T *= 0.5;
304
305 T = T + delta_T;
306 }
307
308 #if 0
309 MSG("Trying to solve using fprops_sat_T");
310
311 /* none of that worked, so use sat_T instead, and iterate. much less efficient */
312 double T1 = 0.999 * d->T_c, T2 = 0.9995 * d->T_c;
313 double p1, p2;
314 res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
315 i = 0;
316 while(1){
317 res = fprops_sat_T(T2, &p1, &rhof, &rhog, d);
318 T = T2;
319 if(fabs(p2 - p) < 1e-6 && !res){
320 MSG("Converged using sat_T");
321 *T_out = T;
322 *rhof_out = rhof;
323 *rhog_out = rhog;
324 return 0;
325 }
326 ++i;
327 /* REFPROP's method */
328 if(fabs(p1 - p1) < 1e-12)break;
329 if(i > 20)break;
330 T = T1 - (T2 - T1) / (p2 - p1) * (p1 - p);
331 if(T > d->T_c && T2 > T1)T = 0.5*(T2 + d->T_c);
332 if(T > d->T_c && T1 > T2)T = 0.5*(T1 + d->T_c);
333 T1 = T2;
334 p1 = p2;
335 T2 = T;
336 }
337 #endif
338
339 MSG("Failed to converge using sat_T approach\n");
340
341 *rhof_out = d->rho_c;
342 *rhog_out = d->rho_c;
343 *T_out = d->T_c;
344 return 99;
345 }
346
347
348
349 /*
350 Iterate to find density as a function of pressure and temperature
351
352 @param use_guess whether (1) or not (0) to use first guess value of *rho.
353 @param rho initial guess (in), and return value (out)
354 @return non-zero on error
355 */
356 int fprops_rho_pT(double p, double T, char phase, char use_guess, const HelmholtzData *d, double *rho){
357
358 double tol = 1e-8;
359 double vlog, Dguess, dpdrho;
360 double vlog_prev;
361 int i;
362
363 #if 0
364 if(p < 1e-11 /*Pa*/){
365 if(T > 0){
366 return /* estimate based on ideal gas equation */
367 }
368 }
369 #endif
370
371 double vc = 1./d->rho_c;
372 if(vc <= 0)vc = 1;
373 double vclog = log(vc);
374
375 char use_liquid_calc = 0;
376
377 if(p < d->p_c * (1. + 7.*(T/d->T_c - 1.)) && T > d->T_c){
378 use_liquid_calc = 0;
379 }else if(p > d->p_c || T > d->T_c){
380 use_liquid_calc = 1;
381 }else if(phase == 'L'){
382 use_liquid_calc = 1;
383 }else{
384 use_liquid_calc = 0;
385 }
386
387 char bad_guess = 0;
388 if(use_guess){
389 if(*rho > 0)vlog = log(1. / *rho);
390 else bad_guess = 1;
391 }
392
393 if(!use_guess || bad_guess){
394 if(p > d->p_c || T > d->T_c){
395 /* supercritical density estimate */
396 vlog = log(0.5 * vc);
397 }else if(use_liquid_calc){
398 /* first guess liquid density */
399 Dguess = fprops_rhof_T_rackett(T, d);
400
401 /* TODO sure first guess is within permitted upper bounds */
402
403 /* if Rackett estimate gives dp/drho < 1, increase guess */
404 int k = 0;
405 while(k++ < 8){
406 double dpdrho = helmholtz_dpdrho_T(T,Dguess,d);
407 if(dpdrho > 0)break;
408 Dguess *= 1.1;
409 };
410 vlog = log(1./Dguess);
411
412 /* upper bound our first guess by log(0.5*v_crit) */
413 double vlog_max = log(0.5 * vc);
414 if(vlog > vlog_max)vlog = vlog_max;
415 }else{
416 /* for gas, first guess is to assume ideal gas */
417 vlog=log(d->R * T / p);
418 }
419
420 if(!use_liquid_calc)vlog=log(d->R*T/p);
421 }
422
423
424 if (use_liquid_calc){
425
426 /* ==== liquid phase iteration ==== */
427
428 for(i=0; i<100; ++i){
429
430 vlog_prev=vlog;
431
432 *rho=1./exp(vlog);
433 double p2 = helmholtz_p_raw(T,*rho,d);
434 double dpdrho = helmholtz_dpdrho_T(T,*rho,d);
435
436 if(dpdrho < 0){
437 // reduce spec volume if we have dp/drho < 0
438 vlog=vlog-0.1;
439 }else{
440 double newt;
441 double dpdlv = -*rho * dpdrho;
442 newt = (p2-p) / dpdlv; // first-order Newton method
443
444 /* loosen the tolerance if we're not getting there */
445 if(i == 10)tol = tol*10;
446 if(i == 15)tol = tol*10;
447
448 if(fabs(newt/p) < 0.0001*tol){
449 // iteration has converged
450 *rho = 1./exp(vlog - newt);
451 return 0;
452 }else{
453 // next guess
454 vlog = vlog_prev - newt;
455
456 // keep step-sizes under control...
457 if(fabs(vlog-vlog_prev) > 0.1 && T < 1.5*d->T_c){
458 vlog = vlog_prev + (vlog-vlog_prev > 0 ? 0.1 : -0.1);
459 }
460 if(vlog > vclog && T < d->T_c){
461 vlog=0.5*(vlog_prev+vclog);
462 }
463 if (vlog < -5.0 && T > d->T_c){
464 use_liquid_calc = 0;
465 MSG(" switching to vapour calc\n");
466 break; /* switch to vapour calc */
467 }
468 }
469 }
470
471 }
472
473 // iteration has not converged
474 *rho = 1./exp(vlog);
475 if (T > d->T_c){
476 use_liquid_calc = 0;
477 }
478
479 if(use_liquid_calc){
480 MSG("liquid iteration did not converge, final value rho = %f\n",*rho);
481 return 1;
482 }
483 }
484
485 /* ==== vapour phase iteration ==== */
486
487 double plog=log(p);
488 for(i=0; i<30; ++i){
489 vlog_prev = vlog;
490 *rho = 1./exp(vlog);
491 double p2 = helmholtz_p_raw(T,*rho,d);
492 double dpdrho = helmholtz_dpdrho_T(T,*rho,d);
493 if(dpdrho < 0. || p2 <= 0.){
494 // increase spec vol if dpdrho < 0 or guess pressure < 0
495 vlog += 0.1;
496 }else{
497 double dpdlv = - *rho * dpdrho;
498 double fvdpl = (log(p2) - plog)*p2/dpdlv; // first-order Newton method
499
500 if(fabs(fvdpl) > 25.)fvdpl=0.1e-5;
501
502 // loosen tol if we're not getting there
503 if(i == 10)tol *= 10;
504 if(i == 15)tol *= 10;
505
506 if(fabs(fvdpl) < 0.0001*tol){
507 // converged, success
508 *rho = 1./exp(vlog);
509 return 0;
510 }else{
511 // next guess
512 vlog=vlog-fvdpl;
513 }
514 }
515 }
516
517 // not converged
518 *rho=1./exp(vlog);
519 MSG(" vapour iteration did not converge\n");
520 return 1;
521 }
522
523

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