Inverse Identification of Elastic Properties of Constituents of Discontinuously Reinforced Composites

This paper is devoted to determination of elastic properties of composite constituents by using an inverse identification procedure. The aim of the developed identification procedure is to compute the elastic constants of individual material phases on the basis of known properties of composite materials. The inverse problem of identification has been solved by combining an evolutionary algorithm with a micromechanical model. The paper also focuses on selection of a suitable micromechanical model for optimization which should ensure a compromise between accuracy and complexity. Two different cases have been studied: composite reinforced with short cylindrical fibers and composite reinforced with cubic particles. Moreover, Monte Carlo simulations have been carried out to expose a difference in outcome of identification which may occur when uncertain input data is considered. Obtained results show that identification is successful only when properties of composite materials with at least two different volume fractions of the reinforcement are known.


Introduction
A prediction of effective properties of composite materials in terms of their microstructural features allows to design new materials in an efficient way. Micromechanical models are useful for investigation of the influence of material properties, morphology and orientation of constituents on the final, effective composite behavior. On the other hand, the mentioned quantities must be known in order to use the micromechanical models successfully. This paper focuses on determination of elastic properties of matrix and reinforcement phases. In some cases, these properties can be determined in a straightforward way by performing standard experimental tests for matrix and reinforcement materials separately. However, some composite materials are fabricated using in situ technologies; thus, quantification of elastic properties of the individual phases may be more difficult. Examples of such materials are: aluminum-aluminum oxide [1], titanium-titanium boride [2], aluminum-titanium carbide [3], cooper-titanium carbide [4] and many others. Moreover, the properties of composite constituents may change during the material processing [5][6][7]. Extraction of specimens of single phase from the composite for ex situ mechanical testing is cumbersome; therefore, in order to measure the local material properties, advanced experimental techniques like for example nanoindentation [8,9] or micropillar compression [10,11] have been developed. However, this study focuses on a different, indirect approach which is connected with the solution of an inverse problem of identification. The inverse problem of identification may be solved by combining an optimization method with a micromechanical model; in this case, microstructural quantities are variables which are adjusted in such a way that the effective material behavior predicted by the micromechanical model fits to the experimental data obtained after testing the composite. This issue has been raised by several researchers who studied different composite materials and used different methods of optimization and micromechanical models. Burczyński and Kuś [12] analyzed composites reinforced with continuous fibers; they combined finite element based homogenization with an evolutionary algorithm. Kaiser and Stommel [13] identified amorphous and crystalline constituent properties of thermoplastic material by using an evolutionary algorithm and Mori-Tanaka micromechanical model. Beluch and Burczyński [14] studied composite reinforced with continuous fibers; they applied finite element based homogenization and evolutionary algorithm as well as artificial immune system. Herrera-Solaz et al. [15] identified properties of a single crystal of AZ31 Mg alloy from polycrystal tests by combining finite element based homogenization with Levenberg-Marquard optimization method. Comellas et al. [16] studied composites reinforced with continuous fibers by using an evolutionary algorithm and mixing theory. The other way of identification of mechanical properties of composite materials presented in the literature is based on fitting the numerical model parameters to the full-field measured displacement data [17][18][19]. This paper is devoted to inverse identification of elastic constants of composite materials with randomly distributed discontinuous reinforcement. Method developed during this study assumes that elastic constants of individual material phases can be reconstructed on the basis of Young modulus and Poisson ratio of composite determined during static tensile tests. For the purpose of this study, virtual tensile tests of composites have been conducted by using the finite element method. Furthermore, the article discusses different micromechanical approaches which can be applied for the solution of the forward problem. The inverse problem has been solved by using an evolutionary algorithm. Two different cases have been studied: composite reinforced with short cylindrical fibers and composite reinforced with cubic particles. Moreover, Monte Carlo simulations have been carried out to expose a difference in outcome of identification which may occur when uncertain input data is considered.

Optimization Problem
The identification procedure developed during this study is based on the solution of optimization problem. The objective function is defined as minimization of the relative difference between elastic constants of the composite predicted by micromechanical model in terms of constituent's elastic constants and elastic constants of the composite determined by experimental testing: where: x 1 , x 2 . . . x k are the elastic constants of the material phases (variables), y i are the elastic constants of the composite predicted by micromechanical model depending on variables, Y i -given elastic constants of the composite (input data), n denotes the number of known elastic constants of the composite. The aim of identification is presented graphically in Figure 1. The investigation is devoted to composites with randomly distributed reinforcement; thus, isotropic effective material behavior is expected. Therefore, after the tensile testing of one coupon two elastic constants (Young modulus and Poisson ratio) can be provided as the input data. However, the input data can be extended by performing tensile tests for composites with different volume fractions of the reinforcement. During the solution of optimization problem formulated in such a way it is important to avoid getting stuck in local optimum. It could be achieved by using global optimization methods like the evolutionary algorithm [20][21][22], artificial immune algorithm [23,24] or particle swarm optimization [25,26]. Other advantages of the global optimization methods over the traditional methods are: no need of computing the objective function gradient and low impact of initial values of the project variables on the optimization results. During this work, the evolutionary algorithm implemented in MATLAB software (MathWorks, Natick, MA, USA), whose simplified scheme has been presented in Figure 2, has been applied. The first step is generation of initial population which has been chosen in random way with respect to the following optimization constraints: if the variable is Young modulus and if the variable is Poisson ratio. The next step is application of the evolutionary operators which are mutation and crossover, then the objective function value is computed for each individual from the population by incorporating a micromechanical model. Finally, a selection procedure is applied, and then computations are continued until the termination condition is fulfilled. In the present paper, population consisting of 50 individuals (chromosomes) has been taken into account and termination condition is fulfilled after 100 iterations (generations).
As mentioned above, it should be pointed out that the evolutionary algorithm requires to calculate the objective function value multiple times in each iteration; therefore, efficiency of the inverse identification procedure strongly depends on efficiency of applied micromechanical model. Therefore, special emphasis must be put on the selection of a suitable micromechanical model which should ensure a compromise between accuracy and complexity. The investigation is devoted to composites with randomly distributed reinforcement; thus, isotropic effective material behavior is expected. Therefore, after the tensile testing of one coupon two elastic constants (Young modulus and Poisson ratio) can be provided as the input data. However, the input data can be extended by performing tensile tests for composites with different volume fractions of the reinforcement. During the solution of optimization problem formulated in such a way it is important to avoid getting stuck in local optimum. It could be achieved by using global optimization methods like the evolutionary algorithm [20][21][22], artificial immune algorithm [23,24] or particle swarm optimization [25,26]. Other advantages of the global optimization methods over the traditional methods are: no need of computing the objective function gradient and low impact of initial values of the project variables on the optimization results. During this work, the evolutionary algorithm implemented in MATLAB software (MathWorks, Natick, MA, USA), whose simplified scheme has been presented in Figure 2, has been applied. The first step is generation of initial population which has been chosen in random way with respect to the following optimization constraints: if the variable is Young modulus and x i ∈ 0.12, 0.4 .
if the variable is Poisson ratio. The next step is application of the evolutionary operators which are mutation and crossover, then the objective function value is computed for each individual from the population by incorporating a micromechanical model. Finally, a selection procedure is applied, and then computations are continued until the termination condition is fulfilled. In the present paper, population consisting of 50 individuals (chromosomes) has been taken into account and termination condition is fulfilled after 100 iterations (generations).
As mentioned above, it should be pointed out that the evolutionary algorithm requires to calculate the objective function value multiple times in each iteration; therefore, efficiency of the inverse identification procedure strongly depends on efficiency of applied micromechanical model. Therefore, special emphasis must be put on the selection of a suitable micromechanical model which should ensure a compromise between accuracy and complexity.

Micromechanical Modeling
Methods of micromechanics allow to predict an overall effective material properties in terms of microstructural features like for example: material properties of phases, their orientation, morphology, etc. One of the most popular methods for solving problems of micromechanics of materials is finite element analysis (FEA) which is typically connected with analysis of representative volume element (RVE) [27,28]. The RVE should contain all the necessary information about the statistical description of the microstructure, the RVE size should be large enough so that the average properties of this volume element are independent of its size and position within the material [29]. In order to find the effective material properties, homogenization procedure can be applied by computing average of field quantities over the RVE's volume. In general, the six analyses should be performed to obtain an effective stiffness tensor by applying boundary conditions which enforce the following unit strains [30] (superscript denotes the number of analysis): After FE solution of the six boundary value problems stresses are averaged in the following way: and afterwards the effective stiffness tensor can be expressed as follows:

Micromechanical Modeling
Methods of micromechanics allow to predict an overall effective material properties in terms of microstructural features like for example: material properties of phases, their orientation, morphology, etc. One of the most popular methods for solving problems of micromechanics of materials is finite element analysis (FEA) which is typically connected with analysis of representative volume element (RVE) [27,28]. The RVE should contain all the necessary information about the statistical description of the microstructure, the RVE size should be large enough so that the average properties of this volume element are independent of its size and position within the material [29]. In order to find the effective material properties, homogenization procedure can be applied by computing average of field quantities over the RVE's volume. In general, the six analyses should be performed to obtain an effective stiffness tensor by applying boundary conditions which enforce the following unit strains [30] (superscript denotes the number of analysis): After FE solution of the six boundary value problems stresses are averaged in the following way: and afterwards the effective stiffness tensor can be expressed as follows: Materials 2018, 11, 2332

of 15
The FE-based homogenization can be applied in modelling microstructures of complex geometry involving arbitrary shapes and orientation distributions of constituents; however, in such cases, it typically requires time-consuming computations [31]. Another group of micromechanical models is based on the Eshelby's fundamental solution [32], here several approaches can be distinguished like for example self-consistent schemes [33], Mori-Tanaka method [34], double inclusion method [35]. These methods provide very efficient solution in comparison with the FE based homogenization. Particularly, the Mori-Tanaka (M-T) method found wide popularity in analysis of composite materials due to good predictive capabilities [36][37][38]. The basic formulation of M-T method provides the solution for two-phase, unidirectionally reinforced materials. In this case, the effective stiffness tensor can be determined in the following way: where C m and C i are isotropic stiffness tensors of matrix and inclusion respectively, I is an identity tensor, f i is volume fraction of the inclusion and I is strain concentration tensor that depends on the Eshelby's tensor S in the following way [34]: The M-T method may be extended for composites with misaligned reinforcement by using the orientation averaging procedure where effective stiffness tensor that describes the behavior of composite with misaligned inclusions C ijkl can be determined in terms of stiffness tensor of unidirectional composite C pqrs as follows: a ip a jq a kr a ls C pqrs ψ(θ, ϕ, β) sin(θ)dθdϕdβ, where ψ(θ,ϕ,β) is the orientation distribution function defined in the Euler coordinates (θ,ϕ,β), a ij is coordinate system transformation matrix [39]. An effectiveness of the orientation-averaging procedure has been presented in numerous works [39][40][41]. The drawback of the M-T method is that it is limited to analysis of spheroidal shape of inclusions only; moreover, error of homogenization increases with increasing volume fraction of reinforcement; thus, only composites with low volume fractions of the reinforcement (approximately up to 0.25) can be successfully analyzed. However, the numerical solution of the equivalent inclusion problem, instead of using the Eshelby's tensor, allowing the M-T method to be extended so as to involve the arbitrary shapes of the inclusions. The equivalent inclusion problem relates to analysis of single inclusion embedded in a large matrix [42]. The medium is typically approximated by a rectangular prism whose finite dimensions are large enough in comparison with the size of the inclusion [42,43]. The strain concentration tensor A defines the relation between the average strain in the single inclusion embedded in infinite matrix ε i and the far field strain (macro strain) ε: Numerical determination of components of A tensor is similar to the direct FE homogenization, the same boundary conditions have to be enforced (Equation (4)) although integration is performed only over the volume of the single inclusion V i as follows: Finally, the strain concentration tensor determined numerically has the following form: and it could be substituted to the Equation (7) in order to find the effective stiffness tensor [44].

Composite Reinforced with Short Fibers
Feasibility and effectiveness of the proposed inverse identification procedure has been tested by investigation of elastic constants of matrix and fiber which are part of composite material reinforced with randomly oriented short fibers. The aim of an inverse identification procedure is to reconstruct the elastic properties of matrix and fiber based on known properties of composite. Virtual tensile tests based on the FE analysis of the RVEs have been conducted for determination of composite elastic constants in terms of known elastic properties of the phases. The RVEs containing three different fiber contents, whose geometrical models are presented in Figure 3 have been considered. Each RVE contains approximately 200 fibers whose orientations are defined randomly, and periodic boundary conditions [27] have been applied. Elastic properties of matrix and fiber, which have been applied for the virtual tensile tests, are shown in Table 1. The effective elastic constants corresponding to composite properties, obtained after direct FE homogenization, are presented in Table 2.      Selection of a suitable micromechanical model for the evolutionary optimization is very important issue as pointed out in Section 2.1. The best accuracy of identification should be provided by direct FE homogenization based on the RVE containing large number of fibers (like presented in Figure 3), although it could lead to prohibitive time of computation. On the other hand, Mori-Tanaka method coupled with orientation averaging should provide relatively short time of computations, but it is not able to consider cylindrical shape of fiber and thus it must be approximated by ellipsoidal shape. Therefore, an influence of fiber shape on effective Young modulus has been tested by computing strain concentration tensor for cylindrical fiber by using FEM and comparing the result with a pure analytical solution. The geometrical model and corresponding finite element mesh of equivalent inclusion problem for cylindrical fiber have been presented in Figure 4. Single inclusion has a volume fraction 0.001 (an influence of the volume fraction on homogenization accuracy is discussed in work [45]). The following elastic constants of the composite constituents have been taken into account: E m = 10 4 MPa, Selection of a suitable micromechanical model for the evolutionary optimization is very important issue as pointed out in Section 2.1. The best accuracy of identification should be provided by direct FE homogenization based on the RVE containing large number of fibers (like presented in Figure 3), although it could lead to prohibitive time of computation. On the other hand, Mori-Tanaka method coupled with orientation averaging should provide relatively short time of computations, but it is not able to consider cylindrical shape of fiber and thus it must be approximated by ellipsoidal shape. Therefore, an influence of fiber shape on effective Young modulus has been tested by computing strain concentration tensor for cylindrical fiber by using FEM and comparing the result with a pure analytical solution. The geometrical model and corresponding finite element mesh of equivalent inclusion problem for cylindrical fiber have been presented in Figure 4. Single inclusion has a volume fraction 0.001 (an influence of the volume fraction on homogenization accuracy is discussed in work [45]). The following elastic constants of the composite constituents have been taken into account: Em = 10 4 MPa, νm = 0.3, Ei = 105 MPa, νi = 0.2.   Strain concentration tensor computed for cylindrical fiber has the following form: and it is in close agreement with the strain concentration tensor for ellipsoidal inclusion determined analytically in terms of Eshelby's tensor: Finally, after the orientation averaging, normalized Young moduli obtained by basing on A NUM_CYLINDER (hybrid solution) and A ELLIPSOID (analytical solution) have been compared ( Figure 5). A minor difference between the results can be noticed and, therefore, usage of the pure analytical Mori-Tanaka method (accounting for ellipsoidal shape of fiber) should provide reasonable accuracy of identification in this case. .0000 0.0000 0.0000 0.0000 0.1386 Finally, after the orientation averaging, normalized Young moduli obtained by basing on A NUM_CYLINDER (hybrid solution) and A ELLIPSOID (analytical solution) have been compared ( Figure 5). A minor difference between the results can be noticed and, therefore, usage of the pure analytical Mori-Tanaka method (accounting for ellipsoidal shape of fiber) should provide reasonable accuracy of identification in this case. As an input data to the identification procedure, different combinations of known composite properties presented in Table 2 have been applied. Due to nondeterministic nature of evolutionary algorithms, three independent simulations for each case have been conducted. Results of the carried out computations and relative errors corresponding to the known data (Table 1) have been collected in Table 3. As an input data to the identification procedure, different combinations of known composite properties presented in Table 2 have been applied. Due to nondeterministic nature of evolutionary algorithms, three independent simulations for each case have been conducted. Results of the carried out computations and relative errors corresponding to the known data (Table 1) have been collected in Table 3. Obtained results show that taking as input data properties of one composite material leads to inaccurate identification results, moreover each of three independent simulations gives completely different results. Taking as input data properties of two composite materials (with two different volume fractions of the reinforcement) lead to obtaining results in good agreement with prescribed, known values. In this case independent simulations give similar results. This same trend was noticed when properties of three composite materials (with three different volume fraction of the reinforcement) have been taken as an input data.
Afterwards, an influence of the input data uncertainty on identification accuracy has been investigated by performing Monte Carlo simulations [46,47]. Gaussian distribution of the input data has been taken into account by considering mean value µ, and two different cases of standard deviation s 1 and s 2 as indicated in Table 4. For each simulation, the input data was selected randomly with probability given by the gaussian distribution (1500 simulations for each standard deviation have been carried out). Monte Carlo simulations exposed a difference in outcome of identification which may occur when uncertain input data is applied. Figures 6 and 7 present the results obtained for the standard deviation s 1 and Figures 8 and 9 present the results obtained for the standard deviation s 2 . Histograms that represents a distribution of identified quantities have been determined on the basis of Monte Carlo simulations ( Figure 10). Table 4. Statistical properties of composites reinforced with three different volume fractions of reinforcement which serves as an input data to Monte Carlo simulations.        The distributions of identified elastic constants of matrix have much smaller widths than the distributions of identified elastic constants of fibers. An irregular shape of the distribution of the Poisson ratio of fibers is caused by reaching the lower optimization constraint during simulations.  The distributions of identified elastic constants of matrix have much smaller widths than the distributions of identified elastic constants of fibers. An irregular shape of the distribution of the Poisson ratio of fibers is caused by reaching the lower optimization constraint during simulations. The distributions of identified elastic constants of matrix have much smaller widths than the distributions of identified elastic constants of fibers. An irregular shape of the distribution of the Poisson ratio of fibers is caused by reaching the lower optimization constraint during simulations.

Composite Reinforced with Cubic Particles
Afterwards, a composite reinforced with cubic particles has been analyzed. Is this case, usage of the analytical Mori-Tanaka method may lead to inaccurate results of identification since substantial difference between Young modulus estimated by using pure analytical method (accounting for the Eshelby's tensor for spherical inclusion) and hybrid method (involving cubic shape of inclusion) has been noticed ( Figure 11). A geometrical model of equivalent inclusion problem for cubic particle is presented in Figure 12a, the following elastic constants of the composite constituents have been assumed: Em = 10 4 MPa, νm = 0.3, Ei = 10 5 MPa, νi = 0.2. Virtual tensile tests based on the FE analysis of the RVEs have been conducted in the same way as previously (Section 3.1). The geometrical model of the RVE representing composite containing 15% of fibers is presented in Figure 12b. Elastic properties for matrix and particle are collected in Table 5. The effective elastic constants corresponding to composite material, obtained after direct FE homogenization, have been presented in Table 6.

Composite Reinforced with Cubic Particles
Afterwards, a composite reinforced with cubic particles has been analyzed. Is this case, usage of the analytical Mori-Tanaka method may lead to inaccurate results of identification since substantial difference between Young modulus estimated by using pure analytical method (accounting for the Eshelby's tensor for spherical inclusion) and hybrid method (involving cubic shape of inclusion) has been noticed ( Figure 11). A geometrical model of equivalent inclusion problem for cubic particle is presented in Figure 12a, the following elastic constants of the composite constituents have been assumed: E m = 10 4 MPa, ν m = 0.3, E i = 10 5 MPa, ν i = 0.2.

Composite Reinforced with Cubic Particles
Afterwards, a composite reinforced with cubic particles has been analyzed. Is this case, usage of the analytical Mori-Tanaka method may lead to inaccurate results of identification since substantial difference between Young modulus estimated by using pure analytical method (accounting for the Eshelby's tensor for spherical inclusion) and hybrid method (involving cubic shape of inclusion) has been noticed ( Figure 11). A geometrical model of equivalent inclusion problem for cubic particle is presented in Figure 12a, the following elastic constants of the composite constituents have been assumed: Em = 10 4 MPa, νm = 0.3, Ei = 10 5 MPa, νi = 0.2. Figure 11. Normalized Young moduli in terms of volume fraction of the reinforcement determined by using pure analytical and hybrid method.
Virtual tensile tests based on the FE analysis of the RVEs have been conducted in the same way as previously (Section 3.1). The geometrical model of the RVE representing composite containing 15% of fibers is presented in Figure 12b. Elastic properties for matrix and particle are collected in Table 5. The effective elastic constants corresponding to composite material, obtained after direct FE homogenization, have been presented in Table 6. Virtual tensile tests based on the FE analysis of the RVEs have been conducted in the same way as previously (Section 3.1). The geometrical model of the RVE representing composite containing 15% of fibers is presented in Figure 12b. Elastic properties for matrix and particle are collected in Table 5. The effective elastic constants corresponding to composite material, obtained after direct FE homogenization, have been presented in Table 6.   Table 7 presents results of identification obtained by using different micromechanical models and errors corresponding to the known data. As was supposed, application of the pure analytical Mori-Tanaka method leads to substantial errors, hybrid M-F/FE approach lead to much better identification accuracy in this case.

Concluding Remarks
This paper focused on determination of elastic properties of composite constituents by using an inverse identification procedure. The aim of the developed identification procedure is to compute the elastic constants of individual material phases on the basis of properties of composite material that may be measured experimentally. The inverse problem of identification has been solved by combining an evolutionary algorithm with a micromechanical model. Obtained results show that identification is successful only when properties of composite materials with at least two different volume fractions of the reinforcement are known, otherwise identification is ambiguous. The paper also focuses on selection of a suitable micromechanical model for optimization which should ensure   Table 7 presents results of identification obtained by using different micromechanical models and errors corresponding to the known data. As was supposed, application of the pure analytical Mori-Tanaka method leads to substantial errors, hybrid M-F/FE approach lead to much better identification accuracy in this case.

Concluding Remarks
This paper focused on determination of elastic properties of composite constituents by using an inverse identification procedure. The aim of the developed identification procedure is to compute the elastic constants of individual material phases on the basis of properties of composite material that may be measured experimentally. The inverse problem of identification has been solved by combining an evolutionary algorithm with a micromechanical model. Obtained results show that identification is successful only when properties of composite materials with at least two different volume fractions of the reinforcement are known, otherwise identification is ambiguous. The paper also focuses on selection of a suitable micromechanical model for optimization which should ensure a compromise between accuracy and complexity. Here, two different materials have been analyzed: in the case of composites reinforced with cylindrical fiber usage of pure analytical Mori-Tanaka method provided reasonable accuracy of identification, in the case of composite reinforced with cubic particle pure analytical method lead to substantial errors and, therefore, a hybrid homogenization method which accounts for actual reinforcement shape has to be applied. An influence of the input data uncertainty on identification accuracy has been investigated by performing Monte Carlo simulations. The Monte Carlo simulations exposed a difference in outcome of identification which may occur when input data is distributed normally. The performed numerical simulations presented the feasibility and effectiveness of proposed inverse identification procedure. On the other hand, investigation exposed identification errors that may occur, especially in the case when input data is uncertain. This study is the basis for further research connected with experimental tests.
Nonetheless, there is still a lot of work to be done to extend the identification procedure for general applications. There are several common features in engineering materials that may cause additional issues connected with the inverse identification procedure: (a) reinforcing fibers may have anisotropic properties; (b) porosity present in the material is not constant and increases with increasing volume fraction of the reinforcement; (c) orientation distribution of the reinforcement is not strictly random but depends on manufacturing process and may vary from unidirectional to random; and (d) the interface between the matrix and the reinforcement is imperfect. Moreover, the directions for future research are connected with improving the efficiency of micromechanical models for non-spheroidal inclusions, for example, by applying boundary element methods [48,49], instead of FEM, during the numerical solution of equivalent inclusion problem or developing metamodels instead of using time-consuming micromechanical models.