Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part II: Sensitivity study on modeling details

Engineering Fracture Mechanics 73 (2006) 786–801 www.elsevier.com/locate/engfracmech

Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part II: Sensitivity study on modeling details De Xie a, Sherrill B. Biggers Jr. a

b,*

Alpha STAR Corporation, 5199 East Pacific Coast Highway, Long Beach, CA 90804, USA b Department of Mechanical Engineering, Clemson University, Clemson, SC 29634, USA Received 28 May 2004; received in revised form 12 April 2005; accepted 29 July 2005 Available online 28 December 2005

Abstract The interface element and VCCT process described in Part I of this two-part paper, developed to compute strain energy release rates of an arbitrary delamination front using non-orthogonal finite element meshes, are further investigated in this paper for robustness and ease of use in tracking delamination growth. Standard 3-D elements are used in conjunction with the interface elements. No special singularity elements are required. Stationary meshes that are independent of the shape of the delamination front can be used. Three cases having different initial delamination shapes are examined. The process is shown to be insensitive to the values used for the interfacial spring stiffness, the orientation of the interface element, or even the mesh pattern if the mesh has a reasonable degree of refinement. Therefore, the method can be used with ease and confidence in general-purpose delamination growth analysis for engineering applications.
2005 Elsevier Ltd. All rights reserved. Keywords: Delamination growth; Virtual crack closure technique; Strain energy release rate; Interface element; Finite element analysis

1. Introduction In delamination growth simulation in composite laminates, determining the location of the delamination front over time is one of the major challenges because the delamination is more likely to propagate in a non-self-similar and mixed-mode manner than in isotropic materials. One way to avoid this difficulty is to develop an approach that does not require the precise shape of the delamination front. A compliant layer approach, which is similar to the Barenblatt-Dugdale (BD) model in fracture mechanics, has been proposed to make use of this idea. The compliant layer describes the delamination growth as progressive damage behavior through an updated constitutive relation across the layer. The constitutive relation can be calibrated from *

Corresponding author. Fax: +1 864 656 4435. E-mail address: [email protected] (S.B. Biggers Jr.).

0013-7944/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2005.07.014

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

787

experimental load–displacement curves [1], estimated by relating to the critical values of fracture energy [2,3], or derived through a thermo-mechanical theory based analysis [4]. Many published works on the topic are confined to this approach. Different authors have proposed models and named them in various ways such as cohesive fracture model [5], cohesive layer [1,6], interfacial decohesion [2,7]. However, at the center of these models is their indirect relations to fracture mechanics, or the introduction of damage parameters on debonding. Therefore, this approach can be referred to as an indirect approach or damage mechanics approach. Another way to study delamination growth is the direct approach or fracture mechanics approach. This approach directly utilizes fracture mechanics criteria in which the strain energy release rates (G’s) are computed and compared to the corresponding critical values. The virtual crack closure technique (VCCT) [8,9] is the most popular and powerful tool to approximately calculate G values. In the VCCT approximation, the nodal forces at the delamination front, the displacement openings just behind the front, and the virtual area just ahead of the front are required. In order to calculate the mixed-mode G’s in a general case, the nodal forces and the displacement openings need to be transformed into an instantaneous local coordinate system that follows the delamination front. Both the local coordinate system and the virtual area can be determined only after the normal direction for the advancing delamination front is known. For a fixed front, the normal direction can be easily determined through using the orthogonal meshes that match the front shape. But this is more difficult for a moving front since its shape and local normal direction change instantaneously during the delamination growth. One method for dealing with this case is to use an automatic adaptive re-meshing technique. In this way, the re-meshing can produce continually changing meshes that are orthogonal to the delamination front as it progresses. However, this approach is difficult to apply to problems that are highly non-linear, especially when global non-linearities accompany or interact with delamination growth, because of the large computational burden. In order to solve the problem in a more computationally efficient way, the present paper proposes an approach to determine the delamination front using two simple vectors that define the delamination front in an approximate way at a given point [10]. The normal vector and tangent vector for the local coordinate system can then be obtained based on the two vectors that define the front. The approach avoids the use of automatic re-meshing techniques. An interfacial element has been developed to implement the approach into the commercial software ABAQUS using the user element subroutine (UEL). Validation of the approach was carried out in Part I of this two-part paper [10] by comparing computed results with the analytical solutions or the FEA calculations from the other sources. However, issues relating to the effects of modeling details of the approach on the results have not been addressed. Therefore, in this paper, the effects of mesh density, mesh pattern, and the orientation of the interfacial element on the ultimate compressive load of delaminated plates are studied. 2. Crack front tracing approach An algorithm to determine the orientation from the instantaneous shape of the delamination front has been proposed [10]. As discussed in Part I of this two-part paper, the delamination front centered on node N (center node) can be located by defining the two of the eight vectors, originating at node N (Fig. 1(a)) and pointing

R8

R7

R6

N

R1

R5

m8 = 0 n b + n e

m e –1 = m1 = 0

me =

R2 (a)

R3

R4

ne

m2 = 0 (b)

m7 = 0

ˆj

m6 = 0

nb mb = N

tˆ

mb +1 =

kˆ

m5 = 0

Rb

Re m e +1 =

mb– 1 =

m3 = 1

m4 = 1

For this example, b=5, e=2

Fig. 1. Two vectors that define a delamination front.

ˆi

nˆ tˆ

kˆ

788

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

toward the surrounding nodes, that divide the bonded nodes from the debonded nodes according to the damage index mi (mi = 0 for perfectly bonded nodes, mi = 1 for separated nodes). Let these two vectors be defined ~b and R ~e where the ‘‘b’’ and ‘‘e’’ subscripts stand for the ‘‘beginning’’ and ‘‘end’’ of the bonded region as R ~b to R ~e in sweeping out the bonded region, as shown in where one moves counterclockwise from R Fig. 1(b). Once these crack front vectors are determined, the local coordinate system at node N can be set to determine the normal and tangent directions at the crack front and these can then be used within the VCCT analysis. After the component G’s are calculated, a mixed-mode fracture criterion is used to predict delamination growth.  Ed ¼

GI GIC

a

 þ

GII GIIC

b

 þ

GIII GIIIC

c P1

where Ed is the delamination growth parameter. GIC, GIIC, and GIIIC are the critical values corresponding to mode I, mode II and mode III loading, respectively, and are assumed to be constant during delamination growth in this study. In this study, the exponents a, b, and c were taken as 1 as suggested by Kutlu and Chang [11]. The spring stiffness, K, was reduced to zero immediately after fracture occurs (i.e., Ed P 1). The calculation of energy release rates and the application of the failure criterion were carried out using a user-defined interfacial element through the ABAQUS user element subroutine UEL. ‘‘Dummy nodes’’ are introduced to extract necessary information including the damage index, nodal coordinates, nodal forces and nodal displacements for those nodes surrounding the center node N. This approach has also applied to the 2D crack growth problems based on VCCT [12]. 3. Case studies To study the sensitivity of the approach to the mesh pattern, mesh density, and the interface element orientation, three cases were examined. Their geometry and dimensions are shown in Fig. 2 and Table 1. The delaminations are embedded in flat laminated plates and have the initial shapes shown in gray in Fig. 2. The three delaminated plates have the same stacking sequence, [0/90/90/D/0], with ‘‘D’’ referring to ‘‘Delamination’’. Therefore, the delamination is located at the interface between the top 0 and 90 plies. The mechanical properties are given in Table 2. The compressive load is applied horizontally (1-direction) on the edge of the plate as a uniform displacement. To initiate the out-of-plane deformation, an imperfection of 2% the plate thickness is prescribed in the center of the plate. The imperfection is imposed upward for the sublaminate above the delamination and downward for the sublaminate below the delamination. This imperfection does not significantly affect the post-buckling behavior of the sublaminates. Geometrically non-linear analysis (NLGEOM) was carried out to include large deformations in the post-buckling stage. The method can just as easily be applied to laminated plates with arbitrary fiber orientations, such as a quasi-isotropic [0/90/ +45/
45]s laminate and the interface element can be located at any interface regardless of fiber orientations. Of course in such cases, symmetry along the central length and width cannot be taken advantage of in the modeling. Each ply of the laminated plate is modeled using ABAQUS brick elements (C3D8) across the 1–2 plane and one element through the thickness. No singular (as used in Ref. [13]) or special elements were used aside from the interface elements and gap elements placed uniformly along the plane of potential delamination growth. Further through-the-thickness refinement would probably be warranted in most cases although the simpler model serves the purpose of this study. Use of singular elements would also improve accuracy at the expense of increasing model complexity. To employ the VCCT, the interface elements are placed between the sublaminates to hold the nodes together on opposite sides in the interlaminar region. To support the general loading, three springs are placed at each node to act in mutually perpendicular directions corresponding to the three global coordinates. The stiffness of the springs in the three directions is chosen to be the same value as ‘‘k’’. Due to the large out-of-plane deformation in the post-buckling stage, contact phenomena may be relevant in the delaminated regions. In order to prevent penetrations between the two contacting surfaces, gap elements (GAPUNI) were used to only allow separation of the node pairs in the out-of-plane direction.

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

789

Fig. 2. Geometrical configurations of delaminations. (a) Through-the-width delamination; (b) octagonal embedded delamination, (c) square embedded delamination.

Table 1 Geometry and dimensions

Case (1) Case (2) Case (3)

t/L

a/L

b/L

0.02 0.0333 0.0333

0.4 0.6 0.8

N/A 0.2 N/A

Table 2 Mechanical properties (carbon/epoxy) Stiffness, Poisson’s ratio E11 (GPa) E22, E33 (GPa) G12, G13 (GPa) G23 (GPa) m12, m13 m23

Critical values for G 109.34 8.82 4.32 3.20 0.342 0.520

GIC (N/mm) GIIC (N/mm) GIIIC (N/mm)

0.306 0.632 0.817

3.1. Through-the-width delamination The first example involves a through-the-width delamination that is investigated to gain some basic understanding of mesh sensitivity in a simple one-dimensional growth case. The geometrical configuration is shown

790

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

in Fig. 2(a). Due to the symmetry of the geometry and loading conditions, only one half of the structure is modeled. To perform the convergence study, four meshes with varying density were used and are shown in Fig. 3. There is only one element across the width so that the delamination is forced to remain straight when it propagates along the 1-direction (local direction). For this specific mesh pattern, the enforced front approach and the proposed front tracing approach using the adaptive normal vector direction predict identical results. Therefore, only one end load (Nx) vs. nominal strain (U/L) is plotted for each mesh density in Fig. 4. Note that the four plots corresponding to four mesh refinements have the same initial slope prior to the delamination buckling point. The two least refined models provide a thin delaminated layer that is overly stiff in bending and thus neither predicts well the buckling that drives the delamination growth in this example. To help see the convergence clearly, the ultimate loads, (Nx)u, taken from the four plots in Fig. 4 are plotted against the number of elements along the length of the model in Fig. 5. It can be seen that the ultimate load decreases rapidly as the meshes are refined and tends to approach a constant after the number of elements reaches a certain degree of refinement, Ne. This indicates that the solution has converged quickly. It is clear in this first example that mesh refinement has no effect on the shape of the delamination front and, therefore, the question of convergence is not as complex as in the following examples.

Fig. 3. Four refined meshes for the through-the-width delamination.

0.7 A05 Nx

Nx, (klb/in)

0.5

A10 0.3

A20 0.1 A40 0.000

0.002

0.004

0.006 U/L

0.008

0.010

Fig. 4. Average running load Nx vs. nominal strain U/L.

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

791

0.7

(Nx)u, (klb/in)

Nx

0.5

0.3

0.1 0

10

30 20 Number of Elements (Ne )

40

Fig. 5. Mesh convergence study.

Another concern is the importance of choosing a proper stiffness value for the spring of the interface element. In order to study this, a series of calculations were performed with varying spring stiffness. The ultimate load (Nx)u is plotted against the spring stiffness (k) normalized with respect to the longitudinal stiffness (E11) in Fig. 6. For each mesh density, there is a wide range for stiffness values that produces the same ultimate load. The calculations will encounter convergence difficulties if the stiffness falls out of the valid range. The finer the mesh, the broader the valid range will be. As expected, the stiffness used has almost no influence on the ultimate load for a wide range of numerical values. In other words, the approach is not sensitive to the interface element stiffness as long as a reasonable value is chosen. For this example, log(k/E11) = 6.0 or k = 109.34 · 106 GPa, which is valid for all of the mesh densities, would be one of many such reasonable values. 3.2. Octagonal embedded delamination Next, consider an embedded delamination having an octagonal shape shown in Fig. 2(b). Due to the symmetry, only a quarter of the plate has been modeled. Two mesh patterns and two interface element orientations were adopted. Investigations were made for combinations of the mesh patterns and the interface element orientations to study the approach’s sensitivities in this regard. Mesh pattern A (MA) and mesh pattern B (MB) are shown in Fig. 7(a) and (b), respectively. Each mesh pattern has four meshes of varying densities.

0.7

(Nx)u, (klb/in)

A05 0.5

A10

0.3

A20 0.1

A40 -2

0

2

4 6 Log(k/E11 )

Fig. 6. Effect of spring stiffness.

8

10

792

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

Fig. 7. Mesh patterns and interface element orientations for octagonal embedded delamination. (a) Mesh pattern A (MA); (b) mesh pattern B (MB).

The orientation of the interface element is also a question that needs to be addressed. Here two sets of orientations are studied to determine how this affects the results. Mesh pattern A with orientation A (OA) for the interface element serves as a baseline model. Four plots of average running load (Nx) vs. nominal strain (U/L) are shown in Fig. 8 for four different mesh densities. Some important features can be observed from the figures. First, the greater the number of elements, the larger the range of loading over which delamination growth can be predicted and followed. For a coarse mesh such as A25, the initial load (Nx)i at which delamination starts to grow, coincides with the ultimate load (Nx)u. As a consequence, there is no delamination growth observed before the ultimate load. For a refined mesh such as A1600, significant delamination growth occurs between the delamination initiation load and the ultimate load. Second, the method converges with respect to mesh refinement for both delamination growth initiation load (Nx)i and the ultimate load (Nx)u. These loads decrease rapidly as the meshes becomes finer but tend to approach constants after the number of elements reaches a certain degree of refinement. To study the sensitivity of the interface element orientation on approximation of the delamination front normal direction, two orientation sets (OA and OB) of the interface element were applied to mesh pattern B. The plots of average running load (Nx) vs. nominal strain (U/L) curves obtained from the present approach are shown in Fig. 9. Again, the finer the mesh, the larger the load range over which the delamination growth can be traced. In addition, it is clear that the approach predicts two nearly identical curves with the two different orientations of the interface elements (OA and OB). This is encouraging since in a general case, the element cannot be oriented in a preferred direction for all locations as delamination growth occurs. Since the delamination starts to grow at the root of the top horizontal edge where the interface elements have the same orientation for OA and OB, the delamination growth initiation loads are hardly affected by the orientation. The same can be said for the ultimate loads, which for model B3200, show a relative difference of less than 2%.

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

1.0

1.0 0.8

Nx

Nx, (klb/in)

Nx, (klb/in)

0.8 0.6

A25 OA

0.4

(Nx)i and (Nx)u

0.2 0.0 0.000

0.6

Nx

(Nx )u A100 OA

0.4

(Nx)i

0.2

0.002

0.004

0.006

0.008

0.0 0.000

0.010

0.002

0.004

U/L

0.8

Nx

Nx, (klb/in)

(Nx )u Nx, (klb/in)

0.008

0.010

0.008

0.010

1.0

0.8

A400 OA

0.4 0.2

0.0 0.000

0.006 U/L

1.0

0.6

793

0.004

0.006 U/L

(Nx)u A1600 OA

0.4 0.2

(Nx)i

0.002

0.6

Nx

(Nx)i 0.008

0.010

0.0 0.000

0.002

0.004

0.006 U/L

Fig. 8. Average running load Nx vs. nominal strain U/L (MA with OA), front tracing approach, octagonal embedded delamination.

There is some difference due to the fact that orientation OB uses interface element orientations that are ‘‘preferential’’ (closer to normal) to the portion of the initial delamination front that is sloped and slightly preferential to the front as is progresses toward the ultimate load configuration (shown latter in Fig. 12). However, this excellent agreement shows that the approach is not sensitive to the orientation of the interface element and the approach works well in detecting the delamination front as it progresses regardless of its orientation relative to the orientation of the interface elements. Comparing Figs. 8 and 9, the approach’s sensitivity to the mesh pattern is clearly seen as a significant difference between the ultimate end load vs. nominal strain curves of two coarse meshes (A25 and B50). However, this is only because the models are not yet converged at these mesh densities. In this problem, the lack of refinement affects the ability of the coarse interface element to accurately predict the G values and the shape of the delamination front as it progresses, and also the ability of the coarse 3D elements forming the delaminated layer to accurately predict the post-buckling behavior of the thin layer in compression. As the mesh is refined to A1600 and B3200, all of these deficiencies diminish and the results become very close for both mesh patterns. The relative difference between the two delamination growth onset loads for B3200 and A1600 is about 3%. For the ultimate loads and comparing to A1600, the values for B3200 with OA and B3200 with OB have relative differences of about 3% and 1%, respectively. The predictions from the two mesh patterns eventually converge to approximately the same value. Thus, the conclusion can be drawn that the present approach is not sensitive to either the mesh pattern or the interface element orientation provided the mesh is reasonably fine in the same sense that any FE mesh must have been determined to have sufficient refinement before confidence can be placed in the results.

794

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

1.0

1.0

0.6

0.4

0.8

(N x) i and (N x) u

Nx

Nx, (klb/in)

Nx, (klb/in)

0.8

B50 OA OB

0.4

0.002

0.004

0.006

0.008

0.0 0.000

0.010

(N x) u

B200 OA OB

0.2

0.2 0.0 0.000

0.6

Nx

(Nx)i

0.002

0.004

0.8

Nx

Nx, (klb/in)

Nx, (klb/in)

0.8

0.4

(Nx) u

B800 OA OB

0.010

0.6

0.4

0.008

0.010

Nx

(Nx)i

(Nx) i 0.002

0.006

0.004 U/L

0.008

0.010

0.0 0.000

(Nx) u

B3200 OA OB

0.2

0.2 0.0 0.000

0.008

1.0

1.0

0.6

0.006

U/L

U/L

0.002

0.004

0.006 U/L

Fig. 9. Average running load Nx vs. nominal strain U/L (MB with OA and OB), front tracing approach, octagonal embedded delamination.

It is now desired to determine how the results would change if the normal vectors are forced to align with the mesh pattern established according to the delamination shape prior to loading. This approach is herein referred to as the ‘‘enforced front’’ approach. By comparing the dashed (OB) and solid (OA) lines in Fig. 10, a significant difference can be seen between the results for the two different element orientation schemes. Furthermore, this difference is not eliminated through mesh refinement. The enforced front algorithm obviously depends on the orientation of the interface elements. Therefore, assuming normal vector directions prior to the analysis may prevent the correct result from being obtained even with very refined meshes depending on what those directions are and how the growth progresses. To summarize the above discussion, the ultimate loads (Nx)u are plotted versus the number of elements (Ne) in Fig. 11. One can see that all three models (MA with OA, MB with OA, and MB with OB) predict very similar results through using the front tracing approach. If the enforced front element directions are selected in an appropriate way (such as MB with OB), a good solution may be obtained. If the enforced front element directions are not selected correctly (such as MB with OA), the enforced front approach may predict a poor result. Of course, if the enforced front approach is used, one cannot be certain that orientations that are the preferred ones prior to delamination growth remain the preferred ones after growth. To help understand this, the delamination shapes at U/L = 0.004 for different combinations of mesh patterns and interface element orientations are shown in Fig. 12. Three combinations using the front tracing approach and the enforced front approach with MB with OB generate a similar shape for the delamination front at this instant, see Fig. 12(a)–(c) and (e). However, MB with OA using the enforced front approach predicts a totally different delamination shape as can be seen in Fig. 12(d). Based on the above studies, the conclusion can be drawn that with the proposed

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

1.0

1.0

(N x) i

Nx

0.6

0.8

(N x) u

B50 OA OB

0.4

Nx, (klb/in)

Nx, (klb/in)

0.8

0.6

0.4

Nx

0.0 0.000

0.002

0.004

0.006 U/L

0.008

0.0 0.000

0.010

(Nx ) u

B200 OA OB (N x) i

0.2

0.2

0.002

0.004

0.006

0.008

0.010

U/L 1.0

1.0 0.8

0.8

Nx

0.6

B800 OA OB

0.4

Nx, (klb/in)

Nx, (klb/in)

795

(Nx)u

0.2

0.002

0.004

0.4

B3200 OA OB

(N x) u

0.2

(Nx)i

0.0 0.000

0.6

Nx

0.006

0.008

0.0 0.000

0.010

(N x) i

0.002

0.004

0.006

0.008

0.010

U/L

U/L

Fig. 10. Average running load Nx vs. nominal strain U/L (MB with OA and OB), enforced front approach, octagonal embedded delamination.

1.0

1.0

0.8

Nx

(Nx) u,(klb/in)

(Nx) u,(klb/in)

0.8 0.6 0.4 Front Tracing

0.2 0.0 0

0.5

1.0

1.5

0.6

0.4

Enforced Front

0.2

MB, OA MB, OB

MB, OA MB, OB

MA, OA

Nx

2.0

2.5

3.0

3.5

0.0 1.0

1.5

2.0

Number of Elements (Ne ) X103

2.5

3.0

3.5

4.0

Log(Ne )

Fig. 11. Summary for octagonal embedded delamination.

front tracing approach all mesh patterns and interface element orientations, even the very simple mesh pattern A with ‘‘OA’’ oriented interface element, achieve very similar results.

796

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

Fig. 12. Front shape for octagonal embedded delamination at U/L = 0.004.

One might ask why the ultimate load curves in Fig. 11 first decrease sharply then very gradually increase as the number of elements increases. The increasing portion of this behavior is not evident in the simpler example of the through-the-width delamination (see Fig. 5). The reason is that there is a rather complex combination of

Fig. 13. Mesh patterns and interface element orientations for square embedded delamination. (a) Mesh pattern A (MA); (b) mesh pattern B (MB).

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

797

effects taking place in the current example as the model is refined and a simpler set of effects exists in the earlier example. The driving force for delamination growth in both examples is the post-buckling response of the delaminated layer. In the earlier example, the delamination shape is not affected by model refinement since it was forced to remain straight. However, the tendency for the thin layer to buckle and therefore drive delamination growth increases as excessive bending stiffness is removed by further mesh refinement. The reduced thin layer bending stiffness, the tendency toward increased post-buckling response, and the associated tendency toward increased delamination size (not shown in figures for this case) all produce a reduced ultimate load as the model is refined. In the current example, the delamination shape is affected by model refinement, becoming more smoothly curved and coming closer to its true shape with increased refinement. As with the earlier example, here the increased refinement brings reduced thin layer bending stiffness, increased tendency toward buckling and post-buckling, and increased delamination size. All of these tend to reduce the ultimate load. However, as the delamination shape grows smoother, the loads are redistributed less abruptly between the thin and thick layers and, perhaps more importantly, around the delamination in the plane of the laminate. The effect is much like reducing stress concentrations by using larger radii in cutouts or fillets. This last effect serves to increase the ultimate load. It is the combination of these conflicting trends that is observed in Fig. 11. Presumably with even further increases in refinement these conflicting trends would continue to reduce in magnitude and an asymptote would eventually be reached. However, for practical purposes, the conservative lower bound of the curves shown can be used in design situations. Indeed, the changes in ultimate loads past the lower bound appear to be quite small as large increases in refinement are made. 1.0

1.0

0.6

0.8

Nx

Nx, (klb/in)

Nx, (klb/in)

0.8

A25 OA (Nx)i and (Nx)u

0.4

(N x) i and (N x)u A100 OA

0.4 0.2

0.2 0.0 0.000

0.6

Nx

0.002

0.004

0.006

0.008

0.0 0.000

0.010

0.002

U/L

Nx

0.010

A400 OA

(N x) u

0.4 0.2

0.6

0.008

0.010

(Nx) u A1600 OA

0.4 0.2

(Nx ) i

0.002

Nx

0.8

Nx, (klb/in)

Nx, (klb/in)

0.8

0.0 0.000

0.008

1.0

1.0

0.6

0.006 0.004 U/L

0.06

0.004 U/L

0.008

0.010

0.0 0.000

(N x) i 0.002

0.006

0.004 U/L

Fig. 14. Average running load Nx vs. nominal strain U/L (MA with OA), front tracing approach, square embedded delamination.

798

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

3.3. Square embedded delamination Finally, the same analysis procedure was also used to study a square shaped delamination oriented on the diagonal and having the geometry shown in Fig. 2(c). The definitions of the mesh patterns (MA and MB) are unchanged and the interface element orientations (OA and OB) are quite similar to those used for the octagonal embedded delamination discussed in the previous section. They are shown in Fig. 13. Again, ‘‘OA’’ oriented interface elements were imposed on mesh pattern A to provide a baseline case and only the front tracing approach was used on this baseline. Plots of the average running load (Nx) vs. nominal strain (U/L) are shown in Fig. 14. Features similar to those observed in Fig. 8 can be seen here. Plots of the average running load (Nx) vs. nominal strain (U/L) are shown in Fig. 15 for mesh pattern B (MB) obtained from the present front tracing approach. It is again observed that the finer meshes allow the delamination growth to be traced over a larger range of loading. For all the mesh refinements, the approach again predicts very similar curves for MB with OA and MB with OB. In this case, the interface elements have different orientations for OA and OB at the place where the delamination starts to grow. This sharp corner in the initial delamination at the growth initiation location presents an extreme test of the method and leads to slightly different delamination growth initiation loads even for mesh B3200. Delamination fronts with sharp corners were also noted as being particularly difficult in Ref. [13]. The relative differences for mesh B3200 with orientations OA and OB for initiation loads and ultimate loads are about 5% and 3%, respectively. These values are a bit higher compared to the previous octagonal delamination example. The slightly larger differences are due to the added sharpness of the initial delamination at the location of growth initiation in this case. In addition, the delamination front remains more linear and does not grow to assume as smoothly curved a shape in this case as in the previous octagonal case (see Figs. 12 and 18). Thus 1.0

1.0

0.6 0.4

(Nx)i and (Nx)u

Nx

Nx

0.8

Nx, (klb/in)

Nx, (klb/in)

0.8

B50 OA OB

0.2

0.6 0.4

(N x) u

B200 OA OB

0.2 (Nx)i

0.0 0.000

0.002

0.004

0.006

0.008

0.0 0.000

0.010

0.002

0.004

U/L

0.010

B800

0.4

OA OB

Nx, (klb/in)

0.6

0.008

0.010

Nx

0.8

Nx

0.8

Nx, (klb/in)

0.008

1.0

1.0

(Nx )u

0.6

0.4

(N x) u

B3200 OA OB

0.2

0.2

(N x) i 0.0 0.000

0.006 U/L

0.002

0.004

(N x) i 0.006 U/L

0.008

0.010

0.0 0.000

0.002

0.004

0.006 U/L

Fig. 15. Average running load Nx vs. nominal strain U/L (MB with OA and OB), front tracing approach, square embedded delamination.

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

799

this example has an added preference for orientation OB compared to orientation OA from growth initiation through attainment of ultimate load. However, despite the differences in initial and final delamination shapes, these percent differences are quite small and indicate there is little sensitivity to the orientation of the interface elements. Finally, note that the effect of interface element orientation is a bit larger for growth initiation than for ultimate load. This is again due to the fact that the sharp corners in the initial delamination shape tend to becomes more rounded, albeit to different extents for the two examples, as the delamination grows. The FEbased VCCT approach will naturally handle sharp corners less well than smoother shapes of the delamination front as pointed out in Ref. [13] and as observed here. Comparison of Figs. 14 and 15 shows the insensitivity of the front tracing approach to the different mesh patterns. In particular, with reference to mesh A1600 and mesh B3200, the curves are very similar. Compared to the delamination growth initiation loads for A1600 with OA, the values for B3200 with OA and B3200 with OB have relative differences of less than 2% and 7%, respectively. Compared to the ultimate loads for A1600 with OA, the values for B3200 with OA and B3200 with OB have relative differences of less than 1% and 3%, respectively. Comments made in the previous paragraph relative to percent differences apply here as well. Good agreement exists among all three combinations (MA with OA, MB with OA, and MB with OB). This implies that the approach is not very sensitive to the mesh pattern, especially for ultimate load predictions, assuming a reasonable mesh density has been used. For this example, computations have also been made for mesh pattern B (MB) with the two different interface element orientations (OA and OB) by using the enforced front approach. The results presented in Fig. 16 show that the enforced front approach proves to be worse in this case than in the previous example. Both the delamination growth initiation loads and ultimate loads are highly dependent on the interface element 1.0

1.0

0.6

0.4

Nx

0.004 0.006 U/L

0.008

0.8

Nx

(Nx)i B800 OA OB

Nx, (klb/in)

Nx, (klb/in)

0.002

0.006

0.004

0.008

0.010

U/L

(Nx)u

0.6

0.4

Nx

(N x) i B3200 OA OB

(Nx) u

0.2

0.2 0.0 0.000

(Nx) u

1.0

0.8

0.4

0.4

(Nx )i B200 OA OB

0.0 0.000

0.010

1.0

0.6

0.6

0.2

(N x) i and (N x) u

0.002

Nx

0.8

B50 OA OB

0.2 0.0 0.000

(N x) i and (N x) u

Nx, (klb/in)

Nx, (klb/in)

0.8

0.002

0.004

0.006 U/L

0.008

0.010

0.0 0.000

0.002

0.004

0.006 U/L

0.008

0.010

Fig. 16. Average running load Nx vs. nominal strain U/L (MB with OA and OB), enforced front approach, square embedded delamination.

800

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

1.0

1.0

0.8

Nx

(Nx) u,(klb/in)

(Nx) u,(klb/in)

0.8 0.6

0.4 0.2

Front Tracing MB, OA MA, OA MB, OB

0.0 0

Enforced Front MB, OA MB, OB

0.5 1.0 1.5 2.0 2.5 3.0 Number of Elements (Ne) X103

Nx

0.6

0.4 0.2

3.5

0.0 1.0

1.5

2.0

2.5 3.0 Log(Ne )

3.5

4.0

Fig. 17. Summary for square embedded delamination.

Fig. 18. Front shape for octagonal embedded delamination at U/L = 0.0045.

orientation. It is also worth noting that for mesh density B200, the enforced front approach even predicts load–displacement curves with different shapes when interface elements are oriented differently. Therefore, assuming normal vector directions prior to the analysis is not reliable and can produce misleading results. For all the cases considered above, the ultimate loads (Nx)u are plotted versus the number of elements (Ne) in Fig. 17. The trends observed in Fig. 11 are also seen in Fig. 17 and the same comments made relative to convergence apply here. The delamination shapes at U/L = 0.0045 are shown in Fig. 18. The same observations made concerning Fig. 12 can also be made relative to Fig. 18. These results once again demonstrate the robust nature of the current front tracing approach in capturing significant aspects of the delamination front shape and location during delamination growth. 4. Conclusions This paper has shown that the proposed algorithm for tracing the front of a delamination growing in a nonself-similar way is not sensitive to the values used for the interfacial spring stiffness, the orientation of the interface element, or the mesh pattern if the mesh is reasonably fine. What constitutes an acceptable refinement must

D. Xie, S.B. Biggers Jr. / Engineering Fracture Mechanics 73 (2006) 786–801

801

be answered in the same way as in other finite element modeling efforts. To study these sensitivities, three cases having different initial delamination shapes were examined. Using a simple through-the width delamination, the spring stiffness value was shown to have almost no influence on the predictions as long as the numerical value is selected within a very wide acceptable range. For the octagonal embedded delamination, the ultimate loads obtained with different interface element orientations and the same mesh patterns, have a relative difference of less than 2%. This excellent agreement shows that the approach is not sensitive to the orientation of the interface element. When the mesh is reasonably fine, the relative difference is less than 3% for delamination growth initiation loads and for ultimate loads between the two very different mesh patterns. This suggests that the front tracing approach is also insensitive to the mesh pattern once the mesh has a reasonable density. The above observations are also true in the difficult case of a square delamination oriented on the diagonal and the approach works well in detecting the progressing delamination front. The proposed approach has been shown to be insensitive to the value chosen for the interface spring stiffness within a wide range of acceptable values, to the orientation of the interface element, and to the mesh pattern once a reasonable refinement has been achieved. This same case shows that an ‘‘enforced front’’ approach can yield very different results depending on interface element orientation even after significant mesh refinement has been achieved. Finally, it should be noted that no numerical difficulties have been encountered in the cases examined. Numerical difficulties that might have been expected include inconsistent behavior for very large values of interface element spring stiffness or for complex delamination shapes with sharp corners. All observations indicate that the approach is robust and easy to apply. The proposed front tracing method provides the capability of using relatively simple, stationary mesh patterns to simulate the moving delamination having an arbitrary shape without adapting the mesh as the delamination shape changes. The approach should prove to be valuable in future studies of delamination growth. Acknowledgements Financial support for this work was provided by NASA, the State of South Carolina, and Clemson University through the EPSCoR grant ‘‘Development and Enhancement of Research Capability for Aircraft Structures and Materials’’. References [1] La Saponara V, Muliana H, Haj-Ali HR, Kardomateas GA. Experimental and numerical analysis of delamination growth in double cantilever laminated beams. Engng Fract Mech 2002;69:687–99. [2] Camanho PP, Davila CG, Ambur DR. Numerical simulation of delamination growth in composite materials. NASA/TP-2001211041. [3] Roudolff F, Ousset Y. Comparison between two approaches for the simulation of delamination growth in a D.C.B. specimen. Aerospace Sci Technol 2002;6:123–30. [4] Chaboche JL, Girard R, Levasseur P. On the interface debonding models. Int J Damage Mech 1997;6:220–57. [5] Fleming DC. Delamination modeling of composites for improved crash analysis. J Compos Mater 2001;35:1777–92. [6] El-Sayed S, Sridharan S. Predicting and tracking interlaminar crack growth in composites using a cohesive layer model. Compos Part B: Engng 2001;32:545–53. [7] Shahwan KW, Waas AM. Non-self-similar decohesion along a finite interface of unilaterally constrained delaminations. Proc R Soc London 1997;453:515–50. [8] Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Engng Fract Mech 1977;9:931–8. [9] Shivakumar KN, Tan PW, Newman Jr JC. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Fract 1988;36:R43–50. [10] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on virtual crack closure technique, Part I: Formulation and validation. Engng Fract Mech, in press, doi:10.1016/j.engfracmech.2005.07.013. [11] Kutlu Z, Chang FK. Modeling compressive failure of laminated composites containing multiple through-the-width delaminations. J Compos Mater 1992;26:350–87. [12] Xie D, Biggers Jr. SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elem Anal Des, submitted for publication. [13] Smith SA, Raju IS. Evaluation of stress-intensity factors using general finite element models. In: Panontin TL, Sheppard SD, editors. Fatigue and Fracture Mechanics. ASTM STP 1332, vol. 29. West Canshohocken (PA): American Society of Testing and Materials; 1999. p. 176–200.