Neural Network Accelerated Process Design of Polycrystalline
  Microstructures
Computational experiments are exploited in finding a well-designed processing
path to optimize material structures for desired properties. This requires
understanding the interplay between the processing-(micro)structure-property
linkages using a multi-scale approach that connects the macro-scale (process
parameters) to meso (homogenized properties) and micro (crystallographic
texture) scales. Due to the nature of the problem's multi-scale modeling setup,
possible processing path choices could grow exponentially as the decision tree
becomes deeper, and the traditional simulators' speed reaches a critical
computational threshold. To lessen the computational burden for predicting
microstructural evolution under given loading conditions, we develop a neural
network (NN)-based method with physics-infused constraints. The NN aims to
learn the evolution of microstructures under each elementary process. Our
method is effective and robust in finding optimal processing paths. In this
study, our NN-based method is applied to maximize the homogenized stiffness of
a Copper microstructure, and it is found to be 686 times faster while achieving
0.053% error in the resulting homogenized stiffness compared to the traditional
finite element simulator on a 10-process experiment.
Junrong Lin, Mahmudul Hasan, Pinar Acar, Jose Blanchet, Vahid Tarokh
# Neural Network Accelerated Process Design of Polycrystalline Microstructures

###### Abstract

Computational experiments are exploited in finding a well-designed processing path to optimize material structures for desired properties. This requires understanding the interplay between the processing-(micro)structure-property linkages using a multi-scale approach that connects the macro-scale (process parameters) to meso (homogenized properties) and micro (crystallographic texture) scales. Due to the nature of the problem's multi-scale modeling setup, possible processing path choices could grow exponentially as the decision tree becomes deeper, and the traditional simulators' speed reaches a critical computational threshold. To lessen the computational burden for predicting microstructural evolution under given loading conditions, we develop a neural network (NN)-based method with physics-infused constraints. The NN aims to learn the evolution of microstructures under each elementary process. Our method is effective and robust in finding optimal processing paths. In this study, our NN-based method is applied to maximize the homogenized stiffness of a Copper microstructure, and it is found to be 686 times faster while achieving 0.053% error in the resulting homogenized stiffness compared to the traditional finite element simulator on a 10-process experiment.
## 1 Introduction

The research on investigating process-structure-property relationships has become more prevalent with the introduction of the Integrated Computational Materials Engineering (ICME) paradigm [1], which aims to solve complex and multi-scale material design problems. As a result of the recent advancements following the introduction of ICME, several aspects of computational materials science and process engineering have significantly improved. For instance, new methodologies are developed for ICME to lower the costs and risks associated with the processing of new materials [2; 3]. Computational experiments are exploited to eliminate the traditional trial-error approach in bridging material features and properties. The investigation of the process-(micro)structure-property linkages requires a multi-scale approach that connects the macro-scale (process parameters) to meso (homogenized properties) and micro (crystallographic texture) scales. Physics-based multi-scale modeling of materials provides an opportunity to achieve the optimum design of materials with enhanced properties. However, their high computational cost prevents these multi-scale approaches from being widely adopted and used by industry in real material design efforts [4]. Nevertheless, the process-structure problem is studied less than the structure-property problem as the physics behind the process-structure problem is generally more complicated. Machine learning tools have demonstrated promise to address this gap by building low-cost process-structure-property surrogate models to replace computationally expensive physics-based models [5].

Control of polycrystalline microstructures is important in material design, processing, and quality control since the orientation-dependent material properties
(e.g., stiffness) could change as the underlying microstructures evolve during a deformation process. Therefore, identification of the optimal processing route to produce a material with desired texture and properties plays an important role in materials research. Traditional experimental approaches used to explore the optimum processing routes for a given texture are usually based on trial and error and, thus, can be tedious and expensive. Hence, computational methods are developed to replace these experiments to accelerate the design of microstructures with desired textures. Fast and accurate prediction of texture evolution during processing can significantly contribute to linking the current material design and manufacturing efforts for polycrystalline materials.

Current research in this field is focused on process-structure [6; 7; 8; 9; 10; 11] or structure-property [12; 13; 14; 15; 16; 17; 18] relationships separately. For example, Sarkar et al. [6] developed a surrogate model for \(ZrO_{2}\)-toughened \(Al_{2}O_{2}\) ceramics to predict the sinter density and grain size from sintering heat treatments process parameters. In another study, Tapia et al. [8] established a Gaussian process regression-based surrogate model for the heat treatment process of NiTi shape-memory alloys where the input parameters are heat treatment temperature and its duration, as well as initial nickel composition, and the output is the final nickel composition after heat treatment. Next, Acar and Sundararaghavan [12] developed a linear solver-based multi-scale approach using reduced-order modeling to design target microstructures that optimize homogenized material properties. In order to learn the reduced basis functions that can adequately represent the crystallographic texture of polycrystalline materials, Achargee et al. [13] and Ganapathysubramanian et al. [14] used the proper orthogonal decomposition (POD) and method of snapshots in Rodrigues orientation space [19] to cre
ate a continuum sensitivity-based optimization technique to compute the material properties that are sensitive to the microstructure. Similarly, Kalidindi et al. [15] designed the microstructure for a thin plate with a circular hole at the center to maximize the uniaxial load-carrying capacity of the plate without plastic deformation. There are few other studies available in the literature that explore the process-structure and structure-property relationships for additively manufactured materials [20; 21; 22]. Very recently, Dornheim et al. [23] developed a model-free deep reinforcement learning algorithm to optimize the processing paths (up to 100 combinations) for a targeted microstructure. Their algorithm does not require prior samples rather it can connect with the processing simulations during optimization. They also extended the method to solve multi-objective optimization problems. In another study, Honarmandi et al. [24] proposed a novel framework based on batch Bayesian optimization to solve the inverse problem to find the material processing specifications using microstructure data. Using both low-fidelity and high-fidelity phase field models, they developed a Gaussian process regression-based surrogate model to replace the computationally expensive process models and integrated it into inverse design optimization.

Recent studies [11; 25] model the texture evolution of polycrystalline materials using a probabilistic representation called orientation distribution function (ODF), which describes the volume density of crystals of different orientations in a microstructure. Predicting the changes of the microstructural texture after applying a particular deformation process (e.g. applying shear force along a particular direction for 1 second) involves solving the ODF conservation equation, which is a differential equation [25; 26]. The traditional physics-based solution is developed using the finite element method to solve the conservation equation
numerically [26]. The change in the ODFs during processing also controls the homogenized (meso-scale) mechanical properties of the microstructures. Significant challenges arise when we carry out a large-scale search task to find out microstructures with desired orientation-dependent properties. Under this context, we expect to optimize some homogenized mechanical properties by sequentially applying different deformation processes. The optimal processing path could be

Figure 1: Schematic of the contribution of this study. Data-driven surrogate model is developed to replace the physics-based simulator on the process-structure-property problem.
found by searching algorithms coupled with a process simulator of microstructure evolution. However, in practice, even a simple task (e.g. given an arbitrary initial texture and 5 different deformation processes and their combinations, finding out the exact optimal path of applying 10 processes sequentially) could be time-consuming to solve using traditional physics-based simulators due to both simulation time requirements and exponentially growing number of possible processing paths. Concurrent multi-process modeling can further increase the computational complexity of the material design. Simulation speed becomes the bottleneck when confronting these tasks. Therefore, this study develops a supervised learning approach together with a local search algorithm to explore and bridge the process-microstructure-property relationships of polycrystalline materials.

The neural network has already been shown to be a revolutionary function approximator in this era of abundant data which is already adopted in the literature [27; 28; 29; 30; 31] to study process-structure-property linkages in order to build low-cost surrogate models. Inspired by physics-informed neural networks (PINN) [32], we will develop a surrogate neural network trained on a small dataset to replace the physics-based simulator in this process design task. The overall contribution of this study is summarized in Fig. 1. The high-speed inference feature of a pre-trained neural network could accelerate the algorithm and therefore enable large-scale searching. Compared to the finite element (FE) simulator, this method is faster and has a little trade-off in prediction accuracy. The organization of this article is as follows: Section 2 and Section 3 discuss the mathematical formulation for modeling of deformation processing and polycrystalline microstructures, respectively. Section 4 and Section 5 describe the physics-based and data-driven modeling of the process-microstructure relationship, respectively.
The performance of the neural network-based surrogate model of process design for improved mechanical performance is reported in Section 6. Finally, Section 7 provides a summary of the paper and a discussion on potential future work.

## 2 Mathematical Modeling of Deformation Process with ODF Approach

Multiple crystals with various crystallographic orientations make up a polycrystalline material, and these orientations determine the microstructural texture, which is mathematically described by the orientation distribution function (ODF). During a deformation process, the texture of polycrystalline microstructure changes under applied loads. Deformation process modeling with the ODF approach is computationally efficient compared to the expensive finite element solver. The ODF, denoted by \(A(\mathbf{r},t)\), indicates the volume density of the crystals in the orientation space, \(\mathbf{r}\) in a certain time \(t\). ODFs are discretized over the Rodrigues orientation space using finite element techniques. The reduced region, called the fundamental region \(\Omega\), is induced based on the crystallographic symmetry of the polycrystal system from the initial Rodrigues orientation space. The Rodrigues angle-axis parameterization method is used to depict the various crystal orientations. In contrast to the Euler angles representation [33; 34], this method uses axis-angle representations to express crystal orientations. The Rodrigues parameterization is described by scaling of the axis of rotation, \(\mathbf{n}\), as \(\mathbf{r}=\mathbf{n}tan(\frac{\theta}{2})\), where \(\theta\) is an angle of crystal rotation. For further information on Rodrigues parameterization of microstructural solution spaces, interested readers are referred to Refs. [11; 26].

The ODF, \(A(\mathbf{r},t)\), could be used to compute homogenized elastic stiffness \(<C>\)
through its volume integration over the fundamental region, \(\Omega\):

\[<C>=\int_{\Omega}C(\mathbf{r})A(\mathbf{r},t)\;dV \tag{1}\]

where \(C(\mathbf{r})\) includes the single-crystal material properties required to compute the homogenized elastic stiffness values given by \(<C>\). By controlling the ODF values, desired homogenized (meso-scale) properties can be obtained as described in Eq. (1). However, the ODF (\(A(\mathbf{r},t)\geq 0\)) must satisfy the volume normalization constraint which is expressed as follows:

\[\int_{\Omega}\;A(\mathbf{r},t)\;dV=1 \tag{2}\]

The homogenized (meso-scale, volume-averaged) properties of the microstructures are obtained using the given expression in Eq. (1). Here, the integration for the homogenized properties is performed over the fundamental region by considering the lattice rotation, \(\mathbf{R}\). Given the Rodrigues orientation parameter, \(\mathbf{r}\), the rotation, \(\mathbf{R}\), can be obtained with the following expression:

\[\mathbf{R}=\frac{1}{1+\mathbf{r}\cdot\mathbf{r}}(I(1-\mathbf{r}\cdot\mathbf{r })+2(\mathbf{r}\otimes\mathbf{r}+I\times\mathbf{r})) \tag{3}\]

The finite element discretization of the microstructural orientation space is exhibited in Fig. 2. Here, each independent nodal point of the finite element mesh represents a unique ODF value for the associated crystallographic orientation. For \(N\) independent nodes with \(N_{elem}\) in the finite element discretization (\(N_{int}\) integration points per element), Eq. (1) can be approximated as follows at a given time,
\[\begin{split}<C>&=\int_{\Omega}C(\mathbf{r})A(\mathbf{r},t )\,dV\\ &=\sum_{n=1}^{N_{elem}}\sum_{m=1}^{N_{int}}C(\mathbf{r_{m}})A(r_{m}) \omega_{m}|J_{n}|\frac{1}{(1+r_{m}.r_{m})^{2}}\end{split} \tag{4}\]

where \(A(r_{m})\) is the ODF value at the \(m^{th}\) integration point with global coordinate \(r_{m}\) (orientation vector) of the \(n^{th}\) element. \(|J_{n}|\) is the Jacobian determinant of the \(n^{th}\) element and \(\omega_{m}\) is the integration weight of the \(m^{th}\) integration point.

Considering the crystallographic symmetry, the homogenized property of Eq. (4) can alternatively be approximated in the linear form as \(<C>\)\(=\)\(P^{T}\mathbf{A}\), where \(P\) is the property matrix that is a product of the single-crystal material properties

Figure 2: Finite element discretization of the orientation space for face-centered cubic (FCC) microstructures. The red-colored nodal points show the independent ODF values while the blue-colored nodes indicate the dependent ODFs as a result of the crystallographic symmetries.
and finite element discretization of the orientation space, and \(\mathbf{A}\) is the column vector of the ODF values for the independent nodes (this work uses 76 independent nodes to model the FCC microstructure) of the finite element mesh (see Fig. 2). Moreover, the ODF should satisfy the normalization constraint (see Eq. 2). The normalization constraint is mathematically equivalent to the fact that the sum of the probabilities for having all possible crystallographic orientations in a microstructure must be one.

During the deformation process, ODFs change due to the reorientation of the grains. It evolves from the initial ODFs (t=0) to the final deformed ODFs (t=t). The evolution of the ODF values is governed by the ODF conservation equation, which satisfies the volume normalization constraint of Eq. 2. Equation (5) shows the Eulerian rate form of the conservation equation in the crystallographic orientation space with a gradient operator given by [25]:

\[\frac{\partial A(\mathbf{r},t)}{\partial t}+\nabla A(\mathbf{r},t)\cdot v( \mathbf{r},t)+A(\mathbf{r},t)\nabla\cdot v(\mathbf{r},t)=0 \tag{5}\]

where \(v(\mathbf{r},t)\) is the reorientation velocity.

The texture evolution can be calculated by the microstructure constitutive model in terms of a velocity gradient (\(\mathbf{L}\)) definition (see Eq. 6 below), which is linked to \(v(\mathbf{r},t)\) by the Taylor macro-micro linking hypothesis. The Taylor hypothesis assumes that the crystal velocity gradient is equal to the macro velocity gradient [25]. To compute the reorientation velocity, a rate-independent constitutive model is adopted. Resulting \(A(\mathbf{r},t)\) (current texture) which evolves from \(A(\mathbf{r},0)\) (initial texture) is solved by utilizing the constitutive model and finite element representation of the Rodrigues orientation space.

Each deformation process (e.g., tension/compression and shear) yields a par
ticular ODF as output after applying the load for a certain time. The macro velocity gradient, \(\mathbf{L}\), for a particular process is provided as input in the crystal plasticity solver to investigate the ODF evolution during that process. While designing a process sequence to obtain the desired texture, the macro velocity gradient describing the processing route (type and sequence of the process and strain rate) is solved as a design variable. The velocity gradient of a crystal with the orientation, \(\mathbf{r}\), can be expressed as

\[\mathbf{L}=S+\mathbf{R}\sum_{\alpha}\dot{\gamma^{\alpha}}\mathbf{\bar{T}^{ \alpha}}\mathbf{R}^{T}, \tag{6}\]

where \(S\) represents the lattice spin, \(\mathbf{R}\) indicates the lattice rotation, \(\dot{\gamma^{\alpha}}\) and \(\mathbf{\bar{T}^{\alpha}}\) indicate the shearing rate and Schmid tensor for the slip system \(\alpha\), respectively. Here, the crystal velocity gradient is defined in terms of the time rate change of the deformation gradient. Equation (6) demonstrates that the velocity gradient is decomposed into two components. The first component, lattice spin \(S\), is related to the elastic deformation gradient by assuming that the deformation gradient (\(\mathbf{F}\)) is decomposed into elastic (\(\mathbf{F}^{e}\)) and plastic (\(\mathbf{F}^{p}\)) parts such that: \(\mathbf{F}=\mathbf{F}^{e}\mathbf{F}^{p}\). In particular, the lattice spin is equal to \(\mathbf{\dot{R}^{e}}\mathbf{R^{e}}\), where \(\mathbf{R^{e}}\) is evaluated through the polar decomposition of the elastic deformation gradient: \(\mathbf{F}^{e}=\mathbf{R^{e}}\mathbf{U}^{e}\), where \(\mathbf{U}^{e}\) is the unitary matrix of the polar decomposition [35]. The second component in Eq. (6) shows the rotated plastic velocity gradient (\(\mathbf{R}\mathbf{L}^{p}\mathbf{R}^{T}\)). Here, \(\mathbf{L}^{p}\) denotes the plastic velocity gradient, which is related to the combined shearing of the slip systems given by \(\sum_{\alpha}\dot{\gamma^{\alpha}}\mathbf{\bar{T}^{\alpha}}\).

The macro velocity gradient expression of Eq. (6) for different deformation processes can also be written in the following matrix form, given in Eq. (7). The detailed derivation from Eq. (6) to Eq. (7) is skipped here for brevity and it can be
found in Ref. [25].

\[\mathbf{L}=\alpha_{1}\begin{bmatrix}1&0&0\\ 0&-0.5&0\\ 0&0&-0.5\end{bmatrix}+\alpha_{2}\begin{bmatrix}0&0&0\\ 0&1&0\\ 0&0&-1\end{bmatrix}+\alpha_{3}\begin{bmatrix}0&1&0\\ 1&0&0\\ 0&0&0\end{bmatrix}\\ +\alpha_{4}\begin{bmatrix}0&0&1\\ 0&0&0\\ 1&0&0\end{bmatrix}+\alpha_{5}\begin{bmatrix}0&0&0\\ 0&0&1\\ 0&1&0\end{bmatrix} \tag{7}\]

Each matrix in Eq. (7) defines a deformation process, e.g. tension/compression (\(\alpha_{1}\)), plane strain compression (\(\alpha_{2}\)), and shear modes (\(\alpha_{3},\alpha_{4},\alpha_{5}\)).

We define a processing path of \(n\) steps with a fixed deformation duration time \(\Delta t\) in an autoregressive pattern below. \(\phi_{F_{i}}(A)\) stands for a deformation process with load \(F_{i}\) (e.g., tension/compression and shear modes defined through the velocity gradient, \(\mathbf{L}\), in the finite element simulations) for a duration \(\Delta t\) on a polycrystalline microstructure described by the ODF, \(A(r,i\Delta t)\), at time \(i\Delta t\):

\[\begin{split} A(\mathbf{r},\Delta t)&=\phi_{F_{0}}(A(\mathbf{ r},0))\\ A(\mathbf{r},2\Delta t)&=\phi_{F_{1}}(A(\mathbf{r},\Delta t))\\ &\vdots\\ A(\mathbf{r},n\Delta t)&=\phi_{F_{n-1}}(A(\mathbf{r},(n-1)\Delta t)) \end{split} \tag{8}\]

Here \(F_{i}\in\mathcal{F}\) is picked from a pre-specified set \(\mathcal{F}\coloneqq\{f_{1},\ldots,f_{k}\}\) (e.g. a set of shear force \(\{xy,xz,yz\}\)) of deformations we can apply in each step. An optimal processing path \(P^{*}\) for an initial ODF \(A(\mathbf{r},0)\) is a sequence of deformation processes that could achieve the best ODF set \(A^{*}(\mathbf{r},n\Delta t)\) at time \(n\Delta t\). \(A(\mathbf{r},n\Delta t)\) is evaluated by the desired orientation-dependent properties which is an integral
of \(A\) and single-crystal material property matrix (weight function) \(C(\mathbf{r})\) over the fundamental region \(\Omega\) (as shown in Eq. (1)).

\[\begin{split}\max_{A}\,F=\sum_{i}^{6}\omega_{ii}<C_{ii}>+\,\sum_{i <j}^{6}\omega_{ij}<C_{ij}>\\ \text{subject to}\int_{\Omega}A(\mathbf{r},i\Delta t)dV=1\\ A(\mathbf{r},i\Delta t)\geq 0\end{split} \tag{9}\]

where the composite objective function (\(F\)) includes the summation of the constants of the \(6\times 6\) anisotropic homogenized elastic stiffness matrix (\(<C_{ii}>\), \(<C_{ij}>\)) that are computed using the formulation given in Eq. (1). In this problem, the diagonal entries of the homogenized stiffness matrix (\(<C_{ii}>\)) are assumed to be more important, and thus multiplied by a weight factor of \(\omega_{ii}=1\), while the cross-diagonal terms (\(<C_{ij}>\)) are multiplied by a weight factor of \(\omega_{ij}=\frac{1}{2}\). Note that the symmetric terms (\(<C_{ij}>\) and \(<C_{ji}>\)) are only counted once (through the \(i<j\) condition in the second term of the objective function) in the summation formula. The underlying idea behind this composite objective function is to improve the homogenized elastic stiffness of the microstructure in different directions rather than improving the stiffness along a particular direction.

To address this task, we need 1) an efficient deformation simulator \(\phi_{F}(A)\) and 2) a path-searching algorithm. For the algorithm, various works have been done on such local search problems. These algorithms work by making a sequence of decisions locally to optimize the objective function, with well-known methods including simulated annealing. In this work, we focus on a novel approach to building up \(\phi_{F}(A)\) by a neural network. As a function approximator, a neural network could learn the behavior of different deformations accurately from data and predict fastly to accelerate process design. It is chosen to address the challenge
arising from the numerically intractable nature of the processing design problem and the speed bottleneck of a traditional simulator.

## 3 End-to-end deformation prediction with Neural Networks

We aim to develop an efficient neural network to replace the FE predictor in computing deformation results. In this data-driven method, we no longer focus on solving the conservative equation presented in Section 2 but exploit some physical constraints to build up a surrogate neural network. The neural network aims to approximate the deformation process \(\phi_{F}\):

\[A(\mathbf{r},\Delta t)=\phi_{F}(A(\mathbf{r},0)),\ \ \mathbf{r}\in\Omega \tag{10}\]

where \(A(\mathbf{r},0)\) denotes the ODFs before process and \(A(\mathbf{r},\Delta t)\) is the ODFs after applying a deformation \(F\) of duration \(\Delta t\).

We define \(NN_{F}(A;\theta)\) to be the surrogate neural network model with \(\theta\) parameters. In this work, we employ multilayer perceptron (MLP) to approximate \(\phi_{F}(A)\). For an MLP with \(L\) hidden layers, we have:

\[\begin{split} a^{(0)}&=x\\ z^{(i+1)}&=M^{(i)}a^{(i)}+b^{(i)}\\ a^{(i+1)}&=\varphi^{(i+1)}(a^{(i)})=\xi^{(i+1)}(z^ {(i+1)})\\ NN_{F}&=\varphi^{(L)}\circ\cdots\circ\varphi^{(1)} \end{split} \tag{11}\]

where \(x\) is the model input, \(\xi^{i+1}\) denotes the activation function of layer \(i\). Considering that the ODF stands for the probability density of the orientation space, the ODF non-negativity and volume normalization constraints must be satisfied. A ReLU followed by a normalization layer which divides each channel with a
material-specific weighted sum between volume fraction constant and output from the previous layer are applied to keep the network outputs physically feasible.

The model parameters can be learned by minimizing the error between exact ODFs and network predictions. One common choice is using squared error to match them. Here we introduce the specific objective value weights to assign different importance to different elastic stiffness constants. Considering that the integral of the homogenized properties in Eq. (1) over a discrete fundamental is well approximated the weighted sum (Eq. (4)), the weighted mean squared error (WMSE) should be a more appropriate choice:

\[WMSE=\frac{1}{B}\frac{\sum_{i=1}^{n}w_{i}\|y_{i}-y_{i}^{*}\|_{2}^{2}}{\sum_{i= 1}^{n}w_{i}} \tag{12}\]

where y and y* stand for true and predicted ODF values, respectively, \(w_{i}\) are the orientation weights to calculate the desired objective material property, and \(B\) is batch size.

In this section, we evaluate our method on a stiffness optimization task. We find a minimal difference in predicted ODFs and objective property value, compared with the traditional FE method. In addition, the neural network predicts hundreds of times faster than traditional methods which enables large-scale searching.

### Model Setup

We consider a stiffness optimization task of finding out the best texture evolution path by applying 31 different possible deformation processing modes (tension, compression, and shear along xy, xz, and yz axes, respectively, and their combinations, except the case when none of them are selected). The fundamental region is discrete with 76 possible outcomes (the FCC microstructure is modeled
with 76 independent ODF values). Only one of the deformation processes can be applied in one step. Each step is 0.1 seconds in duration. An evolution path consists of 10 sub-processes executed in order. The elastic stiffness is selected as the objective mechanical property to be maximized while the material of interest is Copper (Cu). The following single crystal properties are taken for Copper (Cu): \(C_{11}=C_{22}=C_{33}=168\) GPa, \(C_{12}=C_{21}=C_{13}=C_{23}=C_{31}=C_{32}=121.4\) GPa, and \(C_{44}=C_{55}=C_{66}=75.4\) GPa [25]. The objective function is defined as the sum of the homogenized elastic stiffness constants (\(C\)) using a composite function definition with weights of 1 for diagonal entries (\(C_{ii}\)) and weights of 0.5 for off-diagonal entries (\(C_{ij}\), \(i\neq j\)) of the elastic stiffness matrix. The symmetric off-diagonal (\(C_{ij}\), \(i\neq j\)) entries are only considered once. The synthetic dataset of size 5000 is uniformly initialized and normalized to keep the ODFs feasible according to the unit volume normalization constraint (Eq. (2)). The deformation mode results of each initial ODF are generated by the FE solver. In the following experiment, data points used for training and testing are selected randomly. The dataset is also randomly divided by an 80%/20% train-test ratio.

### Searching Algorithm and network structure

Due to the complex nature of the combinational optimization problem, calculating all the \(31^{10}\) results is impractical even with fast NN. Instead, we employ a heuristic algorithm that helps to find out the optimal process path by sampling from a multinominal distribution constructed from the predicted value at each stage.
The neural network has 1 hidden layer of 760 neurons with hyperbolic tangent activation functions except the last ReLU. ReLU and normalization are applied to fit the physical constraint in the last layer. The network is trained with a mini-batch of size 128. ADAM with warm restarts [36] is used for network training. When searching, restart 1000 times with weight base \(\beta=5\) for each initial texture.

### Network performance and searching samples

Figure 3 depicts the comparison between the textures obtained from the neural network prediction and FE simulator at different time steps of processing. The loading conditions are as follows: for Fig. 3 (a) and (d) shear force in the YZ
plane, for Fig. 3 (b) and (e) a combined force of tension, compression, and shear in XZ and YZ planes, and for Fig. 3 (c) and (f) a combined force of tension, compression, and shear in the YZ plane. The performance of the neural network is also evident from Fig. 4 and Table 1, which presents the relative average errors of independent ODFs over the fundamental region for different deformation modes.

Figure 3: The comparison between the predicted textures by the neural network model (a, b and c) and finite element crystal plasticity model (d, e and f) at time steps t=0.3 sec (a and d), t=0.8 sec (b and e), and t=1 sec (c and f).
The computation time for predicting all 31 deformation modes' results from an arbitrary set of initial ODF values in the neural network is 0.2213s while the traditional finite element simulator takes 152s on average when running in a 31 process-based parallelism. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 and 100000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 1000 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 100 training data. The neural network gives predictions 686.85 times faster. The predictions of the neural network were generated on NVIDIA(R) with 10 training data.
Tesla P100-16G with single-core Intel(R) Xeon(R)@2.3Ghz and simulator run on AMD(R) EPYC C2D-highcpu-32@2.45 Ghz. The average difference in the elastic stiffness prediction values of the physics-based model and neural network after 10 processes is 0.4663 GPa.

The neural network-based surrogate model replaced the physics-based finite

Figure 5: Texture evolution prediction through the optimum processing path by the neural network surrogate model. The figure shows the different steps of deformation processing from an initial texture to a final optimum texture which maximizes an objective function defined for the homogenized elastic stiffness constants of a Copper microstructure. T and C stand for tension and compression, respectively, and XY, XZ & YZ represent the corresponding shear processes.
element model to find the optimum processing paths. It is already mentioned that an optimization problem is developed with the goal of maximizing the sum of the elastic stiffness constants of Copper. The solution of the optimization problem provides a maximum objective function value with the corresponding texture. It also suggests an optimum processing route to obtain the optimum texture. We

Figure 6: Texture evolution through the optimum processing path by the physics-based simulator. The figure shows the different steps of deformation processing from an initial texture to a final optimum texture which maximizes an objective function defined for the homogenized elastic stiffness constants of a Copper microstructure. T and C stand for tension and compression, respectively, and XY, XZ & YZ represent the corresponding shear processes.
determine the processing route for five different initial textures. Though the optimum solutions vary with the initial textures, almost equal objective function values are obtained for all cases. The average of the sum of the maximum stiffness constants of Copper is found to be 885.2 GPa from the neural network surrogate model whereas this value is 884.9 GPa for the FE-based process model.

The neural network predicted the optimum processing path with texture evolution after each step to obtain the optimum texture at the end starting from a certain initial texture, as shown in Fig. 5. Similarly, the optimum process route and texture evolution obtained from the physics-based model for another initial texture is displayed in Fig. 6. Here, T and C stand for tension and compression, respectively, and XY, XZ & YZ represent the corresponding shear processes. The strain rate is constant for all steps (i.e., 1 \(s^{-1}\)). The processing is allowed to have

Figure 7: A single crystal optimum texture is obtained using linear programming to maximize the homogenized elastic stiffness constants without considering processing.
only one loading condition or combined loading conditions as reported in Fig. 5 and Fig. 6. Both models are found to provide polycrystalline textures as optimum solutions.

If we solve the same optimization problem using linear programming without considering the processing (only structure-property problem), a single crystal texture is found as an optimum ODF solution (see Fig. 7). In this case, the maximum objective function value is 896.2 GPa, which indicates the theoretical maximum value of the objective function. On the other hand, the neural network surrogate model, which accounts for processing (process-structure-property problem), provides an optimum texture with 885.2 GPa for the objective function value (which is a higher value compared to the randomly oriented texture providing an objective function value of 878.4997 GPa). However, the theoretical solution of the structure-property problem results in a single crystal solution while the solutions are polycrystalline textures for the process-structure-property problems using both the neural network and physics-based model. Even though single crystal textures provide the theoretically possible maximum value for elastic stiffness constants, their manufacturing is difficult. On the other hand, the polycrystalline textures of the process-structure-property problem can easily be manufactured with the presented simple deformation processing modes. Therefore, considering the effects of the processing is significant for manufacturing and bridging materials design and manufacturing [37].

This work can be a valuable addition to the literature to investigate the process-structure-property linkages of polycrystalline microstructures. The neural network surrogate model can be a substitute for the physics-based simulator as it is computationally faster and accurate. Moreover, the processing-produced opti
mum textures predicted by the surrogate model can still improve the mechanical properties (e.g. elastic stiffness) near to the maximum theoretical values.

## 4 Conclusions

In this work, we developed a surrogate neural network model to accelerate the process design of polycrystalline microstructures. The traditional physics-based FE simulator was found inefficient to address the large-scale searching task. With the representation potential of the neural network, a pre-trained network could predict the ODF after deformation significantly fast and accurately. An example design problem was solved to find the optimum processing route maximizing an objective function defined for the homogenized elastic stiffness constants of Copper. The results demonstrate a good match between the predictions of the physics-based simulator and neural network surrogate model. Studies on accumulative error analysis of neural network are reserved in the future. Future work may also involve the integration of a similar data-driven modeling strategy for the concurrent multi-scale modeling of metallic components.

## 5 Acknowledgements

Material in this paper is based upon work supported in part by the Air Force Office of Scientific Research under award number FA9550-20-1-0397. MH and PA also acknowledge the support from the Air Force Office of Scientific Research Young Investigator Program under grant FA9550-21-1-0120 and from the National Science Foundation under award number 2053840.
## 6 Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. The data will be made available upon request.
