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

Tony's fixes for Pitzer, High Pressure is still not right,

Browse Source 
will require more fitting.

git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/trunk@7854 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2013-07-23 14:51:26 +00:00
parent 7da9c7b1e3
commit 79698b010d
4 changed files with 19 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Phreeqc.h
View File

@ -1677,7 +1677,7 @@ protected:
int count_sys, max_sys;
LDBLE sys_tot;
LDBLE V_solutes, rho_0, kappa_0, p_sat/*, ah2o_x0*/;
LDBLE V_solutes, rho_0, rho_0_sat, kappa_0, p_sat/*, ah2o_x0*/;
LDBLE eps_r; // relative dielectric permittivity
LDBLE DH_A, DH_B, DH_Av; // Debye-Hueckel A, B and Av
LDBLE QBrn; // Born function d(ln(eps_r))/dP / eps_r * 41.84004, for supcrt calc'n of molal volume

27
pitzer.cpp
View File

@ -738,43 +738,18 @@ C
C SUBROUTINE TO CALUCLATE TEMPERATURE DEPENDENCE OF PITZER PARAMETER
C
*/
LDBLE DC0;
int i;
LDBLE TR = 298.15;
if (fabs(TK - OTEMP) < 0.001e0 && fabs(patm_x - OPRESS) < 0.1)
if (fabs(TK - OTEMP) < 0.001 && fabs(patm_x - OPRESS) < 0.1)
return OK;
/*
C Set DW0
*/
//DW(TK);
//rho_0 = DW0;
//patm_x = VP;
DW0 = rho_0 = calc_rho_0(TK - 273.15, patm_x);
VP = patm_x;
for (i = 0; i < count_pitz_param; i++)
{
calc_pitz_param(pitz_params[i], TK, TR);
}
//DC0 = DC(TK);
calc_dielectrics(TK - 273.15, patm_x);
DC0 = eps_r;
// Only at 1.0 atm??? if (fabs(TK - TR) < 0.001 && fabs(patm_x - OPRESS) < 0.1)
if (fabs(TK - TR) < 0.001 && fabs(patm_x - 1.0) < 0.1)
{
A0 = 0.392e0;
//Cite Jonathan Toner
//A0 = 86.6836498e0 + (TK)*0.084879594e0 - (0.0000888785e0)*pow(TK,(LDBLE) 2.0e0) + (4.88096e-8)*pow(TK,(LDBLE) 3.0e0) - (1327.31477e0/(TK)) - 17.6460172e0*log(TK);
}
else
{
//DC0 = DC(TK);
calc_dielectrics(TK - 273.15, patm_x);
DC0 = eps_r;
A0 = 1.400684e6 * sqrt(DW0 / (pow((DC0 * TK), (LDBLE) 3.0e0)));
/*A0=1.400684D6*(DW0/(DC0*TK)**3.0D0)**0.5D0 */
}
OTEMP = TK;
OPRESS = patm_x;
return OK;

18
sit.cpp
View File

@ -1331,7 +1331,6 @@ C
C SUBROUTINE TO CALUCLATE TEMPERATURE DEPENDENCE OF PITZER PARAMETER
C
*/
LDBLE DC0;
int i;
LDBLE TR = 298.15;
@ -1339,23 +1338,14 @@ C
/*
C Set DW0
*/
DW(TK);
DW0 = rho_0 = calc_rho_0(TK - 273.15, patm_x);
VP = patm_x;
for (i = 0; i < count_sit_param; i++)
{
calc_sit_param(sit_params[i], TK, TR);
}
DC0 = DC(TK);
// Only at 1.0 atm??? if (fabs(TK - TR) < 0.001 && fabs(patm_x - OPRESS) < 0.1)
if (fabs(TK - TR) < 0.001 && fabs(patm_x - 1.0) < 0.1)
{
sit_A0 = 0.392e0;
}
else
{
DC0 = DC(TK);
sit_A0 = 1.400684e6 * sqrt(DW0 / (pow((DC0 * TK), (LDBLE) 3.0e0)));
/*sit_A0=1.400684D6*(DW0/(DC0*TK)**3.0D0)**0.5D0 */
}
calc_dielectrics(TK - 273.15, patm_x);
sit_A0 = A0;
OTEMP = TK;
OPRESS = patm_x;
return OK;

20
utilities.cpp
View File

@ -108,7 +108,7 @@ calc_rho_0(LDBLE tc, LDBLE pa)
LDBLE Tc = 647.096, th = 1 - T / Tc;
LDBLE b1 = 1.99274064, b2 = 1.09965342, b3 = -0.510839303,
b4 = -1.75493479, b5 = -45.5170352, b6 = -6.7469445e5;
rho_0 = 322.0 * (1.0 + b1 * pow(th, (LDBLE) 1./3.) + b2 * pow(th, (LDBLE) 2./3.) + b3 * pow(th, (LDBLE) 5./3.) +\
rho_0_sat = 322.0 * (1.0 + b1 * pow(th, (LDBLE) 1./3.) + b2 * pow(th, (LDBLE) 2./3.) + b3 * pow(th, (LDBLE) 5./3.) +\
b4 * pow(th, (LDBLE) 16./3.) + b5 * pow(th, (LDBLE) 43./3.) + b6 * pow(th, (LDBLE) 110./3));
//pressure...
LDBLE p0 = 5.1880000E-02 + tc * (-4.1885519E-04 + tc * ( 6.6780748E-06 + tc * (-3.6648699E-08 + tc * 8.3501912E-11)));
@ -128,7 +128,7 @@ calc_rho_0(LDBLE tc, LDBLE pa)
if (!use.Get_gas_phase_in())
patm_x = pa;
pa -= (p_sat - 1e-6);
rho_0 += pa * (p0 + pa * (p1 + pa * (p2 + sqrt(pa) * p3)));
rho_0 = rho_0_sat + pa * (p0 + pa * (p1 + pa * (p2 + sqrt(pa) * p3)));
if (rho_0 < 0.01)
rho_0 = 0.01;
@ -159,8 +159,8 @@ calc_dielectrics(LDBLE tc, LDBLE pa)
LDBLE d1000 = u1 * exp(T * (u2 + T * u3)); // relative dielectric constant at 1000 bar
LDBLE c = u4 + u5 / (u6 + T);
LDBLE b = u7 + u8 / T + u9 * T;
pa *= 1.01325; // pa in bar
eps_r = d1000 + c * log((b + pa) / (b + 1e3)); // relative dielectric constant
LDBLE pb = pa * 1.01325; // pa in bar
eps_r = d1000 + c * log((b + pb) / (b + 1e3)); // relative dielectric constant
/* qe^2 / (eps_r * kB * T) = 4.803204e-10**2 / 1.38065e-16 / (eps_r * T)
= 1.671008e-3 (esu^2 / (erg/K)) / (eps_r * T) */
@ -170,14 +170,20 @@ calc_dielectrics(LDBLE tc, LDBLE pa)
DH_A = DH_B * e2_DkT / (2. * LOG_10); //(mol/kg)^-0.5
/* A0 in pitzer */
if (pitzer_model || sit_model)
{
A0 = DH_B * e2_DkT / 6.0;
}
/* Debye-Hueckel limiting slope = DH_B * e2_DkT * RT * (d(ln(eps_r)) / d(P) - compressibility) */
DH_Av = DH_B * e2_DkT * R_LITER_ATM * 1e3 * T * (c / (b + pa) * 1.01325 / eps_r - kappa_0 / 3.); // (cm3/mol)(mol/kg)^-0.5
DH_Av = DH_B * e2_DkT * R_LITER_ATM * 1e3 * T * (c / (b + pb) * 1.01325 / eps_r - kappa_0 / 3.); // (cm3/mol)(mol/kg)^-0.5
DH_B /= 1e8; // kappa, 1/Angstrom(mol/kg)^-0.5
/* the Born functions, * 41.84 to give molal volumes in cm3/mol... */
ZBrn = (- 1 / eps_r + 1.0) * 41.84004;
QBrn = c / (b + pa) / eps_r / eps_r * 41.84004;
QBrn = c / (b + pb) / eps_r / eps_r * 41.84004;
/* dgdP from subroutine gShok2 in supcrt92, g is neglected here (at tc < 300)...
and, dgdP is small. Better, adapt Wref to experimental Vm's */
dgdP = 0;
@ -192,7 +198,7 @@ calc_dielectrics(LDBLE tc, LDBLE pa)
// dgdP = - sa * sb * pow(1.0 - rho_0, sb - 1.0) * rho_0 * kappa_0 / 1.01325;
// LDBLE ft = pow((tc - 155.0)/300.0, 4.8) + csc[1] * pow((tc - 155.0)/300.0, 16.0);
// LDBLE dfdP = ft * (-3.0 * csc[2] * pow(1000.0 - pa, 2) - 4.0 * csc[3] * pow(1000.0 - pa, 3));
// LDBLE dfdP = ft * (-3.0 * csc[2] * pow(1000.0 - pb, 2) - 4.0 * csc[3] * pow(1000.0 - pb, 3));
// dgdP -= dfdP;
//}