/[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 2680 - (hide annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (9 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 12648 byte(s)
Working on problem with solve_ph. Could be that one of the deriv routines is wrong in the saturation region?
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 jpye 2661 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 jpye 2122 *//** @file
17     Routines to calculate saturation properties using Helmholtz correlation
18     data. We first include some 'generic' saturation equations that make use
19     of the acentric factor and critical point properties to predict saturation
20     properties (pressure, vapour density, liquid density). These correlations
21     seem to be only very rough in some cases, but it is hoped that they will
22     be suitable as first-guess values that can then be fed into an iterative
23     solver to converge to an accurate result.
24 jpye 2114 */
25    
26 jpye 2654 #include "rundata.h"
27 jpye 2114 #include "sat.h"
28 jpye 2654 #include "fprops.h"
29     #include "zeroin.h"
30 jpye 2114
31 jpye 2654 // report lots of stuff
32 jpye 2665 //#define SAT_DEBUG
33 jpye 2654 #define SAT_ERRORS
34 jpye 2283
35 jpye 2654 // assertions for NANs etc
36     //#define SAT_ASSERT
37    
38     #ifdef SAT_ASSERT
39 jpye 2124 # include <assert.h>
40 jpye 2271 #else
41     # define assert(ARGS...)
42     #endif
43    
44 jpye 2114 #include <math.h>
45 jpye 2262 #include <stdio.h>
46 jpye 2114
47     #define SQ(X) ((X)*(X))
48    
49 jpye 2654 #define TCRIT(DATA) (DATA->data->T_c)
50     #define PCRIT(DATA) (DATA->data->p_c)
51     #define RHOCRIT(DATA) (DATA->data->rho_c)
52     #define OMEGA(DATA) (DATA->data->omega)
53     #define RGAS(DATA) (DATA->data->R)
54     #define TTRIP(DATA) (DATA->data->T_t)
55 jpye 2283
56 jpye 2654 #define THROW_FPE
57    
58 jpye 2283 #ifdef THROW_FPE
59 jpye 2262 #define _GNU_SOURCE
60     #include <fenv.h>
61 jpye 2290 int feenableexcept (int excepts);
62     int fedisableexcept (int excepts);
63     int fegetexcept (void);
64 jpye 2262 #endif
65    
66 jpye 2654 # include "color.h"
67    
68 jpye 2288 #ifdef SAT_DEBUG
69 jpye 2654 # define MSG(FMT, ...) \
70     color_on(stderr,ASC_FG_BRIGHTRED);\
71     fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
72     color_on(stderr,ASC_FG_BRIGHTBLUE);\
73     fprintf(stderr,"%s: ",__func__);\
74     color_off(stderr);\
75     fprintf(stderr,FMT "\n",##__VA_ARGS__)
76 jpye 2288 #else
77 jpye 2654 # define MSG(ARGS...) ((void)0)
78 jpye 2288 #endif
79    
80 jpye 2654 #ifdef SAT_ERRORS
81     # define ERRMSG(STR,...) \
82     color_on(stderr,ASC_FG_BRIGHTRED);\
83     fprintf(stderr,"ERROR:");\
84     color_off(stderr);\
85     fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
86     #else
87     # define ERRMSG(ARGS...) ((void)0)
88     #endif
89    
90 jpye 2122 /**
91     Estimate of saturation pressure using H W Xiang ''The new simple extended
92     corresponding-states principle: vapor pressure and second virial
93     coefficient'', Chemical Engineering Science,
94     57 (2002) pp 1439-1449.
95 jpye 2654
96     This seems to be hopeless, or buggy. Tried with water at 373 K, gives
97     525 kPa...
98 jpye 2122 */
99 jpye 2654 double fprops_psat_T_xiang(double T, const FluidData *data){
100     double Zc = data->p_c / (8314. * data->rho_c * data->T_c);
101 jpye 2114
102 jpye 2117 #ifdef TEST
103 jpye 2114 fprintf(stderr,"Zc = %f\n",Zc);
104 jpye 2117 #endif
105 jpye 2114
106     double theta = SQ(Zc - 0.29);
107    
108 jpye 2117 #ifdef TEST
109 jpye 2114 fprintf(stderr,"theta = %f\n",theta);
110 jpye 2117 #endif
111 jpye 2114
112     double aa[] = {5.790206, 4.888195, 33.91196
113     , 6.251894, 15.08591, -315.0248
114     , 11.65859, 46.78273, -1672.179
115     };
116    
117 jpye 2654 double a0 = aa[0] + aa[1]*data->omega + aa[2]*theta;
118     double a1 = aa[3] + aa[4]*data->omega + aa[5]*theta;
119     double a2 = aa[6] + aa[7]*data->omega + aa[8]*theta;
120 jpye 2114
121 jpye 2654 double Tr = T / data->T_c;
122 jpye 2114 double tau = 1 - Tr;
123    
124     double taupow = pow(tau, 1.89);
125     double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow;
126    
127     double logpr = temp * log(Tr);
128     double p_r = exp(logpr);
129    
130     #ifdef TEST
131     fprintf(stderr,"a0 = %f\n", a0);
132     fprintf(stderr,"a1 = %f\n", a1);
133     fprintf(stderr,"a2 = %f\n", a2);
134     fprintf(stderr,"temp = %f\n", temp);
135     fprintf(stderr,"T_r = %f\n", Tr);
136     fprintf(stderr,"p_r = %f\n", p_r);
137     #endif
138    
139 jpye 2654 return p_r * data->p_c;
140 jpye 2114 }
141    
142 jpye 2115 /**
143 jpye 2230 Estimate saturation pressure using acentric factor. This algorithm
144     is used for first estimates for later refinement in the program REFPROP.
145 jpye 2654
146 jpye 2664 This is derived from the definition of the acentric factor,
147     omega = -log10(psat(T1) - 1 where T1/Tc = Tr = 0.7
148    
149     together with the saturation curve obtained if h_fg(T) is assumed constant:
150     ln(psat(T)) = A - B/(T + C)
151    
152     See Sandler 5e, p 320.
153    
154     s Such a curve will pass through (pc,Tc) and (psat(Tr),Tr) where Tr = 0.7,
155     but will probably be inaccurate at Tt. Given additional data, such as the
156     normal boiling point, a better equation like the Antoine equation can be
157     fitted. But we don't use that here.
158    
159 jpye 2230 */
160 jpye 2654 double fprops_psat_T_acentric(double T, const FluidData *data){
161 jpye 2230 /* first guess using acentric factor */
162 jpye 2654 double p = data->p_c * pow(10, -7./3 * (1.+data->omega) * (data->T_c / T - 1.));
163 jpye 2230 return p;
164     }
165    
166    
167     /**
168 jpye 2120 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
169     see http://dx.doi.org/10.1002/aic.690250412
170     */
171 jpye 2654 double fprops_rhof_T_rackett(double T, const FluidData *data){
172     //MSG("RHOCRIT=%f, RGAS=%f, TCRIT=%f, PCRIT=%f",data->rho_c,data->R,data->T_c,data->p_c);
173     double Zc = data->rho_c * data->R * data->T_c / data->p_c;
174     double Tau = 1. - T/data->T_c;
175     //MSG("Zc = %f, Tau = %f",Zc,Tau);
176     double vf = (data->R * data->T_c / data->p_c) * pow(Zc, -1 - pow(Tau, 2./7));
177     //MSG("got vf(T=%f) = %f",T,vf);
178 jpye 2120 return 1./vf;
179     }
180    
181 jpye 2654 /*
182     TODO add Yamada & Gunn sat rhof equation eg from RPP5 eq 4-11.4a, should
183     be more accurate?
184    
185     ^Sean? I think you are referring to http://dx.doi.org/10.1002/aic.690170613,
186     is that correct? That equation requires extra data for V_SC according to
187     the paper. I don't have RPP5, only RPP4 unfortunately -- JP.
188     */
189    
190 jpye 2267 /**
191 jpye 2654 Inverse of fprops_rhof_T_rackett. TODO: this need checking.
192 jpye 2267 */
193 jpye 2654 double fprops_T_rhof_rackett(double rhof, const FluidData *data){
194     double Zc = data->rho_c * data->R * data->T_c / data->p_c;
195     double f1 = data->p_c / data->R / data->T_c / rhof;
196 jpye 2267 double f2 = -log(f1)/log(Zc);
197     return pow(f2 -1, 3./2);
198     }
199 jpye 2120
200 jpye 2654
201 jpye 2120 /**
202     Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
203     see http://dx.doi.org/10.1016/j.tca.2004.05.017
204     */
205 jpye 2654 double fprops_rhog_T_chouaieb(double T, const FluidData *data){
206     double Tau = 1. - T/data->T_c;
207 jpye 2272 #if 0
208 jpye 2654 double Zc = RHOCRIT(d) * RGAS(d) * TCRIT(d) / PCRIT(d);
209 jpye 2120 # define N1 -0.1497547
210 jpye 2649 # define N2 0.6006565
211 jpye 2120 # define P1 -19.348354
212     # define P2 -41.060325
213     # define P3 1.1878726
214     double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */
215 jpye 2122 //MMM = 2.4686277;
216     double PPP = Zc / (P1 + P2*Zc*log(Zc) + P3/Zc);
217     fprintf(stderr,"PPP = %f\n",PPP);
218     //PPP = -0.6240188;
219 jpye 2661 double NNN = PPP + 1./(N1*D->omega + N2);
220 jpye 2120 #else
221 jpye 2122 /* exact values from Chouaieb for CO2 */
222 jpye 2120 # define MMM 2.4686277
223     # define NNN 1.1345838
224     # define PPP -0.6240188
225     #endif
226    
227     double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM));
228 jpye 2654 return data->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
229 jpye 2120 }
230    
231 jpye 2654 void fprops_sat_T(double T, double *psat, double *rhof, double *rhog, const PureFluid *d, FpropsError *err){
232     *psat = d->sat_fn(T,rhof,rhog,d->data,err);
233     }
234    
235 jpye 2283 /**
236 jpye 2654 Calculate the triple point pressure and densities using T_t from the FluidData.
237 jpye 2283 */
238 jpye 2654 void fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const PureFluid *d, FpropsError *err){
239     static const PureFluid *d_last = NULL;
240     static double p_t, rhof_t, rhog_t;
241     if(d == d_last){
242     *p_t_out = p_t;
243     *rhof_t_out = rhof_t;
244     *rhog_t_out = rhog_t;
245     return;
246 jpye 2274 }
247 jpye 2120
248 jpye 2654 if(d->data->T_t == 0){
249     ERRMSG("Note: data for '%s' does not include a valid triple point temperature.",d->name);
250 jpye 2293 }
251    
252 jpye 2664 MSG("Calculating for '%s' (type %d, T_t = %f, T_c = %f, p_c = %f)",d->name, d->type, d->data->T_t, d->data->T_c, d->data->p_c);
253 jpye 2654 fprops_sat_T(d->data->T_t, &p_t, &rhof_t, &rhog_t,d,err);
254     if(*err)return;
255     d_last = d;
256     *p_t_out = p_t;
257     *rhof_t_out = rhof_t;
258     *rhog_t_out = rhog_t;
259 jpye 2664 MSG("p_t = %f, rhof_t = %f, rhog_t = %f", p_t, rhof_t, rhog_t);
260 jpye 2654 }
261 jpye 2119
262 jpye 2122
263 jpye 2654 typedef struct{
264     const PureFluid *P;
265 jpye 2680 double logp;
266 jpye 2654 FpropsError *err;
267 jpye 2662 double Terr;
268 jpye 2680 #ifdef SAT_DEBUG
269     int neval;
270     #endif
271 jpye 2654 } SatPResidData;
272 jpye 2119
273 jpye 2654 static ZeroInSubjectFunction sat_p_resid;
274 jpye 2680 static double sat_p_resid(double rT, void *user_data){
275 jpye 2654 #define D ((SatPResidData *)user_data)
276 jpye 2680 double p, rhof, rhog, resid;
277     fprops_sat_T(1/rT, &p, &rhof, &rhog, D->P, D->err);
278     if(*(D->err))D->Terr = 1/rT;
279     resid = log(p) - D->logp;
280     MSG("T = %f --> p = %f, rhof = %f, rhog = %f, RESID %f", 1/rT, p, rhof, rhog, resid);
281 jpye 2654 //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));
282     //if(*(D->err))return -1;
283 jpye 2680 #ifdef SAT_DEBUG
284     D->neval++;
285     #endif
286     return resid;
287 jpye 2654 #undef D
288     }
289 jpye 2270
290 jpye 2124
291 jpye 2267 /**
292 jpye 2654 Solve saturation conditions as a function of pressure.
293 jpye 2271
294 jpye 2680 Currently this is just a Brent solver. We've tried to improve it slightly
295     by solving for the residual of log(p)-log(p1) as a function of 1/T, which
296     should make the function a bit more linear.
297    
298     TODO Shouldn't(?) be hard at all to improve this to use a Newton solver via
299     the Clapeyron equation?
300    
301     dp/dT = h_fg/(T*v_fg)
302    
303     where
304    
305     h_fg = h(T, rhog) - h(T, rhof)
306     v_fg = 1/rhog - 1/rhof
307     rhof, rhog are evaluated at T.
308    
309     We guess T, calculate saturation conditions at T, then evaluate dp/dT,
310     use Newton solver to converge, while checking that we remain within
311     Tt < T < Tc. It may be faster to iterate using 1/T as the free variable,
312     and log(p) as the residual variable.
313    
314     TODO Better still would be to reexamine the EOS and find a strategy similar to
315     the Akasaka algorithm that works with with rhof, rhog as the free variables,
316     see helmholtz.c...? Method above uses nested iteration on T inside p, so is
317     going to be slooooooow, even if it's fairly reliable.
318     */
319 jpye 2654 void fprops_sat_p(double p, double *T_sat, double *rho_f, double *rho_g, const PureFluid *P, FpropsError *err){
320     if(*err){
321     MSG("ERROR FLAG ALREADY SET");
322 jpye 2271 }
323 jpye 2654 if(p == P->data->p_c){
324     MSG("Requested pressure is critical point pressure, returning CP conditions");
325     *T_sat = P->data->T_c;
326     *rho_f = P->data->rho_c;
327     *rho_g = P->data->rho_c;
328     return;
329 jpye 2271 }
330 jpye 2662 /* FIXME what about checking triple point pressure? */
331    
332    
333 jpye 2680 SatPResidData D = {
334     P, log(p), err, 0
335     #ifdef SAT_DEBUG
336     ,0
337     #endif
338     };
339 jpye 2654 MSG("Solving saturation conditions at p = %f", p);
340 jpye 2680 double p1, rT, T, resid;
341 jpye 2654 int errn;
342     double Tt = P->data->T_t;
343     if(Tt == 0)Tt = 0.2* P->data->T_c;
344 jpye 2680 errn = zeroin_solve(&sat_p_resid, &D, 1./P->data->T_c, 1./Tt, 1e-10, &rT, &resid);
345 jpye 2654 if(*err){
346 jpye 2662 MSG("FPROPS error within zeroin_solve iteration ('%s', p = %f, p_c = %f): %s"
347     , P->name, p, P->data->p_c, fprops_error(*err)
348     );
349 jpye 2271 }
350 jpye 2654 if(errn){
351     ERRMSG("Failed to solve saturation at p = %f.",p);
352     *err = FPROPS_SAT_CVGC_ERROR;
353     return;
354 jpye 2662 }else{
355     if(*err){
356     ERRMSG("Ignoring error inside zeroin_solve iteration at T = %f",D.Terr);
357     }
358     *err = FPROPS_NO_ERROR;
359 jpye 2264 }
360 jpye 2680 T = 1./rT;
361 jpye 2654 fprops_sat_T(T, &p1, rho_f, rho_g, P, err);
362     if(!*err)*T_sat = T;
363 jpye 2680 MSG("Got p1 = %f, p = %f in %d iterations", p1, p, D.neval);
364 jpye 2123 }
365    
366    
367 jpye 2267 /**
368     Calculate Tsat based on a value of hf. This value is useful in setting
369     first guess Temperatures when solving for the coordinates (p,h).
370 jpye 2290 This function uses the secant method for the iterative solution.
371 jpye 2267 */
372 jpye 2654 void fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const PureFluid *P, FpropsError *err){
373     double T1 = 0.4 * P->data->T_t + 0.6 * P->data->T_c;
374     double T2 = P->data->T_t;
375 jpye 2267 double h1, h2, p, rhof, rhog;
376 jpye 2654 fprops_sat_T(T2, &p, &rhof, &rhog, P, err);
377     if(*err){
378     ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,P->name);
379     return;
380 jpye 2269 }
381 jpye 2270 double tol = 1e-6;
382 jpye 2654 h2 = P->h_fn(T2,rhof,P->data, err);
383     if(*err){
384     ERRMSG("Unable to calculate h(T=%f K,rhof=%f kg/m3",T2,rhof);
385     }
386 jpye 2290 if(hf < h2){
387     ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
388 jpye 2654 *err = FPROPS_RANGE_ERROR;
389     return;
390 jpye 2290 }
391    
392 jpye 2267 int i = 0;
393 jpye 2290 while(i++ < 60){
394 jpye 2654 assert(T1 >= P->data->T_t - 1e-4);
395     assert(T1 <= P->data->T_c);
396 jpye 2301 MSG("T1 = %f\n",T1);
397 jpye 2654 fprops_sat_T(T1, &p, &rhof, &rhog, P, err);
398     if(*err){
399     ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,P->name);
400     return;
401 jpye 2269 }
402 jpye 2654 h1 = P->h_fn(T1,rhof, P->data, err);
403     if(*err){
404     ERRMSG("Unable to calculate h");
405     return;
406     }
407 jpye 2270 if(fabs(h1 - hf) < tol){
408 jpye 2267 *Tsat_out = T1;
409     *psat_out = p;
410     *rhof_out = rhof;
411     *rhog_out = rhog;
412 jpye 2654 return; /* SUCCESS */
413 jpye 2267 }
414 jpye 2290 if(h1 == h2){
415 jpye 2654 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",P->name,h1,hf);
416     *err = FPROPS_SAT_CVGC_ERROR;
417     return;
418 jpye 2290 }
419 jpye 2124
420 jpye 2267 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
421     T2 = T1;
422     h2 = h1;
423 jpye 2654 while(T1 + delta_T > P->data->T_c)delta_T *= 0.5;
424 jpye 2267 T1 += delta_T;
425 jpye 2654 if(T1 < P->data->T_t)T1 = P->data->T_t;
426 jpye 2270 if(i==20 || i==30)tol*=100;
427 jpye 2267 }
428 jpye 2270 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
429 jpye 2267 *Tsat_out = T1;
430     *psat_out = p;
431     *rhof_out = rhof;
432     *rhog_out = rhog;
433 jpye 2654 *err = FPROPS_SAT_CVGC_ERROR;
434 jpye 2267 }
435 jpye 2123
436 jpye 2124

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