diff --git a/src/thameslib/ChemicalSystem.h b/src/thameslib/ChemicalSystem.h
index f7a88d6..c817117 100644
--- a/src/thameslib/ChemicalSystem.h
+++ b/src/thameslib/ChemicalSystem.h
@@ -1467,18 +1467,18 @@ unsigned int getDCId (const string &dcname)
if (p != DCIdLookup_.end()) {
return p->second;
} else {
- cout << "WARNING: Could not find DCIdLookup_ match to " << dcname << endl;
- cout << "WARNING: Here are the ones I know about:" << endl;
- cout.flush();
- p = DCIdLookup_.begin();
- while (p != DCIdLookup_.end()) {
- cout << "WARNING: " << p->first << " ("
- << p->second << ")" << endl;
- cout.flush();
- p++;
- }
- cout << "WARNING:" << endl;
- cout.flush();
+ // cout << "WARNING: Could not find DCIdLookup_ match to " << dcname << endl;
+ // cout << "WARNING: Here are the ones I know about:" << endl;
+ // cout.flush();
+ // p = DCIdLookup_.begin();
+ // while (p != DCIdLookup_.end()) {
+ // cout << "WARNING: " << p->first << " ("
+ // << p->second << ")" << endl;
+ // cout.flush();
+ // p++;
+ // }
+ // cout << "WARNING:" << endl;
+ // cout.flush();
return (numDCs_ + 9999);
}
}
@@ -4851,8 +4851,6 @@ double getDCChemicalPotential (const long int dcidx, bool norm = false)
/**
@brief Get the chemical activity of a dependent component (DC).
-@note NOT USED.
-
@param dcidx is the index of the DC being queried
@return the chemical activity of the DC
*/
@@ -4861,10 +4859,25 @@ double getDCActivity (const long int dcidx)
return (node_->Get_aDC(dcidx));
}
+/**
+@brief Get the chemical activity of a dependent component (DC) by name
+
+@param dcname is the name of the DC being queried
+@return the chemical activity of the DC
+*/
+double getDCActivity (const string &dcname)
+{
+ int dcidx = getDCId(dcname);
+ if (dcidx < numDCs_) {
+ return (node_->Get_aDC(dcidx));
+ }
+ return(0.0);
+}
+
/**
@brief Get the concentration of a dependent component (DC).
-The units of the returned value depend on the typd of DC being queried:
+The units of the returned value depend on the type of DC being queried:
- Aqueous species [molal]
- Gas species [partial pressure]
@@ -4881,6 +4894,30 @@ double getDCConcentration (const long int dcidx)
return (node_->Get_cDC(dcidx));
}
+/**
+@brief Get the concentration of a dependent component (DC) by name
+
+The units of the returned value depend on the type of DC being queried:
+
+ - Aqueous species [molal]
+ - Gas species [partial pressure]
+ - Surface complexes [mol/m2]
+ - Species in other phases [mole fraction]
+
+Will return a value of zero if the DC does not exist
+
+@param dcname is the name of the DC component begin queried
+@return the concentration of the DC in appropriate units
+*/
+double getDCConcentration (const string &dcname)
+{
+ int dcidx = getDCId(dcname);
+ if (dcidx < numDCs_) {
+ return (node_->Get_cDC(dcidx));
+ }
+ return(0.0);
+}
+
/**
@brief Set the red channel in the rgb triplet for color of a microstructure phase.
diff --git a/src/thameslib/KineticModel.cc b/src/thameslib/KineticModel.cc
index e5d4682..ad0728f 100644
--- a/src/thameslib/KineticModel.cc
+++ b/src/thameslib/KineticModel.cc
@@ -34,6 +34,7 @@ KineticModel::KineticModel ()
name_.clear();
isKinetic_.clear();
+ isPK_.clear();
isThermo_.clear();
isSoluble_.clear();
microPhaseId_.clear();
@@ -65,6 +66,8 @@ KineticModel::KineticModel ()
sulfateAttackTime_ = 1.0e10;
leachTime_ = 1.0e10;
+ pfk1_ = pfk2_ = pfk3_ = 1.0; // Default PK pozzolanic factors for no pozzolans
+
return;
}
@@ -139,6 +142,7 @@ KineticModel::KineticModel (ChemicalSystem *cs,
name_.clear();
isKinetic_.clear();
+ isPK_.clear();
isThermo_.clear();
isSoluble_.clear();
microPhaseId_.clear();
@@ -230,6 +234,8 @@ KineticModel::KineticModel (ChemicalSystem *cs,
}
chemSys_->setMicroInitVolume(totmicvol);
+ pfk1_ = pfk2_ = pfk3_ = 1.0; // Default PK factors for no pozzolanic material
+
return;
}
@@ -348,6 +354,7 @@ void KineticModel::parsePhase (xmlDocPtr doc,
string testname;
bool kineticfound = false;
bool iskin = false;
+ bool isPK = false;
bool istherm = false;
bool issol = false;
@@ -413,6 +420,12 @@ void KineticModel::parsePhase (xmlDocPtr doc,
}
isKinetic_.push_back(iskin);
+ string sfume("SFUME");
+ if (kineticData.name == sfume) {
+ isPK_.push_back(false);
+ } else {
+ isPK_.push_back(iskin);
+ }
isThermo_.push_back(istherm);
isSoluble_.push_back(issol);
initScaledMass_.push_back(0.0);
@@ -610,7 +623,7 @@ void KineticModel::setInitialPhaseVolumeFractions()
molarMass = chemSys_->getDCMolarMass(DCId); // g/mol
molarVolume = chemSys_->getDCMolarVolume(DCId); // m3/mol
density = 0.0;
- if (molarVolume > 1.0e-12) P
+ if (molarVolume > 1.0e-12) {
density = molarMass / molarVolume / 1.0e6; // g/cm3
}
microPhaseMass[microPhaseId] = volumeFraction * density;
@@ -738,7 +751,6 @@ void KineticModel::calculateKineticStep (const double timestep,
double DOH = 0.0;
double arrhenius = 1.0;
-
double ngrate = 1.0e-10; // Nucleation and growth rate
double hsrate = 1.0e-10; // Hydration shell rate
double diffrate = 1.0e-10; // Diffusion rate
@@ -821,7 +833,7 @@ void KineticModel::calculateKineticStep (const double timestep,
// Set the proper amount of water for the total solid mass
// and the water-solid ratio
- // Determine total intial solid mass
+ // Determine total initial solid mass
double solidMass = 0.0;
vector initSolidMasses = getInitScaledMass();
for (int i = 0; i < initSolidMasses.size(); i++) {
@@ -984,6 +996,39 @@ void KineticModel::calculateKineticStep (const double timestep,
p++;
}
+ // @todo BULLARD PLACEHOLDER
+ // While we are still using PK model for clinker phases, we
+ // must scan for and account for pozzolanic reactive components
+ // that will alter the hydration rate of clinker phases, presumably
+ // by making a denser hydration shell with slower transport
+
+ for (int i = 0; i < microPhaseId_.size(); i++) {
+
+ microPhaseId = microPhaseId_[i];
+
+ // @todo BULLARD PLACEHOLDER
+ // Currently silica fume is the only pozzolanically reactive
+ // component
+
+ if (chemSys_->getMicroPhaseName(microPhaseId) == "SFUME") {
+ double betarea = k1_[i];
+ double area = betarea * scaledMass_[i];
+ double loi = k2_[i];
+ double sio2 = critDOH_[i];
+
+ // Handle silica fume as a special case here:
+ // Only for silica fume, we let k1 = BET surface area in m2/g,
+ // k2 = LOI in percent by solid mass.
+ // critDOH = silica content in PERCENT BY MASS of silica fume
+
+ for (int ii = 0; ii < microPhaseId_.size(); ii++) {
+ if (isPK(ii)) {
+ pfk1_ = pfk2_ = pfk3_ = ((betarea/24.0) * (sio2/94.0) * (sio2/94.0));
+ }
+ }
+ }
+ }
+
} // End of special first-time tasks
if (hyd_time < leachTime_ && hyd_time < sulfateAttackTime_) {
@@ -1017,9 +1062,9 @@ void KineticModel::calculateKineticStep (const double timestep,
cout << "SS kinetic phase " << microPhaseId << endl;
cout << "SS name = " << chemSys_->getMicroPhaseName(microPhaseId)
<< endl;
- cout << "SS k1 = " << k1_[i] << endl;
- cout << "SS k2 = " << k2_[i] << endl;
- cout << "SS k3 = " << k3_[i] << endl;
+ cout << "SS k1 = " << k1_[i] * pfk1_ << endl;
+ cout << "SS k2 = " << k2_[i] * pfk2_ << endl;
+ cout << "SS k3 = " << k3_[i] * pfk3_ << endl;
cout << "SS n1 = " << n1_[i] << endl;
cout << "SS n3 = " << n3_[i] << endl;
cout << "SS Ea = " << activationEnergy_[i] << endl;
@@ -1028,6 +1073,9 @@ void KineticModel::calculateKineticStep (const double timestep,
if (initScaledMass_[i] > 0.0) {
DOH = (initScaledMass_[i] - scaledMass_[i]) /
(initScaledMass_[i]);
+ if (!isPK(i)) {
+ DOH = min(DOH,0.99); // prevents DOH from prematurely stopping PK calculations
+ }
} else {
throw FloatException("KineticModel","calculateKineticStep",
"initScaledMass_ = 0.0");
@@ -1037,70 +1085,154 @@ void KineticModel::calculateKineticStep (const double timestep,
(pow(((critDOH_[i] * wcRatio_) - DOH),4.0)));
arrhenius = exp((activationEnergy_[i]/GASCONSTANT)*((1.0/refT_) - (1.0/T)));
- // if (verbose_) {
- // cout << "QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ" << endl;
- // cout << "Calculating Arrhenius effect:" << endl;
- // cout << " GASCONSTANT = " << GASCONSTANT << endl;
- // cout << " activationEnergy_[" << i << "] = " << activationEnergy_[i] << endl;
- // cout << " refT_ = " << refT_ << endl;
- // cout << " T = " << T << endl;
- // cout << " arrhenius factor = " << arrhenius << endl;
- // cout << "QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ" << endl;
- // }
if (DOH < 1.0 && !doTweak) {
+
+ // @todo BULLARD PLACEHOLDER
+ // Handle silica fume as a special case here:
+ // Only for silica fume, we let k1 = BET surface area in m2/g,
+ // k2 = LOI in percent by solid mass.
+ // critDOH = silica content in PERCENT BY MASS of silica fume
+ //
+ // This is a TOTAL KLUGE right now from here...
+
+ if (chemSys_->getMicroPhaseName(microPhaseId) == "SFUME") {
+
+ // TBD, Dove and Crerar 1990 (Geochim et Cosmochim Acta)
+ // report k+ (aSiO2)(aH2O)^2 (1 - Q/K) for pure water, where
+ // k+ = 5.6e-9 mol m-2 s-1 in 0.05 m NaCl at 100 C and
+ // an activation enthalpy of about 71 kJ/mol
+
+ double baserateconst = 5.6e-3 * arrhenius; // mol m-2 s-1
+ //
+
+ double ca = chemSys_->getDCConcentration("Ca+2");
+ double kca = 4.0e-7; // mol m-2 s-1 ads. rate const for Ca (guess)
+ double Kca = 10.0; // adsorption equilibrium constant is a guess
+ double na = chemSys_->getDCConcentration("Na+");
+ double kna = 6.35e-7; // mol m-2 s-1 ads. rate const from Dove and Crerar
+ double Kna = 58.3; // adsorption equilibrium constant from Dove and Crerar
+ double k = chemSys_->getDCConcentration("K+");
+ double kk = 5.6e-7; // mol m-2 s-1 ads. rate const from Dove and Crerar
+ double Kk = 46.6; // adsorption equilibrium constant from Dove and Crerar
+
+ // Langmuir adsorption isotherms assumed to be additive
+
+ baserateconst += (kca * Kca * ca / (1.0 + (Kca * ca)));
+ baserateconst += (kna * Kna * na / (1.0 + (Kna * na)));
+ baserateconst += (kk * Kk * k / (1.0 + (Kk * k)));
+
+ double ohact = chemSys_->getDCActivity("OH-");
+ double betarea = k1_[i];
+ double area = betarea * scaledMass_[i];
+ double loi = k2_[i];
+ double sio2 = critDOH_[i];
+ arrhenius = exp((activationEnergy_[i]/GASCONSTANT)*((1.0/373.0) - (1.0/T)));
+ double molarmass = chemSys_->getDCMolarMass("Amor-Sl"); // g mol-1
+ vector GEMPhaseIndex;
+ GEMPhaseIndex.clear();
+ GEMPhaseIndex = chemSys_->getMicroPhaseToGEMPhase(microPhaseId);
+
+ // saturation index , but be sure that there is only one GEM Phase
+
+ double satindex = solut_->getSI(GEMPhaseIndex[0]);
+
+ // activity of water
+ double wateractivity = chemSys_->getDCActivity(chemSys_->getDCId("H2O@"));
+
+ // Make sure sio2 is given in percent by mass of silica fume
+ // This equation basically implements the Dove and Crerar rate equation
+ // for quartz. Needs to be calibrated for silica fume, but hopefully
+ // the BET area and LOI will help do that.
+
+ rate = baserateconst * ohact * area * pow(wateractivity,2.0)
+ * (1.0 - (loi/100.0)) * (sio2/100.0) * (1.0 - satindex);
+
+ massDissolved = (rate * timestep) * molarmass;
+
+ cout << "SFUME: baserateconst = " << baserateconst << " mol/m2/s" << endl;
+ cout << "SFUME: betarea = " << betarea << " m2/g" << endl;
+ cout << "SFUME: scaledMass = " << scaledMass_[i] << " g/(100 g)" << endl;
+ cout << "SFUME: area = " << area << " m2" << endl;
+ cout << "SFUME: loi = " << loi << " %" << endl;
+ cout << "SFUME: sio2 = " << sio2 << " %" << endl;
+ cout << "SFUME: arrhenius = " << arrhenius << endl;
+ cout << "SFUME: satindex = " << satindex << endl;
+ cout << "SFUME: wateractivity = " << wateractivity << endl;
+ cout << "SFUME: rate = " << rate << " mol/s" << endl;
+ cout << "SFUME: massDissolved = " << massDissolved << " g/(100 g)" << endl;
+ cout.flush();
+
+ scaledMass_[i] -= massDissolved;
+
+
+ // @todo BULLARD PLACEHOLDER ... to here
+
+ } else {
+
+ wcFactor = 1.0 + (3.333 * (pow(((critDOH_[i] * wcRatio_) - DOH),4.0)));
+
+ // Normal Parrott and Killoh implementation here
+
+ // @todo BULLARD PLACEHOLDER pfk1_, pfk2_, pfk3_ are
+ // pozzolanic modification factors for the Parrott and Killoh
+ // rate constants k1, k2, and k3.
- if (fabs(n1_[i]) > 0.0) {
- ngrate = (k1_[i]/n1_[i]) * (1.0 - DOH)
- * pow((-log(1.0 - DOH)),(1.0 - n1_[i]));
- ngrate *= (blaineFactor_); // only used for the N+G rate
+ if (fabs(n1_[i]) > 0.0) {
+ ngrate = (pfk1_ * k1_[i]/n1_[i]) * (1.0 - DOH)
+ * pow((-log(1.0 - DOH)),(1.0 - n1_[i]));
+ ngrate *= (blaineFactor_); // only used for the N+G rate
- if (ngrate < 1.0e-10) ngrate = 1.0e-10;
- } else {
- throw FloatException("KineticModel","calculateKineticStep",
- "n1_ = 0.0");
- }
+ if (ngrate < 1.0e-10) ngrate = 1.0e-10;
+ } else {
+ throw FloatException("KineticModel","calculateKineticStep",
+ "n1_ = 0.0");
+ }
- hsrate = k3_[i] * pow((1.0 - DOH),n3_[i]);
- if (hsrate < 1.0e-10) hsrate = 1.0e-10;
-
- if (DOH > 0.0) {
- diffrate = (k2_[i] * pow((1.0 - DOH),(2.0/3.0))) /
- (1.0 - pow((1.0 - DOH),(1.0/3.0)));
- if (diffrate < 1.0e-10) diffrate = 1.0e-10;
- } else {
- diffrate = 1.0e9;
- }
+ hsrate = pfk3_ * k3_[i] * pow((1.0 - DOH),n3_[i]);
+ if (hsrate < 1.0e-10) hsrate = 1.0e-10;
+
+ if (DOH > 0.0) {
+ diffrate = (pfk2_ * k2_[i] * pow((1.0 - DOH),(2.0/3.0))) /
+ (1.0 - pow((1.0 - DOH),(1.0/3.0)));
+ if (diffrate < 1.0e-10) diffrate = 1.0e-10;
+ } else {
+ diffrate = 1.0e9;
+ }
- rate = (ngrate < hsrate) ? ngrate : hsrate;
- if (diffrate < rate) rate = diffrate;
- rate *= (wcFactor * rhFactor * arrhenius);
- newDOH = DOH + (rate * timestep);
- if (verbose_) {
- cout << "PK model for " << getName(i)
- << ", ngrate = " << ngrate
- << ", hsrate = " << hsrate
- << ", diffrate = " << diffrate
- << ", rhFactor = " << rhFactor
- << ", wcFactor = " << wcFactor
- << ", timestep " << timestep
- << ", oldDOH = " << DOH << ", new DOH = "
- << newDOH << endl;
- cout.flush();
- }
+ rate = (ngrate < hsrate) ? ngrate : hsrate;
+ if (diffrate < rate) rate = diffrate;
+ rate *= (wcFactor * rhFactor * arrhenius);
+ newDOH = DOH + (rate * timestep);
+ if (verbose_) {
+ cout << "PK model for " << getName(i)
+ << ", ngrate = " << ngrate
+ << ", hsrate = " << hsrate
+ << ", diffrate = " << diffrate
+ << ", rhFactor = " << rhFactor
+ << ", wcFactor = " << wcFactor
+ << ", timestep " << timestep
+ << ", oldDOH = " << DOH << ", new DOH = "
+ << newDOH << endl;
+ cout.flush();
+ }
- /// @note This where we can figure out the volume dissolved
- /// and link it back to the current volume to see how many
- /// voxels need to dissolve
- ///
+ /// @note This where we can figure out the volume dissolved
+ /// and link it back to the current volume to see how many
+ /// voxels need to dissolve
+ ///
- /// @note This all depends on concept of degree of
- /// hydration as defined by the PK model
+ /// @note This all depends on concept of degree of
+ /// hydration as defined by the PK model
- /// @todo Make this independent of PK model
+ /// @todo Make this independent of PK model
- scaledMass_[i] = initScaledMass_[i] * (1.0 - newDOH);
- massDissolved = (newDOH - DOH) * initScaledMass_[i];
+ scaledMass_[i] = initScaledMass_[i] * (1.0 - newDOH);
+ massDissolved = (newDOH - DOH) * initScaledMass_[i];
+
+ // end it here?
+ }
+
chemSys_->setMicroPhaseMass(microPhaseId,scaledMass_[i]);
chemSys_->setMicroPhaseMassDissolved(microPhaseId,massDissolved);
if (verbose_) {
@@ -1114,10 +1246,10 @@ void KineticModel::calculateKineticStep (const double timestep,
/// @note impurityRelease index values are assumed to
/// be uniquely associated with particular chemical
/// elements
-
+
/// @todo Make this more general so that any indexing
/// can be used
-
+
/// @todo Allow any IC element to be an impurity, not
/// necessarily just the ones hard-coded here
@@ -1131,9 +1263,19 @@ void KineticModel::calculateKineticStep (const double timestep,
chemSys_->getSo3(microPhaseId));
for (int ii = 0; ii < ICMoles.size(); ii++) {
- ICMoles[ii] += ((massDissolved
- / chemSys_->getDCMolarMass(DCId_[i]))
- * chemSys_->getDCStoich(DCId_[i],ii));
+
+ /// @todo BULLARD PLACEHOLDER
+ /// Special case for SFUME here to account for SiO2 < 100%
+
+ double purity = 1.0;
+ if (chemSys_->getMicroPhaseName(microPhaseId) == "SFUME") {
+ purity = critDOH_[i] / 100.0;
+ }
+
+ ICMoles[ii] += ((purity * massDissolved
+ / chemSys_->getDCMolarMass(DCId_[i]))
+ * chemSys_->getDCStoich(DCId_[i],ii));
+
if (ICName[ii] == "O") {
// Dissolved K2O in this phase
if (chemSys_->isIC("K")) {
diff --git a/src/thameslib/KineticModel.h b/src/thameslib/KineticModel.h
index 9597973..a93b3e2 100644
--- a/src/thameslib/KineticModel.h
+++ b/src/thameslib/KineticModel.h
@@ -220,6 +220,7 @@ vector SO4Target_;
vector name_; /**< List of names of phases in the kinetic model */
vector isKinetic_; /**< List of ids of phases that are kinetically controlled */
+vector isPK_; /**< List of ids of phases that are PK-model controlled */
vector isThermo_; /**< List of ids of phases that are thermodynamically
controlled */
vector isSoluble_; /**< List of ids of phases that are instantly dissolved */
@@ -241,6 +242,8 @@ vector critDOH_; /**< List of critical degrees of hydration for w
effect in the Parrot and Killoh model */
vector degreeOfHydration_; /**< Degree of hydration of each kinetic (clinker) phase */
+double pfk1_,pfk2_,pfk3_; /**< Pozzolanic factors for the Parrott and Killoh parameters */
+
bool verbose_; /**< Flag for verbose output */
bool warning_; /**< Flag for warnining output */
bool debug_; /**< Flag for debugging output */
@@ -539,6 +542,44 @@ bool isKinetic (const string micname)
}
}
+/**
+@brief Get the Parrott and Killoh model status of a microstructure phase by its id
+
+@param idx is a microstructure id number
+@return true if it is kinetically controlled by PK model
+*/
+bool isPK (const unsigned int idx)
+{
+ try { return isPK_.at(idx); }
+ catch (out_of_range &oor) {
+ EOBException ex("KineticModel","isPK",
+ "isPK_",isPK_.size(),idx);
+ ex.printException();
+ exit(1);
+ }
+}
+
+/**
+@brief Get the Parrott and Killoh model status of a microstructure phase by its name
+
+@param micname is the name of a microstructure phase
+@return true if it is kinetically controlled by PK model
+*/
+bool isPK (const string micname)
+{
+ int idx;
+ try {
+ idx = chemSys_->getMicroPhaseId(micname);
+ return isPK_.at(idx);
+ }
+ catch (out_of_range &oor) {
+ EOBException ex("KineticModel","isPK",
+ "isPK_",isPK_.size(),idx);
+ ex.printException();
+ exit(1);
+ }
+}
+
/**
@brief Get the thermo status of a microstructure phase by its id
@@ -1182,6 +1223,66 @@ double getCritDOH (const unsigned int i) const
}
}
+/**
+@brief Set the "pozzolanic factor" for the Parrott and Killoh k1 parameter
+
+@param pk1 is the value to set for the k1 parameter's pozzolanic factor
+*/
+void setPfk1 (const double pk1)
+{
+ pfk1_ = pk1;
+}
+
+/**
+@brief Get the "pozzolanic factor" for the Parrott and Killoh k1 parameter
+
+@return the "pozzolanic factor" for the Parrott and Killoh k1 parameter
+*/
+double getPfk1 () const
+{
+ return pfk1_;
+}
+
+/**
+@brief Set the "pozzolanic factor" for the Parrott and Killoh k2 parameter
+
+@param pk2 is the value to set for the k2 parameter's pozzolanic factor
+*/
+void setPfk2 (const double pk2)
+{
+ pfk2_ = pk2;
+}
+
+/**
+@brief Get the "pozzolanic factor" for the Parrott and Killoh k2 parameter
+
+@return the "pozzolanic factor" for the Parrott and Killoh k2 parameter
+*/
+double getPfk2 () const
+{
+ return pfk2_;
+}
+
+/**
+@brief Set the "pozzolanic factor" for the Parrott and Killoh k3 parameter
+
+@param pk3 is the value to set for the k3 parameter's pozzolanic factor
+*/
+void setPfk3 (const double pk3)
+{
+ pfk3_ = pk3;
+}
+
+/**
+@brief Get the "pozzolanic factor" for the Parrott and Killoh k3 parameter
+
+@return the "pozzolanic factor" for the Parrott and Killoh k3 parameter
+*/
+double getPfk3 () const
+{
+ return pfk3_;
+}
+
/**
@brief Get the list of degrees of hydration for clinker phases in the kinetic model.