/[ascend]/trunk/models/johnpye/fprops/pengrob.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/pengrob.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations) (download) (as text)
Mon Jan 21 06:54:33 2013 UTC (7 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 24196 byte(s)
comment out debug file creation
1 /* ASCEND modelling environment
2 Copyright (C) 2011-2012 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 Implementation of the Peng-Robinson equation of state.
18
19 For nomenclature see Sandler, 5e and IAPWS-95 and IAPWS Derivatives.
20
21 Authors: John Pye, Richard Towers, Sean Muratet
22 */
23
24 #include "pengrob.h"
25 #include "rundata.h"
26 #include "cubicroots.h"
27 #include "fprops.h"
28 #include "sat.h"
29 #include "ideal_impl.h"
30 #include "cp0.h"
31 #include "zeroin.h"
32 #include <math.h>
33 #include <stdio.h>
34 #include <assert.h>
35
36 #include "helmholtz.h" // for helmholtz_prepare
37
38 /* these are the 'raw' functions, they don't do phase equilibrium. */
39 PropEvalFn pengrob_p;
40 PropEvalFn pengrob_u;
41 PropEvalFn pengrob_h;
42 PropEvalFn pengrob_s;
43 PropEvalFn pengrob_a;
44 PropEvalFn pengrob_g;
45 PropEvalFn pengrob_cp;
46 PropEvalFn pengrob_cv;
47 PropEvalFn pengrob_w;
48 PropEvalFn pengrob_dpdrho_T;
49 PropEvalFn pengrob_alphap;
50 PropEvalFn pengrob_betap;
51 SatEvalFn pengrob_sat;
52 //SatEvalFn pengrob_sat_akasaka;
53
54 static double MidpointPressureCubic(double T, const FluidData *data, FpropsError *err);
55
56 //#define PR_DEBUG
57 #define PR_ERRORS
58
59 #ifdef PR_DEBUG
60 # include "color.h"
61 # define MSG(FMT, ...) \
62 color_on(stderr,ASC_FG_BRIGHTRED);\
63 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
64 color_on(stderr,ASC_FG_BRIGHTBLUE);\
65 fprintf(stderr,"%s: ",__func__);\
66 color_off(stderr);\
67 fprintf(stderr,FMT "\n",##__VA_ARGS__)
68 #else
69 # define MSG(ARGS...) ((void)0)
70 #endif
71
72 /* TODO centralise declaration of our error-reporting function somehow...? */
73 #ifdef PR_ERRORS
74 # include "color.h"
75 # define ERRMSG(STR,...) \
76 color_on(stderr,ASC_FG_BRIGHTRED);\
77 fprintf(stderr,"ERROR:");\
78 color_off(stderr);\
79 fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
80 #else
81 # define ERRMSG(ARGS...) ((void)0)
82 #endif
83
84
85 PureFluid *pengrob_prepare(const EosData *E, const ReferenceState *ref){
86 MSG("Preparing PR fluid...");
87 PureFluid *P = FPROPS_NEW(PureFluid);
88 P->data = FPROPS_NEW(FluidData);
89
90 /* metadata */
91 P->name = E->name;
92 P->source = E->source;
93 P->type = FPROPS_PENGROB;
94
95 #define D P->data
96 /* common data across all correlation types */
97 switch(E->type){
98 case FPROPS_HELMHOLTZ:
99 #define I E->data.helm
100 D->M = I->M;
101 D->R = I->R;
102 D->T_t = I->T_t;
103 D->T_c = I->T_c;
104 D->rho_c = I->rho_c;
105 D->omega = I->omega;
106 D->cp0 = cp0_prepare(I->ideal, D->R, D->T_c);
107
108 // helmholtz data doesn't include p_c so we need to calculate it somehow...
109
110 #ifdef USE_HELMHOLTZ_RHO_C
111 // use Zc = 0.307 to calculate p_c from the T_c and rho_c in the
112 // helmholtz data. This seems to be a bad choice
113 // since eg for CO2 it gives a significantly overestimated p_c.
114 {
115 double Zc = 0.307;
116 D->p_c = Zc * D->R * D->T_c * D->rho_c;
117 MSG("p_c = %f", D->p_c);
118 }
119 #else // USE HELMHOLTZ P_C
120 // Probably the preferred alternative in this case is to
121 // use rho_c and T_c to calculated helmholtz_p(T_c, rho_c), then
122 // use that pressure as the PR p_c, together with updating the
123 // value of rho_c for consistency with PR EOS (from the known
124 // Z_c = 0.307.
125 {
126 FpropsError herr = FPROPS_NO_ERROR;
127 MSG("Preparing helmholtz data '%s'...",E->name);
128 PureFluid *PH = helmholtz_prepare(E,ref);
129 if(!PH){
130 ERRMSG("Failed to create Helmholtz runtime data");
131 return NULL;
132 }
133 D->p_c = PH->p_fn(D->T_c, D->rho_c, PH->data, &herr);
134 MSG("Calculated p_c = %f from Helmholtz data",D->p_c);
135 if(herr){
136 ERRMSG("Failed to calculate critical pressure (%s)",fprops_error(herr));
137 return NULL;
138 }
139 double Zc = 0.307;
140 D->rho_c = D->p_c / (Zc * D->R * D->T_c);
141 }
142 #endif
143 break;
144 #undef I
145 case FPROPS_CUBIC:
146 #define I E->data.cubic
147 D->M = I->M;
148 D->R = R_UNIVERSAL / I->M;
149 D->T_t = I->T_t;
150 D->T_c = I->T_c;
151 D->p_c = I->p_c;
152 D->rho_c = I->rho_c;
153 D->omega = I->omega;
154 D->cp0 = cp0_prepare(I->ideal, D->R, D->T_c);
155 break;
156 default:
157 fprintf(stderr,"Invalid EOS data\n");
158 return NULL;
159 }
160
161 if(D->p_c == 0){
162 ERRMSG("ERROR p_c is zero in this data, need to calculate it here somehow");
163 return NULL;
164 }
165
166 /* FIXME note that in the following paper, the constants in the PR EOS
167 are given with lots more decimal places. Need to figure out if it's
168 preferable to use those extra DPs, or not...?
169 http://www.che.uah.edu/courseware/che641/peng-robinson-derivatives-joule-thompson.pdf
170 */
171
172 /* NOTE: we're using a mass basis for all our property calculations. That
173 means that our 'b' is the usual 'b/M' since our 'R' is 'Rm/M'. Because of
174 this, our p is still OK, since Ru/M / (Vm/M - b/M) is still the same value.*/
175 #define C P->data->corr.pengrob
176 C = FPROPS_NEW(PengrobRunData);
177 C->aTc = 0.45724 * SQ(D->R * D->T_c) / D->p_c;
178 C->b = 0.07780 * D->R * D->T_c / D->p_c;
179 C->kappa = 0.37464 + (1.54226 - 0.26992 * D->omega) * D->omega;
180
181 /* function pointers... more to come still? */
182 #define FN(VAR) P->VAR##_fn = &pengrob_##VAR
183 FN(p); FN(u); FN(h); FN(s); FN(a); FN(g); FN(cp); FN(cv); FN(w);
184 FN(dpdrho_T); FN(alphap); FN(betap);
185 FN(sat);
186 #undef FN
187 #undef I
188 #undef D
189 #undef C
190 //P->sat_fn = &pengrob_sat_akasaka;
191
192 return P;
193 }
194
195
196 /* shortcuts take us straight into the correct data structure */
197 #define PD data->corr.pengrob
198 #define PD_TCRIT data->T_c
199 #define PD_RHOCRIT data->rho_c
200 #define PD_M data->M
201 #define SQRT2 1.4142135623730951
202
203 #define DEFINE_SQRTALPHA \
204 double sqrtalpha = 1 + PD->kappa * (1 - sqrt(T / PD_TCRIT));
205 #define DEFINE_A \
206 double a = PD->aTc * SQ(sqrtalpha);
207
208 #define DEFINE_V double v = 1. / rho;
209
210 /**
211 Maxima code:
212 a(T) := aTc * (1 + kappa * (1 - sqrt(T/Tc)));
213 diff(a(T),T,1);
214 XXX
215 */
216 #define DEFINE_DADT \
217 double dadT = -PD->kappa * PD->aTc * sqrtalpha / sqrt(T * PD_TCRIT)
218
219 /**
220 Maxima code:
221 a(T) := aTc * (1 + kappa * (1 - sqrt(T/Tc)));
222 diff(a(T),T,2);
223 XXX
224 */
225 #define DEFINE_D2ADT2 \
226 double d2adt2 = PD->aTc * PD->kappa * sqrt(PD_TCRIT/T) * (1 + PD->kappa) / (2 * T * PD_TCRIT);
227
228 #define DEFINE_DPDT_RHO \
229 double dpdT_rho = data->R/(v - PD->b) - dadT/(v*(v + PD->b) + PD->b*(v - PD->b))
230
231 double pengrob_p(double T, double rho, const FluidData *data, FpropsError *err){
232 DEFINE_SQRTALPHA;
233 DEFINE_A;
234 DEFINE_V;
235 double b = PD->b;
236 if(rho > 1./b){
237 /* TODO check this: is it because it gives p < 0? */
238 MSG("Density exceeds limit value 1/b = %f",1./b);
239 *err = FPROPS_RANGE_ERROR;
240 }
241 //MSG("v = %f, b = %f",v,b);
242 double p = (data->R * T)/(v - b) - a/(v*(v + b) + b*(v - b));
243 //MSG("p(T = %f, rho = %f) = %f",T,rho,p);
244 return p;
245 }
246
247
248 double pengrob_h(double T, double rho, const FluidData *data, FpropsError *err){
249 DEFINE_SQRTALPHA;
250 DEFINE_A;
251 DEFINE_V;
252 if(rho > 1./PD->b){
253 MSG("Density exceeds limit value 1/b = %f",1./PD->b);
254 *err = FPROPS_RANGE_ERROR;
255 return 0;
256 }
257 double h0 = ideal_h(T,rho,data,err);
258 double p = pengrob_p(T, rho, data, err);
259 double Z = p * v / (data->R * T);
260 double B = p * PD->b / (data->R * T);
261 DEFINE_DADT;
262 double hr = data->R * T * (Z - 1) + (T*dadT - a)/(2*SQRT2 * PD->b) * log((Z + (1+SQRT2)*B) / (Z + (1-SQRT2)*B));
263 return h0 + hr;
264 }
265
266
267 double pengrob_s(double T, double rho, const FluidData *data, FpropsError *err){
268 DEFINE_SQRTALPHA;
269 DEFINE_V;
270 double b = PD->b;
271 if(rho > 1./b){
272 MSG("Density exceeds limit value 1/b = %f",1./b);
273 *err = FPROPS_RANGE_ERROR;
274 return 0;
275 }
276 double s0 = ideal_s(T,rho,data,err);
277 double p = pengrob_p(T, rho, data, err);
278 double Z = p * v / (data->R * T);
279 double B = p * b / (data->R * T);
280 DEFINE_DADT;
281 //MSG("s0 = %f, p = %f, Z = %f, B = %f, dadt = %f", s0, p, Z, B, dadt);
282 //MSG("log(Z-B) = %f, log((Z+(1+SQRT2)*B)/(Z+(1-SQRT2)*B = %f", log(Z-B), log((Z+(1+SQRT2)*B)/(Z+(1-SQRT2)*B)));
283 double sr = data->R * log(Z-B) + dadT/(2*SQRT2*b) * log((Z+(1+SQRT2)*B)/(Z+(1-SQRT2)*B));
284 return s0 + sr;
285 }
286
287
288 double pengrob_a(double T, double rho, const FluidData *data, FpropsError *err){
289 // FIXME maybe we can improve this with more direct maths
290 double h = pengrob_h(T,rho,data,err);
291 double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!
292 double p = pengrob_p(T,rho,data,err); // duplicated calculation of p!
293 MSG("h = %f, p = %f, s = %f, rho = %f, T = %f",h,p,s,rho,T);
294 return (h - p/rho) - T * s;
295 // previous code from Richard, probably fine but need to check
296 // DEFINE_A;
297 // double vm = PD_M / rho;
298 // double p = pengrob_p(T, rho, data, err);
299 // return -log(p*(v - PD->b)/(data->R * T))+a/(2*SQRT2 * PD->b * data->R *T)*log((v + (1-SQRT2)*PD->b)/(v + (1+SQRT2)*PD->b));
300 }
301
302
303 double pengrob_u(double T, double rho, const FluidData *data, FpropsError *err){
304 // FIXME work out a cleaner approach to this...
305 double p = pengrob_p(T, rho, data, err);
306 return pengrob_h(T,rho,data,err) - p/rho; // duplicated calculation of p!
307 }
308
309 double pengrob_g(double T, double rho, const FluidData *data, FpropsError *err){
310 if(rho > 1./PD->b){
311 MSG("Density exceeds limit value 1/b = %f",1./PD->b);
312 *err = FPROPS_RANGE_ERROR;
313 }
314 #if 0
315 double h = pengrob_h(T,rho,data,err);
316 double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!
317 if(isnan(h))MSG("h is nan");
318 if(isnan(s))MSG("s is nan");
319 return h - T*s;
320 #else
321 // previous code from Richard, probably fine but need to check
322 DEFINE_SQRTALPHA;
323 DEFINE_A;
324 DEFINE_V;
325 double p = pengrob_p(T, rho, data, err);
326 double Z = p*v/(data->R * T);
327 double B = p*PD->b/(data->R * T);
328 double A = p * a / SQ(data->R * T);
329 return log(fabs(Z-B))-(A/(sqrt(8)*B))*log(fabs((Z+(1+sqrt(2))*B)/(Z+(1-sqrt(2))*B)))+Z-1;
330 #endif
331 }
332
333 /**
334 \f[ c_p = \left(\frac{\partial u}{\partial T} \right)_v \f]
335
336 See Pratt, 2001. "Thermodynamic Properties Involving Derivatives",
337 ASEE Chemical Engineering Division Newsletter, Malaysia
338 http://www.che.uah.edu/courseware/che641/peng-robinson-derivatives-joule-thompson.pdf
339
340 Maxima code:
341 u(T,x):= (T*diff(a(T),T) - a(T))/b/sqrt(8) * log((Z(T) + B(T)*(1+sqrt(2)))/(Z(T)+B(T)*(1-sqrt(2))));
342 diff(u(T,v),T);
343 subst(B(T)*(alpha_p - 1/T), diff(B(T),T), %);
344 subst(Z(T)*(alpha_p - 1/T), diff(Z(T),T), %);
345 subst(B(T)*v/Z(T),b,%);
346 subst(b*Z(T)/B(T),v,%);
347 ratsimp(%);
348 */
349 double pengrob_cv(double T, double rho, const FluidData *data, FpropsError *err){
350 DEFINE_V;
351 DEFINE_D2ADT2;
352 double cv0 = ideal_cv(T, rho, data, err);
353 MSG("cv0 = %f",cv0);
354 #define DEFINE_CVR \
355 double p = pengrob_p(T, rho, data, err); \
356 double Z = p * v / (data->R * T); \
357 double B = p * PD->b / (data->R * T); \
358 double cvr1 = T * d2adt2/ (PD->b * 2*SQRT2); \
359 double cvr2 = (Z+B*(1+SQRT2))/(Z+B*(1-SQRT2)); \
360 double cvr = cvr1*log(cvr2);
361 DEFINE_CVR;
362 MSG("d2adT2 = %f", d2adt2);
363 MSG("b = %f",PD->b);
364 MSG("cvr1 = %f, cvr2 = %f, log(cvr2) = %f", cvr1, cvr2, log(cvr2));
365 MSG("cvr = %f",cvr);
366 return cv0 + cvr;// J/K*Kg
367 }
368
369 /**
370 Isobaric specific heat capacity
371 \f[ c_p = \left(\frac{\partial h}{\partial T} \right)_p \f]
372
373 See Pratt, 2001. "Thermodynamic Properties Involving Derivatives",
374 ASEE Chemical Engineering Division Newsletter, Malaysia
375 http://www.che.uah.edu/courseware/che641/peng-robinson-derivatives-joule-thompson.pdf
376
377 TODO this function needs to be checked.
378 */
379 double pengrob_cp(double T, double rho, const FluidData *data, FpropsError *err){
380 //these calculations are broken apart intentionally to help provide clarity as the calculations are tedious
381 DEFINE_SQRTALPHA;
382 DEFINE_A;
383 DEFINE_V;
384 DEFINE_DADT;
385 DEFINE_D2ADT2;
386 DEFINE_CVR;
387 double cp0 = ideal_cp(T, rho, data, err);
388 DEFINE_DPDT_RHO;
389
390 #define DEFINE_CPR \
391 double A = p * a / SQ(data->R * T); \
392 double dAdT_p = p / SQ(data->R*T) * (dadT - 2*a/ T); \
393 double dBdT_p = - PD->b * p / (data->R * SQ(T)); \
394 double num = dAdT_p * (B-Z) + dBdT_p * (6*B*Z + 2*Z - 3*SQ(B) - 2*B + A - SQ(Z)); \
395 double den = 3*SQ(Z) + 2*(B-1)*Z + (A - 2*B - 3*SQ(B)); \
396 double dZdT_p = num / den; \
397 double dvdT_p = data->R / p * (T * dZdT_p + Z); \
398 double cpr = cvr + T * dpdT_rho * dvdT_p - data->R
399 DEFINE_CPR;
400 return cp0 + cpr;
401 }
402
403 /**
404 Speed of sound (see ��engel and Boles, 2012. "Thermodynamics, An Engineering Approach", McGraw-Hill, 7SIe):
405 \f[ w = \sqrt{ \left(\frac{\partial p}{\partial \rho} \right)_s }
406 = v \sqrt{ -\frac{c_p}{c_v} \left(\frac{\partial p}{\partial v}\right)_T}
407 \f]
408
409 See Pratt, 2001. "Thermodynamic Properties Involving Derivatives",
410 ASEE Chemical Engineering Division Newsletter, Malaysia
411 http://www.che.uah.edu/courseware/che641/peng-robinson-derivatives-joule-thompson.pdf
412
413 TODO this function needs to be checked.
414 */
415 double pengrob_w(double T, double rho, const FluidData *data, FpropsError *err){
416 DEFINE_SQRTALPHA;
417 DEFINE_V;
418 DEFINE_DADT;
419 DEFINE_D2ADT2;
420 DEFINE_A;
421 DEFINE_DPDT_RHO;
422 double cv0 = ideal_cv(T, rho, data, err);
423 double cp0 = cv0 + data->R;
424 DEFINE_CVR;
425 DEFINE_CPR;
426 double k = (cp0 + cpr) / (cv0 + cvr);
427 double dpdv_T = - SQ(rho) * pengrob_dpdrho_T(T,rho,data,err);
428 return v * sqrt(-k * dpdv_T);
429 }
430
431 double pengrob_dpdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
432 DEFINE_SQRTALPHA;
433 DEFINE_A;
434 DEFINE_V;
435 #define b PD->b
436 // CHECK: I check this...I think the previous expression was for dpdv_T
437 return -SQ(v)*( data->R*T / SQ(v-b) - 2 * a * (v + b) / SQ(v*(v+b) + b*(v-b)));
438 #undef b
439 }
440
441 /**
442 Relative pressure coefficient, as defined in http://www.iapws.org/relguide/Advise3.pdf
443 \f[ \alpha_p = \frac{1}{p} \left( \frac{\partial p}{\partial T} \right)_v \f]
444
445 Maxima code:
446 p(T,v) := R*T/(v-b) - a(T)/(v*(v+b)+b*(v-b))
447 alpha_p = 1/p * diff(p(T,v),T)
448
449 TODO the function is not yet checked/tested.
450 */
451 double pengrob_alphap(double T, double rho, const FluidData *data, FpropsError *err){
452 DEFINE_SQRTALPHA;
453 DEFINE_V;
454 DEFINE_DADT;
455 double p = pengrob_p(T, rho, data, err);
456 DEFINE_DPDT_RHO;
457 return 1/p * dpdT_rho;
458 }
459
460
461 /**
462 Isothermal stress coefficient, as defined in http://www.iapws.org/relguide/Advise3.pdf
463 \f[ \beta_p = - \frac{1}{p} \left( \frac{\partial p}{\partial v} \right)_T \f]
464
465 Maxima code:
466 p(T,v) := R*T/(v-b) - a(T)/(v*(v+b)+b*(v-b))
467 */
468 double pengrob_betap(double T, double rho, const FluidData *data, FpropsError *err){
469 double p = pengrob_p(T, rho, data, err);
470 return -1/p * SQ(rho) * pengrob_dpdrho_T(T,rho,data,err);
471 }
472
473 /**
474 Saturation calculation for a pure Peng-Robinson fluid. Algorithm as
475 outlined in Sandler 5e, sect 7.5. Another source of information is at
476 https://www.e-education.psu.edu/png520/m17.html
477
478 @return psat(T)
479 @param T temperature [K]
480 @param rhof_ret (output) saturated liquid density
481 @param rhog_ret (output) saturated vapour density
482 @param err (output) error flag, if any. should be reset first by caller if needed.
483
484 TODO implement accelerated successive substitution method (ASSM) or the
485 'Minimum Variable Newton Raphson (MVNR) Method' as detailed at the above
486 link.
487
488 Possible reference: Ghanbari & Davoodi Majd, 2004, "Liquid-Vapor Equilibrium
489 Curves for Methane System by Using Peng-Robinson Equation of State",
490 Petroleum & Coal 46(1) 23-27.
491 http://www.vurup.sk/sites/vurup.sk/archivedsite/www.vurup.sk/pc/vol46_2004/issue1/pdf/pc_ghanbari.pdf
492 http://www.doaj.org/doaj?func=abstract&id=251368
493 */
494 double pengrob_sat(double T,double *rhof_ret, double *rhog_ret, const FluidData *data, FpropsError *err){
495 /* rewritten, using GSOC code from Sean Muratet as a starting point -- thanks Sean! */
496
497 if(fabs(T - data->T_c) < 1e-3){
498 MSG("Saturation conditions requested at critical temperature");
499 *rhof_ret = data->rho_c;
500 *rhog_ret = data->rho_c;
501 return data->p_c;
502 }
503
504 //double rhog;
505 double A, B;
506 double sqrt2 = sqrt(2);
507
508 double p = fprops_psat_T_acentric(T, data);
509 MSG("Initial guess: p = %f from acentric factor",p);
510
511 int i = 0;
512 double Zg, Z1, Zf, vg, vf;
513 double oldfratio = 1e9;
514
515 #ifdef PR_DEBUG
516 FILE *F1 = fopen("pf.txt","w");
517 #endif
518
519 // FIXME test upper iteration limit required
520 while(++i < 200){
521 MSG("iter %d: p = %f, rhof = %f, rhog = %f", i, p, 1/vf, 1/vg);
522 // Peng & Robinson eq 17
523 double sqrtalpha = 1. + PD->kappa * (1. - sqrt(T / PD_TCRIT));
524 // Peng & Robinson eq 12
525 double a = PD->aTc * SQ(sqrtalpha);
526 // Peng & Robinson eq 6
527 A = a * p / SQ(data->R*T);
528 B = PD->b * p / (data->R*T);
529
530 // use GSL function to return real roots of polynomial: Peng & Robinson eq 5
531 Zf = 0; Z1 = 0; Zg = 0;
532 if(3 == cubicroots(-(1.-B), A-3.*SQ(B)-2.*B, -(A*B-SQ(B)*(1.+B)), &Zf,&Z1,&Zg)){
533 assert(Zf < Z1);
534 assert(Z1 < Zg);
535
536 //MSG(" roots: Z = %f, %f, %f", Zf, Z1, Zg);
537 //MSG(" Zf = %f, Zg = %f", Zf, Zg);
538 // three real roots in this case
539 vg = Zg*data->R*T / p;
540 vf = Zf*data->R*T / p;
541 if(vf < 0 || vg < 0){
542 // FIXME find out what does this mean; how does it happen?
543 MSG("Got a density root less than 0");
544 *err = FPROPS_SAT_CVGC_ERROR;
545 return 0;
546 }
547 //MSG(" vf = %f, vg = %f",vf,vg);
548 //MSG(" VMf = %f, VMg = %f",vf*data->M,vg*data->M);
549
550 // TODO can we use a function pointer for the fugacity expression, so that we can reuse this function for other cubic EOS?
551 #define FUG(Z,A,B) \
552 exp( (Z-1) - log(Z - B) - A/(2*sqrt2*B)*log( (Z + (1+sqrt2) * B) / (Z + (1-sqrt2) * B)) )
553
554 double ff = FUG(Zf,A,B);
555 double fg = FUG(Zg,A,B);
556 double fratio = ff/fg;
557 MSG(" ff = %f, fg = %f, fratio = %f", ff, fg, fratio);
558
559 //double hf = pengrob_h(T, 1/vf, data, err);
560 //double hg = pengrob_h(T, 1/vg, data, err);
561 //MSG(" HMf = %f, HMg = %f", hf*data->M/1000, hg*data->M/1000);
562
563 if(fabs(fratio - 1) < 1e-7){
564 *rhof_ret = 1 / vf;
565 *rhog_ret = 1 / vg;
566 p = pengrob_p(T, *rhog_ret, data, err);
567 MSG("Solved for T = %f: p = %f, rhof = %f, rhog = %f", T, p, *rhof_ret, *rhog_ret);
568 #ifdef PR_DEBUG
569 fclose(F1);
570 #endif
571 return p;
572 }
573 #ifdef PR_DEBUG
574 fprintf(F1,"%f\t%f\t%f\n",p,ff,fg);
575 #endif
576 if(fratio > oldfratio){
577 MSG("fratio increased!");
578 //fratio = 0.5 + 0.5*fratio;
579 }
580 p *= fratio;
581 if(p < 0){
582 p = 0.5 * p/fratio;
583 }
584 oldfratio = fratio;
585 }else{
586 MSG("Midpoint pressure calculation");
587 /* In this case we need to adjust our guess p(T) such that we get
588 into the narrow range of values that gives multiple solutions. */
589 p = MidpointPressureCubic(T, data, err);
590 if(*err){
591 ERRMSG("Failed to solve for a midpoint pressure");
592 #ifdef PR_DEBUG
593 fclose(F1);
594 #endif
595 return p;
596 }
597 MSG(" single root: Z = %f. new pressure guess: %f", Zf, p);
598 }
599 }
600 MSG("Did not converge");
601 *err = FPROPS_SAT_CVGC_ERROR;
602 #ifdef PR_DEBUG
603 fclose(F1);
604 #endif
605 return 0;
606 }
607
608 /*
609 FIXME can we generalise this to work with other *cubic* EOS as well?
610 Currently not easily done since pointers are kept at a higher level in our
611 data structures than FluidData.
612 */
613
614 typedef struct{
615 const FluidData *data;
616 FpropsError *err;
617 double T;
618 } MidpointSolveData;
619
620 static ZeroInSubjectFunction resid_dpdrho_T;
621
622 /**
623 This function is trying to find the locations of the stationary points
624 of the p(rho,T) function at a specified temperature, assumed to be
625 close to the critical temperature.
626
627 We return the midpoint of the pressure between these two.
628
629 It is assumed that if using this function we are close enough to the
630 critical point that vf < vc < vg, and that the stationary points
631 will be within those interval, ie vf < v1 < vc < v2 < vg.
632 */
633 double MidpointPressureCubic(double T, const FluidData *data, FpropsError *err){
634 MidpointSolveData msd = {data, err, T};
635 double rhomin = 0.5 * data->rho_c;
636 double rhomax = data->rho_c;
637 double rho, resid;
638
639 if(T > data->T_c){
640 ERRMSG("Invalid temperature T > T_c");
641 *err = FPROPS_RANGE_ERROR;
642 return data->p_c;
643 }
644
645 // look for a stationary point in density range less than rho_c
646 int res = zeroin_solve(&resid_dpdrho_T, &msd, rhomin, rhomax, 1e-9, &rho, &resid);
647 if(res){
648 ERRMSG("Failed to solve density for first stationary point");
649 *err = FPROPS_NUMERIC_ERROR;
650 return data->p_c;
651 }
652 double p1 = pengrob_p(T,rho, data, err);
653
654 // look for the other stationary point in density range less than rho_c
655 rhomin = data->rho_c;
656 rhomax = 2 * data->rho_c;
657 if(rhomax + 1e-2 > 1./(data->corr.pengrob->b)) rhomax = 1./(data->corr.pengrob->b) - 1e-3;
658
659 res = zeroin_solve(&resid_dpdrho_T, &msd, rhomin, rhomax, 1e-9, &rho, &resid);
660 if(res){
661 ERRMSG("Failed to solve density for second stationary point");
662 *err = FPROPS_NUMERIC_ERROR;
663 return data->p_c;
664 }
665
666 double p2 = pengrob_p(T,rho, data, err);
667 return 0.5*(p1 + p2);
668 }
669
670 static double resid_dpdrho_T(double rho, void *user_data){
671 #define D ((MidpointSolveData *)user_data)
672 return pengrob_dpdrho_T(D->T,rho,D->data,D->err);
673 #undef D
674 }
675
676
677 #if 0
678
679 /* we would like to try the Akasaka approach and use it with cubic EOS.
680 but at this stage it's not working correctly. Not sure why... */
681
682 /**
683 Solve saturation condition for a specified temperature using approach of
684 Akasaka, but adapted for general use to non-helmholtz property correlations.
685 @param T temperature [K]
686 @param psat_out output, saturation pressure [Pa]
687 @param rhof_out output, saturated liquid density [kg/m^3]
688 @param rhog_out output, saturated vapour density [kg/m^3]
689 @param d helmholtz data object for the fluid in question.
690 @return 0 on success, non-zero on error (eg algorithm failed to converge, T out of range, etc.)
691 */
692 double pengrob_sat_akasaka(double T, double *rhof_out, double * rhog_out, const FluidData *data, FpropsError *err){
693 if(T < data->T_t - 1e-8){
694 ERRMSG("Input temperature %f K is below triple-point temperature %f K",T,data->T_t);
695 return FPROPS_RANGE_ERROR;
696 }
697
698 if(T > data->T_c){
699 ERRMSG("Input temperature is above critical point temperature");
700 return FPROPS_RANGE_ERROR;
701 }
702
703 // we're at the critical point
704 if(fabs(T - data->T_c) < 1e-9){
705 *rhof_out = data->rho_c;
706 *rhog_out = data->rho_c;
707 return data->p_c;
708 }
709
710 // FIXME at present step-length multiplier is set to 0.4 just because of
711 // ONE FLUID, ethanol. Probably our initial guess data isn't good enough,
712 // or maybe there's a problem with the acentric factor or something like
713 // that. This factor 0.4 will be slowing down the whole system, so it's not
714 // good. TODO XXX.
715
716 // initial guesses for liquid and vapour density
717 double rhof = fprops_rhof_T_rackett(T,data);
718 double rhog= fprops_rhog_T_chouaieb(T,data);
719 double R = data->R;
720 double pc = data->p_c;
721
722 MSG("initial guess rho_f = %f, rho_g = %f\n",rhof,rhog);
723 MSG("calculating for T = %.12e",T);
724
725 int i = 0;
726 while(i++ < 70){
727 if(rhof > 1/PD->b){
728 MSG("limit rhof to 1/b");
729 rhof = 1/PD->b;
730 }
731
732 MSG("iter %d: T = %f, rhof = %f, rhog = %f",i,T, rhof, rhog);
733
734 double pf = pengrob_p(T,rhof,data,err);
735 double pg = pengrob_p(T,rhog,data,err);
736 double gf = pengrob_g(T,rhof,data,err);
737 double gg = pengrob_g(T,rhog,data,err);
738 double dpdrf = pengrob_dpdrho_T(T,rhof,data,err);
739 double dpdrg = pengrob_dpdrho_T(T,rhog,data,err);
740 if(*err){
741 ERRMSG("error returned");
742 }
743
744 MSG("gf = %f, gg = %f", gf, gg);
745
746 // jacobian for [F;G](rhof, rhog) --- derivatives wrt rhof and rhog
747 double F = (pf - pg)/pc;
748 double G = (gf - gg)/R/T;
749
750 if(fabs(F) + fabs(G) < 1e-12){
751 //fprintf(stderr,"%s: CONVERGED\n",__func__);
752 *rhof_out = rhof;
753 *rhog_out = rhog;
754 return pengrob_p(T, *rhog_out, data, err);
755 /* SUCCESS */
756 }
757
758 double Ff = dpdrf/pc;
759 double Fg = -dpdrg/pc;
760 MSG("Ff = %e, Fg = %e",Ff,Fg);
761
762 double Gf = dpdrf/rhof/R/T;
763 double Gg = -dpdrg/rhog/R/T;
764 MSG("Gf = %e, Gg = %e",Gf,Gg);
765
766 double DET = Ff*Gg - Fg*Gf;
767 MSG("DET = %f",DET);
768
769 MSG("F = %f, G = %f", F, G);
770
771 // 'gamma' needs to be increased to 0.5 for water to solve correctly (see 'test/sat.c')
772 #define gamma 1
773 rhof += gamma/DET * (Fg*G - Gg*F);
774 rhog += gamma/DET * ( Gf*F - Ff*G);
775 #undef gamma
776
777 if(rhog < 0)rhog = -0.5*rhog;
778 if(rhof < 0)rhof = -0.5*rhof;
779 }
780 *rhof_out = rhof;
781 *rhog_out = rhog;
782 *err = FPROPS_SAT_CVGC_ERROR;
783 ERRMSG("Not converged: with T = %e (rhof=%f, rhog=%f).",T,*rhof_out,*rhog_out);
784 return pengrob_p(T, rhog, data, err);
785 }
786
787 #endif
788
789

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