diff --git a/database/Amm.dat b/database/Amm.dat index 17d8d243..604e8e32 100644 --- a/database/Amm.dat +++ b/database/Amm.dat @@ -1,7 +1,3 @@ -# File 1 = C:\GitPrograms\phreeqc3-1\database\Amm.dat, 22/05/2024 19:38, 1948 lines, 55817 bytes, md5=78b3659799b73ddca128328b6ee7533b -# Created 22 May 2024 19:55:37 -# C:\3rdParty\lsp\lsp.exe -f2 -k=asis -ts Amm.dat - # PHREEQC.DAT for calculating temperature and pressure dependence of reactions, and the specific conductance and viscosity of the solution. Based on: # diffusion coefficients and molal volumina of aqueous species, solubility and volume of minerals, and critical temperatures and pressures of gases in Peng-Robinson's EOS. # Details are given at the end of this file. @@ -67,13 +63,13 @@ SOLUTION_SPECIES H+ = H+ -gamma 9 0 -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.57 # for viscosity parameters see ref. 4 - -dw 9.31e-9 838 16.315 0 2.376 24.01 0 -# Dw(25 C) dw_T a a2 visc a3 a_v_dif + -dw 9.31e-9 838 6.96 -2.285 0.206 24.01 0 +# Dw(25 C) dw_T a a2 visc a3 a_v_dif # Dw(TK) = 9.31e-9 * exp(838 / TK - 838 / 298.15) * viscos_0_25 / viscos_0_tc -# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif +# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif for tracer diffusion. -# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 2.376 for H+) -# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Debye-Onsager eqn. (a2 = Vm = 0 for H+, the reference for Vm) +# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 0.206 for H+) +# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Onsager-Falkenhagen eqn. (For H+, the reference ion, vm = v0 = 0, a *= (1 + mu)^a2.) # a3 = -10 ? ka = DH_B * a * mu^a2 (Define a3 = -10, not used in this database.) (a3 = 24.01 for H+, a flag.) # -3 < a3 < 4 ? ka = DH_B * a2 * mu^0.5 / (1 + mu^a3), Appelo, 2017: Dw(I) = Dw(TK) * exp(-a * DH_A * z * sqrt_mu / (1 + ka)) (Sr+2 in this database) @@ -176,7 +172,7 @@ F- = F- Br- = Br- -gamma 3 0 -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 - -viscosity -1.15e-2 -5.75e-2 5.72e-2 1.46e-2 0.116 0.9295 0.82 + -viscosity -6.98e-2 -0.141 1.78e-2 0.159 7.76e-3 6.25e-2 0.859 -dw 2.09e-9 208 3.5 0 0.5737 Zn+2 = Zn+2 -gamma 5 0 @@ -1327,7 +1323,17 @@ Pb(OH)2 389 Pb(OH)2 + 2 H+ = Pb+2 + 2 H2O -log_k 8.15 -delta_h -13.99 kcal - +GAS_BINARY_PARAMETERS +H2O(g) CO2(g) 0.19 +H2O(g) H2S(g) 0.19 +H2O(g) H2Sg(g) 0.19 +H2O(g) CH4(g) 0.49 +H2O(g) Mtg(g) 0.49 +H2O(g) Methane(g) 0.49 +H2O(g) N2(g) 0.49 +H2O(g) Ntg(g) 0.49 +H2O(g) Ethane(g) 0.49 +H2O(g) Propane(g) 0.55 EXCHANGE_MASTER_SPECIES X X- EXCHANGE_SPECIES @@ -1908,10 +1914,11 @@ END # PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). # Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. # Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. -# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are +# These binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are # hard-coded in calc_PR(): # kij CH4 CO2 H2S N2 # H2O 0.49 0.19 0.19 0.49 +# but are overwritten by the data block GAS_BINARY_PARAMETERS of this file. # ============================================================================================= # The molar volumes of solids are entered with # -Vm vm cm3/mol diff --git a/database/phreeqc.dat b/database/phreeqc.dat index 9d281f84..00ca4920 100644 --- a/database/phreeqc.dat +++ b/database/phreeqc.dat @@ -67,13 +67,13 @@ SOLUTION_SPECIES H+ = H+ -gamma 9 0 -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.57 # for viscosity parameters see ref. 4 - -dw 9.31e-9 838 16.315 0 2.376 24.01 0 -# Dw(25 C) dw_T a a2 visc a3 a_v_dif + -dw 9.31e-9 838 6.96 -2.285 0.206 24.01 0 +# Dw(25 C) dw_T a a2 visc a3 a_v_dif # Dw(TK) = 9.31e-9 * exp(838 / TK - 838 / 298.15) * viscos_0_25 / viscos_0_tc -# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif +# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif for tracer diffusion. -# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 2.376 for H+) -# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Debye-Onsager eqn. (a2 = Vm = 0 for H+, the reference for Vm) +# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 0.206 for H+) +# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Onsager-Falkenhagen eqn. (For H+, the reference ion, vm = v0 = 0, a *= (1 + mu)^a2.) # a3 = -10 ? ka = DH_B * a * mu^a2 (Define a3 = -10, not used in this database.) (a3 = 24.01 for H+, a flag.) # -3 < a3 < 4 ? ka = DH_B * a2 * mu^0.5 / (1 + mu^a3), Appelo, 2017: Dw(I) = Dw(TK) * exp(-a * DH_A * z * sqrt_mu / (1 + ka)) (Sr+2 in this database) @@ -176,7 +176,7 @@ F- = F- Br- = Br- -gamma 3 0 -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 - -viscosity -1.15e-2 -5.75e-2 5.72e-2 1.46e-2 0.116 0.9295 0.82 + -viscosity -6.98e-2 -0.141 1.78e-2 0.159 7.76e-3 6.25e-2 0.859 -dw 2.09e-9 208 3.5 0 0.5737 Zn+2 = Zn+2 -gamma 5 0 @@ -1327,7 +1327,17 @@ Pb(OH)2 389 Pb(OH)2 + 2 H+ = Pb+2 + 2 H2O -log_k 8.15 -delta_h -13.99 kcal - +GAS_BINARY_PARAMETERS +H2O(g) CO2(g) 0.19 +H2O(g) H2S(g) 0.19 +H2O(g) H2Sg(g) 0.19 +H2O(g) CH4(g) 0.49 +H2O(g) Mtg(g) 0.49 +H2O(g) Methane(g) 0.49 +H2O(g) N2(g) 0.49 +H2O(g) Ntg(g) 0.49 +H2O(g) Ethane(g) 0.49 +H2O(g) Propane(g) 0.55 EXCHANGE_MASTER_SPECIES X X- EXCHANGE_SPECIES @@ -1908,10 +1918,11 @@ END # PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). # Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. # Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. -# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are +# These binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are # hard-coded in calc_PR(): # kij CH4 CO2 H2S N2 # H2O 0.49 0.19 0.19 0.49 +# but are overwritten by the data block GAS_BINARY_PARAMETERS of this file. # ============================================================================================= # The molar volumes of solids are entered with # -Vm vm cm3/mol diff --git a/database/phreeqc_rates.dat b/database/phreeqc_rates.dat index aa397b64..9ee0f4d6 100644 --- a/database/phreeqc_rates.dat +++ b/database/phreeqc_rates.dat @@ -1,7 +1,3 @@ -# File 1 = C:\GitPrograms\phreeqc3-1\database\phreeqc_rates.dat, 24/05/2024 01:41, 3147 lines, 110328 bytes, md5=7fc916311a573d0ad7ce880f996a9bbf -# Created 24 May 2024 01:58:45 -# C:\3rdParty\lsp\lsp.exe -f2 -k=asis -ts phreeqc_rates.dat - # PHREEQC.DAT for calculating temperature and pressure dependence of reactions, and the specific conductance and viscosity of the solution. Augmented with kinetic rates for minerals from compilations. Based on: # diffusion coefficients and molal volumina of aqueous species, solubility and volume of minerals, and critical temperatures and pressures of gases in Peng-Robinson's EOS. # Details are given at the end of this file. @@ -67,13 +63,13 @@ SOLUTION_SPECIES H+ = H+ -gamma 9 0 -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.57 # for viscosity parameters see ref. 4 - -dw 9.31e-9 838 16.315 0 2.376 24.01 0 -# Dw(25 C) dw_T a a2 visc a3 a_v_dif + -dw 9.31e-9 838 6.96 -2.285 0.206 24.01 0 +# Dw(25 C) dw_T a a2 visc a3 a_v_dif # Dw(TK) = 9.31e-9 * exp(838 / TK - 838 / 298.15) * viscos_0_25 / viscos_0_tc -# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif +# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif for tracer diffusion. -# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 2.376 for H+) -# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Debye-Onsager eqn. (a2 = Vm = 0 for H+, the reference for Vm) +# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 0.206 for H+) +# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Onsager-Falkenhagen eqn. (For H+, the reference ion, vm = v0 = 0, a *= (1 + mu)^a2.) # a3 = -10 ? ka = DH_B * a * mu^a2 (Define a3 = -10, not used in this database.) (a3 = 24.01 for H+, a flag.) # -3 < a3 < 4 ? ka = DH_B * a2 * mu^0.5 / (1 + mu^a3), Appelo, 2017: Dw(I) = Dw(TK) * exp(-a * DH_A * z * sqrt_mu / (1 + ka)) (Sr+2 in this database) @@ -176,7 +172,7 @@ F- = F- Br- = Br- -gamma 3 0 -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 - -viscosity -1.15e-2 -5.75e-2 5.72e-2 1.46e-2 0.116 0.9295 0.82 + -viscosity -6.98e-2 -0.141 1.78e-2 0.159 7.76e-3 6.25e-2 0.859 -dw 2.09e-9 208 3.5 0 0.5737 Zn+2 = Zn+2 -gamma 5 0 @@ -1327,7 +1323,17 @@ Pb(OH)2 389 Pb(OH)2 + 2 H+ = Pb+2 + 2 H2O -log_k 8.15 -delta_h -13.99 kcal - +GAS_BINARY_PARAMETERS +H2O(g) CO2(g) 0.19 +H2O(g) H2S(g) 0.19 +H2O(g) H2Sg(g) 0.19 +H2O(g) CH4(g) 0.49 +H2O(g) Mtg(g) 0.49 +H2O(g) Methane(g) 0.49 +H2O(g) N2(g) 0.49 +H2O(g) Ntg(g) 0.49 +H2O(g) Ethane(g) 0.49 +H2O(g) Propane(g) 0.55 EXCHANGE_MASTER_SPECIES X X- EXCHANGE_SPECIES @@ -3098,10 +3104,11 @@ Wollastonite -6.97 700 56 0.4 0 0 # PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). # Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. # Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. -# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are +# These binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are # hard-coded in calc_PR(): # kij CH4 CO2 H2S N2 # H2O 0.49 0.19 0.19 0.49 +# but are overwritten by the data block GAS_BINARY_PARAMETERS of this file. # ============================================================================================= # The molar volumes of solids are entered with # -Vm vm cm3/mol diff --git a/database/pitzer.dat b/database/pitzer.dat index 508bd1c3..eb1c87f2 100644 --- a/database/pitzer.dat +++ b/database/pitzer.dat @@ -1,7 +1,3 @@ -# File 1 = C:\GitPrograms\phreeqc3-1\database\pitzer.dat, 22/05/2024 19:46, 1033 lines, 38088 bytes, md5=d70476773ed110a269ebbcaf334f1133 -# Created 22 May 2024 19:49:25 -# C:\3rdParty\lsp\lsp.exe -f2 -k=asis -ts pitzer.dat - # Pitzer.DAT for calculating temperature and pressure dependence of reactions, and the specific conductance and viscosity of the solution, using # diffusion coefficients of species, molal volumina of aqueous species and minerals, and critical temperatures and pressures of gases used in Peng-Robinson's EOS. # Details are given at the end of this file. @@ -40,13 +36,13 @@ Ntg Ntg 0 Ntg 28.0134 # N2 gas SOLUTION_SPECIES H+ = H+ -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.57 # for viscosity parameters see ref. 4 - -dw 9.31e-9 823 5.55 0 3.07 24.01 0 -# Dw(25 C) dw_T a a2 visc a3 a_v_dif + -dw 9.31e-9 838 4.02 -1.836 0.415 24.01 0 +# Dw(25 C) dw_T a a2 visc a3 a_v_dif # Dw(TK) = 9.31e-9 * exp(823 / TK - 823 / 298.15) * viscos_0_25 / viscos_0_tc -# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif +# a = DH ion size, a2 = exponent, visc = viscosity exponent, a3(H+) = 24.01 = new dw calculation from A.D. 2024, a_v_dif = exponent in (viscos_0_tc / viscos)^a_v_dif for tracer diffusion. -# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 3.07 for H+) -# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Debye-Onsager eqn. (a2 = Vm = 0 for H+, the reference for Vm) +# For SC, Dw(TK) *= (viscos_0_tc / viscos)^visc (visc = 0.415 for H+) +# a3 > 5 or a3 = 0 or not defined ? ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5, in Onsager-Falkenhagen eqn. (For H+, the reference ion, vm = v0 = 0, a *= (1 + mu)^a2.) # a3 = -10 ? ka = DH_B * a * mu^a2 in DHO. (Define a3 = -10.) # -5 < a3 < 5 ? ka = DH_B * a2 * mu^0.5 / (1 + mu^a3), Appelo, 2017: Dw(I) = Dw(TK) * exp(-a * DH_A * z * sqrt_mu / (1 + ka)) @@ -107,7 +103,7 @@ B(OH)3 = B(OH)3 -dw 1.1e-9 Br- = Br- -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 - -viscosity -1.16e-2 -5.23e-2 5.54e-2 1.22e-2 0.119 0.9969 0.818 + -viscosity -6.98e-2 -0.141 1.78e-2 0.159 7.76e-3 6.25e-2 0.859 -dw 2.09e-9 208 3.5 0 0.5737 H4SiO4 = H4SiO4 -Vm 10.5 1.7 20 -2.7 0.1291 # supcrt 2*H2O in a1 @@ -787,6 +783,17 @@ PITZER K+ OH- SO4-2 -0.05 Mg+2 Na+ SO4-2 -0.015 Na+ OH- SO4-2 -0.009 +GAS_BINARY_PARAMETERS +H2O(g) CO2(g) 0.19 +H2O(g) H2S(g) 0.19 +H2O(g) H2Sg(g) 0.19 +H2O(g) CH4(g) 0.49 +H2O(g) Mtg(g) 0.49 +H2O(g) Methane(g) 0.49 +H2O(g) N2(g) 0.49 +H2O(g) Ntg(g) 0.49 +H2O(g) Ethane(g) 0.49 +H2O(g) Propane(g) 0.55 EXCHANGE_MASTER_SPECIES X X- EXCHANGE_SPECIES @@ -980,10 +987,11 @@ END # PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). # Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. # Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. -# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are +# These binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are # hard-coded in calc_PR(): # kij CH4 CO2 H2S N2 # H2O 0.49 0.19 0.19 0.49 +# but are overwritten by the data block GAS_BINARY_PARAMETERS of this file. # ============================================================================================= # The molar volumes of solids are entered with # -Vm vm cm3/mol diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.cpp index 672b22d2..c5c892c4 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.cpp @@ -1756,6 +1756,7 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc) b2 = pSrc->b2; b_sum = pSrc->b_sum; R_TK = pSrc->R_TK; + gas_binary_parameters = pSrc->gas_binary_parameters; /* input.cpp ------------------------------- */ check_line_return = 0; reading_db = FALSE; diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.h b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.h index ddff8262..cd542516 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.h +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.h @@ -473,6 +473,7 @@ class Phreeqc #endif int calc_gas_pressures(void); int calc_fixed_volume_gas_pressures(void); + double calc_gas_binary_parameter(std::string name1, std::string name2) const; int calc_ss_fractions(void); int gammas(LDBLE mu); int gammas_a_f(int i); @@ -700,6 +701,7 @@ class Phreeqc int read_rate_parameters_svd(void); int read_rate_parameters_hermanska(void); int read_mean_gammas(void); + int read_gas_binary_parameters(void); int read_mix(void); int read_entity_mix(std::map& mix_map); //int read_solution_mix(void); @@ -1674,6 +1676,7 @@ class Phreeqc std::vector x_arg, res_arg, scratch; /* gases.cpp ------------------------------- */ LDBLE a_aa_sum, b2, b_sum, R_TK; + std::map < std::pair, double > gas_binary_parameters; /* input.cpp ------------------------------- */ int check_line_return; diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.cpp index 8b0e8c00..590fc597 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.cpp @@ -134,6 +134,7 @@ std::map::value_type("rate_parameters_pk" std::map::value_type("rate_parameters_svd", Keywords::KEY_RATE_PARAMETERS_SVD), std::map::value_type("rate_parameters_hermanska", Keywords::KEY_RATE_PARAMETERS_HERMANSKA), std::map::value_type("mean_gammas", Keywords::KEY_MEAN_GAMMAS), +std::map::value_type("gas_binary_parameters", Keywords::KEY_GAS_BINARY_PARAMETERS), std::map::value_type("solution_mix", Keywords::KEY_SOLUTION_MIX), std::map::value_type("mix_solution", Keywords::KEY_SOLUTION_MIX), std::map::value_type("exchange_mix", Keywords::KEY_EXCHANGE_MIX), @@ -228,7 +229,8 @@ std::map::value_type(Keywords::KEY_REACTI std::map::value_type(Keywords::KEY_RATE_PARAMETERS_PK, "RATE_PARAMETERS_PK"), std::map::value_type(Keywords::KEY_RATE_PARAMETERS_SVD, "RATE_PARAMETERS_SVD"), std::map::value_type(Keywords::KEY_RATE_PARAMETERS_HERMANSKA, "RATE_PARAMETERS_HERMANSKA"), -std::map::value_type(Keywords::KEY_MEAN_GAMMAS, "RATE_MEAN_GAMMAS"), +std::map::value_type(Keywords::KEY_MEAN_GAMMAS, "MEAN_GAMMAS"), +std::map::value_type(Keywords::KEY_GAS_BINARY_PARAMETERS, "GAS_BINARY_PARAMETERS"), std::map::value_type(Keywords::KEY_SOLUTION_MIX, "SOLUTION_MIX"), std::map::value_type(Keywords::KEY_EXCHANGE_MIX, "EXCHANGE_MIX"), std::map::value_type(Keywords::KEY_GAS_PHASE_MIX, "GAS_PHASE_MIX"), diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.h b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.h index 06de5596..ae10e6fa 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.h +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.h @@ -80,6 +80,7 @@ class Keywords KEY_RATE_PARAMETERS_SVD, KEY_RATE_PARAMETERS_HERMANSKA, KEY_MEAN_GAMMAS, + KEY_GAS_BINARY_PARAMETERS, KEY_SOLUTION_MIX, KEY_EXCHANGE_MIX, KEY_GAS_PHASE_MIX, diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/basicsubs.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/basicsubs.cpp index 917353da..adac2736 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/basicsubs.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/basicsubs.cpp @@ -420,7 +420,7 @@ calc_SC(void) q = 1 / ((t1 / z_plus + (1 - t1) / z_min) * (z_min + z_plus)); sqrt_q = sqrt(q); - // B1 = relaxtion, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka)) + // B1 = relaxation, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka)) a = 1.60218e-19 * 1.60218e-19 / (6 * pi); B1 = a / (2 * 8.8542e-12 * eps_r * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10) B2 = a * AVOGADRO / viscos_0 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4) @@ -486,8 +486,10 @@ calc_SC(void) //av += 0 * t1; } Dw *= Dw_SC * l_z; - if (!a2 || !strcmp(s_x[i]->name, "H+")) + if (!a2) t1 = 1; + else if (!strcmp(s_x[i]->name, "H+")) + t1 = pow(1 + mu_x, a2); else { v0 = calc_vm0(s_x[i]->name, tc_x, 1, 0); diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/gases.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/gases.cpp index bfcc64d7..dea1e551 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/gases.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/gases.cpp @@ -432,36 +432,37 @@ calc_PR(void) continue; a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * phase_ptr1->pr_a * phase_ptr1->pr_alpha); - if (!strcmp(phase_ptr->name, "H2O(g)")) - { - if (!strcmp(phase_ptr1->name, "CO2(g)")) - a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Propane(g)")) - a_aa *= 0.45; - } - if (!strcmp(phase_ptr1->name, "H2O(g)")) - { - if (!strcmp(phase_ptr->name, "CO2(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Propane(g)")) - a_aa *= 0.45; - } + a_aa *= calc_gas_binary_parameter(phase_ptr->name, phase_ptr1->name); + //if (!strcmp(phase_ptr->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr1->name, "CO2(g)")) + // a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + // else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Propane(g)")) + // a_aa *= 0.45; + //} + //if (!strcmp(phase_ptr1->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr->name, "CO2(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Propane(g)")) + // a_aa *= 0.45; + //} a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; a_aa_sum2 += phase_ptr1->fraction_x * a_aa; } @@ -566,12 +567,14 @@ calc_PR(void) if (ri + rq / 2 <= 0) { V_m = pow(ri - rq / 2, one_3) + pow(- ri - rq / 2, one_3) - r3[1] / 3; - } else + } + else { ri = - pow(ri + rq / 2, one_3); V_m = ri - rp / (3.0 * ri) - r3[1] / 3; } - } else // use complex plane... + } + else // use complex plane... { ri = sqrt(- rp3 / 27); // rp < 0 ri1 = acos(- rq / 2 / ri); @@ -695,3 +698,52 @@ calc_fixed_volume_gas_pressures(void) return (OK); } +/* ---------------------------------------------------------------------- */ +double Phreeqc:: +calc_gas_binary_parameter(std::string name1, std::string name2) const +/* ---------------------------------------------------------------------- */ +{ + double f = 1.0; + std::pair < std::string, std::string > p; + p = { name1, name2 }; + std::map, double>::const_iterator gas_pair_it; + gas_pair_it = gas_binary_parameters.find(p); + if (gas_pair_it != gas_binary_parameters.end()) + { + f = (1.0 - gas_pair_it->second); + } + else + { + if (!strcmp(name1.c_str(), "H2O(g)")) + { + if (!strcmp(name2.c_str(), "CO2(g)")) + f = 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + else if (!strcmp(name2.c_str(), "H2S(g)") || !strcmp(name2.c_str(), "H2Sg(g)")) + f = 0.81; + else if (!strcmp(name2.c_str(), "CH4(g)") || !strcmp(name2.c_str(), "Mtg(g)") || !strcmp(name2.c_str(), "Methane(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "N2(g)") || !strcmp(name2.c_str(), "Ntg(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "Ethane(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "Propane(g)")) + f = 0.45; + } + if (!strcmp(name2.c_str(), "H2O(g)")) + { + if (!strcmp(name1.c_str(), "CO2(g)")) + f = 0.81; + else if (!strcmp(name1.c_str(), "H2S(g)") || !strcmp(name1.c_str(), "H2Sg(g)")) + f = 0.81; + else if (!strcmp(name1.c_str(), "CH4(g)") || !strcmp(name1.c_str(), "Mtg(g)") || !strcmp(name1.c_str(), "Methane(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "N2(g)") || !strcmp(name1.c_str(), "Ntg(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "Ethane(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "Propane(g)")) + f = 0.45; + } + } + return f; +} \ No newline at end of file diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/prep.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/prep.cpp index fb867e08..9ab2668e 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/prep.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/prep.cpp @@ -3843,36 +3843,37 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) continue; a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * phase_ptr1->pr_a * phase_ptr1->pr_alpha); - if (!strcmp(phase_ptr->name, "H2O(g)")) - { - if (!strcmp(phase_ptr1->name, "CO2(g)")) - a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Propane(g)")) - a_aa *= 0.45; - } - if (!strcmp(phase_ptr1->name, "H2O(g)")) - { - if (!strcmp(phase_ptr->name, "CO2(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Propane(g)")) - a_aa *= 0.45; - } + a_aa *= calc_gas_binary_parameter(phase_ptr->name, phase_ptr1->name); + //if (!strcmp(phase_ptr->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr1->name, "CO2(g)")) + // a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + // else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Propane(g)")) + // a_aa *= 0.45; + //} + //if (!strcmp(phase_ptr1->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr->name, "CO2(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Propane(g)")) + // a_aa *= 0.45; + //} a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; a_aa_sum2 += phase_ptr1->fraction_x * a_aa; } @@ -3946,7 +3947,8 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) } if (P <= 0) // iterations = -1 P = 1; - } else + } + else { if (P < 1e-10) P = 1e-10; diff --git a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/read.cpp b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/read.cpp index bb451db6..1b0bb84d 100644 --- a/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/read.cpp +++ b/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/read.cpp @@ -147,6 +147,9 @@ read_input(void) case Keywords::KEY_MEAN_GAMMAS: read_mean_gammas(); break; + case Keywords::KEY_GAS_BINARY_PARAMETERS: + read_gas_binary_parameters(); + break; case Keywords::KEY_SOLUTION_MIX: //read_solution_mix(); read_entity_mix(Rxn_solution_mix_map); @@ -2551,6 +2554,103 @@ read_mean_gammas(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +read_gas_binary_parameters(void) +/* ---------------------------------------------------------------------- */ +{ + /* + * Reads GAS_BINARY_PARAMETERS data + * + * Arguments: + * none + * + * Returns: + * KEYWORD if keyword encountered, input_error may be incremented if + * a keyword is encountered in an unexpected position + * EOF if eof encountered while reading mass balance concentrations + * ERROR if error occurred reading data + * + */ + std::string token; + int return_value, opt; + const char* next_char; + const char* opt_list[] = { + "xxxx", /* 0 */ + }; + int count_opt_list = 0; + /* + * Read rate parameters + */ + return_value = UNKNOWN; + for (;;) + { + opt = get_option(opt_list, count_opt_list, &next_char); + switch (opt) + { + case OPTION_EOF: /* end of file */ + return_value = EOF; + break; + case OPTION_KEYWORD: /* keyword */ + return_value = KEYWORD; + break; + case OPTION_DEFAULT: /* add to gas_binary_parameters map */ + { + bool error = false; + std::string gas1, gas2; + int j = copy_token(token, &next_char); + if (j != EMPTY) + { + gas1 = token; + } + else + { + error = true; + } + j = copy_token(token, &next_char); + if (j != EMPTY) + { + gas2 = token; + } + else + { + error = true; + } + j = copy_token(token, &next_char); + double d; + if (j != EMPTY) + { + j = sscanf(token.c_str(), SCANFORMAT, &d); + } + else + { + error = true; + } + if (!error) + { + std::pair p; + p = { gas1, gas2 }; + gas_binary_parameters[p] = d; + p = { gas2, gas1 }; + gas_binary_parameters[p] = d; + } + else + { + error_msg("Error reading gas binary parameter", CONTINUE); + } + } + break; + case OPTION_ERROR: + input_error++; + error_msg("Unknown input in GAS_BINARY_PARAMETERS keyword.", CONTINUE); + error_msg(line_save, CONTINUE); + break; + } + if (return_value == EOF || return_value == KEYWORD) + break; + } + return (return_value); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: read_rate_parameters_svd(void) /* ---------------------------------------------------------------------- */ {