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

Annotation of /branches/sid/models/johnpye/fprops/sat.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3024 - (hide annotations) (download) (as text)
Thu Jul 23 04:16:35 2015 UTC (3 years, 10 months ago) by sid
File MIME type: text/x-csrc
File size: 12645 byte(s)
first attempt at two phase

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 sid 3024
149 jpye 2664 together with the saturation curve obtained if h_fg(T) is assumed constant:
150     ln(psat(T)) = A - B/(T + C)
151    
152 sid 3024 See Sandler 5e, p 320.
153 jpye 2664
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 sid 3024
233    
234    
235 jpye 2654 *psat = d->sat_fn(T,rhof,rhog,d->data,err);
236     }
237    
238 jpye 2283 /**
239 jpye 2654 Calculate the triple point pressure and densities using T_t from the FluidData.
240 jpye 2283 */
241 jpye 2654 void fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const PureFluid *d, FpropsError *err){
242     static const PureFluid *d_last = NULL;
243     static double p_t, rhof_t, rhog_t;
244     if(d == d_last){
245     *p_t_out = p_t;
246     *rhof_t_out = rhof_t;
247     *rhog_t_out = rhog_t;
248     return;
249 jpye 2274 }
250 jpye 2120
251 jpye 2654 if(d->data->T_t == 0){
252     ERRMSG("Note: data for '%s' does not include a valid triple point temperature.",d->name);
253 jpye 2293 }
254    
255 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);
256 jpye 2654 fprops_sat_T(d->data->T_t, &p_t, &rhof_t, &rhog_t,d,err);
257     if(*err)return;
258     d_last = d;
259     *p_t_out = p_t;
260     *rhof_t_out = rhof_t;
261     *rhog_t_out = rhog_t;
262 jpye 2664 MSG("p_t = %f, rhof_t = %f, rhog_t = %f", p_t, rhof_t, rhog_t);
263 jpye 2654 }
264 jpye 2119
265 jpye 2122
266 jpye 2654 typedef struct{
267     const PureFluid *P;
268 jpye 2680 double logp;
269 jpye 2654 FpropsError *err;
270 jpye 2662 double Terr;
271 jpye 2680 #ifdef SAT_DEBUG
272     int neval;
273     #endif
274 jpye 2654 } SatPResidData;
275 jpye 2119
276 jpye 2654 static ZeroInSubjectFunction sat_p_resid;
277 jpye 2680 static double sat_p_resid(double rT, void *user_data){
278 jpye 2654 #define D ((SatPResidData *)user_data)
279 jpye 2680 double p, rhof, rhog, resid;
280     fprops_sat_T(1/rT, &p, &rhof, &rhog, D->P, D->err);
281     if(*(D->err))D->Terr = 1/rT;
282     resid = log(p) - D->logp;
283     MSG("T = %f --> p = %f, rhof = %f, rhog = %f, RESID %f", 1/rT, p, rhof, rhog, resid);
284 jpye 2654 //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));
285     //if(*(D->err))return -1;
286 jpye 2680 #ifdef SAT_DEBUG
287     D->neval++;
288     #endif
289     return resid;
290 jpye 2654 #undef D
291     }
292 jpye 2270
293 jpye 2124
294 jpye 2267 /**
295 sid 3024 Solve saturation conditions as a function of pressure.
296 jpye 2271
297 jpye 2680 Currently this is just a Brent solver. We've tried to improve it slightly
298     by solving for the residual of log(p)-log(p1) as a function of 1/T, which
299     should make the function a bit more linear.
300    
301     TODO Shouldn't(?) be hard at all to improve this to use a Newton solver via
302     the Clapeyron equation?
303    
304     dp/dT = h_fg/(T*v_fg)
305 sid 3024
306 jpye 2680 where
307    
308     h_fg = h(T, rhog) - h(T, rhof)
309     v_fg = 1/rhog - 1/rhof
310     rhof, rhog are evaluated at T.
311    
312     We guess T, calculate saturation conditions at T, then evaluate dp/dT,
313 sid 3024 use Newton solver to converge, while checking that we remain within
314 jpye 2680 Tt < T < Tc. It may be faster to iterate using 1/T as the free variable,
315     and log(p) as the residual variable.
316    
317     TODO Better still would be to reexamine the EOS and find a strategy similar to
318     the Akasaka algorithm that works with with rhof, rhog as the free variables,
319     see helmholtz.c...? Method above uses nested iteration on T inside p, so is
320     going to be slooooooow, even if it's fairly reliable.
321     */
322 jpye 2654 void fprops_sat_p(double p, double *T_sat, double *rho_f, double *rho_g, const PureFluid *P, FpropsError *err){
323     if(*err){
324     MSG("ERROR FLAG ALREADY SET");
325 jpye 2271 }
326 jpye 2654 if(p == P->data->p_c){
327     MSG("Requested pressure is critical point pressure, returning CP conditions");
328     *T_sat = P->data->T_c;
329     *rho_f = P->data->rho_c;
330     *rho_g = P->data->rho_c;
331     return;
332 jpye 2271 }
333 jpye 2662 /* FIXME what about checking triple point pressure? */
334    
335 sid 3024
336 jpye 2680 SatPResidData D = {
337     P, log(p), err, 0
338     #ifdef SAT_DEBUG
339     ,0
340     #endif
341     };
342 jpye 2654 MSG("Solving saturation conditions at p = %f", p);
343 jpye 2680 double p1, rT, T, resid;
344 jpye 2654 int errn;
345     double Tt = P->data->T_t;
346     if(Tt == 0)Tt = 0.2* P->data->T_c;
347 jpye 2680 errn = zeroin_solve(&sat_p_resid, &D, 1./P->data->T_c, 1./Tt, 1e-10, &rT, &resid);
348 jpye 2654 if(*err){
349 jpye 2662 MSG("FPROPS error within zeroin_solve iteration ('%s', p = %f, p_c = %f): %s"
350     , P->name, p, P->data->p_c, fprops_error(*err)
351     );
352 jpye 2271 }
353 jpye 2654 if(errn){
354     ERRMSG("Failed to solve saturation at p = %f.",p);
355     *err = FPROPS_SAT_CVGC_ERROR;
356     return;
357 jpye 2662 }else{
358     if(*err){
359     ERRMSG("Ignoring error inside zeroin_solve iteration at T = %f",D.Terr);
360     }
361     *err = FPROPS_NO_ERROR;
362 jpye 2264 }
363 jpye 2680 T = 1./rT;
364 jpye 2654 fprops_sat_T(T, &p1, rho_f, rho_g, P, err);
365     if(!*err)*T_sat = T;
366 jpye 2680 MSG("Got p1 = %f, p = %f in %d iterations", p1, p, D.neval);
367 jpye 2123 }
368    
369    
370 jpye 2267 /**
371     Calculate Tsat based on a value of hf. This value is useful in setting
372     first guess Temperatures when solving for the coordinates (p,h).
373 jpye 2290 This function uses the secant method for the iterative solution.
374 jpye 2267 */
375 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){
376     double T1 = 0.4 * P->data->T_t + 0.6 * P->data->T_c;
377     double T2 = P->data->T_t;
378 jpye 2267 double h1, h2, p, rhof, rhog;
379 jpye 2654 fprops_sat_T(T2, &p, &rhof, &rhog, P, err);
380     if(*err){
381     ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,P->name);
382     return;
383 jpye 2269 }
384 jpye 2270 double tol = 1e-6;
385 jpye 2654 h2 = P->h_fn(T2,rhof,P->data, err);
386     if(*err){
387     ERRMSG("Unable to calculate h(T=%f K,rhof=%f kg/m3",T2,rhof);
388     }
389 jpye 2290 if(hf < h2){
390     ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
391 jpye 2654 *err = FPROPS_RANGE_ERROR;
392     return;
393 jpye 2290 }
394    
395 jpye 2267 int i = 0;
396 jpye 2290 while(i++ < 60){
397 jpye 2654 assert(T1 >= P->data->T_t - 1e-4);
398     assert(T1 <= P->data->T_c);
399 jpye 2301 MSG("T1 = %f\n",T1);
400 jpye 2654 fprops_sat_T(T1, &p, &rhof, &rhog, P, err);
401     if(*err){
402     ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,P->name);
403     return;
404 jpye 2269 }
405 jpye 2654 h1 = P->h_fn(T1,rhof, P->data, err);
406     if(*err){
407     ERRMSG("Unable to calculate h");
408     return;
409     }
410 jpye 2270 if(fabs(h1 - hf) < tol){
411 jpye 2267 *Tsat_out = T1;
412     *psat_out = p;
413     *rhof_out = rhof;
414     *rhog_out = rhog;
415 jpye 2654 return; /* SUCCESS */
416 jpye 2267 }
417 jpye 2290 if(h1 == h2){
418 jpye 2654 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",P->name,h1,hf);
419     *err = FPROPS_SAT_CVGC_ERROR;
420     return;
421 jpye 2290 }
422 jpye 2124
423 jpye 2267 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
424     T2 = T1;
425     h2 = h1;
426 jpye 2654 while(T1 + delta_T > P->data->T_c)delta_T *= 0.5;
427 jpye 2267 T1 += delta_T;
428 jpye 2654 if(T1 < P->data->T_t)T1 = P->data->T_t;
429 jpye 2270 if(i==20 || i==30)tol*=100;
430 jpye 2267 }
431 jpye 2270 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
432 jpye 2267 *Tsat_out = T1;
433     *psat_out = p;
434     *rhof_out = rhof;
435     *rhog_out = rhog;
436 jpye 2654 *err = FPROPS_SAT_CVGC_ERROR;
437 jpye 2267 }
438 jpye 2123
439 jpye 2124

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