Updated the shear modulus to depend on T in the EOS (#198).

ORNL-Fusion · Oct 18, 2024 · b2c1377 · b2c1377
1 parent 8c6267d
commit b2c1377
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 4 deletions.
diff --git a/xolotl/core/include/xolotl/core/Constants.h b/xolotl/core/include/xolotl/core/Constants.h
@@ -112,7 +112,9 @@ constexpr IdType psiVThreshold = 150;
 
 // Parameters for W Benedict EOS
 constexpr double gammaEOS = 2.67;
-constexpr double gEOS = 156.1e9;
+constexpr double gEOS1 = 1602.8;
+constexpr double gEOS2 = 0.1456;
+constexpr double gEOS3 = 3.28e-5;
 constexpr double kEOS = 1.380649e-23;
 constexpr double A_11 = -5.5991e-27;
 constexpr double A01 = 1.74e-26;

diff --git a/xolotl/util/include/xolotl/util/MathUtils.h b/xolotl/util/include/xolotl/util/MathUtils.h
@@ -121,18 +121,22 @@ getMaxHePerVLoop(double amtV, double latticeParameter, double temp) noexcept
 	double omega =
 		0.5 * latticeParameter * latticeParameter * latticeParameter * 1.0e-27;
 	double bEOS = latticeParameter * 0.5 * sqrt(3.0) * 1.0e-9;
+	double tempCelsius = temp - 273.0;
+	double gEOS = (xolotl::core::gEOS1 - xolotl::core::gEOS2 * tempCelsius -
+					  xolotl::core::gEOS3 * tempCelsius * tempCelsius) *
+		1.0e6;
 	double nume = amtV * omega;
 	double Fs = pow(3.0 * omega / (4.0 * xolotl::core::pi), 1.0 / 3.0);
 	double term1 = computeBenedictF1(temp) * pow(Fs, 1.0 / 3.0) *
 		pow((double)amtV, 1.0 / 9.0) *
-		pow(2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS,
+		pow(2.0 * xolotl::core::gammaEOS + gEOS * bEOS,
 			-1.0 / 3.0);
 	double term2 = computeBenedictF2(temp) * pow(Fs, 2.0 / 3.0) *
 		pow((double)amtV, 2.0 / 9.0) *
-		pow(2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS,
+		pow(2.0 * xolotl::core::gammaEOS + gEOS * bEOS,
 			-2.0 / 3.0);
 	double term3 = computeBenedictF3(temp) * Fs * pow((double)amtV, 1.0 / 3.0) /
-		(2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS);
+		(2.0 * xolotl::core::gammaEOS + gEOS * bEOS);
 
 	return util::max((nume / (term1 + term2 + term3)),
 		maxHePerV[maxHePerV.size() - 1] + amtV - maxHePerV.size() + 1.0);