- Email: [email protected]
An approach for predicting the elastic modulus of heterogeneous materials
Materials and Design 30 (2009) 2938–2945
Contents lists available at ScienceDirect
Materials and Design journal homepage: www.elsevier.com/locate/matdes
An approach for predicting the elastic modulus of heterogeneous materials S. Yilmaz Department of Mechanical Engineering, Technical University of Istanbul, 34437 Beyoglu, Istanbul, Turkey
a r t i c l e
i n f o
Article history: Received 3 October 2008 Accepted 6 January 2009 Available online 9 January 2009 Keywords: Shape affect Arrangement affect Orientation affect Heterogeneous materials Elastic modulus Void
a b s t r a c t This paper presents a numerical method for the calculation of the elastic modulus of heterogeneous materials with a continuous matrix, including discontinuous inhomogenity. The effects of the shape, orientation, and arrangement of the inhomogenity were estimated. For an existing heterogeneous material, microstructural geometry was taken into consideration using a simple geometrical parameter for the inhomogenity and the representative volume element (RVE) in which the inhomogenity was embedded. Elastic moduli, predicted using a numerical method, in the case of voids and rigid inhomogenity were in very good agreement with results obtained by the Mori–Tanaka method (MTM), which gives the elastic moduli for the case of an ellipsoidal-shaped inhomogenity. The numerically predicted elastic moduli presented here were also in very good agreement with the results of finite element analysis for different cases of inhomogenity, such as the simple cubic (SC), body-centered cubic (BCC), face-centered cubic (FCC), close packing hexagonal (CPH) type regular arrangements and random distributions. There is an important advantage to the numerical method presented here, as it allows systematic- and parametric investigation of effects of microstructural geometry on the elastic modulus of a heterogeneous material. In order to design a heterogeneous material with an optimal microstructural architecture or tailor a specific elastic modulus, the numerical method proposed in this study can quickly provide information and allows easy implementation of calculations. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction In order to design a heterogeneous material with functional orientation, the effects of microstructural parameters on the material properties should be well understood. The overall properties of a heterogeneous material depend not only on the properties of its constituents but also on the architecture of its microstructure, i.e., the shape, size, orientation, and arrangement geometries of its constituents. Therefore, material engineers need a powerful property prediction method that considers all the requisite microstructural parameters for designing a heterogeneous material. Finite element method (FEM) simulations, using a model based on the real microstructural geometry of heterogeneous materials, can predict the local- and overall behaviors, provided there is sufficient computational power available. The computational limitations are overcome by using the simplification assumption that continuous media can be built up using a representative volume element (RVE) [1] and unit cells [2,3]. Obviously, the success of FEM analysis depends on the validity of the RVE used rather than the actual microstructural geometry. There are, however, several analytical methods proposed for special microstructure cases, i.e., homogenization methods [4,5]. The most sophisticated analytical approach, proposed by Eshelby [6], considers the shape of the
E-mail address: [email protected] 0261-3069/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.matdes.2009.01.001
inhomogenity by means of the Eshelby tensor. Mori–Tanaka [7] modified the Eshelby approach in order to take into account inhomogenity interactions. The analytical solutions of the Mori–Tanaka method (MTM) for ellipsoidal-shaped inhomogenities embedded into ellipsoidal RVEs can be found elsewhere [4,5]. For the case of a continuous matrix with a discontinuous inhomogenity, Hashin and Shtrikman proposed upper- and lower bounds for a spherically shaped inhomogenity [8]. A generalized method for calculating elastic moduli must consider not only the elastic properties and volume fractions of the constituent phases but also the shape, orientation, and arrangement of the inhomogenities. In this work, a numerical approach, which takes into account the geometric properties of the microstructure, is proposed as a method to estimate the elastic moduli of heterogeneous materials with a continuous matrix that have a discontinuous inhomogenity bounded to the matrix. 2. Computational approach and material models For a heterogeneous domain with perfectly bounded double constituents, a continuous matrix, and a discontinuous inhomogenity, the elastic modulus lies in the upper- and lower bounds. The upper bound considers the parallel loading case (i.e., equal strain at the components) [9]. The lower bound, however, considers the serial loading case (i.e., equal stress at components) [10]. Depending
2939
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
on the shape, orientation, and arrangement of inhomogenities, the elastic modulus may be close to either the upper- or lower bound value. It is obvious that the difference between upper- and lower bounds increases as the difference between the properties of the constituents increases. However, if the volume fraction of inhomogenity is very close to 1, the difference in the elastic modulus between the upper- and lower bound calculations remains very small. Hence, it is plausible to estimate the elastic modulus of the heterogeneous domain as a weighted average of the upperand lower bounds as the volume fraction takes very high (close to 1) values. Obviously, the value of the weighting ratio determines the correctness of the estimation. The procedure that was applied in the current paper to determine the weighting ratio is discussed below. Apparently, a matrix region including an inhomogenity with a volume fraction Vf can be composed of a collection of a nested sequence of skins, as shown in Fig. 1. In such a case, each skin in-
a i th matrix skin covering the domain Ω i-1 (i -1)th domain Ω i-1
Inhomogenity domain Ω o
cludes a heterogeneous domain with a volume faction very close to 1. It was assumed that elastic modulus of the ith heterogeneous domain, as shown in Fig. 1, is calculated using the following equations: LB
Ei ¼ ð1 Gf ÞEUB i þ GfEi EUB i
¼
ELB i
¼
S Vf i Em
þ ð1 Em Ei1 S
S
ð1 Vf i ÞEm þ Vf i Ei1
L
b
F
ð3Þ
2ATP AS
ð5Þ
where AS is surface area of inhomogenity, ALP is the projection area of inhomogenity in the direction of loading and ATP is the average projection area of inhomogenity in the transverse direction (Fig. 2). The multiplier 2 represents the top- and bottom faces of the projection area. S If Vf i is the volume fraction of the ith skin covering the inner heterogeneous domain and N is the total number of skins built up in the matrix region, then the volume fraction (Vf) of inhomogenity embedded in the heterogeneous domain can be calculated using the following equation:
Vf ¼
!1 N Y S ð1 þ Vf i Þ
ð6Þ
i¼1
APT
Projection area at transverse direction of θi
Elipsoid shaped inhomogenity
ðLower bound modelÞ
ð4Þ
T
the inhomogenity
ð2Þ
2ALP AS
Gf ¼ 1 Gf ¼ 1
Boundary
ðUpper bound modelÞ
LB where EUB i ; Ei , and Ei are the upper bound, lower bound, and weighted average values for the elastic modulus of the ith heterogeneous domain, respectively. Ei1 is the elastic modulus of the inner (i 1)th heterogeneous domain, which is covered by the ith skin. S Vf i is the relative volume fraction of the ith skin covering the inner (i 1)th heterogeneous domain. Gf is the weighted constant for the lower bound value. It is apparent that the elastic modulus of the heterogeneous domain, including an inhomogenity with a volume fraction of Vf, would equal the elastic modulus value calculated by Eq. (1) for the domain covered by the outermost skin. The ratio of the projection area to the surface area of the inhomogenity can be used to determine Gf, the value of the weighted constant. The projection area can be calculated either parallel or transverse to the loading direction (depending on corresponding parallel or serial loading cases). Thus, the Gf value was calculated by one of the following equations:
Gf ¼ Gf ¼
1th matrix skin covering
ð1Þ
S Vf i ÞEi1
100
z
PAA considering Gf
L
PAA considering Gf
T
c
a x
Log ( E 3 / Ematrix )
MTM
θi b
y
10 RIGID
Prolate elps.
INHOM..
Spheroid
1
Oblate elps.
0.1
VOID
L
Projecion area AP at loading direction
0.01
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 1. (a) Schematic representation of the division of a heterogeneous domain into a nested sequence. (b) Schematic representation of projection areas of a ellipsoidal inhomogenity at loading and transverse directions.
Fig. 2. Normalized elastic moduli in the axial direction versus volume fraction. The heterogeneous domains include prolate ellipsoidal (10a = 10b = c), oblate ellipsoidal (a = b = 10c) and spherical-type inhomogenity.
2940
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
Obviously, Ei1 equals the elastic modulus of inhomogenity for the S first step (i is equal to 1). If Vf i is constant for each skin and equals S Vf , N can be calculated using the following equation:
N¼
lnðV f Þ
ð7Þ
S
lnð1 þ Vf Þ
A set of Visual Basic codes was written in order to compute the numerical calculations mentioned above. For simplicity, the proposed numerical approach is abbreviated to projection area approach (PAA) in the rest of the text. The representative volume element (RVE) models were used for modeling the heterogeneous domain and used to examine the effects of shape, orientation, and arrangement of inhomogenities on the elastic modulus. In some models, the shape of the RVE with embedded inhomogenities was assumed to be similar in shape to the inhomogenity in order to isolate the shape effect from the arrangement effect. For this purpose, the simple case of ellipsoidal inhomogenity embedded in ellipsoidal RVE systems was chosen, precisely like those used in the MTM. Thus, the Gf value was chosen to equal the value of inhomogenity and was kept constant during the numerical iterations relating to the analysis of the shape effect. In order to investigate the arrangement effect, material models with an ideal arrangement of inhomogenities were considered. In an ideal arrangement, a heterogeneous domain can be assumed to have periodically distributed inhomogenities in 3D space. Thus, an inhomogenity-matrix packing system can be established, and a space built up using three-dimensional Voronoi tessellations [11]. Voronoi tessellations, which contained a complete inhomogenity, were used as a RVE for the heterogeneous domain [1]. In order to include the arrangement effect in the analysis, the Gf values were assumed equal to the value of inhomogenity at the first step, and they then linearly approached the value of the RVE in the numerical iteration steps mentioned above. 3. Results 3.1. The effect of shape To determine the pure shape effect, voids were considered inhomogenities in order to eliminate their elastic behavior contribution. Likewise, to determine the extreme stiffness effect of inhomogenities, rigid inhomogenities were also considered. A value of 1032 GPa was used as the elastic modulus to model the case of rigid inhomogenity. The different axisymmetric ellipsoidal geometries investigated are listed in Table 1. The weighted constant, Gf, which was based on GfL or GfT values used in the PAA calculations, are also given in Table 1. Empirical Gfe values, with a relative difference of 0.1% compared to the MTM results, are also tabulated in Table 1 for comparison. The empirical Gf values were determined by trial and error. The results of PAA associated with selected geometries, such as the prolate- and oblate ellipsoid and the spherical cases, are plotted in Fig. 2. The results of the MTM method are also given in Table 1 Gf weighting constant for different ellipsoid geometries and empirical Gfe values, which yield elastic moduli values similar to those obtained from MTM calculations. Ellipsoid geometry
Gf = GfL
Gf = 1 GfT
Gfe
a = b = 10c a = b = 5c a = b = 2c a=b=c 2a = 2b = c 5a = 5b = c 10a = 10b = c
0.971 0.914 0.725 0.500 0.293 0.125 0.063
0.903 0.817 0.638 0.500 0.415 0.374 0.367
0.920 0.850 0.685 0.501 0.303 0.113 0.044
Fig. 2 for comparison. As expected, elastic moduli calculated for domains, decrease with increasing void volume fraction in the case of voids and increase with increasing void volume fraction in the case of rigid inhomogenities. The elastic moduli predictions vary directly with the value of Gf in Eq. (1). As shown in Fig. 2, PAA calculations yielded different elastic moduli curves depending on the Gf parameter, which varies with GfL or GfT. PAA as a function of GfL estimates higher elastic moduli for prolate ellipsoidally shaped (a > b = c) voids and rigid inhomogenities. In contrast, PAA as a function of GfT estimates higher elastic moduli for oblate ellipsoidally shaped (a = b > c) voids or rigid inhomogenity (Fig. 2). The Gf values calculated using GfL were lower than those obtained using GfT for prolate ellipsoids, and Gf values calculated using GfT were lower for oblate ellipsoids, as seen in Table 1. Since Gf determines the contribution of the lower bound value on the elastic modulus in Eq. (1), PAA calculations using Gf were calculated using GfL estimated higher elastic modulus in the case of prolate ellipsoids. Similarly, PAA calculations using Gf calculated using GfT estimated higher elastic modulus in the case of oblate ellipsoid cases. Furthermore, the discrepancy between GfL- and GfT -based estimations were higher for the prolate ellipsoidal geometry than for the oblate one in the rigid inhomogenity case, as seen in Fig. 2. Since there was a large difference between GfL- and (1 GfT) values for prolate ellipsoids and a small difference for oblate ellipsoids (Table 1), this evidence was comprehensible for the rigid inhomogenity case. The Gf parameter dominated the elastic moduli calculations using PAA in the case of rigid inhomogenities. However, there was a large difference between the elastic moduli curves for the oblate ellipsoidal case using voids and a small difference for the prolate ellipsoidal case (Fig. 2), despite small differences in the Gf values calculated using GfL and GfT for the oblate ellipsoidal case, as seen in Table 1. Actually, the (1 GfT) values representing the lower bound contribution in Eq. (1) were almost three times larger than the corresponding parameter, GfL, for the oblate ellipsoidal cases, as indicated in Table 1. The (1 Gf) value sensitivity of PAA was higher than that of Gf in the case of voids, although the Gf value was higher than the (1 Gf) value (Table 1) in the case of oblate ellipsoid type inhomogenities. Therefore, one can say that the upper bound value used in Eq. (1) was dominant in the elastic moduli estimations using PAA in the case of voids. Since GfL and GfT are equal for a spherical inhomogenity, PAA gives similar results for the elastic moduli of spherical inhomogenities (Fig. 2). In addition, the elastic moduli curve of PAA calculated using GfL was closer to the MTM results in the prolate ellipsoidal case. Similarly, the elastic moduli curve of PAA calculated using GfT was closer to the MTM results for the oblate ellipsoidal case (Fig. 2). Furthermore, the agreement between the trend of those closer to the PAA and MTM curves was also excellent, as seen in Fig. 2. The lower values of Gf given in Table 1 were closer to Gfe values, which indicates similar PAA and MTM results, as mentioned above. The lower values of Gf were calculated using GfT values for oblate ellipsoids and GfL values for prolate ellipsoids. Furthermore, the difference between the lower Gf values and Gfe values diminishes if the aspect ratio of inhomogenity approaches 1, as seen in Table 1. 3.2. The effect of orientation In the aforementioned axisymmetric ellipsoidal geometries, the calculated elastic moduli in the transverse direction were plotted in Fig. 3a in order to determine the effects of orientation. The Gf values for transverse direction were also calculated using the data from Table 1. Furthermore, the PAA predictions for the axial directions of a generalized ellipsoidal geometry (c = 2b = 4a) were plotted in Fig. 3b. For clarity, only PAA estimations calculated using small Gf values were plotted in Fig. 3b. Again, the results using
2941
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
a
4
100 L
PAA considering Gf MTM
FEM
HA p /PLLA [12] Al / SiC p [13] PAA
3
Prolate eps.
1
EComposite / EMatrix
Log (E / Ematrix )
10
Oblate elps.
0.1
2 0.01
0
0.2
0.4
0.6
0.8
1
Volume Fraction
b
100
PAA considering Gf PAA considering Gf MTM
L
1
T
Log (E / Ematrix )
0.1
0.2
0.3
0.4
0.5
0.6
Volume Fraction
10
Fig. 4. Normalized elastic moduli versus volume fraction. Various heterogeneous systems with BCC- and FCC arrangements, and comparison with FEM results.
E
1
1
E
2
E
3
Gf 0.1
0.204 0.408 0.693
0.01
0
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 3. (a) Normalized elastic moduli in the transverse direction versus volume fraction. The heterogeneous domains include prolate ellipsoidal (10a = 10b = c) and oblate ellipsoidal type (a = b = 10c) inhomogenity. (b) Normalized elastic moduli in the axial directions versus volume fraction. The heterogeneous domain includes general ellipsoidal type inhomogenity (4a = 2b = c).
the MTM method are also plotted in Fig. 3 for comparison. As seen from Fig. 3, elastic moduli as a function of increasing void volume fractions varied as expected for all loading orientations. The PAA predictions and MTM results were well matched for each of the loading orientations shown in Fig. 3. Therefore, the proposed PAA can represent the effects of orientation. 3.3. The effect of arrangement A space formed by three-dimensional Voronoi tessellations can also be formed by using appropriate unit cell models [11]. For example, typical unit cells, which can be used to form heterogeneous domains with periodic arrangements, may be body-centered cubic (BCC) or face-centered cubic (FCC) as in the conventional space lattice [2]. The elastic moduli of HA/PEEK [12] and Al/SiC
[13] composites calculated using the finite element method with unit cell models for BCC and FCC are compared with the results of the proposed PAA in Fig. 4. The shape of the Voronoi tessellation and, therefore, the shape of the RVE for domains represented by the BCC unit cell is a tetrakaidecahedron. The shape of the RVE for an FCC unit cell is a rhombic dodecahedron. The Gf values for the Voroni tessellation of a BCC and FCC in the [1 0 0] direction are 0.523 and 0.471, respectively. The BCC and FCC unit cells and corresponding RVEs are shown in Fig. 4. For the results shown in Fig. 4, numerically calculated Gf values varied linearly from the spherical inhomogenity (Gf = 0.5) to RVE values. As seen from Fig. 4, the agreement between the PAA results and FEM simulations was very good except for the high volume fractions. The narrow gaps between the inhomogenities at high volume fractions produce Poisson contraction effects. Therefore, the FEM analysis calculated higher values than the PAA predictions at high volume fractions (Fig. 4), as expected. Obviously, the degree of constrained deformations of constituents increases if the difference between the elastic moduli and Poisson ratios of components increases. Since the difference between elastic behaviors of the matrix and inhomogenity was lower in calculations using the BCC model than those using the FCC model, the agreement between the PAA and FEM results for the BCC model remained consistent even at higher volume fractions. In the proposed PAA, only the elastic moduli of the constituents were taken into consideration. However, the agreement between the PAA results and FEM simulations for low- and moderately high volume fractions reveals the competence of the proposed PAA. 4. Discussion 4.1. The effect of inhomogenity stiffness In order to determine the contribution of inhomogenity stiffness, the cases of rigid inhomogenity and Gfe values belonging to
2942
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
the case of voids reflecting only the shape effects (Table 1) were considered. The elastic moduli calculated were normalized by elastic modulus found using the MTM approach, and are plotted with respect to volume fractions in Fig. 5. As shown in Fig. 5a, there were differences between the PAA and MTM results, particularly for the case of an oblate ellipsoid geometry that has more flat surfaces. The constrained deformation of constituents depends not only on the differences in stiffness but also on the mutual Poisson contraction interactions between the constituents. Calculation methods differ in the way they account for the elastic interaction between inhomogenities. MTM calculations consider the contribution of inhomogenity stiffness using the elastic modulus and Poisson ratios of the matrix material only [4,5]. The differences between the results of the two approaches are seen in Fig. 5a, clarifying the effect of Poisson contraction (i.e., contribution of constrained deformation). Therefore, one can say that the Poisson contraction effect due to stiffer inhomogenities is compounded with the pure geometric effect of the inhomogenity and makes the heterogeneous domain stiffer, as expected. The Poisson contraction contribution to the stiffness of the heterogeneous domain is higher if the value of Gf increases. In other words, Poisson contraction effects at the narrow gaps between the matrix and inhomogenities should be considered in the PAA calculations,
a
1.1 1.05
Prolate elps 10a=10b=c
.
E 3 PAA / E MTM
1 Spheroid 0.95 Oblate elps. a = b = 5c
0.9
a= b=10c
0.85 0.8
0
0.2
0.4
0.6
0.8
1
Volume Fraction
b
1.2
E 3 PAA / E MTM
Eem ¼ EPS ¼
E ð1 v 2 Þ
ð8Þ
where Eem is the equivalent elastic modulus and EPS is the plane strain elastic modulus. Similarly, full Poisson contraction of an elastic matrix due to the anchored boundary condition of rigid inhomogenity can be modeled by the following equation:
Eem ¼ EAS ¼
E ð1 2v 2 Þ
ð9Þ
where EAS is the elastic modulus in the axial (one-dimensional) strain case. Therefore, the equivalent stiffness values of the matrix given by Eqs. (8) and (9) can be used in the PAA approach referred above for simulating the Poisson contraction effect. The elastic moduli calculated using the values for Gfe and Eem were normalized using the elastic modulus calculated using the MTM approach, and the results are plotted in Fig. 5b. Apparently, constrained deformation of the matrix is valid at high volume fractions. As a result, the PAA estimations were higher than the MTM results for lower volume fractions. The agreement between both methods increased with increasing volume fractions for equivalent elastic moduli, simulating constrained deformation of the matrix, as shown in Fig. 5b. For the oblate ellipsoid case of (a = b = 5c), results of the plain strain elastic modulus yields are closer to those from the MTM for high volume fractions. However, for the oblate ellipsoid case of (a = b = 10c), the fully constrained matrix model simulated by Eq. (9) yields results are closer to those obtained from the MTM for high volume fractions. Obviously, increasing the flatness of inhomogenity increases the Poisson contraction effect. Since the ellipsoid with radius ratio a = b = 5c has more concave surfaces than the ellipsoid with a radius ratio of a = b = 10c, a plain strain elastic modulus was sufficient to simulate the Poisson contraction effect for high volume fractions in the oblate ellipsoid case (a = b = 5c). The PAA method estimations, using both the regular- and effective elastic modulus of the matrix mentioned above, also give the property band for the elastic modulus of heterogeneous materials with a continuous matrix and discontinuous inhomogenity. 4.2. Simple cubic (SC) arrangement
Oblate elps. ------ a = b = 5c —— a = b =10c
Eme=EAS 1
Eme=EPS
Eme=Em 0.8
particularly if the rigid inhomogenity has flat surfaces (i.e., platelet-type inhomogenities). In order to clarify this situation, the equivalent stiffness of a matrix under constrained deformation due to the Poisson contraction effect can be modeled using the plain strain elastic modulus:
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 5. (a) The differences between the results of the PAA using the equivalent value of Gfe and results from the MTM approaches for the case of rigid inhomogenity. (b) Comparison of PAA with MTM results. In the PAA calculations, an equivalent value of Gfe and a different equivalent value of Eem , which considers Poisson contraction effects on the matrix were used.
A transversely aligned arrangement is the inhomogenity arrangement in which similarly signed stress fields around the nearest-neighbor inhomogenities overlap [13]. The unit cell for the transversely aligned arrangement is similar to the simple cubic (SC) space lattice. The RVE for the SC unit cell that includes spherical inhomogenity has cubic geometry. Since the heterogeneous domain was assumed to be a 3D infinitely periodic array of transversely aligned particles embedded in the matrix, the stress–strain distributions in the unit cells are symmetric. Due to the shape and symmetric distribution of inhomogenity, the boundary conditions are defined so that the angles between the surfaces of the RVE are preserved; thus, deformed unit cell surfaces remain parallel. Therefore, the outermost layer of the RVE is under a parallel loading condition. The Gf value for the Voronoi cell of an SC unit cell becomes zero (0) since it has parallelepiped walls and the angles between walls non-deform as if they are under parallel loading. Therefore, the Gf values were initially set to the value of the inhomogenity and then reduced linearly down to 0 during the numerical iteration carried out for an SC arrangement. In Fig. 6a, FEM calculations and PAA estimations were given for the different het-
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
a
4.3. The effect of orientation for SC-, BCC-, FCC-, and CPH arrangements
4
PAA FEM HAp /PLLA [12] Glass /Epoxy [13] Al /SiC p [13]
E Composite / EMatrix
3
2
1 0
0.1
0.2
0.3
0.4
0.5
Volume Fraction
b
2943
4
FEM U&L bounds[14]
PAA
E Composite / EMatrix
3
If the Gf values of inhomogenity and of the Voronoi tessellation RVE are calculated with respect to the loading direction, the effects of orientation and arrangement can be considered. In Table 2, the elastic moduli for an Al/SiCp system calculated by the finite element method using different unit cell models and different particulate shapes are compared with the results of PAA. The Gf values at the [0 0 1] and [1 1 0] directions for RVE Voroni cells of SC-, BCC-, and FCC unit cells are given in Table 2. The Gf values for spherical-, cubic-, and cylindrical particulates are also given in Table 2 for the same directions. As seen from Table 2, there was very good agreement between the FEM results and PAA estimations. The Voronoi tessellations for the CPH and FCC arrangements are identical. The Gf values for the Voronoi CPH unit cells are 0.450 in the [0 0 1] direction and 0.491 in the [1 2 0] direction. For an epoxy– 30% glass particulate heterogeneous system (EEpoxy = 3.01 GPa, vEpoxy = 0.394, EGlass = 76 GPa, vGlass = 0.23), the elastic moduli were calculated as 5.810 and 5.681 GPa for the [0 0 1] and [1 2 0] orientations, using an FEM analysis for a solid model with CPH unit cells and spherical particulates [1]. The proposed PAA estimated elastic moduli of 5.470 and 5.344 GPa for the [0 0 1] and [1 2 0] orientations, respectively. The PAA values were less than, but closely approached, the FEM results, which was to be expected since the PAA does not consider the Poisson contraction effect. Although the PAA predictions yielded lower values than those obtained using the FEM analysis, it should be noted that the proposed PAA approach was able to predict higher elastic modulus in the [0 0 1] direction than in the [1 2 0] direction, precisely like those obtained using the FEM analysis. It should also be remembered that Voronoi tessellation of the CPH arrangement has some parallelepiped surfaces retained during deformation such as the SC arrangement mentioned above. If a Gf value of 0.29 was used for the RVE, the elastic modulus was calculated to be 5.621 GPa in the [0 0 1] direction, which was closer to the FEM predicted value of 5.810 GPa. This was because 44% of the total RVE surfaces remain parallelepiped during deformation in the [0 0 1] direction (i.e., (0.44 0 + 0.66 0.45 = 0.29)).
2
4.4. The effect of randomly arranged spherical inhomogenity
1
0
0.1
0.2
0.3
0.4
0.5
Volume Fraction Fig. 6. (a) Normalized elastic moduli of different heterogeneous systems having an SC arrangement and comparison with FEM results. (b) Normalized elastic moduli of a glass–epoxy system with an SC arrangement and comparison with FEM results for iso-displacement and iso-stress loading conditions.
erogeneous systems. As seen in Fig. 6, there is a good correlation between the PAA estimations and FEM results. The expected discrepancy between the FEM results and PAA estimations seen for high volume fractions is due to the constrained deformation of the matrix [13]. The FEM results for the elastic moduli of a glass– epoxy heterogeneous system is also shown in Fig. 6b. For the FEM results [14] given in Fig. 6b, the SC unit cell model is loaded with iso-displacement from one face, as mentioned above, and with iso-stress, in order to reduce the Poisson contraction effect. The PAA curve lies between those lower- and upper bound values as expected (Fig. 6b).
To obtain an effective response for a periodic arrangement, one can use either the smallest RVE or the smallest unit cell with a smaller number of inhomogenities [15]. For a random arrangement, however, the RVE must include a sufficient number of inhomogenities. Thus, the behavior of the RVE represents the overall behavior of the heterogeneous domain. Since the RVE in the proposed PAA calculations includes only a single inhomogenity, the overall effect of arrangement was reflected in the Gf parameter. Obviously, each of the sub-volumes with a single inhomogenity in a random arrangement may have different shapes (irregular or regular). Further, each individual sub-volume may have parallelepiped surfaces retained during deformation as mentioned above for RVE of the CPH arrangement. Hence, an imaginary RVE may be considered to represent the random arrangement case so that it would have an appropriate Gf value that reflects the overall effect of random arrangement. The elastic moduli of a continuous-, isotropic-, and elastic domain with statically isotropic-, randomly arranged spherical inhomogenities was calculated for extreme cases (with voids and with spheres assumed to be rigid), analytically using Torquato’s third order approximation (TOA) [16] and numerically by FEM simulations [17]. A cubic RVE cell containing 30 spheres randomly located with periodic boundary conditions was used in the FEM analysis. Although the Gf can take a value between 0 and 1, a value of
2944
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
Table 2 Elastic moduli of a heterogeneous domain for different arrangements and loading orientations calculated using the finite element method with unit cell models [2] and by the proposed PAA. System
Al/20%SiCp EAl = 67.2 GPa, vAl = 0.35 ESiCp = 429 GPa vSiCp = 0.17
Unit cell
Inhmgn. shape
FEM (GPa)
SC
Sphere Cube Cylinder Sphere Cube Cylinder Sphere Cube Cylinder
96.4 101.0 102.6 89.9 – – 90.0 94.4 94.5
BCC
FCC
E
001
PAA (GPa) E
110
90.2 92.1 91.7 91.1 – – 91.6 93.8 91.4
0.333 was found adequate to reflect the overall behavior of a random arrangement, as determined by trial and error based on the data obtained by the TOA and FEM results shown in Fig. 7. There was very good agreement between the FEM and TOA estimations at low volume fractions. The Gf value of 0.333, which is smaller than 0.5, implies that the parallel loading condition is crucial for the arrangement effects of random distribution. The proposed PAA estimated harder behavior in the high volume fractions for the void case. Matrix regions, which were under localized stress– strain intensities, expanded with increasing volume fractions; thus, the FEM calculations showed softer elastic behavior in the high volume fractions, as expected for the void case. In contrast, the proposed PAA estimated softer behavior for the rigid inhomogenity case with higher volume fractions, as seen in Fig. 7. Obviously, the constrained deformation of the matrix was pronounced with increasing volume fractions in the case of rigid inhomogenity. Therefore, increasing the stiffness of the matrix due to the elastic interaction between inhomogenity and the matrix resulted in harder FEM calculated elastic behavior, as was expected for the rigid inhomogenity case. However, general agreement of the PAA curve with the curves of the analytical TOA and FEM simulations indicates that the proposed PAA can follow the arrangement effects of a random distribution. It should be remembered that more valid Gf values that reflect real situa-
4
PAA TOA [16] FEM [17]
E / E matrix
3
2
Rigid inhomogenity 1 Void
0 0
0.1
0.2
0.3
0.4
E
001
E
110
Gf001 for RVE
Gf110 for RVE
Gf001 for Inhm.
Gf110 for Inhm.
0.500 0.333 0.333 0.500 0.333 0.333 0.500 0.333 0.333
0.500 0.471 0.424 0.500 0.471 0.424 0.500 0.471 0.424
95.3 102.5 102.5 90.3
90.7 91.5 92.7 91.1
0.000
0.471
0.523
0.422
90.7 95.6 95.6
90.5 91.2 92.4
0.471
0.500
tions could be determined by comparison with experimental results. 5. Concluding remarks
1. The elastic modulus of heterogeneous domains with a continuous matrix and discontinuous inhomogenities was estimated using a numerical calculation method. This method can consider the shape, orientation and arrangement effects of inhomogenity via a simple parameter, the geometrical factor (Gf). Gf values vary from 0 to 1 and indicate stress–strain distribution conditions of the constituents in a heterogeneous domain. The Gf value approaches 0 and 1 as the stress–strain distributions approach the parallel- and serial loading conditions, respectively. 2. The shape and orientation effects of the inhomogenity were considered with the Gf value of the inhomogenity. Arrangement effect of the inhomogenity was taken into account by using the Gf value belonging to a representative volume element (RVE) in which the inhomogenity was embedded. 3. In general, the Gf was determined using the ratios (i) between the projection area in the loading direction and the surface area of the inhomogenity and (ii) between the projection area in the loading direction and the surface area of RVE. If the RVE, with the inhomogenity, had parallelepiped surfaces during deformation, such as that in an SC arrangement, Gf values were chosen to be 0 (i.e., the parallel loading condition is dominant in the heterogeneous domain). 4. There was very good agreement between the estimations of the proposed numerical approach (PAA) and results of analytical calculations obtained from MTM and TOA. There was also very good agreement between the estimations of the PAA and results of FEM simulations taken from the literature for SC-, BCC-, FCC-, and CPH packing systems and for unit cells belonging to randomly distributed inhomogenities. 5. The PAA can be used in systematic- and parametric investigations into the effects of microstructural geometry on the elastic modulus of a heterogeneous domain. The numerical method proposed can provide quick information and is easily implemented by engineers. 6. The PAA still needs to consider Poisson contraction effects appearing at high volume fractions. However, the PAA can provide property bands for elastic moduli of heterogeneous domains (for a matrix which is under plane strain condition and for a matrix which is free from Poisson contraction).
0.5
Volume Fraction Fig. 7. Normalized elastic moduli of a heterogeneous domain with randomly arranged voids and rigid inhomogenities and comparison with both FEM and TOA results.
Acknowledgment The author gratefully acknowledges financial support of the Technical University of Istanbul, BAP Program.
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
References [1] Li S, Wongsto A. Unit cells for micromechanical analyses of particle-reinforced heterogeneous. Mech Mater 2004;26:543–72. [2] Weissenbek E, Bohm HJ, Rammerstorfer FG. Micromechanical investigation of arrangement effects in particle reinforced metal matrix heterogeneous. Comput Mater Sci 1994;3:263–78. [3] Levy A, Papazian JM. Elasto-plastic finite element analysis of short fiber reinforced Sic/Al heterogeneous: effect of thermal treatment. Acta Metall Mater 1991;39:2255–66. [4] Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials. The Netherlands: Elsevier; 1999. [5] Mura T. Micromechanics of defect in solids. The Netherlands: Kluwer; 1998. [6] Eshelby CC. The determination of the elastic field of an ellipsoidal inclusion, and related field. Proc R Soc Lond 1957;A241:376–96. [7] Mori T, Tanaka K. Average stress in the matrix and average elastic energy of materials with misfitting inclusions. Acta Metall Mater 1973;21:571–4. [8] Hashin Z, Shtrikman S. A variational approach to the theory of elastic behavior of multiphase materials. J Mech Phys Solid 1963;11:127–40.
2945
[9] Voight W. Lehrbuch Der Kristallphysik. Berlin-Leipzig: Teubner Verlag; 1928. [10] Reuss A. Berechung der fliessgrenzen von mischkristallen aurf grund der plastiizitatsbedingug fur einkristalle. Z Angew Math Mech 1929;9:49–58. [11] Ahuja N, Schachter BJ. Pattern models. New York: Wiley; 1983. [12] Balac I, Milovancevic M, Tang C, Uskokovic PS, Uskokovic DP. Estimation of elastic properties of a particulate polymer composite using a face-centered cubic FE model. Mater Lett 2004;58:2437–41. [13] Yilmaz S, Aran A. Finite element analysis of deformation behavior in ductile matrix containing hard particles. Mater Sci Tech 1998;14:1154–62. [14] Sun C, Saffari P, Ranade R, Sadeghipour K, Baran G. Finite element analysis of elastic property bounds of a composite with randomly distributed particles. Composites 2007;A38:80–6. [15] Jiang M, Jasiuk I, Ostoja-Starzewski M. Apparent elastic and elastoplastic behavior of periodic composites. Int J Solids Struct 2002;39:199–212. [16] Torquato S. Effective stiffness tensor of heterogeneous media. II: Application to isotropic dispersions. J Mech Phys Solids 1998;46:1411–40. [17] Segurado J, Llorca J. A numerical approximation to the elastic properties of sphere reinforced heterogeneous. J Mech Phys Solids 2002;50:2107–21.
Contents lists available at ScienceDirect
Materials and Design journal homepage: www.elsevier.com/locate/matdes
An approach for predicting the elastic modulus of heterogeneous materials S. Yilmaz Department of Mechanical Engineering, Technical University of Istanbul, 34437 Beyoglu, Istanbul, Turkey
a r t i c l e
i n f o
Article history: Received 3 October 2008 Accepted 6 January 2009 Available online 9 January 2009 Keywords: Shape affect Arrangement affect Orientation affect Heterogeneous materials Elastic modulus Void
a b s t r a c t This paper presents a numerical method for the calculation of the elastic modulus of heterogeneous materials with a continuous matrix, including discontinuous inhomogenity. The effects of the shape, orientation, and arrangement of the inhomogenity were estimated. For an existing heterogeneous material, microstructural geometry was taken into consideration using a simple geometrical parameter for the inhomogenity and the representative volume element (RVE) in which the inhomogenity was embedded. Elastic moduli, predicted using a numerical method, in the case of voids and rigid inhomogenity were in very good agreement with results obtained by the Mori–Tanaka method (MTM), which gives the elastic moduli for the case of an ellipsoidal-shaped inhomogenity. The numerically predicted elastic moduli presented here were also in very good agreement with the results of finite element analysis for different cases of inhomogenity, such as the simple cubic (SC), body-centered cubic (BCC), face-centered cubic (FCC), close packing hexagonal (CPH) type regular arrangements and random distributions. There is an important advantage to the numerical method presented here, as it allows systematic- and parametric investigation of effects of microstructural geometry on the elastic modulus of a heterogeneous material. In order to design a heterogeneous material with an optimal microstructural architecture or tailor a specific elastic modulus, the numerical method proposed in this study can quickly provide information and allows easy implementation of calculations. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction In order to design a heterogeneous material with functional orientation, the effects of microstructural parameters on the material properties should be well understood. The overall properties of a heterogeneous material depend not only on the properties of its constituents but also on the architecture of its microstructure, i.e., the shape, size, orientation, and arrangement geometries of its constituents. Therefore, material engineers need a powerful property prediction method that considers all the requisite microstructural parameters for designing a heterogeneous material. Finite element method (FEM) simulations, using a model based on the real microstructural geometry of heterogeneous materials, can predict the local- and overall behaviors, provided there is sufficient computational power available. The computational limitations are overcome by using the simplification assumption that continuous media can be built up using a representative volume element (RVE) [1] and unit cells [2,3]. Obviously, the success of FEM analysis depends on the validity of the RVE used rather than the actual microstructural geometry. There are, however, several analytical methods proposed for special microstructure cases, i.e., homogenization methods [4,5]. The most sophisticated analytical approach, proposed by Eshelby [6], considers the shape of the
E-mail address: [email protected] 0261-3069/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.matdes.2009.01.001
inhomogenity by means of the Eshelby tensor. Mori–Tanaka [7] modified the Eshelby approach in order to take into account inhomogenity interactions. The analytical solutions of the Mori–Tanaka method (MTM) for ellipsoidal-shaped inhomogenities embedded into ellipsoidal RVEs can be found elsewhere [4,5]. For the case of a continuous matrix with a discontinuous inhomogenity, Hashin and Shtrikman proposed upper- and lower bounds for a spherically shaped inhomogenity [8]. A generalized method for calculating elastic moduli must consider not only the elastic properties and volume fractions of the constituent phases but also the shape, orientation, and arrangement of the inhomogenities. In this work, a numerical approach, which takes into account the geometric properties of the microstructure, is proposed as a method to estimate the elastic moduli of heterogeneous materials with a continuous matrix that have a discontinuous inhomogenity bounded to the matrix. 2. Computational approach and material models For a heterogeneous domain with perfectly bounded double constituents, a continuous matrix, and a discontinuous inhomogenity, the elastic modulus lies in the upper- and lower bounds. The upper bound considers the parallel loading case (i.e., equal strain at the components) [9]. The lower bound, however, considers the serial loading case (i.e., equal stress at components) [10]. Depending
2939
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
on the shape, orientation, and arrangement of inhomogenities, the elastic modulus may be close to either the upper- or lower bound value. It is obvious that the difference between upper- and lower bounds increases as the difference between the properties of the constituents increases. However, if the volume fraction of inhomogenity is very close to 1, the difference in the elastic modulus between the upper- and lower bound calculations remains very small. Hence, it is plausible to estimate the elastic modulus of the heterogeneous domain as a weighted average of the upperand lower bounds as the volume fraction takes very high (close to 1) values. Obviously, the value of the weighting ratio determines the correctness of the estimation. The procedure that was applied in the current paper to determine the weighting ratio is discussed below. Apparently, a matrix region including an inhomogenity with a volume fraction Vf can be composed of a collection of a nested sequence of skins, as shown in Fig. 1. In such a case, each skin in-
a i th matrix skin covering the domain Ω i-1 (i -1)th domain Ω i-1
Inhomogenity domain Ω o
cludes a heterogeneous domain with a volume faction very close to 1. It was assumed that elastic modulus of the ith heterogeneous domain, as shown in Fig. 1, is calculated using the following equations: LB
Ei ¼ ð1 Gf ÞEUB i þ GfEi EUB i
¼
ELB i
¼
S Vf i Em
þ ð1 Em Ei1 S
S
ð1 Vf i ÞEm þ Vf i Ei1
L
b
F
ð3Þ
2ATP AS
ð5Þ
where AS is surface area of inhomogenity, ALP is the projection area of inhomogenity in the direction of loading and ATP is the average projection area of inhomogenity in the transverse direction (Fig. 2). The multiplier 2 represents the top- and bottom faces of the projection area. S If Vf i is the volume fraction of the ith skin covering the inner heterogeneous domain and N is the total number of skins built up in the matrix region, then the volume fraction (Vf) of inhomogenity embedded in the heterogeneous domain can be calculated using the following equation:
Vf ¼
!1 N Y S ð1 þ Vf i Þ
ð6Þ
i¼1
APT
Projection area at transverse direction of θi
Elipsoid shaped inhomogenity
ðLower bound modelÞ
ð4Þ
T
the inhomogenity
ð2Þ
2ALP AS
Gf ¼ 1 Gf ¼ 1
Boundary
ðUpper bound modelÞ
LB where EUB i ; Ei , and Ei are the upper bound, lower bound, and weighted average values for the elastic modulus of the ith heterogeneous domain, respectively. Ei1 is the elastic modulus of the inner (i 1)th heterogeneous domain, which is covered by the ith skin. S Vf i is the relative volume fraction of the ith skin covering the inner (i 1)th heterogeneous domain. Gf is the weighted constant for the lower bound value. It is apparent that the elastic modulus of the heterogeneous domain, including an inhomogenity with a volume fraction of Vf, would equal the elastic modulus value calculated by Eq. (1) for the domain covered by the outermost skin. The ratio of the projection area to the surface area of the inhomogenity can be used to determine Gf, the value of the weighted constant. The projection area can be calculated either parallel or transverse to the loading direction (depending on corresponding parallel or serial loading cases). Thus, the Gf value was calculated by one of the following equations:
Gf ¼ Gf ¼
1th matrix skin covering
ð1Þ
S Vf i ÞEi1
100
z
PAA considering Gf
L
PAA considering Gf
T
c
a x
Log ( E 3 / Ematrix )
MTM
θi b
y
10 RIGID
Prolate elps.
INHOM..
Spheroid
1
Oblate elps.
0.1
VOID
L
Projecion area AP at loading direction
0.01
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 1. (a) Schematic representation of the division of a heterogeneous domain into a nested sequence. (b) Schematic representation of projection areas of a ellipsoidal inhomogenity at loading and transverse directions.
Fig. 2. Normalized elastic moduli in the axial direction versus volume fraction. The heterogeneous domains include prolate ellipsoidal (10a = 10b = c), oblate ellipsoidal (a = b = 10c) and spherical-type inhomogenity.
2940
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
Obviously, Ei1 equals the elastic modulus of inhomogenity for the S first step (i is equal to 1). If Vf i is constant for each skin and equals S Vf , N can be calculated using the following equation:
N¼
lnðV f Þ
ð7Þ
S
lnð1 þ Vf Þ
A set of Visual Basic codes was written in order to compute the numerical calculations mentioned above. For simplicity, the proposed numerical approach is abbreviated to projection area approach (PAA) in the rest of the text. The representative volume element (RVE) models were used for modeling the heterogeneous domain and used to examine the effects of shape, orientation, and arrangement of inhomogenities on the elastic modulus. In some models, the shape of the RVE with embedded inhomogenities was assumed to be similar in shape to the inhomogenity in order to isolate the shape effect from the arrangement effect. For this purpose, the simple case of ellipsoidal inhomogenity embedded in ellipsoidal RVE systems was chosen, precisely like those used in the MTM. Thus, the Gf value was chosen to equal the value of inhomogenity and was kept constant during the numerical iterations relating to the analysis of the shape effect. In order to investigate the arrangement effect, material models with an ideal arrangement of inhomogenities were considered. In an ideal arrangement, a heterogeneous domain can be assumed to have periodically distributed inhomogenities in 3D space. Thus, an inhomogenity-matrix packing system can be established, and a space built up using three-dimensional Voronoi tessellations [11]. Voronoi tessellations, which contained a complete inhomogenity, were used as a RVE for the heterogeneous domain [1]. In order to include the arrangement effect in the analysis, the Gf values were assumed equal to the value of inhomogenity at the first step, and they then linearly approached the value of the RVE in the numerical iteration steps mentioned above. 3. Results 3.1. The effect of shape To determine the pure shape effect, voids were considered inhomogenities in order to eliminate their elastic behavior contribution. Likewise, to determine the extreme stiffness effect of inhomogenities, rigid inhomogenities were also considered. A value of 1032 GPa was used as the elastic modulus to model the case of rigid inhomogenity. The different axisymmetric ellipsoidal geometries investigated are listed in Table 1. The weighted constant, Gf, which was based on GfL or GfT values used in the PAA calculations, are also given in Table 1. Empirical Gfe values, with a relative difference of 0.1% compared to the MTM results, are also tabulated in Table 1 for comparison. The empirical Gf values were determined by trial and error. The results of PAA associated with selected geometries, such as the prolate- and oblate ellipsoid and the spherical cases, are plotted in Fig. 2. The results of the MTM method are also given in Table 1 Gf weighting constant for different ellipsoid geometries and empirical Gfe values, which yield elastic moduli values similar to those obtained from MTM calculations. Ellipsoid geometry
Gf = GfL
Gf = 1 GfT
Gfe
a = b = 10c a = b = 5c a = b = 2c a=b=c 2a = 2b = c 5a = 5b = c 10a = 10b = c
0.971 0.914 0.725 0.500 0.293 0.125 0.063
0.903 0.817 0.638 0.500 0.415 0.374 0.367
0.920 0.850 0.685 0.501 0.303 0.113 0.044
Fig. 2 for comparison. As expected, elastic moduli calculated for domains, decrease with increasing void volume fraction in the case of voids and increase with increasing void volume fraction in the case of rigid inhomogenities. The elastic moduli predictions vary directly with the value of Gf in Eq. (1). As shown in Fig. 2, PAA calculations yielded different elastic moduli curves depending on the Gf parameter, which varies with GfL or GfT. PAA as a function of GfL estimates higher elastic moduli for prolate ellipsoidally shaped (a > b = c) voids and rigid inhomogenities. In contrast, PAA as a function of GfT estimates higher elastic moduli for oblate ellipsoidally shaped (a = b > c) voids or rigid inhomogenity (Fig. 2). The Gf values calculated using GfL were lower than those obtained using GfT for prolate ellipsoids, and Gf values calculated using GfT were lower for oblate ellipsoids, as seen in Table 1. Since Gf determines the contribution of the lower bound value on the elastic modulus in Eq. (1), PAA calculations using Gf were calculated using GfL estimated higher elastic modulus in the case of prolate ellipsoids. Similarly, PAA calculations using Gf calculated using GfT estimated higher elastic modulus in the case of oblate ellipsoid cases. Furthermore, the discrepancy between GfL- and GfT -based estimations were higher for the prolate ellipsoidal geometry than for the oblate one in the rigid inhomogenity case, as seen in Fig. 2. Since there was a large difference between GfL- and (1 GfT) values for prolate ellipsoids and a small difference for oblate ellipsoids (Table 1), this evidence was comprehensible for the rigid inhomogenity case. The Gf parameter dominated the elastic moduli calculations using PAA in the case of rigid inhomogenities. However, there was a large difference between the elastic moduli curves for the oblate ellipsoidal case using voids and a small difference for the prolate ellipsoidal case (Fig. 2), despite small differences in the Gf values calculated using GfL and GfT for the oblate ellipsoidal case, as seen in Table 1. Actually, the (1 GfT) values representing the lower bound contribution in Eq. (1) were almost three times larger than the corresponding parameter, GfL, for the oblate ellipsoidal cases, as indicated in Table 1. The (1 Gf) value sensitivity of PAA was higher than that of Gf in the case of voids, although the Gf value was higher than the (1 Gf) value (Table 1) in the case of oblate ellipsoid type inhomogenities. Therefore, one can say that the upper bound value used in Eq. (1) was dominant in the elastic moduli estimations using PAA in the case of voids. Since GfL and GfT are equal for a spherical inhomogenity, PAA gives similar results for the elastic moduli of spherical inhomogenities (Fig. 2). In addition, the elastic moduli curve of PAA calculated using GfL was closer to the MTM results in the prolate ellipsoidal case. Similarly, the elastic moduli curve of PAA calculated using GfT was closer to the MTM results for the oblate ellipsoidal case (Fig. 2). Furthermore, the agreement between the trend of those closer to the PAA and MTM curves was also excellent, as seen in Fig. 2. The lower values of Gf given in Table 1 were closer to Gfe values, which indicates similar PAA and MTM results, as mentioned above. The lower values of Gf were calculated using GfT values for oblate ellipsoids and GfL values for prolate ellipsoids. Furthermore, the difference between the lower Gf values and Gfe values diminishes if the aspect ratio of inhomogenity approaches 1, as seen in Table 1. 3.2. The effect of orientation In the aforementioned axisymmetric ellipsoidal geometries, the calculated elastic moduli in the transverse direction were plotted in Fig. 3a in order to determine the effects of orientation. The Gf values for transverse direction were also calculated using the data from Table 1. Furthermore, the PAA predictions for the axial directions of a generalized ellipsoidal geometry (c = 2b = 4a) were plotted in Fig. 3b. For clarity, only PAA estimations calculated using small Gf values were plotted in Fig. 3b. Again, the results using
2941
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
a
4
100 L
PAA considering Gf MTM
FEM
HA p /PLLA [12] Al / SiC p [13] PAA
3
Prolate eps.
1
EComposite / EMatrix
Log (E / Ematrix )
10
Oblate elps.
0.1
2 0.01
0
0.2
0.4
0.6
0.8
1
Volume Fraction
b
100
PAA considering Gf PAA considering Gf MTM
L
1
T
Log (E / Ematrix )
0.1
0.2
0.3
0.4
0.5
0.6
Volume Fraction
10
Fig. 4. Normalized elastic moduli versus volume fraction. Various heterogeneous systems with BCC- and FCC arrangements, and comparison with FEM results.
E
1
1
E
2
E
3
Gf 0.1
0.204 0.408 0.693
0.01
0
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 3. (a) Normalized elastic moduli in the transverse direction versus volume fraction. The heterogeneous domains include prolate ellipsoidal (10a = 10b = c) and oblate ellipsoidal type (a = b = 10c) inhomogenity. (b) Normalized elastic moduli in the axial directions versus volume fraction. The heterogeneous domain includes general ellipsoidal type inhomogenity (4a = 2b = c).
the MTM method are also plotted in Fig. 3 for comparison. As seen from Fig. 3, elastic moduli as a function of increasing void volume fractions varied as expected for all loading orientations. The PAA predictions and MTM results were well matched for each of the loading orientations shown in Fig. 3. Therefore, the proposed PAA can represent the effects of orientation. 3.3. The effect of arrangement A space formed by three-dimensional Voronoi tessellations can also be formed by using appropriate unit cell models [11]. For example, typical unit cells, which can be used to form heterogeneous domains with periodic arrangements, may be body-centered cubic (BCC) or face-centered cubic (FCC) as in the conventional space lattice [2]. The elastic moduli of HA/PEEK [12] and Al/SiC
[13] composites calculated using the finite element method with unit cell models for BCC and FCC are compared with the results of the proposed PAA in Fig. 4. The shape of the Voronoi tessellation and, therefore, the shape of the RVE for domains represented by the BCC unit cell is a tetrakaidecahedron. The shape of the RVE for an FCC unit cell is a rhombic dodecahedron. The Gf values for the Voroni tessellation of a BCC and FCC in the [1 0 0] direction are 0.523 and 0.471, respectively. The BCC and FCC unit cells and corresponding RVEs are shown in Fig. 4. For the results shown in Fig. 4, numerically calculated Gf values varied linearly from the spherical inhomogenity (Gf = 0.5) to RVE values. As seen from Fig. 4, the agreement between the PAA results and FEM simulations was very good except for the high volume fractions. The narrow gaps between the inhomogenities at high volume fractions produce Poisson contraction effects. Therefore, the FEM analysis calculated higher values than the PAA predictions at high volume fractions (Fig. 4), as expected. Obviously, the degree of constrained deformations of constituents increases if the difference between the elastic moduli and Poisson ratios of components increases. Since the difference between elastic behaviors of the matrix and inhomogenity was lower in calculations using the BCC model than those using the FCC model, the agreement between the PAA and FEM results for the BCC model remained consistent even at higher volume fractions. In the proposed PAA, only the elastic moduli of the constituents were taken into consideration. However, the agreement between the PAA results and FEM simulations for low- and moderately high volume fractions reveals the competence of the proposed PAA. 4. Discussion 4.1. The effect of inhomogenity stiffness In order to determine the contribution of inhomogenity stiffness, the cases of rigid inhomogenity and Gfe values belonging to
2942
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
the case of voids reflecting only the shape effects (Table 1) were considered. The elastic moduli calculated were normalized by elastic modulus found using the MTM approach, and are plotted with respect to volume fractions in Fig. 5. As shown in Fig. 5a, there were differences between the PAA and MTM results, particularly for the case of an oblate ellipsoid geometry that has more flat surfaces. The constrained deformation of constituents depends not only on the differences in stiffness but also on the mutual Poisson contraction interactions between the constituents. Calculation methods differ in the way they account for the elastic interaction between inhomogenities. MTM calculations consider the contribution of inhomogenity stiffness using the elastic modulus and Poisson ratios of the matrix material only [4,5]. The differences between the results of the two approaches are seen in Fig. 5a, clarifying the effect of Poisson contraction (i.e., contribution of constrained deformation). Therefore, one can say that the Poisson contraction effect due to stiffer inhomogenities is compounded with the pure geometric effect of the inhomogenity and makes the heterogeneous domain stiffer, as expected. The Poisson contraction contribution to the stiffness of the heterogeneous domain is higher if the value of Gf increases. In other words, Poisson contraction effects at the narrow gaps between the matrix and inhomogenities should be considered in the PAA calculations,
a
1.1 1.05
Prolate elps 10a=10b=c
.
E 3 PAA / E MTM
1 Spheroid 0.95 Oblate elps. a = b = 5c
0.9
a= b=10c
0.85 0.8
0
0.2
0.4
0.6
0.8
1
Volume Fraction
b
1.2
E 3 PAA / E MTM
Eem ¼ EPS ¼
E ð1 v 2 Þ
ð8Þ
where Eem is the equivalent elastic modulus and EPS is the plane strain elastic modulus. Similarly, full Poisson contraction of an elastic matrix due to the anchored boundary condition of rigid inhomogenity can be modeled by the following equation:
Eem ¼ EAS ¼
E ð1 2v 2 Þ
ð9Þ
where EAS is the elastic modulus in the axial (one-dimensional) strain case. Therefore, the equivalent stiffness values of the matrix given by Eqs. (8) and (9) can be used in the PAA approach referred above for simulating the Poisson contraction effect. The elastic moduli calculated using the values for Gfe and Eem were normalized using the elastic modulus calculated using the MTM approach, and the results are plotted in Fig. 5b. Apparently, constrained deformation of the matrix is valid at high volume fractions. As a result, the PAA estimations were higher than the MTM results for lower volume fractions. The agreement between both methods increased with increasing volume fractions for equivalent elastic moduli, simulating constrained deformation of the matrix, as shown in Fig. 5b. For the oblate ellipsoid case of (a = b = 5c), results of the plain strain elastic modulus yields are closer to those from the MTM for high volume fractions. However, for the oblate ellipsoid case of (a = b = 10c), the fully constrained matrix model simulated by Eq. (9) yields results are closer to those obtained from the MTM for high volume fractions. Obviously, increasing the flatness of inhomogenity increases the Poisson contraction effect. Since the ellipsoid with radius ratio a = b = 5c has more concave surfaces than the ellipsoid with a radius ratio of a = b = 10c, a plain strain elastic modulus was sufficient to simulate the Poisson contraction effect for high volume fractions in the oblate ellipsoid case (a = b = 5c). The PAA method estimations, using both the regular- and effective elastic modulus of the matrix mentioned above, also give the property band for the elastic modulus of heterogeneous materials with a continuous matrix and discontinuous inhomogenity. 4.2. Simple cubic (SC) arrangement
Oblate elps. ------ a = b = 5c —— a = b =10c
Eme=EAS 1
Eme=EPS
Eme=Em 0.8
particularly if the rigid inhomogenity has flat surfaces (i.e., platelet-type inhomogenities). In order to clarify this situation, the equivalent stiffness of a matrix under constrained deformation due to the Poisson contraction effect can be modeled using the plain strain elastic modulus:
0
0.2
0.4
0.6
0.8
1
Volume Fraction Fig. 5. (a) The differences between the results of the PAA using the equivalent value of Gfe and results from the MTM approaches for the case of rigid inhomogenity. (b) Comparison of PAA with MTM results. In the PAA calculations, an equivalent value of Gfe and a different equivalent value of Eem , which considers Poisson contraction effects on the matrix were used.
A transversely aligned arrangement is the inhomogenity arrangement in which similarly signed stress fields around the nearest-neighbor inhomogenities overlap [13]. The unit cell for the transversely aligned arrangement is similar to the simple cubic (SC) space lattice. The RVE for the SC unit cell that includes spherical inhomogenity has cubic geometry. Since the heterogeneous domain was assumed to be a 3D infinitely periodic array of transversely aligned particles embedded in the matrix, the stress–strain distributions in the unit cells are symmetric. Due to the shape and symmetric distribution of inhomogenity, the boundary conditions are defined so that the angles between the surfaces of the RVE are preserved; thus, deformed unit cell surfaces remain parallel. Therefore, the outermost layer of the RVE is under a parallel loading condition. The Gf value for the Voronoi cell of an SC unit cell becomes zero (0) since it has parallelepiped walls and the angles between walls non-deform as if they are under parallel loading. Therefore, the Gf values were initially set to the value of the inhomogenity and then reduced linearly down to 0 during the numerical iteration carried out for an SC arrangement. In Fig. 6a, FEM calculations and PAA estimations were given for the different het-
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
a
4.3. The effect of orientation for SC-, BCC-, FCC-, and CPH arrangements
4
PAA FEM HAp /PLLA [12] Glass /Epoxy [13] Al /SiC p [13]
E Composite / EMatrix
3
2
1 0
0.1
0.2
0.3
0.4
0.5
Volume Fraction
b
2943
4
FEM U&L bounds[14]
PAA
E Composite / EMatrix
3
If the Gf values of inhomogenity and of the Voronoi tessellation RVE are calculated with respect to the loading direction, the effects of orientation and arrangement can be considered. In Table 2, the elastic moduli for an Al/SiCp system calculated by the finite element method using different unit cell models and different particulate shapes are compared with the results of PAA. The Gf values at the [0 0 1] and [1 1 0] directions for RVE Voroni cells of SC-, BCC-, and FCC unit cells are given in Table 2. The Gf values for spherical-, cubic-, and cylindrical particulates are also given in Table 2 for the same directions. As seen from Table 2, there was very good agreement between the FEM results and PAA estimations. The Voronoi tessellations for the CPH and FCC arrangements are identical. The Gf values for the Voronoi CPH unit cells are 0.450 in the [0 0 1] direction and 0.491 in the [1 2 0] direction. For an epoxy– 30% glass particulate heterogeneous system (EEpoxy = 3.01 GPa, vEpoxy = 0.394, EGlass = 76 GPa, vGlass = 0.23), the elastic moduli were calculated as 5.810 and 5.681 GPa for the [0 0 1] and [1 2 0] orientations, using an FEM analysis for a solid model with CPH unit cells and spherical particulates [1]. The proposed PAA estimated elastic moduli of 5.470 and 5.344 GPa for the [0 0 1] and [1 2 0] orientations, respectively. The PAA values were less than, but closely approached, the FEM results, which was to be expected since the PAA does not consider the Poisson contraction effect. Although the PAA predictions yielded lower values than those obtained using the FEM analysis, it should be noted that the proposed PAA approach was able to predict higher elastic modulus in the [0 0 1] direction than in the [1 2 0] direction, precisely like those obtained using the FEM analysis. It should also be remembered that Voronoi tessellation of the CPH arrangement has some parallelepiped surfaces retained during deformation such as the SC arrangement mentioned above. If a Gf value of 0.29 was used for the RVE, the elastic modulus was calculated to be 5.621 GPa in the [0 0 1] direction, which was closer to the FEM predicted value of 5.810 GPa. This was because 44% of the total RVE surfaces remain parallelepiped during deformation in the [0 0 1] direction (i.e., (0.44 0 + 0.66 0.45 = 0.29)).
2
4.4. The effect of randomly arranged spherical inhomogenity
1
0
0.1
0.2
0.3
0.4
0.5
Volume Fraction Fig. 6. (a) Normalized elastic moduli of different heterogeneous systems having an SC arrangement and comparison with FEM results. (b) Normalized elastic moduli of a glass–epoxy system with an SC arrangement and comparison with FEM results for iso-displacement and iso-stress loading conditions.
erogeneous systems. As seen in Fig. 6, there is a good correlation between the PAA estimations and FEM results. The expected discrepancy between the FEM results and PAA estimations seen for high volume fractions is due to the constrained deformation of the matrix [13]. The FEM results for the elastic moduli of a glass– epoxy heterogeneous system is also shown in Fig. 6b. For the FEM results [14] given in Fig. 6b, the SC unit cell model is loaded with iso-displacement from one face, as mentioned above, and with iso-stress, in order to reduce the Poisson contraction effect. The PAA curve lies between those lower- and upper bound values as expected (Fig. 6b).
To obtain an effective response for a periodic arrangement, one can use either the smallest RVE or the smallest unit cell with a smaller number of inhomogenities [15]. For a random arrangement, however, the RVE must include a sufficient number of inhomogenities. Thus, the behavior of the RVE represents the overall behavior of the heterogeneous domain. Since the RVE in the proposed PAA calculations includes only a single inhomogenity, the overall effect of arrangement was reflected in the Gf parameter. Obviously, each of the sub-volumes with a single inhomogenity in a random arrangement may have different shapes (irregular or regular). Further, each individual sub-volume may have parallelepiped surfaces retained during deformation as mentioned above for RVE of the CPH arrangement. Hence, an imaginary RVE may be considered to represent the random arrangement case so that it would have an appropriate Gf value that reflects the overall effect of random arrangement. The elastic moduli of a continuous-, isotropic-, and elastic domain with statically isotropic-, randomly arranged spherical inhomogenities was calculated for extreme cases (with voids and with spheres assumed to be rigid), analytically using Torquato’s third order approximation (TOA) [16] and numerically by FEM simulations [17]. A cubic RVE cell containing 30 spheres randomly located with periodic boundary conditions was used in the FEM analysis. Although the Gf can take a value between 0 and 1, a value of
2944
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
Table 2 Elastic moduli of a heterogeneous domain for different arrangements and loading orientations calculated using the finite element method with unit cell models [2] and by the proposed PAA. System
Al/20%SiCp EAl = 67.2 GPa, vAl = 0.35 ESiCp = 429 GPa vSiCp = 0.17
Unit cell
Inhmgn. shape
FEM (GPa)
SC
Sphere Cube Cylinder Sphere Cube Cylinder Sphere Cube Cylinder
96.4 101.0 102.6 89.9 – – 90.0 94.4 94.5
BCC
FCC
E
001
PAA (GPa) E
110
90.2 92.1 91.7 91.1 – – 91.6 93.8 91.4
0.333 was found adequate to reflect the overall behavior of a random arrangement, as determined by trial and error based on the data obtained by the TOA and FEM results shown in Fig. 7. There was very good agreement between the FEM and TOA estimations at low volume fractions. The Gf value of 0.333, which is smaller than 0.5, implies that the parallel loading condition is crucial for the arrangement effects of random distribution. The proposed PAA estimated harder behavior in the high volume fractions for the void case. Matrix regions, which were under localized stress– strain intensities, expanded with increasing volume fractions; thus, the FEM calculations showed softer elastic behavior in the high volume fractions, as expected for the void case. In contrast, the proposed PAA estimated softer behavior for the rigid inhomogenity case with higher volume fractions, as seen in Fig. 7. Obviously, the constrained deformation of the matrix was pronounced with increasing volume fractions in the case of rigid inhomogenity. Therefore, increasing the stiffness of the matrix due to the elastic interaction between inhomogenity and the matrix resulted in harder FEM calculated elastic behavior, as was expected for the rigid inhomogenity case. However, general agreement of the PAA curve with the curves of the analytical TOA and FEM simulations indicates that the proposed PAA can follow the arrangement effects of a random distribution. It should be remembered that more valid Gf values that reflect real situa-
4
PAA TOA [16] FEM [17]
E / E matrix
3
2
Rigid inhomogenity 1 Void
0 0
0.1
0.2
0.3
0.4
E
001
E
110
Gf001 for RVE
Gf110 for RVE
Gf001 for Inhm.
Gf110 for Inhm.
0.500 0.333 0.333 0.500 0.333 0.333 0.500 0.333 0.333
0.500 0.471 0.424 0.500 0.471 0.424 0.500 0.471 0.424
95.3 102.5 102.5 90.3
90.7 91.5 92.7 91.1
0.000
0.471
0.523
0.422
90.7 95.6 95.6
90.5 91.2 92.4
0.471
0.500
tions could be determined by comparison with experimental results. 5. Concluding remarks
1. The elastic modulus of heterogeneous domains with a continuous matrix and discontinuous inhomogenities was estimated using a numerical calculation method. This method can consider the shape, orientation and arrangement effects of inhomogenity via a simple parameter, the geometrical factor (Gf). Gf values vary from 0 to 1 and indicate stress–strain distribution conditions of the constituents in a heterogeneous domain. The Gf value approaches 0 and 1 as the stress–strain distributions approach the parallel- and serial loading conditions, respectively. 2. The shape and orientation effects of the inhomogenity were considered with the Gf value of the inhomogenity. Arrangement effect of the inhomogenity was taken into account by using the Gf value belonging to a representative volume element (RVE) in which the inhomogenity was embedded. 3. In general, the Gf was determined using the ratios (i) between the projection area in the loading direction and the surface area of the inhomogenity and (ii) between the projection area in the loading direction and the surface area of RVE. If the RVE, with the inhomogenity, had parallelepiped surfaces during deformation, such as that in an SC arrangement, Gf values were chosen to be 0 (i.e., the parallel loading condition is dominant in the heterogeneous domain). 4. There was very good agreement between the estimations of the proposed numerical approach (PAA) and results of analytical calculations obtained from MTM and TOA. There was also very good agreement between the estimations of the PAA and results of FEM simulations taken from the literature for SC-, BCC-, FCC-, and CPH packing systems and for unit cells belonging to randomly distributed inhomogenities. 5. The PAA can be used in systematic- and parametric investigations into the effects of microstructural geometry on the elastic modulus of a heterogeneous domain. The numerical method proposed can provide quick information and is easily implemented by engineers. 6. The PAA still needs to consider Poisson contraction effects appearing at high volume fractions. However, the PAA can provide property bands for elastic moduli of heterogeneous domains (for a matrix which is under plane strain condition and for a matrix which is free from Poisson contraction).
0.5
Volume Fraction Fig. 7. Normalized elastic moduli of a heterogeneous domain with randomly arranged voids and rigid inhomogenities and comparison with both FEM and TOA results.
Acknowledgment The author gratefully acknowledges financial support of the Technical University of Istanbul, BAP Program.
S. Yilmaz / Materials and Design 30 (2009) 2938–2945
References [1] Li S, Wongsto A. Unit cells for micromechanical analyses of particle-reinforced heterogeneous. Mech Mater 2004;26:543–72. [2] Weissenbek E, Bohm HJ, Rammerstorfer FG. Micromechanical investigation of arrangement effects in particle reinforced metal matrix heterogeneous. Comput Mater Sci 1994;3:263–78. [3] Levy A, Papazian JM. Elasto-plastic finite element analysis of short fiber reinforced Sic/Al heterogeneous: effect of thermal treatment. Acta Metall Mater 1991;39:2255–66. [4] Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials. The Netherlands: Elsevier; 1999. [5] Mura T. Micromechanics of defect in solids. The Netherlands: Kluwer; 1998. [6] Eshelby CC. The determination of the elastic field of an ellipsoidal inclusion, and related field. Proc R Soc Lond 1957;A241:376–96. [7] Mori T, Tanaka K. Average stress in the matrix and average elastic energy of materials with misfitting inclusions. Acta Metall Mater 1973;21:571–4. [8] Hashin Z, Shtrikman S. A variational approach to the theory of elastic behavior of multiphase materials. J Mech Phys Solid 1963;11:127–40.
2945
[9] Voight W. Lehrbuch Der Kristallphysik. Berlin-Leipzig: Teubner Verlag; 1928. [10] Reuss A. Berechung der fliessgrenzen von mischkristallen aurf grund der plastiizitatsbedingug fur einkristalle. Z Angew Math Mech 1929;9:49–58. [11] Ahuja N, Schachter BJ. Pattern models. New York: Wiley; 1983. [12] Balac I, Milovancevic M, Tang C, Uskokovic PS, Uskokovic DP. Estimation of elastic properties of a particulate polymer composite using a face-centered cubic FE model. Mater Lett 2004;58:2437–41. [13] Yilmaz S, Aran A. Finite element analysis of deformation behavior in ductile matrix containing hard particles. Mater Sci Tech 1998;14:1154–62. [14] Sun C, Saffari P, Ranade R, Sadeghipour K, Baran G. Finite element analysis of elastic property bounds of a composite with randomly distributed particles. Composites 2007;A38:80–6. [15] Jiang M, Jasiuk I, Ostoja-Starzewski M. Apparent elastic and elastoplastic behavior of periodic composites. Int J Solids Struct 2002;39:199–212. [16] Torquato S. Effective stiffness tensor of heterogeneous media. II: Application to isotropic dispersions. J Mech Phys Solids 1998;46:1411–40. [17] Segurado J, Llorca J. A numerical approximation to the elastic properties of sphere reinforced heterogeneous. J Mech Phys Solids 2002;50:2107–21.