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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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
233
234
235 *psat = d->sat_fn(T,rhof,rhog,d->data,err);
236 }
237
238 /**
239 Calculate the triple point pressure and densities using T_t from the FluidData.
240 */
241 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 }
250
251 if(d->data->T_t == 0){
252 ERRMSG("Note: data for '%s' does not include a valid triple point temperature.",d->name);
253 }
254
255 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 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 MSG("p_t = %f, rhof_t = %f, rhog_t = %f", p_t, rhof_t, rhog_t);
263 }
264
265
266 typedef struct{
267 const PureFluid *P;
268 double logp;
269 FpropsError *err;
270 double Terr;
271 #ifdef SAT_DEBUG
272 int neval;
273 #endif
274 } SatPResidData;
275
276 static ZeroInSubjectFunction sat_p_resid;
277 static double sat_p_resid(double rT, void *user_data){
278 #define D ((SatPResidData *)user_data)
279 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 //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));
285 //if(*(D->err))return -1;
286 #ifdef SAT_DEBUG
287 D->neval++;
288 #endif
289 return resid;
290 #undef D
291 }
292
293
294 /**
295 Solve saturation conditions as a function of pressure.
296
297 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
306 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 use Newton solver to converge, while checking that we remain within
314 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 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 }
326 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 }
333 /* FIXME what about checking triple point pressure? */
334
335
336 SatPResidData D = {
337 P, log(p), err, 0
338 #ifdef SAT_DEBUG
339 ,0
340 #endif
341 };
342 MSG("Solving saturation conditions at p = %f", p);
343 double p1, rT, T, resid;
344 int errn;
345 double Tt = P->data->T_t;
346 if(Tt == 0)Tt = 0.2* P->data->T_c;
347 errn = zeroin_solve(&sat_p_resid, &D, 1./P->data->T_c, 1./Tt, 1e-10, &rT, &resid);
348 if(*err){
349 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 }
353 if(errn){
354 ERRMSG("Failed to solve saturation at p = %f.",p);
355 *err = FPROPS_SAT_CVGC_ERROR;
356 return;
357 }else{
358 if(*err){
359 ERRMSG("Ignoring error inside zeroin_solve iteration at T = %f",D.Terr);
360 }
361 *err = FPROPS_NO_ERROR;
362 }
363 T = 1./rT;
364 fprops_sat_T(T, &p1, rho_f, rho_g, P, err);
365 if(!*err)*T_sat = T;
366 MSG("Got p1 = %f, p = %f in %d iterations", p1, p, D.neval);
367 }
368
369
370 /**
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 This function uses the secant method for the iterative solution.
374 */
375 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 double h1, h2, p, rhof, rhog;
379 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 }
384 double tol = 1e-6;
385 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 if(hf < h2){
390 ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
391 *err = FPROPS_RANGE_ERROR;
392 return;
393 }
394
395 int i = 0;
396 while(i++ < 60){
397 assert(T1 >= P->data->T_t - 1e-4);
398 assert(T1 <= P->data->T_c);
399 MSG("T1 = %f\n",T1);
400 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 }
405 h1 = P->h_fn(T1,rhof, P->data, err);
406 if(*err){
407 ERRMSG("Unable to calculate h");
408 return;
409 }
410 if(fabs(h1 - hf) < tol){
411 *Tsat_out = T1;
412 *psat_out = p;
413 *rhof_out = rhof;
414 *rhog_out = rhog;
415 return; /* SUCCESS */
416 }
417 if(h1 == h2){
418 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",P->name,h1,hf);
419 *err = FPROPS_SAT_CVGC_ERROR;
420 return;
421 }
422
423 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
424 T2 = T1;
425 h2 = h1;
426 while(T1 + delta_T > P->data->T_c)delta_T *= 0.5;
427 T1 += delta_T;
428 if(T1 < P->data->T_t)T1 = P->data->T_t;
429 if(i==20 || i==30)tol*=100;
430 }
431 fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
432 *Tsat_out = T1;
433 *psat_out = p;
434 *rhof_out = rhof;
435 *rhog_out = rhog;
436 *err = FPROPS_SAT_CVGC_ERROR;
437 }
438
439

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