Overall properties in fibrous elastic composite with imperfect contact condition

International Journal of Engineering Science 61 (2012) 142–155

Contents lists available at SciVerse ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Overall properties in fibrous elastic composite with imperfect contact condition F.J. Sabina a,⇑, R. Guinovart-Díaz b, R. Rodriguez-Ramos b, J.C. López-Realpozo b, J. Bravo-Castillero b a Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas, Universidad Nacional Autónoma de México, Apartado Postal 20-726, Delegación de Álvaro Obregón, 01000 México D.F., Mexico b Facultad de Matemática y Computación, Universidad de la Habana, San Lázaro y L, Vedado, Habana 4, CP10400, Cuba

a r t i c l e

i n f o

Article history: Received 20 February 2012 Accepted 2 April 2012 Available online 9 August 2012 Keywords: Effective properties Imperfect contact Periodic composites Asymptotic homogenization

a b s t r a c t In this contribution, the complete set of effective elastic moduli are obtained by means of the asymptotic homogenization method (AHM), for two-phase fibrous periodic composites with imperfect contact conditions of linear spring type. This work is an extension of previous reported results, where only perfect contact for elastic composite with square cells were considered. The constituents of the composites exhibit transversely isotropic properties. As validation of the present method, some numerical examples and comparisons with theoretical and experimental results verified that the present model is efficient for the analysis of composites with presence of imperfect interface. The present method can provide benchmark results for other numerical and approximate methods. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Many multi-scale approaches are associated with composites where perfect interface conditions were assumed. Most of the models devoted to the micromechanical analysis of fibre-reinforced composites consider a perfect fibre–matrix bond and suppose that stress and displacement continuity conditions are satisfied along the interface. However, experiments show that local or partial debonding at interfaces is a rule rather than the exception in materials such as reinforced composites Hui-Zu and Tsu-Wei (1995), Rokhlin, Chu, Gosz, and Achenbach (1994). The existence of a stiff region is only an idealization of the complex phenomenon that occurs at the interface where a transition zone (interphase) between the fibre and the matrix is omnipresent. This third phase results from the material fabrication process (chemical treatments of fibre surfaces, resin crystallization), Achenbach and Zhu (1989). In most fibre-reinforced composites, the fibre–matrix adhesion is imperfect; the continuity conditions for stresses and displacements are not satisfied. Thus various approaches have been used, in which the bond between the reinforcements and the matrix is modeled by an interphase with specified thickness and elastic properties different from those of the fibres and the matrix. The research works of Jasiuk et al. (1989), Jasiuk, Chen, and Thorpe (1992) and Hashin (1990) are classified in this group. On the other hand, we will distinguish works supposing the interphase as a boundary layer (no thickness specified) through which it is assumed that the stresses are continuous but not the displacements. This approach has been employed by several authors, like, Hashin (1991a, 1991b) and Hashin (2002).

⇑ Corresponding author. Tel.: +52 55 5622 3544; fax: +52 55 5622 4564 E-mail address: [email protected] (F.J. Sabina). 0020-7225/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2012.06.017

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

143

According to preceding observations the selection of an appropriate model to represent the interphase is important. The micromechanical models are able to answer these requirements. Their purpose is to determine the response of fibre-reinforced composites by the knowledge of the fibre and the matrix behaviors. Some more recently researches are studied the influence of the interphase/interface in the composites. Caporale, Luciano, and Sacco (2006) investigated the behavior of unidirectional fiber reinforced composites with imperfect interfacial bonding by using the finite element method. In this work, an interfacial failure model is implemented by connecting the fibers and the matrix at the finite element nodes by normal and tangential brittle-elastic springs. Furthermore (Artioli, Bisegna, & Maceri, 2010; Bisegna & Caselli, 2008) obtained closed-form expression for the homogenized longitudinal shear moduli of a linear elastic composite material reinforced by long, parallel, circular fibres with a periodic arrangement. An imperfect linear elastic fibre–matrix interface is allowed. The aim of this paper is to point out the effects of the fibre–matrix interface on the mechanical properties of elastic two phase composite materials based on the analytical expressions obtained by micromechanical model of cell assembly, in particular, the two scale asymptotic homogenization method, proposed by Bakhvalov and Panasenko (1989), Pobedrya (1984), Sanchez-Palencia (1980). This work is an extension of previous reported results (see, Guinovart-Díaz, Bravo-Castillero, Rodríguez-Ramos, &Sabina,Sabi Rodriguez-Ramos, Sabina, Guinovart-Díaz, & Bravo-Castillero, 2001 among others) where only perfect contact for elastic composite were considered. In this contribution, the complete set of effective elastic moduli by means of the asymptotic homogenization method (AHM), for two-phase fibrous periodic composites with imperfect contact conditions of linear spring type are obtained. As validation of the present method, the Hills’s universal relations are satisfied and some numerical examples, comparisons with theoretical and experimental results verified that the present model is efficient for the analysis of composites with presence of imperfect interface. 2. Statement of the problem. Local problems based on asymptotic homogenization method The mechanical behavior of imperfect interface is modeled via a layer of mechanical springs of zero thickness. The spring constants are the measures for the magnitude of the associated continuities on C, where Kn, Ks and Kt are spring constant type material parameters which have dimension of stress divided by length. We shall call these constants the interface parameters. It is seen that infinite values of the parameters imply vanishing of displacement jumps and therefore perfect interface conditions. At the other extremity zero values of the parameters implies vanishing of interface tractions and therefore disband. Any finite position values of the interface parameters define an imperfect interface. This may be due to the presence of an interphase but also due to interface bond deterioration caused for reasons such as fatigue damage or environmental and chemical effects. Within this approach the composite is modeled as a two phase material with imperfect interface conditions. The vanishing value of Kn, Ks and Kt corresponds to pure debonding (normal perfect debonding), in-plane pure sliding, and out-of-plane pure sliding, respectively. The status of the mechanical bonding is completely determined by appropriate values of these constants. For large enough values of the constants, the perfect bonding interface is achieved. Using the vector notation and defining the spring stiffness matrix, the mechanic displacement and the traction vectors in the following manner

0

1 un B C u ¼ @ ut A; us

0

1 Tn B C T ¼ @ T t A; Ts

0

en K B K¼@ 0

0 et K

0

0

0

1

C 0 A; es K

ð1Þ

the mechanical imperfect condition (Hashin, 1990; Shodja, Tabatabaei, & Kamali, 2006) in general may be expressed as

Tð1Þ þ Tð2Þ ¼ 0;

TðcÞ ¼ ð1Þcþ1 Kkuk;

on C:

ð2Þ

In these relations kk indicates the jump in the quantity at the common interface between the fiber and the matrix denoted by C; un, ut, us are the tangential and normal components of the mechanic displacement vector; Tn, Tt, Ts are the tangential and normal components of the traction vector T (Ti = rijnj) and n is the outward unit normal on C. The superscripts (c), c = 1,2 denote the matrix and fiber respectively. We can consider this representation to be preferred over the three phase description Guinovart-Díaz, Rodriguez-Ramos, Bravo-Castillero, Sabina, and Maugin (2005), for some reasons, for instance: (a) It is more general in the sense that the three phase description is a special case of (2) since, the interface parameters can be expressed in terms of interphase properties. Furthermore, it is applicable also when an interphase cannot be defined or identified; (b) it characterizes imperfect interface in terms of a small number of parameters which may hopefully be determined experimentally; (c) it is mathematically simpler and more transparent. A two-phase uniaxial reinforced material is considered here in which fibers and matrix have transversely isotropic elastic properties; the axis of transverse symmetry coincides with the fiber direction, which is taken as the Ox3 axis. The fibers cross-section is circular. Moreover, the fibers are periodically distributed without overlapping in directions parallel to the Ox1- and Ox2-axis, w1 = 1/2 and w2 = i/2 are two complex numbers which define the square periodic cell of the two-phase composite. (see, Fig. 1). Therefore the composite X consists of a square array of identical circular cylinders embedded in a homogeneous medium (Fig. 1).

144

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

Fig. 1. The cross section of a quadratic and periodic array of circular fibers.

The overall properties of the above periodic medium are sought using the well-known asymptotic homogenization method. Then, it follows that in terms of the fast variable y, the appropriate periodic unit cell S is taken as a regular square in the y1y2-plane so that S = S1 [ S2 with S = S1 \ S2 = ;, where the domain S2 is occupied by the matrix and its complement S1 (fiber) is considered by a circle of radius R and center at the origin O (Fig. 1). The common interface between the fiber and the matrix is denoted by C. The fiber and matrix associated quantities are also referred below by means of super-indices in brackets (1) and (2), respectively. Using the conventional indicial notation in which repeated subscripts are summed over the range of i, j, k, l = 1, 2, 3, the constitutive equation is

rij ¼ C ijkl ekl ;

ð3Þ

where rij,eij are the stress and strain tensors, respectively. The quantities Cijkl are components of the elastic stiffness tensor. The elastic equilibrium equation as the body forces are absent is,

rij;j ¼ 0; in X;

ð4Þ

where the subscript comma denotes partial differentiation. The gradient equations, which are the strain–displacement equations are

ekl ¼

  1 @uk @ul ; þ 2 @xl @xk

ð5Þ

where ui are the components of the mechanic displacement. Substituting Eqs. (3) and (5) into (4) we obtain a system of partial differential equations with coefficients rapidly oscillating

  C ijkl ðyÞuek;l ðxÞ ¼ 0;

in X:

ð6Þ

;j

Eq. (6) represents a system of equations for finding ui. For a complete solution, it is necessary to assign suitable boundary and interface conditions, for instance

uei ¼ u0i ; u0i ;

reij nj ¼ S0i ; on @ X;

ð7Þ

S0i

where are the prescribed displacement, force on the boundary of the composite. In order to study the imperfect contact conditions (2) we consider the relations between the mechanic displacement and the traction vectors (1) with Cartesian mechanic displacement (ui) and the traction (Ti) vectors respectively by the following expressions

0

un

1

0

cos u

sinu

0

10

u1

1

CB C B C B @ ut A ¼ @ sinu cos u 0 A@ u2 A; us u3 0 0 1

0

Tn

1

0

cos u

sinu

0

10

T1

1

CB C B C B @ T t A ¼ @ sinu cos u 0 A@ T 2 A: Ts T3 0 0 1

ð8Þ

e n ! 1; it is On the interface C, the functions u(x) satisfies a particular case of non-ideal contact conditions (2), taking K e t < 1; 0 < K e s < 1. Thus the equivalent to the continuity of the displacement in the normal direction, kunk = 0 and 0 < K expression (2) on C, can be rewritten in the following indicial form,

Tð1Þ þ Tð2Þ ¼ 0;

ð9Þ

145

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

kun k ¼ 0;

ð10Þ

e t kut k; T tðcÞ ¼ ð1Þcþ1 K

ð11Þ

e s kus k: T sðcÞ ¼ ð1Þcþ1 K

ð12Þ

By means of the AHM it is possible to obtain an asymptotic solution of the boundary-value problem (6)–(12) as e ? 0. The solution is sought in the form of a series in powers of e with coefficients depending on both the variables x and y = x/e treated as independent; they are referred as the slow or macroscopic and fast or microscopic variables, respectively, where e = l/L is a small dimensionless parameter and L is a linear dimension of the body. Here, the solution is explicitly posed as

ue ðxÞ ¼ uð0Þ ðxÞ þ euð1Þ ðx; yÞ þ Oðe2 Þ; (0)

where u

ð13Þ

satisfy the homogenized system of differential equations ð0Þ

C ijlk ul;kj ¼ 0: on X;

ð14Þ (1)

(0)

(1)

and the asterisk denotes the overall elastic properties. The term u represent a correction of the u . The function u is found in combination of the function M(y) (solution of the local problems) and the partial derivatives of the function u(0). Then, the original constitutive relations with rapidly oscillating material coefficients Eq. (6) are transformed into equivalent system Eq. (14) with constant coefficients C⁄ which represent the elastic properties of an equivalent homogeneous medium. They are called effective coefficients of the composite X. The main problem to obtain average coefficients is to find the periodic solutions of six pq L local problems on S in terms of the fast variable y, where p, q = 1, 2, 3. Each local problem uncouples into sets of equations, i.e. plane-strain and antiplanestrain systems. In the following Table 1, the correspondence between the effective properties and the local problems is shown. In the present work, the pq L problem consist to find the displacements pq M(c)(y) in Sc, c = 1,2 (double periodic functions with half periods w1 = 1/2, w2 = i/2) as solution of the following system of partial differential equations pq

ðcÞ rid;d ¼ 0 in Sc ;

ð15Þ

pq

ðcÞ ðcÞ ridðcÞ ¼ C idkk pq M k;k ;

ð16Þ

where

the comma notation denotes a partial derivative relative to the yd component, i.e., U,d  @ U/@yd; the summation convention is also understood for Greek indices, which run from 1 to 2; no summation is carried out over upper case indices, whether Latin on Greek. Thus the expression (9)–(12) on C for the pq L problem can be expressed in the following indicial form, ð1Þ pq T

þ pq Tð2Þ ¼ 0;

ð17Þ

kpq Mn k ¼ 0; ðcÞ pq T t ðcÞ pq T s

ð18Þ

e t kpq Mt k; ¼ ð1Þcþ1 K cþ1 e ¼ ð1Þ K s kpq Ms k;

ð19Þ ð20Þ

where Mt, Ms, Mn and Tt, Ts, Tn have the same meaning as (1) but adequate to pq L problems. To assure the only one solution of the pq L problems, the functions also satisfy that hpq Mki = 0, where the angular brackets define the volume average per unit R length over the unit periodic cell, that is hFi ¼ jSj FðyÞdy. The symmetry between the indices p and q shows right away that at most six problems need to be considered. Once the local problems are solved, the homogenized moduli C ijpq may be determined by using the following formulae:

C ijpq ¼ hC ijpq þ C ijklpq M k;l i:

ð21Þ

The potential method of complex variables z = y1 + iy2, (y1, y2) 2 S and the properties of doubly periodic Weierstrass o P n }ðzÞ ¼ z12 þ 0s;t ðzP1 Þ2  P12 Pst ¼ s þ it; s; t ¼ 0; 1; 2; . . ., and related functions (Z- function f(z) = }0 (z) and Natanzon’s st

st

Table 1 Effective properties related to the local problems. 11L

22L

33L

23L

13L

12L

C 1111 C 2211 C 3311 0 0 0

C 1122 C 2222 C 3322 0 0 0

C 1133 C 2233 C 3333 0 0 0

0 0 0 C 1313 0 0

0 0 0 0 C 2323 0

0 0 0 0 0 C 1212

146

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155 ðaÞ

function Q(z)) are used for the solution of the local problems (15)–(20). Hence, the non-zero solution pq M k in Sa of the problem (15)–(20) must be found among doubly periodic functions of half periods w1 = 1/2, w2 = i/2 (see Fig. 1). Each local problem (15)–(20) uncouples into sets of equations. An in-plane strain system for

ðaÞ pq M k ;

k ¼ 1; 2 and an out-of-plane strain

ðaÞ pq M 3

Laplace’s equation for has to be solved. Then the solution of the in-plane (out-of-plane) strain problems involves the determination of the in-plane (out-of-plane) displacements, strains and stresses over each phase Sa of the composite. ðaÞ

Due to the non-vanishing components of the elastic tensors C ijpq , the only non-homogeneous problems, that have a non-zero solution, correspond to the four in-plane strain problems jj L and 12 L, and the two out-of-plane strain ones 23L and 13L. In that way the solutions of both (in-plane and out-of-plane) local problems (Guinovart-Díaz et al., 2001, Rodriguez-Ramos et al., 2001) (six pq L problems need to be considered due to the symmetry between p and q) lead to obtain the average coefficients of the composite given in Fig. 1. 3. Solution of plane problem

jj L

and the effective coefficients C 1111 ; C 1122 ; C 1133 and C 3333

For the sake of clarity, in this section let U(c) = jj M(c) and the jj presubindices are dropped from all relevant quantities. From the (15)–(20) using the relation Eq. (8) we have Þ ðcÞ rðacd;dÞ ¼ 0 in Sc ; where rðacdÞ ¼ C ðacdbk U b;k ;

ð22Þ

   ðcÞ  U a na ¼ 0 on C;

ð23Þ

   ðcÞ  rad nd  ¼ kC adjj k1 nd

on C;

1

ð24Þ

   i  i h   ðcÞ  ðcÞ ðcÞ ðcÞ ðcÞ ðcÞ  ðcÞ  et  r11  r22 þ C 11jj  C 22jj n1 n2 þ r12 n22  n21 ¼ ð1Þ1þc K on C; U 2 n1  U 1 n2

h

ð25Þ

where ðcÞ ðcÞ ðcÞ r11 ¼ ðkc þ mc ÞU 1;1 þ ðkc  mc ÞU 2;2 ; ðcÞ ðcÞ ðcÞ r22 ¼ ðkc  mc ÞU 1;1 þ ðkc þ mc ÞU 2;2 ; ðcÞ ðcÞ ðcÞ r12 ¼ mc ðU 1;2 þ U 2;1 Þ;

ð26Þ

where k = (C1111 + C1122)/2 is the plane-strain bulk modulus for lateral dilatation without longitudinal extension; m = C1212 = (C1111  C1122)/2 is the rigid modulus for longitudinal uniaxial straining. It can be recognized that the structure of (22) and (26) is one of plane-strain isotropic elasticity except that the usual Lame constants k and l are here identified with kc  mc and mc, respectively. The method of complex variables in terms of two harmonic functions uc(z) and wc(z) and the Kolosov–Muskhelishvili complex potentials are applicable. The potentials are related to the displacement and stress components by means of the classical formulae

  ðcÞ ðcÞ  c ðzÞ;  c ðzÞ  w 2mc U 1 þ iU 2 ¼ jc uc ðzÞ  zu ðcÞ 11

ðcÞ 22

0

ð27Þ

0

r þ r ¼ 2ðuc ðzÞ þ uc ðzÞÞ; ðcÞ ðcÞ ðcÞ r22  r11 þ 2ir12 ¼ 2ðzu00c ðzÞ þ w0c ðzÞÞ;

ð28Þ ð29Þ

where the prime denotes a derivative with respect to z, the over bar a complex conjugate and jc ¼ 3  4mTc , here mTc is the transverse Poisson’s ratio, The non-ideal contact conditions (23)–(25) are transformed according to the Kolosov–Muskhelishvili complex potentials by the formulae (28) and (29) in the following form

h

i

h

i

vm ½j1 /1 ðzÞ  z/=1 ðzÞ  w1 ðzÞ  j2 /2 ðzÞ  z/=2 ðzÞ  w2 ðzÞ þ vm j1 /1 ðzÞ  z/=1 ðzÞ  w1 ðzÞ e2ih  =

h

i

j2 /2 ðzÞ  z/=2 ðzÞ  w2 ðzÞ e2ih ¼ 0;

ð30Þ =

z/1 ðzÞ þ w1 ðzÞ þ /1 ðzÞ þ zc1j þ zc2j ¼ z/2 ðzÞ þ w2 ðzÞ þ /2 ðzÞ; h i h i = = ih 3ih  2Rvm z/==  2c3j Rvm ðe3ih  eih Þ 2Rvm z/== 2 ðzÞ þ w2 ðzÞ e 2 ðzÞ þ w2 ðzÞ e 8 h 9 i h i > > < vm j1 /1 ðzÞ  z/=1 ðzÞ  w1 ðzÞ  j2 /2 ðzÞ  z/=2 ðzÞ  w2 ðzÞ  = h i h i ¼ K t ; > : vm j1 /1 ðzÞ  z/=1 ðzÞ  w1 ðzÞ e2ih  j2 /2 ðzÞ  z/=2 ðzÞ  w2 ðzÞ e2ih > ;

ð31Þ

ð32Þ

147

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

where ð1Þ ð2Þ ð2Þ ð1Þ ð1Þ ð2Þ ð1Þ ð2Þ ð2Þ ð2 Þ vm ¼ m2 =m1 ; 2c1j ¼ C 22jj  C 22jj þ C 11jj  C 11jj ; 2c2j ¼ C 11jj  C 11jj þ C 22jj  C 22jj ; 2c3j ¼ C 22jj  C 11jj

ð33Þ

e t R=m1 is a dimensionless parameter. The complex potentials uc(z) and wc (z) are looked for the periodic cell that and K t ¼ K contains the origin of coordinates in the following form

u1 ðzÞ ¼ a0 z þ w1 ðzÞ ¼ b0 z þ

1 ðk1Þ X o ak f ðzÞ

ðk  1Þ!

k¼1 1 X o

; ð34Þ

o

1 bk fðk1Þ ðzÞ X ak Q ðk1Þ ðzÞ þ ; ðk  1Þ! ðk  1Þ! k¼1

k¼1

for the matrix phase and

u2 ðzÞ ¼

1 X o

c k zk ;

w2 ðzÞ ¼

k¼1

1 X o

dk zk ;

ð35Þ

k¼1

for the fiber phase. In the foregoing definitions f(k1)(z) denotes the doubly–periodic (k  1)th derivative of the Weierstrass Zeta quasi–periodic function f(z), Q(k1)(z) is the doubly–periodic (k  1)th derivative of Natanzon’s quasi–periodic function, and the symbol Po indicates a sum over the indices k = 1, 3, 5, . . ., meanwhile coefficients a0,b0,ak, bk, ck and dk (k = 1, 3, 5, . . .) are all real and undefined. Next, by taking into consideration the double periodicity of matrix displacements U(1), by means of Legendre’s relation and the properties of the doubly periodic functions, it is found that coefficients a0 and b0 are

a0 ¼

V2 b ; j1  1 1

b0 ¼ P

where V2 = pR2 and S4 ¼



j1 þ

5S4

p

2

 V 2 a1 ;

ð36Þ

1 p;q ðpþiqÞ4 .

Replacing (34) and (35) into (30)–(32) the following relations between the unknown constants of the above expansions are obtained

bpþ2 ¼ pap 

1 0 X

gkðpþ2Þ ak þ cpþ2 ;

ð37Þ

k¼1

ðp þ 2Þcpþ2 þ dp ¼ ap þ b0 d1p þ Rc1j d1p þ g1p b1 þ

1 o X

pþk

ðp þ 2Þgkðpþ2Þ ak þ gkþ2p bkþ2 þ kR

C ppþk T pþk ak ;

ð38Þ

k¼1

cpþ2 ¼  b1 ¼ F

1 0 X Dp Bp K t vm ðj1 þ 1Þap  K t vm ðj1 þ 1Þ gkðpþ2Þ ak ; Ep Ep k¼1

1 o X

ð39Þ

gk1 ak  PRc2j ;

ð40Þ

k¼1

"

#

1 o X vm ðj1 þ 1Þ gk1 ak þ ð1  V 2 ÞRc2j ; 2a0 k¼1   1 þ vm j1  vm  j2 j2  1 V2 1 ; F¼ ; P¼ ; a0 ¼ vm V 1 þ ðj2  1Þ þ 2a0 a0 j1  1 2

c1 ¼

ð41Þ

where dik is the Kronecker’s delta function, C nk is the binomial coefficient, V1,V2 are the area of the matrix and fiber respectively, V1 + V2 = 1. The parameters involved in this expression are listed as follows

gkp ¼ k

ðk þ p  1Þ! Skþp Rkþp ; k!p!

Snþk ¼

X

nþk

p;q

p2 þ q2 – 0; p; q  integer numbers;

Ap ¼ vm ðp þ 1Þ  K t ðvm  1Þ ðp þ 2Þ; C p ¼ K t ðvm þ j2 Þ þ v ðp þ 2Þðp þ 1Þ;

1 ðp þ iqÞ

;

T nþk ¼

X

p  iq

p;q

ðp þ iqÞnþkþ1

;

Bp ¼ K t ð1  vm Þ þ vm p; Dp ¼ vm p;

Ep ¼ Ap Dp  C p Bp :

The above expression Eqs. (37)–(41) depends on the unknown parameter ak which can be calculated from the following system of algebraic equations

ap þ B11 V 2



j1 þ

1 1p Bp PR 2j

¼g

5S4

p

2



a1 d1p þ

c þ Ec1j Rd1p ;

1 o X

 1 B1p Gkp  B1p C 0k r kp þ B1p F g1p gk1 þ B0p A0p gkpþ2  B1p D0p gkþ2p ak

k¼1

ð42Þ

148

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

where p = 1, 3, 5, . . .,

 Bp B0p ¼ vm ð1 þ j1 Þ 1 þ ðvm þ j2 Þ K t ; Ep  1  1 D p D0p ¼ K t vm ð1 þ j1 Þ; B1p ¼ ð1  vm Þ A0p ; E ¼  A01 ; Ep

A0p ¼ 1 þ vm j1 þ K t vm ð1 þ j1 Þðvm þ j2 Þ C 0p ¼ 1 þ K t vm ð1 þ j1 Þ r kp ¼

1 o X

Bp ; Ep

Dp ; Ep

gki gip ; Gkp ¼ ðp þ 2Þgkðpþ2Þ þ kgkþ2p þ kRpþk C ppþk T pþk :

i¼3

It should be noticed that Eq. (42) represents an infinite linear system from which we can calculate the sought coefficients ap. The unknown coefficients ap can be obtained as in Rodriguez-Ramos et al. (2001). We present here only the most important results in order to obtain the effective properties. The residue a1 of the function u1(z) is

Ec1j R

a1 ¼

1 þ Hþ  U1 U1 U2

where Hþ ¼ B11 V 2



ð43Þ

;



j1 þ 5S4 =p2 þ B11 G11  B11 C 01 r11 þ B01 g13 =A01  B11 D01 g31 and U; U1 and U2 denotes a matrix and two vectors

of infinite order respectively and U1 is the inverse of the matrix U, given by the following expressions

h U½uts  ¼ dts þ B14sþ1 G4tþ14sþ1  B14sþ1 C 04tþ1 r4tþ14sþ1 þ B04sþ1 g4tþ1 4sþ3 =A04sþ1  B14sþ1 D04sþ1 g4tþ3 h i U1 ½us  ¼ B14sþ1 G14sþ1  B14sþ1 C 01 r 14sþ1 þ B04sþ1 g1 4sþ3 =A04sþ1  B14sþ1 D04sþ1 g3 4sþ1 ; h i U2 ½ut  ¼ B11 G4tþ11  B11 C 04tþ1 r 4tþ11 þ B01 g4tþ1 3 =A04sþ1  B11 D01 g4tþ31 :

i 4sþ1

;

Once the solution of the previous system (42) is obtained, the coefficient c1 in (41) can be written in the following form

vm c2j R Pðj1 þ 1ÞN1 Z1 N2 þ V 1 ; 2a0

c1 ¼

ð44Þ

where again the capital letters Z; N1 and N2 denotes a matrix and two vectors of infinite order respectively, with the following components

  ZðZts Þ ¼ dts þ B14s1 G4t14s1  B14s1 C 04t1 r4t14s1 þ B14s1 F g4t11 g14s1 þ B04s1 g4t14sþ1 A04s1  B14s1 D04s1 g4tþ14s1 ;   N1 ðNs Þ ¼ B14s1 g14s1 ;

N2 ðNt Þ ¼ ðg4t11 Þ;

t; s ¼ 1; 2; 3; . . .

In order to obtain analytic expressions of the effective properties C 1111 ; C 1122 ; C 1133 and C 3333 , is necessary to apply the Green’s theorem to (21) on the interface C and using the condition (27) and the orthogonality of the trigonometric functions we get

C 11jj þ C 22jj ¼ hC 11jj þ C 22jj i  2pR C 22jj C 33jj

kkkðj2  1Þc1 ; vm m1

ð46Þ

ð47Þ

¼ hC 22jj  C 11jj i þ 2pR½ðj1 þ 1Þa1 þ c1j ;  kC 1133 kðj2  1Þc1 ¼ C 33jj  pR : vm m1



ð45Þ

C 11jj

The following analytical expressions of the effective properties are obtained replacing the previous equalities Eqs. (33), (43) and (44) into Eqs. (45)–(47)

C 1111 ¼ hC 1111 i þ V 2 ðkkk2 K1 þ kmkK2 Þ;

ð48Þ

C 1122 ¼ hC 1122 i þ V 2 ðkkk2 K1  kmkK2 Þ;

ð49Þ

C 1133 ¼ C 2233 ¼ hC 1133 i þ V 2 kC 1133 kkkkK1 ;

ð50Þ

2

C 3333

¼ hC 3333 i þ V 2 kC 1133 k K1 ;

K1 ¼

  ð1  j2 Þ j2  1 Eðj1 þ 1Þ V1 þ ðj1 þ 1ÞN1 Z 1 N2 ; K2 ¼ þ1 : þ 1 2m1 a0 2a0 1 þ H  U1 U U2

where

4. The plane problem

12L

Let U(c) = 12M(c), from the (15)–(20) and using the relation (8), the local problem is stated as follows

ð51Þ

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155 Þ ðcÞ rðacd;dÞ ¼ 0 in Sc ; where rðacdÞ ¼ C ðacdbk U b;k ;

   ðcÞ  U a na ¼ 0 on C;    ðcÞ  rad nd  ¼ kC ad12 knd on C;    i h i   h   ðcÞ  ðcÞ ðcÞ ðcÞ ðcÞ ðcÞ  et  r11  r22 n1 n2 þ r12 þ C 1212 n22  n21 ¼ ð1Þ1þc K on C: U 2 n1  U 1 n2

149

ð52Þ ð53Þ ð54Þ ð55Þ

Once the plane problems jj L have been solved, it is easier to solve the plane problem 12L. This can be realized by means of the next relationship between the Kolosov–Muskhelishvili potentials for plane local problems, 12

uc ðzÞ ¼ iðjj uc ðzÞÞ;

12 wc ðzÞ

¼ iðjj wc ðzÞÞ;

ð56Þ

pffiffiffiffiffiffiffi where i ¼ 1 is the imaginary unit and jjuc(z) are the potential functions given in (34) and (35) (Pobedrya, 1984). The unknown coefficients ap that appear now in the relation (56) can be obtained as solution of following system of equations

ap þ B11 V 2



j1 

5S4



p2

a1 d1p 

1 o X k¼1

 1 1 B1p Gkp þ B1p C 0k r kp þ B0p A0p gkpþ2  B1p D0p gkþ2p ak ¼ 0 kmkRd1p : Ap

ð57Þ

Solving the above system (Rodriguez-Ramos 2001), we have

a1 ¼

EkmkR ; 1 þ H  V1 V1 V2

ð58Þ

where

  1 H ¼ B11 V 2 ðj1  5S4 =p2 Þ  B11 G11 þ B11 C 01 r 11 þ B01 A01 g13  B11 D01 g31 ;   1 g14sþ3  B14sþ1 D04sþ1 g34sþ1 ; V1 ½v1s  ¼ B14sþ1 G14sþ1 þ B14sþ1 C 01 r 14sþ1 þ B04sþ1 A04sþ1   1 g4tþ13  B11 D01 g4tþ31 ; V2 ðv12 Þ ¼ B11 G4tþ11 þ B11 C 04tþ1 r 4tþ11 þ B01 A01   1 g4tþ14sþ3 þ B14sþ1 D04sþ1 g4tþ34sþ1 : V½vts  ¼ dts  B14sþ1 G4tþ14sþ1  B14sþ1 C 04tþ1 r4tþ14sþ1  B04sþ1 A04sþ1 The effective coefficient C 1212 is obtained from (21). Again it is necessary to apply the Green’s theorem on the interface C and taking into consideration the expression (58), we can write the following expression for the shear module, ð1Þ

C 1212 ¼ C 1212 þ kC 1212 kV 2 M;

ð59Þ

j1 þ1ÞE where M ¼ 1þHðþV . V1 V 1

2

5. The antiplane problem

13L

and the effective coefficient C1313

We have previously reported in Lopez-Realpozo, Rodriguez-Ramos, Guinovart-Díaz, Bravo-Castillero, and Sabina (2011) the axial effective properties of composites with rhombic periodic cell. From formula (36) and (39) of this abovementioned work, the expression of C 44 is computed as a particular case, i.e. when the basic cell is square. The final expression in this case is given by ð1Þ

C 1313 ¼ C 2323 ¼ C 1313 ð1  V 2 K3 Þ;

ð60Þ

where

K3 ¼

2b1 ; 1 þ b1 V 2  b1 N1 Y1 N2

N1 ðns Þ ¼ ðb4s1 g14s1 Þ;

bp ¼

ð1  vÞK s þ vp ; ð1 þ vÞK s þ vp

ð1Þ v ¼ C ð2Þ 1313 =C 1313 ;

N2 ðnt Þ ¼ ðg4t11 Þ and Yðnts Þ ¼ dts  b24s1 r 4t14s1 :

Finally, as a summary, we have calculated the exact analytic formulae for the effective coefficients C 1111 ; C 1122 ; C 1133 ; C 3333 ; C 1313 and C 1212 given by the formulae (48)–(51), **Eqs. (58) and (59) for fibrous composites with square unit cell under spring imperfect contact at the interface. These closed–form expressions obtained show an explicit dependence on: (a) The mechanical properties of the constituents, (b) the volumetric fraction of each material phase and, (c) the dimensions of the representative volume of analysis and (d) the tangential imperfect interface parameters Kt and Ks. Besides, the following connections between the effective properties for two phase elastic fibrous composites with transversely isotropic constituents and parallelogram periodic cell are obtained in a similar way as Guinovart-Díaz et al. (2001) and Rodriguez-Ramos et al. (2001),

150

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

   C 1111 þ C 1122 =2  hðC 1111 þ C 1122 Þ=2i C 1133  hC 1133 i kðC 1111 þ C 1122 Þ=2k ¼  ¼ : C 1133  hC 1133 i C 3333  hC 3333 i kC 1133 k

ð61Þ

These universal relations are the Hill’s universal relations (see Hill 1964). Furthermore, (61) are valid for any shape of the interface C, they are independent of the fiber volume fraction and the imperfect parameters. 6. Analysis of numerical results In this section, some numerical results will be shown as a consequence of the obtained analytic expressions (48)–(51), (59) and (60). Moreover, comparisons with other theoretical results will be given in order to validate the aformentioned expressions. For simplicity, the effective properties are written in the short index notation. The computations by AHM were made for N0 = 10, where N0 denotes the number of equations considered in the solution of the infinite algebraic system of Eqs. (42) and (57). In general, for low volume fraction of fiber (V2 < 0.4) the accuracy and convergence of the results are good for small value of N0 (N0 6 2). Only for a high volume fraction the value of N0 requires to be changed for getting better accuracy. The effect of interface in the overall properties of composites is studied in Figs. 2–4. The considered material constituents ð2Þ ð1Þ in these figures are isotropic with C 66 =C 66 ¼ 10 and Poisson’s ratio taken as m1 = m2 = 0.3. Fig. 2 shows the ratio of the effecð1Þ  tive transverse shear moduli C 66 =C 66 versus fiber volume fraction V2 for different values of imperfect parameter Kt. Note that the imperfect parámeter Kt has a significant effect on C 66 . For instance, Kt approaches to infinity corresponds to perfect case bonding. Moreover, as Kt tends to zero the coefficient C 66 asymptotically approaches to the solutions for pure sliding case. The effective shear modulus increases as the fiber volume fraction increases whereas this effective modulus decreases as Kt decreases as well. Fig. 3 displays the ratio of the bulk effective modulus k⁄/k1 versus fiber volume fraction V2 for differents values of imperfect parameter Kt. This figure illustrates that k⁄ is independent of the imperfect parameter Kt, due to the fact of continuity of the displacement in the normal direction are considered, in contrast with the previous figure. Small differences are observed only for high-volume fraction. Similar behavior was reported by Jasiuk and Tong (1989). ð1Þ In Fig. 4 the change of the effective shear modulus C 66 =C 66 with respect to the reciprocal imperfect parameter Kt for a fiber volume fraction V2 = 0.5 is sketched. The present spring model is compared with two approaches, finite element method (FEM) reported in Otero et al. (2012) and three phase model published by Guinovart-Díaz et al. (2005). The local in-plane and out-plane problems in FEM are solved using the semi-analytical method combining asymptotic homogenization and finite element methods. Three phase model find analytically the solutions of the local problems, considering the existence of

ð1Þ

Fig. 2. Effect of imperfect parameter Kt on effective transverse shear modulus C 66 =C 66 versus fiber volume fraction V2.

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

151

   Fig. 3. The plane-strain bulk modulus for lateral dilatation without longitudinal extension k ¼ C 11 þ C 12 =2 versus fiber volume fraction V2 for different values of imperfection parameter.

an isotropic thin layer interphase between the matrix and fibers. In this model, the properties of this intermediate layer are considered as a function of the imperfect parameters Kn, Kt and Ks, and the relation (19) reported by Hashin (2002) is used. ð1Þ In the present case K t ¼ K s ¼ C I66 =V I C 66 and Kn = 1012, where the capital superscript I denotes the intherphase property and the interphase volume fraction is VI = 2.51  104. Fig. 4. exhibits good match between the spring, FEM and three phase models in all range of variation of 1/Kt. Notice that C 66 corresponds to the perfect case bonding for Kt ? 1 whereas C 66 approaches asymptotically to the solutions of fiber porous case as Kt ? 0. All effective properties of fibrous composites with empty fibers are presented in Table 2. Young’s modulus E1 = 70 GPa and Poisson’s ratio m1 = 0.3 are the material parameters used in the computation. The set of all effective elastic coefficients are given for different values of porous volume fraction V2. A comparison between the spring model and finite element method (FEM) Otero et al. (2012) is listed. The parameters Kn and Kt are zero for FEM and for the spring models is considered Kt = 0.01 and the fibers properties are zero. The effective properties decrease as the porosity increase. It can be noticed a very good coincidence between these two models. The following numerical calculations are conducted for three composites: glass/epoxy, carbon/epoxy and FP (Al2O3,)/Al reported in Hui-Zu and Tsu-Wei (1995). The carbon fiber is transversely isotropic, and all of the other constituents of the composites mentioned above are isotropic. The elastic properties of the fiber and matrix materials are listed in Table 3. To investigate the consistency of the calculation with experiments results, Figs. 5, 6 show a comparison between AHM spring model for differents values of imperfect parameter Kt and experimental data reported by Tsai (1980) and Zhu, Qi, Wu, Leung, and Choy (1986) which can be indirectly taken from Figs. 9 and 10 by Hui-Zu and Tsu-Wei (1995). Fig. 5 illustrates the analytical transverse Young’s modulus Et of the glass/epoxy composite. Fig. 6 displays the analytical axial shear modulus Ga of the carbon/epoxy composite. In both figures, the corresponding experimental results are also plotted, showing good agreement between analytical and experimental results. Fig. 5 shows that experimental results constitute a dispersed cloud, which might be motivated by the effect of the contact between the phases. Observe that for different values of the imperfect parameter Kt the experiments are reached to the theoretical model. In Fig. 6 the experimental results are closed to the theoretical prediction for Kt = 1010 which corresponds with perfect bonding. In Figs. 5 and 6, lines with symbols ‘‘+’’ represent the analytical results assuming a composite with empty fiber reported by Sabina, Bravo-Castillero, Guinovart-Díaz, Rodriguez-Ramos, and Valdiviezo-Mijangos, (2002). The curves corresponding to the perfect bonding (Kt = 1 1010) and empty fiber are upper and lower bounds respectively for all the curves obtained as the imperfect parameterKt belongs to the interval (0,1). Fig. 7 gives a comparison between the spring and three phase models for the six analytical engineering constants (axial and transverse Young, Poisson’s ratio and shear moduli) for FP (Al2O3,)/Al transversely isotropic composite. In order to study

152

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

ð1Þ

Fig. 4. The effective transverse shear moduli C 66 =C 66 versus reciprocal imperfect parameter Kt for a fiber volume fraction V2 = 0.5. The spring model is compared with finite element method (FEM) and three phase model.

Table 2 Comparison between the analytical results obtained by spring model and the finite element method (FEM) for porous fiber composites.

C 11 C 12 C 13 C 33 C 44 C 66

V2

0.05

0.2

0.35

0.55

0.75

FEM Spring FEM Spring FEM Spring FEM Spring FEM Spring FEM Spring

80.525367 80.525340 33.149807 33.149803 34.102552 34.102543 86.961532 86.961526 24.358970 24.358970 23.209442 23.209425

53.390014 53.390011 18.365994 18.365993 21.526802 21.526801 68.916082 68.916081 17.945057 17.945057 13.422078 13.422076

36.536036 36.536035 9.781059 9.781059 13.895129 13.895128 53.837077 53.837077 12.915297 12.915297 6.616327 6.616326

20.498584 20.498584 3.378678 3.378679 7.163179 7.163179 35.797907 35.797907 7.459089 7.459089 1.808771 1.808767

5.849942 5.849969 0.289122 0.289123 1.841719 1.841728 18.60504 18.60504 2.111431 2.111428 0.078901 0.078818

Table 3 Elastic constants of constituents of composite glass/epoxy, carbon/epoxy and FP/Al.

Glass fiber Epoxy resin FP alumina fiber Aluminium Epoxy resin Carbon Axial Fiber Transverse

Young’s modulus E (GPa)

Poisson’s ratio m

Shear modulus G(GPa)

73.1 3.45 379 68.9 5.35 232 15

0.22 0.35 0.2 0.345 0.354 0.279 0.490

30 1.28 157.9 25.6 1.98 5.03 –

the influence of the imperfect parameter, in the spring model, different values for the constant Kt are considered. The properties of the interphase are estimated assuming it as isotropic material in the three phase model. The symbols without color denote the spring model and with red color the three phase model. The coefficients are showing good agreement between

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

153

Fig. 5. The transverse Young’s modulus Et of glass/epoxy composite for various imperfect parameters are compared with experimental results. Perfect bonding (⁄) and empty fiber property (+) are upper and lower bound for different imperfect cases 0 < Kt < 1.

Fig. 6. The transverse shear axial modulus Ga of carbon/epoxy composite for various imperfect parameters are compared with experimental results.

154

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

Fig. 7. Comparison between the spring and three phase models for FP (Al2O3,)/Al transversely isotropic composite. The axial and transverse Young, Poisson’s ratio and shear moduli are calculate for different imperfect parameters. Empty fiber is denoted by the symbol ‘‘+’’.

the two analytical results. The lines with symbols ‘‘+’’ represent the analytical results assuming a composite with empty fiber. 7. Conclusions In the paper analytical expressions are obtained using asymptotic homogenization method for all elastic coefficients considering imperfect contact conditions at the interface between the matrix and fiber. The formulae are functions of imperfect parameters. The Hill universal relations are satisfied for these formulae. The following items are obtained from the previous numerical analysis where the interface effects are modeled as distributed mechanical springs. (a) The mechanical behaviors of fibrous composite are influenced by the stiffness of the spring parameter. The present numerical results of the effective shear and Young modulus versus the fiber volume ratio are correlated very well with other results for perfect and imperfect bonding case, and by the three phase model. (b) At the fixed fiber volume ratio, the effective shear modulus of composite increases as the stiffness of the interphase increases. (c) In the asymptotic limit, we can simulate different degrees of the interface’s response, the case of the perfect bonding as Kt, Ks ? 1 and complete separation of the matrix and inclusions as Kt, Ks ? 0. (d) The analytical expressions are validated by experimental results.

Acknowledgments The funding of the CoNaCyT project 129658 is acknowledged. Thanks are due to Ana Pérez Arteaga and Rámiro Chávez Tovar for computational help. References Achenbach, J. D., & Zhu, H. (1989). Effect of interfacial zone on mechanical behavior and failure of fiber-reinforced composites. Journal of the Mechanics and Physics of Solids, 37(3), 381–393.

F.J. Sabina et al. / International Journal of Engineering Science 61 (2012) 142–155

155

Artioli, E., Bisegna, P., & Maceri, F. (2010). Effective longitudinal shear moduli of periodic fibre reinforced composites with radially-graded fibres. International Journal of Solids and Structures, 383–397. Bakhvalov, N. S., & Panasenko, G. P. (1989). homogenization averaging processes in periodic media. Dordrecht: Kluwer. Bisegna, P., & Caselli, F. (2008). A simple formula for the effective complex conductivity of periodic fibrous composites with interfacial impedance and applications to biological tissues. Journal of Physics D: Applied Physics, 41, 115506. Caporale, A., Luciano, R., & Sacco, E. (2006). Micromechanical analysis of interfacial debonding in unidirectional fiber-reinforced composites. Computers & Structures, 84, 2200–2211. Guinovart-Díaz, R., Bravo-Castillero, J., Rodríguez-Ramos, R., & Sabina, F. J. (2001). Closed-form expressions for the effective coefficients of fibre-reinforced composite with transversely isotropic constituents – I. Elastic and square symmetry. Journal of the Mechanics and Physics of Solids, 49 1445–1462. Guinovart-Díaz, R., Rodriguez-Ramos, R., Bravo-Castillero, J., Sabina, F. J., & Maugin, G. A. (2005). Closed-form thermoelastic moduli of a periodic three-phase fiber reinforced composite. Journal of Thermal Stresses, 28, 1067–1093. Hashin, Z. (1990). Thermoelastic properties of fiber composites with imperfect interface. Mechanics of Materials, 8, 333–338. Hashin, Z. (1991a). The spherical inclusion with imperfect interface. Journal of Applied Mechanics, 58, 444–449. Hashin, Z. (1991b). Thermoelastic properties of particulate composites with imperfect interface. Journal of the Mechanics and Physics of Solids, 39, 745–762. Hashin, Z. (2002). Thin interphase/imperfect interface in elasticity with application to coated fiber composites. Journal of the Mechanics and Physics of Solids, 50, 2509–2537. Hill, R. (1964). Theory of mechanical properties of fibre-strengthened materials: I. Elastic behavior. Journal of the Mechanics and Physics of Solids, 12, 199–212. Hui-Zu, S., & Tsu-Wei, C. (1995). Transverse elastic moduli of unidirectional fiber composite with fiber/matrix interfacial debonding. Composites Science and Technology, 5, 383–391. Jasiuk, I., & Tong, Y. (1989). The effect of interface on the elastic stiffness of composites. Mechanics of Composite Materials, 100, 49–54. Jasiuk, I., Chen, J., & Thorpe, M. F. (1992). Elastic moduli of composites with rigid sliding inclusions. Journal of the Mechanics and Physics of Solids, 40, 373–391. Lopez-Realpozo, J. C., Rodriguez-Ramos, R., Guinovart-Díaz, R., Bravo-Castillero, J., & Sabina, F. J. (2011). Transport properties in fibrous elastic rhombic composite with imperfect contact condition. International Journal of Mechanical Sciences, 53, 98–107. Pobedrya, B. E. (1984). Mechanics of composite materials. Moscow: Moscow State University Press. in Russian. Rodriguez-Ramos, R., Sabina, F. J., Guinovart-Díaz, R., & Bravo-Castillero, J. (2001). Closed-form expressions for the effective coefficients of fibre-reinforced composite with transversely isotropic constituents-I. Elastic and square symmetry. Mechanics of Materials, 33, 223–235. Rokhlin, S. I., Chu, Y. C., Gosz, M., & Achenbach, D. J. (1994). Determination of interphasial elastic properties in fiber composites from ultrasonic bulk wave velocities. In D. O. Thompson & D. E. Chimenti (Eds.). Review of progress in QNDE (Vol. 13, pp. 1461–1468). New York: Plenum Press. Sabina, F. J., Bravo-Castillero, J., Guinovart-Díaz, R., Rodriguez-Ramos, R., & Valdiviezo-Mijangos, O. C. (2002). Overall behavior of two-dimensional periodic composites. Int. J. Solids Struct, 39(2), 483–497. Sanchez-Palencia, E. (1980). Non homogeneous media and vibration theory. Lecture notes in physics (Vol. 127). Berlin: Springer. Shodja, H. M., Tabatabaei, S. M., & Kamali, M. T. (2006). A piezoelectric inhomogeneity system with imperfect interface. International Journal of Engineering Science, 44, 291–311. Otero, J. A., Rodriguez-Ramos, R., Bravo-Castillero, J., Guinovart-Díaz, R., Sabina, F.J., Monsivais, G. (2012). Semi-analytical method for the computing effective properties in elastic composite under imperfect contact. International Journal of Solids and Structures In press, http://dx.doi.org/10.1016/ j.ijengsci.2012.06.017. Tsai, S. W. (1980). Introduction to composite materials. Westport, CT: Technomic Publishing Co., Inc.. Zhu, X., Qi, Z., Wu, R., Leung, W. P., & Choy, C. L. (1986). Elastic moduli of unidirectional carbon fiber composite. In T. T. Loo & C. T. Sun (Eds.), Proceedings of the international symposium on composite materials and structures (pp. 5–40). Westport, CT: Technomic Publishing Co., Inc..