Mini-constitutive finite element model for plastic response of unidirectional fiber composites

Computers & S~rucrures Vol. 55. No 3. pp. 463470. 1995 Elsewer Science Ltd Printed tn Great Britain 0045-794919s 169.50+ 0.00

00457949(94)00477-3

MINI-CONSTITUTIVE FINITE ELEMENT MODEL FOR PLASTIC RESPONSE OF UNIDIRECTIONAL FIBER COMPOSITES L. Parietti, S.-Y. Hsu and 0. H. Griffin, Jr Department

of Engineering

Science and Mechanics, Virginia Polytechnic Blacksburg, VA 24061, U.S.A. (Receiced 18 December

Institute

and State University,

1993)

Abstract-A micromechanical finite element model to compute the overall instantaneous stiffness of fiber-reinforced composites in elastoplastic response is presented. This model is applicable to a periodic diamond array of elastic circular fibers embedded in an elastic--plastic matrix subjected to a plane stress loading. This model enforces symmetry and anti-symmetry conditions isolating the smallest unit cell and micromechanics within a larger finite element should greatly increase the speed of doing “built-in” program, because of the small number of degrees of freedom (12-14 DOF). At this stage of development, effective properties of a composite with a nearly incompressible matrix are presented to assess the model’s plasticity performance. Comparison with a fine grid finite element solution shows very good results and demonstrates the effectiveness of the mini-model presented.

I. INTRODUCTION

AND

variational method. In this model, the composite microstructure is modeled as a system of cylindrical fibers embedded in a continuous matrix phase. Later, Dvorak and Rao [IO] used the CCA model to derive closed-form solutions of the elastic-plastic response of fiber composites under axisymmetric loads. Comparison with finite element plastic solutions shows the accuracy of their results, but their model is restrictive because of the nature of the loads applied. Another popular model which was also first used for the elastic range and then applied to elasticplastic analysis is the self-consistent model (SCM). First devised by Hershey [I I] as a means to model the behavior of polycrystalline materials, the SCM, which is based on the solution of an elastic ellipsoidal inclusion problem [12], was extended later to multiphase media [l3, 141. This model is very simple. but it takes great liberties with the geometry and gives inaccurate results for some quantities. When a variant of the self-consistent scheme is used to study elastic-plastic response of fibrous composites. it overestimates the initial yield stress [I 51 and underestimates plastic strains in the early stages of deformation [ 161. Dvorak and Bahei-El-Din [I71 used a micromechanical model consisting of elastic filaments of vanishingly small diameters in an elastic-plastic matrix. Analytical expressions for overall compliances of the composite in terms of the phase properties and the volume fraction were derived. However, at higher fiber volume fractions, this model provides a rather low estimate of the flow stress in transverse tension response.

BACKGROUND

Even though overall elastic behavior of fiber-reinforced composites is reasonably well understood, elasticPplastic behavior is still under investigation. Fibrous composites exhibit significant non-linear stress-strain behavior under off-axis loading. This non-linearity is mainly due to the matrix, while the fibers are in general linearly elastic up to rupture. Studying this non-linear behavior is crucial since the strain at which many common matrix materials deform plastically is well below fiber failure. Therefore, a good design must exploit the fact that a composite can carry loads even after yielding has occurred. Several researchers have investigated the inelastic behavior of unidirectional fiber composites. Two different approaches have been used: a macromechanical and a micromechanical approach. In the macromechanical approach, the composite is regarded as a “mixture” of two different materials and approximated as a macroscopically homogeneous medium. Although this approach is usually simpler, it often relies heavily on experimental data. Examples of such approaches are given in Refs [l-8]. On the other hand, the micromechanical approach regards the composite as a non-homogeneous structure. It allows more complicated composite behaviors to be modeled and reduces testing needs by taking internal iteraction into account. Several analytical models were developed first. The composite cylinder assemblage (CCA) model was introduced by Hashin and Rosen [9], who derived bounds and expressions for the effective elastic moduli of fiber-reinforced composites using a 463

464

L. Parietti er rrl.

Aboudi [IS] analyzed a square unit cell consisting of four square subcells (a fiber subcell surrounded by three matrix subcells),in each of which uniform stress and strain are presumed; the fibers are assumed to be arranged throughout a periodic square model (PSM). This model can actually be viewed as an assembly of four non-conforming linear triangular finite elements [ 191, and thus much more concise formulation and programming than those found in the literature can be achieved. Moreover. Aboudi’s model leads to an indeterminate displacement field when the applied load involves transverse shear [20], which is ascribable to the non-conforming elements. Aboudi further derived a system of closed-form constitutivc equations [2l]. which has been applied by Pindcra (~1 (11.[22] to simulate inelastic behavior of a boron/aluminium composite. and by Arenburg and Reddy [23] to analyze metal-matrix laminates. However, the use of the closed-form relations is not numerically more advantageous than the usual finite element procedure for the simple triangular elements. In addition. the system of closed-form equations is not suitable for the tangential stiffness method, either at the micromcchanical Icvel or at the macromechanical level. More recently, Sun and Chen [24] developed a micromechanical model consisting of three rectangular subcells (a fiber subcell and two matrix subcells) to calculate off-axis responses of fiber composites. This model is very simple and can be implemented easily, but. like Aboudi’s model, it does not allow modeling of matrix property gradients. Besides. both models are only applicable to a periodic square array of fibers, and it is known that significant differences in inelastic predictions exist bctwcen the periodic square (PSM) and the hexagonal (PHM) models [l9]. The PSM seems too compliant for off-axis responses at small fiber angles and too stiff at large fiber angles. Therefore, no matter how the matrix material parameters are adjusted, it is impossible to obtain reliable results for both cases. This suggests that the PHM is more suitable for modcling internal constraint hardening than is the PSM. Finally. some authors [25-271 have used the finite element method in a conventional way to analyze unit cells of unidirectional composites in the elastoplastic range. If effective behavior of composites is of primary interest and computational eficiency is a serious concern. modification of the conventional finite elcmcnt implementation is necessary to meet the following requirements: (I) the smallest number of degrees of freedom (DOF) used for a smallest unit cell; (2) efficient imposition of micromcchanical boundary conditions: (3) cficient imposition or computation of macromechanical variables such as average strains and stresses. The periodic hexagonal array (PHA) model developed by Tcply and Dvorak [28. 291 seems to bc the only mini-finite element model that satisfies the above requisites. The PHA model is capable of predicting the three-dimensional

instantaneous stiffness of a composite. and uses 18 DOF for three-dimensional strain-controlled simulation. However, the model approximates the fibers as hexagonal cylinders: we will show that the geometry of the fibers has a great influence on effective inelastic behavior of composites. In particular, significant differences in effective propcrtics exist between models using hexagonal and circular fibers. Moreover. further reduction for the PHA model is needed to promote its computational performance for the important plane stress loading. but the dcrivation to the current form of the model is already lengthy [28]. The objective of this paper is to present an cfticicnt composite model for prediction of eff‘cctive behavior of unidirectional fibrous composites under plane stress loading. Recently. a micromechanical finite element model involving 318 DOF has been used by Hsu [I91 for elastoplastic analysis. The finite clement analysis is performed on the smallest unit ccl1 for a periodic diamond array of circular fibers by cnl‘orcing symmetry and anti-symmetry boundary conditions under plane stress loading. The periodic diamond model with circular fibers is used for better modeling of overall inelastic deformation, which is influenced by the internal constraint. Thus. the finite element method is adopted for its geometric modeling capability. Good agrccmcnt has been demonstrated bctwcen the predictions of Hsu’s model and experimental results provided bq Bcckcr rt rll. [30] for various elastoplastic characteristics. Thel-cfore. Hsu’s model will be used as a reference, and the aim of this work is to develop I mini-model involving considerably fewer degrees of freedom (almost 30 times less), but capable of predicting the same cffcctive properties. The number of degrees of freedom of the mini-model is expected to bc convcnicntlh adjusted with respect to composite deformation characteristics, which arc determined by factors such as differences m constituent stiffncsscs, inconlprcssiblc matrix plasticity. and radial property gradients 111the matrix. To meet the aforcmcntioncd rcquircmcnts fog mini-models. a specialized tinitc clement Implcmcntation is ncedcd. Nevertheless. we attempt to keep our formulation simple and its computer implcmcntation convenient.

As mentioned earlier, our goal is to dc\elop a mim-model to predict the effcctikc propcrtics of a unidirectional fiber-reinforced lamina under plane stress loading. As shown in Fig. I. the lamina III the plant .Y,-.vz is subjected to plane stress loading. The fibers arc aligned in the .I-, direction and assumed to be arranged in a periodic diamond array. Both the fibers and the matrix are assumed orthotropic; the fibers are elastic. while the matrix is elastoplastic. Pcrfcct bonding is assumed between fiber and matrix.

Mini-constitutive finite element model

465

mogeneity at the microscopic field u is thus:

level. The displacement

u=ii+Xc”,

Fig. 1. Plane stress loading conditions and axis system.

where X is a matrix containing the spatial coordinates (.Y,, x2, x3) of the point where the displacement field is evaluated, i.e.

2.1. Unit cell The geometry of the array is defined by the fiber diameter df, and the fiber spacings h and e, as shown in Fig. 2. Because of the assumed periodicity, a typical repeating hexagonal volume can be isolated as indicated by the solid lines in the figure. Furthermore, because of the symmetry about AB and BC, only one quadrant ABCDM of the repeating unit volume needs to be analyzed. This quadrant constitutes our unit cell. It represents by itself the entire composite if proper boundary conditions are applied (see Section 2.3). With M defined as the center of the rectangle BCB’C’, the sloping edge AMD is defined by the angle x, which can take any value as long as AMD lies entirely in the matrix phase. 2.2. Displacement .$eld decomposition At the macroscopic level, the composite is considered as a homogeneous orthotropic material undergoing deformations represented by a macroscopic strain field,

C”=jf;),,~r~.tEi,2C,?,,T . 0

(1)

consistent with the plane stress loading conditions. The gradients of that macroscopic deformation field are assumed to be small so that L” can be considered constant across several unit cells. With this assumption, the micromechanical displacement field u in the unit cell can be seen as the superposition of two fields. One field corresponds to the macroscopic strain field co, and the other represents a perturbation fi from that state of uniform deformation due to the nonho-

(2)

x=

s,

0

0

.Y?

0

_Y2

0

0.

0

0

xi

0

I

(3)

The perturbation field ii will be referred to as the fhctuuting displacement,field. Since the fibers are very long, the fluctuating displacement field does not depend on .Y,, the direction of the fibers. 2.3. Boundar.:r conditionA Boundary conditions are imposed on the fluctuating displacements fi and the tractions t on the sides of the unit cell to account for the periodicity, symmetries and antisymmetries of the geometry and deformations. They are discussed in detail in Ref. [19]. We only give here a brief description. The fluctuating displacement components ti, and zi? are anti-symmetric about the Y,~_Y~planes at x2 = 0 and .Y?= e. so ~7,= & = 0 on the sides AB and DC of the unit cell (see Fig. 3). The traction component t, is also zero along these sides. The .Y,-.Y: plane at .Y?= 0 is a plane of anti-symmetry for L&, so 6, is equal to zero on BC. The tractions t, and tz along this side are also equal to zero. Finally, the fluctuating displacements are anti-symmetric with respect to the midpoint M of the sloping edge of the unit cell and therefore vanish at M. The tractions t are symmetric with respect to M on the sloping edge. 2.4. Finite element ,fbrmulation Using a displacement formulation, the fluctuating displacement field ii within element number e is interpolated from element nodal degrees of freedom

Fig. 2. Unit cell.

L. Parietti

466

et al

which 4

___I

i

t _AM=t,p

)

to be zero. eqn (6) can be rewritten

After

these

___.

where

{oy,. IT;', aTz}' of macroscopic

fO =

It;), C(;?, We can now eliminate ments the effective as: CT”=

t, = t*= 0 Fig. 3.

displacerelation

(8)

‘

The matrix C, can be identified to the stiffness matrix of an orthotropic material and engineering moduli of the composite can be determined. 3. REFERENCE

MODEL

As mentioned earlier, the 318 DOF finite element grid introduced by Hsu [I91 will be used to evaluate the performance of our mini-model. The two-dimensional finite element mesh used in Ref. [I91 is shown in Fig. 4. It consists of 216 constant-strain triangles and 126 nodes. Three variables which represent the fluctuating displacement field are allocated at each node. By applying the boundary conditions described in Section 2.3, the number of degrees of freedom reduces to 318. 4. INFLL’ENCE

OF INTERNAL

GEOMETRY

In previous works, simplified geometries such as a square geometry [ 18,20-241 or an hexagonal geometry [28,29] have been used to represent the fibers. However, significant differences in effective properties exist between models using square, hexagonal and

where {e;‘, scopic plane strain field

is the vector associated total cell, of nodal associated displacement field matrix nodal K: is the coupling between and macrofields. of the other edge, off associated components On the sloping edge, displacement field and the traction about M. Therefore, of the numbers and degrees in eqn (6) can be performed nodes about edge. reduction for nodal displace-

Fig. 4. Reference grid.

Mini-constitutive finite element X3

0 hierarchic

t

DOF

Fig. 5. Mini-grid.

circular fibers, especially when modeling plastic behavior. As an illustration, effective properties of a boron/aluminium composite have been computed with two different fine grids, one using hexagonal fibers and the other circular fibers (reference grid). The cross-sectional area of the fiber is the same for both grids. In order to simulate a plastic behavior, the properties of the matrix were modified. The matrix was made almost incompressible by adjusting its Poisson’s ratio to 0.495, while keeping its bulk modulus unchanged. The modification of the Poisson’s ratio results in a change in the ratio of the bulk modulus to the shear modulus from 2.6 to 100. A difference of 42% was observed between the two geometries on the prediction of the composite Young’s moduli in the direction x1 and xj defined in Fig. 4 (/I$ and ,I&). This difference indicates that internal geometry affects the effective properties of the composite when plastic response is analyzed. The mini-model developed here uses the more realistic circular geometry to represent the fibers. 5. MINI MODEL

Our goal is to develop a mini-model capable of predicting the same material properties as the reference grid presented in Section 3. The mini-model uses the same unit cell, the same superposition method for the displacement field, and the same boundary conditions as the reference model. The mini-grid is shown in Fig. 5; it contains four elements. The matrix is divided into three hierarchic [3l] elements I, II and III, as shown in Fig. 5. Hierarchic

/ , id

i,” = 10, ip, iilr 0, 0, 2i,?JT.

Ul

,0’

Nl ----____

Fig. 6. Linear

elements

,,/ a

By applying the boundary Section 2.3 and fluctuating

UJ

,_.:‘.

conditions described in displacement continuity

a3 ,c’ _.--*

/ //

u2

with standard

(9)

The remaining non-zero components are the 3 DOF allotted to the fiber, i.e. t;?, E,, and 2i,? (see Fig. 7). The associated fluctuating displacement field 8,” inside the fiber is a linear function of x2 and _yj and is written as:

"'NZ

----I-/

467

elements are isoparametric bilinear elements with quadratic corrections added along the sides. The resulting elements are either subparametric or isoparametric (distorted elements). The additional degrees of freedom associated with the quadratic hierarchic shape functions measure the departure from linearity at the midpoint of each side of the element. Examples of simple, one-dimensional linear elements with standard shape functions and with quadratic hierarchic shape functions are shown in Fig. 6. For the linear element with standard shape functions, the displacement field u = N, U, + N,uz is linear, while the displacement field for the linear element with quadratic hierarchic shape function u = N, u, + NZuz + N,a, is quadratic. The degree of freedom a, measures the amplitude of the quadratic correction. Note that the use of hierarchic shape functions instead of regular serendipity shape functions allows us to switch from quadratic to linear approximation in a particular direction by simply setting the associated quadratic correction to zero. For this reason, the use of “hierarchic elements” proved very valuable in the process of minimizing the number of DOF in our model. The fiber is represented by element IV, which has a circular geometry. The,fluctuating strain field inside the fiber C,, = {i,, , iz2, ix3, 2i2), 2C,, , 2C,2}Tis assumed to be uniform. The boundary conditions described in Section 2.3 along the sides x2 = 0 and .xj = 0 of the fiber lead to i,, = CZ3= 0. Also, i,, is equal to zero since the fluctuating displacement field does not depend on xi. Thus the fluctuating strain field inside the fiber reduces to

,’

I

model

.,..-.. ‘. . .N3

“2

‘.

3

and hierarchic

shape

2

functions.

L. Parietti

468

er ul Table

I. Elastic effective boron/aluminum

properties material composite

3

Fine grid 318 DOF

Mini-grid 18 DOF

E, (GPa)

223.29 139.97 139.97 54.19 0.2637 0.2637 0.3557

223.48 140.48 140.40 54.19 0.2637 0.263X 0.3539

NJ-

for

ti

&(GPa) E,(GPa) G,,(GPa) “I? VII \‘i-

_)

X2

0 Fig. 7. Degrees

hierarchic

of freedom

V, = 46%. E, = 400 GPa.

\lr= 0.2, E,,, = 72.4 Gpa, v,,, = 0.33.

6.2.

rurzg~

DOF

of the mini-grid

(IX DOF).

at the seven nodes of the fiber-matrix interface, the number of DOF reduces to 18. The remaining DOF are shown in Fig. 7 and are as follows: one regular DOF represented by a bullet, 13 hierarchic DOF represented by circles, three normalized fluctuating strains in the fiber and the unknown macrostrain co13 When plasticity dominates, the matrix behaves nearly incompressibly. If too few degrees of freedom are used to model the deformation of the matrix. it may not deform at all. This phenomenon is called plastic locking and it is usually avoided by adding numerous DOF and using a reduced numerical integration rule. Our first concern was to develop a mini-model immune to plastic locking; that is the reason why quadratic corrections were added along every side of the elements discretizing the matrix. However. as already mentioned. minimizing the number of DOF is also desirable. WC will see in Section 6.3 that not all of the 13 quadratic corrections of Fig. 7 are needed to prevent plastic locking and linear displacement approximations in some directions are sufficient along some sides of the elements representing the matrix. conditions

6. RESULTS

The mini-model is used to predict effective properties of a boron/aluminium composite. The geometry of the composite is defined by a volume fraction ~‘t= 0.46. a fiber diameter dr = 0.142 mm and a fiber spacing 11= 0.173 mm as defined in Fig. 2 (from which c = 0.0995 mm). This geometry is very close to an hexagonal array. The angle z defined in Fig. 2 is set to 30 The elastic properties of the fibers and the matrix are: Etlher = 400 GPa, tlllhcr= 0.2, En,,,,,, = 72.4 GPa, and vnlJlrly= 0.33 The effective properties predicted by the mini-model are given in Table I and compared to those predicted by the reference model. Very good agreement is observed. The largest error is about 0.5% in the value of tlj2

Sirnuluted

plustic

The mini-model is intended to provide reliable numerical results when a certain plasticity theory representing the matrix deformation is incorporated into the composite model. At this stage of development, elastoplastic characteristics arc not yet available. But as mentioned in Section 5, difficulties such as Poisson’s ratio locking may arise in the plastic range. In our case, this may be crucial because of the very small number of DOF used. In order to test the immunity of our model to plastic locking without involving any plasticity theory or computation, we chose to simulate the effective behavior of the boron/aluminum composite with the elastic properties of the matrix modified as described in Section 4. while the fiber properties remain unchanged. The modified Young’s modulus is E,,,,,,,, = 2.12 GPa, and the modified Poisson’s ratio is v,,,,,,~~, = 0.495. which corresponds to a ratio of the bulk modulus to the shear modulus equal to 100. The geometry of the periodic diamond array and the fiber properties are as defined in Section 6.1. The effective properties predicted by the mini models are computed using a 2-by-2 Gauss rule for the elements I, II and III of Fig. 7 and are compared once again with predictions of the reference model (see Table 2). Very good agreement is observed. The largest error is about 3% in the value of Ez and E,. The results obtained indicate that locking has been avoided, dcspitc the

Table 2. Simulated plastic effective material boron/aluminum composite 3

J--

Fine grid 318 DOF

properties

for

Mini-grid 18 DOF

25

E, (GPa)

E,(GPa) E; (GPa) G‘,I(GPa) I’ll “I 1

I’:,

\‘, = 46%. I’,,,= 0.495.

184.97 8.67 8.67 1.91 0.3589 0.3589 0.9478 E, = 400 GPa.

V,= 0.2.

185.21 8.93 x.93 I.00 0.3585 lJ.3590 0.9464 E,,, = 2. I2 GPa.

Mini-constitutive

n

Set

finite element

properties

E,,

,

, V 32

s.t. c(QQF)is MINIMUM

Fig. 8. Optimization

very small number of DOF. Note that the ratio of the bulk modulus to the shear modulus which is characteristic of inelastic deformation of matrix materials is usually lower than 100, so that our mini-model should be immune to Poisson’s ratio locking in most practical cases. 6.3. Further minimization

469

DOF to zero

Compute equivalent

Select MF

model

of model size

As previously discussed, performing a plastic analysis of a composite can be very costly, and even a small reduction of DOF can be significant for minimizing computational cost. This can be achieved by removing one by one the DOF which is least needed to mode1 the fluctuating displacements (i.e. produces the smallest deterioration in simulated plastic effective properties when removed from the model). The algorithm is described in Fig. 8. As an illustration, the number of DOF for a boron/aluminum composite with the characteristics defined in Section 6.2 can be reduced from I8 to I2 if a tolerance of 8% on the effective properties is chosen. During the first stage of the optimization process, when going from 18 to 14 DOF, the errors tend to become evenly distributed among the effective

process.

properties. All quadratic corrections in the X, direction (number 4, 7, IO and I3 in Fig. 7) are removed during this phase and decrease the accuracy of G,?, while the errors on other engineering moduli remain essentially unaffected. Then, when the number of DOF is further reduced from I4 to II, the quadratic corrections I. 6 and 3 (see Fig. 7) are successively eliminated and the error on E2 increases up to 9.25%, while the error on Glz remains unchanged. The evolution of the relative errors on E2 and G,, during the optimization is shown in Fig. 9. In comparison to typical errors introduced by existing plasticity models, an error of 5% seems acceptable and the I2 DOF mini-model will be used for future elastoplastic analyses. The effective properties computed by the 14 and the I2 DOF models are given in Table 3. The maximum errors with respect to fine grid are about 3.1% on G,z for the model containing I4 DOF and 4.6% on Ez. E, for the I2 DOF model. Note however that the optimal reduction of DOF depends on the nature of the composite (both its geometry and the properties of its constituents). Therefore, this simple optimization process could be performed for each composite system before attempting any macro+micro integrated analysis. Table 3. Simulated plastic effective material properties for boron,/aluminum

II--+--_-

a

2

’\ .

-

E, (GPa)

.L.

\c -1

0 IO

I

I

11

12

1

1

I

13

14

15

Number

(

16

.\& 17

18

of DOF

Fig. 9. Relative errors on simulated plastic effective properties vs the number of DOF for boron/aluminum.

E,(GPa) & (GPa) G,:(GPa) “I? 1’1 1

viz 11,= 46%. I’,”= 0.495.

composite

predicted erids

by reduced

mini-

Fine grid 318 DOF

Mini-grid 14 DOF

Mini-grid 12 DOF

184.97 8.67 8.67 1.91 0.3589 0.3589 0.9478

185.21 8.93 8.93 1.97 0.3585 0.3589 0.9464

185.24 9.07 9.07 1.97 0.3606 0.3564 0.9452

E, = 400 GPa.

c, = 0.2,

E,,, = 2.12 GPa.

L. Parietti el al

470 7. CONCLUDING

REMARKS

A mini-model involving very few DOF (18) has been presented. It is capable of predicting effective elastic properties of unidirectional fiber-reinforced composites under plane stress loading with excellent accuracy (0.5%). The small number of DOF makes the new model a valuable tool for performing macro-micro integrated plastic analysis. Its immunity to plastic locking was demonstrated by obtaining effective properties with incompressible matrix material. Its number of DOF can be further reduced by a simple optimization process. In view of current results, good behavior in the plastic range can be expected

and

is currently

under

investigation.

Acknowledgement-S.-Y. Hsu is supported Virginia Institute for Material Systems.

in part by the

REFERENCES

1. P. H. Petit and M. E. Waddoups. A method of predicting the nonlinear behavior of laminated composites. J. Composite Mater. 3, 2-19 (1969). 2. H. T. Hahn, Nonhnear behavior of laminated composites. J. Composite Mater. 7, 257~-21 I (1973). 3. 2. D. Hashin, D. Bagchi and B. W. Rosen, Nonlinear behavior of fiber composite laminates. NASA CR-2 I3 (1974). 4. R. S. Sandhu, Nonhnear behavior of unidirectional and angle-ply laminates. J. Aircrqft 13, 10441 I1 (1974). 5. R. M. Jones and H. S. Morgan, Analysis of nonlinear stress-strain behavior of fiber-reinforced composite materials. AIAA J. 15. 1669-1676 (1977). 6. S. Amijima and T. Adachi. Nonlinear stress-strain response of laminated composites. J. Composite Mater. 13, 2066218 (1979). 7. 0. H. Griffin Jr, M. P. Kamat and C. T. Herakovich, Three-dimensional inelastic finite element analysis of laminated composites. J. Compo.yite Ma/w. 15, 543 560 (1981). 8. M. N. Nahas, Analysts of nonlinear stress-strain response of laminated fiber-reinforced composites. Fiber Sci. Technol. 20, 297 3 13 (I 984). 9. Z. Hashin and B. W. Rosen, The elastic moduli of fiber-reinforced material. J. Appl. Mrch. 31, 223-232 (1964). IO. G. J. Dvorak and M. S. Rao, Axisymmetric plasticity theory of fibrous composites. In/. J. Engng Sci. 14, 361-373 (1976). I I. A. V. Hershey. The elasticity of an isotropic aggregate of anisotropic cubic crystals. J. Appl. Mech. 21,236- 240 (1954). 12. J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. R. Sot. Land. 241, 376 396 (1957).

13. R. Hill, Continuum micro-mechanics of elastoplastic polycrystal. J. Mech. Phys. Sol. 13, 89-101 (1965). 14. B. Budiansky, On the elastic moduli of some heterogeneous materials, J. M&I. Phys. Sol. 27, 5 I 72 ( 1979). 15. G. J. Dvorak and Y. A. Bahei-El-Din, Elastic-plastic behavior of fibrous composites. J. Me&. fhys. Sol. 27, 51-72 (1979). 16. J. W. Hutchinson, Elastic-plastic behavior of polycrystalline metals and composites. Proc. R. Sot. Land. 319, 2477272 (1970). 17. G. J. Dvorak and Y. A. Bahei-El-Din, Plasticity analysis of fibrous composites. J. Appl. Mech. 49, 327 355 (1982). 18. J. Aboudi, Elastoplasticity theory for composite materials. So/id Mech. Arch. 11, 141 -183 (1986). 19. S.-Y. Hsu, Finite element micromechanics modeling of inelastic deformation of unidirectionally fiber-reinforced composites. Ph.D. Dissertation, Virginia Polytechnic Institute and State University (1992). 20. M.-J. Pindera and J. Aboudi, Micromechanical analysis of yielding of metal matrix composites. Inr. J. Plus. 4, 195-214 (1988). 21. J. Aboudi, Closed form constitutive equations for metal matrix composites. In/. J. Engng Sci. 25, 1229%1240 (1987). 22. M. J. Pindera, C. T. Herakovich. W. Becker and J. Aboudi, Nonlinear response of unidirectional boronialuminum. J. Compos. Mater. 24, 2-~21 (1990). 23. R. T. Arenburg and J. N. Reddy. Analysis of composite metal -matrix structures--II. Laminate analyses. Comput. Struct. 40, 1369 -1385 (1991). 24. C. T. Sun and J. L. Chen, A micromechanical model for plastic behavior of fibrous composites. J. Compos. Sci. Technol. 40, 115-129 (1991). 25. D. E. Adams, Inelastic analysis of a unidirectional composite subjected to transverse normal loading. J. Compos. Mater. 4, 310-328 (1970). 26. R. L. Foye, Theoretical post-yielding behavior of composite laminates-I. Inelastic micromechanics. J. Compas. Mater. 7, 178-193 (1973). 27. S. Jansson, Mechanical characterization and modeling of non-linear deformation and fracture of a fiber reinforced metal matrix composite. J. Mech. Mater. 12, 47-62 (1991). 28. J. L. Teply and G. J. Dvorak, Dual estimates of overall instantaneous properties of elasticcplastic composites. Proc. 5th Int. Symposium on Continuum Mod& of Discrete Slstems (edited by A. J. M. Spencer). Nottingham, pp. 205-216 (1985): 29. J. L. Teolv . _ and G. J. Dvorak. Bounds on overall instantaneous properties of elastic plastic composites. J. Mec,h. Phys. Sol. 36, 29-58 (1988). 30. W. Becker, M. J. Pindera and C. T. Herakovich. unidirectional boron; Mechanical response of aluminum under combmed loading. Technical Report VPI-E-81-17, College of Engmeering. Virginia Polytechnic Institute and State University (1987). 3 I. 0. C. Zienkiewicz and R. L. Taylor, The Finite &men/ Method. 4th edn. Section 7.2. McGraw-Hill. Maidenhead ( 1989).