diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index 86c60359..058b63ed 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -3380,20 +3380,22 @@ factor(struct LOC_exec * LINK) case toksetdiff_c: { - double d; + double d, d_v_d = 0; require(toklp, LINK); const char* str = stringfactor(STR1, LINK); require(tokcomma, LINK); - // double arugument d = realexpr(LINK); + if (LINK->t != NULL && LINK->t->kind == tokcomma) + { + LINK->t = LINK->t->next; + d_v_d = realexpr(LINK); + } require(tokrp, LINK); - n.UU.val = (parse_all) ? 1 : PhreeqcPtr->setdiff_c(str, d); - - //PhreeqcPtr->PHRQ_free((void *) str); + n.UU.val = (parse_all) ? 1 : PhreeqcPtr->setdiff_c(str, d, d_v_d); } break; diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 7a1e7e9f..fec912d7 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -110,7 +110,7 @@ public: LDBLE aqueous_vm(const char* species_name); LDBLE phase_vm(const char* phase_name); LDBLE diff_c(const char* species_name); - LDBLE setdiff_c(const char* species_name, double d); + LDBLE setdiff_c(const char * species_name, double d, double d_v_d); LDBLE flux_mcd(const char* species_name, int option); LDBLE sa_declercq(double type, double sa, double d, double m, double m0, double gfw); LDBLE calc_SC(void); diff --git a/phreeqcpp/basicsubs.cpp b/phreeqcpp/basicsubs.cpp index 9f4116b8..3ef3d6c5 100644 --- a/phreeqcpp/basicsubs.cpp +++ b/phreeqcpp/basicsubs.cpp @@ -178,7 +178,6 @@ sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw error_msg(error_string, CONTINUE); input_error++; return (MISSING); - } /* ---------------------------------------------------------------------- */ @@ -187,97 +186,58 @@ diff_c(const char *species_name) /* ---------------------------------------------------------------------- */ { class species *s_ptr; - LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av; - sqrt_mu = sqrt(mu_x); + LDBLE Dw; s_ptr = s_search(species_name); if (s_ptr == NULL) return(0); if ((Dw = s_ptr->dw) == 0) - { - if (correct_Dw) - Dw = default_Dw; - else - return(0); - } + return(0); if (correct_Dw) { - if ((l_z = fabs(s_ptr->z)) == 0) - { - //l_z = 1; // only a 1st approximation for correct_Dw in electrical field - } - else - { - 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; - } + calc_SC(); + Dw = s_ptr->dw_corr; } - if (tk_x != 298.15 && s_ptr->dw_t) - Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + else + { + 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 * viscos_0_25 / viscos_0; - return s_ptr->dw_corr; + Dw *= viscos_0_25 / viscos_0; + } + if (s_ptr->dw_a_v_dif) + Dw *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif); + return Dw; } + /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -setdiff_c(const char *species_name, double d) +setdiff_c(const char *species_name, double d, double d_v_d) /* ---------------------------------------------------------------------- */ { class species *s_ptr; - LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av; - sqrt_mu = sqrt(mu_x); + LDBLE Dw; s_ptr = s_search(species_name); if (s_ptr == NULL) return(0); Dw = s_ptr->dw = d; + s_ptr->dw_a_v_dif = d_v_d; if (correct_Dw) { - if ((l_z = fabs(s_ptr->z)) == 0) - { - //l_z = 1; // only a 1st approximation for correct_Dw in electrical field - } - else - { - 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; - } + calc_SC(); + Dw = s_ptr->dw_corr; } - if (tk_x != 298.15 && s_ptr->dw_t) + else + { + 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 * viscos_0_25 / viscos_0; - return s_ptr->dw_corr; + Dw *= viscos_0_25 / viscos_0; + } + if (d_v_d) + Dw *= pow(viscos_0 / viscos, d_v_d); + return Dw; } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: @@ -287,6 +247,7 @@ calc_SC(void) class species *s_ptr; int i; LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av, v_Cl = 1; + SC = 0; sqrt_mu = sqrt(mu_x); bool Falk = false; s_ptr = s_search("H+"); @@ -294,7 +255,6 @@ calc_SC(void) return(0); else if (s_ptr->dw_a3 > 24) Falk = true; - SC = 0; //LDBLE ta1, ta2, ta3, ta4; //for (i = 0; i < (int)this->s_x.size(); i++) //{ @@ -320,7 +280,10 @@ calc_SC(void) if (correct_Dw) Dw = default_Dw; else + { + s_x[i]->dw_corr = 0; continue; + } } if (s_x[i]->lm < min_dif_LM) continue; @@ -343,27 +306,6 @@ calc_SC(void) } else { - //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) { @@ -384,6 +326,8 @@ calc_SC(void) } Dw *= ff; + if (correct_Dw) + 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; } @@ -416,8 +360,13 @@ calc_SC(void) 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 (correct_Dw) + Dw = default_Dw; // or charge based...Dw = l_z > 0 ? 1.6e-9 / l_z : 2e-9 / -l_z; + else + { + s_x[i]->dw_corr = 0; + 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); @@ -496,6 +445,8 @@ calc_SC(void) ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka)); Dw *= ff; + if (correct_Dw) + 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 * 1e3 * Dw_SC; } @@ -540,7 +491,9 @@ calc_SC(void) //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); + t1 *= pow(viscos_0 / viscos, av); + if (correct_Dw) + s_x[i]->dw_corr *= t1 / Dw; // fractional contribution in mu, and correct for charge imbalance a2 = 2 / (eq_plus + eq_min); diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index e9d82f9d..50aae97a 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -14,7 +14,7 @@ # define NAN nan("1") # endif #endif -#define MISSING -9999.999 +#define MISSING -9999.999 #include "NA.h" /* NA = not available */ #define F_C_MOL 96493.5 /* C/mol or joule/volt-eq */ @@ -712,18 +712,18 @@ public: secondary = NULL; gfw = 0; // gram formula wt of species z = 0; // charge of species - // tracer diffusion coefficient in water at 25oC, m2/s - dw = 0; - // correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15) - dw_t = 0; + dw = 0; // tracer diffusion coefficient in water at 25oC, m2/s + dw_t = 0; // correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15) // 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)) + // with DHO: ka = DH_B * dw_a * (1 + DD(V_apparent)^dw_a2 * sqrt_mu, dw_a3 is a switch, see calc_SC in PBasic dw_a = 0; dw_a2 = 0; dw_a3 = 0; - dw_a_visc = 0; // viscosity correction of SC + dw_a_visc = 0; // exponent in viscosity correction of SC + dw_a_v_dif = 0; // exponent in viscosity correction of D, the diffusion coefficient of the species dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC dw_t_visc = 0; // contribution to viscosity - dw_corr = 0; // dw corrected for TK and mu + dw_corr = 0; // dw corrected for mu and TK erm_ddl = 0; // enrichment factor in DDL equiv = 0; // equivalents in exchange species alk = 0; // alkalinity of species, used for cec in exchange @@ -784,6 +784,7 @@ public: LDBLE dw_a2; LDBLE dw_a3; LDBLE dw_a_visc; + LDBLE dw_a_v_dif; LDBLE dw_t_SC; LDBLE dw_t_visc; LDBLE dw_corr; @@ -1489,6 +1490,8 @@ public: Dwt = 0; // temperature factor for Dw dw_t = 0; + // viscosity factor for Dw + dw_a_v_dif = 0; // enrichment factor in ddl erm_ddl = 0; } @@ -1502,6 +1505,7 @@ public: LDBLE z; LDBLE Dwt; LDBLE dw_t; + LDBLE dw_a_v_dif; LDBLE erm_ddl; }; @@ -1517,7 +1521,9 @@ public: count_exch_spec = 0; // total moles of X-, max X- in transport step in sol_D[1], tk exch_total = 0, x_max = 0, tk_x = 0; - // (tk_x * viscos_0_25) / (298 * viscos) + // (tk_x * viscos_0_25) / (298 * viscos_0) + viscos_f0 = 0; + // (viscos_0) / (298 * viscos) viscos_f = 0; spec = NULL; spec_size = 0; @@ -1525,7 +1531,7 @@ public: int count_spec; int count_exch_spec; LDBLE exch_total, x_max, tk_x; - LDBLE viscos_f; + LDBLE viscos_f0, viscos_f; class spec* spec; int spec_size; }; diff --git a/phreeqcpp/read.cpp b/phreeqcpp/read.cpp index a5eed844..24f6c8c5 100644 --- a/phreeqcpp/read.cpp +++ b/phreeqcpp/read.cpp @@ -5499,9 +5499,9 @@ read_species(void) input_error++; break; } - 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); + 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; s_ptr->dw_a_v_dif = 0; + i = sscanf(next_char, SCANFORMAT 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, &s_ptr->dw_a_v_dif); if (i < 1) { input_error++; diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 31860db9..d6b0535c 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -115,6 +115,7 @@ transport(void) sol_D[i].count_exch_spec = 0; sol_D[i].exch_total = 0; sol_D[i].x_max = 0; + sol_D[i].viscos_f0 = 1.0; sol_D[i].viscos_f = 1.0; sol_D[i].tk_x = 298.15; sol_D[i].spec = NULL; @@ -1640,14 +1641,14 @@ init_heat_mix(int l_nmix) { if (implicit) { - LDBLE viscos_f; + LDBLE viscos_f0; l_heat_nmix = l_nmix; for (i = 1; i <= count_cells + 1; i++) { heat_mix_array[i - 1] = heat_mix_array[i] / l_heat_nmix; /* for implicit, m[i] has mixf with higher cell */ - viscos_f = sol_D[i - 1].viscos_f * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15); - viscos_f += sol_D[i].viscos_f * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15); - heat_mix_array[i - 1] *= (viscos_f / 2); + viscos_f0 = sol_D[i - 1].viscos_f0 * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15); + viscos_f0 += sol_D[i].viscos_f0 * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15); + heat_mix_array[i - 1] *= (viscos_f0 / 2); } } else @@ -1821,7 +1822,7 @@ fill_spec(int l_cell_no, int ref_cell) class master *master_ptr; LDBLE dum, dum2; LDBLE lm; - LDBLE por, por_il, viscos_f, viscos_il_f, viscos; + LDBLE por, por_il, viscos_f0, viscos_f, viscos_il_f0, viscos_il_f, viscos; bool x_max_done = false; std::set loc_spec_names; @@ -1860,13 +1861,14 @@ fill_spec(int l_cell_no, int ref_cell) sol_D[l_cell_no].spec[i].z = 0.0; sol_D[l_cell_no].spec[i].Dwt = 0.0; sol_D[l_cell_no].spec[i].dw_t = 0.0; + sol_D[l_cell_no].spec[i].dw_a_v_dif = 0.0; sol_D[l_cell_no].spec[i].erm_ddl = 0.0; sol_D[l_cell_no].count_exch_spec = sol_D[l_cell_no].count_spec = 0; } sol_D[l_cell_no].tk_x = tk_x; - viscos_f = viscos_il_f = 1.0; + viscos_f0 = viscos_il_f0 = viscos_f = viscos_il_f = 1.0; if (l_cell_no == 0) { por = cell_data[1].por; @@ -1883,10 +1885,10 @@ fill_spec(int l_cell_no, int ref_cell) por_il = cell_data[l_cell_no].por_il; } if (por < multi_Dpor_lim) - por = viscos_f = 0.0; + por = viscos_f0 = viscos_f = 0.0; if (por_il < interlayer_Dpor_lim) - por_il = viscos_il_f = 0.0; + por_il = viscos_il_f0 = viscos_il_f = 0.0; /* * 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 @@ -1898,10 +1900,14 @@ fill_spec(int l_cell_no, int ref_cell) /* * put temperature factor in por_factor which corrects for porous medium... */ - dum = viscos_0_25 / viscos; - viscos_f *= dum; - viscos_il_f *= dum; - sol_D[l_cell_no].viscos_f = dum; + dum = viscos_0_25 / viscos_0; + dum2 = viscos_0 / viscos; + viscos_f0 *= dum; + viscos_il_f0 *= dum; + sol_D[l_cell_no].viscos_f0 = dum; + viscos_f *= dum2; + viscos_il_f *= dum2; + sol_D[l_cell_no].viscos_f = dum2; count_spec = count_exch_spec = 0; /* @@ -2002,6 +2008,10 @@ fill_spec(int l_cell_no, int ref_cell) else sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr2->dw * viscos_il_f; } + if (s_ptr2->dw_a_v_dif) + sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = s_ptr2->dw_a_v_dif; + else + sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = 0.0; //if (implicit) // && name_ret.second && (l_cell_no > 1 || (l_cell_no == 1 && bcon_first != 2))) //{ @@ -2089,22 +2099,29 @@ fill_spec(int l_cell_no, int ref_cell) sol_D[l_cell_no].spec[count_spec].lg = s_ptr->lg; sol_D[l_cell_no].spec[count_spec].z = s_ptr->z; if (s_ptr->dw == 0) - sol_D[l_cell_no].spec[count_spec].Dwt = default_Dw * viscos_f; + sol_D[l_cell_no].spec[count_spec].Dwt = default_Dw; else { if (s_ptr->dw_t) { sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw * - exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15) * viscos_f; + exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); sol_D[l_cell_no].spec[count_spec].dw_t = s_ptr->dw_t; } else - sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw * viscos_f; + sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw; } + if (s_ptr->dw_a_v_dif) + { + sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = s_ptr->dw_a_v_dif; + sol_D[l_cell_no].spec[count_spec].Dwt *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif); + } + else + sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = 0.0; if (correct_Dw) { //calc_SC(); // removed that neutral species are corrected as if z = 1, but is viscosity-dependent - sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr * viscos_f; + sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr; } if (l_cell_no <= count_cells + 1 && sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn) > diffc_max) diffc_max = sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn); @@ -4272,7 +4289,11 @@ 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 / ct[icell].visc1; + dum = ct[icell].visc1; + if (sol_D[icell].spec[i].dw_a_v_dif) + dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif); + g_i *= sol_D[icell].spec[i].erm_ddl / dum; + //g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4295,7 +4316,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } } } - g_j *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc2; + dum = ct[icell].visc2; + if (sol_D[icell].spec[i].dw_a_v_dif) + dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif); + g_j *= sol_D[icell].spec[i].erm_ddl / dum; + //g_j *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc2; } } @@ -4305,11 +4330,13 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) b_j *= sol_D[icell].spec[i].Dwt; else { - dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f; + dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f0; dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x); - dum2 *= sol_D[jcell].viscos_f; + dum2 *= sol_D[jcell].viscos_f0; b_j *= dum2; } + if (sol_D[icell].spec[i].dw_a_v_dif) + b_j *= pow(sol_D[jcell].viscos_f / sol_D[icell].viscos_f, sol_D[icell].spec[i].dw_a_v_dif); calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); k++; @@ -4388,7 +4415,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } } } - g_i *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc1; + dum = ct[icell].visc1; + if (sol_D[jcell].spec[j].dw_a_v_dif) + dum = pow(dum, sol_D[icell].spec[j].dw_a_v_dif); + g_i *= sol_D[jcell].spec[j].erm_ddl / dum; + //g_i *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4396,7 +4427,12 @@ 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 / ct[icell].visc2; + + dum = ct[icell].visc2; + if (sol_D[jcell].spec[j].dw_a_v_dif) + dum = pow(dum, sol_D[jcell].spec[j].dw_a_v_dif); + g_j *= sol_D[jcell].spec[j].erm_ddl / dum; + //g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2; } } b_i = A1; @@ -4405,11 +4441,14 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) b_i *= sol_D[jcell].spec[j].Dwt; else { - dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f; + dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f0; dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x); - dum2 *= sol_D[icell].viscos_f; + dum2 *= sol_D[icell].viscos_f0; b_i *= dum2; } + if (sol_D[icell].spec[i].dw_a_v_dif) + b_i *= pow(sol_D[icell].viscos_f / sol_D[jcell].viscos_f, sol_D[icell].spec[i].dw_a_v_dif); + calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); k++; @@ -4482,7 +4521,11 @@ 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 / ct[icell].visc1; + dum = ct[icell].visc1; + if (sol_D[icell].spec[i].dw_a_v_dif) + dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif); + g_i *= sol_D[icell].spec[i].erm_ddl / dum; + //g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1; } if (dl_aq2) { @@ -4490,7 +4533,11 @@ 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 / ct[icell].visc2; + dum = ct[icell].visc2; + if (sol_D[jcell].spec[j].dw_a_v_dif) + dum = pow(dum, sol_D[jcell].spec[j].dw_a_v_dif); + g_j *= sol_D[jcell].spec[j].erm_ddl / dum; + //g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2; } } b_i = A1 * sol_D[icell].spec[i].Dwt; @@ -4874,12 +4921,11 @@ Step (from cell 1 to count_cells + 1): LDBLE lav, A_ij, por, Dp1, Dp2; cxxMix * mix_ptr; cxxSurface *surface_ptr1, *surface_ptr2; - LDBLE viscos_f; + LDBLE viscos_f0; /* * temperature and viscosity correction for MCD coefficient, D_T = D_298 * viscos_298 / viscos */ - viscos_f = viscos_0; - viscos_f = viscos_0_25 / viscos_f; + viscos_f0 = viscos_0_25 / viscos_0; //n1 = 0; //n2 = n1 + 1; @@ -5089,7 +5135,7 @@ Step (from cell 1 to count_cells + 1): Dp2 = (Dp2 + Dp1) / 2; /* and the mixing factor... */ - mcd_mixf = Dp2 * viscos_f * A_ij * DDt / lav; + mcd_mixf = Dp2 * viscos_f0 * A_ij * DDt / lav; } mixf = mixf_store + mcd_mixf; @@ -5165,7 +5211,7 @@ Step (from cell 1 to count_cells + 1): Dp1 = (Dp1 + Dp2) / 2; /* and the mixing factor... */ - mcd_mixf = Dp1 * viscos_f * A_ij * DDt / lav; + mcd_mixf = Dp1 * viscos_f0 * A_ij * DDt / lav; } mixf = mixf_store + mcd_mixf; @@ -5545,12 +5591,11 @@ diff_stag_surf(int mobile_cell) LDBLE Dp1, Dp2; cxxMix *mix_ptr; cxxSurface *surface_ptr1, *surface_ptr2; - LDBLE viscos_f; + LDBLE viscos_f0; /* * temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos) */ - viscos_f = viscos_0; - viscos_f = viscos_0_25 / viscos_f; + viscos_f0 = viscos_0_25 / viscos_0; cxxSurface surface_n1, surface_n2; cxxSurface *surface_n1_ptr = &surface_n1; @@ -5698,7 +5743,7 @@ diff_stag_surf(int mobile_cell) if (multi_Dflag) { Dp2 = comp_k_ptr->Get_Dw() * - pow(cell_data[i2].por, multi_Dn) * viscos_f; + pow(cell_data[i2].por, multi_Dn) * viscos_f0; Dp1 = 0; if (surf1) { @@ -5708,7 +5753,7 @@ diff_stag_surf(int mobile_cell) if (strcmp(comp_k1_ptr->Get_formula().c_str(), comp_k_ptr->Get_formula().c_str()) != 0) continue; - Dp1 = comp_k1_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f; + Dp1 = comp_k1_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f0; break; } } @@ -5757,7 +5802,7 @@ diff_stag_surf(int mobile_cell) /* find diffusion coefficients of surfaces... */ if (multi_Dflag) { - Dp1 = comp_k_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f; + Dp1 = comp_k_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f0; Dp2 = 0; if (surf2) @@ -5768,7 +5813,7 @@ diff_stag_surf(int mobile_cell) if (strcmp(comp_k1_ptr->Get_formula().c_str(), comp_k_ptr->Get_formula().c_str()) != 0) continue; - Dp2 = comp_k1_ptr->Get_Dw() * pow(cell_data[i2].por, multi_Dn) * viscos_f; + Dp2 = comp_k1_ptr->Get_Dw() * pow(cell_data[i2].por, multi_Dn) * viscos_f0; break; } } @@ -6121,7 +6166,7 @@ viscosity(cxxSurface *surf_ptr) 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 + if (t3 < -1e-5) t3 = 0; Dc += t3; if (!surf_ptr) s_x[i]->dw_t_visc = dw_t_visc + t3; @@ -6186,10 +6231,10 @@ viscosity(cxxSurface *surf_ptr) V_an /= m_an; if (!V_Cl) V_Cl = calc_vm_Cl(); - if (V_an) + if (V_an && V_Cl) fan = 2 - V_an / V_Cl; - //else - // fan = 1; + else + fan = 1; if (Dc < 0) Dc = 0; // provisional... viscos += viscos_0 * fan * (Bc + Dc);