A Supervised Machine Learning Approach for Accelerating the Design of Particulate Composites: Application to Thermal Conductivity
Corresponding author: Azadeh Sheidaei, PhD, Aerospace Engineering Department, Iowa State University, Ames, IA 50011, United States, Tel: 515-294-2956 (O)/Fax: 515-294-3262 (O), Email: [email protected]
A Supervised Machine Learning Approach for Accelerating the Design of Particulate Composites: Application to Thermal Conductivity
Mohammad Saber Hashemi , Masoud Safdari , Azadeh Sheidaei
1* 1
Aerospace Engineering Department, Iowa State University, Ames, IA 50011, United States Aerospace Engineering Department, University of Illinois, Urbana-Champaign, IL 61820, United States
Abstract
A supervised machine learning (ML) based computational methodology for the design of particulate multifunctional composite materials with desired thermal conductivity (TC) is presented. The design variables are physical descriptors of the material microstructure that directly link microstructure to the material’s properties. A sufficiently large and uniformly sampled database was generated based on the Sobol sequence. Microstructures were realized using an efficient dense packing algorithm, and the TCs were obtained using our previously developed Fast Fourier Transform (FFT) homogenization method. Our optimized ML method is trained over the generated database and establishes the complex relationship between the structure and properties. Finally, the application of the trained ML model in the inverse design of a new class of composite materials, liquid metal (LM) elastomer, with desired TC is discussed. The results show that the surrogate model is accurate in predicting the microstructure behavior with respect to high-fidelity FFT simulations, and inverse design is robust in finding microstructure parameters according to case studies.
Keywords:
Thermal properties, Particle-reinforced composites, Flexible composites, Computational mechanics, Machine learning Introduction
Computational material design, an emerging field of study, is a powerful technique in developing advanced multifunctional materials. Accomplishing the goal of these studies depends on the appropriate representation of the material microstructure as the design variables. Microstructure characterization and reconstruction (MCR) techniques, which are generally considered to represent the microstructure, can be categorized into (1) correlation function-based methods [1–4], (2) physical descriptor-based methods [5–8], (3) spectral density function-based characterization and reconstruction by level-cutting a random field [9] or by disk-packing [10], (4) ML-based methods such as convolutional deep neural networks [11], instance-based learning [12], and encoding/decoding methods [13], and (5) texture synthesis-based methods [14–16]. Categories 1, 4, and 5 cannot be used for material design since they do not provide specific or physical design variables. Others may involve dimensional reduction due to high-dimensional representations [17], which should be cautiously studied to avoid significant information loss and decrease the structural variability. All in all, the most convenient yet capable category for material design is the physical descriptor-based method. ML methods have been used to learn the complex relationship between microstructure descriptors and their homogenized response when dealing with a massive database. For instance, Hashemi et al. [18] recently developed a novel ML-based computational framework for homogenizing heterogeneous soft materials. Furthermore, Bessa et al. [19] have proposed a framework for data-driven analysis and material systems design under uncertainty. In such frameworks, the computational cost of high-fidelity analyses is reported as the main hurdle in the data generation phase as material analyses are inherently complex for several reasons, e.g., complexities of resolving heterogeneities of the material, non-linearity of material’s response and boundary conditions, and excessive dimensionality of the design space. Reduced-order models (ROMs) could be utilized to accelerate the data generation phase. Several research works have been devoted to such developments. For example, Liu et al. [20] have developed a self-consistent clustering analysis to predict irreversible processes accurately. In this study, we focus on designing particulate composites with LM elastomer as our case study. LM composites constitute a new class of multifunctional materials with concurrently tuned thermal, dielectric, and mechanical properties. LM composite is shown to have promising applications in areas such as wearable devices, electronics, robotics, and biomedical. Carbon-based fibers in micron or nano-size limits the flexibility of the polymeric materials due to the huge difference in the properties with host polymer that results in a fracture within a few percent strain. However, LM elastomer stretches for 500% with 70% VF without fracture while maintaining high thermal conductivity [21]. Over the past few years, most research work has been focused on developing methods to synthesize LM droplets and their suspension within various matrix materials; for example, methods are developed for precise controlling of the size distribution and the volume fraction of LM droplets with a wide variety of surfactants, polymer coatings, and dispersion media [21–26]. As these methods are further refined, there will be an increasing need for computational tools to design LM composites with target material properties. The determination of the effective properties of composite materials given their specific constituents has been widely explored in the past decades. High-fidelity finite element (FE) simulations of composite materials’ response yield accurate predictions, but the associated computational time limits their applicability in the design phase. Thus, we developed a computational framework to obtain accurate and inexpensive predictions of the TC of LM composites and understand their dependence on the microstructural geometry based on optimal ML algorithms. The overview of this material discovery framework is shown in Fig. 1. Phase 1, data generation, is necessary to train ML models over a labeled dataset. Phase 2 involves finding the complex relationship between the structure and properties using appropriate ML algorithms. The first objective of this phase is finding the effective properties of a material system given its microstructural parameters, such as volume fraction (VF), size distribution, and aspect ratio (AR). The second objective is the inverse design, i.e., finding microstructure parameters given its properties as inputs. The microstructure could be realized, given its parameters or 3-D visualization through computational packing algorithms in the studied material system. Thus, Phase 3 encompasses the inference of direct structure-property relationships and generating and visualizing candidate microstructures from the inverse design framework. The computational times reported in this study are based on what has been achieved on a machine with AMD Ryzen 1950x CPU and 32 GB DDR4 RAM.
Fig. 1.
Overview of ML predictive modeling and inverse design of composites. Methods
Data generation
To have a sufficiently large dataset for advanced ML algorithms, we cannot only rely on the experimental results of the material. Even best designed experimental procedures cannot cover all feature vectors, material system characteristics in this study, required for the ML training as diverse and representative as their computational counterparts. Therefore, a robust and efficient design scheme of virtual experiments, i.e., computational simulations, was necessary to have a representative dataset of the random variables, which were the physical microstructure descriptors in this study, affecting the performance of a trained ML predictive model. The output of this phase is an open-source dataset that could open new horizons for research on the material discovery, mainly when it can be applied to similar material systems of particulate composites with different material constituents.
Design of experiment
The target variables to be predicted are an effective TC tensor of a composite given a specified set of materials for the composite constituents. The properties will be computationally measured using a homogenization technique. Since the TC is governed by a linear Laplacian equation, assuming constant material coefficients, a set of constant periodic boundary conditions can be prescribed to perform homogenization. Thus, the only remaining parameter that affects the property is the microstructural morphology. Since the studied composite was insulated LM elastomer, MCR type we chose was based on the physical descriptors. It has been shown [21–26] that LM particulates are roughly encapsulated in an ellipsoidal shape with varying aspect ratios ( 𝐴𝑅 s), following a normal distribution. To account for different 𝐴𝑅 s, two 𝐴𝑅 s for particles were considered: i) 𝐴𝑅 = 1 for spherical particles, and ii)
𝐴𝑅 ≠ 1 for all other ellipsoidal particles. Another important geometrical factor in composites is the volume fraction ( 𝑉𝐹 ) of the constituents. The physical descriptors and their bounds, as well as the numbers of particles, which are necessary for packing algorithm performance, are summarized in Table 1 where 𝑉𝐹(%) denotes the volume fraction, 𝐴𝑅 denotes the ellipsoidal aspect ratio, 𝐴𝑣𝑔(𝜇𝑚) is the average particle size,
𝑆𝑡𝑑(𝜇𝑚) is the standard deviation of the particle sizes, is the number of ellipsoidal particles, and is the number of spherical particles. The bounds of variables were selected based on an experimental work [24].
Table 1.
The bounds on the physical descriptors of the microstructure and the numbers of particles inside a pack.
𝑽𝑭(%) 𝑨𝑹 𝑨𝒗𝒈(𝝁𝒎)
𝑺𝒕𝒅(𝝁𝒎)
Lower Bounds 1 0.5 1 0.1 0.1 1 Upper Bounds 60 3.0 500 100.0 100.0 500
After identifying and limiting the microstructure parameters affecting the properties, a DoE method was used to explore the design or input variables’ domain for training machine learners efficiently. Since there was no prior knowledge of the conditional probabilities of the microstructural inputs and the property output, space-filling designs that equally cover different regions of the design space were chosen. Different optimum Latin Hypercube Samplings [27] and the Sobol sequence [28], a deterministic low discrepancy sequence, have shown a better balance between more regular distribution or randomness and being closer to a regular grid or better coverage of the input variables space [29]. Thus, we chose the Sobol sequence, which is also very fast in generating the experiment points.
Microstructure generation
Rocpack [30], a derivative of the Lubachevsky-Stillinger (LS) algorithm [31] for packing disks and spheres, was used to generate a microstructure for each point of DoE. It allows the specification of material characterization parameters presented in Section 2.1.1 and different realizations of a microstructure with unique parameters by a random initial seeding of particles. Not all experiment points generated by the Sobol sequence are consistent or physically meaningful. For instance, the continuous normal size distribution of the particles would be discretized according to the total number of particles resulting in different growth rates and final sizes for the particles. The minimum and maximum sizes of particles can be consequently determined by
𝐼𝑓 𝐸𝑥𝑝𝑟 ≥ 0 ⟹ (cid:4682)𝑟 (cid:3040)(cid:3036)(cid:3041) = 𝐴𝑣𝑔 − (cid:3493)𝐸𝑥𝑝𝑟𝑟 (cid:3040)(cid:3028)(cid:3051) = 𝐴𝑣𝑔 + (cid:3493)𝐸𝑥𝑝𝑟 𝑃𝐷𝐹(𝑟) = (cid:2869)(cid:3020)(cid:3047)(cid:3031)√(cid:2870)(cid:3095) exp (− (cid:2869)(cid:2870) (cid:4672) (cid:3045)(cid:2879)(cid:3002)(cid:3049)(cid:3034)(cid:3020)(cid:3047)(cid:3031) (cid:4673) (cid:2870) ) (cid:3017)(cid:3005)(cid:3007)(cid:2880) (cid:3117)(cid:3263) (cid:4659)(cid:4621)(cid:4621)(cid:4621)(cid:4656) (𝑟 − 𝐴𝑣𝑔) (cid:2870) = −2𝑆𝑡𝑑 (cid:2870) ln (cid:4672) (cid:3020)(cid:3047)(cid:3031)(cid:3015) √2𝜋(cid:4673) = 𝐸𝑥𝑝𝑟 (1) where 𝑃𝐷𝐹 and 𝑁 are the Gaussian probability distribution function and the number of particles, respectively. In this equation, if the minimum size 𝑟 (cid:3040)(cid:3036)(cid:3041) is lower than zero or 𝐸𝑥𝑝𝑟 is lower than zero, the parameters are physically inconsistent. The numbers of particles also determine the dimension of the microstructure in a periodic cube. This can be inferred from Eq. (2), which elucidates the implicit relationship between the physical descriptors of the microstructure given by
𝐷𝑜𝑚𝑎𝑖𝑛 (cid:3020)(cid:3036)(cid:3053)(cid:3032) = 𝑉𝑜𝑙 (cid:3021)(cid:3042)(cid:3047)(cid:3028)(cid:3039) (cid:3117)(cid:3119) = (cid:4672) (cid:3023)(cid:3042)(cid:3039) (cid:3258)(cid:3289)(cid:3278)(cid:3287)(cid:3296)(cid:3294)(cid:3284)(cid:3290)(cid:3289)(cid:3294) (cid:3023)(cid:3007) (cid:4673) (cid:3117)(cid:3119) = (cid:4672) (cid:3023)(cid:3042)(cid:3039) (cid:3254)(cid:3287)(cid:3287)(cid:3284)(cid:3291)(cid:3294)(cid:3290)(cid:3284)(cid:3279)(cid:3294) (cid:2878)(cid:3023)(cid:3042)(cid:3039) (cid:3268)(cid:3291)(cid:3283)(cid:3280)(cid:3293)(cid:3280)(cid:3294) (cid:3023)(cid:3007) (cid:4673) (cid:3117)(cid:3119) =(cid:4676) (cid:2869)(cid:3023)(cid:3007) (cid:4672)∫ (cid:3427)4𝜋 × 𝐴𝑅 × 𝑟 (cid:2870) × 𝑁 (cid:3006)(cid:3039)(cid:3039)(cid:3036)(cid:3043)(cid:3046)(cid:3042)(cid:3036)(cid:3031)(cid:3046) × 𝑃𝐷𝐹(𝑟)(cid:3431)𝑑𝑟 (cid:3019) (cid:3288)(cid:3276)(cid:3299) (cid:3019) (cid:3288)(cid:3284)(cid:3289) + ∫ (cid:3427)4𝜋 × 𝑟 (cid:2870) × (cid:3019) (cid:3288)(cid:3276)(cid:3299) (cid:3019) (cid:3288)(cid:3284)(cid:3289) 𝑁 (cid:3020)(cid:3043)(cid:3035)(cid:3032)(cid:3045)(cid:3032)(cid:3046) × 𝑃𝐷𝐹(𝑟)(cid:3431)𝑑𝑟(cid:4673)(cid:4677) (cid:3117)(cid:3119) (2) The outputs of the code, i.e., 3D realization of the packs, were given as 2-D images of the sliced 3-D microstructures in one arbitrary direction due to the isotropy. The resolution of slicing limits voxelization, and it cannot be arbitrarily increased since the FFT homogenization cost scales super-linearly with the number of voxels used. Therefore, we set the number of pixels in all directions to 300. This setting may be coarse for packs with tiny LM particles to capture the exact geometrical shape of the particles, but it resulted in a reasonable FFT computation time, about 3 hours on average for each pack. FFT homogenization to calculate effective TC
TC values of 0.29W/mK and 26.4W/mK were selected for the silicone elastomer matrix and eutectic gallium-indium (EGaIn), respectively. Conventional numerical methods such as Finite element for finding the effective properties of random heterogeneous materials suffer from their dependency on very fine and high-quality mesh conforming to intricate geometries of phases. The FFT method is shown to be efficient with voxelized representative volume elements (RVE) as no conformal meshing is required [32]. It is also superior to other numerical methods in terms of scalability
𝑂(𝑁𝑙𝑜𝑔𝑁) in complexity vs.
𝑂(𝑁 (cid:2871) ) . In a separate study, we have validated this homogenization method with the experimental results of the LM composites [33], and the reader is referred to this work for more details. ML model training
Phase 2 aims to find an efficient and optimal ML model to replace the time-consuming homogenization process. Neural networks are versatile and robust as a regression tool for modeling complex functions. Each neuron can be a nonlinear function and using different network architectures, number of neurons, number of layers, and links between neurons may arbitrarily increase their complexity. Therefore, we considered different architectures and trained them on the dataset according to the n-fold cross-validation technique. Although ReLU units generally perform better when using data ranging outside the regular interval [-1, 1], the perceptron function was Sigmoid. Therefore, the input data have been linearly normalized into [-1, 1]. The inputs were vectors of physical descriptors and other packing parameters needed for microstructure reconstruction. Simultaneously, the only output was the homogenized TC, which was the average value of the diagonal elements of the TC tensor, assuming the material system under study is almost isotropic. The whole available dataset of homogenized packs was randomly divided into five equally sized sections. The neural networks were trained five times by using a section of data, which has not been considered previously, as a test set and the rest as a training set each time. After five training processes, the average training accuracy and its standard deviation were calculated so that the performance of different architectures could be compared. The best performing network with the highest average training accuracy was chosen for the final training on the whole dataset.
Inferring complex structure-property relationship
Based on Section 2.2, A fast and reliable ML model for material properties prediction can be found to act as a surrogate of relatively expensive direct numerical solvers and to establish the direct relationship between the structure of the studied particulate material and its effective homogenized properties. However, the more demanding problem is the inverse design, which has been challenging due to inefficient and expensive methods of finding the material properties for a given microstructure, especially when dealing with the complex characterization of microstructure images with too many features. Since our studied material system could be characterized by only six features (or physical descriptors), and we have already established a reliable yet fast surrogate model for the direct structure-property relationship discussed in Section 2.2, an elitist genetic algorithm (GA) was utilized to optimize the structure according to limits imposed in experimental studies. The main prohibiting factor in evolutionary algorithms is the computational complexity due to the fitness calculation of many design points in each iteration [34]; however, the objective function in our study was calculated by the trained neural network, which is very fast in inference. The single target of such an optimization is the isotropic TC of the composite, and the design variables are the physical descriptors of Section 2.1. The objective function was selected to be the absolute difference between the ML prediction and the desired property to cast the problem into a minimization form. The flowchart of the optimization method of inverse material design is shown in Fig. 2. A green color distinguishes the evaluation step with the trained surrogate model.
Fig. 2 . The optimization flowchart for the inverse material design. Results and Discussion
Generated database
A subset of the first DoE results of almost 10,000 points based on the Sobol sequence is given in Table 2. In this table,
SID denotes the Sobol ID or the position of the parameters in the sequence, and columns 2-7 are needed for microstructure generation. The last columns contain minimum and maximum radii of the ellipsoid and spherical particles, and the domain size of the microstructure, which is required for the packing algorithm, respectively. The unit of all dimensions is micrometer ( 𝜇𝑚 ). Naturally, Sobol IDs should be 1, 2, 3, …, and the absent Sobol IDs are due to physically inconsistent sets of parameters or others for which packs were not generated under an hour time limit as discussed in Section 2.1.2. Table 2.
The first Sobol-based inputs of the material structure database; columns 2-7 are the features of the neural network, while the last column can be calculated from them by Eq. (2).
𝐒𝐈𝐃 𝐕𝐅 𝐀𝐑 𝐀𝐯𝐠
𝐒𝐭𝐝 𝐑 𝐦𝐢𝐧𝐄𝐥 𝐑 𝐦𝐚𝐱𝐄𝐥 𝐑 𝐦𝐢𝐧𝐒𝐩𝐡 𝐑 𝐦𝐚𝐱𝐒𝐩𝐡 𝐃𝐒
10 7.3 2.1 180 86.5 30.7 211 46.5 126.6 42.9 130.2 866.2 19 40.3 2.6 364 94.2 17.0 50 58.9 129.4 84.5 103.8 646.4 30 10.9 2.8 115 90.2 40.7 160 70.3 110.1 51.6 128.8 602.6
Fig. 3.
The first 1000 feasible DoE points generated by the Sobol sequence are projected on different 2D planes. Space-filling designs should cover the design space almost homogeneously while they need to maintain non-collapsing constraints. From Fig. 3 of the generated packs, it can be inferred that the criteria are met for our Sobol DoE although some generated sets of parameters were not used in the final simulations due to physical inconsistencies or some long times needed for packing. This design has an advantage of successive coverage of space along with the sequence generation so that the dataset can be successively improved, i.e., the design space can be further explored by continuing the previous number sequences. For instance, the first 50 DoE points of generated packs were shown by a red color, then the next 200 and the next 750 points were plotted by blue and green colors, respectively. Additionally, the projection of 6D points on different 2D planes, VF-Mean Size and VF-AR, did not overlap each other. The 3D visualization of a microstructure pack (purple spherical particles and yellow ellipsoidal ones) and a sample of its 2D slices (a black-and-white image) are shown in the top part of Fig. 4. The respective FFT simulation results (thermal gradient) due to thermal loading of ∆T=[0,10,0] are shown at the bottom of Fig. 4. The gradient is larger in high-concentration regions.
Fig. 4 . The FFT simulation results of a pack under a thermal loading of ∆T=[0,10,0].
Optimized surrogate model of direct structure-property relationship
As stated in Section 2, several fully connected neural network architectures were considered, trained, and compared to find a network with high expected prediction accuracy. Their characteristics are described in Table 3.
Table 3 . Neural network architectures grid-searched for ML hyperparameter optimization. Network type Layer Number of neurons 1-layer Layer 1 3 6 10 20 50 100 2-layer Layer 1 3 6 10 20 50 100 Layer 2 1 6 10 20 50 100 According to the cross-validation technique, the best network with the lowest mean squared error (MSE) is shown in Fig. 5(a). Its regression plot over the whole dataset, Fig. 5(b), shows that the accuracy is lower for large conductivity composites due to fewer DoE points covering regions of design space with higher VFs. Furthermore, the error histogram in the subplot of Fig. 5 (b) indicates that most errors with respect to the homogenized packs are relatively small. It is worth mentioning that the speed of the surrogate model in terms of a trained neural network is an order of 0.1 seconds compared with the conventional method of packing plus homogenization, which took an average time of 1+3 hours for each microstructure in our developed database.
Fig. 5. (a) The optimized fully connected neural network architecture, (b) the regression plot of the network prediction (ordinate) vs. reference values from FFT homogenization (abscissa); the subplot is the error histogram of the optimized network over the whole dataset. Following our objective of inferring the direct relationship between microstructure and its properties, several response surfaces of the studied LM composite were plotted using the trained neural network. For each surface, all microstructure features, network inputs, were fixed except for two of them, warm colors show high conductivity composites, and (a) (b) black lines on the surface are constant TC contours. Fig. 6(a) is the response surface of TC-AR-VF.
Fig. 6.
Inferring the direct Structure-Property relationship of the studied LM composite via plotting the response surfaces: (a) TC-VF-AR, (b) TC-Avg-Std projection on Avg-Std plane, (c) TC-VF-Avg, and (d) TC-VF-Number of Particles. As expected, VF has a prominent effect on the property although very high ARs, which were not studied in this work based on experimental results, may play an important role in determining material property. Since some design points are not feasible, as pointed out in Eq. (1), the projection of the TC-Avg-Std surface on the Avg-Std plane is empty in some regions of Fig. 6(b). Thus, it is deduced that the TC has a direct relationship with the Avg and Std in particle sizes; however, their effects are much less significant than that of VF. The effect of VF and Avg is illustrated in Fig. 6(c). Again, VF is shown (a) (b) (c) (d) to be the most critical factor in TC, and Avg has a negligible effect on TC. Fig. 6(d) shows that TC has been almost constant with an increasing number of particles. It is satisfactory in this study in that the standard deviation of predicted property due to the variation in microstructure size was low. In other words, the calculated domain sizes based on Eq. (2) were sufficient to define RVE sizes.
Inverse design via GA optimization
A case study was done to show how the proposed method in Section 2.3 performs and to verify its results through the data generation process discussed in Section 2.1. The goal of optimization was set to 3 J/mK heat conductivity. The best design point among the last population as well as the predicted property value and the FFT calculated one are [0.522, 1.61, 223, 61.04, 69.37, 238], 2.9, and 2.7, respectively. The error may be due to the limited resolution of reconstruction and smaller surrogate model prediction errors. Following multiple tests, it can be concluded that the inverse design method is efficient and accurate enough. Additionally, the total inverse design optimization took 1 min on average. Conclusions
This paper proposed a new supervised machine learning approach for accelerating the prediction of the thermal conductivity of the particulate composites, as well as designing a composite with the desired property. This framework has the advantages of superior computational speed compared to conventional optimization techniques. A comprehensive database for particulate composites has been generated covering the whole design space. The microstructure reconstructions based on this study’s DoE can also be used for similar heterogeneous particulate materials with different constituents. Additionally, a surrogate ML model was trained on the database to establish the direct links between the microstructure and the conductivity property and visualize them with various response surfaces in minutes compared with days for the traditional method of microstructure reconstruction and direct numerical solution. For the studied material system, the VF is far more important in determining the conductivity; however, greater particle sizes and higher ARs slightly improve TC. The smart and physically aware choice of the specified physical descriptors for MCR not only provided less-complicated modeling of structure-property links with respect to the image-based convolution neural networks which require many more training data but also connected the results of this study directly to the process phase, which is readily prepared for material synthesization. Finally, the low number of characterization features, the target TC, and a fully connected neural network as the fast surrogate model trained on our generated database enabled us to use an evolutionary optimization, GA, to explore the design space and find the physical descriptors of an LM composite which will have a given TC in about a minute.
Data availability : Supplementary data to this article can be found online at https://github.com/ms-hashemi/Insulated-LM-elastomer-conductivity
Declaration of Competing Interest : The authors declare that there are no conflicts of interest.
Acknowledgement:
This research has been funded by Iowa State University