Full ranged ZrCoHx PCT modelling capabilities

doc/content/source/interfacekernels/ADMatInterfaceReactionZrCoHxPCT.md

@@ -4,52 +4,45 @@
 
 ## Description
 
-The composition of metal hydrides like zirconium-cobalt hydride (ZrCoHx) is described by pressure-composition-temperature (PCT) data.
+The composition of metal hydrides like zirconium-cobalt hydride (ZrCoH$_{x}$) is described by pressure-composition-temperature (PCT) data.
 In TMAP8, the PCT curve can be imposed by an interface kernel that then dictates the material composition.
-Hence, the `ADMatInterfaceReactionZrCoHxPCT` interface kernel imposes the surface concentration of H in ZrCoHx based on the input pressure (Pa) and temperature (K). `ADMatInterfaceReactionZrCoHxPCT` is related to [ADMatInterfaceReaction.md], but includes the ZrCoHx PCT curves.
+Hence, the `ADMatInterfaceReactionZrCoHxPCT` interface kernel imposes the surface concentration of H in ZrCoH$_{x}$ based on the input pressure (Pa) and temperature (K). `ADMatInterfaceReactionZrCoHxPCT` is related to [ADMatInterfaceReaction.md], but includes the ZrCoH$_{x}$ PCT curves.
  At the interface between a solid (main) and a gas (neighbor), it imposes:
 \begin{equation} \label{eq:test_interfacereactionZr2FeHx}
 \frac{d C_s}{dt} = 0 = K_b f_{at}(T,P) \rho -  K_f C_s ,
 \end{equation}
 where $C_s$ is the surface H concentration in mol/m$^3$,
 $K_b$ and $K_f$ are the backward and forward surface reaction rate in 1/s, respectively,
-$f_{at}(T,P)$ is the composition in atomic fraction of H in ZrCoHx given a gas temperature $T$ and gas pressure $P$,
-and $\rho$ is the zirconium-cobalt atomic density in mol/m$^3$.
+$f_{at}(T,P)$ is the composition in atomic fraction of H in ZrCoH$_{x}$ given a gas temperature $T$ and gas pressure $P$,
+and $\rho$ is the zirconium-cobalt molar density in mol/m$^3$.
 
-[ZrCoHx_PCT_Data] shows the data used in this interface kernel. The experimental data was selected from several authors [!cite](jat2013hydrogen) for (433.15 K) and (604.15 K) and [!cite](nagasaki1986zirconium) for (573.15 K) since a clean PCT curve is required for accurate modelling.
+[ZrCoHx_PCT_combined] shows the experimental data, fitted equation, and TMAP8 simulation results. The experimental data was selected from several authors for 433.15 K [!cite](nagasaki1986zirconium) and [!cite](jat2013hydrogen) for 524.15 K, 544.15 K, 584.15 K, 604.15 K, and 624.15 K.
 
-!media comparison_ZrCoHx_PCT.py
-       image_name=ZrCoHx_PCT_Data.png
-       style=width:80%;margin-bottom:2%;margin-left:auto;margin-right:auto
-       id=ZrCoHx_PCT_Data
-       caption=PCT data for ZrCoHx from [!cite](jat2013hydrogen) and [!cite](nagasaki1986zirconium).
 
 To include this PCT data in TMAP8 modelling capabilities the high and low pressure regions were extracted and regressed for the resulting equations.
 
 The low pressure is captured as:
-\begin{equation}\label{eq:low_pressure}
-f_{at}(T,P) = 0.5-\left[1.00 \times 10^{-3}+\exp(-4.29+ 0.020 T + (-1.07+ 5.69 \times 10^{-4} T) (\log\left(P_{lim}(T) - P\right)))\right]^{-1},
+\begin{equation}\label{eq:atomic_fraction_LP}
+f_{at,LP}(T,P) = 0.7-1\left[5 \times 10^{-3}+\exp(-4.37 + 1.34 \times 10^{-2} T + ( -8.22 \times 10^{-2} - 3.97 \times 10^{-4} T) (\log\left(P_{lim}(T) - P \right)))\right]^{-1}.
 \end{equation}
 
 The high pressure is captured as:
-\begin{equation}\label{eq:high_pressure}
-f_{at}(T,P) = 2.5-3.42\left[1.4+\exp(7.97 - 0.012 T + (-0.17 + 1.19 \times 10^{-3} T) (\log\left(P - P_{lim}(T)\right)))\right]^{-1},
+\begin{equation}\label{eq:atomic_fraction_HP}
+f_{at,HP}(T,P) = 2.7-1.45\left[1.00+\exp(6.57 - 2.21 \times 10^{-2} T + ( 6.52 \times 10^{-1} - 1.17 \times 10^{-5} T) (\log\left(P - P_{lim}(T)\right)))\right]^{-1}.
 \end{equation}
 
-[ZrCoHx_PCT_fit_2D] shows the fitting and TMAP8 results used in this interface kernel.
-
 !media comparison_ZrCoHx_PCT.py
-       image_name=ZrCoHx_PCT_fit_2D.png
+       image_name=ZrCoHx_PCT_combined.png
        style=width:80%;margin-bottom:2%;margin-left:auto;margin-right:auto
-       id=ZrCoHx_PCT_fit_2D
-       caption=PCT data for ZrCoHx from [!cite](jat2013hydrogen) and [!cite](nagasaki1986zirconium) implemented in TMAP8, and test results.
+       id=ZrCoHx_PCT_combined
+       caption=PCT data for $ZrCoH_{x}$ from [!cite](jat2013hydrogen) and [!cite](nagasaki1986zirconium) with fits implemented in TMAP8 and test results.
 
 
 The plateau representing phase transition is captured as:
 \begin{equation} \label{eq:pressure_plateau}
-P_{lim} = \exp\left(12.43-4.84 \times 10^{-2} T +7.15 \times 10^{-5} T^2 \right),
+P_{lim} = \exp\left(-9.41 + 3.32 \times 10^{-2} T - 3.30 \times 10^{-6} T^2 \right).
 \end{equation}
-with $P_{lim}$ being the hydrogen partial pressure limit delineating the plateau in Pa and $T$ being the temperature in K.
+
 This fit is shown in [ZrCoHx_PCT_plateau_pressure_fit].
 
 !media comparison_ZrCoHx_PCT.py
@@ -58,7 +51,77 @@ This fit is shown in [ZrCoHx_PCT_plateau_pressure_fit].
        id=ZrCoHx_PCT_plateau_pressure_fit
        caption=Fit phase transition region as pressure as a function of temperature.
 
-The jump between the low low high pressure occurs if the atomic fraction is equal to 0.50, while from high to low pressure occurs if the atomic fraction is equal to 1.4. These end points follow the [eq:low_pressure] and [eq:high_pressure] maximum and minimum,respectively.
+The subsequent set of equations are bridged using a simoid blending methodology that is present below.
+
+First we express a linear logarithmic interpolation:
+\begin{equation}\label{eq:mid_interp}
+f_{mid}(P)= m{0} + m_{1}\ln(P),
+\end{equation}
+
+where the slope of [eq:mid_interp] are given as
+
+\begin{equation}\label{eq:m0_definition}
+m_{0} = f_{LP}(T,\alpha P_{lim}) - \frac{f_{HP}(T,\beta P{lim}) - f_{LP}(T,\alpha P{lim})}{L_{b} - L_{a}} L_{a},
+\end{equation}
+and
+\begin{equation}\label{eq:m1_definition}
+m_{1} = \frac{f_{HP}(T,\beta P_{lim}) - f_{LP}(T,\alpha P{lim})}{L_{b} - L_{a}},
+\end{equation}
+with
+\begin{equation}
+\label{eq:log_bounds}
+L_{a} = \ln(\alpha P_{lim}), \qquad
+L_{b} = \ln(\beta P_{lim}),
+\end{equation}
+
+where $\beta$ and $\alpha$ represent the high- and low-pressure transition as a ratio of high-to-plateau pressure and plateau-to-low pressure regions. For the ratio of plateau-to-low pressure regions a constant value of 1.008 was found to be suitable, while for the high-to-plateau pressure a fitted function, $\beta$, was required.
+
+$\beta$ and $\alpha$ are defined as
+\begin{equation}
+\label{eq:alpha_beta_defs}
+\beta = 2.39 - 5.10 \times 10^{-3} T + 5.42 \times 10^{-6} T^{2}, \qquad
+\alpha = 1.008,
+\end{equation}
+
+with the $f_{mid}(P)$ fully described a weight sigmoid blending function can be applied with the following equations:
+
+\begin{equation}
+\label{eq:scaled_variables}
+r = \frac{P}{P_{lim}}, \qquad
+x = \ln(r), \qquad
+x_{\alpha} = \ln(\alpha), \qquad
+x_{\beta} = \ln(\beta)
+\end{equation}
+and
+\begin{equation}
+\label{eq:switching_functions}
+s_{(LP\to mid)} = \frac{1}{1 + e^{-(x - x_{\alpha})/\Delta_{\alpha}}}, \qquad
+s_{(mid\to HP)} = \frac{1}{1 + e^{-(x - x_{\beta})/\Delta_{\beta}}},
+\end{equation}
+
+where $\Delta_{\beta}$ is the tunable base widths of the smooth blending function set to a default 0.08. Subsequent of the sigmoid blending function, a normalization of their weights to induce a smooth transition is written as such:
+
+\begin{equation}
+\label{eq:weight_normalization}
+W = (1 - s_{(LP\to mid)}) + s_{(LP\to mid)}(1 - s_{(mid\to HP)}) + s_{(mid\to HP)},
+\end{equation}
+
+\begin{equation}
+\label{eq:weight_definitions}
+w_{LP} = \frac{1 - s_{(LP\to mid)}}{W}, \qquad
+w_{mid} = \frac{s_{(LP\to mid)}(1 - s_{(mid\to HP)})}{W}, \qquad
+w_{HP} = \frac{s_{(mid\to HP)}}{W}
+\end{equation}
+
+\begin{equation}
+\label{eq:overall_interpolation}
+f_{overall}(T,P) = w_{LP}f_{LP}(T,P) + w_{mid}f_{mid}(T,P) + w_{HP}f_{HP}(T,P).
+\end{equation}
+
+Lastly, the validity of this present fit is between:
+\begin{equation} \label{eq:bounds}
+20 < P\,\text{[Pa]} < 2 \times 10^{5}.
+\end{equation}
 
 ## Test
 
@@ -69,9 +132,30 @@ To model the interface, the input file employs the [InterfaceDiffusion.md] objec
 \begin{equation} \label{eq:test_interfacereaction}
 \frac{d C_s}{dt} = 0 = K (f_{at}(T,P) \rho - C_s),
 \end{equation}
-where $\rho$ is the zirconium-cobalt atomic density.
 
-The results of the high pressure test for ($T$, $P$) = (433.15 K, $3 \times 10^{4}$ Pa), (433.15 K, $1 \times 10^{4}$ Pa), (573.15 K, $1 \times 10^{4}$ Pa), and (604.15 K, $5 \times 10^{4}$ Pa). The results of the low pressure test for ($T$, $P$) = (433.15 K, $1 \times 10^{2}$ Pa), (573.15 K, $1 \times 10^{3}$ Pa), (604.15 K, $1 \times 10^{4}$ Pa), and (604.15 K, $3 \times 10^{3}$ Pa) are shown in [ZrCoHx_PCT_fit_2D]. Both show good agreement.
+
+The results of the high pressure test for ($T$, $P$) = (433.15 K, $3 \times 10^{4}$ Pa), (433.15 K, $1 \times 10^{4}$ Pa), (573.15 K, $1 \times 10^{4}$ Pa), and (604.15 K, $5 \times 10^{4}$ Pa). The results of the low pressure test for ($T$, $P$) = (433.15 K, $1 \times 10^{2}$ Pa), (573.15 K, $1 \times 10^{3}$ Pa), (604.15 K, $1 \times 10^{4}$ Pa), and (604.15 K, $3 \times 10^{3}$ Pa) are shown in [ZrCoHx_PCT_combined]. Overall, the analytical fits show reasonable agreement with the measurements, yielding average root‑mean‑square errors of 0.056 for the low‑pressure region and 0.089 for the high‑pressure region.
+
+
+### Overall PCT testing
+
+[!file](/ZrCoHx_PCT.i) also tests the entire PCT modelling curves in TMAP8. The simulation follows the same structure as before, but has an arbitrarily high diffusion value to acheive steady-state quickly. The model also includes a linear pressure increase to cover the entire PCT curve. The [ZrCoHx_PCT_combined] shows the PCT fit against the experimental data. Evidently, the fit exhibits moderate average error of 24% , but the deviations remain within an acceptable range for modeling purposes
+
+
+\begin{equation}
+P = P_{\text{initial}} + t \frac{P_{\text{max}} - P_{\text{initial}}}{t_{\text{end}}}.
+\end{equation}
+
+The testing conditions include
+($T$, $P_{initial}$) = (423.15 K, $2.5 \times 10^{1}$ Pa),
+(524.15 K, $2.5 \times 10^{1}$ Pa),
+(544.15 K, $2.5 \times 10^{1}$ Pa),
+(564.15 K, $2.5 \times 10^{1}$ Pa),
+(584.15 K, $2.5 \times 10^{1}$ Pa),
+(604.15 K, $2.5 \times 10^{1}$ Pa),
+and (624.15 K, $2.5 \times 10^{1}$ Pa).
+
+
 
 ## Example Input File Syntax
 

src/interfacekernels/ADMatInterfaceReactionZrCoHxPCT.C

@@ -45,61 +45,149 @@ ADMatInterfaceReactionZrCoHxPCT::ADMatInterfaceReactionZrCoHxPCT(const InputPara
 ADReal
 ADMatInterfaceReactionZrCoHxPCT::computeQpResidual(Moose::DGResidualType type)
 {
+
+  // Variables
+  ADReal r = 0.;
+  ADReal m1 = 0.0;
+  ADReal m0 = 0.0;
+
+  // Local conveniences
   using std::exp;
   using std::log;
   using std::max;
-  using std::pow;
-  ADReal r = 0;
-  // tolerance for the pressure being closed to the plateau region
-  const Real tolerance = 10; // Pa
-  // Calculate the equilibrium concentration value based on PCT curve
-  // (/2 because two atoms for a molecule) (pressure in Pa)V3_projects/TMAP8/src/bcs
-  auto neighbor_pressure =
-      PhysicalConstants::ideal_gas_constant * _neighbor_temperature[_qp] * _neighbor_value[_qp] / 2;
-
-  // Calculate the value of the pressures for the phase transition plateau (pressure in Pa)
-  auto limit_pressure = exp(12.43 - 4.84e-2 * _neighbor_temperature[_qp] +
-                            7.15e-5 * Utility::pow<2>(_neighbor_temperature[_qp]));
-
-  // define atomic fraction variable
-  ADReal atomic_fraction = 0.0;
+
+  // Fitting constants of High-Pressure (HP)
+  const Real HP_A = 2.7;
+  const Real HP_B = 1.45;
+  const Real HP_C = 1.00;
+  const Real HP_D = 6.57;
+  const Real HP_E = 2.21e-02;
+  const Real HP_F = 6.52e-01;
+  const Real HP_G = 1.17e-05;
+  // Fitting constants of  Low-Pressure (LP)
+  const Real LP_A = 0.7;
+  const Real LP_B = 1.0;
+  const Real LP_C = 5.e-03;
+  const Real LP_D = -4.37;
+  const Real LP_E = -1.34e-02;
+  const Real LP_F = -8.22e-02;
+  const Real LP_G = 3.97e-04;
+
+  // -------------------------------
+  // Shared PCT correlation shape
+  // -------------------------------
+  // Both the LP and HP branches (and their values at the blending boundaries
+  // P_a and P_b) evaluate the exact same functional form:
+  //
+  //   af(P) = A - B / (C + exp(D - E*T + (F - G*T) * log(max(pressure_diff, 1e-10))))
+  //
+  // Only the fitting-constant set (A..G) and the "distance from the plateau
+  // pressure" term (pressure_diff) differ between calls
+
+  auto pctCorrelation = [&](const auto & A,
+                            const auto & B,
+                            const auto & C,
+                            const auto & D,
+                            const auto & E,
+                            const auto & F,
+                            const auto & G,
+                            const ADReal & pressure_diff) -> ADReal
+  {
+    return A - B / (C + exp(D - E * _neighbor_temperature[_qp] +
+                            (F - G * _neighbor_temperature[_qp]) * log(max(pressure_diff, 1e-10))));
+  };
+
+  // Gas pressure (Pa): R * T * c / 2 (two atoms per molecule)
+  const ADReal neighbor_pressure = PhysicalConstants::ideal_gas_constant *
+                                   _neighbor_temperature[_qp] * _neighbor_value[_qp] / 2.0;
 
   // Give a warning if the initial or computed neighbor pressure is out of the analytical model
-  if (((neighbor_pressure > 9.e06) || (neighbor_pressure < 0.011)))
+  if (((neighbor_pressure < 20) || (neighbor_pressure > 2.e5)))
     mooseDoOnce(mooseWarning("In ZrCoHxPCT: pressure ",
                              neighbor_pressure,
                              "Pa and temperature ",
                              _neighbor_temperature[_qp],
                              "K are outside the bounds of the atomic fraction correlation. See "
                              "documentation for ZrCoHxPCT material."));
 
-  if (neighbor_pressure > limit_pressure && abs(neighbor_pressure - limit_pressure) < tolerance)
-  {
-    // High pressure region, near limit.
-    atomic_fraction = 0.50;
-  }
-  else if (neighbor_pressure > limit_pressure)
-  {
-    // High pressure region
-    atomic_fraction =
-        2.5 - 3.42 / (1.40 + exp(7.97 - 1.99e-02 * _neighbor_temperature[_qp] +
-                                 (-1.69e-01 + 1.19e-03 * _neighbor_temperature[_qp]) *
-                                     log(max(neighbor_pressure - limit_pressure, 1.e-10))));
-  }
-  else if (neighbor_pressure < limit_pressure &&
-           abs(neighbor_pressure - limit_pressure) < tolerance)
-  {
-    // Low pressure region, near limit
-    atomic_fraction = 1.4;
-  }
-  else if (neighbor_pressure < limit_pressure)
-  {
-    // Low pressure region
-    atomic_fraction =
-        0.5 - 1. / (0.001 + exp(-4.29 + 1.98e-02 * _neighbor_temperature[_qp] +
-                                (-1.07 + 5.69e-04 * _neighbor_temperature[_qp]) *
-                                    log(max(limit_pressure - neighbor_pressure, 1.e-10))));
-  }
+  // Plateau / limit pressure (Pa)
+  const ADReal PLim = exp(-9.41 + 3.32e-02 * _neighbor_temperature[_qp] -
+                          3.30e-06 * Utility::pow<2>(_neighbor_temperature[_qp]));
+
+  // Transition fitted parameters, beta -> high pressure; alpha -> low pressure
+  const ADReal beta_corr = 2.39 - 5.1e-03 * _neighbor_temperature[_qp] +
+                           5.42e-06 * Utility::pow<2>(_neighbor_temperature[_qp]);
+  const Real alpha = 1.008;
+
+  // -------------------------------
+  // Ratio r = P / P_limit   (used for blending thresholds in log-space)
+  // -------------------------------
+  const ADReal ratio_r = neighbor_pressure / PLim;
+
+  // Low pressure (LP) branch: distance below the plateau pressure
+  const ADReal f_LP =
+      pctCorrelation(LP_A, LP_B, LP_C, LP_D, LP_E, LP_F, LP_G, PLim - neighbor_pressure);
+
+  // High pressure (HP) branch: distance above the plateau pressure
+  const ADReal f_HP =
+      pctCorrelation(HP_A, HP_B, HP_C, HP_D, HP_E, HP_F, HP_G, neighbor_pressure - PLim);
+
+  // -------------------------------
+  // Mid branch (exact continuity at alpha*Plim and beta*Plim)
+  // -------------------------------
+  // Boundaries in absolute pressure
+  const ADReal alpha_Plim = alpha * PLim;    // LP
+  const ADReal beta_Plim = beta_corr * PLim; // HP
+
+  // LP value at alpha_Plim
+  const ADReal f_LP_alpha_Plim =
+      pctCorrelation(LP_A, LP_B, LP_C, LP_D, LP_E, LP_F, LP_G, PLim - alpha_Plim);
+
+  // HP value at beta_Plim
+  const ADReal f_HP_beta_Plim =
+      pctCorrelation(HP_A, HP_B, HP_C, HP_D, HP_E, HP_F, HP_G, beta_Plim - PLim);
+
+  // Solve for af_mid(P) = m0 + m1 * log(P)
+  const ADReal L_a = log(max(alpha_Plim, 1e-10));
+  const ADReal L_b = log(max(beta_Plim, 1e-10));
+
+  // Determine af_mid(P) slopes
+  m1 = (f_HP_beta_Plim - f_LP_alpha_Plim) / (L_b - L_a);
+  m0 = f_LP_alpha_Plim - (f_HP_beta_Plim - f_LP_alpha_Plim) / (L_b - L_a) * L_a;
+
+  // Mid branch at current pressure using limited slope
+  const ADReal f_mid = m0 + m1 * log(max(neighbor_pressure, 1e-10));
+
+  // -------------------------------
+  // Smooth blending in log-space (LP ↔ mid ↔ HP)
+  // -------------------------------
+  // Base widths (tunable)
+  const Real base_delta_alpha_log = 0.08; // typical: 0.05–0.12
+  const Real base_delta_beta_log = 0.08;  // typical: 0.05–0.12
+
+  const ADReal x = log(max(ratio_r, 1e-10));      // log(ratio)
+  const Real x_alpha = log(alpha);                // constant threshold (Real)
+  const ADReal x_beta = log(max(beta_corr, 1.0)); // log(beta(T))
+
+  // Sigmoid steps (AD-safe)
+  const ADReal s_LP_to_mid = 1.0 / (1.0 + exp(-(x - x_alpha) / base_delta_alpha_log)); // LP→mid
+  const ADReal s_mid_to_HP = 1.0 / (1.0 + exp(-(x - x_beta) / base_delta_beta_log));   // mid→HP
+
+  // Weights
+  ADReal w_LP = 1.0 - s_LP_to_mid;                  // Low pressure weights
+  ADReal w_mid = s_LP_to_mid * (1.0 - s_mid_to_HP); // Mid pressure weights
+  ADReal w_HP = s_mid_to_HP;                        // High pressure weights
+
+  // Normalization
+  const ADReal w_sum = w_LP + w_mid + w_HP; // Sum of weights
+  w_LP /= w_sum;                            // Ratio of weight for low pressure
+  w_mid /= w_sum;                           // Ratio of weight for mid pressure
+  w_HP /= w_sum;                            // Ratio of weight for high pressure
+
+  // -------------------------------
+  // Final atomic fraction
+  // -------------------------------
+  ADReal atomic_fraction = w_LP * f_LP + w_mid * f_mid + w_HP * f_HP;
 
   // Convert to concentration
   auto _surface_equilibrium_concentration = atomic_fraction * _density[_qp];

test/tests/ZrCo_hydrogen_system/PCT_data/423.csv

@@ -0,0 +1,14 @@
+Atom Ratio,Partial Pressure
+0.08738286,31.22088945
+0.182491238,50.86914399
+0.30880347,54.78024527
+0.471377355,56.41037766
+0.736350578,58.31820093
+1.061739771,59.38947601
+1.458867648,67.29307926
+1.926949587,133.7730298
+2.401767212,2026.015385
+2.648476825,22534.45061
+2.648476825,22534.45061
+2.670760119,45063.02163
+2.6923795,116145.6631