Merge pull request idaholab#21165 from dschwen/slkks_5862

Implement Sublattice KKS phase field model

modules/phase_field/doc/content/bib/phase_field.bib

@@ -157,3 +157,27 @@ @article{Kobayashi1993
 author = "Ryo Kobayashi",
 abstract = "A simple phase field model for one component melt growth is presented, which includes anisotropy in a certain form. The formation of various dendritic patterns can be shown by a series of numerical simulations of this model. Qualitative relations between the shapes of crystals and some physical parameters are discussed. Also it is shown that noises give a crucial influence on the side branch structure of dendrites in some situations."
 }
+
+@article{Schwen2021,
+title = {A sublattice phase-field model for direct CALPHAD database coupling},
+journal = {Computational Materials Science},
+volume = {195},
+pages = {110466},
+year = {2021},
+issn = {0927-0256},
+doi = {https://doi.org/10.1016/j.commatsci.2021.110466},
+url = {https://www.sciencedirect.com/science/article/pii/S0927025621001919},
+author = {D. Schwen and C. Jiang and L.K. Aagesen},
+keywords = {Phase-field, CALPHAD, Automatic differentiation},
+abstract = {The phase-field method has been established as a de facto standard for simulating the microstructural evolution of materials. In quantitative modeling the assessment and compilation of thermodynamic/kinetic data is largely dominated by the CALPHAD approach, which has produced a large set of experimentally and computationally generated Gibbs free energy and atomic mobility data in a standardized format: the thermodynamic database (TDB) file format. Harnessing this data for the purpose of phase-field modeling is an ongoing effort encompassing a wide variety of approaches. In this paper, we aim to directly link CALPHAD data to the phase-field method, without intermediate fitting or interpolation steps. We introduce a model based on the Kim-Kim-Suzuki (KKS) approach. This model includes sublattice site fractions and can directly utilize data from TDB files. Using this approach, we demonstrate the model on the U-Zr and Mo-Ni-Re systems.}
+}
+
+@article{jacob2018revised,
+  title={Revised thermodynamic description of the Fe-Cr system based on an improved sublattice model of the $\sigma$ phase},
+  author={Jacob, Aur{\'e}lie and Povoden-Karadeniz, Erwin and Kozeschnik, Ernst},
+  journal={Calphad},
+  volume={60},
+  pages={16--28},
+  year={2018},
+  publisher={Elsevier}
+}

modules/phase_field/doc/content/modules/phase_field/CALPHAD.md

@@ -2,27 +2,45 @@
 
 This is a work in progress.
 
-A script to export MOOSE readable free energy expressions from `*.tdb` thermodynamic database files can be found at
+A script to export MOOSE readable free energy expressions from `*.tdb`
+thermodynamic database files can be found at
 
 ```text
 moose/python/calphad/free_energy.py
 ```
 
-The `free_energy.py` tool allows users to extract free energy expressions from `*.tdb`
-thermodynamic database files (ThermoCalc format). The tool exports a list of MOOSE Material blocks
-for each phase (or a user specified subset of phases). The CALPHAD functional expressions are
-implemented using the [`DerivativeParsedMaterial`](/DerivativeParsedMaterial.md) class.
+The `free_energy.py` tool allows users to extract free energy expressions from
+`*.tdb` thermodynamic database files (ThermoCalc format). The tool exports a
+list of MOOSE Material blocks for each phase (or a user specified subset of
+phases). The CALPHAD functional expressions are implemented using the
+[`DerivativeParsedMaterial`](/DerivativeParsedMaterial.md) class.
 
-It is up to the user to construct a full MOOSE input file around these material blocks and to rename
-the variables used in the exported form to more suitable names.
+It is up to the user to construct a full MOOSE input file around these material
+blocks and to rename the variables used in the exported form to more suitable
+names.
 
 ## Installation
 
-The [pycalphad](https://github.com/richardotis/pycalphad) Python module is required to perform the
-parsing of the `*.tdb` files. [SymPy](https://github.com/sympy/sympy) is required to build the
-functional forms of the calphad expressions.
+The [pycalphad](https://github.com/richardotis/pycalphad) Python module is
+required to perform the parsing of the `*.tdb` files.
+[SymPy](https://github.com/sympy/sympy) is required to build the functional
+forms of the calphad expressions.
 
-To install and/or upgrade these prerequisites use pip:
+### conda installation
+
+If you are using the MOOSE conda/mamba environment it is best to create a
+dedicated pycalphad environment:
+
+```
+mamba create -n pycalphad pycalphad sympy jupyterlab
+mamba activate pycalphad
+```
+
+`jupyterlab` is optional for viewing the example notebook under `moose/modules/phase_field/examples/slkks`.
+
+### pip Installation
+
+To install and/or upgrade these prerequisites using pip:
 
 ```text
 pip install --upgrade sympy
@@ -34,3 +52,4 @@ pip install --upgrade pycalphad
 Database files can be obtained online at
 
 * [Computational Phase Diagram Database](http://cpddb.nims.go.jp/index_en.html) (CPDDB)
+* [Thermodynamic Database Database](https://avdwgroup.engin.brown.edu/) (TDBDB)

modules/phase_field/doc/content/modules/phase_field/MultiPhase/SLKKS.md

@@ -0,0 +1,134 @@
+# Sublattice KKS model
+
+The sublattice Kim-Kim-Suzuki (SLKKS) [!cite](Schwen2021) model is an extension
+of the original KKS model incorporating an additional equal chemical potential
+constraint among the sublattice concentration in each phase.
+
+As such each component in the system will have corresponding phase
+concentrations, which are split up into per-sublattice concentrations.
+
+\begin{equation}
+c_i = \sum_j h(\eta_j)\sum_k a_{jk} c_{ijk},
+\end{equation}
+
+## Nomenclature
+
+$i$ indexes a component, $j$ a phase, and $k$ a sublattice in the given
+phase. $a_{jk}$ is the fraction of $k$ sublattice sites in phase $j$.
+
+## Sublattice equilibrium
+
+All sublattice pairs $k$ and $k'$ in a given phase are assumed to be always in
+equilibrium, thus
+
+\begin{equation}
+\frac 1{a_{jk}} \frac{\partial F_j}{\partial c_{ijk}} = \frac 1{a_{jk'}} \frac{\partial F_j}{\partial c_{ijk'}}
+\end{equation}
+
+This condition is enforced by the [`SLKKSChemicalPotential`](SLKKSChemicalPotential.md) kernel.
+
+```style=background-color:gray
+[chempot1a1b]
+  type = SLKKSChemicalPotential
+  variable = Cijk
+  k        = ajk
+  cs       = Cijk'
+  ks       = ajk'
+  F  = Fj
+[]
+```
+
+With $N$ sublattices in a phase $N-1$ such kernels are required. That leaves one
+sublattice concentration variable in a given phase to be covered by a kernel.
+
+# Phase equilibrium
+
+At the same time the original KKS model chemical potential equilibria still hold
+
+\begin{equation}
+\frac{\partial F}{\partial c_i} = \frac{\partial F_j}{\partial c_{ijk}} \quad ,\quad \forall \{j,k\}.
+\end{equation}
+
+meaning that sublattice chemical potentials (corrected for the sublattice site
+fractions) must be in equilibrium across phases.
+
+\begin{equation}
+\frac 1{a_{jk}} \frac{\partial F_j}{\partial c_{ijk}} = \frac 1{a_{j'k}} \frac{\partial F_{j'}}{\partial c_{ij'k}}
+\end{equation}
+
+That remaining sublattice concentration will couple to the next phase using a
+[`KKSPhaseChemicalPotential`](/KKSPhaseChemicalPotential.md) kernel.
+
+```style=background-color:gray
+[chempot1c2a]
+  type = KKSPhaseChemicalPotential
+  variable = Cijk
+  ka = ajk
+  fa_name = Fj
+  cb = Cij'k
+  kb = aj'k
+  fb_name = Fj'
+[]
+```
+
+Note that in this kernel the `ajk` and `aj'k` must be true site fractions ranging from 0 to 1.
+
+## Mass transport
+
+The global concentration variables $c_i$ govern the mass transport along
+chemical potential gradients.
+
+\begin{equation}
+   \frac{\partial c_i}{\partial t} = \nabla\cdot D_i\sum_j h_j\sum_k a_{jk}\nabla c_{ijk}. \label{eq:chdiff}
+\end{equation}
+
+This equation can be implemented using multiple `MatDiffusion` kernels in MOOSE.
+To illustrate this we can re-order the summation to yield
+
+\begin{equation}
+   \frac{\partial c_i}{\partial t} = \sum_j\sum_k \nabla\cdot \underbrace{D_i h_j a_{jk}}_{D'_i}\nabla c_{ijk}. \label{eq:chdiff2}
+\end{equation}
+
+This means we need to add a `MatDiffusion` kernel for sublattice in all phases,
+each operating on the variable $c_i$. We use the
+[!param](/Kernels/MatDiffusion/v) parameter to  specify $c_{ijk}$ as the
+variable to take the gradient of (rather than $c_i$). The effective diffusivity
+$D'_i$ passed into the kernel is $D\cdot h_j \cdot a_{jk}$, which can be provided
+by a [`DerivativeParsedMaterial`](DerivativeParsedMaterial.md)
+
+```style=background-color:gray
+[D'i]
+  type = DerivativeParsedMaterial
+  f_name = D'i
+  function = Di*hi*ajk
+  material_property_names = 'Di hi'
+  constant_names = ajk
+  constant_expressions = 1/3
+[]
+```
+
+Here we assume the site fraction of the current sublattice to be $1/3$.
+
+!alert note title=Extension to phase/sublattice dependent diffusivities
+The derivation should be check to see if a phase dependent concentration $D_{ij}$ or even sublattice dependent concentration $D_{ijk}$ is feasible.
+
+
+## Example
+
+An open TDB file for the Fe-Cr system was graciously provided by Aurélie Jacob
+at TU Wien, based on the work in [!cite](jacob2018revised).
+
+The TDB file can be converted into MOOSE [`DerivativeParsedMaterial`](DerivativeParsedMaterial.md)
+syntax using the [`free_energy.py` script](/CALPHAD.md).
+
+```
+mamba activate pycalphad
+cd $MOOSE_DIR/modules/phase_field/examples/slkks
+../../../../python/calphad/free_energy.py CrFe_Jacob.tdb BCC_A2 SIGMA
+```
+
+The example directory also contains an Jupyter notebook to visualize the phase
+diagram for the supplied Fe-Cr thermodynamic database file and the simulation
+results,
+
+!listing modules/phase_field/examples/slkks/CrFe.i

modules/phase_field/doc/content/modules/phase_field/index.md

@@ -29,6 +29,7 @@ MOOSE provides capabilities that enable the easy development of multiphase field
 
 - [phase_field/MultiPhase/WBMTwoPhase.md]: Two phases, one phase order parameter
 - [phase_field/MultiPhase/KKS.md]: per-phase concentrations, two phases
+- [phase_field/MultiPhase/SLKKS.md]: per-phase concentrations, per-sublattice concentrations, multiple phases
 - [phase_field/MultiPhase/WBM.md]: $N$ phases, $N$ phase order parameters
 - [Grand Potential Model](/GrandPotentialKernelAction.md): solving a Legendre transform of the phase field equations, where the independent variable is the chemical potential
 

modules/phase_field/doc/content/source/kernels/SLKKSChemicalPotential.md

@@ -0,0 +1,16 @@
+# SLKKSChemicalPotential
+
+!syntax description /Kernels/SLKKSChemicalPotential
+
+Implements
+\begin{equation}
+\frac 1{a_{jk}} \frac{\partial F_j}{\partial c_{ijk}} = \frac 1{a_{jk'}} \frac{\partial F_j}{\partial c_{ijk'}}
+\end{equation}
+
+where $c_{ijk}$ and $c_{ijk}$ are two sublattice concentrations (for sublattices $k$ and $k'$) in the same phase $j$, for component $i$.
+
+!syntax parameters /Kernels/SLKKSChemicalPotential
+
+!syntax inputs /Kernels/SLKKSChemicalPotential
+
+!syntax children /Kernels/SLKKSChemicalPotential

modules/phase_field/doc/content/source/kernels/SLKKSMultiACBulkC.md

@@ -0,0 +1,17 @@
+# SLKKSMultiACBulkC
+
+!syntax description /Kernels/SLKKSMultiACBulkC
+
+Implements the weak form
+\begin{equation}
+\left( -M\frac1{a_{jk}}\frac{\partial F_j}{\partial c_ijk} \sum_j \frac{\partial h_j}{\partial \eta_j}\sum_k c_{ijk}a_{jk},\psi\right)
+\end{equation}
+
+where $c_i$ is the phase concentration for phase $i$ and $h_i$ is the interpolation
+function for phase $i$ defined in [!cite](Folch05) (referred to as $g_i$ there, but we use $h_i$ to maintain consistency with other interpolation functions in MOOSE). Since in the KKS model, chemical potentials are constrained to be equal at each position, $\frac{\partial F_1}{\partial c_1} = \frac{\partial F_2}{\partial c_2} = \frac{\partial F_3}{\partial c_3}$.
+
+!syntax parameters /Kernels/SLKKSMultiACBulkC
+
+!syntax inputs /Kernels/SLKKSMultiACBulkC
+
+!syntax children /Kernels/SLKKSMultiACBulkC

modules/phase_field/doc/content/source/kernels/SLKKSMultiPhaseConcentration.md

@@ -0,0 +1,22 @@
+# SLKKSMultiPhaseConcentration
+
+!syntax description /Kernels/SLKKSMultiPhaseConcentration
+
+For a sub lattice KKS (SLKKS) model with $n$ phases, the phase concentration
+constraint equation is
+\begin{equation}
+ c_i = \sum_j h_j\sum_k a_{jk} c_{ijk},
+\end{equation}
+where $i$ indexes a component, $j$ a phase, and $k$ a sublattice in the given
+phase. $a_{jk}$ is the fraction of $k$ sublattice sites in phase $j$. Thus $c_i$
+is the global concentration of component $i$ and $h_j$ is the interpolation
+function for phase $j$.
+
+The version of this class for a two phase model with a single switching function
+is [SLKKSPhaseConcentration](/SLKKSPhaseConcentration.md).
+
+!syntax parameters /Kernels/SLKKSMultiPhaseConcentration
+
+!syntax inputs /Kernels/SLKKSMultiPhaseConcentration
+
+!syntax children /Kernels/SLKKSMultiPhaseConcentration

modules/phase_field/doc/content/source/kernels/SLKKSPhaseConcentration.md

@@ -0,0 +1,22 @@
+# SLKKSPhaseConcentration
+
+!syntax description /Kernels/SLKKSPhaseConcentration
+
+This class is the two-phase version of
+[SLKKSMultiPhaseConcentration](/SLKKSMultiPhaseConcentration.md).
+
+For a sub lattice KKS (SLKKS) model with two phases $\alpha$ and $\beta$, the
+phase concentration constraint equation is
+\begin{equation}
+ c_i = h\sum_k a_{\alpha k} c_{i\alpha k} + (1-h)\sum_k a_{\beta k} c_{i\beta k}
+\end{equation}
+where $i$ indexes a component and $k$ a sublattice in the given phase.
+$a_{\alpha k}$ and $a_{\beta k}$ is the fraction of $k$ sublattice sites in
+phases $\alpha$ and $\beta$. Thus $c_i$ is the global concentration of component
+$i$ and $h$ is the phase switching function.
+
+!syntax parameters /Kernels/SLKKSPhaseConcentration
+
+!syntax inputs /Kernels/SLKKSPhaseConcentration
+
+!syntax children /Kernels/SLKKSPhaseConcentration

modules/phase_field/doc/content/source/kernels/SLKKSSum.md

@@ -0,0 +1,29 @@
+# SLKKSSum
+
+!syntax description /Kernels/SLKKSSum
+
+## Overview
+
+This kernel is used if the simulation only contains a single phase with multiple
+sublattices. Otherwise
+[`SLKKSMultiPhaseConcentration`](SLKKSMultiPhaseConcentration.md) needs to be used.
+
+One potential application is a static solve for the sublattice concentrations
+(and free energy) of a single phase. This information can be used to tabulate
+sublattice concentrations as functions of total concentrations to inform
+initial conditions of full solve simulations.
+
+### See also
+
+- `modules/phase_field/examples/slkks/CrFe_sigma`  for an example input that tabulates the sublattice concentrations of the Fe-Cr sigma phase
+- `modules/phase_field/examples/slkks/CrFe`  for an example input for a full solve that uses the tabulated sublattice concentrations to set the initial conditions
+
+## Example Input File Syntax
+
+!! Describe and include an example of how to use the SLKKSSum object.
+
+!syntax parameters /Kernels/SLKKSSum
+
+!syntax inputs /Kernels/SLKKSSum
+
+!syntax children /Kernels/SLKKSSum

modules/phase_field/doc/content/source/postprocessors/WeightedVariableAverage.md

@@ -0,0 +1,20 @@
+# WeightedVariableAverage
+
+!syntax description /Postprocessors/WeightedVariableAverage
+
+Compute the ratio $a$ of volume integrals
+\begin{equation}
+a=\frac{\int_Omega w\cdot v dr}{\int_Omega w dr},
+\end{equation}
+where $v$ (`v`) is a coupled variable and $w$ (`weight`) a material property.
+
+For constant weight values $w$ this object is equivalent to
+[ElementAverageValue](/ElementAverageValue.md).
+
+!syntax parameters /Postprocessors/WeightedVariableAverage
+
+!syntax inputs /Postprocessors/WeightedVariableAverage
+
+!syntax children /Postprocessors/WeightedVariableAverage
+
+!bibtex bibliography