Simulation tool for damage in composite laminates

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.

BACKGROUND OF THE INVENTION

Damage processes in laminates can be inherently complex. As a result, prior methods to simulate damage in composite laminates rely on high fidelity numerical (i.e., finite element) models. A primary disadvantage of prior high fidelity simulation models is that they can be cost inhibitive in terms of required user expertise, model development time, and computational resources. Often existing tools are not feasible for use outside of research or academic type environments because, even for consideration of a simple damage process in small structures, their use requires a high level of training and time. Furthermore, if a model becomes too high fidelity and involves too many interacting damage processes, there are more sources for error to appear and magnify. There are ways to mitigate these disadvantages such as dividing models of complex structures into component level models, or locally refining fidelity around the expected area of damage. However, in doing this, additional time and expertise is required, knowledge of expected behavior beforehand is required, and the model results can be influenced by the model creator. Difficulties associated with high fidelity models could be overcome if an accurate lower-fidelity simulation tool was available.

Accordingly, there is a need for a tool to simulate damage processes in composite laminates that is rapid and less demanding in terms of user expertise, model development time and computational resources, but of equal or similar accuracy and as a higher fidelity model.

BRIEF SUMMARY OF THE INVENTION

One aspect of the present invention is a numerical simulation tool for progressive damage in laminates. The simulation tool may utilize several numerical techniques (i.e., Floating Node Method (FNM), Virtual Crack Closure Technique (VCCT), finite element analysis) in connection with a developmental damage simulation theory in a manner necessary to fully capture the formation of a three dimensional internal crack network in a laminate using a model composed of a low fidelity mesh of planar type finite elements (i.e., a shell element, a plate element, or similar finite element). Shell/plate elements are general in nature, easily usable, and computationally efficient. The model mesh remains low fidelity throughout an analysis and increases in fidelity only locally as needed to suit a damage process occurring throughout a solution procedure.

The tool can be used to simulate three dimensional laminate damage processes and is in the form of an enriched shell finite element. One component of the tool includes numerical representation in a finite element mesh of material and structural discontinuities that can change/evolve according to an ongoing damage process using a discrete type modeling approach. One embodiment of this type of modeling approach includes a technique disclosed in “The Floating Node Method (FNM),” as described in Chen, B Y, S. T. Pinho, N. V. De Carvalho, P. M. Baiz, T. E. Tay, 2014, “A floating node method for the modeling of discontinuities in composites,” Engineering Fracture Mechanics 127:104-134, which is hereby incorporated by reference in its entirety.

The tool includes a criteria to predict delamination growth. This may utilize the Virtual Crack Closure Technique (VCCT). A summary of this approach is described in Krueger, R. 2004, “Virtual crack closure technique: History, approach, and applications, Applied Mechanics Review,” (57) 2:109-143, which is hereby incorporated by reference in its entirety.

The tool may also include a criteria to predict delamination migration in the form of matrix cracks internal to a laminate. In one embodiment, no details of this damage feature are simulated other than its location and occurrence. Details regarding investigation of a delamination migration criteria are disclosed in Ratcliffe, J. G., M. W. Czabaj, T. K. Obrien, 2013, “A test for characterizing delamination migration in carbon/epoxy tape laminates,” NASA/TM-2013-218028, 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 and Greenhalgh, E. S., Rogers, C., Robinson, P. “Fractographic observations on delamination growth and the subsequent migration through the laminate,” Composites Science and Technology, 2009, 69:2345-2351, which are hereby incorporated by reference in their entirety.

The tool may also include a feature that numerically represents a geometric material discontinuity in the form of a transverse matrix crack that is adjacent to a delamination. A transverse matrix crack adjacent to a delamination may be represented in the finite element mesh as a discontinuity in mesh stiffness along an element boundary (i.e., in the case where a delamination migrates via a transverse matrix crack the laminate material on either side of the delamination changes in thickness at the migration location).

The enriched element is based on the formulation of a shear deformable shell/plate element. Suitable elements may be developed as user defined subroutines with commercial software ABAQUS, but the element could be developed using other suitable finite element software tools.

The invention may be utilized to simulate damage in composite laminate structures in various applications, including the fields of aerospace, automobile and marine industries. The invention may be useful for a structure design and certification by simulating impact damage, compression after impact damage, or any similar delamination driven damage process in a laminate. For example, the simulation tool of the present invention may be utilized to evaluate various composite laminates. A structural component for an aircraft, vehicle, space craft, ship, etc. may be designed based, at least in part, on the results provided by the simulation tool. More specifically, a composite laminate having acceptable resistance to damage (delamination matrix cracking, etc.) may be selected for a structural component of an aircraft or other vehicle using the simulation tool of the present invention, and a composite laminate having a ply orientation etc. as evaluated/selected using the simulation tool may be utilized in fabricating the component. The simulation tool of the present invention may reduce the need to fabricate and test various laminates during the design process, and also utilizes significantly less computer resources than prior multi-layer finite element methods utilized for damage simulation/evaluation of composite materials.

DETAILED DESCRIPTION OF THE INVENTION

FIG. 1depicts a process1according to one aspect the present invention. Following the start of process1, a shell finite element model is created at step2and 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.

Referring again toFIG. 1, at a next step4, a load increment is solved. The process1may then proceed directly to step14(discussed below), or the process may proceed to optional steps or features shown in the dashed region12. As discussed in more detail below, the optional process features relate to determining if a crack or flaw forms, and includes evaluating damage initiation criterion6throughout the model. This is followed by determining if a crack or flaw is formed at step8. 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 step10, and the process proceeds to step16. If no crack or flaw is formed at step8, the optional processes proceeds to step14. At14, the process involves determining if there is an existing crack or flaw. If not, the process returns to4, and the next load increment is solved. However, if it is determined at step14that there is an existing crack or flaw, the process proceeds to step16. 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 to4, 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 step18. If the delamination crack does migrate, the next elements in the direction of the growth are split at a new interface as indicated at step20, and the process then continues with the next load increment/solution4. 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.

If a delamination does not migrate at18, the process continues as shown at22. At22, the next elements are split in the direction of growth at the current interface, and the process continues with the next load increment/solution4. It will be understood that the process includes checking if the solution is complete after20and22, and ending the process if the solution is complete. If the solution is not complete after20or22, the process returns to4as shown inFIG. 1. It will be understood that the present invention is not necessarily limited to the steps and sequences described above in connection withFIG. 1. In general, there may be additional steps not shown inFIG. 1following steps20and22for 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.

With further reference toFIG. 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:

The laminate theory constitutive material matrices are as follows:

The laminate shell element stiffness integration is as follows:

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, ΩAand ΩB, where region ΩBcorresponds to the DOF of the floating nodes. The stiffness matrix for a split element is given by:
K(e)=KΩA(e)+KΩB(e)(4.0)

With further reference toFIG. 3, an undamaged shell element24initially includes nodes N1-N4that 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.

The stiffness matrix for the undamaged shell element24is given by:

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.

With further reference toFIG. 4, the shell element24can be split into first and second overlapping shell elements or components24A and24B. Shell element24A includes real nodes RN1-RN4, and the shell element24B includes floating nodes FN1-FN4. The stiffness matrix for the split elements24A and24B ofFIG. 4is as follows:

In general, the undamaged shell element24(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 element24may be split to form the shell elements or components24A and24B (FIG. 4) to simulate a delamination between the shell elements24A and24B. The stiffnesses of the shell elements24A and24B correspond to the portions of a composite laminate disposed on either side of a delamination crack. Shell elements or components24A and24B may be considered to be separate elements or they may be considered to be separate parts or components of a single element.

With further reference toFIG. 5, a finite element model50may 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 Mechanics54.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 inFIG. 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:

G1(+x)=-12⁢Δ⁢⁢A(e)[Fz(i,j)⁡(w(i-1,j)′-w(i-1,j))+Mx(i,j)(θx(i-1,j)′-θx(i-1,j)+My(i,j)⁡(θy(i-1,j)′-θy(i-1,j))](7.0)GII(+x)=-12⁢⁢Δ⁢⁢A(e)⁡[Fx(i,j)⁡(u(i-1,j)′-u(i-1,j))](8.0)GT(+x)=GI(+x)+GII(+x)(9.0)
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(e)=b2, is the area of an element for a square regular mesh. It will be understood that other mesh shapes may be utilized.

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, Gcmay be calculated using the Benzeggagh-Kenane criterion as follows:
Gc(+x)=GIc+(GIIc−GIc)(GII(+x)/GT(+x))ηBK(10.0)

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.

With further reference toFIG. 6, the enriched shell element24can be used in a mesh to represent four different damage states. These are illustrated inFIG. 6and 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,” 12thAustralian International Aerospace Congress, Australia). Thus, the Virtual Crack Closure Technique (VCCT) may be used at a delamination crack tip or delamination front30utilizing 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 elements24as labeled inFIG. 6, preceding the delamination crack tip or front30include nodes N that utilize real nodes RN only. Split elements32immediately preceding crack tip or front30include nodes TN at each corner of the split elements32. 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.

The split elements include lower and upper element components34A and34B, respectively. The lower components34A and upper components34B adjacent to the “open” side of delamination crack tip or front30include tied nodes TN at the crack tip30. The components34A include real nodes RN along boundaries36A away from the crack tip30, and the components34B include floating nodes FN along boundaries36B away from the crack tip30. The split components38A and38B away from the crack tip30include floating nodes FN and real nodes RN that are all free/not tied. In general, an offset is applied to the elements24A,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 as34A and34B, are actually the same element with two separate regions or components.

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.

With further reference toFIG. 7, a finite element model50similar to the model50ofFIG. 6but with a staggered delamination front is shown in plan view. Delamination propagation is evaluated, as show inFIG. 7, along directions X and Y aligned with the mesh. The total energy release rate (GT) is calculated along the mesh lines or boundaries40in four directions at each tied node TN on a delamination front. If the total energy release rate (GT) in any of the four propagation directions is greater than the critical energy release rate (GC), the tie is released.

Although the propagation of a delamination crack may be determined utilizing the VCCT approach as described previously and as shown inFIGS. 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.

With further reference toFIG. 8, a laminate45may comprise a plurality of layers of fibers52-56. The fibers may be oriented at 0° and 90° as shown inFIG. 8, or in other orientations as required for a particular application. A delamination crack60may include a first portion60A with layers52-54disposed above the first crack portion60A, and layers55-56disposed below the first crack portion60A, forming a first laminate portion62having a thickness t1, and a second laminate portion64having a thickness t2. The delamination60may migrate via a matrix crack66from the first portion or interface60A to a second delamination crack portion or interface60B. Upper laminate portion68(layers52and53) has a thickness t3above the new crack portion or interface60B, and a lower portion70(layers54-56) below the new crack or interface60B has a thickness t4.

As shown by the arrow A, nodes FN2and RN2of finite element model50correspond to the matrix crack66location. The matrix crack66is represented by a discontinuity corresponding to the integration of the stiffness matrix (equation 3.0, supra) across different thicknesses domains as follows:
KΩA(e)=∫∫A(e)HΩAT[∫−h/Zz′C(A,B,D,G)ΩAdz]HΩAdA(11.0)
KΩB(e)=∫∫A(e)HΩBT[∫z′−h/ZC(A,B,D,G)ΩBdz]HΩBdA(12.0)
where z′ is the location of a delamination along the z-axis in a laminate.

Thus, the element24B1(thickness t1) will have a stiffness corresponding to the layers52-54. The element24A1(thickness t2) has a stiffness matrix corresponding to the layers55and56. However, the element24B2(formed after crack60migrates to the new interface60B) has a thickness t3corresponding to layers52and53. Element24A2has a thickness t4corresponding to layers54-56.

With further reference toFIG. 9, in the illustrated example delamination60propagates via a matrix crack from a first interface76(crack60A) between layers54and55to a second interface78(crack60B) between layers53and54. However, the matrix crack66may also propagate in the opposite direction from interface76to a third interface80between layers55and56. The direction of the propagation of matrix crack66may be determined utilizing the direction (sign, positive or negative) of the tie shear forces Fxat the crack tip or delamination front30. The shear force sign indicates the orientation of microcracks, or cusps, that precede the delamination60. 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 A71:192-203.

With further reference toFIG. 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 crack60may have a three dimensional crack tip or delamination front30with a physical representation shown in dashed lines inFIG. 10superimposed over the staggered approximation made in the mesh. At a given node84located on delamination front30, the angle of the force vector Fnis calculated from the individual tie shear forces Fxand Fyacting on node84. The cusp orientation about the Z-axis may be assumed to be orthogonal to Fn. The relative angle, α, between Fnand 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).

With further reference toFIG. 11, energy criteria can be utilized to determine migration and/or delamination at a delamination front30. At delamination front30, the total energy release rate (GT) is calculated. The critical energy release rate (GC) is calculated for the first interface76. The total energy release rate and critical energy release rate (designated G′Tand G′CinFIG. 11at interface78) 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).

GICis the Mode I critical energy release rate, wherein Mode I is a crack that is “opening.” GIICis the Mode II critical energy release rate, wherein Mode II is a shear or “sliding” type crack. GIIis the Mode TI energy release rate. The following assumptions are utilized:
GIC(fiber)=∞
ΔUmig<<ΔUdelam
Were ΔUmigis the energy dissipated due a matrix crack (i.e. migration), and ΔUdelamis the energy dissipation due to growth of the delamination crack60.
The energy release rates and critical energy release rates can be utilized to predict one of three possibilities at a delamination front node86. In the following, GIcrefers to matrix crack (i.e., cusp) critical energy release rate and Gcrefers to delamination critical energy release rate. GTis 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 A71:192-203. Specifically, migration and delamination occur if:
GT>GIc
&
GT>Gc

No growth occurs if:
GT<Gc

Referring again toFIG. 1, the steps designated16and18may utilize the delamination/migration criteria discussed above in connection withFIG. 11. The cusp orientation is determined using the tie shear forces as discussed above in connection withFIGS. 10 and 11. As discussed above in connection withFIG. 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.

Referring again toFIG. 1, the splitting of elements either at a new migrated interface (designated20inFIG. 1) and/or splitting the next element in the direction of growth at the current interface (designated22inFIG. 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.

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.

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.

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