From 14b0fe80375cf4692d8bcbc746150684807b8261 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 9 Nov 2024 10:27:38 -0700 Subject: [PATCH] GAS_BINARY_PARAMETERS implemented. gas_binary_parameters test case --- database/phreeqc.dat | 12 +- mytest/gas_binary_parameters | 15 ++ mytest/gas_binary_parameters.out | 192 ++++++++++++++++ mytest/gas_binary_parameters_101.sel | 3 + src/Phreeqc.cpp | 1 + src/Phreeqc.h | 2 + src/PhreeqcKeywords/Keywords.cpp | 4 +- src/PhreeqcKeywords/Keywords.h | 1 + src/gases.cpp | 325 +++++++++++++++++++++++++++ src/prep.cpp | 312 +++++++++++++++++++++++++ src/read.cpp | 98 ++++++++ 11 files changed, 963 insertions(+), 2 deletions(-) create mode 100644 mytest/gas_binary_parameters create mode 100644 mytest/gas_binary_parameters.out create mode 100644 mytest/gas_binary_parameters_101.sel diff --git a/database/phreeqc.dat b/database/phreeqc.dat index dbb3d0814..9ead9b549 100644 --- a/database/phreeqc.dat +++ b/database/phreeqc.dat @@ -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 diff --git a/mytest/gas_binary_parameters b/mytest/gas_binary_parameters new file mode 100644 index 000000000..4a9ee8c54 --- /dev/null +++ b/mytest/gas_binary_parameters @@ -0,0 +1,15 @@ +DATABASE ../database/phreeqc.dat +SELECTED_OUTPUT 101 + -file gas_binary_parameters_101.sel +USER_PUNCH 101 + -headings Mu SC + -start +10 PUNCH STR_F$(MU, 20, 12) +20 PUNCH STR_F$(SC, 20, 10) + -end +SOLUTION +GAS_PHASE +H2O(g) 0 +CO2(g) 1 +Hdg(g) 1 +END \ No newline at end of file diff --git a/mytest/gas_binary_parameters.out b/mytest/gas_binary_parameters.out new file mode 100644 index 000000000..83e20e126 --- /dev/null +++ b/mytest/gas_binary_parameters.out @@ -0,0 +1,192 @@ + Input file: gas_binary_parameters + Output file: gas_binary_parameters.out +Database file: ../database/phreeqc.dat + +------------------ +Reading data base. +------------------ + + SOLUTION_MASTER_SPECIES + SOLUTION_SPECIES + PHASES + GAS_BINARY_PARAMETERS + EXCHANGE_MASTER_SPECIES + EXCHANGE_SPECIES + SURFACE_MASTER_SPECIES + SURFACE_SPECIES + MEAN_GAMMAS + RATES + END +------------------------------------ +Reading input data for simulation 1. +------------------------------------ + + DATABASE ../database/phreeqc.dat + SELECTED_OUTPUT 101 + file gas_binary_parameters_101.sel + USER_PUNCH 101 + headings Mu SC + start + 10 PUNCH STR_F$(MU, 20, 12) + 20 PUNCH STR_F$(SC, 20, 10) + end + SOLUTION + GAS_PHASE + H2O(g) 0 + CO2(g) 1 + Hdg(g) 1 + END +------------------------------------------- +Beginning of initial solution calculations. +------------------------------------------- + +Initial solution 1. + +-----------------------------Solution composition------------------------------ + + Elements Molality Moles + + Pure water + +----------------------------Description of solution---------------------------- + + pH = 7.000 + pe = 4.000 + Specific Conductance (µS/cm, 25°C) = 0 + Density (g/cm³) = 0.99704 + Volume (L) = 1.00297 + Viscosity (mPa s) = 0.89002 + Activity of water = 1.000 + Ionic strength (mol/kgw) = 1.007e-07 + Mass of water (kg) = 1.000e+00 + Total alkalinity (eq/kg) = 1.217e-09 + Temperature (°C) = 25.00 + Electrical balance (eq) = -1.217e-09 + Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.60 + Iterations = 0 + Total H = 1.110124e+02 + Total O = 5.550622e+01 + +----------------------------Distribution of species---------------------------- + + Log Log Log mole V + Species Molality Activity Molality Activity Gamma cm³/mol + + OH- 1.013e-07 1.012e-07 -6.995 -6.995 -0.000 -4.14 + H+ 1.001e-07 1.000e-07 -7.000 -7.000 -0.000 0.00 + H2O 5.551e+01 1.000e+00 1.744 0.000 0.000 18.07 +H(0) 1.416e-25 + H2 7.079e-26 7.079e-26 -25.150 -25.150 0.000 28.61 +O(0) 0.000e+00 + O2 0.000e+00 0.000e+00 -42.080 -42.080 0.000 30.40 + +------------------------------Saturation indices------------------------------- + + Phase SI** log IAP log K(298 K, 1 atm) + + H2(g) -22.05 -25.15 -3.10 H2 + H2O(g) -1.50 0.00 1.50 H2O + O2(g) -39.19 -42.08 -2.89 O2 + +**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm. + For ideal gases, phi = 1. + +----------------------------------------- +Beginning of batch-reaction calculations. +----------------------------------------- + +Reaction step 1. + +Using solution 1. +Using gas phase 1. + +-----------------------------------Gas phase----------------------------------- + +Total pressure: 1.00 atmospheres (Peng-Robinson calculation) + Gas volume: 1.72e+00 liters + Molar volume: 2.44e+01 liters/mole + P * Vm / RT: 0.99861 (Compressibility Factor Z) + + Moles in gas + ---------------------------------- +Component log P P phi Initial Final Delta + +CO2(g) -0.41 3.924e-01 0.996 4.101e-02 2.763e-02 -1.338e-02 +H2O(g) -1.50 3.165e-02 0.992 0.000e+00 2.229e-03 2.229e-03 +Hdg(g) -0.24 5.759e-01 1.001 4.101e-02 4.056e-02 -4.568e-04 + +-----------------------------Solution composition------------------------------ + + Elements Molality Moles + + C 1.338e-02 1.338e-02 + Hdg 4.568e-04 4.568e-04 + +----------------------------Description of solution---------------------------- + + pH = 4.114 Charge balance + pe = 13.747 Adjusted to redox equilibrium + Specific Conductance (µS/cm, 25°C) = 30 + Density (g/cm³) = 0.99716 + Volume (L) = 1.00340 + Viscosity (mPa s) = 0.89111 + Activity of water = 1.000 + Ionic strength (mol/kgw) = 7.767e-05 + Mass of water (kg) = 1.000e+00 + Total alkalinity (eq/kg) = 1.217e-09 + Total CO2 (mol/kg) = 1.338e-02 + Temperature (°C) = 25.00 + Electrical balance (eq) = -1.217e-09 + Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.00 + Iterations = 15 + Total H = 1.110080e+02 + Total O = 5.553075e+01 + +----------------------------Distribution of species---------------------------- + + Log Log Log mole V + Species Molality Activity Molality Activity Gamma cm³/mol + + H+ 7.767e-05 7.689e-05 -4.110 -4.114 -0.004 0.00 + OH- 1.330e-10 1.316e-10 -9.876 -9.881 -0.004 -4.13 + H2O 5.551e+01 9.998e-01 1.744 -0.000 0.000 18.07 +C(-4) 0.000e+00 + CH4 0.000e+00 0.000e+00 -120.374 -120.374 0.000 35.46 +C(4) 1.338e-02 + CO2 1.330e-02 1.330e-02 -1.876 -1.876 0.000 34.43 + HCO3- 7.767e-05 7.689e-05 -4.110 -4.114 -0.004 24.56 + (CO2)2 3.245e-06 3.245e-06 -5.489 -5.489 0.000 68.87 + CO3-2 4.884e-11 4.689e-11 -10.311 -10.329 -0.018 -4.02 +H(0) 2.685e-39 + H2 1.343e-39 1.343e-39 -38.872 -38.872 0.000 28.61 +Hdg 4.568e-04 + Hdg 4.568e-04 4.568e-04 -3.340 -3.340 0.000 28.61 +O(0) 4.623e-15 + O2 2.312e-15 2.312e-15 -14.636 -14.636 0.000 30.40 + +------------------------------Saturation indices------------------------------- + + Phase SI** log IAP log K(298 K, 1 atm) + + CH4(g) -117.57 -120.37 -2.80 CH4 + CO2(g) -0.41 -1.88 -1.47 CO2 Pressure 0.4 atm, phi 0.996 + H2(g) -35.77 -38.87 -3.10 H2 + H2O(g) -1.50 -0.00 1.50 H2O Pressure 0.0 atm, phi 0.992 + Hdg(g) -0.24 -3.34 -3.10 Hdg Pressure 0.6 atm, phi 1.001 + O2(g) -11.74 -14.64 -2.89 O2 + +**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm. + For ideal gases, phi = 1. + +------------------ +End of simulation. +------------------ + +------------------------------------ +Reading input data for simulation 2. +------------------------------------ + +------------------------------- +End of Run after 0.063 Seconds. +------------------------------- + diff --git a/mytest/gas_binary_parameters_101.sel b/mytest/gas_binary_parameters_101.sel new file mode 100644 index 000000000..8871eecae --- /dev/null +++ b/mytest/gas_binary_parameters_101.sel @@ -0,0 +1,3 @@ + Mu SC + 0.000000100661 0.0546997633 + 0.000077673739 30.2878526752 diff --git a/src/Phreeqc.cpp b/src/Phreeqc.cpp index 672b22d2c..c5c892c4b 100644 --- a/src/Phreeqc.cpp +++ b/src/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/Phreeqc.h b/src/Phreeqc.h index ddff8262b..83b01260b 100644 --- a/src/Phreeqc.h +++ b/src/Phreeqc.h @@ -700,6 +700,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 +1675,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/PhreeqcKeywords/Keywords.cpp b/src/PhreeqcKeywords/Keywords.cpp index 8b0e8c000..590fc597c 100644 --- a/src/PhreeqcKeywords/Keywords.cpp +++ b/src/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/PhreeqcKeywords/Keywords.h b/src/PhreeqcKeywords/Keywords.h index 06de5596e..ae10e6fa0 100644 --- a/src/PhreeqcKeywords/Keywords.h +++ b/src/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/gases.cpp b/src/gases.cpp index bfcc64d74..da0cd715b 100644 --- a/src/gases.cpp +++ b/src/gases.cpp @@ -325,6 +325,8 @@ build_fixed_volume_gas(void) } return (OK); } +//#define ORIGINAL +#ifdef ORIGINAL /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: calc_PR(void) @@ -618,6 +620,329 @@ calc_PR(void) } return (V_m); } +#else +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +calc_PR(void) +/* ---------------------------------------------------------------------- */ +/* Calculate fugacity and fugacity coefficient for gas pressures if critical T and P + are defined. + 1) Solve molar volume V_m or total pressure P from Peng-Robinson's EOS: + P = R * T / (V_m - b) - a * aa / (V_m^2 + 2 * b * V_m - b^2) + a = 0.457235 * (R * T_c)^2 / P_c + b = 0.077796 * R * T_c / P_c + aa = (1 + kk * (1 - T_r^0.5))^2 + kk = 0.37464 + 1.54226 * omega - 0.26992 * omega^2 + T_r = T / T_c + multicomponent gas phase: + use: b_sum = Sum(x_i * b), x_i is mole-fraction + a_aa_sum = Sum_i( Sum_j(x_i * x_j * (a_i * aa_i * a_j * aa_j)^0.5) ) + 2) Find the fugacity coefficient phi for gas i: + log(phi_i) = B_ratio * (z - 1) - log(z - B) + A / (2.8284 * B) * (B_ratio - 2 / a_aa_sum * a_aa_sum2) *\ + log((z + 2.4142 * B) / (z - 0.4142 * B)) + B_ratio = b_i / b_sum + A = a_aa_sum * P / R_TK^2 + B = b_sum * P / R_TK + a_aa_sum2 = Sum_j(x_j * (a_aa_i * a_aa_j)^0.5 + 3) correct the solubility of gas i with: + pr_si_f = log10(phi_i) - Delta_V_i * (P - 1) / (2.303 * R * TK); +*/ +{ + LDBLE T_c, P_c; + LDBLE A, B, B_r, /*b2,*/ kk, oo, a_aa, T_r; + LDBLE m_sum, /*b_sum, a_aa_sum,*/ a_aa_sum2; + LDBLE phi; + LDBLE /*R_TK,*/ R = R_LITER_ATM; /* L atm / (K mol) */ + LDBLE r3[4], r3_12, rp, rp3, rq, rz, ri, ri1, one_3 = 0.33333333333333333; + LDBLE disct, vinit, v1, ddp, dp_dv, dp_dv2; + int it; + class phase* phase_ptr; + LDBLE V_m = 0, P = 0; + + LDBLE TK = tk_x; + bool halved; + R_TK = R * TK; + m_sum = b_sum = a_aa_sum = 0.0; + size_t i; + //std::map < std::pair, double > gas_binary_parameters; + //gas_binary_parameters[p = { "H2O(g)", "CO2(g)" }] = 0.19; + //gas_binary_parameters[p = { "H2O(g)", "H2S(g)" }] = 0.19; + //gas_binary_parameters[p = { "H2O(g)", "H2Sg(g)" }] = 0.19; + //gas_binary_parameters[p = { "H2O(g)", "CH4(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "Mtg(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "Methane(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "N2(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "Ntg(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "Ethane(g)" }] = 0.49; + //gas_binary_parameters[p = { "H2O(g)", "Propane(g)" }] = 0.55; + + if (gas_unknowns.size() == 0) + { + error_msg("No gas unknowns.", STOP); + } + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + + for (i = 0; i < gas_unknowns.size(); i++) + { + m_sum += gas_unknowns[i]->moles; + phase_ptr = gas_unknowns[i]->phase; + if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0) + error_msg("Cannot calculate a mixture of ideal and Peng_Robinson gases,\n please define Tc and Pc for the active gases in PHASES.", STOP); + //continue; + if (!phase_ptr->pr_a) + { + T_c = phase_ptr->t_c; + P_c = phase_ptr->p_c; + phase_ptr->pr_a = 0.457235 * R * R * T_c * T_c / P_c; + phase_ptr->pr_b = 0.077796 * R * T_c / P_c; + T_r = TK / T_c; + oo = phase_ptr->omega; + kk = 0.37464 + oo * (1.54226 - 0.26992 * oo); + phase_ptr->pr_alpha = pow(1 + kk * (1 - sqrt(T_r)), 2); + phase_ptr->pr_tk = TK; + phase_ptr->pr_in = true; + } + if (phase_ptr->pr_tk != TK) + { + T_r = TK / phase_ptr->t_c; + oo = phase_ptr->omega; + kk = 0.37464 + oo * (1.54226 - 0.26992 * oo); + phase_ptr->pr_alpha = pow(1 + kk * (1 - sqrt(T_r)), 2); + phase_ptr->pr_tk = TK; + phase_ptr->pr_in = true; + } + } + if (m_sum == 0) + return (OK); + gas_phase_ptr->Set_v_m(gas_phase_ptr->Get_volume() / m_sum); + for (i = 0; i < gas_unknowns.size(); i++) + { + phase_ptr = gas_unknowns[i]->phase; + phase_ptr->fraction_x = gas_unknowns[i]->moles / m_sum; // phase_ptr->fraction_x updated + } + + for (i = 0; i < gas_unknowns.size(); i++) + { + a_aa_sum2 = 0.0; + phase_ptr = gas_unknowns[i]->phase; + //if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0) + // continue; + b_sum += phase_ptr->fraction_x * phase_ptr->pr_b; + size_t i1; + class phase* phase_ptr1; + + for (i1 = 0; i1 < gas_unknowns.size(); i1++) + { + phase_ptr1 = gas_unknowns[i1]->phase; + //if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0) + // continue; + if (phase_ptr1->fraction_x == 0) + continue; + a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * + phase_ptr1->pr_a * phase_ptr1->pr_alpha); + std::pair < std::string, std::string > p; + p = { phase_ptr->name, phase_ptr1->name }; + std::map, double>::const_iterator gas_pair_it; + gas_pair_it = gas_binary_parameters.find(p); + if (gas_pair_it == gas_binary_parameters.end()) + { + p = { phase_ptr1->name, phase_ptr->name }; + gas_pair_it = gas_binary_parameters.find(p); + } + if (gas_pair_it != gas_binary_parameters.end()) + { + a_aa *= (1.0 - gas_pair_it->second); + } + + //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; + } + phase_ptr->pr_aa_sum2 = a_aa_sum2; + } + b2 = b_sum * b_sum; + + if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME) + { + V_m = gas_phase_ptr->Get_volume() / m_sum; + P = 0.0; + while (P <= 0) + { + P = R_TK / (V_m - b_sum) - a_aa_sum / (V_m * (V_m + 2 * b_sum) - b2); + if (P <= 0.0) + { + V_m *= 2.0; + //a_aa_sum /= 2.0; + } + } + if (iterations > 0 && P < 150 && V_m < 1.01) + { + // check for 3-roots... + r3[1] = b_sum - R_TK / P; + r3[2] = -3.0 * b2 + (a_aa_sum - R_TK * 2.0 * b_sum) / P; + r3[3] = b2 * b_sum + (R_TK * b2 - b_sum * a_aa_sum) / P; + // the discriminant of the cubic eqn... + disct = 18. * r3[1] * r3[2] * r3[3] - + 4. * pow(r3[1], 3) * r3[3] + + r3[1] * r3[1] * r3[2] * r3[2] - + 4. * pow(r3[2], 3) - + 27. * r3[3] * r3[3]; + if (disct > 0) + { + // 3-roots, find the largest P... + it = 0; + halved = false; + ddp = 1e-9; + v1 = vinit = 0.729; + dp_dv = f_Vm(v1, this); + while (fabs(dp_dv) > 1e-11 && it < 40) + { + it += 1; + dp_dv2 = f_Vm(v1 - ddp, this); + v1 -= (dp_dv * ddp / (dp_dv - dp_dv2)); + if (!halved && (v1 > vinit || v1 < 0.03)) + { + if (vinit > 0.329) + vinit -= 0.1; + else + vinit -= 0.05; + if (vinit < 0.03) + { + vinit = halve(f_Vm, 0.03, 1.0, 1e-3); + if (f_Vm(vinit - 2e-3, this) < 0) + vinit = halve(f_Vm, vinit + 2e-3, 1.0, 1e-3); + halved = true; + } + v1 = vinit; + } + dp_dv = f_Vm(v1, this); + if (fabs(dp_dv) < 1e-11) + { + if (f_Vm(v1 - 1e-4, this) < 0) + { + v1 = halve(f_Vm, v1 + 1e-4, 1.0, 1e-3); + dp_dv = f_Vm(v1, this); + } + } + } + if (it == 40) + { + // accept a (possible) whobble in the curve... + // error_msg("No convergence when calculating P in Peng-Robinson.", STOP); + } + if (V_m < v1 && it < 40) + P = R_TK / (v1 - b_sum) - a_aa_sum / (v1 * (v1 + 2 * b_sum) - b2); + } + } + if (P <= 0) // iterations = -1 + P = 1.; + gas_phase_ptr->Set_total_p(P); // phase_ptr->total_p updated + gas_phase_ptr->Set_v_m(V_m); + } + else + { + assert(false); + P = gas_phase_ptr->Get_total_p(); + r3[1] = b_sum - R_TK / P; + r3_12 = r3[1] * r3[1]; + r3[2] = -3.0 * b2 + (a_aa_sum - R_TK * 2.0 * b_sum) / P; + r3[3] = b2 * b_sum + (R_TK * b2 - b_sum * a_aa_sum) / P; + // solve t^3 + rp*t + rq = 0. + // molar volume V_m = t - r3[1] / 3... + rp = r3[2] - r3_12 / 3; + rp3 = rp * rp * rp; + rq = (2.0 * r3_12 * r3[1] - 9.0 * r3[1] * r3[2]) / 27 + r3[3]; + rz = rq * rq / 4 + rp3 / 27; + if (rz >= 0) // Cardono's method... + { + ri = sqrt(rz); + if (ri + rq / 2 <= 0) + { + V_m = pow(ri - rq / 2, one_3) + pow(-ri - rq / 2, one_3) - r3[1] / 3; + } + else + { + ri = -pow(ri + rq / 2, one_3); + V_m = ri - rp / (3.0 * ri) - r3[1] / 3; + } + } + else // use complex plane... + { + ri = sqrt(-rp3 / 27); // rp < 0 + ri1 = acos(-rq / 2 / ri); + V_m = 2.0 * pow(ri, one_3) * cos(ri1 / 3) - r3[1] / 3; + } + gas_phase_ptr->Set_v_m(V_m); // phase_ptr->fraction_x updated + } + // calculate the fugacity coefficients... + for (i = 0; i < gas_unknowns.size(); i++) + { + phase_ptr = gas_unknowns[i]->phase; + if (phase_ptr->fraction_x == 0.0) + { + phase_ptr->pr_p = 0; + phase_ptr->pr_phi = 1; + phase_ptr->pr_si_f = 0.0; + //phase_ptr->logk[vm0] = 0.0; // *** + continue; + } + phase_ptr->pr_p = phase_ptr->fraction_x * P; + + //if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0) + //{ + // phase_ptr->pr_phi = 1; + // continue; + //} + rz = P * V_m / R_TK; + A = a_aa_sum * P / (R_TK * R_TK); + B = b_sum * P / R_TK; + B_r = phase_ptr->pr_b / b_sum; + if (rz > B) + { + phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * + log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); + phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); + } + else + phi = -4.6; // fugacity coefficient = 0.01 + phase_ptr->pr_phi = exp(phi); + phase_ptr->pr_si_f = phi / LOG_10; // pr_si_f updated + // **** + //phase_ptr->logk[vm0] = V_m * 1e3 * phase_ptr->fraction_x; + phase_ptr->pr_in = true; + } + return (V_m); +} +#endif /* ---------------------------------------------------------------------- */ int Phreeqc:: diff --git a/src/prep.cpp b/src/prep.cpp index fb867e087..4906cd4fa 100644 --- a/src/prep.cpp +++ b/src/prep.cpp @@ -3738,6 +3738,8 @@ setup_master_rxn(const std::vector &master_ptr_list, const std:: } return (OK); } +//#define ORIGINAL +#ifdef ORIGINAL /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) @@ -4033,7 +4035,317 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) } return (V_m); } +#else +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) +/* ---------------------------------------------------------------------- */ +/* Calculate fugacity and fugacity coefficient for gas pressures if critical T and P + are defined. + 1) Solve molar volume V_m or total pressure P from Peng-Robinson's EOS: + P = R * T / (V_m - b) - a * aa / (V_m^2 + 2 * b * V_m - b^2) + a = 0.457235 * (R * T_c)^2 / P_c + b = 0.077796 * R * T_c / P_c + aa = (1 + kk * (1 - T_r^0.5))^2 + kk = 0.37464 + 1.54226 * omega - 0.26992 * omega^2 + T_r = T / T_c + multicomponent gas phase: + use: b_sum = Sum(x_i * b), x_i is mole-fraction + a_aa_sum = Sum_i( Sum_j(x_i * x_j * (a_i * aa_i * a_j * aa_j)^0.5) ) + 2) Find the fugacity coefficient phi for gas i: + log(phi_i) = B_ratio * (z - 1) - log(z - B) + A / (2.8284 * B) * (B_ratio - 2 / a_aa_sum * a_aa_sum2) *\ + log((z + 2.4142 * B) / (z - 0.4142 * B)) + B_ratio = b_i / b_sum + A = a_aa_sum * P / R_TK^2 + B = b_sum * P / R_TK + a_aa_sum2 = Sum_j(x_j * (a_aa_i * a_aa_j)^0.5 + 3) correct the solubility of gas i with: + pr_si_f = log10(phi_i) - Delta_V_i * (P - 1) / (2.303 * R * TK); +*/ +{ + int i, i1, n_g = (int)phase_ptrs.size(); + LDBLE T_c, P_c; + LDBLE A, B, B_r, /*b2,*/ kk, oo, a_aa, T_r; + LDBLE m_sum, /*b_sum, a_aa_sum,*/ a_aa_sum2; + LDBLE phi; + LDBLE /*R_TK,*/ R = R_LITER_ATM; /* L atm / (K mol) */ + LDBLE r3[4], r3_12, rp, rp3, rq, rz, ri, ri1, one_3 = 0.33333333333333333; + LDBLE disct, vinit, v1, ddp, dp_dv, dp_dv2; + int it; + class phase* phase_ptr, * phase_ptr1; + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + bool halved; + R_TK = R * TK; + m_sum = b_sum = a_aa_sum = 0.0; + for (i = 0; i < n_g; i++) + { + phase_ptr = phase_ptrs[i]; + if (n_g > 1) + { + if (phase_ptr->moles_x == 0) + continue; + m_sum += phase_ptr->moles_x; + } + if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0) + error_msg("Cannot calculate a mixture of ideal and Peng_Robinson gases,\n please define Tc and Pc for the active gases in PHASES.", STOP); + //continue; + if (!phase_ptr->pr_a) + { + T_c = phase_ptr->t_c; + P_c = phase_ptr->p_c; + phase_ptr->pr_a = 0.457235 * R * R * T_c * T_c / P_c; + phase_ptr->pr_b = 0.077796 * R * T_c / P_c; + T_r = TK / T_c; + oo = phase_ptr->omega; + kk = 0.37464 + oo * (1.54226 - 0.26992 * oo); + phase_ptr->pr_alpha = pow(1 + kk * (1 - sqrt(T_r)), 2); + phase_ptr->pr_tk = TK; + // phase_ptr->pr_in = true; + } + if (phase_ptr->pr_tk != TK) + { + T_r = TK / phase_ptr->t_c; + oo = phase_ptr->omega; + kk = 0.37464 + oo * (1.54226 - 0.26992 * oo); + phase_ptr->pr_alpha = pow(1 + kk * (1 - sqrt(T_r)), 2); + phase_ptr->pr_tk = TK; + // phase_ptr->pr_in = true; + } + } + for (i = 0; i < n_g; i++) + { + phase_ptr = phase_ptrs[i]; + if (n_g == 1) + { + phase_ptr->fraction_x = 1.0; + break; + } + if (m_sum == 0) + return (OK); + phase_ptr->fraction_x = phase_ptr->moles_x / m_sum; + } + for (i = 0; i < n_g; i++) + { + a_aa_sum2 = 0.0; + phase_ptr = phase_ptrs[i]; + //if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0) + // continue; + b_sum += phase_ptr->fraction_x * phase_ptr->pr_b; + for (i1 = 0; i1 < n_g; i1++) + { + phase_ptr1 = phase_ptrs[i1]; + //if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0) + // continue; + if (phase_ptr1->fraction_x == 0) + continue; + a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * + phase_ptr1->pr_a * phase_ptr1->pr_alpha); + std::pair < std::string, std::string > p; + p = { phase_ptr->name, phase_ptr1->name }; + std::map, double>::const_iterator gas_pair_it; + gas_pair_it = gas_binary_parameters.find(p); + if (gas_pair_it == gas_binary_parameters.end()) + { + p = { phase_ptr1->name, phase_ptr->name }; + gas_pair_it = gas_binary_parameters.find(p); + } + if (gas_pair_it != gas_binary_parameters.end()) + { + a_aa *= (1.0 - gas_pair_it->second); + } + //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; + } + phase_ptr->pr_aa_sum2 = a_aa_sum2; + } + b2 = b_sum * b_sum; + + if (V_m) + { + P = R_TK / (V_m - b_sum) - a_aa_sum / (V_m * (V_m + 2 * b_sum) - b2); + if (iterations > 0 && P < 150 && V_m < 1.01) + { + // check for 3-roots... + r3[1] = b_sum - R_TK / P; + r3[2] = -3.0 * b2 + (a_aa_sum - R_TK * 2.0 * b_sum) / P; + r3[3] = b2 * b_sum + (R_TK * b2 - b_sum * a_aa_sum) / P; + // the discriminant of the cubic eqn... + disct = 18. * r3[1] * r3[2] * r3[3] - + 4. * pow(r3[1], 3) * r3[3] + + r3[1] * r3[1] * r3[2] * r3[2] - + 4. * pow(r3[2], 3) - + 27. * r3[3] * r3[3]; + //if (iterations > 50) + // it = 0; // debug + if (disct > 0) + { + // 3-roots, find the largest P... + it = 0; + halved = false; + ddp = 1e-9; + v1 = vinit = 0.729; + dp_dv = f_Vm(v1, this); + while (fabs(dp_dv) > 1e-11 && it < 40) + { + it += 1; + dp_dv2 = f_Vm(v1 - ddp, this); + v1 -= (dp_dv * ddp / (dp_dv - dp_dv2)); + if (!halved && (v1 > vinit || v1 < 0.03)) + { + if (vinit > 0.329) + vinit -= 0.1; + else + vinit -= 0.05; + if (vinit < 0.03) + { + vinit = halve(f_Vm, 0.03, 1.0, 1e-3); + if (f_Vm(vinit - 2e-3, this) < 0) + vinit = halve(f_Vm, vinit + 2e-3, 1.0, 1e-3); + halved = true; + } + v1 = vinit; + } + dp_dv = f_Vm(v1, this); + if (fabs(dp_dv) < 1e-11) + { + if (f_Vm(v1 - 1e-4, this) < 0) + { + v1 = halve(f_Vm, v1 + 1e-4, 1.0, 1e-3); + dp_dv = f_Vm(v1, this); + } + } + } + if (it == 40) + { + // accept a (possible) whobble in the curve... + // error_msg("No convergence when calculating P in Peng-Robinson.", STOP); + } + if (V_m < v1 && it < 40) + P = R_TK / (v1 - b_sum) - a_aa_sum / (v1 * (v1 + 2 * b_sum) - b2); + } + } + if (P <= 0) // iterations = -1 + P = 1; + } + else + { + if (P < 1e-10) + P = 1e-10; + r3[1] = b_sum - R_TK / P; + r3_12 = r3[1] * r3[1]; + r3[2] = -3.0 * b2 + (a_aa_sum - R_TK * 2.0 * b_sum) / P; + r3[3] = b2 * b_sum + (R_TK * b2 - b_sum * a_aa_sum) / P; + // solve t^3 + rp*t + rq = 0. + // molar volume V_m = t - r3[1] / 3... + rp = r3[2] - r3_12 / 3; + rp3 = rp * rp * rp; + rq = (2.0 * r3_12 * r3[1] - 9.0 * r3[1] * r3[2]) / 27 + r3[3]; + rz = rq * rq / 4 + rp3 / 27; + if (rz >= 0) // Cardono's method... + { + ri = sqrt(rz); + if (ri + rq / 2 <= 0) + { + V_m = pow(ri - rq / 2, one_3) + pow(-ri - rq / 2, one_3) - r3[1] / 3; + } + else + { + ri = -pow(ri + rq / 2, one_3); + V_m = ri - rp / (3.0 * ri) - r3[1] / 3; + } + } + else // use complex plane... + { + ri = sqrt(-rp3 / 27); // rp < 0 + ri1 = acos(-rq / 2 / ri); + V_m = 2.0 * pow(ri, one_3) * cos(ri1 / 3) - r3[1] / 3; + } + } + // calculate the fugacity coefficients... + for (i = 0; i < n_g; i++) + { + phase_ptr = phase_ptrs[i]; + if (phase_ptr->fraction_x == 0.0) + { + phase_ptr->pr_p = 0; + phase_ptr->pr_phi = 1; + phase_ptr->pr_si_f = 0.0; + continue; + } + phase_ptr->pr_p = phase_ptr->fraction_x * P; + rz = P * V_m / R_TK; + A = a_aa_sum * P / (R_TK * R_TK); + B = b_sum * P / R_TK; + B_r = phase_ptr->pr_b / b_sum; + if (rz > B) + { + phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * + log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); + phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); + //if (phi > 4.44) + // phi = 4.44; + } + else + phi = -4.6; // fugacity coefficient = 0.01 + //if (!strcmp(phase_ptr->name, "H2O(g)") && phi < -4.6) + //{ + //// avoid such phi... + // phi = -4.6; + //} + phase_ptr->pr_phi = exp(phi); + phase_ptr->pr_si_f = phi / LOG_10; + // for initial equilibrations, adapt log_k of the gas phase... + if (state < REACTION) + { + rho_0 = calc_rho_0(TK - 273.15, P); + calc_dielectrics(TK - 273.15, P); + phase_ptr->lk = calc_lk_phase(phase_ptr, TK, P); + } + phase_ptr->pr_in = true; + } + if (gas_phase_ptr && iterations > 2) + { + if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME) + { + gas_phase_ptr->Set_total_p(P); + } + gas_phase_ptr->Set_v_m(V_m); + return (OK); + } + return (V_m); +} +#endif LDBLE Phreeqc:: f_Vm(LDBLE v1, void *cookie) /* ---------------------------------------------------------------------- */ diff --git a/src/read.cpp b/src/read.cpp index bb451db68..5317c5ded 100644 --- a/src/read.cpp +++ b/src/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,101 @@ 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; + } + 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) /* ---------------------------------------------------------------------- */ {