Simulation Tool for Damage in Composite Laminates

Abstract

A numerical simulation tool for progressive failure in laminates utilizes a low fidelity approach. The numerical model includes an enriched element that is initially in a low fidelity form. The enriched elements may increase fidelity by splitting locally to simulate an ongoing damage process such as delamination.

Claims

1. A computerized method of simulating damage growth in composite laminates having a plurality of interfaces between adjacent layers of fibers disposed in a matrix material and at least one delamination crack having a delamination front at a current interface, the method comprising: forming a finite element model of a composite material by providing a plurality of nodes and an associated single layer of elements comprising shell or plate finite elements, wherein the elements comprise stiffness matrices corresponding to a plurality of layers of fibers disposed in a matrix material; storing the finite element model in computer memory; using a computer to determine a finite element solution to the finite element model subjected to a load; providing separate current stiffness matrices corresponding to laminate layers on each side of a delamination crack at the current interface; utilizing predefined criteria to determine if the delamination crack grows as a result of a most recent incremental finite element load step solution as part of a finite element model solution for an applied load; if the delamination crack does grow, predicting the spatial orientation of microcracks that precede the delamination front and evaluating if an adjacent layer of fibers arrests the microcracks to thereby determine if the delamination front remains at the current interface or migrates in a direction that is transverse to the elements through at least one of the adjacent layers of fibers via a matrix crack to a new interface; if the delamination crack grows and remains at the current interface, split one or more adjacent elements in the direction of growth into new separate stiffness matrices that are the same as the current stiffness matrices whereby the new separate stiffness matrices correspond to laminate layers on each side of the delamination crack at its current interface; if the delamination grows and migrates to a new interface that is offset from the current interface, split one or more adjacent elements in the direction of growth into new separate stiffness matrices that are modified relative to the current stiffness matrices and correspond to the new interface location that is offset from the current interface to thereby account for migration of the delamination to the new interface.

2. The method of claim 1, wherein: forming a finite element model includes forming a plurality of floating nodes for each single layer shell or plate element.

3. The method of claim 2, wherein: the separate stiffness matrices are formed utilizing the floating nodes.

4. The method of claim 1, including; forming tied nodes at the delamination front; predicting the spatial orientation of microcracks that precede the delamination front utilizing tie forces at the tied nodes.

5. The method of claim 4, including: utilizing tie forces of the tied nodes to predict if an adjacent layer of fibers arrests the microcracks; and splitting one or more adjacent elements in the direction of growth at the new interface if an adjacent layer of fibers does not arrest the microcracks.

6. The method of claim 1, wherein: if the delamination crack grows according to a first delamination growth prediction at the current interface or a new interface, the interface location of the adjacent element that is determined to be split is superseded by a second delamination growth prediction that has also determined the adjacent element to be split but at a different interface if the second delamination growth prediction is energetically dominant.

7. The method of claim 1, including: using a computer to determine at each load step if a new delamination crack forms according to a predefined damage initiation criteria.

8. The method of claim 1, wherein: the finite element model includes a plurality of delamination cracks; and including: forming separate stiffness matrices corresponding to the laminate layers on each side of each delamination crack.

9. The method of claim 8, including: utilizing predefined criteria to determine if the delamination cracks grow at a current interface corresponding to each delamination crack as a result of a most recent finite element solution to an applied load.

10. The method of claim 1, wherein: the differences in the stiffness matrices define a stiffness discontinuity at a boundary between adjacent shell elements and thereby model a matrix crack.

11. A computerized method of simulating damage growth in composite laminates having a plurality of interfaces between adjacent layers of fibers disposed in a matrix material, the method comprising: forming a finite element model of a composite material, the finite element model comprising a plurality of nodes and having a region having a single layer of shell or plate elements comprising a composite stiffness matrix corresponding to a plurality of layers of fibers disposed in a matrix material; storing the finite element model in computer memory; using a computer to determine a finite element solution to the finite element model subjected to a load; if the finite element model has a delamination crack, providing separate stiffness matrices corresponding to laminate layers on each side of the delamination crack; predicting the spatial orientation of microcracks that precede the delamination front to thereby determine if the delamination front migrates via a matrix crack to a new interface; if the delamination migrates to a new interface, split one or more adjacent shell elements in the direction of growth to form new shell elements, the new shell elements having stiffness matrices that are modified from those of the current interface split location to correspond to a new interface location that is offset from the current interface to thereby account for migration of the delamination to the new interface by providing a stiffness discontinuity at a boundary between adjacent shell elements thereby modeling the matrix crack.

12. The method of claim 11, including; forming tied nodes at the delamination front; predicting the spatial orientation of microcracks that precede the delamination front utilizing tie forces of the tied nodes.

13. The method of claim 12, including: utilizing tie forces of the tied nodes to predict if an adjacent layer of fibers arrests the microcracks.

14. The method of claim 11, including: utilizing predefined criteria to determine if the delamination crack grows at the current interface between adjacent layers of fibers as a result of a most recent incremental finite element load step solution as part of a finite element model solution for an applied load.

15. The method of claim 14, wherein: if the delamination crack grows, with or without migration, the interface location of an adjacent element that is determined to be split may be superseded by another nearby element.

16. A computerized method of simulating damage growth in composite laminates having a plurality of interfaces between adjacent layers of fibers disposed in a matrix material, the method comprising: forming a finite element model of a composite material, the finite element model comprising a plurality of nodes and having a region having a single layer of shell or plate elements comprising a composite stiffness matrix corresponding to a plurality of layers of fibers disposed in a matrix material; storing the finite element model in computer memory; using a computer to determine a finite element solution to the finite element model subjected to a load; using a computer to evaluate at each load step in a numerical finite element solution procedure if a delamination crack exists in the finite element mesh and/or if a delamination crack forms according to a predefined damage initiation criteria; if a delamination crack forms or is already present, determining the location in the laminate including a present interface; if a delamination crack does not form, or if a delamination crack does not exist, repeat the step of determining a finite element solution; if a delamination crack forms or already exists, providing separate stiffness matrices corresponding to laminate layers on each side of the delamination crack; predicting the spatial orientation of microcracks that precede the delamination front and evaluating if an adjacent layer of fibers arrests the microcracks to thereby determine if the delamination front remains in the current interface or migrates in a direction that is traverse to the elements through at least one of the adjacent layers of fibers via a matrix crack to a new interface; if the delamination grows and migrates to a new interface, split one or more adjacent shell elements to form new shell elements having stiffness matrices to thereby account for migration of the delamination to the new interface.

17. The method of claim 16, including: providing a stiffness discontinuity at a boundary between adjacent shell elements thereby modeling the matrix crack.

18. The method of claim 16, including: utilizing predefined criteria to determine if the delamination crack grows at the current interface between adjacent layers of fibers as a result of a most recent incremental finite element load step solution as part of a finite element model solution for an applied load.

19. The method of claim 18, including: if the delamination crack grows and remains at the current interface, split one or more adjacent elements in the direction of growth into separate stiffness matrices corresponding to laminate layers on each side of the delamination crack at its current interface.

20. The method of claim 16, including, forming tied nodes at the delamination front; predicting the spatial orientation of microcracks that precede the delamination front utilizing tie forces of the tied nodes at a delamination front.

Description

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

[0012] FIG. 1 is a flow chart showing a process of simulating damage in composite laminates according to one aspect of the present invention;

[0013] FIG. 2 is a schematic view showing a Mindlin composite shell element;

[0014] FIG. 3 is a schematic view of an undamaged shell element;

[0015] FIG. 4 is a schematic view of a split shell element corresponding to delaminations in a composite material;

[0016] FIG. 5 is a schematic view showing coordinates and variables utilized in a Virtual Crack Closure Technique (VCCT);

[0017] FIG. 6 is a schematic view showing implementation of a Virtual Crack Closure Technique (VCCT) showing opposing nodes tied together according to element state in relation to a crack tip or delamination front;

[0018] FIG. 7 is a schematic view showing an evaluation of delamination propagation in a finite element model;

[0019] FIG. 8 is a schematic view including a physical schematic of delamination migration that has occurred and a corresponding representation in a shell element model where a stiffness discontinuity at an element boundary corresponds to a matrix crack in the physical schematic;

[0020] FIG. 9 is a schematic view showing delamination that has migrated in a composite laminate;

[0021] FIG. 10 is a schematic view showing cusp orientation and adjacent fiber orientation; and

[0022] FIG. 11 is a physical schematic view showing one potential delamination-migration energy criteria for use in a shell element model.

DETAILED DESCRIPTION OF THE INVENTION

[0023] For purposes of description herein, the terms “upper,” “lower,” “right,” “left,” “rear,” “front,” “vertical,” “horizontal,” “positive,” “negative,” and derivatives thereof shall relate to the invention as oriented in FIG. 2. A summary of the process is shown in FIG. 1. However, it is to be understood that the invention may assume various alternative orientations and step sequences, except where expressly specified to the contrary. It is also to be understood that the specific devices and processes illustrated in the attached drawings, and described in the following specification, are simply exemplary embodiments of the inventive concepts defined in the appended claims. Hence, specific dimensions and other physical characteristics relating to the embodiments disclosed herein are not to be considered as limiting, unless the claims expressly state otherwise.

[0024] FIG. 1 depicts a process 1 according to one aspect the present invention. Following the start of process 1, a shell finite element model is created at step 2 and a load is applied. As discussed in more detail below, the shell finite element model is initially a model with a low-fidelity mesh including a plurality of adjacent (side-by-side, not “stacked”) shell elements. The shell elements include a stiffness matrix corresponding to a composite laminate having a plurality of layers of fibers and a matrix material. The shell finite element model may initially include a single layer mesh composed of a plurality of single-layer shell elements with “floating” nodes that permit the shell element to be split as required to simulate delamination and delamination-matrix crack interaction. As used herein, the term “shell” generally refers to any finite element (e.g. shell, plate, etc.) capable of numerically modeling damage propagation in composite laminates as described herein. It will be understood that the present invention is not limited to shell elements.

[0025] Referring again to FIG. 1, at a next step 4, a load increment is solved. The process 1 may then proceed directly to step 14 (discussed below), or the process may proceed to optional steps or features shown in the dashed region 12. As discussed in more detail below, the optional process features relate to determining if a crack or flaw forms, and includes evaluating damage initiation criterion 6 throughout the model. This is followed by determining if a crack or flaw is formed at step 8. If a crack or flaw does form the shell elements are split as needed (preferably only as needed) to insert the crack or flaw at step 10, and the process proceeds to step 16. If no crack or flaw is formed at step 8, the optional processes proceeds to step 14. At 14, the process involves determining if there is an existing crack or flaw. If not, the process returns to 4, and the next load increment is solved. However, if it is determined at step 14 that there is an existing crack or flaw, the process proceeds to step 16. As discussed in more detail below, predefined criteria may be utilized to determine if a delamination grows at the crack or flaw boundary. If not, the process returns to 4, and the next load increment is solved. If the delamination does grow, the process then involves determining if the delamination migrates via a matrix crack at step 18. If the delamination crack does migrate, the next elements in the direction of the growth are split at a new interface as indicated at step 20, and the process then continues with the next load increment/solution 4. As discussed in more detail below, the stiffness matrices of the laminate at the migrated interface are not equal to the stiffness matrices at the prior interface to thereby account for the migration of the delamination crack via a matrix crack in the mesh at the element boundary.

[0026] If a delamination does not migrate at 18, the process continues as shown at 22. At 22, the next elements are split in the direction of growth at the current interface, and the process continues with the next load increment/solution 4. It will be understood that the process includes checking if the solution is complete after 20 and 22, and ending the process if the solution is complete. If the solution is not complete after 20 or 22, the process returns to 4 as shown in FIG. 1. It will be understood that the present invention is not necessarily limited to the steps and sequences described above in connection with FIG. 1. In general, there may be additional steps not shown in FIG. 1 following steps 20 and 22 for scenarios including: 1) when multiple delaminations interact with one another and may or may not link together via matrix cracks; and 2) an element that has migration predicted in one region but not in another in the same load increment.

[0027] With further reference to FIG. 2, a Mindlin composite shell element includes 4 nodes N1-N4. The stiffness matrix integration can be expressed utilizing the following general shell element stiffness integration equation:

[00001]$\begin{array}{cc}{K}^{\left(e\right)}={\int }_x.Math.{\int }_y.Math.H^T\left[{\int }_z.Math..Math.\mathrm{Cdz}\right].Math.H.Math..Math.\mathrm{dydx}={\int }_b.Math.{\int }_b.Math.H^T\left[{\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.\mathrm{Cdz}\right].Math.H.Math..Math.\mathrm{dydx}& \left(1.0\right)\end{array}$

where:
H=strain-displacement matrix
C=constitutive material matrix
b=edge dimension of an element
h=total thickness of the laminate

[0028] The laminate theory constitutive material matrices are as follows:

[00002]$\begin{array}{cc}\begin{array}{cc}A={\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.C_p.Math.\mathrm{dz},& B={\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.{\mathrm{zC}}_p.Math.\mathrm{dz}\\ D={\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.z^2.Math.C_p.Math.\mathrm{dz},& G={\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.C_s.Math.\mathrm{dz}\end{array}& \left(2.0\right)\end{array}$

wherein C.sub.p and C.sub.s=planar and shear constitutive material matrices.

[0029] The laminate shell element stiffness integration is as follows:

[00003]$\begin{array}{cc}{K}^{\left(e\right)}={\int }_b.Math.{\int }_b.Math.H^T\left[{\int }_{-\frac{h}{2}}^{\frac{h}{2}}.Math.C\left(A,B,D,G\right).Math.\mathrm{dz}\right].Math.H.Math..Math.\mathrm{dydx}& \left(3.0\right)\end{array}$

[0030] When a discontinuity is introduced at a single uniform z-coordinate, as in the case of a delamination, the element material is split into two regions, Ω.sub.A and Ω.sub.B, where region Ω.sub.B corresponds to the DOF of the floating nodes. The stiffness matrix for a split element is given by:

K.sup.(e)=K.sub.Ω.sub.A.sup.(e)+K.sub.Ω.sub.B.sup.(e)  (4.0)

[0031] With further reference to FIG. 3, an undamaged shell element 24 initially includes nodes N1-N4 that comprise real nodes (RN) and unused floating nodes (FN). Floating nodes and real nodes are generally known in the art (see, Chen, “The Floating Node Method (FNM)” supra) Thus, a detailed description of these terms is not believed to be required.

[0032] The stiffness matrix for the undamaged shell element 24 is given by:

[00004]$\begin{array}{cc}{K}^{\left(e\right)}={.Math.\begin{array}{cc}{\left[K_{{\Omega }_A}^{\left(e\right)}\right]}_{24\times 24}& \left[0\right]\\ \left[0\right]& {\left[0\right]}_{24\times 24}\end{array}.Math.}_{48\times 48}& \left(5.0\right)\end{array}$

[0033] The configuration of equation 5.0 allows for one delamination. However, it will be understood that the same approach may be utilized to allow for additional delaminations if required.

[0034] With further reference to FIG. 4, the shell element 24 can be split into first and second overlapping shell elements or components 24A and 24B. Shell element 24A includes real nodes RN1-RN4, and the shell element 24B includes floating nodes FN1-FN4. The stiffness matrix for the split elements 24A and 24B of FIG. 4 is as follows:

[00005]$\begin{array}{cc}{K}^{\left(e\right)}={.Math.\begin{array}{cc}{\left[K_{{\Omega }_A}^{\left(e\right)}\right]}_{24\times 24}& \left[0\right]\\ \left[0\right]& {\left[K_{{\Omega }_A}^{\left(e\right)}\right]}_{24\times 24}\end{array}.Math.}_{48\times 48}& \left(6.0\right)\end{array}$

[0035] In general, the undamaged shell element 24 (FIG. 3) may have a stiffness matrix corresponding to a composite laminate including a plurality of layers of reinforcing fibers disposed at different orientations. The shell element 24 may be split to form the shell elements or components 24A and 24B (FIG. 4) to simulate a delamination between the shell elements 24A and 24B. The stiffnesses of the shell elements 24A and 24B correspond to the portions of a composite laminate disposed on either side of a delamination crack. Shell elements or components 24A and 24B may be considered to be separate elements or they may be considered to be separate parts or components of a single element.

[0036] With further reference to FIG. 5, a finite element model 50 may utilize the Virtual Crack Closure Technique (VCCT). Energy release rate equations originally given by Wang, J T, and Raju, I S. “Strain energy release rate formulas for skin-stiffener debond modeled with plate elements.” Engineering Fracture Mechanics 54.2 (1996): 211-228, may be utilized. These equations are shown here corresponding to a delamination growth direction along the positive x direction as described in FIG. 5. Using delamination growth in the positive x-direction as an example, Mode I, Mode II, and total energy release rates at node (i,j) are determined as follows:

[00006]$\begin{array}{cc}{G}_1^{\left(+x\right)}=\frac{-1}{2.Math.\Delta .Math..Math.A^{\left(e\right)}}[F_z^{\left(i,j\right)}\left(w^{{\left(i-1,j\right)}^{\prime }}-w^{\left(i-1,j\right)}\right)+M_x^{\left(i,j\right)}\left({\theta }_x^{{\left(i-1,j\right)}^{\prime }}-{\theta }_x^{\left(i-1,j\right)}+M_y^{\left(i,j\right)}\left({\theta }_y^{{\left(i-1,j\right)}^{\prime }}-{\theta }_y^{\left(i-1,j\right)}\right)\right]& \left(7.0\right)\end{array}$

where F, M, w, u, and θ are nodal force, moment, z-displacement, x-displacement, and rotation, respectively. Nodal designations with a prime superscript refer to the “upper” set of elements (i.e., floating nodes) and crack extension area, ΔA.sup.(e)=b.sup.2, is the area of an element for a square regular mesh. It will be understood that other mesh shapes may be utilized.

[0037] As described in Benzeggagh, M. L., M. Kenane. 1996. “Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus,” Composites Science and Technology, 56(4):439-449, the mixed-mode critical energy release rate, G.sub.c may be calculated using the Benzeggagh-Kenane criterion as follows:

G.sub.c.sup.(+x)=G.sub.Ic+(G.sub.IIc−G.sub.Ic)(G.sub.II.sup.(+x)/G.sub.T.sup.(+x)).sup.η.sup.BK   (10.0)

[0038] A summary of the VCCT is disclosed in Kreuger, “Virtual crack closure technique: History, approach, and applications, Applied Mechanics Review,” supra. Thus, a detailed description of the VCCT is not believed to be required.

[0039] With further reference to FIG. 6, the enriched shell element 24 can be used in a mesh to represent four different damage states. These are illustrated in FIG. 6 and consist of 1) undamaged element, 2) split element ahead of the crack tip with all nodes tied, 3) split element at the crack tip with one or more nodes free, and 4) split element in the crack wake with all nodes free. At the end of every solution increment, VCCT is used on each delamination front tied node to determine if the tie should remain in place or be released resulting in a change in the adjacent element damage states (see Orifici, A. C., R. S. Thomson, R. Degenhardt, S. Busing, J. Bayandor. 2007. “Development of a Finite Element Methodology for Modelling Mixed-mode Delamination Growth in Composite Structures,” 12.sup.th Australian International Aerospace Congress, Australia). Thus, the Virtual Crack Closure Technique (VCCT) may be used at a delamination crack tip or delamination front 30 utilizing tied nodes TN. The tied nodes TN comprise real nodes RN that are tied to floating nodes FN. One method to tie nodes together is by inserting high stiffness terms into the stiffness matrix coupling translational degrees of freedom. The undamaged elements 24 as labeled in FIG. 6, preceding the delamination crack tip or front 30 include nodes N that utilize real nodes RN only. Split elements 32 immediately preceding crack tip or front 30 include nodes TN at each corner of the split elements 32. As discussed above, each TN comprises a real node RN that is tied to a floating node FN. Alternatively, more than one row of split elements with all nodes tied may be defined ahead of a delamination to improve model accuracy and mesh independence.

[0040] The split elements include lower and upper element components 34A and 34B, respectively. The lower components 34A and upper components 34B adjacent to the “open” side of delamination crack tip or front 30 include tied nodes TN at the crack tip 30. The components 34A include real nodes RN along boundaries 36A away from the crack tip 30, and the components 34B include floating nodes FN along boundaries 36B away from the crack tip 30. The split components 38A and 38B away from the crack tip 30 include floating nodes FN and real nodes RN that are all free/not tied. In general, an offset is applied to the elements 24A, 24B, etc. that have been split coupling certain rotational and membrane degrees of freedom to account for the offset of the material on each side of a delamination from the original undamaged element's neutral axis. When using the FNM, opposing components such as 34A and 34B, are actually the same element with two separate regions or components.

[0041] An example of a VCCT is disclosed in Krueger, “Virtual crack closure technique: History, approach, and applications, Applied Mechanics Review,” supra. Thus, a detailed description of VCCT is not believed to be required.

[0042] With further reference to FIG. 7, a finite element model 50 similar to the model 50 of FIG. 6 but with a staggered delamination front is shown in plan view. Delamination propagation is evaluated, as show in FIG. 7, along directions X and Y aligned with the mesh. The total energy release rate (G.sub.T) is calculated along the mesh lines or boundaries 40 in four directions at each tied node TN on a delamination front. If the total energy release rate (G.sub.T) in any of the four propagation directions is greater than the critical energy release rate (G.sub.C), the tie is released.

[0043] Although the propagation of a delamination crack may be determined utilizing the VCCT approach as described previously and as shown in FIGS. 6 and 7, it will be understood that other VCCT approaches such as described in Xie, D., S. B. Biggers Jr. “Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part I: Formulation and validation,” Engineering Fracture Mechanics, 2006, 73:771-785 or a Cohesive Zone (CZ) approach as described in Wisnom, M R, “Modelling discrete failure in composites with interface elements,” Composites: Part A, 2010, 41:795-805 may also be utilized.

[0044] With further reference to FIG. 8, a laminate 45 may comprise a plurality of layers of fibers 52-56. The fibers may be oriented at 0° and 90° as shown in FIG. 8, or in other orientations as required for a particular application. A delamination crack 60 may include a first portion 60A with layers 52-54 disposed above the first crack portion 60A, and layers 55-56 disposed below the first crack portion 60A, forming a first laminate portion 62 having a thickness t.sub.1, and a second laminate portion 64 having a thickness t.sub.2. The delamination 60 may migrate via a matrix crack 66 from the first portion or interface 60A to a second delamination crack portion or interface 60B. Upper laminate portion 68 (layers 52 and 53) has a thickness t.sub.3 above the new crack portion or interface 60B, and a lower portion 70 (layers 54-56) below the new crack or interface 60B has a thickness t.sub.4.

[0045] As shown by the arrow A, nodes FN2 and RN2 of finite element model 50 correspond to the matrix crack 66 location. The matrix crack 66 is represented by a discontinuity corresponding to the integration of the stiffness matrix (equation 3.0, supra) across different thicknesses domains as follows:

K.sub.Ω.sub.A.sup.(e)=∫∫.sub.A.sub.(e)H.sub.Ω.sub.A.sup.T[∫.sub.−h/Z.sup.z′C(A,B,D,G).sub.Ω.sub.Adz]H.sub.Ω.sub.AdA   (11.0)

K.sub.Ω.sub.B.sup.(e)=∫∫.sub.A.sub.(e)H.sub.Ω.sub.B.sup.T[∫.sub.z′.sup.−h/ZC(A,B,D,G).sub.Ω.sub.Bdz]H.sub.Ω.sub.BdA   (12.0)

where z′ is the location of a delamination along the z-axis in a laminate.

[0046] Thus, the element 24B1 (thickness t.sub.1) will have a stiffness corresponding to the layers 52-54. The element 24A1 (thickness t.sub.2) has a stiffness matrix corresponding to the layers 55 and 56. However, the element 24B2 (formed after crack 60 migrates to the new interface 60B) has a thickness t.sub.3 corresponding to layers 52 and 53. Element 24A2 has a thickness t.sub.4 corresponding to layers 54-56.

[0047] With further reference to FIG. 9, in the illustrated example delamination 60 propagates via a matrix crack from a first interface 76 (crack 60A) between layers 54 and 55 to a second interface 78 (crack 60B) between layers 53 and 54. However, the matrix crack 66 may also propagate in the opposite direction from interface 76 to a third interface 80 between layers 55 and 56. The direction of the propagation of matrix crack 66 may be determined utilizing the direction (sign, positive or negative) of the tie shear forces F.sub.x at the crack tip or delamination front 30. The shear force sign indicates the orientation of microcracks, or cusps, that precede the delamination 60. The tendency of propagation in the Z coordinate (FIG. 9) transverse to the present interface, i.e., up or down, may be determined by determining if the fiber orientation in the adjacent ply in the direction the cusps are pointing towards (i.e., the bounding fiber orientation) is such that the cusps are arrested or allowed to pass through as described in more detail in De Carvalho, N V, Chen, B Y, Pinho, S T, Ratcliffe, J G, Baiz, P M, Tay, T E. 2015, “Modeling delamination migration in cross-ply tape laminates,” Composite: Part A 71:192-203.

[0048] With further reference to FIG. 10, another aspect of the present invention involves determining the cusp orientation about the Z axis as discussed in Canturri, C., E. S., Greenhalgh, S. T. Pinho, J. Ankersen, 2013, “Delamination growth directionality and the subsequent migration process: The key to damage tolerant design, Composites Part A: Applied Science and Manufacturing,” (54):79-87 relative to the bounding ply fibers. The crack 60 may have a three dimensional crack tip or delamination front 30 with a physical representation shown in dashed lines in FIG. 10 superimposed over the staggered approximation made in the mesh. At a given node 84 located on delamination front 30, the angle of the force vector F.sub.n is calculated from the individual tie shear forces F.sub.x and F.sub.y acting on node 84. The cusp orientation about the Z-axis may be assumed to be orthogonal to F.sub.n. The relative angle, α, between F.sub.n and the bounding fiber orientation may be used as the basis for a criterion determining if cusps are allowed to pass through and result in a delamination migration. If the conditions for described above are met, migration still may not occur if the cusps are greatly misaligned with the direction of delamination growth and/or the tie shear forces are small (i.e., possible in a Mode I dominated delamination).

[0049] With further reference to FIG. 11, energy criteria can be utilized to determine migration and/or delamination at a delamination front 30. At delamination front 30, the total energy release rate (G.sub.T) is calculated. The critical energy release rate (G.sub.C) is calculated for the first interface 76. The total energy release rate and critical energy release rate (designated G′.sub.T and G′.sub.C in FIG. 11 at interface 78) are assumed to be very similar (or identical) to those of the current interface. A critical energy release rate is given by the equation 10.0 (supra).

[0050] G.sub.IC is the Mode I critical energy release rate, wherein Mode I is a crack that is “opening.” G.sub.IIC is the Mode II critical energy release rate, wherein Mode II is a shear or “sliding” type crack. G.sub.II is the Mode TI energy release rate. The following assumptions are utilized:

G.sub.I.sub.C.sup.(fiber)=∞

ΔU.sub.mig<<ΔU.sub.delam

Were ΔU.sub.mig is the energy dissipated due a matrix crack (i.e. migration), and ΔU.sub.delam is the energy dissipation due to growth of the delamination crack 60.
The energy release rates and critical energy release rates can be utilized to predict one of three possibilities at a delamination front node 86. In the following, G.sub.Ic refers to matrix crack (i.e., cusp) critical energy release rate and G.sub.c refers to delamination critical energy release rate. G.sub.T is compared to both toughness quantities in a manner similar to De Carvalho, N V, Chen, B Y, Pinho, S T, Ratcliffe, J G, Baiz, P M, Tay, T E. 2015, “Modeling delamination migration in cross-ply tape laminates,” Composite: Part A 71:192-203. Specifically, migration and delamination occur if:

G.sub.T>G.sub.Ic

&

G.sub.T>G.sub.c

[0051] Delamination occurs if:

G.sub.T<G.sub.Ic

&

G.sub.T>G.sub.c

[0052] No growth occurs if:

G.sub.T<G.sub.c

[0053] Referring again to FIG. 1, the steps designated 16 and 18 may utilize the delamination/migration criteria discussed above in connection with FIG. 11. The cusp orientation is determined using the tie shear forces as discussed above in connection with FIGS. 10 and 11. As discussed above in connection with FIG. 8, the matrix crack is modeled utilizing a discontinuity in the mesh stiffness at opposing element boundaries. Alternatively, matrix cracks may be inserted into the mesh by using floating nodes to break nodal connectivity between adjacent elements that represent a single ply.

[0054] Referring again to FIG. 1, the splitting of elements either at a new migrated interface (designated 20 in FIG. 1) and/or splitting the next element in the direction of growth at the current interface (designated 22 in FIG. 1) utilize adaptive fidelity (i.e. splitting elements on demand as needed, but not otherwise). It will be understood that this approach permits modeling of crack growth without the need to provide a large number of parallel/overlapping elements corresponding to the individual layers of the composite material at the time the initial finite element model is created.

[0055] Initially separate delaminations in a finite element model may grow independently but at some point in a solution procedure the initially separate delaminations may reach a common nodal location in the mesh. Or, similarly, separate deliminations may be adjacent to one another during growth. The migration criteria as described previously or a similar variation may be applied in these instances to determine if the delaminations link together via a matrix crack and whether a TN is released or the tie is maintained. Furthermore, if separate delaminations exist at different interfaces, the delaminations may grow to a common location in the mesh or they may be adjacent to one another during growth. The migration criteria as described above can be used to determine how elements are split and which ties are released, if any, in the region where the two delaminations interact.

[0056] The simulation tool may optionally include a fiber failure simulation capability. One method of doing this is use of continuum damage mechanics as described in Matzenmiller, A. J. Lubliner, R. L. Taylor, 1995, “A constitutive model for anisotropic damage in fiber-composites,” Mechanics of Materials, (20)2:125-152.

[0057] All references contained herein are hereby incorporated by reference in their entirety.