Mobility
Both precipitation and diffusion modeling in kawin depend on knowing the mobility of a component in a phase. For precipitation, this enables precipitate growth rate calculations, and for diffusion, diffusive fluxes can be computed. kawin supports the usage of Calphad-based mobility models to determine both the mobility and diffusivity. This is handled by an extension of pycalphad’s Model class. The following equation shows how the mobility is determined by two terms: MF and MQ, where both are treated as Redlich-Kister polynomials. These two terms can be stored in a thermodynamic database (.tdb) file as MF and MQ.
$$ M = \frac{\theta}{RT} \exp{\left(-\frac{Q}{RT}\right)} = \frac{1}{RT} \exp{\left(M_F + M_Q\right)} $$
In addition, diffusivity can be directly stored in .tdb files as DF and DQ. Note that kawin only supports cross-diffusion effects if the mobility parameters (MF and MQ) are used.
$$ D = D_0 \exp{\left(-\frac{Q}{RT}\right)} = \exp{\left(D_F + D_Q\right)} $$
Interdiffusivity
If using mobility terms, the fluxes in a lattice fixed frame ($J_k$) is related to the atomic mobility and the chemical potential gradient. Here, we’ll assume all components are substitional such that $ u_k = x_k $ (see Diffusion Overview for more on u-fractions).
$$ J_k = -x_k M_k \frac{\partial \mu_k}{\partial z} $$
In volume fixed frame, the sum of fluxes must be 0 or $ \sum_{k}{J_k^v} = 0 $. The fluxes in a volume fixed frame is related to that of the lattice fixed frame by:
$$ J_k^v = J_k - x_k \sum_{i}{J_i} = -\sum_{i}{(\delta_{ki} - x_k) x_i M_i \frac{\partial \mu_i}{\partial z}} $$
Expanding the chemical potential gradient using chain-rule derivatives, we obtain the following where we can define $ L_{ki} $ and $ D_{kj} $ to simplify this equation.
$$ J_k^v = -\sum_{i}{(\delta_{ki} - x_k) x_i M_i \sum_{j}{\frac{\partial \mu_i}{\partial x_j} \frac{\partial x_j}{\partial z}}} $$
$$ L_{ki} = (\delta_{ki} - x_k) x_i M_i $$
$$ D_{kj} = \sum_{i}{L_{ki} \frac{\partial \mu_i}{\partial x_j}} $$
$$ J_k^v = -\sum_{j}{D_{kj} \frac{\partial x_j}{\partial z}} $$
Due to the summation constraint, there are only n-1 independent components. Thus, a reference element (n) can be defined and the interdiffusivities ($D_{kj}^n$) can be determined. Only n-1 components need to be tracked if interdiffusivities are used. Note that the definition for $D_{kj}^n$ is different whether k is a substitutional or interstitial component. This is due to the assumption that the interstitials do not contribute to the volume of the system.
$$ \sum_{j}{\frac{\partial x_j}{\partial z}} = 0 $$
$$ \frac{\partial x_n}{\partial z} = -\sum_{j}^{n-1}{\frac{\partial x_j}{\partial z}} $$
$$ J_k^v = -\left( \sum_{j}^{n-1}{D_{kj} \frac{\partial x_j}{\partial z}} \right) - D_{kn} \frac{\partial x_n}{\partial z} = \left( -\sum_{j}^{n-1}{D_{kj} \frac{\partial x_j}{\partial z}} \right) - \left( -\sum_{j}^{n-1}{D_{kn} \frac{\partial x_j}{\partial z}} \right) $$
$$ J_k^v = -\sum_{j}^{n-1}{D_{kj}^n \frac{\partial x_j}{\partial z}} $$
Where
$$ D_{kj}^n = D_{kj} - D_{kn} = \sum_{i}{L_{ki} \left(\frac{\partial \mu_i}{\partial x_j} - \frac{\partial \mu_i}{\partial x_n}\right)} $$
When interstitials are present, it is assumed that interstitials don’t contribute to the total volume. Working in u-fractions, only the substitional components ($S$) have to fulfill the summation constraints for the composition and fluxes.
$$ \sum_{k \epsilon S}{u_k} = 1 $$
$$ \sum_{k \epsilon S}{J_k^v} = 0 $$
This changes the definition of the fluxes and interdiffusivities. For substitutional components, no change is needed for the interdiffusivities apart from excluding the interstitials:
$$ J_k^v = J_k - u_k \sum_{j \epsilon S}{J_j} $$
$$ L_{ki} = (\delta_{ki} - u_k) u_i M_i \text { for } i \epsilon S $$
$$ D_{kj} = \sum_{i \epsilon S}{L_{ki} \frac{\partial \mu_i}{\partial u_j}} $$
$$ D_{kj}^n = D_{kj} - D_{kn} $$
For interstitials (c), the u-fraction is not bound by the summation constraint. The vacancy site fraction in the interstitial sublattice ($y_{va}$) is also included to account for the average number of available jump sites of an interstitial increasing with increasing vacancy content. For most phases, $y_{va} \approx 1 $ on the interstitial sublattice. However, for phases such as oxides or carbides, the vacancy content is near 0 and $ y_{va} M_c $ can be lumped into a single term. Note that this is implicitly done for the substitional species since the vacancy site fraction in the substitutional sublattice is near 0 so the effects of vacancies is already included in $ M_k $.
$$ J_c^v = J_c = -y_{va} u_c M_c \frac{\partial \mu_c}{\partial z} $$
$$ D_{cj}^n = D_{cj} = y_{va} u_c M_c \frac{\partial \mu_c}{\partial u_c} $$
In kawin, lumping the $ y_{va} $ and $ M_c $ can be done for specific phases by GeneralThermodynamics.vacancyPoorInterstitialSublattice['phase'] = True (which will ignore $ y_{va}$ and assume that vacancies are accounted for in $ M_c $). Note that the mobility parameters for these phases must also be assessed with this assumption.
Chemical potential gradient (free energy curvature)
In the above equations, calculating the interdiffusivity from atomic mobility requires the chemical potential gradient ($\frac{\partial \mu_i}{\partial x_j}$) which can be determined from the curvature of the free energy surface.
$$ \frac{\partial \mu_i}{\partial x_j} = \frac{\partial^2 G_m^\alpha}{\partial x_i \partial x_j} $$
For models representing solid solutions with a single sublattice, the curvature of the Gibbs free energy can be easily found. However, for most other models, the chemical potential gradient must be found numerically. This is done by first solving for equilibrium, then using chain rule differentiation to get the derivative with respect to any parameter.
Solving equilibrium for a single phase system is done by minimizing the Gibbs free energy subjected to mass balance and summation constraints ($\boldsymbol g$).
$$ \mathcal{L} = n^\alpha G_m^\alpha + \sum_{i}{\Lambda_i g_i} $$
For a single phase system, assuming no charged species, the Langragian becomes the following:
$$ \mathcal{L} = n^\alpha G_m^\alpha + \sum_{A}{\mu_A \left(x_A - n^\alpha \frac{M_A^\alpha}{M}\right)} + \sum_{s}{\lambda^s \left(1 - \sum_{i \epsilon s}{y_i^s} \right)} $$
The terms in the above equation can be divided into two groups: the fixed parameters (z) and the degrees of freedom ($\boldsymbol\theta$). The fixed parameters include temperature (T), pressure (P) and composition ($x_A$). The degrees of freedom include the site fractions ($y_i^s$), phase fraction ($n^\alpha$), chemical potential ($\mu_A$) and summation constraint multiplier ($\lambda^s$).
The $\boldsymbol\theta$ group will be further divided into $\boldsymbol{\mathcal{v}}$ (which includes $y_i^s$ and $n^\alpha$) and $\boldsymbol \Lambda$ (which includes the Langragian multipliers, $\mu_A$ and $\lambda^s$).
Solving for equilibrium involves finding the solution to the following equations.
$$ \frac{\partial \mathcal{L}}{\partial \mathcal{v}_i} = \frac{\partial}{\partial \mathcal{v}_i} (n^\alpha G_m^\alpha) + \sum_i{\Lambda_j \frac{\partial g_j}{\partial \mathcal{v}_i}} = 0 $$
$$ \frac{\partial \mathcal{L}}{\partial \Lambda_j} = g_j = 0 $$
Once equilibrium is found, a system of equations can be constructed using chain rule derivative to find the derivatives of $\boldsymbol\theta$ with respect to z.
$$ \sum_{j}{\frac{\partial^2 \mathcal{L}}{\partial \theta_i \partial \theta_j} \frac{\partial \theta_j}{\partial z_k}} = \frac{\partial^2 \mathcal{L}}{\partial \theta_i \partial z_k} $$
$$ \left[ \begin{array}{cc} \frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{\mathcal{v}}^2} & \frac{\partial \boldsymbol{g}}{\partial \boldsymbol{\mathcal{v}}} \newline \left(\frac{\partial \boldsymbol{g}}{\partial \boldsymbol{\mathcal{v}}} \right)^T & 0 \end{array} \right] \left[ \begin{array}{c} \frac{\partial \boldsymbol{\mathcal{v}}}{\partial z} \newline \frac{\partial \boldsymbol{\Lambda}}{\partial z} \end{array} \right] = \left[ \begin{array}{c} -\sum{\Lambda_i \frac{\partial^2 g_i}{\partial \boldsymbol{v} \partial z}} \newline \frac{\partial \boldsymbol{g}}{\partial z} \end{array} \right] $$
Single phase systems
For a single phase system, the above system of equations become
$$ \left[ \begin{array}{cccc} \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y}^2} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right) \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right)^T & 0 & 0 & \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right) \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right)^T & 0 & 0 & 0 \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right)^T & \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right)^T & 0 & 0 \end{array} \right] \left[ \begin{array}{c} \frac{\partial \boldsymbol{y}}{\partial x_i} \newline \frac{\partial n^\alpha}{\partial x_i} \newline \frac{\partial \boldsymbol{\lambda}}{\partial x_i} \newline \frac{\partial \boldsymbol{\mu}}{\partial x_i} \end{array} \right] = \left[ \begin{array}{c} 0 \newline \vdots \newline 0 \newline \frac{\partial x_A}{\partial x_i} \newline \vdots \newline \frac{\partial x_N}{\partial x_i} \end{array} \right] $$
Where
$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y}^2} \right)_{ij} = n^\alpha \frac{\partial^2 G_m^\alpha}{\partial y_i \partial y_j} $$
$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right)_i = \frac{\partial G_m^\alpha}{\partial y_i} - \sum_A{\mu_A \frac{\partial M_A}{\partial y_i}} $$
$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right)_{ij} = -1 \text{ if } y_i \epsilon \text{ sublattice j else } 0 $$
$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right)_{ij} = -n^\alpha \frac{\partial M_A}{\partial y_i} $$
$$ \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right)_j = -M_j^\alpha $$
$$ \frac{\partial x_A}{\partial x_i} = \delta_{Ai} $$
It will be noted that this is to be done using terms normalized to the moles of formula units instead of moles of atoms.
Defining the chemical potential gradient
The FreeEnergyHessian module in kawin provides two methods: dMudX and partialdMudX. These methods give different definitions of the derivative of the chemical potential. Additionally, it will give different values than if the chemical potential gradient is determined using a finite difference approach using thermodynamic software. The partialdMudX method outputs a ($\text{n x n}$) matrix of the partial derivatives of the chemical potential with respect to the composition, treating all components as independent. Using a finite difference method to calculate the chemical potential gradient will calculate the total derivative (with reference element n) since most thermodynamic software only require n-1 components to be defined. This will give a ($\text{n x n-1}$) matrix. The dMudX method gives a ($\text{n-1 x n-1}$) matrix, where each index is the total derivative of the chemical potential with reference element n subtracted by the total derivative of the chemical potential of the reference element. This is more equivalent to the curvature (Laplacian) of the Gibbs free energy surface with the reference element. This method is required for finding $M^{-1}$ for the multicomponent growth rate.
$$ \text{partialdMudX} = \frac{\partial \mu_i}{\partial x_j} $$
$$ \text{dMudX} = \frac{\partial}{\partial x_j} (\mu_i - \mu_n) = (\frac{\partial \mu_i}{\partial x_j} - \frac{\partial \mu_i}{\partial x_n}) - (\frac{\partial \mu_n}{\partial x_j} - \frac{\partial \mu_n}{\partial x_n}) $$
$$ \text{finite difference} = \frac{\partial \mu_i}{\partial x_j} - \frac{\partial \mu_i}{\partial x_n} $$
References
A. Borgenstam, A. Engstrom, L. Hoglund and J. Agren, “DICTRA, a tool for simulation of diffusional transformations in alloys” Journal of Phase Equilibria 21 (2000) p. 269
U. R. Kattner and C. E. Campbell, “Modelling of thermodynamics and diffusion in multicomponent systems” Materials Science and Technology 25 (2009) p. 443
H. Larsson and M. Jansson, “Rate of change at equilibrium” Calphad: Computer Coupling of Phase Diagrams and Thermochemistry 51 (2015) p. 220