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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2301 - (hide annotations) (download) (as text)
Sat Aug 21 13:28:35 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 12783 byte(s)
Regen toluene model working, next water.
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 2283 //#define SAT_DEBUG
32    
33     #ifdef SAT_DEBUG
34 jpye 2124 # include <assert.h>
35 jpye 2271 #else
36     # define assert(ARGS...)
37     #endif
38    
39 jpye 2114 #include <math.h>
40 jpye 2262 #include <stdio.h>
41 jpye 2114
42     #define SQ(X) ((X)*(X))
43    
44 jpye 2301 //#define THROW_FPE
45 jpye 2283
46     #ifdef THROW_FPE
47 jpye 2262 #define _GNU_SOURCE
48     #include <fenv.h>
49 jpye 2290 int feenableexcept (int excepts);
50     int fedisableexcept (int excepts);
51     int fegetexcept (void);
52 jpye 2262 #endif
53    
54 jpye 2288 #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 jpye 2290 #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
60 jpye 2288
61 jpye 2122 /**
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 jpye 2114 double fprops_psat_T_xiang(double T, const HelmholtzData *d){
68 jpye 2271 double p_c = fprops_pc(d);
69     double Zc = p_c / (8314. * d->rho_c * d->T_c);
70 jpye 2114
71 jpye 2117 #ifdef TEST
72 jpye 2114 fprintf(stderr,"Zc = %f\n",Zc);
73 jpye 2117 #endif
74 jpye 2114
75     double theta = SQ(Zc - 0.29);
76    
77 jpye 2117 #ifdef TEST
78 jpye 2114 fprintf(stderr,"theta = %f\n",theta);
79 jpye 2117 #endif
80 jpye 2114
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 jpye 2271 return p_r * p_c;
109 jpye 2114 }
110    
111 jpye 2115 /**
112 jpye 2230 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 jpye 2271 double p_c = fprops_pc(d);
118     double p = p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
119 jpye 2230 return p;
120     }
121    
122    
123     /**
124 jpye 2120 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
125     see http://dx.doi.org/10.1002/aic.690250412
126     */
127 jpye 2121 double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
128 jpye 2271 double p_c = fprops_pc(D);
129     double Zc = D->rho_c * D->R * D->T_c / p_c;
130 jpye 2120 double Tau = 1. - T/D->T_c;
131 jpye 2271 double vf = (D->R * D->T_c / p_c) * pow(Zc, -1 - pow(Tau, 2./7));
132 jpye 2120
133     return 1./vf;
134     }
135    
136 jpye 2267 /**
137     Inverse of fprops_rhof_T_rackett. FIXME this need checking.
138     */
139     double fprops_T_rhof_rackett(double rhof, const HelmholtzData *D){
140 jpye 2271 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 jpye 2267 double f2 = -log(f1)/log(Zc);
144     return pow(f2 -1, 3./2);
145     }
146 jpye 2120
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 jpye 2121 double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
152 jpye 2272 double Tau = 1. - T/D->T_c;
153     #if 0
154 jpye 2271 double p_c = fprops_pc(D);
155     double Zc = D->rho_c * D->R * D->T_c / p_c;
156 jpye 2120 # 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 jpye 2122 //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 jpye 2120 double NNN = PPP + 1./(N1*D->omega + N2);
167     #else
168 jpye 2122 /* exact values from Chouaieb for CO2 */
169 jpye 2120 # 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 jpye 2121 return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
176 jpye 2120 }
177    
178 jpye 2283 /**
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 jpye 2262 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 jpye 2283 #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 jpye 2289 if(T < d->T_t - 1e-8){
199 jpye 2292 ERRMSG("Input temperature is below triple-point temperature");
200 jpye 2274 return 1;
201     }
202 jpye 2120
203 jpye 2293 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 jpye 2262 int i = 0;
211 jpye 2270 while(i++ < 20){
212     assert(!__isnan(delg));
213 jpye 2283 #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 jpye 2262 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 jpye 2270 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 jpye 2119
229 jpye 2262 #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 jpye 2119
242 jpye 2262 double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
243 jpye 2270 assert(!__isnan(DELTA));
244 jpye 2122
245 jpye 2270 #define gamma 1.0
246 jpye 2262 delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
247     delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);
248 jpye 2119
249 jpye 2270 assert(!__isnan(delg));
250     assert(!__isnan(delf));
251    
252 jpye 2262 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 jpye 2270 if(__isnan(*rhog_out)){
257     fprintf(stderr,"%s: T = %.12e\n",__func__,T);
258     }
259 jpye 2262 *psat_out = helmholtz_p_raw(T, *rhog_out, d);
260     return 0;
261     }
262 jpye 2283 if(delg < 0)delg = -0.5*delg;
263     if(delf < 0)delf = -0.5*delf;
264 jpye 2124 }
265 jpye 2262 *rhof_out = delf * d->rho_c;
266     *rhog_out = delg * d->rho_c;
267     *psat_out = helmholtz_p_raw(T, *rhog_out, d);
268 jpye 2292 ERRMSG("Not converged: '%s' with T = %e (rhof=%f, rhog=%f).",d->name,T,*rhof_out,*rhog_out);
269 jpye 2262 return 1;
270 jpye 2124
271 jpye 2122 }
272    
273 jpye 2267 /**
274 jpye 2271 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 jpye 2288 MSG("Calculating saturation for T = %f",d->T_t);
300 jpye 2271 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 jpye 2267 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 jpye 2262 int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
320 jpye 2289 MSG("Calculating for %s at p = %.12e Pa",d->name,p);
321 jpye 2271 double T1;
322     double p_c = fprops_pc(d);
323     if(fabs(p - p_c)/p_c < 1e-6){
324 jpye 2289 MSG("Very close to critical pressure: using critical temperature without iteration.");
325 jpye 2271 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 jpye 2289 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 jpye 2271 }
346 jpye 2264 double p1, rhof, rhog;
347     int i = 0;
348 jpye 2271 while(i++ < 50){
349 jpye 2264 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
350 jpye 2289 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 jpye 2264 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 jpye 2289 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 jpye 2264 else T1 += delta_T;
371     }
372 jpye 2289 MSG("Exceeded iteration limit, returning last guess with error code");
373 jpye 2264 *Tsat_out = T1;
374     *rhof_out = rhof;
375     *rhog_out = rhog;
376 jpye 2262 return 1;
377 jpye 2123 }
378    
379    
380 jpye 2124
381 jpye 2267 /**
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 jpye 2290 This function uses the secant method for the iterative solution.
385 jpye 2267 */
386     int fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
387 jpye 2270 double T1 = 0.4 * d->T_t + 0.6 * d->T_c;
388 jpye 2267 double T2 = d->T_t;
389     double h1, h2, p, rhof, rhog;
390     int res = fprops_sat_T(T2, &p, &rhof, &rhog, d);
391 jpye 2269 if(res){
392 jpye 2290 ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,d->name);
393 jpye 2269 return 1;
394     }
395 jpye 2270 double tol = 1e-6;
396 jpye 2267 h2 = helmholtz_h(T2,rhof,d);
397 jpye 2290 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 jpye 2267 int i = 0;
403 jpye 2290 while(i++ < 60){
404     assert(T1 >= d->T_t - 1e-4);
405 jpye 2270 assert(T1 <= d->T_c);
406 jpye 2301 MSG("T1 = %f\n",T1);
407 jpye 2267 res = fprops_sat_T(T1, &p, &rhof, &rhog, d);
408 jpye 2269 if(res){
409 jpye 2290 ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,d->name);
410 jpye 2269 return 1;
411     }
412 jpye 2267 h1 = helmholtz_h(T1,rhof, d);
413 jpye 2270 if(fabs(h1 - hf) < tol){
414 jpye 2267 *Tsat_out = T1;
415     *psat_out = p;
416     *rhof_out = rhof;
417     *rhog_out = rhog;
418     return 0;
419     }
420 jpye 2290 if(h1 == h2){
421     MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",d->name,h1,hf);
422     return 2;
423     }
424 jpye 2124
425 jpye 2267 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
426     T2 = T1;
427     h2 = h1;
428 jpye 2269 while(T1 + delta_T > d->T_c)delta_T *= 0.5;
429 jpye 2267 T1 += delta_T;
430 jpye 2290 if(T1 < d->T_t)T1 = d->T_t;
431 jpye 2270 if(i==20 || i==30)tol*=100;
432 jpye 2267 }
433 jpye 2270 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
434 jpye 2267 *Tsat_out = T1;
435     *psat_out = p;
436     *rhof_out = rhof;
437     *rhog_out = rhog;
438     return 1;
439     }
440 jpye 2123
441 jpye 2124
442 jpye 2271

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