Micromechanics of stress transfer through the interphase in fiber-reinforced composites

Mechanics of Materials 89 (2015) 190–201

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Micromechanics of stress transfer through the interphase in fiber-reinforced composites Priyank Upadhyaya a, S. Kumar a,b,⇑ a Institute Center for Energy (iEnergy), Department of Mechanical and Materials Engineering, Masdar Institute of Science and Technology, PO Box 54224, Abu Dhabi, United Arab Emirates b Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139-4307, United States

a r t i c l e

i n f o

Article history: Received 8 February 2015 Received in revised form 13 May 2015 Available online 19 June 2015 Keywords: Polymer matrix composites (PMCs) Micromechanics Interphase Interface Stress transfer Eigenfunction expansion method

a b s t r a c t A micromechanical model to predict the interphasial/interfacial stress transfer in a three-phase fiber-reinforced composite is presented. The axisymmetric system consists of a fiber embedded in a compliant matrix having an interphase between them. Each constituent of the composite is regarded as a linear elastic continuum. The matrix is treated as an isotropic material while the fiber and interphase are considered as a transversely isotropic material. Traction-free boundary conditions are strictly enforced. It is assumed that the interfaces are perfect and strong. A pair of uncoupled governing partial differential equations is obtained in terms of unknown displacements. Furthermore, assuming that the Eigenvalues exist for this system of equations, Eigenfunction expansion method is employed to derive an exact solution in terms of the Bessel functions. Analytical solutions are obtained for free boundary conditions at the external surface of the matrix cylinder to model a single fiber pull-out problem, and for fixed boundary conditions to approximately model a hexagonal array of fibers in the matrix material. This formulation provides an analytical framework for the analysis of interphasial and interfacial stresses as well as displacements in the entire 3D axisymmetric system. Finite element (FE) analysis was also performed to simulate stress transfer from the fiber to the matrix through the interphase. Analytically obtained stress fields are verified with FE results. Shear and radial interphasial stresses provide insight into the design of engineered interfaces/interphases. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Fiber-reinforced composites (FRCs) are widely used in advanced engineering applications due to their low specific weight and superior thermo-mechanical stability. Furthermore, the bi-material nature of FRCs can be exploited to advantageously tailor the properties according ⇑ Corresponding author at: Institute Center for Energy (iEnergy), Department of Mechanical and Materials Engineering, Masdar Institute of Science and Technology, PO Box 54224, Abu Dhabi, UAE. Tel.: +971 2 810 9239; fax: +971 2 698 8121. E-mail addresses: [email protected], [email protected], [email protected] (S. Kumar). http://dx.doi.org/10.1016/j.mechmat.2015.06.012 0167-6636/Ó 2015 Elsevier Ltd. All rights reserved.

to application requirements. One of the most important phenomena in FRCs is the stress transfer between the fiber and the matrix across the interphase/interface. When composites are subjected to various loading conditions, the efficiency of load transfer across the interface plays an important role in overall performance of the composites (Kim and Mai, 1998). Several researchers have investigated the influence of interfacial strength and the quality of adhesion on overall performance of composites. As a result of physical and chemical interactions between the fiber and the matrix, a nano-meter length scale thin layer (Wu et al., 2014) forms between them during processing of composites. This is referred to as an interphase. The

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

chemical composition and thermo-mechanical properties of the interphase differ from both the reinforcing fiber and the matrix material (Sottos et al., 1992; Duek, 1986). The formation of interphase depends on multiple factors such as specimen geometry, fiber volume fraction and the bonding mechanism between the matrix and the reinforcement (Atkins, 1975; Hughes, 1991; Naslain, 1998; Gao and Mäder, 2002; Liu et al., 2008). Within the interphase, the elastic properties vary in the radial direction (Sottos et al., 1992). To characterize the influence of interface on the load bearing capacity of composites, a number of experimental techniques, such as pull-out, push-out, microbond or fragmentation tests are commonly performed. Among these popular techniques, single fiber pull-out test is widely used mainly due to its simple geometry. In spite of its simple geometry, a complete closed form solution of a fiber pull-out problem is hopelessly complicated. In the past few decades, a number of analytical solutions have been proposed in order to better understand the stress transfer mechanism across the interfaces between the fiber and the matrix. Generally, these solutions are obtained by adopting a few simplistic assumptions. One of the commonly employed assumptions is the state of plane-strain (Cox, 1952) in the assembly in order to relate radial and shear stresses in the equilibrium equations. For example, Hsueh (1990, 1991,, 1992) published a series of papers on different aspects of fiber pull-out/push-out problems, adopting plane-strain conditions. However, the plane-strain assumption is not a good approximation for perfectly bonded interfaces where the load transfer length is very small (only a few fiber radii) (Hsueh, 1991). Another such simplification is to ignore the radial dependence of axial stresses, either in the fiber or matrix, or in both (Hsueh, 1991; Cox, 1952), which is justifiable only when embedded fibers are sufficiently long. Owing to the mathematically involved nature of this boundary value problem, a few investigators have chosen to impose boundary conditions in an approximate sense. For instance, Hsueh (1990, 1991) developed an analytical solution based on shear-lag theory. In this model, shear stress in the matrix region was assumed to be inversely proportional to the radial distance. Consequently, a stress free boundary condition was satisfied only in a limited sense which is appropriate only for problems where the outer radius of the matrix is considerably large compared to the fiber radius. Nairn (1997) showed that the shear-lag analysis is best suited for composites with high fiber to matrix stiffness ratio. Moreover, a shear-lag analysis predicts inaccurate results for composites with low fiber volume fractions. However, the total energy prediction through shear-lag analysis is reasonably accurate and therefore shear-lag theory has been utilized to study fracture mechanics of embedded fibers (Nairn, 1992). To develop tractable models, many researchers have modeled the interphase region as a homogeneous material (Christensen and Lo, 1979; Hashin, 1990; Hayes et al., 2001; Qiu and Weng, 1991; Rjafiallah et al., 2010; Tsai et al., 1990). However, a few studies considered the inhomogeneous nature of interphase adopting a stair-case variation of material properties across the thickness of the

191

interphase layer (Jiang et al., 2008; Mogilevskaya, 2001; Wang et al., 2006). Alternatively, a few investigators proposed an effective interphase model (EIM) and uniform replacement model (URM) to replace the fiber and the surrounding interphase by an effective homogeneous fiber in order to convert a three-phase composite into a two-phase composite (see, for e.g. Shen and Li, 2005). For mathematical convenience and to better describe the variation of properties within the interphase region, several researchers treated the interphase as an inhomogeneous material by smoothly varying the material properties as a function of radius. Usually in such models, the material properties are graded by adopting an empirical law (Huang and Young, 1996; Jayaraman and Reifsnider, 1992, 1993; Kiritsi and Anifantis, 2001; Low et al., 1994; Romanowicz, 2010; Shen and Li, 2003). Since, it is rather difficult to arrive at an explicit closed form 3D solution for such problems, many researchers have adopted numerical methods, particularly FE analysis (Kovalev et al., 1998; Qing, 2013; Kiritsi and Anifantis, 2001; Needleman et al., 2010). Most of the analytical studies that exist in the literature do not account for the anisotropy of fibers. Nevertheless, Wu et al. (2000) presented a generic solution that is valid for both isotropic (ISO) and transversely isotropic (TISO) materials for a two-phase fiber-reinforced composite. In this study, following Wu et al. (2000), a 3D axisymmetric theoretical model is presented to predict the stress transfer in three-phase composite material. A thin explicit interphase having transversely isotropic material behavior is introduced to model the interphase region. The main goal of this study is to develop a 3D axisymmetric micromechanical model without simplifying assumption on the stress- and/or strain-state in order to characterize the stress transfer at the interfaces in a three phase composite while treating the interfaces to be strong (both stresses and displacements are continuous).

2. Problem statement The first step towards analyzing fiber-reinforced composites is to choose a representative volume element (RVE) based on the fiber packing arrangement. The most commonly preferred arrangements are square packed and hexagonal packed arrays of fibers in the matrix. In the current study, the hexagonal fiber packing is idealized as a coaxial cylinder composed of a cylindrical fiber surrounded by the matrix material. Within this RVE, an explicit interphase is considered. A schematic of such an RVE is shown in Fig. 1. The task is to develop a tractable micromechanical model so as to obtain the stress fields in a three-phase system under mechanical loading conditions. The coordinate system ðr; h; zÞ with origin (denoted by letter O) at center of the three-phase system is shown in Fig. 1(a). L is the half-length of the RVE. r f and r m are the radii of the fiber and the matrix outer surface respectively. r0 is the stress applied at the top end of the fiber. Fig. 1(b) shows the simplified model with boundary conditions, in which the outer surface of the matrix is traction-free which is an idealization of a single fiber pull-out problem. Fig. 1(c) shows the

192

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

(a)

(b)

(c)

Fig. 1. Three-phase RVE with boundary and loading conditions: (a) axisymmetric geometric model, with applied traction r0 , (b) traction-free outer matrix surface, (c) restrained outer matrix surface.

geometric model in which the outer surface of matrix is fixed so as to approximately model a hexagonal array of fibers in the matrix material. Note that the three phases are constrained in the z-direction at z ¼ 0 and the matrix alone is constrained in z-direction at the top end as shown in Fig. 1(b) and (c).

3. Analytical formulation 3.1. Model geometry and materials As discussed in Section 2, a 3D axisymmetric RVE consists of three material phases. A geometric model of a three-phase axisymmetric RVE is shown in Fig. 1. In the current analysis, the fiber radius rf is 5 lm and the RVE radius is 12:5 lm. In between the fiber and the matrix a thin layer with thickness t ¼ 0:1rf is introduced to model an explicit interphase region. The fiber half-length L is taken to be 10r f (see Fig. 1). The material used in this study is a carbon fiber reinforced epoxy composite. In reality, the interlayer between the fiber and the matrix is inhomogeneous with a smooth transition in mechanical properties. At this scale, it becomes extremely difficult to conduct

measurements, although properties of interphases can be determined through molecular dynamics simulations. However, researchers have demonstrated the use of the SIEM (Speckle Interferometry with Electron Microscopy) technique to measure the elastic properties of interphases (Wang and Chiang, 1996). Elastic properties of the fiber, matrix and the interphase are given in Table 1. Most of the analytical studies treat both fiber and the matrix as isotropic materials. However, experimental observations indicate that the fiber exhibits transverse isotropy and therefore, the fiber is regarded as a transversely isotropic material in this study. As shown by Wu et al. (2000), thickness of interphase region ranges from 100 to 200 nm for carbon fiber-epoxy composites. However, thickness and properties of the interphase depend upon local changes of chemical and physical processes such as crystallization and cross-linking in the vicinity of the fiber. For instance, for S-glass epoxy composites, interphase thickness varies from 10 to 1000 nm (Tsai et al., 1990). An interphase thickness of 500 nm is considered in this study to better visualize the results. Since, the micromechanical model presented in this study is based on an assumption that all constituents of the composite are continua, the choice of interphase thickness does not

Table 1 Mechanical properties of fiber, matrix and the interphase (Wu et al., 2000). Materials

Er (GPa)

Eh (GPa)

Ez (GPa)

Grz (GPa)

mrh

mrz

Fiber Matrix Interphase

20 2.7 11.0

20 2.7 11.0

210 2.7 105

23 1.0 12.0

0.25 0.35 0.3

0.025 0.35 0.16

193

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

matter as the formulation does not account for any bonded/non-bonded atomic interactions which usually occur at nanometer scale. The interphase is either compliant or brittle depending upon the composite system. For S-glass epoxy and graphite-epoxy composites, the interphase is much compliant than the matrix (Tsai et al., 1990). On the other hand, the interphase is stiffer than the matrix material for some uncoated carbon-epoxy composites or for GF/polyester composites (Cech et al., 2013). In this study, carbon-epoxy composite is considered and therefore the interphase is assumed to be stiffer than the matrix. Therefore, average properties are assigned to the interphase region to model the stiff response of the interphase assuming a single step variation in elastic properties across the thickness of the interphase. However, the mathematical formulation presented here can be easily extended to model a multi-step or smooth variation in properties across the thickness of the interphase. Micrographic image analysis of the interphase region suggests that the interphase has a gradient in material properties such that it matches the properties of the fiber and the matrix at the boundaries (Cech et al., 2013). The relative stiffness of the interphase region also depends on the thickness of the interphase layer. A thicker interlayer translates into a stiffer interphase whereas a thin interlayer exhibits a softer behavior (Tsai et al., 1990). In the current study, the case of a stiffer interlayer is studied.

The components of the strain vector in a cylindrical coordinate system are expressed in terms of displacements u and w as

fg ¼

"

ð1aÞ

1 @ @ rz ðr srz Þ þ ¼0 r @r @z

ð1bÞ

The linear elastic constitutive relationship in cylindrical coordinates for an axisymmetric system is given by

frg ¼ ½Cfg where,

ð2Þ

stress

and T



strain

T

vectors

are

frg ¼ frr ; rh ; rz ; srz g ; fg ¼ r ; h ; z ; crz , and stiffness matrices ½C for TISO and ISO materials are given by

2 C TISO

6 6 C 12 6 ¼6 6 C 13 4 2

C ISO

C 11

C 12

C 13

C 11

C 13

0

3

7 0 7 7 7 0 7 5

C 13

C 33

0

0

C 11

C 12

C 12

0

C 11

C 12

0

C 12

C 11

0

0

0

ðC 11  C 12 Þ=2

0

0

@z

0

1=r

@z

@r

T 

u



ð4Þ

w

ð5Þ where @ r ; @ z are partial derivatives with respect to r and z respectively and C 55 ¼ C 13 þ C 44 is a material constant. After a few algebraic manipulations of operators on Eq. (5), a pair of following fourth order governing differential equations in u and w can be derived.

ðr21 þ n1 @ 2z Þðr21 þ n2 @ 2z Þu ¼ 0 ðr22 þ n1 @ 2z Þðr22 þ n2 @ 2z Þw ¼ 0



r21 ¼ @ 2r þ



@r 1  r r2

r22 ¼ @ 2r þ

@r r

 ð6aÞ

 ð6bÞ

where n1;2 is given by

1 ¼ 2 ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! C j  j2  4 33 here; C 11

C 33 C 44 ðC 13 þ C 44 Þ2 þ  C 44 C 11 C 11 C 44

j ð7Þ

Closely observing the differential operators r21 and r22 , it can be seen that @ z r21 @ r ¼ @ r r22 @ z . It follows from this relation that the displacements u and w are related through a common function / such that u ¼ @ r / and w ¼ k@ z /. Substituting stresses in terms of the function / in equilibrium equation (Eq. (1)) gives

 @ r C 11 r22 þ ðC 44 þ kC 55 Þ@ 2z / ¼ 0

ð8aÞ

 @ z ðC 55 þ kC 44 Þr22 þ kC 33 @ 2z / ¼ 0

ð8bÞ

here, the constant k can be determined by combining Eqs. (6) and (8) as

ki ¼

C 11 ni  C 44 C 55 ni ¼ C 55 C 33  C 44 ni

ð9Þ

The Eigenvalue solution for the function / can be written as

0 6 6 C 12 6 ¼6 6 C 12 4

0

#     C 11 @ 2r þ @rr  r12 þ C 44 @ 2z C 55 @ r @ z u 0 ¼   2 @r 2 1 w 0 C 55 @ r þ r @ z C 44 @ r þ r þ C 33 @ z

n1;2

@ rr rr  rh @ srz þ þ ¼0 @r r @z

@r

Stress components can be expressed in terms of displacements by substituting Eqs. (3) and (4) in Eq. (2). Using resulting expressions for stresses, equilibrium equations (Eq. (1)) can be rewritten as

3.2. 3D axisymmetric solution To maintain the axial symmetry of the RVE, displacement (v) and variation (@ h ) along the hoop direction have to be zero. Therefore, non-zero stresses and displacements are rr ; rh ; rz ; srz and u; w respectively. As a consequence of these model assumptions, the non-trivial equilibrium equations in cylindrical coordinates are



/i ¼ ½a1 sinhðzi Þ þ a2 coshðzi Þ½b1 J 01 þ b2 Y 01 þ b3 J 02 þ b4 Y 02 

C 44

ð10Þ

3 7 7 7 7 7 5

ð3Þ

where J pq and Y pq are the Bessel functions of first (J p ) and pffiffiffiffiffiffiffiffi second (Y p ) kind evaluated at rq ¼ nq ki r respectively, pffiffiffiffi and zi ¼ ki z, for each Eigenvalue ki . The Eigenvalue solution for displacements u and w can be derived from / as

194

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

pffiffiffiffi ui ¼  ki ½a1 sinhðzi Þ h pffiffiffiffiffi pffiffiffiffiffii þa2 coshðzi Þ ðb1 J 11 þ b2 Y 11 Þ n1 þ ðb3 J 12 þ b4 Y 12 Þ n2 ð11aÞ

pffiffiffiffi wi ¼ ki ½a1 coshðzi Þ þ a2 sinhðzi Þ½ðb1 J 01 þ b2 Y 01 Þk1 þðb3 J 02 þ b4 Y 02 Þk2 

ð11bÞ

where, a1 ; a2 ; b1 ; b2 ; b3 and b4 are integration constants. For zero Eigenvalue (k ¼ 0) solution, limit case of Eq. (6b) can be solved to get w0 . Then, plugging w0 in Eq. (1a), u0 can be determined. Two linearly independent solutions for u0 and w0 are

u0 ¼ b30 r þ

b40 C 55  b20 r log r r 2C 11

w0 ¼ ða10 þ zÞðb10 þ b20 log rÞ

ð12aÞ ð12bÞ

Combining the non-zero and zero Eigenvalue solutions, complete displacement solutions can be written as

u¼

1 X

1 X

i¼1

i¼1

ui þ u0 and w ¼

wi þ w0

ð13Þ

2

ðr22 þ @ 2z Þ u ¼ 0 and ðr22 þ @ 2z Þ w ¼ 0

ð14Þ

Corresponding / solution for isotropic case can be written as h i pffiffiffiffi /i ¼ ½a1 sinhðzi Þ þ a2 coshðzi Þ b1 J 0 þ b2 Y 0 þ ki rðb3 J 1 þ b4 Y 1 Þ

ð15Þ With k ¼ 1, the general Eigenvalue solution for displacements u and w can be derived as ui ¼

The displacement expressions presented above involve many unknown integration constants. Displacement and stress boundary conditions are imposed to determine these constants. For clarity, all the boundary conditions are grouped into four categories. Continuity of displacements (u and w) and stresses (rr and srz ) is enforced at the interfaces. 3.3.1. Continuity at the fiber-interphase (FI) interface

u f ðr f ; zÞ ¼ ui ðr f ; zÞ;

w f ðr f ; zÞ ¼ wi ðrf ; zÞ ð0 6 z < LÞ ð17aÞ

rrf ðrf ; zÞ ¼ rir ðrf ; zÞ; srzf ðrf ; zÞ ¼ sirz ðrf ; zÞ ð0 6 z < LÞ ð17bÞ 3.3.2. Continuity at the interphase-matrix (IM) interface

ui ðr i ; zÞ ¼ um ðri ; zÞ;

wi ðri ; zÞ ¼ wm ðr i ; zÞ ð0 6 z < LÞ ð18aÞ

For isotropic materials, material constants k1 ¼ k2 ¼ 1 and similarly n1 ¼ n2 ¼ 1. Therefore, the solution for non-zero Eigenvalues (ki – 0) obtained for TISO material cannot be used for ISO materials. However, it can be proved that the zero Eigenvalue solution given in Eq. (12) for TISO material can be used for ISO material as well. To determine the nonzero Eigenvalue solution, the uncoupled equations for ISO material take the following biharmonic form 2

3.3. Boundary conditions

h i pffiffiffiffi pffiffiffiffi ki ½a1 sinhðzi Þ þ a2 coshðzi Þ b1 J 1  b2 Y 1 þ ki rðb3 J0 þ b4 Y0Þ

rir ðri ; zÞ ¼ rmr ðri ; zÞ; sirz ðri ; zÞ ¼ smrz ðri ; zÞ ð0 6 z < LÞ ð18bÞ 3.3.3. Outer surface conditions Both stress free (Case-I) and radially constrained (Case-II) boundary conditions on the outer surface of the matrix are explored in the current study. Case-I represents the pull out of single fiber embedded in the matrix whereas Case-II approximately models the pull-out of fibers from hexagonally arranged fibers within the matrix. In the latter case, to maintain the periodicity of the fiber arrangement the outer surface is not allowed to displace radially. In both cases, the shear stress srz is assumed to be zero at the outer surface of the matrix.

Case  I :

rmrr ðrm ; zÞ ¼ 0 or

Case  II :

um ðr m ; zÞ ¼ 0 ð0 6 z < LÞ

smrz ðrm ; zÞ ¼ 0 ð0 6 z < LÞ

ð19bÞ

ð16aÞ

wi ¼

h i pffiffiffiffi pffiffiffiffi ki ½a1 coshðzi Þ þ a2 sinhðzi Þ b1 J0 þ b2 Y 0 þ ki rðb3 J1 þ b4 Y 1 Þ

3.3.4. Axial end conditions

ð16bÞ

w f ðr; 0Þ ¼ 0;

srzf ðr; 0Þ ¼ 0 ð0 < r < rf Þ

ð20aÞ

The complete solution will be obtained by combining the non-zero Eigenvalue solution and the zero Eigenvalue solution. Integration constants a1 ; a2 ; a3 ; a4 ; b1 ; b2 ; b3 ; b4 ; a10 ; b10 ; b20 ; b30 and b40 can be determined using appropriate boundary conditions for the boundary value problem which are discussed in next section. Generic solutions for both transversely isotropic and isotropic material systems are formulated. A complete solution is obtained by utilizing the respective formulation based on the material symmetry of each constituent within the three-phase composite system.

wi ðr; 0Þ ¼ 0;

sirz ðr; 0Þ ¼ 0 ðrf < r < ri Þ

ð20bÞ

wm ðr; 0Þ ¼ 0;

smrz ðr; 0Þ ¼ 0 ðri < r < rm Þ

ð20cÞ

sirz ðr; LÞ ¼ 0 ðrf < r < ri Þ

ð20dÞ

wm ðr; LÞ ¼ 0;

ð20eÞ

Z 0

rf

smrz ðr; LÞ ¼ 0 ðri < r < rm Þ

rzf ðr; LÞ ¼ r0 pr2f ; srzf ðr; LÞ ¼ 0 ð0 < r < rf Þ

ð20fÞ

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

In addition to aforementioned boundary conditions, finite displacements and stresses at the center of the fiber are imposed to avoid singular stresses and displacements. All homogeneous boundary conditions are satisfied for non-zero and zero Eigenvalue solutions separately. Eigenvalues for the system can be obtained through the end conditions in the matrix region using Eqs. (20c) and (20e) as

2 ip ki ¼  L

ð21Þ

Applying other boundary conditions, a system of linear equations can be formed which then can be solved to determine the unknown integration constants. There are totally six integration constants associated with each nonzero Eigenvalue solution. All homogeneous boundary conditions are satisfied separately for each zero/nonzero Eigenvalue solution. Even with an increase in number of terms in the series solution, the system remains deterministic with adequate number of equations to solve for unknown integration constants. 4. Computational model To verify the predictions of the analytical model, a 3D axisymmetric FE analysis of stress transfer in a three phase composite was performed using a commercial FE code Abaqus FEA version 6.12. The geometric model and material constants used in the FE analysis are identical to those

195

of analytical study. A mesh convergence study was conducted to see the effect of element size on the stress distribution in the interphase. Stress field at singular point (at the top where the fiber and the interphase are bonded) increased with increase in mesh density. However, the shear stress at the mid-surface of the interlayer does not vary noticeably with change in element size as the stresses at the midsurface of interlayer are finite. Therefore, the stresses at the mid-surface of the interphase were used to ascertain the accuracy of the analytical model. The focus of the study is not to capture the bi-material singular stress field that exists at the free surface of assembly. FE mesh shown in Fig. 2 consists of 1899 bilinear axisymmetric (CAX4R) elements with two degrees of freedom (displacements in r and z-directions) per node. Identical to the analytical solution, bottom of the RVE (z ¼ 0) is constrained in z-direction. On the top surface (z ¼ L), the matrix is constrained in z-direction while the interphase is free to move. A static stress analysis was conducted under a uniform tensile stress r0 at the free fiber end as shown in Fig. 2. The outer surface of the matrix is either constrained in the radial direction or left free, as discussed in Section 3.3.3. Results from numerical and analytical studies are discussed in the next section. 5. Results and discussion With known material and geometric parameters, a linear system of equations was solved to obtain the

Fig. 2. Finite element model with appropriate boundary conditions.

r0 is the applied traction. (Case-I).

196

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

Eigenvalue solutions for unknown displacements u and w. Since the complete displacement solution is obtained by summing up all the nonzero Eigenvalue solutions to the zero Eigenvalue solution, the value of displacement increases with increase in number of Eigenvalues and eventually stabilizes. To minimize the computational demand, only the first twenty Eigenvalue solutions were determined. Therefore, the contribution of first few Eigenvalues is significantly higher than the later ones. This truncation in series expansion also results in a smooth transition of shear stress at the interface at the loading end (Wu et al., 2000), where singular boundary conditions exist. Analyses were conducted for two cases of boundary

conditions at the outer surface (r ¼ r m ) in matrix region. Case-I (rr ¼ 0) represents the pull out of single fiber embedded in matrix whereas Case-II (u ¼ 0) approximately models the pull-out of fibers from hexagonally arranged fibers within matrix. No significant change in stress distribution around the interphase region between Case-I and Case-II was observed and therefore only Case-I results are presented here for brevity. Fig. 3 depicts the shear stress srz distribution for Case-I along z at three different radial locations (r ¼ r f ; rmid and ri ) and along r at different z locations, normalized by applied traction r0 . Here, r ¼ rf represents the fiber-interphase (FI) interface, r ¼ ri represents the interphase-matrix (IM)

(a)

(b)

Fig. 3. Variation of normalized shear stress srrz0 for Case-I: (a) Along the axial direction (z) in the vicinity of loaded fiber end. (b) Along the axial direction (z) for the entire fiber length.

Fig. 4. Surface plot showing the variation of normalized shear stress srrz0 in the vicinity of the loaded fiber end for Case-I.

197

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

(a)

(b)

Fig. 5. Variation of normalized radial stress rr0r for Case-I: (a) Along the axial direction (z) in the vicinity of loaded fiber end. (b) Along the axial direction (z) for the entire fiber length.

(a)

(b)

Fig. 6. Variation of normalized axial stress rr0z for Case-I: (a) Along the axial direction (z) in the vicinity of loaded fiber end. (b) Along the axial direction (z) for the entire fiber length.

interface and r ¼ rmid represents the radial mid-surface of the interphase. Fig. 3(a) shows the axial variation of srz in the vicinity of the loaded end of the fiber. It is very clear

from Fig. 3 that the shear stress srz peaks towards the loaded end of the fiber (Lz 61) and decays rapidly as we r f

move towards the fixed fiber end. Within

Lz rf

¼ 1, shear

198

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

stress reduces down to 10% of its peak value. Radial variation of shear stress at six different locations (Lz ¼ 0; 0:1; 0:2 . . . 0:5) is shown in Fig. 3(b). It is evident r f

from Fig. 3(b) that shear stress is concentrated around the interphase region and major portion of the fiber and matrix have negligible shear stress. Also, it is worth noting in Fig. 3 that the shear stress continuity at both interfaces is perfectly enforced. Shear stress in rz-plane is shown as a 3D surface in Fig. 4, in the vicinity of loaded fiber end (0 < Lz 6 1) for better visualization. Such a concentrated r f

stress distribution close to stress singular point makes the interphase region very susceptible to mixed mode fracture (mode I and II), depending upon the elastic mismatch and the fiber length (long/short fiber). This underlines the need for a better interphase design towards improving overall performance of composite materials. Closely observing the slope of shear stress in different regions (see, Fig. 3(b)), it is seen that the slope in the interphase region is slightly lower than that in the fiber region which in effect mitigates the shear stress. This is because of reduction in mismatch of elastic properties due to the presence of an interphase between the stiff fiber and compliant matrix. This idea can be extended to continuously grade the properties within the interphase, in order to minimize the shear stress where actual transition in elastic

properties from the fiber to matrix is achieved in infinitely many number of steps. Axial variation of radial stress rr at different radial locations (r ¼ rf ; rmid and ri ) is presented in Fig. 5(a) for Case-I. Maximum radial stress along z occurs right at the loaded end of the fiber and diminishes rapidly (within half fiber radius). The magnitude of radial stress is relatively lower compared to the magnitude of shear stress. However, radial tension at FI interface could initiate debonding (mode-I) at the loaded fiber end if the radial stress exceeds the interfacial tensile strength. Consequently, such micro-scale damages can grow rapidly leading to macro scale failure of the component. In Fig. 5(b), radial variation of radial stress rr shows that the FI interface is under tensile stress which makes it a relatively weaker bonding surface in the RVE. The radial stress continuity can be clearly seen from Fig. 5(b). Radial stress profiles are very similar for Case-I and Case-II. Only a small offset in Case-I is observed which forces the radial stress to vanish at the outer matrix surface (r ¼ r m ) as prescribed in Case-I. The interphasial shear, tensile and compressive strengths for carbon fiber reinforced epoxy composites as reported by Wang et al. (2011) are equal to 50 MPa. Analytical model presented in this study gives maximum normalized shear stress of 0.175 and maximum normalized tensile stress of 0.1 in the interphase

(a)

(b)

(c)

Fig. 7. Comparison of normalized stress fields at

Lz rf

¼ 0:5 for Case-I: (a) Shear stress srrz0 . (b) Radial stress rrrr0 . (c) Axial stress rrzz0 .

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

layer. Therefore, for a carbon-fiber reinforced epoxy composites shear induced damage would occur first within the interphase when the applied tensile stress on the fiber reaches 285 MPa. Axial and radial variation of axial stress

199

is presented in Fig. 6. Similar to the radial stress, rzz reaches its peak value at the loaded fiber end and diminishes steeply. In the radial direction (Fig. 6(b)), rzz is almost constant and is approximately equal to  0:9r0 (lower

(a)

(b)

33 Fig. 8. Variation of normalized radial displacement uC r0 L for Case-I: (a) Along the axial direction (z) in the vicinity of loaded fiber end. (b) Along the axial direction (z) for the entire fiber length.

(a)

(b)

33 Fig. 9. Variation of normalized axial displacement wC r0 L for Case-I: (a) Along the axial direction (z) in the vicinity of loaded fiber end. (b) Along the axial direction (z) for the entire fiber length.

200

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201

than the applied stress r0 ) in the major part of the fiber and steeply rises close to the FI interface beyond the applied stress ( 1:5r0 ). rzz distribution follows the global force balance i.e. at the loaded end the average stress in the fiber is still r0 . Analytically estimated stresses as a function

 of dimensionless radial distance rr at a cross section f

where

Lz rf

¼ 0:5 are compared with FE predictions in

Fig. 7. In Fig. 7(a), the shear stress results from the present model match reasonably well with FE calculations. In Fig. 7(b) radial stress comparison is shown. In the interphase and matrix region a perfect match between analytical and FE predictions is observed, unlike in the fiber region. The possible explanation for this deviation could be the manner in which force boundary condition in the fiber region is enforced. The analytical solution was obtained by satisfying axial stress condition in a global sense which ensures the equilibrium but fails to capture the point-wise variation of stresses in the fiber region. Similarly, axial stress obtained from the analytical model was compared with that of the computational model in Fig. 7(c). Analytically determined axial stress solution agrees reasonably well with FE predictions. Here again, the effect of satisfying boundary forces macroscopically is clearly seen. For completeness, distribution of

 33 non-dimensional radial displacement uC for Case-I is r0 L shown in Fig. 8. In the absence of any radial forces, radial displacements are very small in the fiber and interphase regions compared to the displacements seen in matrix part. A perfect continuity of radial displacement at both interfaces is evident from Fig. 8(b). Displacement profiles for Case-I and II are also very similar to each other. Radial displacement u vary linearly along r-direction in the fiber and in the interphase region, and it increases steeply after IM interface and reaches its maximum value in the matrix region. Beyond the peak location, radial dis  33 placements plateaux uC r0 L ¼ 0:25 . Finally, distribution

 33 for case-I is of dimensionless axial displacement wC r0 L presented in Fig. 9. Along the axial direction z, with peak value at the loaded fiber end (z ¼ L), w reaches 0 at z ¼ 0 to satisfy the constrained boundary conditions. From the radial variation of w shown in Fig. 9(b) it is clear that the displacement is constant throughout the fiber and the interphase but reduces significantly in the matrix region. 6. Conclusion A theoretical framework to analyze the micromechanics of stress transfer in a three-phase axisymmetric composite system based on an Eigenfunction expansion method is presented. The interface is modeled as a finite width zone with distinct material properties. In contrast to many solutions available in literature where fiber is modeled as an isotropic material, the proposed model takes into account the transverse symmetry of the fiber. In this displacement based formulation, a pair of uncoupled fourth order differential equations in terms of non-zero displacements (u and w) were obtained. Eigenvalue displacement solutions were obtained in terms of the Bessel functions. Displacements

obtained were subsequently used to calculate strains and stresses in the system. Elastic strong interface conditions and appropriate homogeneous stress boundary conditions were enforced. Global axial force equilibrium (in an integral sense) was imposed at the loaded end of the fiber. On the outer surface of the matrix, two cases of boundary conditions were studied. A FE study with identical geometry, material properties and loading conditions was carried out to verify the correctness of the analytical solution. Analytically obtained stresses are in reasonably good agreement with FE predictions. It was observed that the shear stress peak occurs close to loaded end of the fiber at IM interface and it reduces to 10% of its peak value within one fiber radius in the axial direction from the loaded end. Shear stress distribution is not influenced by the boundary condition imposed on the outer surface of the matrix. It was also observed that the peak radial stress occur at the loaded end of the fiber at the FI interface. Due to stress singularity at the FI interface at the free surface and mismatch in elastic properties, the FI interface is susceptible to shear induced mixed mode (mode-I and II) failure. The peak stresses are reduced due to the presence of an explicit interlayer. Extension of this idea to additively manufacture composites using multi-material 3D printers at micron- and nano-length scales will enable design of engineered interphases with graded properties. Such an engineered design of interphases with gradation in properties can lead to optimized mechanical performance and structural integrity of composites. Acknowledgements This research was sponsored by the MIT & Masdar Institute Cooperative Program (Project Code: 12MAMA1). Authors would like to acknowledge the financial support from Masdar Institute of Science and Technology for this research work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.mechmat.2015.06.012. References Atkins, A.G., 1975. Intermittent bonding for high toughness/high strength composites. J. Mater. Sci. 10 (5), 819–832. Cech, V., Palesch, E., Lukes, J., 2013. The glass fiber–polymer matrix interface/interphase characterized by nanoscale imaging techniques. Compos. Sci. Technol. 83, 22–26. Christensen, R., Lo, K., 1979. Solutions for effective shear properties in three phase sphere and cylinder models. J. Mech. Phys. Solids 27 (4), 315–330. Cox, H.L., 1952. The elasticity and strength of paper and other fibrous materials. Br. J. Appl. Phys. 3, 72–79. Drzal L. The interphase in epoxy composites. In: Duek K, editor. Epoxy Resins and Composites II. vol. 75 of Advances in Polymer Science. Springer, Berlin Heidelberg; 1986. p. 1–32.. Gao, S.L., Mäder, E., 2002. Characterisation of interphase nanoscale property variations in glass fibre reinforced polypropylene and epoxy resin composites. Compos. Part A Appl. Sci. Manuf. 33 (4), 559–576. Hashin, Z., 1990. Thermoelastic properties of fiber composites with imperfect interface. Mech. Mater. 8 (4), 333–348.

P. Upadhyaya, S. Kumar / Mechanics of Materials 89 (2015) 190–201 Hayes, S., Lane, R., Jones, F., 2001. Fibre/matrix stress transfer through a discrete interphase. Part 1: single-fibre model composites. Compos. Part A Appl. Sci. Manuf. 32 (3), 379–389. Hsueh, C.H., 1990. Interfacial debonding and fiber pull-out stresses of fiber-reinforced composites. Mater. Sci. Eng. A 123 (1), 1–11. Hsueh, C.H., 1991. Interfacial debonding and fiber pull-out stresses of fiber-reinforced composites III: with residual radial and axial stresses. Mater. Sci. Eng. A 145 (2), 135–142. Hsueh, C.H., 1991. Interfacial debonding and fiber pull-out stresses of fiber-reinforced composites IV: sliding due to residual stresses. Mater. Sci. Eng. A 145 (2), 143–150. Hsueh, C.H., 1991. Interfacial debonding and fiber pull-out stresses of fiber-reinforced composites Part V. With a viscous interface. Mater. Sci. Eng. A 149 (1), 1–9. Hsueh, C.H., 1992. Interfacial debonding and fiber pull-out stresses of fiber-reinforced composites VII: improved analyses for bonded interfaces. Mater. Sci. Eng. A 154 (2), 125–132. Huang, Y., Young, R.J., 1996. Interfacial micromechanics in thermoplastic and thermosetting matrix carbon fibre composites. Compos. Part A Appl. Sci. Manuf. 27 (10), 973–980. Hughes, J., 1991. The carbon fibre/epoxy interface–a review. Compos. Sci. Technol. 41 (1), 13–45. Jayaraman, K., Reifsnider, K.L., 1992. Residual stresses in a composite with continuously varying Young’s modulus in the fiber/matrix interphase. J. Compos. Mater. 26 (6), 770–791. Jayaraman, K., Reifsnider, K.L., 1993. The interphase in unidirectional fiber-reinforced epoxies: effect on residual thermal stresses. Compos. Sci. Technol. 47 (2), 119–129. Jiang, Y., Guo, W., Yang, H., 2008. Numerical studies on the effective shear modulus of particle reinforced composites with an inhomogeneous inter-phase. Comput. Mater. Sci. 43 (4), 724–731. Kim, J.K., Mai, Y.W., 1998. Engineered Interfaces in Fiber Reinforced Composites. Elsevier Science. Kiritsi, C., Anifantis, N., 2001. Load carrying characteristics of short fiber composites containing a heterogeneous interphase region. Comput. Mater. Sci. 20 (1), 86–97. Kovalev, S.P., Miranzo, P., Osendi, M.I., 1998. Finite element simulation of thermal residual stresses in joining ceramics with thin metal interlayers. J. Am. Ceram. Soc. 81 (9), 2342–2348. Liu, L., Song, Y., Fu, H., Jiang, Z., Zhang, X., Wu, L., et al., 2008. The effect of interphase modification on carbon fiber/polyarylacetylene resin composites. Appl. Surf. Sci. 254 (17), 5342–5347. Low, B.Y., Gardner, S.D., Pittman Jr., C.U., Hackett, R.M., 1994. A micromechanical characterization of graphite-fiber/epoxy composites containing a heterogeneous interphase region. Compos. Sci. Technol. 52 (4), 589–606. Mogilevskaya, S., Crouch, S., 2001. A Galerkin boundary integral method for multiple circular elastic inclusions. Int. J. Numer. Methods Eng. 52 (10), 1069–1106. Nairn, J.A., 1992. A variational mechanics analysis of the stresses around breaks in embedded fibers. Mech. Mater. 13 (2), 131–154.

201

Nairn, J.A., 1997. On the use of shear-lag methods for analysis of stress transfer in unidirectional composites. Mech. Mater. 26 (2), 63–80. Naslain, R.R., 1998. The design of the fibre-matrix interfacial zone in ceramic matrix composites. Compos. Part A Appl. Sci. Manuf. 29 (9), 1145–1155. Needleman, A., Borders, T., Brinson, L., Flores, V., Schadler, L., 2010. Effect of an interphase region on debonding of a CNT reinforced polymer composite. Compos. Sci. Technol. 70 (15), 2207–2215. Qing, H., 2013. A new theoretical model of the quasistatic single-fiber pullout problem: analysis of stress field. Mech. Mater. 60 (0), 66–79. Qiu, Y., Weng, G., 1991. Elastic moduli of thickly coated particle and fiberreinforced composites. J. Appl. Mech. 58 (2), 388–398. Rjafiallah, S., Guessasma, S., Bizot, H., 2010. Effect of surface etching on interphase and elastic properties of a biocomposite reinforced using glass–silica particles. Compos. Sci. Technol. 70 (8), 1272– 1279. Romanowicz, M., 2010. Progressive failure analysis of unidirectional fiber-reinforced polymers with inhomogeneous interphase and randomly distributed fibers under transverse tensile loading. Compos. Part A Appl. Sci. Manuf. 41 (12), 1829–1838. Shen, L., Li, J., 2003. Effective elastic moduli of composites reinforced by particle or fiber with an inhomogeneous interphase. Int. J. Solids Struct. 40 (6), 1393–1409. Shen, L., Li, J., 2005. Homogenization of a fibre/sphere with an inhomogeneous interphase for the effective elastic moduli of composites. Proc. R. Soc. A Math. Phys. Eng. Sci. 461 (2057), 1475– 1504. Sottos, N.R., McCullough, R.L., Scott, W.R., 1992. The influence of interphase regions on local thermal displacements in composites. Compos. Sci. Technol. 44 (4), 319–332. Tsai, H.C., Arocho, A.M., Gause, L.W., 1990. Prediction of fiber-matrix interphase properties and their influence on interface stress, displacement and fracture toughness of composite material. Mater. Sci. Eng. A 126 (1), 295–304. Wang, Q., Chiang, F.P., 1996. Experimental characterization of interphase mechanical properties of composites. Compos. Part B Eng. 27 (2), 123–128. Wang, J., Crouch, S.L., Mogilevskaya, S.G., 2006. Numerical modeling of the elastic behavior of fiber-reinforced composites with inhomogeneous interphases. Compos. Sci. Technol. 66 (1), 1–18. Wang, X., Zhang, J., Wang, Z., Zhou, S., Sun, X., 2011. Effects of interphase properties in unidirectional fiber reinforced composite materials. Mater. Des. 32 (6), 3486–3492. Wu, Z.J., Ye, J.Q., Cabrera, J.G., 2000. 3D analysis of stress transfer in the micromechanics of fiber reinforced composites by using an eigen-function expansion method. J. Mech. Phys. Solids 48 (5), 1037–1063. Wu, Q., Li, M., Gu, Y., Li, Y., Zhang, Z., 2014. Nano-analysis on the structure and chemical composition of the interphase region in carbon fiber composite. Compos. Part A Appl. Sci. Manuf. 56, 143–149.