1 //-------------------------------------------------------------------
2 // $Id: s_fgl.cpp 724 2012-10-02 14:25:25Z kulik $
3 //
4 /// \file s_solmod2.cpp
5 /// Implementation of TSolMod derived classes for fluid phase models
6 /// (TPRSVcalc, TCGFcalc, TSRKcalc, TPR78calc, TCORKcalc and TSTPcalc classes)
7 //
8 // Copyright (c) 2004-2012 T.Wagner, D.Kulik, S. Dmitrieva, S.Churakov
9 // <GEMS Development Team, mailto:gems2.support@psi.ch>
10 //
11 // This file is part of the GEMS4K code for thermodynamic modelling
12 // by Gibbs energy minimization <http://gems.web.psi.ch/GEMS4K/>
13 //
14 // GEMS4K is free software: you can redistribute it and/or modify
15 // it under the terms of the GNU Lesser General Public License as
16 // published by the Free Software Foundation, either version 3 of
17 // the License, or (at your option) any later version.
18
19 // GEMS4K is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU Lesser General Public License for more details.
23
24 // You should have received a copy of the GNU General Public License
25 // along with GEMS4K code. If not, see <http://www.gnu.org/licenses/>.
26 //-------------------------------------------------------------------
27
28 #include <cmath>
29 #include <cstdio>
30 #include <iostream>
31 #include <iomanip>
32 #include <fstream>
33 using namespace std;
34 #include "s_solmod_.h"
35 #include "verror.h"
36
37 namespace solmod
38 {
39 //=======================================================================================================
40 // Peng-Robinson-Stryjek-Vera (PRSV) model for fluid mixtures
41 // References: Stryjek and Vera (1986), Proust and Vera (1989)
42 // (c) TW July 2006
43 //=======================================================================================================
44
45
46 // generic constructor
TPRSVcalc(long int NCmp,double Pp,double Tkp)47 TPRSVcalc::TPRSVcalc( long int NCmp, double Pp, double Tkp ):
48 TSolMod( NCmp, 'P', Tkp, Pp )
49
50 {
51 // aGEX = 0;
52 // aVol = 0;
53 Pparc = 0;
54 alloc_internal();
55 }
56
57
58
TPRSVcalc(SolutionData * sd)59 TPRSVcalc::TPRSVcalc( SolutionData *sd ):
60 TSolMod( sd )
61 {
62 // aGEX = arGEX;
63 // aVol = arVol;
64 Pparc = aPparc;
65 alloc_internal();
66 }
67
68
69
~TPRSVcalc()70 TPRSVcalc::~TPRSVcalc()
71 {
72 free_internal();
73 }
74
75
76
77 /// allocate work arrays for pure fluid and fluid mixture properties
alloc_internal()78 void TPRSVcalc::alloc_internal()
79 {
80 Eosparm = new double [NComp][6];
81 Pureparm = new double [NComp][4];
82 Fugpure = new double [NComp][6];
83 Fugci = new double [NComp][4];
84 a = new double *[NComp];
85 b = new double *[NComp];
86 KK = new double *[NComp];
87 dKK = new double *[NComp];
88 d2KK = new double *[NComp];
89 AA = new double *[NComp];
90
91 for (long int i=0; i<NComp; i++)
92 {
93 a[i] = new double[NComp];
94 b[i] = new double[NComp];
95 KK[i] = new double[NComp];
96 dKK[i] = new double[NComp];
97 d2KK[i] = new double[NComp];
98 AA[i] = new double[NComp];
99 }
100 }
101
102
free_internal()103 void TPRSVcalc::free_internal()
104 {
105 long int i;
106
107 for (i=0; i<NComp; i++)
108 {
109 delete[]a[i];
110 delete[]b[i];
111 delete[]KK[i];
112 delete[]dKK[i];
113 delete[]d2KK[i];
114 delete[]AA[i];
115 }
116
117 delete[]Eosparm;
118 delete[]Pureparm;
119 delete[]Fugpure;
120 delete[]Fugci;
121 delete[]a;
122 delete[]b;
123 delete[]KK;
124 delete[]dKK;
125 delete[]d2KK;
126 delete[]AA;
127 }
128
129
130
131 /// high-level method to retrieve pure fluid fugacities
PureSpecies()132 long int TPRSVcalc::PureSpecies()
133 {
134 long int j, retCode = 0;
135
136 for( j=0; j<NComp; j++)
137 {
138 // Calling PRSV EoS for pure fugacity
139 retCode = FugacityPT( j, aDCc+j*NP_DC );
140 aGEX[j] = log( Fugpure[j][0] );
141 Pparc[j] = Fugpure[j][0]*Pbar; // fure fluid fugacity (required for performance)
142 aVol[j] = Fugpure[j][4]*10.; // molar volume of pure fluid component, J/bar to cm3
143 } // j
144
145 if ( retCode )
146 {
147 char buf[150];
148 sprintf(buf, "PRSV fluid: calculation of pure fugacity failed");
149 Error( "E71IPM IPMgamma: ", buf );
150 }
151 return 0;
152 }
153
154
155
156 /// high-level method to calculate T,P corrected binary interaction parameters
PTparam()157 long int TPRSVcalc::PTparam()
158 {
159 long int j, i;
160
161 PureSpecies();
162
163 // set all interaction parameters zero
164 for( j=0; j<NComp; j++ )
165 {
166 for( i=0; i<NComp; i++ )
167 {
168 KK[j][i] = 0.;
169 dKK[j][i] = 0.;
170 d2KK[j][i] = 0.;
171 }
172 }
173
174 switch ( MixCode )
175 {
176 case MR_UNDEF_:
177 case MR_WAAL_:
178 MixingWaals();
179 break;
180 case MR_CONST_:
181 MixingConst();
182 break;
183 case MR_TEMP_:
184 MixingTemp();
185 break;
186 default:
187 break;
188 }
189
190 return 0;
191 }
192
193
194
195 /// high-level method to retrieve activity coefficients of the fluid mixture
MixMod()196 long int TPRSVcalc::MixMod()
197 {
198 long int j, iRet;
199
200 iRet = FugacitySpec( aPparc );
201
202 phVOL[0] = PhVol * 10.;
203
204 for(j=0; j<NComp; j++)
205 {
206 if( Fugci[j][3] > 1e-23 )
207 lnGamma[j] = log( Fugci[j][3] );
208 else
209 lnGamma[j] = 0;
210 }
211 if ( iRet )
212 {
213 char buf[150];
214 sprintf(buf, "PRSV fluid: calculation failed");
215 Error( "E71IPM IPMgamma: ", buf );
216 }
217 return iRet;
218 }
219
220
221
222 /// high-level method to retrieve residual functions of the fluid mixture
ExcessProp(double * Zex)223 long int TPRSVcalc::ExcessProp( double *Zex )
224 {
225 long int iRet;
226
227 iRet = ResidualFunct( aPparc );
228
229 if ( iRet )
230 {
231 char buf[150];
232 sprintf(buf, "PRSV fluid: calculation failed");
233 Error( "E71IPM IPMgamma: ", buf );
234 }
235
236 Ars = Grs - Vrs*Pbar;
237 Urs = Hrs - Vrs*Pbar;
238
239 // assignments (residual functions)
240 Zex[0] = Grs;
241 Zex[1] = Hrs;
242 Zex[2] = Srs;
243 Zex[3] = CPrs;
244 Zex[4] = Vrs;
245 Zex[5] = Ars;
246 Zex[6] = Urs;
247
248 return iRet;
249 }
250
251
252
253 /// calculates ideal mixing properties
IdealProp(double * Zid)254 long int TPRSVcalc::IdealProp( double *Zid )
255 {
256 long int j;
257 double s, sc, sp;
258
259 s = 0.0;
260 for (j=0; j<NComp; j++)
261 {
262 if ( x[j] > 1.0e-32 )
263 s += x[j]*log(x[j]);
264 }
265 sc = (-1.)*R_CONST*s;
266 sp = (-1.)*R_CONST*log(Pbar);
267 Hid = 0.0;
268 CPid = 0.0;
269 Vid = 0.0;
270 Sid = sc + sp;
271 Gid = Hid - Sid*Tk;
272 Aid = Gid - Vid*Pbar;
273 Uid = Hid - Vid*Pbar;
274
275 // assignments (ideal mixing properties)
276 Zid[0] = Gid;
277 Zid[1] = Hid;
278 Zid[2] = Sid;
279 Zid[3] = CPid;
280 Zid[4] = Vid;
281 Zid[5] = Aid;
282 Zid[6] = Uid;
283
284 return 0;
285 }
286
287
288
289 /// basic van der waals mixing rule
MixingWaals()290 long int TPRSVcalc::MixingWaals()
291 {
292 // currently no calculations
293
294 return 0;
295 }
296
297
298
299 /// constant one-term interaction parameter
MixingConst()300 long int TPRSVcalc::MixingConst()
301 {
302 long int ip, i1, i2;
303 double k, dk, d2k;
304
305 if( NPcoef > 0 )
306 {
307 // transfer interaction parameters that have non-standard value
308 for ( ip=0; ip<NPar; ip++ )
309 {
310 i1 = aIPx[MaxOrd*ip];
311 i2 = aIPx[MaxOrd*ip+1];
312 k = aIPc[NPcoef*ip];
313 dk = 0.;
314 d2k = 0.;
315 KK[i1][i2] = k;
316 dKK[i1][i2] = dk;
317 d2KK[i1][i2] = d2k;
318 KK[i2][i1] = k; // symmetric case
319 dKK[i2][i1] = dk;
320 d2KK[i2][i1] = d2k;
321 }
322 }
323
324 return 0;
325 }
326
327
328
329 /// temperature dependent one-term interaction parameter
MixingTemp()330 long int TPRSVcalc::MixingTemp()
331 {
332 long int i, j, ip, i1, i2;
333 double ai, aj, bi, bj, di, dj, dai, daj, d2ai, d2aj, ddi, ddj, d2di, d2dj,
334 U, V, dU, dV, d2U, d2V, tmp, k, dk, d2k, C;
335
336 // set model specific interaction parameters zero
337 for( j=0; j<NComp; j++ )
338 {
339 for( i=0; i<NComp; i++ )
340 {
341 a[j][i] = 0.;
342 b[j][i] = 0.;
343 }
344 }
345
346 if( NPcoef > 0 )
347 {
348 // transfer parameters that have non-standard value
349 for ( ip=0; ip<NPar; ip++ )
350 {
351 i1 = aIPx[MaxOrd*ip];
352 i2 = aIPx[MaxOrd*ip+1];
353 a[i1][i2] = aIPc[NPcoef*ip];
354 b[i1][i2] = aIPc[NPcoef*ip+1];
355 a[i2][i1] = aIPc[NPcoef*ip]; // symmetric case
356 b[i2][i1] = aIPc[NPcoef*ip+1];
357 }
358 }
359
360 // calculate binary interaction parameters
361 for (i=0; i<NComp; i++)
362 {
363 for (j=0; j<NComp; j++)
364 {
365 if ( a[i][j] == 0.0 )
366 tmp = 1.0;
367 else
368 tmp = a[i][j];
369
370 // read a, b, da, d2a
371 ai = Pureparm[i][0];
372 aj = Pureparm[j][0];
373 bi = Pureparm[i][1];
374 bj = Pureparm[j][1];
375 dai = Pureparm[i][2];
376 daj = Pureparm[j][2];
377 d2ai = Pureparm[i][3];
378 d2aj = Pureparm[j][3];
379
380 // calculate k and derivatives
381 di = sqrt(ai)/bi;
382 dj = sqrt(aj)/bj;
383 ddi = (0.5/bi) * pow(ai,-0.5) * dai;
384 ddj = (0.5/bj) * pow(aj,-0.5) * daj;
385 d2di = (0.5/bi) * ( (-0.5)*pow(ai,-1.5)*dai*dai + pow(ai,-0.5)*d2ai );
386 d2dj = (0.5/bj) * ( (-0.5)*pow(aj,-1.5)*daj*daj + pow(aj,-0.5)*d2aj );
387
388 C = ( b[i][j]/tmp - 1. );
389 U = a[i][j]*pow((298.15/Tk),C) - pow((di-dj),2.);
390 V = (2.*di*dj);
391 dU = - ( a[i][j]*C*pow((298.15/Tk),C) ) / Tk - 2.*(di-dj)*(ddi-ddj);
392 dV = 2.*( ddi*dj + di*ddj );
393 d2U = ( a[i][j]*pow(C,2.)*pow((298.15/Tk),C) ) / pow(Tk,2.)
394 + ( a[i][j]*C*pow((298.15/Tk),C) ) / pow(Tk,2.)
395 - 2.*( pow((ddi-ddj),2.) + (di-dj)*(d2di-d2dj) );
396 d2V = 2.*( d2di*dj + 2.*ddi*ddj + di*d2dj );
397 k = U/V;
398 dk = (dU*V-U*dV)/pow(V,2.);
399 d2k = (d2U*V+dU*dV)*pow(V,2.)/pow(V,4.) - (dU*V)*(2.*V*dV)/pow(V,4.)
400 - (dU*dV+U*d2V)*pow(V,2.)/pow(V,4.) + (U*dV)*(2.*V*dV)/pow(V,4.);
401
402 // assignments
403 KK[i][j] = k;
404 dKK[i][j] = dk;
405 d2KK[i][j] = d2k;
406 }
407 }
408
409 return 0;
410 }
411
412
413
414 /// retrieve pure fluid properties
FugacityPT(long int i,double * EoSparam)415 long int TPRSVcalc::FugacityPT( long int i, double *EoSparam )
416 {
417 long int iRet = 0;
418 double Tcrit, Pcrit, omg, k1, k2, k3, apure, bpure, da, d2a;
419
420 // reads EoS parameters from database into work array
421 if( !EoSparam )
422 return -1; // Memory alloc error
423
424 Eosparm[i][0] = EoSparam[0]; // critical temperature in K
425 Eosparm[i][1] = EoSparam[1]; // critical pressure in bar
426 Eosparm[i][2] = EoSparam[2]; // Pitzer acentric factor omega
427 Eosparm[i][3] = EoSparam[3]; // empirical EoS parameter k1
428 Eosparm[i][4] = EoSparam[4]; // empirical EoS parameter k2
429 Eosparm[i][5] = EoSparam[5]; // empirical EoS parameter k3
430 Tcrit = Eosparm[i][0];
431 Pcrit = Eosparm[i][1];
432 omg = Eosparm[i][2];
433 k1 = Eosparm[i][3];
434 k2 = Eosparm[i][4];
435 k3 = Eosparm[i][5];
436
437 AB( Tcrit, Pcrit, omg, k1, k2, k3, apure, bpure, da, d2a );
438
439 Pureparm[i][0] = apure; // a parameter
440 Pureparm[i][1] = bpure; // b parameter
441 Pureparm[i][2] = da; // da/dT
442 Pureparm[i][3] = d2a; // d2a/dT2
443
444 iRet = FugacityPure( i );
445 if( iRet)
446 return iRet;
447
448 return iRet;
449 }
450
451
452
453 /// calculates attractive (a) and repulsive (b) parameter of PRSV equation of state
454 /// and partial derivatives of alpha function
AB(double Tcrit,double Pcrit,double omg,double k1,double k2,double k3,double & apure,double & bpure,double & da,double & d2a)455 long int TPRSVcalc::AB( double Tcrit, double Pcrit, double omg, double k1, double k2, double k3,
456 double &apure, double &bpure, double &da, double &d2a )
457 {
458 double Tred, k0, k, alph, ac, sqa, dsqa, d2sqa;
459
460 Tred = Tk/Tcrit;
461 k0 = 0.378893 + 1.4897153*omg - 0.17131848*pow(omg,2.) + 0.0196554*pow(omg,3.);
462 if(Tk >= Tcrit)
463 {
464 k1 = 0.0;
465 k2 = 0.0;
466 k3 = 0.0;
467 }
468 k = k0 + (k1 + k2*(k3-Tred)*(1.-sqrt(Tred))) * (1.+sqrt(Tred)) * (0.7-Tred);
469 alph = pow(1. + k*(1.-sqrt(Tred)), 2.);
470 ac = (0.457235)*pow(R_CONST,2.)*pow(Tcrit,2.) / Pcrit;
471 apure = alph*ac;
472 bpure = (0.077796)*R_CONST*Tcrit/Pcrit;
473 sqa = 1.+k*(1.-sqrt(Tred));
474 // dsqa = (-1.)*k0/(2.*sqrt(Tk*Tcrit)) - 1.7*k1/Tcrit + 2.*k1*Tk/(pow(Tcrit,2.)); // extend dA/dT for k2, k3
475 dsqa = ( k1*(0.7-Tred)/(2.*sqrt(Tred)*Tcrit) - k1*(1.+sqrt(Tred))/Tcrit ) * (1.-sqrt(Tred))
476 - (k0 + k1*(1.+sqrt(Tred))*(0.7-Tred))/(2*sqrt(Tred)*Tcrit);
477 da = 2.*ac*(sqa*dsqa);
478 d2sqa = ( - (k1*(0.7-Tred))/(4.*pow(Tred,1.5)*pow(Tcrit,2.)) - k1/(sqrt(Tred)*pow(Tcrit,2.)) ) * (1.-sqrt(Tred))
479 + ( - k1*(0.7-Tred)/(2.*sqrt(Tred)*Tcrit) - k1*(1.+sqrt(Tred))/Tcrit ) / (sqrt(Tred)*Tcrit)
480 + ( k0 + k1*(1.+sqrt(Tred))*(0.7-Tred) ) / (4.*pow(Tred,1.5)*pow(Tcrit,2.));
481 d2a = 2.*ac*(dsqa*dsqa + sqa*d2sqa);
482
483 return 0;
484 }
485
486
487
488 /// calculates fugacities and residual functions of pure fluid species
FugacityPure(long int i)489 long int TPRSVcalc::FugacityPure( long int i )
490 {
491 double Tcrit, Pcrit, Tred, aprsv, bprsv, alph, da, d2a, k, A, B, a2, a1, a0,
492 z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, z, vol, lnf;
493 double /*gig, hig, sig, cpig,*/ fugpure, grs, hrs, srs, cprs,
494 cv, dPdT, dPdV, dVdT;
495
496 // ideal gas changes from 1 bar to P at T of interest
497 //hig = 0.;
498 //sig = (-1.)*R_CONST*log(Pbar);
499 //gig = hig - Tk*sig;
500 //cpig = 0.;
501
502 // retrieve a and b terms of cubic EoS
503 Tcrit = Eosparm[i][0];
504 Pcrit = Eosparm[i][1];
505 Tred = Tk/Tcrit;
506 aprsv = Pureparm[i][0];
507 bprsv = Pureparm[i][1];
508 da = Pureparm[i][2];
509 d2a = Pureparm[i][3];
510
511 // solve cubic equation
512 A = aprsv*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
513 B = bprsv*Pbar/(R_CONST*Tk);
514 a2 = B - 1.;
515 a1 = A - 3.*pow(B,2.) - 2.*B;
516 a0 = pow(B,3.) + pow(B,2.) - A*B;
517 Cardano(a2, a1, a0, z1, z2, z3);
518
519 // find stable roots
520 vol1 = z1*R_CONST*Tk/Pbar;
521 vol2 = z2*R_CONST*Tk/Pbar;
522 vol3 = z3*R_CONST*Tk/Pbar;
523 if (z1 > B)
524 lnf1 = (-1.)*log(z1-B)
525 - A/(B*sqrt(8.))*log((z1+(1.+sqrt(2.))*B)/(z1+(1.-sqrt(2.))*B))+z1-1.;
526 else
527 lnf1 = 1000.;
528 if (z2 > B)
529 lnf2 = (-1.)*log(z2-B)
530 - A/(B*sqrt(8.))*log((z2+(1.+sqrt(2.))*B)/(z2+(1.-sqrt(2.))*B))+z2-1.;
531 else
532 lnf2 = 1000.;
533 if (z3 > B)
534 lnf3 = (-1.)*log(z3-B)
535 - A/(B*sqrt(8.))*log((z3+(1.+sqrt(2.))*B)/(z3+(1.-sqrt(2.))*B))+z3-1.;
536 else
537 lnf3 = 1000.;
538
539 if (lnf2 < lnf1)
540 {
541 z = z2; vol = vol2; lnf = lnf2;
542 }
543 else
544 {
545 z = z1; vol = vol1; lnf = lnf1;
546 }
547 if (lnf3 < lnf)
548 {
549 z = z3; vol = vol3; lnf = lnf3;
550 }
551 else
552 {
553 z = z; vol = vol; lnf = lnf;
554 }
555
556 // calculate thermodynamic properties
557 alph = aprsv/((0.457235)*pow(R_CONST,2.)*pow(Tcrit,2.) / Pcrit);
558 k = (sqrt(alph)-1.)/(1.-sqrt(Tred));
559 grs = R_CONST*Tk*(z-1.-log(z-B)-A/(B*sqrt(8.))
560 *log((z+(1+sqrt(2.))*B)/(z+(1-sqrt(2.))*B)));
561 hrs = R_CONST*Tk*(z-1.-log((z+(1+sqrt(2.))*B)/(z+(1-sqrt(2.))*B))
562 *A/(B*sqrt(8.))*(1+k*sqrt(Tred)/sqrt(alph)));
563 srs = (hrs-grs)/Tk;
564
565 // heat capacity part
566 cv = Tk*d2a/(bprsv*sqrt(8.))
567 * log( (z+B*(1.+sqrt(2.)))/(z+B*(1.-sqrt(2.))) );
568 dPdT = R_CONST/(vol-bprsv) - da/( vol*(vol+bprsv) + bprsv*(vol-bprsv) );
569 dPdV = - R_CONST*Tk/pow((vol-bprsv),2.) + 2.*aprsv*(vol+bprsv)/pow((vol*(vol+bprsv)+bprsv*(vol-bprsv)),2.);
570 dVdT = (-1.)*(1./dPdV)*dPdT;
571 cprs = cv + Tk*dPdT*dVdT - R_CONST;
572
573 // increment thermodynamic properties
574 fugpure = exp(lnf);
575 Fugpure[i][0] = fugpure;
576 Fugpure[i][1] = grs; // changed to residual functions, 31.05.2008 (TW)
577 Fugpure[i][2] = hrs;
578 Fugpure[i][3] = srs;
579 Fugpure[i][4] = vol;
580 Fugpure[i][5] = cprs;
581
582 return 0;
583 }
584
585
586
587 /// cubic equation root solver based on Cardanos method
Cardano(double a2,double a1,double a0,double & z1,double & z2,double & z3)588 long int TPRSVcalc::Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 )
589 {
590 double q, rc, q3, rc2, theta, ac, bc;
591
592 q = (pow(a2,2.) - 3.*a1)/9.;
593 rc = (2.*pow(a2,3.) - 9.*a2*a1 + 27.*a0)/54.;
594 q3 = pow(q,3.);
595 rc2 = pow(rc,2.);
596 if (rc2 < q3) // three real roots
597 {
598 theta = acos(rc/sqrt(q3));
599 z1 = (-2.)*sqrt(q)*cos(theta/3.)-a2/3.;
600 z2 = (-2.)*sqrt(q)*cos(theta/3.+2./3.*3.1415927)-a2/3.;
601 z3 = (-2.)*sqrt(q)*cos(theta/3.-2./3.*3.1415927)-a2/3.;
602 }
603 else // one real root
604 {
605 ac = (-1.)*rc/fabs(rc)*pow(fabs(rc)+sqrt(rc2-q3), 1./3.);
606 if (ac != 0.)
607 bc = q/ac;
608 else
609 bc = 0.;
610 z1 = ac+bc-a2/3.;
611 z2 = ac+bc-a2/3.;
612 z3 = ac+bc-a2/3.;
613 }
614 return 0;
615 }
616
617
618
619 /// calculates mixing properties of the fluid mixture
MixParam(double & amix,double & bmix)620 long int TPRSVcalc::MixParam( double &amix, double &bmix )
621 {
622 long int i, j;
623 double K;
624 amix = 0.;
625 bmix = 0.;
626
627 // calculate binary aij parameters
628 for (i=0; i<NComp; i++)
629 {
630 for (j=0; j<NComp; j++)
631 {
632 K = KK[i][j];
633 AA[i][j] = sqrt(Pureparm[i][0]*Pureparm[j][0])*(1.-K);
634 }
635 }
636 // find a and b of the mixture
637 for (i=0; i<NComp; i++)
638 {
639 for (j=0; j<NComp; j++)
640 {
641 amix = amix + x[i]*x[j]*AA[i][j];
642 }
643 }
644 for (i=0; i<NComp; i++)
645 {
646 bmix = bmix + x[i]*Pureparm[i][1];
647 }
648 return 0;
649 }
650
651
652
653 /// calculates fugacity of the bulk fluid mixture
FugacityMix(double amix,double bmix,double & fugmix,double & zmix,double & vmix)654 long int TPRSVcalc::FugacityMix( double amix, double bmix, double &fugmix, double &zmix,
655 double &vmix )
656 {
657 double A, B, a2, a1, a0, z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, lnf;
658
659 // solve cubic equation
660 A = amix*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
661 B = bmix*Pbar/(R_CONST*Tk);
662 a2 = B - 1.;
663 a1 = A - 3.*pow(B,2.) - 2.*B;
664 a0 = pow(B,3.) + pow(B,2.) - A*B;
665 Cardano( a2, a1, a0, z1, z2, z3 );
666
667 // find stable roots
668 vol1 = z1*R_CONST*Tk/Pbar;
669 vol2 = z2*R_CONST*Tk/Pbar;
670 vol3 = z3*R_CONST*Tk/Pbar;
671 if (z1 > B)
672 lnf1 = (-1.)*log(z1-B)
673 - A/(B*sqrt(8.))*log((z1+(1.+sqrt(2.))*B)/(z1+(1.-sqrt(2.))*B))+z1-1.;
674 else
675 lnf1 = 1000.;
676 if (z2 > B)
677 lnf2 = (-1.)*log(z2-B)
678 - A/(B*sqrt(8.))*log((z2+(1.+sqrt(2.))*B)/(z2+(1.-sqrt(2.))*B))+z2-1.;
679 else
680 lnf2 = 1000.;
681 if (z3 > B)
682 lnf3 = (-1.)*log(z3-B)
683 - A/(B*sqrt(8.))*log((z3+(1.+sqrt(2.))*B)/(z3+(1.-sqrt(2.))*B))+z3-1.;
684 else
685 lnf3 = 1000.;
686
687 if (lnf2 < lnf1)
688 {
689 zmix = z2; vmix = vol2; lnf = lnf2;
690 }
691 else
692 {
693 zmix = z1; vmix = vol1; lnf = lnf1;
694 }
695 if (lnf3 < lnf)
696 {
697 zmix = z3; vmix = vol3; lnf = lnf3;
698 }
699 else
700 {
701 zmix = zmix; vmix = vmix; lnf = lnf;
702 }
703 fugmix = exp(lnf);
704 PhVol = vmix;
705 return 0;
706 }
707
708
709
710 /// calculates fugacities and activities of fluid species in the mixture,
FugacitySpec(double * fugpure)711 long int TPRSVcalc::FugacitySpec( double *fugpure )
712 {
713 long int i, j, iRet=0;
714 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0., sum=0.;
715 double A, B, lnfci, fci;
716
717 // Reload params to Pureparm
718 for( j=0; j<NComp; j++ )
719 {
720 Fugpure[j][0] = fugpure[j]/Pbar;
721 }
722
723 // retrieve properties of the mixture
724 iRet = MixParam( amix, bmix );
725 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix );
726 A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
727 B = bmix*Pbar/(R_CONST*Tk);
728
729 // calculate fugacity coefficient, fugacity and activity of species i
730 for (i=0; i<NComp; i++)
731 {
732 sum = 0.;
733 for (j=0; j<NComp; j++)
734 {
735 sum = sum + x[j]*AA[i][j];
736 }
737 lnfci = Pureparm[i][1]/bmix*(zmix-1.) - log(zmix-B)
738 + A/(sqrt(8.)*B)*(2.*sum/amix-Pureparm[i][1]/bmix)
739 * log((zmix+B*(1.-sqrt(2.)))/(zmix+B*(1.+sqrt(2.))));
740 fci = exp(lnfci);
741 Fugci[i][0] = fci; // fugacity coefficient using engineering convention
742 Fugci[i][1] = x[i]*fci; // fugacity coefficient using geology convention
743 Fugci[i][2] = Fugci[i][1]/Fugpure[i][0]; // activity of species
744 if (x[i]>1.0e-20)
745 Fugci[i][3] = Fugci[i][2]/x[i]; // activity coefficient of species
746 else
747 Fugci[i][3] = 1.0;
748 }
749
750 return iRet;
751 }
752
753
754
755 /// calculates residual functions in the mixture
ResidualFunct(double * fugpure)756 long int TPRSVcalc::ResidualFunct( double *fugpure )
757 {
758 long int i, j, iRet=0;
759 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0.;
760 double /*A,*/ B, K, dK, d2K, Q, dQ, d2Q, damix, d2amix, ai, aj, dai, daj, d2ai, d2aj,
761 cv, dPdT, dPdV, dVdT;
762
763 // Reload params to Pureparm (probably now obsolete?)
764 for( j=0; j<NComp; j++ )
765 {
766 Fugpure[j][0] = fugpure[j]/Pbar;
767 }
768
769 // retrieve properties of the mixture
770 iRet = MixParam( amix, bmix );
771 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix );
772 //A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
773 B = bmix*Pbar/(R_CONST*Tk);
774
775 // calculate total state functions of the mixture
776 damix = 0.;
777 d2amix = 0.;
778 for (i=0; i<NComp; i++)
779 {
780 for (j=0; j<NComp; j++)
781 {
782 // pull parameters
783 ai = Pureparm[i][0];
784 aj = Pureparm[j][0];
785 dai = Pureparm[i][2];
786 daj = Pureparm[j][2];
787 d2ai = Pureparm[i][3];
788 d2aj = Pureparm[j][3];
789 K = KK[i][j];
790 dK = dKK[i][j];
791 d2K = d2KK[i][j];
792
793 // increments to derivatives
794 Q = sqrt(ai*aj);
795 dQ = 0.5*( sqrt(aj/ai)*dai + sqrt(ai/aj)*daj );
796 d2Q = 0.5*( dai*daj/sqrt(ai*aj) + d2ai*sqrt(aj)/sqrt(ai) + d2aj*sqrt(ai)/sqrt(aj)
797 - 0.5*( pow(dai,2.)*sqrt(aj)/sqrt(pow(ai,3.))
798 + pow(daj,2.)*sqrt(ai)/sqrt(pow(aj,3.)) ) );
799 damix = damix + x[i]*x[j] * ( dQ*(1.-K) - Q*dK );
800 d2amix = d2amix + x[i]*x[j] * ( d2Q*(1.-K) - 2.*dQ*dK - Q*d2K );
801 }
802 }
803
804 // calculate thermodynamic properties
805 Grs = (amix/(R_CONST*Tk*sqrt(8.)*bmix) * log((vmix+(1.-sqrt(2.))*bmix)
806 / (vmix+(1.+sqrt(2.))*bmix))-log(zmix*(1.-bmix/vmix))+zmix-1.)*R_CONST*Tk;
807 Hrs = ((amix-Tk*damix)/(R_CONST*Tk*sqrt(8.)*bmix)*log((vmix+(1.-sqrt(2.))
808 *bmix)/(vmix+(1.+sqrt(2.))*bmix))+zmix-1.)*R_CONST*Tk;
809 Srs = (Hrs - Grs)/Tk;
810
811 // heat capacity part
812 cv = Tk*d2amix/(bmix*sqrt(8.))
813 * log( (zmix+B*(1.+sqrt(2.)))/(zmix+B*(1.-sqrt(2.))) );
814 dPdT = R_CONST/(vmix-bmix) - damix/( vmix*(vmix+bmix) + bmix*(vmix-bmix) );
815 dPdV = - R_CONST*Tk/pow((vmix-bmix),2.) + 2.*amix*(vmix+bmix)/pow((vmix*(vmix+bmix)+bmix*(vmix-bmix)),2.);
816 dVdT = (-1.)*(1./dPdV)*dPdT;
817 CPrs = cv + Tk*dPdT*dVdT - R_CONST;
818 Vrs = vmix;
819
820 return iRet;
821 }
822
823
824
825 #ifndef IPMGEMPLUGIN_
826
827 /// calculates properties of pure fluids when called from DCthermo
PRSVCalcFugPure(double Tmin,double * Cpg,double * FugProps)828 long int TPRSVcalc::PRSVCalcFugPure( double Tmin, double *Cpg, double *FugProps )
829 {
830 long int retCode = 0;
831 double Coeff[7];
832
833 for( int ii=0; ii<7; ii++ )
834 Coeff[ii] = Cpg[ii];
835
836 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
837 {
838 retCode = FugacityPT( 0, Coeff );
839 for( int i=0; i<6; i++ )
840 FugProps[i] = Fugpure[0][i];
841 return retCode;
842 }
843
844 else
845 {
846 for( int i=1; i<6; i++ )
847 FugProps[i] = 0.;
848 FugProps[0] = 1.;
849 FugProps[4] = 8.31451*Tk/Pbar;
850 return -1;
851 }
852 }
853
854 #endif
855
856
857
858
859
860 //=======================================================================================================
861 // Churakov-Gottschalk (CG) model for fluid mixtures
862 // References: Churakov and Gottschalk (2003a, 2003b)
863 // (c) SC June 2003
864 //=======================================================================================================
865
866
867 /// Generic constructor
TCGFcalc(long int NCmp,double Pp,double Tkp)868 TCGFcalc::TCGFcalc( long int NCmp, double Pp, double Tkp ):
869 TSolMod( NCmp, 'F', Tkp, Pp )
870 {
871 // Pparc = 0;
872 phWGT = 0;
873 aX = 0;
874 // aGEX = 0;
875 // aVol = 0;
876
877 set_internal();
878 alloc_internal();
879 }
880
881
882
TCGFcalc(SolutionData * sd,double * arphWGT,double * arX)883 TCGFcalc::TCGFcalc( SolutionData *sd, double *arphWGT, double *arX ):
884 TSolMod( sd )
885 {
886 Pparc = aPparc;
887 phWGT = arphWGT;
888 aX = arX;
889 // aGEX = arGEX;
890 // aVol = arVol;
891 set_internal();
892 alloc_internal();
893 }
894
895
896
897 // destructor
~TCGFcalc()898 TCGFcalc::~TCGFcalc()
899 {
900 free_internal();
901 }
902
903
904
905 /// set internally used parameters
set_internal()906 void TCGFcalc::set_internal()
907 {
908 PI_1 = 3.141592653589793120; // pi
909 TWOPI = 6.283185307179586230; // 2.*pi
910 PISIX = 0.523598775598298927; // pi/6.
911 TWOPOW1SIX = 1.12246204830937302; // 2^ = 1/6)
912 DELTA = 0.00001;
913 DELTAMOLLIM = 0.0000001;
914 R = 8.31439; // R constant
915 NA = 0.6023;
916 P1 = 1.186892378996;
917 PP2 = -0.4721963005527;
918 P3 = 3.259515855283;
919 P4 = 3.055229342609;
920 P5 = 1.095409321023;
921 P6 = 1.282306659774E-2;
922 P7 = 9.55712461425E-2;
923 P8 = 13.67807693107;
924 P9 = 35.75464856619;
925 P10 = 16.04724381643;
926 AA1 = -0.120078459237;
927 AA2 = -.808712488307;
928 AA3 = .321543801337;
929 A4 = 1.16965477132;
930 A5 = -.410564939543;
931 A6 = -.516834310691;
932 BB1 = -2.18839961483;
933 BB2 = 1.59897428009;
934 BB3 = -.392578806128;
935 B4 = -.189396607904;
936 B5 = -.576898496254;
937 B6 = -0.0185167641359;
938 A00 = .9985937977069455;
939 A01 = .5079834224407451;
940 A10 = 1.021887697885469;
941 A11 = -5.136619463333883;
942 A12 = -5.196188074016755;
943 A21 = -6.049240839050804;
944 A22 = 18.67848155616692;
945 A23 = 20.10652684217768;
946 A31 = 9.896491419756988;
947 A32 = 14.6738380473899;
948 A33 = -77.44825116542995;
949 A34 = -4.82871082941229;
950 }
951
952
953
alloc_internal()954 void TCGFcalc::alloc_internal()
955 {
956 paar = 0;
957 paar1 = 0;
958 FugCoefs = 0;
959 EoSparam = 0;
960 EoSparam1 = 0;
961 Cf = new double [NComp][8];
962 }
963
964
965
free_internal()966 void TCGFcalc::free_internal()
967 {
968 if( paar ) delete paar;
969 paar = 0;
970 if( paar1 ) delete paar1;
971 paar1 = 0;
972 if( FugCoefs ) delete[]FugCoefs;
973 if( EoSparam ) delete[]EoSparam;
974 if( EoSparam1 ) delete[]EoSparam1;
975 delete[]Cf;
976 }
977
978
979
980 /// high-level method to retrieve pure fluid fugacities
PureSpecies()981 long int TCGFcalc::PureSpecies()
982 {
983 double Fugacity = 0.1, Volume = 0.0;
984 double X[1] = {1.};
985 double Eos4parPT[4] = { 0.0, 0.0, 0.0, 0.0 },
986 Eos4parPT1[4] = { 0.0, 0.0, 0.0, 0.0 } ;
987 double roro; // added, 21.06.2008 (TW)
988 long int j, retCode = 0;
989
990 for( j=0; j<NComp; j++)
991 {
992 // Calling CG EoS for pure fugacity
993 if( Tk >= 273.15 && Tk < 1e4 && Pbar >= 1e-6 && Pbar < 1e5 )
994 retCode = CGFugacityPT( aDCc+j*NP_DC, Eos4parPT, Fugacity, Volume, Pbar, Tk, roro );
995 else {
996 Fugacity = Pbar;
997 Volume = 8.31451*Tk/Pbar;
998 Cf[j][0] = aDCc[j*NP_DC];
999 if( Cf[j][0] < 1. || Cf[j][0] > 10. )
1000 Cf[j][0] = 1.; // foolproof temporary
1001 Cf[j][1] = aDCc[j*NP_DC+1];
1002 Cf[j][2] = aDCc[j*NP_DC+2];
1003 Cf[j][3] = aDCc[j*NP_DC+3];
1004 Cf[j][4] = 0.;
1005 Cf[j][5] = 0.;
1006 Cf[j][6] = 0.;
1007 Cf[j][7] = 0.;
1008 continue;
1009 }
1010
1011 aGEX[j] = log( Fugacity / Pbar ); // now here (since 26.02.2008) DK
1012 Pparc[j] = Fugacity; // fure fluid fugacity (required for performance)
1013 aVol[j] = Volume*10.; // molar volume of pure fluid component, J/bar to cm3
1014
1015 // passing corrected EoS coeffs to calculation of fluid mixtures
1016 Cf[j][0] = Eos4parPT[0];
1017 if( Cf[j][0] < 1. || Cf[j][0] > 10. )
1018 Cf[j][0] = 1.; // foolproof temporary
1019 Cf[j][1] = Eos4parPT[1];
1020 Cf[j][2] = Eos4parPT[2];
1021 Cf[j][3] = Eos4parPT[3];
1022
1023 CGFugacityPT( aDCc+j*NP_DC, Eos4parPT1, Fugacity, Volume, Pbar, Tk+Tk*GetDELTA(), roro );
1024
1025 // passing corrected EoS coeffs for T+T*DELTA
1026 Cf[j][4] = Eos4parPT1[0];
1027 if( Cf[j][4] < 1. || Cf[j][4] > 10. )
1028 Cf[j][4] = 1.; // foolproof temporary
1029 Cf[j][5] = Eos4parPT1[1];
1030 Cf[j][6] = Eos4parPT1[2];
1031 Cf[j][7] = Eos4parPT1[3];
1032
1033 // Calculation of departure functions
1034 CGResidualFunct( X, Eos4parPT, Eos4parPT1, 1, roro, Tk ); // changed, 21.06.2008 (TW)
1035 } // j
1036
1037 if ( retCode )
1038 {
1039 char buf[150];
1040 sprintf(buf, "CG fluid: calculation of pure fugacity failed");
1041 Error( "E71IPM IPMgamma: ", buf );
1042 }
1043 return 0;
1044 }
1045
1046
1047
1048 /// calculates T,P corrected binary interaction parameters
PTparam()1049 long int TCGFcalc::PTparam()
1050 {
1051 long int i, j;
1052
1053 if( FugCoefs ) delete[]FugCoefs;
1054 if( EoSparam ) delete[]EoSparam;
1055 if( EoSparam1 ) delete[]EoSparam1;
1056
1057 FugCoefs = new double[ NComp ];
1058 EoSparam = new double[ NComp*4 ];
1059 EoSparam1 = new double[ NComp*4 ];
1060
1061 PureSpecies();
1062
1063 // Copying T,P corrected coefficients
1064 for( j=0; j<NComp; j++)
1065 {
1066 for( i=0; i<4; i++)
1067 EoSparam[j*4+i] = Cf[j][i];
1068 for( i=0; i<4; i++)
1069 EoSparam1[j*4+i] = Cf[j][i+4];
1070 }
1071 return 0;
1072 }
1073
1074
1075
1076 /// high-level method to retrieve activity coefficients in the fluid mixture
MixMod()1077 long int TCGFcalc::MixMod()
1078 {
1079 long int j;
1080 double roro; // changed, 21.06.2008 (TW)
1081
1082 if( Tk >= 273.15 && Tk < 1e4 && Pbar >= 1e-6 && Pbar < 1e5 )
1083 {
1084 CGActivCoefPT( aX, EoSparam, FugCoefs, NComp, Pbar, Tk, roro ); // changed, 21.06.2008 (TW)
1085 if (roro <= 0. )
1086 {
1087 char buf[150];
1088 sprintf(buf, "CG fluid: bad calculation of density ro= %lg", roro);
1089 Error( "E71IPM IPMgamma: ", buf );
1090 }
1091
1092 // Phase volume of the fluid in cm3 (not needed any more?)
1093 phVOL[0] = phWGT[0] / roro;
1094
1095 }
1096
1097 else // Setting Fugcoefs to 0 outside TP interval
1098 for( j=0; j<NComp; j++ )
1099 FugCoefs[ j ] = 0.0;
1100
1101 for( j=0; j<NComp; j++ )
1102 {
1103 if( FugCoefs[j] > 1e-23 )
1104 lnGamma[j] = log(FugCoefs[j]/Pparc[j]);
1105 else
1106 lnGamma[j] = 0;
1107 } // j
1108 return 0;
1109 }
1110
1111
1112
1113 /// high-level method to calculate residual functions
ExcessProp(double * Zex)1114 long int TCGFcalc::ExcessProp( double *Zex )
1115 {
1116 double roro; // changed, 21.06.2008 (TW)
1117
1118 if( Tk >= 273.15 && Tk < 1e4 && Pbar >= 1e-6 && Pbar < 1e5 )
1119 {
1120 CGActivCoefPT( aX, EoSparam, FugCoefs, NComp, Pbar, Tk, roro ); // changed, 21.06.2008 (TW)
1121 if (roro <= 0. )
1122 {
1123 char buf[150];
1124 sprintf(buf, "CG fluid: bad calculation of density ro= %lg", roro);
1125 Error( "E71IPM IPMgamma: ", buf );
1126 }
1127
1128 // calculate residual functions
1129 CGResidualFunct( aX, EoSparam, EoSparam1, NComp, roro, Tk );
1130
1131 }
1132
1133 else // setting residual functions to 0 outside TP interval
1134 {
1135 Grs = 0.;
1136 Srs = 0.;
1137 Hrs = 0.;
1138 CPrs = 0.;
1139 Vrs = 0.;
1140 }
1141
1142 Ars = Grs - Vrs*Pbar;
1143 Urs = Hrs - Vrs*Pbar;
1144
1145 // assignments (departure functions)
1146 Zex[0] = Grs;
1147 Zex[1] = Hrs;
1148 Zex[2] = Srs;
1149 Zex[3] = CPrs;
1150 Zex[4] = Vrs;
1151 Zex[5] = Ars;
1152 Zex[6] = Urs;
1153
1154 return 0;
1155 }
1156
1157
1158
1159 /// calculates ideal mixing properties
IdealProp(double * Zid)1160 long int TCGFcalc::IdealProp( double *Zid )
1161 {
1162 long int j;
1163 double s, sc, sp;
1164
1165 s = 0.0;
1166 for (j=0; j<NComp; j++)
1167 {
1168 if ( x[j] > 1.0e-32 )
1169 s += x[j]*log(x[j]);
1170 }
1171 sc = (-1.)*R_CONST*s;
1172 sp = (-1.)*R_CONST*log(Pbar);
1173 Hid = 0.0;
1174 CPid = 0.0;
1175 Vid = 0.0;
1176 Sid = sc + sp;
1177 Gid = Hid - Sid*Tk;
1178 Aid = Gid - Vid*Pbar;
1179 Uid = Hid - Vid*Pbar;
1180
1181 // assignments (ideal mixing properties)
1182 Zid[0] = Gid;
1183 Zid[1] = Hid;
1184 Zid[2] = Sid;
1185 Zid[3] = CPid;
1186 Zid[4] = Vid;
1187 Zid[5] = Aid;
1188 Zid[6] = Uid;
1189
1190 return 0;
1191 }
1192
1193
1194
1195 /// high-level method to retrieve pure fluid properties
CGFugacityPT(double * EoSparam,double * EoSparPT,double & Fugacity,double & Volume,double P,double T,double & roro)1196 long int TCGFcalc::CGFugacityPT( double *EoSparam, double *EoSparPT, double &Fugacity,
1197 double &Volume, double P, double T, double &roro )
1198 {
1199 long int iRet = 0;
1200 // double ro;
1201 double X[1] = {1.};
1202 double FugPure[1];
1203
1204 // modification to simplify CG database structure, 20.03.2007 (TW)
1205 EoSparPT[0] = EoSparam[0] + EoSparam[4]*exp(T*EoSparam[5]);
1206 EoSparPT[1] = EoSparam[1] + EoSparam[6]*exp(T*EoSparam[7]);
1207 EoSparPT[2] = EoSparam[2] + EoSparam[8]/(T+EoSparam[9]);
1208 EoSparPT[3] = EoSparam[3] + EoSparam[10]/(T+EoSparam[11]);
1209
1210 // returns density
1211 CGActivCoefPT( X, EoSparPT, FugPure, 1, P, T, roro ); // changed, 21.06.2008 (TW)
1212 if( roro < 0. )
1213 {
1214 return -1;
1215 };
1216
1217 Fugacity = FugPure[0];
1218 roro = DENSITY( X, EoSparPT, 1, P, T );
1219
1220 if( roro < 0 )
1221 { // error - density could not be calculated
1222 iRet = -2;
1223 roro = 1.0;
1224 }
1225
1226 Volume = 0.1/roro; // in J/bar
1227 // roro = ro; // added, 21.06.2008 (TW)
1228
1229 return iRet;
1230 }
1231
1232
1233
CGActivCoefPT(double * X,double * param,double * act,unsigned long int NN,double Pbar,double T,double & roro)1234 long int TCGFcalc::CGActivCoefPT( double *X,double *param, double *act,
1235 unsigned long int NN, double Pbar, double T, double &roro )
1236 {
1237 double *xtmp,*Fx;
1238 double P = Pbar/10.;
1239 xtmp = new double [NN];
1240 Fx = new double [NN];
1241
1242 if(!paar)
1243 paar = new EOSPARAM(X, param, NN);
1244 else
1245 paar->init( X, param, NN );
1246
1247 double F0,Z,F1,fideal;
1248 double ro,delta = DELTA,ax,dx /*,tmp*/;
1249 long int i;
1250
1251 norm(paar->XX0,paar->NCmp());
1252 copy(paar->XX0,xtmp,paar->NCmp());
1253
1254 paar->ParamMix(xtmp);
1255
1256 ro = ROTOTALMIX(P,T,paar);
1257
1258 if( ro < 0.0 ) // Too low pressure, no corrections will be done
1259 return ( -1 );
1260
1261 Z = P/(R*T*ro);
1262 F0 = FTOTALMIX(T,ro,paar);
1263
1264 // fideal=log(R*T*ro/BARMPA);
1265 fideal = log(R*T*ro/0.1);
1266 ax = Z - 1.+fideal;
1267
1268 for ( i=0;i<paar->NCmp();i++)
1269 {
1270 if ( xtmp[i]>0. )
1271 {
1272 copy(paar->XX0,xtmp,paar->NCmp());
1273 dx = xtmp[i]*delta;
1274 xtmp[i] += dx;
1275 norm(xtmp,paar->NCmp());
1276 paar->ParamMix(xtmp);
1277 F1 = FTOTALMIX(T,ro,paar)*(1.+dx);
1278 Fx[i] = (F1-F0)/(dx);
1279 }
1280 else Fx[i] = 0.;
1281 };
1282
1283 // GMix=0.;
1284 for ( i=0;i<paar->NCmp();i++)
1285 {
1286 if ( xtmp[i]>0. && Fx[i]< 100. )
1287 {
1288 // tmp=log(paar.XX0[i]);
1289 // GMix+=tmp*paar.XX0[i];
1290 act[i] = exp(ax+Fx[i]);
1291 }
1292 else
1293 {
1294 act[i] = 0.;
1295 }
1296 };
1297
1298 // GMix+=F0 + ax;
1299 // MLPutRealList(stdlink,act,paar.NCmp());
1300 delete[]xtmp;
1301 delete[]Fx;
1302 roro = ro; // added, 21.06.2008 (TW)
1303
1304 return 0; // changed, 21.06.2008 (TW)
1305 }
1306
1307
1308
1309 /// calculate residual functions through numerical derivative
CGResidualFunct(double * X,double * param,double * param1,unsigned long int NN,double ro,double T)1310 long int TCGFcalc::CGResidualFunct( double *X, double *param, double *param1, unsigned long int NN,
1311 double ro, double T )
1312 {
1313 double F0, Z, F1, vmix;
1314 double delta = DELTA;
1315 double *xtmp = new double [NN];
1316
1317 if(!paar)
1318 paar = new EOSPARAM(X, param, NN);
1319 else
1320 paar->init( X, param, NN );
1321
1322 if(!paar1)
1323 paar1 = new EOSPARAM(X, param1, NN);
1324 else
1325 paar1->init( X, param1, NN );
1326
1327 norm(paar->XX0,paar->NCmp());
1328 norm(paar1->XX0,paar1->NCmp());
1329 copy(paar->XX0,xtmp,paar->NCmp());
1330 paar->ParamMix(xtmp);
1331 paar1->ParamMix(xtmp);
1332 Z = ZTOTALMIX(T,ro,paar);
1333
1334 F0 = FTOTALMIX(T,ro,paar);
1335 // recalculate param1 for T+T*delta
1336 F1 = FTOTALMIX(T+T*delta,ro,paar1);
1337 // F1 = FTOTALMIX(T+T*delta,ro,paar);
1338
1339 Srs = - ( (F1-F0)/(delta*Tk)*Tk + F0 ) * R_CONST; // corrected, 20.06.2008 (TW)
1340 Hrs = (F0*Tk*R_CONST + Tk*Srs) + Z*R_CONST*Tk;
1341 Grs = Hrs - Tk*Srs;
1342 CPrs = 0.;
1343 vmix = Z*R_CONST*Tk/Pbar;
1344 Vrs = vmix;
1345
1346 delete[]xtmp;
1347 return 0;
1348
1349 }
1350
1351
1352
1353 // void ACTDENS(double *data,long nn, double *act )
CGActivCoefRhoT(double * X,double * param,double * act,unsigned long int NN,double ro,double T)1354 long int TCGFcalc::CGActivCoefRhoT( double *X, double *param, double *act,
1355 unsigned long int NN, double ro, double T )
1356 {
1357 double F0,Z,F1,GMix,fideal;
1358 double delta = DELTA,ax,dx,tmp;
1359 long int i;
1360 double *Fx,*xtmp;
1361 xtmp = new double [NN];
1362 Fx = new double [NN];
1363
1364 if(!paar)
1365 paar = new EOSPARAM(X, param, NN);
1366 else
1367 paar->init( X, param, NN );
1368
1369 norm(paar->XX0,paar->NCmp());
1370 copy(paar->XX0,xtmp,paar->NCmp());
1371 paar->ParamMix(xtmp);
1372 Z = ZTOTALMIX(T,ro,paar);
1373 F0 = FTOTALMIX(T,ro,paar);
1374 fideal = log(R*T*ro/0.1);
1375 ax = Z - 1.+fideal;
1376
1377 for ( i=0;i<paar->NCmp();i++)
1378 {
1379 if ( xtmp[i]>0. )
1380 {
1381 copy(paar->XX0,xtmp,NN);
1382 if ( xtmp[i]>DELTAMOLLIM )
1383 {
1384 dx = xtmp[i]*delta;
1385 }
1386 else
1387 {
1388 dx = DELTAMOLLIM*delta;
1389 }
1390
1391 xtmp[i] += dx;
1392 norm(xtmp,paar->NCmp());
1393 paar->ParamMix(xtmp);
1394 F1 = FTOTALMIX(T,ro,paar)*(1.+dx);
1395 Fx[i] = (F1-F0)/(dx);
1396
1397 }
1398 else Fx[i] = 0.;
1399 };
1400
1401 GMix = 0.;
1402 for ( i=0;i<paar->NCmp();i++)
1403 {
1404 if ( xtmp[i]>0. )
1405 {
1406 tmp = log(paar->XX0[i]);
1407 GMix += tmp*paar->XX0[i];
1408 act[i] = exp(ax+Fx[i]);
1409 }
1410 else
1411 {
1412 act[i] = 0.;
1413 }
1414 };
1415
1416 delete[]xtmp;
1417 delete[]Fx;
1418 return 0;
1419 // MLPutRealList(stdlink,act,paar.NCmp());
1420 }
1421
1422
1423
DIntegral(double T,double ro,unsigned long int IType)1424 double TCGFcalc::DIntegral( double T, double ro, unsigned long int IType )
1425 {
1426 static double TOld,roOld;
1427 static double a,b,c,d,e;
1428 static double data[][6]=
1429 {{-0.257431, 0.439229, 0.414783, -0.457019, -0.145520, 0.299666},
1430 {-0.396724, 0.690721, 0.628935, -0.652622, -0.201462, -0.23163 },
1431 {-0.488498, 0.863195, 0.761344, -0.750086, -0.218562, -0.538463},
1432 {-0.556600, 0.995172, 0.852903, -0.804710, -0.214736, -0.761700},
1433 {-0.611295, 1.103390, 0.921359, -0.838804, -0.197999, -0.940714},
1434 {-0.657866, 1.196189, 0.975721, -0.862346, -0.172526, -1.091678},
1435 {-0.698790, 1.278054, 1.020604, -0.880027, -0.140749, -1.222733},
1436 {-0.735855, 1.351533, 1.058986, -0.894024, -0.104174, -1.338626},
1437 {-0.769504, 1.418223, 1.092052, -0.905347, -0.063730, -1.442391},
1438 {-0.800934, 1.479538, 1.121453, -0.914864, -0.020150, -1.536070},
1439 {-0.829779, 1.535822, 1.147161, -0.922381, 0.026157 , -1.621183},
1440 {-0.856655, 1.587957, 1.169885, -0.928269, 0.074849 , -1.698853},
1441 {-0.881757, 1.636402, 1.190082, -0.932668, 0.125590 , -1.769898},
1442 {-0.904998, 1.681421, 1.207610, -0.935419, 0.178283 , -1.835070},
1443 {-0.926828, 1.723393, 1.223088, -0.936667, 0.232649 , -1.894899},
1444 {-0.946773, 1.762571, 1.236007, -0.936403, 0.288687 , -1.949858},
1445 {-0.965248, 1.799170, 1.246887, -0.934650, 0.346207 , -2.000344}};
1446
1447 // static double dt12[]=
1448 // {-2.139734,1.971553, 0.945513, -1.901492,-0.588630,-5.390941};
1449 // {-0.637684, 0.708107, 0.222086, -0.481116, -0.332141, -3.492213};
1450
1451 unsigned long int n;
1452 double *dtmp,rez;
1453
1454 if ( (T!=TOld) || (ro!=roOld) )
1455 {
1456 TOld = T;
1457 roOld = ro;
1458 e = log(T);
1459 b = ro*ro;
1460 d = ro;
1461 c = ro*e;
1462 a = b*e;
1463 }
1464
1465 // special case
1466 /*
1467 if ( IType==12 )
1468 {
1469 rez=(dt12[0]*T + dt12[1])*b +
1470 (dt12[2]*T + dt12[3])*ro + dt12[4]*T + dt12[5];
1471 return exp(rez);
1472 }
1473 */
1474
1475 n = IType-4;
1476 dtmp = data[n];
1477 rez = dtmp[0]*a + dtmp[1]*b + dtmp[2]*c + dtmp[3]*d + dtmp[4]*e + dtmp[5];
1478 return exp(rez);
1479 }
1480
1481
1482
LIntegral(double T,double ro,unsigned long int IType)1483 double TCGFcalc::LIntegral( double T, double ro,unsigned long int IType )
1484 {
1485 static double TOld,roOld;
1486 static double a,b,c,d,e;
1487 static double data[][6]=
1488 {{ -1.010391, 1.628552, 2.077476, -2.30162 , -0.689931, -2.688117},
1489 { -1.228611, 2.060090, 2.463396, -2.453303, -0.573894, -3.350638},
1490 { -1.354004, 2.402034, 2.718124, -2.462814, -0.412252, -4.018632}};
1491
1492 double *dtmp,rez;
1493
1494 if ( (T!=TOld) || (ro!=roOld) )
1495 {
1496 TOld = T;
1497 roOld = ro;
1498 a = ro*ro*log(T);
1499 b = ro*ro;
1500 c = ro*log(T);
1501 d = ro;
1502 e = log(T);
1503 }
1504
1505 switch ( IType )
1506 {
1507 case 662:
1508 dtmp = data[0];
1509 break;
1510 case 1262:
1511 dtmp = data[1];
1512 break;
1513 case 12122:
1514 dtmp = data[2];
1515 break;
1516 default:
1517 return 0;
1518 }
1519 rez = dtmp[0]*a + dtmp[1]*b + dtmp[2]*c + dtmp[3]*d + dtmp[4]*e + dtmp[5];
1520 return -exp(rez);
1521 }
1522
1523
1524
KIntegral(double T,double ro,unsigned long int IType)1525 double TCGFcalc::KIntegral( double T, double ro,unsigned long int IType )
1526 {
1527 static double TOld,roOld;
1528 static double a,b,c,d,e;
1529 static double data[][6]=
1530 {{ -1.050534, 1.747476, 1.749366, -1.999227, -0.661046, -3.028720},
1531 { -1.309550, 2.249120, 2.135877, -2.278530, -0.773166, -3.704690},
1532 { -1.490116, 2.619997, 2.404319, -2.420706, -0.829466, -3.930928},
1533 { -1.616385, 2.881007, 2.577600, -2.484990, -0.828596, -4.175589},
1534 { -1.940503, 3.552034, 2.940925, -2.593808, -0.724353, -4.899975}};
1535
1536 double *dtmp,rez;
1537
1538 if ( (T!=TOld) || (ro!=roOld) )
1539 {
1540 TOld = T;
1541 roOld = ro;
1542 a = ro*ro*log(T);
1543 b = ro*ro;
1544 c = ro*log(T);
1545 d = ro;
1546 e = log(T);
1547 }
1548
1549 switch ( IType )
1550 {
1551 case 222333:
1552 dtmp = data[0];
1553 break;
1554 case 233344:
1555 dtmp = data[1];
1556 break;
1557 case 334445:
1558 dtmp = data[2];
1559 rez = dtmp[0]*a + dtmp[1]*b + dtmp[2]*c + dtmp[3]*d + dtmp[4]*e + dtmp[5];
1560 return -exp(rez);
1561 case 444555:
1562 dtmp = data[3];
1563 break;
1564 case 666777:
1565 dtmp = data[4];
1566 break;
1567 default:
1568 return 0;
1569
1570 }
1571
1572 rez = dtmp[0]*a + dtmp[1]*b + dtmp[2]*c + dtmp[3]*d + dtmp[4]*e + dtmp[5];
1573 return exp(rez);
1574 }
1575
1576
1577
K23_13(double T,double ro)1578 double TCGFcalc::K23_13( double T, double ro )
1579 {
1580 static double TOld,roOld,KOLD;
1581 static double a,b,c,d,e;
1582 static double dtmp[]=
1583 { -1.050534, 1.747476, 1.749366, -1.999227, -0.661046, -3.028720};
1584
1585 if ( (T!=TOld) || (ro!=roOld) )
1586 {
1587 TOld = T;
1588 roOld = ro;
1589 a = ro*ro*log(T);
1590 b = ro*ro;
1591 c = ro*log(T);
1592 d = ro;
1593 e = log(T);
1594 }
1595 else return KOLD;
1596
1597 KOLD = dtmp[0]*a + dtmp[1]*b + dtmp[2]*c + dtmp[3]*d + dtmp[4]*e + dtmp[5];
1598 KOLD = exp(KOLD/3.);
1599 return KOLD;
1600 }
1601
1602
1603
DENSITY(double * X,double * param,unsigned long NN,double Pbar,double T)1604 double TCGFcalc::DENSITY( double *X, double *param, unsigned long NN, double Pbar, double T )
1605 {
1606 double P = Pbar * 0.1;
1607 double *xtmp;
1608 double ro;
1609
1610 xtmp = new double [NN];
1611 if( !paar1 )
1612 paar1 = new EOSPARAM(X,param,NN);
1613 else
1614 paar1->init( X, param, NN );
1615
1616 norm(paar1->XX0,paar1->NCmp());
1617 copy(paar1->XX0,xtmp,paar1->NCmp());
1618 paar1->ParamMix(xtmp);
1619 ro = ROTOTALMIX(P,T,paar1);
1620
1621 delete [] xtmp;
1622 if( ro < 0. )
1623 Error( ""," Error - density cannot be found at this T,P" );
1624 return ro;
1625 }
1626
1627
1628
PRESSURE(double * X,double * param,unsigned long int NN,double ro,double T)1629 double TCGFcalc::PRESSURE( double *X,double *param,
1630 unsigned long int NN,double ro, double T )
1631 {
1632 double *xtmp;
1633 xtmp = new double [NN];
1634
1635 if( !paar1 )
1636 paar1 = new EOSPARAM(X,param,NN);
1637 else
1638 paar1->init( X, param, NN );
1639
1640 norm(paar1->XX0,paar1->NCmp());
1641 copy(paar1->XX0,xtmp,paar1->NCmp());
1642 paar1->ParamMix(xtmp);
1643 double P = PTOTALMIX(T,ro,paar1);
1644 delete [] xtmp;
1645 return P*10.;
1646 }
1647
1648
1649
copy(double * sours,double * dest,unsigned long int num)1650 void TCGFcalc::copy( double* sours,double *dest,unsigned long int num )
1651 {
1652 unsigned long int i;
1653 for ( i=0; i<num; i++)
1654 {
1655 dest[i]=sours[i];
1656 };
1657 }
1658
1659
1660
norm(double * X,unsigned long int mNum)1661 void TCGFcalc::norm( double *X,unsigned long int mNum )
1662 {
1663 double tmp=0.;
1664 unsigned long int i;
1665 for ( i=0; i<mNum; i++ )
1666 {
1667 tmp += X[i];
1668 }
1669 tmp = 1./tmp;
1670 for ( i=0; i<mNum; i++ )
1671 {
1672 X[i] *= tmp;
1673 }
1674 }
1675
1676
1677
RPA(double beta,double nuw)1678 double TCGFcalc::RPA( double beta,double nuw )
1679 {
1680 double fi1,fi2;
1681 fi1 = (1.20110+(0.064890+(-76.860+(562.686+(-2280.090+(6266.840+(-11753.40+(14053.8
1682 +(-9491.490 +2731.030*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw;
1683 fi2 = (0.588890+(-7.455360+(40.57590+(-104.8970+(60.25470+(390.6310+(-1193.080
1684 +(1576.350+(-1045.910+283.7580*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw)*nuw*nuw;
1685 return (-12.*fi1 + 192.*fi2*beta)*beta*beta/PI_1;
1686 }
1687
1688
1689
dHS(double beta,double ro)1690 double TCGFcalc::dHS( double beta,double ro )
1691 {
1692 // service constants
1693 double DV112 = 1./12.;
1694 double DV712 = 7./12.;
1695 // local variables
1696 double T12, T112, T712, B13, dB, delta, d;
1697 double a0, a1, a6, a3, a4, a7, a9, a12;
1698 double p0, p2, p6, p3, p5, p8, p11;
1699 double dbdl, ri6ro, ri6ro2, d3, d2, dnew, F0, F1;
1700 unsigned long int i;
1701
1702 T12 = sqrt(beta);
1703 T112 = exp(DV112*log(beta));
1704 T712 = exp(DV712*log(beta));
1705 B13 = (1+beta);
1706 B13 = B13*B13*B13;
1707
1708 dB = (P1*T112+PP2*T712+(P3+(P4+P5*beta)*beta)*beta)/B13;
1709 delta = (P6+P7*T12)/(1.+(P8+(P9+P10*T12)*T12)*T12);
1710
1711 dbdl = dB*delta;
1712 ri6ro = PISIX*ro;
1713 ri6ro2 = ri6ro*ri6ro;
1714
1715 a0 = dB+dbdl;
1716 a1 = -1.;
1717 a3 = (-1.5*dB -3.75*dbdl)*ri6ro;
1718 a4 = (1.5*ri6ro);
1719 a6 = (2.*dB + dbdl)*0.25*ri6ro2;
1720 a7 = -0.5*ri6ro2;
1721 a9 = -2.89325*ri6ro2*ri6ro*dbdl;
1722 a12 = -0.755*ri6ro2*ri6ro2*dbdl;
1723
1724 p0 = -1.;
1725 p2 = a3*3.;
1726 p3 = a4*4.;
1727 p5 = a6*6.;
1728 p6 = a7*7.;
1729 p8 = a9*9.;
1730 p11 = a12*12.;
1731
1732 d = dB;
1733 i = 0;
1734
1735 while ( i++<21 )
1736 {
1737 d2 = d*d;
1738 d3 = d*d*d;
1739 F0 = a0+(a1+(a3+(a4+(a6+(a7+(a9+a12*d3)*d2)*d)*d2)*d)*d2)*d;
1740 F1 = p0+(p2+(p3+(p5+(p6+(p8+p11*d3)*d2)*d)*d2)*d)*d2;
1741 dnew = d-F0/F1;
1742 if ( fabs(dnew-d)<1.E-7 )
1743 {
1744 return dnew;
1745 }
1746 d = dnew;
1747 }
1748
1749 if ( i>=20 )
1750 {
1751 return dB;
1752 }
1753 return dnew;
1754 }
1755
1756
1757
FWCA(double T,double ro)1758 double TCGFcalc::FWCA( double T,double ro )
1759 {
1760 static double TOld,roOld,F;
1761 double d,beta,nu,nuw;
1762 double nu1w1,nu1w2,nu1w3,nu1w4,nu1w5;
1763 double a0,a1,a2,a3;
1764 double I2;
1765 double I1_6,I1_12;
1766 double dW,dW12,dW6;
1767 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7;
1768 double F0,F1,FA;
1769 double rm,rmdw1,rmdw2,rmdw3,rmdw4,rmdw5;
1770
1771 if ((T==TOld) && (ro==roOld))
1772 {
1773 return F;
1774 }
1775 else
1776 {
1777 TOld = T;
1778 roOld = ro;
1779 }
1780
1781 rm = TWOPOW1SIX;
1782 beta = 1./T;
1783 d = dHS( beta, ro );
1784 tmp2 = PISIX*d*d*d;
1785 nu = tmp2*ro;
1786 tmp1 = (1. - nu/16.);
1787 nuw = nu*tmp1;
1788 dW = d*exp(1./3.*log(tmp1));
1789
1790 nu1w1 = (1.-nuw);
1791 nu1w2 = nu1w1*nu1w1;
1792 nu1w3 = nu1w2*nu1w1;
1793 nu1w4 = nu1w2*nu1w2;
1794 nu1w5 = nu1w2*nu1w3;
1795
1796 tmp1 = (1-nu);
1797 tmp1 = tmp1*tmp1;
1798 F0 = ((4.-3.*nu)*nu)/tmp1;
1799
1800 a0 = fa0( nuw , nu1w2);
1801 a1 = fa1( nuw , nu1w3);
1802 a2 = fa2( nuw , nu1w4);
1803 a3 = fa3( nuw , nu1w5);
1804
1805 I1_6 = fI1_6( nuw );
1806 I1_12 = fI1_12( nuw );
1807
1808 rmdw1 = rm/dW;
1809 rmdw2 = rmdw1*rmdw1;
1810 rmdw3 = rmdw1*rmdw2;
1811 rmdw4 = rmdw2*rmdw2;
1812 rmdw5 = rmdw3*rmdw2;
1813
1814 dW6 = dW*dW*dW;
1815 dW6 = 1./(dW6*dW6);
1816 dW12 = dW6*dW6;
1817
1818 tmp1 = (a0/4.+ a1/12. + a2/24. + a3/24.)*dW6;
1819 tmp2 = (a0/10.+ a1/90. + a2/720. + a3/5040.)*(-dW12);
1820 tmp3 = (a0 - a1/3. + a2/12 - a3/60)/8.;
1821 tmp4 = (a0 - a1 + a2/2. - a3/6.)*rmdw2*(-9.)/40.;
1822 tmp5 = (a1 - a2 + a3/2)*rmdw3*(-2.)/9.;
1823 tmp6 = (a2 - a3)*rmdw4*(-9.)/64.;
1824 tmp7 = a3*(-3.)/35.*rmdw5;
1825
1826 I2 = tmp1+tmp2+tmp3+tmp4+tmp5+tmp6+tmp7;
1827
1828 F1 = 48.*nuw*(I1_12*dW12-I1_6*dW6 + I2)*beta;
1829 FA = RPA(beta,nuw);
1830
1831 F = F0+F1+FA;
1832
1833 return F;
1834 }
1835
1836
1837
ZWCANum(double T,double ro)1838 double TCGFcalc::ZWCANum( double T,double ro )
1839 {
1840 double delta = DELTA;
1841 double a0,a1;
1842 a1 = FWCA(T,ro*(1.+delta));
1843 a0 = FWCA(T,ro);
1844 return 1.+(a1-a0)/delta;
1845 }
1846
1847
1848
UWCANum(double T,double ro)1849 double TCGFcalc::UWCANum( double T,double ro )
1850 {
1851 double delta = DELTA;
1852 double a0,a1,beta0,beta1;
1853 beta0 = 1./T;
1854 beta1 = beta0*(1.+delta);
1855 a1 = FWCA(1./beta1,ro);
1856 a0 = FWCA(T,ro);
1857 return (a1-a0)/(beta1-beta0);
1858 }
1859
1860
1861
FDipPair(double T,double ro,double m2)1862 double TCGFcalc::FDipPair( double T,double ro,double m2 )
1863 {
1864 double kappa,Z,U,beta,F;
1865 kappa = m2*m2/(24.*T);
1866 beta = 1./T;
1867 Z = ZWCANum(T,ro);
1868 U = UWCANum(T,ro);
1869 F = kappa*(4.*beta*U-Z+1.);
1870 return F;
1871 }
1872
1873
1874
J6LJ(double T,double ro)1875 double TCGFcalc::J6LJ( double T,double ro )
1876 {
1877 double kappa,Z,U,beta,F;
1878 beta = 1./T;
1879 Z = ZWCANum(T,ro);
1880 kappa = -16.*PI_1*ro*beta;
1881 U = UWCANum(T,ro);
1882 F = (4.*beta*U-Z+1.)/kappa;
1883 return F;
1884 }
1885
1886
1887
FTOTALMIX(double T_Real,double ro_Real,EOSPARAM * param)1888 double TCGFcalc::FTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param )
1889 {
1890 double FF,A0,A2,A3,AP,A1;
1891 // unsigned iall,inopol;
1892 double emix,s3mix,rotmp,T2R;
1893 double Jind,Jdp;
1894 long int /*itmp,jtmp,ktmp,*/ i,j,k;
1895 double s3tmp,mtmp,IK /*,atmp*/;
1896 double imtmp,jmtmp,iatmp,jatmp;
1897 double m2i,m2j,m2k;
1898 double s3tmpij,s3tmpik,s3tmpjk;
1899 double IKtmpij,IKtmpik,IKtmpjk;
1900
1901 // iall=param.inonzero();
1902 // inopol=param.inonpolar();
1903 emix = param->EMIX();
1904 s3mix = param->S3MIX();
1905
1906 rotmp = NA*ro_Real;
1907 T2R = T_Real*T_Real;
1908
1909 A0 = FWCA(T_Real/emix,s3mix*rotmp);
1910 // if ( inopol< iall )
1911 {
1912 // dipole part
1913 A2 = 0.;
1914 for ( i=0; i<param->NCmp()-1; i++ )
1915 {
1916 for ( j=i+1; j<param->NCmp(); j++ )
1917 {
1918 s3tmp = param->MIXS3(i,j);
1919 Jdp = J6LJ(T_Real*s3tmp/param->MIXES3(i,j) , s3tmp*rotmp);
1920 A2 += param->M2R(i)*param->M2R(j)*Jdp*
1921 param->X(i)*param->X(j)/s3tmp;
1922 }
1923 }
1924 A2 *= 2.;
1925 for ( i=0; i<param->NCmp(); i++ )
1926 {
1927 // itmp=param->ind(i);
1928 mtmp = param->M2R(i);
1929 s3tmp = param->SIG3(i);
1930 Jdp = J6LJ(T_Real/param->EPS(i),s3tmp*rotmp);
1931 A2 += mtmp*mtmp*Jdp*param->X(i)*param->X(i)/s3tmp;
1932 }
1933
1934 A2 = -A2*TWOPI*rotmp/(3.*T2R);
1935 // A2 done
1936
1937 if ( A2!=0. )
1938 {
1939
1940 A3 = 0.;
1941 for ( i=0; i<param->NCmp(); i++ )
1942 {
1943 // itmp=param->ind(i);
1944 m2i = param->M2R(i);
1945
1946 for ( j=0; j<param->NCmp(); j++ )
1947 {
1948 // jtmp=param->ind(j);
1949 m2j = param->M2R(j);
1950
1951 s3tmpij = param->MIXS3(i,j);
1952 IKtmpij = K23_13(T_Real*s3tmpij/param->MIXES3(i,j),
1953 s3tmpij*rotmp);
1954 for ( k=0; k<param->NCmp(); k++ )
1955 {
1956 // ktmp=param->ind(k);
1957 m2k = param->M2R(k);
1958
1959 s3tmpik = param->MIXS3(i,k);
1960 s3tmpjk = param->MIXS3(j,k);
1961
1962 IKtmpik = K23_13(T_Real*s3tmpik/param->MIXES3(i,k),
1963 s3tmpik*rotmp);
1964 IKtmpjk = K23_13(T_Real*s3tmpjk/param->MIXES3(j,k),
1965 s3tmpjk*rotmp);
1966
1967 IK = IKtmpij*IKtmpik*IKtmpjk;
1968 A3 += m2i*m2j*m2k*IK*pow(s3tmpij*s3tmpik*s3tmpjk,-1./3.)*
1969 param->X(i)*param->X(j)*param->X(k);
1970 }
1971 }
1972 }
1973 A3 = A3*32.*sqrt(14.*PI_1/5.)*
1974 rotmp*rotmp*PI_1*PI_1*PI_1/(135.*T_Real*T2R);
1975 AP = A2/(1. - A3/A2);
1976 }
1977 else AP = 0.;
1978
1979 // induced interaction
1980 A1 = 0.;
1981 for ( i=0; i<param->NCmp(); i++ )
1982 {
1983 // itmp=param->ind(i);
1984 iatmp = param->A(i);
1985 imtmp = param->M2R(i);
1986 for ( j=0; j<param->NCmp(); j++ )
1987 {
1988 // jtmp=param->ind(j);
1989 jatmp = param->A(j);
1990 jmtmp = param->M2R(j);
1991
1992 s3tmp = param->MIXS3(i,j);
1993 Jind = J6LJ(T_Real*s3tmp/param->MIXES3(i,j),s3tmp*rotmp);
1994
1995 A1 += (iatmp*jmtmp + jatmp*imtmp)
1996 *Jind*param->X(i)*param->X(j)/s3tmp;
1997 }
1998 }
1999 A1 = -A1*TWOPI*rotmp/T_Real;
2000 // A1=-A1*FOURPI*rotmp/T_Real;
2001 // A1=0.;
2002
2003 } // end of polar contribution
2004
2005 FF = A0 + A1 + AP;
2006 //printf("%g %g %g %g %g",A0,A1,A2,A3,AP);
2007 //exit(1) ;
2008
2009 return FF;
2010 }
2011
2012
2013
UTOTALMIX(double T_Real,double ro_Real,EOSPARAM * param)2014 double TCGFcalc::UTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param )
2015 {
2016 double T /*,ro,s3 */;
2017 double delta = DELTA;
2018 double a0,a1,beta0,beta1,eps;
2019 eps = param->EMIX();
2020 T = T_Real/eps;
2021
2022 beta0 = 1./T;
2023 beta1 = beta0*(1.+delta);
2024 a1 = FTOTALMIX((1./beta1)*eps,ro_Real,param);
2025 a0 = FTOTALMIX(T_Real,ro_Real,param);
2026 return (a1-a0)/(beta1-beta0);
2027 }
2028
2029
2030
ZTOTALMIX(double T_Real,double ro_Real,EOSPARAM * param)2031 double TCGFcalc::ZTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param )
2032 {
2033 double delta = DELTA;
2034 double a0,a1;
2035 a1 = FTOTALMIX(T_Real,ro_Real*(1.+delta),param);
2036 a0 = FTOTALMIX(T_Real,ro_Real,param);
2037
2038 return 1.+(a1-a0)/delta;
2039 }
2040
2041
2042
PTOTALMIX(double T_Real,double ro_Real,EOSPARAM * param)2043 double TCGFcalc::PTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param )
2044 {
2045 double Z;
2046 Z = ZTOTALMIX(T_Real,ro_Real,param);
2047 return Z*R*T_Real*ro_Real;
2048 }
2049
2050
2051
2052 /// melting density
Melt(double T)2053 double TCGFcalc::Melt( double T )
2054 {
2055
2056 return T*0.+.9;
2057
2058 }
2059
2060
2061
Melt2(double T)2062 double TCGFcalc::Melt2(double T)
2063 {
2064 return T*0.+3.;
2065 }
2066
2067
2068 #define FIRSTSEED (15)
2069 #define ROMIN (1.E-2)
2070 #define NPOINT (5)
2071
2072
2073
choose(double * pres,double P,unsigned long int & x1,unsigned long int & x2)2074 void TCGFcalc::choose( double *pres, double P,unsigned long int &x1,unsigned long int &x2 )
2075 {
2076 unsigned long int i;
2077 double deltam = -10000000.,tmp;
2078 double deltap = 10000000.;
2079
2080 for ( i=0; i<NPOINT; i++ )
2081 {
2082 tmp = P-pres[i];
2083
2084 if ( tmp>0. )
2085 {
2086 if ( tmp<deltap )
2087 {
2088 deltap=tmp;
2089 x1 = i;
2090 }
2091 }
2092 else
2093 {
2094 if ( tmp>deltam )
2095 {
2096 deltam=tmp;
2097 x2 = i;
2098 }
2099 }
2100 }
2101 return ;
2102 }
2103
2104
2105
ROTOTALMIX(double P,double TT,EOSPARAM * param)2106 double TCGFcalc::ROTOTALMIX( double P,double TT,EOSPARAM* param )
2107 {
2108 unsigned long int i;
2109 double T /*,ro*/;
2110 double fact, fact0, romax, dro, roarr[FIRSTSEED];
2111 double Ptmp[FIRSTSEED], ro0, ro1, rotest, PP0, PP1 /* ,Ptest */;
2112 double a,b;
2113 double inttofloat;
2114 double f[4],x[4],ff,dens[5],pres[5];
2115 unsigned long int x1=0L,x2=0L;
2116 // double ptmp;
2117
2118 T = TT/param->EMIX();
2119 fact0 = 1./(param->S3MIX()*NA);
2120 fact = R*TT*fact0;
2121
2122 romax = Melt(T);
2123 inttofloat = FIRSTSEED-1;
2124 dro = (romax-ROMIN)/inttofloat;
2125 roarr[0] = ROMIN;
2126 roarr[1] = 2.*ROMIN;
2127
2128 for ( i=2; i<FIRSTSEED; i++)
2129 {
2130 inttofloat = i;
2131 roarr[i] = ROMIN+inttofloat*dro;
2132 }
2133
2134 for ( i=0; i<FIRSTSEED; i++)
2135 {
2136 Ptmp[i] = ZTOTALMIX(TT,roarr[i]*fact0,param);
2137 Ptmp[i] *= roarr[i] * fact;
2138 if ( Ptmp[i] > P )
2139 {
2140 break;
2141 }
2142 }
2143
2144 if ( i==0 ) // Uses aproximation of ideal gas
2145 {
2146 return P/(R*TT);
2147 }
2148
2149 // additional high pressure inteval
2150 if ( i==FIRSTSEED )
2151 {
2152
2153 // roarr[0]=romax-0.0001;
2154 roarr[0] = roarr[FIRSTSEED-1];
2155 Ptmp[0] = Ptmp[FIRSTSEED-1];
2156
2157 romax = Melt2(T);
2158 inttofloat = FIRSTSEED-1;
2159 dro = (romax-ROMIN)/inttofloat;
2160 for ( i=1; i<FIRSTSEED; i++)
2161 {
2162 inttofloat = i;
2163 roarr[i] = ROMIN+inttofloat*dro;
2164 }
2165
2166 for ( i=1; i<FIRSTSEED; i++)
2167 {
2168 Ptmp[i] = ZTOTALMIX(TT,roarr[i]*fact0,param)*roarr[i]*fact;
2169 if ( Ptmp[i]>P )
2170 {
2171 break;
2172 }
2173 }
2174
2175 if ( i==FIRSTSEED || i==0 )
2176 {
2177 printf( "Input pressure is too high!\n" );
2178 // exit(1);
2179 return (-1.0);
2180 }
2181 }
2182
2183 ro0 = roarr[i-1];
2184 ro1 = roarr[i];
2185 PP0 = Ptmp[i-1];
2186 PP1 = Ptmp[i];
2187 i = 0;
2188
2189 while ( i++<20 )
2190 {
2191 // Start interp
2192 ff = ro0;
2193 dens[0] = ro0;
2194 dens[1] = ro1;
2195 pres[0] = PP0;
2196 pres[1] = PP1;
2197
2198 // first order
2199 x[0] = P-pres[0];
2200 f[0] = (dens[1]-dens[0])/(pres[1]-pres[0]);
2201 ff += f[0]*x[0];
2202
2203 // second order
2204 dens[2] = ff;
2205 pres[2] = ZTOTALMIX(TT,ff*fact0,param)*ff*fact;
2206
2207 if ( fabs(pres[2]-P)<1E-5 )
2208 {
2209 return ff*fact0;
2210 }
2211
2212 x[1] = x[0]*(P-pres[1]);
2213 f[1] = (dens[2]-dens[1])/(pres[2]-pres[1]);
2214
2215 f[0] = (f[1]-f[0])/(pres[2]-pres[0]);
2216 ff += f[0]*x[1];
2217
2218 // third order
2219 dens[3] = ff;
2220 pres[3] = ZTOTALMIX(TT,ff*fact0,param)*ff*fact;
2221 if ( fabs(pres[3]-P)<1E-6 )
2222 {
2223 return ff*fact0;
2224 }
2225 x[2] = x[1]*(P-pres[2]);
2226 f[2] = (dens[3]-dens[2])/(pres[3]-pres[2]);
2227 f[1] = (f[2]-f[1])/(pres[3]-pres[1]);
2228 f[0] = (f[1]-f[0])/(pres[3]-pres[0]);
2229 ff += f[0]*x[2];
2230 dens[4] = ff;
2231 pres[4] = ZTOTALMIX(TT,ff*fact0,param)*ff*fact;
2232 if ( fabs(pres[4]-P)<1e-6 )
2233 {
2234 return ff*fact0;
2235 }
2236
2237 choose(pres,P,x1,x2);
2238
2239 ro0 = dens[x1];
2240 ro1 = dens[x2];
2241 PP0 = pres[x1];
2242 PP1 = pres[x2];
2243
2244 if ( fabs((ro1-ro0))<0.001 )
2245 {
2246 a = (PP1-PP0)/(ro1-ro0);
2247 b = PP1-a*ro1;
2248 rotest = (P-b)/a;
2249 return rotest*(fact0);
2250 }
2251 }
2252 //return 10.;
2253 /// bad result
2254
2255 a = (PP1-PP0)/(ro1-ro0);
2256 b = PP1-a*ro1;
2257 rotest = (P-b)/a;
2258 return rotest*(fact0);
2259
2260 }
2261
2262 #ifndef IPMGEMPLUGIN_
2263
2264
2265
2266 /// calculates properties of pure fluids when called from DCthermo
CGcalcFugPure(double Tmin,double * Cemp,double * FugProps)2267 long int TCGFcalc::CGcalcFugPure(double Tmin, double *Cemp, double *FugProps )
2268 {
2269 long int retCode = 0;
2270 double T, P, Fugacity = 0.1, Volume = 0.0;
2271 double X[1] = {1.};
2272 double roro = 1.; // added 21.06.2008 (TW)
2273 double Coeff[12]; // MAXEOSPARAM = 20;
2274 double Eos4parPT[4] = { 0.0, 0.0, 0.0, 0.0 },
2275 Eos4parPT1[4] = { 0.0, 0.0, 0.0, 0.0 } ;
2276
2277 T = Tk;
2278 P = Pbar;
2279
2280 for(long int ii=0; ii<12; ii++ )
2281 Coeff[ii] = (double)Cemp[ii];
2282
2283 // Calling CG EoS functions here
2284 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
2285 {
2286 retCode = CGFugacityPT( Coeff, Eos4parPT, Fugacity, Volume, P, T, roro );
2287 FugProps[0] = Fugacity/Pbar;
2288 FugProps[1] = 8.31451 * Tk * log( Fugacity / P );
2289 FugProps[4] = Volume;
2290 retCode = CGFugacityPT( Coeff, Eos4parPT1, Fugacity, Volume, P, T+T*DELTA, roro );
2291 CGResidualFunct( X, Eos4parPT, Eos4parPT1, 1, roro, T );
2292 FugProps[2] = Hrs;
2293 FugProps[3] = Srs;
2294 return retCode;
2295 }
2296
2297 else
2298 {
2299 for( int i=1; i<6; i++ )
2300 FugProps[i] = 0.;
2301 FugProps[0] = 1.;
2302 FugProps[4] = 8.31451*Tk/Pbar;
2303 return -1;
2304 }
2305 }
2306
2307 #endif
2308
2309
2310
2311
2312 //=======================================================================================================
2313 // Implementation of EOSPARAM class (used by TCGFcalc class)
2314 //=======================================================================================================
2315
2316
free()2317 void EOSPARAM::free()
2318 {
2319 long int i;
2320
2321 if ( NComp > 0)
2322 {
2323 for ( i=0;i<NComp;i++ )
2324 delete[]mixpar[i];
2325 delete[]mixpar;
2326
2327 delete[]epspar;
2328 delete[]sig3par;
2329 delete[]XX;
2330 delete[]eps;
2331 delete[]eps05;
2332 delete[]sigpar;
2333 delete[]mpar;
2334 delete[]apar;
2335 delete[]aredpar;
2336 delete[]m2par;
2337 delete[]XX0;
2338 NComp = 0;
2339 }
2340 }
2341
2342
2343
allocate()2344 void EOSPARAM::allocate()
2345 {
2346 long int i;
2347
2348 mixpar = new double*[NComp];
2349 for ( i=0; i<NComp; i++ )
2350 mixpar[i] = new double[NComp];
2351
2352 epspar = new double[NComp];
2353 sig3par = new double[NComp];
2354 XX = new double[NComp];
2355 eps = new double[NComp];
2356 eps05 = new double[NComp];
2357 sigpar = new double[NComp];
2358 mpar = new double[NComp];
2359 apar = new double[NComp];
2360 aredpar = new double[NComp];
2361 m2par = new double[NComp];
2362 XX0 = new double[NComp];
2363 }
2364
2365
2366
init(double * Xinp,double * data,long int nn)2367 void EOSPARAM::init( double *Xinp, double * data, long int nn )
2368 {
2369 long int i,j;
2370 double tmp;
2371
2372 if( nn != NComp )
2373 { // or error message
2374 free();
2375 NComp = nn;
2376 allocate();
2377 }
2378
2379 for ( i=0;i<NComp;i++ )
2380 {
2381 XX0[i] = Xinp[i];
2382
2383 sigpar[i] = data[i*4 ];
2384 eps[i] = data[i*4 + 1];
2385 mpar[i] = data[i*4 + 2];
2386 apar[i] = data[i*4 + 3];
2387 }
2388 for ( i=0; i<NComp; i++ )
2389 {
2390 tmp = sigpar[i];
2391 tmp = tmp*tmp*tmp;
2392 sig3par[i] = tmp;
2393 eps05[i] = sqrt(eps[i]);
2394 epspar[i] = tmp*eps[i];
2395 m2par[i] = mpar[i]*mpar[i]/(1.38048E-4);
2396 aredpar[i] = apar[i]/tmp;
2397 }
2398
2399 // calculation of mixing properties
2400 for ( i=0; i<NComp-1; i++ )
2401 {
2402 for ( j=i+1; j<NComp; j++ )
2403 {
2404 tmp = (sigpar[i]+sigpar[j])*0.5;
2405 tmp = tmp*tmp*tmp;
2406 mixpar[i][j] = tmp;
2407 mixpar[j][i] = tmp*eps05[i]*eps05[j];
2408 }
2409 }
2410 }
2411
2412
ParamMix(double * Xin)2413 long int EOSPARAM::ParamMix( double *Xin )
2414 {
2415 long int j,i;
2416 double tmp,tmp1,tmp2;
2417
2418 for ( i=0; i<NComp; i++ )
2419 XX[i] = Xin[i];
2420
2421 emix = 0.;
2422 s3mix = 0.;
2423 for ( i=0; i<NComp-1; i++ )
2424 {
2425 for ( j=i+1; j<NComp; j++ )
2426 {
2427 tmp = XX[i]*XX[j];
2428 tmp2 = mixpar[j][i]; //eps
2429 tmp1 = mixpar[i][j]; //signa
2430 s3mix += tmp1*tmp;
2431 emix += tmp2*tmp;
2432 }
2433 }
2434 s3mix *= 2.;
2435 emix *= 2.;
2436 for ( i=0; i<NComp; i++ )
2437 {
2438 tmp = XX[i]*XX[i];
2439
2440 s3mix += sig3par[i]*tmp;
2441 emix += epspar[i]*tmp;
2442 }
2443 emix = emix/s3mix;
2444 return NComp;
2445 }
2446
2447
2448
2449 //=======================================================================================================
2450 // Soave-Redlich-Kwong (SRK) model for fluid mixtures
2451 // References: Soave (1972), Soave (1993)
2452 // (c) TW December 2008
2453 //=======================================================================================================
2454
2455
2456 // constructor
TSRKcalc(long int NCmp,double Pp,double Tkp)2457 TSRKcalc::TSRKcalc( long int NCmp, double Pp, double Tkp ):
2458 TSolMod( NCmp, 'E', Tkp, Pp )
2459 {
2460 // aGEX = 0;
2461 // aVol = 0;
2462 Pparc = 0;
2463 alloc_internal();
2464 }
2465
2466
2467
TSRKcalc(SolutionData * sd)2468 TSRKcalc::TSRKcalc( SolutionData *sd ):
2469 TSolMod( sd )
2470 {
2471 Pparc = aPparc;
2472 // aGEX = arGEX;
2473 // aVol = arVol;
2474 alloc_internal();
2475 }
2476
2477
2478
~TSRKcalc()2479 TSRKcalc::~TSRKcalc()
2480 {
2481 free_internal();
2482 }
2483
2484
2485
2486 /// allocate work arrays for pure fluid and fluid mixture properties
alloc_internal()2487 void TSRKcalc::alloc_internal()
2488 {
2489 Eosparm = new double [NComp][4];
2490 Pureparm = new double [NComp][4];
2491 Fugpure = new double [NComp][6];
2492 Fugci = new double [NComp][4];
2493 a = new double *[NComp];
2494 b = new double *[NComp];
2495 KK = new double *[NComp];
2496 dKK = new double *[NComp];
2497 d2KK = new double *[NComp];
2498 AA = new double *[NComp];
2499
2500 for (long int i=0; i<NComp; i++)
2501 {
2502 a[i] = new double[NComp];
2503 b[i] = new double[NComp];
2504 KK[i] = new double[NComp];
2505 dKK[i] = new double[NComp];
2506 d2KK[i] = new double[NComp];
2507 AA[i] = new double[NComp];
2508 }
2509 }
2510
2511
2512
free_internal()2513 void TSRKcalc::free_internal()
2514 {
2515 long int i;
2516
2517 for (i=0; i<NComp; i++)
2518 {
2519 delete[]a[i];
2520 delete[]b[i];
2521 delete[]KK[i];
2522 delete[]dKK[i];
2523 delete[]d2KK[i];
2524 delete[]AA[i];
2525 }
2526
2527 delete[]Eosparm;
2528 delete[]Pureparm;
2529 delete[]Fugpure;
2530 delete[]Fugci;
2531 delete[]a;
2532 delete[]b;
2533 delete[]KK;
2534 delete[]dKK;
2535 delete[]d2KK;
2536 delete[]AA;
2537
2538 }
2539
2540
2541
2542 /// high-level method to retrieve pure fluid fugacities
PureSpecies()2543 long int TSRKcalc::PureSpecies()
2544 {
2545 long int j, retCode = 0;
2546
2547 for( j=0; j<NComp; j++)
2548 {
2549 // Calling SRK EoS for pure fugacity
2550 retCode = FugacityPT( j, aDCc+j*NP_DC );
2551 aGEX[j] = log( Fugpure[j][0] );
2552 Pparc[j] = Fugpure[j][0]*Pbar; // fure fluid fugacity (required for performance)
2553 aVol[j] = Fugpure[j][4]*10.; // molar volume of pure fluid component, J/bar to cm3
2554 } // j
2555
2556 if ( retCode )
2557 {
2558 char buf[150];
2559 sprintf(buf, "SRK fluid: calculation of pure fugacity failed");
2560 Error( "E71IPM IPMgamma: ", buf );
2561 }
2562
2563 return 0;
2564 }
2565
2566
2567
2568 /// high-level method to calculate T,P corrected binary interaction parameters
PTparam()2569 long int TSRKcalc::PTparam()
2570 {
2571 long int j, i;
2572
2573 PureSpecies();
2574
2575 // set all interaction parameters zero
2576 for( j=0; j<NComp; j++ )
2577 {
2578 for( i=0; i<NComp; i++ )
2579 {
2580 KK[j][i] = 0.;
2581 dKK[j][i] = 0.;
2582 d2KK[j][i] = 0.;
2583 }
2584 }
2585
2586 switch ( MixCode )
2587 {
2588 case MR_UNDEF_:
2589 case MR_WAAL_:
2590 MixingWaals();
2591 break;
2592 case MR_CONST_:
2593 MixingConst();
2594 break;
2595 case MR_TEMP_:
2596 MixingTemp();
2597 break;
2598 default:
2599 break;
2600 }
2601
2602 return 0;
2603 }
2604
2605
2606
2607 /// high-level method to retrieve activity coefficients of the fluid mixture
MixMod()2608 long int TSRKcalc::MixMod()
2609 {
2610 long int j, iRet;
2611
2612 iRet = FugacitySpec( Pparc );
2613 phVOL[0] = PhVol * 10.;
2614
2615 for(j=0; j<NComp; j++)
2616 {
2617 if( Fugci[j][3] > 1e-23 )
2618 lnGamma[j] = log( Fugci[j][3] );
2619 else
2620 lnGamma[j] = 0;
2621 }
2622 if ( iRet )
2623 {
2624 char buf[150];
2625 sprintf(buf, "SRK fluid: calculation failed");
2626 Error( "E71IPM IPMgamma: ", buf );
2627 }
2628 return iRet;
2629 }
2630
2631
2632
2633 /// high-level method to retrieve residual functions of the fluid mixture
ExcessProp(double * Zex)2634 long int TSRKcalc::ExcessProp( double *Zex )
2635 {
2636 long int iRet;
2637
2638 iRet = ResidualFunct( Pparc );
2639
2640 if ( iRet )
2641 {
2642 char buf[150];
2643 sprintf(buf, "SRK fluid: calculation failed");
2644 Error( "E71IPM IPMgamma: ", buf );
2645 }
2646
2647 Ars = Grs - Vrs*Pbar;
2648 Urs = Hrs - Vrs*Pbar;
2649
2650 // assignments (residual functions)
2651 Zex[0] = Grs;
2652 Zex[1] = Hrs;
2653 Zex[2] = Srs;
2654 Zex[3] = CPrs;
2655 Zex[4] = Vrs;
2656 Zex[5] = Ars;
2657 Zex[6] = Urs;
2658
2659 return iRet;
2660
2661 }
2662
2663
2664
2665 /// basic van der waals mixing rule
MixingWaals()2666 long int TSRKcalc::MixingWaals()
2667 {
2668 // currently no calculations
2669
2670 return 0;
2671 }
2672
2673
2674
2675 /// constant one-term interaction parameter
MixingConst()2676 long int TSRKcalc::MixingConst()
2677 {
2678 long int ip, i1, i2;
2679 double k, dk, d2k;
2680
2681 if( NPcoef > 0 )
2682 {
2683 // transfer interaction parameters that have non-standard value
2684 for ( ip=0; ip<NPar; ip++ )
2685 {
2686 i1 = aIPx[MaxOrd*ip];
2687 i2 = aIPx[MaxOrd*ip+1];
2688 k = aIPc[NPcoef*ip];
2689 dk = 0.;
2690 d2k = 0.;
2691 KK[i1][i2] = k;
2692 dKK[i1][i2] = dk;
2693 d2KK[i1][i2] = d2k;
2694 KK[i2][i1] = k; // symmetric case
2695 dKK[i2][i1] = dk;
2696 d2KK[i2][i1] = d2k;
2697 }
2698 }
2699
2700 return 0;
2701 }
2702
2703
2704
2705 /// temperature dependent one-term interaction parameter
MixingTemp()2706 long int TSRKcalc::MixingTemp()
2707 {
2708 long int i, j, ip, i1, i2;
2709 double ai, aj, bi, bj, di, dj, dai, daj, d2ai, d2aj, ddi, ddj, d2di, d2dj,
2710 U, V, dU, dV, d2U, d2V, tmp, k, dk, d2k, C;
2711
2712 // set model specific interaction parameters zero
2713 for( j=0; j<NComp; j++ )
2714 {
2715 for( i=0; i<NComp; i++ )
2716 {
2717 a[j][i] = 0.;
2718 b[j][i] = 0.;
2719 }
2720 }
2721
2722 if( NPcoef > 0 )
2723 {
2724 // transfer parameters that have non-standard value
2725 for ( ip=0; ip<NPar; ip++ )
2726 {
2727 i1 = aIPx[MaxOrd*ip];
2728 i2 = aIPx[MaxOrd*ip+1];
2729 a[i1][i2] = aIPc[NPcoef*ip];
2730 b[i1][i2] = aIPc[NPcoef*ip+1];
2731 a[i2][i1] = aIPc[NPcoef*ip]; // symmetric case
2732 b[i2][i1] = aIPc[NPcoef*ip+1];
2733 }
2734 }
2735
2736 // calculate binary interaction parameters
2737 for (i=0; i<NComp; i++)
2738 {
2739 for (j=0; j<NComp; j++)
2740 {
2741 if ( a[i][j] == 0.0 )
2742 tmp = 1.0;
2743 else
2744 tmp = a[i][j];
2745
2746 // read a, b, da, d2a
2747 ai = Pureparm[i][0];
2748 aj = Pureparm[j][0];
2749 bi = Pureparm[i][1];
2750 bj = Pureparm[j][1];
2751 dai = Pureparm[i][2];
2752 daj = Pureparm[j][2];
2753 d2ai = Pureparm[i][3];
2754 d2aj = Pureparm[j][3];
2755
2756 // calculate k and derivatives
2757 di = sqrt(ai)/bi;
2758 dj = sqrt(aj)/bj;
2759 ddi = (0.5/bi) * pow(ai,-0.5) * dai;
2760 ddj = (0.5/bj) * pow(aj,-0.5) * daj;
2761 d2di = (0.5/bi) * ( (-0.5)*pow(ai,-1.5)*dai*dai + pow(ai,-0.5)*d2ai );
2762 d2dj = (0.5/bj) * ( (-0.5)*pow(aj,-1.5)*daj*daj + pow(aj,-0.5)*d2aj );
2763
2764 C = ( b[i][j]/tmp - 1. );
2765 U = a[i][j]*pow((298.15/Tk),C) - pow((di-dj),2.);
2766 V = (2.*di*dj);
2767 dU = - ( a[i][j]*C*pow((298.15/Tk),C) ) / Tk - 2.*(di-dj)*(ddi-ddj);
2768 dV = 2.*( ddi*dj + di*ddj );
2769 d2U = ( a[i][j]*pow(C,2.)*pow((298.15/Tk),C) ) / pow(Tk,2.)
2770 + ( a[i][j]*C*pow((298.15/Tk),C) ) / pow(Tk,2.)
2771 - 2.*( pow((ddi-ddj),2.) + (di-dj)*(d2di-d2dj) );
2772 d2V = 2.*( d2di*dj + 2.*ddi*ddj + di*d2dj );
2773 k = U/V;
2774 dk = (dU*V-U*dV)/pow(V,2.);
2775 d2k = (d2U*V+dU*dV)*pow(V,2.)/pow(V,4.) - (dU*V)*(2.*V*dV)/pow(V,4.)
2776 - (dU*dV+U*d2V)*pow(V,2.)/pow(V,4.) + (U*dV)*(2.*V*dV)/pow(V,4.);
2777
2778 // assignments
2779 KK[i][j] = k;
2780 dKK[i][j] = dk;
2781 d2KK[i][j] = d2k;
2782 }
2783 }
2784
2785 return 0;
2786 }
2787
2788
2789
2790 /// calculates ideal mixing properties
IdealProp(double * Zid)2791 long int TSRKcalc::IdealProp( double *Zid )
2792 {
2793 long int j;
2794 double s, sc, sp;
2795
2796 s = 0.0;
2797 for (j=0; j<NComp; j++)
2798 {
2799 if ( x[j] > 1.0e-32 )
2800 s += x[j]*log(x[j]);
2801 }
2802 sc = (-1.)*R_CONST*s;
2803 sp = (-1.)*R_CONST*log(Pbar);
2804 Hid = 0.0;
2805 CPid = 0.0;
2806 Vid = 0.0;
2807 Sid = sc + sp;
2808 Gid = Hid - Sid*Tk;
2809 Aid = Gid - Vid*Pbar;
2810 Uid = Hid - Vid*Pbar;
2811
2812 // assignments (ideal mixing properties)
2813 Zid[0] = Gid;
2814 Zid[1] = Hid;
2815 Zid[2] = Sid;
2816 Zid[3] = CPid;
2817 Zid[4] = Vid;
2818 Zid[5] = Aid;
2819 Zid[6] = Uid;
2820
2821 return 0;
2822 }
2823
2824
2825
2826 /// High-level method to retrieve pure fluid properties
FugacityPT(long int i,double * EoSparam)2827 long int TSRKcalc::FugacityPT( long int i, double *EoSparam )
2828 {
2829 long int iRet = 0;
2830 double Tcrit, Pcrit, omg, N, apure, bpure, da, d2a;
2831
2832 // reads EoS parameters from database into work array
2833 if( !EoSparam )
2834 return -1; // Memory alloc error
2835
2836 Eosparm[i][0] = EoSparam[0]; // critical temperature in K
2837 Eosparm[i][1] = EoSparam[1]; // critical pressure in bar
2838 Eosparm[i][2] = EoSparam[2]; // Pitzer acentric factor omega
2839 Eosparm[i][3] = EoSparam[3]; // empirical EoS parameter N
2840 Tcrit = Eosparm[i][0];
2841 Pcrit = Eosparm[i][1];
2842 omg = Eosparm[i][2];
2843 N = Eosparm[i][3];
2844
2845 AB( Tcrit, Pcrit, omg, N, apure, bpure, da, d2a );
2846
2847 Pureparm[i][0] = apure; // a parameter
2848 Pureparm[i][1] = bpure; // b parameter
2849 Pureparm[i][2] = da; // da/dT
2850 Pureparm[i][3] = d2a; // d2a/dT2
2851
2852 iRet = FugacityPure( i );
2853 if( iRet)
2854 return iRet;
2855
2856 return iRet;
2857 }
2858
2859
2860
2861 /// Calculates attractive (a) and repulsive (b) parameter of SRK equation of state
2862 /// and partial derivatives of alpha function
AB(double Tcrit,double Pcrit,double omg,double,double & apure,double & bpure,double & da,double & d2a)2863 long int TSRKcalc::AB( double Tcrit, double Pcrit, double omg, double /*N*/,
2864 double &apure, double &bpure, double &da, double &d2a )
2865 {
2866 double Tred, m, alph, ac, sqa, dsqa, d2sqa;
2867
2868 Tred = Tk/Tcrit;
2869 m = 0.48 + 1.574*omg - 0.176*pow(omg,2.);
2870 alph = pow(1. + m*(1-sqrt(Tred)), 2.);
2871 ac = (0.42747)*pow(R_CONST,2.)*pow(Tcrit,2.) / Pcrit;
2872 apure = alph*ac;
2873 bpure = (0.08664)*R_CONST*Tcrit/Pcrit;
2874 sqa = 1. + m*(1-sqrt(Tred));
2875 dsqa = -0.5*m/(sqrt(Tred)*Tcrit);
2876 da = 2.*ac*(sqa*dsqa);
2877 d2sqa = 0.25*m/(pow(Tred,1.5)*pow(Tcrit,2.));
2878 d2a = 2.*ac*(dsqa*dsqa + sqa*d2sqa);
2879
2880 return 0;
2881 }
2882
2883
2884
2885 /// Calculates fugacities and residual functions of pure fluid species
FugacityPure(long int i)2886 long int TSRKcalc::FugacityPure( long int i )
2887 {
2888 double /*Tcrit, Pcrit, Tred,*/ asrk, bsrk, da, d2a, A, B, a2, a1, a0,
2889 z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, z, vol, lnf;
2890 double /*gig, hig, sig, cpig,*/ fugpure, grs, hrs, srs, cprs,
2891 cv, dPdT, dPdV, dVdT;
2892
2893 // ideal gas changes from 1 bar to P at T of interest
2894 //hig = 0.;
2895 //sig = (-1.)*R_CONST*log(Pbar);
2896 //gig = hig - Tk*sig;
2897 //cpig = 0.;
2898
2899 // retrieve a and b terms of cubic EoS
2900 //Tcrit = Eosparm[i][0];
2901 //Pcrit = Eosparm[i][1];
2902 //Tred = Tk/Tcrit;
2903 asrk = Pureparm[i][0];
2904 bsrk = Pureparm[i][1];
2905 da = Pureparm[i][2];
2906 d2a = Pureparm[i][3];
2907
2908 // solve cubic equation
2909 A = asrk*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
2910 B = bsrk*Pbar/(R_CONST*Tk);
2911 a2 = (-1.);
2912 a1 = A-B-pow(B,2.);
2913 a0 = -A*B;
2914 Cardano(a2, a1, a0, z1, z2, z3);
2915
2916 // find stable roots
2917 vol1 = z1*R_CONST*Tk/Pbar;
2918 vol2 = z2*R_CONST*Tk/Pbar;
2919 vol3 = z3*R_CONST*Tk/Pbar;
2920
2921 if (z1 > B)
2922 lnf1 = z1 - 1 - log(z1-B) - A/B*log(1.+B/z1);
2923 else
2924 lnf1 = 1000.;
2925 if (z2 > B)
2926 lnf2 = z2 - 1 - log(z2-B) - A/B*log(1.+B/z2);
2927 else
2928 lnf2 = 1000.;
2929 if (z3 > B)
2930 lnf3 = z3 - 1 - log(z3-B) - A/B*log(1.+B/z3);
2931 else
2932 lnf3 = 1000.;
2933
2934 if (lnf2 < lnf1)
2935 {
2936 z = z2; vol = vol2; lnf = lnf2;
2937 }
2938 else
2939 {
2940 z = z1; vol = vol1; lnf = lnf1;
2941 }
2942 if (lnf3 < lnf)
2943 {
2944 z = z3; vol = vol3; lnf = lnf3;
2945 }
2946 else
2947 {
2948 z = z; vol = vol; lnf = lnf;
2949 }
2950
2951 // calculate thermodynamic properties
2952 hrs = - ( 1 - z + 1./(bsrk*R_CONST*Tk) * (asrk-Tk*da) * log(1.+ bsrk/vol) )*R_CONST*Tk;
2953 srs = ( log(z*(1.-bsrk/vol)) + 1./(bsrk*R_CONST) * da * log(1.+bsrk/vol) )*R_CONST;
2954 grs = hrs - Tk*srs;
2955
2956 // heat capacity part
2957 cv = Tk*d2a/bsrk * log(1.+B/z);
2958 dPdT = R_CONST/(vol-bsrk) - da/(vol*(vol+bsrk));
2959 dPdV = - R_CONST*Tk/pow((vol-bsrk),2.) + asrk*(2.*vol+bsrk)/pow((vol*(vol+bsrk)),2.);
2960 dVdT = (-1.)*(1./dPdV)*dPdT;
2961 cprs = cv + Tk*dPdT*dVdT - R_CONST;
2962
2963 // increment thermodynamic properties
2964 fugpure = exp(lnf);
2965 Fugpure[i][0] = fugpure;
2966 Fugpure[i][1] = grs;
2967 Fugpure[i][2] = hrs;
2968 Fugpure[i][3] = srs;
2969 Fugpure[i][4] = vol;
2970 Fugpure[i][5] = cprs;
2971
2972 return 0;
2973 }
2974
2975
2976
2977 /// Cubic equation root solver based on Cardanos method
Cardano(double a2,double a1,double a0,double & z1,double & z2,double & z3)2978 long int TSRKcalc::Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 )
2979 {
2980 double q, rc, q3, rc2, theta, ac, bc;
2981
2982 q = (pow(a2,2.) - 3.*a1)/9.;
2983 rc = (2.*pow(a2,3.) - 9.*a2*a1 + 27.*a0)/54.;
2984 q3 = pow(q,3.);
2985 rc2 = pow(rc,2.);
2986
2987 if (rc2 < q3) // three real roots
2988 {
2989 theta = acos(rc/sqrt(q3));
2990 z1 = (-2.)*sqrt(q)*cos(theta/3.)-a2/3.;
2991 z2 = (-2.)*sqrt(q)*cos(theta/3.+2./3.*3.1415927)-a2/3.;
2992 z3 = (-2.)*sqrt(q)*cos(theta/3.-2./3.*3.1415927)-a2/3.;
2993 }
2994
2995 else // one real root
2996 {
2997 ac = (-1.)*rc/fabs(rc)*pow(fabs(rc)+sqrt(rc2-q3), 1./3.);
2998 if (ac != 0.)
2999 bc = q/ac;
3000 else
3001 bc = 0.;
3002 z1 = ac+bc-a2/3.;
3003 z2 = ac+bc-a2/3.;
3004 z3 = ac+bc-a2/3.;
3005 }
3006 return 0;
3007 }
3008
3009
3010
3011 /// Calculates mixing properties of the fluid mixture
MixParam(double & amix,double & bmix)3012 long int TSRKcalc::MixParam( double &amix, double &bmix )
3013 {
3014 long int i, j;
3015 double K;
3016 amix = 0.;
3017 bmix = 0.;
3018
3019 // calculate binary aij parameters
3020 for (i=0; i<NComp; i++)
3021 {
3022 for (j=0; j<NComp; j++)
3023 {
3024 K = KK[i][j];
3025 AA[i][j] = sqrt(Pureparm[i][0]*Pureparm[j][0])*(1.-K);
3026 }
3027 }
3028
3029 // find a and b of the mixture
3030 for (i=0; i<NComp; i++)
3031 {
3032 for (j=0; j<NComp; j++)
3033 {
3034 amix = amix + x[i]*x[j]*AA[i][j];
3035 }
3036 }
3037 for (i=0; i<NComp; i++)
3038 {
3039 bmix = bmix + x[i]*Pureparm[i][1];
3040 }
3041 return 0;
3042 }
3043
3044
3045
3046 /// Calculates fugacity of the bulk fluid mixture
FugacityMix(double amix,double bmix,double & fugmix,double & zmix,double & vmix)3047 long int TSRKcalc::FugacityMix( double amix, double bmix,
3048 double &fugmix, double &zmix, double &vmix )
3049 {
3050 double A, B, a2, a1, a0, z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, lnf;
3051
3052 // solve cubic equation
3053 A = amix*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
3054 B = bmix*Pbar/(R_CONST*Tk);
3055 a2 = (-1.);
3056 a1 = A-B-pow(B,2.);
3057 a0 = -A*B;
3058 Cardano(a2, a1, a0, z1, z2, z3);
3059
3060 // find stable roots
3061 vol1 = z1*R_CONST*Tk/Pbar;
3062 vol2 = z2*R_CONST*Tk/Pbar;
3063 vol3 = z3*R_CONST*Tk/Pbar;
3064
3065 if (z1 > B)
3066 lnf1 = z1 - 1 - log(z1-B) - A/B*log(1.+B/z1);
3067 else
3068 lnf1 = 1000.;
3069 if (z2 > B)
3070 lnf2 = z2 - 1 - log(z2-B) - A/B*log(1.+B/z2);
3071 else
3072 lnf2 = 1000.;
3073 if (z3 > B)
3074 lnf3 = z3 - 1 - log(z3-B) - A/B*log(1.+B/z3);
3075 else
3076 lnf3 = 1000.;
3077
3078 if (lnf2 < lnf1)
3079 {
3080 zmix = z2; vmix = vol2; lnf = lnf2;
3081 }
3082 else
3083 {
3084 zmix = z1; vmix = vol1; lnf = lnf1;
3085 }
3086 if (lnf3 < lnf)
3087 {
3088 zmix = z3; vmix = vol3; lnf = lnf3;
3089 }
3090 else
3091 {
3092 zmix = zmix; vmix = vmix; lnf = lnf;
3093 }
3094
3095 fugmix = exp(lnf);
3096 PhVol = vmix;
3097 return 0;
3098 }
3099
3100
3101
3102 /// Calculates fugacities and activities of fluid species in the mixture,
FugacitySpec(double * fugpure)3103 long int TSRKcalc::FugacitySpec( double *fugpure )
3104 {
3105 long int i, j, iRet=0;
3106 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0., sum=0.;
3107 double A, B, bi, Bi, lnfci, fci;
3108
3109 // Reload params to Pureparm (possibly not required any more)
3110 for( j=0; j<NComp; j++ )
3111 {
3112 Fugpure[j][0] = fugpure[j]/Pbar;
3113 }
3114
3115 // calculate properties of the mixture
3116 iRet = MixParam( amix, bmix);
3117 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix);
3118 A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
3119 B = bmix*Pbar/(R_CONST*Tk);
3120
3121 // calculate fugacity coefficient, fugacity and activity of species i
3122 for (i=0; i<NComp; i++)
3123 {
3124 bi = Pureparm[i][1];
3125 Bi = bi*Pbar/(R_CONST*Tk);
3126
3127 sum = 0.;
3128 for (j=0; j<NComp; j++)
3129 {
3130 sum = sum + x[j]*AA[i][j];
3131 }
3132
3133 lnfci = Bi/B*(zmix-1.) - log(zmix-B)
3134 + A/B * ( Bi/B - 2./amix*sum ) * log(1.+B/zmix);
3135 fci = exp(lnfci);
3136 Fugci[i][0] = fci; // fugacity coefficient using engineering convention
3137 Fugci[i][1] = x[i]*fci; // fugacity coefficient using geology convention
3138 Fugci[i][2] = Fugci[i][1]/Fugpure[i][0]; // activity of species
3139 if (x[i]>1.0e-20)
3140 Fugci[i][3] = Fugci[i][2]/x[i]; // activity coefficient of species
3141 else
3142 Fugci[i][3] = 1.0;
3143 }
3144
3145 return iRet;
3146 }
3147
3148
3149
3150 /// calculates residual functions in the mixture
ResidualFunct(double * fugpure)3151 long int TSRKcalc::ResidualFunct( double *fugpure )
3152 {
3153 long int i, j, iRet=0;
3154 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0.;
3155 double /*A,*/ B, K, dK, d2K, Q, dQ, d2Q, damix, d2amix, ai, aj, dai, daj, d2ai, d2aj,
3156 cv, dPdT, dPdV, dVdT;
3157
3158 // Reload params to Pureparm (possibly not required any more)
3159 for( j=0; j<NComp; j++ )
3160 {
3161 Fugpure[j][0] = fugpure[j]/Pbar;
3162 }
3163
3164 // calculate properties of the mixture
3165 iRet = MixParam( amix, bmix);
3166 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix);
3167 //A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
3168 B = bmix*Pbar/(R_CONST*Tk);
3169
3170 // calculate total state functions of the mixture
3171 damix = 0.;
3172 d2amix = 0.;
3173 for (i=0; i<NComp; i++)
3174 {
3175 for (j=0; j<NComp; j++)
3176 {
3177 // pull parameters
3178 ai = Pureparm[i][0];
3179 aj = Pureparm[j][0];
3180 dai = Pureparm[i][2];
3181 daj = Pureparm[j][2];
3182 d2ai = Pureparm[i][3];
3183 d2aj = Pureparm[j][3];
3184 K = KK[i][j];
3185 dK = dKK[i][j];
3186 d2K = d2KK[i][j];
3187
3188 // increments to derivatives
3189 Q = sqrt(ai*aj);
3190 dQ = 0.5*( sqrt(aj/ai)*dai + sqrt(ai/aj)*daj );
3191 d2Q = 0.5*( dai*daj/sqrt(ai*aj) + d2ai*sqrt(aj)/sqrt(ai) + d2aj*sqrt(ai)/sqrt(aj)
3192 - 0.5*( pow(dai,2.)*sqrt(aj)/sqrt(pow(ai,3.))
3193 + pow(daj,2.)*sqrt(ai)/sqrt(pow(aj,3.)) ) );
3194 damix = damix + x[i]*x[j] * ( dQ*(1.-K) - Q*dK );
3195 d2amix = d2amix + x[i]*x[j] * ( d2Q*(1.-K) - 2.*dQ*dK - Q*d2K );
3196 }
3197 }
3198
3199 // calculate thermodynamic properties
3200 Hrs = - ( 1. - zmix + 1./(bmix*R_CONST*Tk) * (amix-Tk*damix )
3201 * log(1.+bmix/vmix) )*R_CONST*Tk;
3202 Srs = ( log(zmix*(1.-bmix/vmix)) + 1./(bmix*R_CONST)*damix
3203 * log(1.+bmix/vmix) )*R_CONST;
3204 Grs = Hrs - Tk*Srs;
3205
3206 // heat capacity part
3207 cv = Tk*d2amix/bmix * log(1.+B/zmix);
3208 dPdT = R_CONST/(vmix-bmix) - damix/(vmix*(vmix+bmix));
3209 dPdV = - R_CONST*Tk/pow((vmix-bmix),2.) + amix*(2.*vmix+bmix)/pow((vmix*(vmix+bmix)),2.);
3210 dVdT = (-1.)*(1./dPdV)*dPdT;
3211 CPrs = cv + Tk*dPdT*dVdT - R_CONST;
3212 Vrs = vmix;
3213
3214 return iRet;
3215 }
3216
3217
3218
3219 #ifndef IPMGEMPLUGIN_
3220
3221 /// Calculates properties of pure fluids when called from DCthermo
SRKCalcFugPure(double Tmin,double * Cpg,double * FugProps)3222 long int TSRKcalc::SRKCalcFugPure(double Tmin, double *Cpg, double *FugProps )
3223 {
3224 long int retCode = 0;
3225 double Coeff[7];
3226
3227 for( int ii=0; ii<7; ii++ )
3228 Coeff[ii] = (double)Cpg[ii];
3229
3230 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
3231 {
3232 retCode = FugacityPT( 0, Coeff );
3233 for( int i=0; i<6; i++ )
3234 FugProps[i] = Fugpure[0][i];
3235 return retCode;
3236 }
3237
3238 else
3239 {
3240 for( int i=1; i<6; i++ )
3241 FugProps[i] = 0.;
3242 FugProps[0] = 1.;
3243 FugProps[4] = 8.31451*Tk/Pbar;
3244 return -1;
3245 }
3246 }
3247
3248 #endif
3249
3250
3251
3252
3253
3254 //=======================================================================================================
3255 // Peng-Robinson (PR78) model for fluid mixtures
3256 // References: Peng and Robinson (1976), Peng and Robinson (1978)
3257 // (c) TW July 2009
3258 //=======================================================================================================
3259
3260
3261 // Constructor
TPR78calc(long int NCmp,double Pp,double Tkp)3262 TPR78calc::TPR78calc( long int NCmp, double Pp, double Tkp ):
3263 TSolMod( NCmp, '7', Tkp, Pp )
3264 {
3265 // aGEX = 0;
3266 // aVol = 0;
3267 Pparc = 0;
3268 alloc_internal();
3269 }
3270
3271
3272
TPR78calc(SolutionData * sd)3273 TPR78calc::TPR78calc( SolutionData *sd ):
3274 TSolMod( sd )
3275 {
3276 Pparc = aPparc;
3277 // aGEX = arGEX;
3278 // aVol = arVol;
3279 alloc_internal();
3280 }
3281
3282
3283
~TPR78calc()3284 TPR78calc::~TPR78calc()
3285 {
3286 free_internal();
3287 }
3288
3289
3290
3291 /// allocate work arrays for pure fluid and fluid mixture properties
alloc_internal()3292 void TPR78calc::alloc_internal()
3293 {
3294 Eosparm = new double [NComp][4];
3295 Pureparm = new double [NComp][4];
3296 Fugpure = new double [NComp][6];
3297 Fugci = new double [NComp][4];
3298 a = new double *[NComp];
3299 b = new double *[NComp];
3300 KK = new double *[NComp];
3301 dKK = new double *[NComp];
3302 d2KK = new double *[NComp];
3303 AA = new double *[NComp];
3304
3305 for (long int i=0; i<NComp; i++)
3306 {
3307 a[i] = new double[NComp];
3308 b[i] = new double[NComp];
3309 KK[i] = new double[NComp];
3310 dKK[i] = new double[NComp];
3311 d2KK[i] = new double[NComp];
3312 AA[i] = new double[NComp];
3313 }
3314 }
3315
3316
3317
free_internal()3318 void TPR78calc::free_internal()
3319 {
3320 long int i;
3321
3322 for (i=0; i<NComp; i++)
3323 {
3324 delete[]a[i];
3325 delete[]b[i];
3326 delete[]KK[i];
3327 delete[]dKK[i];
3328 delete[]d2KK[i];
3329 delete[]AA[i];
3330 }
3331
3332 delete[]Eosparm;
3333 delete[]Pureparm;
3334 delete[]Fugpure;
3335 delete[]Fugci;
3336 delete[]a;
3337 delete[]b;
3338 delete[]KK;
3339 delete[]dKK;
3340 delete[]d2KK;
3341 delete[]AA;
3342
3343 }
3344
3345
3346
3347 /// High-level method to retrieve pure fluid fugacities
PureSpecies()3348 long int TPR78calc::PureSpecies()
3349 {
3350 long int j, retCode = 0;
3351
3352 for( j=0; j<NComp; j++)
3353 {
3354 // Calling SRK EoS for pure fugacity
3355 retCode = FugacityPT( j, aDCc+j*NP_DC );
3356 aGEX[j] = log( Fugpure[j][0] );
3357 Pparc[j] = Fugpure[j][0]*Pbar; // fure fluid fugacity (required for performance)
3358 aVol[j] = Fugpure[j][4]*10.; // molar volume of pure fluid component, J/bar to cm3
3359 } // j
3360
3361 if ( retCode )
3362 {
3363 char buf[150];
3364 sprintf(buf, "PR78 fluid: calculation of pure fugacity failed");
3365 Error( "E71IPM IPMgamma: ", buf );
3366 }
3367
3368 return 0;
3369 }
3370
3371
3372 /// High-level method to calculate T,P corrected binary interaction parameters
PTparam()3373 long int TPR78calc::PTparam()
3374 {
3375 long int j, i;
3376
3377 PureSpecies();
3378
3379 // set all interaction parameters zero
3380 for( j=0; j<NComp; j++ )
3381 {
3382 for( i=0; i<NComp; i++ )
3383 {
3384 KK[j][i] = 0.;
3385 dKK[j][i] = 0.;
3386 d2KK[j][i] = 0.;
3387 }
3388 }
3389
3390 switch ( MixCode )
3391 {
3392 case MR_UNDEF_:
3393 case MR_WAAL_:
3394 MixingWaals();
3395 break;
3396 case MR_CONST_:
3397 MixingConst();
3398 break;
3399 case MR_TEMP_:
3400 MixingTemp();
3401 break;
3402 default:
3403 break;
3404 }
3405
3406 return 0;
3407 }
3408
3409
3410
3411 /// High-level method to retrieve activity coefficients of the fluid mixture
MixMod()3412 long int TPR78calc::MixMod()
3413 {
3414 long int j, iRet;
3415
3416 iRet = FugacitySpec( Pparc );
3417 phVOL[0] = PhVol * 10.;
3418
3419 for(j=0; j<NComp; j++)
3420 {
3421 if( Fugci[j][3] > 1e-23 )
3422 lnGamma[j] = log( Fugci[j][3] );
3423 else
3424 lnGamma[j] = 0;
3425 }
3426 if ( iRet )
3427 {
3428 char buf[150];
3429 sprintf(buf, "PR78 fluid: calculation failed");
3430 Error( "E71IPM IPMgamma: ", buf );
3431 }
3432 return iRet;
3433 }
3434
3435
3436
3437 /// High-level method to retrieve residual functions of the fluid mixture
ExcessProp(double * Zex)3438 long int TPR78calc::ExcessProp( double *Zex )
3439 {
3440 long int iRet;
3441
3442 iRet = ResidualFunct( Pparc );
3443
3444 if ( iRet )
3445 {
3446 char buf[150];
3447 sprintf(buf, "PR78 fluid: calculation failed");
3448 Error( "E71IPM IPMgamma: ", buf );
3449 }
3450
3451 Ars = Grs - Vrs*Pbar;
3452 Urs = Hrs - Vrs*Pbar;
3453
3454 // assignments (residual functions)
3455 Zex[0] = Grs;
3456 Zex[1] = Hrs;
3457 Zex[2] = Srs;
3458 Zex[3] = CPrs;
3459 Zex[4] = Vrs;
3460 Zex[5] = Ars;
3461 Zex[6] = Urs;
3462
3463 return iRet;
3464
3465 }
3466
3467
3468
3469 /// basic van der waals mixing rule
MixingWaals()3470 long int TPR78calc::MixingWaals()
3471 {
3472 // currently no calculations
3473
3474 return 0;
3475 }
3476
3477
3478
3479 /// constant one-term interaction parameter
MixingConst()3480 long int TPR78calc::MixingConst()
3481 {
3482 long int ip, i1, i2;
3483 double k, dk, d2k;
3484
3485 if( NPcoef > 0 )
3486 {
3487 // transfer interaction parameters that have non-standard value
3488 for ( ip=0; ip<NPar; ip++ )
3489 {
3490 i1 = aIPx[MaxOrd*ip];
3491 i2 = aIPx[MaxOrd*ip+1];
3492 k = aIPc[NPcoef*ip];
3493 dk = 0.;
3494 d2k = 0.;
3495 KK[i1][i2] = k;
3496 dKK[i1][i2] = dk;
3497 d2KK[i1][i2] = d2k;
3498 KK[i2][i1] = k; // symmetric case
3499 dKK[i2][i1] = dk;
3500 d2KK[i2][i1] = d2k;
3501 }
3502 }
3503
3504 return 0;
3505 }
3506
3507
3508
3509 /// temperature dependent one-term interaction parameter
MixingTemp()3510 long int TPR78calc::MixingTemp()
3511 {
3512 long int i, j, ip, i1, i2;
3513 double ai, aj, bi, bj, di, dj, dai, daj, d2ai, d2aj, ddi, ddj, d2di, d2dj,
3514 U, V, dU, dV, d2U, d2V, tmp, k, dk, d2k, C;
3515
3516 // set model specific interaction parameters zero
3517 for( j=0; j<NComp; j++ )
3518 {
3519 for( i=0; i<NComp; i++ )
3520 {
3521 a[j][i] = 0.;
3522 b[j][i] = 0.;
3523 }
3524 }
3525
3526 if( NPcoef > 0 )
3527 {
3528 // transfer parameters that have non-standard value
3529 for ( ip=0; ip<NPar; ip++ )
3530 {
3531 i1 = aIPx[MaxOrd*ip];
3532 i2 = aIPx[MaxOrd*ip+1];
3533 a[i1][i2] = aIPc[NPcoef*ip];
3534 b[i1][i2] = aIPc[NPcoef*ip+1];
3535 a[i2][i1] = aIPc[NPcoef*ip]; // symmetric case
3536 b[i2][i1] = aIPc[NPcoef*ip+1];
3537 }
3538 }
3539
3540 // calculate binary interaction parameters
3541 for (i=0; i<NComp; i++)
3542 {
3543 for (j=0; j<NComp; j++)
3544 {
3545 if ( a[i][j] == 0.0 )
3546 tmp = 1.0;
3547 else
3548 tmp = a[i][j];
3549
3550 // read a, b, da, d2a
3551 ai = Pureparm[i][0];
3552 aj = Pureparm[j][0];
3553 bi = Pureparm[i][1];
3554 bj = Pureparm[j][1];
3555 dai = Pureparm[i][2];
3556 daj = Pureparm[j][2];
3557 d2ai = Pureparm[i][3];
3558 d2aj = Pureparm[j][3];
3559
3560 // calculate k and derivatives
3561 di = sqrt(ai)/bi;
3562 dj = sqrt(aj)/bj;
3563 ddi = (0.5/bi) * pow(ai,-0.5) * dai;
3564 ddj = (0.5/bj) * pow(aj,-0.5) * daj;
3565 d2di = (0.5/bi) * ( (-0.5)*pow(ai,-1.5)*dai*dai + pow(ai,-0.5)*d2ai );
3566 d2dj = (0.5/bj) * ( (-0.5)*pow(aj,-1.5)*daj*daj + pow(aj,-0.5)*d2aj );
3567
3568 C = ( b[i][j]/tmp - 1. );
3569 U = a[i][j]*pow((298.15/Tk),C) - pow((di-dj),2.);
3570 V = (2.*di*dj);
3571 dU = - ( a[i][j]*C*pow((298.15/Tk),C) ) / Tk - 2.*(di-dj)*(ddi-ddj);
3572 dV = 2.*( ddi*dj + di*ddj );
3573 d2U = ( a[i][j]*pow(C,2.)*pow((298.15/Tk),C) ) / pow(Tk,2.)
3574 + ( a[i][j]*C*pow((298.15/Tk),C) ) / pow(Tk,2.)
3575 - 2.*( pow((ddi-ddj),2.) + (di-dj)*(d2di-d2dj) );
3576 d2V = 2.*( d2di*dj + 2.*ddi*ddj + di*d2dj );
3577 k = U/V;
3578 dk = (dU*V-U*dV)/pow(V,2.);
3579 d2k = (d2U*V+dU*dV)*pow(V,2.)/pow(V,4.) - (dU*V)*(2.*V*dV)/pow(V,4.)
3580 - (dU*dV+U*d2V)*pow(V,2.)/pow(V,4.) + (U*dV)*(2.*V*dV)/pow(V,4.);
3581
3582 // assignments
3583 KK[i][j] = k;
3584 dKK[i][j] = dk;
3585 d2KK[i][j] = d2k;
3586 }
3587 }
3588
3589 return 0;
3590 }
3591
3592
3593
3594 /// calculates ideal mixing properties
IdealProp(double * Zid)3595 long int TPR78calc::IdealProp( double *Zid )
3596 {
3597 long int j;
3598 double s, sc, sp;
3599
3600 s = 0.0;
3601 for (j=0; j<NComp; j++)
3602 {
3603 if ( x[j] > 1.0e-32 )
3604 s += x[j]*log(x[j]);
3605 }
3606 sc = (-1.)*R_CONST*s;
3607 sp = (-1.)*R_CONST*log(Pbar);
3608 Hid = 0.0;
3609 CPid = 0.0;
3610 Vid = 0.0;
3611 Sid = sc + sp;
3612 Gid = Hid - Sid*Tk;
3613 Aid = Gid - Vid*Pbar;
3614 Uid = Hid - Vid*Pbar;
3615
3616 // assignments (ideal mixing properties)
3617 Zid[0] = Gid;
3618 Zid[1] = Hid;
3619 Zid[2] = Sid;
3620 Zid[3] = CPid;
3621 Zid[4] = Vid;
3622 Zid[5] = Aid;
3623 Zid[6] = Uid;
3624
3625 return 0;
3626 }
3627
3628
3629
3630 /// High-level method to retrieve pure fluid properties
FugacityPT(long int i,double * EoSparam)3631 long int TPR78calc::FugacityPT( long int i, double *EoSparam )
3632 {
3633 long int iRet = 0;
3634 double Tcrit, Pcrit, omg, N, apure, bpure, da, d2a;
3635
3636 // reads EoS parameters from database into work array
3637 if( !EoSparam )
3638 return -1; // Memory alloc error
3639
3640 Eosparm[i][0] = EoSparam[0]; // critical temperature in K
3641 Eosparm[i][1] = EoSparam[1]; // critical pressure in bar
3642 Eosparm[i][2] = EoSparam[2]; // Pitzer acentric factor omega
3643 Eosparm[i][3] = EoSparam[3]; // empirical EoS parameter N
3644 Tcrit = Eosparm[i][0];
3645 Pcrit = Eosparm[i][1];
3646 omg = Eosparm[i][2];
3647 N = Eosparm[i][3];
3648
3649 AB(Tcrit, Pcrit, omg, N, apure, bpure, da, d2a);
3650
3651 Pureparm[i][0] = apure; // a parameter
3652 Pureparm[i][1] = bpure; // b parameter
3653 Pureparm[i][2] = da; // da/dT
3654 Pureparm[i][3] = d2a; // d2a/dT2
3655
3656 iRet = FugacityPure( i );
3657 if( iRet)
3658 return iRet;
3659
3660 return iRet;
3661 }
3662
3663
3664
3665 /// Calculates attractive (a) and repulsive (b) parameter of SRK equation of state
3666 /// and partial derivatives of alpha function
AB(double Tcrit,double Pcrit,double omg,double,double & apure,double & bpure,double & da,double & d2a)3667 long int TPR78calc::AB( double Tcrit, double Pcrit, double omg, double /*N*/,
3668 double &apure, double &bpure, double &da, double &d2a )
3669 {
3670 double Tred, k, alph, ac, sqa, dsqa, d2sqa;
3671
3672 Tred = Tk/Tcrit;
3673 if (omg <= 0.491)
3674 k = 0.37464 + 1.54226*omg - 0.26992*pow(omg,2.);
3675 else
3676 k = 0.379642 + 1.48503*omg - 0.164423*pow(omg,2.) + 0.0166666*pow(omg,3.);
3677
3678 alph = pow(1.+k*(1.-sqrt(Tred)),2.);
3679 ac = (0.457235529)*pow(R_CONST,2.)*pow(Tcrit,2.) / Pcrit;
3680 apure = alph*ac;
3681 bpure = (0.0777960739)*R_CONST*Tcrit/Pcrit;
3682 sqa = 1.+k*(1.-sqrt(Tred));
3683 dsqa = -0.5*k/(sqrt(Tred)*Tcrit);
3684 da = 2.*ac*(sqa*dsqa);
3685 d2sqa = 0.25*k/(pow(Tred,1.5)*pow(Tcrit,2.));
3686 d2a = 2.*ac*(dsqa*dsqa + sqa*d2sqa);
3687
3688 return 0;
3689 }
3690
3691
3692
3693 /// Calculates fugacities and residual functions of pure fluid species
FugacityPure(long int i)3694 long int TPR78calc::FugacityPure( long int i )
3695 {
3696 double Tcrit, Pcrit, Tred, apr, bpr, alph, da, d2a, k, A, B, a2, a1, a0,
3697 z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, z, vol, lnf;
3698 double /*gig, hig, sig, cpig,*/ fugpure, grs, hrs, srs, cprs,
3699 cv, dPdT, dPdV, dVdT;
3700
3701 // ideal gas changes from 1 bar to P at T of interest
3702 //hig = 0.;
3703 //sig = (-1.)*R_CONST*log(Pbar);
3704 //gig = hig - Tk*sig;
3705 //cpig = 0.;
3706
3707 // retrieve a and b terms of cubic EoS
3708 Tcrit = Eosparm[i][0];
3709 Pcrit = Eosparm[i][1];
3710 Tred = Tk/Tcrit;
3711 apr = Pureparm[i][0];
3712 bpr = Pureparm[i][1];
3713 da = Pureparm[i][2];
3714 d2a = Pureparm[i][3];
3715
3716 // solve cubic equation
3717 A = apr*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
3718 B = bpr*Pbar/(R_CONST*Tk);
3719 a2 = B - 1.;
3720 a1 = A - 3.*pow(B,2.) - 2.*B;
3721 a0 = pow(B,3.) + pow(B,2.) - A*B;
3722 Cardano(a2, a1, a0, z1, z2, z3);
3723
3724 // find stable roots
3725 vol1 = z1*R_CONST*Tk/Pbar;
3726 vol2 = z2*R_CONST*Tk/Pbar;
3727 vol3 = z3*R_CONST*Tk/Pbar;
3728 if (z1 > B)
3729 lnf1 = (-1.)*log(z1-B)
3730 - A/(B*sqrt(8.))*log((z1+(1.+sqrt(2.))*B)/(z1+(1.-sqrt(2.))*B))+z1-1.;
3731 else
3732 lnf1 = 1000.;
3733 if (z2 > B)
3734 lnf2 = (-1.)*log(z2-B)
3735 - A/(B*sqrt(8.))*log((z2+(1.+sqrt(2.))*B)/(z2+(1.-sqrt(2.))*B))+z2-1.;
3736 else
3737 lnf2 = 1000.;
3738 if (z3 > B)
3739 lnf3 = (-1.)*log(z3-B)
3740 - A/(B*sqrt(8.))*log((z3+(1.+sqrt(2.))*B)/(z3+(1.-sqrt(2.))*B))+z3-1.;
3741 else
3742 lnf3 = 1000.;
3743
3744 if (lnf2 < lnf1)
3745 {
3746 z = z2; vol = vol2; lnf = lnf2;
3747 }
3748 else
3749 {
3750 z = z1; vol = vol1; lnf = lnf1;
3751 }
3752 if (lnf3 < lnf)
3753 {
3754 z = z3; vol = vol3; lnf = lnf3;
3755 }
3756 else
3757 {
3758 z = z; vol = vol; lnf = lnf;
3759 }
3760
3761 // calculate thermodynamic properties
3762 alph = apr/((0.457235529)*pow(R_CONST,2.)*pow(Tcrit,2.) / Pcrit);
3763 k = (sqrt(alph)-1.)/(1.-sqrt(Tred));
3764 grs = R_CONST*Tk*(z-1.-log(z-B)-A/(B*sqrt(8.))
3765 *log((z+(1+sqrt(2.))*B)/(z+(1-sqrt(2.))*B)));
3766 hrs = R_CONST*Tk*(z-1.-log((z+(1+sqrt(2.))*B)/(z+(1-sqrt(2.))*B))
3767 *A/(B*sqrt(8.))*(1+k*sqrt(Tred)/sqrt(alph)));
3768 srs = (hrs-grs)/Tk;
3769
3770 // heat capacity part
3771 cv = Tk*d2a/(bpr*sqrt(8.))
3772 * log( (z+B*(1.+sqrt(2.)))/(z+B*(1.-sqrt(2.))) );
3773 dPdT = R_CONST/(vol-bpr) - da/( vol*(vol+bpr) + bpr*(vol-bpr) );
3774 dPdV = - R_CONST*Tk/pow((vol-bpr),2.) + 2.*apr*(vol+bpr)/pow((vol*(vol+bpr)+bpr*(vol-bpr)),2.);
3775 dVdT = (-1.)*(1./dPdV)*dPdT;
3776 cprs = cv + Tk*dPdT*dVdT - R_CONST;
3777
3778 // increment thermodynamic properties
3779 fugpure = exp(lnf);
3780 Fugpure[i][0] = fugpure;
3781 Fugpure[i][1] = grs;
3782 Fugpure[i][2] = hrs;
3783 Fugpure[i][3] = srs;
3784 Fugpure[i][4] = vol;
3785 Fugpure[i][5] = cprs;
3786
3787 return 0;
3788 }
3789
3790
3791
3792 /// Cubic equation root solver based on Cardanos method
Cardano(double a2,double a1,double a0,double & z1,double & z2,double & z3)3793 long int TPR78calc::Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 )
3794 {
3795 double q, rc, q3, rc2, theta, ac, bc;
3796
3797 q = (pow(a2,2.) - 3.*a1)/9.;
3798 rc = (2.*pow(a2,3.) - 9.*a2*a1 + 27.*a0)/54.;
3799 q3 = pow(q,3.);
3800 rc2 = pow(rc,2.);
3801 if (rc2 < q3) // three real roots
3802 {
3803 theta = acos(rc/sqrt(q3));
3804 z1 = (-2.)*sqrt(q)*cos(theta/3.)-a2/3.;
3805 z2 = (-2.)*sqrt(q)*cos(theta/3.+2./3.*3.1415927)-a2/3.;
3806 z3 = (-2.)*sqrt(q)*cos(theta/3.-2./3.*3.1415927)-a2/3.;
3807 }
3808 else // one real root
3809 {
3810 ac = (-1.)*rc/fabs(rc)*pow(fabs(rc)+sqrt(rc2-q3), 1./3.);
3811 if (ac != 0.)
3812 bc = q/ac;
3813 else
3814 bc = 0.;
3815 z1 = ac+bc-a2/3.;
3816 z2 = ac+bc-a2/3.;
3817 z3 = ac+bc-a2/3.;
3818 }
3819 return 0;
3820 }
3821
3822
3823
3824 /// Calculates mixing properties of the fluid mixture
MixParam(double & amix,double & bmix)3825 long int TPR78calc::MixParam( double &amix, double &bmix )
3826 {
3827 long int i, j;
3828 double K;
3829 amix = 0.;
3830 bmix = 0.;
3831
3832 // calculate binary aij parameters
3833 for (i=0; i<NComp; i++)
3834 {
3835 for (j=0; j<NComp; j++)
3836 {
3837 K = KK[i][j];
3838 AA[i][j] = sqrt(Pureparm[i][0]*Pureparm[j][0])*(1.-K);
3839 }
3840 }
3841
3842 // find a and b of the mixture
3843 for (i=0; i<NComp; i++)
3844 {
3845 for (j=0; j<NComp; j++)
3846 {
3847 amix = amix + x[i]*x[j]*AA[i][j];
3848 }
3849 }
3850 for (i=0; i<NComp; i++)
3851 {
3852 bmix = bmix + x[i]*Pureparm[i][1];
3853 }
3854
3855 return 0;
3856 }
3857
3858
3859
3860 /// Calculates fugacity of the bulk fluid mixture
FugacityMix(double amix,double bmix,double & fugmix,double & zmix,double & vmix)3861 long int TPR78calc::FugacityMix( double amix, double bmix,
3862 double &fugmix, double &zmix, double &vmix )
3863 {
3864 double A, B, a2, a1, a0, z1, z2, z3, vol1, vol2, vol3, lnf1, lnf2, lnf3, lnf;
3865
3866 // solve cubic equation
3867 A = amix*Pbar/(pow(R_CONST,2.)*pow(Tk,2.));
3868 B = bmix*Pbar/(R_CONST*Tk);
3869 a2 = B - 1.;
3870 a1 = A - 3.*pow(B,2.) - 2.*B;
3871 a0 = pow(B,3.) + pow(B,2.) - A*B;
3872 Cardano( a2, a1, a0, z1, z2, z3 );
3873
3874 // find stable roots
3875 vol1 = z1*R_CONST*Tk/Pbar;
3876 vol2 = z2*R_CONST*Tk/Pbar;
3877 vol3 = z3*R_CONST*Tk/Pbar;
3878 if (z1 > B)
3879 lnf1 = (-1.)*log(z1-B)
3880 - A/(B*sqrt(8.))*log((z1+(1.+sqrt(2.))*B)/(z1+(1.-sqrt(2.))*B))+z1-1.;
3881 else
3882 lnf1 = 1000.;
3883 if (z2 > B)
3884 lnf2 = (-1.)*log(z2-B)
3885 - A/(B*sqrt(8.))*log((z2+(1.+sqrt(2.))*B)/(z2+(1.-sqrt(2.))*B))+z2-1.;
3886 else
3887 lnf2 = 1000.;
3888 if (z3 > B)
3889 lnf3 = (-1.)*log(z3-B)
3890 - A/(B*sqrt(8.))*log((z3+(1.+sqrt(2.))*B)/(z3+(1.-sqrt(2.))*B))+z3-1.;
3891 else
3892 lnf3 = 1000.;
3893
3894 if (lnf2 < lnf1)
3895 {
3896 zmix = z2; vmix = vol2; lnf = lnf2;
3897 }
3898 else
3899 {
3900 zmix = z1; vmix = vol1; lnf = lnf1;
3901 }
3902 if (lnf3 < lnf)
3903 {
3904 zmix = z3; vmix = vol3; lnf = lnf3;
3905 }
3906 else
3907 {
3908 zmix = zmix; vmix = vmix; lnf = lnf;
3909 }
3910 fugmix = exp(lnf);
3911 PhVol = vmix;
3912
3913 return 0;
3914 }
3915
3916
3917
3918 /// Calculates fugacities and activities of fluid species in the mixture,
FugacitySpec(double * fugpure)3919 long int TPR78calc::FugacitySpec( double *fugpure )
3920 {
3921 long int i, j, iRet=0;
3922 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0., sum=0.;
3923 double A, B, lnfci, fci;
3924
3925 // Reload params to Pureparm
3926 for( j=0; j<NComp; j++ )
3927 {
3928 Fugpure[j][0] = fugpure[j]/Pbar;
3929 }
3930
3931 // retrieve properties of the mixture
3932 iRet = MixParam( amix, bmix );
3933 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix );
3934 A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
3935 B = bmix*Pbar/(R_CONST*Tk);
3936
3937 // calculate fugacity coefficient, fugacity and activity of species i
3938 for (i=0; i<NComp; i++)
3939 {
3940 sum = 0.;
3941 for (j=0; j<NComp; j++)
3942 {
3943 sum = sum + x[j]*AA[i][j];
3944 }
3945 lnfci = Pureparm[i][1]/bmix*(zmix-1.) - log(zmix-B)
3946 + A/(sqrt(8.)*B)*(2.*sum/amix-Pureparm[i][1]/bmix)
3947 * log((zmix+B*(1.-sqrt(2.)))/(zmix+B*(1.+sqrt(2.))));
3948 fci = exp(lnfci);
3949 Fugci[i][0] = fci; // fugacity coefficient using engineering convention
3950 Fugci[i][1] = x[i]*fci; // fugacity coefficient using geology convention
3951 Fugci[i][2] = Fugci[i][1]/Fugpure[i][0]; // activity of species
3952 if (x[i]>1.0e-20)
3953 Fugci[i][3] = Fugci[i][2]/x[i]; // activity coefficient of species
3954 else
3955 Fugci[i][3] = 1.0;
3956 }
3957
3958 return iRet;
3959 }
3960
3961
3962
3963 /// calculates residual functions in the mixture
ResidualFunct(double * fugpure)3964 long int TPR78calc::ResidualFunct( double *fugpure )
3965 {
3966 long int i, j, iRet=0;
3967 double fugmix=0., zmix=0., vmix=0., amix=0., bmix=0.;
3968 double /*A,*/ B, K, dK, d2K, Q, dQ, d2Q, damix, d2amix, ai, aj, dai, daj, d2ai, d2aj,
3969 cv, dPdT, dPdV, dVdT;
3970
3971 // Reload params to Pureparm (probably now obsolete?)
3972 for( j=0; j<NComp; j++ )
3973 {
3974 Fugpure[j][0] = fugpure[j]/Pbar;
3975 }
3976
3977 // retrieve properties of the mixture
3978 iRet = MixParam( amix, bmix );
3979 iRet = FugacityMix( amix, bmix, fugmix, zmix, vmix );
3980 //A = amix*Pbar/(pow(R_CONST, 2.)*pow(Tk, 2.));
3981 B = bmix*Pbar/(R_CONST*Tk);
3982
3983 // calculate total state functions of the mixture
3984 damix = 0.;
3985 d2amix = 0.;
3986 for (i=0; i<NComp; i++)
3987 {
3988 for (j=0; j<NComp; j++)
3989 {
3990 // pull parameters
3991 ai = Pureparm[i][0];
3992 aj = Pureparm[j][0];
3993 dai = Pureparm[i][2];
3994 daj = Pureparm[j][2];
3995 d2ai = Pureparm[i][3];
3996 d2aj = Pureparm[j][3];
3997 K = KK[i][j];
3998 dK = dKK[i][j];
3999 d2K = d2KK[i][j];
4000
4001 // increments to derivatives
4002 Q = sqrt(ai*aj);
4003 dQ = 0.5*( sqrt(aj/ai)*dai + sqrt(ai/aj)*daj );
4004 d2Q = 0.5*( dai*daj/sqrt(ai*aj) + d2ai*sqrt(aj)/sqrt(ai) + d2aj*sqrt(ai)/sqrt(aj)
4005 - 0.5*( pow(dai,2.)*sqrt(aj)/sqrt(pow(ai,3.))
4006 + pow(daj,2.)*sqrt(ai)/sqrt(pow(aj,3.)) ) );
4007 damix = damix + x[i]*x[j] * ( dQ*(1.-K) - Q*dK );
4008 d2amix = d2amix + x[i]*x[j] * ( d2Q*(1.-K) - 2.*dQ*dK - Q*d2K );
4009 }
4010 }
4011
4012 // calculate thermodynamic properties
4013 Grs = (amix/(R_CONST*Tk*sqrt(8.)*bmix) * log((vmix+(1.-sqrt(2.))*bmix)
4014 / (vmix+(1.+sqrt(2.))*bmix))-log(zmix*(1.-bmix/vmix))+zmix-1.)*R_CONST*Tk;
4015 Hrs = ((amix-Tk*damix)/(R_CONST*Tk*sqrt(8.)*bmix)*log((vmix+(1.-sqrt(2.))
4016 *bmix)/(vmix+(1.+sqrt(2.))*bmix))+zmix-1.)*R_CONST*Tk;
4017 Srs = (Hrs - Grs)/Tk;
4018
4019 // heat capacity part
4020 cv = Tk*d2amix/(bmix*sqrt(8.))
4021 * log( (zmix+B*(1.+sqrt(2.)))/(zmix+B*(1.-sqrt(2.))) );
4022 dPdT = R_CONST/(vmix-bmix) - damix/( vmix*(vmix+bmix) + bmix*(vmix-bmix) );
4023 dPdV = - R_CONST*Tk/pow((vmix-bmix),2.) + 2.*amix*(vmix+bmix)/pow((vmix*(vmix+bmix)+bmix*(vmix-bmix)),2.);
4024 dVdT = (-1.)*(1./dPdV)*dPdT;
4025 CPrs = cv + Tk*dPdT*dVdT - R_CONST;
4026 Vrs = vmix;
4027
4028 return iRet;
4029 }
4030
4031
4032
4033 #ifndef IPMGEMPLUGIN_
4034
4035 /// Calculates properties of pure fluids when called from DCthermo
PR78CalcFugPure(double Tmin,double * Cpg,double * FugProps)4036 long int TPR78calc::PR78CalcFugPure( double Tmin, double *Cpg, double *FugProps )
4037 {
4038 long int retCode = 0;
4039 double Coeff[7];
4040
4041 for( int ii=0; ii<7; ii++ )
4042 Coeff[ii] = (double)Cpg[ii];
4043
4044 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
4045 {
4046 retCode = FugacityPT( 0, Coeff );
4047 for( int i=0; i<6; i++ )
4048 FugProps[i] = Fugpure[0][i];
4049 return retCode;
4050 }
4051
4052 else
4053 {
4054 for( int i=1; i<6; i++ )
4055 FugProps[i] = 0.;
4056 FugProps[0] = 1.;
4057 FugProps[4] = 8.31451*Tk/Pbar;
4058 return -1;
4059 }
4060 }
4061
4062 #endif
4063
4064
4065
4066
4067
4068 //=======================================================================================================
4069 // Compensated Redlich-Kwong (CORK) model for fluid mixtures
4070 // References: Holland and Powell (1991)
4071 // (c) TW May 2010
4072 //=======================================================================================================
4073
4074
4075 // Constructor
TCORKcalc(long int NCmp,double Pp,double Tkp,char Eos_Code)4076 TCORKcalc::TCORKcalc( long int NCmp, double Pp, double Tkp, char Eos_Code ):
4077 TSolMod( NCmp, '8', Tkp, Pp )
4078 {
4079 RR = 0.00831451; // gas constant in kbar
4080 Pkb = Pbar/1000.; // pressure in kbar
4081 // aGEX = 0;
4082 // aVol = 0;
4083 Pparc = 0;
4084 alloc_internal();
4085 EosCode[0] = Eos_Code; // must be changed in m_dcomp.cpp line 818
4086
4087 }
4088
4089
4090
TCORKcalc(SolutionData * sd)4091 TCORKcalc::TCORKcalc( SolutionData *sd ):
4092 TSolMod( sd )
4093 {
4094 RR = 0.00831451; // gas constant in kbar
4095 Pkb = Pbar/1000.; // pressure in kbar
4096 Pparc = aPparc;
4097 // aGEX = arGEX;
4098 // aVol = arVol;
4099 alloc_internal();
4100
4101 for( long int j=0; j<NComp; j++ )
4102 // EosCode[j] = sd->DC_Codes[j]; // may be change to EosCode[j] = sd->TP_Code[j][3]; SD
4103 EosCode[j] = sd->TP_Code[j][3];
4104 }
4105
4106
4107
~TCORKcalc()4108 TCORKcalc::~TCORKcalc()
4109 {
4110 free_internal();
4111 }
4112
4113
4114
4115 /// allocate work arrays for pure fluid and fluid mixture properties
alloc_internal()4116 void TCORKcalc::alloc_internal()
4117 {
4118 EosCode = new char[NComp];
4119 phi = new double [NComp];
4120 dphi = new double [NComp];
4121 d2phi = new double [NComp];
4122 dphip = new double [NComp];
4123 Eosparm = new double [NComp][2];
4124 Fugpure = new double [NComp][6];
4125 Fugci = new double [NComp][4];
4126 Rho = new double [NComp][11];
4127 A = new double *[NComp];
4128 W = new double *[NComp];
4129 B = new double *[NComp];
4130 dB = new double *[NComp];
4131 d2B = new double *[NComp];
4132 dBp = new double *[NComp];
4133
4134 for (long int i=0; i<NComp; i++)
4135 {
4136 A[i] = new double [NComp];
4137 W[i] = new double [NComp];
4138 B[i] = new double [NComp];
4139 dB[i] = new double [NComp];
4140 d2B[i] = new double [NComp];
4141 dBp[i] = new double [NComp];
4142 }
4143 }
4144
4145
4146
free_internal()4147 void TCORKcalc::free_internal()
4148 {
4149 long int i;
4150
4151 for (i=0; i<NComp; i++)
4152 {
4153 delete[]A[i];
4154 delete[]W[i];
4155 delete[]B[i];
4156 delete[]dB[i];
4157 delete[]d2B[i];
4158 delete[]dBp[i];
4159 }
4160
4161 delete[]EosCode;
4162 delete[]phi;
4163 delete[]dphi;
4164 delete[]d2phi;
4165 delete[]dphip;
4166 delete[]Eosparm;
4167 delete[]Fugpure;
4168 delete[]Fugci;
4169 delete[]Rho;
4170 delete[]A;
4171 delete[]W;
4172 delete[]B;
4173 delete[]dB;
4174 delete[]d2B;
4175 delete[]dBp;
4176 }
4177
4178
4179
4180 /// high-level method to retrieve pure fluid fugacities
PureSpecies()4181 long int TCORKcalc::PureSpecies()
4182 {
4183 long int j, retCode = 0;
4184
4185 for( j=0; j<NComp; j++ )
4186 {
4187 // Calling CORK EoS for pure fugacity
4188 retCode = FugacityPT( j, aDCc+j*NP_DC );
4189 aGEX[j] = log( Fugpure[j][0] );
4190 Pparc[j] = Fugpure[j][0]*Pbar; // fure fluid fugacity (required for performance)
4191 aVol[j] = Fugpure[j][4]*10.; // molar volume of pure fluid component, J/bar to cm3
4192 } // j
4193
4194 if ( retCode )
4195 {
4196 char buf[150];
4197 sprintf(buf, "CORK fluid: calculation of pure fluid fugacity failed");
4198 Error( "E71IPM IPMgamma: ", buf );
4199 }
4200
4201 return 0;
4202 }
4203
4204
4205
4206 /// high-level method to calculate T,P corrected binary interaction parameters
PTparam()4207 long int TCORKcalc::PTparam()
4208 {
4209 long int j, i, ip, i1, i2;
4210 double a;
4211
4212 Pkb = Pbar/1000.;
4213
4214 PureSpecies();
4215
4216 // set all interaction parameters zero
4217 for( j=0; j<NComp; j++ )
4218 {
4219 for( i=0; i<NComp; i++ )
4220 {
4221 A[j][i] = 0.;
4222 W[j][i] = 0.;
4223 B[j][i] = 0.;
4224 dB[j][i] = 0.;
4225 d2B[j][i] = 0.;
4226 dBp[j][i] = 0.;
4227 }
4228 }
4229
4230 // transfer interaction parameters that have non-standard value
4231 if( NPcoef > 0 )
4232 {
4233 for ( ip=0; ip<NPar; ip++ )
4234 {
4235 i1 = aIPx[MaxOrd*ip];
4236 i2 = aIPx[MaxOrd*ip+1];
4237 a = aIPc[NPcoef*ip];
4238 A[i1][i2] = a;
4239 A[i2][i1] = a; // symmetric case
4240 }
4241 }
4242
4243 return 0;
4244 }
4245
4246
4247
4248 /// high-level method to retrieve activity coefficients of the fluid mixture
MixMod()4249 long int TCORKcalc::MixMod()
4250 {
4251 long int i, j, k;
4252 double dj, dk, sumphi, lnGam, Gam, vi, vj, vk;
4253
4254 // calculate phi values
4255 sumphi = 0.;
4256 for (i=0; i<NComp; i++)
4257 {
4258 vi = Fugpure[i][4];
4259 sumphi = sumphi + x[i]*vi;
4260 }
4261 for (i=0; i<NComp; i++)
4262 {
4263 vi = Fugpure[i][4];
4264 phi[i] = x[i]*vi/sumphi;
4265 }
4266
4267 // interaction parameters
4268 for (i=0; i<NComp; i++)
4269 {
4270 for (j=i+1; j<NComp; j++)
4271 {
4272 vi = Fugpure[i][4];
4273 vj = Fugpure[j][4];
4274 W[i][j] = A[i][j]*(vi+vj)/(vi*vj);
4275 }
4276 }
4277
4278 // activity coefficients
4279 for (i=0; i<NComp; i++)
4280 {
4281 lnGam = 0.;
4282 for (j=0; j<NComp; j++)
4283 {
4284 for (k=j+1; k<NComp; k++)
4285 {
4286 vi = Fugpure[i][4];
4287 vj = Fugpure[j][4];
4288 vk = Fugpure[k][4];
4289 if (i==j)
4290 dj = 1.;
4291 else
4292 dj = 0.;
4293 if (i==k)
4294 dk = 1.;
4295 else
4296 dk = 0.;
4297 lnGam = lnGam - (dj-phi[j])*(dk-phi[k])*W[j][k]*2.*vi/(vj+vk);
4298 }
4299 }
4300 Gam = exp(lnGam/(R_CONST*Tk));
4301 Fugci[i][0] = Gam;
4302 if( Fugci[i][0] > 1e-23 )
4303 lnGamma[i] = log( Fugci[i][0] );
4304 else
4305 lnGamma[i] = 0.;
4306
4307 } // i
4308
4309 return 0;
4310 }
4311
4312
4313
4314 /// high-level method to retrieve residual functions of the fluid mixture
ExcessProp(double * Zex)4315 long int TCORKcalc::ExcessProp( double *Zex )
4316 {
4317 long int iRet;
4318
4319 iRet = ResidualFunct();
4320
4321 if ( iRet )
4322 {
4323 char buf[150];
4324 sprintf(buf, "CORK fluid: calculation failed");
4325 Error( "E71IPM IPMgamma: ", buf );
4326 }
4327
4328 Ars = Grs - Vrs*Pbar;
4329 Urs = Hrs - Vrs*Pbar;
4330
4331 // assignments (residual functions)
4332 Zex[0] = Grs;
4333 Zex[1] = Hrs;
4334 Zex[2] = Srs;
4335 Zex[3] = CPrs;
4336 Zex[4] = Vrs;
4337 Zex[5] = Ars;
4338 Zex[6] = Urs;
4339
4340 return iRet;
4341 }
4342
4343
4344
4345 /// calculates ideal mixing properties
IdealProp(double * Zid)4346 long int TCORKcalc::IdealProp( double *Zid )
4347 {
4348 long int j;
4349 double s, sc, sp;
4350
4351 s = 0.0;
4352 for ( j=0; j<NComp; j++ )
4353 {
4354 if ( x[j] > 1.0e-32 )
4355 s += x[j]*log(x[j]);
4356 }
4357 sc = (-1.)*R_CONST*s;
4358 sp = (-1.)*R_CONST*log(Pbar);
4359 Hid = 0.0;
4360 CPid = 0.0;
4361 Vid = 0.0;
4362 Sid = sc + sp;
4363 Gid = Hid - Sid*Tk;
4364 Aid = Gid - Vid*Pbar;
4365 Uid = Hid - Vid*Pbar;
4366
4367 // assignments (ideal mixing properties)
4368 Zid[0] = Gid;
4369 Zid[1] = Hid;
4370 Zid[2] = Sid;
4371 Zid[3] = CPid;
4372 Zid[4] = Vid;
4373 Zid[5] = Aid;
4374 Zid[6] = Uid;
4375
4376 return 0;
4377 }
4378
4379
4380
4381 /// high-level method to retrieve pure fluid properties
FugacityPT(long int j,double * EoSparam)4382 long int TCORKcalc::FugacityPT( long int j, double *EoSparam )
4383 {
4384 long int iErr = 0;
4385
4386 // reads EoS parameters from database into work array
4387 if( !EoSparam )
4388 return -1; // Memory alloc error
4389 Eosparm[j][0] = EoSparam[0]; // critical temperature in K
4390 Eosparm[j][1] = EoSparam[1]; // critical pressure in bar
4391
4392 // select subroutines for different fluid types
4393 switch ( EosCode[j] )
4394 {
4395 // case DC_GAS_H2O_: // H2O
4396 case CEM_H2O_: // H2O
4397 iErr = FugacityH2O( j );
4398 break;
4399 // case DC_GAS_CO2_: // CO2
4400 case CEM_CO2_: // CO2
4401 iErr = FugacityCO2( j );
4402 break;
4403 // case DC_GAS_COMP_: // other fluids
4404 // case DC_GAS_H2_:
4405 // case DC_GAS_N2_:
4406 case CEM_GAS_: // other fluids
4407 case CEM_CH4_:
4408 case CEM_N2_:
4409 case CEM_H2_:
4410 case CEM_O2_:
4411 case CEM_AR_:
4412 case CEM_PO_:
4413 case CEM_NP_:
4414 iErr = FugacityCorresponding( j );
4415 break;
4416 default:
4417 iErr = 3;
4418 break;
4419 }
4420
4421 return iErr;
4422 }
4423
4424
4425
4426 /// calculates fugacity and state functions of H2O
FugacityH2O(long int j)4427 long int TCORKcalc::FugacityH2O( long int j )
4428 {
4429 long int phState; // 1: vapor, 2: liquid
4430 double a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a, da, d2a, b, c0, c1, c,
4431 d0, d1, d, e0, e, /*dc, dd,*/ p0, Psat, vc, fc, v, g1, g2, g3, rho, z;
4432 double grs, hrs, srs, cprs, cvrs, dpr, dprr, dpt, dptt, dprt,
4433 drt, drtt, drp, drpp, drtp, dpv, dvt;
4434
4435 a0 = 1113.4;
4436 a1 = -0.88517;
4437 a2 = 4.53e-3;
4438 a3 = -1.3183e-5;
4439 a4 = -0.22291;
4440 a5 = -3.8022e-4;
4441 a6 = 1.7791e-7;
4442 a7 = 5.8487;
4443 a8 = -2.1370e-2;
4444 a9 = 6.8133e-5;
4445 b = 1.465;
4446 c0 = -8.909e-2;
4447 c1 = 0.;
4448 d0 = 1.9853e-3;
4449 d1 = 0.;
4450 e0 = 8.0331e-2;
4451 c = c0 + c1*Tk;
4452 d = d0 + d1*Tk;
4453 e = e0;
4454 //dc = c1;
4455 //dd = d1;
4456 p0 = 2.;
4457
4458 // molar volume and fugacity coefficient
4459 if (Tk > 695.) // supercritical fluid phase
4460 {
4461 phState = 1;
4462 a = a0 + a4*(Tk-673.) + a5*pow((Tk-673.),2.) + a6*pow((Tk-673.),3.);
4463 da = -a4 - 2.*a5*(673.-Tk) - 3.*a6*pow((673.-Tk),2.);
4464 d2a = 2.*a5 + 6.*a6*(673.-Tk);
4465 VolumeFugacity(phState, Pkb, p0, a, b, c, d, e, vc, fc);
4466 }
4467
4468 else // subcritical region
4469 {
4470 Psat = (-13.627e-3) + (7.29395e-7)*pow(Tk,2.) - (2.34622e-9)*pow(Tk,3.) + (4.83607e-15)*pow(Tk,5.);
4471 if (Pkb < Psat)
4472 {
4473 phState = 1; // vapor
4474 if (Tk < 673.)
4475 {
4476 a = a0 + a7*(673.-Tk) + a8*pow((673.-Tk),2.) + a9*pow((673.-Tk),3.);
4477 da = -a7 - 2.*a8*(673.-Tk) - 3.*a9*pow((673.-Tk),2.);
4478 d2a = 2.*a8 + 6.*a9*(673.-Tk);
4479 }
4480 else
4481 {
4482 a = a0 + a4*(673.-Tk) + a5*pow((673.-Tk),2.) + a6*pow((673.-Tk),3.);
4483 da = -a4 - 2.*a5*(673.-Tk) - 3.*a6*pow((673.-Tk),2.);
4484 d2a = 2.*a5 + 6.*a6*(673.-Tk);
4485 }
4486 VolumeFugacity(phState, Pkb, p0, a, b, c, d, e, vc, fc);
4487 }
4488
4489 else
4490 {
4491 phState = 1; // gaseous phase at Psat
4492 if (Tk < 673.)
4493 {
4494 a = a0 + a7*(673.-Tk) + a8*pow((673.-Tk),2.) + a9*pow((673.-Tk),3.);
4495 da = -a7 - 2.*a8*(673.-Tk) - 3.*a9*pow((673.-Tk),2.);
4496 d2a = 2.*a8 + 6.*a9*(673.-Tk);
4497 }
4498 else
4499 {
4500 a = a0 + a4*(673.-Tk) + a5*pow((673.-Tk),2.) + a6*pow((673.-Tk),3.);
4501 da = -a4 - 2.*a5*(673.-Tk) - 3.*a6*pow((673.-Tk),2.);
4502 d2a = 2.*a5 + 6.*a6*(673.-Tk);
4503 }
4504 VolumeFugacity(phState, Psat, p0, a, b, c, d, e, v, g1);
4505
4506 phState = 2; // liquid phase at Psat
4507 if (Tk < 673.)
4508 {
4509 a = a0 + a1*(673.-Tk) + a2*pow((673.-Tk),2.) + a3*pow((673.-Tk),3.);
4510 da = -a1 - 2.*a2*(673.-Tk) - 3.*a3*pow((673.-Tk),2.);
4511 d2a = 2.*a2 + 6.*a3*(673.-Tk);
4512 }
4513 else
4514 {
4515 a = a0 + a4*(673.-Tk) + a5*pow((673.-Tk),2.) + a6*pow((673.-Tk),3.);
4516 da = -a4 - 2.*a5*(673.-Tk) - 3.*a6*pow((673.-Tk),2.);
4517 d2a = 2.*a5 + 6.*a6*(673.-Tk);
4518 }
4519 VolumeFugacity(phState, Psat, p0, a, b, c, d, e, v, g2);
4520
4521 phState = 2; // fluid phase at Pkb
4522 if (Tk < 673.)
4523 {
4524 a = a0 + a1*(673.-Tk) + a2*pow((673.-Tk),2.) + a3*pow((673.-Tk),3.);
4525 da = -a1 - 2.*a2*(673.-Tk) - 3.*a3*pow((673.-Tk),2.);
4526 d2a = 2.*a2 + 6.*a3*(673.-Tk);
4527 }
4528 else
4529 {
4530 a = a0 + a4*(673.-Tk) + a5*pow((673.-Tk),2.) + a6*pow((673.-Tk),3.);
4531 da = -a4 - 2.*a5*(673.-Tk) - 3.*a6*pow((673.-Tk),2.);
4532 d2a = 2.*a5 + 6.*a6*(673.-Tk);
4533 }
4534 VolumeFugacity(phState, Pkb, p0, a, b, c, d, e, vc, g3); // corrected from (phState, Pkb, psat ...)
4535 fc = g1/g2*g3; // this formula needs checking
4536 }
4537 }
4538
4539 // pressure and density derivatives (MRK part)
4540 rho = (1./vc);
4541 dpr = (RR*Tk)/(pow((1./rho-b),2.)*pow(rho,2.)) - a/((1./rho+b)*pow(Tk,0.5))
4542 - a/(rho*pow((1./rho+b),2.)*pow(Tk,0.5));
4543 dprr = (2.*RR*Tk)/(pow((1./rho-b),3.)*pow(rho,4.)) - (2.*RR*Tk)/(pow((1./rho-b),2.)*pow(rho,3.))
4544 - (2.*a)/(pow(rho,3.)*pow((1./rho+b),3.)*pow(Tk,0.5));
4545 dpt = RR/(1./rho-b) - (da*rho)/((1./rho+b)*pow(Tk,0.5)) + (0.5*a*rho)/((1./rho+b)*pow(Tk,1.5));
4546 dptt = - (d2a*rho)/((1./rho+b)*pow(Tk,0.5)) + (da*rho)/((1./rho+b)*pow(Tk,1.5))
4547 - (0.75*a*rho)/((1./rho+b)*pow(Tk,2.5));
4548 dprt = RR/(pow((1./rho-b),2.)*pow(rho,2.)) - da/((1./rho+b)*pow(Tk,0.5))
4549 - da/(rho*pow((1./rho+b),2.)*pow(Tk,0.5)) + (0.5*a)/((1./rho+b)*pow(Tk,1.5))
4550 + (0.5*a)/(rho*pow((1./rho+b),2.)*pow(Tk,1.5));
4551 drp = (1./dpr);
4552 drpp = (-1.)*dprr*pow(dpr,-3.);
4553 drt = (-1.)*(1./dpr)*dpt;
4554 drtt = (dpt - dptt + dpt*(dprt-dpr)*pow(dpr,-1.))*pow(dpr,-1.)
4555 + (dprt - dpt*dprr*pow(dpr,-1.))*pow(dpr,-2.)*dpt;
4556 drtp = - pow(dpr,-1.) - (dprt-dpr)*pow(dpr,-2.) + dprr*dpt*pow(dpr,-3.);
4557
4558 // thermodynamic properties (MRK part)
4559 z = (Pkb*vc)/(RR*Tk);
4560 grs = R_CONST*Tk*log(fc);
4561 hrs = ( (z-1.) - 1./(b*RR*Tk)*(pow(Tk,0.5)*da-(3./2.)*a/pow(Tk,0.5))*log(vc/(vc+b)) )*R_CONST*Tk;
4562 srs = (hrs-grs)/Tk;
4563 cvrs = ( -1./(b*RR)*(-da/pow(Tk,0.5)+pow(Tk,0.5)*d2a+(3.*a)/(4.*pow(Tk,1.5))))*log(vc/(vc+b) )*R_CONST;
4564 dpv = (- RR*Tk/pow((vc-b),2.) + a/(pow(vc,2.)*(vc+b)*pow(Tk,0.5))
4565 + a/(vc*pow((vc+b),2.)*pow(Tk,0.5)))*(1000.);
4566 dvt = (-1.)*(1./dpv)*(dpt*1000.);
4567 cprs = cvrs + Tk*(dpt*1000.)*dvt - R_CONST;
4568
4569 // copy results
4570 Fugpure[j][0] = fc; // fugacity coefficient
4571 Fugpure[j][1] = grs;
4572 Fugpure[j][2] = hrs;
4573 Fugpure[j][3] = srs;
4574 Fugpure[j][4] = vc;
4575 Fugpure[j][5] = cprs;
4576 Rho[j][0] = rho*(0.1); // mol cm-3
4577 Rho[j][1] = drt*(0.1); // mol cm-3 K-1
4578 Rho[j][2] = drtt*(0.1); // mol cm-3 K-2
4579 Rho[j][3] = drp*(0.001); // mol cm-3 MPa-1
4580 Rho[j][4] = drpp*(0.00001); // mol cm-3 MPa-2
4581 Rho[j][5] = drtp*(0.001); // mol cm-3 K-1 MPa-1
4582 Rho[j][6] = dpr*(1000.); // MPa cm3 mol-1
4583 Rho[j][7] = dprr*(10000.); // MPa cm6 mol-2
4584 Rho[j][8] = dpt*(100.); // MPa K-1
4585 Rho[j][9] = dptt*(100.); // MPa K-2
4586 Rho[j][10] = dprt*(1000.); // MPa cm3 mol-1 K-1
4587
4588 return 0;
4589 }
4590
4591
4592
4593 /// calculates fugacity and state functions of CO2
FugacityCO2(long int j)4594 long int TCORKcalc::FugacityCO2( long int j )
4595 {
4596 long int phState;
4597 double a0, a1, a2, a, b, c0, c1, c, d0, d1, d, e0, e, da, d2a, p0, vc, fc, rho, z;
4598 double grs, hrs, srs, cprs, cvrs, dpr, dprr, dpt, dptt, dprt,
4599 drt, drtt, drp, drpp, drtp, dpv, dvt;
4600
4601 a0 = 659.8;
4602 a1 = 0.21078;
4603 a2 = -6.3976e-4;
4604 b = 3.057;
4605 c0 = -1.78198e-1;
4606 c1 = 2.45317e-5;
4607 d0 = 5.40776e-3;
4608 d1 = -1.59046e-6;
4609 e0 = 0.;
4610 a = a0 + a1*Tk + a2*pow(Tk,2.);
4611 da = a1 + 2.*a2*Tk;
4612 d2a = 2.*a2;
4613 c = c0 + c1*Tk;
4614 d = d0 + d1*Tk;
4615 e = e0;
4616 p0 = 5.;
4617
4618 // molar volume and fugacity coefficient
4619 phState = 1;
4620 VolumeFugacity(phState, Pkb, p0, a, b, c, d, e, vc, fc);
4621
4622 // pressure and density derivatives (MRK part)
4623 rho = (1./vc);
4624 dpr = (RR*Tk)/(pow((1./rho-b),2.)*pow(rho,2.)) - a/((1./rho+b)*pow(Tk,0.5))
4625 - a/(rho*pow((1./rho+b),2.)*pow(Tk,0.5));
4626 dprr = (2.*RR*Tk)/(pow((1./rho-b),3.)*pow(rho,4.)) - (2.*RR*Tk)/(pow((1./rho-b),2.)*pow(rho,3.))
4627 - (2.*a)/(pow(rho,3.)*pow((1./rho+b),3.)*pow(Tk,0.5));
4628 dpt = RR/(1./rho-b) - (da*rho)/((1./rho+b)*pow(Tk,0.5)) + (0.5*a*rho)/((1./rho+b)*pow(Tk,1.5));
4629 dptt = - (d2a*rho)/((1./rho+b)*pow(Tk,0.5)) + (da*rho)/((1./rho+b)*pow(Tk,1.5))
4630 - (0.75*a*rho)/((1./rho+b)*pow(Tk,2.5));
4631 dprt = RR/(pow((1./rho-b),2.)*pow(rho,2.)) - da/((1./rho+b)*pow(Tk,0.5))
4632 - da/(rho*pow((1./rho+b),2.)*pow(Tk,0.5)) + (0.5*a)/((1./rho+b)*pow(Tk,1.5))
4633 + (0.5*a)/(rho*pow((1./rho+b),2.)*pow(Tk,1.5));
4634 drp = (1./dpr);
4635 drpp = (-1.)*dprr*pow(dpr,-3.);
4636 drt = (-1.)*(1./dpr)*dpt;
4637 drtt = (dpt - dptt + dpt*(dprt-dpr)*pow(dpr,-1.))*pow(dpr,-1.)
4638 + (dprt - dpt*dprr*pow(dpr,-1.))*pow(dpr,-2.)*dpt;
4639 drtp = - pow(dpr,-1.) - (dprt-dpr)*pow(dpr,-2.) + dprr*dpt*pow(dpr,-3.);
4640
4641 // thermodynamic properties (MRK part)
4642 z = (Pkb*vc)/(RR*Tk);
4643 grs = R_CONST*Tk*log(fc);
4644 hrs = ( (z-1.) - 1./(b*RR*Tk)*(pow(Tk,0.5)*da-(3./2.)*a/pow(Tk,0.5))*log(vc/(vc+b)) )*R_CONST*Tk;
4645 srs = (hrs-grs)/Tk;
4646 cvrs = ( -1./(b*RR)*(-da/pow(Tk,0.5)+pow(Tk,0.5)*d2a+(3.*a)/(4.*pow(Tk,1.5))))*log(vc/(vc+b) )*R_CONST;
4647 dpv = (- RR*Tk/pow((vc-b),2.) + a/(pow(vc,2.)*(vc+b)*pow(Tk,0.5))
4648 + a/(vc*pow((vc+b),2.)*pow(Tk,0.5)))*(1000.);
4649 dvt = (-1.)*(1./dpv)*(dpt*1000.);
4650 cprs = cvrs + Tk*(dpt*1000.)*dvt - R_CONST;
4651
4652 // copy results
4653 Fugpure[j][0] = fc; // fugacity coefficient
4654 Fugpure[j][1] = grs;
4655 Fugpure[j][2] = hrs;
4656 Fugpure[j][3] = srs;
4657 Fugpure[j][4] = vc;
4658 Fugpure[j][5] = cprs;
4659 Rho[j][0] = rho*(0.1); // mol cm-3
4660 Rho[j][1] = drt*(0.1); // mol cm-3 K-1
4661 Rho[j][2] = drtt*(0.1); // mol cm-3 K-2
4662 Rho[j][3] = drp*(0.001); // mol cm-3 MPa-1
4663 Rho[j][4] = drpp*(0.00001); // mol cm-3 MPa-2
4664 Rho[j][5] = drtp*(0.001); // mol cm-3 K-1 MPa-1
4665 Rho[j][6] = dpr*(1000.); // MPa cm3 mol-1
4666 Rho[j][7] = dprr*(10000.); // MPa cm6 mol-2
4667 Rho[j][8] = dpt*(100.); // MPa K-1
4668 Rho[j][9] = dptt*(100.); // MPa K-2
4669 Rho[j][10] = dprt*(1000.); // MPa cm3 mol-1 K-1
4670
4671 return 0;
4672 }
4673
4674
4675
4676 /// calculates fugacity and state functions of fluids other than H2O and CO2
FugacityCorresponding(long int j)4677 long int TCORKcalc::FugacityCorresponding( long int j )
4678 {
4679 double a0, a1, a, b0, b, c0, c1, c, d0, d1, d, tcr, pcr, da, dc, dd, vc, fc, rtlnf, rho;
4680 double grs, hrs, srs, cprs, drt, drtt, drp, drpp, drtp,
4681 dvt, dvtt, dvp, dvpp, dvtp;
4682
4683 a0 = 5.45963e-5;
4684 a1 = -8.6392e-6;
4685 b0 = 9.18301e-4;
4686 c0 = -3.30558e-5;
4687 c1 = 2.30524e-6;
4688 d0 = 6.93054e-7;
4689 d1 = -8.38293e-8;
4690 tcr = Eosparm[j][0];
4691 pcr = Eosparm[j][1]/1000.; // kbar
4692 a = a0*pow(tcr,2.)*sqrt(tcr)/pcr + a1*tcr*sqrt(tcr)/pcr*Tk;
4693 b = b0*tcr/pcr;
4694 c = c0*tcr/(pcr*sqrt(pcr)) + c1/(pcr*sqrt(pcr))*Tk;
4695 d = d0*tcr/pow(pcr,2.) + d1/pow(pcr,2.)*Tk;
4696 da = a1*pow(tcr,1.5)/pcr;
4697 dc = c1/pow(pcr,1.5);
4698 dd = d1/pow(pcr,2.);
4699
4700 // molar volume and fugacity coefficient
4701 vc = RR*Tk/Pkb + b - a*RR*sqrt(Tk)/((RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb)) + c*sqrt(Pkb) + d*Pkb;
4702 rtlnf = RR*Tk*log(1000*Pkb) + b*Pkb + a/(b*sqrt(Tk))*(log(RR*Tk+b*Pkb) - log(RR*Tk+2*b*Pkb))
4703 + (2./3.)*c*Pkb*sqrt(Pkb) + (d/2.)*pow(Pkb,2.);
4704 fc = exp(rtlnf/(RR*Tk))/(1000*Pkb);
4705
4706 // volume and density derivatives
4707 dvt = ( RR/Pkb - (da*RR*pow(Tk,0.5))/((RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
4708 - (a*RR)/(2.*pow(Tk,0.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
4709 + (a*pow(RR,2.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4710 + (a*pow(RR,2.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4711 + dc*pow(Pkb,0.5) + dd*Pkb );
4712 dvtt = ( - (da*RR)/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
4713 + (2.*da*pow(RR,2.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4714 + (2.*da*pow(RR,2.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4715 + (a*RR)/(4.*pow(Tk,1.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
4716 + (a*pow(RR,2.))/(pow(Tk,0.5)*pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4717 + (a*pow(RR,2.))/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4718 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
4719 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
4720 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.)) );
4721 dvp = ( - RR*Tk/pow(Pkb,2.) + (a*RR*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4722 + (2.*a*RR*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4723 + c/(2.*pow(Pkb,0.5)) + d ) / (1000.);
4724 dvpp = ( (2.*RR*Tk)/pow(Pkb,3.) - (2.*a*RR*pow(Tk,0.5)*pow(b,2.))/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
4725 - (4.*a*RR*pow(Tk,0.5)*pow(b,2.))/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
4726 - (8.*a*RR*pow(Tk,0.5)*pow(b,2.))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.))
4727 - c/(4.*pow(Pkb,1.5)) ) / (1.0e6);
4728 dvtp = ( - RR/pow(Pkb,2.) + (da*RR*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4729 + (a*RR*b)/(2.*pow(Tk,0.5)*pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
4730 - (2.*a*pow(RR,2.)*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
4731 - (3.*a*pow(RR,2.)*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
4732 + (2.*da*RR*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4733 + (a*RR*b)/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
4734 - (4.*a*pow(RR,2.)*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.))
4735 + dc/(2.*pow(Pkb,2.)) + dd ) / (1000.);
4736 rho = (1./vc);
4737 drt = - dvt*pow(rho,2.);
4738 drp = - dvp*pow(rho,2.);
4739 drtt = - (dvtt*pow(rho,3.)-2.*pow(drt,2.))/rho;
4740 drpp = - (dvpp*pow(rho,3.)-2.*pow(drp,2.))/rho;
4741 drtp = - (dvtp*pow(rho,3.)-2.*drt*drp)/rho;
4742
4743 // thermodynamic properties
4744 grs = R_CONST*Tk*log(fc);
4745 hrs = ( (pow(Tk,0.5)*da*log(RR*Tk+2.*b*Pkb))/b - (pow(Tk,0.5)*da*log(RR*Tk+b*Pkb))/b
4746 - (3.*a*log(RR*Tk+2.*b*Pkb))/(2.*pow(Tk,0.5)*b) + (3.*a*log(RR*Tk+b*Pkb))/(2.*pow(Tk,0.5)*b)
4747 - (pow(Tk,0.5)*a*RR)/(b*(RR*Tk+b*Pkb)) + (pow(Tk,0.5)*a*RR)/(b*(RR*Tk+2.*b*Pkb))
4748 - (2.*Tk*dc*pow(Pkb,1.5))/3. - (Tk*dd*pow(Pkb,2.))/2. + b*Pkb
4749 + (2.*c*pow(Pkb,1.5))/3. + (d*pow(Pkb,2.))/2. )*1000.;
4750 srs = (hrs-grs)/Tk;
4751 cprs = ( - (2.*pow(Tk,0.5)*da*RR)/(b*(RR*Tk+b*Pkb)) + (pow(Tk,0.5)*a*pow(RR,2.))/(b*pow((RR*Tk+b*Pkb),2.))
4752 + (2.*pow(Tk,0.5)*da*RR)/(b*(RR*Tk+2.*b*Pkb)) - (a*RR)/(pow(Tk,0.5)*b*(RR*Tk+2.*b*Pkb))
4753 - (da*log(RR*Tk+2.*b*Pkb))/(pow(Tk,0.5)*b) + (da*log(RR*Tk+b*Pkb))/(pow(Tk,0.5)*b)
4754 + (3.*a*log(RR*Tk+2.*b*Pkb))/(4.*pow(Tk,1.5)*b) - (3.*a*log(RR*Tk+b*Pkb))/(4.*pow(Tk,1.5)*b)
4755 - (pow(Tk,0.5)*a*pow(RR,2.))/(b*pow((RR*Tk+2.*b*Pkb),2.)) + (a*RR)/(pow(Tk,0.5)*b*(RR*Tk+b*Pkb)) )*1000.;
4756
4757 // copy results
4758 Fugpure[j][0] = fc; // fugacity coefficient
4759 Fugpure[j][1] = grs;
4760 Fugpure[j][2] = hrs;
4761 Fugpure[j][3] = srs;
4762 Fugpure[j][4] = vc;
4763 Fugpure[j][5] = cprs;
4764 Rho[j][0] = rho*(0.1); // mol cm-3
4765 Rho[j][1] = drt*(0.1); // mol cm-3 K-1
4766 Rho[j][2] = drtt*(0.1); // mol cm-3 K-2
4767 Rho[j][3] = drp; // mol cm-3 MPa-1
4768 Rho[j][4] = drpp*(10.); // mol cm-3 MPa-2
4769 Rho[j][5] = drtp; // mol cm-3 K-1 MPa-1
4770 Rho[j][6] = 1.0; // MPa cm3 mol-1
4771 Rho[j][7] = 1.0; // MPa cm6 mol-2
4772 Rho[j][8] = 1.0; // MPa K-1
4773 Rho[j][9] = 1.0; // MPa K-2
4774 Rho[j][10] = 1.0; // MPa cm3 mol-1 K-1
4775
4776 return 0;
4777 }
4778
4779
4780
4781 /// calculate volume and fugacity coefficient
VolumeFugacity(long int phState,double pp,double p0,double a,double b,double c,double d,double e,double & vol,double & fc)4782 long int TCORKcalc::VolumeFugacity( long int phState, double pp, double p0, double a, double b, double c,
4783 double d, double e, double &vol, double &fc )
4784 {
4785 double cb, cc, cd, v1, v2, v3, vmrk, vvir, lng, lnvir;
4786
4787 cb = (-1.)*RR*Tk/pp;
4788 cc = (-1.)*(b*RR*Tk+pow(b,2.)*pp-a/sqrt(Tk))/pp;
4789 cd = (-1.)*a*b/(sqrt(Tk)*pp);
4790 Cardano(cb, cc, cd, v1, v2, v3);
4791
4792 if ( phState == 1 ) // vapor root
4793 {
4794 if (v1>0.)
4795 vmrk = v1;
4796 else if (v2>0.)
4797 vmrk = v2;
4798 else
4799 vmrk = v3;
4800 if (v2>vmrk)
4801 vmrk = v2;
4802 if (v3>vmrk)
4803 vmrk=v3;
4804 }
4805
4806 else // liquid root
4807 {
4808 if (v1>0.)
4809 vmrk = v1;
4810 else if (v2>0.)
4811 vmrk = v2;
4812 else
4813 vmrk = v3;
4814 if ( (v2<vmrk) && (v2>0.) )
4815 vmrk = v2;
4816 if ( (v3<vmrk) && (v2>0.) )
4817 vmrk = v3;
4818 }
4819
4820 // calculate fugacity coefficient
4821 lng = pp*vmrk/(RR*Tk) - 1. - log((vmrk-b)*pp/(RR*Tk)) - a/(b*RR*Tk*sqrt(Tk))*log(1.+b/vmrk);
4822 if (pp>p0)
4823 {
4824 vvir = c*sqrt(pp-p0) + d*(pp-p0) + e*sqrt(sqrt(pp-p0));
4825 lnvir = ((2./3.)*c*(pp-p0)*sqrt(pp-p0) + d/2.*pow((pp-p0),2.) + (4./5.)*e*(pp-p0)*sqrt(sqrt(pp-p0)))/(RR*Tk);
4826 }
4827 else
4828 {
4829 vvir = 0.;
4830 lnvir = 0.;
4831 }
4832
4833 vol = vmrk + vvir;
4834 lng = lng + lnvir;
4835 fc = exp(lng);
4836 return 0;
4837 }
4838
4839
4840
4841 /// finds roots of cubic equation
Cardano(double cb,double cc,double cd,double & v1,double & v2,double & v3)4842 long int TCORKcalc::Cardano(double cb, double cc, double cd, double &v1, double &v2, double &v3)
4843 {
4844 double co, cp, cq, cr, cp2, cq3;
4845
4846 cp = (2.*pow(cb,3.)-9.*cb*cc+27.*cd)/54.;
4847 cq = (pow(cb,2.)-3.*cc)/9.;
4848 cp2 = pow(cp,2.);
4849 cq3 = pow(cq,3.);
4850
4851 if ( (cp2-cq3)>0. ) // one real root
4852 {
4853 cr = -cp/fabs(cp)*pow(fabs(cp)+sqrt(cp2-cq3),1./3.);
4854 if (cr!=0)
4855 {
4856 v1 = cr + cq/cr - cb/3.;
4857 v2 = cr + cq/cr - cb/3.;
4858 v3 = cr + cq/cr - cb/3.;
4859 }
4860 else
4861 {
4862 v1 = -cb/3.;
4863 v2 = -cb/3.;
4864 v3 = -cb/3.;
4865 }
4866 }
4867
4868 else //three real roots
4869 {
4870 co = atan(sqrt(1.-cp2/(cq3))/(cp/sqrt(cq3)));
4871 if (co<0.)
4872 co = co + 3.1415927;
4873 v1 = (-2.)*sqrt(cq)*cos(co/3.) - cb/3.;
4874 v2 = (-2.)*sqrt(cq)*cos((co+2.*3.1415927)/3.) - cb/3.;
4875 v3 = (-2.)*sqrt(cq)*cos((co-2.*3.1415927)/3.) - cb/3.;
4876 }
4877
4878 return 0;
4879 }
4880
4881
4882
4883 /// calculates excess state functions of the fluid mixture
ResidualFunct()4884 long int TCORKcalc::ResidualFunct()
4885 {
4886 long int i, j;
4887 double sumphi, vi, vj, dvi, dvj, d2vi, d2vj, dvip, dvjp, Y, dY, d2Y, dYp,
4888 rhoi, drti, drtti, drpi, rhoj, drtj, drttj, drpj;
4889 double gex, dgex, d2gex, dgexp, sex, hex, cpex, vex, gph, sph, hph, cpph, vph;
4890 gex = 0.; dgex = 0.; d2gex = 0.; dgexp = 0.;
4891 sex = 0.; hex = 0.; vex = 0.; cpex = 0.;
4892 gph = 0.; sph = 0.; hph = 0.; cpph = 0.; vph = 0.;
4893 sumphi = 0.; Y = 0.; dY = 0.; d2Y = 0.; dYp = 0.;
4894
4895 // phi values (and derivatives)
4896 for (i=0; i<NComp; i++)
4897 {
4898 rhoi = Rho[i][0];
4899 drti = Rho[i][1];
4900 drtti = Rho[i][2];
4901 drpi = Rho[i][3];
4902 vi = 1./rhoi;
4903 dvi = - pow(rhoi,-2.)*drti;
4904 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
4905 dvip = - pow(rhoi,-2.)*drpi;
4906 sumphi += x[i]*vi;
4907 Y += x[i]*vi;
4908 dY += x[i]*dvi;
4909 d2Y += x[i]*d2vi;
4910 dYp += x[i]*dvip;
4911 }
4912
4913 for (i=0; i<NComp; i++)
4914 {
4915 rhoi = Rho[i][0];
4916 drti = Rho[i][1];
4917 drtti = Rho[i][2];
4918 drpi = Rho[i][3];
4919 vi = 1./rhoi;
4920 dvi = - pow(rhoi,-2.)*drti;
4921 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
4922 dvip = - pow(rhoi,-2.)*drpi;
4923 phi[i] = x[i]*vi/sumphi;
4924 dphi[i] = x[i]*(dvi*Y-vi*dY)/pow(Y,2.);
4925 dphip[i] = x[i]*(dvip*Y-vi*dYp)/pow(Y,2.);
4926 d2phi[i] = x[i]*(d2vi*Y+dvi*dY)/pow(Y,2.) - x[i]*(dvi*Y)*(2.*dY)/pow(Y,3.)
4927 - x[i]*(dvi*dY+vi*d2Y)/pow(Y,2.) + x[i]*(vi*dY)*(2.*dY)/pow(Y,3.);
4928 }
4929
4930 // interaction parameters
4931 for (i=0; i<NComp; i++)
4932 {
4933 for (j=i+1; j<NComp; j++)
4934 {
4935 rhoi = Rho[i][0];
4936 drti = Rho[i][1];
4937 drtti = Rho[i][2];
4938 drpi = Rho[i][3];
4939 rhoj = Rho[j][0];
4940 drtj = Rho[j][1];
4941 drttj = Rho[j][2];
4942 drpj = Rho[j][3];
4943
4944 vi = 1./rhoi;
4945 dvi = - pow(rhoi,-2.)*drti;
4946 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
4947 dvip = - pow(rhoi,-2.)*drpi;
4948 vj = 1./rhoj;
4949 dvj = - pow(rhoj,-2.)*drtj;
4950 d2vj = 2.*pow(rhoj,-3.)*pow(drtj,2.) - pow(rhoj,-2.)*drttj;
4951 dvjp = - pow(rhoj,-2.)*drpj;
4952
4953 B[i][j] = A[i][j]*2*Y/(vi*vj);
4954 dB[i][j] = A[i][j]*(2.*dY*vi*vj-2.*Y*(dvi*vj+vi*dvj))/(pow(vi,2.)*pow(vj,2.));
4955 d2B[i][j] = A[i][j]*((2.*d2Y)*(vi*vj) + (2.*dY)*(dvi*vj+vi*dvj))/pow((vi*vj),2.)
4956 - A[i][j]*((2.*dY)*(vi*vj))*(2.*(dvi*vj+vi*dvj))/pow((vi*vj),3.)
4957 - A[i][j]*((2.*dY)*(dvi*vj+vi*dvj) + (2.*Y)*(d2vi*vj+2.*dvi*dvj+vi*d2vj))/pow((vi*vj),2.)
4958 + A[i][j]*((2.*Y)*(dvi*vj+vi*dvj))*(2.*(dvi*vj+vi*dvj))/pow((vi*vj),3.);
4959 dBp[i][j] = A[i][j]*(2.*dYp*vi*vj-2.*Y*(dvip*vj+vi*dvjp))/(pow(vi,2.)*pow(vj,2.));
4960 }
4961 }
4962
4963 // excess properties
4964 for (i=0; i<NComp; i++)
4965 {
4966 for (j=i+1; j<NComp; j++)
4967 {
4968 gex += phi[i]*phi[j]*B[i][j];
4969 dgex += dphi[i]*phi[j]*B[i][j] + phi[i]*dphi[j]*B[i][j] + phi[i]*phi[j]*dB[i][j];
4970 d2gex +=d2phi[i]*phi[j]*B[i][j] + 2.*dphi[i]*dphi[j]*B[i][j] + 2.*dphi[i]*phi[j]*dB[i][j]
4971 + phi[i]*d2phi[j]*B[i][j] + 2.*phi[i]*dphi[j]*dB[i][j] + phi[i]*phi[j]*d2B[i][j];
4972 dgexp += dphip[i]*phi[j]*B[i][j] + phi[i]*dphip[j]*B[i][j] + phi[i]*phi[j]*dBp[i][j];
4973 }
4974 }
4975
4976 // thermodynamic properties
4977 sex = - dgex;
4978 hex = gex + Tk*sex;
4979 cpex = - d2gex*Tk;
4980 vex = dgexp;
4981 Gex = gex*(10.);
4982 Sex = sex*(10.);
4983 Hex = hex*(10.);
4984 CPex = cpex*(10.);
4985 Vex = vex;
4986 for (j=0; j<NComp; j++)
4987 {
4988 gph += x[j]*Fugpure[j][1];
4989 hph += x[j]*Fugpure[j][2];
4990 sph += x[j]*Fugpure[j][3];
4991 vph += x[j]*Fugpure[j][4];
4992 cpph += x[j]*Fugpure[j][5];
4993 }
4994 Grs = gph + Gex;
4995 Hrs = hph + Hex;
4996 Srs = sph + Sex;
4997 Vrs = vph + Vex;
4998 CPrs = cpph + CPex;
4999
5000 return 0;
5001 }
5002
5003
5004
5005 //#ifndef IPMGEMPLUGIN_
5006
5007 /// Calculates properties of pure fluids when called from DCthermo
CORKCalcFugPure(double Tmin,double * Cpg,double * FugProps)5008 long int TCORKcalc::CORKCalcFugPure( double Tmin, double *Cpg, double *FugProps )
5009 {
5010 long int retCode = 0;
5011 double Coeff[7];
5012
5013 for( int ii=0; ii<7; ii++ )
5014 Coeff[ii] = (double)Cpg[ii];
5015
5016 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
5017 {
5018 retCode = FugacityPT( 0, Coeff );
5019 for( int i=0; i<6; i++ )
5020 FugProps[i] = Fugpure[0][i];
5021 return retCode;
5022 }
5023
5024 else
5025 {
5026 for( int i=1; i<6; i++ )
5027 FugProps[i] = 0.;
5028 FugProps[0] = 1.;
5029 FugProps[4] = 8.31451*Tk/Pbar;
5030 return -1;
5031 }
5032 }
5033
5034 //#endif
5035
5036
5037
5038
5039
5040 //=======================================================================================================
5041 // Sterner-Pitzer (STP) model for fluid mixtures
5042 // References: Sterner and Pitzer (1994)
5043 // (c) TW December 2010
5044 //=======================================================================================================
5045
5046
5047 // Constructor
TSTPcalc(long int NCmp,double Pp,double Tkp,char Eos_Code)5048 TSTPcalc::TSTPcalc( long int NCmp, double Pp, double Tkp, char Eos_Code ):
5049 TSolMod( NCmp, '6', Tkp, Pp )
5050 {
5051
5052 UpdateTauP();
5053 Pparc = 0;
5054 alloc_internal();
5055 set_internal();
5056 EosCode[0] = Eos_Code;
5057
5058 // provisional
5059 if ( EosCode[0] == CEM_H2O_ ) // H2O
5060 Mw[0] = 18.015268;
5061 if ( EosCode[0] == CEM_CO2_ ) // CO2
5062 Mw[0] = 44.0098;
5063 }
5064
5065
5066
TSTPcalc(SolutionData * sd)5067 TSTPcalc::TSTPcalc( SolutionData *sd ):
5068 TSolMod( sd )
5069 {
5070 UpdateTauP();
5071 Pparc = aPparc;
5072 alloc_internal();
5073 set_internal();
5074 for( long int j=0; j<NComp; j++ )
5075 EosCode[j] = sd->TP_Code[j][3];
5076
5077 // provisional
5078 for ( long int j=0; j<NComp; j++ )
5079 {
5080 if ( EosCode[j] == CEM_H2O_ ) // H2O
5081 Mw[j] = 18.015268;
5082 if ( EosCode[j] == CEM_CO2_ ) // CO2
5083 Mw[j] = 44.0098;
5084 }
5085
5086 }
5087
5088
5089
~TSTPcalc()5090 TSTPcalc::~TSTPcalc()
5091 {
5092 free_internal();
5093 }
5094
5095
5096
5097 /// allocate work arrays for pure fluid and fluid mixture properties
alloc_internal()5098 void TSTPcalc::alloc_internal()
5099 {
5100 long int k;
5101
5102 EosCode = new char[NComp];
5103 Tc = new double [NComp];
5104 Pc = new double [NComp];
5105 Psat = new double [NComp];
5106 Rhol = new double [NComp];
5107 Rhov = new double [NComp];
5108 Mw = new double [NComp];
5109 Phi = new double [NComp];
5110 dPhiD = new double [NComp];
5111 dPhiDD = new double [NComp];
5112 dPhiT = new double [NComp];
5113 dPhiTT = new double [NComp];
5114 dPhiDT = new double [NComp];
5115 dPhiDDD = new double [NComp];
5116 dPhiDDT = new double [NComp];
5117 dPhiDTT = new double [NComp];
5118 phi = new double [NComp];
5119 dphi = new double [NComp];
5120 d2phi = new double [NComp];
5121 dphip = new double [NComp];
5122 lng = new double [NComp];
5123 Fugpure = new double [NComp][7];
5124 Rho = new double [NComp][11];
5125 cfh = new double *[10];
5126 cfc = new double *[10];
5127 A = new double *[NComp];
5128 W = new double *[NComp];
5129 B = new double *[NComp];
5130 dB = new double *[NComp];
5131 d2B = new double *[NComp];
5132 dBp = new double *[NComp];
5133
5134 for (k=0; k<10; k++)
5135 {
5136 cfh[k] = new double [6];
5137 cfc[k] = new double [6];
5138 }
5139
5140 for (k=0; k<NComp; k++)
5141 {
5142 A[k] = new double [NComp];
5143 W[k] = new double [NComp];
5144 B[k] = new double [NComp];
5145 dB[k] = new double [NComp];
5146 d2B[k] = new double [NComp];
5147 dBp[k] = new double [NComp];
5148 }
5149
5150 }
5151
5152
5153
free_internal()5154 void TSTPcalc::free_internal()
5155 {
5156 long int k;
5157
5158 for (k=0; k<10; k++)
5159 {
5160 delete[]cfh[k];
5161 delete[]cfc[k];
5162 }
5163
5164 for (k=0; k<NComp; k++)
5165 {
5166 delete[]A[k];
5167 delete[]W[k];
5168 delete[]B[k];
5169 delete[]dB[k];
5170 delete[]d2B[k];
5171 delete[]dBp[k];
5172 }
5173
5174 delete[]EosCode;
5175 delete[]Tc;
5176 delete[]Pc;
5177 delete[]Psat;
5178 delete[]Rhol;
5179 delete[]Rhov;
5180 delete[]Mw;
5181 delete[]Phi;
5182 delete[]dPhiD;
5183 delete[]dPhiDD;
5184 delete[]dPhiT;
5185 delete[]dPhiTT;
5186 delete[]dPhiDT;
5187 delete[]dPhiDDD;
5188 delete[]dPhiDDT;
5189 delete[]dPhiDTT;
5190 delete[]phi;
5191 delete[]dphi;
5192 delete[]d2phi;
5193 delete[]dphip;
5194 delete[]lng;
5195 delete[]Fugpure;
5196 delete[]Rho;
5197 delete[]cfh;
5198 delete[]cfc;
5199 delete[]A;
5200 delete[]W;
5201 delete[]B;
5202 delete[]dB;
5203 delete[]d2B;
5204 delete[]dBp;
5205 }
5206
5207
5208
5209 /// set model specific parameters
set_internal()5210 void TSTPcalc::set_internal()
5211 {
5212 long int j, k;
5213
5214 RC = R_CONST; // gas constant (bar)
5215 RR = 0.00831451; // gas constant (kbar)
5216 TMIN = 200.;
5217 TMAX = 2000.;
5218 PMIN = 1.0e-20;
5219 PMAX = 10000.;
5220
5221 // EoS coefficients (temperature dependence)
5222 double cfh_[10][6] = { { 0.0, 0.0, 0.24657688e6, 0.51359951e2, 0.0, 0.0 },
5223 { 0.0, 0.0, 0.58638965, -0.28646939e-2, 0.31375577e-4, 0.0 },
5224 { 0.0, 0.0, -0.62783840e1, 0.14791599e-1, 0.35779579e-3, 0.15432925e-7 },
5225 { 0.0, 0.0, 0.0, -0.42719875, -0.16325155e-4, 0.0 },
5226 { 0.0, 0.0, 0.56654978e4, -0.16580167e2, 0.76560762e-1, 0.0 },
5227 { 0.0, 0.0, 0.0, 0.10917883, 0.0, 0.0 },
5228 { 0.38878656e13, -0.13494878e9, 0.30916564e6, 0.75591105e1, 0.0, 0.0 },
5229 { 0.0, 0.0, -0.65537898e5, 0.18810675e3, 0.0, 0.0 },
5230 { -0.14182435e14, 0.18165390e9, -0.19769068e6, -0.23530318e2, 0.0, 0.0 },
5231 { 0.0, 0.0, 0.92093375e5, 0.12246777e3, 0.0, 0.0 } };
5232
5233
5234 double cfc_[10][6] = { { 0.0, 0.0, 0.18261340e7, 0.79224365e2, 0.0, 0.0 },
5235 { 0.0, 0.0, 0.0, 0.66560660e-4, 0.57152798e-5, 0.30222363e-9 },
5236 { 0.0, 0.0, 0.0, 0.59957845e-2, 0.71669631e-4, 0.62416103e-8 },
5237 { 0.0, 0.0, -0.13270279e1, -0.15210731, 0.53654244e-3, -0.71115142e-7 },
5238 { 0.0, 0.0, 0.12456776, 0.49045367e1, 0.98220560e-2, 0.55962121e-5 },
5239 { 0.0, 0.0, 0.0, 0.75522299, 0.0, 0.0 },
5240 { -0.39344644e12, 0.90918237e8, 0.42776716e6, -0.22347856e2, 0.0, 0.0 },
5241 { 0.0, 0.0, 0.40282608e3, 0.11971627e3, 0.0, 0.0 },
5242 { 0.0, 0.22995650e8, -0.78971817e5, -0.63376456e2, 0.0, 0.0 },
5243 { 0.0, 0.0, 0.95029765e5, 0.18038071e2, 0.0, 0.0 } };
5244
5245 // transfer coefficients
5246 for (j=0; j<10; j++)
5247 {
5248 for (k=0; k<6; k++)
5249 {
5250 cfh[j][k] = cfh_[j][k];
5251 cfc[j][k] = cfc_[j][k];
5252 }
5253 }
5254 }
5255
5256
5257
UpdateTauP()5258 long int TSTPcalc::UpdateTauP()
5259 {
5260 Pkbar = Pbar/1000.;
5261 Pkb = Pbar/1000.;
5262 Pmpa = Pbar/10.;
5263
5264 return 0;
5265 }
5266
5267
5268
5269 /// Calculates pure species properties (pure fugacities)
PureSpecies()5270 long int TSTPcalc::PureSpecies()
5271 {
5272 long int j, iErr = 0;
5273
5274 // error conditions
5275 // 1: input temperature/pressure outside valid range
5276 // 2: density iteration not converged
5277 // 3: wrong EoS subroutine code
5278
5279 for( j=0; j<NComp; j++ )
5280 {
5281 // Calling STP EoS for pure fugacity
5282 iErr = FugacityPT( j, aDCc+j*NP_DC );
5283 aGEX[j] = log( Fugpure[j][0] );
5284 Pparc[j] = Fugpure[j][0]*Pbar; // fure fluid fugacity (required for performance)
5285 aVol[j] = Fugpure[j][4]*10.; // molar volume of pure fluid component, J/bar to cm3
5286 } // j
5287
5288 if ( iErr )
5289 {
5290 char buf[150];
5291 sprintf(buf, "STP fluid: calculation of pure fluid fugacity failed");
5292 Error( "E71IPM IPMgamma: ", buf );
5293 }
5294
5295 return 0;
5296 }
5297
5298
5299
5300 /// Calculates T,P corrected interaction parameters
PTparam()5301 long int TSTPcalc::PTparam()
5302 {
5303 long int j, i, ip, i1, i2;
5304 double a;
5305
5306 UpdateTauP();
5307
5308 PureSpecies();
5309
5310 // set all interaction parameters zero
5311 for( j=0; j<NComp; j++ )
5312 {
5313 for( i=0; i<NComp; i++ )
5314 {
5315 A[j][i] = 0.;
5316 W[j][i] = 0.;
5317 B[j][i] = 0.;
5318 dB[j][i] = 0.;
5319 d2B[j][i] = 0.;
5320 dBp[j][i] = 0.;
5321 }
5322 }
5323
5324 // transfer interaction parameters that have non-standard value
5325 if( NPcoef > 0 )
5326 {
5327 for ( ip=0; ip<NPar; ip++ )
5328 {
5329 i1 = aIPx[MaxOrd*ip];
5330 i2 = aIPx[MaxOrd*ip+1];
5331 a = aIPc[NPcoef*ip];
5332 A[i1][i2] = a;
5333 A[i2][i1] = a; // symmetric case
5334 }
5335 }
5336
5337 return 0;
5338 }
5339
5340
5341
5342 /// Calculates activity coefficients
MixMod()5343 long int TSTPcalc::MixMod()
5344 {
5345 long int i, j, k;
5346 double dj, dk;
5347 double sumphi, lnGam, Gam, vi, vj, vk;
5348
5349 // phi values
5350 sumphi = 0.;
5351 for (i=0; i<NComp; i++)
5352 {
5353 vi = Fugpure[i][4];
5354 sumphi = sumphi + x[i]*vi;
5355 }
5356 for (i=0; i<NComp; i++)
5357 {
5358 phi[i] = x[i]*Fugpure[i][4]/sumphi;
5359 }
5360
5361 // interaction parameters
5362 for (i=0; i<NComp; i++)
5363 {
5364 for (j=i+1; j<NComp; j++)
5365 {
5366 vi = Fugpure[i][4];
5367 vj = Fugpure[j][4];
5368 W[i][j] = A[i][j]*(vi+vj)/(vi*vj);
5369 }
5370 }
5371
5372 // activity coefficients
5373 for (i=0; i<NComp; i++)
5374 {
5375 lnGam = 0.;
5376 for (j=0; j<NComp; j++)
5377 {
5378 for (k=j+1; k<NComp; k++)
5379 {
5380 vi = Fugpure[i][4];
5381 vj = Fugpure[j][4];
5382 vk = Fugpure[k][4];
5383 if (i==j)
5384 dj = 1.;
5385 else
5386 dj = 0.;
5387 if (i==k)
5388 dk = 1.;
5389 else
5390 dk = 0.;
5391 lnGam = lnGam - (dj-phi[j])*(dk-phi[k])*W[j][k]*2.*vi/(vj+vk);
5392 }
5393 }
5394 Gam = exp(lnGam/(RC*Tk));
5395 lnGamma[i] = log(Gam);
5396 } // i
5397
5398 return 0;
5399 }
5400
5401
5402
5403 /// calculates excess properties
ExcessProp(double * Zex)5404 long int TSTPcalc::ExcessProp( double *Zex )
5405 {
5406 long int iErr;
5407
5408 iErr = ResidualFunct();
5409
5410 if ( iErr )
5411 {
5412 char buf[150];
5413 sprintf(buf, "STP fluid: calculation failed");
5414 Error( "E71IPM IPMgamma: ", buf );
5415 }
5416
5417 Ars = Grs - Vrs*Pbar;
5418 Urs = Hrs - Vrs*Pbar;
5419
5420 // assignments (residual functions)
5421 Zex[0] = Grs;
5422 Zex[1] = Hrs;
5423 Zex[2] = Srs;
5424 Zex[3] = CPrs;
5425 Zex[4] = Vrs;
5426 Zex[5] = Ars;
5427 Zex[6] = Urs;
5428
5429 return iErr;
5430 }
5431
5432
5433
5434 /// calculates ideal mixing properties
IdealProp(double * Zid)5435 long int TSTPcalc::IdealProp( double *Zid )
5436 {
5437 long int j;
5438 double s, sc, sp;
5439
5440 s = 0.0;
5441 for ( j=0; j<NComp; j++ )
5442 {
5443 if ( x[j] > 1.0e-32 )
5444 s += x[j]*log(x[j]);
5445 }
5446 sc = (-1.)*R_CONST*s;
5447 sp = (-1.)*R_CONST*log(Pbar);
5448 Hid = 0.0;
5449 CPid = 0.0;
5450 Vid = 0.0;
5451 Sid = sc + sp;
5452 Gid = Hid - Sid*Tk;
5453 Aid = Gid - Vid*Pbar;
5454 Uid = Hid - Vid*Pbar;
5455
5456 // assignments (ideal mixing properties)
5457 Zid[0] = Gid;
5458 Zid[1] = Hid;
5459 Zid[2] = Sid;
5460 Zid[3] = CPid;
5461 Zid[4] = Vid;
5462 Zid[5] = Aid;
5463 Zid[6] = Uid;
5464
5465 return 0;
5466 }
5467
5468
5469
5470 /// high-level method to retrieve pure fluid properties
FugacityPT(long int j,double * EoSparam)5471 long int TSTPcalc::FugacityPT( long int j, double *EoSparam )
5472 {
5473 long int iErr = 0;
5474
5475 // reads EoS parameters from database into work array
5476 if( !EoSparam )
5477 return -1; // Memory alloc error
5478 Tc[j] = EoSparam[0]; // critical temperature (K)
5479 Pc[j] = EoSparam[1]/(10.); // critical pressure (MPa)
5480
5481 // check for temperature/pressure limits
5482 if ( (Tk<TMIN) || (Tk>TMAX) || (Pmpa<PMIN) || (Pmpa>PMAX) )
5483 {
5484 iErr = 1;
5485 return iErr;
5486 }
5487
5488 // select subroutines for different fluid types
5489 switch ( EosCode[j] )
5490 {
5491 case CEM_H2O_: // H2O
5492 iErr = FugacityH2O( j );
5493 break;
5494 case CEM_CO2_: // CO2
5495 iErr = FugacityCO2( j );
5496 break;
5497 case CEM_GAS_: // other fluids
5498 case CEM_CH4_:
5499 case CEM_N2_:
5500 case CEM_H2_:
5501 case CEM_O2_:
5502 case CEM_AR_:
5503 case CEM_PO_:
5504 case CEM_NP_:
5505 iErr = FugacityCorresponding( j );
5506 break;
5507 default:
5508 iErr = 3;
5509 break;
5510 }
5511
5512 return iErr;
5513 }
5514
5515
5516
5517 /// calculates fugacity and state functions of H2O
FugacityH2O(long int j)5518 long int TSTPcalc::FugacityH2O( long int j )
5519 {
5520 long int k, maxit, iErr = 0;
5521 double pmpa, tol, rhoguess, rho, pnew, pgrad, rhoold, rhonew, errnew, /*errold,*/ step,
5522 rhomin, rhomax, vol, lnfug, fug, fugc;
5523 double /*ar,*/ dard, dardd, darddd, dart, dartt, dardt, darddt, dardtt, pig,
5524 g, h, s, cp, cv, /*u, a,*/ dpr, dprr, dpt, dptt, dprt, drt, drtt, drp, drpp, drtp;
5525
5526 pmpa = Pbar/10.; // MPa
5527 tol = 1.e-10;
5528 maxit = 1000;
5529 rhomin = (1.e-20/1000.)/Mw[j]; // 1.e-20 kg m-3
5530 rhomax = (1800./1000.)/Mw[j]; // 1800 kg m-3
5531
5532 // find initial density guess
5533 DensityGuess( j, rhoguess );
5534 Pressure( rhoguess, pnew, pgrad, cfh );
5535
5536 rhoold = rhoguess;
5537 rhonew = rhoguess;
5538 errnew = pnew - pmpa;
5539 //errold = pnew - pmpa;
5540
5541 // find density by Newton iteration
5542 k = 0;
5543 do
5544 {
5545 k++;
5546
5547 if ( k>= maxit )
5548 {
5549 iErr = 2;
5550 return iErr;
5551 }
5552
5553 rhoold = rhonew;
5554 //errold = errnew;
5555 Pressure( rhonew, pnew, pgrad, cfh );
5556 errnew = pnew - pmpa;
5557 step = - errnew/pgrad;
5558 rhonew = rhoold + step;
5559
5560 if ( rhonew < rhomin )
5561 rhonew = rhomin ;
5562 if ( rhonew > rhomax )
5563 rhonew = rhomax;
5564
5565 } while ( fabs(1.-pnew/pmpa) > tol );
5566
5567 // calculate thermodynamic properties
5568 rho = rhonew;
5569 Helmholtz( j, rho, cfh );
5570 //ar = RC*Tk*Phi[j];
5571 dard = RC*Tk*dPhiD[j];
5572 dardd = RC*Tk*dPhiDD[j];
5573 darddd = RC*Tk*dPhiDDD[j];
5574 dart = RC*(dPhiT[j]*Tk+Phi[j]);
5575 dartt = RC*(dPhiTT[j]*Tk+2.*dPhiT[j]);
5576 dardt = RC*(dPhiDT[j]*Tk+dPhiD[j]);
5577 darddt = RC*(dPhiDDT[j]*Tk+dPhiDD[j]);
5578 dardtt = RC*(dPhiDTT[j]*Tk+2.*dPhiDT[j]);
5579
5580 vol = (1./rho)/10.; // J bar-1
5581 lnfug = ( log(rho) + Phi[j] + pmpa/(rho*RC*Tk) ) + log(RC*Tk) - 1.;
5582 fug = exp(lnfug);
5583 fugc = fug/pmpa;
5584 pig = rho*RC*Tk;
5585 g = log(fugc)*RC*Tk;
5586 s = - dart + RC*log(pmpa/pig);
5587 h = g + Tk*s;
5588 //a = ar - RC*Tk*log(pmpa/pig);
5589 cv = - Tk*dartt - RC;
5590 cp = cv + RC*pow((1.+(rho/RC)*dardt),2.) / (1.+(2.*rho)/(RC*Tk)*dard + pow(rho,2.)/(RC*Tk)*dardd);
5591 dpr = RC*Tk + 2.*rho*dard + pow(rho,2.)*dardd;
5592 dprr = 2.*dard + 4.*rho*dardd + pow(rho,2.)*darddd;
5593 dpt = rho*RC + pow(rho,2.)*dardt;
5594 dptt = pow(rho,2.)*dardtt;
5595 dprt = RC + 2.*rho*dardt + pow(rho,2.)*darddt;
5596 drp = (1./dpr);
5597 drpp = (-1.)*dprr*pow(dpr,-3.);
5598 drt = (-1.)*(1./dpr)*dpt;
5599 drtt = (dpt - dptt + dpt*(dprt-dpr)*pow(dpr,-1.))*pow(dpr,-1.)
5600 + (dprt - dpt*dprr*pow(dpr,-1.))*pow(dpr,-2.)*dpt;
5601 drtp = - pow(dpr,-1.) - (dprt-dpr)*pow(dpr,-2.) + dprr*dpt*pow(dpr,-3.);
5602
5603 // copy results
5604 Fugpure[j][0] = fugc; // fugacity coef.
5605 Fugpure[j][1] = g; // Gres
5606 Fugpure[j][2] = h; // Hres
5607 Fugpure[j][3] = s; // Sres
5608 Fugpure[j][4] = vol; // V
5609 Fugpure[j][5] = cp; // CPres
5610 Fugpure[j][6] = cv; // CVres
5611 Rho[j][0] = rho; // mol cm-3
5612 Rho[j][1] = drt; // mol cm-3 K-1
5613 Rho[j][2] = drtt; // mol cm-3 K-2
5614 Rho[j][3] = drp; // mol cm-3 MPa-1
5615 Rho[j][4] = drpp; // mol cm-3 MPa-2
5616 Rho[j][5] = drtp; // mol cm-3 K-1 MPa-1
5617 Rho[j][6] = dpr; // MPa cm3 mol-1
5618 Rho[j][7] = dprr; // MPa cm6 mol-2
5619 Rho[j][8] = dpt; // MPa K-1
5620 Rho[j][9] = dptt; // MPa K-2
5621 Rho[j][10] = dprt; // MPa cm3 mol-1 K-1
5622
5623 return iErr;
5624 }
5625
5626
5627
5628 /// calculates fugacity and state functions of CO2
FugacityCO2(long int j)5629 long int TSTPcalc::FugacityCO2( long int j )
5630 {
5631 long int k, maxit, iErr = 0;
5632 double pmpa, tol, rhoguess, rho, pnew, pgrad, rhoold, rhonew, errnew, /*errold,*/ step,
5633 rhomin, rhomax, vol, lnfug, fug, fugc;
5634 double /*ar,*/ dard, dardd, darddd, dart, dartt, dardt, darddt, dardtt, pig,
5635 g, h, s, cp, cv, /*u, a,*/ dpr, dprr, dpt, dptt, dprt, drt, drtt, drp, drpp, drtp;
5636
5637 pmpa = Pbar/10.; // MPa
5638 tol = 1.e-10;
5639 maxit = 1000;
5640 rhomin = (1.e-20)/1000./Mw[j]; // 1.e-20 kg m-3
5641 rhomax = (2400.)/1000./Mw[j]; // 2400 kg m-3
5642
5643 // find initial density guess
5644 DensityGuess( j, rhoguess );
5645 Pressure( rhoguess, pnew, pgrad, cfc );
5646
5647 rhoold = rhoguess;
5648 rhonew = rhoguess;
5649 errnew = pnew - pmpa;
5650 //errold = errnew;
5651
5652 k = 0;
5653 do
5654 {
5655 k++;
5656
5657 if ( k>= maxit )
5658 {
5659 iErr = 2;
5660 return iErr;
5661 }
5662
5663 rhoold = rhonew;
5664 //errold = errnew;
5665 Pressure( rhonew, pnew, pgrad, cfc );
5666 errnew = pnew - pmpa;
5667 step = - errnew/pgrad;
5668 rhonew = rhoold + step;
5669
5670 if ( rhonew < rhomin )
5671 rhonew = rhomin ;
5672 if ( rhonew > rhomax )
5673 rhonew = rhomax;
5674
5675 } while ( fabs(1.-pnew/pmpa) > tol );
5676
5677 // calculate thermodynamic properties
5678 rho = rhonew;
5679 Helmholtz( j, rho, cfc );
5680 //ar = RC*Tk*Phi[j];
5681 dard = RC*Tk*dPhiD[j];
5682 dardd = RC*Tk*dPhiDD[j];
5683 darddd = RC*Tk*dPhiDDD[j];
5684 dart = RC*(dPhiT[j]*Tk+Phi[j]);
5685 dartt = RC*(dPhiTT[j]*Tk+2.*dPhiT[j]);
5686 dardt = RC*(dPhiDT[j]*Tk+dPhiD[j]);
5687 darddt = RC*(dPhiDDT[j]*Tk+dPhiDD[j]);
5688 dardtt = RC*(dPhiDTT[j]*Tk+2.*dPhiDT[j]);
5689
5690 vol = (1./rho)/10.; // J bar-1
5691 lnfug = ( log(rho) + Phi[j] + pmpa/(rho*RC*Tk) ) + log(RC*Tk) - 1.;
5692 fug = exp(lnfug);
5693 fugc = fug/pmpa;
5694 pig = rho*RC*Tk;
5695 g = log(fugc)*RC*Tk;
5696 s = - dart + RC*log(pmpa/pig);
5697 h = g + Tk*s;
5698 //a = ar - RC*Tk*log(pmpa/pig);
5699 cv = - Tk*dartt - RC;
5700 cp = cv + RC*pow((1.+(rho/RC)*dardt),2.) / (1.+(2.*rho)/(RC*Tk)*dard + pow(rho,2.)/(RC*Tk)*dardd);
5701 dpr = RC*Tk + 2.*rho*dard + pow(rho,2.)*dardd;
5702 dprr = 2.*dard + 4.*rho*dardd + pow(rho,2.)*darddd;
5703 dpt = rho*RC + pow(rho,2.)*dardt;
5704 dptt = pow(rho,2.)*dardtt;
5705 dprt = RC + 2.*rho*dardt + pow(rho,2.)*darddt;
5706 drp = (1./dpr);
5707 drpp = (-1.)*dprr*pow(dpr,-3.);
5708 drt = (-1.)*(1./dpr)*dpt;
5709 drtt = (dpt - dptt + dpt*(dprt-dpr)*pow(dpr,-1.))*pow(dpr,-1.)
5710 + (dprt - dpt*dprr*pow(dpr,-1.))*pow(dpr,-2.)*dpt;
5711 drtp = - pow(dpr,-1.) - (dprt-dpr)*pow(dpr,-2.) + dprr*dpt*pow(dpr,-3.);
5712
5713 // copy results
5714 Fugpure[j][0] = fugc; // fugacity coef.
5715 Fugpure[j][1] = g; // Gres
5716 Fugpure[j][2] = h; // Hres
5717 Fugpure[j][3] = s; // Sres
5718 Fugpure[j][4] = vol; // V
5719 Fugpure[j][5] = cp; // CPres
5720 Fugpure[j][6] = cv; // CVres
5721 Rho[j][0] = rho; // mol cm-3
5722 Rho[j][1] = drt; // mol cm-3 K-1
5723 Rho[j][2] = drtt; // mol cm-3 K-2
5724 Rho[j][3] = drp; // mol cm-3 MPa-1
5725 Rho[j][4] = drpp; // mol cm-3 MPa-2
5726 Rho[j][5] = drtp; // mol cm-3 K-1 MPa-1
5727 Rho[j][6] = dpr; // MPa cm3 mol-1
5728 Rho[j][7] = dprr; // MPa cm6 mol-2
5729 Rho[j][8] = dpt; // MPa K-1
5730 Rho[j][9] = dptt; // MPa K-2
5731 Rho[j][10] = dprt; // MPa cm3 mol-1 K-1
5732
5733 return iErr;
5734 }
5735
5736
5737
5738 /// calculates fugacity and state functions of fluids other than H2O and CO2
5739 /// adapted from CORK fluid model for consistency with Thermocalc
FugacityCorresponding(long int j)5740 long int TSTPcalc::FugacityCorresponding( long int j )
5741 {
5742 double a0, a1, a, b0, b, c0, c1, c, d0, d1, d, tcr, pcr, da, dc, dd, vc, fc, rtlnf, rho;
5743 double grs, hrs, srs, cprs, drt, drtt, drp, drpp, drtp,
5744 dvt, dvtt, dvp, dvpp, dvtp;
5745
5746 a0 = 5.45963e-5;
5747 a1 = -8.6392e-6;
5748 b0 = 9.18301e-4;
5749 c0 = -3.30558e-5;
5750 c1 = 2.30524e-6;
5751 d0 = 6.93054e-7;
5752 d1 = -8.38293e-8;
5753 tcr = Tc[j];
5754 // pcr = Eosparm[j][1]/1000.; // kbar
5755 pcr = Pc[j]*(10.)/(1000.); // kbar (convert from MPa)
5756 a = a0*pow(tcr,2.)*sqrt(tcr)/pcr + a1*tcr*sqrt(tcr)/pcr*Tk;
5757 b = b0*tcr/pcr;
5758 c = c0*tcr/(pcr*sqrt(pcr)) + c1/(pcr*sqrt(pcr))*Tk;
5759 d = d0*tcr/pow(pcr,2.) + d1/pow(pcr,2.)*Tk;
5760 da = a1*pow(tcr,1.5)/pcr;
5761 dc = c1/pow(pcr,1.5);
5762 dd = d1/pow(pcr,2.);
5763
5764 // molar volume and fugacity coefficient
5765 vc = RR*Tk/Pkb + b - a*RR*sqrt(Tk)/((RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb)) + c*sqrt(Pkb) + d*Pkb;
5766 rtlnf = RR*Tk*log(1000*Pkb) + b*Pkb + a/(b*sqrt(Tk))*(log(RR*Tk+b*Pkb) - log(RR*Tk+2*b*Pkb))
5767 + (2./3.)*c*Pkb*sqrt(Pkb) + (d/2.)*pow(Pkb,2.);
5768 fc = exp(rtlnf/(RR*Tk))/(1000*Pkb);
5769
5770 // volume and density derivatives
5771 dvt = ( RR/Pkb - (da*RR*pow(Tk,0.5))/((RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
5772 - (a*RR)/(2.*pow(Tk,0.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
5773 + (a*pow(RR,2.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5774 + (a*pow(RR,2.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5775 + dc*pow(Pkb,0.5) + dd*Pkb );
5776 dvtt = ( - (da*RR)/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
5777 + (2.*da*pow(RR,2.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5778 + (2.*da*pow(RR,2.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5779 + (a*RR)/(4.*pow(Tk,1.5)*(RR*Tk+b*Pkb)*(RR*Tk+2.*b*Pkb))
5780 + (a*pow(RR,2.))/(pow(Tk,0.5)*pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5781 + (a*pow(RR,2.))/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5782 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
5783 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
5784 - (2.*a*pow(RR,3.)*pow(Tk,0.5))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.)) );
5785 dvp = ( - RR*Tk/pow(Pkb,2.) + (a*RR*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5786 + (2.*a*RR*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5787 + c/(2.*pow(Pkb,0.5)) + d ) / (1000.);
5788 dvpp = ( (2.*RR*Tk)/pow(Pkb,3.) - (2.*a*RR*pow(Tk,0.5)*pow(b,2.))/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
5789 - (4.*a*RR*pow(Tk,0.5)*pow(b,2.))/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
5790 - (8.*a*RR*pow(Tk,0.5)*pow(b,2.))/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.))
5791 - c/(4.*pow(Pkb,1.5)) ) / (1.0e6);
5792 dvtp = ( - RR/pow(Pkb,2.) + (da*RR*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5793 + (a*RR*b)/(2.*pow(Tk,0.5)*pow((RR*Tk+b*Pkb),2.)*(RR*Tk+2.*b*Pkb))
5794 - (2.*a*pow(RR,2.)*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),3.)*(RR*Tk+2.*b*Pkb))
5795 - (3.*a*pow(RR,2.)*pow(Tk,0.5)*b)/(pow((RR*Tk+b*Pkb),2.)*pow((RR*Tk+2.*b*Pkb),2.))
5796 + (2.*da*RR*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5797 + (a*RR*b)/(pow(Tk,0.5)*(RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),2.))
5798 - (4.*a*pow(RR,2.)*pow(Tk,0.5)*b)/((RR*Tk+b*Pkb)*pow((RR*Tk+2.*b*Pkb),3.))
5799 + dc/(2.*pow(Pkb,2.)) + dd ) / (1000.);
5800 rho = (1./vc);
5801 drt = - dvt*pow(rho,2.);
5802 drp = - dvp*pow(rho,2.);
5803 drtt = - (dvtt*pow(rho,3.)-2.*pow(drt,2.))/rho;
5804 drpp = - (dvpp*pow(rho,3.)-2.*pow(drp,2.))/rho;
5805 drtp = - (dvtp*pow(rho,3.)-2.*drt*drp)/rho;
5806
5807 // thermodynamic properties
5808 grs = R_CONST*Tk*log(fc);
5809 hrs = ( (pow(Tk,0.5)*da*log(RR*Tk+2.*b*Pkb))/b - (pow(Tk,0.5)*da*log(RR*Tk+b*Pkb))/b
5810 - (3.*a*log(RR*Tk+2.*b*Pkb))/(2.*pow(Tk,0.5)*b) + (3.*a*log(RR*Tk+b*Pkb))/(2.*pow(Tk,0.5)*b)
5811 - (pow(Tk,0.5)*a*RR)/(b*(RR*Tk+b*Pkb)) + (pow(Tk,0.5)*a*RR)/(b*(RR*Tk+2.*b*Pkb))
5812 - (2.*Tk*dc*pow(Pkb,1.5))/3. - (Tk*dd*pow(Pkb,2.))/2. + b*Pkb
5813 + (2.*c*pow(Pkb,1.5))/3. + (d*pow(Pkb,2.))/2. )*1000.;
5814 srs = (hrs-grs)/Tk;
5815 cprs = ( - (2.*pow(Tk,0.5)*da*RR)/(b*(RR*Tk+b*Pkb)) + (pow(Tk,0.5)*a*pow(RR,2.))/(b*pow((RR*Tk+b*Pkb),2.))
5816 + (2.*pow(Tk,0.5)*da*RR)/(b*(RR*Tk+2.*b*Pkb)) - (a*RR)/(pow(Tk,0.5)*b*(RR*Tk+2.*b*Pkb))
5817 - (da*log(RR*Tk+2.*b*Pkb))/(pow(Tk,0.5)*b) + (da*log(RR*Tk+b*Pkb))/(pow(Tk,0.5)*b)
5818 + (3.*a*log(RR*Tk+2.*b*Pkb))/(4.*pow(Tk,1.5)*b) - (3.*a*log(RR*Tk+b*Pkb))/(4.*pow(Tk,1.5)*b)
5819 - (pow(Tk,0.5)*a*pow(RR,2.))/(b*pow((RR*Tk+2.*b*Pkb),2.)) + (a*RR)/(pow(Tk,0.5)*b*(RR*Tk+b*Pkb)) )*1000.;
5820
5821 // copy results
5822 Fugpure[j][0] = fc; // fugacity coefficient
5823 Fugpure[j][1] = grs;
5824 Fugpure[j][2] = hrs;
5825 Fugpure[j][3] = srs;
5826 Fugpure[j][4] = vc;
5827 Fugpure[j][5] = cprs;
5828 Rho[j][0] = rho*(0.1); // mol cm-3
5829 Rho[j][1] = drt*(0.1); // mol cm-3 K-1
5830 Rho[j][2] = drtt*(0.1); // mol cm-3 K-2
5831 Rho[j][3] = drp; // mol cm-3 MPa-1
5832 Rho[j][4] = drpp*(10.); // mol cm-3 MPa-2
5833 Rho[j][5] = drtp; // mol cm-3 K-1 MPa-1
5834 Rho[j][6] = 1.0; // MPa cm3 mol-1
5835 Rho[j][7] = 1.0; // MPa cm6 mol-2
5836 Rho[j][8] = 1.0; // MPa K-1
5837 Rho[j][9] = 1.0; // MPa K-2
5838 Rho[j][10] = 1.0; // MPa cm3 mol-1 K-1
5839
5840 return 0;
5841 }
5842
5843
5844
5845 /// calculates and returns density guess for pure fluids
DensityGuess(long int j,double & Rhoguess)5846 long int TSTPcalc::DensityGuess( long int j, double &Rhoguess )
5847 {
5848 double tred, pred, rhomin, rhomax, rhoguess;
5849
5850 tred = Tk/Tc[j];
5851 pred = Pmpa/Pc[j];
5852
5853 if ( EosCode[j] == CEM_H2O_ ) // H2O
5854 {
5855 rhomin = (1.e-20/1000.)/Mw[j]; // 1.e-20 kg m-3
5856 rhomax = (1800./1000.)/Mw[j]; // 1800 kg m-3
5857
5858 if ( Tk < Tc[j] )
5859 {
5860 PsatH2O( j );
5861 if ( Pmpa < Psat[j] )
5862 rhoguess = Rhov[j];
5863 else
5864 rhoguess = Rhol[j];
5865 }
5866
5867 else
5868 {
5869 if ( Pmpa < 400. )
5870 rhoguess = 1000./(Tk-273.15)*2.*Pmpa;
5871 else
5872 rhoguess = 800.;
5873 }
5874
5875 Rhoguess = (rhoguess/1000.)/Mw[j]; // mol cm-3
5876
5877 if ( Rhoguess < rhomin )
5878 Rhoguess = rhomin;
5879 if ( Rhoguess > rhomax )
5880 Rhoguess = rhomax;
5881 }
5882
5883 else // CO2 (other fluids)
5884 {
5885 rhomin = (1.e-20)/1000./Mw[j]; // 1.e-20 kg m-3
5886 rhomax = (2400.)/1000./Mw[j]; // 2400 kg m-3
5887
5888 if ( Tk < Tc[j] )
5889 {
5890 PsatCO2( j );
5891 if ( Pmpa < Psat[j] )
5892 rhoguess = Rhov[j]; // kg m-3
5893 else
5894 rhoguess = Rhol[j];
5895 }
5896
5897 else
5898 {
5899 rhoguess = (pred*Pc[j]*1000.)/((0.1889241)*tred*Tc[j]); // kg m-3
5900 }
5901
5902 Rhoguess = (rhoguess/1000.)/Mw[j]; // mol cm-3
5903
5904 if ( Rhoguess < rhomin )
5905 Rhoguess = rhomin;
5906 if ( Rhoguess > rhomax )
5907 Rhoguess = rhomax;
5908 }
5909
5910 return 0;
5911 }
5912
5913
5914
5915 /// calculates pressure (P) and first density derivative (dP/dRho)
Pressure(double rho,double & p,double & dpdrho,double ** cf)5916 long int TSTPcalc::Pressure ( double rho, double &p, double &dpdrho, double **cf )
5917 {
5918 long int k;
5919 double pred, dpred;
5920 double c[10];
5921
5922 for (k=0; k<10; k++)
5923 {
5924 c[k] = cf[k][0]*pow(Tk,-4.) + cf[k][1]*pow(Tk,-2.) + cf[k][2]*pow(Tk,-1.) + cf[k][3]
5925 + cf[k][4]*Tk + cf[k][5]*pow(Tk,2.);
5926 }
5927
5928 pred = rho + c[0]*pow(rho,2.) - pow(rho,2.) * ( (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
5929 / pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.) )
5930 + c[6]*pow(rho,2.)*exp(-c[7]*rho) + c[8]*pow(rho,2.)*exp(-c[9]*rho);
5931 dpred = 1. + 2.*c[0]*rho - 2.*rho*(c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
5932 / pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5933 - pow(rho,2.)*(2.*c[3]+6.*c[4]*rho+12.*c[5]*pow(rho,2.))
5934 / pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5935 + 2.*pow(rho,2.)*pow((c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.)),2.)
5936 / pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
5937 + 2.*c[6]*rho*exp(-c[7]*rho) - c[6]*pow(rho,2.)*c[7]*exp(-c[7]*rho)
5938 + 2.*c[8]*rho*exp(-c[9]*rho) - c[8]*pow(rho,2.)*c[9]*exp(-c[9]*rho);
5939
5940 p = pred*(RC*Tk);
5941 dpdrho = dpred*(RC*Tk);
5942
5943 return 0;
5944 }
5945
5946
5947
5948 /// calculates reduced Helmholtz energy and derivatives
Helmholtz(long int j,double rho,double ** cf)5949 long int TSTPcalc::Helmholtz( long int j, double rho, double **cf )
5950 {
5951 long int k;
5952 double c[10], dc[10], d2c[10];
5953
5954 for (k=0; k<10; k++)
5955 {
5956 c[k] = cf[k][0]*pow(Tk,-4.) + cf[k][1]*pow(Tk,-2.) + cf[k][2]*pow(Tk,-1.) + cf[k][3]
5957 + cf[k][4]*Tk + cf[k][5]*pow(Tk,2.);
5958 dc[k] = - 4.*cf[k][0]*pow(Tk,-5.) - 2.*cf[k][1]*pow(Tk,-3.) - cf[k][2]*pow(Tk,-2.)
5959 + cf[k][4] + 2.*cf[k][5]*Tk;
5960 d2c[k] = 20.*cf[k][0]*pow(Tk,-6.) + 6.*cf[k][1]*pow(Tk,-4.) + 2.*cf[k][2]*pow(Tk,-3.)
5961 + 2.*cf[k][5];
5962 }
5963
5964 Phi[j] = c[0]*rho + ( 1./(c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)) - 1./c[1] )
5965 - (c[6]/c[7])*(exp(-c[7]*rho)-1.) - (c[8]/c[9])*(exp(-c[9]*rho)-1.);
5966 dPhiD[j] = c[0] - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5967 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
5968 + c[6]*exp(-c[7]*rho) + c[8]*exp(-c[9]*rho);
5969 dPhiDD[j] = 2./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
5970 * pow((c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.)),2.)
5971 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5972 * (2.*c[3]+6.*c[4]*rho+12.*c[5]*pow(rho,2.))
5973 - c[6]*c[7]*exp(-c[7]*rho) - c[8]*c[9]*exp(-c[9]*rho);
5974 dPhiT[j] = dc[0]*rho - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5975 * (dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.)) + 1./pow(c[1],2.)*dc[1]
5976 - dc[6]/c[7]*(exp(-c[7]*rho)-1.) + c[6]/pow(c[7],2.)*(exp(-c[7]*rho)-1.)*dc[7]
5977 + c[6]/c[7]*dc[7]*rho*exp(-c[7]*rho) - dc[8]/c[9]*(exp(-c[9]*rho)-1)
5978 + c[8]/pow(c[9],2.)*(exp(-c[9]*rho)-1.)*dc[9] + c[8]/c[9]*dc[9]*rho*exp(-c[9]*rho);
5979 dPhiTT[j] = - 2./pow(c[1],3.)*pow(dc[1],2.)
5980 + 2./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
5981 * pow((dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.)),2.)
5982 - d2c[8]/c[9]*(exp(-c[9]*rho)-1.) - d2c[6]/c[7]*(exp(-c[7]*rho)-1.)
5983 + 2.*dc[6]/pow(c[7],2.)*(exp(-c[7]*rho)-1)*dc[7] + d2c[0]*rho
5984 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5985 * (d2c[1]+d2c[2]*rho+d2c[3]*pow(rho,2.)+d2c[4]*pow(rho,3.)+d2c[5]*pow(rho,4.))
5986 - 2.*c[8]/pow(c[9],2.)*pow(dc[9],2.)*rho*exp(-c[9]*rho) + c[8]/pow(c[9],2.)*(exp(-c[9]*rho)-1.)*d2c[9]
5987 - 2.*c[8]/pow(c[9],3.)*(exp(-c[9]*rho)-1.)*pow(dc[9],2.) + 2.*dc[8]/c[9]*dc[9]*rho*exp(-c[9]*rho)
5988 + 2.*dc[8]/pow(c[9],2.)*(exp(-c[9]*rho)-1.)*dc[9] - c[6]/c[7]*pow(dc[7],2.)*pow(rho,2.)*exp(-c[7]*rho)
5989 + c[6]/c[7]*d2c[7]*rho*exp(-c[7]*rho) - 2.*c[6]/pow(c[7],2.)*pow(dc[7],2.)*rho*exp(-c[7]*rho)
5990 + c[6]/pow(c[7],2.)*(exp(-c[7]*rho)-1.)*d2c[7] - 2.*c[6]/pow(c[7],3.)*(exp(-c[7]*rho)-1.)*pow(dc[7],2.)
5991 + 2.*dc[6]/c[7]*dc[7]*rho*exp(-c[7]*rho) + 1./pow(c[1],2.)*d2c[1]
5992 - c[8]/c[9]*pow(dc[9],2.)*pow(rho,2.)*exp(-c[9]*rho) + c[8]/c[9]*d2c[9]*rho*exp(-c[9]*rho);
5993 dPhiDT[j] = dc[0] + 2./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
5994 * (dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.))
5995 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
5996 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
5997 * (dc[2]+2.*dc[3]*rho+3.*dc[4]*pow(rho,2.)+4.*dc[5]*pow(rho,3.))
5998 + dc[6]*exp(-c[7]*rho) - c[6]*dc[7]*rho*exp(-c[7]*rho)
5999 + dc[8]*exp(-c[9]*rho) - c[8]*dc[9]*rho*exp(-c[9]*rho);
6000 dPhiDDD[j] = - 6./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),4.)
6001 * pow((c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.)),3.)
6002 + 6./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
6003 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))*(2.*c[3]+6.*c[4]*rho+12.*c[5]*pow(rho,2.))
6004 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.) * (6.*c[4]+24.*c[5]*rho)
6005 + c[6]*pow(c[7],2.)*exp(-c[7]*rho)+c[8]*pow(c[9],2.)*exp(-c[9]*rho);
6006 dPhiDDT[j] = - 6./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),4.)
6007 * (dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.))
6008 * pow((c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.)),2.)
6009 + 4./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
6010 * (dc[2]+2.*dc[3]*rho+3.*dc[4]*pow(rho,2.)+4.*dc[5]*pow(rho,3.))
6011 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
6012 + 2./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
6013 * (dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.))
6014 * (2.*c[3]+6.*c[4]*rho+12.*c[5]*pow(rho,2.))
6015 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
6016 * (2.*dc[3]+6.*dc[4]*rho+12.*dc[5]*pow(rho,2.))
6017 - dc[6]*c[7]*exp(-c[7]*rho) - c[6]*dc[7]*exp(-c[7]*rho)
6018 + c[6]*dc[7]*rho*c[7]*exp(-c[7]*rho) - dc[8]*c[9]*exp(-c[9]*rho)
6019 - c[8]*dc[9]*exp(-c[9]*rho) + c[8]*dc[9]*rho*c[9]*exp(-c[9]*rho);
6020 dPhiDTT[j] = 2./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
6021 * (d2c[1]+d2c[2]*rho+d2c[3]*pow(rho,2.)+d2c[4]*pow(rho,3.)+d2c[5]*pow(rho,4.))
6022 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.)) + d2c[0]
6023 - 2.*dc[8]*dc[9]*rho*exp(-c[9]*rho) - c[6]*d2c[7]*rho*exp(-c[7]*rho)
6024 + c[6]*pow(dc[7],2.)*pow(rho,2.)*exp(-c[7]*rho)
6025 - 6./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),4.)
6026 * pow((dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.)),2.)
6027 * (c[2]+2.*c[3]*rho+3.*c[4]*pow(rho,2.)+4.*c[5]*pow(rho,3.))
6028 + 4./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),3.)
6029 * (dc[1]+dc[2]*rho+dc[3]*pow(rho,2.)+dc[4]*pow(rho,3.)+dc[5]*pow(rho,4.))
6030 * (dc[2]+2.*dc[3]*rho+3.*dc[4]*pow(rho,2.)+4.*dc[5]*pow(rho,3.)) + d2c[6]*exp(-c[7]*rho)
6031 - 1./pow((c[1]+c[2]*rho+c[3]*pow(rho,2.)+c[4]*pow(rho,3.)+c[5]*pow(rho,4.)),2.)
6032 * (d2c[2]+2.*d2c[3]*rho+3.*d2c[4]*pow(rho,2.)+4.*d2c[5]*pow(rho,3.))
6033 - 2.*dc[6]*dc[7]*rho*exp(-c[7]*rho) + c[8]*pow(dc[9],2.)*pow(rho,2.)*exp(-c[9]*rho)
6034 - c[8]*d2c[9]*rho*exp(-c[9]*rho) + d2c[8]*exp(-c[9]*rho);
6035
6036 return 0;
6037 }
6038
6039
6040
6041 /// calculates residual state functions of the fluid mixture
ResidualFunct()6042 long int TSTPcalc::ResidualFunct()
6043 {
6044 long int i, j;
6045 double sumphi, vi, vj, dvi, dvj, d2vi, d2vj, dvip, dvjp, Y, dY, d2Y, dYp,
6046 rhoi, drti, drtti, drpi, rhoj, drtj, drttj, drpj;
6047 double gex, dgex, d2gex, dgexp, sex, hex, cpex, vex, gph, sph, hph, cpph, vph;
6048 gex = 0.; dgex = 0.; d2gex = 0.; dgexp = 0.;
6049 sex = 0.; hex = 0.; vex = 0.; cpex = 0.;
6050 gph = 0.; sph = 0.; hph = 0.; cpph = 0.; vph = 0.;
6051 sumphi = 0.; Y = 0.; dY = 0.; d2Y = 0.; dYp = 0.;
6052
6053 // phi values (and derivatives)
6054 for (i=0; i<NComp; i++)
6055 {
6056 rhoi = Rho[i][0];
6057 drti = Rho[i][1];
6058 drtti = Rho[i][2];
6059 drpi = Rho[i][3];
6060 vi = 1./rhoi;
6061 dvi = - pow(rhoi,-2.)*drti;
6062 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
6063 dvip = - pow(rhoi,-2.)*drpi;
6064 sumphi += x[i]*vi;
6065 Y += x[i]*vi;
6066 dY += x[i]*dvi;
6067 d2Y += x[i]*d2vi;
6068 dYp += x[i]*dvip;
6069 }
6070
6071 for (i=0; i<NComp; i++)
6072 {
6073 rhoi = Rho[i][0];
6074 drti = Rho[i][1];
6075 drtti = Rho[i][2];
6076 drpi = Rho[i][3];
6077 vi = 1./rhoi;
6078 dvi = - pow(rhoi,-2.)*drti;
6079 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
6080 dvip = - pow(rhoi,-2.)*drpi;
6081 phi[i] = x[i]*vi/sumphi;
6082 dphi[i] = x[i]*(dvi*Y-vi*dY)/pow(Y,2.);
6083 dphip[i] = x[i]*(dvip*Y-vi*dYp)/pow(Y,2.);
6084 d2phi[i] = x[i]*(d2vi*Y+dvi*dY)/pow(Y,2.) - x[i]*(dvi*Y)*(2.*dY)/pow(Y,3.)
6085 - x[i]*(dvi*dY+vi*d2Y)/pow(Y,2.) + x[i]*(vi*dY)*(2.*dY)/pow(Y,3.);
6086 }
6087
6088 // interaction parameters
6089 for (i=0; i<NComp; i++)
6090 {
6091 for (j=i+1; j<NComp; j++)
6092 {
6093 rhoi = Rho[i][0];
6094 drti = Rho[i][1];
6095 drtti = Rho[i][2];
6096 drpi = Rho[i][3];
6097 rhoj = Rho[j][0];
6098 drtj = Rho[j][1];
6099 drttj = Rho[j][2];
6100 drpj = Rho[j][3];
6101
6102 vi = 1./rhoi;
6103 dvi = - pow(rhoi,-2.)*drti;
6104 d2vi = 2.*pow(rhoi,-3.)*pow(drti,2.) - pow(rhoi,-2.)*drtti;
6105 dvip = - pow(rhoi,-2.)*drpi;
6106 vj = 1./rhoj;
6107 dvj = - pow(rhoj,-2.)*drtj;
6108 d2vj = 2.*pow(rhoj,-3.)*pow(drtj,2.) - pow(rhoj,-2.)*drttj;
6109 dvjp = - pow(rhoj,-2.)*drpj;
6110
6111 B[i][j] = A[i][j]*2*Y/(vi*vj);
6112 dB[i][j] = A[i][j]*(2.*dY*vi*vj-2.*Y*(dvi*vj+vi*dvj))/(pow(vi,2.)*pow(vj,2.));
6113 d2B[i][j] = A[i][j]*((2.*d2Y)*(vi*vj) + (2.*dY)*(dvi*vj+vi*dvj))/pow((vi*vj),2.)
6114 - A[i][j]*((2.*dY)*(vi*vj))*(2.*(dvi*vj+vi*dvj))/pow((vi*vj),3.)
6115 - A[i][j]*((2.*dY)*(dvi*vj+vi*dvj) + (2.*Y)*(d2vi*vj+2.*dvi*dvj+vi*d2vj))/pow((vi*vj),2.)
6116 + A[i][j]*((2.*Y)*(dvi*vj+vi*dvj))*(2.*(dvi*vj+vi*dvj))/pow((vi*vj),3.);
6117 dBp[i][j] = A[i][j]*(2.*dYp*vi*vj-2.*Y*(dvip*vj+vi*dvjp))/(pow(vi,2.)*pow(vj,2.));
6118 }
6119 }
6120
6121 // excess properties
6122 for (i=0; i<NComp; i++)
6123 {
6124 for (j=i+1; j<NComp; j++)
6125 {
6126 gex += phi[i]*phi[j]*B[i][j];
6127 dgex += dphi[i]*phi[j]*B[i][j] + phi[i]*dphi[j]*B[i][j] + phi[i]*phi[j]*dB[i][j];
6128 d2gex +=d2phi[i]*phi[j]*B[i][j] + 2.*dphi[i]*dphi[j]*B[i][j] + 2.*dphi[i]*phi[j]*dB[i][j]
6129 + phi[i]*d2phi[j]*B[i][j] + 2.*phi[i]*dphi[j]*dB[i][j] + phi[i]*phi[j]*d2B[i][j];
6130 dgexp += dphip[i]*phi[j]*B[i][j] + phi[i]*dphip[j]*B[i][j] + phi[i]*phi[j]*dBp[i][j];
6131 }
6132 }
6133
6134 // thermodynamic properties
6135 sex = - dgex;
6136 hex = gex + Tk*sex;
6137 cpex = - d2gex*Tk;
6138 vex = dgexp;
6139 Gex = gex*(10.);
6140 Sex = sex*(10.);
6141 Hex = hex*(10.);
6142 CPex = cpex*(10.);
6143 Vex = vex;
6144
6145 // incrementing residual properties
6146 for (j=0; j<NComp; j++)
6147 {
6148 gph += x[j]*Fugpure[j][1];
6149 hph += x[j]*Fugpure[j][2];
6150 sph += x[j]*Fugpure[j][3];
6151 vph += x[j]*Fugpure[j][4];
6152 cpph += x[j]*Fugpure[j][5];
6153 }
6154 Grs = gph + Gex;
6155 Hrs = hph + Hex;
6156 Srs = sph + Sex;
6157 Vrs = vph + Vex;
6158 CPrs = cpph + CPex;
6159
6160 return 0;
6161 }
6162
6163
6164
6165 /// calculates saturation pressure of H2O (required only for initial density guess)
PsatH2O(long int j)6166 long int TSTPcalc::PsatH2O( long int j )
6167 {
6168 double Rhoc, a1, a2, a3, a4, a5, a6, b1, b2, b3, b4, b5, b6, c1, c2, c3, c4, c5, c6,
6169 tau, psat, rhol, rhov, ppc, rrc;
6170
6171 Rhoc = 322.0; // kg m-3
6172 a1 = -7.85951783;
6173 a2 = 1.84408259;
6174 a3 = -11.7866497;
6175 a4 = 22.6807411;
6176 a5 = -15.9618719;
6177 a6 = 1.80122502;
6178 b1 = 1.99274064;
6179 b2 = 1.09965342;
6180 b3 = -0.510839303;
6181 b4 = -1.75493479;
6182 b5 = -45.5170352;
6183 b6 = -6.74694450e5;
6184 c1 = -2.03150240;
6185 c2 = -2.68302940;
6186 c3 = -5.38626492;
6187 c4 = -17.2991605;
6188 c5 = -44.7586581;
6189 c6 = -63.9201063;
6190
6191 tau = 1 - Tk/Tc[j];
6192 ppc = (Tc[j]/Tk)*( a1*tau + a2*pow(tau,1.5) + a3*pow(tau,3.) + a4*pow(tau,3.5)
6193 + a5*pow(tau,4.) + a6*pow(tau,7.5) );
6194 psat = exp(ppc) * Pc[j];
6195 rhol = ( 1. + b1*pow(tau,(1./3.)) + b2*pow(tau,(2./3.)) + b3*pow(tau,(5./3.)) + b4*pow(tau,(16./3.))
6196 + b5*pow(tau,(43./3.)) + b6*pow(tau,(110./3.)) ) * Rhoc;
6197 rrc = c1*pow(tau,2./6.) + c2*pow(tau,4./6.) + c3*pow(tau,8./6.) + c4*pow(tau,18./6.)
6198 + c5*pow(tau,37./6.) + c6*pow(tau,71./6.);
6199 rhov = exp(rrc) * Rhoc;
6200
6201 // copy results
6202 Psat[j] = psat; // MPa
6203 Rhol[j] = rhol;
6204 Rhov[j] = rhov;
6205
6206 return 0;
6207 }
6208
6209
6210
6211 /// calculates saturation pressure of CO2 (required only for initial density guess)
PsatCO2(long int j)6212 long int TSTPcalc::PsatCO2( long int j )
6213 {
6214 double Rhoc, a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, c5, tau, ppc, psat,
6215 rhol, rhov, rholc, rhovc;
6216
6217 Rhoc = 467.6; // kg m-3
6218 a1 = -7.0602087;
6219 a2 = 1.9391218;
6220 a3 = -1.6463597;
6221 a4 = -3.2995634;
6222 b1 = 1.9245108;
6223 b2 = -0.62385555;
6224 b3 = -0.32731127;
6225 b4 = 0.39245142;
6226 c1 = -1.7074879;
6227 c2 = -0.82274670;
6228 c3 = -4.6008549;
6229 c4 = -10.111178;
6230 c5 = -29.742252;
6231
6232 tau = 1. - Tk/Tc[j];
6233 ppc = (Tc[j]/Tk)*( a1*pow(tau,1.) + a2*pow(tau,1.5) + a3*pow(tau,2.) + a4*pow(tau,4.) );
6234 psat = exp(ppc)*Pc[j];
6235 rholc = b1*pow(tau,0.34)+ b2*pow(tau,0.5) + b3*pow(tau,(10./6.)) + b4*pow(tau,(11./6.));
6236 rhol = exp(rholc)*Rhoc;
6237 rhovc = c1*pow(tau,0.34) + c2*pow(tau,0.5) + c3*pow(tau,1.) + c4*pow(tau,(7./3.)) + c5*pow(tau,(14./3.));
6238 rhov = exp(rhovc)*Rhoc;
6239
6240 // copy results
6241 Psat[j] = psat; // MPa
6242 Rhol[j] = rhol; // kg m-3
6243 Rhov[j] = rhov; // kg m-3
6244
6245 return 0;
6246 }
6247
6248
6249
6250 #ifndef IPMGEMPLUGIN_
6251
6252 /// Calculates pure species properties (called from DCthermo)
STPCalcFugPure(double Tmin,double * Cpg,double * FugProps)6253 long int TSTPcalc::STPCalcFugPure(double Tmin, double *Cpg, double *FugProps )
6254 {
6255 long int iErr = 0;
6256 double Coeff[7];
6257
6258 for( int ii=0; ii<7; ii++ )
6259 Coeff[ii] = (double)Cpg[ii];
6260
6261 if( (Tk >= Tmin) && (Tk < 1e4) && (Pbar >= 1e-5) && (Pbar < 1e5) )
6262 {
6263 iErr = FugacityPT( 0, Coeff );
6264 for( int i=0; i<6; i++ )
6265 FugProps[i] = Fugpure[0][i];
6266 return iErr;
6267 }
6268
6269 else
6270 {
6271 for( int i=1; i<6; i++ )
6272 FugProps[i] = 0.;
6273 FugProps[0] = 1.;
6274 FugProps[4] = 8.31451*Tk/Pbar;
6275 return -1;
6276 }
6277
6278 return 0;
6279 }
6280
6281
6282 #endif
6283
6284 }
6285
6286
6287
6288
6289
6290 //--------------------- End of s_solmod2.cpp ---------------------------
6291
6292
6293