Grain boundary structure–property model inference using polycrystals: the overdetermined case

Abstract

Efforts to construct predictive grain boundary (GB) structure–property models have historically relied on property measurements or calculations made on bicrystals. Experimental bicrystals can be difficult or expensive to fabricate, and computational constraints limit atomistic bicrystal simulations to high-symmetry GBs (i.e., those with small enough GB periodicity). Although the use of bicrystal property data to construct GB structure–property models is more direct, in many experimental situations the only type of data available may be measurements of the effective properties of polycrystals. In this work, we investigate the possibility of inferring GB structure–property models from measurements of the homogenized effective properties of polycrystals when the form of the structure–property model is unknown. We present an idealized case study in which GB structure–property models for diffusivity are inferred from noisy simulation results of two-dimensional microstructures, under the assumption that the number of polycrystal measurements available is larger than the number of parameters in the inferred model. We also demonstrate how uncertainty quantification for the inferred structure–property models is easily performed within this framework.

Appendices

Appendix A: Derivation of diffusivity homogenization equation

The derivation of Eq. 1 is related to the idea of [40] that one can obtain the effective diffusivity from Fick’s First law by considering the flux across some “slice” through the microstructure. However, in contrast to [40], we are interested in the effective diffusivity of the entire GBN, rather than the average diffusivity of individual GBs in the network; the two ideas are, nevertheless, intimately related (see footnote 1 in [26]). Our approach to calculate the effective diffusivity of the entire GBN is equivalent to an adaptation of the finite volume method (see [27]). For convenience, we set the concentration at the source node (node a) to be $c_0$ and consider the total flux ($J_b = Q_b/A$) arriving at the sink node (node b), whose concentration we fix at 0 (Neumann boundary conditions). The effective diffusivity is then

$$\begin{aligned} \bar{D}^{\text{pred}} = \frac{L}{c_0} \frac{Q_b}{A} \end{aligned}$$

(17)

and the mass flow rate arriving at the sink is simply

$$\begin{aligned} Q_b = \sum _{i \sim b} -\frac{D_{ib} A_{ib}}{L_{ib}} \left( c_b - c_i\right) \end{aligned}$$

(18)

where $D_{ib}$, $A_{ib}$, and $L_{ib}$ are, respectively, the diffusivity, cross-sectional area, and length of the edge connecting node i to node b, $c_k$ is the concentration at node k, and the sum is over all edges in the GBN that are incident to node b ($i \sim b$ denotes that there is an edge between nodes i and b). Making use of the definition of the GBN Laplacian matrix ($\mathcal{L}$) [27], this can be rewritten as

$$\begin{aligned} Q_b = -\mathcal{L}_b^{\intercal } c \end{aligned}$$

(19)

where $\mathcal{L}_b$ is the bth column of $\mathcal{L}$ and c is a vector containing the concentration at every node, which is calculated by solving the linear system of equations defined by

$$\begin{aligned} \widehat{\mathcal{L}} c = c_0 e_b \end{aligned}$$

(20)

where $\widehat{\mathcal{L}}$ is a matrix formed by replacing the ath and bth rows of $\mathcal{L}$ with $e_a^{\intercal }$ and $e_b^{\intercal },$ respectively, and where $e_k$ is a vector having the kth element equal to one and all others equal to zero. We note that Eq. 20 is simply a statement of conservation of mass at every node. Substituting the solution of Eq. 20 into Eq. 19 and substituting the result into Eq. 17, we arrive at Eq. 1. We note also that the direct inversion of $\widehat{\mathcal{L}}$ is not necessary for the solution of Eq. 20, and is not computationally efficient; rather, standard methods for the solution of linear systems of equations will be preferred. In our implementation of Eq. 1, we employ MATLAB’s mldivide operator.

Appendix B: Inverse problem theory

Here we give a brief overview of relevant aspects of Tarantola’s approach to Bayesian inverse problem theory, and then demonstrate its application to the problem of GB property localization. Inverse problem theory is a method of inferring model parameters $(\mathbf{x})$ that characterize a system using the results of some measurements/observations of the system $(\mathbf{y})$ [42,43,44,45,46]. In a given system, ${\mathbf{x}}=\left\{ x_1,x_2,\dots \right\}$ is a set containing the independent parameters and ${\mathbf{y}}=\left\{ y_1,y_2,\dots \right\}$ is a set containing the dependent parameters, both of which we may only hope to know with some imperfect degree of certainty. Tarantola⁵ proposed that our state of information (what we know about $\mathbf{y}$ and $\mathbf{x}$) can be described by a PDF, called the a posteriori state of information, $\sigma (\mathbf{y},\mathbf{x})$, which is equal to the conjunction of the a priori state of information, $\rho (\mathbf{y},\mathbf{x})$, and the theoretical state of information, $\Theta (\mathbf{y},\mathbf{x})$ [44]. The a priori state of information is what we know before ever making any observations and may represent some known physical constraints. The theoretical state of information encodes correlations between $\mathbf{x}$ and $\mathbf{y}$ resulting from a homogenization or other physical theory and corresponding uncertainty. Using the Kolmogorov axioms, Tarantola and Valette showed that the a posteriori state of information is given by [42]:

$$\begin{aligned} \sigma ({\mathbf{y},\mathbf{x}})=k\frac{\rho ({\mathbf{y},\mathbf{x}})\Theta (\mathbf{y},\mathbf{x})}{\mu (\mathbf{y},\mathbf{x})} \end{aligned}$$

(21)

Here, k is a normalization constant, and $\mu (\mathbf{y},\mathbf{x})$ is the homogeneous state of information, which is the PDF that assigns a probability to each region of the parameter space that is equal to the volume of that region [43]. Functionally, $\mu (\mathbf{y},\mathbf{x})$ represents a state of complete ignorance (i.e., the absence of any a priori information), so that the ratio in Eq. 21 can be interpreted as quantifying how much is known relative to a state of complete ignorance. In the present context, Eq. 1 represents the forward problem, where the dependent parameter is the observed effective diffusivity $\mathbf{y}=\{\bar{D}^{\text{obs}}\}$ and the independent parameters are the structure–property model and sample microstructure $\mathbf{x}=\{\mathscr{D},M\}$. As described in the “GB property localization” section, GB property localization will typically leverage information from multiple samples, so that we have ${\mathbf{M}}=\{M_1,M_2,\ldots ,M_N\}$ and ${\bar{\mathbf{D}^{\text{obs}}}}=\{\bar{D}_1^{\text{obs}},\bar{D}_2^{\text{obs}},\ldots ,\bar{D}_N^{\text{obs}}\}$. Ignoring the normalization constant, we can then rewrite Eq. 21 as:

$$\begin{aligned} \sigma ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})\propto \frac{\rho ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})\Theta ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})}{\mu ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})} \end{aligned}$$

(22)

The resolution of the inverse problem consists in identifying the structure–property model, $\mathscr{D}$, that is most probable given our observations of ${\bar{\mathbf{D}^{\text{obs}}}}$ and ${\mathbf{M}}$. This is accomplished by integrating Eq. 22 to compute the a posteriori state of information over the space of independent parameters:

$$\begin{aligned} \sigma ({\mathscr{D},\mathbf{M}})=\int{\sigma ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}}) d(\bar{\mathbf{D}^{\text{obs}}})} \end{aligned}$$

(23)

The evaluation of this integral is facilitated by considering relevant simplifications. Because the a priori information about $\mathbf{M}$ and $\mathscr{D}$ is not obtained from measurements of ${\bar{\mathbf{D}^{\text{obs}}}}$, their states of information are independent [44], which implies that:

$$\begin{aligned}&\rho ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})=\rho _{\bar{\mathbf{D}^{\text{obs}}}}({\bar{\mathbf{D}^{\text{obs}}}})\rho _{\{\mathscr{D},\mathbf{M}\}}\left({\mathscr{D},\mathbf{M}}\right) \end{aligned}$$

(24)

$$\begin{aligned}&\mu ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})=\mu _{\bar{\mathbf{D}^{\text{obs}}}}({\bar{\mathbf{D}^{\text{obs}}}})\mu _{\{\mathscr{D},\mathbf{M}\}}\left({\mathscr{D},\mathbf{M}}\right) \end{aligned}$$

(25)

It should be noted here that subscripts are used to distinguish the functions.

Let us assume that Eq. 1, as the physical theory relating the independent and dependent parameters, is at most mildly⁶ nonlinear. Combining this assumption with the Kolmogorov definition for conditional probability [55], and taking the homogeneous probability of the independent parameters as their marginal probability, $\Theta ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})$ can be written, according to the treatment of Tarantola and Vallete [42, 44], as:

$$\begin{aligned} \Theta ({\bar{\mathbf{D}^{\text{obs}}},\mathscr{D},\mathbf{M}})=\theta \left({\bar{\mathbf{D}^{\text{obs}}}}\mid{\{\mathscr{D},\mathbf{M}\}}\right) \mu _{\{\mathscr{D},\mathbf{M}\}}\left({\mathscr{D},\mathbf{M}}\right) \end{aligned}$$

(26)

Substituting Eqs. 22 and 24–26 into Eq. 23, we obtain:

$$\begin{aligned} \sigma \left({\mathscr{D},\mathbf{M}}\right) \propto \rho _{\{\mathscr{D},\mathbf{M}\}}\left({\mathscr{D},\mathbf{M}}\right) \int{\frac{\rho _{\bar{\mathbf{D}^{\text{obs}}}}({\bar{\mathbf{D}^{\text{obs}}}})\theta \left({\bar{\mathbf{D}^{\text{obs}}}}\mid{\{\mathscr{D},\mathbf{M}\}}\right) }{\mu _{\bar{\mathbf{D}^{\text{obs}}}}({\bar{\mathbf{D}^{\text{obs}}}})} d({\bar{\mathbf{D}^{\text{obs}}}})} \end{aligned}$$

(27)

Because the manifold that ${\bar{\mathbf{D}^{\text{obs}}}}$ inhabits is a linear space, using the definition of homogeneous probability distribution presented by Mosegaard and Tarantola [43], $\mu ({\bar{\mathbf{D}^{\text{obs}}}})$ is constant. We also make the simplifying approximation that any uncertainty in Eq. 1 is negligible, which implies that $\theta \left({\bar{\mathbf{D}^{\text{obs}}}}\mid{\{\mathscr{D},\mathbf{M}\}}\right) =\delta ({\bar{\mathbf{D}^{\text{obs}}}}-{\bar{\mathbf{D}^{\text{pred}}}}(\mathscr{D},{\mathbf{M}}))$. Under these conditions, the integration operation in Eq. 27 results in:

$$\begin{aligned} \sigma ({\mathscr{D},\mathbf{M}})\propto \rho _{\{\mathscr{D},\mathbf{M}\}}\left({\mathscr{D},\mathbf{M}}\right) \rho _{\bar{\mathbf{D}^{\text{obs}}}}\left({\bar{\mathbf{D}^{\text{pred}}}}(\mathscr{D},{\mathbf{M}})\right) \end{aligned}$$

(28)

which is the result given in Eq. 3. The last term, $\rho _{\bar{\mathbf{D}^{\text{obs}}}}\left({\bar{\mathbf{D}^{\text{pred}}}}(\mathscr{D},{\mathbf{M}})\right)$, is a likelihood function, which quantifies how well the model explains the data. In other words, $\rho _{\bar{\mathbf{D}^{\text{obs}}}}\left({\bar{\mathbf{D}^{\text{pred}}}}(\mathscr{D},{\mathbf{M}})\right)$ quantifies how well the independent parameters explain the dependent parameters. Thus, for the conditions presented here, the a posteriori state of information about the independent parameters is proportional to the product of the likelihood function and the a priori state of information about the independent parameters.