- Email: [email protected]
epoxy composites using homogenization method and Monte Carlo simulation
Renewable Energy 65 (2014) 219e226
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Probabilistic analysis for mechanical properties of glass/epoxy composites using homogenization method and Monte Carlo simulation Seung-Pyo Lee a, Ji-Won Jin b, Ki-Weon Kang b, * a b
R&D Center, ILJIN Global, 128-5, Samsung-dong, Kangnam-go, Seoul 135-875, Republic of Korea Department of Mechanical Engineering, Kunsan National University, Kunsan, Jeonbuk 573-701, Republic of Korea
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 February 2013 Accepted 12 September 2013 Available online 1 October 2013
In this paper, we used a homogenization method to estimate the equivalent properties of composite material, which depend on a microstructure and behavior of constituents (fiber and matrix). The estimated results conform well to the results of the rule of mixtures, which is a conventional method for estimating equivalent properties of composite materials. In addition, to assess the uncertainty of the equivalent properties of the composite materials according to the variability in the basic properties of the constituents, a probabilistic analysis using a homogenization-based Monte Carlo simulation was carried out. In this way, the variation of basic properties of the constituents was identified to have a significant effect on the uncertainty of the equivalent properties of composite materials. Moreover, a sensitivity of the properties of the constituents on the equivalent properties was assessed through a correlation and regression analysis. 2013 Elsevier Ltd. All rights reserved.
Keywords: Glass/epoxy composites Homogenization Mechanical properties Microstructure Monte Carlo simulation
1. Introduction Fiber-reinforced composites have superior mechanical properties such as specific stiffness and specific strength, and they have enabled us a stiffness design through changes in the stacking sequence and fiber volume fraction. Then, such composites have been used extensively in aerospace, automotive, and wind energy conversion industries [1]. Composite materials have a structure wherein reinforcing fiber and a matrix, each with unique characteristics, are physically coupled on a microscopic scale. The complicated microstructure and the behaviors of the constituents lead to the inherent characteristics of composite materials, such as strong heterogeneity and anisotropy [2,3]. And, because the basic properties of their constituents must necessarily be dispersed, there is a considerable uncertainty in the mechanical properties of composite materials [4]. Therefore, it is very important to assess the mechanical properties in a manner that takes into consideration the behavior of the microstructure and constituents of a composite material, as well as the dispersion of their basic properties. The homogenization method is one of the most useful methods to assess the mechanical properties of composite materials, which
* Corresponding author. Tel.: þ82 63 469 4872; fax: þ82 63 469 4727. E-mail address: [email protected] (K.-W. Kang). 0960-1481/$ e see front matter 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.renene.2013.09.012
takes into consideration the behavior of the constituents and microstructure [5]. The homogenization method is a procedure whereby a heterogeneous media is converted into an equivalent material model that is energetically equivalent to the heterogeneous media. On the basis of the fundamental concepts of the homogenization method, Orlik [6] applied the theory of the two-scale convergence to homogenization of initial flow stresses and hardening constants in exponential hardening laws for elasto-plastic composites with a periodic microstructure. Ji-Wei et al. [7] applied the asymptotic expansion homogenization (AEH) method to simulate the non-linear behavior of the cumulative damage to various types of unit cell. This allowed them to research the tensile strength of 3D braided composites from a microscopic point of view. Li [8] also developed a non-linear algorithm together with the finite element analysis and homogenization method to implement kinematic limit analysis for a microstructure and the macroscopic strength of a composite material with anisotropic constituents. And others have also carried out researches that coupled a stochastic approach to the homogenization method, which made it possible to consider the behavior of the constituents and microstructures together with the uncertainty in the basic properties of the constituents. Kaminski et al. [9] have proposed the stochastic second order and second moment analysis to analyze the first two probabilistic moments of the effective elasticity tensor components for homogenization of the two-phase periodic composite structure.
220
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
ckl i
Nomenclature
sεij xεj
fi uεi nj Fi eij Dεijkl ε ~ ð1Þ u i
stress generalized coordinate system for the macroscopic scale body force displacement normal vector traction strain constitutive matrix scale parameter for asymptotic series approximation arbitrary additive constant in yi
Sakata et al. [10] have analyzed the influence of microscopic uncertainty on a homogenized macroscopic elastic properties of composite materials using higher order perturbation-based homogenization method. Also, Kaminski et al. [11] have proposed the generalized stochastic perturbation method to analyze stochastic aging process of the metallic fiber reinforced composites. They have used the computational homogenization procedure based on finite element analysis. Despite offering a useful depiction of the uncertainty of the macroscopic mechanical properties due to the dispersion of the microscopic properties, these studies are of limited applicability and the difficulties in their mathematical formulations. The present study aims to probabilistically estimate the equivalent properties of glass fiber-reinforced composite materials, by using a homogenization method and Monte Carlo simulation. This first made it possible to estimate the homogenized equivalent properties, which take into consideration the microstructure and the behavior of constituents. The results were then assessed in comparison to the properties estimated from an ROM (rule of mixtures). The homogenization method-based Monte Carlo simulation also allowed us to estimate the variability of the equivalent properties of composite materials, caused by the dispersion in the properties of the constituent, as well as to assess how this variability is impacted by properties of the constituents. 2. Theoretical background and analysis procedure 2.1. Homogenization method To compare the dimensions of the composite structure composed of periodic heterogeneous unit cell such as wind turbine composite blade shown in Fig. 1, let us assume that the dimensions of the unit cell are relatively very small, with a value of ε. Given an appropriate load and boundary conditions for such a structure, there will be rapid changes in the stress and strain as a result of the
DH ijmn yi E11,eq E22,eq G12,eq Ef Em Gf Gm
nf nm
characteristic function equivalent material property tensor second coordinate system for the microscopic scale longitudinal elastic modulus of unit cell transverse elastic modulus of unit cell shear modulus of unit cell elastic modulus of fiber elastic modulus of matrix shear modulus of fiber shear modulus of matrix Poisson’s ratio of fiber Poisson’s ratio of matrix
repetition of the unit cell. In other words, the stress and strain undergo very rapid changes in the very small circumference, ε, of a point, x, in the structure. To portray this effectively, the macroscopic and microscopic levels can be defined in terms of x and y ¼ x/ε, respectively [12]. In conventional elasticity, a body subject to body force fi should satisfy the following equilibrium equation [12e14]
vsεij vxεj
þ fi ¼ 0 in U:
(1)
And the prescribed displacement uεi on v1 U and traction Fi on the boundary v2 U can be stated as
uεi ¼ 0 on v1 U;
(2)
sεij nj ¼ Fi on v2 U:
(3)
The strainedisplacement and the stressestrain relationships can be written as ε 1 vuεi vuj eij ðu Þ ¼ þ 2 vxεj vxεi
!
ε
(4)
sεij ¼ Dεijkl ekl ðuε Þ
(5)
The displacements can be approximated by an asymptotic series representation in ε given by Guedes and Kikuchi [12] ð0Þ
ð1Þ
uεi ðxÞ ¼ ui ðx; yÞ þ εui ðx; yÞ þ /:
(6)
Substituting eq. (6) into eq. (4) and then substituting the result into eq. (5) enables us to obtain an equation that is dependent on powers of ε. Each term of the same power of ε is equal to zero to ensure that the asymptotic series approximation is valid as ε approaches zero. Hence, ð0Þ
vuk v D ¼ 0 vyj ijkl vyl
(7)
ð0Þ
ð0Þ
ð1Þ
vuk vuk vu v v D þ D þ k vxj ijkl vyl vyj ijkl vxl vyl
ð0Þ
Fig. 1. Wind turbine composite blade.
ð1Þ
vuk vu v D þ k vxj ijkl vxl vyl
!
! ¼ 0;
ð1Þ
(8)
ð2Þ
vuk vu v D þ k þ vyj ijkl vxl vyl
! þ fi ¼ 0
(9)
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
ð1Þ
ui
¼ ckl i
221
ð0Þ
vuk ~ ið1Þ ðxÞ; þu vxl
(10)
ð1Þ
~ i ðxÞ is constant and ckl in which u i is a characteristic function that is the solution of the auxiliary equation given by
Z Dijkl Y
vcmn k vyi dy ¼ vyl vyj
Z Dijmn Y
vyi dy cyi ˛VY : vyj
(11)
ð1Þ
Let us assume that the solution of ui takes the form as following [12,13]. Substituting eq. (10) into eq. (9) results in eq. (12). ð0Þ
v H vum D þ fi ¼ 0 vxj ijmn vxn
(12)
The equivalent material property tensor DH ijmn can be written as [12,13]
DH ijmn
1 ¼ jYj
# Z " vckl p Dijkl Dijpq dy vyq
(13)
Y
2.2. Analysis procedure In composite materials, the distribution of fiber and matrix is likely to be random across the cross-section as shown in Fig. 2. Since it is difficult to model random fiber arrangement, many micromechanical models [15,16] assume a periodic arrangement of constituents for which unit cell can be isolated. The periodic fiber and matrix arrangements are the square and the hexagonal arrays. The unit cell has the identical elastic properties and volume fraction of fiber as the composite materials. Thus, a micromechanics-based finite element analysis for a unit cell can lead to the mechanical properties of a whole structure, which comprises the repeated rectangular unit cell in Fig. 3. In this paper, accordingly, the unit cell is defined as the square-packed array in Fig. 4 and the homogenization method is applied to the unit cell to estimate the equivalent mechanical properties of the E-glass/epoxy composite.
Fig. 3. Idealized array cells of unidirectional composites; square-packed array.
In order to calculate the equivalent mechanical properties of the unit cell, the present paper utilizes an ANSYS SOLID45 element (an eight-node structural solid element) with eight nodes in one element and three degrees of freedom in one node. The total number of nodes herein is 386 and the total number of elements is 180 (refer to Fig. 5). Periodic boundary conditions are applied as eqs. (14)e(16) to the boundaries of the unit cell in order to apply the homogenization method. The mutually facing sides of the unit cell have the same behavior, on the basis of these boundary conditions [13]. jk cjk ðy1 ; y2 ; y3 Þ ¼ ci ðy1 þ Y1 ; y2 ; y3 Þ i
(14)
cjk ðy1 ; y2 ; y3 Þ ¼ cjk ðy1 ; y2 þ Y2 ; y3 Þ i i
(15)
jk cjk ðy1 ; y2 ; y3 Þ ¼ ci ðy1 ; y2 ; y3 þ Y3 Þ i
(16)
The characteristic function ckl i given in eq. (11) is the matrix. Each column of this matrix matches a column of a constitutive matrix, Dijkl. That is, eq. (11) has the columns comprising the
Fig. 2. Microscopic configuration of multiaxial glass/epoxy composites.
Fig. 4. Geometric models of square unit cells for glass/epoxy composites.
222
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
h¼
g1 gþ2
(22)
In addition, a symbol g is a mechanical property ratio between the fiber and matrix and is represented by the following equation:
g¼
Pf Pm
(23)
And 2 is the variable for the partitioning effect of the fiber and matrix. When 2 ¼ 0 in eq. (22), eq. (21) is expressed as follows:
1 V Vm ¼ fþ P Pf Pm
(24)
When 2 ¼ N and h ¼ 0, eq. (21) is expressed as follows:
P ¼ Pm 1 Vf þ Pf Vf
(25)
3. Results and discussion 3.1. Equivalent mechanical properties using homogenization method
Fig. 5. Finite element model of square unit cell for glass/epoxy composites.
solution of a constitutive matrix, in a three-dimensional elasticity problem. The solution for these six columns can be obtained from the six columns in eq. (11). 2.3. Rule of mixture (ROM) The mathematical equations for the simple rule-of-mixtures (SROM), which is another method for the equivalent mechanical properties, are as follows.
E11;eq ¼ Ef Vf þ Em 1 Vf
n12;eq ¼ nf Vf þ nm 1 Vf
(17) (18)
1 V 1 Vf ¼ fþ E22;eq Ef Em
(19)
1 V 1 Vf ¼ fþ G12;eq Gf Gm
(20)
where the subscripts f and m represent the fibers and matrix, respectively, and Vf signifies the volume fraction of the fibers. Hyer et al. [17] determined the mechanical properties using a structural analysis of the unit cell and a simple rule-of-mixtures method. In a comparison of the results, the values for E22,eq and G12,eq showed a large difference. Thus, in order to correct for this, the modified rule of mixtures (MROM) proposed by Halphin [18] was introduced. The equation in the modified rule of mixtures is as follows.
P 1 þ 2hVf ¼ Pm 1 hVf
The equivalent mechanical properties of the glass/epoxy composites can be calculated by means of the homogenization method mentioned in Section 2.1 and Section 2.2 for the unit cell, which comprise a matrix and fiber. For this, the characteristic function in eq. (11) was first solved, the result of which was substituted into eq. (13) to calculate the equivalent mechanical properties. The geometry of the unit cell that was used to calculate the equivalent mechanical properties is shown in Fig. 4, and the boundary conditions were as shown in eqs. (14)e(16). Respective CE commands in the ANSYS, were applied to mutually facing sides to express the boundary conditions in eqs. (14)e(16). And the basic mechanical properties of the E-glass fiber and epoxy resin were used, as indicated in Table 1. The load conditions used to solve the characteristic function were based on the right-side operand in eq. (11). Because the above equations include an element strain matrix and material properties matrix, it is not possible to express these using a general method for load input. Thus, in the present study, this was resolved by taking the load terms as the initial stress terms and utilizing the user subroutine USTRESS in ANSYS to represent the initial stress terms as per eq. (11). In addition, APDL [19] in ANSYS was used to calculate the equivalent mechanical properties from eq. (13). The equivalent mechanical properties calculated using the above procedure of the homogenization method is summarized in Table 2. To verify the homogenization method, the mechanical properties were determined using the rule of mixtures and compared with the results by the homogenization method. The volume fraction of the fiber was varied from 0.1 to 0.7. As shown in Fig. 6, E11,eq and n12,eq are consistent with the SROM results, and E22,eq and G12,eq are consistent with the MROM results presented by Hyer [17]. Herein, the coefficient z can be calculated using a finite element analysis of the unit cell. In the present study, when z ¼ 1.8 and z ¼ 1.2, the equivalent properties E22,eq and G12,eq were found to be consistent with the homogenized results. Based on the above
Table 1 Material properties of constituents.
(21)
In eq. (21), P signifies E22,eq or G12,eq, which are the mechanical properties of the composites, and Pm is a mechanical property of the matrix. Herein, h is defined as follows:
Elastic modulus, E [GPa] Shear modulus, G [GPa] Poisson’s ratio, n
Fiber; E-glass
Matrix; epoxy
72.40 28.96 0.25
2.45 0.91 0.35
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
223
Vf
E11,eq [GPa]
E22,eq [GPa]
G12,eq [GPa]
G12,eq [GPa]
n12,eq
n23,eq
materials according to the variation of the basic mechanical properties. The homogenization based Monte Carlo simulation has three main steps:
0.1 0.2 0.3 0.4 0.5 0.6 0.7
9.369 16.286 23.203 30.119 37.035 43.951 50.868
3.170 3.889 4.855 6.214 8.196 11.278 16.838
1.095 1.326 1.617 2.001 2.541 3.395 5.090
1.051 1.203 1.377 1.595 1.897 2.382 3.406
0.338 0.326 0.315 0.304 0.294 0.284 0.273
0.460 0.453 0.423 0.380 0.326 0.266 0.212
1. Generate a randomly distributed set of basic properties of the constituents. 2. Perform calculations based on the homogenization method using the generated set of basic properties of the constituents. 3. Determine probability from a large number of repetitions.
results, it can be stated that the homogenization method can lead to the equivalent mechanical properties using the basic mechanical properties of the constituents. 3.2. Homogenization based Monte Carlo simulation
Normalized Longitudinal Modulus, E11,eq/Ef
The above-presented homogenization method can estimate the equivalent mechanical properties of composites using the basic properties of the constituents. However, since the mechanical properties of constituents are random in nature, one can observe great differences in the macroscopic mechanical properties of composites among the samples in the same condition. The mechanical properties of composites should, therefore, be treated in a probabilistic, rather than a deterministic way [1,20]. In this study, a homogenization based Monte Carlo simulation was applied to assess the uncertainty of the macroscopic mechanical properties, i.e., the equivalent properties of the composite 1.0
Equivalent longitudinal modulus (E11,eq) 0.8
Results by SROM Results by Homogenization
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
For the sampling of basic properties and assessment of uncertainty of the equivalent properties, the commercial code HyperStudy [21] was used. The sampling method used here was the Hamersley method [21], which is one type of quasi-Monte Carlo method and is proven to be effective for a small number of samples; the sample number was 100. The selected design variables (basic properties of the constituents) were the elastic modulus (Ef and Em) and Poisson’s ratio (nf and nm) of the fiber and matrix. Table 3 shows the design variables for the homogenization based Monte Carlo simulation. Herein, the term of CoV is defined as the coefficient of variation and can be calculated from eq. (26).
CoV ¼
standard deviation mean
10
Equivalent transverse modulus (E22,eq) Results by SROM Results by MROM Results by Homogenization
8
6
4
2
0 0.0
0.6
0.8
2.0
10
Equivalent shear modulus (G12,eq)
Normalized Poisson's Ratio, ν12,eq/ν12,f
Normalized Shear Modulus, G12,eq/Gm
0.4
(b) Equivalent transverse modulus
(a) Equivalent longitudinal modulus
Results by SROM Results by MROM Results by Homogenization
6
4
2
0 0.0
0.2
Volume Fraction of Fiber, Vf
Volume Fraction of Fiber, Vf
8
(26)
Among the equivalent mechanical properties obtained from the homogenization method, the selected response variables were E11,eq, E22,eq, G12,eq, G23,eq, n12,eq and n23,eq, respectively. The results
Normalized Transverse Modulus, E22,eq/Em
Table 2 Equivalent material properties of unit cell.
0.2
0.4
0.6
Volume Fraction of Fiber, Vf
(c) Equivalent shear modulus
0.8
Equivalent Poisson's ratio (ν12,eq) 1.5
1.0
0.5
0.0 0.0
Results by SROM Results by Homogenization
0.2
0.4
0.6
Volume Fraction of Fiber, Vf
(d) Equivalent Poisson's ratio
Fig. 6. Comparison of results by SROM, MROM and homogenization technique.
0.8
224
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
Table 3 Design variables for Monte Carlo simulation. Design variable
Type
Mode
Distribution
Mean
CoV
Ef [GPa] Em [GPa]
Real Real Real Real
Continuous Continuous Continuous Continuous
Normal Normal Normal Normal
72.4 2.45 0.25 0.35
0.1 0.1 0.1 0.1
20
20
Equivalent transverse modulus (E22,eq)
V =0.265
V =0.465
V =0.365
V =0.565
15
V =0.265
Number of Occurence
Number of Occurence
Equivalent longitudinal modulus (E11,eq)
V =0.465 V =0.265 V =0.365
V =0.565
10
5
0 15
20
25
30
35
40
45
50
V =0.365
15
V =0.365
V =0.565
5
0 2
4
6
8
10
12
14
Equivalent Transverse Modulus (E22,eq)
(b) Equivalent transverse modulus 25
25
Equivalent shear modulus (G23,eq)
Equivalent shear modulus (G12,eq) V =0.265
V =0.465
V =0.365
V =0.565
V =0.365
20
Number of Occurence
20
Number of Occurence
V =0.565 V =0.465
10
55
(a) Equivalent longitudinal modulus
V =0.465
V =0.265
Equivalent Longitudinal Modulus (E11,eq)
V =0.565
V =0.265 V =0.465
15
10
5
0 1.0
1.5
2.0
2.5
3.0
3.5
V =0.265
V =0.465
V =0.365
V =0.565 V =0.465
15
V =0.565
V =0.365 V =0.265
10
5
0 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75
4.0
Equivalent Shear Modulus (G23,eq)
Equivalent Shear Modulus (G12,eq)
(c) Equivalent shear modulus (G12,eq) (d) Equivalent shear modulus (G23,eq) 25
Equivalent Poisson's ratio (ν12,eq)
20
15
V =0.265
V =0.465
V =0.365
V =0.565
V =0.565
V =0.465
V =0.265 V =0.365
10
5
0 0.20
0.25
Equivalent Poisson's ratio (ν23,eq)
20
0.30
0.35 0.40 Equivalent Poisson's Ratio (ν12,eq)
(e) Equivalent Poisson's raio (ν12,eq )
Number of Occurence
25
Number of Occurence
nf nm
were derived from four different volume fractions: 0.265, 0.365, 0.465, and 0.565. The results of the stochastic analyses calculated for various fiber volume fractions (Vf) are shown in Fig. 7. As can be seen from the figure, the best fitted distribution to the response variables is the normal distribution, which is the same distribution as the design variables. The standard deviation of E11,eq is larger than those for the other equivalent properties; next, the sequence went E22,eq, G12,eq and G23,eq. In particular, the standard deviations of these
V =0.265
V =0.465
V =0.365
V =0.565
V =0.465
15
V =0.365 V =0.265
V =0.565
10
5
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
Equivalent Poisson's Ratio (ν23,eq)
(f) Equivalent Poisson's ratio (ν 23,eq)
Fig. 7. Distribution of material properties calculated by homogenization technique and Monte-Carlo simulation.
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
response variables show a trend where they all increased as the fiber volume fraction increased from 0.265 to 0.565. The standard deviations of n12,eq and n23,eq were about 0.2e0.3 and 0.5e0.6, depending on the change in the volume fraction, respectively. To summarize, E11,eq, E22,eq, G12,eq and G23,eq increased as the fiber volume fraction increased, but n12,eq and n23,eq showed the reverse trend. Also, while E11,eq for Vf ¼ 0.265 is 20.8 GPa in a deterministic value, it could be expressed as 20.8 GPa þ4.22/e4.38 GPa with a reliability of 90% in a probabilistic way. From the above results, it can be stated that the proposed homogenization based Monte Carlo simulation leads to a more reliable results for design of composite structure. From the results of a linear correlation analysis that was conducted to investigate the correlation between the design and response variables, the results corresponding to a volume fraction of 0.365 are shown in Fig. 8. Correlation is a bivariate measure of association (strength) of the relationship between two variables. It varies from 0 (random relationship) to 1 (perfect linear relationship) or 1 (perfect negative linear relationship). HyperStudy uses the pairing procedure developed by Iman and Conover to induce specified correlations between random variables [21]. Fig. 8 shows that Ef has a strong correlation with E11,eq and has a weak correlation with E22,eq, G12,eq, G23,eq, n12,eq and n23,eq. That is, it is found that E11,eq is mainly governed by the properties of the fiber. In addition, Em is shown to correlate with E22,eq, G12,eq and G23,eq, but not to
Fig. 8. Linear correlation analysis results for Vf ¼ 0.365.
40
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 50
60
70
80
90
Equivalent Longitudinal Modulus, E11,eq
Equivalent Longitudinal Modulus, E11,eq
40
100
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 1.8
2.0
2.2
2.4
2.6
2.8
3.0
3.2
Elastic Modulus of Matrix, Em[GPa]
Elastic Modulus of Fiber, EF[GPa]
(b)Effect of Em
(a) Effect of E f 40
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32
Poisson's Ratio of Fiber, νf
(c) Effect of
f
Equivalent Longitudinal Modulus, E11,eq
40
Equivalent Longitudinal Modulus, E11,eq
225
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 0.24
0.28
0.32
0.36
0.40
Poisson's Ratio of Matrix, νm
(d)Effect of m
Fig. 9. Effect plot of design variables on E11,eq.
0.44
226
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
have any correlation with E11,eq, n12,eq and n23,eq, signifying that Em plays a dominant role for E22,eq, G12,eq and G23,eq. It was found that nf does not correlate to any of the mechanical properties, while nm is found to correlate to n12,eq and n23,eq. Fig. 9 shows the effect of the design variables on E11,eq through a regression analysis. In Fig. 9, the dotted lines indicate the minimal, mean, and maximum values of the equivalent properties. As can be seen in Fig. 9, the design variable having the greatest effect on E11,eq is Ef, followed by Em. On the other hand, nm and nf have no substantial effect. 4. Conclusions In the present study, we used a homogenization method and Monte Carlo simulation to probabilistically estimate the equivalent properties of glass fiber-reinforced composite materials and to assess how these are affected by the constituents of the composite materials. The obtained conclusions are following. 1. The homogenization method was implemented to estimate the equivalent properties of composite material, which depend on a microstructure and behavior of constituents. The homogenized mechanical properties conform well to the results from the rule of mixtures ranging from Vf ¼ 0.1 to Vf ¼ 0.7 of volume fraction of fiber. 2. To assess the uncertainty of the homogenized equivalent properties of the composite materials, Monte-Carlo simulation coupled with the homogenization method was carried out according to the variability in the basic properties of the constituents. The equivalent properties follow the normal distribution and are simultaneously affected by the basic properties of the constituents. 3. From the sensitivity analysis on the basic properties of the constituents, it was found that the equivalent properties of composite materials are mainly governed by the modulus of fiber and matrix, not their Poisson’s ratio. Acknowledgment This work was financially supported by the Human Resources Development of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea government Ministry of Knowledge Economy (No. 20124010203240) and by
Basic Science Research Program through the NRF funded by the Ministry of Education, Science and Technology (2011-0007012).
References [1] Kang KW, Lim DM, Kim JK. Probabilistic analysis for the fatigue life of carbon/ epoxy laminates. Compos Struct 2008;85:258e64. [2] Hashin Z, Rosen BW. The elastic moduli of fiber-reinforced materials. J Appl Mech 1964;21:223e32. [3] Tanaka S, Ichikawa M, Akita S. A probabilistic investigation of fatigue life and cumulative cycle ratio. Eng Fract Mech 1984;20:501e13. [4] Kam TY, Tsai SY, Chu KH. Fatigue reliability analysis of composite laminates under spectrum stress. Intern J Solids Struct 1997;34:1441e61. [5] Hassani B, Hinton E. A review of homogenization and topology optimization II e analytical and numerical solution of homogenization equations. Comput Struct 1998;69:719e38. [6] Orlik J. Asymptotic homogenization algorithm for reinforced metal-matrix elasto-plastic composites. Compos Struct 2010;92:1581e90. [7] Ji-Wei D, Miao-Lin F. Asymptotic expansion homogenization for simulating progressive damage of 3D braided composites. Compos Struct 2010;92: 873e82. [8] Li HX. Limit analysis of composite materials with anisotropic microstructures: a homogenization approach. Mech Mater 2011;43:574e85. [9] Kaminski M, Kleiber M. Perturbation based stochastic finite element method for homogenization of two-phase elastic composites. Comput Struct 2000;78: 811e26. [10] Sakata S, Ashida F, Kojima T, Zako M. Three-dimensional stochastic analysis using a perturbation based homogenization method for elastic properties of composite material considering microscopic uncertainty. Intern J Solids Struct 2008;45:894e907. [11] Kaminski M, Lesniak M. Homogenization of metallic fiber-reinforced composites under stochastic ageing. Compos Struct 2012;94:386e93. [12] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Methods Appl Mech Eng 1990;83:143e98. [13] Lee SP, Kim JH, Kim DJ, Ha SK. Numerical analysis of dielectric ceramic and void containing Ni electrode layers of multilayer ceramic capacitors. J Compos Mater 2012. http://dx.doi.org/10.1177/0021998312449882 [under publication]. [14] Hsueh CH, Ferber MK. Apparent coefficient of thermal expansion and residual stresses in multilayer capacitors. Composites Part A 2002;33:1115e21. [15] Aghdam MM, Smith DJ, Pavier MJ. Finite element micromechanical modeling of yield and collapse behavior of metal matrix composites. J Mech Phys Solids 2000;48:499e528. [16] Sun CT, Vaidya RS. Prediction of composite properties from a representative volume element. Compos Sci Technol 1996;56:171e9. [17] Hyer MW, Waas AM. Micromechanics of linear elastic continuous fiber composites. New York: WCB/McGraw-Hill; 1999. [18] Halphin JC. Effect of environmental factors on composite materials. AFML-TR67e423; 1969. [19] ANSYSANSYS APDL programmer’s guide ver. 10.0. PA, USA: ANSYS. Inc; 2005. [20] Kang KW, Koh SK, Kim DK, Kim KJ. Assessment of the statistical distribution of flexural strength of woven-fabric laminates with impact-induced damage. Compos Struct 2009;90:60e6. [21] HyperStudyUser’s guide ver. 11.0. Korea: Altair Inc; 2012.
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Probabilistic analysis for mechanical properties of glass/epoxy composites using homogenization method and Monte Carlo simulation Seung-Pyo Lee a, Ji-Won Jin b, Ki-Weon Kang b, * a b
R&D Center, ILJIN Global, 128-5, Samsung-dong, Kangnam-go, Seoul 135-875, Republic of Korea Department of Mechanical Engineering, Kunsan National University, Kunsan, Jeonbuk 573-701, Republic of Korea
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 February 2013 Accepted 12 September 2013 Available online 1 October 2013
In this paper, we used a homogenization method to estimate the equivalent properties of composite material, which depend on a microstructure and behavior of constituents (fiber and matrix). The estimated results conform well to the results of the rule of mixtures, which is a conventional method for estimating equivalent properties of composite materials. In addition, to assess the uncertainty of the equivalent properties of the composite materials according to the variability in the basic properties of the constituents, a probabilistic analysis using a homogenization-based Monte Carlo simulation was carried out. In this way, the variation of basic properties of the constituents was identified to have a significant effect on the uncertainty of the equivalent properties of composite materials. Moreover, a sensitivity of the properties of the constituents on the equivalent properties was assessed through a correlation and regression analysis. 2013 Elsevier Ltd. All rights reserved.
Keywords: Glass/epoxy composites Homogenization Mechanical properties Microstructure Monte Carlo simulation
1. Introduction Fiber-reinforced composites have superior mechanical properties such as specific stiffness and specific strength, and they have enabled us a stiffness design through changes in the stacking sequence and fiber volume fraction. Then, such composites have been used extensively in aerospace, automotive, and wind energy conversion industries [1]. Composite materials have a structure wherein reinforcing fiber and a matrix, each with unique characteristics, are physically coupled on a microscopic scale. The complicated microstructure and the behaviors of the constituents lead to the inherent characteristics of composite materials, such as strong heterogeneity and anisotropy [2,3]. And, because the basic properties of their constituents must necessarily be dispersed, there is a considerable uncertainty in the mechanical properties of composite materials [4]. Therefore, it is very important to assess the mechanical properties in a manner that takes into consideration the behavior of the microstructure and constituents of a composite material, as well as the dispersion of their basic properties. The homogenization method is one of the most useful methods to assess the mechanical properties of composite materials, which
* Corresponding author. Tel.: þ82 63 469 4872; fax: þ82 63 469 4727. E-mail address: [email protected] (K.-W. Kang). 0960-1481/$ e see front matter 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.renene.2013.09.012
takes into consideration the behavior of the constituents and microstructure [5]. The homogenization method is a procedure whereby a heterogeneous media is converted into an equivalent material model that is energetically equivalent to the heterogeneous media. On the basis of the fundamental concepts of the homogenization method, Orlik [6] applied the theory of the two-scale convergence to homogenization of initial flow stresses and hardening constants in exponential hardening laws for elasto-plastic composites with a periodic microstructure. Ji-Wei et al. [7] applied the asymptotic expansion homogenization (AEH) method to simulate the non-linear behavior of the cumulative damage to various types of unit cell. This allowed them to research the tensile strength of 3D braided composites from a microscopic point of view. Li [8] also developed a non-linear algorithm together with the finite element analysis and homogenization method to implement kinematic limit analysis for a microstructure and the macroscopic strength of a composite material with anisotropic constituents. And others have also carried out researches that coupled a stochastic approach to the homogenization method, which made it possible to consider the behavior of the constituents and microstructures together with the uncertainty in the basic properties of the constituents. Kaminski et al. [9] have proposed the stochastic second order and second moment analysis to analyze the first two probabilistic moments of the effective elasticity tensor components for homogenization of the two-phase periodic composite structure.
220
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
ckl i
Nomenclature
sεij xεj
fi uεi nj Fi eij Dεijkl ε ~ ð1Þ u i
stress generalized coordinate system for the macroscopic scale body force displacement normal vector traction strain constitutive matrix scale parameter for asymptotic series approximation arbitrary additive constant in yi
Sakata et al. [10] have analyzed the influence of microscopic uncertainty on a homogenized macroscopic elastic properties of composite materials using higher order perturbation-based homogenization method. Also, Kaminski et al. [11] have proposed the generalized stochastic perturbation method to analyze stochastic aging process of the metallic fiber reinforced composites. They have used the computational homogenization procedure based on finite element analysis. Despite offering a useful depiction of the uncertainty of the macroscopic mechanical properties due to the dispersion of the microscopic properties, these studies are of limited applicability and the difficulties in their mathematical formulations. The present study aims to probabilistically estimate the equivalent properties of glass fiber-reinforced composite materials, by using a homogenization method and Monte Carlo simulation. This first made it possible to estimate the homogenized equivalent properties, which take into consideration the microstructure and the behavior of constituents. The results were then assessed in comparison to the properties estimated from an ROM (rule of mixtures). The homogenization method-based Monte Carlo simulation also allowed us to estimate the variability of the equivalent properties of composite materials, caused by the dispersion in the properties of the constituent, as well as to assess how this variability is impacted by properties of the constituents. 2. Theoretical background and analysis procedure 2.1. Homogenization method To compare the dimensions of the composite structure composed of periodic heterogeneous unit cell such as wind turbine composite blade shown in Fig. 1, let us assume that the dimensions of the unit cell are relatively very small, with a value of ε. Given an appropriate load and boundary conditions for such a structure, there will be rapid changes in the stress and strain as a result of the
DH ijmn yi E11,eq E22,eq G12,eq Ef Em Gf Gm
nf nm
characteristic function equivalent material property tensor second coordinate system for the microscopic scale longitudinal elastic modulus of unit cell transverse elastic modulus of unit cell shear modulus of unit cell elastic modulus of fiber elastic modulus of matrix shear modulus of fiber shear modulus of matrix Poisson’s ratio of fiber Poisson’s ratio of matrix
repetition of the unit cell. In other words, the stress and strain undergo very rapid changes in the very small circumference, ε, of a point, x, in the structure. To portray this effectively, the macroscopic and microscopic levels can be defined in terms of x and y ¼ x/ε, respectively [12]. In conventional elasticity, a body subject to body force fi should satisfy the following equilibrium equation [12e14]
vsεij vxεj
þ fi ¼ 0 in U:
(1)
And the prescribed displacement uεi on v1 U and traction Fi on the boundary v2 U can be stated as
uεi ¼ 0 on v1 U;
(2)
sεij nj ¼ Fi on v2 U:
(3)
The strainedisplacement and the stressestrain relationships can be written as ε 1 vuεi vuj eij ðu Þ ¼ þ 2 vxεj vxεi
!
ε
(4)
sεij ¼ Dεijkl ekl ðuε Þ
(5)
The displacements can be approximated by an asymptotic series representation in ε given by Guedes and Kikuchi [12] ð0Þ
ð1Þ
uεi ðxÞ ¼ ui ðx; yÞ þ εui ðx; yÞ þ /:
(6)
Substituting eq. (6) into eq. (4) and then substituting the result into eq. (5) enables us to obtain an equation that is dependent on powers of ε. Each term of the same power of ε is equal to zero to ensure that the asymptotic series approximation is valid as ε approaches zero. Hence, ð0Þ
vuk v D ¼ 0 vyj ijkl vyl
(7)
ð0Þ
ð0Þ
ð1Þ
vuk vuk vu v v D þ D þ k vxj ijkl vyl vyj ijkl vxl vyl
ð0Þ
Fig. 1. Wind turbine composite blade.
ð1Þ
vuk vu v D þ k vxj ijkl vxl vyl
!
! ¼ 0;
ð1Þ
(8)
ð2Þ
vuk vu v D þ k þ vyj ijkl vxl vyl
! þ fi ¼ 0
(9)
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
ð1Þ
ui
¼ ckl i
221
ð0Þ
vuk ~ ið1Þ ðxÞ; þu vxl
(10)
ð1Þ
~ i ðxÞ is constant and ckl in which u i is a characteristic function that is the solution of the auxiliary equation given by
Z Dijkl Y
vcmn k vyi dy ¼ vyl vyj
Z Dijmn Y
vyi dy cyi ˛VY : vyj
(11)
ð1Þ
Let us assume that the solution of ui takes the form as following [12,13]. Substituting eq. (10) into eq. (9) results in eq. (12). ð0Þ
v H vum D þ fi ¼ 0 vxj ijmn vxn
(12)
The equivalent material property tensor DH ijmn can be written as [12,13]
DH ijmn
1 ¼ jYj
# Z " vckl p Dijkl Dijpq dy vyq
(13)
Y
2.2. Analysis procedure In composite materials, the distribution of fiber and matrix is likely to be random across the cross-section as shown in Fig. 2. Since it is difficult to model random fiber arrangement, many micromechanical models [15,16] assume a periodic arrangement of constituents for which unit cell can be isolated. The periodic fiber and matrix arrangements are the square and the hexagonal arrays. The unit cell has the identical elastic properties and volume fraction of fiber as the composite materials. Thus, a micromechanics-based finite element analysis for a unit cell can lead to the mechanical properties of a whole structure, which comprises the repeated rectangular unit cell in Fig. 3. In this paper, accordingly, the unit cell is defined as the square-packed array in Fig. 4 and the homogenization method is applied to the unit cell to estimate the equivalent mechanical properties of the E-glass/epoxy composite.
Fig. 3. Idealized array cells of unidirectional composites; square-packed array.
In order to calculate the equivalent mechanical properties of the unit cell, the present paper utilizes an ANSYS SOLID45 element (an eight-node structural solid element) with eight nodes in one element and three degrees of freedom in one node. The total number of nodes herein is 386 and the total number of elements is 180 (refer to Fig. 5). Periodic boundary conditions are applied as eqs. (14)e(16) to the boundaries of the unit cell in order to apply the homogenization method. The mutually facing sides of the unit cell have the same behavior, on the basis of these boundary conditions [13]. jk cjk ðy1 ; y2 ; y3 Þ ¼ ci ðy1 þ Y1 ; y2 ; y3 Þ i
(14)
cjk ðy1 ; y2 ; y3 Þ ¼ cjk ðy1 ; y2 þ Y2 ; y3 Þ i i
(15)
jk cjk ðy1 ; y2 ; y3 Þ ¼ ci ðy1 ; y2 ; y3 þ Y3 Þ i
(16)
The characteristic function ckl i given in eq. (11) is the matrix. Each column of this matrix matches a column of a constitutive matrix, Dijkl. That is, eq. (11) has the columns comprising the
Fig. 2. Microscopic configuration of multiaxial glass/epoxy composites.
Fig. 4. Geometric models of square unit cells for glass/epoxy composites.
222
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
h¼
g1 gþ2
(22)
In addition, a symbol g is a mechanical property ratio between the fiber and matrix and is represented by the following equation:
g¼
Pf Pm
(23)
And 2 is the variable for the partitioning effect of the fiber and matrix. When 2 ¼ 0 in eq. (22), eq. (21) is expressed as follows:
1 V Vm ¼ fþ P Pf Pm
(24)
When 2 ¼ N and h ¼ 0, eq. (21) is expressed as follows:
P ¼ Pm 1 Vf þ Pf Vf
(25)
3. Results and discussion 3.1. Equivalent mechanical properties using homogenization method
Fig. 5. Finite element model of square unit cell for glass/epoxy composites.
solution of a constitutive matrix, in a three-dimensional elasticity problem. The solution for these six columns can be obtained from the six columns in eq. (11). 2.3. Rule of mixture (ROM) The mathematical equations for the simple rule-of-mixtures (SROM), which is another method for the equivalent mechanical properties, are as follows.
E11;eq ¼ Ef Vf þ Em 1 Vf
n12;eq ¼ nf Vf þ nm 1 Vf
(17) (18)
1 V 1 Vf ¼ fþ E22;eq Ef Em
(19)
1 V 1 Vf ¼ fþ G12;eq Gf Gm
(20)
where the subscripts f and m represent the fibers and matrix, respectively, and Vf signifies the volume fraction of the fibers. Hyer et al. [17] determined the mechanical properties using a structural analysis of the unit cell and a simple rule-of-mixtures method. In a comparison of the results, the values for E22,eq and G12,eq showed a large difference. Thus, in order to correct for this, the modified rule of mixtures (MROM) proposed by Halphin [18] was introduced. The equation in the modified rule of mixtures is as follows.
P 1 þ 2hVf ¼ Pm 1 hVf
The equivalent mechanical properties of the glass/epoxy composites can be calculated by means of the homogenization method mentioned in Section 2.1 and Section 2.2 for the unit cell, which comprise a matrix and fiber. For this, the characteristic function in eq. (11) was first solved, the result of which was substituted into eq. (13) to calculate the equivalent mechanical properties. The geometry of the unit cell that was used to calculate the equivalent mechanical properties is shown in Fig. 4, and the boundary conditions were as shown in eqs. (14)e(16). Respective CE commands in the ANSYS, were applied to mutually facing sides to express the boundary conditions in eqs. (14)e(16). And the basic mechanical properties of the E-glass fiber and epoxy resin were used, as indicated in Table 1. The load conditions used to solve the characteristic function were based on the right-side operand in eq. (11). Because the above equations include an element strain matrix and material properties matrix, it is not possible to express these using a general method for load input. Thus, in the present study, this was resolved by taking the load terms as the initial stress terms and utilizing the user subroutine USTRESS in ANSYS to represent the initial stress terms as per eq. (11). In addition, APDL [19] in ANSYS was used to calculate the equivalent mechanical properties from eq. (13). The equivalent mechanical properties calculated using the above procedure of the homogenization method is summarized in Table 2. To verify the homogenization method, the mechanical properties were determined using the rule of mixtures and compared with the results by the homogenization method. The volume fraction of the fiber was varied from 0.1 to 0.7. As shown in Fig. 6, E11,eq and n12,eq are consistent with the SROM results, and E22,eq and G12,eq are consistent with the MROM results presented by Hyer [17]. Herein, the coefficient z can be calculated using a finite element analysis of the unit cell. In the present study, when z ¼ 1.8 and z ¼ 1.2, the equivalent properties E22,eq and G12,eq were found to be consistent with the homogenized results. Based on the above
Table 1 Material properties of constituents.
(21)
In eq. (21), P signifies E22,eq or G12,eq, which are the mechanical properties of the composites, and Pm is a mechanical property of the matrix. Herein, h is defined as follows:
Elastic modulus, E [GPa] Shear modulus, G [GPa] Poisson’s ratio, n
Fiber; E-glass
Matrix; epoxy
72.40 28.96 0.25
2.45 0.91 0.35
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
223
Vf
E11,eq [GPa]
E22,eq [GPa]
G12,eq [GPa]
G12,eq [GPa]
n12,eq
n23,eq
materials according to the variation of the basic mechanical properties. The homogenization based Monte Carlo simulation has three main steps:
0.1 0.2 0.3 0.4 0.5 0.6 0.7
9.369 16.286 23.203 30.119 37.035 43.951 50.868
3.170 3.889 4.855 6.214 8.196 11.278 16.838
1.095 1.326 1.617 2.001 2.541 3.395 5.090
1.051 1.203 1.377 1.595 1.897 2.382 3.406
0.338 0.326 0.315 0.304 0.294 0.284 0.273
0.460 0.453 0.423 0.380 0.326 0.266 0.212
1. Generate a randomly distributed set of basic properties of the constituents. 2. Perform calculations based on the homogenization method using the generated set of basic properties of the constituents. 3. Determine probability from a large number of repetitions.
results, it can be stated that the homogenization method can lead to the equivalent mechanical properties using the basic mechanical properties of the constituents. 3.2. Homogenization based Monte Carlo simulation
Normalized Longitudinal Modulus, E11,eq/Ef
The above-presented homogenization method can estimate the equivalent mechanical properties of composites using the basic properties of the constituents. However, since the mechanical properties of constituents are random in nature, one can observe great differences in the macroscopic mechanical properties of composites among the samples in the same condition. The mechanical properties of composites should, therefore, be treated in a probabilistic, rather than a deterministic way [1,20]. In this study, a homogenization based Monte Carlo simulation was applied to assess the uncertainty of the macroscopic mechanical properties, i.e., the equivalent properties of the composite 1.0
Equivalent longitudinal modulus (E11,eq) 0.8
Results by SROM Results by Homogenization
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
For the sampling of basic properties and assessment of uncertainty of the equivalent properties, the commercial code HyperStudy [21] was used. The sampling method used here was the Hamersley method [21], which is one type of quasi-Monte Carlo method and is proven to be effective for a small number of samples; the sample number was 100. The selected design variables (basic properties of the constituents) were the elastic modulus (Ef and Em) and Poisson’s ratio (nf and nm) of the fiber and matrix. Table 3 shows the design variables for the homogenization based Monte Carlo simulation. Herein, the term of CoV is defined as the coefficient of variation and can be calculated from eq. (26).
CoV ¼
standard deviation mean
10
Equivalent transverse modulus (E22,eq) Results by SROM Results by MROM Results by Homogenization
8
6
4
2
0 0.0
0.6
0.8
2.0
10
Equivalent shear modulus (G12,eq)
Normalized Poisson's Ratio, ν12,eq/ν12,f
Normalized Shear Modulus, G12,eq/Gm
0.4
(b) Equivalent transverse modulus
(a) Equivalent longitudinal modulus
Results by SROM Results by MROM Results by Homogenization
6
4
2
0 0.0
0.2
Volume Fraction of Fiber, Vf
Volume Fraction of Fiber, Vf
8
(26)
Among the equivalent mechanical properties obtained from the homogenization method, the selected response variables were E11,eq, E22,eq, G12,eq, G23,eq, n12,eq and n23,eq, respectively. The results
Normalized Transverse Modulus, E22,eq/Em
Table 2 Equivalent material properties of unit cell.
0.2
0.4
0.6
Volume Fraction of Fiber, Vf
(c) Equivalent shear modulus
0.8
Equivalent Poisson's ratio (ν12,eq) 1.5
1.0
0.5
0.0 0.0
Results by SROM Results by Homogenization
0.2
0.4
0.6
Volume Fraction of Fiber, Vf
(d) Equivalent Poisson's ratio
Fig. 6. Comparison of results by SROM, MROM and homogenization technique.
0.8
224
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
Table 3 Design variables for Monte Carlo simulation. Design variable
Type
Mode
Distribution
Mean
CoV
Ef [GPa] Em [GPa]
Real Real Real Real
Continuous Continuous Continuous Continuous
Normal Normal Normal Normal
72.4 2.45 0.25 0.35
0.1 0.1 0.1 0.1
20
20
Equivalent transverse modulus (E22,eq)
V =0.265
V =0.465
V =0.365
V =0.565
15
V =0.265
Number of Occurence
Number of Occurence
Equivalent longitudinal modulus (E11,eq)
V =0.465 V =0.265 V =0.365
V =0.565
10
5
0 15
20
25
30
35
40
45
50
V =0.365
15
V =0.365
V =0.565
5
0 2
4
6
8
10
12
14
Equivalent Transverse Modulus (E22,eq)
(b) Equivalent transverse modulus 25
25
Equivalent shear modulus (G23,eq)
Equivalent shear modulus (G12,eq) V =0.265
V =0.465
V =0.365
V =0.565
V =0.365
20
Number of Occurence
20
Number of Occurence
V =0.565 V =0.465
10
55
(a) Equivalent longitudinal modulus
V =0.465
V =0.265
Equivalent Longitudinal Modulus (E11,eq)
V =0.565
V =0.265 V =0.465
15
10
5
0 1.0
1.5
2.0
2.5
3.0
3.5
V =0.265
V =0.465
V =0.365
V =0.565 V =0.465
15
V =0.565
V =0.365 V =0.265
10
5
0 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75
4.0
Equivalent Shear Modulus (G23,eq)
Equivalent Shear Modulus (G12,eq)
(c) Equivalent shear modulus (G12,eq) (d) Equivalent shear modulus (G23,eq) 25
Equivalent Poisson's ratio (ν12,eq)
20
15
V =0.265
V =0.465
V =0.365
V =0.565
V =0.565
V =0.465
V =0.265 V =0.365
10
5
0 0.20
0.25
Equivalent Poisson's ratio (ν23,eq)
20
0.30
0.35 0.40 Equivalent Poisson's Ratio (ν12,eq)
(e) Equivalent Poisson's raio (ν12,eq )
Number of Occurence
25
Number of Occurence
nf nm
were derived from four different volume fractions: 0.265, 0.365, 0.465, and 0.565. The results of the stochastic analyses calculated for various fiber volume fractions (Vf) are shown in Fig. 7. As can be seen from the figure, the best fitted distribution to the response variables is the normal distribution, which is the same distribution as the design variables. The standard deviation of E11,eq is larger than those for the other equivalent properties; next, the sequence went E22,eq, G12,eq and G23,eq. In particular, the standard deviations of these
V =0.265
V =0.465
V =0.365
V =0.565
V =0.465
15
V =0.365 V =0.265
V =0.565
10
5
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
Equivalent Poisson's Ratio (ν23,eq)
(f) Equivalent Poisson's ratio (ν 23,eq)
Fig. 7. Distribution of material properties calculated by homogenization technique and Monte-Carlo simulation.
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
response variables show a trend where they all increased as the fiber volume fraction increased from 0.265 to 0.565. The standard deviations of n12,eq and n23,eq were about 0.2e0.3 and 0.5e0.6, depending on the change in the volume fraction, respectively. To summarize, E11,eq, E22,eq, G12,eq and G23,eq increased as the fiber volume fraction increased, but n12,eq and n23,eq showed the reverse trend. Also, while E11,eq for Vf ¼ 0.265 is 20.8 GPa in a deterministic value, it could be expressed as 20.8 GPa þ4.22/e4.38 GPa with a reliability of 90% in a probabilistic way. From the above results, it can be stated that the proposed homogenization based Monte Carlo simulation leads to a more reliable results for design of composite structure. From the results of a linear correlation analysis that was conducted to investigate the correlation between the design and response variables, the results corresponding to a volume fraction of 0.365 are shown in Fig. 8. Correlation is a bivariate measure of association (strength) of the relationship between two variables. It varies from 0 (random relationship) to 1 (perfect linear relationship) or 1 (perfect negative linear relationship). HyperStudy uses the pairing procedure developed by Iman and Conover to induce specified correlations between random variables [21]. Fig. 8 shows that Ef has a strong correlation with E11,eq and has a weak correlation with E22,eq, G12,eq, G23,eq, n12,eq and n23,eq. That is, it is found that E11,eq is mainly governed by the properties of the fiber. In addition, Em is shown to correlate with E22,eq, G12,eq and G23,eq, but not to
Fig. 8. Linear correlation analysis results for Vf ¼ 0.365.
40
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 50
60
70
80
90
Equivalent Longitudinal Modulus, E11,eq
Equivalent Longitudinal Modulus, E11,eq
40
100
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 1.8
2.0
2.2
2.4
2.6
2.8
3.0
3.2
Elastic Modulus of Matrix, Em[GPa]
Elastic Modulus of Fiber, EF[GPa]
(b)Effect of Em
(a) Effect of E f 40
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32
Poisson's Ratio of Fiber, νf
(c) Effect of
f
Equivalent Longitudinal Modulus, E11,eq
40
Equivalent Longitudinal Modulus, E11,eq
225
Equivalent longitudinal modulus (E11,eq) Estimated Results
35
30
25
20 0.24
0.28
0.32
0.36
0.40
Poisson's Ratio of Matrix, νm
(d)Effect of m
Fig. 9. Effect plot of design variables on E11,eq.
0.44
226
S.-P. Lee et al. / Renewable Energy 65 (2014) 219e226
have any correlation with E11,eq, n12,eq and n23,eq, signifying that Em plays a dominant role for E22,eq, G12,eq and G23,eq. It was found that nf does not correlate to any of the mechanical properties, while nm is found to correlate to n12,eq and n23,eq. Fig. 9 shows the effect of the design variables on E11,eq through a regression analysis. In Fig. 9, the dotted lines indicate the minimal, mean, and maximum values of the equivalent properties. As can be seen in Fig. 9, the design variable having the greatest effect on E11,eq is Ef, followed by Em. On the other hand, nm and nf have no substantial effect. 4. Conclusions In the present study, we used a homogenization method and Monte Carlo simulation to probabilistically estimate the equivalent properties of glass fiber-reinforced composite materials and to assess how these are affected by the constituents of the composite materials. The obtained conclusions are following. 1. The homogenization method was implemented to estimate the equivalent properties of composite material, which depend on a microstructure and behavior of constituents. The homogenized mechanical properties conform well to the results from the rule of mixtures ranging from Vf ¼ 0.1 to Vf ¼ 0.7 of volume fraction of fiber. 2. To assess the uncertainty of the homogenized equivalent properties of the composite materials, Monte-Carlo simulation coupled with the homogenization method was carried out according to the variability in the basic properties of the constituents. The equivalent properties follow the normal distribution and are simultaneously affected by the basic properties of the constituents. 3. From the sensitivity analysis on the basic properties of the constituents, it was found that the equivalent properties of composite materials are mainly governed by the modulus of fiber and matrix, not their Poisson’s ratio. Acknowledgment This work was financially supported by the Human Resources Development of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea government Ministry of Knowledge Economy (No. 20124010203240) and by
Basic Science Research Program through the NRF funded by the Ministry of Education, Science and Technology (2011-0007012).
References [1] Kang KW, Lim DM, Kim JK. Probabilistic analysis for the fatigue life of carbon/ epoxy laminates. Compos Struct 2008;85:258e64. [2] Hashin Z, Rosen BW. The elastic moduli of fiber-reinforced materials. J Appl Mech 1964;21:223e32. [3] Tanaka S, Ichikawa M, Akita S. A probabilistic investigation of fatigue life and cumulative cycle ratio. Eng Fract Mech 1984;20:501e13. [4] Kam TY, Tsai SY, Chu KH. Fatigue reliability analysis of composite laminates under spectrum stress. Intern J Solids Struct 1997;34:1441e61. [5] Hassani B, Hinton E. A review of homogenization and topology optimization II e analytical and numerical solution of homogenization equations. Comput Struct 1998;69:719e38. [6] Orlik J. Asymptotic homogenization algorithm for reinforced metal-matrix elasto-plastic composites. Compos Struct 2010;92:1581e90. [7] Ji-Wei D, Miao-Lin F. Asymptotic expansion homogenization for simulating progressive damage of 3D braided composites. Compos Struct 2010;92: 873e82. [8] Li HX. Limit analysis of composite materials with anisotropic microstructures: a homogenization approach. Mech Mater 2011;43:574e85. [9] Kaminski M, Kleiber M. Perturbation based stochastic finite element method for homogenization of two-phase elastic composites. Comput Struct 2000;78: 811e26. [10] Sakata S, Ashida F, Kojima T, Zako M. Three-dimensional stochastic analysis using a perturbation based homogenization method for elastic properties of composite material considering microscopic uncertainty. Intern J Solids Struct 2008;45:894e907. [11] Kaminski M, Lesniak M. Homogenization of metallic fiber-reinforced composites under stochastic ageing. Compos Struct 2012;94:386e93. [12] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Methods Appl Mech Eng 1990;83:143e98. [13] Lee SP, Kim JH, Kim DJ, Ha SK. Numerical analysis of dielectric ceramic and void containing Ni electrode layers of multilayer ceramic capacitors. J Compos Mater 2012. http://dx.doi.org/10.1177/0021998312449882 [under publication]. [14] Hsueh CH, Ferber MK. Apparent coefficient of thermal expansion and residual stresses in multilayer capacitors. Composites Part A 2002;33:1115e21. [15] Aghdam MM, Smith DJ, Pavier MJ. Finite element micromechanical modeling of yield and collapse behavior of metal matrix composites. J Mech Phys Solids 2000;48:499e528. [16] Sun CT, Vaidya RS. Prediction of composite properties from a representative volume element. Compos Sci Technol 1996;56:171e9. [17] Hyer MW, Waas AM. Micromechanics of linear elastic continuous fiber composites. New York: WCB/McGraw-Hill; 1999. [18] Halphin JC. Effect of environmental factors on composite materials. AFML-TR67e423; 1969. [19] ANSYSANSYS APDL programmer’s guide ver. 10.0. PA, USA: ANSYS. Inc; 2005. [20] Kang KW, Koh SK, Kim DK, Kim KJ. Assessment of the statistical distribution of flexural strength of woven-fabric laminates with impact-induced damage. Compos Struct 2009;90:60e6. [21] HyperStudyUser’s guide ver. 11.0. Korea: Altair Inc; 2012.