max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 08:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge commit 'ccdbd0ce2040eb00006775fe0cf2ab2faf4138b5'

Browse Source
This commit is contained in:
Darth Vader 2024-04-17 00:17:30 +00:00
commit a9dc9b84f8
23 changed files with 1145 additions and 790 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/phreeqcpp/Phreeqc.cpp
View File

@ -583,7 +583,6 @@ void Phreeqc::init(void)
solution_pe_x = 0; solution_pe_x = 0;
mu_x = 0; mu_x = 0;
ah2o_x = 1.0; ah2o_x = 1.0;
density_x = 0;
total_h_x = 0; total_h_x = 0;
total_o_x = 0; total_o_x = 0;
cb_x = 0; cb_x = 0;
@ -898,6 +897,7 @@ void Phreeqc::init(void)
viscos = 0.0; viscos = 0.0;
viscos_0 = 0.0; viscos_0 = 0.0;
viscos_0_25 = 0.0; viscos_0_25 = 0.0;
density_x = 0.0;
rho_0 = 0.0; rho_0 = 0.0;
kappa_0 = 0.0; kappa_0 = 0.0;
p_sat = 0.0; p_sat = 0.0;
@ -1714,6 +1714,7 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc)
viscos = pSrc->viscos; viscos = pSrc->viscos;
viscos_0 = pSrc->viscos_0; 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 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_pore_volume = pSrc->cell_pore_volume;;
cell_porosity = pSrc->cell_porosity; cell_porosity = pSrc->cell_porosity;
cell_volume = pSrc->cell_volume; cell_volume = pSrc->cell_volume;
@ -1722,9 +1723,6 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc)
sys_tot = pSrc->sys_tot; sys_tot = pSrc->sys_tot;
// solution properties // solution properties
V_solutes = pSrc->V_solutes; 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; rho_0 = pSrc->rho_0;
kappa_0 = pSrc->kappa_0; kappa_0 = pSrc->kappa_0;
p_sat = pSrc->p_sat; p_sat = pSrc->p_sat;

7
src/phreeqcpp/Phreeqc.h
View File

@ -283,7 +283,7 @@ public:
int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1); int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1);
int calc_all_donnan(void); int calc_all_donnan(void);
int calc_init_donnan(void); int calc_init_donnan(void);
LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr); LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, LDBLE f_free, std::vector<LDBLE> &zcorr);
LDBLE g_function(LDBLE x_value); LDBLE g_function(LDBLE x_value);
LDBLE midpnt(LDBLE x1, LDBLE x2, int n); LDBLE midpnt(LDBLE x1, LDBLE x2, int n);
void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv, void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv,
@ -555,6 +555,7 @@ public:
LDBLE calc_PR(std::vector<class phase*> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m); LDBLE calc_PR(std::vector<class phase*> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m);
LDBLE calc_PR(); LDBLE calc_PR();
int calc_vm(LDBLE tc, LDBLE pa); int calc_vm(LDBLE tc, LDBLE pa);
LDBLE calc_vm0(const char *species_name, LDBLE tc, LDBLE pa, LDBLE mu);
int clear(void); int clear(void);
int convert_units(cxxSolution* solution_ptr); int convert_units(cxxSolution* solution_ptr);
class unknown* find_surface_charge_unknown(std::string& str_ptr, int plane); class unknown* find_surface_charge_unknown(std::string& str_ptr, int plane);
@ -995,7 +996,7 @@ public:
LDBLE new_Dw); LDBLE new_Dw);
int reformat_surf(const char* comp_name, LDBLE fraction, const char* new_comp_name, int reformat_surf(const char* comp_name, LDBLE fraction, const char* new_comp_name,
LDBLE new_Dw, int cell); LDBLE new_Dw, int cell);
LDBLE viscosity(void); LDBLE viscosity(cxxSurface *surf_ptr);
LDBLE calc_f_visc(const char *name); LDBLE calc_f_visc(const char *name);
LDBLE calc_vm_Cl(void); LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant); int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
@ -1274,7 +1275,6 @@ protected:
LDBLE solution_pe_x; LDBLE solution_pe_x;
LDBLE mu_x; LDBLE mu_x;
LDBLE ah2o_x; LDBLE ah2o_x;
LDBLE density_x;
LDBLE total_h_x; LDBLE total_h_x;
LDBLE total_o_x; LDBLE total_o_x;
LDBLE cb_x; LDBLE cb_x;
@ -1614,6 +1614,7 @@ protected:
int print_viscosity; int print_viscosity;
LDBLE viscos, viscos_0, viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C 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_pore_volume;
LDBLE cell_porosity; LDBLE cell_porosity;
LDBLE cell_volume; LDBLE cell_volume;

34
src/phreeqcpp/Surface.cxx
View File

@ -36,9 +36,10 @@ cxxSurface::cxxSurface(PHRQ_io *io)
dl_type = NO_DL; dl_type = NO_DL;
sites_units = SITES_ABSOLUTE; sites_units = SITES_ABSOLUTE;
only_counter_ions = false; only_counter_ions = false;
correct_GC = false; correct_D = false;
thickness = 1e-8; thickness = 1e-8;
debye_lengths = 0.0; debye_lengths = 0.0;
calc_DDL_viscosity = false;
DDL_viscosity = 1.0; DDL_viscosity = 1.0;
DDL_limit = 0.8; DDL_limit = 0.8;
transport = false; transport = false;
@ -56,9 +57,10 @@ cxxNumKeyword(io)
dl_type = NO_DL; dl_type = NO_DL;
sites_units = SITES_ABSOLUTE; sites_units = SITES_ABSOLUTE;
only_counter_ions = false; only_counter_ions = false;
correct_GC = false; correct_D = false;
thickness = 1e-8; thickness = 1e-8;
debye_lengths = 0.0; debye_lengths = 0.0;
calc_DDL_viscosity = false;
DDL_viscosity = 1.0; DDL_viscosity = 1.0;
DDL_limit = 0.8; DDL_limit = 0.8;
transport = false; 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 << indent1;
s_oss << "-only_counter_ions " << this->only_counter_ions << "\n"; s_oss << "-only_counter_ions " << this->only_counter_ions << "\n";
s_oss << indent1; s_oss << indent1;
s_oss << "-correct_GC " << this->correct_GC << "\n"; s_oss << "-correct_D " << this->correct_D << "\n";
s_oss << indent1; s_oss << indent1;
s_oss << "-thickness " << this->thickness << "\n"; s_oss << "-thickness " << this->thickness << "\n";
s_oss << indent1; s_oss << indent1;
@ -193,7 +195,7 @@ cxxSurface::read_raw(CParser & parser, bool check)
this->Set_tidied(true); this->Set_tidied(true);
bool only_counter_ions_defined(false); bool only_counter_ions_defined(false);
//bool correct_GC_defined(false); //bool correct_D_defined(false);
bool thickness_defined(false); bool thickness_defined(false);
bool type_defined(false); bool type_defined(false);
bool dl_type_defined(false); bool dl_type_defined(false);
@ -395,7 +397,7 @@ cxxSurface::read_raw(CParser & parser, bool check)
case 11: // DDL_viscosity case 11: // DDL_viscosity
if (!(parser.get_iss() >> this->DDL_viscosity)) if (!(parser.get_iss() >> this->DDL_viscosity))
{ {
this->DDL_viscosity = 0.0; this->DDL_viscosity = 1.0;
parser.incr_input_error(); parser.incr_input_error();
parser.error_msg("Expected numeric value for DDL_viscosity.", parser.error_msg("Expected numeric value for DDL_viscosity.",
PHRQ_io::OT_CONTINUE); PHRQ_io::OT_CONTINUE);
@ -473,16 +475,16 @@ cxxSurface::read_raw(CParser & parser, bool check)
PHRQ_io::OT_CONTINUE); PHRQ_io::OT_CONTINUE);
} }
break; break;
case 19: // correct_GC case 19: // correct_D
if (!(parser.get_iss() >> this->correct_GC)) if (!(parser.get_iss() >> this->correct_D))
{ {
this->correct_GC = false; this->correct_D = false;
parser.incr_input_error(); parser.incr_input_error();
parser. parser.
error_msg("Expected boolean value for correct_GC.", error_msg("Expected boolean value for correct_D.",
PHRQ_io::OT_CONTINUE); PHRQ_io::OT_CONTINUE);
} }
//correct_GC_defined = true; //correct_D_defined = true;
break; break;
} }
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) 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.", error_msg("Only_counter_ions not defined for SURFACE_RAW input.",
PHRQ_io::OT_CONTINUE); PHRQ_io::OT_CONTINUE);
} }
//if (correct_GC_defined == false) //if (correct_D_defined == false)
//{ //{
// parser.incr_input_error(); // parser.incr_input_error();
// parser. // 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); // PHRQ_io::OT_CONTINUE);
//} //}
if (thickness_defined == false) if (thickness_defined == false)
@ -582,7 +584,7 @@ cxxSurface::add(const cxxSurface & addee_in, LDBLE extensive)
if (this->surface_comps.size() == 0) if (this->surface_comps.size() == 0)
{ {
this->only_counter_ions = addee.only_counter_ions; 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->dl_type = addee.dl_type;
this->type = addee.type; this->type = addee.type;
this->sites_units = addee.sites_units; 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->debye_lengths);
doubles.push_back(this->DDL_viscosity); doubles.push_back(this->DDL_viscosity);
doubles.push_back(this->DDL_limit); 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); ints.push_back(this->transport ? 1 : 0);
this->totals.Serialize(dictionary, ints, doubles); this->totals.Serialize(dictionary, ints, doubles);
ints.push_back(this->solution_equilibria ? 1 : 0); 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->debye_lengths = doubles[dd++];
this->DDL_viscosity = doubles[dd++]; this->DDL_viscosity = doubles[dd++];
this->DDL_limit = 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->transport = (ints[ii++] != 0);
this->totals.Deserialize(dictionary, ints, doubles, ii, dd); this->totals.Deserialize(dictionary, ints, doubles, ii, dd);
this->solution_equilibria = (ints[ii++] != 0); 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("n_solution"), // 16
std::vector< std::string >::value_type("totals"), // 17 std::vector< std::string >::value_type("totals"), // 17
std::vector< std::string >::value_type("tidied"), // 18 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]); const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

9
src/phreeqcpp/Surface.h
View File

@ -65,10 +65,12 @@ public:
void Set_debye_lengths(LDBLE t) {debye_lengths = t;} void Set_debye_lengths(LDBLE t) {debye_lengths = t;}
LDBLE Get_DDL_viscosity(void) const {return DDL_viscosity;} LDBLE Get_DDL_viscosity(void) const {return DDL_viscosity;}
void Set_DDL_viscosity(LDBLE t) {DDL_viscosity = t;} 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;} LDBLE Get_DDL_limit(void) const {return DDL_limit;}
void Set_DDL_limit(LDBLE t) {DDL_limit = t;} void Set_DDL_limit(LDBLE t) {DDL_limit = t;}
bool Get_correct_GC(void) const { return correct_GC; } bool Get_correct_D(void) const { return correct_D; }
void Set_correct_GC(bool tf) { correct_GC = tf; } void Set_correct_D(bool tf) { correct_D = tf; }
std::vector<LDBLE> Donnan_factors; std::vector<LDBLE> Donnan_factors;
bool Get_transport(void) const {return transport;} bool Get_transport(void) const {return transport;}
void Set_transport(bool tf) {transport = tf;} void Set_transport(bool tf) {transport = tf;}
@ -93,8 +95,9 @@ protected:
LDBLE thickness; LDBLE thickness;
LDBLE debye_lengths; LDBLE debye_lengths;
LDBLE DDL_viscosity; LDBLE DDL_viscosity;
bool calc_DDL_viscosity;
LDBLE DDL_limit; LDBLE DDL_limit;
bool correct_GC; bool correct_D;
bool transport; bool transport;
cxxNameDouble totals; cxxNameDouble totals;
bool solution_equilibria; bool solution_equilibria;

45
src/phreeqcpp/SurfaceCharge.cxx
View File

@ -36,6 +36,8 @@ PHRQ_base(io)
grams = 0.0; grams = 0.0;
charge_balance = 0.0; charge_balance = 0.0;
mass_water = 0.0; mass_water = 0.0;
DDL_viscosity = 0.0;
f_free = 0.0;
la_psi = 0.0; la_psi = 0.0;
capacitance[0] = 1.0; capacitance[0] = 1.0;
capacitance[1] = 5.0; capacitance[1] = 5.0;
@ -68,6 +70,7 @@ cxxSurfaceCharge::dump_xml(std::ostream & s_oss, unsigned int indent) const
charge_balance << "\"" << "\n"; charge_balance << "\"" << "\n";
s_oss << indent0 << "mass_water=\"" << this-> s_oss << indent0 << "mass_water=\"" << this->
mass_water << "\"" << "\n"; mass_water << "\"" << "\n";
s_oss << indent0 << "f_free=\"" << this->f_free << "\"" << "\n";
s_oss << indent0 << "la_psi=\"" << this->la_psi << "\"" << "\n"; s_oss << indent0 << "la_psi=\"" << this->la_psi << "\"" << "\n";
s_oss << indent0 << "capacitance=\"" << this-> s_oss << indent0 << "capacitance=\"" << this->
capacitance[0] << " " << this->capacitance[0] << "\"" << "\n"; 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 << "-grams " << this->grams << "\n";
s_oss << indent0 << "-charge_balance " << this->charge_balance << "\n"; s_oss << indent0 << "-charge_balance " << this->charge_balance << "\n";
s_oss << indent0 << "-mass_water " << this->mass_water << "\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 << "-la_psi " << this->la_psi << "\n";
s_oss << indent0 << "-capacitance0 " << this->capacitance[0] << "\n"; s_oss << indent0 << "-capacitance0 " << this->capacitance[0] << "\n";
s_oss << indent0 << "-capacitance1 " << this->capacitance[1] << "\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 capacitance0_defined(false);
bool capacitance1_defined(false); bool capacitance1_defined(false);
bool g_map_first(true); bool g_map_first(true);
bool DDL_viscosity_defined(false);
for (;;) for (;;)
{ {
@ -225,7 +231,6 @@ cxxSurfaceCharge::read_raw(CParser & parser, bool check)
mass_water_defined = true; mass_water_defined = true;
break; break;
case 5: // la_psi case 5: // la_psi
if (!(parser.get_iss() >> this->la_psi)) if (!(parser.get_iss() >> this->la_psi))
{ {
@ -366,10 +371,27 @@ cxxSurfaceCharge::read_raw(CParser & parser, bool check)
} }
} }
opt_save = 16; opt_save = 16;
break; 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) if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break; break;
@ -454,9 +476,11 @@ cxxSurfaceCharge::add(const cxxSurfaceCharge & addee, LDBLE extensive)
this->mass_water += addee.mass_water * extensive; this->mass_water += addee.mass_water * extensive;
this->la_psi = this->la_psi * f1 + addee.la_psi * f2; this->la_psi = this->la_psi * f1 + addee.la_psi * f2;
this->capacitance[0] = 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] =
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); 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->sigma1);
doubles.push_back(this->sigma2); doubles.push_back(this->sigma2);
doubles.push_back(this->sigmaddl); 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()); ints.push_back((int) this->g_map.size());
{ {
std::map<LDBLE, cxxSurfDL>::iterator it; std::map<LDBLE, cxxSurfDL>::iterator it;
@ -523,6 +549,8 @@ cxxSurfaceCharge::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->sigma1 = doubles[dd++]; this->sigma1 = doubles[dd++];
this->sigma2 = doubles[dd++]; this->sigma2 = doubles[dd++];
this->sigmaddl = doubles[dd++]; this->sigmaddl = doubles[dd++];
this->f_free = doubles[dd++];
this->DDL_viscosity = doubles[dd++];
{ {
this->g_map.clear(); this->g_map.clear();
int count = ints[ii++]; 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("sigma2"), // 13
std::vector< std::string >::value_type("sigmaddl"), // 14 std::vector< std::string >::value_type("sigmaddl"), // 14
std::vector< std::string >::value_type("g_map"), // 15 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]); const std::vector< std::string > cxxSurfaceCharge::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

6
src/phreeqcpp/SurfaceCharge.h
View File

@ -87,6 +87,10 @@ public:
void Set_charge_balance(LDBLE d) {this->charge_balance = d;} void Set_charge_balance(LDBLE d) {this->charge_balance = d;}
LDBLE Get_mass_water() const {return this->mass_water;} LDBLE Get_mass_water() const {return this->mass_water;}
void Set_mass_water(LDBLE d) {this->mass_water = d;} 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;} LDBLE Get_la_psi() const {return this->la_psi;}
void Set_la_psi(LDBLE d) {this->la_psi = d;} void Set_la_psi(LDBLE d) {this->la_psi = d;}
LDBLE Get_capacitance0() const {return this->capacitance[0];} LDBLE Get_capacitance0() const {return this->capacitance[0];}
@ -117,6 +121,8 @@ protected:
LDBLE grams; LDBLE grams;
LDBLE charge_balance; LDBLE charge_balance;
LDBLE mass_water; LDBLE mass_water;
LDBLE DDL_viscosity;
LDBLE f_free;
LDBLE la_psi; LDBLE la_psi;
LDBLE capacitance[2]; LDBLE capacitance[2];
cxxNameDouble diffuse_layer_totals; cxxNameDouble diffuse_layer_totals;

737
src/phreeqcpp/basicsubs.cpp
View File

File diff suppressed because it is too large Load Diff

6
src/phreeqcpp/class_main.cpp
View File

@ -299,7 +299,7 @@ write_banner(void)
/* version */ /* version */
#ifdef NPP #ifdef NPP
len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "3.7.1"); len = sprintf(buffer, "* PHREEQC-%s *", "3.7.3");
#else #else
len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "@VERSION@"); len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "@VERSION@");
#endif #endif
@ -323,7 +323,7 @@ write_banner(void)
/* date */ /* date */
#ifdef NPP #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 #else
len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@"); len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@");
#endif #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(" Input file: %s\n", in_file.c_str()));
output_msg(sformatf(" Output file: %s\n", out_file.c_str())); output_msg(sformatf(" Output file: %s\n", out_file.c_str()));
#ifdef NPP #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 #endif
output_msg(sformatf("Database file: %s\n\n", token.c_str())); output_msg(sformatf("Database file: %s\n\n", token.c_str()));
#ifdef NPP #ifdef NPP

3
src/phreeqcpp/gases.cpp
View File

@ -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) * 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)); 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)); phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi));
} }
else else
phi = -3.0; // fugacity coefficient = 0.05 phi = -4.6; // fugacity coefficient = 0.01
phase_ptr->pr_phi = exp(phi); phase_ptr->pr_phi = exp(phi);
phase_ptr->pr_si_f = phi / LOG_10; // pr_si_f updated phase_ptr->pr_si_f = phi / LOG_10; // pr_si_f updated
// **** // ****

4
src/phreeqcpp/global_structures.h
View File

@ -716,9 +716,10 @@ public:
dw = 0; dw = 0;
// correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15) // correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15)
dw_t = 0; 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_a = 0;
dw_a2 = 0; dw_a2 = 0;
dw_a3 = 0;
dw_a_visc = 0; // viscosity correction of SC dw_a_visc = 0; // viscosity correction of SC
dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC
dw_t_visc = 0; // contribution to viscosity dw_t_visc = 0; // contribution to viscosity
@ -781,6 +782,7 @@ public:
LDBLE dw_t; LDBLE dw_t;
LDBLE dw_a; LDBLE dw_a;
LDBLE dw_a2; LDBLE dw_a2;
LDBLE dw_a3;
LDBLE dw_a_visc; LDBLE dw_a_visc;
LDBLE dw_t_SC; LDBLE dw_t_SC;
LDBLE dw_t_visc; LDBLE dw_t_visc;

225
src/phreeqcpp/integrate.cpp
View File

@ -733,40 +733,25 @@ calc_all_donnan(void)
{ {
bool converge; bool converge;
int cd_m; 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_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; 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); return (OK);
f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
tk_x * mu_x); tk_x * mu_x);
bool only_count = use.Get_surface_ptr()->Get_only_counter_ions(); bool only_count = surf_ptr->Get_only_counter_ions();
bool correct_GC = use.Get_surface_ptr()->Get_correct_GC(); bool correct_D = surf_ptr->Get_correct_D();
/* calculate g for each surface... /* calculate g for each surface...
*/ */
if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths() || if (!calculating_deriv || surf_ptr->Get_debye_lengths() ||
correct_GC) // DL_pitz && correct_GC correct_D) // DL_pitz && correct_D
initial_surface_water(); initial_surface_water();
LDBLE nDbl = 1; // z1, z2, fr_cat2 are the counter-ions, z_1, z_2, fr_ani2 are for co-ions.
if (correct_GC) LDBLE nDbl = 1, db_lim = 2, f_free, fr_cat2, fr_ani2;
{ LDBLE z1, z2, z_1, z_2;
if ((nDbl = use.Get_surface_ptr()->Get_debye_lengths()) == 0) z1 = z2 = z_1 = z_2 = f_free = fr_cat2 = fr_ani2 = 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);
}
}
converge = TRUE;
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);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf("Calc_all_g, X[%d]\n", j));
/* /*
* sum eq of each charge number in solution... * sum eq of each charge number in solution...
*/ */
@ -782,9 +767,44 @@ calc_all_donnan(void)
continue; continue;
charge_group_map[s_x[i]->z] += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl; 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 = surf_ptr->Get_debye_lengths()) == 0)
{
LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3);
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 = surf_ptr->Find_charge(x[j]->surface_charge);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf("Calc_all_g, X[%d]\n", j));
/* find surface charge from potential... */ /* find surface charge from potential... */
A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams(); 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 = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2; 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; surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
LDBLE lim_seq = 5e3; LDBLE lim_seq = 5e3;
if (correct_GC) lim_seq = 5e1; if (correct_D) lim_seq = 5e3;
if (fabs(surf_chrg_eq) > lim_seq) if (fabs(surf_chrg_eq) > lim_seq)
{ {
surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq); surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq);
@ -816,23 +836,24 @@ calc_all_donnan(void)
std::vector<LDBLE> zcorr(charge_group_map.size()); std::vector<LDBLE> zcorr(charge_group_map.size());
std::vector<LDBLE> zcorr2(charge_group_map.size()); std::vector<LDBLE> zcorr2(charge_group_map.size());
//LDBLE fD = 0; //LDBLE fD = 0;
psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, zcorr); 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, zcorr2); 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)); */ /*output_msg(sformatf( "psi's %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */
/* fill in g's */ /* fill in g's */
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; 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; int z_iter = 0;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{ {
LDBLE z = it->first, z1 = z; LDBLE z = it->first, z1 = z;
co_ion = surf_chrg_eq * z; co_ion = surf_chrg_eq * z;
if (correct_GC) if (correct_D)
z1 = zcorr[z_iter]; z1 = zcorr[z_iter];
//z1 *= cgc[0] * pow(z_factor, abs(z));
if (!ratio_aq) if (!ratio_aq)
{ {
@ -878,18 +899,18 @@ calc_all_donnan(void)
/* save Boltzmann factor * water fraction for MCD calc's in transport */ /* save Boltzmann factor * water fraction for MCD calc's in transport */
if (converge) if (converge)
{ {
if (only_count) if (only_count && co_ion > 0) // co-ions are not in the DL
{
if (co_ion > 0) // co-ions are not in the DL
charge_ptr->Get_z_gMCD_map()[z] = 0; 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 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++; z_iter++;
} }
charge_ptr->Set_f_free(f_free);
if (debug_diffuse_layer == TRUE) if (debug_diffuse_layer == TRUE)
{ {
std::string name = x[j]->master[0]->elt->name; 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; 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); return (OK);
f_sinh = f_sinh =
//sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * //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) if (x[j]->type != SURFACE_CB)
continue; 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(); charge_ptr->Get_g_map().clear();
/* find surface charge from potential... */ /* find surface charge from potential... */
A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams(); 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 = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2; 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; surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
/* find psi_avg that matches surface charge... */ /* 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 */ /* fill in g's */
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; 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)); 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)
|| (surf_chrg_eq > 0 && z > 0))) || (surf_chrg_eq > 0 && z > 0)))
charge_ptr->Get_g_map()[z].Set_g(-ratio_aq); charge_ptr->Get_g_map()[z].Set_g(-ratio_aq);
@ -1038,45 +1060,84 @@ calc_init_donnan(void)
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
LDBLE Phreeqc:: LDBLE Phreeqc::
calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr) calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, LDBLE f_free, std::vector<LDBLE> &zcorr)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
/* /*
* calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge * 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; ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
p = 0; p = 0;
if (surf_chrg_eq == 0 || ratio_aq == 0) if (surf_chrg_eq == 0 || ratio_aq == 0)
return (0.0); return (0.0);
else if (surf_chrg_eq < 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) 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 * 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: * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
* g(p) = (exp(-p *z_i) - 1) * ratio_aq * g(p) = (exp(-p *z_i) - 1) * ratio_aq
* with correct_GC true: * with correct_D true and f_free > 0:
* correct ions to better match the integrated concentrations: c_edl = c_free * (f_free + (1 - f_free) * exp(-p * z_i))
z == 1? z *= 0.285 cgc[6] * with correct_D true and f_free == 0:
z == 2? z *= 0.372 cgc[7] * correct ions to better match the integrated PB concentrations:
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]) Gamma = abs(surf_chrg_eq / A_surf / 1e-6)
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]) 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(); cxxSurface *surf_ptr = use.Get_surface_ptr();
bool correct_GC = surf_p->Get_correct_GC(), local_correct_GC = correct_GC; bool correct_D = surf_ptr->Get_correct_D(), local_correct_D = correct_D;
bool only_count = surf_p->Get_only_counter_ions(); bool only_count = surf_ptr->Get_only_counter_ions();
LDBLE Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6, LDBLE Gamma, cgc[12] = { 0.3805, -0.0106, 1.96, 0.812,
cgc[10] = { 0.36, 0.1721, 0.798, 0.287, 0.1457, 1.2, 0.285, 0.372 }; 0.395, 2.13,
0.380, 0.0408, 0.799, 0.594,
if (!surf_p->Donnan_factors.empty()) 0.373, 1.181 };
std::copy(surf_p->Donnan_factors.begin(), surf_p->Donnan_factors.end(), cgc); if (correct_D)
{
cgc[1] *= pow(nDbl, cgc[2]) * pow(Gamma, cgc[3]) * pow(mu_x, cgc[4]); 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; int l_iter = 0, z_iter;
sum_co = sum_counter = 0; 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; fd = surf_chrg_eq;
fd1 = 0.0; fd1 = 0.0;
z_iter = 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; l_iter = 0;
} }
std::map<LDBLE, LDBLE>::iterator it; std::map<LDBLE, LDBLE>::iterator it;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{ {
z = it->first; z = it->first;
z1 = z; Z1 = z;
if (l_iter == 0) zcorr[z_iter] = z; if (l_iter == 0) zcorr[z_iter] = z;
co_ion = surf_chrg_eq * z; co_ion = surf_chrg_eq * z;
if (!z || (only_count && co_ion > 0)) 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++; z_iter++;
continue; continue;
} }
if (nDbl && local_correct_GC) if (nDbl && local_correct_D)
{ {
/*psi_DL = fabs(p * p_psi);*/ /*psi_DL = fabs(p * p_psi);*/
if (co_ion < 0) if (co_ion < 0)
{//counter-ion {//counter-ion
if (fabs(z) > 1) temp = cgc[7]; if (fabs(z) > 1.5) temp = z2;
else temp = cgc[6]; else temp = z1;
sum_counter += z * temp; sum_counter += z * temp;
} }
else else
{// co-ion {// co-ion
if (fabs(z) > 1) temp = cgc[0] * pow(mu_x, cgc[1] * cgc[5]); if (fabs(z) > 1.5) temp = z_2;
else temp = cgc[0] * pow(mu_x, cgc[1]); else temp = z_1;
sum_co += z * temp; sum_co += z * temp;
} }
zcorr[z_iter] = z * temp; zcorr[z_iter] = z * temp;
} }
z1 = zcorr[z_iter]; Z1 = zcorr[z_iter];
eq = it->second; eq = it->second;
temp = exp(-z1 * p) * ratio_aq; temp = exp(-Z1 * p) * ratio_aq * (1 - f_free);
fd += eq * temp; fd += eq * temp;
fd1 -= z1 * eq * temp; fd1 -= Z1 * eq * temp;
if (z == 1) z1_c = z1; if (z == 1) Z1_c = Z1;
z_iter++; z_iter++;
} }
fd /= -fd1; fd /= -fd1;
@ -1152,7 +1213,7 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::
if (debug_diffuse_layer == TRUE) if (debug_diffuse_layer == TRUE)
output_msg(sformatf( output_msg(sformatf(
"iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n", "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); return (p);
} }

5
src/phreeqcpp/kinetics.cpp
View File

@ -1579,7 +1579,10 @@ set_and_run(int i, int use_mix, int use_kinetics, int nsaver,
converge = model(); converge = model();
} }
sum_species(); 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); return (converge);
} }

14
src/phreeqcpp/mainsubs.cpp
View File

@ -434,7 +434,10 @@ initial_solutions(int print)
diagonal_scale = (diag) ? TRUE : FALSE; diagonal_scale = (diag) ? TRUE : FALSE;
converge1 = check_residuals(); converge1 = check_residuals();
sum_species(); 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); add_isotopes(solution_ref);
punch_all(); punch_all();
print_all(); print_all();
@ -537,7 +540,6 @@ initial_exchangers(int print)
converge = model(); converge = model();
converge1 = check_residuals(); converge1 = check_residuals();
sum_species(); sum_species();
viscosity();
species_list_sort(); species_list_sort();
print_exchange(); print_exchange();
if (pr.user_print) if (pr.user_print)
@ -1260,12 +1262,12 @@ xsolution_save(int n_user)
temp_solution.Set_pe(solution_pe_x); temp_solution.Set_pe(solution_pe_x);
temp_solution.Set_mu(mu_x); temp_solution.Set_mu(mu_x);
temp_solution.Set_ah2o(ah2o_x); temp_solution.Set_ah2o(ah2o_x);
//temp_solution.Set_density(density_x); // 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(calc_dens()); temp_solution.Set_density(density_x);
temp_solution.Set_viscosity(this->viscos); temp_solution.Set_viscosity(viscos);
temp_solution.Set_total_h(total_h_x); temp_solution.Set_total_h(total_h_x);
temp_solution.Set_total_o(total_o_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 */ /* does include surface charge after step */
temp_solution.Set_mass_water(mass_water_aq_x); temp_solution.Set_mass_water(mass_water_aq_x);
temp_solution.Set_total_alkalinity(total_alkalinity); temp_solution.Set_total_alkalinity(total_alkalinity);

1
src/phreeqcpp/model.cpp
View File

@ -4891,7 +4891,6 @@ sum_species(void)
solution_pe_x = -s_eminus->la; solution_pe_x = -s_eminus->la;
ah2o_x = exp(s_h2o->la * LOG_10); ah2o_x = exp(s_h2o->la * LOG_10);
density_x = 1.0;
if (s_o2 != NULL) if (s_o2 != NULL)
s_o2->moles = under(s_o2->lm) * mass_water_aq_x; s_o2->moles = under(s_o2->lm) * mass_water_aq_x;
if (s_h2 != NULL) if (s_h2 != NULL)

54
src/phreeqcpp/prep.cpp
View File

@ -3996,17 +3996,16 @@ calc_PR(std::vector<class phase *> 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) * 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)); 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)); phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi));
//if (phi > 4.44) //if (phi > 4.44)
// phi = 4.44; // phi = 4.44;
} }
else else
phi = -3.0; // fugacity coefficient = 0.05 phi = -4.6; // fugacity coefficient = 0.01
//if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3) //if (!strcmp(phase_ptr->name, "H2O(g)") && phi < -4.6)
//{ //{
//// avoid such phi... //// avoid such phi...
// phi = -3; // phi = -4.6;
//} //}
phase_ptr->pr_phi = exp(phi); phase_ptr->pr_phi = exp(phi);
phase_ptr->pr_si_f = phi / LOG_10; phase_ptr->pr_si_f = phi / LOG_10;
@ -5404,6 +5403,53 @@ calc_vm(LDBLE tc, LDBLE pa)
return OK; 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:: int Phreeqc::
k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */

32
src/phreeqcpp/print.cpp
View File

@ -270,6 +270,16 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
output_msg(sformatf( output_msg(sformatf(
"\tWater in diffuse layer: %8.3e kg, %4.1f%% of total DDL-water.\n", "\tWater in diffuse layer: %8.3e kg, %4.1f%% of total DDL-water.\n",
(double) charge_ptr->Get_mass_water(), (double) d)); (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) if (use.Get_surface_ptr()->Get_debye_lengths() > 0 && d > 0)
{ {
sum_surfs = 0.0; sum_surfs = 0.0;
@ -279,8 +289,7 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
continue; continue;
cxxSurfaceCharge * charge_ptr_search = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); cxxSurfaceCharge * charge_ptr_search = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
sum_surfs += sum_surfs +=
charge_ptr_search->Get_specific_area() * charge_ptr_search->Get_specific_area() * charge_ptr_search->Get_grams();
charge_ptr_search->Get_grams();
} }
r = 0.002 * mass_water_bulk_x / sum_surfs; r = 0.002 * mass_water_bulk_x / sum_surfs;
output_msg(sformatf( output_msg(sformatf(
@ -304,10 +313,8 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
if (s_x[j]->type > HPLUS) if (s_x[j]->type > HPLUS)
continue; continue;
molality = under(s_x[j]->lm); molality = under(s_x[j]->lm);
moles_excess = mass_water_aq_x * molality * (charge_ptr->Get_g_map()[s_x[j]->z].Get_g() * moles_excess = mass_water_aq_x * molality * (charge_ptr->Get_g_map()[s_x[j]->z].Get_g() * s_x[j]->erm_ddl +
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; moles_surface = mass_water_surface * molality + moles_excess;
if (debug_diffuse_layer == TRUE) if (debug_diffuse_layer == TRUE)
{ {
@ -336,18 +343,27 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
} }
else 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; 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( output_msg(sformatf(
"\n\tTotal moles in diffuse layer (excluding water), Donnan corrected to match Poisson-Boltzmann.")); "\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 else
{
output_msg(sformatf( output_msg(sformatf(
"\n\tTotal moles in diffuse layer (excluding water), Donnan calculation.")); "\n\tTotal moles in diffuse layer (excluding water), Donnan calculation."));
output_msg(sformatf( 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", "\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)); psi_DL, exp_g));
} }
}
output_msg(sformatf("\tElement \t Moles\n")); output_msg(sformatf("\tElement \t Moles\n"));
for (j = 0; j < count_elts; j++) for (j = 0; j < count_elts; j++)
{ {

34
src/phreeqcpp/read.cpp
View File

@ -5499,9 +5499,17 @@ read_species(void)
input_error++; input_error++;
break; break;
} }
s_ptr->dw_t = 0; s_ptr->dw_a = 0; s_ptr->dw_a2 = 0; s_ptr->dw_a_visc = 0; 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, &s_ptr->dw, &s_ptr->dw_t, i = sscanf(next_char, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT,
&s_ptr->dw_a, &s_ptr->dw_a2, &s_ptr->dw_a_visc); &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; s_ptr->dw_corr = s_ptr->dw;
opt_save = OPTION_DEFAULT; opt_save = OPTION_DEFAULT;
break; break;
@ -6414,7 +6422,7 @@ read_surface(void)
if (thickness != 0) if (thickness != 0)
{ {
error_msg 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); CONTINUE);
error_msg(line_save, CONTINUE); error_msg(line_save, CONTINUE);
input_error++; input_error++;
@ -6442,12 +6450,12 @@ read_surface(void)
copy_token(token1, &next_char); copy_token(token1, &next_char);
if (token1[0] == 'T' || token1[0] == 't' || token1[0] == 'F' || token1[0] == 'f') 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; continue;
} else } else
{ {
error_msg 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); CONTINUE);
error_msg(line_save, CONTINUE); error_msg(line_save, CONTINUE);
input_error++; input_error++;
@ -6460,9 +6468,19 @@ read_surface(void)
if (j == DIGIT) if (j == DIGIT)
{ {
(void)sscanf(token1.c_str(), SCANFORMAT, &dummy); (void)sscanf(token1.c_str(), SCANFORMAT, &dummy);
if(dummy == 0)
{
dummy = 1; temp_surface.Calc_DDL_viscosity(true);
}
temp_surface.Set_DDL_viscosity(dummy); temp_surface.Set_DDL_viscosity(dummy);
continue; 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) else if (j != EMPTY)
{ {
error_msg error_msg
@ -6601,10 +6619,10 @@ read_surface(void)
i1++; i1++;
continue; continue;
} }
else if (i != EMPTY || i1 > 8) else if (i != EMPTY || i1 > 4)
{ {
error_msg 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); CONTINUE);
error_msg(line_save, CONTINUE); error_msg(line_save, CONTINUE);
input_error++; input_error++;

30
src/phreeqcpp/spread.cpp
View File

@ -162,8 +162,6 @@ read_solution_spread(void)
{ {
case 0: /* temp */ case 0: /* temp */
case 1: /* temperature */ case 1: /* temperature */
case 2: /* dens */
case 3: /* density */
case 10: /* water */ case 10: /* water */
if ((count == 2 || count == 3) && num == TRUE) if ((count == 2 || count == 3) && num == TRUE)
{ {
@ -174,6 +172,18 @@ read_solution_spread(void)
opt = OPTION_DEFAULT; opt = OPTION_DEFAULT;
} }
break; 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 6: /* ph */
case 7: /* pe */ case 7: /* pe */
if ((count == 2 || count == 3 || count == 4) if ((count == 2 || count == 3 || count == 4)
@ -285,9 +295,10 @@ read_solution_spread(void)
break; break;
case 2: /* density */ case 2: /* density */
case 3: case 3:
//sscanf(next_char, SCANFORMAT, &(soln_defaults.density));
{ {
copy_token(token, &next_char); int j = copy_token(token, &next_char);
if (j == DIGIT)
{
if (sscanf(token.c_str(), SCANFORMAT, &dummy) != 1) if (sscanf(token.c_str(), SCANFORMAT, &dummy) != 1)
{ {
error_msg("Expecting numeric value for density.", PHRQ_io::OT_CONTINUE); error_msg("Expecting numeric value for density.", PHRQ_io::OT_CONTINUE);
@ -297,22 +308,23 @@ read_solution_spread(void)
else else
{ {
soln_defaults.density = dummy; soln_defaults.density = dummy;
copy_token(token, &next_char);
if (token[0] == 'c' || token[0] == 'C')
soln_defaults.calc_density = true;
} }
int j = copy_token(token, &next_char); }
if (j != EMPTY) else if (j != EMPTY)
{ {
if (token[0] != 'c' && token[0] != 'C') 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); error_msg(line_save, PHRQ_io::OT_CONTINUE);
input_error++; input_error++;
} }
else else
{
soln_defaults.calc_density = true; soln_defaults.calc_density = true;
} }
} }
}
break; break;
case 4: /* units */ case 4: /* units */
case 8: /* unit */ case 8: /* unit */

2
src/phreeqcpp/step.cpp
View File

@ -327,6 +327,7 @@ xsolution_zero(void)
solution_pe_x = 0.0; solution_pe_x = 0.0;
mu_x = 0.0; mu_x = 0.0;
ah2o_x = 0.0; ah2o_x = 0.0;
viscos = 0.0;
density_x = 0.0; density_x = 0.0;
total_h_x = 0.0; total_h_x = 0.0;
total_o_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; solution_pe_x += solution_ptr->Get_pe() * intensive;
mu_x += solution_ptr->Get_mu() * intensive; mu_x += solution_ptr->Get_mu() * intensive;
ah2o_x += solution_ptr->Get_ah2o() * intensive; ah2o_x += solution_ptr->Get_ah2o() * intensive;
viscos += solution_ptr->Get_viscosity() * intensive;
density_x += solution_ptr->Get_density() * intensive; density_x += solution_ptr->Get_density() * intensive;
total_h_x += solution_ptr->Get_total_h() * extensive; total_h_x += solution_ptr->Get_total_h() * extensive;

1
src/phreeqcpp/structures.cpp
View File

@ -1561,6 +1561,7 @@ s_init(class species *s_ptr)
s_ptr->dw_t = 0.0; s_ptr->dw_t = 0.0;
s_ptr->dw_a = 0.0; s_ptr->dw_a = 0.0;
s_ptr->dw_a2 = 0.0; s_ptr->dw_a2 = 0.0;
s_ptr->dw_a3 = 0.0;
s_ptr->erm_ddl = 1.0; s_ptr->erm_ddl = 1.0;
s_ptr->equiv = 0; s_ptr->equiv = 0;
s_ptr->alk = 0.0; s_ptr->alk = 0.0;

226
src/phreeqcpp/transport.cpp
View File

@ -1891,6 +1891,9 @@ fill_spec(int l_cell_no, int ref_cell)
* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos * 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 * modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959
*/ */
if (print_viscosity)
viscos = Utilities::Rxn_find(Rxn_solution_map, l_cell_no)->Get_viscosity();
else
viscos = viscos_0; viscos = viscos_0;
/* /*
* put temperature factor in por_factor which corrects for porous medium... * 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. // 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. // 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, // 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; int i, icell, cp, comp;
// ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells); // ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells);
int ifirst, ilast; int ifirst, ilast;
@ -3726,36 +3730,73 @@ fill_m_s(class J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant)
void Phreeqc:: 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) 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)); // Oct. 2023, with g_i,j = exp(g*z) * SS (charge_ptr-water / aq_x)
// At filterends, concentrations of ions change step-wise to the DL. LDBLE fg_i = (1 - free_i) * g_i,
// We take the harmonic mean for f_free, the average for the DL. fg_j = (1 - free_j) * g_j;
if (ct[icell].v_m[k].z) 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) if (!g_i && g_j)
{ ct[icell].v_m[k].b_ij = b_i * free_j * b_j / (b_i + b_j) +
ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) + (b_i * (1 - free_j) + b_j * fg_j) / 4;
b_i * (1 - free_j) / 4 + b_j * g_j / 4; if (g_i && !g_j)
}
else if (g_i && !g_j)
ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_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) if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell, { /* 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... */ and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j) 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) 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 else
{ {
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1)) 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) 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) 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; 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)}. 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 22/2/18: now calculates diffusion through EDL's of multiple, differently charged surfaces
* stagnant TRUE: * 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_ij = mixf / (Dw / init_tort_f) / new_tort_f * new_por / init_por
* mixf is defined in MIX; Dw is default multicomponent diffusion coefficient; * mixf is defined in MIX; Dw is default multicomponent diffusion coefficient;
* init_tort_f equals multi_Dpor^(-multi_Dn); new_pf = new tortuosity factor. * 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; only_counter = TRUE;
ct[icell].visc1 = s_ptr1->Get_DDL_viscosity(); 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... */ /* find the immobile surface charges with DL... */
for (i = 0; i < (int)s_charge_p.size(); i++) 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; only_counter = TRUE;
ct[icell].visc2 = s_ptr2->Get_DDL_viscosity(); 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++) 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 += 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) 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) 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 += 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; 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 += 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) 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 += 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; b_i = A1 * sol_D[icell].spec[i].Dwt;
@ -5881,9 +5926,12 @@ LDBLE new_Dw, int l_cell)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
LDBLE Phreeqc:: 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. */ /* from Atkins, 1994. Physical Chemistry, 5th ed. */
//viscos = //viscos =
@ -5892,7 +5940,7 @@ viscosity(void)
// 0.000836 * (tc_x - 20) * (tc_x - 20)) / (109 + tc_x)); // 0.000836 * (tc_x - 20) * (tc_x - 20)) / (109 + tc_x));
/* Huber et al., 2009, J. Phys. Chem. Ref. Data, Vol. 38, 101-125 */ /* 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 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; int i, j, i1;
for (i = 1; i < 4; i++) for (i = 1; i < 4; i++)
denom += H[i] / pow(Tb, i); denom += H[i] / pow(Tb, i);
@ -5964,11 +6012,73 @@ viscosity(void)
both weighted by the equivalent concentration. 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 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; 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; tc = (tc_x > 200) ? 200 : tc_x;
l_water = mass_water_aq_x;
l_mu_x = mu_x;
int i1_last;
if (surf_ptr == NULL)
i1_last = 1;
else
{
i1_last = (int)surf_ptr->Get_surface_charges().size();
if (i1_last > 1)
i1_last += 1;
}
std::map <LDBLE, LDBLE> 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)
{
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
{
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;
}
for (i = 0; i < (int)this->s_x.size(); i++) for (i = 0; i < (int)this->s_x.size(); i++)
{ {
@ -5976,15 +6086,20 @@ viscosity(void)
continue; continue;
if ((lm = s_x[i]->lm) < -9) if ((lm = s_x[i]->lm) < -9)
continue; continue;
l_moles = s_x[i]->moles;
if (surf_ptr != NULL)
{
l_moles *= z_g_map.find(s_x[i]->z)->second;
}
if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3]) 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 = 0; dw_t_visc = 0;
t1 = s_x[i]->moles / mass_water_aq_x; t1 = l_moles / l_water;
l_z = fabs(s_x[i]->z); l_z = fabs(s_x[i]->z);
if (l_z) if (l_z)
f_z = (l_z * l_z + l_z) / 2; f_z = (l_z * l_z + l_z) / 2;
else else
f_z = mu_x / t1; 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 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]) if (s_x[i]->Jones_Dole[7] || s_x[i]->Jones_Dole[8])
{ {
@ -5994,20 +6109,22 @@ viscosity(void)
s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0); s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0);
} }
// find B * m and D * m * mu^d3 // find B * m and D * m * mu^d3
s_x[i]->dw_t_visc = (s_x[i]->Jones_Dole[0] + dw_t_visc = (s_x[i]->Jones_Dole[0] +
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1; s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
Bc += s_x[i]->dw_t_visc; Bc += dw_t_visc;
// define f_I from the exponent of the D * m^d3 term... // define f_I from the exponent of the D * m^d3 term...
if (s_x[i]->Jones_Dole[5] >= 1) if (s_x[i]->Jones_Dole[5] >= 1)
t2 = mu_x / 3 / s_x[i]->Jones_Dole[5]; t2 = l_mu_x / 3 / s_x[i]->Jones_Dole[5];
else if (s_x[i]->Jones_Dole[5] > 0.4) else if (s_x[i]->Jones_Dole[5] > 0.4)
t2 = -0.8 / s_x[i]->Jones_Dole[5]; t2 = -0.8 / s_x[i]->Jones_Dole[5];
else else
t2 = -1; t2 = -1;
t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) * 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); 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);
s_x[i]->dw_t_visc += t3; if (t3 < -1e-5) // add this check
t3 = 0;
Dc += t3; 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 )); //output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc ));
} }
// parms for A... // parms for A...
@ -6026,19 +6143,19 @@ viscosity(void)
// volumina for f_an... // volumina for f_an...
{ {
V_Cl = s_x[i]->logk[vm_tc]; V_Cl = s_x[i]->logk[vm_tc];
V_an += V_Cl * s_x[i]->moles; V_an += V_Cl * l_moles;
m_an += s_x[i]->moles; m_an += l_moles;
} }
else// if (s_x[i]->Jones_Dole[6]) else// if (s_x[i]->Jones_Dole[6])
{ {
V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * s_x[i]->moles; V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * l_moles;
m_an += s_x[i]->moles; m_an += l_moles;
} }
if (Dw) if (Dw)
{ {
// anions for A... // anions for A...
m_min += s_x[i]->moles; m_min += l_moles;
t1 = s_x[i]->moles * l_z; t1 = l_moles * l_z;
eq_min -= t1; eq_min -= t1;
eq_dw_min -= t1 / Dw; eq_dw_min -= t1 / Dw;
} }
@ -6046,8 +6163,8 @@ viscosity(void)
else if (Dw) else if (Dw)
{ {
// cations for A... // cations for A...
m_plus += s_x[i]->moles; m_plus += l_moles;
t1 = s_x[i]->moles * l_z; t1 = l_moles * l_z;
eq_plus += t1; eq_plus += t1;
eq_dw_plus += t1 / Dw; eq_dw_plus += t1 / Dw;
} }
@ -6064,7 +6181,7 @@ viscosity(void)
} }
else else
A = 0; A = 0;
viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / mass_water_aq_x); viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / l_water);
if (m_an) if (m_an)
V_an /= m_an; V_an /= m_an;
if (!V_Cl) if (!V_Cl)
@ -6073,16 +6190,39 @@ viscosity(void)
fan = 2 - V_an / V_Cl; fan = 2 - V_an / V_Cl;
//else //else
// fan = 1; // fan = 1;
if (Dc < 0)
Dc = 0; // provisional...
viscos += viscos_0 * fan * (Bc + Dc); viscos += viscos_0 * fan * (Bc + Dc);
if (viscos < 0) if (viscos < 0)
{ {
viscos = viscos_0; viscos = viscos_0;
warning_msg("viscosity < 0, reset to viscosity of pure water"); warning_msg("viscosity < 0, reset to viscosity of pure water");
} }
if (!surf_ptr)
{
for (i = 0; i < (int)this->s_x.size(); i++) 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); 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; return viscos;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */