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

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