FINITE ELEMENTS IN ANALYSIS A N D DESIGN
ELSEVIER
Finite Elements in Analysis and Design 23 (1996) 333347
The influence of imperfections on the creep behavior of woven polymer composites at elevated temperatures S. Govindarajan*, N.A. Langrana, G.J. Weng Department of Mechanical and Aerospace Engineering, Rutgers University, New Brunswick, NJ 08855, USA
Abstract
The dynamic behavior of compression molded polymer/woven graphite fiber composites at elevated temperatures is investigated analytically. This is performed with the objective of predicting the initiation of catastrophic failure that may occur after prolonged usage of the material at these temperatures. Special attention is paid to the behavior of the voids present in them where the failure may occur. The polymer matrix is modeled as a 4parameter model (MaxwellVoigt combination) (Govindarajan et al., in: Advances in ComputerAided Engineering, ASME, 1994) while the composite structure is modeled using the fiber undulation model (Ishikawa and Chou, J. Mater. Sci. 17, pp. 32113220, 1982). The relation between the polymer properties and the ambient temperature is modeled after Arhenius' relation (Govindarajan et al., 1994; Ferry, Viscoelastic Properties of Polymers, Wiley, New York, 1961). The multiple phases in the matrix are taken into account through Eshelby's theory (Proc. Royal Soc. London A 241, pp. 376396, 1957) and its extension for multiple occurrences of the same phase (Tanden and Weng, Polymer Composites 5, pp. 327333, 1984; Weng, Internat. J. Eng. Sci. 22 (7), pp. 845856, 1984) which assumes an ellipsoidal shape for inclusions. The resulting elastic equations are transformed into the time domain using Laplace transformation and the correspondence principle (Govindarajan et al., 1994; Wang and Weng, ASME J. Appl. Mech., 1992). All the voids are considered to be prolate ellipsoids with the 1axis being the axis of symmetry. The distribution of voids is assumed to be of a Gaussian form with respect to the aspect ratio. The response of the composite under creep condition (constant load) has been simulated. Relations between the applied stress and the stresses in the matrix/void phase are also supplied, so that the influence of the voids may be characterized. The model is then applied to simulate the behavior of an epoxy/woven graphite composite to obtain the numerical results.
Keywords." Woven polymer composites; Viscoelasticity; Elevated temperature; Creep; Ellipsoidal voids; Gaussian distribution; Laplace transform
1. Introduction Polymer composites have assumed importance in the current engineering industry owing to their ease of processing and a higher strengthtoweight ratio than most metals and alloys. Current research in polymer science has resulted in the development of new heatresistant polymers that also * Corresponding author. 0168874X/96/$15.00 (~) 1996 Elsevier Science B.V. All rights reserved
PII S 0 1 6 8  8 7 4 X ( 9 6 ) 0 0 0 3 52
334
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
possess better longterm properties. These developments and the inherent advantages of the composite materials themselves have increased the potential for usage of these materials in the aircraft and engineering industries. However, the variation of their material properties with temperature and the presence of multiple phases (inclusions, resin and voids) make analytical modeling of their behavior complicated. As a result, most of the current research undertakes to explore their fundamental behavior traits through simple models. It is now understood that the constituent ratio in a polymer composite can be controlled by varying the process parameters. In a compression molded composite, this translates to varying the temperature and applied pressure cycles [1]. Thus, different fiber, resin and void volume fractions can be obtained through different process cycles. While the volume fraction of the voids may be controlled to a certain extent during manufacturing, it is impossible to restrict their aspect ratio (and the shape) within a certain value or range. Thus, in our modeling, the void volume density is formulated through a Gaussian distribution of the aspect ratio. The upper and lower bounds of the aspect ratio are obtained through observation of sections of the specimen [2]. The presence of voids affects the dynamic behavior of the composite significantly, particularly at elevated temperatures. It has been observed that the void fraction in the polymer matrix increases under continuous loading which results in increased stresses in other phases. Thus, in addition to weakening the material, the voids may serve as starting points for crack initiation and propagation. This is especially true in high temperature creep and fatigue, where the strength of the material is considerably lowered by prolonged usage. In this study, we aim to investigate these phenomenon and use it as a means to predict the onset of catastrophic failure under adverse conditions. As far as we know, this approach of considering the different phases and their volume distribution and geometry has not been used before to simulate the composite's behavior.
2. Analytical modeling 2.1. M a t r i x modeling
The matrix is modeled as an isotropic viscoelastic material in which tiny pockets of volatiles are trapped. These volatiles are produced during the manufacturing process and are that part that have not yet escaped from the matrix during the curing process. The thermal and mechanical behavior of the matrix are separated to make the modeling simpler. The thermal behavior is modeled using Arhenius'equation [3, 4] while the elastic behavior is obtained from classical elastic theories. In this section, the elastic behavior of the matrixvoid phase is discussed. The thermal behavior which is related to the viscoelastic behavior of the polymer will be discussed in latter sections. It has been observed that the voids possess a shape that roughly resembles a prolate ellipsoid symmetric about the 1axis [2]. The aspect ratio, 7, of the ellipsoid is defined as the ratio of axis length of the ellipsoid in direction 1 to that along either of the two axes. It is assumed that the voids in the composite have aspect ratios ranging from (~Minto aMax with a mean aspect ratio of ~Ave. The
S. Govindarajan et al. / Finite Elements in Analysis and Desion 23 (1996) 333347
335
voidvolume density at aspect ratio ~ is obtained as: ~gv(a) = c v e(~'~'A'°)' cNo,~ ,
(1)
where c N°~m is the normalizing factor which ensures that the volume fraction density over all the aspect ratio equals c v, i.e. c N°rm =
f
~Max
e ~~A'°)2 d~.
(2)
~Min
The modulus of the matrix, with prolate voids with constant aspect ratio ~ as the inclusions at a concentration of c v, is derived in Appendix B using Eshelby's S tensors [5, 6]. To obtain the modulus of the matrix with the Gaussian void distribution, an incremental doping method is used. This method, as the name suggests, calculates the properties of the matrix phase through an iterative method. During each iteration, voids with aspect ratio ~ and with a total volume computed from Eq. (1) is added to the existing matrix. This process is undertaken until the matrix possesses the necessary distribution. Thus, at the end of i iterations, the modulus of the matrix containing voids with i different aspect ratios is: EMV
V MV = E iMV(~i,c~;i,Ei_),
1 + cgv"
(3)
It is assumed that the matrix at each iteration before doping is homogeneous and has the same properties as that of a matrix with an elastic modulus E ~rv. Therefore, it can be substituted with a material that has the same modulus. 2.2. Composite behavior The woven fiber composite has, in addition to the void and the matrix phases, two more constituents, namely, the weft and the warp fibers. These fibers are woven and their axes of symmetry lie in mutually perpendicular directions. The fiber undulation model [3, 7] has been used to model the geometry of the composite. A representation of the unit cell of this model is shown in Fig. 1. Under the constant strain condition, a stress ~ll, applied to the composite, produces a strain ~ll(X) at any section x equal to ~II(X) = cf(x)Efl(x ) k
cW(x)E'~l + cM(x)E MV"
(4)
The quantities cPhase(x) represents the volume fraction of the phase at that section and can be obtained from geometry. The equations for the weft and warp fibers are derived as [7],
hi(x) =
ht[1 + s i n { ( x  51a ) 7~ } ] I, ht/2
(O<~x<~ao) (a0 <~x~a2) (a2 <~x <<,nga/2)
(5)
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
336
¢ h
1 0
I
s/a
BE
Warp
(ng.O a/2
Fill fibers
fibers
[
MatrixVoid phase
Fig. 1. Unit cell of the fiber undulation model used for the geometry of the composite.
for the weft fiber, while for the warp fiber,
ht/2 ¼hi[1  sin{(x  ½a)~,,}] h2(x)= ]¼hit1  s i n { ( x ½a)~,,}]
( h,/2
( O<<.x<<.ao ) (ao ~ x <~a/2 ) (a/2 ~ x ~ a2)
(6)
(a2 ~ x <.nga/2).
In this case, the matrix phase contains the voids and in Eq. (4), E My with the value obtained from Eq. (3) in the previous section should be substituted.
2.3. Time dependent behavior To model the dynamic behavior of the composites, the individual behavior of each constituent has been considered separately. The matrix has been modeled as a 4parameter model (series combination of the Maxwell and Voigt models). Under creep conditions, it has been observed in many polymers that there is an initial elastic strain, followed by a transient state. This state is characterized by a progressively decreasing strain rate that ends in an almost constant value. Thus, towards the end, the matrix attains a steady state which results in a constant strain rate. The 4parameter model has
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
337
O
rl m Ev
Fig. 2. The 4parameter model as a combination of the Maxwell and Voigt models.
been effectively used to model all the above states [3, 8, 9] and Fig. 2 shows a representation of this model. In the model, the Maxwell elastic element E m determines the initial elastic deformation of the matrix, the Voigt elements E v and qv determine the transient creep state and the Maxwell viscous element r/m gives rise to the steady state creep. From earlier work [3, 8], the following relations are obtained for the matrix alone under the creep condition in the Laplace transformed domain. 1
1
1
~M  E TM
+__
s~m
1
+ _ _ E v + st/V,
fM = vM.
(7)
To obtain the behavior of the composite under creep, Laplace transformation of the matrixvoid phase has to be obtained first. This is achieved by substituting Eq. (7) in Eq. (3). Following these steps, the following relations are obtained. =
1+
(8)
where the initial/~uv is calculated by substituting the transformed values in Eq. (13). Subsequent values are obtained through iteration of the same equation. Eq. (8) gives the Ell for the matrixvoid phase in the transformed domain. The fibers are assumed to possess negligible viscoelastic properties [10] and are therefore considered purely elastic. To obtain the behavior of the composite, the correspondence principle [8,11] is used. From Eq. (4), the following relations are obtained. O'll ~II(X) ~
cf(x)I~I(X) ~ cW(x)Ell
811(X, t) ~  l / l l ( X ) ,
~ cM(x)/~ MV'
811(/) =
where l is the length of the unit cell.
~ l l ( X , t ) d X,
(9)
s. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
338
Eq. (9) can be used to derive the dynamic strain under a time varying stress ~11 in terms of the elastic moduli of the weft and warp fibers and the combined matrixvoid phase and the geometry of the woven composite. For simple loading conditions like creep or fatigue, the Laplace transformation of the stress is readily available. For more complex loading patterns, the convolution integral can be used in conjunction with the impulse response of the composite to determine its strain behavior [12]. Among the constituents, the matrix is significantly affected by temperature variations. The other phases, i.e. the fibers, retain their properties in the range of our interest. Thus, the matrix phase must be studied closely and a suitable model must be adopted to simulate its behavior under different temperatures.
2.4. Thermal behavior The polymer matrix softens with increasing temperature until it fails at its glass temperature (Tg). Within its operating range, Arhenius' equation gives a relation between the ambient temperature and the physical properties of the polymer matrix. Since the elastic and viscoelastic moduli of the matrix represents the resistance of its molecules to the change required by the applied force, they possess an inverse exponential relation with the temperature [3]. Thus,
~M(T) oc e Q/Rr,
EV(T) = E~e Q/Rr,
qV(T)
=
qVoeQ/Rr,
qm(T) = q~e O/gr,
(10)
where Q refers to the activation energy of the polymer resin, T is the temperature in Kelvin, R is the gas constant and E~, r/~ and r/~' are material constants. The fibers are assumed to possess properties that remain relatively constant in the operating range, and hence demonstrate a purely elastic behavior as mentioned in the previous section.
2.5. Dynamic stress in the voidmatrix phase Under constant strain condition, the voidmatrix phase bears a portion of the load, that contributes to its and the composite's viscoelastic behavior, which varies with time. It can be seen that the applied stress is borne by the three phases such that the sum of the loads they bear add up to the applied load. From Eq. 4, the following relation may be obtained for trMv. 1 o'1My 1 (x, t)  cMV(x) [~,1  {cf(x)E[l (x) + cW(x)E~ }g,,(x, t)],
( 11 )
o'~v(t) =
(12)
~l~v(x, t) dx.
Eq. (12) gives the stress in the matrixvoid phase as a function of time. It may be noted that the the stress in both the composites tends down to a zero value at the end, implying that all the applied load is borne by the fiber. This follows from the matrix model and constant strain case. 3. Results and discussion
The above model was used to simulate the behavior of epoxy/woven graphite composite under an applied creep load of 132 MPa (~ll = 132 MPa). The material parameters used for the simulation
S. Govindarajan et al./ Finite Elements in Analysis and Design 23 (1996) 333347
339
are derived from a previous study [1] and are shown in Table 1. For the epoxy matrix, the glass temperature (Tg) is specified as 250 F. Research in molded composites indicate that it is possible to vary the void content in the composite between as much as 1% and 12% and the fiber content in the composite between 50% and 65% [8]. This implies that the void volume content in the matrix varies between 2% and 35%. It has been observed that a typical composite specimen possesses between 2% and 5% void volume fraction, and between 50% and 60% fiber content. The results that follow were determined keeping these variations in consideration. The creep simulations considered the two volume concentrations of voids to observe the effects of its variations. Future simulations would incorporate much larger void volume concentrations. For molded polymer composites, it has been determined that 0~Min and ~Max are 0.05 and 10.0 corresponding to a needle and disc shapes, respectively [9]. The mean value aAve is 5.025. The above values were substituted into the analytical solution and the creep behavior of the composite was determined. Fig. 3(a) shows the variation of EI~v with varying void volume fraction in the matrix (CvM). The initial value (at 0% c v) was obtained from the manufacturer at different temperatures. It is seen that there is an inverse relationship between the void volume fraction and the modulus. It is also seen at all the three temperatures, the modulus tends to a zero value as the void volume increases. This is to be expected since voids have no elastic moduli and thus, a matrix with 100% void content (no resin) will have its modulus equal to zero. The overall modulus of the composite E~I is plotted as a function of its void content at different temperatures in Fig. 3(b). This shows a similar behavior to Fig. 3(a) and has a similar explanation. It is however noticed that the change in the modulus with temperature is not as much as in the previous case, nor does it approach a zero value at higher void content. This is so since the moduli of the graphite fiber is much higher than that of the matrix and dominates the overall modulus of the composite. It remains constant in the operating temperature range, accounting for the only slight variation in the modulus at any void volume fraction. Thus, the composite's modulus would have a nonzero value even when the entire matrix section is removed. The creep strains at two different fiber volume contents 40% and 60% are compared in Figs. 4(a) and (b). It can be seen that the fiber content affects the creep behavior significantly (decrease by nearly 18% with increase in the fiber content). This is to be expected since the fiber is the main load bearing phase in the composite initially and the only load bearing member in the end and any change in its content would be significant throughout the phenomenon. Also, the difference between the two strains also increases as time progresses. It is noted that the difference at the end of 4 h is Table 1 The MaxwellVoigt constants used for simulation Temp.
Fiber
Matrix
EF
Ew
(° C)
(GPa)
(GPa)
22 66 93
240.0 240.0 240.0
22.0 22.0 22.0
vF
0.38 0.38 0.38
E~
?]m
EV
?~v
(GPa)
(GPa s)
(GPa)
(GPa s)
4.51 4.00 3.68
1.17 x 105 4.60 × 104 2.66 × 10 4
12.72 5.00 2.90
3.52 x 104 1.38 × l04 8.00 × 10 3
vM
0.311 0.311 0.311
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
340
E" (GPa) 5.0
4.5
72 F (22 C) ...........
4,0
3.5
"....
150 F (66 C)
......................... 200 F (93 C)
" ,.
3.0
2.5
2.0
1.5 (a)
~ 5
o
Ell
15
10
20
c.V (%)
(GPa) 6O
58 57 56
........... " :" .
55
T=72
54
F(22C)
T=lSOF(66C)
53
...........
52
......................... T = 2 0 0
F(93C)
51
cV(%)
50 (b)
o
5
10
15
2O
Fig. 3. Variation o f (a) matrix's elastic modulus, EllMv and (b) the composite's elastic modulus, EH, as a function of its void content at different temperatures.
S. Govindarajan et aL /Finite Elements in Analysis and Design 23 (1996) 333347
341
~11 I
'
'
'
'
I
'
'
'
'
I
'
'
'
'
I
0.0034
C F = 40% 0.0032
C F = 60%
0.0030
0.0028
0,0026
0.0024
0.0022
I
0,0020
*
L
I
I
L
I
,
I
]
l
l
l
l
0
l
l
L
l
l
l
2
(a)
'
I
'
'
'
'
I
'
'
'
'
3
'Time
I
I
(Hrs)

0.0034
C F = 40% 0,0032
C F : 60%
.......... 0.0030
0,0028
0.0026
0.0024 /
0.0022
0.0020
,
0
=
,
,
I
1
,
,
J
,
I
2
,
~
,
,
I
3
,
~
,
~
I
4Time
(Hrs)
Fig. 4. Creep strain of the composite with varying fiber volume fraction (40% and 60%) and constant c v = 0% and (a) T = 22°C and (b) T = 93°C.
S. Govindarajan et al./ Finite Elements in Analysis and Desion 23 (1996) 333347
342 D
I
'
'
'
'
i
'
'
'
'
I
'
'
'
'
i
,
,
I
0.0034
0.0032
... ..............................  ...............
0.0030
..~
~±,'~L~LL'S',,SL3L,'~,LX±~
3t
0,0028
0.0026
T = 72 F (22 C) ...........
0,0024
T=
150 F(66C)
........................... T = 2 0 0
0.0022
0.0020
,
J
,
~
I
0
(a)
J
,
,
,
I
1
,
F(93C) ,
,
=
I
2
,
,
3
' T i m e (Hrs)
E1 I
'
'
'
'
[
'
'
'
'
I
I
0.0034
0.0032
...... ~ : ~ _ ~ _
0.0030
.......
0.0028
0,0026
T=72F(22C)
0.0024
...........
T=
150 F(66C)
...........................
T=200
F(93C)
0,0022
0.0020
(b)

0
'
~
~
~
I
1
~
~
~
~
I
2
,
~
*
~
I
3
~
~
~
~
I
4 T i m e (Hrs)
Fig. 5. Creep strain o f the composite at different temperatures w i t h c F = 4 0 % and ( a ) ¢v = 0 % and ( b ) c v = 5%.
S. Govindarajan et al./ Finite Elements in Analysis and Design 23 (1996) 333347
343
almost 20%. Since the matrix has a higher volume content in the composite with lesser fiber content, it affects the creep behavior to a greater extent. This manifests itself as an increased creep strain and strain rate at any time. A similar behavior is observed in the simulation at 93°C. Figs. 5(a) and (b) shows the creep behavior of the composite at different temperatures and with c F = 40% and c v = 0% and 5%, respectively. Figs. 6(a) and (b) simulate the behavior at c F = 60% for similar void volume fractions. As the matrix becomes weaker, it loses both its viscous and elastic properties. Thus, the transition state become shorter. It can, however, be noted that the three strains converge to the same value implying that the fiber bears all the stress in the end  a condition specified by the geometric model and the constant strain condition. A similar behavior is encountered in all the cases whatever be the fiber and void volume fractions. Figs. 7(a) and (b) compare the creep strain at two different void volume fractions and at two different temperatures for the same fiber volume content. It can be noted that the lower void content implies a higher modulus (due to more matrix material) and so enables lesser elastic and creep strains at both the temperatures. The transition state is more prominent for the case with more matrix (lesser void content). However, as before, the end strain would tend to the same value in both cases.
m
~11 I
'
'
'
'
I
'
'
'
'
I
'
'
'
'
I
0.0034
T=72
0.0032
F (22C)
...........
T=
150 F(66C)
.........................
T =200
0.0030
F(93C)
0,0028
0.0026
0.0024
0.0022
0.0020
(a)
~
0
,
,
~
I
1
~
~
~
~
I
2
~
,
~
~
I
3
~
~
~
~
L
4Tim e (H rs)
Fig. 6. Creep strain o f the composite at different temperatures with c F = 60% and (a) c v = 0% and ( b ) c v = 5%.
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
344 D
Ell I
'
'
'
'
I
'
'
'
'
I
I

0.0034
T = 72 F (22 C)
0.0032
...........
T=
150F(66C)
.........................
T=200
0.0030
F(93C)
0.0028
0.0026
0.0024
..~.~.2~2'2_
............
0.0022
,
0.0020
J
,
,
I
0
,
,
,
,
I
,
1
,
,
,
I
,
,
t
2
,
3
I
4Time (Hrs)
(b) Fig. 6. Continued
I
'
'
'
'
1
'
'
'
'
1
'
,
,
'
'
'
1
0.0034
0.0032
0,0030
0.0028
0.0026
cV= 0% 0.0024
V
c =5%
0.0022
0.0020
(a)
,
0
,
,
,
k
1
,
,
,
,
I
2
[
,
I
3
,
,
,
,
i
4Time (Hrs)
Fig. 7. Creep strain of the composite with varying void volume fraction (0% and 5%) and constant c F = 40% and (a) T = 22°C and (b) T = 93°C.
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
345
~11 I
0.0034
0.0032
0.0030
0.0028
0.0026
cV= 0% 0,0024
Cv .
.
.
.
.
:5%
0.0022
0.0020
(b)
,
o
,
,
,
I
1
L
,
,
,
I
~
,
,
2
,
I
3
,
,
,
,
[
4Time (Hrs)
Fig. 7. Continued
4. Future work The proposed work extending this model to predict the void growth phenomenon in the matrixvoid phase using the dynamic stress variation at different ambient temperatures extending up to the Tg of the matrix. Experimental investigation will also be carried on composites with different constituent ratios to obtain the experimental validation of the analytical model. This model will then be extended to predict crack initiation and growth in the composite under dynamic loading conditions.
5. Conclusions A thermomechanical model to simulate the creep behavior of a polymer woven composite has been presented. This model considers all the phases that are present in a typical composite and considers the different behavior of each constituent to obtain the overall response of the composite. Also, a method to derive the stress field in the polymer matrix has been presented. The stress values thus obtained would aid in predicting the behavior of the resin matrix near the voids and be able to predict the crack initiation phenomenon in the matrix.
Acknowledgements This research was made possible by a grant from the FAA Center of Excellence through the Center for Computational Methods in Aircraft Structures at the Rutgers University.
346
S.
Govindarajan et al./ Finite Elements in Analysis and Design 23 (1996) 333347
Appendix A. Nomenclature A.1. Symbol convention
In modeling the time dependent characteristics of the composite, the following convention has been assumed. For a first order vector, a subscript denotes the axis of interest. For a tensor, the subscripts denote the indices. An uppercase superscript denotes a particular phase, while a lower case superscript denotes a component or constituent of the phase. An overbar denotes that the quantity considered is for the entire composite and not for any specific phase. A hat (^) denotes the quantity in Laplace transformed domain.
A.2. Symbols
a e q a c c~ c Mv ¢g E S
aspect ratio of the ellipsoid strain viscosity stress volume concentration/fraction volume fraction of the void in the matrix phase volume fraction of the matrixvoid phase volume fraction density elastic modulus Eshelby's S tensor
Superscripts
f F m
M V
V W
weft fiber phase fiber phase (combined weft and warp phase) Maxwell model component matrix phase Voigt model component void phase warp fiber phase
Appendix B: Modulus of a matrix with eilipsoidal voids The following relations are obtained from literature [6]. These serve to derive the modulus Ell of an isotropic matrix material containing prolate ellipsoidal voids of constant aspect ratio ~ at a volume fraction of c v. EM
E ~ ( ~ ' e V ' E M ) = 1 + cV(Al + 2vMA2)/A '
(13)
S. Govindarajan et al./Finite Elements in Analysis and Design 23 (1996) 333347
347
where A1 = DI(B4 + Bs)  2B2,
A2 = (1 + Dr)B2  (B4 q Bs).
(14)
The various geometric and elastic terms in Eq. (14) are obtained as B3 c v q D3 + (1  cV)[s1 Ill + (1 + D1 )$2211]
B4 = cVD1 + 02 q (1  cv)($1122 } DIS2222 + $2233) B5  c v q D3 q (1  cV)(Sl122 k $2222 ] D1S2233)
(15)
where Si#t is the Eshelby's S tensor and is a function of ~ and the material properties. For a matrix with inclusions with zero Eij, we obtain D1 = 1 + 2#M/2 M,
D2 =  D I ,
D3 =  1 .
(16)
References [1] A. Farouk and T.H. Kwon, "Effect of processing parameters on compression molded pmr15/c3k composites", Polymer Composites 11(6), pp. 379386, 1990. [2] A. Farouk, N.A. Langrana and G.J. Weng, "Modulus prediction of crossply fiber reinforced fabric composite with voids", Polymer Composites 13(4), pp. 285294, 1992. [3] S. Govindamjan, N.A. Langrana and G.J. Weng, "Modeling creep behavior in polymeric woven composites", in: Advances in ComputerAided Engineering (CAE) of Polymer Processing, pp. 218296, ASME, 1994. [4] J.D. Ferry, Viscoelastic Properties of Polymers, Wiley, New York, 1961. [5] J.D. Eshelby, "The determination of the elastic field of an ellipsoidal inclusion, and related problems", Proc. Royal Soc. London A 241, pp. 376396, 1957. [6] G.P. Tandon and G.J. Weng, "The effect of aspect ratio of inclusions on the elastic properties of unidirectionally aligned composites", Polymer Composites 5, pp. 327333, 1984. [7] T. Ishikawa and T.W. Chou, "Stiffness and strength behavior of woven fabric composites", J. Mater. Sci. 17, pp. 32113220, 1982. [8] Y.M. Wang and G.J. Weng, "The influence of inclusion shape on the overall viscoelastic behavior of composites", ASME J. Appl. Mech. pp. 510518, 1992. [9] R.N. Yancey and M.J. Pindera, "Micomechanical analysis of the creep response of unidirectional composites", J. Eng. Mater. Technol. 112, pp. 157163, 1990. [10] W.N. Reynolds, Physical Properties of Graphite, Elsevier, Amsterdam, 1968. [11] R. Bellman, R.E. Kalaba and J.A. LockeR, Numerical Inversion of Laplace Transform: Applications to Biology, Economics, Engineering and Physics, Elsevier, New York, 1966. [12] N.A. Langrana and S. Govindarajan, "The application of the convolution integral to simulate a woven composite's behavior under random loading", Technical Report 957, FAA Center of Excellence for Computational Modeling of Aircraft Structures, Rutgers University, 1995. [13] G.J. Weng, "Some elastic properties of reinforced solids, with special reference to isotropic ones containing spherical inclusions", Internat. J. Eng. Sci. 22 (7), pp. 845856, 1984.