Multi Scale Approach For Composites

A Multi-scale Approach to Modeling Progressive Damage in Composite Structures Tay T.E.*, Liu G., Yudhanto A. and Tan V.B.C. Department of Mechanical Engineering National University of Singapore 9 Engineering Drive 1 Singapore 117576

Abstract: Modeling progressive damage in composite materials and structures poses considerable challenges because the damage is in general complex and involves multiple modes such as delamination, transverse cracking, fiber breakage, fiber pullout, etc. Clearly, damage in composites can be investigated at different length scales, ranging from the micromechanical to the macromechanical specimen and structural scales. In this paper, a simple but novel finite element-based method for modeling progressive damage in fiber-reinforced composites is presented. The element-failure method (EFM) is based on the simple idea that the nodal forces of an element of a damaged composite material can be modified to reflect the general state of damage and loading. This has an advantage over the usual material property degradation approaches in that because the stiffness matrix of the element is not changed, computational convergence is theoretically guaranteed, resulting in a robust modeling tool. The EFM, when employed with suitable micromechanics-based failure criteria, may be a practical method for mapping damage initiation and propagation in composite structures. In this paper, we present a micromechanical analysis for a new failure criterion called the Strain Invariant Failure Theory (SIFT) and the application of the EFM in the modeling of open-hole tension specimens. The micromechanical analysis yields a set of amplification factors which are used to establish a set of micromechanically enhanced strain invariants for the failure criterion. The effects of material properties and volume fraction on the amplification factors are discussed. Keywords: Progressive Damage, Multi-scale Damage, Element-Failure Method, Strain Invariant Failure Theory, Micromechanical Amplification. 1. Introduction Damage progression in composite materials and structures is in general very complicated and involves multiple failure modes, such as fiber breakage, fiber pullout, delamination between plies, matrix cracking, fiber-matrix debonding, etc. Modeling the damage process accurately poses a very difficult problem because these mechanisms clearly operate at various length scales [1]. However, if more rational design and damage tolerance approaches are to be developed for composite structures, it is necessary to develop engineering tools that will enable analysts to model damage and its propagation. At the micromechanics scale, the fibers may be modeled individually and microcracks introduced [2]. However, the large variations in local fiber *

Corresponding Author’s email: [email protected]

distribution and possible flaw sizes quickly renders such models impractical, and invariably, simplifications and idealizations are introduced in the form of representative volume elements (RVE) to make the problem more tractable. The goal of many RVE analyses is to obtain some form of weighted or effective composite property or behavior so that it may be used in constitutive stress-strain relations. Furthermore, it is also not practical nor desirable to model explicitly the very numerous microcracks often observed in damaged composite materials under the microscope. At the macromechanical level of the effective constitutive relations of the damaged material, there have been significant strides in the development of damage mechanics [3, 4]. In most damage mechanics approaches, the crucial component is the damage variables, which are sometimes obtained phenomenologically. In a sense, the traditional continuum damage mechanics approach may be viewed as part of the class of methodology whereby damage is modeled by degradation of selective elastic or mechanical properties. This material property degradation method (MPDM) [5 – 7] essentially involves directly reducing the values of certain material properties in the constitutive relations when and if damage or failure is determined. The severity of damage may be characterized by damage evolution laws which are typically determined through damage mechanics. However, a drawback of this approach is that by reducing the material properties, there is a possibility that the stiffness matrix of the finite element (FE) may inadvertently become ill-conditioned, so that convergence to a solution is not always assured. For a conservative analysis, it is usual to set the values of certain material properties to zero when damage or failure is assumed. For example, if failure is determined to have occurred in the fiber direction (as in the breaking of fibers in tension), the fiber-direction Young’s modulus E11 may be set to zero. However, while not explicitly stated in most published literature, in practice, the material stiffness properties (typically the elastic moduli) are sometimes not completely set to zero but to a certain (often arbitrarily) small percentage of the original value in order to overcome computational difficulties. Another issue involves deciding which of the many material properties to degrade. For example, if transverse cracking or failure is predicted in a composite, it is not clear if only the transverse Young’s modulus E22 or both it and the in-plane shear modulus G12 should be set to zero; additionally, the effect on the Poisson’s ratios ν12 and ν21 is not obvious nor easily determined. The MPDM is popularly used in laminated plate and shell FE codes with failure criteria such as the Tsai-Wu tensor theory because it is convenient to zero material property values in laminate theory, hence the so-called ply discount method. The application of MPDM in three-dimensional problems however, where delaminations and interactions with numerous damage mechanisms are significant, may be cumbersome. At the moment, it is still very difficult to translate information from a micromechanical analysis to guide the development of modified constitutive relations for damaged composite materials. If this can be achieved, the modeling of damage at the components and structures level, supported by the micromechanics of failure and damage, should result in greater accuracy. In this paper, a novel FE-based element-failure method (EFM) is proposed for the modeling of damage in composite structures, in which only the nodal forces are changed to reflect changes in the stress-bearing capability of the damaged material. Because the stiffness matrix remains unaltered, there are no similar computational problems associated with the MPDM. Another consequence is that there should be savings in computational effort since no reformulation of the stiffness matrix with damage progression is involved; each change in the damage state is modeled by appropriately modifying the nodal forces only. Moreover, a new micromechanics-based failure theory, recently proposed by Gosse et al. [8, 9], is used together with the EFM to determine damage initiation and propagation. Called the strain invariant failure theory (SIFT), a unique feature of this failure theory is that information at the micromechanical scale is extracted from micromechanical FE block analyses to amplify strain invariant quantities which are in turn used as criteria to determine failure. In combination with a simple nodal force modification scheme,

the SIFT-EFM approach enables the modeling of damage progression in composite structures that takes into account mechanisms that bridge micro- and macro-length scales. The paper first describes the concept and implementation of the EFM, followed by the strain invariant failure theory (SIFT). It should be noted that while we have chosen SIFT for the work reported herein, the EFM may in general be used with any appropriate failure theory. Indeed, in an earlier work by the authors [10], the Tsai-Wu failure criterion [11] was used with EFM in a plane strain analysis to predict the pattern of delamination growth by low-velocity impact on composite laminates. However, the Tsai-Wu failure criterion (and some of the other more established failure criteria) is best known for application in essentially plane stress analyses involving the classical lamination theory (CLT) of laminated plates. Its extension to the more general three-dimensional cases is considerably more complex and it is difficult to determine the parameters which will ensure closed or admissible failure envelopes in the six-component stress or strain space. While keeping to plane stress analyses significantly simplifies matters and may be sufficient for many design purposes, there is growing realization that damage and delamination in composite structures are truly three-dimensional events with the need for appropriate analytical and modeling tools. SIFT is selected partly because it offers a fully three-dimensional theory with inherent micromechanical features. There are other new micromechanics-based failure theories; an example is the so-called multi-continuum failure theory developed by Mayes et al. [12]. Finally, the paper describes two illustrations of the application of EFM and SIFT in the prediction of damage propagation. The first involves the implementation of EFM in three-dimensional explicit FE, and the problem shown is the evolution of damage in low-velocity impact of composite laminates. The second concerns the EFM’s implementation in three-dimensional implicit FE code. The problem of an open hole in composite laminates is analyzed, and the damage map of each layer is shown. 2. Concept of the EFM The idea and assumption of the EFM is that the effects of damage on the mechanical behaviour can be essentially described by the effective nodal forces of a finite element. It was first developed for dynamic fracture in metals [13], but the modified EFM was used for impact damage in fiber-reinforced composites [10]. The manner by which these effects due to damage translate to the effective nodal forces will in general depend upon the damage evolution law appropriate to the local mode of damage experienced by the composite material, as well as the FE formulation. Traditionally, when damage is assumed to have occurred within an element of a material, the stiffness matrix of the element is altered to reflect the damaged state. One is able to develop explicit relations between the nodal forces and the elastic stiffnesses of an FE: The force-stiffness relation for a finite element is given by

Ku = f

(1)

where u is the vector of nodal displacements, f the vector of nodal forces and K the elemental stiffness matrix of undamaged material, integrated over the domain Ω,

K = ∫ B T CB ∂Ω . Ω

(2)

For two-dimensional (2-D) plane strain or plane stress problems, elements of the material stiffness matrix Cij is defined for i,j = 1,2,6. Although equations for 2-D FE are used for illustrative purpose, the extension to 3-D FE is similar and straightforward. The matrix B is defined as

⎡ N 1,x ⎢ B=⎢ 0 ⎢ N 1,y ⎣

0

N 2 ,x

0

L N m,x

N 1, y N 1, x

0 N 2 ,y

N 2 ,y N 2 ,x

0 L L N m,y

0 ⎤ ⎥ N m,y ⎥ N m,x ⎥⎦

where m is the number of nodes for the element, and N i,x =

(3)

∂N i ∂N i and N i,y = are the ∂x ∂y

derivatives of the shape functions N i with respect to the global x and y coordinates, respectively. The elemental stiffness matrix K can be written as

K 12 L K 1m ⎤ ⎡ K 11 ⎢ K 22 L K 2 m ⎥⎥ ⎢ K= ⎢ O M ⎥ ⎢ ⎥ K mm ⎦ ⎣symm .

(4)

where Κ ιj is a matrix of the form

⎡C11 N i,x N j,x ∂Ω + C 66 N i,y N j,y ∂Ω C12 N i,x N j,y ∂Ω + C 66 N i,y N j,x ∂Ω ⎤ ∫ ∫ ∫ ⎢ ∫ ⎥ Ω Ω Ω K ij = ⎢ Ω ⎥ , (5) ∂ + ∂ ∂ + ∂ C N N Ω C N N Ω C N N Ω C N N Ω i,x j,y i,y j,y i,x j,x 66 ∫ 22 ∫ 66 ∫ ⎢ 12 ∫ i,y j,x ⎥ Ω Ω Ω ⎣ Ω ⎦ for i,j = 1,2,K, m . Substituting Eqns. (2) through (5) into Eqn. (1), the x and y components of the nodal force for node i can be explicitly written as

⎛ m ⎞ f xi = C11 ∫ N i,x ⎜⎜ ∑ N j,x u xj ⎟⎟ ∂Ω Ω ⎝ j ⎠ ⎛ m ⎞ + C12 ∫ N i,x ⎜⎜ ∑ N j,y u yj ⎟⎟ ∂Ω Ω ⎝ j ⎠ m ⎛ m ⎞ + C 66 ∫ N i,y ⎜⎜ ∑ N j,y u xj + ∑ N j,x u yj ⎟⎟ ∂Ω j Ω ⎝ j ⎠

(6)

and

⎞ ⎛ m f yi = C12 ∫ N i , y ⎜⎜ ∑ N j,x u xj ⎟⎟ ∂Ω Ω ⎝ j ⎠ m ⎛ ⎞ + C 22 ∫ N i,y ⎜⎜ ∑ N j,y u yj ⎟⎟ ∂Ω Ω ⎝ j ⎠ m m ⎛ ⎞ + C 66 ∫ N i,x ⎜⎜ ∑ N j,y u xj + ∑ N j,x u yj ⎟⎟ ∂Ω. j Ω ⎝ j ⎠

(7)

Eqns. (6) and (7) show the relationship between nodal forces and the material stiffness properties. It is clear that for every change in the stiffness coefficients Cij, there are corresponding changes in the nodal forces. The EFM, however, works through the direct manipulation of the nodal forces. While it is possible to work out the required nodal force modifications with any given set of changes in material stiffness properties, the reverse is not generally true. In other words, a set of changes to the nodal forces may not translate to a set of identifiable changes in material stiffness properties. For this reason, the EFM is a more general method than the MPDM. Consider an FE of an undamaged composite material (Fig.1 (a)), experiencing a set of nodal forces, which have been obtained from the FE solution of the problem. On the other hand, in an FE containing damaged material, the load-carrying capacity of the FE will be compromised, very likely in a directionally and spatially dependent manner. If much of the damage consists of transverse matrix microcracks, it is reasonable to assume that the FE will have reduced loadbearing capacity in the direction transverse to the fibers (Fig.1 (b)). In conventional MPDM’s, this reduction is achieved by reducing or zeroing certain pertinent material stiffness properties of the damaged finite element. In the EFM however, the reduction is effected by applying a set of external nodal forces such that the nett internal nodal forces of elements adjacent to the damaged element are reduced or zeroed (the latter if complete failure or fracture is implied (Fig.1 (c)). The decision whether to fail an element is guided by a suitable failure theory and in each step, only one or two elements are failed at a time. The required set of applied nodal forces to achieve the reduction within each step is determined by successive iterations until the nett internal nodal forces (residuals) of the adjacent elements converge to the desired values. Typically, less than 200 iterations are required and convergence is guaranteed. Note that it is not the internal nodal forces of the damaged element that is zeroed (for the case of complete failure (Fig.1 (c)), but the nett internal nodal forces of adjacent elements. Thus the “stresses” within the failed element no longer have physical meaning, although compatibility may be preserved. This process leaves the original (undamaged) material stiffness properties unchanged, and is thus computationally efficient as every step and iteration is simply an analysis with the updated set of applied nodal forces. For this reason, it may also be called the nodal force modification method. Hence, no reformulation of the FE stiffness matrix is necessary. 3. Implementation of the EFM The EFM has been implemented in an in-house three-dimensional implicit FE code. Since the method merely involves nodal force modification, it can also be readily implemented in commercial FE codes. However, a commercial general-purpose code was not used to implement EFM because it was not possible to over-ride the internal computational housekeeping associated

with each step of nodal force modification, which would unnecessarily increase the computational time. Hence it was necessary to validate the in-house codes by comparing their solutions for simple problems without damage, with those of the general purpose FE code Abaqus. Fig.2 is the implementation flow-chart of the EFM. It begins with a straightforward linear elastic FE analysis of the problem. The solution in terms of strains is fed into the portion of the program which deals with determining failure; in this case, the program uses SIFT. Within this subroutine, the strains are amplified through amplification factors pre-determined from micromechanical block analyses, and the three strain invariants are calculated. These are in turn compared with their critical values, and failure is determined when any one of them reaches its critical value. A more complete description of SIFT and method of strain amplification is reserved for a subsequent section. The analysis then proceeds to the part of the program that deals with nodal force modification, where the appropriate scheme is applied to all the nodes of the failed element. For the moment, a simple nodal force modification scheme is employed; it zeroes the components of the internal nodal forces that are transverse to the fiber direction. This assumes that the dominant damage mode at the micromechanics level is transverse micro-cracking, which could in reality, be various combinations of micro-fractures within the matrix material and fiber-matrix debonding. This assumption is considered reasonable and conservative for initial-stage damage or damage that is not yet very severe or extensive. The authors are developing more comprehensive nodal force modification schemes that will incorporate fiber failure, so that more extensive damage may be modeled more accurately in the future. The external forces are applied until the desired residual internal nodal force values (in this case, zero) are reached. This process takes typically a few iterations, but convergence is always assured. After convergence has been achieved, the program, guided by the failure theory, searches out the next element that has failed, and performs the nodal modification for that element. If and when no elements are deemed to have failed, the applied load or displacement for the structure is increased until failure is again predicted.

4. The strain invariant failure theory (SIFT) In this section, a description of the new Strain Invariant Failure Theory (SIFT) is given. The EFM, of course, may be used with other failure theories, but SIFT is chosen because it is fully three-dimensional. Recently proposed by Gosse et al. [8, 9], the theory determines if failure has occurred by considering the criticality of three strain invariant values, which have been “amplified” through micromechanical analysis. The procedure is more fully described later, but essentially, the strain components from the homogenous finite element solution are amplified with thermo-mechanical amplification factors extracted from unit cell finite element micromechanical models, before the invariants are calculated. The first of the invariants is J1, defined by

J 1 = ε xx + ε yy + ε zz

(8)

where εxx, εyy, εzz, εxy, εyz and εxz are the six components of the strain vector in general Cartesian coordinates. Where distortional deformation is significant, a criterion based on the second deviatoric strain invariant may be more useful.

The second deviatoric strain invariant J 2′ is defined by

J 2′ =

[

] (

)

1 (ε xx − ε yy )2 + (ε yy − ε zz )2 + (ε xx − ε zz )2 − 1 ε xy2 + ε yz2 + ε xz2 . 6 4

(9)

SIFT employs the von Mises (or equivalent) strain, which is related to the second deviatoric strain invariant by

ε vm = 3J 2′

(10)

These strain invariants are amplified through the use of representative micromechanical blocks, whereby individual fiber and matrix are modeled by three-dimensional finite elements (Fig. 3). Three fiber arrangements or arrays are considered: square, hexagonal and diamond. The diamond arrangement is in fact the same as the square, but rotated through a 45° angle. These representative micromechanical blocks are given prescribed unit displacements in three cases of normal and three cases of shear deformations. For example, in order to obtain strain amplification factors for prescribed displacement in the fiber (or 1- ) direction for one of the faces, the model is constrained in the other five faces (Fig. 4 (a)). The procedure is repeated each time in order to obtain strain amplification factors for displacements in the other two orthogonal (2- and 3- ) directions. Similarly, for shear deformations, the prescribed shear strain is applied in each of the three directions (Fig. 4 (b)). The local micromechanical strains are extracted from various positions within the model (Fig. 5) and normalized with respect to the prescribed strain. In addition to the above mechanical amplification factors, so-called thermo-mechanical amplification factors may be obtained by constraining all the faces from expansion and performing a thermo-mechanical analysis by prescribing a unit temperature differential above the strain-free temperature. In the subsequent analysis, thermal effects are included by prescribing a ∆T of 155°C, between the cure or stress-free temperature of 180°C and the use temperature of 25°C. The locations chosen for the extraction of local amplification factors are shown in Fig. 5(a) for the square array model and Fig. 5(b) for the hexagonal array model. The points F1 through F8 are located at the fiber-matrix interface for values of amplification factors computed using fiber properties, while points M1 through M8, also located at the same positions along the fiber-matrix interface, are for values of amplification factors computed with matrix properties. The point F9 is located at the center of the (assumed circular) fiber, IF1 and IF2 are inter-fiber positions and IS corresponds to the interstitial position. The total number of locations is therefore 20. There are 6 mechanical and 6 thermo-mechanical strain amplification factors for each position; since there are 20 positions and 3 fiber arrangements, the total number of amplification factors is 720 (i.e. 12 × 20 × 3). It should be noted that for a given matrix and fiber material system, the suite of micromechanical block analyses need only be performed once; the resulting amplification factors are stored in a look-up table or subroutine. The output of strains from a macro-finite element analysis is efficiently amplified through this look-up subroutine before the strain invariant values are calculated and compared with the corresponding critical values. The amplification factors for carbon fiber-epoxy system used in all the analyses reported in this paper were obtained from Gosse [9], and coded in a look-up subroutine. However, the amplification factors have been independently verified by the authors, who performed the micromechanical FE block analyses. The matrix and carbon fiber (IM7) material properties are given in Table 1. The subscripts m and f refer to matrix and fiber respectively; the subscript 1 indicates the axial fiber direction, the

subscripts 2 and 3 the transverse directions. The material properties for IM7 and epoxy are obtained from Ha [14] and Gosse et al [8, 9]. The properties for S-glass fiber in the table, used in the subsequent section on the effect of changing materials on the amplification factors, are derived from Gibson [15]. The first strain invariant J1 (Eqn (8)) is calculated with strains amplified only at the IF1, IF2 and IS positions within the matrix material in the micromechanical block. It is generally believed that J1-driven failure is dominated by volumetric changes in the matrix material. On the other hand, the von Mises strain (Eqn (10)) may be amplified with factors not only within the matrix material (IF1, IF2, IS and M1 through M8), but also the fiber and fiber-matrix interface (F1 through F9). m We designate the superscript m for the former case to denote “matrix” (i.e. ε vm ), and the f superscript f for the latter case to denote “fiber” (i.e. ε vm ). SIFT states that failure occurs when

m either of the three strain invariant values reaches its respective critical values (i.e. J 1Crit , ε vm Crit f ), which are determined from analysis of coupon tests of composite laminates with and ε vm Crit various lay-ups [8, 9]. It is clear from the foregoing that for each failure prediction, it is possible not only to determine the invariant that has become critical, but also the position within the micromechanical block model where this has occurred. f It should also be noted that the third critical SIFT criterion ε vm is an effective property, Crit o o obtained through the testing of unnotched 0 unidirectional and [10 /-10o]ns laminates in tension [9]. The critical values of the strain invariants used in the work reported here are derived from the experimental data of Gosse et al. [9] and given in Table 2. In fact, all the SIFT criteria are effective properties of the composite and not only of dry fibers or pure polymer. Analysis of extensive test data have shown that the failure of fibers within a 0o laminate (unnotched tension) and the failure of fibers in a [10o /-10o]ns laminate (also unnotched tension) fail near a common ε vmf more than a common strain along the fiber length. Since the extraction of the SIFT invariants

f is done through the micromechanical block models, we hypothesize that the criterion ε vm incorporates not just the behavior of the fiber phase but also the interfacial behavior (or at least f is not just a fiber-only strength the matrix phase very near the fiber). Viewed in this way, ε vm value. This third invariant can become critical if the fiber material or if the matrix material very close to the fiber is compromised or both. Admittedly, at this present time, it is very difficult to test this hypothesis directly.

5. Effects of changing volume fractions and materials The maximum amplification factors computed from the micromechanical block models with three different volume fractions are shown in Table 3 for the square array and Table 4 for the hexagonal array. In each table, the three orthogonal directions of tensile loading are denoted by “Dir-1”, “Dir-2” and “Dir-3”, while the shear loadings are given by “Dir-12”, “Dir-13” and “Dir23”. Clearly, due to rotational symmetry, the case for Dir-2 yields the same results for Dir-3, and for the square array model, Dir-12 produces the same results for Dir-13. The maximum or highest values of amplification factors are shown in bold fonts, while the next highest are shown below the highest in italics. Similarly, the locations of the occurrence of the maximum amplification factors are denoted by bold fonts, and the locations of the next highest by italic fonts.

We see that in loading modes normal to the direction of the fibers, the maximum amplification factors appear in the inter-fiber regions (IF1 and IF2) suggesting possible failure in the matrix material, although the next highest values occur at the fiber-matrix interface (M1, M5, M3 and M7). For shear cases in the 1-2 and 1-3 planes, amplification factors for the highest and next highest values are extremely close, especially for the fiber volume Vf = 60 % case. This suggests that failure in the case of pure shear is almost equally likely to occur in the matrix (IF1 and IF2) as in the fiber-matrix interface (M1, M5, M3 and M7). For the case of shear across the fibers in the 2-3 plane, failure in the matrix is more likely to be in the interstitial position (IS) although failure in the fiber-matrix interface may still occur. At the high volume fraction of Vf = 70 %, the preferred failure site appears to switch to the interface region from the interstitial position. Generally, maximum amplification factors increase with increasing fiber volume fraction. However, this does not mean that resin-rich composites are necessarily more resistant to damage progression, because the macroscopic composite elastic properties will also change with volume fraction and the critical invariants will also very likely be functions of volume fraction. The effect of changing the fiber material from that of carbon IM7 to S-glass is shown in Table 5 (for square array) and Table 6 (for hexagonal array). The fiber volume in each table is kept at 60 % for comparison. The properties of the fibers and epoxy are given in Table 2. 6. Three-dimensional implicit FE analysis of damage in open-hole tension laminates A 3-D implicit EFM code (with SIFT) was written and applied to the analysis of open-hole composite laminates under uniform remote tension loading. Damage progression for a carbonepoxy cross-ply laminate with a stacking sequence [03/904]s is shown in Figs. 6 (a) and 6 (b). The diameter of the hole is 12.7 mm and the width of the plate is 76.2 mm. The macroscopic elastic properties of the composite are given in Table 2. The mesh consists of three-dimensional 20noded isoparametric brick elements, with one element per ply in the thickness direction. One element per ply in the thickness direction is probably sufficient for modeling in-plane damage, as the following results will show, but is very likely to be inadequate for modeling interlaminar stresses which lead to delamination. However, due to limits in computational power, for the moment, delamination is not modeled, although in principle, with sufficient mesh refinement in the thickness direction, this may be achieved. The nodal force modification employed is simply to zero the component of the force transverse to the fiber direction, once failure is determined by SIFT. This rather simple scheme is reasonable for initial-stage damage. The authors are developing more comprehensive schemes that will include fiber failure. Fig. 6 (a) shows the development of damage emanating from the hole edge in the innermost 90° ply. The dominant f strain invariant which governs the failure of this ply is J1, although some elements did fail by ε vm . The pattern indicates that failure in the rather thick group of 90° plies is predominantly due to a large transverse crack, although other smaller transverse matrix cracks may also arise. Fig. 6 (b) shows the damage progression in the outermost 0° ply. Damage in the 0° plies appears much later than damage in the 90° plies, since the latter are weaker. The pattern suggests matrix cracking in the longitudinal direction, and a picture of a tested specimen with typical longitudinal cracks along the 0° outer ply is also shown for comparison. An alternate representation of damage, with the dominant microcracks drawn in, is shown in Fig. 6 (c), where the damage from the 0° and 90° plies is superposed to form one diagram. This type of representation has the advantage that the damage can be more easily compared with X-ray images. A failed cross-ply specimen with typical longitudinal cracks in the surface 0° ply is shown in Fig. 6 (d).

Damage maps for a laminate with a more general stacking sequence [-45/+45/03/90/0] s under remote uniform tensile load are shown in Fig. 7. The results for the outermost -45° ply are shown in Fig. 7 (a); when compared to the results for the +45° ply (Fig. 7 (b)), the former appears to suffer less damage. This suggests that stacking sequence has a significant effect on damage f propagation. Both ε vm and J1 appear to have contributed to damage in these plies. Damage in the f fifth layer is typical of a 0° ply (Fig. 7 (c)) and is dominated by ε vm . The only 90° ply in the layup (Fig. 7 (d)) exhibits J1 dominated failure; the interesting pattern suggests a distributed system of parallel transverse cracks originating from the edge of the hole. This is in contrast to the single transverse crack found in the 90° ply of the [03/904]s laminate described earlier (Fig. 6 (a)). Clearly, this difference is due to the level of constraint experienced by the 90° ply in the two laminates. It is reasonable to assume that when the level of constraint is high, as in the case of the single ply in the [-45/+45/03/90/0]s laminate, the damage becomes more diffused and less likely to form large single cracks. A similar analysis but with a mesh more refined around the hole is also performed and the results shown in Fig. 8. An encouraging feature is that the results do not show great sensitivity to mesh design and density. This increases the confidence that the predicted damage pattern is indeed a good reflection of the mechanics and physics of damage development. Finally, the matrix cracks representation of damage is shown in Fig. 13. It shows a reasonable distribution of damage that may be expected from an OHT specimen.

7. Conclusion The element-failure or nodal force modification method is introduced for the modeling of damage in composite laminates. It assumes that the deleterious effects of damage on mechanical properties can be effectively modeled through the manipulation of nodal forces, which are then applied to the FE as damage progresses. The method has the advantage that no elements are removed from the FE mesh and therefore has a mechanism to ensure that interpenetration of crack surfaces does not occur. Furthermore, the material stiffness properties are not altered, ensuring that no recalculation of FE stiffness matrices is necessary. This is in contrast to conventional material property degradation techniques, which may result in computational problems because convergence is not always assured and results may be highly mesh-dependent. The EFM has been applied in this paper with a recent micromechanics-based composite strain invariant failure theory (SIFT) to predict damage pattern evolution for open-hole tension problems. This combined SIFT-EFM approach appears to predict reasonable damage patterns emanating from a hole in laminates loaded under remote tension. Furthermore, the results are independent of mesh design and density. Acknowledgement Partial support for this research was provided by the US Air Force Office for Scientific Research. References: 1. Tay T.E., “Characterization and analysis of delamination fracture in composites – an overview of developments from 1990 to 2001”, Appl. Mechs. Revs., 2003, Vol 56, No 1, 1 – 32. 2. Mahishi J.M., “An integrated micromechanical and macromechanical approach to fracture behavior of fiber-reinforced composites”, Engineering Fracture Mechanics, 1986, vol 25, no

2, 197 – 228. 3. Talreja R., “Damage characterization by internal variables”, in Damage Mechanics of Composite Materials, 1994, Elsevier Science, Amsterdam, The Netherlands, (ed. R. Talreja), 53 – 78. 4. Chow C.L., Yang F. and Asundi A., “A three-dimensional analysis of symmetric composite laminates with damage”, Int. Jour. Damage Mechanics, 1993, vol 2, 229 – 245. 5. Tserpes K.I., Papanikos P. and Kermanidis T., “A three-dimensional progressive damage model for bolted joints in composite laminates subjected to tensile loading”, Fatigue Fract Engng Mater Struct, 2001, vol 24, 663 – 675. 6. Camanho P.P. and Matthews F.L., “A progressive damage model for mechanically fastened joints in composite laminates”, Jour. Composite Mater., 1999, Vol 33, No 24, 2248 – 2279. 7. Shokrieh M.M. and Lessard L.B., “Progressive fatigue damage modeling of composite materials, part I: modeling”, Jour. Composite Mater., 2000, Vol 34, No 13, 1056 – 1079. 8. Gosse J.H. and Christensen S., “Strain invariant failure criteria for polymers in composite materials”, AIAA-2001-1184. 2001. 9. Gosse J.H., Christensen S., Wollschlager J.A. and Llanos A.S., “A strain invariant failure theory (SIFT) for composite materials”, Jour. of Composite Mater., (to appear). 10. Tay T.E., Tan V.B.C. and Deng M., “Element-failure concepts for dynamic fracture and delamination in low-velocity impact of composites”, Int. Jour. of Solids & Structures, 2003, Vol 40, No 3, 555 – 571. 11. Tsai S.W., Theory of Composites Design, 1992, Think Composites, Dayton, USA. 12. Mayes J.S. and Hansen A.C., “Multicontinuum failure analysis of composite structural laminates”, Mechanics of Composite Materials and Structures, 2001, Vol 8, 249 – 262. 13. Beissel S.R., Johnson G.R. and Popelar C.H., “An element-failure algorithm for dynamic crack propagation in general directions”, Engineering Fracture Mechanics, 1998, Vol 61, No 3-4, 407-425. 14. Ha S.K., “Micromechanics in the analysis of composite structures”, 6th Composites Durability Workshop (CDW-6), Tokyo, Japan, 14-15 November 2002. 15. Gibson R.F., Principles of Composite Material Mechanics, 1994, McGraw-Hill, Inc., New York, USA.

Table 1 Fiber and epoxy elastic properties used in micromechanical block models. Ef1 Ef2 Ef3 νf12 νf13 νf23 Gf12 Gf13 Gf23

αf1 αf2 αf3

IM7 carbon fiber 303 GPa 15.2 GPa 15.2 GPa 0.2 0.2 0.2 9.65 GPa 9.65 GPa 6.32 GPa 0 8.28 × 10-6 /°C 8.28 × 10-6 /°C

S-glass fiber 85.50 GPa 85.50 GPa 85.50 GPa 0.2 0.2 0.2 35.65 GPa 35.65 GPa 35.65 GPa 5.04 × 10-6 /°C 5.04 × 10-6 /°C 5.04 × 10-6 /°C

Epoxy Em

3.31 GPa

νm

0.35

αm

5.76 × 10-5 /°C

Table 2 Material properties of carbon-epoxy composite used in OHT FE models. Modulus in fiber direction E1 (GPa) Transverse moduli E2 = E3 (GPa) Shear moduli G12 = G13 (GPa) Shear modulus G23 (GPa) Thermal expansion coefficient in fiber direction α1 (/°C) Thermal expansion coefficients in transverse direction α2 = α3 (/°C) Critical invariant J 1 Crit

161.3 8.3 5.16 3.38 0.01 × 10-6 32.7 × 10-6 0.0274

m Critical invariant ε vm Crit

0.103

f Critical invariant ε vm Crit

0.0182

Table 3 Effect of fiber volume fraction Vf on amplification factors in square array model (figures in bold are maximum values; figures in italics for next highest values. Positions are marked in Fig. 5). Fiber volume fraction

Vf = 50%

Maximum amplification factor Position

Vf = 60%

Maximum amplification factor Position

Vf = 70%

Maximum amplification factor Position

Dir-1

Dir-2

Dir-3

Dir-12

Dir-13

Dir-23

1.0

2.494 2.012

2.494 2.012

3.308 3.049

3.308 3.049

2.280 2.041

All points

IF1 M1, M5

IF2 M3, M7

IF1 M1,M5

IF2 M3,M7

IS M1,M3,M5,M7

1.0

2.897 2.383

2.897 2.383

4.662 4.639

4.662 4.639

2.623 2.575

All points

IF1 M1, M5

IF2 M3, M7

M1, M5 IF1

M3, M7 IF2

IS M1,M3,M5,M7

1.0

3.156 2.771

3.156 2.771

7.502 7.347

7.502 7.347

3.904 3.747

All points

IF1 M1,M5

IF2 M3,M7

IF1 M1,M5

IF2 M3,M7

M1,M3,M5,M7 IF1,IF2

Table 4 Effect of fiber volume fraction Vf on amplification factors in hexagonal array model (figures in bold are maximum values; figures in italics for next highest values. Positions are marked in Fig. 5). Fiber volume fraction

Vf = 50%

Maximum amplification factor Position

Vf = 60%

Maximum amplification factor Position

Vf = 70%

Maximum amplification factor Position

Dir-1

Dir-2

Dir-3

Dir-12

Dir-13

Dir-23

1.0

1.843 1.696

2.346 1.869

2.691 2.061

3.023 2.710

2.437 1.948

All points

M2,M4,M6,M8 IF2

IF1 M3,M7

M1,M5 IF2

M3,M7 IF1

M1,M5 IF2

1.0

2.050 1.833

2.786 2.160

3.360 2.428

3.529 3.524

2.913 1.899

All points

M2,M4,M6,M8 IF2

IF1 M3,M7

M1,M5 IF2

M3,M7 IF1

M1,M5 M3,M7

1.0

2.177 2.101

2.880 2.482

3.357 3.199

3.767 3.594

2.712 2.085

All points

IF2 M2,M4,M6,M8

IF1 M3,M7

IF2 M1,M5

M3,M7 IF1

M1,M5 M3,M7

Table 5 Effect of materials on amplification factors in square array model (figures in bold are maximum values; figures in italics for next highest values. Positions are marked in Fig. 5).

IM7epoxy

Position

S-glassepoxy

Dir-1

Dir-2

Dir-3

Dir-12

Dir-13

Dir-23

1

2.897 2.383

2.897 2.383

4.662 4.639

4.662 4.639

2.623 2.575

All points

IF1 M1, M5

IF2 M3, M7

M1, M5 IF1

M3, M7 IF2

IS M1,M3,M5,M7

1

6.952 5.442

6.952 5.442

6.878 6.754

6.878 6.754

3.892 3.470

All points

IF1 M1, M5

IF2 M3, M7

M1, M5 IF1

M3, M7 IF2

M1,M3,M5,M7 IS

Maximum amplification factor

Maximum amplification factor Position

Table 6 Effect of materials on amplification factors in hexagonal array model (figures in bold are maximum values; figures in italics for next highest values. Positions are marked in Fig. 5).

IM7epoxy

Maximum amplification factor Position

S-glassepoxy

Maximum amplification factor Position

Dir-1

Dir-2

Dir-3

Dir-12

Dir-13

Dir-23

1

2.050 1.833

2.786 2.160

3.360 2.428

3.529 3.524

2.913 1.899

All points

M2,M4,M6,M8 IF2

IF1 M3,M7

M1,M5 IF2

M3,M7 IF1

M1,M5 M3,M7

1

3.157 3.134

5.248 3.921

4.330 4.309

4.778 4.731

4.403 3.387

All points

M2,M4,M6,M8 IF

IF1 M3,M7

M1,M5 IF2

M3,M7 IF1

M1,M5 IF2

List of figures Fig. 1. (a) FE of undamaged composite with internal nodal forces (b) FE of composite with transverse matrix cracks. Components of internal nodal forces transverse to fiber direction are modified. (c) Completely failed or fractured element. All nett internal nodal forces of adjacent elements are zeroed. Fig. 2. Damage and delamination in composite specimen subjected to 3-point bend, and comparison with EFM predicted damage pattern. Fig. 3. FE models of micromechanical blocks with (a) square (b) hexagonal and (c) diamond packing arrays. Fig. 4. (a) Prescribed normal displacements, (b) prescribed shear deformations. Fig. 5. Locations for the extraction of amplification factors within the micromechanical block models: (a) square array (b) hexagonal array. Fig. 6. Damage maps for cross-ply [03/904]s laminate: (a) Innermost 90° ply; (b) Outermost 0° ply; (c) Representation by matrix cracks; (d) Damaged cross-ply specimen. Fig. 7. Damage maps for [-45/+45/03/90/0]s laminate: (a) -45° ply; (b) 45° ply; (c) 0° ply; (d) 90° ply. Fig. 8. Damage maps for [-45/+45/03/90/0]s laminate (refined mesh around hole): (a) -45° ply; (b) 45° ply; (c) 0° ply; (d) 90° ply. Fig. 9. Matrix cracks representation for the [-45/+45/03/90/0]s laminate.

Fiber direction

(a) FE of undamaged material and nodal force components

Fiber direction

(b) Partially failed FE with damage and modified nodal forces

(c) Completely failed FE with extensive damage

Fig. 1. (a) FE of undamaged composite with internal nodal forces (b) FE of composite with transverse matrix cracks. Components of internal nodal forces transverse to fiber direction are modified. (c) Completely failed or fractured element. All nett internal nodal forces of adjacent elements are zeroed.

Perform FEA Recover strains

Apply SIFT and determine next failed element

Find nodal forces of failed elements Proceed to fail additional elements until the desired number is reached

Apply external nodal forces at nodes Continue iterations Solve for new displacements and new nodal forces

Yes Check Convergence

Fig. 2. Implementation flowchart of EFM.

No

90˚ 60˚

(a)

(b)

45˚

(c)

Fig. 3. Micromechanical blocks with (a) square (b) hexagonal and (c) diamond packing arrays.

3

2

1 (a)

(b) Fig. 4. (a) Prescribed normal displacements, (b) prescribed shear deformations.

IF1

Matrix 3

IS

M1

M8

M2

F1 F2

F8 F9

2 M7

F3

F7

M3

IF2

Fiber F6 M6

1

IF1

F4

F5

M4

M5

(a)

IF1 IS M1

3 M8 2 M7

IF2 M2 F1 F8 F9 F2 F3 M3 F7 F6

M6 1

F5

F4 M4

M5

(b)

Fig. 5. Locations for the extraction of amplification factors within the micromechanical block models: (a) square array (b) hexagonal array.

J1 m ε vm ε vmf

(a) 90° ply (6th ply)

(c) Matrix cracks representation

(b) 0° ply (surface ply)

(d) A failed cross-ply specimen

Fig. 6. Damage maps for cross-ply [03/904]s laminate: (a) Innermost 90° ply; (b) Outermost 0° ply; (c) Representation by matrix cracks; (d) Damaged cross-ply specimen.

J1 m ε vm ε vmf

(a) -45° ply

(c) 0° (fifth) ply

(b) 45° ply

(d) 90° ply

Fig. 7. Damage maps for [-45/+45/03/90/0]s laminate: (a) -45° ply; (b) 45° ply; (c) 0° ply; (d) 90° ply.

J1 m ε vm

ε vmf

(a) -45° ply

(c) 0° (fifth) ply

(b) 45° ply

(d) 90° ply

Fig. 8. Damage maps for [-45/+45/03/90/0]s laminate (refined mesh around hole): (a) -45° ply; (b) 45° ply; (c) 0° ply; (d) 90° ply.

Fig. 9. Matrix cracks representation for the [-45/+45/03/90/0]s laminate.