/[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 2680 - (show annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (9 years, 4 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 /* 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, see <http://www.gnu.org/licenses/>.
16 *//** @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 */
25
26 #include "rundata.h"
27 #include "sat.h"
28 #include "fprops.h"
29 #include "zeroin.h"
30
31 // report lots of stuff
32 //#define SAT_DEBUG
33 #define SAT_ERRORS
34
35 // assertions for NANs etc
36 //#define SAT_ASSERT
37
38 #ifdef SAT_ASSERT
39 # include <assert.h>
40 #else
41 # define assert(ARGS...)
42 #endif
43
44 #include <math.h>
45 #include <stdio.h>
46
47 #define SQ(X) ((X)*(X))
48
49 #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
56 #define THROW_FPE
57
58 #ifdef THROW_FPE
59 #define _GNU_SOURCE
60 #include <fenv.h>
61 int feenableexcept (int excepts);
62 int fedisableexcept (int excepts);
63 int fegetexcept (void);
64 #endif
65
66 # include "color.h"
67
68 #ifdef SAT_DEBUG
69 # 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 #else
77 # define MSG(ARGS...) ((void)0)
78 #endif
79
80 #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 /**
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
96 This seems to be hopeless, or buggy. Tried with water at 373 K, gives
97 525 kPa...
98 */
99 double fprops_psat_T_xiang(double T, const FluidData *data){
100 double Zc = data->p_c / (8314. * data->rho_c * data->T_c);
101
102 #ifdef TEST
103 fprintf(stderr,"Zc = %f\n",Zc);
104 #endif
105
106 double theta = SQ(Zc - 0.29);
107
108 #ifdef TEST
109 fprintf(stderr,"theta = %f\n",theta);
110 #endif
111
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 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
121 double Tr = T / data->T_c;
122 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 return p_r * data->p_c;
140 }
141
142 /**
143 Estimate saturation pressure using acentric factor. This algorithm
144 is used for first estimates for later refinement in the program REFPROP.
145
146 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 */
160 double fprops_psat_T_acentric(double T, const FluidData *data){
161 /* first guess using acentric factor */
162 double p = data->p_c * pow(10, -7./3 * (1.+data->omega) * (data->T_c / T - 1.));
163 return p;
164 }
165
166
167 /**
168 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
169 see http://dx.doi.org/10.1002/aic.690250412
170 */
171 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 return 1./vf;
179 }
180
181 /*
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 /**
191 Inverse of fprops_rhof_T_rackett. TODO: this need checking.
192 */
193 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 double f2 = -log(f1)/log(Zc);
197 return pow(f2 -1, 3./2);
198 }
199
200
201 /**
202 Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
203 see http://dx.doi.org/10.1016/j.tca.2004.05.017
204 */
205 double fprops_rhog_T_chouaieb(double T, const FluidData *data){
206 double Tau = 1. - T/data->T_c;
207 #if 0
208 double Zc = RHOCRIT(d) * RGAS(d) * TCRIT(d) / PCRIT(d);
209 # define N1 -0.1497547
210 # define N2 0.6006565
211 # 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 //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 double NNN = PPP + 1./(N1*D->omega + N2);
220 #else
221 /* exact values from Chouaieb for CO2 */
222 # 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 return data->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
229 }
230
231 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 /**
236 Calculate the triple point pressure and densities using T_t from the FluidData.
237 */
238 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 }
247
248 if(d->data->T_t == 0){
249 ERRMSG("Note: data for '%s' does not include a valid triple point temperature.",d->name);
250 }
251
252 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 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 MSG("p_t = %f, rhof_t = %f, rhog_t = %f", p_t, rhof_t, rhog_t);
260 }
261
262
263 typedef struct{
264 const PureFluid *P;
265 double logp;
266 FpropsError *err;
267 double Terr;
268 #ifdef SAT_DEBUG
269 int neval;
270 #endif
271 } SatPResidData;
272
273 static ZeroInSubjectFunction sat_p_resid;
274 static double sat_p_resid(double rT, void *user_data){
275 #define D ((SatPResidData *)user_data)
276 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 //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));
282 //if(*(D->err))return -1;
283 #ifdef SAT_DEBUG
284 D->neval++;
285 #endif
286 return resid;
287 #undef D
288 }
289
290
291 /**
292 Solve saturation conditions as a function of pressure.
293
294 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 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 }
323 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 }
330 /* FIXME what about checking triple point pressure? */
331
332
333 SatPResidData D = {
334 P, log(p), err, 0
335 #ifdef SAT_DEBUG
336 ,0
337 #endif
338 };
339 MSG("Solving saturation conditions at p = %f", p);
340 double p1, rT, T, resid;
341 int errn;
342 double Tt = P->data->T_t;
343 if(Tt == 0)Tt = 0.2* P->data->T_c;
344 errn = zeroin_solve(&sat_p_resid, &D, 1./P->data->T_c, 1./Tt, 1e-10, &rT, &resid);
345 if(*err){
346 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 }
350 if(errn){
351 ERRMSG("Failed to solve saturation at p = %f.",p);
352 *err = FPROPS_SAT_CVGC_ERROR;
353 return;
354 }else{
355 if(*err){
356 ERRMSG("Ignoring error inside zeroin_solve iteration at T = %f",D.Terr);
357 }
358 *err = FPROPS_NO_ERROR;
359 }
360 T = 1./rT;
361 fprops_sat_T(T, &p1, rho_f, rho_g, P, err);
362 if(!*err)*T_sat = T;
363 MSG("Got p1 = %f, p = %f in %d iterations", p1, p, D.neval);
364 }
365
366
367 /**
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 This function uses the secant method for the iterative solution.
371 */
372 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 double h1, h2, p, rhof, rhog;
376 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 }
381 double tol = 1e-6;
382 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 if(hf < h2){
387 ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
388 *err = FPROPS_RANGE_ERROR;
389 return;
390 }
391
392 int i = 0;
393 while(i++ < 60){
394 assert(T1 >= P->data->T_t - 1e-4);
395 assert(T1 <= P->data->T_c);
396 MSG("T1 = %f\n",T1);
397 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 }
402 h1 = P->h_fn(T1,rhof, P->data, err);
403 if(*err){
404 ERRMSG("Unable to calculate h");
405 return;
406 }
407 if(fabs(h1 - hf) < tol){
408 *Tsat_out = T1;
409 *psat_out = p;
410 *rhof_out = rhof;
411 *rhog_out = rhog;
412 return; /* SUCCESS */
413 }
414 if(h1 == h2){
415 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",P->name,h1,hf);
416 *err = FPROPS_SAT_CVGC_ERROR;
417 return;
418 }
419
420 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
421 T2 = T1;
422 h2 = h1;
423 while(T1 + delta_T > P->data->T_c)delta_T *= 0.5;
424 T1 += delta_T;
425 if(T1 < P->data->T_t)T1 = P->data->T_t;
426 if(i==20 || i==30)tol*=100;
427 }
428 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
429 *Tsat_out = T1;
430 *psat_out = p;
431 *rhof_out = rhof;
432 *rhog_out = rhog;
433 *err = FPROPS_SAT_CVGC_ERROR;
434 }
435
436

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