In designing a microstructural materials system, there are several key questions associated with design representation, design evaluation, and design synthesis: how to quantitatively represent the design space of a heterogeneous microstructure system using a small set of design variables, how to efficiently reconstruct statistically equivalent microstructures for design evaluation, and how to quickly search for the optimal microstructure design to achieve the desired material properties. This paper proposes a new descriptor-based methodology for designing microstructural materials systems. It is proposed to use a small set of microstructure descriptors to represent material morphology features quantitatively. The descriptor set should be able to cover microstructure features at different levels, including composition, dispersion status, and phase geometry. A descriptor-based multiphase microstructure reconstruction algorithm is developed accordingly that allows efficient stochastic reconstructions of microstructures in both 2D and 3D spaces for finite element analysis (FEA) of material behavior. Finally, the descriptor-based representation allows the use of parametric optimization approach to search the optimal microstructure design that meets the target material properties. To improve the search efficiency, this paper integrates state-of-the-art computational design methods such as design of experiment (DOE), metamodeling, statistical sensitivity analysis, and multi-objective optimization, into one design optimization framework to automate the microstructure design process. The proposed methodology is demonstrated using the design of a polymer nanocomposites system. The choice of descriptors for polymer nanocomposites is verified by establishing a mapping between the finite set of descriptors and the infinite dimensional correlation function.

## Introduction

Design of microstructural materials systems is an emerging topic in material design for accelerating advanced material development. As brought up by Olson and McDowell, superior performance of material systems can be achieved through the “microstructural mediated design of materials” [1,2] by changing the material compositions, phases, and microstructure morphologies. However, great challenges exist in (1) how to model the complex interactions among multiple material constituents to predict the properties of interests; and (2) how to automate the material design optimization at the microstructure level.

In the existing literature of materials design, combinatorial search has been used to choose different material compositions from materials databases [3,4], that store information of individual compositions' properties and their interactions. However, this composition-based material analysis and design approach only involves choosing constituents and determining their percentages, but no longer suffices in designing complex microstructural materials system. The morphology of microstructure (i.e., the spatial arrangements of local microstructural features) has a large impact on the overall properties and performances of a materials system. Taking polymer nanocomposites as an example, the electrical conductivity of such material is highly associated with microstructure percolation (phase dispersion); the morphology of interphases between nanoparticles (fillers) and polymer matrix is crucial for predicting damping properties [5]. Furthermore, heterogeneity in microstructure is the root cause of material randomness at multiple length scales. This paper considers “microstructure” as the material structure features, including composition, dispersion status, and geometry of constituents/phases.

Quantitative representation of microstructures plays a key role as “the taxonomy basis for describing microstructural materials” under the paradigm of integrated computational design of Materials [2,6]. Following the “process-structure-property” relation [7,8] that defines the whole path line in material design, this paper is focused on utilizing the microstructure-property relation to optimize microstructures for desired properties. To automate materials design, there are three major research questions involved: (1) **design representation**: how to describe a heterogeneous microstructure quantitatively using a small set of descriptors; (2) **design analysis**: how to efficiently evaluate material property based on microstructure information; (3) **design synthesis**: how to efficiently explore the vast, high dimensional design space to search for optimal material designs. This paper focuses on the quantitative representation of microstructure, as well as establishment of a computational framework to search the optimal microstructures based on structure-property simulations.

Microstructure characterization and reconstruction are important for both microstructure design representation and design evaluation. Design representation requires a finite set of quantitative descriptors of the infinite dimensional microstructure image as design variables, which can be obtained via microstructure characterization. In design evaluation, statistical FEA analysis of material properties [9,10] requires microstructure reconstruction to generate multiple random but statistically equivalent microstructures. Microstructure characterization describes the heterogeneous microstructure quantitatively with a small set of microstructure parameters based on image analysis. For given microstructure descriptions, microstructure reconstruction generates statistically equivalent digital microstructures for stochastic analysis of material behavior [11]. Three basic categories of microstructure characterization/reconstruction approaches are widely adopted in material analysis and design: (1) the physical descriptor-based approach, (2) the correlation function-based approach, and (3) the random field-based approach. In the physical descriptor-based approach, important structural parameters (descriptors) highly related to the interested material properties are chosen as microstructure representation. Typical microstructure descriptors fall into three categories: composition, dispersion status, and geometry of the inclusions. **Composition descriptors** indicate the percentage of each material constituent, such as volume fraction of filler. **Dispersion descriptors** depict the inclusions' spatial relation and their neighbor status, such as the nearest neighbor distance, the ranked neighbor distances, pair correlation [12–16], and spatial pattern descriptors obtained from the infinite dimensional *F*-, *G*-, and *K*-functions [17]. **Geometry descriptors** provide information of the shapes of inclusions such as the size/radius distribution, surface area, surface-to-area/volume fraction, roundness, eccentricity, elongation, rectangularity, tortuosity, aspect ratio, etc. [12,16,18–23]. The major strength of physical descriptors is the clear physical meanings they offer and the meaningful mappings to processing parameters [24,25].

*N*-point correlation function [18] is another category of widely used characterization/reconstruction techniques. Torquato et al. have shown that the microstructure of heterogeneous materials can be characterized statistically via various types of *N*-point correlation functions [13,18–20]. With a good balance between computational cost and accuracy, 2-point correlation functions (*N* = 2) and 2-point cluster correlation functions [11] are widely adopted in practice. Despite its high accuracy, the correlation function approach often involves a high computational cost in reconstruction (especially for 3D reconstruction) due to the pixel-moving optimization strategy and the cost in evaluating 2-point cluster correlations. The correlation function is infinite dimensional, causing difficulties for design representation, and the function parameters do not have any physical meaning.

Random field-based approaches are emerging techniques which model microstructures as random fields. Compromised in accuracy, it requires much less computational power in reconstruction by using a level-cut strategy based on a Gaussian random field [18,26,27] or a Poisson random field [28]. The random field approach has a similar disadvantage of being infinite dimensional as the correlation function approach.

Despite the large amount of researches done on microstructure characterization and reconstruction there still exists the need of developing an efficient microstructure reconstruction algorithm without lowering the accuracy. Our research goal includes bridging the information gap between the descriptor-based approach and the correlation function-based approach to facilitate microstructure design.

For design synthesis, combinatorial search has been used to choose material compositions from databases [3,4] without considering the complex microstructure morphology and the interactions among constituents. Data mining has been applied to analyzing and designing materials based on archived microstructure-property information [29,30], but the predictive capability of the approach is limited. Optimal microstructure design has also been explored using topology optimization [31], which optimizes the distribution of different phases of materials. However, this type of works focuses primarily on the periodic cell structure, neglecting the inherent heterogeneity and randomness in structure and property.

This paper proposes a new descriptor-based methodology for designing microstructural materials systems. We develop a descriptor-based characterization method to represent material morphology using a small set of microstructure descriptors that enable fast microstructure reconstruction (Sec. 2). The descriptor-based representation also allows the use of parametric optimization approach to search the optimal microstructure design for achieving the desired material properties. Integration of state-of-the-art computational design methods such as DOE, dimension reduction, metamodeling, and multi-objective optimization improves the search efficiency (Sec. 3). We demonstrate the proposed methodology using the design of polymer nanocomposites (Sec. 4).

## Descriptor-Based Microstructure Characterization and Reconstruction

In Sec. 2.1, we present the proposed multiphase microstructure characterization and reconstruction algorithm using the example of nanoparticle polymer composite. The heterogeneous microstructure is represented quantitatively by a finite set of descriptors covering composition, dispersion status, and shape information of the inclusions in matrix. A hybrid optimization and random assignment approach for microstructure reconstruction is then presented and demonstrated on multiple microstructure samples with distinctive features. In Sec. 2.2, the completeness of the descriptor set is verified by establishing a mapping between the finite chosen descriptors and the infinite dimensional correlation function.

### Descriptor-Based Microstructure Characterization and Multiphase Reconstruction Algorithm.

The main issue of microstructure representation is to determine a small set of descriptors that are equivalent to the high dimensional representation of morphology. Preprocessing of the microstructure image (e.g., from techniques like scanning electron microscope (SEM) or transmission electron microscope [11]) distinguishes different phases in a digital image. The phase information of each pixel is marked by different greyscale values. Figure 1 illustrates the transformation of the bi-phase microstructure (left) to a binary image (right), where black pixels stand for nanoparticle filler and white pixels represent polymer matrix.

As introduced in Sec. 1, three levels of microstructure features: composition, dispersion, and geometry, are important to fully describe microstructure morphology. At the highest coarse level, *composition* is the lowest order of microstructure information as it only captures the homogenized response of an entire material system. One level down, *dispersion* is capable to capture material properties induced by local features. For example, the nearest distance between inclusions and voids in alloy determines the fracture behaviors, and the orientation of the fillers has an impact on polycrystal plasticity behavior. At the lowest level, the impact of *geometry* should not be neglected [32] in cases such as determining the yield point or the critical load in microvoid-induced microbuckling. Descriptors for each level can be either collected based on experience or identified using statistical learning of the structure-property relation. It should be noted that descriptors used at each level are statistical, quantified by mean and higher order moments such as variance, skewness etc.

Figure 2 shows the process of descriptor-based multiphase microstructure reconstruction to match the prespecified descriptors. The reconstructed microstructures are stochastic but following the same statistical microstructure characteristics. Randomness is inherent in the heterogeneous material microstructures. The goal of statistical reconstruction is to generate statistically equivalent random digital microstructure images to simulate the heterogeneity and randomness of material behavior [33]. Detailed steps are provided as follows:

- (1)
Prespecification of target microstructure characteristics (composition, dispersion, and geometry descriptors). The target microstructure characteristics are obtained from microstructure characterization or determined by new material designs.

- (2)
Dispersion reconstruction: this step reconstructs particles' dispersion status. The location of each particle's center is determined by matching dispersion descriptors using the simulated annealing (SA) algorithm. The amount of particles is determined by both composition and key geometry descriptors such as the mean particle radius.

- (3)
Geometry reconstruction: at the locations of particle centers, the algorithm reconstructs each particle based on the geometry descriptors. The geometry descriptors' values are generated randomly following predefined distributions (e.g., exponential distribution of nanoparticle's diameter). Steps (1)–(3) are applied to each material phase (except matrix) in multiphase microstructure reconstruction. If a pixel reconstruction is required, the pixelated particle profile will be generated according to the geometry descriptors.

- (4)
Adjustment: this is the final step of reconstruction that fine tunes the digital image to match the composition descriptors. The first part of adjustment is an optional step which eliminates the overlap between particles. For example, for “dispersed hard disk” structure (e.g., long fiber composites' cross-section), the algorithm will switch the positions of fillers randomly until the overlap area is minimized (SA). This adjustment is also needed to eliminate the overlap between different phases in multiphase reconstruction. For some other materials in which the particles are overlapping or connecting with each other to form the porous network, the algorithm skips the overlap-elimination step and directly starts the

*composition compensation*step. Discrepancy between the actual composition (volume fraction, (VF)) of the image generated from step (3) and the target composition exists due to the overlap between particles and the incomplete particles on the image boundary. For pixelated images, pixels are added or subtracted on particle boundaries until the composition is matched exactly. For smooth-profile images, the radius values of randomly chosen particles are perturbed to match the VF.

When applying the proposed methodology to nanoparticle reinforced polymer composites in this study, we choose four descriptors based on the existing literature [12,13]:

- (1)
Composition descriptor: volume fraction ($VF\xaf$). It is the volume percentage of the inclusion (particle) phase. Though it is actually “area fraction” in 2D, we stick to the widely-accepted term “volume fraction”. In this paper, having bar on the top of microstructure descriptor symbol indicates the variable is deterministic; otherwise it is a random variable.

- (2)
Dispersion descriptor: average nearest cluster center distances ($r\xafd$), which is the center distance between each particle and its nearest neighbor particle.

- (3)Geometry descriptor: distribution of equivalent radius (
*r*). We define the equivalent radius of a single particle as_{c}$rc=rl\xd7rs$(1)in that

*r*is radius along the long axis of the particle, and_{l}*r*is the radius along the short axis of the particle. The average of_{s}*r*is equivalent to particle number $N\xaf$ when $VF\xaf$ is known. $N\xaf=(l\xaf2\xb7VF\xaf/\pi r\xafc2)$, where $l\xaf$ is the side length of the square reconstruction window._{c} - (4)
Geometry descriptor: distribution of elongation ratio (

*e*). Also known as aspect ratio, elongation ratio equals to the ratio of the major diameter to the minor diameter. When we simplify a particle's geometry into an ellipsoid, elongation ratio, and average radius together can give a complete description of a particle's shape._{l}

The final microstructure image can be either pixelated image or smooth-profile image in order to adapt to different mesh techniques. The pixelated reconstruction output is in the format of a 0-1 matrix, that is suitable for mapped mesh (one pixel is one element). Detailed information of FE modeling using mapped mesh will be provided in Sec. 4.2 and Fig. 11. The smooth-profile reconstruction output is represented by a matrix that lists each particle's center location, major radius, minor radius, and the orientation. Figure 3 provides the visualization of a smooth-profile reconstruction in abaqus and the free mesh on it. Kanit et al. [34] reported the influence of different mesh techniques on simulated properties and claim a small difference between mapped mesh and free mesh. This paper chooses the mapped mesh, and leaves the free mesh for future work.

One reconstruction case is shown in Fig. 4 to demonstrate the effectiveness of the proposed algorithm. Figure 4(a) is the target microstructure image with the size of 300 × 300 pixel. Figures 4(b)–4(d) are three realizations using the descriptor-based reconstruction algorithm that takes 110 s, respectively, on Intel Core i5 CPU, with a single 2.40 GHz processor. This is a significant saving compared to ∼10 h when using the correlation function-based reconstruction [11]. The latter approach is computationally infeasible for iterative search in the microstructure optimization process.

This paper evaluates the accuracy of reconstructions quantitatively via the matching between descriptor-based representations with the correlation functions, which are higher-order representations of microstructures. We compare both 2-point correlation and 2-point cluster correlation functions (defined in Sec. 2.2) between the target microstructure and the reconstructed microstructures. In Fig. 5, we can see that the target and reconstructed correlation functions match very well, that indicates that the descriptor-based characterization and reconstruction are both accurate. We can further refine the descriptor-based reconstruction by applying the correlation function-based matching algorithm [11] if higher accuracy is desired.

Both 2D and 3D microstructures can be generated using this algorithm, depending on the type of inputs. If 3D microstructure descriptors are available, 3D microstructures can be reconstructed following the same procedure. One demonstration is shown in Fig. 6. The 3D reconstructed structure matches the listed 3D descriptor values. It takes less than 0.5 h to generate this 300^{3} voxel size image on Intel^{®} Xeon^{®} CPU 3.20 GHz using the descriptor-based reconstruction. In contrast, it takes more than 100 h for reconstructing a much smaller structure (90^{3} voxels) using the correlation function-based algorithm. Reconstruction of a cubic structure of 300^{3} voxel size using correlation function-based voxel moving techniques is infeasible with the conventional computing resource due to the extremely large amount of memory required for restoring the voxel-voxel distance information and the voxel-by-voxel moving strategy in optimization search. In the remaining paper, we illustrate our design methodologies using 2D descriptor-based reconstructions which also reduce the cost of FEA simulations.

### Linking the Descriptor-Based Representation and 2-Point Correlation Function

where *x*_{1} and *x*_{2} are two arbitrary points in the system. We can interpret 2-point correlation as a normalized probability of finding two points at locations *x*_{1} and *x*_{2} that are both in the phase of inclusion. Another form of *S*_{2}(*x*_{1}, *x*_{2}) is *S*(*r*), a function of the distance *r* between any two points randomly chosen from the microstructure. In this case, 2-point correlation function is interpreted as the normalized probability that both of the two ends of a randomly tossed line section with length *r* are in the same phase.

The infinite dimension of 2-point correlation function limits its usage in material analysis and design. To represent the 2-point correlation function with a small set of parameters, several analytical forms were developed in previous research to fit the 2-point correlation function (Table 1). A common feature of fitting functions (I)–(IV) is the use of an exponential function to capture the decreasing trend of the 2-point correlation function (see Fig. 5 for a 2-point correlation function from a real microstructure image). These four functions all consider 2-point correlation functions as smooth and monotone decreasing functions, that is an inaccurate representation of the fluctuations observed in the long tail of the true correlation functions. Information contained in the “tail with fluctuation” is not trivial [40], because it captures the dispersion features of the microstructure, such as the nearest distances between inclusion clusters. This tail fluctuation can also be demonstrated using a toy example shown in Fig. 7, where the black phase is distributed in square grids with interval of 20 pixels. X coordinate of the first peak on 2-point correlation function's tail is correspondent to the nearest center distance between clusters, that is 20 pixels. The fluctuation in the 2-point correlation function of this example is very prominent due to the periodicity of the inclusion's distribution pattern. For real heterogeneous microstructure with great randomness, smaller oscillation amplitudes have been often observed (see Fig. 5(a)). The only fitting function that captures the fluctuation feature is function (V) listed in Table 1. Function (V) introduces a fluctuation term by using a sine function. However, a common issue of fitting functions (I)–(V) is that they all require some artificial fitting parameters with no or vague physical meanings. Function (VI) does not use any fitting parameters, but it is only applicable to microstructures dispersed with equal size disk-shape inclusions. Again, it is incapable of capturing the tail fluctuation either.

All parameters in the above equation are the microstructure descriptors used in characterization (Sec. 2.1). $VF\xaf$ is the volume fraction; $e\xafl$ is the elongation ratio of the particles; $r\xafc$ is the average radius of particles; and $r\xafd$ is the average nearest center distance.

There are two terms in this analytical model, the exponential term and the fluctuation term. Similar to previous models, we use an exponential function as the foundation to describe the decreasing trend of 2-point correlation function. Elongation ratio included in the first exponential term captures the geometry of inclusions, that is approximated by ellipsoids. In the second term, a sine function mimics the fluctuations in the 2-point correlation function, that implies the dispersion status of inclusions. X-coordinate of the first peak of fluctuation equals to $r\xafd$, the nearest center distances. The proposed analytical 2-point correlation function establishes a mapping relation bridging two different microstructure characterization approaches, and verifies the capability of the four-descriptor-set capability of capturing microstructure features across all levels. It reduces the dimensionality of microstructure design from infinite to four design variables in designing polymer nanocomposites.

Accuracy of the analytical 2-point correlation function is verified in two ways: (1) confirming boundary conditions and (2) comparing with the real 2-point correlation functions evaluated from microstructure image samples. Three boundary conditions should be satisfied for any 2-point correlation functions:

Boundary Condition 1: if

*r*= 0,*S*(0) equals to the volume fraction;Boundary Condition 2: if

*r*→ ∞,*S*(∞) equals to the square of volume fraction;Boundary Condition 3: |

*S*'(0)| equals to the interphase area amount per unit volume in 3D, or the interphase length per area in 2D.

that equals to the length of interphase per area for 2D microstructures, when the inclusion geometries are approximated by ellipsoids.

The proposed analytical function form is tested on three sample microstructures (Fig. 8). The first sample is featured by dispersion of a massive number of small round shape particles. The second microstructure is featured by the elongated shape particles (larger *e _{l}*) induced by material extrusion in processing. The third microstructure has a relatively higher volume fraction. It is observed that the descriptor-predicted correlation function (Eq. (3)) matches very well with the real 2-point correlation function for low volume fraction microstructures (microstructures 1 and 2). However, increasing volume fraction will bring larger prediction errors due to the approximations indicated by function (4). We recommend using the analytical function for microstructures with VF ≤ 15% based on our empirical studies.

## Parametric Optimization Framework for Designing Microstructural Materials

The descriptor-based microstructure representation enables the parametric design of material microstructures. Each combination of four descriptor values is one microstructure design. The first challenge in design optimization lies in how to efficiently explore the vast design space to establish the microstructure-property relation. The high computational cost of FEA of structure-property relation makes it impractical to integrate the FEA directly with an optimization search engine. The second challenge is how to search for material designs that meet multiple targets of material properties.

To overcome the above challenges, we establish a parametric design framework (Fig. 9), that is marked by two modules: the **Analysis Module** that efficiently assesses the microstructure descriptor-property relation and the **Design Module** that searches the optimal values of microstructure descriptors:

- (1)
**Analysis Module**. This module integrates DOE [41], microstructure reconstruction algorithm, advanced microstructure-property simulation, and metamodeling techniques [42] for rapid assessment of microstructure-property relations. DOE explores the design space formed by microstructure descriptors. Each sample point represents one microstructure design, for which we reconstruct one or multiple statistically equivalent microstructures and simulate the properties using FEA. The number of reconstructions depends on the window size of statistical volume element in FEA. Only one reconstruction is needed when the window size reaches representative volume element (RVE) [43,44]. After all simulations, We apply sensitivity analysis [45] to identify key descriptors that highly impact the properties to reduce the dimension of metamodels. Surrogate models, also known as metamodels [42], are created to replace the computationally expensive FEA models based on the “key descriptors” that are highly correlated with properties of interests. The Design Module uses the metamodels with reduced dimensionality for iterative optimization search. - (2)
**Design Module**. This module formulates and solves microstructure design as a design optimization problem. The design space is defined by the reduced set of microstructure descriptors. In the case that microstructure descriptors are dependent, the dependency can be modeled as design constraints. This module obtains the Pareto frontier for multiple design criteria associated with multiple material properties using heuristic search algorithms such as genetic algorithm (GA), simulated annealing, etc. based on the fast-to-compute metamodels that capture the microstructure descriptor-property relation.

## Design of Polymer Nanocomposites' Microstructure

The addition of reinforcing particles to polymer nanocomposites' matrix can lead to significant improvements in homogenized properties even at a very low filler concentration [46,47]. High impact of the quantity and morphologies of nanoparticle fillers on damping property makes it an interesting design problem for microstructure optimization.

### Design Requirements.

where G″ is the shear loss modulus (GPa), G′ is the shear storage modulus (GPa). tan δ, G′, and G″ are all functions of ω, the frequency of excitation (Hz).

The viscoelastic property is calculated by applying a sinusoidal load (stress *σ*) to a material and then measuring the resulting displacement (strain). A phase lag δ between stress and the resulting strain will be observed in materials with viscosity. For example, *δ* = 0 for a perfectly elastic solid (stress and strain are perfectly in phase); *δ* = 90 deg for a purely viscous fluid. Viscoelastic materials will show a phase lag *δ* between 0 and 90 deg. The elastic portion of viscoelastic materials is represented by the shear storage modulus that measures the stored energy; and the viscous portion is measured by the shear loss modulus that measures the energy dissipated as heat. The shear storage modulus and shear dynamic modulus is calculated from *δ* [48]. The range of loading frequency in this study is set between 6.3 × 10^{−4} and 6.3 × 10^{6 }Hz.

Corresponding to the aforementioned three design requirements, the microstructure of nanocomposites is desired to have a tan *δ* curve with low value in the low frequency domain (smaller than 1 × 10^{−1 }Hz), high value in the normal (from 1 × 10^{−1} to 1 × 10^{3 }Hz) and high frequency domains (larger than 1 × 10^{3 }Hz). Shown in Fig. 10, we choose three property characteristics as the design criteria: value of the first point on tan δ (*L*, representing damping property in low frequency domain), value of the peak on tan δ (*P*, representing normal frequency domain), and value of the last point on tan δ (*H*, representing high frequency domain). *L*, *P*, and *H* represent three different desired properties of tire material, respectively: low *L* leads to low rolling resistance; high *P* leads to high wet traction; high *H* leads to low wear. The goal is to achieve good performances on all three properties simultaneously. Thus, the design problem is formulated as a multiobjective optimization problem:

*Find a set of microstructure descriptors, s.t.: Min L; Max P (Min -P); Max H (Min -H);*

### Finite Element Modeling for Material-Property Simulation.

The 2D finite element model takes the reconstructed microstructures (binary digital image) as inputs. Each pixel becomes one element in the FEM model. Influences on the filled elastomers' properties mainly come from three aspects: the filler, the polymer matrix, and the interphase. The fillers are considered as elastic with Young's modulus of 73 (GPa), Poisson's ratio 0.19, and density 2.2 g/cm^{3}. The filled particles determine the strength and infinitesimal strain stiffness. The matrix is assigned linear viscoelastic material properties whose frequency response is given by the property curves (solid line) in Fig. 11(b), that is obtained from physical experiments. The polymer matrix provides energy dissipation and resistance of large deformation.

Due to either physical or chemical interactions, the material properties of interphase are different from those in the bulk matrix. Reinforcing fillers modify the composite's property because they result in an alteration of the matrix polymer properties at the polymer/filler interface [49]. Here, we use a double layer gradient interphase model with two-pixel thick interphase as shown in Fig. 11(a) where we assume the interphase properties decay to the bulk behavior as the distance from the fillers increases. It is assumed that the interphase viscoelastic responses are related to those of matrix by a simple shifting in the frequency domain [50] toward lower frequencies as in Fig. 11(b). The amount of shifting decays in a gradient fashion outwards the particle/polymer interface. We choose the relaxation time of the outer and inner interphase layers 10 times and 100 times that of the bulk polymer, respectively. In Fig. 11(b), the interphase moduli (dashed curves) are obtained by shifting the matrix's property curves toward low frequency. Frequency-domain analysis is run on the reconstructed microstructure with a harmonic displacement loading. In each simulation, the loss and storage moduli are assessed for a spectrum of 50 frequencies between 6.3 × 10^{−4 }Hz and 6.3 × 10^{−4 }Hz.

This paper chooses 300 × 300 as the pixel size of the microstructure images. The corresponding physical size is 10.93 *μ*m, that is considered as RVE because it is orders of magnitude larger than the carbon black particles (typically much less than 1 *μ*m in diameter). FEA simulations take the pixelated binary images as inputs to calculate the composite's frequency-dependent moduli that serve as the material properties of interest.

### Metamodel of Microstructure-Property Relation.

A single run of the 2D model on 300 × 300 pixel size structure takes roughly 15 min on 4 Quad-core Intel 5550 CPUs, 2 processors, 2.66 GHz. The computation time will increase superlinearly if the resolution of microstructure image is increased or the high fidelity 3D models are used. For example, it takes 60 h for simulating an 80^{3} voxel size cube on 2 Quad-core Intel 5550 CPUs, 16 processors, 2.66 GHz. In this paper, we demonstrate the use of metamodels for replacing high cost simulations so that microstructure optimization is computationally feasible. We present a systematic way of exploring a microstructure design space and establishing the microstructure-property relation using metamodels.

#### Exploration of the Design Variable Space Based on DOE.

DOE explores the four-dimensional design variables space (volume fraction $VF\xaf$, cluster number $N\xaf$, average elongation ratio $e\xafl$, and average nearest distances $r\xafd$). Table 2 defines the range of each design variable. It should be noted that the analytical correlation function is not used in the descriptor-based design process, so the recommended upper bound of $VF\xaf$ in analytical correlation function is not a constraint in the design formulation. The descriptors' ranges used here are slightly wider than the ranges observed from multiple SEM images of polymer nanocomposites used for tire materials.

where $l\xaf$ is the side length of the microstructure (in pixel).

With the above constraint on $N\xaf$ and $r\xafd$ values, a constrained latin hypercube sampling (cLHS) algorithm [51] is operated to spread 50 sample designs evenly throughout the four-dimensional design space (Fig. 12). The algorithm of cLHS has two steps. First, optimal latin hypercube sampling ([52]) is employed to spread all the points over the entire design space. Second, the design sequences of $N\xaf$ and $r\xafd$ are randomly switched, respectively, until no DOE point violates the constraint (using SA optimization). Each DOE point represents one material design.

#### Dimension Reduction and Metamodels of Microstructure-Property Relation.

The responses of design evaluation are the three key characteristics in property (*L*, *P*, *H*), and their relationship with four microstructure descriptors ($VF\xaf$, $N\xaf$, $r\xafd$, $e\xafl$) are modeled via Kriging model. Among the various metamodeling techniques [42], Kriging model [53] is the most widely used one due to its flexibility in fitting highly nonlinear functions and the capability of predicting the uncertainties at the locations without sample points.

In the situation of a large number of design variables (microstructure descriptors), statistical sensitivity analysis [54] should be used to identify the key microstructure descriptors and reduce the dimension of a design space. Variational sensitivity analysis such as Sobol index [55] indicates how the variation in response is apportioned to the influences of different inputs. The Sobol method is adopted here to demonstrate the dimension reduction process, even though the dimensionality of the problem is low. For the three material property responses, Fig. 13 plots and compares the total sensitivity index (TSI) for all four microstructure descriptors. High TSI means a high impact of the input variable on the response. In the dimension reduction process, when TSI > 0.2 is chosen as the threshold, $VF\xaf$, $N\xaf$, and $r\xafd$ are identified as key descriptors and treated as design variables. The influence of $e\xafl$ is neglected and its value is kept at its mean in design optimization. When TSI > 0.7 is chosen, only $VF\xaf$ and $N\xaf$ are kept as design variables. The influences of $r\xafd$ and $e\xafl$ are neglected. The reduced descriptor sets lead to lower dimensional descriptor-property relations. Taken the TSI > 0.2 case for instance, the *L*-[$VF\xaf$, $N\xaf$, $r\xafd$] relation, *P*-[$VF\xaf$, $N\xaf$, $r\xafd$] relation, and *H*-[$VF\xaf$, $N\xaf$, $r\xafd$] relation are used as objective functions in design optimization, while the value of $e\xafl$ is set to the mean of its range.

It should be noted that reducing the number of design variables based on sensitivity index values in this paper is for demonstration only. The use of such technique is more meaningful when the number of descriptors is large. In practical applications, threshold needs to be determined carefully considering the tradeoff between the reduced computational cost and the accuracy sacrificed.

#### Optimization Results and Verification.

A multi-objective optimization formulation (*MIN* L*, MIN –*P, and *MIN -*H) is established to achieve multiple targeted material properties simultaneously. We solve this problem using the multi-objective GA [56]. This design problem is solved based on two different descriptor-property models: metamodel built on the full descriptor set, and the reduced dimensional metamodel (TSI > 0.2, leaving out $e\xafl$). The Pareto solutions from both optimizations are plotted in the response space shown in Fig. 14. The optimal microstructures on the Pareto frontier vary a lot in terms of morphology. Figure 14 shows three distinctively different optimal microstructures at different locations on the Pareto frontier. The design at the upper left corner of the frontier achieves the best in low rolling resistance (low *L*); the design at the lower right corner archives the best in low wear (high *H*); the design in the middle has a high wet traction (high *P*). To select the final microstructure design, a tradeoff function among the three desired properties needs to be created.

To verify the effectiveness of Pareto solutions, we compare the property of the optimal designs on the Pareto frontier with two existing materials (called baseline designs). The microstructure images and descriptors of the two materials are shown in Fig. 15 together with the baseline designs. Figure 16 compares the properties between the baseline designs and the optimal designs. For both optimization runs using the full descriptor set and the reduced set, three optimal designs are chosen respectively from different locations on the Pareto front. It is observed that the optimal designs on the Pareto frontier outperform the baseline designs, that means the proposed design procedure successfully finds the optimal microstructure designs with superior performances. The results obtained using the full descriptor set and the reduced descriptor set are close with each other, that implies the choice of TSI > 0.2 as threshold is appropriate. We also tested TSI > 0.7 (leaving out $r\xafd$ and $e\xafl$), that shows the result is unsatisfactory due to the information loss on *P*-[$r\xafd$, $e\xafl$] and *H*-[$r\xafd$, $e\xafl$] relations.

## Conclusion

This paper presents a descriptor-based microstructure design framework. It facilitates microstructure parameterization and optimization to obtain optimal material properties. The descriptor-based microstructure characterization approach reduces an infinite dimensional microstructure design space into a small set of physically meaningful microstructure descriptors. The multiphase reconstruction algorithm enables fast and accurate reconstruction of statistically equivalent microstructures for design evaluation. We verify the accuracy of the descriptor-based characterization/reconstruction approach by analytically mapping the microstructure descriptors to the 2-point correlation function. We also establish a parametric design framework via integrating advanced FEA simulations, metamodeling, statistical sensitivity analysis, and optimization search algorithms. This framework enables efficient microstructure optimization based on microstructure-property relations.

Our proposed methodology for designing microstructural materials system has several advantages. First, descriptor-based approach shows several strengths in design representation. A reduced set of descriptors defines the microstructure design space, that can be further linked to the processing variables with good physical meanings. The proposed approach also presents a new viewpoint for microstructure analysis. Microstructure features are sorted into three different levels: composition, dispersion, and geometry. The multiphase reconstruction algorithm is featured by its strength in accuracy, efficiency, and stochasticity. Even though this algorithm is only demonstrated by 2D case studies in this paper, it is also applicable to 3D reconstruction if 3D microstructure descriptors are provided as inputs. Second, we develop a mapping relation between microstructure descriptors and the 2-point correlation function to verify the accuracy of descriptor-based characterization and reconstruction approach. The proposed analytical correlation function establishes a mathematical relation between two different microstructure characterization approaches. Third, a microstructure design framework is built via integrating state-of-the-art computational design methods for efficient exploration of the design space in the microstructure domain to obtain optimal material properties.

In future work, statistical learning techniques will be adopted to replace the empirical approach in identifying the full set of key microstructure descriptors by taking into account the impact of descriptors on material properties and the trade-off between the amount of descriptors and the computational efficiency. In addition, process-microstructure relations need to be included into the proposed material design framework to cover the whole spectrum of material design and ensure that the obtained optimal microstructure can be processed using the right processing conditions. This remains to be a future research as creating such relationship requires comprehensive empirical and theoretical studies, and a large number of test samples corresponding to different processing conditions.

## Acknowledgment

The work is supported by NSF (CMMI-0928320) and ONR, Materials Informatics Applied to Nanocomposites for Advanced Technology (A12161/N000014-10-1-0244). The authors also sincerely thank Professor Linda Schadler and Bharath Natarajan from Rensselaer Polytechnic Institute for providing the SEM images, and Dr. Amin Ajdari and Dr. Anqi Hu from Northwestern for providing the image of meshed microstructure.