Computation of the Effective Properties of Materials with Microstructure

Effective material properties are computed for materials with microstructure. Two types of materials are considered. The first type of material studied is a discontinuously reinforced composite material. This type of material has a microstructure that consists of stiff, brittle particles surrounded by a compliant matrix. A micromechanical model is utilized to estimate the elastic-plastic deformation characteristics of these composite materials. The model consists of an elastic-plastic matrix material reinforced with a dilute concentration of axisymmetric elastic inclusions. The model is based on the dilute approximation which models the composite material as an inclusion of finite size embedded in an infinite material. The matrix material is assumed to obey J2 flow theory with isotropic hardening. The finite element method is utilized to solve the nonlinear boundary value problem. Results are presented which illustrate the effect of reinforcement shape on the evolution of the local matrix state variables. A volume averaging scheme is then employed to obtain dilute estimates for the macroscopic response. Results showing the effect of reinforcement shape and volume fraction on the predicted stress-strain response of the composite are presented. It is shown that the reinforcement shape can have a significant effect on the effective properties of the composite material. The second type of material studied is a material which has a granular and/or fiberous microstructure. The microstructures of these types of materials are often modelled by a latticelike arrangement of load carrying members. Within the framework of the finite element method, a general technique is presented to compute the effective properties of materials which are modelled by these latticeli.ke microstructures. The equivalency between the continuum and microstructural stiffness matrices is utilized to produce an over-determined system of equations which is solved using the Moore-Penrose generalized inverse procedure. Although the resulting solution is not exact, it is unique in the least squares sense. The effective properties are reported in the form of a continuum constitutive matrix which approximates the behavior of the latticelike microstructure. Several specific examples are given to demonstrate the effectiveness and accuracy of the method.

25. Effect of inclusion aspect ratio on the macroscopic composite response for a volume fraction of . 10. 48 x Figure  Introduction A wide range of metal matrix composites are currently being developed for numerous structural applications. These advanced composite materials may contain either continuous fiber reinforcement or discontinuous particle or whisker reinforcements. These composites offer the advantages of improved strength, stiffness and dimensional stability as compared to the unreinforced metals. Their primary limitation in many applications is their relatively low ductility. As a result, an improved understanding of the plastic deformation mechanisms is essential to their successful implementation in structural applications. In this paper, the elastic-plastic response of metal matrix composites with discontinuous reinforcements of varying geometry is considered.
The elastic deformation characteristics of composite materials have been the subject of numerous research efforts in the past [ 1,2,3,4 ]. Several of these studies have provided methcxis which can be utilized to predict the effective elastic properties of composite materials. These studies have also provided a great deal of insight into how the details of the microstructure effect the elastic deformation characteristics of the composite material. The elastic-plastic deformation behavior of composite materials, however, is not well understood. Since the low ductility exhibited by particle reinforced metal matrix composites is their primary drawback, an improved understanding of their elastic-plastic response is required. Since analytical formulations become intractable when applied to this type of elastic-plastic analysis, the finite element methcxi has been utilized extensively in several recent studies which examine the elastic-plastic deformation characteristics of metal matrix composites. Several of these papers utilize the "unit cell" formulation to predict the elastic-plastic stress-strain response of the composite [5,6,7]. The unit cell typically consists of a single particle surrounded by a finite sized matrix shell. The geometry and boundary conditions for the unit cell problem arc deduced by assuming certain perfectly periodic and uniform distributions. The predicted macroscopic stressstra.in response of the composite is based upon the response of this unit cell. Although a considerable amount of useful information has been obtained with this approach, it suffers from the drawback that the macroscopic stress-strain response is very strongly dependent on the assumed boundary conditions. A similar method, known as the method of cells, has been utilized by Aboudi[8] to compute the effective properties of inelastic composites.
Another recent study [9] predicted the elastic-plastic response of a metal matrix composite reinforced with spherical inclusions by applying the dilute, self-consistent and generalized self-consistent (three-phase) approximations. These models, originally developed to estimate the composite effective elastic properties, were extended to the plastic range by an incremental application of the volume averaging scheme.
The purpose of the present work is to systematically examine the evolution of the local state variables in the matrix of the composite, to determine the effect of reinforcement shape on the evolution of these variables, and to determine the effect that the reinforcement shape has on the macroscopic stress-strain response of the composite.
The following two sections of this paper describe, respectively, the algorithm utilized to compute the macroscopic stress-strain response, and the boundary value problem selected to model the behavior of the composite. Typical results of this work are presented in sections three, four, and five. Section three consists of a systematic examination of the evolution of local matrix state variables for several reinforcement shapes. Section four consists of an examination of the effect of matrix yielding on the distribution of axial stress in the inclusion, for each of the reinforcement geometries considered. Section five 3 consists of an examination of the effect of reinforcement shape and volume fraction on the predicted macroscopic stress-strain response of the composite.

Macroscopic Composite Response Formulation
The macroscopic constitutive behavior of a composite subjected to a uniform far-field loading is given by the relation between the average stresses and strains within each phase.
Expressions for the macroscopic tangent moduli for a two-phase composite under uniaxial loading are developed. The analysis extends the work of Willis [ 10] for the effective elastic properties to give an incremental elastic-plastic formulation. The incremental macroscopic constitutive behavior of the composite can be represented by (1) where o-and tare the macroscopic stress and strain increments which are related through L -, which is a fourth-order tensor of effective composite tangent mcxluli, and L :E denotes the product Li.ikJEkl. Equilibrium requires that the macroscopic stress increment equals the volume average of the local stress fields (2) where O'(x) is the distribution of local stress increments and x is the position vector.
Similarly, the macroscopic strain increment equals the volume average of the local strain increments, t(x) 4 (3) where .V represents the total volume of the composite.
In a similar manner, the volume average inclusion and matrix stress and strain increments (O'i• om•~. and Eui. respectively) may be defined by (7) where Vi and V mare the total inclusion and matrix volumes, respectively. For the elastic inclusion, the relation between the average stress increment and the average strain increment is given by (8) where L 1 is the inclusion elastic mcxiuli tensor. The incremental stress-strain relation for the elastic-plastic matrix is given, approximately. by the relation a =L (~):t • • • • (9) where L.(£!) is the matrix tangent mcxiuli tensor and £!=[(2/3)€f!ft~ is the average matrix effective plastic strain.
For a two-phase composite, since Vi+Vm=V. equations (2-7) may be rewritten as 5 (10) (11) where vi = V fV and vm = V JV. Substituting equations (8) and (9) into (10) yields (12) Rewriting equation (11) as (13) and substituting into equation (12) yields (14) If the inclusion average strain increment can be related to the macroscopic stain increment through a relation of the form (15) where the fourth-order tensor A is the incremental strain concentration tensor, then equation (14)  . (16) Comparing equation (1) with equation (16) gives an expression for the the macroscopic tangent moduli tensor (17) Note that except for the approximation involving Lm(e:), this expression is exact.
Hence, estimates of the incremental strain concentration tensor, A, gives an estimate of the effective tangent moduli. In general, the macroscopic behavior is influenced by the 6 reinforcement volume fraction, shape, spatial distribution and orientation. The effect of reinforcement shape, distribution and orientation are all contained in the incremental strain concc~tration tensor.
The formulation presented above is now specialized to the case of uniaxial tension.
The elastic-plastic matrix material is assumed to obey 1 2 flow theory. If it also assumed that the composite obeys 1 2 flow theory, the composite uniaxial response may be used to determine a generalized, multiaxial response. The composite uniaxial stress-strain curve is determi n~ by applying a far-field loading (or straining) which is consistent with a uniform uniaxial stress. For the case of axisymmetric stressing (a'; 2 =a; 3 , ah=a~3, ~2=d;1 3 , and a ii=() for i;~j) the inclusion and matrix constitutive behavior can be written in matrix form as (18) and (19) where C~ and ~i are inclusion elastic moduli and matrix tangent moduli which are most conveniently expressed in terms of compliance components (where [C]=[Ml 1 ) as ,.. 2v,.. 1 ,.. ( v,.. 1 ,.. 1-v,.. 1 ,.. h E,.. 2h E,.. 2h (21) where p 1 , Bui. v, and vm are the inclusion and matrix Young's mcxiuli and Poisson's ratios, respectively and h(E~) is the matrix hardening parameter (=dcr/d£P in uniaxial tension).
For the case of uniaxial composite loading (<J7 1 =cr-, 02 2 =aj 3 =0), we can define the incremental strain concentration factors, A 1 =t; 1 1t7 1 , and A~t7 1 , a tangent mcxiulus, equation (17) gives (22) and (23) With appropriate estimates of the incremental strain concentration factors, A 1 and A 2 , equations (22) and (23) give an estimate of the effective composite tangent mcxiulus. The incremental uniaxial stress-strain relation, 0'7 1 =E-e7 1 , may then be integrated to yield the uniaxial composite response. In this study, the incremental strain concentration factors are estimated by neglecting interactions between neighboring inclusions. This approximation leads to the boundary value problem described below.

Boundary Value Problem
The micromechanical mcxiel utilized in this work is based on the dilute approximation. This methcxi assumes that, for low to moderate inclusion volume fractions, the interaction between stress inhomogeneities created by adjacent particles is 8 negligible. Thus, the dilute model assumes that the local state variables in an inclusion and its surrounding matrix shell can be obtained by solving the probl~m of a finite size inclusion embedded in an infinite matrix material. The boundary conditions utilized in this incremental boundary value problem are the known far-field boundary conditions imposed on the composite material; in this case, boundary conditions consistent with an applied uniaxial far-field stress. Perfect bonding at the inclusion/matrix interface is assumed. The solution of this incremental boundary value problem provides the details of the evolution of the local matrix and inclusion state variables, as well as the strain concentration factor tensor and volume average state variables required in the effective properties formula.
The material model chosen for the infinite matrix material is an initially isotropic elastic-plastic material. The infinite matrix material is assumed to obey the J 2 flow theory of plasticity with isotropic hardening with a uniaxial power law hardening description (24) where cr 0 and E 0 are the matrix yield stress and yield strain, respectively, CEm=crjE 0 ) and n is the strain-hardening exponent (n~l). The value of the strain-hardening exponent utilized in the present study is 10. The inclusion is modeled as an isotropic elastic material. All reported stress and strain values are normalized with respect to the matrix uniaxial yield stress and strain, respectively. As discussed by Tvergaard [7], and Christman et al. [5], for a typical 2124 Al-SiC whisker composite the assumption of rigid fibers does not result in any qualitative difference in the predicted behavior of the composite. Thus, since most discontinuously reinforced composite materials utilize a 9 reinforcing material that is much stiffer than the matrix material, it is appropriate to assume that the inclusions are nearly rigid. To mcxlel the case of rigid inclusions, the inclusipn mcxluli are taken to be three orders of magnitude larger than the corresponding matrix moduli.
The incremental boundary value problem described in the previous paragraphs was solved by the finite element methcxl. Due to the axial symmetry of the inclusions and the assumption of axisymmetric deformations, two-dimensional axisymmetric finite element mcxlels were utilized to mcxlel the boundary value problems under consideration. Also, due to mid-plane symmetry, it is only neccessary to mcxlel the upper half of the domain.
A finite element ccxle written for this project was utilized to solve the incremental boundary value problems under consideration. A four-ncxled quadrilateral isoparametric element belonging to the serendipity family is utilized in this ccxle [l 1,12]. The element as implemented utilizes a selective reduced integration procedure which avoids the problem of mesh-locking, and thus makes the element suitable for the analysis of elasticplastic problems [13]. The elastic-plastic constitutive equations are integrated by a forward Euler methcxl in conjunction with a consistency condition correction algorithm and a subincrementation strategy [14,15]. The incremental boundary value problem is solved using a displacement incrementation without partitioning continuation algorithm in conjunction with the piecewise linear self-correcting methcxl [16]. Since the strain-to-failure of typical metal matrix composites is small, a small strain formulation is utiliz.ed. Further details of the finite element formulation used are presented in appendix C.

Evolution of Local Matrix State Variables
It can be shown through volume averaging arguments that the macroscopic stiffness 10 of a composite is largely determined by the stress level in the inclusion. Since the stress level in the inclusions is largely determined by the matrix/inclusion load transfer process it is very important to gain insight into the details of this process. In this section the effect of reinforcement shape on the evolution of local matrix state variables, and thus the load transfer characteristics of the matrix/inclusion system, is examined. Five reinforcement geometries, one spherical and four ellipsoidal, are considered. Each of the ellipsoidal geometries considered has an axis of axial symmetry aligned with the direction of the applied uniaxial load(see Figure 1). The inclusion geometry is identified by its aspect ratio, i.e. the ratio of its axis length in the direction of the applied uniaxial load to the length of a perpendicular axis. Thus, an ellipsoidal inclusion with an aspect ratio greater than 1 is a fiber-like inclusion with its major axis aligned with the direction of the applied uniaxial load. Whereas an ellipsoidal inclusion with an aspect ratio less than 1 is a platelet-like inclusion with its minor axis aligned with the direction of the applied uniaxial load. The five inclusion geometries examined in this work have aspect ratios of 1 (sphere}, 2, 5, 1/2 and 1/5. The contours displayed at the smallest of the four values of far-field stress describe the load transfer characteristics of the inclusion geometry for the elastic case, as this value of far-field stress is reached immediately prior to any local matrix yielding. Thus, the smallest value of far-field stress at which the stress contours are presented is the value of far-_ field stress which must be surpassed to initiate local yielding in the matrix.
Effective plastic strain contours are not shown for the smallest value of far-field stress because no local yielding is present at that value of applied load. Contours shown for the three larger values of far-field stress are presented to illustrate the growth of the plastic zones and the effect these inhomogeneities have on the load transfer characteristics of the inclusion geometries. Contours are not displayed for far-field normalized stress values greater than 1 because at that value of applied load the entire matrix yields and no significant changes in load transfer behavior subsequently occur.
The evolution of matrix state variables for the case of a spherical inclusion is presented in Figures 2-4. Examination of the Mises stress contours in Figure 2 reveals that the majority of the load transferred to the inclusion is transmitted at or near the pole of the inclusion. It is interesting to note that the maximum Mises stress does not occur at the pole but rather at a polar angle of about 40 degrees. This is due to the large hydrostatic stresses which develop at the pole of the spherical inclusion. This result is in agreement with previous experimental [l 7, 18], analytical[l9] and nurnerical [20,9] results for the case of spherical inclusions. As illustrated in Figure 4 local yielding initially occurs at a polar angle of about 40 degrees. Early local yielding is also observed directly above the pole of the inclusion; although in the initial stages of yielding no plastic deformation is observed at the pole. This is also due to the large hydrostatic stresses which develop at the pole.
By comparing  it is seen that decreasing the aspect ratio from 1 to 1/2 has a detrimental effect on the efficiency of the load transfer process. This is revealed by the increase in the value of applied load required to initiate yielding, by the lower values of Mises and hydrostatic stress, and by the significant decrease in the magni~de of plastic strains which occur inspite of the lower values of hydrostatic stress.
Examination of Figure 7 indicates that the location of initial yielding occurs at a polar angle of approximately 70 degrees, and that early local yielding is again observed a shon distance above the pole. As in the case of the spherical inclusion, the absence of early local yielding at the poles is due to the large hydrostatic stresses present at that location.
Examination of Figures 8-10 indicates that a funher decrease in the aspect ratio of the inclusion to a value of 1/5 does not result in a noticeable change in the efficiency of the load transfer process. There is no change in the magnitude of the applied stress level required to initiate local yielding. Also, a comparison of the state variable contours for the 1/2 and 1/5 ellipsoids does not provide any data that makes it possible to conclude which platelet geometry provides the more effective load tranfer process. Although, it is interesting to note that the load transfer characteristics of the two ellipsoidal inclusions are quite different, e.g. the 1/5 ellipsoid has smaller stress and plastic strain levels at its pole but higher stress and plastic strain values along its lateral boundary.  indicate that increasing the aspect ratio of the inclusion to a value of 2 dramatically increases the efficiency of the load transfer process. This is clearly indicated by the significantly lower value of applied stress required to initiate local yielding, and by the large values of Mises stress that occur at that low applied stress level. An examination of  reveals that the majority of the load transferred to the inclusion is transmitted at or near the poles of the inclusion. This is indicated by the high levels of Mises and hydrostatic stress which exist near the poles, and by the fact that the initial yielding in the matrix occurs at a polar angle of about 20 degress. As observed in the previously considered geometries, large hydrostatic stresses located at the poles cause early yielding to occur at a distance above the poles rather than at the poles.
Althou_gh, the hydrostatic stresses located at the poles of the 2/1 ellipsoid are higher and more localized than observed in the spherical, and platelet-like inclusions.  indicate that further increasing the aspect ratio of the inclusion to a value of 5 dramatically increases the efficiency of the load tranfer process. This is again clearly indicated by the significantly lower value of applied stress required to initiate local yielding, and by the large values of Mises stress that occur at that low applied stress level. An examination  reveals that the majority of the load transferred to the inclusion is transmitted at or near the poles of the inclusion. This is indicated by the high levels of Mises and hydrostatic stress which exist near the poles, and by the fact that the initial yielding in the matrix occurs at a polar angle that is close to 0 degress. As observed in previously considered geometries large hydrostatic stresses exist at the poles of the 5/1 inclusion. However, in this case the region of large hydrostatic stresses is very small which allows initial yielding to occur very close to the poles.
An examination of the Mises stress contours presented for the different inclusion geometries reveals that, as the aspect ratio of the inclusion is increased the far-field stress value required to initiate local yielding decreases. Thus, as the inclusion aspect ratio increases the stress inhomogeneity becomes more severe. This indicates that as the inclusion becomes more fiber-like the load transfer process becomes more efficient, and the inclusion becomes a more effective load carrying component of the composite.
Comparing Figures 2,5 and 8 it is seen that as the aspect ratio is decreased from a value of 1 the decrease in load carrying effectiveness of the inclusion is modest. Whereas by comparing Figures 2, 11 and 14 it is seen that increasing the aspect ratio from a value of 14 1 results in a significant increase in the load carrying effectiveness of the inclusion.
These results indicate that the value of the aspect ratio of a fiber-like inclusion can have a significant effect on its effectiveness as a load carrying component of the composite, whereas the aspect ratio of a platelet-like inclusion has a relatively minor effect on its effectiveness as a load carrier. to the inclusion by shear along its lateral boundary becomes a significant percentage of the total load carried by the inclusion. Thus, as the aspect ratio of the discontinuous inclusion is increased its load transfer characteristics approach those of a continuous fiber.
It is also interesting to note that fiber-like inclusions( 2/1 and 5/1) develop significantly larger hydrostatic stresses than the platelet-like inclusions. These large hydrostatic stresses that develop at the poles of the fiber-like inclusions make interfacial decohesion and/or void nucleation likely. Among the likely consequences of interfacial decohesion and/or void formation are a decrease in load transfer efficiency and a tendency for the inclusion to behave more like a fiber, i.e. the majority of the load would be transferred to the inclusion by shear along its lateral boundary.

Inclusion Axial Str~ Distribution
Due to the work of Eshelby[21], it is well known that prior to any matrix yielding the stress distribution inside the inclusion is uniform. In this section contours of axial stress arc presented for the 1/1, 2/1, and 5/1 inclusion geometries; contours are not presented for the platelet-like geometries because only a slight deviation from the initial uniform stress state was observed. The contours are displayed for a far-field normalized stress value equal to 1. These contour plots ilustrate the effect of matrix yielding on the stress distribution inside the inclusion. An examination of Figure 17 reveals that the axial stress distribution inside the inclusions is no longer uniform. That figure also reveals that the deviation from a uniform stress state is more pronounced for the fiber-like inclusions than for the platelet-like inclusions. It is also interesting to note that the distribution of axial stress in the 5/1 inclusion is similar to the distribution expected for a continuous fiber, this supports the assertion that as the aspect ratio is increased the load transfer process more closely approximates that observed for continuous fibers.
It is important to note that the levels of axial stress observed in the fiber-like inclusions are significantly larger than those observed in the platelet-like inclusions; this observation is in agreement with the results of the previous section that indicated that the load transfer process is more efficient for the fiber-like inclusions. To further investigate the efficiency and characteristics of the load transfer process, the values of the volume average inclusion stresses are shown in Table 1. The results shown in Table 1 were computed for a far-field normalized stress value equal to 1. An examination of Table 1 indicates that the fiber-like inclusions strengthen the composite by carrying a considerable percentage of the applied axial load, and that decreasing the inclusion aspect ratio decreases the percentage of axial load carried by the inclusion. An examination of Table  1 also indicates that the levels of transverse normal stress which occur inside the inclusion increase with decreasing aspect ratio. Since large values of transverse normal stress cause the vo~ume averaging equations to predict large levels of matrix hydrostatic stress, the platelet geometries result in higher levels of matrix triaxiality. Thus, it appears that for fiber-like inclusions the dominant strengthening mechanism is direct load transfer, whereas for platelet-like geometries the dominant strengthening mechanism is constraint of the matrix.

Effective Composite Response
In this section estimates of the effective composite response, in the form of uniaxial stress-strain curves, is presented. These curves illustrate the effect of the inclusion aspect ratio and volume fraction on the macroscopic properties of the composite. These anticipated trends, which are in agreement with previous experimental [5,6] and numerical [5,6,7 ,9] results, are observed for each of the inclusion shapes considered.
An examination of  also indicates that neither the inclusion volume fraction, or the inclusion geometry, have a sigi;iificant effect on the amount of strainhardening predicted for the composite. This indicates that the efficiency of the load transfer process has virtually no effect on the amount of composite strain-hardening.
Thus, a strengthening mechanism other than direct load transfer must be responsible for the increase in composite strain-hardening observed in experimental studies [5,6].
Experimental results discussed by Christman et al. [6] indicate that at an inclusion volume fraction of 13.2 vol.% the strain-hardening exponent of a typical whisker reinforced composite decreased by approximately a factor of three, as compared to the unreinforced matrix material; whereas for a typical paniculate reinforced composite the decrease in strain hardening exponent was negligible even at inclusion volume fractions as large as 20 vol.%. Christman et al. [5,6] concluded, based on the results of a unit cell finite element analysis, that the increase in composite strain-hardening is due to the existence of large hydrostatic stresses caused by the geometric constraint of neighboring inclusions.
Since the dilute approximation utilized in the present work does not account for interaction between neighboring inclusions, the large hydrostatic stresses identified by Therefore, it is expected that neglecting the geometric constraint of neighboring inclusions, and the induced hydrostatic stresses, is a more severe approximation for the fiber-like inclusion geometries. This fact is also supported by comparing the negligible increase in composite strain-hardening observed in Figure 22 to the significant increase observed experimentally for the whisker reinforced composite, and also by comparing the negligi'ble increase in composite strain-hardening observed in Figure 18 to the negligible increase observed experimentally for the particulate reinforced composite. That conclusion leads to the following related conclusions: (1) the accuracy of the dilute model predictions will depend upon the inclusion aspect ratio, as well as the inclusion volume fraction, of the specific composite under consideration, and (2) that the dilute model is better suited to the analysis of particulate reinforced composites than whisker reinforced composites. Thus, for a given volume fraction, results computed by the method utilized in this work will tend to be more accurate for particulate inclusion geometries than for fiber-like inclusion geometries. It is also interesting to note that the 1/5 inclusion geometry results in the most significant increase in composite strain-hardening. This is explained by noting that nearly the entire boundary of the 5/1 inclusion is subjected to significant levels of hydrostatic stress which result in large levels of transverse normal stress( crb & cr~3) inside the inclusion. Due to these large inclusion transverse normal stresses the volume averaging equations predict large levels of matrix hydrostatic stresses which cause the observed increase in composite strain-hardening.

An examination of Figures 23-26 reveals that as the inclusion volume fraction is
reduced the curves converge to a single curve, i.e. the uniaxial stress-strain curve of the matrix material. It is also observed that the same trends are observed in each of these figures; although the trends are much more pronounced for the higher inclusion volume fractions. An examination of Figure 26 reveals that the composite elastic moduli computed for the 1/2 and 1/5 inclusion geometries are virtually identical, and that the composite elastic moduli computed for the spherical(l/1) inclusion geometry is only slightly larger. Whereas, the computed value of the composite elastic moduli increases dramatically as the inclusion aspect ratio is increased to the 2/1 and 5/1 values considered.
It is iII?portant to note that the inclusion geometries which result in the largest values of composite elastic modulus were identified in the previous section as the most efficient load carrying geometries. This result leads to the following conclusions: (1) in the case of fiber-like inclusions direct load transfer between the matrix and the inclusion is an important strengthening mechanism which significantly effects the elastic modulus of the composite, (2) in platelet-like inclusions constraint of the matrix is an important strengthening mechanism, (3) since the aspect ratio of the inclusion significantly effects the efficiency and characteristics of the load transfer process, it also significantly effects the ability of the inclusion to increase the elastic modulus of the composite, and (4) (25) where for the 1/5, 1/2, 1/1, 2/1 and 5/1 inclusion aspect ratios the value of the constant C is equal to 1.7,1.6,2.0,3.3 and 8.8, respectively. The composite elastic modulus computed for the case of the 1/5 inclusion is larger than for the 1/2 inclusion, in spite of the fact that the 1/2 inclusion results in a more efficient direct load transfer process. This is due to the large levels of transverse normal stress which occur inside the 1/5 inclusion. inclusions. This indicates that the yield strength of the composite is also strongly effected by the efficiency and characteristics of the load transfer process, and that the fiber-like inclusipns provide greater increases in yield strength than the platelet-like geometries.

Concluding Remarks
In this study, the dilute approximation was utilized to examine the effect of reinforcement geometry on the macroscopic response of discontinuously reinforced composites. A systematic examination of the evolution of local matrix state variables revealed that the aspect ratio of the inclusion has a significant effect on the efficiency and characteristics of the load transfer process. It was determined that fiber-like inclusion geometries resulted in a significantly more efficient direct load transfer process than the platelet-like geometries, but that the platelet-like geometries provide more efficient constraint of the matrix. Those conclusions were further witnessed by the higher levels of axial stress present inside the fiber-like inclusions, and by the higher levels of transverse normal stress observed inside the platelet-like inclusion geometries.
An examination of the computed uniaxial stress-strain response of the composite revealed that the aspect ratio of the reinforcement phase also has a significant effect on the uniaxial response of the composite. It was shown that the fiber-like inclusion geometries resulted in significantly larger values of composite yield strength and elastic modulus, as compared to the platelet-like inclusion geometries. That result led to the conclusion that, for dilute concentrations of reinforcement, the fiber-like inclusion geometries are considerably more effective as reinforcement than the platelet-like geometries. It was also determined that the inclusion aspect ratio and volume fraction had only a very small effect on the amount of composite strain-hardening predicted. This led 21 to the conclusion that the amount of composite strain-hardening is nearly independent of the efficiency and characteristics of the load transfer process.     [7] Tvergaard, V., "Analysis of Tensile Properties for a Whisker-Reinforced Metal-Matrix Composite", Acta Metall., Vol. 38, pp. 185-194, 1990.
[9] Taggart, D. G., "Modeling of the Elastic-Plastic Response of Metal Matrix Composites with Spherical Reinforcements," in preparation.

Introduction
The mechanical behavior of materials composed of granular and/or fiberous microstructures, is inherently involved with the transmission of loadings along discrete paths within the material. This behavior is fundamentally different from that predicted by classical theories of continuum mechanics. It has been observed, see for example Oda [l], that the particular distribution of internal load transfer depends strongly on the material's microstructure. A considerable amount of contemporary research has been 52 conducted in order to understand and explain this complex microstructural behavior. Bun and Dougill [2] proposed a simple planar pin-jointed truss network to simulate the stressstrain . behavior of heterogeneous materials. Based on the fact that granular materials transmit loads only through contact mechanisms, they have been modeled with network theories, e.g. Trollope and Burman [3], Bagster and Kirk [4], Bideau et.al. [5], Thornton and Barnes [6] , and Sadd, Qiu and Boardman [7]. The connection between various lattice gridworks and an equivalent micropolar continuum has been investigated by Banks and Sokolowski [8], Bazant and Christensen [9] and Sun and Yang (10] . The relationship between the micromechanics and the overall macro material behavior is a very imponant and fundamental issue in current materials research. For example, this relationship is needed to understand how localized microstrutural failure will lead to global failure of the entire body, and to predict the effective global properties of a material knowing its i i microstructure.
In addition to this type of research, there has also been considerable interest by the structural mechanics community in the continuum modeling of large repetitive lattice structures. The relationship between a continuum and a gridwork of discrete elements was examined in some early work which involved analyzing a continuum by replacing it with an equivalent elastic grid work [ 11,12]. More recently the inverse problem of continuum modelling of repetitive latticelike structures has received considerable attention (13)(14)(15)(16)(17), see the excellent review by Noor (18]. The ability to replace a large repetitive gridwork with an equivalent continuum provides a means to simplify cenain calculations for large space structures. Based on finite element procedures, several methods of computing the properties of an equivalent continuum have been developed. The majority 11 of these methods were developed for beamlike and platelike lattices, where a reduction 53 in dimensionality was an integral part of the equivalency problem. Current methods of establishing the equivalency between structural and continuum models have been primarily based on stiffness or energy methods [18]. These techniques have proposed continuum finite elements with stiffness (constitutive) and strain and/or kinetic energy properties that approximate to some degree of accuracy, the structural system. The present state of research indicates that while satisfactory equivalent continuum models for linear behavior of beamlike and platelike lattices can be developed, methods for general continua still need further refinement and study.
The current article addresses this general problem of developing continuum models for materials with latticelike microstructure. The main objective is to construct a somewhat general technique to compute the constitutive matrix of an equivalent continuum given the microstructural properties. The study is limited to two-dimensional plane problems; however, the technique can easily be extended to three dimensions. After envoking equivalency between the continuum and microstructural stiffness matrices, the resulting over-determined system is solved using the Moore-Penrose generalized inverse procedure. This provides a unique solution in the least-squares sense. Several specific examples are then presented to demonstrate the effectiveness and accuracy of the method.

Formulation of Equivalent Continuum Properties
In the method presented here it is assumed that the stiffness matrix of a repeating element of the latticelike structure is known. This stiffness matrix may be computed by utilizing the direct stiffness method to assemble each of the discrete clement stiffness matrices in the repeating element into a global stiffness matrix. For the type of problems under consideration here the repeating element of the lanicclike structure and the element of the equivalent continuum are to have the same degrees of freedom. Therefore, the two elements will exhibit identical behavior if their stiffness matrices are identical. Thus, the stiffne~s matrix of the element of the equivalent continuum is required to satisfy the equation ( 1) where [Klt and [K]c are, respectively, the stiffness maaices of the repeating cell and the equivalent continuum.
Using standard procedures from finite element analysis, the continuum stiffness maaix can be written as (2) where  (2), the equivalency equation (I) can be rewritten as (3) where lo· 3u· and IJul are, respectively, the thickness, the weighting factor, and the determinant of the Jacobian, all evaluated at the integration point r 1 , sJ. The values of the weighting factors and the location of the integration points (which are specified in natural coordinates r and s) are determined by the order of Gaussian quadrature chosen [ 19]. The values of the :rain-displacement maaix and the Jacobian maaix are determined by the continuum finite element selected to model the equivalent continuum. Thus, once the interpolation scheme has been chosen, the value of [D] is the only unknown in equation (3). Therefore, equation (3) can be used to solve for the constitutive matrix [D] of the equiva_lent continuum.
At this point it should be noted that equation (3) is only valid for a homogeneous equivalent continuum, i.e. the value of [DJ is assumed to be the same at each sampling point in the continuum element. It should also be noted that it has been implicitly assumed that a strain energy function exists for the equivalent continuum, i.e. that [DJ = [DJT. Therefore, it is only necessary to solve for the upper diagonal terms of the constitutive matrix.
As mentioned previously, the solution represented by equation (8) is not an exact solution.
However, provided that [A) has linearly independent columns, it is a unique solution in the least-squares sense.

Examples and Discussion
The method presented m the previous section was utilized, with a bi-linear This alternate method involves relating the force and deformation characteristics of a small segment of the grid to those of a small segment of the equivalent continuum as each is subjected to a homogeneous deformation [14].  Figure 2 and then computing the percent difference in the stored strain energy.
The displacement vector specified by the selected nonhomogeneous deformation state contains only non-zero terms. Thus, each term in the stiffness matrices will be involved in the strain energy computation. This avoids the possibility of inadvenently accounting for the accurate terms and neglecting the inaccurate terms, or vice versa. It should be noted that the error measure listed in Table 1 is not a function of the magnitude of the chosen displacement vector; however, it is a function of the "direction" of the displacement vector. Thus, since the error in stored strain energy is a function of the rotations. That type of orthotropic material behavior seems to be more reasonable than the isotropic behavior predicted by the alternate method.
The example shown in Figure Figure 3 is the same cell which is shown in Figure la; thus, the equivalent continuum properties computed for case 1-2 in Table 1 can be utilized in each of the examples shown in Figure 3. Table 2 shows the percent difference in stored strain energy between the patches of repeating cells and the continuum element models, when each is subjected to the nonhomogeneous deformation state shown in Figure 2.
The lattice shown in Figure 3a consists of four repeating cells of the type examined in Figure la. The equivalent continuum model utilized to approximate the behavior of 61 the lattice shown in Figure 3a consists of four continuum elements. As listed for case 2-1 in Table 2 the percent difference in stored strain energy computed for the lattice and the patch ~f continuum elements is 9.0 percent. This seems to indicate that the equivalent continuum model of an entire latticelike structure would be nearly as accurate as the equivalent continuum model of a single repeating cell.
The lattice shown in Figure 3b also consists of four repeating cells of the type shown in Figure 1 a. The equivalent continuum model utilized to approximate the behavior of the lattice shown in Figure 3b consists of one continuum element. As listed for case 2-2 in Table 2 the percent difference in stored strain energy for this case is 9.2 percent. The results of this example seem to indicate that an entire latticelike structure could be modelled by a relatively small number of continuum elements, with a degree of accuracy that is largely determined by the accuracy of the equivalent continuum properties utilized.
Thus, this example illustrates the benefit of the formulation presented in the previous section, as well as the power of equivalent continuum modelling.

Concluding Remarks
The results presented in Table 1 indicate that the method under consideration yields equivalent continuum properties which are a good approximation to a repeating cell, e.g.
the mean error computed for the repeating cells shown in Figure 1 is 3%. The results presented for the repeating cells in Figures la and 1 b seem to indicate that for cases which involve nonhomogeneous deformation states the method presented here is significantly more accurate than the alternate method examined.
The results presented in Table 2 seem to indicate that by utilizing the present method to compute the equivalent continuum properties of a repeating cell an entire latticelike 62 structure could be modelled, with a reasonable degree of accuracy, by using a relatively small number of continuum elements.  [14]. 2Equivalent continuum properties computed by method described in this work. 3 Equivalent continuum properties were normalized by dividing them by the modulus of elasticity utilized for the discrete members in the repeating cell. Table 2. Computed difference in stored strain energy between a patch of repeating cells and the equivalent continuum models shown in Figure 3.

Limitations Imposed on the Range of Applicability
The micromechanical model utilized for the work presented in manuscript 1 is an idealized representation of a discontinuously reinforced composite. That idealized representation is based upon the assumption that the distance between neighbooring inclusions is sufficient to insure that the interaction between stress inhomogeneities created by adjacent inclusions is negligible. Since this assumption is only satisfied for low to moderate inclusion volume fractions, the results presented in manuscript 1 are only applicable to composites with a low to moderate inclusion volume fraction.

71
The purpose of the research presented in manuscript 1 is to predict the effective properties of a composite which consists of stiff inclusions in a ductile, strain-hardening elastic~plastic matrix, and to determine the effect of matrix yielding on these effective properties. It is imponant to note that no attempt is made in this work to determine the strength of the composite, or its mode of failure. With these considerations in mind, the following additional assumptions are made in the development of the numerical model: (1) perfect bonding at the inclusion/matrix interface is assumed, (2) it is assumed that the modulus of elasticity and yield strength of the inclusion are three orders of magnitude larger than the corresponding matrix values. Thus, effects such as interfacial decohesion, fiber elasticity, and fiber breakage are not accounted for. Therefore, the results presented in manuscript 1 are applicable to composites consisting of a low to moderate inclusion volume fraction of stiff, high strength inclusions surrounded by (and perfectly bonded to) a ductile, strain-hardening elastic-plastic matrix. Also, since the nucleation and growth of voids is not considered (and a small strain formulation is utilized), the results presented in manuscript 1 are not valid for large strains.
The method of computing effective propenies presented in manuscript 2 assumes that the degrees of freedom of the repeating cell are identical to those of the equivalent continuum element. That assumption places restrictions on the type of discrete elements that can be utilized in the repeating cell. For example, if a conventional continuum is utilized to model the equivalent continuum, the discrete elements can not contain rotational degrees of freedom.

Recommendations for Future Research
The proposed research topics are designed to remove some of the restrictions which 72 currently apply to the research presented in manuscript 1. The propoposed topics are: (1). Utilize a strain-softening constitutive model with an appropriate yield criteria . and finite element formulation to study interfacial decohesion.
(2). Utilize a strain-softening constitutive model with an appropriate yield criteria and finite element formulation to study fiber breakage.
(3). Implement appropriate advanced constitutive models in the framework of a large strain finite element formulation to study modes of failure, and the effect of the composite characteristics on ductility.
The following research topic is proposed to allow the use of more general discrete elements in the effective property formulation presented in manuscript 2. The proposed topic is: (4). Utilize a micropolar continuum to model the equivalent continuum. This would allow for a more sophisticated micromechanical model to be constructed enabling concentrated moments to exist in the media.

APPENDIX B
Comprehensive Literature Review of the Research Area · The way a material behaves when it is under the influence of external loads is a ramification of the microstructure of the material. However, the scale of the structure to be analyzed is usually many orders of magnitude larger than the microstructural details of the material. For that reason, it is not possible to directly include the details of the material's microstructure in the model of the structure. Therefore, the material must be treated as a continuum and its behavior modelled by a continuum constitutive equation. which Jail by progressive damage is reviewed. In addition to those topics, which fall under the heading of micromechanics, an extensive review of the numerical methods commonly utilized to analyze micromechanics problems has also been completed. A review of the literature which penains to these numerical methods is contained in the following appendix.

Survey of Literature Pertaining to Reinforced Composites
Composite materials, which usually consist of stiff, brittle reinforcements embedded in a compliant matrix, constitute a large percentage of the advanced new materials which are currently being developed for many structural applications. Thus, the problem of computing effective properties for these materials is currently receiving alot of attention.
A common type of advanced composite material is the metal matrix composite. A wide range of metal matrix composites are currently being developed for numerous structural applications. These advanced composite materials may contain either continuous, or discontinuous, reinforcements. Metals reinforced with continuous brittle fibers are often known as fiber reinforced, or continuously reinforced, metal matrix composites. Metals reinforced with discontinuous brittle particles are often referred to as particle reinforced, or discontinuously reinforced. Typically, the reinforcing particles(or whiskers) have dimensions of a few micrometers. The work presented in this dissertation concerns the discontinuously reinforced metal matrix composites.
The elastic deformation characteristics of composite materials have been the subject of numerous research efforts in the past (Christensen 1979, Dewey 1947, Halpin 1984, Mura 1987, Walpole 1984, Willis 1977. Several of these studies have provided methods which can be utilized to predict the effective elastic properties of composite materials.
These _studies have also provided a great deal of insight into how the details of the microstructure effect the elastic deformation characteristics of the composite material.
However, the elastic-plastic deformation characteristics of composite materials is not well understood. Since the low ductility exhibited by panicle reinforced metal matrix composites is their primary drawback, their elastic-plastic deformation characteristics need to be more fully understood. Due to the fact that analytical formulations become intractable when applied to this type of elastic-plastic analysis, the finite element method has been utilized extensively in several recent studies which examine the elastic-plastic deformation characteristics of metal matrix composites. Several of these papers utilize the "unit cell" formulation to predict the elastic-plastic stress-strain response of the composite ). The unit cell typically consists of a single brittle panicle surrounded by a finite sized matrix shell. The geometry and boundary conditions for the unit cell are deduced by assuming certain perfectly periodic and uniform distributions. The predicted macroscopic stress-strain response of the composite is based upon the response of this unit cell. Although a considerable amount of useful information has been obtained with this approach, it suffers from the drawback that the macroscopic stress-strain response is very strongly dependent on the assumed boundary conditions.
Another recent study (Taggart, in preparation) predicted the elastic-plastic response of a metal matrix composite reinforced with spherical inclusions by utilizing the dilute, self-consistent and generalized self-consistent approximations. These models, originally developed to estimate the composite effective elastic properties, were extended to the plastic range by an incremental application of the volume averaging scheme. Materials reinforced with brittle spherical particles have also been the focus of other numerical (Wilner 1986), analytical (Gcxxiier 1933, and experimental Heikens, 1984, 1985) studies.

Survey of Literature Pertaining to Granular Materials
The second type of material studied in this dissertation is the type of material which contains a granular and/or fiberous microstructure. The mechanical behavior of materials composed of granular and/or fiberous microstructures, is inherently involved with the transmission of loadings along discrete paths within the material. This behavior is fundamentally different from that predicted by classical theories of continuum mechanics.
It has been observed, see for example Oda , that the particular distribution of internal load transfer depends strongly on the material ' s microstructure. A considerable amount of contemporary research has been conducted in order to understand and explain this complex microstructural behavior. One such project utilized a simple planar pin-jointed truss network to simulate the stress-strain behavior of heterogeneous materials . Based on the fact that granular materials transmit loads only through contact mechanisms, they have been modeled with network theories , Sadd, Qiu and Boardman in preparation, Thornton and Barnes 1986, Trollope and Burman 1980. The connection between various lattice gridworks and an equivalent micropolar continuum has been investigated by several researchers .
In addition to this type of research, there has also been considerable interest by the strUctural mechanics community in the continuum modeling of large repetitive lattice strUctures. Since the work presented in this dissertation involves computing effective propeqies for materials with microstructure that is often modelled as a latticelike network of elements, the aforementioned body of research regarding repetitive lattice structures is applicable to the micromechanics problem under consideration. The relationship between a continuum and a gridwork of discrete elements was examined in some early work which involved analyzing a continuum by replacing it with an equivalent elastic gridwork (Hrennikoff 1941. More recently the inverse problem of continuum modelling of repetitive latticelike structures has received considerable attention , see the excellent review by Noor . The ability to replace a large repetitive gridwork with an equivalent continuum provides a means to simplify certain calculations for large space structures. Based on finite element procedures, several methods of computing the properties of an equivalent continuum have been developed. The majority of these methods were developed for beamlike and platelike lattices, where a reduction in dimensionality was an integral part of the equivalency problem.
Current methods of establishing the equivalency between structural and continuum mcx:lels have been primarily based on stiffness or energy methods

survey of Literature Regarding Failure by Progr~ive Damage
Many common materials fail by progressive damage. In brittle materials the damage is oftep in the form of a system of densely distributed microcracks. Materials which exhibit microcracking include: concrete and rock, stiff clays, ice, sea ice, wood, particle board, paper, two-phase ceramic composites, fiber-reinforced polymers and fiberreinforced concrete, see the excellent review by Bazant (Bazant 1986). In ductile materials the damage is often characterized by localized bands of high strain, i.e. shear bands. This behavior is typical of structural metals and saturated clays (Leroy et al. 1987) In mechanical testing of materials which fail by progressive damage, the test specimen does not fail when the maximum stress is reached; instead a gradual decrease of stress with increasing strain is observed. This phenomenon is known as strainsoftening. In order to analyze strain-softening materials, within the framework of continuum mechanics, it is neccessary to utilize constitutive relations that exhibit strainsoftening behavior. It is well known that any elastic, or elastic-plastic material becomes unstable when the matrix of its tangent moduli is no longer positive definite. Since the tangent modulus of the declining branch of the stress-strain diagram is negative, this indicates that the test specimen must lose stability and fail as soon as the maximum stress is reached. However, the existence of strain-softening is an experimental fact. Bazant (Bazant 1976) examined this apparent inconsistency by formulating a stability analysis of a uniaxial test specimen. The conclusion drawn from the stability analysis is that strain-softening is not possible in a continuum, i.e. strain-softening is only possible in a heterogeneous material. In view of that conclusion, the problem may be stated as follows .
In order to formulate a phenomenological model of strain-softening, the material must be treated as a continuum. However, strain-softening is not possible in a continuum. This 79 apparent contradiction is resolved by assuming that a certain size exists such that the material clement is large enough to be approximately treated as a continuum, but small enough to prevent instability. By vinue of that assumption, it is seen that utilizing a strain-softening constitutive equation is not an unreasonable approach to modelling progressive damage.
A mathematical model of the phenomenon of progressive damage is obtained by combining a strain-softening constitutive equation with the relevant field equations of continuum mechanics. Since the resulting governing equations are quite complicated, the primary method of solution used has been the finite element method. Additionally, some analytical solutions have been obtained (Belytschko et al. 1986). As pointed out by several researchers (Bazant 1976, there are certain mathematical difficulties inherent to the finite element formulation of a strain-softening problem. Consequently, if this model is to be of any use, these difficulties must be handled properly. A discussion of the problems which arise in finite element analyses of strain-softening materials follows.
A common problem encountered in finite element analyses of strain-softening materials is that of spurious mesh size sensitivity (Bazant 1976, 1985a, 1985b, Bazant and Cedolin 1979. This problem received alot of attention after it was discovered that several widely used finite element codes with strain-softening models computed results that were a function of the mesh size chosen by the analyst. In the case of localized cracking, the cause of spurious mesh size sensitivity can be determined by examining the fracture characteristics of the material. Fracture energy(G,) is defined as the energy consumed per unit extension( and per unit thickness) of the crack surface. Fracture energy may be computed as (B-1) where "Ne is the width of the fracture process zone(the width of the finite element which lies along the crack band), _ and Wr is the total work of the tensile stress per unit volume, which is equivalent to the total area under the stress-strain curve. Since both Wr and G, are material constants, it is seen from equation (B-1) that w c is also a material constant.
The value of wc can be determined empirically. The value of wc can also be determined theoretically by utilizing a micromechanics analysis (Zubelewicz and Bazant 1987).
The cause of the spurious mesh size sensitivity is revealed by examining equation (Bl). If the same stress-strain curve is utilized, without regard to element size, the fracture energy will converge to zero as the size of the element is decreased. Since the fracture energy of the material is a constant, this demonstrates that utilizing the same stress-strain curve, regardless of element size, is not a rational approach.
The problem of spurious mesh size sensitivity can be solved in many ways. A few of the more common approaches, which have proven effective in the case of localized cracking, will be described. Instead of utilizing the strength criterion which yields incorrect results, the problem of spurio " mesh sensitivity can be alleviated by utilizing an energy criterion (Bazant 1985b) to determine crack band propagation. The energy criterion consists of modifying the stress-strain curve to maintain a constant value for the fracture energy. Thus, the stress-strain curve utilized depends on the size of the mesh.
Another effective approach involves embedding a shear band(or crack band), of a prescribed size, into a finite element (Belytschko, Fish andEngelmann 1988, Pietruszczak andMroz 1981). Hillerborg (Hillerborg, Modeer and Petersson 1976) solved the problem by utilizing a strain-softening, stress-displacement relation, instead of the stress-strain 81 curve.
Another difficulty which arises in finite element analyses of strain-softening materials, is that the solution can be dependent upon the orientation of the mesh. This problem arises because large strains can accumulate inside the strain-softening region without substantially affecting the strains in the material adjacent to this region. Since most conventional finite element formulations are only able to accurately describe such a localization of strain if the element boundaries are parallel to the shear band(or crack band), this results in a directional dependence. A few of the more common methcxls of dealing with this difficulty are described below.
In many analyses the direction of the crack band is known, or may be determined, in advance. In those cases the finite element mesh may be aligned with the crack band to help ensure an accurate description of the strain field. That approach avoids the problem of directional dependence but does not solve it. A more general approach involves conducting a localization analysis, at the element level, to determine the orientation of the shear band (Belytschko, Fish and Engelmann 1988, Leroy et al. 1987, Ortiz, Leroy and Needleman 1987. Once the orientation is determined, appropriate additional terms are appended to the interpolating polynomial of the element. These additional terms allow the element to accurately describe the strain field under consideration.
Another consideration that must be addressed is the stability of strain-softening finite element analyses. An instability, in the finite element model, is manifested as a nonpositive definite, incremental stiffness matrix. Therefore, stability requires that the mcremental stiffness matrix be positive definite. That criterion can be utilized to detect instabilities. For example, in the case of nonlocalized strain-softening, a nonpositive definite, incremental stiffness matrix may indicate a strain localization instability, i.e. the trUe solution may actually be localized strain-softening.
If t}ie aforementioned difficulties are correctly resolved, the finite element method can be utilized to compute correct solutions to problems involving strain-softening. The results obtained, however, are only useful if the constitutive equation is an accurate model of the real material. For that reason, the determination of accurate constitutive equations has been the subject of a great deal of research. Many attempts have been made to develop accurate strain-softening constitutive relations based upon a macroscopic, phenomenologic viewpoint. This work has produced various sophisticated strain-softening models, e.g. plastic, fracturing, plastic-fracturing, see the excellent review by Bazant (Bazant 1986). Since the macroscopic, phenomenologic approach has not been completely successful, many researchers are now attempting to derive strain-softening constitutive equations based upon physics and micromechanics analysis (Bazant and Oh 1985, Budiansky and O'Connell 1981, Frantziskonis and Desai 1987a, 1987b, Hoenig 1978, Horii and Nemat-Nasser 1983, Kachanov 1987, Montagut and Kachanov 1987, Simo and Ju 1987a, 1987b, Sumarac and Krajcinovic 1987.
Several finite element formulations have been developed to analyze strain-softening materials. A few of the more effective methods are described in the following paragraphs. The models to be discussed are shown in Figures B-1 and B-2.
Bazant and coworkers (Bazant 1976, 1985a, 1985b, Bazant and Cedolin 1979 analyzed the localized fracture of heterogeneous materials by utilizing an approach known as the blunt crack band model. In the analysis of localized fracture, a fracture process zone of finite size is assumed to exist( see Figure B-1 (a)). The fracture process zone is a region of material adjacent to a macrocrack that contains a system of densely distributed 83 microcracks. The blunt crack band model involves mcxielling the fracture process zone as a single band of strain-softening finite elements. The blunt crack band mcxiel has also been utilized successfully in cases involving nonlocalized cracking, i.e. cases where many neighbooring elements have entered the strain-softening region of the stress-strain curve.
It has been reported that this approach yields results, which are independent of mesh size, for localized and nonlocalized cracking(strain-softening) situations. The blunt crack band model, however, is less effective in the analysis of cracks that are not parallel to the sides of the finite elements. Pietruszczak and Mroz (Pietruszczak and Mroz 1981) employed an approach which involves embedding a shear band, of a prescribed size, in a finite element(see Figure B-1 (c)). The width of the shear band is not computed in this model, and thus, must be inferred from other sources. This method appears to have mitigated the problem of spurious mesh size sensitivity. However, in this approach the shape functions were not modified. Therefore, the elements ability to capture the strain localization will depend on their orientation. A more recent model developed by Belytschko and coworkers (Belytschko, Fish and Engelmann 1988) utilizes a finite element with embedded localization zones and modified shape functions. It appears that Belytschko's model alleviates the problem of directional dependence, as well as the problem of spurious mesh size sensitivity. The size of the shear band utilized in Belytschko's model, however, must be determined from other sources, e.g. from experimental measurements. The model by Bclytschko utilizes several features of. an earlier model developed by Ortiz and coworkers (Leroy ct al. 1987, Ortiz, Leroy andNeedleman 1987).
Bazant and coworkers (Bazant 1984, Bazant, Belytschko and Chang 1984, Bazant and Chang 1987) also developed a model based upon the idea of a nonlocal continuum. It 85 is known as the imbricate(overlapping) finite element model(see fig. B -2). This model can be utilized in cases where it is neccesary to use finite elements which are smaller than the strain-softening zone. For element sizes which are equal to, or greater than, the width of the strain-softening zone, the imbricate model is equivalent to the crack band model.

Computer Code Outline and Verification
This appendix contains a description of the computer code written for the thesis research and a discussion of the numerical methods that were implemented in the code.
The appendix is organized in the following manner. First the organization of the computer code is discussed and presented in a flowchart. Next, each major step in the flowchart is discussed. Finally the solutions to several verification problems are presented and compared to available solutions.

Computer Code Outline
The computer code written for this project is a nonlinear, elastic-plastic, twodimensional finite element code known as NEPAC(Nonlinear Elastic-Plastic Analysis Code). The code has an element library that consists of three isoparametric twodimensional continuum elements: a three-node constant strain triangular element, a fournode quadrilateral element, and an eight-node quadrilateral element. The elements implemented in NEPAC are members of the serendipity family of isoparametric elements (Ergatoudis, Irons and Zienkiewicz 1968

Read Input Data
Subroutines in this section of the code read and store the input data, compute the reference load vector, and control the restart feature. The input data consists of: control data(e.g. number of nodes and elements, etc.), nodal numbers and locations, element connectivity and material type, load data, and material property data. In addition to storing the input data, the subroutines in this section of the code also compute parameters utilized in the storage and solution schemes. The input data storage scheme utilized in NEPAC is described, and implemented in a 1-d linear truss code, in Bathe .
The applied load input data is utilized in this section of the code to compute the reference load vector.
The code NEPAC has a restart feature. This allows the user to compute the solution for a certain load level, check the results, then restart the code to continue the analysis. Write plot and restart data to unformatted files, and output data to a formatted file. As reported by Nagtegaal (Nagtegaal, Parks and Rice 1973), analyses which combine an incompressible material model(e.g. isochoric perfect plasticity) with a situation where all three direct strain components are computed from the interpolation(i.c. plane strain, axisymmetric and three-dimensional problems) are subject to the problem of meshlocking. Thus, in this type of problem the finite element model exhibits spurious stiffness, i.e. the mesh locks. Several methods of alleviating this difficulty are described in the literature (Belytschko and Tsay 1983, Hughes 1980, Malkus and Hughes 1978, Nagtegaal, Parks and Rice 1973, Taylor, Beresford and Wilson 1976. The method utilized in this code is a selective reduced integration procedure known as the meandilatation formulation (Nagtegaal, Parks and Rice 1973 where [B]dev and [8]du denote the deviatoric and dilatational parts of the straindisplacement matrix, and [B]ave dil is denoted by (8]ave dil = lo. (8] d!l / lo. d!l (C-4) and a is a small stabilization parameter which can be utilized to suppress zero energy dcfonnation modes (Belytschko and Tsay 1983). Thus, if selective reduced integration is utilized the stiffness matrix is computed as (C-5) The three-node triangular element and the four-node quadrilateral clement utilize the mean-dilatation formulation described above. The eight-node quadrilateral utilizes fully reduced integration. Thus, an integration scheme of one order lower than that required to integrate the element stiffness matrix exactly is used. In the case of the eight-node quadrilateral, a 2 by 2 Gauss rule is utilized instead of a 3 by 3 Gauss rule. The main drawback to a fully reduced integration scheme is that spurious zero-energy modes can be intr¢uced into the solution. However, this is not a problem for meshes of reduced integration eight-node quadrilateral elements (ABAQUS Examples Manual 1987).
The constitutive matrix [Dlep is also computed by subroutines in this section of the code. The constitutive equation can be computed by utilizing one of three different material models: an isotropic elastic model, a transversely isotropic elastic model, or an initially isotropic elastic-plastic material model. The elastic-plastic material model assumes J 2 flow theory with isotropic hardening (Chen and Han 1988). A power law hardening decription is utilized to describe the uniaxial stress-strain response of the elastic-plastic material.

Assemble Global Stiffness Matrix
In this section of the code the global stiffness matrix is assembled by the direct stiffness method. The assembled stiffness matrix is then stored in a column vector according to the skyline storage method . Those assembly and storage algorithms are described, and implemented in Fortran in a 1-d linear truss code, in Bathe . This section of the code also contains a subroutine which modifies the global stiffness matrix and load vectors to account for fixed non-zero displacement boundary conditions.

Increment the Applied Load
Since the problems to be analyzed in this project are nonlinear, the analysis must be performed incrementally and/or iteratively. Algorithms presently available to solve nonlinear finite element equations include: Newton-Raphson method, modified Newton-Raphson method, BFGS method, piecewise linear incremental method and the piecewise linear incremental self-correcting method. The following references contain descriptions of some of the more common solution algorithms for nonlinear finite element cquations , Bathe and Cimento 1980, Bathe and Dvorkin 1983, Chen and Han 1988, Kardestuncer 1987, Matthies and Strang 1979. The piecewise linear self-correcting algorithm is utilized in program NEPAC. That algorithm solves nonlinear boundary value problems incrementally.
Several methods are available to increment the load in nonlinear analyses. Probably the most common method is to simply increase the applied loads by equal amounts at the beginning of each load step. However, this approach is not very efficient or versatile.
The drawbacks to this method include the fact that the same size load increment is utilized for each load step; thus leading to a load increment that may be unneccessarily small in the initial stages of the analysis, but too large as the analysis progresses and the structure softens. Also, this algorithm will fail in analyses where it is neccessary to trace the response of the structure into the post-limit point range, i.e. past the maximum load level. There are several so-called continuation algorithm methods available which allow analyses to be continued past the maximum load level. Continuation algorithms presently available include: displacement incrementation via panitioning, displacement incrementation without partitioning,combined finite element/Rayleigh-Ritz method and the fictitious spring approach. A description of several of the more popular continuation algorithms can be found in the following references (Batoz and Dhatt 1979, Bergan et al. 1978, Haisler and Stricklin 1977, Kardestuncer 1987, Noor 1982, Powell and Simons 1981, Ramm 1981, Riks 1972, Zienkiewicz 1971. A description of the combined finite element/Rayleigh-Ritz method can be found in the article by Needleman and Tvergaard (Needleman and Tvergaard, 1983).
Th_ e continuation algorithm implemented in NEPAC is known as the displacement incrementation without partitioning method (Batoz and Dhatt 1979). scheme is known as the piecewise linear method. The vector {R} is computed as (C-8) where {F}y is the total applied load vector, the integral represents the total equilibrated load vector, and Vis the volume of the structure. Thus, {R} is defined as the vector of unequilibrated nodal point forces. Combining equations (C-6) and (C-7), the problem represented by equation (C-6) can be decomposed into the following two problems: Thus, the displacement vector {U} can be written as where the nodal displacement for the qth degree of freedom can be written as After equations (C-9) and (C-10) have been solved, the value of A. can be determined by specifying the incremental displacement at the qth degree of freedom(Uq). Equation (C-11) is then utilized to determine the magnitudes of the remaining incremental nodal point displacements. Thus, the load is incremented, by suitably varying amounts, by incrementing a single control degree of freedom by the same specified value for each load step. The control degree of freedom is selected based upon the characteristics of the specific problem under consideration. It should be noted that this method can be utilized to solve problems involving applied displacements, as well as applied loads. Also, since the only additional computational effort required by this algorithm is the reduction and back-substitution of an additional load vector, the method is quite efficient.
The load incrementation algorithm implemented in NEPAC contains the following improvement. The value of the load parameter A. computed for the 2nd load step is 95 JI I I calculated by determining the load required to cause initial yielding at some location in the structure. Thus, regardless of the size of the control incremental displacement, the linear range of the analysis is always traversed in two load steps. This modification increases the computational efficiency of the algorithm.
Solve System of Simultaneous Equations The system of simultaneous equations is solved by a direct algorithm based on Gauss elimination known as the skyline reduction method. The algorithm is described, and implemented in Fortran, in the finite element text by Bathe . In this solution method the global stiffness matrix is factorized and stored, then the load vector(s) is reduced and back-substituted to obtain the nodal displacements. This algorithm is very computationally efficient because no operations are performed on zero elements outside the skyline. Also, since no elements above the skyline are required in the solution process only elements below the skyline need to be stored. Thus, this algorithm allows the very efficient skyline storage scheme to be utilized to store the global stiffness matrix.

Compute Element Incremental Stresses and Plastic Strains
In this section of the code the element state variables are computed, at the element Gauss points, for each element in the model. The plasticity formulation utilized in NEPAC is the von Mises model with isotropic hardening (Chen and Han 1982). where for the isotropic-hardening von Mises model F=(3Ji.} 112 , where J 2 is the 2nd invariant of the deviatoric stress tensor, and the isotropic hardening parameter k is a 96 function of the plastic work WT'. The strain increment is decomposed into two parts (C-14) where the 1st pan is the elastic strain increment, and the 2nd part is the plastic strain increment. The stresses are related to the elastic strain increments by the constitutive equation (C-15) where the constitutive matrix [D] contains the elastic moduli. Substituting equation (C-14) into equation (C-15) yields where the increment in plastic strains is given by the associated flow rule (C-17) where dA is a scalar multiplier that is determined by the consistency condition df = 0.
By examining equations (C-16) and (C-17) it is seen that once dA is known, the constitutive equation (C-16) is completely defined. The consistency condition yields the following expression for d.A, -18) where H is known as the plastic modulus, and is defined as the slope of the uniaxial stress vs. plastic strain curve. The uniaxial stress vs. strain curve is defined by the power law hardening description (Chen and Han 1988) 97 l I (C-19) where 0 0 and E 0 are the matrix yield stress and yield strain, respectively, (E=CJjE 0 ) and n is the strain-hardening exponent (~1). Due to the form of the uniaxial stress vs. strain curve the expression for the plastic modulus is given by (C-20) where E is the modulus of elasticity, and Oerr is the effective stress.
The constitutive equation represented by equation (C-16) relates an infinitesimal strain increment to an infinitesimal stress increment. In a finite element analysis, however, a finite size load increment is applied which results in a finite size strain increment. Thus, the constitutive equation must be integrated numerically. Several effective algorithms for this purpose have been presented and evaluated in the literature (Chen and Han 1988, Hibbitt 1985, Krieg and Krieg 1977, Ortiz and Popov 1985, Owen and Hinton 1980, Runesson, Sture and William 1988, Schreyer, Kulak and Kramer 1979. The forward Euler method in conjunction with a consistency condition correction algorithm and a subincrementation strategy is utilized in program NEPAC. That algorithm is described and implemented in Fortran, for a linear-hardening von Mises material, in the following texts Hinton 1980, Owen andFawkes 1983). It should be noted, however, that the text (Owen and Hinton 1980)  is then computed. This trial stress increment is added to the total stress value computed at the previous load increment to yield a trial total stress value. This trial total stress value is then substituted into the yield criteria to determine whether the assumption of elastic material behavior was correct. If the assumption of elastic material behavior is found to be correct, the total stress value is set equal to the trial total stress value and, the process is complete. For yielded Gauss points the ponion of the trial stress increment (1r){AOJ which satisfies the yield criteria is determined, where the scaling factor r is defined as (C-23) and where Oy is the current value of the uniaxial yield stress. Equation (C-23) is valid if a previously elastic Gauss point yields during the current load step. Otherwise the value of the scaling factor r is determined by the previous and current values of the yield function. For example, for a Gauss point that was previously plastic and that remains plastic, r is equal to 1; whereas for a Gauss point that was previously plastic and that unloads during the current load step r is equal to 0. It should be noted that other methods of computing the scaling factor are also available (Chen and Han 1988). Once the value of the scaling factor is known the strain increment is divided into two parts, where the 1st part corresponds to a pure elastic response, and the 2nd part corresponds to an elastic-plastic response. The portion of the strain increment corresponding to the elastic-plastic response is then divided into m equal steps; this is known as a subincrementation strategy. At this point the forward Euler method is utilized to numerically integr~te equation (C-16) which yields the stress increment corresponding to the m subincrements of strain. This is accomplished by evaluating the scalar multiplier dA and the incremental plastic strains, which correspond to the total stress value that exists at the beginning of a strain sub-increment This value of incremental plastic strain is then uitlized to evaluate equation (C-16) for the current subincrement, which yields a subincrement of stress. Upon completion of the numerical integration process a consistency condition correction is performed to compensate for any numerical errors which occured in the integration process. This consistency condition correction consists of scaling the updated total stress vector to insure that the yield function is satisfied for the updated value of effective plastic strain. The correction is represented by where the scaling factor a is represented by a = cr/crerr (C-26) where cry represents the updated value of the uniaxial yield stress, and crerr is the effective stress. This correction insures that the updated stress state lies on the yield surf ace. It should be noted that there are other methods of implementing this consistency condition correction (Chen and Han 1988).

Update State Variables
In whether appropriate step sizes have been selected, and in alerting the user to potential problems such as ill-conditioned stiffness matrices (Bergan et al. 1978).

Computer Code Verification
The solutions to several common verification problems were computed with NEPAC.
Since the accuracy of the finite element formulations and numerical algorithms 101 I I implemented in NEPAC have been well documented (as witnessed by the abundance of literature referenced in this appendix), the aim of this section is merely to demonstrate the reliability of the code. Thus, the verification problems selected are designed to show that the selected algorithms are implemented correctly. The results reported pertain to the four-node quadrilateral element implemented with selective reduced integration. Since the eight-node quadrilateral element was not utilized for the thesis research, verification results for this element will not be presented. Verification results are not reponed separately for the triangular element because it is developed by simply collapsing one side of the four-node quadrilateral element. Thus, the only difference between the four-node and three-node elements is the value of the nodal point coordinates input for the element.
Since no modification to the code is required to implement the three-node triangular element, separate verification tests for the three-node element are not required. For the sake of completeness, however, it should be mentioned that the three-node and eight-node elements implemented in NEPAC were utilized in several verification tests; in each case satisfactory results were obtained.
Verification problem number 1 is a simple loading test. This problem verifies the ability of the element to accurately model constant strain states. It also insures that distributed element loads are handled correctly. This problem is described on page 105 for the plane stress and plane strain elements, and on page 106 for the axisymmetric element. The source of this example problem and its analytical solution was the ABAQUS examples manual (ABAQUS Examples Manual 1987).
Verification problem number 2 is the problem of a stretched plate containing a small central hole. This example is used to verify that the element performs well in a problem that involves a severely non-uniform stress field. This problem is solved using a plane stress element. This verification problem is described on page 107. The source of this problem and its analytical solution was the elasticity text by Timoshenko and Goodier (Timoshenko and Goodier 1970).
Verification problem number 3 is the plane strain problem of a thick cylinder subjected to a gradually increasing internal pressure. This problem verifies the plasticity formulation and the ability of the selected element to accurately predict the response of an incompressible material. This example problem is described on page 108. The source of this example problem and its analytical solution was the plasticity text by Prager and Hodge (Prager and Hodge 1951), although it should be noted that the text by Owen and Hinton was the source of some very helpful numerical results for this problem (Owen and Hinton 1980).
Verification problem number 4 is the problem of a single element subjected to a gradually increasing uniaxial load. Since verification problem number 3 involved a perfectly plastic material, the current example problem was required to verify the subroutine utilized to compute the plastic modulus for a nonlinear hardening law. This example problem is described on page 109. The exact solution of this simple example was computed by the author.
The 4-node isoparametric element utilized in this work is a conforming element, provided that full Gaussian quadrature is utilized. Therefore, the interpolating polynomials used for the element satisfy the conditions of completeness and compatibility, and thus, monotonic convergence is guaranteed. However, the use of selective reduced integration renders the element non-conforming(incomparible). Consequently, monotonic convergence is no longer guaranteed. This does not represent a serious problem because the element will still be effective as long as nonmonotonic convergence is assured.

103
I Nonmonotonic convergence is assured if an assemblage of nonconforming elements can represent constant strain states, i.e. if the assemblage of elements satisfies the compl~teness condition. Satisfaction of this condition guarantees that the incompatibilities between elements do not prohibit the assemblage from representing constant strain states.
The patch test is designed to determine whether or not the assemblage of elements satisfies the completeness condition. The patch test involves applying ncxlal point displacements, that correspond to constant states of strain, to a patch of elements. If the computed element strains and stresses correspond to the constant strain states, the element satisfies the patch test. Verification problem 5 is the patch test. As shown in the ABAQUS examples manual, the 4-node element under consideration implemented with the mean-dilatational selective reduced integration scheme is a very effective element.
The purpose of verification problem 5 is not to re-prove that point, but rather to insure that the element is implemented correctly in NEPAC. This example problem is described Analytical Solution: Stresses of CJYY = 300. at point A which corresponds to a stress concentration factor of 3.
Numerical Solution: The stress computed by NEPAC is CJ" = 297.4 at point A, which corresponds to a stress concentration factor of 2.97. This result corresponds to an error of 0.86%.
Comments: This problem was solved with two different meshes, a coarse mesh, and a finer mesh. The convergence properties of the element appear to be satisfactory. Comments: The results computed by program NEPAC agree exactly with the analytical solution. This patch test was also performed for different element aspect ratios, and in all cases excellent agreement with the exact solution was obtained.
111 ! I I I 11