/[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 2290 - (show annotations) (download) (as text)
Sun Aug 15 10:11:54 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 12695 byte(s)
Trying to catch errors associated with pressure below triple point for CO2.
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, 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 */
27
28 #include "sat.h"
29 #include "helmholtz_impl.h"
30
31 //#define SAT_DEBUG
32
33 #ifdef SAT_DEBUG
34 # include <assert.h>
35 #else
36 # define assert(ARGS...)
37 #endif
38
39 #include <math.h>
40 #include <stdio.h>
41
42 #define SQ(X) ((X)*(X))
43
44 #define THROW_FPE
45
46 #ifdef THROW_FPE
47 #define _GNU_SOURCE
48 #include <fenv.h>
49 int feenableexcept (int excepts);
50 int fedisableexcept (int excepts);
51 int fegetexcept (void);
52 #endif
53
54 #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 #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
60
61 /**
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 double fprops_psat_T_xiang(double T, const HelmholtzData *d){
68 double p_c = fprops_pc(d);
69 double Zc = p_c / (8314. * d->rho_c * d->T_c);
70
71 #ifdef TEST
72 fprintf(stderr,"Zc = %f\n",Zc);
73 #endif
74
75 double theta = SQ(Zc - 0.29);
76
77 #ifdef TEST
78 fprintf(stderr,"theta = %f\n",theta);
79 #endif
80
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 return p_r * p_c;
109 }
110
111 /**
112 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 double p_c = fprops_pc(d);
118 double p = p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
119 return p;
120 }
121
122
123 /**
124 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
125 see http://dx.doi.org/10.1002/aic.690250412
126 */
127 double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
128 double p_c = fprops_pc(D);
129 double Zc = D->rho_c * D->R * D->T_c / p_c;
130 double Tau = 1. - T/D->T_c;
131 double vf = (D->R * D->T_c / p_c) * pow(Zc, -1 - pow(Tau, 2./7));
132
133 return 1./vf;
134 }
135
136 /**
137 Inverse of fprops_rhof_T_rackett. FIXME this need checking.
138 */
139 double fprops_T_rhof_rackett(double rhof, const HelmholtzData *D){
140 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 double f2 = -log(f1)/log(Zc);
144 return pow(f2 -1, 3./2);
145 }
146
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 double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
152 double Tau = 1. - T/D->T_c;
153 #if 0
154 double p_c = fprops_pc(D);
155 double Zc = D->rho_c * D->R * D->T_c / p_c;
156 # 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 //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 double NNN = PPP + 1./(N1*D->omega + N2);
167 #else
168 /* exact values from Chouaieb for CO2 */
169 # 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 return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
176 }
177
178 /**
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 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 #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 if(T < d->T_t - 1e-8){
199 MSG("Input temperature is below triple-point temperature");
200 return 1;
201 }
202
203 int i = 0;
204 while(i++ < 20){
205 assert(!__isnan(delg));
206 #ifdef SAT_DEBUG
207 fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delf*d->rho_c, delg*d->rho_c);
208 #endif
209 double phirf = helm_resid(tau,delf,d);
210 double phirf_d = helm_resid_del(tau,delf,d);
211 double phirf_dd = helm_resid_deldel(tau,delf,d);
212 double phirg = helm_resid(tau,delg,d);
213 double phirg_d = helm_resid_del(tau,delg,d);
214 double phirg_dd = helm_resid_deldel(tau,delg,d);
215 assert(!__isnan(phirf));
216 assert(!__isnan(phirf_d));
217 assert(!__isnan(phirf_dd));
218 assert(!__isnan(phirg));
219 assert(!__isnan(phirg_d));
220 assert(!__isnan(phirg_dd));
221
222 #define J(FG) (del##FG * (1. + del##FG * phir##FG##_d))
223 #define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG))
224 #define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd)
225 #define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG)
226 double Jf = J(f);
227 double Jg = J(g);
228 double Kf = K(f);
229 double Kg = K(g);
230 double Jf_del = J_del(f);
231 double Jg_del = J_del(g);
232 double Kf_del = K_del(f);
233 double Kg_del = K_del(g);
234
235 double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
236 assert(!__isnan(DELTA));
237
238 #define gamma 1.0
239 delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
240 delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);
241
242 assert(!__isnan(delg));
243 assert(!__isnan(delf));
244
245 if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){
246 //fprintf(stderr,"%s: CONVERGED\n",__func__);
247 *rhof_out = delf * d->rho_c;
248 *rhog_out = delg * d->rho_c;
249 if(__isnan(*rhog_out)){
250 fprintf(stderr,"%s: T = %.12e\n",__func__,T);
251 }
252 *psat_out = helmholtz_p_raw(T, *rhog_out, d);
253 return 0;
254 }
255 if(delg < 0)delg = -0.5*delg;
256 if(delf < 0)delf = -0.5*delf;
257 //if(delg > 1)delg = 1.001;
258 //if(delf > 1)delf = 0.999;
259
260 }
261 *rhof_out = delf * d->rho_c;
262 *rhog_out = delg * d->rho_c;
263 *psat_out = helmholtz_p_raw(T, *rhog_out, d);
264 MSG("NOT CONVERGED for '%s' with T = %e (rhof=%f, rhog=%f)\n",d->name,T,*rhof_out,*rhog_out);
265 return 1;
266
267 }
268
269 /**
270 Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.
271 */
272 double fprops_pc(const HelmholtzData *d){
273 static const HelmholtzData *d_last = NULL;
274 static double p_c = 0;
275 if(d == d_last){
276 return p_c;
277 }
278 p_c = helmholtz_p_raw(d->T_c, d->rho_c,d);
279 d_last = d;
280 return p_c;
281 }
282
283 /**
284 Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.
285 */
286 int fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const HelmholtzData *d){
287 static const HelmholtzData *d_last = NULL;
288 static double p_t, rhof_t, rhog_t;
289 if(d == d_last){
290 *p_t_out = p_t;
291 *rhof_t_out = rhof_t;
292 *rhog_t_out = rhog_t;
293 return 0;
294 }
295 MSG("Calculating saturation for T = %f",d->T_t);
296 int res = fprops_sat_T(d->T_t, &p_t, &rhof_t, &rhog_t,d);
297 if(res)return res;
298 else{
299 d_last = d;
300 *p_t_out = p_t;
301 *rhof_t_out = rhof_t;
302 *rhog_t_out = rhog_t;
303 return 0;
304 }
305 }
306
307
308 /**
309 Solve saturation properties in terms of pressure.
310 This function makes calls to fprops_sat_T, and solves for temperature using
311 a Newton solver algorith. Derivatives dp/dT are calculated using the
312 Clapeyron equation.
313 @return 0 on success.
314 */
315 int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
316 MSG("Calculating for %s at p = %.12e Pa",d->name,p);
317 double T1;
318 double p_c = fprops_pc(d);
319 if(fabs(p - p_c)/p_c < 1e-6){
320 MSG("Very close to critical pressure: using critical temperature without iteration.");
321 T1 = d->T_c;
322 double p1, rhof, rhog;
323 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
324 *Tsat_out = T1;
325 *rhof_out = rhof;
326 *rhog_out = rhog;
327 return res;
328 }else{
329 /*
330 Estimate of saturation temperature using definition of acentric factor and
331 the assumed p(T) relationship:
332 log10(p)=A + B/T
333 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
334 */
335 T1 = d->T_c / (1. - 3./7. / (1.+d->omega) * log10(p / p_c));
336 MSG("Estimated using acentric factor: T = %f",T1);
337 if(T1 < d->T_t){
338 T1 = d->T_t;
339 MSG("Estimate moved up to T_t = %f",T1);
340 }
341 }
342 double p1, rhof, rhog;
343 int i = 0;
344 while(i++ < 50){
345 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
346 if(res){
347 MSG("Got error %d from fprops_sat_T at T = %.12e", res,T1);
348 return 1;
349 }
350 MSG("T1 = %f ������> p = %f bar\trhof = %f\trhog = %f",T1, p1/1e5, rhof, rhog);
351 if(fabs(p1 - p) < 1e-5){
352 *Tsat_out = T1;
353 *rhof_out = rhof;
354 *rhog_out = rhog;
355 return 0;
356 }
357 double hf = helmholtz_h_raw(T1, rhof, d);
358 double hg = helmholtz_h_raw(T1, rhog, d);
359 double dpdT_sat = (hg - hf) / T1 / (1./rhog - 1./rhof);
360 //fprintf(stderr,"\t\tdpdT_sat = %f bar/K\n",dpdT_sat/1e5);
361 double delta_T = -(p1 - p)/dpdT_sat;
362 if(T1 + delta_T < d->T_t - 1e-2){
363 MSG("Correcting sub-triple-point temperature guess");
364 T1 = 0.5 * (d->T_t + T1);
365 }
366 else T1 += delta_T;
367 }
368 MSG("Exceeded iteration limit, returning last guess with error code");
369 *Tsat_out = T1;
370 *rhof_out = rhof;
371 *rhog_out = rhog;
372 return 1;
373 }
374
375
376
377 /**
378 Calculate Tsat based on a value of hf. This value is useful in setting
379 first guess Temperatures when solving for the coordinates (p,h).
380 This function uses the secant method for the iterative solution.
381 */
382 int fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
383 double T1 = 0.4 * d->T_t + 0.6 * d->T_c;
384 double T2 = d->T_t;
385 double h1, h2, p, rhof, rhog;
386 int res = fprops_sat_T(T2, &p, &rhof, &rhog, d);
387 if(res){
388 ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,d->name);
389 return 1;
390 }
391 double tol = 1e-6;
392 h2 = helmholtz_h(T2,rhof,d);
393 if(hf < h2){
394 ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
395 return 2;
396 }
397
398 int i = 0;
399 while(i++ < 60){
400 assert(T1 >= d->T_t - 1e-4);
401 assert(T1 <= d->T_c);
402 res = fprops_sat_T(T1, &p, &rhof, &rhog, d);
403 if(res){
404 ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,d->name);
405 return 1;
406 }
407 h1 = helmholtz_h(T1,rhof, d);
408 if(fabs(h1 - hf) < tol){
409 *Tsat_out = T1;
410 *psat_out = p;
411 *rhof_out = rhof;
412 *rhog_out = rhog;
413 return 0;
414 }
415 if(h1 == h2){
416 MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",d->name,h1,hf);
417 return 2;
418 }
419
420 double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
421 T2 = T1;
422 h2 = h1;
423 while(T1 + delta_T > d->T_c)delta_T *= 0.5;
424 T1 += delta_T;
425 if(T1 < d->T_t)T1 = d->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 return 1;
434 }
435
436
437

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