Optimizing Thermal-Elastic Properties of C/C–SiC Composites Using a Hybrid Approach and PSO Algorithm

Carbon fiber-reinforced multi-layered pyrocarbon–silicon carbide matrix (C/C–SiC) composites are widely used in aerospace structures. The complicated spatial architecture and material heterogeneity of C/C–SiC composites constitute the challenge for tailoring their properties. Thus, discovering the intrinsic relations between the properties and the microstructures and sequentially optimizing the microstructures to obtain composites with the best performances becomes the key for practical applications. The objective of this work is to optimize the thermal-elastic properties of unidirectional C/C–SiC composites by controlling the multi-layered matrix thicknesses. A hybrid approach based on micromechanical modeling and back propagation (BP) neural network is proposed to predict the thermal-elastic properties of composites. Then, a particle swarm optimization (PSO) algorithm is interfaced with this hybrid model to achieve the optimal design for minimizing the coefficient of thermal expansion (CTE) of composites with the constraint of elastic modulus. Numerical examples demonstrate the effectiveness of the proposed hybrid model and optimization method.


Introduction
Carbon fiber-reinforced multi-layered pyrocarbon-silicon carbide matrix (C/C-SiC) composites exhibit attractive properties for thermal-structural applications, including low density, high strength, and high oxidation resistance. The multi-layered matrices consist in alternating sub-layers of pyrocarbon (PyC) and silicon carbide (SiC) [1,2]. The complicated spatial architecture and material heterogeneity of C/C-SiC composites constitute the challenge to understand their properties. Thus, discovering the intrinsic relations between the properties and the microstructures and sequentially optimizing the microstructures to obtain composites with the best possible performances becomes the key for practical applications of C/C-SiC composites.
The multi-layered matrices can be obtained by using the chemical vapor infiltration (CVI) process [3,4]. The controllable parameter in the CVI process is the layer thickness of each material. The layer thicknesses have to be properly controlled since the thickness variation of each layer affects the material microstructure and the effective properties of the composite as well. Many typical C/C-SiC composite components in aerospace engineering would be loaded in high temperature environments over hundreds even thousands of hours [5]. In such environments, the primary concern is to use materials with low thermal expansion behaviors and large elastic modulus. Motivated by this situation, optimization of the thicknesses of matrix layers within microstructure to minimize the coefficient of thermal expansion (CTE) of the unidirectional C/C-SiC composites is proposed in this paper. The constraint is imposed on the allowable elastic modulus according to real applications.
The micromechanical modeling approach, which provides the overall behavior of the composite through a finite element analysis of a unit cell model [6,7], is applied to obtain the CTE [8] and elastic modulus [9] of composites. The advantage of this approach is not only in obtaining global properties of the composite but also the behaviors that can be related to the composite microstructure. However, due to the complex multi-layered microstructure and large heterogeneity of multi-phase materials, a detailed unit cell finite element model of the unidirectional C/C-SiC composites involves a large number of elements. Generalization of the relationship between the microstructure and the overall properties of the composites using this finite element procedure is extremely difficult, especially in an optimization procedure. A new finite element mesh has to be rebuilt for each new situation and an iterative finite element analysis has to be carried out. This is extremely time consuming and computationally expensive. Thus, a hybrid approach is proposed in this paper by integrating the micromechanical model and artificial neural network for the identification of CTE and elastic modulus of the C/C-SiC composites. The artificial neural network has been extensively used in modeling composite material properties [10][11][12][13], especially for composite designing [14][15][16], as the relationship of the properties of the designed composite with its design parameters is very difficult to be represented as an explicit mathematical model. However, there is a lack of studies on applying the neural network in predicting the properties of C/C-SiC composites.
The nonlinear and non-differentiable nature of the presented optimization problem induces difficulty in using classical deterministic approaches for solutions. To solve this nonlinear optimization problem, a particle swarm optimization (PSO) algorithm [17,18] is used. The PSO algorithm belongs to the category of swarm intelligence techniques. It has only a small number of parameters that need to be adjusted, and is easy to implement. Although the PSO algorithm has been applied to a wide range of engineering problems in the literature [19][20][21][22], few applications to C/C-SiC composites are known.
In this study, a hybrid approach integrating the micromechanical model and artificial neural network is firstly proposed for the identification of CTE and elastic modulus of the unidirectional C/C-SiC composites. Predictions are compared with the results of a micromechanical model to assess the predictive capability of the proposed hybrid approach. The comparison shows that the forecast errors of the hybrid approach are inside the range of the relative fluctuations of testing samples. Although the neural network predictions partly agree with the micromechanical model, it is essential to improve the current neural network model in the future for an enhanced predicting capability. Then, a modified PSO algorithm is interfaced with the hybrid predictive model to minimize the CTE of a unidirectional C/C-SiC composite with six layers of alternating PyC and SiC matrix. The design variables are the thicknesses of matrix layers within the microstructure, and a constraint is imposed on the allowable elastic modulus. The classical PSO algorithm is modified to satisfy the constraints and the variable limits. The multi-stage penalty function method is adopted within PSO to satisfy the constraints, and the Harmony Search algorithm is used to deal with the particles that fly outside the variable boundaries.

Unit Cell Model
The architecture of the preform of unidirectional C/C-SiC composite consists of closely arranged fibers. The multi-layered PyC and SiC matrices are infiltrated within the porous fiber preforms by CVI process. Figure 1 shows the scanning electron microscope (SEM) photograph of a C/C-SiC composite [23], with the matrix consisting of alternate SiC (white color) and PyC (black color). It is clearly observed that the multi-layered matrices are distributed around the fibers. In addition, the pores are usually generated between adjacent fibers due to incomplete infiltration. For the CVI-processed composites, the research of Chateau and Gélébart [24] has indicated that the residual pores after CVI process have an important influence on the mechanical behavior of composites. Thus, for an accurate simulation of the material behavior, one must carefully introduce these manufacturing flaws in a computing scheme. However, in the present study, a highly idealized unit cell model is employed. Our purpose is to use this idealized model to develop a numerical scheme in an efficient manner for optimizing the thermal-elastic properties of composites. Thus, the presented research in this paper puts more emphasis on creating a validated and expandable optimization scheme. However, it should be noted that for an accurate microstructure design of the C/C-SiC composites, the pores and fiber positions must be carefully captured and modeled. To do this, the X-ray micro-computed tomography is an effective tool that has shown its applicability in the work of Chateau and Gélébart [24]. In our future study, for a high-quality optimization of C/C-SiC composites, a more real unit cell including the heterogeneous pore and fiber distributions would be carefully modeled.
In this paper, a unidirectional C/C-SiC composite with six layers of alternating PyC and SiC is considered as a case study. A geometrical model of the unit cell is displayed in Figure 2a. Characteristic geometric parameters of the unit cell are given: φ f is fiber diameter and d1-d6 are thicknesses of the matrix layers. The six layers of matrices are alternating PyC and SiC material layers (denoted as PyC/SiC/PyC/SiC/PyC/SiC). The unit cell model is then meshed using the 3-D twenty-node, thermal-structural coupled element (SOLID 96) of ANSYS finite element software, as depicted in Figure 2b.  For the CVI-processed composites, the research of Chateau and Gélébart [24] has indicated that the residual pores after CVI process have an important influence on the mechanical behavior of composites. Thus, for an accurate simulation of the material behavior, one must carefully introduce these manufacturing flaws in a computing scheme. However, in the present study, a highly idealized unit cell model is employed. Our purpose is to use this idealized model to develop a numerical scheme in an efficient manner for optimizing the thermal-elastic properties of composites. Thus, the presented research in this paper puts more emphasis on creating a validated and expandable optimization scheme. However, it should be noted that for an accurate microstructure design of the C/C-SiC composites, the pores and fiber positions must be carefully captured and modeled. To do this, the X-ray micro-computed tomography is an effective tool that has shown its applicability in the work of Chateau and Gélébart [24]. In our future study, for a high-quality optimization of C/C-SiC composites, a more real unit cell including the heterogeneous pore and fiber distributions would be carefully modeled.
In this paper, a unidirectional C/C-SiC composite with six layers of alternating PyC and SiC is considered as a case study. A geometrical model of the unit cell is displayed in Figure 2a. Characteristic geometric parameters of the unit cell are given: ϕ f is fiber diameter and d 1 -d 6 are thicknesses of the matrix layers. The six layers of matrices are alternating PyC and SiC material layers (denoted as PyC/SiC/PyC/SiC/PyC/SiC). The unit cell model is then meshed using the 3-D twenty-node, thermal-structural coupled element (SOLID 96) of ANSYS finite element software, as depicted in Figure 2b. For the CVI-processed composites, the research of Chateau and Gélébart [24] has indicated that the residual pores after CVI process have an important influence on the mechanical behavior of composites. Thus, for an accurate simulation of the material behavior, one must carefully introduce these manufacturing flaws in a computing scheme. However, in the present study, a highly idealized unit cell model is employed. Our purpose is to use this idealized model to develop a numerical scheme in an efficient manner for optimizing the thermal-elastic properties of composites. Thus, the presented research in this paper puts more emphasis on creating a validated and expandable optimization scheme. However, it should be noted that for an accurate microstructure design of the C/C-SiC composites, the pores and fiber positions must be carefully captured and modeled. To do this, the X-ray micro-computed tomography is an effective tool that has shown its applicability in the work of Chateau and Gélébart [24]. In our future study, for a high-quality optimization of C/C-SiC composites, a more real unit cell including the heterogeneous pore and fiber distributions would be carefully modeled.
In this paper, a unidirectional C/C-SiC composite with six layers of alternating PyC and SiC is considered as a case study. A geometrical model of the unit cell is displayed in Figure 2a. Characteristic geometric parameters of the unit cell are given: φ f is fiber diameter and d1-d6 are thicknesses of the matrix layers. The six layers of matrices are alternating PyC and SiC material layers (denoted as PyC/SiC/PyC/SiC/PyC/SiC). The unit cell model is then meshed using the 3-D twenty-node, thermal-structural coupled element (SOLID 96) of ANSYS finite element software, as depicted in Figure 2b.

Computation of the Elastic Modulus
In this study, a strain energy-based finite element approach is applied to evaluate effective elastic properties. It is assumed that each unit cell in the composites has the same deformation mode and that there is no separation or overlap between neighboring unit cells. Therefore, the periodic boundary conditions (PBC) [25][26][27] must be applied to the unit cell model. The PBC can be applied using nodes coupling (CP) and constraint equations (CE) defining in ANSYS. Then, the explicit formulations between the stiffness matrix coefficients and the strain energy of the unit cell model under specific loadings are derived. The detailed description of this method can be found in [28,29]. Here, only a basic introduction is presented.
In the elastic regime, the macroscopic behaviors of the unit cell can be characterized by the effective stress tensor σ and strain tensor ε over the homogeneous equivalent model. They are interrelated by the effective, also termed homogenized, stiffness matrix C H .
ş Ω εdΩ, and V is the volume of the unit cell. Consider the case of 3-D orthotropic materials; Equation (1) corresponds to The strain energy related to the microstructure is equal to: ż Ω 1 2 pσ 11 ε 11`σ22 ε 22`σ33 ε 33`σ12 ε 12`σ23 ε 23`σ31 ε 31 q dΩ " 1 2 pσ 11 ε 11`σ22 ε 22`σ33 ε 33`σ12 ε 12`σ23 ε 23`σ31 ε 31 q V With the help of specific loadings, the combination of Equation (2) and Equation (3) can be used to deduce the effective stiffness matrix C H for the unit cell. Suppose a unit initial strain is imposed in the y1 direction; i.e.,´p 1q ε " p1 0 0 0 0 0q T . Note that the superscript (1) represents the first load case. The corresponding average stress is then obtained by Equation (2): By replacing´p 1q σ and´p 1q ε into Equation (3), one obtains the following expression of the strain energy: The matrix coefficient C H 1111 can be derived: In the same way, demonstrations can be made for other coefficients, and all the results are summarized in Table 1. The elastic properties can be derived by inversing the elastic matrix. In practice, the considered unit cell will be discretized into a finite element model on which the initial strain will be imposed to evaluate the strain energy. Table 1. Different loadings and the coefficients of elastic matrix.

Computation of the CTE
Here, the CTE of the composites is determined by finite element computation of the unit cell with specific structural and thermal loadings [30]. As shown in Figure 1a, along the planes x 1 = 0, x 2 = 0, and x 3 = 0, the model is restricted to move in the x 1 , x 2 , and x 3 directions. Planes x 1 = l 1 , x 2 = l 2 , and x 3 = l 3 are free to move but have to remain planar in a parallel way to preserve the compatibility with adjacent cells. Suppose the deformation of the unit cell is caused by a temperature rise of ∆T. During the deformation, x i = l i becomes x i = l i + ∆l i , and the displacement, ∆l i , can be determined from the finite element analysis. The CTE in direction i then corresponds to

Experiment
To test the accuracy of the micromechanical model, three samples of the unidirectional C/C-SiC composites with different layer thicknesses were fabricated. The fiber preforms were close-packed 1K T-300 carbon fiber yearns from Nippon Toray (Tokyo, Japan). The multi-layered PyC and SiC matrix were deposited by CVI process using butane and methyltrichlorosilane (MTS) as the reactive materials in School of Material, Northwestern Polytechnical University, Xi'an, PR China. The infiltration condition of PyC was: temperature 960˝C, pressure 5 KPa, Ar flow 200 mL/min, butane flow 15 mL/min. The infiltration condition of SiC was: temperature 1000˝C, pressure 5 KPa, H 2 flow 350 mL/min, Ar flow 350 mL/min, and the molar ratio of H 2 and MTS of 10. Different layer thicknesses are obtained by controlling the infiltration time. The detailed illustration for the above process can be found in [3]. Elastic constants and CTEs of each material phase were taken from [3] and listed in Table 2. In this study, three groups of layers thicknesses were designed by controlling the infiltration time and are listed in Table 3. The fiber volume fractions of three samples were 19.7%, 19.7%, and 19.1%, respectively. Here, it should be noted that there indeed exists discrepancy between the designed thickness and the measured thickness, due to the complexity of the CVI process. However, since the purpose of this experiment is to verify the micromechanical model in an efficient manner, the discrepancy between the designed and true thickness values is neglected to simplify the experiment implementation (i.e., avoid the complex measurement of layer thicknesses in the SEM microphotographs). Uniaxial tensile tests were conducted at room temperature to obtain the longitudinal tensile modulus of the composites. Quasi-static tension tests were performed on a DNS-100 electronic universal testing system (CIMACH, Changchun, China). To measure the longitudinal and transverse CTEs of composites, DIL 402C dilatometer made by NETZSCH Company (Selb, Germany) was employed. Predictions based on the micromechanical model were compared with experimental data. The diameter of T-300 carbon fiber is 7.0 µm. The modulus and CTEs obtained in the previous tests for composite samples (denoted as A, B, and C) with different layer thicknesses, listed in Table 3, were chosen for comparison. Table 3 shows the comparison of measured and predicted modulus and CTEs for the various composite samples. It can be seen that the predicted results coincide well with the experimental data. The comparative results may highlight the predictive capacity of the proposed micromechanical model for predicting the elastic modulus and CTEs of unidirectional C/C-SiC composite. However, it should be pointed out that because rudimental cavities generated within the composite for the emission of large infiltrated by-products during the infiltration are not considered in this paper; thus, the modulus and CTEs computed numerically by the present model are larger than the experimental results.

Optimization Problem
In this paper, a unidirectional C/C-SiC composite consisting of six layers of matrices made up of alternate PyC and SiC is used as a case study. In high temperature environments, one of the common requirements is to use C/C-SiC composites with low thermal expansion behaviors and high elastic modulus. Thus, the objective of this study is to minimize the CTE of composites with elastic modulus constraints. The optimization problems given in the present study include two cases, which include, respectively, the minimization of the longitudinal and transverse CTEs. A constant fiber volume fraction of 30% is defined for the composites. Therefore, an equality constraint is imposed on the sum of the thicknesses of matrix layers. Note that in order to simplify the programming implementation, this equality constraint is transferred to inequality constraints, as illustrated in Equations (8) and (9) (∆ = 0.01). D 0 is a constant derived from the fiber volume fraction. In this study, D 0 is equal to 5.168 for a 30% fiber volume fraction. In addition, since the load-bearing capability of C/C-SiC composite structures in industrial applications is primarily related to the tensile modulus E 33 , another constraint is imposed on the allowable value of E 33 according to the real applications.
Mathematically, the optimization problems can be formulated as follows: The design variables are the thicknesses of the matrix layers. The fiber diameter ϕ f is 7.0 µm. Elastic constants and CTEs of each material phase are listed in Table 2. The upper bound of matrix layer thickness is 1.0 µm. The lower bound of thickness of matrix is set to 0.2 µm for reducing the complexity of the fabrication process.

Back Propagation (BP) Neural Network Model
The proposed micromechanical modeling approach has shown favorable predicting capability for the thermal-elastic properties of C/C-SiC composites. However, a high computational cost would be induced due to the large number of elements of the complex multi-layered microstructure. Especially in an optimization procedure, for each new situation, a finite element mesh has to be rebuilt and an iterative finite element analysis has to be carried out. This is extremely time consuming and computationally expensive. The most important benefit of an artificial neural network is the high computing efficiency. Therefore, in this study a hybrid approach integrating the micromechanical model and an artificial neural network is proposed for the determination of the CTE and elastic modulus of the C/C-SiC composites.
The BP neural network has the powerful ability of non-linear interpolation to obtain the mathematical mapping reflecting the internal law of the experimental data. In this study, a four-layer BP neural network containing one input layer, two hidden layers and one output layer is developed to construct the mapping between the layer thicknesses and the thermal-elastic properties of unidirectional C/C-SiC composites. Every neural network has exactly one input layer and one output layer. So we only need to determine the number of hidden layers. Heaton [31] summarized the capabilities of neural network architectures with various hidden layers: the hidden layer is not needed if the function is linearly separable; one hidden layer can approximate any function that contains a continuous mapping from one finite space to another; two hidden layers can represent an arbitrary decision boundary to arbitrary accuracy with rational activation functions and can approximate any smooth mapping to any accuracy. Therefore, two hidden layers are used in this study. Although two hidden layers increases the computational cost, their capability of representing functions with any kind of shape provides a promising tool for our further study. Figure 3 illustrates the BP neural network architecture used in this study. The network contains three parts: one input layer having six neurons related to the layer thicknesses; two hidden layers with 20 neurons for each one, and one output layer having three neurons representing the transverse and longitudinal CTEs α 11 and α 33 and the modulus E 33 . There are many general methods for determining the number of neurons in the hidden layers. However, these rules just provide a starting point for users to consider. For the problem considered in this study, one of the commonly-used rules gives an approximate range for the number of neurons in a hidden layer: where n is the number of neurons in the hidden layer, n i is the number of neurons in the input layer, n o is the number of neurons in the output layer, α is a constant from 1 to 10. According to the empirical equation, in this study the number of neurons in the hidden layer is between 4 and 13. Considering a BP neural network with more neurons in the hidden layer can give a higher precision solution [32]; the number of neurons in the hidden layer is finally selected as 20 in this paper. In the network, the total input, inj, received by the jth neuron in the hidden layer from all of the neurons in the preceding layer is: where N is the number of inputs to the jth neuron in the hidden layer, xi is the input from the ith neuron in the preceding layer, and wij is the connection weight from the ith neuron in the forward layer to the jth neuron in the hidden layer. The neuron then processes the input through a transfer function fs as below to produces its output outj: Before the above BP neural network system can be used to predict the thermal-elastic properties of the composite, it must be trained by the data obtained from the micromechanical model. The connection weight wij will be calculated by minimizing the error between the predictive value and the actual value during the training process. Details about the training process will be discussed in the following section.

Generation of Training Data
In this study, the training data are obtained from the micromechanical computations. In order to reflect the inner relationship between the thermal-elastic properties and the matrix layer thicknesses, a full factorial experimental design is no doubt an excellent idea. However, the full factorial experimental design means that a large number of computations (15,625 in this study) should be taken, which will obviously consume much time. Therefore, the Taguchi orthogonal array [33] which suggests using less simulation to find out the relationship between parameters is employed in this study. 25 samples designed by the L25 orthogonal array as well as another 35 samples randomly In the network, the total input, in j , received by the jth neuron in the hidden layer from all of the neurons in the preceding layer is: where N is the number of inputs to the jth neuron in the hidden layer, x i is the input from the ith neuron in the preceding layer, and w ij is the connection weight from the ith neuron in the forward layer to the jth neuron in the hidden layer. The neuron then processes the input through a transfer function f s as below to produces its output out j : Before the above BP neural network system can be used to predict the thermal-elastic properties of the composite, it must be trained by the data obtained from the micromechanical model. The connection weight w ij will be calculated by minimizing the error between the predictive value and the actual value during the training process. Details about the training process will be discussed in the following section.

Generation of Training Data
In this study, the training data are obtained from the micromechanical computations. In order to reflect the inner relationship between the thermal-elastic properties and the matrix layer thicknesses, a full factorial experimental design is no doubt an excellent idea. However, the full factorial experimental design means that a large number of computations (15,625 in this study) should be taken, which will obviously consume much time. Therefore, the Taguchi orthogonal array [33] which suggests using less simulation to find out the relationship between parameters is employed in this study. 25 samples designed by the L25 orthogonal array as well as another 35 samples randomly generated by computer, 60 samples in sum, are used to train the designed network.

Neural Network Training
During the training process, the connection weights should be calculated by minimizing the mean square error between network predictions and training data. Equation (13) [31] is used to update the connection weights iteratively. At the beginning of the training, the weights are given at random, and then they are iteratively updated until convergence to the certain values by using the gradient descent method.
where E is the mean square error and is set as 1ˆ10´4. η is the learning rate parameter controlling the stability and rate of convergence of the network, which is usually a constant between 0 and 1 and is chosen to be 0.01 in this study. Totally, the number of the connection weights to be identified is 580.
The training process takes about 1200 s of CPU time on HP personal workstation for 4.0ˆ10 5 training iterations. Figure 4 gives the variation curve of the mean square error with the iteration of connection weights (according to Equation (13)). It can be observed that with the updating of the connection weights, the mean square error has been gradually declined and converged to 1ˆ10´4. The mathematic mapping between the layer thicknesses and the CTE and elastic modulus is then stored in the trained net. The mathematic function can be expressed as: S piq " f l´ÿ w 3 f s´ÿ w 2 f s´ÿ w 1 X¯¯¯ (14) where S(i) (i = 1, 2, 3) represents the CTE and elastic modulus; X = [x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ] is the vector consisting of the thickness values of six matrix layers; f l is the linear transfer function between hidden layer 2 and the output layer; f s is the transfer function between the input layer and hidden layer 1, as well as hidden layers 1 and 2; w 1 , w 2 , w 3 represent the connection weights between the input layer and hidden layer 1, hidden layer 1 and hidden layer 2, and hidden layer 2 and the output layer, respectively. Materials 2016, 9,222 random, and then they are iteratively updated until convergence to the certain values by using the gradient descent method. (13) where E is the mean square error and is set as 1 × 10 −4 . η is the learning rate parameter controlling the stability and rate of convergence of the network, which is usually a constant between 0 and 1 and is chosen to be 0.01 in this study. Totally, the number of the connection weights to be identified is 580.
The training process takes about 1200 s of CPU time on HP personal workstation for 4.0 × 10 5 training iterations. Figure 4 gives the variation curve of the mean square error with the iteration of connection weights (according to Equation (13)). It can be observed that with the updating of the connection weights, the mean square error has been gradually declined and converged to 1 × 10 −4 . The mathematic mapping between the layer thicknesses and the CTE and elastic modulus is then stored in the trained net. The mathematic function can be expressed as: (14) where S(i) (i = 1, 2, 3) represents the CTE and elastic modulus; X = [x1, x2, x3, x4, x5, x6] is the vector consisting of the thickness values of six matrix layers; fl is the linear transfer function between hidden layer 2 and the output layer; fs is the transfer function between the input layer and hidden layer 1, as well as hidden layers 1 and 2; w1, w2, w3 represent the connection weights between the input layer and hidden layer 1, hidden layer 1 and hidden layer 2, and hidden layer 2 and the output layer, respectively.

Neural Network Testing
In order to demonstrate the ability of a neural network system to generalize the training data, the neural network is used to estimate the modulus and CTEs of the input design parameter combination. Twenty groups of layer thicknesses randomly generated by computer (not used in the training process) are used in the testing and are listed in Table 4.

Neural Network Testing
In order to demonstrate the ability of a neural network system to generalize the training data, the neural network is used to estimate the modulus and CTEs of the input design parameter combination. Twenty groups of layer thicknesses randomly generated by computer (not used in the training process) are used in the testing and are listed in Table 4.  Table 5 shows the comparison between the neural network prediction and micromechanical computation. The forecast error is defined as: where F e represents the forecast error of the prediction system. S p is the result of neural network prediction, and S m stands for the result of micromechanical computation. 2 ), we can easily obtain the relative fluctuations of α 11 and α 33 within the 20 groups of testing samples: (´10%, 10%) and (´11%, 11%). The above results clearly indicate that both the forecast errors of α 11 and α 33 fall in the range of fluctuations. Therefore we want to emphasize that although the neural network predictions partly agree with the micromechanical model, the improvement of the current neural network is still needed in future for an enhanced predicting capability. Possible refining approaches include supplementing adequate training samples with a wider range and optimizing the neural network. The running time of the prediction system is sharply decreased compared to that of the micromechanical analysis. The average running time of one micromechanical computation (including 12 finite element analyses) is about 2500 s and that of neural network prediction system is only about 0.001 s.

Particle Swarm Optimization Algorithm
In a PSO algorithm, each individual of the swarm is considered as a flying particle in the design space that has a position and a velocity. The particles remember the best position that they have seen during the flying. Members of a swarm remember the location where they had their best success and communicate good positions to each other, then update their own position and velocity based on these good positions as follows: where V i and X i represent the velocity and the position of the ith particle, respectively (the subscripts k and k + 1 refer to the recent and the next iterations, respectively). P i is the best previous position of the ith particle and P g is the best global position among all the particles in the swarm. ω is the inertia weight controlling the impact of the previous history of velocities on the current velocity and is set to 0.875 in this study. c 1 and c 2 are acceleration constants indicating the stochastic acceleration terms which pull each particle towards the best position attained by the particle or the best position attained by the swarm. In this work, c 1 = 2 and c 2 = 2 are chosen. r 1 and r 2 are two random numbers between 0 and 1. Most optimization problems include the problem-specific constraints and the variable limits. For the present optimization, the problem-specific constraint is the elastic modulus, and the variable limits are design bounds of the thicknesses of matrix layers. If the particle flies out of the variable boundaries, the solution cannot be used even if the problem-specific constraint is satisfied, so it is essential to make sure that all of the particles fly inside the variable boundaries, and then to check whether they violate the problem-specific constraint.

Harmony Search Scheme: Handling the Variable Limits
A method introduced by Li et al. [34] dealing with the particles that fly outside the variable boundaries is used in the present study. This method is derived from the harmony search (HS) algorithm. In the HS algorithm, the harmony memory (HM) stores the feasible vectors, which are all in the feasible space. The harmony memory size determines how many vectors can be stored. A new vector is generated by selecting the components of different vectors randomly in the harmony memory. Undoubtedly, the new vector does not violate the variables boundaries. When it is generated, the harmony memory will be updated by accepting this new vector if it gets a better solution and deleting the worst vector. Similarly, the PSO algorithm stores the feasible and "good" vectors (particles) in the p best swarm, as does the harmony memory in the HS algorithm. Hence, the vector (particle) violating the variable boundaries can be generated randomly again by such a strategy: if any component of the current particle violates its corresponding boundary, then it will be replaced by selecting the corresponding component of the particle from p best swarm randomly. To highlight the presentation, a schematic diagram is given in Figure 5 to illustrate this strategy.
Materials 2016, 9,222 12 replaced by selecting the corresponding component of the particle from pbest swarm randomly. To highlight the presentation, a schematic diagram is given in Figure 5 to illustrate this strategy.

Penalty Functions Method: Handling the Problem-Specific Constraints
The most common method of handling the constraints is the use of a penalty function. The constrained problem is transformed to an unconstrained one, by penalizing the constraints and building a single objective function. Hence, the optimization problem becomes one of minimizing the objective function and the penalty together. In this paper, a non-stationary, multi-stage penalty function method implemented by Parsopoulos and Vrahatis [35] is adopted for constraint handling with PSO. The penalty function is where f(X) is the original objective function to be optimized. h(k) is a penalty value which is modified according to the algorithm's current iteration number k and is usually set to    h k k . H(X) is a penalty factor defined as: where qi(X) is a relative violated function of the constraints, which is defined as qi(X) = max {0, gi(X)}

Results
The optimization problems illustrated in Section 3 are implemented by using the proposed hybrid approach and PSO algorithm. For all the optimization problems, a population of 50 particles is used. The stopping criterion can be defined based on the number of iterations without an update

Penalty Functions Method: Handling the Problem-Specific Constraints
The most common method of handling the constraints is the use of a penalty function. The constrained problem is transformed to an unconstrained one, by penalizing the constraints and building a single objective function. Hence, the optimization problem becomes one of minimizing the objective function and the penalty together. In this paper, a non-stationary, multi-stage penalty function method implemented by Parsopoulos and Vrahatis [35] is adopted for constraint handling with PSO. The penalty function is where f(X) is the original objective function to be optimized. h(k) is a penalty value which is modified according to the algorithm's current iteration number k and is usually set to h pkq " ? k. H(X) is a penalty factor defined as: where q i (X) is a relative violated function of the constraints, which is defined as q i (X) = max {0, g i (X)} (Note that g i (X) is the constraint), θ pq i pXqq is an assignment function, and γ pq i pXqq is the power of the penalty function. Regarding the [10], the following values are used for the penalty function: If q i pXq ă 1, then γ pq i pXqq " 1; otherwise γ pq i pXqq " 2; If q i pXq ă 0.001, then θ pq i pXqq " 10; else if q i pXq ď 0.1, then θ pq i pXqq " 20; else if q i pXq ď 1, then θ pq i pXqq " 100; otherwise θ pq i pXqq " 300.

Results
The optimization problems illustrated in Section 3 are implemented by using the proposed hybrid approach and PSO algorithm. For all the optimization problems, a population of 50 particles is used. The stopping criterion can be defined based on the number of iterations without an update in the best values of the swarm or the number of iterations the algorithm executes. Although the latter is not a real physical stopping criterion, it is quite easy in programming implementation and hence is widely used in PSO algorithms. In this work, the maximum number of iterations is limited to 100 and is adopted as the stopping criterion. Figure 6 provides the convergence rates of the optimization procedure for minimizing the longitudinal CTE. The algorithm achieves the best solution after about 50 iterations. The longitudinal CTE has been effectively reduced to 2.89ˆ10´6/˝C. The convergence of the design variables during the iterations is shown in Figure 7. It is observed that the thicknesses of the first (PyC) matrix layer increases to its higher bound. The third (PyC) and last (PyC) matrix layers both reach median values between the lower and higher bounds. The second (SiC), fourth (SiC), and fifth (PyC) matrix layers all iterate to values near the lower bound. The final optimized thicknesses are 0.999/0.259/0.557/0.215/0.276/ 0.525 µm.
Materials 2016, 9,222 13 in the best values of the swarm or the number of iterations the algorithm executes. Although the latter is not a real physical stopping criterion, it is quite easy in programming implementation and hence is widely used in PSO algorithms. In this work, the maximum number of iterations is limited to 100 and is adopted as the stopping criterion. Figure 6 provides the convergence rates of the optimization procedure for minimizing the longitudinal CTE. The algorithm achieves the best solution after about 50 iterations. The longitudinal CTE has been effectively reduced to 2.89 × 10 −6 /°C. The convergence of the design variables during the iterations is shown in Figure 7. It is observed that the thicknesses of the first (PyC) matrix layer increases to its higher bound. The third (PyC) and last (PyC) matrix layers both reach median values between the lower and higher bounds. The second (SiC), fourth (SiC), and fifth (PyC) matrix layers all iterate to values near the lower bound. The final optimized thicknesses are 0.999/0.259/0.557/0.215/0.276/0.525 μm.    Figure 9. It is observed that the thicknesses of the first (PyC) and third (PyC) matrix layers both iterate to median values between the lower and higher bounds. The second (SiC) Materials 2016, 9,222 13 in the best values of the swarm or the number of iterations the algorithm executes. Although the latter is not a real physical stopping criterion, it is quite easy in programming implementation and hence is widely used in PSO algorithms. In this work, the maximum number of iterations is limited to 100 and is adopted as the stopping criterion. Figure 6 provides the convergence rates of the optimization procedure for minimizing the longitudinal CTE. The algorithm achieves the best solution after about 50 iterations. The longitudinal CTE has been effectively reduced to 2.89 × 10 −6 /°C. The convergence of the design variables during the iterations is shown in Figure 7. It is observed that the thicknesses of the first (PyC) matrix layer increases to its higher bound. The third (PyC) and last (PyC) matrix layers both reach median values between the lower and higher bounds. The second (SiC), fourth (SiC), and fifth (PyC) matrix layers all iterate to values near the lower bound. The final optimized thicknesses are 0.999/0.259/0.557/0.215/0.276/0.525 μm.    Figure 9. It is observed that the thicknesses of the first (PyC) and third (PyC) matrix layers both iterate to median values between the lower and higher bounds. The second (SiC)

Conclusions
In this study, optimal design of unidirectional C/C-SiC composites with respect to the thermalelastic properties is obtained by the use of the non-gradient PSO algorithm. A hybrid methodology using a micromechanical model and a BP neural network is firstly proposed for predicting the elastic modulus and CTEs. Numerical results demonstrate its ability to find out the highly non-linear relationship between the multi-layers thicknesses and the CTEs and elastic moduli. However, it should be mentioned that the forecast errors of the presented neural network model are in the range of the relative fluctuations of testing samples. Therefore, although the neural network predictions partly agree with the micromechanical model, the improvement of the current neural network model is still needed in the future for an enhanced predicting capability.
Then, an optimization scheme which combines a PSO algorithm and the hybrid methodology was used to minimize the CTE of composites with the constraint of elastic modulus by designing the thicknesses of matrix layers. The minimization procedures of the longitudinal and transverse CTEs generate quite different thicknesses of matrix layers. The final optimized thicknesses are

Conclusions
In this study, optimal design of unidirectional C/C-SiC composites with respect to the thermalelastic properties is obtained by the use of the non-gradient PSO algorithm. A hybrid methodology using a micromechanical model and a BP neural network is firstly proposed for predicting the elastic modulus and CTEs. Numerical results demonstrate its ability to find out the highly non-linear relationship between the multi-layers thicknesses and the CTEs and elastic moduli. However, it should be mentioned that the forecast errors of the presented neural network model are in the range of the relative fluctuations of testing samples. Therefore, although the neural network predictions partly agree with the micromechanical model, the improvement of the current neural network model is still needed in the future for an enhanced predicting capability.
Then, an optimization scheme which combines a PSO algorithm and the hybrid methodology was used to minimize the CTE of composites with the constraint of elastic modulus by designing the thicknesses of matrix layers. The minimization procedures of the longitudinal and transverse CTEs generate quite different thicknesses of matrix layers. The final optimized thicknesses are

Conclusions
In this study, optimal design of unidirectional C/C-SiC composites with respect to the thermal-elastic properties is obtained by the use of the non-gradient PSO algorithm. A hybrid methodology using a micromechanical model and a BP neural network is firstly proposed for predicting the elastic modulus and CTEs. Numerical results demonstrate its ability to find out the highly non-linear relationship between the multi-layers thicknesses and the CTEs and elastic moduli. However, it should be mentioned that the forecast errors of the presented neural network model are in the range of the relative fluctuations of testing samples. Therefore, although the neural network predictions partly agree with the micromechanical model, the improvement of the current neural network model is still needed in the future for an enhanced predicting capability.
Then, an optimization scheme which combines a PSO algorithm and the hybrid methodology was used to minimize the CTE of composites with the constraint of elastic modulus by designing the thicknesses of matrix layers. The minimization procedures of the longitudinal and transverse CTEs generate quite different thicknesses of matrix layers. The final optimized thicknesses are 0.999/0.259/0.557/0.215/0.276/0.525 µm for minimizing the longitudinal CTE, while for minimization of the transverse CTE the final optimized thicknesses are 0.589/0.301/0.510/0.349/0.703/0.279 µm.
Here, we want to mention that the emphasis of this work was to develop an effective optimization scheme for optimizing the thermal-elastic properties of unidirectional C/C-SiC composites. Numerical examples have shown the effectiveness of the proposed method. However, to obtain a variation law of the multi-layer thicknesses for the thermal-elastic design of C/C-SiC composites, more optimization cases would need to be investigated in our further studies.