1 //------------------------------------------------------------------- 2 // $Id: s_solmod.h 725 2012-10-02 15:43:37Z kulik $ 3 // 4 /// \file s_solmod.h 5 /// Declarations of TSolMod and derived classes implementing built-in models 6 /// of mixing in fluid, liquid, aqueous and solid-solution phases 7 8 // Copyright (C) 2003-2014 T.Wagner, D.Kulik, S.Dmitrieva, F.Hingerl, 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 29 #ifndef _s_solmod_h_ 30 #define _s_solmod_h_ 31 32 #include <vector> 33 #include <cstring> 34 #include <iostream> 35 using namespace std; 36 37 namespace solmod 38 { 39 40 // re-declaration of enums below required for GEMS4K 41 // dc_class_codes for fluids will be replaced by tp_codes 42 enum fluid_mix_rules { /// codes for mixing rules in EoS models (see m_phase.h) 43 MR_UNDEF_ = 'N', 44 MR_WAAL_ = 'W', 45 MR_CONST_ = 'C', 46 MR_TEMP_ = 'T', 47 MR_LJ_ = 'J', 48 MR_KW1_ = 'K', 49 MR_PITZ5_ = '5', 50 MR_PITZ6_ = '6', 51 MR_PITZ8_ = '8', 52 MR_B_RCPT_ = 'R' 53 }; 54 55 enum dc_class_codes { /// codes for fluid types in EoS models (see v_mod.h) 56 DC_GAS_H2O_ = 'V', 57 DC_GAS_CO2_ = 'C', 58 DC_GAS_H2_ = 'H', 59 DC_GAS_N2_ = 'N', 60 DC_GAS_COMP_ = 'G' 61 }; 62 63 enum tp_codes { /// codes for fluid subroutines in EoS models (see v_mod.h) 64 CEM_OFF_ = 'N', 65 CEM_GAS_ = 'G', 66 CEM_H2O_ = 'V', 67 CEM_CO2_ = 'C', 68 CEM_CH4_ = 'M', 69 CEM_N2_ = 'T', 70 CEM_H2_ = 'H', 71 CEM_O2_ = 'O', 72 CEM_AR_ = 'A', 73 CEM_PO_ = 'P', 74 CEM_NP_ = 'Q' 75 }; 76 77 78 // ------------------------------------------------------------------ 79 80 #define MAXPHASENAME 16 81 82 /// Base class for subclasses of built-in mixing models. 83 /// (c) March 2007 DK/TW 84 struct SolutionData { 85 long int NSpecies; ///< Number of species (end members) in the phase 86 long int NParams; ///< Total number of non-zero interaction parameters 87 long int NPcoefs; ///< Number of coefficients per interaction parameter 88 long int MaxOrder; ///< Maximum order of interaction parameters 89 long int NPperDC; ///< Number of parameters per species (DC) 90 long int NSublat; ///< number of sublattices nS 91 long int NMoiet; ///< number of moieties nM 92 93 // long int NlPhs; ///< new: Number of linked phases 94 // long int NlPhC; ///< new: Number of linked phase parameter coefficient per link (default 0) 95 long int NDQFpDC; ///< new: Number of DQF parameters per species (end member) 96 // long int NrcPpDC; ///< new: Number of reciprocal parameters per species (end member) 97 98 char Mod_Code; ///< Code of the mixing model 99 char Mix_Code; ///< Code for specific EoS mixing rule 100 char *DC_Codes; ///< DC class codes for species -> NSpecies 101 char (*TP_Code)[6]; ///< Codes for TP correction methods for species ->NSpecies 102 long int *arIPx; ///< Pointer to list of indexes of non-zero interaction parameters 103 104 // long int *arPhLin; ///< new: indexes of linked phase and link type codes [NlPhs*2] read-only 105 106 double *arIPc; ///< Table of interaction parameter coefficients 107 double *arDCc; ///< End-member properties coefficients 108 double *arMoiSN; ///< End member moiety- site multiplicity number tables -> NSpecies x NSublat x NMoiet 109 double *arSitFr; ///< Tables of sublattice site fractions for moieties -> NSublat x NMoiet 110 double *arSitFj; ///< new: Table of end member sublattice activity coefficients -> NSpecies x NSublat 111 double *arGEX; ///< Pure-species fugacities, G0 increment terms -> NSpecies 112 113 // double *lPhc; ///< new: array of phase link parameters -> NlPhs x NlPhC (read-only) 114 double *DQFc; ///< new: array of DQF parameters for DCs in phases -> NSpecies x NDQFpDC; (read-only) 115 // double *rcpc; ///< new: array of reciprocal parameters for DCs in phases -> NSpecies x NrcPpDC; (read-only) 116 117 double *arPparc; ///< Partial pressures -> NSpecies 118 double *arWx; ///< Species (end member) mole fractions ->NSpecies 119 double *arlnGam; ///< Output: activity coefficients of species (end members) 120 121 // Detailed output on terms of partial end-member properties, allocated in MULTI 122 double *arlnDQFt; ///< new: DQF terms adding to overall activity coefficients [Ls_] 123 double *arlnRcpt; ///< new: reciprocal terms adding to overall activity coefficients [Ls_] 124 double *arlnExet; ///< new: excess energy terms adding to overall activity coefficients [Ls_] 125 double *arlnCnft; ///< new: configurational terms adding to overall activity [Ls_] 126 127 double *arVol; ///< molar volumes of end-members (species) cm3/mol ->NSpecies 128 double *aphVOL; ///< phase volumes, cm3/mol (now obsolete) !!!!!!! check usage! 129 double T_k; ///< Temperature, K (initial) 130 double P_bar; ///< Pressure, bar (initial) 131 132 void set_def(); 133 }; 134 135 136 class TSolMod 137 { 138 protected: 139 char ModCode; ///< Code of the mixing model 140 char MixCode; ///< Code for specific EoS mixing rules 141 char *DC_Codes; ///< Class codes of end members (species) ->NComp 142 143 char PhaseName[MAXPHASENAME+1]; ///< Phase name (for specific built-in models) 144 145 long int NComp; ///< Number of components in the solution phase 146 long int NPar; ///< Number of non-zero interaction parameters 147 long int NPcoef; ///< Number of coeffs per parameter (columns in the aIPc table) 148 long int MaxOrd; ///< max. parameter order (or number of columns in aIPx) 149 long int NP_DC; ///< Number of coeffs per one DC in the phase (columns in aDCc) 150 long int NSub; ///< number of sublattices nS 151 long int NMoi; ///< number of moieties nM 152 153 // long int NlPh; ///< new: Number of linked phases 154 // long int NlPc; ///< new: Number of linked phase parameter coefficient per link (default 0) 155 long int NDQFpc; ///< new: Number of DQF parameters per species (end member), 0 or 4 156 // long int NrcPpc; ///< new: Number of reciprocal parameters per species (end member) 157 158 // long int NPTP_DC; // Number of properties per one DC at T,P of interest (columns in aDC) !!!! Move to CG EOS subclass 159 long int *aIPx; // Pointer to list of indexes of non-zero interaction parameters 160 // long int (*PhLin)[2]; ///< new: indexes of linked phase and link type codes [NlPhs][2] read-only 161 162 double R_CONST; ///< R constant 163 double Tk; ///< Temperature, K 164 double Pbar; ///< Pressure, bar 165 166 double *aIPc; ///< Table of interaction parameter coefficients 167 double *aIP; ///< Vector of interaction parameters corrected to T,P of interest 168 double *aDCc; ///< End-member properties coefficients 169 double *aGEX; ///< Reciprocal energies, DQF terms, pure fugacities of DC (corrected to TP) 170 double *aPparc; ///< Output partial pressures (activities, fugacities) -> NComp 171 double **aDC; ///< Table of corrected end member properties at T,P of interest !!!!!! Move to GC EOS subclass! 172 double *aMoiSN; ///< End member moiety- site multiplicity number tables -> NComp x NSub x NMoi 173 double *aSitFR; ///< Table of sublattice site fractions for moieties -> NSub x NMoi 174 175 // double *lPhcf; ///< new: array of phase link parameters -> NlPh x NlPc (read-only) 176 double *DQFcf; ///< new: array of DQF parameters for DCs in phases -> NComp x NDQFpc; (read-only) 177 ///< x_DQF[j]: mole fraction at transition; a, b, c - coefficients of T,P correction 178 ///< according to the equation aGEX[j] = A + B*T + C*P (so far only binary Margules) 179 // double *rcpcf; ///< new: array of reciprocal parameters for DCs in phases -> NComp x NrcPpc; (read-only) 180 181 double *x; ///< Pointer to mole fractions of end members (provided) 182 double *aVol; ///< molar volumes of species (end members) 183 double *phVOL; ///< phase volume, cm3/mol (now obsolete) !!!!!!!!!!!! Check usage! 184 185 // Results 186 // double Gam; ///< work cell for activity coefficient of end member 187 // double lnGamRT; 188 // double lnGam; 189 double Gex, Hex, Sex, CPex, Vex, Aex, Uex; ///< molar excess properties of the phase 190 double Gid, Hid, Sid, CPid, Vid, Aid, Uid; ///< molar ideal mixing properties 191 double Gdq, Hdq, Sdq, CPdq, Vdq, Adq, Udq; ///< molar Darken quadratic terms 192 double Grs, Hrs, Srs, CPrs, Vrs, Ars, Urs; ///< molar residual functions (fluids) 193 double *lnGamConf, *lnGamRecip, *lnGamEx, *lnGamDQF; ///< Work pointers for lnGamma components 194 double *lnGamma; ///< Pointer to ln activity coefficients of end members (check that it is collected from three above arrays) 195 196 double **y; ///< table of moiety site fractions [NSub][NMoi] 197 double ***mn; ///< array of end member moiety-site multiplicity numbers [NComp][NSub][NMoi] 198 double *mns; ///< array of total site multiplicities [NSub] 199 double **fjs; ///< array of site activity coefficients [NComp][NSub] 200 double *aSitFj; ///< new: pointer to return table of site activity coefficients NComp x NSub 201 202 // functions for calculation of configurational term for multisite ideal mixing 203 void alloc_multisite(); 204 long int init_multisite(); 205 void free_multisite(); 206 void free_sdata(); 207 208 209 /// Functions for calculation of configurational term for multisite ideal mixing 210 long int IdealMixing(); 211 double ideal_conf_entropy(); 212 void return_sitefr(); 213 void retrieve_sitefr(); 214 215 216 public: 217 218 /// Generic constructor 219 TSolMod( SolutionData *sd ); 220 221 /// Generic constructor for DComp/DCthermo 222 TSolMod( long int NSpecies, char Mod_Code, double T_k, double P_bar ); 223 224 /// Destructor 225 virtual ~TSolMod(); 226 PureSpecies()227 virtual long int PureSpecies() 228 { 229 return 0; 230 } 231 PTparam()232 virtual long int PTparam() 233 { 234 return 0; 235 } 236 MixMod()237 virtual long int MixMod() 238 { 239 return 0; 240 } 241 ExcessProp(double *)242 virtual long int ExcessProp( double */*Zex*/ ) 243 { 244 return 0; 245 } 246 IdealProp(double *)247 virtual long int IdealProp( double */*Zid*/ ) 248 { 249 return 0; 250 } 251 Set_Felect_bc(long int,double,double)252 virtual long int Set_Felect_bc (long int /*Flagelect*/, double /*Bc*/, double /*Ac*/) 253 { 254 return 0; 255 } 256 257 258 /// Set new system state 259 long int UpdatePT ( double T_k, double P_bar ); 260 261 bool testSizes( SolutionData *sd ); 262 263 /// Getting phase name 264 void GetPhaseName( const char *PhName ); 265 266 267 // copy activity coefficients into provided array lngamma Get_lnGamma(double * lngamma)268 inline void Get_lnGamma( double* lngamma ) 269 { 270 for( int i=0; i<NComp; i++ ) 271 lngamma[i] = lnGamma[i]; 272 } 273 274 void getSolutionData( SolutionData *sd ); 275 276 // access from node 277 void Set_aIPc( const vector<double> aIPc_ ); 278 void Get_aIPc ( vector<double> &aIPc_ ); 279 280 void Set_aDCc( const vector<double> aDCc_ ); 281 void Get_aDCc( vector<double> &aDCc_ ); 282 283 void Get_aIPx( vector<long int> &aIPx_ ); 284 285 void Get_NPar_NPcoef_MaxOrd_NComp_NP_DC ( long int &NPar_, long int &NPcoef_, 286 long int &MaxOrd_, long int &NComp_, long int &NP_DC_ ); 287 288 }; 289 290 291 /// Subclass for the ideal model (both simple and multi-site) 292 class TIdeal: public TSolMod 293 { 294 private: 295 296 public: 297 298 /// Constructor 299 TIdeal( SolutionData *sd ); 300 301 /// Destructor 302 ~TIdeal(); 303 304 /// Calculates T,P corrected interaction parameters 305 long int PTparam(); 306 307 /// Calculates (fictive) activity coefficients 308 long int MixMod(); 309 310 /// Calculates excess properties 311 long int ExcessProp( double *Zex ); 312 313 /// Calculates ideal mixing properties 314 long int IdealProp( double *Zid ); 315 316 }; 317 318 319 320 /// Churakov & Gottschalk (2003) EOS calculations 321 /// declaration of EOSPARAM class (used by the TCGFcalc class) 322 class EOSPARAM 323 { 324 private: 325 //static double Told; 326 // unsigned long int isize; // int isize = NComp; 327 long int NComp; 328 double emix, s3mix; 329 double *epspar,*sig3par; 330 double *XX; 331 double *eps; 332 double *eps05; 333 double *sigpar; 334 double *mpar; 335 double *apar; 336 double *aredpar; 337 double *m2par; 338 double **mixpar; 339 340 void allocate(); 341 void free(); 342 343 public: 344 345 double *XX0; 346 347 //EOSPARAM():isize(0),emix(0),s3mix(0),NComp(0){}; 348 //EOSPARAM(double*data, unsigned nn):isize(0){allocate(nn);init(data,nn);}; 349 EOSPARAM(double * Xtmp,double * data,long int nn)350 EOSPARAM( double *Xtmp, double *data, long int nn ) 351 :NComp(nn), emix(0),s3mix(0) 352 { allocate(); init(Xtmp,data,nn); } 353 ~EOSPARAM()354 ~EOSPARAM() 355 { free(); } 356 357 void init( double*,double *, long int ); NCmp()358 long int NCmp() {return NComp; } 359 EPS05(long int i)360 double EPS05( long int i) 361 { return eps05[i]; } X(long int i)362 double X( long int i) 363 { return XX[i]; } EPS(long int i)364 double EPS( long int i) 365 { return eps[i]; } EMIX(void)366 double EMIX(void) 367 { return emix; } S3MIX(void)368 double S3MIX(void) 369 { return s3mix; } 370 MIXS3(long int i,long int j)371 double MIXS3( long int i, long int j) 372 { 373 if (i==j) return sig3par[i]; 374 if (i<j) return mixpar[i][j]; 375 else return mixpar[j][i]; 376 } 377 MIXES3(long int i,long int j)378 double MIXES3( long int i, long int j) 379 { 380 if ( i==j ) return epspar[i]; 381 if (i<j) return mixpar[j][i]; 382 else return mixpar[i][j]; 383 } 384 SIG3(long int i)385 double SIG3( long int i){ return sig3par[i]; } M2R(long int i)386 double M2R( long int i) { return m2par[i]; } A(long int i)387 double A( long int i) { return apar[i]; } 388 389 long int ParamMix( double *Xin); 390 }; 391 392 393 394 // ------------------------------------------------------------------------------------- 395 /// Churakov and Gottschalk (2003) EOS calculations. 396 /// Added 09 May 2003 397 /// Declaration of a class for CG EOS calculations for fluids 398 /// Incorporates a C++ program written by Sergey Churakov (CSCS ETHZ) 399 /// implementing papers by Churakov and Gottschalk (2003a, 2003b) 400 class TCGFcalc: public TSolMod 401 { 402 private: 403 404 double 405 PI_1, ///< pi 406 TWOPI, ///< 2.*pi 407 PISIX, ///< pi/6. 408 TWOPOW1SIX, ///< 2^(1/6) 409 DELTA, 410 DELTAMOLLIM, 411 R, NA, P1, 412 PP2, P3, P4, 413 P5, P6, P7, 414 P8, P9, P10, 415 AA1, AA2, AA3, 416 A4, A5, A6, 417 BB1, BB2, BB3, 418 B4, B5, B6, 419 A00, A01, A10, 420 A11, A12, A21, 421 A22, A23, A31, 422 A32, A33, A34; 423 424 // double PhVol; // phase volume in cm3 425 double *Pparc; ///< DC partial pressures (pure fugacities) 426 double *phWGT; 427 double *aX; ///< DC quantities at eqstate x_j (moles) 428 // double *aGEX; // Increments to molar G0 values of DCs from pure fugacities 429 // double *aVol; // DC molar volumes, cm3/mol [L] 430 431 // main work arrays 432 EOSPARAM *paar; 433 EOSPARAM *paar1; 434 double *FugCoefs; 435 double *EoSparam; 436 double *EoSparam1; 437 double (*Cf)[8]; ///< corrected EoS coefficients 438 439 // internal functions 440 void alloc_internal(); 441 void free_internal(); 442 void set_internal(); 443 444 void choose( double *pres, double P,unsigned long int &x1,unsigned long int &x2 ); 445 double Melt2( double T ); 446 double Melt( double T ); 447 void copy( double* sours,double *dest,unsigned long int num ); 448 void norm( double *X,unsigned long int mNum ); 449 double RPA( double beta,double nuw ); 450 double dHS( double beta,double ro ); 451 fI1_6(double nuw)452 inline double fI1_6( double nuw ) 453 { 454 return (1.+(A4+(A5+A6*nuw)*nuw)*nuw)/ 455 ((1.+(AA1+(AA2+AA3*nuw)*nuw)*nuw)*3.); 456 } 457 fI1_12(double nuw)458 inline double fI1_12( double nuw ) 459 { 460 return (1.+(B4+(B5+B6*nuw)*nuw)*nuw)/ 461 ((1.+(BB1+(BB2+BB3*nuw)*nuw)*nuw)*9.); 462 } 463 fa0(double nuw,double nu1w2)464 inline double fa0( double nuw ,double nu1w2 ) 465 { 466 return (A00 + A01*nuw)/nu1w2; 467 } 468 fa1(double nuw,double nu1w3)469 inline double fa1( double nuw ,double nu1w3 ) 470 { 471 return (A10+(A11+A12*nuw)*nuw)/nu1w3; 472 } 473 fa2(double nuw,double nu1w4)474 inline double fa2( double nuw ,double nu1w4 ) 475 { 476 return ((A21+(A22+A23*nuw)*nuw)*nuw)/nu1w4; 477 } 478 fa3(double nuw,double nu1w5)479 inline double fa3( double nuw ,double nu1w5 ) 480 { 481 return ((A31+(A32+(A33+A34*nuw)*nuw)*nuw)*nuw)/nu1w5; 482 } 483 484 double DIntegral( double T, double ro, unsigned long int IType ); // not used 485 double LIntegral( double T, double ro, unsigned long int IType ); // not used 486 double KIntegral( double T, double ro, unsigned long int IType ); // not used 487 double K23_13( double T, double ro ); 488 double J6LJ( double T,double ro ); 489 double FDipPair( double T,double ro,double m2 ); // not used 490 double UWCANum( double T,double ro ); 491 double ZWCANum( double T,double ro ); 492 493 double FWCA( double T,double ro ); 494 double FTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param ); 495 double UTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param ); // not used 496 double ZTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param ); 497 double PTOTALMIX( double T_Real,double ro_Real,EOSPARAM* param ); 498 double ROTOTALMIX( double P,double TT,EOSPARAM* param ); 499 500 double PRESSURE( double *X, double *param, unsigned long int NN, double ro, double T ); // not used 501 double DENSITY( double *X,double *param, unsigned long int NN ,double Pbar, double T ); 502 long int CGActivCoefRhoT( double *X,double *param, double *act, unsigned long int NN, 503 double ro, double T ); // not used 504 505 long int CGActivCoefPT(double *X,double *param,double *act, unsigned long int NN, 506 double Pbar, double T, double &roro ); 507 508 public: 509 510 /// Constructor 511 TCGFcalc( long int NCmp, double Pp, double Tkp ); 512 TCGFcalc( SolutionData *sd, double *aphWGT, double *arX ); 513 514 /// Destructor 515 ~TCGFcalc(); 516 517 /// Calculates of pure species properties (pure fugacities) 518 long int PureSpecies( ); 519 520 /// Calculates T,P corrected interaction parameters 521 long int PTparam(); 522 523 /// Calculates activity coefficients 524 long int MixMod(); 525 526 /// Calculates excess properties 527 long int ExcessProp( double *Zex ); 528 529 ///< calculates ideal mixing properties 530 long int IdealProp( double *Zid ); 531 532 /// CGofPureGases, calculates fugacity for 1 species at (X=1) 533 long int CGcalcFugPure( double Tmin, double *Cemp, double *FugProps ); // called from DCthermo 534 long int CGFugacityPT( double *EoSparam, double *EoSparPT, double &Fugacity, 535 double &Volume, double P, double T, double &roro ); 536 537 /// Calculates departure functions 538 long int CGResidualFunct( double *X, double *param, double *param1, unsigned long int NN, 539 double ro, double T ); 540 GetDELTA(void)541 double GetDELTA( void ) 542 { 543 return DELTA; 544 } 545 }; 546 547 548 549 // ------------------------------------------------------------------------------------- 550 /// Peng-Robinson-Stryjek-Vera (PRSV) model for fluid mixtures. 551 /// References: Stryjek and Vera (1986) 552 /// (c) TW July 2006 553 class TPRSVcalc: public TSolMod 554 555 { 556 private: 557 558 double PhVol; ///< phase volume in cm3 559 double *Pparc; ///< DC partial pressures (pure fugacities) 560 // double *aGEX; // Increments to molar G0 values of DCs from pure fugacities 561 // double *aVol; // DC molar volumes, cm3/mol [L] 562 563 // main work arrays 564 double (*Eosparm)[6]; ///< EoS parameters 565 double (*Pureparm)[4]; ///< Parameters a, b, da/dT, d2a/dT2 for cubic EoS 566 double (*Fugpure)[6]; ///< fugacity parameters of pure gas species 567 double (*Fugci)[4]; ///< fugacity parameters of species in the mixture 568 569 double **a; ///< arrays of generic parameters 570 double **b; 571 double **KK; ///< binary interaction parameter 572 double **dKK; ///< derivative of interaction parameter 573 double **d2KK; ///< second derivative 574 double **AA; ///< binary a terms in the mixture 575 576 577 578 // internal functions 579 void alloc_internal(); 580 void free_internal(); 581 long int AB( double Tcrit, double Pcrit, double omg, double k1, double k2, double k3, 582 double &apure, double &bpure, double &da, double &d2a ); 583 long int FugacityPT( long int i, double *EoSparam ); 584 long int FugacityPure( long int j ); // Calculates the fugacity of pure species 585 long int Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 ); 586 long int MixParam( double &amix, double &bmix ); 587 long int FugacityMix( double amix, double bmix, double &fugmix, double &zmix, double &vmix ); 588 long int FugacitySpec( double *fugpure ); 589 long int ResidualFunct( double *fugpure ); 590 long int MixingWaals(); 591 long int MixingConst(); 592 long int MixingTemp(); 593 594 public: 595 596 /// Constructor 597 TPRSVcalc( long int NCmp, double Pp, double Tkp ); 598 TPRSVcalc( SolutionData *sd ); 599 600 /// Destructor 601 ~TPRSVcalc(); 602 603 /// Calculates pure species properties (pure fugacities) 604 long int PureSpecies(); 605 606 /// Calculates T,P corrected interaction parameters 607 long int PTparam(); 608 609 /// Calculates activity coefficients 610 long int MixMod(); 611 612 /// Calculates excess properties 613 long int ExcessProp( double *Zex ); 614 615 /// Calculates ideal mixing properties 616 long int IdealProp( double *Zid ); 617 618 /// Calculates pure species properties (called from DCthermo) 619 long int PRSVCalcFugPure( double Tmin, double *Cpg, double *FugProps ); 620 621 }; 622 623 624 625 // ------------------------------------------------------------------------------------- 626 /// Soave-Redlich-Kwong (SRK) model for fluid mixtures. 627 /// References: Soave (1972); Soave (1993) 628 /// (c) TW December 2008 629 class TSRKcalc: public TSolMod 630 631 { 632 private: 633 634 double PhVol; ///< phase volume in cm3 635 double *Pparc; ///< DC partial pressures (pure fugacities) 636 // double *aGEX; // Increments to molar G0 values of DCs from pure fugacities 637 // double *aVol; // DC molar volumes, cm3/mol [L] 638 639 // main work arrays 640 double (*Eosparm)[4]; ///< EoS parameters 641 double (*Pureparm)[4]; ///< Parameters a, b, da/dT, d2a/dT2 for cubic EoS 642 double (*Fugpure)[6]; ///< Fugacity parameters of pure gas species 643 double (*Fugci)[4]; ///< Fugacity parameters of species in the mixture 644 645 double **a; ///< arrays of generic parameters 646 double **b; 647 double **KK; ///< binary interaction parameter 648 double **dKK; ///< derivative of interaction parameter 649 double **d2KK; ///< second derivative 650 double **AA; ///< binary a terms in the mixture 651 652 // internal functions 653 void alloc_internal(); 654 void free_internal(); 655 long int AB( double Tcrit, double Pcrit, double omg, double N, 656 double &apure, double &bpure, double &da, double &d2a ); 657 long int FugacityPT( long int i, double *EoSparam ); 658 long int FugacityPure( long int j ); // Calculates the fugacity of pure species 659 long int Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 ); 660 long int MixParam( double &amix, double &bmix ); 661 long int FugacityMix( double amix, double bmix, double &fugmix, double &zmix, double &vmix ); 662 long int FugacitySpec( double *fugpure ); 663 long int ResidualFunct( double *fugpure ); 664 long int MixingWaals(); 665 long int MixingConst(); 666 long int MixingTemp(); 667 668 public: 669 670 /// Constructor 671 TSRKcalc( long int NCmp, double Pp, double Tkp ); 672 TSRKcalc( SolutionData *sd ); 673 674 /// Destructor 675 ~TSRKcalc(); 676 677 /// Calculates pure species properties (pure fugacities) 678 long int PureSpecies(); 679 680 /// Calculates T,P corrected interaction parameters 681 long int PTparam(); 682 683 /// Calculates activity coefficients 684 long int MixMod(); 685 686 /// Calculates excess properties 687 long int ExcessProp( double *Zex ); 688 689 /// Calculates ideal mixing properties 690 long int IdealProp( double *Zid ); 691 692 /// Calculates pure species properties (called from DCthermo) 693 long int SRKCalcFugPure( double Tmin, double *Cpg, double *FugProps ); 694 695 }; 696 697 698 699 // ------------------------------------------------------------------------------------- 700 /// Peng-Robinson (PR78) model for fluid mixtures. 701 /// References: Peng and Robinson (1976); Peng and Robinson (1978) 702 /// (c) TW July 2009 703 class TPR78calc: public TSolMod 704 705 { 706 private: 707 708 double PhVol; ///< phase volume in cm3 709 double *Pparc; ///< DC partial pressures (pure fugacities) 710 // double *aGEX; // Increments to molar G0 values of DCs from pure fugacities 711 // double *aVol; // DC molar volumes, cm3/mol [L] 712 713 // main work arrays 714 double (*Eosparm)[4]; ///< EoS parameters 715 double (*Pureparm)[4]; ///< Parameters a, b, da/dT, d2a/dT2 for cubic EoS 716 double (*Fugpure)[6]; ///< Fugacity parameters of pure gas species 717 double (*Fugci)[4]; ///< Fugacity parameters of species in the mixture 718 719 double **a; ///< arrays of generic parameters 720 double **b; 721 double **KK; ///< binary interaction parameter 722 double **dKK; ///< derivative of interaction parameter 723 double **d2KK; ///< second derivative 724 double **AA; ///< binary a terms in the mixture 725 726 // internal functions 727 void alloc_internal(); 728 void free_internal(); 729 long int AB( double Tcrit, double Pcrit, double omg, double N, 730 double &apure, double &bpure, double &da, double &d2a ); 731 long int FugacityPT( long int i, double *EoSparam ); 732 long int FugacityPure( long int j ); // Calculates the fugacity of pure species 733 long int Cardano( double a2, double a1, double a0, double &z1, double &z2, double &z3 ); 734 long int MixParam( double &amix, double &bmix ); 735 long int FugacityMix( double amix, double bmix, double &fugmix, double &zmix, double &vmix ); 736 long int FugacitySpec( double *fugpure ); 737 long int ResidualFunct( double *fugpure ); 738 long int MixingWaals(); 739 long int MixingConst(); 740 long int MixingTemp(); 741 742 public: 743 744 /// Constructor 745 TPR78calc( long int NCmp, double Pp, double Tkp ); 746 TPR78calc( SolutionData *sd ); 747 748 /// Destructor 749 ~TPR78calc(); 750 751 /// Calculates pure species properties (pure fugacities) 752 long int PureSpecies(); 753 754 /// Calculates T,P corrected interaction parameters 755 long int PTparam(); 756 757 /// Calculates activity coefficients 758 long int MixMod(); 759 760 /// Calculates excess properties 761 long int ExcessProp( double *Zex ); 762 763 /// Calculates ideal mixing properties 764 long int IdealProp( double *Zid ); 765 766 /// Calculates pure species properties (called from DCthermo) 767 long int PR78CalcFugPure( double Tmin, double *Cpg, double *FugProps ); 768 769 }; 770 771 772 773 // ------------------------------------------------------------------------------------- 774 /// Compensated Redlich-Kwong (CORK) model for fluid mixtures. 775 /// References: Holland and Powell (1991) 776 /// (c) TW May 2010 777 class TCORKcalc: public TSolMod 778 779 { 780 private: 781 782 // constants and external parameters 783 double RR; ///< gas constant in kbar 784 double Pkb; ///< pressure in kbar 785 double PhVol; ///< phase volume in cm3 786 double *Pparc; ///< DC partial pressures (pure fugacities) 787 // double *aGEX; // Increments to molar G0 values of DCs from pure fugacities 788 // double *aVol; // DC molar volumes, cm3/mol [L] 789 790 // internal work data 791 double (*Eosparm)[2]; ///< EoS parameters 792 double (*Fugpure)[6]; ///< Fugacity parameters of pure gas species 793 double (*Fugci)[4]; ///< Fugacity parameters of species in the mixture 794 double (*Rho)[11]; ///< density parameters 795 char *EosCode; ///< identifier of EoS routine 796 double *phi; 797 double *dphi; 798 double *d2phi; 799 double *dphip; 800 double **A; ///< binary interaction parameters 801 double **W; ///< volume scaled interaction parameters (derivatives) 802 double **B; 803 double **dB; 804 double **d2B; 805 double **dBp; 806 807 // internal functions 808 void alloc_internal(); 809 void free_internal(); 810 long int FugacityPT( long int j, double *EoSparam ); 811 long int FugacityH2O( long int j ); 812 long int FugacityCO2( long int j ); 813 long int FugacityCorresponding( long int j ); 814 long int VolumeFugacity( long int phState, double pp, double p0, double a, double b, double c, 815 double d, double e, double &vol, double &fc ); 816 long int Cardano( double cb, double cc, double cd, double &v1, double &v2, double &v3 ); 817 long int FugacityMix(); 818 long int ResidualFunct(); 819 820 public: 821 822 /// Constructor 823 TCORKcalc( long int NCmp, double Pp, double Tkp, char Eos_Code ); 824 TCORKcalc( SolutionData *sd ); 825 826 /// Destructor 827 ~TCORKcalc(); 828 829 /// Calculates pure species properties (pure fugacities) 830 long int PureSpecies(); 831 832 /// Calculates T,P corrected interaction parameters 833 long int PTparam(); 834 835 /// Calculates activity coefficients 836 long int MixMod(); 837 838 /// Calculates excess properties 839 long int ExcessProp( double *Zex ); 840 841 /// Calculates ideal mixing properties 842 long int IdealProp( double *Zid ); 843 844 /// Calculates pure species properties (called from DCthermo) 845 long int CORKCalcFugPure( double Tmin, /*float*/ double *Cpg, double *FugProps ); 846 847 }; 848 849 850 851 // ------------------------------------------------------------------------------------- 852 /// Sterner-Pitzer (STP) model for fluid mixtures. 853 /// References: Sterner and Pitzer (1994) 854 /// (c) TW December 2010 855 class TSTPcalc: public TSolMod 856 857 { 858 private: 859 860 // constants and external parameters 861 double RC, RR, TMIN, TMAX, PMIN, PMAX; 862 double Pkbar, Pkb, Pmpa; 863 double PhVol; ///< phase volume in cm3 864 double *Pparc; ///< DC partial pressures (pure fugacities) 865 866 // internal work data 867 char *EosCode; 868 double *Tc; 869 double *Pc; 870 double *Psat; 871 double *Rhol; 872 double *Rhov; 873 double *Mw; 874 double *Phi; 875 double *dPhiD; 876 double *dPhiDD; 877 double *dPhiT; 878 double *dPhiTT; 879 double *dPhiDT; 880 double *dPhiDDD; 881 double *dPhiDDT; 882 double *dPhiDTT; 883 double (*Fugpure)[7]; 884 double (*Rho)[11]; 885 double *phi; 886 double *dphi; 887 double *d2phi; 888 double *dphip; 889 double *lng; 890 double **cfh; 891 double **cfc; 892 double **A; 893 double **W; 894 double **B; 895 double **dB; 896 double **d2B; 897 double **dBp; 898 899 // internal functions 900 void alloc_internal(); 901 void free_internal(); 902 void set_internal(); 903 long int UpdateTauP(); 904 long int FugacityPT( long int j, double *EoSparam ); 905 long int FugacityH2O( long int j ); 906 long int FugacityCO2( long int j ); 907 long int FugacityCorresponding( long int j ); 908 long int DensityGuess( long int j, double &Delguess ); 909 long int PsatH2O( long int j ); 910 long int PsatCO2( long int j ); 911 long int Pressure( double rho, double &p, double &dpdrho, double **cf ); 912 long int Helmholtz( long int j, double rho, double **cf ); 913 long int ResidualFunct(); 914 915 public: 916 917 /// Constructor 918 TSTPcalc ( long int NCmp, double Pp, double Tkp, char Eos_Code ); 919 TSTPcalc ( SolutionData *sd ); 920 921 /// Destructor 922 ~TSTPcalc(); 923 924 /// Calculates pure species properties (pure fugacities) 925 long int PureSpecies(); 926 927 /// Calculates T,P corrected interaction parameters 928 long int PTparam(); 929 930 /// Calculates activity coefficients 931 long int MixMod(); 932 933 /// Calculates excess properties 934 long int ExcessProp( double *Zex ); 935 936 /// Calculates ideal mixing properties 937 long int IdealProp( double *Zid ); 938 939 /// Calculates pure species properties (called from DCthermo) 940 long int STPCalcFugPure( double Tmin, double *Cpg, double *FugProps ); 941 942 }; 943 944 945 946 // ------------------------------------------------------------------------------------- 947 /// Van Laar model for solid solutions. 948 /// References: Holland and Powell (2003) 949 /// (c) TW March 2007 950 951 class TVanLaar: public TSolMod 952 { 953 private: 954 double *Wu; 955 double *Ws; 956 double *Wv; 957 double *Wpt; ///< Interaction coeffs at P-T 958 double *Phi; ///< Mixing terms 959 double *PsVol; ///< End member volume parameters 960 961 void alloc_internal(); 962 void free_internal(); 963 964 public: 965 966 /// Constructor 967 TVanLaar( SolutionData *sd ); 968 969 /// Destructor 970 ~TVanLaar(); 971 972 /// Calculates T,P corrected interaction parameters 973 long int PTparam(); 974 975 /// Calculates of activity coefficients 976 long int MixMod(); 977 978 /// Calculates excess properties 979 long int ExcessProp( double *Zex ); 980 981 /// Calculates ideal mixing properties 982 long int IdealProp( double *Zid ); 983 984 }; 985 986 987 988 // ------------------------------------------------------------------------------------- 989 /// Regular model for multicomponent solid solutions. 990 /// References: Holland and Powell (1993) 991 /// (c) TW March 2007 992 class TRegular: public TSolMod 993 { 994 private: 995 double *Wu; 996 double *Ws; 997 double *Wv; 998 double *Wpt; ///< Interaction coeffs at P-T 999 1000 void alloc_internal(); 1001 void free_internal(); 1002 1003 public: 1004 1005 /// Constructor 1006 TRegular( SolutionData *sd ); 1007 1008 /// Destructor 1009 ~TRegular(); 1010 1011 /// Calculates T,P corrected interaction parameters 1012 long int PTparam( ); 1013 1014 /// Calculates of activity coefficients 1015 long int MixMod(); 1016 1017 /// Calculates excess properties 1018 long int ExcessProp( double *Zex ); 1019 1020 /// Calculates ideal mixing properties 1021 long int IdealProp( double *Zid ); 1022 1023 }; 1024 1025 1026 1027 // ------------------------------------------------------------------------------------- 1028 /// Redlich-Kister model for multicomponent solid solutions. 1029 /// References: Hillert (1998) 1030 /// (c) TW March 2007 1031 class TRedlichKister: public TSolMod 1032 { 1033 private: 1034 double (*Lu)[4]; 1035 double (*Ls)[4]; 1036 double (*Lcp)[4]; 1037 double (*Lv)[4]; 1038 double (*Lpt)[4]; 1039 1040 void alloc_internal(); 1041 void free_internal(); 1042 1043 public: 1044 1045 /// Constructor 1046 TRedlichKister( SolutionData *sd ); 1047 1048 /// Destructor 1049 ~TRedlichKister(); 1050 1051 /// Calculates T,P corrected interaction parameters 1052 long int PTparam(); 1053 1054 /// Calculates activity coefficients 1055 long int MixMod(); 1056 1057 /// Calculates excess properties 1058 long int ExcessProp( double *Zex ); 1059 1060 /// Calculates ideal mixing properties 1061 long int IdealProp( double *Zid ); 1062 1063 }; 1064 1065 1066 1067 // ------------------------------------------------------------------------------------- 1068 /// Non-random two liquid (NRTL) model for liquid solutions. 1069 /// References: Renon and Prausnitz (1968), Prausnitz et al. (1997) 1070 /// (c) TW June 2008 1071 class TNRTL: public TSolMod 1072 { 1073 private: 1074 double **Tau; 1075 double **dTau; 1076 double **d2Tau; 1077 double **Alp; 1078 double **dAlp; 1079 double **d2Alp; 1080 double **G; 1081 double **dG; 1082 double **d2G; 1083 1084 void alloc_internal(); 1085 void free_internal(); 1086 1087 public: 1088 1089 /// Constructor 1090 TNRTL( SolutionData *sd ); 1091 1092 /// Destructor 1093 ~TNRTL(); 1094 1095 /// Calculates T,P corrected interaction parameters 1096 long int PTparam(); 1097 1098 /// Calculates activity coefficients 1099 long int MixMod(); 1100 1101 /// Calculates excess properties 1102 long int ExcessProp( double *Zex ); 1103 1104 /// Calculates ideal mixing properties 1105 long int IdealProp( double *Zid ); 1106 1107 }; 1108 1109 1110 1111 // ------------------------------------------------------------------------------------- 1112 /// Wilson model for liquid solutions. 1113 /// References: Prausnitz et al. (1997) 1114 /// (c) TW June 2008 1115 class TWilson: public TSolMod 1116 { 1117 private: 1118 double **Lam; 1119 double **dLam; 1120 double **d2Lam; 1121 1122 void alloc_internal(); 1123 void free_internal(); 1124 1125 public: 1126 1127 /// Constructor 1128 TWilson( SolutionData *sd ); 1129 1130 /// Destructor 1131 ~TWilson(); 1132 1133 /// Calculates T,P corrected interaction parameters 1134 long int PTparam(); 1135 1136 /// Calculates activity coefficients 1137 long int MixMod(); 1138 1139 /// Calculates excess properties 1140 long int ExcessProp( double *Zex ); 1141 1142 /// Calculates ideal mixing properties 1143 long int IdealProp( double *Zid ); 1144 1145 }; 1146 1147 1148 1149 // ------------------------------------------------------------------------------------- 1150 /// Berman model for multi-component sublattice solid solutions. 1151 /// To be extended with reciprocal terms. 1152 /// References: Wood and Nicholls (1978); Berman and Brown (1993) 1153 /// (c) DK/TW December 2010, June 2011 1154 class TBerman: public TSolMod 1155 { 1156 private: 1157 long int NrcR; ///< max. possible number of reciprocal reactions (allocated) 1158 long int Nrc; ///< number of reciprocal reactions (actual) 1159 long int *NmoS; ///< number of different moieties (in end members) on each sublattice 1160 long int ***XrcM; ///< Table of indexes of end members, sublattices and moieties involved in 1161 ///< reciprocal reactions [NrecR][4][2], two left and two right side. 1162 ///< for each of 4 reaction components: j, mark, // s1, m1, s2, m2. 1163 1164 double *Wu; ///< Interaction parameter coefficients a 1165 double *Ws; ///< Interaction parameter coefficients b (f(T)) 1166 double *Wv; ///< Interaction parameter coefficients c (f(P)) 1167 double *Wpt; ///< Interaction parameters corrected at P-T of interest 1168 double **fjs; ///< array of site activity coefficients for end members [NComp][NSub] 1169 1170 double *Grc; ///< standard molar reciprocal energies (constant) 1171 double *oGf; ///< molar Gibbs energies of end-member compounds 1172 double *G0f; ///< standard molar Gibbs energies of end members (constant) 1173 double *DGrc; ///< molar effects of reciprocal reactions [NrecR] 1174 double *pyp; ///< Products of site fractions for end members (CEF mod.) [NComp] 1175 // double *pyn; // Products of site fractions for sites not in the end member [NComp] 1176 void alloc_internal(); 1177 void free_internal(); 1178 long int choose( const long int n, const long int k ); 1179 bool CheckThisReciprocalReaction( const long int r, const long int j, long int *xm ); 1180 long int CollectReciprocalReactions2( void ); 1181 // long int CollectReciprocalReactions3( void ); 1182 long int FindIdenticalSublatticeRow(const long int si, const long int ji, const long jp, 1183 const long int jb, const long int je ); 1184 // long int &nsx, long int *sx, long int *mx ); 1185 long int ExcessPart(); 1186 ///< Arrays for ideal conf part must exist in base TSolMod instance 1187 double PYproduct( const long int j ); 1188 long int em_which(const long int s, const long int m , const long int jb, const long int je); 1189 long int em_howmany( long int s, long int m ); 1190 double ysigma( const long int j, const long int s ); 1191 double KronDelta( const long int j, const long int s, const long int m ); 1192 double dGref_dysigma(const long int l, const long int s, const long int ex_j ); 1193 double dGref_dysm( const long int s, const long m, const long int ex_j ); 1194 double RefFrameTerm( const long int j, double G_ref ); 1195 long int ReciprocalPart(); ///< Calculation of reciprocal contributions to activity coefficients 1196 1197 public: 1198 1199 /// Constructor 1200 TBerman( SolutionData *sd, double *G0 ); 1201 1202 /// Destructor 1203 ~TBerman(); 1204 1205 /// Calculates T,P corrected interaction parameters 1206 long int PTparam(); 1207 1208 /// Calculates activity coefficients 1209 long int MixMod(); 1210 1211 /// Calculates excess properties 1212 long int ExcessProp( double *Zex ); 1213 1214 /// Calculates ideal mixing properties 1215 long int IdealProp( double *Zid ); 1216 1217 }; 1218 1219 1220 // ------------------------------------------------------------------------------------- 1221 /// CEF (Calphad) model for multi-component sublattice solid solutions with reciprocal terms 1222 /// References: Sundman & Agren (1981); Lucas et al. (2006); Hillert (1998). 1223 /// (c) DK/SN since August 2014 (still to change the excess Gibbs energy terms). 1224 class TCEFmod: public TSolMod 1225 { 1226 private: 1227 long int *NmoS; ///< number of different moieties (in end members) on each sublattic 1228 1229 double *Wu; ///< Interaction parameter coefficients a 1230 double *Ws; ///< Interaction parameter coefficients b (f(T)) 1231 double *Wc; ///< Interaction parameter coefficients b (f(TlnT)) 1232 double *Wv; ///< Interaction parameter coefficients c (f(P)) 1233 double *Wpt; ///< Interaction parameters corrected at P-T of interest 1234 double **fjs; ///< array of site activity coefficients for end members [NComp][NSub] 1235 1236 double *Grc; ///< standard molar reciprocal energies (constant) 1237 double *oGf; ///< molar Gibbs energies of end-member compounds 1238 double *G0f; ///< standard molar Gibbs energies of end members (constant) 1239 double *pyp; ///< Products of site fractions for end members (CEF mod.) [NComp] 1240 // double *pyn; // Products of site fractions for sites not in the end member [NComp] 1241 void alloc_internal(); 1242 void free_internal(); 1243 long int ExcessPart(); 1244 ///< Arrays for ideal conf part must exist in base TSolMod instance 1245 double PYproduct( const long int j ); 1246 long int em_which(const long int s, const long int m , const long int jb, const long int je); 1247 long int em_howmany( long int s, long int m ); 1248 double ysm( const long int j, const long int s ); 1249 double KronDelta( const long int j, const long int s, const long int m ); 1250 double dGref_dysigma(const long int l, const long int s ); 1251 double dGref_dysm(const long int s, const long m ); 1252 double dGm_dysm(const long int s, const long m ); // added by Nichenko 1253 double RefFrameTerm( const long int j, double G_ref ); 1254 long int ReciprocalPart(); ///< Calculation of reciprocal contributions to activity coefficients 1255 1256 long int IdealMixing(); // NSergii: added by Nichenko to rewrite the ideal part contribution 1257 long int CalcSiteFractions(); // NSergii: 1258 double dGrefdnNum(const long int i); // NSergii: 1259 double dGidmixdnNum(const long int i); // NSergii: 1260 double dGexcdnNum(const long int i); // NSergii: 1261 double Gmix(); // NSergii: 1262 double Gexc(); 1263 double Gref(); 1264 double Gidmix(); 1265 public: 1266 1267 /// Constructor 1268 TCEFmod( SolutionData *sd, double *G0 ); 1269 1270 /// Destructor 1271 ~TCEFmod(); 1272 1273 /// Calculates T,P corrected interaction parameters 1274 long int PTparam(); 1275 1276 /// Calculates activity coefficients 1277 long int MixMod(); 1278 1279 /// Calculates excess properties 1280 long int ExcessProp( double *Zex ); 1281 1282 /// Calculates ideal mixing properties 1283 long int IdealProp( double *Zid ); 1284 1285 }; 1286 1287 1288 // ------------------------------------------------------------------------------------- 1289 /// SIT model reimplementation for aqueous electrolyte solutions. 1290 /// (c) DK/TW June 2009 1291 class TSIT: public TSolMod 1292 { 1293 private: 1294 1295 // data objects copied from MULTI 1296 double *z; ///< vector of species charges (for aqueous models) 1297 double *m; ///< vector of species molalities (for aqueous models) 1298 double *RhoW; ///< water density properties 1299 double *EpsW; ///< water dielectrical properties 1300 1301 // internal work objects 1302 double I; ///< ionic strength 1303 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1304 double *LnG; ///< activity coefficient 1305 double *dLnGdT; ///< derivatives 1306 double *d2LnGdT2; 1307 double *dLnGdP; 1308 double **E0; ///< interaction parameter 1309 double **E1; 1310 double **dE0; 1311 double **dE1; 1312 double **d2E0; 1313 double **d2E1; 1314 1315 // internal functions 1316 double IonicStrength(); 1317 void alloc_internal(); 1318 void free_internal(); 1319 1320 public: 1321 1322 /// Constructor 1323 TSIT( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1324 1325 /// Destructor 1326 ~TSIT(); 1327 1328 /// Calculates activity coefficients 1329 long int MixMod(); 1330 1331 /// Calculates excess properties 1332 long int ExcessProp( double *Zex ); 1333 1334 /// Calculates ideal mixing properties 1335 long int IdealProp( double *Zid ); 1336 1337 /// Calculation of internal tables (at each GEM iteration) 1338 long int PTparam(); 1339 1340 }; 1341 1342 1343 1344 // ------------------------------------------------------------------------------------- 1345 /// Pitzer model, Harvie-Moller-Weare (HMW) version, with explicit temperature dependence. 1346 /// References: 1347 /// (c) SD/FH February 2009 1348 class TPitzer: public TSolMod 1349 { 1350 1351 private: 1352 long int Nc; ///< Number of cations 1353 long int Na; ///< Number of anions 1354 long int Nn; ///< Number of neutral species 1355 long int Ns; ///< Total number of aqueous species (without H2O); index of H2O in aq phase 1356 ///< Conversion of species indexes between aq phase and Pitzer parameter tables 1357 long int *xcx; ///< list of indexes of Nc cations in aqueous phase 1358 long int *xax; ///< list of indexes of Na anions in aq phase 1359 long int *xnx; ///< list of indexes of Nn neutral species in aq phase 1360 double *aZ; ///< Vector of species charges (for aqueous models) 1361 double *zc; 1362 double *za; 1363 double *aM; ///< Vector of species molality (for aqueous models) 1364 double *mc; 1365 double *ma; 1366 double *mn; 1367 double *RhoW; ///< water density properties 1368 double *EpsW; ///< water dielectrical properties 1369 1370 double Aphi, dAphidT, d2AphidT2, dAphidP; ///< Computing A-Factor 1371 double I; ///< Ionic Strength 1372 double Is; ///< Ionic Strength square root 1373 double Ffac; ///< F-Factor 1374 double Zfac; ///< Z-Term 1375 1376 // Input parameter arrays 1377 //for Gex and activity coefficient calculation 1378 double **Bet0; ///< Beta0 table for cation-anion interactions [Nc][Na] 1379 double **Bet1; ///< Beta1 table for cation-anion interactions [Nc][Na] 1380 double **Bet2; ///< Beta2 table for cation-anion interactions [Nc][Na] 1381 double **Cphi; ///< Cphi table for cation-anion interactions [Nc][Na] 1382 double **Lam; ///< Lam table for neutral-cation interactions [Nn][Nc] 1383 double **Lam1; ///< Lam1 table for neutral-anion interactions [Nn][Na] 1384 double **Theta; ///< Theta table for cation-cation interactions [Nc][Nc] 1385 double **Theta1; ///< Theta1 table for anion-anion interactions [Na][Na] 1386 double ***Psi; ///< Psi array for cation-cation-anion interactions [Nc][Nc][Na] 1387 double ***Psi1; ///< Psi1 array for anion-anion-cation interactions [Na][Na][Nc] 1388 double ***Zeta; ///< Zeta array for neutral-cation-anion interactions [Nn][Nc][Na] 1389 1390 1391 // Work parameter arrays 1392 // double *B1; /// B' table for cation-anion interactions corrected for IS [Nc][Na] 1393 // double *B2; /// B table for cation-anion interactions corrected for IS [Nc][Na] 1394 // double *B3; /// B_phi table for cation-anion interactions corrected for IS [Nc][Na] 1395 // double *Phi1; /// Phi' table for anion-anion interactions corrected for IS [Na][Na] 1396 // double *Phi2; /// Phi table for cation-cation interactions corrected for IS [Nc][Nc] 1397 // double *Phi3; /// PhiPhi table for anion-anion interactions corrected for IS [Na][Na] 1398 // double *C; /// C table for cation-anion interactions corrected for charge [Nc][Na] 1399 // double *Etheta; /// Etheta table for cation-cation interactions [Nc][Nc] 1400 // double *Ethetap; /// Etheta' table for anion-anion interactions [Na][Na] 1401 // double bk[21]; /// work space 1402 // double dk[21]; /// work space 1403 1404 /// McInnes parameter array and gamma values 1405 double *McI_PT_array; 1406 double *GammaMcI; 1407 1408 enum eTableType 1409 { 1410 bet0_ = -10, bet1_ = -11, bet2_ = -12, Cphi_ = -20, Lam_ = -30, Lam1_ = -31, 1411 Theta_ = -40, Theta1_ = -41, Psi_ = -50, Psi1_ = -51, Zeta_ = -60 1412 }; 1413 1414 // internal setup 1415 void calcSizes(); 1416 void alloc_internal(); 1417 void free_internal(); 1418 1419 /// build conversion of species indexes between aq phase and Pitzer parameter tables 1420 void setIndexes(); 1421 double setvalue(long int ii, int Gex_or_Sex); 1422 1423 // internal calculations 1424 /// Calculation of Etheta and Ethetap values 1425 void Ecalc( double z, double z1, double I, double DH_term, 1426 double& Etheta, double& Ethetap ); getN()1427 inline long int getN() const 1428 { 1429 return Nc+Na+Nn; 1430 } 1431 1432 double Z_Term( ); 1433 double IonicStr( double& I ); 1434 void getAlp( long int c, long int a, double& alp, double& alp1 ); 1435 double get_g( double x_alp ); 1436 double get_gp( double x_alp ); 1437 double G_ex_par5( long int ii ); 1438 double G_ex_par8( long int ii ); 1439 double S_ex_par5( long int ii ); 1440 double S_ex_par8( long int ii ); 1441 double CP_ex_par5( long int ii ); 1442 double CP_ex_par8( long int ii ); 1443 double F_Factor( double DH_term ); 1444 double lnGammaN( long int N ); 1445 double lnGammaM( long int M, double DH_term ); 1446 double lnGammaX( long int X, double DH_term ); 1447 double lnGammaH2O( double DH_term ); 1448 1449 /// Calc vector of interaction parameters corrected to T,P of interest 1450 void PTcalc( int Gex_or_Sex ); 1451 1452 /// Calculation KCl activity coefficients for McInnes scaling 1453 double McInnes_KCl(); 1454 getIc(long int jj)1455 inline long int getIc( long int jj ) 1456 { 1457 for( long int ic=0; ic<Nc; ic++ ) 1458 if( xcx[ic] == jj ) 1459 return ic; 1460 return -1; 1461 } 1462 getIa(long int jj)1463 inline long int getIa( long int jj ) 1464 { 1465 for( long int ia=0; ia<Na; ia++ ) 1466 if( xax[ia] == jj ) 1467 return ia; 1468 return -1; 1469 } 1470 getIn(long int jj)1471 inline long int getIn( long int jj ) 1472 { 1473 for( long int in=0; in<Nn; in++ ) 1474 if( xnx[in] == jj ) 1475 return in; 1476 return -1; 1477 } 1478 p_sum(double * arr,long int * xx,long int Narr)1479 inline double p_sum( double* arr, long int *xx, long int Narr ) 1480 { 1481 double sum_ =0.; 1482 for( long int i=0; i<Narr; i++ ) 1483 sum_ += arr[xx[i]]; 1484 return sum_; 1485 } 1486 1487 public: 1488 1489 /// Constructor 1490 TPitzer( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1491 1492 /// Destructor 1493 ~TPitzer(); 1494 1495 /// Calculation of T,P corrected interaction parameters 1496 long int PTparam(); 1497 1498 1499 long int MixMod(); 1500 1501 /// Calculates activity coefficients 1502 long int Pitzer_calc_Gamma(); 1503 long int Pitzer_McInnes_KCl(); 1504 1505 /// Calculates excess properties 1506 long int ExcessProp( double *Zex ); 1507 1508 /// Calculates ideal mixing properties 1509 long int IdealProp( double *Zid ); 1510 1511 void Pitzer_test_out( const char *path, double Y ); 1512 1513 }; 1514 1515 1516 1517 // ------------------------------------------------------------------------------------- 1518 /// Extended universal quasi-chemical (EUNIQUAC) model for aqueous electrolyte solutions. 1519 /// References: Nicolaisen et al. (1993), Thomsen et al. (1996), Thomsen (2005) 1520 /// (c) TW/FH May 2009 1521 class TEUNIQUAC: public TSolMod 1522 { 1523 private: 1524 1525 // data objects copied from MULTI 1526 double *z; ///< species charges 1527 double *m; ///< species molalities 1528 double *RhoW; ///< water density properties 1529 double *EpsW; ///< water dielectrical properties 1530 1531 // internal work objects 1532 double *R; ///< volume parameter 1533 double *Q; ///< surface parameter 1534 double *Phi; 1535 double *Theta; 1536 double **U; ///< interaction energies 1537 double **dU; ///< first derivative 1538 double **d2U; ///< second derivative 1539 double **Psi; 1540 double **dPsi; 1541 double **d2Psi; 1542 double IS; ///< ionic strength 1543 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1544 1545 ///< objects needed for debugging output 1546 double gammaDH[200]; 1547 double gammaC[200]; 1548 double gammaR[200]; 1549 1550 // internal functions 1551 void alloc_internal(); 1552 void free_internal(); 1553 long int IonicStrength(); 1554 1555 public: 1556 1557 /// Constructor 1558 TEUNIQUAC( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1559 1560 /// Destructor 1561 ~TEUNIQUAC(); 1562 1563 /// Calculates T,P corrected interaction parameters 1564 long int PTparam(); 1565 1566 /// Calculates activity coefficients 1567 long int MixMod(); 1568 1569 /// Calculates excess properties 1570 long int ExcessProp( double *Zex ); 1571 1572 /// Calculates ideal mixing properties 1573 long int IdealProp( double *Zid ); 1574 1575 void Euniquac_test_out( const char *path ); 1576 1577 }; 1578 1579 // ------------------------------------------------------------------------------------- 1580 // ELVIS activity model for aqueous electrolyte solutions 1581 // (c) FFH Aug 2011 1582 1583 class TELVIS: public TSolMod 1584 { 1585 private: 1586 // data objects copied from MULTI 1587 double *z; // species charges 1588 double *m; // species molalities 1589 double *RhoW; // water density properties 1590 double *EpsW; // water dielectrical properties 1591 double aDH; // averaged ion size term in DH term 1592 double A, dAdT, d2AdT2, dAdP; // A term of DH equation (and derivatives) 1593 double B, dBdT, d2BdT2, dBdP; // B term of DH equation (and derivatives) 1594 1595 double **beta0; 1596 double **beta1; 1597 double **alpha; 1598 1599 double **coord; // coordinaiton number parameter 1600 1601 double **RA; 1602 double **RC; 1603 double **QA; 1604 double **QC; 1605 1606 double CN; 1607 1608 #ifdef ELVIS_SPEED 1609 #define ELVIS_NCOMP 10 1610 // internal work objects 1611 double R[ELVIS_NCOMP]; 1612 double Q[ELVIS_NCOMP]; 1613 double Phi[ELVIS_NCOMP]; 1614 double Theta[ELVIS_NCOMP]; 1615 1616 double EffRad[ELVIS_NCOMP]; 1617 1618 double dRdP[ELVIS_NCOMP]; 1619 double dRdT[ELVIS_NCOMP]; 1620 double d2RdT2[ELVIS_NCOMP]; 1621 double dQdP[ELVIS_NCOMP]; 1622 double dQdT[ELVIS_NCOMP]; 1623 double d2QdT2[ELVIS_NCOMP]; 1624 1625 double WEps[ELVIS_NCOMP][ELVIS_NCOMP]; 1626 double U[ELVIS_NCOMP][ELVIS_NCOMP]; 1627 double dU[ELVIS_NCOMP][ELVIS_NCOMP]; 1628 double d2U[ELVIS_NCOMP][ELVIS_NCOMP]; 1629 double Psi[ELVIS_NCOMP][ELVIS_NCOMP]; 1630 double dPsi[ELVIS_NCOMP][ELVIS_NCOMP]; 1631 double d2Psi[ELVIS_NCOMP][ELVIS_NCOMP]; 1632 double TR[ELVIS_NCOMP][4]; 1633 1634 double U[ELVIS_NCOMP][ELVIS_NCOMP]; 1635 double dUdP[ELVIS_NCOMP][ELVIS_NCOMP]; 1636 double dUdT[ELVIS_NCOMP][ELVIS_NCOMP]; 1637 double d2UdT2[ELVIS_NCOMP][ELVIS_NCOMP]; 1638 1639 double ELVIS_lnGam_DH[ELVIS_NCOMP]; 1640 double ELVIS_lnGam_Born[ELVIS_NCOMP]; 1641 double ELVIS_OsmCoeff_DH[ELVIS_NCOMP]; 1642 double ELVIS_lnGam_UNIQUAC[ELVIS_NCOMP]; 1643 1644 #endif 1645 1646 #ifndef ELVIS_SPEED 1647 // internal work objects 1648 double *R; // volume parameter 1649 double *Q; // surface parameter 1650 double *Phi; 1651 double *Theta; 1652 double *EffRad; // effective ionic radii 1653 double **U; // interaction energies 1654 double **dU; // first derivative 1655 double **d2U; // second derivative 1656 double **Psi; 1657 double **dPsi; 1658 double **d2Psi; 1659 double **TR; // TR interpolation parameter array 1660 double **WEps; // indices for electrolyte specific permittivity calculation 1661 1662 double* dRdP; 1663 double* dRdT; 1664 double* d2RdT2; 1665 double* dQdP; 1666 double* dQdT; 1667 double* d2QdT2; 1668 1669 double** dUdP; 1670 double** dUdT; 1671 double** d2UdT2; 1672 1673 double* ELVIS_lnGam_DH; 1674 double* ELVIS_lnGam_Born; 1675 double* ELVIS_OsmCoeff_DH; 1676 double* ELVIS_lnGam_UNIQUAC; 1677 #endif 1678 1679 double IS; // ionic strength 1680 double molT; // total molality of aqueous species (except water solvent) 1681 double molZ; // total molality of charged species 1682 1683 1684 // objects needed for debugging output 1685 double gammaDH[200]; 1686 double gammaBorn[200]; 1687 double gammaQUAC[200]; 1688 double gammaC[200]; 1689 double gammaR[200]; 1690 1691 // internal functions 1692 void alloc_internal(); 1693 void free_internal(); 1694 1695 long int IonicStrength(); 1696 1697 // activity coefficient contributions 1698 void ELVIS_DH(double* ELVIS_lnGam_DH, double* ELVIS_OsmCoeff_DH); 1699 void ELVIS_Born(double* ELVIS_lnGam_Born); 1700 void ELVIS_UNIQUAC(double* ELVIS_lnGam_UNIQUAC); 1701 1702 // Osmotic coefficient 1703 double Int_OsmCoeff(); 1704 void molfrac_update(); 1705 double FinDiff( double m_j, int j ); // Finite Difference of lnGam of electrolyte 'j' with respect to its molality 'm[j]'; 1706 double CalcWaterAct(); 1707 1708 // Apparent molar volume 1709 // double FinDiffVol( double m_j, void* params ); // Finite differences of mean lnGam with respect to pressure 1710 1711 double trapzd( const double lower_bound, const double upper_bound, int& n, long int& species, int select ); // from Numerical Recipes in C, 2nd Ed. 1712 double qsimp( const double lower_bound, const double upper_bound, long int& species, int select ); // from Numerical Recipes in C, 2nd Ed. 1713 1714 1715 public: 1716 // Constructor 1717 TELVIS( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1718 1719 // Destructor 1720 ~TELVIS(); 1721 1722 // calculates T,P corrected interaction parameters 1723 long int PTparam(); 1724 1725 // calculates activity coefficients and osmotic coefficient by 1726 // numerical Integration of Bjerrum Relation 1727 long int MixMod(); 1728 long int CalcAct(); 1729 1730 // Compute apparent molar volume of electrolyte 1731 double App_molar_volume(); 1732 1733 // calculates excess properties 1734 long int ExcessProp( double *Zex ); 1735 1736 // calculates ideal mixing properties 1737 long int IdealProp( double *Zid ); 1738 1739 // plot debug results 1740 void TELVIS_test_out( const char *path, const double M ) const; 1741 1742 // ELVIS_FIT: get lnGamma array 1743 void get_lnGamma( std::vector<double>& ln_gamma ); 1744 1745 double FinDiffVol( double m_j, int j ); // Finite differences of mean lnGam with respect to pressure 1746 1747 }; 1748 1749 1750 // ------------------------------------------------------------------------------------- 1751 /// Extended Debye-Hueckel (EDH) model for aqueous electrolyte solutions, Helgesons variant. 1752 /// References: Helgeson et al. (1981); Oelkers and Helgeson (1990); Pokrovskii and Helgeson (1995; 1997a; 1997b) 1753 /// (c) TW July 2009 1754 class THelgeson: public TSolMod 1755 { 1756 private: 1757 1758 // status flags copied from MULTI 1759 long int flagH2O; ///< flag for water 1760 long int flagNeut; ///< flag for neutral species 1761 long int flagElect; ///< flag for selection of background electrolyte model 1762 1763 // data objects copied from MULTI 1764 double *z; ///< species charges 1765 double *m; ///< species molalities 1766 double *RhoW; ///< water density properties 1767 double *EpsW; ///< water dielectrical properties 1768 double *an; ///< individual ion size-parameters 1769 double *bg; ///< individual extended-term parameters 1770 double ac; ///< common ion size parameters 1771 double bc; ///< common extended-term parameter 1772 1773 // internal work objects 1774 double ao, daodT, d2aodT2, daodP; ///< ion-size parameter (TP corrected) 1775 double bgam, dbgdT, d2bgdT2, dbgdP; ///< extended-term parameter (TP corrected) 1776 double *LnG; ///< activity coefficient 1777 double *dLnGdT; ///< derivatives 1778 double *d2LnGdT2; 1779 double *dLnGdP; 1780 double IS; ///< ionic strength 1781 double molT; ///< total molality of aqueous species (except water solvent) 1782 double molZ; ///< total molality of charged species 1783 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1784 double B, dBdT, d2BdT2, dBdP; ///< B term of DH equation (and derivatives) 1785 double Gf, dGfdT, d2GfdT2, dGfdP; ///< g function (and derivatives) 1786 1787 // internal functions 1788 void alloc_internal(); 1789 void free_internal(); 1790 long int IonicStrength(); 1791 long int BgammaTP(); 1792 long int IonsizeTP(); 1793 long int Gfunction(); 1794 long int GShok2( double T, double P, double D, double beta, 1795 double alpha, double daldT, double &g, double &dgdP, 1796 double &dgdT, double &d2gdT2 ); 1797 1798 public: 1799 1800 /// Constructor 1801 THelgeson( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1802 1803 /// Destructor 1804 ~THelgeson(); 1805 1806 /// Calculates T,P corrected interaction parameters 1807 long int PTparam(); 1808 1809 /// Calculates activity coefficients 1810 long int MixMod(); 1811 1812 /// Calculates excess properties 1813 long int ExcessProp( double *Zex ); 1814 1815 /// Calculates ideal mixing properties 1816 long int IdealProp( double *Zid ); 1817 1818 /// Set function used in GEMSFITS for mixed electrolites DM 15.08.2014 1819 long int Set_Felect_bc (long int Flagelect, double Bc, double Ac); 1820 1821 }; 1822 1823 1824 1825 // ------------------------------------------------------------------------------------- 1826 /// Extended Debye-Hueckel (EDH) model for aqueous electrolyte solutions, Davies variant. 1827 /// References: Langmuir (1997) 1828 /// (c) TW July 2009 1829 class TDavies: public TSolMod 1830 { 1831 private: 1832 1833 // status flags copied from MULTI 1834 long int flagH2O; ///< flag for water 1835 long int flagNeut; ///< flag for neutral species 1836 long int flagMol; ///< flag for molality correction 1837 1838 // data objects copied from MULTI 1839 double *z; ///< species charges 1840 double *m; ///< species molalities 1841 double *RhoW; ///< water density properties 1842 double *EpsW; ///< water dielectrical properties 1843 1844 // internal work objects 1845 double *LnG; ///< activity coefficient 1846 double *dLnGdT; ///< derivatives 1847 double *d2LnGdT2; 1848 double *dLnGdP; 1849 double IS; ///< ionic strength 1850 double molT; ///< total molality of aqueous species (except water solvent) 1851 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1852 1853 // internal functions 1854 void alloc_internal(); 1855 void free_internal(); 1856 long int IonicStrength(); 1857 1858 public: 1859 1860 /// Constructor 1861 TDavies( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1862 1863 /// Destructor 1864 ~TDavies(); 1865 1866 /// Calculates T,P corrected interaction parameters 1867 long int PTparam(); 1868 1869 /// Calculates activity coefficients 1870 long int MixMod(); 1871 1872 /// Calculates excess properties 1873 long int ExcessProp( double *Zex ); 1874 1875 /// Calculates ideal mixing properties 1876 long int IdealProp( double *Zid ); 1877 1878 }; 1879 1880 1881 1882 // ------------------------------------------------------------------------------------- 1883 /// Debye-Hueckel (DH) limiting law for aqueous electrolyte solutions. 1884 /// References: Langmuir (1997) 1885 /// (c) TW July 2009 1886 class TLimitingLaw: public TSolMod 1887 { 1888 private: 1889 1890 // status flags copied from MULTI 1891 long int flagH2O; ///< flag for water 1892 long int flagNeut; ///< flag for neutral species 1893 1894 // data objects copied from MULTI 1895 double *z; ///< species charges 1896 double *m; ///< species molalities 1897 double *RhoW; ///< water density properties 1898 double *EpsW; ///< water dielectrical properties 1899 1900 // internal work objects 1901 double *LnG; ///< activity coefficient 1902 double *dLnGdT; ///< derivatives 1903 double *d2LnGdT2; 1904 double *dLnGdP; 1905 double IS; ///< ionic strength 1906 double molT; ///< total molality of aqueous species (except water solvent) 1907 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1908 1909 // internal functions 1910 void alloc_internal(); 1911 void free_internal(); 1912 long int IonicStrength(); 1913 1914 public: 1915 1916 /// Constructor 1917 TLimitingLaw( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1918 1919 /// Destructor 1920 ~TLimitingLaw(); 1921 1922 /// calculates T,P corrected interaction parameters 1923 long int PTparam(); 1924 1925 /// Calculates activity coefficients 1926 long int MixMod(); 1927 1928 /// Calculates excess properties 1929 long int ExcessProp( double *Zex ); 1930 1931 /// Calculates ideal mixing properties 1932 long int IdealProp( double *Zid ); 1933 1934 }; 1935 1936 1937 1938 // ------------------------------------------------------------------------------------- 1939 /// Two-term Debye-Hueckel (DH) model for aqueous electrolyte solutions. 1940 /// References: Helgeson et al. (1981) 1941 /// uses individual ion-size parameters, optionally individual salting-out coefficients 1942 /// (c) TW July 2009 1943 class TDebyeHueckel: public TSolMod 1944 { 1945 private: 1946 1947 // status flags copied from MULTI 1948 long int flagH2O; ///< flag for water 1949 long int flagNeut; ///< flag for neutral species 1950 1951 // data objects copied from MULTI 1952 double *z; ///< species charges 1953 double *m; ///< species molalities 1954 double *RhoW; ///< water density properties 1955 double *EpsW; ///< water dielectrical properties 1956 double *an; ///< individual ion size-parameters 1957 double *bg; ///< individual extended-term parameters 1958 double ac; ///< common ion size parameters 1959 double bc; ///< common extended-term parameter 1960 1961 // internal work objects 1962 double ao; ///< average ion-size parameter 1963 double *LnG; ///< activity coefficient 1964 double *dLnGdT; ///< derivatives 1965 double *d2LnGdT2; 1966 double *dLnGdP; 1967 double IS; ///< ionic strength 1968 double molT; ///< total molality of aqueous species (except water solvent) 1969 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 1970 double B, dBdT, d2BdT2, dBdP; ///< B term of DH equation (and derivatives) 1971 1972 // internal functions 1973 void alloc_internal(); 1974 void free_internal(); 1975 long int IonicStrength(); 1976 1977 public: 1978 1979 /// Constructor 1980 TDebyeHueckel( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 1981 1982 /// Destructor 1983 ~TDebyeHueckel(); 1984 1985 /// Calculates T,P corrected interaction parameters 1986 long int PTparam(); 1987 1988 /// Calculates activity coefficients 1989 long int MixMod(); 1990 1991 /// Calculates excess properties 1992 long int ExcessProp( double *Zex ); 1993 1994 /// Calculates ideal mixing properties 1995 long int IdealProp( double *Zid ); 1996 1997 }; 1998 1999 2000 2001 // ------------------------------------------------------------------------------------- 2002 /// Extended Debye-Hueckel (EDH) model for aqueous electrolyte solutions, Karpovs variant. 2003 /// References: Karpov et al. (1997); Helgeson et al. (1981); Oelkers and Helgeson (1990); 2004 /// Pokrovskii and Helgeson (1995; 1997a; 1997b) 2005 /// (c) TW July 2009 2006 class TKarpov: public TSolMod 2007 { 2008 private: 2009 2010 // status flags copied from MULTI 2011 long int flagH2O; ///< flag for water 2012 long int flagNeut; ///< flag for neutral species 2013 long int flagElect; ///< flag for selection of background electrolyte model 2014 2015 // data objects copied from MULTI 2016 double *z; ///< species charges 2017 double *m; ///< species molalities 2018 double *RhoW; ///< water density properties 2019 double *EpsW; ///< water dielectrical properties 2020 double *an; ///< individual ion size-parameters at T,P 2021 double *bg; ///< individual extended-term parameters 2022 double ac; ///< common ion size parameters 2023 double bc; ///< common extended-term parameter 2024 2025 // internal work objects 2026 double ao; ///< average ion-size parameter 2027 double bgam, dbgdT, d2bgdT2, dbgdP; ///< extended-term parameter (TP corrected) 2028 double *LnG; ///< activity coefficient 2029 double *dLnGdT; ///< derivatives 2030 double *d2LnGdT2; 2031 double *dLnGdP; 2032 double IS; ///< ionic strength 2033 double molT; ///< total molality of aqueous species (except water solvent) 2034 double molZ; ///< total molality of charged species 2035 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 2036 double B, dBdT, d2BdT2, dBdP; ///< B term of DH equation (and derivatives) 2037 double Gf, dGfdT, d2GfdT2, dGfdP; ///< g function (and derivatives) 2038 2039 // internal functions 2040 void alloc_internal(); 2041 void free_internal(); 2042 long int IonicStrength(); 2043 long int BgammaTP(); 2044 long int IonsizeTP(); 2045 long int Gfunction(); 2046 long int GShok2( double T, double P, double D, double beta, 2047 double alpha, double daldT, double &g, double &dgdP, 2048 double &dgdT, double &d2gdT2 ); 2049 2050 public: 2051 2052 /// Constructor 2053 TKarpov( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 2054 2055 /// Destructor 2056 ~TKarpov(); 2057 2058 /// Calculates T,P corrected interaction parameters 2059 long int PTparam(); 2060 2061 /// Calculates activity coefficients 2062 long int MixMod(); 2063 2064 /// Calculates excess properties 2065 long int ExcessProp( double *Zex ); 2066 2067 /// Calculates ideal mixing properties 2068 long int IdealProp( double *Zid ); 2069 2070 }; 2071 2072 2073 2074 // ------------------------------------------------------------------------------------- 2075 /// Extended Debye-Hueckel (EDH) model for aqueous electrolyte solutions, Shvarov variant. 2076 /// References: Shvarov (2007); Oelkers and Helgeson (1990); 2077 /// Pokrovskii and Helgeson (1995; 1997a; 1997b) 2078 /// (c) TW July 2009 2079 class TShvarov: public TSolMod 2080 { 2081 private: 2082 2083 // status flags copied from MULTI 2084 long int flagH2O; ///< new flag for water 2085 long int flagNeut; ///< new flag for neutral species 2086 long int flagElect; ///< flag for selection of background electrolyte model 2087 2088 // data objects copied from MULTI 2089 double *z; ///< species charges 2090 double *m; ///< species molalities 2091 double *RhoW; ///< water density properties 2092 double *EpsW; ///< water dielectrical properties 2093 double *bj; ///< individual ion parameters 2094 double ac; ///< common ion size parameters 2095 double bc; ///< common extended-term parameter 2096 2097 // internal work objects 2098 double ao, daodT, d2aodT2, daodP; ///< ion-size parameter (TP corrected) 2099 double bgam, dbgdT, d2bgdT2, dbgdP; ///< extended-term parameter (TP corrected) 2100 double *LnG; ///< activity coefficient 2101 double *dLnGdT; ///< derivatives 2102 double *d2LnGdT2; 2103 double *dLnGdP; 2104 double IS; ///< ionic strength 2105 double molT; ///< total molality of aqueous species (except water solvent) 2106 double A, dAdT, d2AdT2, dAdP; ///< A term of DH equation (and derivatives) 2107 double B, dBdT, d2BdT2, dBdP; ///< B term of DH equation (and derivatives) 2108 double Gf, dGfdT, d2GfdT2, dGfdP; ///< g function (and derivatives) 2109 2110 // internal functions 2111 void alloc_internal(); 2112 void free_internal(); 2113 long int IonicStrength(); 2114 long int BgammaTP(); 2115 long int IonsizeTP(); 2116 long int Gfunction(); 2117 long int GShok2( double T, double P, double D, double beta, 2118 double alpha, double daldT, double &g, double &dgdP, 2119 double &dgdT, double &d2gdT2 ); 2120 2121 public: 2122 2123 /// Constructor 2124 TShvarov( SolutionData *sd, double *arM, double *arZ, double *dW, double *eW ); 2125 2126 /// Destructor 2127 ~TShvarov(); 2128 2129 /// Calculates T,P corrected interaction parameters 2130 long int PTparam(); 2131 2132 /// Calculates activity coefficients 2133 long int MixMod(); 2134 2135 /// Calculates excess properties 2136 long int ExcessProp( double *Zex ); 2137 2138 /// Calculates ideal mixing properties 2139 long int IdealProp( double *Zid ); 2140 2141 }; 2142 2143 2144 2145 // ------------------------------------------------------------------------------------- 2146 /// Class for hardcoded models for solid solutions. 2147 /// (c) TW January 2009 2148 class TModOther: public TSolMod 2149 { 2150 private: 2151 double PhVol; ///< phase volume in cm3 2152 // double *Pparc; /// DC partial pressures/ pure fugacities, bar (Pc by default) [0:L-1] 2153 // double *aGEX; /// Increments to molar G0 values of DCs from pure fugacities or DQF terms, normalized [L] 2154 // double *aVol; /// DC molar volumes, cm3/mol [L] 2155 double *Gdqf; ///< DQF correction terms 2156 double *Hdqf; 2157 double *Sdqf; 2158 double *CPdqf; 2159 double *Vdqf; 2160 2161 void alloc_internal(); 2162 void free_internal(); 2163 2164 public: 2165 2166 /// Constructor 2167 TModOther( SolutionData *sd, double *dW, double *eW ); 2168 2169 /// Destructor 2170 ~TModOther(); 2171 2172 /// Calculates pure species properties (pure fugacities, DQF corrections) 2173 long int PureSpecies(); 2174 2175 /// Calculates T,P corrected interaction parameters 2176 long int PTparam(); 2177 2178 /// Calculates activity coefficients 2179 long int MixMod(); 2180 2181 /// Calculates excess properties 2182 long int ExcessProp( double *Zex ); 2183 2184 /// Calculates ideal mixing properties 2185 long int IdealProp( double *Zid ); 2186 2187 // functions for individual models (under construction) 2188 long int Amphibole1(); 2189 long int Biotite1(); 2190 long int Chlorite1(); 2191 long int Clinopyroxene1(); 2192 long int Feldspar1(); 2193 long int Feldspar2(); 2194 long int Garnet1(); 2195 long int Muscovite1(); 2196 long int Orthopyroxene1(); 2197 long int Staurolite1(); 2198 long int Talc(); 2199 2200 }; 2201 2202 2203 2204 // ------------------------------------------------------------------------------------- 2205 /// Ternary Margules (regular) model for solid solutions. 2206 /// References: Anderson and Crerar (1993); Anderson (2006) 2207 /// (c) TW/DK June 2009 2208 class TMargules: public TSolMod 2209 { 2210 private: 2211 2212 double WU12, WS12, WV12, WG12; 2213 double WU13, WS13, WV13, WG13; 2214 double WU23, WS23, WV23, WG23; 2215 double WU123, WS123, WV123, WG123; 2216 2217 public: 2218 2219 /// Constructor 2220 TMargules( SolutionData *sd ); 2221 2222 /// Destructor 2223 ~TMargules(); 2224 2225 /// Calculates T,P corrected interaction parameters 2226 long int PTparam( ); 2227 2228 /// Calculates of activity coefficients 2229 long int MixMod(); 2230 2231 /// Calculates excess properties 2232 long int ExcessProp( double *Zex ); 2233 2234 /// Calculates ideal mixing properties 2235 long int IdealProp( double *Zid ); 2236 2237 }; 2238 2239 2240 2241 // ------------------------------------------------------------------------------------- 2242 /// Binary Margules (subregular) model for solid solutions. 2243 /// References: Anderson and Crerar (1993); Anderson (2006) 2244 /// (c) TW/DK June 2009, DQF part added by DK on April 5, 2015 2245 class TSubregular: public TSolMod 2246 { 2247 private: 2248 2249 double WU12, WS12, WV12, WG12; 2250 double WU21, WS21, WV21, WG21; 2251 double DQFX1, DQFG1, DQFS1, DQFV1, DQF1; 2252 double DQFX2, DQFG2, DQFS2, DQFV2, DQF2; 2253 2254 public: 2255 2256 /// Constructor 2257 TSubregular( SolutionData *sd ); 2258 2259 /// Destructor 2260 ~TSubregular(); 2261 2262 /// Calculates T,P corrected interaction parameters 2263 long int PTparam( ); 2264 2265 /// Calculates of activity coefficients 2266 long int MixMod(); 2267 2268 /// Calculates excess properties 2269 long int ExcessProp( double *Zex ); 2270 2271 /// Calculates ideal mixing properties 2272 long int IdealProp( double *Zid ); 2273 2274 }; 2275 2276 2277 2278 // ------------------------------------------------------------------------------------- 2279 /// Binary Guggenheim (Redlich-Kister) model for solid solutions. 2280 /// References: Anderson and Crerar (1993); Anderson (2006) 2281 /// uses normalized (by RT) interaction parameters 2282 /// (c) TW/DK June 2009 2283 class TGuggenheim: public TSolMod 2284 { 2285 private: 2286 2287 double a0, a1, a2; 2288 2289 public: 2290 2291 /// Constructor 2292 TGuggenheim( SolutionData *sd ); 2293 2294 /// Destructor 2295 ~TGuggenheim(); 2296 2297 /// Calculates T,P corrected interaction parameters 2298 long int PTparam( ); 2299 2300 /// Calculates of activity coefficients 2301 long int MixMod(); 2302 2303 /// Calculates excess properties 2304 long int ExcessProp( double *Zex ); 2305 2306 /// Calculates ideal mixing properties 2307 long int IdealProp( double *Zid ); 2308 2309 }; 2310 2311 #endif 2312 } 2313 2314 2315 2316 /// _s_solmod_h 2317