diff --git a/Phreeqc.cpp b/Phreeqc.cpp index 817bfa262..f67a2be5b 100644 --- a/Phreeqc.cpp +++ b/Phreeqc.cpp @@ -583,7 +583,6 @@ void Phreeqc::init(void) solution_pe_x = 0; mu_x = 0; ah2o_x = 1.0; - density_x = 0; total_h_x = 0; total_o_x = 0; cb_x = 0; @@ -898,6 +897,7 @@ void Phreeqc::init(void) viscos = 0.0; viscos_0 = 0.0; viscos_0_25 = 0.0; + density_x = 0.0; rho_0 = 0.0; kappa_0 = 0.0; p_sat = 0.0; @@ -1714,6 +1714,7 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc) viscos = pSrc->viscos; viscos_0 = pSrc->viscos_0; viscos_0_25 = pSrc->viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C + density_x = pSrc->density_x; cell_pore_volume = pSrc->cell_pore_volume;; cell_porosity = pSrc->cell_porosity; cell_volume = pSrc->cell_volume; @@ -1722,9 +1723,6 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc) sys_tot = pSrc->sys_tot; // solution properties V_solutes = pSrc->V_solutes; - viscos = pSrc->viscos; - viscos_0 = pSrc->viscos_0; - viscos_0_25 = pSrc->viscos_0_25; rho_0 = pSrc->rho_0; kappa_0 = pSrc->kappa_0; p_sat = pSrc->p_sat; diff --git a/Phreeqc.h b/Phreeqc.h index 5510379f7..7a1e7e9fe 100644 --- a/Phreeqc.h +++ b/Phreeqc.h @@ -283,7 +283,7 @@ class Phreeqc int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1); int calc_all_donnan(void); int calc_init_donnan(void); - LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector &zcorr); + LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, LDBLE f_free, std::vector &zcorr); LDBLE g_function(LDBLE x_value); LDBLE midpnt(LDBLE x1, LDBLE x2, int n); void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv, @@ -555,6 +555,7 @@ class Phreeqc LDBLE calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m); LDBLE calc_PR(); int calc_vm(LDBLE tc, LDBLE pa); + LDBLE calc_vm0(const char *species_name, LDBLE tc, LDBLE pa, LDBLE mu); int clear(void); int convert_units(cxxSolution* solution_ptr); class unknown* find_surface_charge_unknown(std::string& str_ptr, int plane); @@ -995,7 +996,7 @@ class Phreeqc LDBLE new_Dw); int reformat_surf(const char* comp_name, LDBLE fraction, const char* new_comp_name, LDBLE new_Dw, int cell); - LDBLE viscosity(void); + LDBLE viscosity(cxxSurface *surf_ptr); LDBLE calc_f_visc(const char *name); LDBLE calc_vm_Cl(void); int multi_D(LDBLE DDt, int mobile_cell, int stagnant); @@ -1274,7 +1275,6 @@ class Phreeqc LDBLE solution_pe_x; LDBLE mu_x; LDBLE ah2o_x; - LDBLE density_x; LDBLE total_h_x; LDBLE total_o_x; LDBLE cb_x; @@ -1614,6 +1614,7 @@ class Phreeqc int print_viscosity; LDBLE viscos, viscos_0, viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C + LDBLE density_x; LDBLE cell_pore_volume; LDBLE cell_porosity; LDBLE cell_volume; diff --git a/Solution.cxx b/Solution.cxx index bc1bd80eb..a326adf56 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -81,7 +81,7 @@ cxxSolution::operator =(const cxxSolution &rhs) this->total_h = rhs.total_h; this->total_o = rhs.total_o; this->density = rhs.density; - this->viscosity = rhs.viscosity; + this->viscosity = rhs.viscosity; this->cb = rhs.cb; this->mass_water = rhs.mass_water; this->soln_vol = rhs.soln_vol; diff --git a/Solution.h b/Solution.h index eef728d46..c4326a8a9 100644 --- a/Solution.h +++ b/Solution.h @@ -49,8 +49,8 @@ class cxxSolution:public cxxNumKeyword void Set_cb(LDBLE l_cb) {this->cb = l_cb;} LDBLE Get_density() const {return this->density;} void Set_density(LDBLE l_density) {this->density = l_density;} - LDBLE Get_viscosity() const { return this->viscosity; } - void Set_viscosity(LDBLE l_viscos) { this->viscosity = l_viscos; } + LDBLE Get_viscosity() const { return this->viscosity; } + void Set_viscosity(LDBLE l_viscos) { this->viscosity = l_viscos; } LDBLE Get_mass_water() const {return this->mass_water;} void Set_mass_water(LDBLE l_mass_water) {this->mass_water = l_mass_water;} LDBLE Get_total_alkalinity() const {return this->total_alkalinity;} diff --git a/Surface.cxx b/Surface.cxx index ed55bb2d6..f3f932a40 100644 --- a/Surface.cxx +++ b/Surface.cxx @@ -36,9 +36,10 @@ cxxSurface::cxxSurface(PHRQ_io *io) dl_type = NO_DL; sites_units = SITES_ABSOLUTE; only_counter_ions = false; - correct_GC = false; + correct_D = false; thickness = 1e-8; debye_lengths = 0.0; + calc_DDL_viscosity = false; DDL_viscosity = 1.0; DDL_limit = 0.8; transport = false; @@ -56,9 +57,10 @@ cxxNumKeyword(io) dl_type = NO_DL; sites_units = SITES_ABSOLUTE; only_counter_ions = false; - correct_GC = false; + correct_D = false; thickness = 1e-8; debye_lengths = 0.0; + calc_DDL_viscosity = false; DDL_viscosity = 1.0; DDL_limit = 0.8; transport = false; @@ -130,7 +132,7 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons s_oss << indent1; s_oss << "-only_counter_ions " << this->only_counter_ions << "\n"; s_oss << indent1; - s_oss << "-correct_GC " << this->correct_GC << "\n"; + s_oss << "-correct_D " << this->correct_D << "\n"; s_oss << indent1; s_oss << "-thickness " << this->thickness << "\n"; s_oss << indent1; @@ -193,7 +195,7 @@ cxxSurface::read_raw(CParser & parser, bool check) this->Set_tidied(true); bool only_counter_ions_defined(false); - //bool correct_GC_defined(false); + //bool correct_D_defined(false); bool thickness_defined(false); bool type_defined(false); bool dl_type_defined(false); @@ -395,7 +397,7 @@ cxxSurface::read_raw(CParser & parser, bool check) case 11: // DDL_viscosity if (!(parser.get_iss() >> this->DDL_viscosity)) { - this->DDL_viscosity = 0.0; + this->DDL_viscosity = 1.0; parser.incr_input_error(); parser.error_msg("Expected numeric value for DDL_viscosity.", PHRQ_io::OT_CONTINUE); @@ -473,16 +475,16 @@ cxxSurface::read_raw(CParser & parser, bool check) PHRQ_io::OT_CONTINUE); } break; - case 19: // correct_GC - if (!(parser.get_iss() >> this->correct_GC)) + case 19: // correct_D + if (!(parser.get_iss() >> this->correct_D)) { - this->correct_GC = false; + this->correct_D = false; parser.incr_input_error(); parser. - error_msg("Expected boolean value for correct_GC.", + error_msg("Expected boolean value for correct_D.", PHRQ_io::OT_CONTINUE); } - //correct_GC_defined = true; + //correct_D_defined = true; break; } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) @@ -498,11 +500,11 @@ cxxSurface::read_raw(CParser & parser, bool check) error_msg("Only_counter_ions not defined for SURFACE_RAW input.", PHRQ_io::OT_CONTINUE); } - //if (correct_GC_defined == false) + //if (correct_D_defined == false) //{ // parser.incr_input_error(); // parser. - // error_msg("correct_GC not defined for SURFACE_RAW input.", + // error_msg("correct_D not defined for SURFACE_RAW input.", // PHRQ_io::OT_CONTINUE); //} if (thickness_defined == false) @@ -582,7 +584,7 @@ cxxSurface::add(const cxxSurface & addee_in, LDBLE extensive) if (this->surface_comps.size() == 0) { this->only_counter_ions = addee.only_counter_ions; - this->correct_GC = addee.correct_GC; + this->correct_D = addee.correct_D; this->dl_type = addee.dl_type; this->type = addee.type; this->sites_units = addee.sites_units; @@ -754,7 +756,7 @@ cxxSurface::Serialize(Dictionary & dictionary, std::vector < int >&ints, doubles.push_back(this->debye_lengths); doubles.push_back(this->DDL_viscosity); doubles.push_back(this->DDL_limit); - ints.push_back(this->correct_GC ? 1 : 0); + ints.push_back(this->correct_D ? 1 : 0); ints.push_back(this->transport ? 1 : 0); this->totals.Serialize(dictionary, ints, doubles); ints.push_back(this->solution_equilibria ? 1 : 0); @@ -801,7 +803,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->debye_lengths = doubles[dd++]; this->DDL_viscosity = doubles[dd++]; this->DDL_limit = doubles[dd++]; - this->correct_GC = (ints[ii++] != 0); + this->correct_D = (ints[ii++] != 0); this->transport = (ints[ii++] != 0); this->totals.Deserialize(dictionary, ints, doubles, ii, dd); this->solution_equilibria = (ints[ii++] != 0); @@ -830,6 +832,6 @@ const std::vector< std::string >::value_type temp_vopts[] = { std::vector< std::string >::value_type("n_solution"), // 16 std::vector< std::string >::value_type("totals"), // 17 std::vector< std::string >::value_type("tidied"), // 18 - std::vector< std::string >::value_type("correct_gc") // 19 + std::vector< std::string >::value_type("correct_d") // 19 }; const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/Surface.h b/Surface.h index f01c716c4..b9ba60dd0 100644 --- a/Surface.h +++ b/Surface.h @@ -65,10 +65,12 @@ class cxxSurface:public cxxNumKeyword void Set_debye_lengths(LDBLE t) {debye_lengths = t;} LDBLE Get_DDL_viscosity(void) const {return DDL_viscosity;} void Set_DDL_viscosity(LDBLE t) {DDL_viscosity = t;} + void Calc_DDL_viscosity(bool tf) {calc_DDL_viscosity = tf;} + bool Get_calc_viscosity(void) const { return calc_DDL_viscosity; } LDBLE Get_DDL_limit(void) const {return DDL_limit;} void Set_DDL_limit(LDBLE t) {DDL_limit = t;} - bool Get_correct_GC(void) const { return correct_GC; } - void Set_correct_GC(bool tf) { correct_GC = tf; } + bool Get_correct_D(void) const { return correct_D; } + void Set_correct_D(bool tf) { correct_D = tf; } std::vector Donnan_factors; bool Get_transport(void) const {return transport;} void Set_transport(bool tf) {transport = tf;} @@ -93,8 +95,9 @@ class cxxSurface:public cxxNumKeyword LDBLE thickness; LDBLE debye_lengths; LDBLE DDL_viscosity; + bool calc_DDL_viscosity; LDBLE DDL_limit; - bool correct_GC; + bool correct_D; bool transport; cxxNameDouble totals; bool solution_equilibria; diff --git a/SurfaceCharge.cxx b/SurfaceCharge.cxx index 5d2857d85..226243692 100644 --- a/SurfaceCharge.cxx +++ b/SurfaceCharge.cxx @@ -36,6 +36,8 @@ PHRQ_base(io) grams = 0.0; charge_balance = 0.0; mass_water = 0.0; + DDL_viscosity = 0.0; + f_free = 0.0; la_psi = 0.0; capacitance[0] = 1.0; capacitance[1] = 5.0; @@ -68,6 +70,7 @@ cxxSurfaceCharge::dump_xml(std::ostream & s_oss, unsigned int indent) const charge_balance << "\"" << "\n"; s_oss << indent0 << "mass_water=\"" << this-> mass_water << "\"" << "\n"; + s_oss << indent0 << "f_free=\"" << this->f_free << "\"" << "\n"; s_oss << indent0 << "la_psi=\"" << this->la_psi << "\"" << "\n"; s_oss << indent0 << "capacitance=\"" << this-> capacitance[0] << " " << this->capacitance[0] << "\"" << "\n"; @@ -98,6 +101,8 @@ cxxSurfaceCharge::dump_raw(std::ostream & s_oss, unsigned int indent) const s_oss << indent0 << "-grams " << this->grams << "\n"; s_oss << indent0 << "-charge_balance " << this->charge_balance << "\n"; s_oss << indent0 << "-mass_water " << this->mass_water << "\n"; + s_oss << indent0 << "-f_free " << this->f_free << "\n"; + s_oss << indent0 << "-ddl_viscosity " << this->DDL_viscosity << "\n"; s_oss << indent0 << "-la_psi " << this->la_psi << "\n"; s_oss << indent0 << "-capacitance0 " << this->capacitance[0] << "\n"; s_oss << indent0 << "-capacitance1 " << this->capacitance[1] << "\n"; @@ -155,6 +160,7 @@ cxxSurfaceCharge::read_raw(CParser & parser, bool check) bool capacitance0_defined(false); bool capacitance1_defined(false); bool g_map_first(true); + bool DDL_viscosity_defined(false); for (;;) { @@ -225,7 +231,6 @@ cxxSurfaceCharge::read_raw(CParser & parser, bool check) mass_water_defined = true; break; - case 5: // la_psi if (!(parser.get_iss() >> this->la_psi)) { @@ -366,10 +371,27 @@ cxxSurfaceCharge::read_raw(CParser & parser, bool check) } } opt_save = 16; - - - break; + case 17: // f_free of water + if (!(parser.get_iss() >> this->f_free)) + { + this->f_free = 0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for f_free of mass_water.", + PHRQ_io::OT_CONTINUE); + } + break; + case 18: // DDL_viscosity + if (!(parser.get_iss() >> this->DDL_viscosity)) + { + this->DDL_viscosity = 1.0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for DDL_viscosity.", + PHRQ_io::OT_CONTINUE); + } + DDL_viscosity_defined = true; + break; + } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -454,9 +476,11 @@ cxxSurfaceCharge::add(const cxxSurfaceCharge & addee, LDBLE extensive) this->mass_water += addee.mass_water * extensive; this->la_psi = this->la_psi * f1 + addee.la_psi * f2; this->capacitance[0] = - this->capacitance[0] * f1 + this->capacitance[0] * f2; + this->capacitance[0] * f1 + addee.capacitance[0] * f2; this->capacitance[1] = - this->capacitance[1] * f1 + this->capacitance[1] * f2; + this->capacitance[1] * f1 + addee.capacitance[1] * f2; + this->f_free = this->f_free * f1 + addee.f_free * f2; + this->DDL_viscosity = this->DDL_viscosity * f1 + addee.DDL_viscosity * f2; this->diffuse_layer_totals.add_extensive(addee.diffuse_layer_totals, extensive); } @@ -486,6 +510,8 @@ cxxSurfaceCharge::Serialize(Dictionary & dictionary, std::vector < int >&ints, doubles.push_back(this->sigma1); doubles.push_back(this->sigma2); doubles.push_back(this->sigmaddl); + doubles.push_back(this->f_free); + doubles.push_back(this->DDL_viscosity); ints.push_back((int) this->g_map.size()); { std::map::iterator it; @@ -523,6 +549,8 @@ cxxSurfaceCharge::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->sigma1 = doubles[dd++]; this->sigma2 = doubles[dd++]; this->sigmaddl = doubles[dd++]; + this->f_free = doubles[dd++]; + this->DDL_viscosity = doubles[dd++]; { this->g_map.clear(); int count = ints[ii++]; @@ -581,6 +609,9 @@ const std::vector< std::string >::value_type temp_vopts[] = { std::vector< std::string >::value_type("sigma2"), // 13 std::vector< std::string >::value_type("sigmaddl"), // 14 std::vector< std::string >::value_type("g_map"), // 15 - std::vector< std::string >::value_type("diffuse_layer_species") // 16 + std::vector< std::string >::value_type("diffuse_layer_species"),// 16 + std::vector< std::string >::value_type("f_free"), // 17 + std::vector< std::string >::value_type("ddl_viscosity") // 18 + }; const std::vector< std::string > cxxSurfaceCharge::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/SurfaceCharge.h b/SurfaceCharge.h index 091b26b76..2333a331b 100644 --- a/SurfaceCharge.h +++ b/SurfaceCharge.h @@ -87,6 +87,10 @@ class cxxSurfaceCharge: public PHRQ_base void Set_charge_balance(LDBLE d) {this->charge_balance = d;} LDBLE Get_mass_water() const {return this->mass_water;} void Set_mass_water(LDBLE d) {this->mass_water = d;} + LDBLE Get_DDL_viscosity(void) const { return DDL_viscosity; } + void Set_DDL_viscosity(LDBLE t) { DDL_viscosity = t; } + LDBLE Get_f_free() const {return this->f_free;} + void Set_f_free(LDBLE d) {this->f_free = d;} LDBLE Get_la_psi() const {return this->la_psi;} void Set_la_psi(LDBLE d) {this->la_psi = d;} LDBLE Get_capacitance0() const {return this->capacitance[0];} @@ -117,6 +121,8 @@ class cxxSurfaceCharge: public PHRQ_base LDBLE grams; LDBLE charge_balance; LDBLE mass_water; + LDBLE DDL_viscosity; + LDBLE f_free; LDBLE la_psi; LDBLE capacitance[2]; cxxNameDouble diffuse_layer_totals; diff --git a/basicsubs.cpp b/basicsubs.cpp index bf545b4fb..9f4116b80 100644 --- a/basicsubs.cpp +++ b/basicsubs.cpp @@ -22,20 +22,20 @@ static char THIS_FILE[] = __FILE__; /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -activity(const char* species_name) +activity(const char *species_name) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; + class species *s_ptr; LDBLE a; s_ptr = s_search(species_name); if (s_ptr == s_h2o) { - a = pow((LDBLE)10., s_h2o->la); + a = pow((LDBLE) 10., s_h2o->la); } else if (s_ptr == s_eminus) { - a = pow((LDBLE)10., s_eminus->la); + a = pow((LDBLE) 10., s_eminus->la); } else if (s_ptr == NULL || s_ptr->in == FALSE) { @@ -43,17 +43,17 @@ activity(const char* species_name) } else { - a = pow((LDBLE)10., s_ptr->lm + s_ptr->lg); + a = pow((LDBLE) 10., s_ptr->lm + s_ptr->lg); } return (a); } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -activity_coefficient(const char* species_name) +activity_coefficient(const char *species_name) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; + class species *s_ptr; LDBLE g, dum = 0.0; s_ptr = s_search(species_name); @@ -61,7 +61,7 @@ activity_coefficient(const char* species_name) { if (s_ptr->type == EX && s_ptr->equiv && s_ptr->alk) dum = log10(s_ptr->equiv / s_ptr->alk); - g = pow((LDBLE)10., s_ptr->lg - dum); + g = pow((LDBLE) 10., s_ptr->lg - dum); } else { @@ -72,10 +72,10 @@ activity_coefficient(const char* species_name) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -log_activity_coefficient(const char* species_name) +log_activity_coefficient(const char *species_name) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; + class species *s_ptr; LDBLE g, dum = 0.0; s_ptr = s_search(species_name); @@ -94,10 +94,10 @@ log_activity_coefficient(const char* species_name) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -aqueous_vm(const char* species_name) +aqueous_vm(const char *species_name) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; + class species *s_ptr; LDBLE g; s_ptr = s_search(species_name); @@ -113,10 +113,10 @@ aqueous_vm(const char* species_name) } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -phase_vm(const char* phase_name) +phase_vm(const char *phase_name) /* ---------------------------------------------------------------------- */ { - class phase* phase_ptr; + class phase *phase_ptr; int l; LDBLE g; @@ -136,7 +136,7 @@ LDBLE Phreeqc:: sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw) /* ---------------------------------------------------------------------- */ { - if (sa_type == 0) + if (sa_type == 0) { // surface-area-calculation-Fixed_Surface return Sa; @@ -147,14 +147,14 @@ sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw double mass0 = m0 * gfw; double V0 = mass0 / d; double St0 = mass0 * Sa; // total surface - double a0 = pow(V0, 1.0 / 3.0); // side length - double Sp0 = 6.0 * a0 * a0; // surface particle + double a0 = pow(V0, 1.0/3.0); // side length + double Sp0 = 6.0 * a0*a0; // surface particle double np = St0 / Sp0; // number of particles - double RATS = Sa / St0; + double RATS = Sa / St0; double mass = m * gfw; double V = mass / d; - double a = pow(V, 1.0 / 3.0); - double St = 6.0 * a * a * np; + double a = pow(V, 1.0/3.0); + double St = 6.0 * a*a*np; return St * RATS; // total current surface } else if (sa_type == 2) @@ -163,47 +163,34 @@ sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw double mass0 = m0 * gfw; double V0 = mass0 / d; // volume double St0 = mass0 * Sa; // total surface - double a0 = pow(3.0 * V0 / (4.0 * pi), 1.0 / 3.0); // ((3*V0)/(4 * 3.14159265359))^(1/3) + double a0 = pow(3.0 * V0/(4.0 * pi), 1.0/3.0); // ((3*V0)/(4 * 3.14159265359))^(1/3) double Sp0 = (4.0 * pi) * a0 * a0; // surface particle double np = St0 / Sp0; // number of particles double RATS = Sa / St0; - + double mass = m * gfw; double V = mass / d; - double a = pow(3.0 * V / (4.0 * pi), 1.0 / 3.0); //((3*V)/(4 * 3.14159265359))^(1/3) + double a = pow(3.0 * V/(4.0 * pi), 1.0/3.0); //((3*V)/(4 * 3.14159265359))^(1/3) double St = 4.0 * pi * a * a * np; return St * RATS; // total current surface } - error_string = sformatf("Unknown surface area type in SA_DECLERCQ %d.", (int)sa_type); + error_string = sformatf( "Unknown surface area type in SA_DECLERCQ %d.", (int) sa_type); error_msg(error_string, CONTINUE); input_error++; return (MISSING); - + } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -diff_c(const char* species_name) +diff_c(const char *species_name) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; - LDBLE ka, l_z, Dw, ff, sqrt_mu; + class species *s_ptr; + LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av; sqrt_mu = sqrt(mu_x); s_ptr = s_search(species_name); - //LDBLE g; - //if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/) - //{ - // g = s_ptr->dw; - // if (s_ptr->dw_t) - // g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - // g *= viscos_0_25 / viscos * tk_x / 298.15; - //} - //else - //{ - // g = 0; - //} - //return (g); if (s_ptr == NULL) return(0); if ((Dw = s_ptr->dw) == 0) @@ -213,179 +200,99 @@ diff_c(const char* species_name) else return(0); } - if ((l_z = fabs(s_ptr->z)) == 0) + if (correct_Dw) { - //l_z = 1; // only a 1st approximation for correct_Dw in electrical field - } - else - { - if (s_ptr->dw_a2) - ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); - else - ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); - if (s_ptr->dw_a) + if ((l_z = fabs(s_ptr->z)) == 0) { - ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); - //if (print_viscosity && s_ptr->dw_a_visc) - // ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc); + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field } else { - ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + if (print_viscosity) + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1); + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : s_ptr->dw_a_visc ? 1 : pow(mu_x, 0.75)); + } + else + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = 1.0; + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : pow(mu_x, 0.75)); + } + ka = DH_B * a2 * sqrt_mu / (1 + a3); + ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka)); + Dw *= ff; } - Dw *= ff; } - if (tk_x != 298.15 && s_ptr->dw_t) Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - s_ptr->dw_corr = Dw; - return (Dw * viscos_0_25 / viscos_0); + s_ptr->dw_corr = Dw * viscos_0_25 / viscos_0; + return s_ptr->dw_corr; } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -setdiff_c(const char* species_name, double d) +setdiff_c(const char *species_name, double d) /* ---------------------------------------------------------------------- */ { - class species* s_ptr; - LDBLE ka, l_z, Dw, ff, sqrt_mu; + class species *s_ptr; + LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av; sqrt_mu = sqrt(mu_x); s_ptr = s_search(species_name); - - //LDBLE g; - //s_ptr = s_search(species_name); - //if (s_ptr != NULL) - //{ - - // s_ptr->dw = d; - // g = s_ptr->dw; - // if (s_ptr->dw_t) - // g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - // g *= viscos_0_25 / viscos * tk_x / 298.15;; - //} - //else - //{ - // g = 0; - //} - //return (g); if (s_ptr == NULL) return(0); Dw = s_ptr->dw = d; - if ((l_z = fabs(s_ptr->z)) == 0) + if (correct_Dw) { - //l_z = 1; // only a 1st approximation for correct_Dw in electrical field - } - else - { - if (s_ptr->dw_a2) - ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); - else - ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); - if (s_ptr->dw_a) + if ((l_z = fabs(s_ptr->z)) == 0) { - ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); - //if (print_viscosity && s_ptr->dw_a_visc) - // ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc); + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field } else { - ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + if (print_viscosity) + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1); + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : 1); + } + else + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = 1.0; + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : s_ptr->dw_a_visc ? 1 : pow(mu_x, 0.75)); + } + ka = DH_B * a2 * sqrt_mu / (1 + a3); + ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka)); + Dw *= ff; } - Dw *= ff; } - if (tk_x != 298.15 && s_ptr->dw_t) - Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - s_ptr->dw_corr = Dw; - return (Dw * viscos_0_25 / viscos_0); + s_ptr->dw_corr = Dw * viscos_0_25 / viscos_0; + return s_ptr->dw_corr; } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: calc_SC(void) /* ---------------------------------------------------------------------- */ { - //int i; - //LDBLE lm, a, l_z, Dw, SC, ff; - - //SC = 0; -# ifdef SKIP - for (i = 0; i < count_species_list; i++) - { - if (species_list[i].s->type == EX) - continue; - if (species_list[i].s->type == SURF) - continue; - if (i > 0 - && strcmp(species_list[i].s->name, - species_list[i - 1].s->name) == 0) - continue; - if (species_list[i].s == s_h2o) - continue; - if ((Dw = species_list[i].s->dw) == 0) - continue; - if ((l_z = fabs(species_list[i].s->z)) == 0) - continue; - - lm = species_list[i].s->lm; - if (lm > -9) - { - ff = (mu_x < .36 * l_z ? 0.6 / sqrt(l_z) : sqrt(mu_x) / l_z); - - ff *= species_list[i].s->lg; - if (ff > 0) ff = 0; - a = under(lm + ff); - if (species_list[i].s->dw_t) - Dw *= exp(species_list[i].s->dw_t / tk_x - species_list[i].s->dw_t / 298.15); // the viscosity multiplier is done in SC - SC += a * l_z * l_z * Dw; - } - } - SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298160.0); - /* correct for temperature dependency... - SC_T = SC_298 * (Dw_T / T) * (298 / Dw_298) and - Dw_T = Dw_298 * (T / 298) * (viscos_298 / viscos_T) give: - SC_T = SC_298 * (viscos_298 / viscos_T) - */ - SC *= viscos_0_25 / viscos; - - return (SC); - //# endif - for (i = 0; i < (int)this->s_x.size(); i++) - { - if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) - continue; - if ((Dw = s_x[i]->dw) == 0) - continue; - if ((l_z = fabs(s_x[i]->z)) == 0) - continue; - - lm = s_x[i]->lm; - if (lm > -9) - { - ff = (mu_x < .36 * l_z ? 0.6 / sqrt(l_z) : sqrt(mu_x) / l_z); - - ff *= s_x[i]->lg; - if (ff > 0) ff = 0; - a = under(lm + ff); - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity multiplier is done in SC - SC += a * l_z * l_z * Dw; - } - } - SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298160.0); - /* correct for temperature dependency... - SC_T = SC_298 * (Dw_T / T) * (298 / Dw_298) and - Dw_T = Dw_298 * (T / 298) * (viscos_298 / viscos_T) give: - SC_T = SC_298 * (viscos_298 / viscos_T) - */ - SC *= viscos_0_25 / viscos; - - return (SC); -# endif + class species *s_ptr; int i; - LDBLE ka, l_z, Dw, ff, sqrt_mu; + LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av, v_Cl = 1; sqrt_mu = sqrt(mu_x); + bool Falk = false; + s_ptr = s_search("H+"); + if (s_ptr == NULL) + return(0); + else if (s_ptr->dw_a3 > 24) Falk = true; SC = 0; //LDBLE ta1, ta2, ta3, ta4; @@ -401,147 +308,252 @@ calc_SC(void) // break; // } //} - for (i = 0; i < (int)this->s_x.size(); i++) + av = 0; + if (!Falk) { - if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) - continue; - if ((Dw = s_x[i]->dw) == 0) + for (i = 0; i < (int)this->s_x.size(); i++) { - if (correct_Dw) - Dw = default_Dw; - else + if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) continue; - } - if (s_x[i]->lm < min_dif_LM) - continue; - if ((l_z = fabs(s_x[i]->z)) == 0) - { - //l_z = 1; // only a 1st approximation for correct_Dw in electrical field - } - else - { - if (s_x[i]->dw_a2) - ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); - else + if ((Dw = s_x[i]->dw) == 0) { - ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); - //ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2)); - //ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2); + if (correct_Dw) + Dw = default_Dw; + else + continue; } - if (s_x[i]->dw_a) + if (s_x[i]->lm < min_dif_LM) + continue; + if (tk_x != 298.15) + { + if (s_x[i]->dw_t) + Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); + //else + //{ + // Dw *= exp(ta1 / tk_x - ta1 / 298.15); + //} + } + // correct for temperature dependent viscosity of pure water... + Dw *= viscos_0_25 / viscos_0; + s_x[i]->dw_corr = Dw; + if ((l_z = fabs(s_x[i]->z)) == 0) { - ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); - //if (print_viscosity && s_x[i]->dw_a_visc) - // ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc); + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field + continue; } else { - ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); - //ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka)); + //if (s_x[i]->dw_a2) + // ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.259)); + // //ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); + //else + //{ + // ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.259)); + // //ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2)); + // //ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2); + //} + //if (s_x[i]->dw_a) + //{ + // ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); + // if (print_viscosity && s_x[i]->dw_a_visc) + // ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc); + //} + //else + //{ + // ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + // //ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka)); + //} + //Dw *= ff; + s_ptr = s_x[i]; + if (print_viscosity) + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1); + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : s_ptr->dw_a_visc ? 1 : pow(mu_x, 0.75)); + } + else + { + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + av = 1.0; + a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : pow(mu_x, 0.75)); + } + ka = DH_B * a2 * sqrt_mu / (1 + a3); + ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka)); } + Dw *= ff; + s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw; + SC += s_x[i]->dw_t_SC; } - if (tk_x != 298.15) - { - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); - //else - //{ - // Dw *= exp(ta1 / tk_x - ta1 / 298.15); - //} - } - s_x[i]->dw_corr = Dw; - - s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw; - SC += s_x[i]->dw_t_SC; - } - SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0); - /* correct for viscosity dependency... - SC_T = SC_298 * (viscos_298 / viscos_T) - */ - SC *= viscos_0_25 / viscos_0; - return (SC); -} -#ifdef SKIP -/*Debye-Onsager according to Robinson and Stokes, 1954, JACS 75, 1991, - but with sqrt charge multiplier for B2 and mu^ff dependent ka */ -LDBLE q, B1, B2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, Sum_m_dw, z_plus, z_min, t1, t2, Dw_SC; - -m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = Sum_m_dw = z_plus = z_min = 0; -SC = 0; -for (i = 0; i < (int)this->s_x.size(); i++) -{ - if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) - continue; - if ((l_z = s_x[i]->z) == 0) - continue; - if ((lm = s_x[i]->lm) < -9) - continue; - if ((Dw = s_x[i]->dw) == 0) - Dw = 1e-9; - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity multiplier cancels in q... - if (l_z > 0) - { - m_plus += s_x[i]->moles; - t1 = s_x[i]->moles * l_z; - eq_plus += t1; - eq_dw_plus += t1 * Dw; - Sum_m_dw += s_x[i]->moles * Dw; + SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0); + return (SC); } else { - m_min += s_x[i]->moles; - t1 = s_x[i]->moles * l_z; - eq_min -= t1; - eq_dw_min -= t1 * Dw; - Sum_m_dw += s_x[i]->moles * Dw; - } -} -// Falkenhagen, q = (Sum(z1 * m1*Dw1) + Sum(z2 *m2*Dw2)) / ((Sum(m1*Dw1) + Sum(m2*Dw2))(av_z1 + av_z2)) -z_plus = eq_plus / m_plus; // |av_z1| -z_min = eq_min / m_min; // |av_z2| -q = (eq_dw_plus + eq_dw_min) / (Sum_m_dw * (z_min + z_plus)); -t1 = 1.60218e-19 * 1.60218e-19 / (6 * pi); -B1 = t1 / (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) -t2 = viscos_0; // (1 - 0.5) * viscos_0 + 0.5 * viscos; -B2 = t1 * AVOGADRO / t2 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4) - -Dw_SC = viscos_0_25 / t2 * 1e4 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298160.0); -for (i = 0; i < (int)this->s_x.size(); i++) -{ - if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) - continue; - if ((l_z = fabs(s_x[i]->z)) == 0) - continue; - if ((lm = s_x[i]->lm) < -9) - continue; - if ((Dw = s_x[i]->dw) == 0) - Dw = 1e-9; - - Dw *= Dw_SC; - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity factor is in Dw_SC - a = (s_x[i]->dw_a ? s_x[i]->dw_a : 3.5); - ka = DH_B * a; - if (s_x[i]->dw_a2) - ka *= pow((double)mu_x, s_x[i]->dw_a2); - else - ka *= sqrt_mu; + /* the phreeqc equation from Appelo, 2017, CCR 101, 102 with viscosity correction, e.g. for SO4-2 and its complexes: + Dw dw_t a a2 visc -5< a3 <5 + -dw 1.07e-9 236 0.7281 3.452 -0.1515 -3.043 # obsolete + or + Debye-Onsager with (1 + ka)^2 in the denominator, + for the individual ions according to their contribution to mu, with sqrt charge multiplier for B2 and + a in ka corrected by volume (or mu^a2, if a3 = -10), and * (viscos_0 / viscos)^av + Dw dw_t a a2 visc a3 = (0) or >5 + -dw 1.03e-9 -14 4.03 0.8341 1.679 # Li+, ka = DH_B * a * (1 + (vm - v0))^a2 * mu^0.5 + */ + LDBLE q, sqrt_q, B1, B2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, z_plus, z_min, t1, Dw_SC; - // Falkenhagen... - //SC += under(lm) * l_z * l_z * (Dw - B2 * l_z * sqrt_mu / (1 + ka)) * (1 - B1 * sqrt_mu / - // ((1 + ka) * (1 + ka * sqrt(q) + ka * ka / 6))); + m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = z_plus = z_min = 0; + SC = 0; + LDBLE eps_c = eps_r; // Cl concentration corrected eps_r + // average z and Dw for transport numbers t1_0 and t2_0 at zero conc's, and q of the solution... + for (i = 0; i < (int)this->s_x.size(); i++) + { + if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) + continue; + if (s_x[i]->lm < min_dif_LM) + continue; + if ((Dw = s_x[i]->dw) == 0) + { + if (correct_Dw) Dw = default_Dw; // or charge based...Dw = l_z > 0 ? 1.6e-9 / l_z : 2e-9 / -l_z; + else continue; + } + if (tk_x != 298.15 && s_x[i]->dw_t) + Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); + Dw *= viscos_0_25 / viscos_0; + s_x[i]->dw_corr = Dw; + + if ((l_z = s_x[i]->z) == 0) + continue; + + if (l_z > 0) + { + m_plus += s_x[i]->moles; + a = s_x[i]->moles * l_z; + eq_plus += a; + eq_dw_plus += a * Dw; + } + else + { + m_min += s_x[i]->moles; + a = s_x[i]->moles * l_z; + eq_min -= a; + eq_dw_min -= a * Dw; + } + } + // q = z1 * z2 / ((z2 * t1_0 + z1 * t2_0) * (z1 + z2)); + // z1 = z_plus, z2 = z_min, t1_0 = (eq_dw_plus / m_plus) / (eq_dw_plus / m_plus + eq_dw_min / m_min) + //a = (eq_plus - eq_min) / 2; + //eq_min += a; + //eq_plus -= a; + z_plus = eq_plus / m_plus; // |av_z1| + z_min = eq_min / m_min; // |av_z2| + t1 = (eq_dw_plus / m_plus) / (eq_dw_plus / m_plus + eq_dw_min / m_min); + 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)) + 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) + //B1 = a / (2 * 8.8542e-12 * eps_c * 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 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4) + + LDBLE mu_plus, mu_min, lz, ll_SC, v0; + // the + and - contributions to mu, assuming -2, -1, 1, 2 charge numbers... + mu_min = 3 * m_min * (z_min - 1) + m_min; + mu_plus = 3 * m_plus * (z_plus - 1) + m_plus; + + Dw_SC = 1e4 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298.15e3); // for recalculating Dw to ll0 + t1 = calc_solution_volume(); + ll_SC = 0.5e3 * (eq_plus + eq_min) / t1 * mass_water_aq_x / t1; // recalculates ll to SC in uS/cm, with mu in mol/kgw + + for (i = 0; i < (int)this->s_x.size(); i++) + { + if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) + continue; + if ((lz = s_x[i]->z) == 0) + continue; + if (s_x[i]->lm < min_dif_LM) + continue; + if ((Dw = s_x[i]->dw_corr) == 0) + continue; + l_z = fabs(lz); + a3 = s_x[i]->dw_a3; + if (a3 != 0 && a3 > -5.01 && a3 < 4.99) + { + // with the phreeqc equation... + s_ptr = s_x[i]; + a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6); + a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73); + a3 = pow(mu_x, a3); + if (print_viscosity) + av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1); + else + av = 1.0; + ka = DH_B * a2 * sqrt_mu / (1 + a3); + ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka)); - t1 = (Dw - (B1 * Dw + B2) * sqrt_mu / (1 + ka)); - //t1 = (Dw - (B1 * Dw + B2 * sqrt(l_z)) * sqrt_mu / (1 + ka)); - //t1 = (Dw - (B1 * Dw + B2 * l_z * l_z) * sqrt_mu / (1 + ka)); - if (t1 > 0) - SC += under(lm) * l_z * l_z * t1; + Dw *= ff; + s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw; + SC += s_x[i]->dw_t_SC * 1e3 * Dw_SC; + } + else + { + // Falkenhagen... + if (!a3) a3 = 10; + a = (s_x[i]->dw_a ? s_x[i]->dw_a : 3.5); + a2 = (s_x[i]->dw_a2);// ? s_x[i]->dw_a2 : 0.5); + av = (s_x[i]->dw_a_visc);// ? s_x[i]->dw_a_visc : 0.8); + if (lz < -0.5/* && lz > -1.5*/) + { + // Mg and Ca give different numbers than H+, Li+, Na+ and K for anions... + t1 = (z_plus - 1); + //a -= a3 / 1000 * l_z * t1; + //a2 += a3 / 1000 * l_z * t1; + a -= 0.126 * l_z * t1; + a2 += 0.126 * l_z * t1; + //av += 0 * t1; + } + Dw *= Dw_SC * l_z; + if (!a2 || !strcmp(s_x[i]->name, "H+")) + t1 = 1; + else + { + v0 = calc_vm0(s_x[i]->name, tc_x, 1, 0); + t1 = 1 + (s_x[i]->rxn_x.logk[vm_tc] - v0); + if (a2 && t1 > 0) + t1 = pow(t1, a2); + //t1 = (t1 <= a3 / 10 ? a3 / 10 : t1); + t1 = (t1 < 1 ? 1 : t1); + } + if (a3 >= 5) + // with the vm correction of a in ka.. + ka = DH_B * a * t1 * sqrt_mu; + else + // with a mu^dw_a2 correction of a.. + ka = DH_B * a * pow((double)mu_x, a2); + + t1 = (Dw - B2 * l_z * sqrt_mu / (1 + ka)) * + (1 - B1 * sqrt_mu / ((1 + ka) * (1 + ka)));// +ka * ka / 6))); // S.cm2/eq / (kgw/L) + //t1 = (Dw - B2 * l_z * sqrt_mu / (1 + ka)) * + // (1 - B1 * sqrt_mu / ((1 + ka) *(1 + ka * sqrt_q + ka * ka / 6))); // S.cm2/eq / (kgw/L) + if (av) + t1 *= pow(viscos_0 / viscos, av); + + // fractional contribution in mu, and correct for charge imbalance + a2 = 2 / (eq_plus + eq_min); + a = (lz > 0 ? mu_plus / (eq_plus * a2) : mu_min / (eq_min * a2)); + t1 *= s_x[i]->moles * l_z * l_z / a; + t1 *= ll_SC; + s_x[i]->dw_t_SC = t1 / (1e3 * Dw_SC); + SC += t1; + } + } + return SC; + } } -return (SC * 1e3); -#endif /* VP: Density Start */ /* ---------------------------------------------------------------------- */ @@ -622,11 +634,12 @@ calc_dens(void) V_solutes += s_x[i]->moles * s_x[i]->logk[vm_tc]; } /* If pure water then return rho_0 */ - if (M_T == 0) - return rho_0; - else - return rho_0 * (1e3 + M_T / mass_water_aq_x) / (rho_0 * V_solutes / mass_water_aq_x + 1e3); - + density_x = rho_0; + if (M_T > 0.0) + { + density_x = rho_0 * (1e3 + M_T / mass_water_aq_x) / (rho_0 * V_solutes / mass_water_aq_x + 1e3); + } + return density_x; //M_T /= 1e3; //solution_mass = mass_water_aq_x + M_T; //V_solutes = M_T - rho_0 * V_solutes / 1e3; @@ -769,7 +782,7 @@ calc_logk_s(const char* name) class species* s_ptr; LDBLE lk, l_logk[MAX_LOG_K_INDICES]; - Utilities::strcpy_safe(token, MAX_LENGTH, name); + Utilities::strcpy_safe(token, MAX_LENGTH, name); s_ptr = s_search(token); if (s_ptr != NULL) { @@ -1246,7 +1259,7 @@ calc_t_sc(const char* name) calc_SC(); if (!SC) return (0); - LDBLE t = s_ptr->dw_t_SC * 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0) * viscos_0_25 / viscos_0; + LDBLE t = s_ptr->dw_t_SC * 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0); return (t / SC); } return (0); @@ -1387,7 +1400,7 @@ equi_phase_delta(const char* phase_name) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -equivalent_fraction(const char* name, LDBLE * eq, std::string & elt_name) +equivalent_fraction(const char* name, LDBLE* eq, std::string& elt_name) /* ---------------------------------------------------------------------- */ { class species* s_ptr = s_search(name); @@ -1888,7 +1901,7 @@ saturation_ratio(const char* phase_name) /* ---------------------------------------------------------------------- */ int Phreeqc:: -saturation_index(const char* phase_name, LDBLE * iap, LDBLE * si) +saturation_index(const char* phase_name, LDBLE* iap, LDBLE* si) /* ---------------------------------------------------------------------- */ { class rxn_token* rxn_ptr; @@ -2062,7 +2075,7 @@ sum_match_ss(const char* mytemplate, const char* name) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -list_ss(std::string ss_name, cxxNameDouble & composition) +list_ss(std::string ss_name, cxxNameDouble& composition) /* ---------------------------------------------------------------------- */ { LDBLE tot = 0; @@ -2711,7 +2724,7 @@ total_mole(const char* total_name) /* ---------------------------------------------------------------------- */ int Phreeqc:: -get_edl_species(cxxSurfaceCharge & charge_ref) +get_edl_species(cxxSurfaceCharge& charge_ref) /* ---------------------------------------------------------------------- */ { @@ -2749,7 +2762,7 @@ get_edl_species(cxxSurfaceCharge & charge_ref) } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -edl_species(const char* surf_name, LDBLE * count, char*** names, LDBLE * *moles, LDBLE * area, LDBLE * thickness) +edl_species(const char* surf_name, LDBLE* count, char*** names, LDBLE** moles, LDBLE* area, LDBLE* thickness) /* ---------------------------------------------------------------------- */ { /* @@ -2806,8 +2819,8 @@ edl_species(const char* surf_name, LDBLE * count, char*** names, LDBLE * *moles, } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -system_total(const char* total_name, LDBLE * count, char*** names, - char*** types, LDBLE * *moles, int isort) +system_total(const char* total_name, LDBLE* count, char*** names, + char*** types, LDBLE** moles, int isort) /* ---------------------------------------------------------------------- */ { /* @@ -2922,7 +2935,7 @@ system_total(const char* total_name, LDBLE * count, char*** names, /* ---------------------------------------------------------------------- */ std::string Phreeqc:: -kinetics_formula(std::string kin_name, cxxNameDouble & stoichiometry) +kinetics_formula(std::string kin_name, cxxNameDouble& stoichiometry) /* ---------------------------------------------------------------------- */ { /* @@ -2973,7 +2986,7 @@ kinetics_formula(std::string kin_name, cxxNameDouble & stoichiometry) } /* ---------------------------------------------------------------------- */ std::string Phreeqc:: -phase_formula(std::string phase_name, cxxNameDouble & stoichiometry) +phase_formula(std::string phase_name, cxxNameDouble& stoichiometry) /* ---------------------------------------------------------------------- */ { /* @@ -2996,7 +3009,7 @@ phase_formula(std::string phase_name, cxxNameDouble & stoichiometry) } /* ---------------------------------------------------------------------- */ std::string Phreeqc:: -species_formula(std::string phase_name, cxxNameDouble & stoichiometry) +species_formula(std::string phase_name, cxxNameDouble& stoichiometry) /* ---------------------------------------------------------------------- */ { /* @@ -3186,13 +3199,13 @@ int Phreeqc:: system_total_aq(void) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ for (i = 0; i < (int)this->s_x.size(); i++) { //if (s_x[i]->type != AQ) @@ -3213,13 +3226,13 @@ int Phreeqc:: system_total_ex(void) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ for (i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type != EX) @@ -3241,13 +3254,13 @@ int Phreeqc:: system_total_surf(void) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ for (i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type != SURF) @@ -3266,20 +3279,20 @@ int Phreeqc:: system_total_gas(void) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i; - /* - * find total in gas phase - */ +/* + * find total in gas phase + */ if (use.Get_gas_phase_ptr() == NULL) return (OK); - cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr(); for (size_t j = 0; j < gas_phase_ptr->Get_gas_comps().size(); j++) { - class phase* phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j].Get_phase_name().c_str(), + class phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j].Get_phase_name().c_str(), &i, FALSE); assert(phase_ptr); size_t count_sys = sys.size(); @@ -3296,24 +3309,24 @@ int Phreeqc:: system_total_equi(void) /* ---------------------------------------------------------------------- */ { - /* - * Equilibrium phases - */ +/* + * Equilibrium phases + */ if (use.Get_pp_assemblage_ptr() == NULL) return (OK); std::map comps = use.Get_pp_assemblage_ptr()->Get_pp_assemblage_comps(); std::map ::iterator it = comps.begin(); - for (; it != comps.end(); it++) + for ( ; it != comps.end(); it++) { - cxxPPassemblageComp* comp_ptr = &(it->second); - int l; - class phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); - size_t count_sys = sys.size(); - sys.resize(count_sys + 1); - sys[count_sys].name = string_duplicate(phase_ptr->name); - sys[count_sys].moles = equi_phase(sys[count_sys].name); - sys_tot += sys[count_sys].moles; - sys[count_sys].type = string_duplicate("equi"); + cxxPPassemblageComp *comp_ptr = &(it->second); + int l; + class phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); + size_t count_sys = sys.size(); + sys.resize(count_sys + 1); + sys[count_sys].name = string_duplicate(phase_ptr->name); + sys[count_sys].moles = equi_phase(sys[count_sys].name); + sys_tot += sys[count_sys].moles; + sys[count_sys].type = string_duplicate("equi"); } return (OK); } @@ -3322,21 +3335,21 @@ int Phreeqc:: system_total_kin(void) /* ---------------------------------------------------------------------- */ { - /* - * Equilibrium phases - */ +/* + * Equilibrium phases + */ if (use.Get_kinetics_ptr() == NULL) return (OK); std::vector comps = use.Get_kinetics_ptr()->Get_kinetics_comps(); - for (size_t i = 0; i < comps.size(); i++) + for (size_t i=0 ; i < comps.size(); i++) { - cxxKineticsComp* comp_ptr = &comps[i]; - size_t count_sys = sys.size(); - sys.resize(count_sys + 1); - sys[count_sys].name = string_duplicate(comp_ptr->Get_rate_name().c_str()); - sys[count_sys].moles = comp_ptr->Get_m(); - sys_tot += sys[count_sys].moles; - sys[count_sys].type = string_duplicate("kin"); + cxxKineticsComp *comp_ptr = &comps[i]; + size_t count_sys = sys.size(); + sys.resize(count_sys + 1); + sys[count_sys].name = string_duplicate(comp_ptr->Get_rate_name().c_str()); + sys[count_sys].moles = comp_ptr->Get_m(); + sys_tot += sys[count_sys].moles; + sys[count_sys].type = string_duplicate("kin"); } return (OK); } @@ -3345,24 +3358,24 @@ int Phreeqc:: system_total_ss(void) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ - /* - * Solid solutions - */ +/* + * Solid solutions + */ if (use.Get_ss_assemblage_ptr() == NULL) return (OK); - std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); + std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); for (size_t k = 0; k < ss_ptrs.size(); k++) { - cxxSS* ss_ptr = ss_ptrs[k]; + cxxSS *ss_ptr = ss_ptrs[k]; for (size_t i = 0; i < ss_ptr->Get_ss_comps().size(); i++) { - cxxSScomp* comp_ptr = &(ss_ptr->Get_ss_comps()[i]); + cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[i]); int l; - class phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); + class phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); size_t count_sys = sys.size(); sys.resize(count_sys + 1); sys[count_sys].name = string_duplicate(phase_ptr->name); @@ -3375,19 +3388,19 @@ system_total_ss(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: -system_total_elt(const char* total_name) +system_total_elt(const char *total_name) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i, j, k; LDBLE molality, moles_excess, moles_surface, mass_water_surface; char name[MAX_LENGTH]; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ for (i = 0; i < (int)this->s_x.size(); i++) { count_elts = 0; @@ -3452,7 +3465,7 @@ system_total_elt(const char* total_name) { if (x[k]->type != SURFACE_CB) continue; - cxxSurfaceCharge* charge_ptr = use.Get_surface_ptr()->Find_charge(x[k]->surface_charge); + cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[k]->surface_charge); i++; /* * Loop through all surface components, calculate each H2O surface (diffuse layer), @@ -3469,8 +3482,8 @@ system_total_elt(const char* total_name) moles_excess = mass_water_aq_x * molality * (charge_ptr->Get_g_map()[s_x[j]->z].Get_g() * s_x[j]->erm_ddl + - mass_water_surface / mass_water_aq_x * (s_x[j]->erm_ddl - - 1)); + mass_water_surface / mass_water_aq_x * (s_x[j]->erm_ddl - + 1)); moles_surface = mass_water_surface * molality + moles_excess; /* * Accumulate elements in diffuse layer @@ -3498,9 +3511,9 @@ system_total_elt(const char* total_name) } } } - /* - * find total moles in mineral phases - */ +/* + * find total moles in mineral phases + */ if (use.Get_pp_assemblage_in() == TRUE && use.Get_pp_assemblage_ptr() != NULL) { for (i = 0; i < count_unknowns; i++) @@ -3509,7 +3522,7 @@ system_total_elt(const char* total_name) continue; //std::map::iterator it; //it = pp_assemblage_ptr->Get_pp_assemblage_comps().find(x[i]->pp_assemblage_comp_name); - cxxPPassemblageComp* comp_ptr = (cxxPPassemblageComp*)x[i]->pp_assemblage_comp_ptr; + cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[i]->pp_assemblage_comp_ptr; //if (it->second.Get_add_formula().size() > 0) if (comp_ptr->Get_add_formula().size() > 0) continue; @@ -3517,7 +3530,7 @@ system_total_elt(const char* total_name) paren_count = 0; int j; //class phase * phase_ptr = phase_bsearch(x[i]->pp_assemblage_comp_name, &j, FALSE); - class phase* phase_ptr = x[i]->phase; + class phase * phase_ptr = x[i]->phase; add_elt_list(phase_ptr->next_elt, x[i]->moles); elt_list_combine(); for (j = 0; j < count_elts; j++) @@ -3535,26 +3548,26 @@ system_total_elt(const char* total_name) } } } - /* - * Solid solutions - */ +/* + * Solid solutions + */ if (use.Get_ss_assemblage_ptr() != NULL) { - std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); + std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); for (size_t k = 0; k < ss_ptrs.size(); k++) { - cxxSS* ss_ptr = ss_ptrs[k]; + cxxSS *ss_ptr = ss_ptrs[k]; if (ss_ptr->Get_ss_in()) { for (size_t i = 0; i < ss_ptr->Get_ss_comps().size(); i++) { - cxxSScomp* comp_ptr = &(ss_ptr->Get_ss_comps()[i]); + cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[i]); int l; - class phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); + class phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); count_elts = 0; paren_count = 0; add_elt_list(phase_ptr->next_elt, - comp_ptr->Get_moles()); + comp_ptr->Get_moles()); elt_list_combine(); for (j = 0; j < count_elts; j++) { @@ -3574,15 +3587,15 @@ system_total_elt(const char* total_name) } } } - /* - * find total in gas phase - */ +/* + * find total in gas phase + */ if (use.Get_gas_phase_ptr() != NULL) { - cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr(); for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++) { - class phase* phase_ptr = + class phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[i].Get_phase_name().c_str(), &k, FALSE); assert(phase_ptr); if (phase_ptr->in == TRUE) @@ -3615,19 +3628,19 @@ system_total_elt(const char* total_name) /* ---------------------------------------------------------------------- */ int Phreeqc:: -system_total_elt_secondary(const char* total_name) +system_total_elt_secondary(const char *total_name) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i, j, k, l; LDBLE molality, moles_excess, moles_surface, mass_water_surface, sum, coef; char name[MAX_LENGTH]; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ for (i = 0; i < (int)this->s_x.size(); i++) { count_elts = 0; @@ -3693,7 +3706,7 @@ system_total_elt_secondary(const char* total_name) { if (x[k]->type != SURFACE_CB) continue; - cxxSurfaceCharge* charge_ptr = use.Get_surface_ptr()->Find_charge(x[k]->surface_charge); + cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[k]->surface_charge); i++; /* * Loop through all surface components, calculate each H2O surface (diffuse layer), @@ -3744,9 +3757,9 @@ system_total_elt_secondary(const char* total_name) } } } - /* - * find total moles in mineral phases - */ +/* + * find total moles in mineral phases + */ if (use.Get_pp_assemblage_in() == TRUE && use.Get_pp_assemblage_ptr() != NULL) { for (i = 0; i < count_unknowns; i++) @@ -3755,7 +3768,7 @@ system_total_elt_secondary(const char* total_name) continue; //std::map::iterator it; //it = pp_assemblage_ptr->Get_pp_assemblage_comps().find(x[i]->pp_assemblage_comp_name); - cxxPPassemblageComp* comp_ptr = (cxxPPassemblageComp*)x[i]->pp_assemblage_comp_ptr; + cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[i]->pp_assemblage_comp_ptr; //if (it->second.Get_add_formula().size() > 0) if (comp_ptr->Get_add_formula().size() > 0) continue; @@ -3763,8 +3776,8 @@ system_total_elt_secondary(const char* total_name) paren_count = 0; int j; //class phase * phase_ptr = phase_bsearch(x[i]->pp_assemblage_comp_name, &j, FALSE); - class phase* phase_ptr = x[i]->phase; - add_elt_list(phase_ptr->next_sys_total, x[i]->moles); + class phase * phase_ptr = x[i]->phase; + add_elt_list(phase_ptr->next_sys_total, x[i]->moles); elt_list_combine(); for (j = 0; j < count_elts; j++) { @@ -3782,26 +3795,26 @@ system_total_elt_secondary(const char* total_name) } } } - /* - * Solid solutions - */ +/* + * Solid solutions + */ if (use.Get_ss_assemblage_ptr() != NULL) { - std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); + std::vector ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize(); for (size_t i = 0; i < ss_ptrs.size(); i++) { - cxxSS* ss_ptr = ss_ptrs[i]; + cxxSS *ss_ptr = ss_ptrs[i]; if (ss_ptr->Get_ss_in()) { for (size_t k = 0; k < ss_ptr->Get_ss_comps().size(); k++) { - cxxSScomp* comp_ptr = &(ss_ptr->Get_ss_comps()[k]); + cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[k]); int l; - class phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); + class phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); count_elts = 0; paren_count = 0; add_elt_list(phase_ptr->next_sys_total, - comp_ptr->Get_moles()); + comp_ptr->Get_moles()); elt_list_combine(); for (j = 0; j < count_elts; j++) { @@ -3821,15 +3834,15 @@ system_total_elt_secondary(const char* total_name) } } } - /* - * find total in gas phase - */ +/* + * find total in gas phase + */ if (use.Get_gas_phase_ptr() != NULL) { - cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); - for (size_t j = 0; j < gas_phase_ptr->Get_gas_comps().size(); j++) + cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr(); + for (size_t j = 0; j < gas_phase_ptr->Get_gas_comps().size(); j++) { - class phase* phase_ptr = + class phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j].Get_phase_name().c_str(), &i, FALSE); assert(phase_ptr); if (phase_ptr->in == TRUE) @@ -3837,13 +3850,13 @@ system_total_elt_secondary(const char* total_name) count_elts = 0; paren_count = 0; add_elt_list(phase_ptr->next_sys_total, - phase_ptr->moles_x); + phase_ptr->moles_x); elt_list_combine(); /* * Look for element */ - for (size_t j1 = 0; j1 < (size_t)count_elts; j1++) + for (size_t j1 = 0; j1 < (size_t) count_elts; j1++) { if (strcmp(elt_list[j1].elt->name, total_name) == 0) { @@ -3867,7 +3880,7 @@ int Phreeqc:: solution_number(void) /* ---------------------------------------------------------------------- */ { - Phreeqc* PhreeqcPtr = this; + Phreeqc * PhreeqcPtr = this; int soln_no = -999; if (PhreeqcPtr->state == TRANSPORT) { @@ -3900,17 +3913,17 @@ solution_number(void) } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -solution_sum_secondary(const char* total_name) +solution_sum_secondary(const char *total_name) /* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in system and lists of species/phases in sort order - */ +/* + * Provides total moles in system and lists of species/phases in sort order + */ int i, j; LDBLE sum; - /* - * find total moles in aq, surface, and exchange - */ +/* + * find total moles in aq, surface, and exchange + */ sum = 0; for (i = 0; i < (int)this->s_x.size(); i++) { @@ -3944,13 +3957,13 @@ solution_sum_secondary(const char* total_name) /* ---------------------------------------------------------------------- */ int Phreeqc:: -system_species_compare(const void* ptr1, const void* ptr2) +system_species_compare(const void *ptr1, const void *ptr2) /* ---------------------------------------------------------------------- */ { - const class system_species* a, * b; + const class system_species *a, *b; - a = (const class system_species*)ptr1; - b = (const class system_species*)ptr2; + a = (const class system_species *) ptr1; + b = (const class system_species *) ptr2; if (a->moles < b->moles) return (1); if (a->moles > b->moles) @@ -3970,21 +3983,21 @@ system_species_compare_name(const void* ptr1, const void* ptr2) /* ---------------------------------------------------------------------- */ int Phreeqc:: -system_total_solids(cxxExchange * exchange_ptr, - cxxPPassemblage * pp_assemblage_ptr, - cxxGasPhase * gas_phase_ptr, - cxxSSassemblage * ss_assemblage_ptr, - cxxSurface * surface_ptr) - /* ---------------------------------------------------------------------- */ +system_total_solids(cxxExchange *exchange_ptr, + cxxPPassemblage *pp_assemblage_ptr, + cxxGasPhase *gas_phase_ptr, + cxxSSassemblage *ss_assemblage_ptr, + cxxSurface *surface_ptr) +/* ---------------------------------------------------------------------- */ { - /* - * Provides total moles in solid phases - */ +/* + * Provides total moles in solid phases + */ count_elts = 0; paren_count = 0; - /* - * find total moles in exchanger - */ +/* + * find total moles in exchanger + */ if (exchange_ptr != NULL) { for (size_t i = 0; i < exchange_ptr->Get_exchange_comps().size(); i++) @@ -4001,17 +4014,17 @@ system_total_solids(cxxExchange * exchange_ptr, } if (ss_assemblage_ptr != NULL) { - std::vector ss_ptrs = ss_assemblage_ptr->Vectorize(); + std::vector ss_ptrs = ss_assemblage_ptr->Vectorize(); for (size_t i = 0; i < ss_ptrs.size(); i++) { - cxxSS* ss_ptr = ss_ptrs[i]; + cxxSS *ss_ptr = ss_ptrs[i]; for (size_t j = 0; j < ss_ptr->Get_ss_comps().size(); j++) { - cxxSScomp* comp_ptr = &(ss_ptr->Get_ss_comps()[j]); + cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[j]); int l; - class phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); + class phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE); add_elt_list(phase_ptr->next_elt, - comp_ptr->Get_moles()); + comp_ptr->Get_moles()); } } } @@ -4020,7 +4033,7 @@ system_total_solids(cxxExchange * exchange_ptr, for (size_t j = 0; j < gas_phase_ptr->Get_gas_comps().size(); j++) { int i; - class phase* phase_ptr = + class phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j].Get_phase_name().c_str(), &i, FALSE); add_elt_list(phase_ptr->next_elt, gas_phase_ptr->Get_gas_comps()[j].Get_moles()); } @@ -4028,13 +4041,13 @@ system_total_solids(cxxExchange * exchange_ptr, if (pp_assemblage_ptr != NULL) { std::map::iterator it; - it = pp_assemblage_ptr->Get_pp_assemblage_comps().begin(); - for (; it != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); it++) + it = pp_assemblage_ptr->Get_pp_assemblage_comps().begin(); + for ( ; it != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); it++) { int j; - class phase* phase_ptr = phase_bsearch(it->first.c_str(), &j, FALSE); + class phase * phase_ptr = phase_bsearch(it->first.c_str(), &j, FALSE); add_elt_list(phase_ptr->next_elt, - it->second.Get_moles()); + it->second.Get_moles()); } } elt_list_combine(); @@ -4042,14 +4055,14 @@ system_total_solids(cxxExchange * exchange_ptr, } LDBLE Phreeqc:: -iso_value(const char* total_name) +iso_value(const char *total_name) { int j; char token[MAX_LENGTH]; char my_total_name[MAX_LENGTH]; Utilities::strcpy_safe(token, MAX_LENGTH, ""); Utilities::strcpy_safe(my_total_name, MAX_LENGTH, total_name); - while (replace(" ", "_", my_total_name)); + while (replace(" ","_",my_total_name)); for (j = 0; j < (int)isotope_ratio.size(); j++) { if (isotope_ratio[j]->ratio == MISSING) @@ -4059,8 +4072,8 @@ iso_value(const char* total_name) return (isotope_ratio[j]->converted_ratio); } Utilities::strcpy_safe(my_total_name, MAX_LENGTH, total_name); - while (replace("[", "", my_total_name)); - while (replace("]", "", my_total_name)); + while (replace("[","",my_total_name)); + while (replace("]","",my_total_name)); Utilities::strcat_safe(token, MAX_LENGTH, "R("); Utilities::strcat_safe(token, MAX_LENGTH, my_total_name); Utilities::strcat_safe(token, MAX_LENGTH, ")"); @@ -4075,16 +4088,16 @@ iso_value(const char* total_name) return -1000.; } -char* Phreeqc:: -iso_unit(const char* total_name) +char * Phreeqc:: +iso_unit(const char *total_name) { int j; char token[MAX_LENGTH], unit[MAX_LENGTH]; - class master_isotope* master_isotope_ptr; + class master_isotope *master_isotope_ptr; char my_total_name[MAX_LENGTH]; Utilities::strcpy_safe(token, MAX_LENGTH, ""); Utilities::strcpy_safe(my_total_name, MAX_LENGTH, total_name); - while (replace(" ", "_", my_total_name)); + while (replace(" ","_",my_total_name)); Utilities::strcpy_safe(unit, MAX_LENGTH, "unknown"); for (j = 0; j < (int)isotope_ratio.size(); j++) { @@ -4100,8 +4113,8 @@ iso_unit(const char* total_name) return string_duplicate(unit); } Utilities::strcpy_safe(my_total_name, MAX_LENGTH, total_name); - while (replace("[", "", my_total_name)); - while (replace("]", "", my_total_name)); + while (replace("[","",my_total_name)); + while (replace("]","",my_total_name)); Utilities::strcat_safe(token, MAX_LENGTH, "R("); Utilities::strcat_safe(token, MAX_LENGTH, my_total_name); Utilities::strcat_safe(token, MAX_LENGTH, ")"); @@ -4122,13 +4135,13 @@ iso_unit(const char* total_name) } int Phreeqc:: -basic_compile(const char* commands, void** lnbase, void** vbase, void** lpbase) +basic_compile(const char *commands, void **lnbase, void **vbase, void **lpbase) { return this->basic_interpreter->basic_compile(commands, lnbase, vbase, lpbase); } int Phreeqc:: -basic_run(char* commands, void* lnbase, void* vbase, void* lpbase) +basic_run(char *commands, void *lnbase, void *vbase, void *lpbase) { return this->basic_interpreter->basic_run(commands, lnbase, vbase, lpbase); } @@ -4145,7 +4158,7 @@ basic_free(void) #include "BasicCallback.h" double Phreeqc:: -basic_callback(double x1, double x2, const char* str) +basic_callback(double x1, double x2, const char * str) { if (this->basicCallback) { @@ -4158,10 +4171,10 @@ basic_callback(double x1, double x2, const char* str) #ifdef IPHREEQC_NO_FORTRAN_MODULE double Phreeqc:: -basic_callback(double x1, double x2, char* str) +basic_callback(double x1, double x2, char * str) #else double Phreeqc:: -basic_callback(double x1, double x2, const char* str) +basic_callback(double x1, double x2, const char * str) #endif { double local_x1 = x1; @@ -4169,35 +4182,35 @@ basic_callback(double x1, double x2, const char* str) if (basic_callback_ptr != NULL) { - return (*basic_callback_ptr) (x1, x2, (const char*)str, basic_callback_cookie); + return (*basic_callback_ptr) (x1, x2, (const char *) str, basic_callback_cookie); } if (basic_fortran_callback_ptr != NULL) { #ifdef IPHREEQC_NO_FORTRAN_MODULE - return (*basic_fortran_callback_ptr) (&local_x1, &local_x2, str, (int)strlen(str)); + return (*basic_fortran_callback_ptr) (&local_x1, &local_x2, str, (int) strlen(str)); #else - return (*basic_fortran_callback_ptr) (&local_x1, &local_x2, str, (int)strlen(str)); + return (*basic_fortran_callback_ptr) (&local_x1, &local_x2, str, (int) strlen(str)); #endif } return 0; } -void -Phreeqc::register_basic_callback(double (*fcn)(double x1, double x2, const char* str, void* cookie), void* cookie1) +void +Phreeqc::register_basic_callback(double (*fcn)(double x1, double x2, const char *str, void *cookie), void *cookie1) { this->basic_callback_ptr = fcn; this->basic_callback_cookie = cookie1; } #ifdef IPHREEQC_NO_FORTRAN_MODULE -void -Phreeqc::register_fortran_basic_callback(double (*fcn)(double* x1, double* x2, char* str, size_t l)) +void +Phreeqc::register_fortran_basic_callback(double ( *fcn)(double *x1, double *x2, char *str, size_t l)) { this->basic_fortran_callback_ptr = fcn; } #else -void -Phreeqc::register_fortran_basic_callback(double (*fcn)(double* x1, double* x2, const char* str, int l)) +void +Phreeqc::register_fortran_basic_callback(double ( *fcn)(double *x1, double *x2, const char *str, int l)) { this->basic_fortran_callback_ptr = fcn; } diff --git a/class_main.cpp b/class_main.cpp index 03fb3e401..731f02274 100644 --- a/class_main.cpp +++ b/class_main.cpp @@ -299,7 +299,7 @@ write_banner(void) /* version */ #ifdef NPP - len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "3.7.1"); + len = sprintf(buffer, "* PHREEQC-%s *", "3.7.3"); #else len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "@VERSION@"); #endif @@ -323,7 +323,7 @@ write_banner(void) /* date */ #ifdef NPP - len = snprintf(buffer, sizeof(buffer), "%s", "March 20, 2023"); + len = sprintf(buffer, "%s", "March 14, 2024, with bug-fixes and new items"); #else len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@"); #endif @@ -507,7 +507,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, output_msg(sformatf(" Input file: %s\n", in_file.c_str())); output_msg(sformatf(" Output file: %s\n", out_file.c_str())); #ifdef NPP - output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 20, 2023\n")); + output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 14, 2024, with bug-fixes and new items\n")); #endif output_msg(sformatf("Database file: %s\n\n", token.c_str())); #ifdef NPP diff --git a/gases.cpp b/gases.cpp index 9f71089f4..bfcc64d74 100644 --- a/gases.cpp +++ b/gases.cpp @@ -606,11 +606,10 @@ calc_PR(void) { 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 < -3 ? -3 : phi)); phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); } else - phi = -3.0; // fugacity coefficient = 0.05 + 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 // **** diff --git a/global_structures.h b/global_structures.h index c4672cfe0..e9d82f9d3 100644 --- a/global_structures.h +++ b/global_structures.h @@ -716,9 +716,10 @@ class species dw = 0; // correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15) dw_t = 0; - // parms for calc'ng SC = SC0 * exp(-dw_a * z * mu^0.5 / (1 + DH_B * dw_a2 * mu^0.5)) + // parms for calc'ng SC = SC0 * exp(-dw_a * z * mu^0.5 / (1 + DH_B * dw_a2 * mu^0.5) / (1 + mu^dw_a3)) dw_a = 0; dw_a2 = 0; + dw_a3 = 0; dw_a_visc = 0; // viscosity correction of SC dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC dw_t_visc = 0; // contribution to viscosity @@ -781,6 +782,7 @@ class species LDBLE dw_t; LDBLE dw_a; LDBLE dw_a2; + LDBLE dw_a3; LDBLE dw_a_visc; LDBLE dw_t_SC; LDBLE dw_t_visc; diff --git a/integrate.cpp b/integrate.cpp index 8f20cd277..f21236961 100644 --- a/integrate.cpp +++ b/integrate.cpp @@ -733,58 +733,78 @@ calc_all_donnan(void) { bool converge; int cd_m; - LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot, co_ion; - LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1; + LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_surf_aq, co_ion; + LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1, viscos; + cxxSurface *surf_ptr = use.Get_surface_ptr(); - if (use.Get_surface_ptr() == NULL) + if (surf_ptr == NULL) return (OK); f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * tk_x * mu_x); - bool only_count = use.Get_surface_ptr()->Get_only_counter_ions(); - bool correct_GC = use.Get_surface_ptr()->Get_correct_GC(); + bool only_count = surf_ptr->Get_only_counter_ions(); + bool correct_D = surf_ptr->Get_correct_D(); /* calculate g for each surface... */ - if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths() || - correct_GC) // DL_pitz && correct_GC + if (!calculating_deriv || surf_ptr->Get_debye_lengths() || + correct_D) // DL_pitz && correct_D initial_surface_water(); - LDBLE nDbl = 1; - if (correct_GC) + // z1, z2, fr_cat2 are the counter-ions, z_1, z_2, fr_ani2 are for co-ions. + LDBLE nDbl = 1, db_lim = 2, f_free, fr_cat2, fr_ani2; + LDBLE z1, z2, z_1, z_2; + z1 = z2 = z_1 = z_2 = f_free = fr_cat2 = fr_ani2 = 0; + /* + * sum eq of each charge number in solution... + */ + std::map::iterator it; + for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) + { + it->second = 0.0; + } + charge_group_map.clear(); + for (int i = 0; i < (int)this->s_x.size(); i++) + { + if (s_x[i]->type > HPLUS) + continue; + charge_group_map[s_x[i]->z] += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl; + } + for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) + { + if (it->first < -1.5) { z_2 += it->second; continue; } + else if (it->first < 0) { z_1 += it->second; continue; } + else if (it->first < 1.5) { z1 += it->second; continue; } + else { z2 += it->second; continue; } + } + if (correct_D) { - if ((nDbl = use.Get_surface_ptr()->Get_debye_lengths()) == 0) + if ((nDbl = surf_ptr->Get_debye_lengths()) == 0) { LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3); - nDbl = use.Get_surface_ptr()->Get_thickness() / debye_length; - //use.Get_surface_ptr()->Set_debye_lengths(nDbl); + nDbl = surf_ptr->Get_thickness() / debye_length; + } + fr_ani2 = z_2 / (z_1 + z_2); + fr_cat2 = z2 / (z1 + z2); + db_lim = 2 - 0.5 * (fr_cat2 + fr_ani2); + if (nDbl > db_lim) + { + f_free = 1 - db_lim / nDbl; + if (f_free < 0) f_free = 0; } } converge = TRUE; + viscos = 0; for (int j = 0; j < count_unknowns; j++) { if (x[j]->type != SURFACE_CB) continue; - cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); + cxxSurfaceCharge *charge_ptr = surf_ptr->Find_charge(x[j]->surface_charge); if (debug_diffuse_layer == TRUE) output_msg(sformatf("Calc_all_g, X[%d]\n", j)); - /* - * sum eq of each charge number in solution... - */ - std::map::iterator it; - for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) - { - it->second = 0.0; - } - charge_group_map.clear(); - for (int i = 0; i < (int)this->s_x.size(); i++) - { - if (s_x[i]->type > HPLUS) - continue; - charge_group_map[s_x[i]->z] += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl; - } + /* find surface charge from potential... */ A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams(); - if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) + if (surf_ptr->Get_type() == cxxSurface::CD_MUSIC) { f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */ f_psi = f_psi / 2; @@ -797,7 +817,7 @@ calc_all_donnan(void) } surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL; LDBLE lim_seq = 5e3; - if (correct_GC) lim_seq = 5e1; + if (correct_D) lim_seq = 5e3; if (fabs(surf_chrg_eq) > lim_seq) { surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq); @@ -816,23 +836,24 @@ calc_all_donnan(void) std::vector zcorr(charge_group_map.size()); std::vector zcorr2(charge_group_map.size()); //LDBLE fD = 0; - psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, zcorr); - psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, zcorr2); + psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, f_free, zcorr); + psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, f_free, zcorr2); /*output_msg(sformatf( "psi's %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */ /* fill in g's */ ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; - ratio_aq_tot = charge_ptr->Get_mass_water() / mass_water_bulk_x; - + ratio_surf_aq = charge_ptr->Get_mass_water() / mass_water_surfaces_x; + //ratio_surf_aq = charge_ptr->Get_mass_water() / mass_water_bulk_x; + if (correct_D) + ratio_aq *= (1 - f_free); int z_iter = 0; for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) { LDBLE z = it->first, z1 = z; co_ion = surf_chrg_eq * z; - if (correct_GC) + if (correct_D) z1 = zcorr[z_iter]; - //z1 *= cgc[0] * pow(z_factor, abs(z)); if (!ratio_aq) { @@ -878,18 +899,18 @@ calc_all_donnan(void) /* save Boltzmann factor * water fraction for MCD calc's in transport */ if (converge) { - if (only_count) - { - if (co_ion > 0) // co-ions are not in the DL + if (only_count && co_ion > 0) // co-ions are not in the DL charge_ptr->Get_z_gMCD_map()[z] = 0; - else // assume that counter-ions have the free water conc for diffusion - charge_ptr->Get_z_gMCD_map()[z] = ratio_aq_tot; - } else - charge_ptr->Get_z_gMCD_map()[z] = (new_g / ratio_aq + 1) * ratio_aq_tot; + { + charge_ptr->Get_z_gMCD_map()[z] = (exp(cd_m * z1 * psi_avg) * (1 - f_free) + f_free) * + ratio_surf_aq;// * s_x[]->moles == mol_DL in charge_ptr + } } z_iter++; } + + charge_ptr->Set_f_free(f_free); if (debug_diffuse_layer == TRUE) { std::string name = x[j]->master[0]->elt->name; @@ -920,8 +941,9 @@ calc_init_donnan(void) /* ---------------------------------------------------------------------- */ { LDBLE f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq; + cxxSurface *surf_ptr = use.Get_surface_ptr(); - if (use.Get_surface_ptr() == NULL) + if (surf_ptr == NULL) return (OK); f_sinh = //sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * @@ -963,12 +985,12 @@ calc_init_donnan(void) { if (x[j]->type != SURFACE_CB) continue; - cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); + cxxSurfaceCharge *charge_ptr = surf_ptr->Find_charge(x[j]->surface_charge); charge_ptr->Get_g_map().clear(); /* find surface charge from potential... */ A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams(); - if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) + if (surf_ptr->Get_type() == cxxSurface::CD_MUSIC) { f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */ f_psi = f_psi / 2; @@ -978,7 +1000,7 @@ calc_init_donnan(void) surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL; /* find psi_avg that matches surface charge... */ - psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, zcorr); + psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, 0, zcorr); /* fill in g's */ ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; @@ -991,7 +1013,7 @@ calc_init_donnan(void) charge_ptr->Get_g_map()[z].Set_g(ratio_aq * (exp(-z * psi_avg) - 1)); - if (use.Get_surface_ptr()->Get_only_counter_ions() + if (surf_ptr->Get_only_counter_ions() && ((surf_chrg_eq < 0 && z < 0) || (surf_chrg_eq > 0 && z > 0))) charge_ptr->Get_g_map()[z].Set_g(-ratio_aq); @@ -1038,45 +1060,84 @@ calc_init_donnan(void) } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector &zcorr) +calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, LDBLE f_free, std::vector &zcorr) /* ---------------------------------------------------------------------- */ { /* * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge */ - LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, z1, z1_c, eq, co_ion, sum_counter, sum_co; - + LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, Z1, Z1_c, eq, co_ion, sum_counter, sum_co; + LDBLE z1, z2, z_1, z_2; ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; p = 0; if (surf_chrg_eq == 0 || ratio_aq == 0) return (0.0); else if (surf_chrg_eq < 0) - p = -0.5 * log(-surf_chrg_eq * ratio_aq / mu_x + 1); + p = -0.5 * log(-surf_chrg_eq * ratio_aq * (1 - f_free) / mu_x + 1); else if (surf_chrg_eq > 0) - p = 0.5 * log(surf_chrg_eq * ratio_aq / mu_x + 1); + p = 0.5 * log(surf_chrg_eq * ratio_aq * (1 - f_free) / mu_x + 1); /* * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq - * g(p) = exp(-p * z_i) * ratio_aq + * g(p) = exp(-p * z_i) * ratio_aq * (1 - f_free) * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0: * g(p) = (exp(-p *z_i) - 1) * ratio_aq - * with correct_GC true: - * correct ions to better match the integrated concentrations: - z == 1? z *= 0.285 cgc[6] - z == 2? z *= 0.372 cgc[7] - z == -1? z *= cgc[0] * (mu_x**( cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4]) - z == -2? z *= cgc[0] * (mu_x**(cgc[5] * cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4]) + * with correct_D true and f_free > 0: + c_edl = c_free * (f_free + (1 - f_free) * exp(-p * z_i)) + * with correct_D true and f_free == 0: + * correct ions to better match the integrated PB concentrations: + Gamma = abs(surf_chrg_eq / A_surf / 1e-6) + a = cgc[1] * nDbl**cgc[2] + b = Gamma**cgc[3] / abs(log10(I)) + counter_ions... + z == 1? z1 = cgc[0] * I**(a * b) + z == 2? z2 = 2 * cgc[4] * I**(cgc[5] * a * b) + co_ions... + c = cgc[7] * nDbl**cgc[8] * Gamma**cgc[9] + z == -1? z_1 = -cgc[6] * I**(c) + z == -2? z_2 = -2 * cgc[10] * I**(c * cgc[11]) + c_edl = c_free * exp(-p * z_i) */ - cxxSurface *surf_p = use.Get_surface_ptr(); - bool correct_GC = surf_p->Get_correct_GC(), local_correct_GC = correct_GC; - bool only_count = surf_p->Get_only_counter_ions(); - LDBLE Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6, - cgc[10] = { 0.36, 0.1721, 0.798, 0.287, 0.1457, 1.2, 0.285, 0.372 }; - - if (!surf_p->Donnan_factors.empty()) - std::copy(surf_p->Donnan_factors.begin(), surf_p->Donnan_factors.end(), cgc); - - cgc[1] *= pow(nDbl, cgc[2]) * pow(Gamma, cgc[3]) * pow(mu_x, cgc[4]); + cxxSurface *surf_ptr = use.Get_surface_ptr(); + bool correct_D = surf_ptr->Get_correct_D(), local_correct_D = correct_D; + bool only_count = surf_ptr->Get_only_counter_ions(); + LDBLE Gamma, cgc[12] = { 0.3805, -0.0106, 1.96, 0.812, + 0.395, 2.13, + 0.380, 0.0408, 0.799, 0.594, + 0.373, 1.181 }; + if (correct_D) + { + if (f_free) + { + z1 = z2 = z_1 = z_2 = 1; + } + else + { + if (!surf_ptr->Donnan_factors.empty()) + { + std::copy(std::begin(surf_ptr->Donnan_factors), std::end(surf_ptr->Donnan_factors), cgc); + z1 = cgc[0]; + z2 = cgc[1]; + z_1 = cgc[2]; + z_2 = cgc[3]; + } + else + { + Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6; + LDBLE a = cgc[1] * pow(nDbl, cgc[2]), + b = pow(Gamma, cgc[3]) / abs(log10(mu_x)); + // counter_ions... + z1 = cgc[0] * pow(mu_x, (a * b)); + z2 = cgc[4] * pow(mu_x, (cgc[5] * a * b)); + if (z1 > 1) z1 = 1; + if (z2 > 1) z2 = 1; + // co_ions... + LDBLE c = cgc[7] * pow(nDbl, cgc[8]) * pow(Gamma, cgc[9]); + z_1 = cgc[6] * pow(mu_x, c); + z_2 = cgc[10] * pow(mu_x, (c * cgc[11])); + } + } + } int l_iter = 0, z_iter; sum_co = sum_counter = 0; @@ -1085,16 +1146,16 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std:: fd = surf_chrg_eq; fd1 = 0.0; z_iter = 0; - if (l_iter == 1 && local_correct_GC && fabs(sum_counter) < fabs(sum_co)) + if (l_iter == 1 && local_correct_D && fabs(sum_counter) < fabs(sum_co)) { - local_correct_GC = false; + local_correct_D = false; l_iter = 0; } std::map::iterator it; for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) { z = it->first; - z1 = z; + Z1 = z; if (l_iter == 0) zcorr[z_iter] = z; co_ion = surf_chrg_eq * z; if (!z || (only_count && co_ion > 0)) @@ -1102,30 +1163,30 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std:: z_iter++; continue; } - if (nDbl && local_correct_GC) + if (nDbl && local_correct_D) { /*psi_DL = fabs(p * p_psi);*/ if (co_ion < 0) {//counter-ion - if (fabs(z) > 1) temp = cgc[7]; - else temp = cgc[6]; + if (fabs(z) > 1.5) temp = z2; + else temp = z1; sum_counter += z * temp; } else {// co-ion - if (fabs(z) > 1) temp = cgc[0] * pow(mu_x, cgc[1] * cgc[5]); - else temp = cgc[0] * pow(mu_x, cgc[1]); + if (fabs(z) > 1.5) temp = z_2; + else temp = z_1; sum_co += z * temp; } zcorr[z_iter] = z * temp; } - z1 = zcorr[z_iter]; + Z1 = zcorr[z_iter]; eq = it->second; - temp = exp(-z1 * p) * ratio_aq; + temp = exp(-Z1 * p) * ratio_aq * (1 - f_free); fd += eq * temp; - fd1 -= z1 * eq * temp; - if (z == 1) z1_c = z1; + fd1 -= Z1 * eq * temp; + if (z == 1) Z1_c = Z1; z_iter++; } fd /= -fd1; @@ -1152,7 +1213,7 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std:: if (debug_diffuse_layer == TRUE) output_msg(sformatf( "iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n", - l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * z1_c * p_psi))); + l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * Z1_c * p_psi))); return (p); } diff --git a/kinetics.cpp b/kinetics.cpp index 0472265c1..d1158e2f5 100644 --- a/kinetics.cpp +++ b/kinetics.cpp @@ -1579,7 +1579,10 @@ set_and_run(int i, int use_mix, int use_kinetics, int nsaver, converge = model(); } sum_species(); - viscosity(); + viscos = viscosity(NULL); + use.Get_solution_ptr()->Set_viscosity(viscos); + if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL) + use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr())); return (converge); } diff --git a/mainsubs.cpp b/mainsubs.cpp index baa20e718..e57c6389c 100644 --- a/mainsubs.cpp +++ b/mainsubs.cpp @@ -434,7 +434,10 @@ initial_solutions(int print) diagonal_scale = (diag) ? TRUE : FALSE; converge1 = check_residuals(); sum_species(); - viscosity(); + viscos = viscosity(NULL); + use.Get_solution_ptr()->Set_viscosity(viscos); + if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL) + use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr())); add_isotopes(solution_ref); punch_all(); print_all(); @@ -537,7 +540,6 @@ initial_exchangers(int print) converge = model(); converge1 = check_residuals(); sum_species(); - viscosity(); species_list_sort(); print_exchange(); if (pr.user_print) @@ -1260,12 +1262,12 @@ xsolution_save(int n_user) temp_solution.Set_pe(solution_pe_x); temp_solution.Set_mu(mu_x); temp_solution.Set_ah2o(ah2o_x); - //temp_solution.Set_density(density_x); - temp_solution.Set_density(calc_dens()); - temp_solution.Set_viscosity(this->viscos); + // the subroutine is called at the start of a new simulation, and the following 2 go wrong since s_x is not updated + temp_solution.Set_density(density_x); + temp_solution.Set_viscosity(viscos); temp_solution.Set_total_h(total_h_x); temp_solution.Set_total_o(total_o_x); - temp_solution.Set_cb(cb_x); /* cb_x does not include surface charge sfter sum_species */ + temp_solution.Set_cb(cb_x); /* cb_x does not include surface charge after sum_species */ /* does include surface charge after step */ temp_solution.Set_mass_water(mass_water_aq_x); temp_solution.Set_total_alkalinity(total_alkalinity); diff --git a/model.cpp b/model.cpp index 60e077fb8..863ff8d05 100644 --- a/model.cpp +++ b/model.cpp @@ -4891,7 +4891,6 @@ sum_species(void) solution_pe_x = -s_eminus->la; ah2o_x = exp(s_h2o->la * LOG_10); - density_x = 1.0; if (s_o2 != NULL) s_o2->moles = under(s_o2->lm) * mass_water_aq_x; if (s_h2 != NULL) diff --git a/prep.cpp b/prep.cpp index a8349d3e1..e08cd5f8f 100644 --- a/prep.cpp +++ b/prep.cpp @@ -3996,17 +3996,16 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) { 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 < -3 ? -3 : phi)); phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); //if (phi > 4.44) // phi = 4.44; } else - phi = -3.0; // fugacity coefficient = 0.05 - //if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3) + phi = -4.6; // fugacity coefficient = 0.01 + //if (!strcmp(phase_ptr->name, "H2O(g)") && phi < -4.6) //{ //// avoid such phi... - // phi = -3; + // phi = -4.6; //} phase_ptr->pr_phi = exp(phi); phase_ptr->pr_si_f = phi / LOG_10; @@ -5404,6 +5403,53 @@ calc_vm(LDBLE tc, LDBLE pa) return OK; } +LDBLE Phreeqc::calc_vm0(const char * species_name, LDBLE tc, LDBLE pa, LDBLE mu) +{ + /* + * Calculate molar volume of an aqueous species at tc, pa and mu + */ + if (llnl_temp.size() > 0) return OK; + class species *s_ptr; + LDBLE g = 0; + s_ptr = s_search(species_name); + if (s_ptr == s_h2o) + return 18.016 / rho_0; + if (s_ptr != NULL && s_ptr->in != FALSE && s_ptr->type < EMINUS && s_ptr->logk[vma1]) + { + LDBLE pb_s = 2600. + pa * 1.01325, TK_s = tc + 45.15, sqrt_mu = sqrt(mu); + /* supcrt volume at I = 0... */ + g = s_ptr->logk[vma1] + s_ptr->logk[vma2] / pb_s + + (s_ptr->logk[vma3] + s_ptr->logk[vma4] / pb_s) / TK_s - + s_ptr->logk[wref] * QBrn; + if (s_ptr->z) + { + /* the ionic strength term * I^0.5... */ + if (s_ptr->logk[b_Av] < 1e-5) + g += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu; + else + { + /* limit the Debye-Hueckel slope by b... */ + /* pitzer... */ + //s_ptr->rxn_x.logk[vm_tc] += s_ptr->z * s_ptr->z * 0.5 * DH_Av * + // log(1 + s_ptr->logk[b_Av] * sqrt(mu_x)) / s_ptr->logk[b_Av]; + /* extended DH... */ + g += s_ptr->z * s_ptr->z * 0.5 * DH_Av * + sqrt_mu / (1 + s_ptr->logk[b_Av] * DH_B * sqrt_mu); + } + /* plus the volume terms * I... */ + if (s_ptr->logk[vmi1] != 0.0 || s_ptr->logk[vmi2] != 0.0 || s_ptr->logk[vmi3] != 0.0) + { + LDBLE bi = s_ptr->logk[vmi1] + s_ptr->logk[vmi2] / TK_s + s_ptr->logk[vmi3] * TK_s; + if (s_ptr->logk[vmi4] == 1.0) + g += bi * mu; + else + g += bi * pow(mu, s_ptr->logk[vmi4]); + } + } + } + return g; +} + /* ---------------------------------------------------------------------- */ int Phreeqc:: k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */ @@ -5419,7 +5465,7 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */ if (pa != current_pa) goto proceed; if (fabs(mu_x - current_mu) > 1e-3 * mu_x) goto proceed; if (mu_terms_in_logk) goto proceed; - return OK; + return OK; proceed: diff --git a/print.cpp b/print.cpp index e372db9bd..325e12430 100644 --- a/print.cpp +++ b/print.cpp @@ -270,6 +270,16 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr) output_msg(sformatf( "\tWater in diffuse layer: %8.3e kg, %4.1f%% of total DDL-water.\n", (double) charge_ptr->Get_mass_water(), (double) d)); + if (charge_ptr->Get_DDL_viscosity()) + { + if (d == 100) + output_msg(sformatf( + "\t\t viscosity: %7.5f mPa s.\n", (double)charge_ptr->Get_DDL_viscosity())); + else + output_msg(sformatf( + "\t\t viscosity: %7.5f mPa s for this DDL water. (%7.5f mPa s for total DDL-water.)\n", (double)charge_ptr->Get_DDL_viscosity(), (double)use.Get_surface_ptr()->Get_DDL_viscosity())); + } + if (use.Get_surface_ptr()->Get_debye_lengths() > 0 && d > 0) { sum_surfs = 0.0; @@ -279,8 +289,7 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr) continue; cxxSurfaceCharge * charge_ptr_search = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); sum_surfs += - charge_ptr_search->Get_specific_area() * - charge_ptr_search->Get_grams(); + charge_ptr_search->Get_specific_area() * charge_ptr_search->Get_grams(); } r = 0.002 * mass_water_bulk_x / sum_surfs; output_msg(sformatf( @@ -304,10 +313,8 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr) if (s_x[j]->type > HPLUS) continue; molality = under(s_x[j]->lm); - moles_excess = mass_water_aq_x * molality * (charge_ptr->Get_g_map()[s_x[j]->z].Get_g() * - s_x[j]->erm_ddl + - mass_water_surface / - mass_water_aq_x * (s_x[j]->erm_ddl - 1)); + moles_excess = mass_water_aq_x * molality * (charge_ptr->Get_g_map()[s_x[j]->z].Get_g() * s_x[j]->erm_ddl + + mass_water_surface / mass_water_aq_x * (s_x[j]->erm_ddl - 1)); moles_surface = mass_water_surface * molality + moles_excess; if (debug_diffuse_layer == TRUE) { @@ -336,17 +343,26 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr) } else { - LDBLE exp_g = charge_ptr->Get_g_map()[1].Get_g() * mass_water_aq_x / mass_water_surface + 1; + LDBLE exp_g = charge_ptr->Get_g_map()[1].Get_g() * mass_water_aq_x / ((1 - charge_ptr->Get_f_free()) * mass_water_surface) + 1; LDBLE psi_DL = -log(exp_g) * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ; - if (use.Get_surface_ptr()->Get_correct_GC()) + if (use.Get_surface_ptr()->Get_correct_D()) + { output_msg(sformatf( "\n\tTotal moles in diffuse layer (excluding water), Donnan corrected to match Poisson-Boltzmann.")); + output_msg(sformatf( + "\n\tDonnan Layer potential, psi_DL = %10.3e V, for (1 - f_free) of DL water = %10.3e kg (f_free = %5.3f).\n\tBoltzmann factor, exp(-psi_DL * z * z_corr * F / RT) = %9.3e (= c_DL / c_free if z is +1)", + psi_DL, (1 - charge_ptr->Get_f_free()) * mass_water_surface, charge_ptr->Get_f_free(), exp_g)); + output_msg(sformatf( + "\n\t\tThus: Moles of Na+ = (c_DL * (1 - f_free) + f_free) * c_free * kg DDL-water\n\n")); + } else + { output_msg(sformatf( "\n\tTotal moles in diffuse layer (excluding water), Donnan calculation.")); - output_msg(sformatf( - "\n\tDonnan Layer potential, psi_DL = %10.3e V.\n\tBoltzmann factor, exp(-psi_DL * F / RT) = %9.3e (= c_DL / c_free if z is +1).\n\n", - psi_DL, exp_g)); + output_msg(sformatf( + "\n\tDonnan Layer potential, psi_DL = %10.3e V.\n\tBoltzmann factor, exp(-psi_DL * F / RT) = %9.3e (= c_DL / c_free if z is +1).\n\n", + psi_DL, exp_g)); + } } output_msg(sformatf("\tElement \t Moles\n")); for (j = 0; j < count_elts; j++) diff --git a/read.cpp b/read.cpp index 2ac1bea75..a5eed844b 100644 --- a/read.cpp +++ b/read.cpp @@ -5499,9 +5499,17 @@ read_species(void) input_error++; break; } - s_ptr->dw_t = 0; s_ptr->dw_a = 0; s_ptr->dw_a2 = 0; s_ptr->dw_a_visc = 0; - i = sscanf(next_char, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT, &s_ptr->dw, &s_ptr->dw_t, - &s_ptr->dw_a, &s_ptr->dw_a2, &s_ptr->dw_a_visc); + s_ptr->dw_t = 0; s_ptr->dw_a = 0; s_ptr->dw_a2 = 0; s_ptr->dw_a3 = 0; s_ptr->dw_a_visc = 0; + i = sscanf(next_char, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT, + &s_ptr->dw, &s_ptr->dw_t, &s_ptr->dw_a, &s_ptr->dw_a2, &s_ptr->dw_a_visc, &s_ptr->dw_a3); + if (i < 1) + { + input_error++; + error_msg("Expecting numeric values for the diffusion coefficient, its temperature dependence, and coefficients for the SC calculation.", + CONTINUE); + return (ERROR); + } + s_ptr->dw_corr = s_ptr->dw; opt_save = OPTION_DEFAULT; break; @@ -6414,7 +6422,7 @@ read_surface(void) if (thickness != 0) { error_msg - ("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5 limit 0.9 correct_GC true\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8 correct_GC false", + ("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5 limit 0.9 correct_D true\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8 correct_D false", CONTINUE); error_msg(line_save, CONTINUE); input_error++; @@ -6442,12 +6450,12 @@ read_surface(void) copy_token(token1, &next_char); if (token1[0] == 'T' || token1[0] == 't' || token1[0] == 'F' || token1[0] == 'f') { - temp_surface.Set_correct_GC(get_true_false(token1.c_str(), TRUE) == TRUE); + temp_surface.Set_correct_D(get_true_false(token1.c_str(), TRUE) == TRUE); continue; } else { error_msg - ("Expected True or False for correct_GC (which brings co-ion concentrations closer to their integrated double layer value).", + ("Expected True or False for correct_D (which brings co-ion concentrations closer to their integrated double layer value).", CONTINUE); error_msg(line_save, CONTINUE); input_error++; @@ -6460,9 +6468,19 @@ read_surface(void) if (j == DIGIT) { (void)sscanf(token1.c_str(), SCANFORMAT, &dummy); + if(dummy == 0) + { + dummy = 1; temp_surface.Calc_DDL_viscosity(true); + } temp_surface.Set_DDL_viscosity(dummy); continue; } + else if (token1[0] == 'C' || token1[0] == 'c' ) + { + temp_surface.Calc_DDL_viscosity(true); + temp_surface.Set_DDL_viscosity(1.0); + continue; + } else if (j != EMPTY) { error_msg @@ -6601,10 +6619,10 @@ read_surface(void) i1++; continue; } - else if (i != EMPTY || i1 > 8) + else if (i != EMPTY || i1 > 4) { error_msg - ("Expected at most 8 numbers for the Donnan_factors for co- and counter-ions,\n z *= cgc[0] * (mu_x**(cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * mu_x**(cgc[4])", + ("Expected 4 numbers for the Donnan_factors of single and double-charged coounter- and co-ions,\n z1, z2, z_1, z_2", CONTINUE); error_msg(line_save, CONTINUE); input_error++; diff --git a/spread.cpp b/spread.cpp index 1f4421381..d2a385c35 100644 --- a/spread.cpp +++ b/spread.cpp @@ -162,8 +162,6 @@ read_solution_spread(void) { case 0: /* temp */ case 1: /* temperature */ - case 2: /* dens */ - case 3: /* density */ case 10: /* water */ if ((count == 2 || count == 3) && num == TRUE) { @@ -174,6 +172,18 @@ read_solution_spread(void) opt = OPTION_DEFAULT; } break; + case 2: /* dens */ + case 3: /* density */ + copy_token(token, &cptr); + if (count == 2 || count == 3 && (num == TRUE || token[0] == 'c' || token[0] == 'C')) + { + /* opt = opt; */ + } + else + { + opt = OPTION_DEFAULT; + } + break; case 6: /* ph */ case 7: /* pe */ if ((count == 2 || count == 3 || count == 4) @@ -285,32 +295,34 @@ read_solution_spread(void) break; case 2: /* density */ case 3: - //sscanf(next_char, SCANFORMAT, &(soln_defaults.density)); { - copy_token(token, &next_char); - if (sscanf(token.c_str(), SCANFORMAT, &dummy) != 1) + int j = copy_token(token, &next_char); + if (j == DIGIT) { + if (sscanf(token.c_str(), SCANFORMAT, &dummy) != 1) + { error_msg("Expecting numeric value for density.", PHRQ_io::OT_CONTINUE); error_msg(line_save, PHRQ_io::OT_CONTINUE); input_error++; + } + else + { + soln_defaults.density = dummy; + copy_token(token, &next_char); + if (token[0] == 'c' || token[0] == 'C') + soln_defaults.calc_density = true; + } } - else - { - soln_defaults.density = dummy; - } - int j = copy_token(token, &next_char); - if (j != EMPTY) + else if (j != EMPTY) { if (token[0] != 'c' && token[0] != 'C') { - error_msg("Only option following density is c[alculate].", PHRQ_io::OT_CONTINUE); + error_msg("Options following density are numeric value or c[alculate].", PHRQ_io::OT_CONTINUE); error_msg(line_save, PHRQ_io::OT_CONTINUE); input_error++; } else - { soln_defaults.calc_density = true; - } } } break; diff --git a/step.cpp b/step.cpp index bd960b36f..be02a1f43 100644 --- a/step.cpp +++ b/step.cpp @@ -327,6 +327,7 @@ xsolution_zero(void) solution_pe_x = 0.0; mu_x = 0.0; ah2o_x = 0.0; + viscos = 0.0; density_x = 0.0; total_h_x = 0.0; total_o_x = 0.0; @@ -379,6 +380,7 @@ add_solution(cxxSolution *solution_ptr, LDBLE extensive, LDBLE intensive) solution_pe_x += solution_ptr->Get_pe() * intensive; mu_x += solution_ptr->Get_mu() * intensive; ah2o_x += solution_ptr->Get_ah2o() * intensive; + viscos += solution_ptr->Get_viscosity() * intensive; density_x += solution_ptr->Get_density() * intensive; total_h_x += solution_ptr->Get_total_h() * extensive; diff --git a/structures.cpp b/structures.cpp index 1580f0ae2..618df1255 100644 --- a/structures.cpp +++ b/structures.cpp @@ -1561,6 +1561,7 @@ s_init(class species *s_ptr) s_ptr->dw_t = 0.0; s_ptr->dw_a = 0.0; s_ptr->dw_a2 = 0.0; + s_ptr->dw_a3 = 0.0; s_ptr->erm_ddl = 1.0; s_ptr->equiv = 0; s_ptr->alk = 0.0; diff --git a/transport.cpp b/transport.cpp index 0bd3d6491..31860db95 100644 --- a/transport.cpp +++ b/transport.cpp @@ -1891,7 +1891,10 @@ fill_spec(int l_cell_no, int ref_cell) * correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos * modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959 */ - viscos = viscos_0; + if (print_viscosity) + viscos = Utilities::Rxn_find(Rxn_solution_map, l_cell_no)->Get_viscosity(); + else + viscos = viscos_0; /* * put temperature factor in por_factor which corrects for porous medium... */ @@ -2221,7 +2224,8 @@ diffuse_implicit(LDBLE DDt, int stagnant) // Transport of aqueous species is summarized into master species. // With electro-migration, transport of anions and cations is calculated in opposite directions since (sign) J = - z * dV. // Only available moles are transported, thus are > 0, but if concentrations oscillate, - // change max_mixf in input file: -implicit true 1 # max_mixf = 1 (default). + // decrease max_mixf in input file: -implicit true 0.7 # max_mixf = 0.7 (default = 1). + // or define time substeps: -time 0.5e6 year 20 # substeps = 20 (default = 1). int i, icell, cp, comp; // ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells); int ifirst, ilast; @@ -3725,37 +3729,74 @@ fill_m_s(class J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant) /* ---------------------------------------------------------------------- */ void Phreeqc:: calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant) -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ +//{ +// ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j)); +// // At filterends, concentrations of ions change step-wise to the DL. +// // We take the harmonic mean for f_free, the average for the DL. +// if (ct[icell].v_m[k].z) +// { +// if (!g_i && g_j) +// { +// ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) + +// b_i * (1 - free_j) / 4 + b_j * g_j / 4; +// } +// else if (g_i && !g_j) +// ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) + +// b_j * (1 - free_i) / 4 + b_i * g_i / 4; +// } +// // for boundary cells... +// if (stagnant > 1) +// { /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell, +// and with the mixf * 2 for the boundary cells in the input... */ +// if (icell == 3 && !g_i && g_j) +// ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2; +// else if (jcell == all_cells - 1 && !g_j && g_i) +// ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2; +// } +// else +// { +// if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1)) +// ct[icell].v_m[k].b_ij = b_j * (free_j + g_j); +// else if (icell == count_cells && jcell == count_cells + 1) +// ct[icell].v_m[k].b_ij = b_i * (free_i + g_i); +// } +// if (ct[icell].v_m[k].z) +// ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z; +// return; +//} { - ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j)); - // At filterends, concentrations of ions change step-wise to the DL. - // We take the harmonic mean for f_free, the average for the DL. - if (ct[icell].v_m[k].z) + // Oct. 2023, with g_i,j = exp(g*z) * SS (charge_ptr-water / aq_x) + LDBLE fg_i = (1 - free_i) * g_i, + fg_j = (1 - free_j) * g_j; + ct[icell].v_m[k].b_ij = b_i * (free_i + fg_i) * b_j * (free_j + fg_j) / (b_i * (free_i + fg_i) + b_j * (free_j + fg_j)); + // At filterends and boundary cells, concentrations of ions change step-wise to the DL. + // filter cells, harmonic mean for f_free, the average for the DL. + if (icell != 0 && icell != count_cells && ct[icell].v_m[k].z) { if (!g_i && g_j) - { - ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) + - b_i * (1 - free_j) / 4 + b_j * g_j / 4; - } - else if (g_i && !g_j) + ct[icell].v_m[k].b_ij = b_i * free_j * b_j / (b_i + b_j) + + (b_i * (1 - free_j) + b_j * fg_j) / 4; + if (g_i && !g_j) ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) + - b_j * (1 - free_i) / 4 + b_i * g_i / 4; + (b_j * (1 - free_i) + b_i * fg_i) / 4; } - // for boundary cells... + // for boundary cells, all z... if (stagnant > 1) { /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell, and with the mixf * 2 for the boundary cells in the input... */ if (icell == 3 && !g_i && g_j) - ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2; + ct[icell].v_m[k].b_ij = b_j * (free_j + fg_j) / 2; else if (jcell == all_cells - 1 && !g_j && g_i) - ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2; + ct[icell].v_m[k].b_ij = b_i * (free_i + fg_i) / 2; } + // regular column... else { if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1)) - ct[icell].v_m[k].b_ij = b_j * (free_j + g_j); + ct[icell].v_m[k].b_ij = b_j * (free_j + fg_j); else if (icell == count_cells && jcell == count_cells + 1) - ct[icell].v_m[k].b_ij = b_i * (free_i + g_i); + ct[icell].v_m[k].b_ij = b_i * (free_i + fg_i); } if (ct[icell].v_m[k].z) ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z; @@ -3788,7 +3829,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) b_i_ani = A1 / (G_i * h_i / 2) * Dw * {f_free + (1 - f_free) / Bm)}. 22/2/18: now calculates diffusion through EDL's of multiple, differently charged surfaces * stagnant TRUE: - * same eqn for J_ij, but multplies with 2 * mixf. (times 2, because mixf = A / (G_i * h_i)) + * same eqn for J_ij, but multiplies with 2 * mixf. (times 2, because mixf = A / (G_i * h_i)) * mixf_ij = mixf / (Dw / init_tort_f) / new_tort_f * new_por / init_por * mixf is defined in MIX; Dw is default multicomponent diffusion coefficient; * init_tort_f equals multi_Dpor^(-multi_Dn); new_pf = new tortuosity factor. @@ -3886,6 +3927,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) only_counter = TRUE; ct[icell].visc1 = s_ptr1->Get_DDL_viscosity(); + if (s_ptr1->Get_calc_viscosity()) + ct[icell].visc1 /= Utilities::Rxn_find(Rxn_solution_map, icell)->Get_viscosity(); /* find the immobile surface charges with DL... */ for (i = 0; i < (int)s_charge_p.size(); i++) { @@ -3913,6 +3956,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) only_counter = TRUE; ct[icell].visc2 = s_ptr2->Get_DDL_viscosity(); + if (s_ptr2->Get_calc_viscosity()) + ct[icell].visc2 /= Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_viscosity(); for (i = 0; i < (int)s_charge_p.size(); i++) { @@ -4227,7 +4272,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) { g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; } - g_i *= sol_D[icell].spec[i].erm_ddl; + g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4250,7 +4295,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } } } - g_j *= sol_D[icell].spec[i].erm_ddl; + g_j *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc2; } } @@ -4343,7 +4388,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } } } - g_i *= sol_D[jcell].spec[j].erm_ddl; + g_i *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4351,7 +4396,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) { g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; } - g_j *= sol_D[jcell].spec[j].erm_ddl; + g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2; } } b_i = A1; @@ -4437,7 +4482,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) { g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; } - g_i *= sol_D[icell].spec[i].erm_ddl; + g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4445,7 +4490,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) { g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; } - g_j *= sol_D[jcell].spec[j].erm_ddl; + g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2; } } b_i = A1 * sol_D[icell].spec[i].Dwt; @@ -5881,9 +5926,12 @@ LDBLE new_Dw, int l_cell) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -viscosity(void) +viscosity(cxxSurface *surf_ptr) /* ---------------------------------------------------------------------- */ { + if (surf_ptr && !surf_ptr->Get_calc_viscosity()) + return surf_ptr->Get_DDL_viscosity(); + /* from Atkins, 1994. Physical Chemistry, 5th ed. */ //viscos = @@ -5892,7 +5940,7 @@ viscosity(void) // 0.000836 * (tc_x - 20) * (tc_x - 20)) / (109 + tc_x)); /* Huber et al., 2009, J. Phys. Chem. Ref. Data, Vol. 38, 101-125 */ LDBLE H[4] = { 1.67752, 2.20462, 0.6366564, -0.241605 }; - LDBLE Tb = tk_x / 647.096, denom = H[0], mu0; + LDBLE Tb = (tc_x + 273.15) / 647.096, denom = H[0], mu0; int i, j, i1; for (i = 1; i < 4; i++) denom += H[i] / pow(Tb, i); @@ -5964,124 +6012,216 @@ viscosity(void) both weighted by the equivalent concentration. */ LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, t3, fan = 1; - LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, tc; + LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, tc, l_moles, l_water, l_mu_x, dw_t_visc; m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = 0; tc = (tc_x > 200) ? 200 : tc_x; - - for (i = 0; i < (int)this->s_x.size(); i++) + l_water = mass_water_aq_x; + l_mu_x = mu_x; + + + int i1_last; + if (surf_ptr == NULL) + i1_last = 1; + else { - if (s_x[i]->type != AQ && s_x[i]->type > HPLUS) - continue; - if ((lm = s_x[i]->lm) < -9) - continue; - if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3]) + i1_last = (int)surf_ptr->Get_surface_charges().size(); + if (i1_last > 1) + i1_last += 1; + } + std::map z_g_map; + cxxSurfaceCharge s_charge_p; + LDBLE ratio_surf_aq = mass_water_surfaces_x / mass_water_aq_x; + for (i1 = 0; i1 < i1_last; i1++) + { + Bc = Dc = Dw = m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = 0; + if (surf_ptr) { - s_x[i]->dw_t_visc = 0; - t1 = s_x[i]->moles / mass_water_aq_x; - l_z = fabs(s_x[i]->z); - if (l_z) - f_z = (l_z * l_z + l_z) / 2; - else - f_z = mu_x / t1; - //if data at tc's other than 25 are scarce, put the values found for 25 C in [7] and [8], optimize [1], [2], and [4]... - if (s_x[i]->Jones_Dole[7] || s_x[i]->Jones_Dole[8]) - { - s_x[i]->Jones_Dole[0] = s_x[i]->Jones_Dole[7] - - s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * 25.0); - s_x[i]->Jones_Dole[3] = - s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0); - } - // find B * m and D * m * mu^d3 - s_x[i]->dw_t_visc = (s_x[i]->Jones_Dole[0] + - s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1; - Bc += s_x[i]->dw_t_visc; - // define f_I from the exponent of the D * m^d3 term... - if (s_x[i]->Jones_Dole[5] >= 1) - t2 = mu_x / 3 / s_x[i]->Jones_Dole[5]; - else if (s_x[i]->Jones_Dole[5] > 0.4) - t2 = -0.8 / s_x[i]->Jones_Dole[5]; + z_g_map.clear(); + if (i1_last == 1 || i1 < i1_last - 1) + { + s_charge_p = surf_ptr->Get_surface_charges()[i1]; + l_water = s_charge_p.Get_mass_water(); + z_g_map.insert(s_charge_p.Get_z_gMCD_map().begin(), s_charge_p.Get_z_gMCD_map().end()); + for (auto& x : z_g_map) + x.second *= ratio_surf_aq; + } else - t2 = -1; - t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) * - t1 * (pow(mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2); - s_x[i]->dw_t_visc += t3; - Dc += t3; - //output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc )); - } - // parms for A... - if ((l_z = s_x[i]->z) == 0) - continue; - Dw = s_x[i]->dw; - if (Dw) - { - Dw *= (0.89 / viscos_0 * tk_x / 298.15); - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); + { + s_charge_p = surf_ptr->Get_surface_charges()[0]; + z_g_map.insert(s_charge_p.Get_z_gMCD_map().begin(), s_charge_p.Get_z_gMCD_map().end()); + for (i = 1; i < i1_last - 1; i++) + { + s_charge_p = surf_ptr->Get_surface_charges()[i]; + for (auto& x : z_g_map) + x.second += s_charge_p.Get_z_gMCD_map()[x.first]; + } + for (auto& x : z_g_map) + x.second *= ratio_surf_aq; + l_water = mass_water_surfaces_x; + } + l_mu_x = eq_plus = eq_min = 0; + for (i = 0; i < (int)this->s_x.size(); i++) + { + if (s_x[i]->type != AQ && s_x[i]->type > HPLUS) + continue; + if (s_x[i]->lm < -9 || s_x[i] == 0) + continue; + l_moles = s_x[i]->moles * s_x[i]->z * z_g_map.find(s_x[i]->z)->second; + if (s_x[i]->z < 0) + eq_min += l_moles; + else + eq_plus += l_moles; + l_mu_x += l_moles * s_x[i]->z; + } + l_mu_x += fabs(eq_plus + eq_min); + l_mu_x /= (2 * l_water); + eq_plus = eq_min = 0; } - if (l_z < 0) + + for (i = 0; i < (int)this->s_x.size(); i++) { - if (!strcmp(s_x[i]->name, "Cl-")) - // volumina for f_an... + if (s_x[i]->type != AQ && s_x[i]->type > HPLUS) + continue; + if ((lm = s_x[i]->lm) < -9) + continue; + l_moles = s_x[i]->moles; + if (surf_ptr != NULL) { - V_Cl = s_x[i]->logk[vm_tc]; - V_an += V_Cl * s_x[i]->moles; - m_an += s_x[i]->moles; + l_moles *= z_g_map.find(s_x[i]->z)->second; } - else// if (s_x[i]->Jones_Dole[6]) + if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3]) { - V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * s_x[i]->moles; - m_an += s_x[i]->moles; - } + dw_t_visc = 0; + t1 = l_moles / l_water; + l_z = fabs(s_x[i]->z); + if (l_z) + f_z = (l_z * l_z + l_z) / 2; + else + f_z = l_mu_x / t1; + //if data at tc's other than 25 are scarce, put the values found for 25 C in [7] and [8], optimize [1], [2], and [4]... + if (s_x[i]->Jones_Dole[7] || s_x[i]->Jones_Dole[8]) + { + s_x[i]->Jones_Dole[0] = s_x[i]->Jones_Dole[7] - + s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * 25.0); + s_x[i]->Jones_Dole[3] = + s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0); + } + // find B * m and D * m * mu^d3 + dw_t_visc = (s_x[i]->Jones_Dole[0] + + s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1; + Bc += dw_t_visc; + // define f_I from the exponent of the D * m^d3 term... + if (s_x[i]->Jones_Dole[5] >= 1) + t2 = l_mu_x / 3 / s_x[i]->Jones_Dole[5]; + else if (s_x[i]->Jones_Dole[5] > 0.4) + t2 = -0.8 / s_x[i]->Jones_Dole[5]; + else + t2 = -1; + t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) * + t1 * (pow(l_mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2); + if (t3 < -1e-5) // add this check + t3 = 0; + Dc += t3; + if (!surf_ptr) s_x[i]->dw_t_visc = dw_t_visc + t3; + //output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc )); + } + // parms for A... + if ((l_z = s_x[i]->z) == 0) + continue; + Dw = s_x[i]->dw; if (Dw) { - // anions for A... - m_min += s_x[i]->moles; - t1 = s_x[i]->moles * l_z; - eq_min -= t1; - eq_dw_min -= t1 / Dw; + Dw *= (0.89 / viscos_0 * tk_x / 298.15); + if (s_x[i]->dw_t) + Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); + } + if (l_z < 0) + { + if (!strcmp(s_x[i]->name, "Cl-")) + // volumina for f_an... + { + V_Cl = s_x[i]->logk[vm_tc]; + V_an += V_Cl * l_moles; + m_an += l_moles; + } + else// if (s_x[i]->Jones_Dole[6]) + { + V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * l_moles; + m_an += l_moles; + } + if (Dw) + { + // anions for A... + m_min += l_moles; + t1 = l_moles * l_z; + eq_min -= t1; + eq_dw_min -= t1 / Dw; + } + } + else if (Dw) + { + // cations for A... + m_plus += l_moles; + t1 = l_moles * l_z; + eq_plus += t1; + eq_dw_plus += t1 / Dw; } } - else if (Dw) + if (m_plus && m_min && eq_dw_plus && eq_dw_min) { - // cations for A... - m_plus += s_x[i]->moles; - t1 = s_x[i]->moles * l_z; - eq_plus += t1; - eq_dw_plus += t1 / Dw; - } - } - if (m_plus && m_min && eq_dw_plus && eq_dw_min) - { - z1 = eq_plus / m_plus; z2 = eq_min / m_min; - D1 = eq_plus / eq_dw_plus; D2 = eq_min / eq_dw_min; + z1 = eq_plus / m_plus; z2 = eq_min / m_min; + D1 = eq_plus / eq_dw_plus; D2 = eq_min / eq_dw_min; - t1 = (D1 - D2) / (sqrt(D1 * z1 + D2 * z2) + sqrt((D1 + D2) * (z1 + z2))); - psi = (D1 * z2 + D2 * z1) / 4.0 - z1 * z2 * t1 * t1; - // Here A is A * viscos_0, avoids multiplication later on... - A = 4.3787e-14 * pow(tk_x, 1.5) / (sqrt(eps_r * (z1 + z2) / ((z1 > z2) ? z1 : z2)) * (D1 * D2)) * psi; - } - else - A = 0; - viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / mass_water_aq_x); - if (m_an) - V_an /= m_an; - if (!V_Cl) - V_Cl = calc_vm_Cl(); - if (V_an) - fan = 2 - V_an / V_Cl; - //else - // fan = 1; - viscos += viscos_0 * fan * (Bc + Dc); - if (viscos < 0) - { - viscos = viscos_0; - warning_msg("viscosity < 0, reset to viscosity of pure water"); - } - for (i = 0; i < (int)this->s_x.size(); i++) - { - s_x[i]->dw_t_visc /= (Bc + Dc); + t1 = (D1 - D2) / (sqrt(D1 * z1 + D2 * z2) + sqrt((D1 + D2) * (z1 + z2))); + psi = (D1 * z2 + D2 * z1) / 4.0 - z1 * z2 * t1 * t1; + // Here A is A * viscos_0, avoids multiplication later on... + A = 4.3787e-14 * pow(tk_x, 1.5) / (sqrt(eps_r * (z1 + z2) / ((z1 > z2) ? z1 : z2)) * (D1 * D2)) * psi; + } + else + A = 0; + viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / l_water); + if (m_an) + V_an /= m_an; + if (!V_Cl) + V_Cl = calc_vm_Cl(); + if (V_an) + fan = 2 - V_an / V_Cl; + //else + // fan = 1; + if (Dc < 0) + Dc = 0; // provisional... + viscos += viscos_0 * fan * (Bc + Dc); + if (viscos < 0) + { + viscos = viscos_0; + warning_msg("viscosity < 0, reset to viscosity of pure water"); + } + if (!surf_ptr) + { + for (i = 0; i < (int)this->s_x.size(); i++) + { + if (s_x[i]->type != AQ && s_x[i]->type > HPLUS) + continue; + if ((lm = s_x[i]->lm) < -9) + continue; + if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3]) + s_x[i]->dw_t_visc /= (Bc + Dc); + } + } + else //if (surf_ptr->Get_calc_viscosity()) + { + if (i1_last == 1) + { + surf_ptr->Get_surface_charges()[i1].Set_DDL_viscosity(viscos); + surf_ptr->Set_DDL_viscosity(viscos); + } + else if (i1 < i1_last - 1) + surf_ptr->Get_surface_charges()[i1].Set_DDL_viscosity(viscos); + else if (i1 == i1_last - 1) + surf_ptr->Set_DDL_viscosity(viscos); + } } return viscos; }