An approach for dealing with high local stresses in finite element analyses

Composites: Part A 41 (2010) 1156–1163

Contents lists available at ScienceDirect

Composites: Part A journal homepage: www.elsevier.com/locate/compositesa

An approach for dealing with high local stresses in finite element analyses F. Hélénon *, M.R. Wisnom, S.R. Hallett, G. Allegri Advanced Composites Centre for Innovation and Science, University of Bristol, Queen’s Building, University Walk, Bristol BS8 1TR, UK

a r t i c l e

i n f o

Article history: Received 2 November 2009 Received in revised form 5 March 2010 Accepted 16 April 2010

Keywords: C. Finite element analysis (FEA) A. Laminates B. Stress concentrations B. Delamination

a b s t r a c t Highly localised through-thickness stress concentrations, higher than the strength of the material, may occur when a linear elastic finite element analysis of a composite structure is performed. Such stresses may be caused by real geometrical or material discontinuities or by artefacts in the model. The objective of this paper is to present a validated approach to determine when these high stresses will not lead to failure by delamination or matrix cracking and can therefore be ignored. Named as the High Stress Concentration (HSC) method, the approach presented in this paper is found to provide good results when applied to several finite element analyses, and is also in agreement with experimental data. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction This work addresses the interpretation of finite element analysis results for complex structural components. One of the fundamental issues in post-processing data from finite element analyses is dealing with very high localised stress concentrations that sometimes arise, see Fig. 1. Are they real, potentially leading to failure, or do they represent a modelling artefact due, for example, to simplifications introduced at the meshing stage? This problem is particularly relevant for complex geometries and material arrangements in large models. Whenever it is difficult to represent fine details at a global level, simplifying assumptions must be adopted. These simplifications may lead to highly localised stresses even for relatively coarse meshes. Examples of stress concentrations due to simplified geometrical models include those arising when local curvature radii are neglected and step changes are introduced into the mesh [1–3]. Material property induced stress concentrations may also occur when the representative scale adopted in the model is much larger than the one featuring the change in mechanical properties so that smooth transitions in the values of the elastic properties are modelled by abrupt changes. A typical case is represented by fibre orientations in the corners of composite parts [4]. Such material property induced stress concentrations also occur at the free lateral surface of laminated composites due to the well-known free edge or free corner effect [5–7]. They are very important as the through thickness stresses are strongly influenced by the stacking sequence [8] and the material property mismatch (coefficient of mutual * Corresponding author. Tel.: +44 117 331 5325; fax: +44 117 927 2771. E-mail address: [email protected] (F. Hélénon). 1359-835X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesa.2010.04.014

influence, Poisson’s ratios and coefficient of thermal expansion) [9]. If high localised stresses are larger than the material strength, they may lead to the identification of spurious failure modes. Several authors have defined rational criteria for assessing their criticality. Whitney and Nuismer [10] proposed both the point stress criterion and the average stress criterion for assessing the criticality of in-plane stress concentrations in laminated composites containing a notch. Their work has been later extended to free edge stresses by Kim and Soni [11]. Although easy to use, their approach assumed the transverse tensile stress to be predominant, the interlaminar shear stresses to be zero and a characteristic length which depends on the material, the overall stacking sequence and the ply orientation. With the prospect of accounting for the other interlaminar stress components and their interaction with the through-thickness stress, authors like Kim and Soni [12], Brewer and Lagace [13], Sun and Zhou [14], and Lecuyer and Engrand [15] have formulated different criteria derived from the Tsai-Wu criterion [16]. However all of them still require a characteristic distance for averaging the stresses that has to be determined empirically. The post-processing approach presented in this paper circumvents the use of any empirically based parameter. Denoted as the ‘‘High Stress Concentration” (HSC) method, it is inspired by linear elastic fracture mechanics theory with the scope to help stress engineers make a rational judgment on the hazard represented by large localised stress values, regardless of whether they are real or just modelling artefacts. The HSC methodology presented in this paper was derived for engineering design purposes and does not determine any of the stress intensity factor components from a given finite element

1157

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

120

Stress [MPa]

110

High localised stress concentration

Finite element stress distribution Material strength

100 90

A

B

80 70 60 0

0.2

0.4

0.6

0.8

1

Selected blocks of materials from a finite element analysis

Distance [mm] Fig. 1. Example of finite element stress distribution with high localised stresses (stress values picked-up along distance A–B).

(b) The material is assumed to have a linear elastic behaviour and the configuration allows the use of linear elastic fracture mechanics. The abstracted inter-laminar crack is therefore considered as propagating in pure mode I, in mode II or in mixed-mode loading, see Fig. 2. Mode III is neglected here, although could also be included in principle. (c) The mathematical equations involved consider an abstracted crack on the point of propagating within a single material. (d) The fracture criterion accounting for the mode-mixity is the power-law [17,18]. (e) The modes are defined using the primary stress tensor components. Indeed, there is a priori no crack physically introduced in the finite element model allowing the usual definition of the modes according to the crack surface displacements [3,19,20].

model. It consists of comparing the finite element stress field with the theoretical values that would occur due to a crack on the point of propagating at the same location. This represents the worst case scenario in terms of stress distribution. If the stress levels predicted by the finite element analysis are below those attained when the presence of the crack is assumed, then the structural component is safe and the stress concentration can be ignored. If the finite element results give stress values above those associated with a critical crack, then the stress concentration marks the onset of a real potential failure mode (failure initiation plus instantaneous crack propagation) which cannot be ignored in structural integrity assessments. The comparison between the two stress distributions is based on plotting their respective fields as functions of the distance from the tip of the assumed crack. This paper is structured as follows. First, the background to the HSC method extracted from linear elastic fracture mechanics is given. Then, the methodology for pure mode I/II and for mixed-mode is presented with validation cases. A unidirectional specimen with cut central plies and a curved beam specimen with cut plies are used respectively for pure mode II and for mixed-mode. Finally, application to a high localised stress concentration occurring at an interface between dissimilar materials is described, with stresses at the free edge of a quasi-isotropic plate used for validation. In all of the cases, the approach is shown to give good results which are in most cases conservative.

Therefore, pure mode I loading is defined locally when the inter-laminar shear stress, s, is zero and the normal through-thickness stress r is positive; pure mode II loading will be defined locally when the normal through-thickness stress, r, is zero or compressive. A stress tensor providing both r > 0 and s – 0 will give a local mixed-mode loading condition. 2.2. Stress field ahead of a crack tip According to linear elastic fracture mechanics theory and based on the local axes drawn in Fig. 2, the relevant singular stress field ahead of a crack tip within an orthotropic material is given by [3,20–22]:

2. Background theory 2.1. Assumptions

K

b

c

z

z

x

z

ð1Þ

where r is the distance from the inter-laminar crack tip, KI and KII are respectively the mode I and mode II stress intensity factors defining the crack tip singularity, and r and s are the normal stress and shear stress components in the x–z plane. In the terms RI(h, l1, l2) and RII(h, l1, l2) which are defined explicitly in [19,21–23], the angle h is given from the x-axis in the x –z plane and the dimensionless parameters l1 and l2 are the complex

(a) The material is orthotropic i.e. composed of one material only or is considered as being a block of several plies with arbitrary orientations but entirely homogenised.

a

K

I II ffi RI ðh; l1 ; l2 Þ and s ¼ pffiffiffiffiffiffiffiffiffi RII ðh; l1 ; l2 Þ r ¼ pffiffiffiffiffiffiffiffi 2p r 2pr

The HSC method aims to assess high localised stresses in linear elastic finite element analyses responsible for inter-laminar failure (delamination) or for matrix cracking if the level of detail in the model is sufficient. Before introducing the theory with which the method is defined, it is important to note some assumptions:

x

y x

Fig. 2. (a) Crack included into a layer. (b) Mode I surface displacements. (c) Mode II surface displacements.

1158

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

conjugate pair roots of the characteristic equation that may be written for plane stress as:

  1 4 1 2mxz 2 1 l þ  l þ ¼0 Ex Gxz Ez Ex

ð2Þ

It is useful to notice that the angular functions RIr ðh; l1 ; l2 Þ and II Rs ðh; l1 ; l2 Þ equal unity along the interface ahead of the crack tip where h = 0. Moreover, when the critical configuration is reached, i.e. when the crack is just on the point of propagating, the stress intensity factors are equal to their critical values KIc and KIIc that are material properties. As a consequence, Eq. (1) become:

K

K

Ic IIc ffi and sc ¼ pffiffiffiffiffiffiffiffiffi rc ¼ pffiffiffiffiffiffiffiffi 2pr 2pr

ð3Þ

where respectively rc and sc are the critical normal and shear stresses. Both equations will be used for defining the critical curves in the HSC method. 2.3. Relationship between the energy release rate and the stress intensity factor

and GII ¼ aII K 2II

ð4Þ

The coefficients aI and aII are defined according to the mode, the configuration and the material. For plane stress conditions they are:

rffiffiffiffiffiffiffiffiffiffiffiffiffirffiffiffiffiffiffiffi 1=2 a11 a33 a33 2a13 þ a55 þ 2 a11 2a11 rffiffiffiffiffiffiffi 1=2 a a33 2a13 þ a55 aII ¼ p11ffiffiffi þ a11 2a11 2

aI ¼

and ð5Þ

where the aij coefficients are the components of the compliance matrix [a] that relates strains {e} to stresses {r} via {e} = [a]{r} i.e. a11 = 1/Ex, a22 = 1/Ey, a33 = 1/Ez, a12 =  mxy/Ex, a32 =  mzy/Ez, a13 =  mxz/Ex and a55 = 1/Gxz. 2.4. Definition of the equivalent stress A pure mode I or mode II condition that enables a crack to propagate usually only occurs in particular loading conditions at the boundary of the structure. Generally, a crack is subject to mixedmode loading conditions. It is therefore necessary to account for this mode-mixity to predict the delamination propagation. Several empirical laws allow accounting for this combination [3,25]. One of them is the simple and widely-used power-law energetic mixedmode fracture criterion [17,18]:



GI GIc

n

 þ

GII GIIc

n P1

2n n 2n K 2n I þ g K II P K Ic

with g ¼

ð6Þ

in which GI and GII are the energy release rates for mode I and mode II respectively and GIc and GIIc are their corresponding critical values

ð7Þ

Substituting Eqs. (1), (3), and (5) into the stress intensity factors from Eq. (7), one obtains: 1

req ¼ ðr2n þ g n s2n Þ2n P rc with g ¼

GIc GIIc

sffiffiffiffiffi Ez Ex

ð8Þ

The above defined quantity req is the ‘‘mode I equivalent stress”. This Eq. (8) means that failure occurs once this mode I equivalent stress is equal or greater than the pure mode I critical stress defined in Eq. (3). As a consequence, to analyse a mixed-mode loading case, it is sufficient to compare this equivalent stress with the mode I critical stress only.

The HSC method defined so far works straightforwardly when assessing stress concentrations in 2D models. For 3D models it is still straightforwardly applicable provided that the finite element stresses are picked-up along a direction of orthotropy e.g. along the local x-axis (b = 0) or y-axis (b = p/2) in Fig. 3. Although not used in the rest of this paper, the developments in the current section are presented here for completeness. In the case of inter-laminar stresses along an off-axis direction with an angle b different from any direction of orthotropy, it is nec c ðr; bÞ or essary to know that angle to plot the critical curve r sc ðr; bÞ, as well as the equivalent stress distribution r eq ðbÞ for mixed-mode. Since determining b may be complicated, the approach chosen for dealing with an off-axis direction is to consider the most conservative configuration among all the possible values  c ðr; bÞ, s c ðr; bÞ and r  eq ðbÞ against b. This is done by finding the of r  c ðr; bÞ or s c ðr; bÞ and by finding the angle bM angle bm minimising r  eq ðbÞ. maximising r The lowest critical curve against b defined in Eq. (3), is obtained by minimising the critical stress intensity factor which may be linked to its corresponding fracture toughness value via Eq. (4). In the rest of the paper, the inter-laminar fracture toughness is assumed constant against b. Such a hypothesis is reasonable at this stage even though it has been reported in the literature that the inter-laminar fracture toughness may vary against b [26]. Under such conditions, the angle bm is obtained  I ðbÞ or a  II ðbÞ respectively for mode by maximising the coefficient a I or mode II against b. Since for any rotation b around the z-axis the characteristic Eq. (2) holds i.e. the material stays 2D orthotropic within the  x  z plane even though anisotropic from the  I ðbÞ and a  II ðbÞ are still similar ; zg axes, the expressions of a f x; y to Eq. (5) but with the compliance matrix components varying against b [27,28]:

11 ðbÞ ¼ a

  4 1 cos4 b 1 2mxy sin b 2 sin b cos2 b þ ¼ þ  Ex Gxy Ey Ex Ex

z=z

z

σ crack

GIc aII GIIc aI

2.5. Considerations for stresses along an off-axis direction

According to the crack closure analysis initially provided by Irwin [24], the energy release rate can be given as a function of the stress intensity factor using the following relations for mode I and mode II [3,20–22]:

GI ¼ aI K 2I

which are material properties measured experimentally [20]. Reorganising the different terms and introducing Eq. (4) into this Eq. (6), it follows:

x

crack

τ

Fig. 3. Details of the local coordinate systems at the crack tip.

β

x x

ð9aÞ

1159

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

13 ðbÞ ¼  a

33 ðbÞ ¼ a

55 ðbÞ ¼ a

mxz Ex

 ¼

mxz Ex

cos2 b þ

myz Ey

2

sin b

 ð9bÞ

1 1 ¼ Ez Ez 1 Gxz

ð9cÞ 2

¼

cos2 b sin b þ Gxz Gyz

ð9dÞ

The highest mode I equivalent stress curve defined in Eq. (8) is  eq ðbÞ against b. Since the stress tensor obtained by maximising r components are directly picked-up from the finite element mesh, the only way to define the most conservative configuration is to determine the highest value of the parameter g against b which is the only one dependent on the material properties, see Eq. (8). By assuming that the fracture toughness and through thickness Young’s modulus do not vary, the required value of g is obtained by minimising the longitudinal in-plane Young’s modulus Ex ðbÞ defined as [27]:

Ex ðbÞ ¼

cos4 b Ex

þ



1 Gxy



2mxy Ex



1 cos2

2

b sin b þ

ð10Þ

sin4 b Ey

Finally, in Eq. (8), the square of the inter-laminar shear stress is defined as being s2 ¼ s2xz þ s2yz . 3. Method for pure mode I/II loading 3.1. Principle For the case of a pure mode loading analysis, the HSC method is based on a comparison between the finite element stress distribution in the neighbourhood of the high localised stress concentration and a critical curve when the finite element throughthickness stress component rF.E. or the shear stress component sF.E. are above the strengths of the material, rmax or smax, see Fig. 4. For a pure mode I loading problem, the method compares the through thickness tensile stress distribution ahead of the localised

Case 2 : Stresses may NOT be ignored

Case 1 :

Stress

Stresses may be ignored

Strength of the material

Crit ical cur ve

Distance Fig. 4. Principle of the HSC method for a pure mode I or mode II loading problem.

stress concentration with the mode I critical stress curve rc(r) defined in Eq. (3), see Fig. 2. In the vicinity of the high localised stress concentration, failure occurs if:

rF:E: P rc when rF:E: P rmax

ksF:E: k P sc

when ksF:E: k P smax

ð12Þ

Any other case that does not meet the criterion given in Eqs. (11) or (12) leads to a finite element localised stress concentration that may not be safely ignored. In such a condition, this means that a crack may initiate and propagate unstably from where the high localised stress concentration has been observed. 3.2. Validation using a unidirectional ply specimen with cut central plies To illustrate and validate the applicability of the method in pure Mode II loading, a unidirectional glass fibre/epoxy test specimen with cut central plies across its full width in pure tension is considered, see Fig. 5 [29–31]. The design with discontinuous plies leads to an inter-laminar shear stress concentration near the cut plies. Very high values above the strength of the material are generated even at low loads, which could in fact be safely ignored since they do not lead to immediate failure. However these stresses do lead to delamination once a critical load is reached. The normal stresses near the cut are compressive except at the extremity of the cut plies. Therefore this test can be considered to be a pure mode II case. The typical material properties are given in Table 1. To simulate this specimen by the finite element method, a 2D analysis of a slice through the thickness was carried out in-plane stress. The mesh was built with Abaqus/CAE using fully integrated quadratic elements of type CPS8. Due to the presence of planes of symmetry, only a quarter-model was used. Four different meshes were built to assess the influence of refinement when applying the HSC method close to the cut where high stress gradients were expected. For each mesh, respectively 1, 2, 4 and 8 elements per ply thickness, being 0.127 mm thick, were used. In the longitudinal direction, the region near the cut was subdivided into elements 125, 62.5,

Table 1 Thermo-elastic properties of E-Glass/913 epoxy. E1 = 43.9 GPa E2 = 15.4 GPa E3 = 15.4 GPa

m23 = 0.45 m13 = 0.30 m12 = 0.30

G23 = 5.31 GPa G13 = 4.34 GPa G12 = 4.34 GPa

20 0

10

ð11Þ

In the case of a pure mode II loading problem, the method is based on comparing the distribution of the absolute value of the inter-laminar shear stresses in the vicinity of the localised stress concentration with the mode II critical stress curve sc(r) defined in Eq. (3). In the vicinity of the high localised stress concentration, failure occurs if:

50

4 continuous plies 2 cut plies 4 continuous plies

Fig. 5. Unidirectional test specimen with central cut plies.

a1 = 6.6  106 °C1 a2 = 30  106 °C1 a3 = 30  106 °C1

1160

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

31.25 and 15.63 lm long in each mesh respectively. A coarser mesh was then used further away from the cut. Using a 744 MPa gross section applied tensile stress which is below the 765 MPa stress value measured at failure, the finite element results give the four shear stress curves close to the ply discontinuity, Nos. 1–4 shown in Fig. 6. Curves 1–4 come from the coarser mesh to the finer one respectively. Each of them was obtained by picking-up the finite element stress values along the interface between the continuous and cut plies. The observed peak of stress is much higher than the inter-laminar shear strength, smax, estimated at 75 MPa [32]. Applying the HSC method using GIIc = 1.08 N/mm [32,33], the result is shown for the 744 MPa applied gross section tensile stress on the test specimen, see Fig. 6. All the finite element curves are under the critical curve close to the stress concentration where the stresses are higher than the shear strength of the material. Moreover, it was found that the degree of mesh refinement is not important since the finite element curves all follow the critical curve curvature in the region above the material strength line. Under these conditions, the level of stresses may be ignored according to the principle of the HSC method. This is consistent with the experimental results which showed that the delamination did not propagate until a 765 MPa average gross stress value, see [29].

4. Method for mixed-mode loading 4.1. Principle The cases of pure mode I or mode II loading are encountered for very simple tests only. A general situation may involve both modes

Shear stress (MPa)

500 Mesh No.1 Mesh No.2 Mesh No.3

400

Mesh No.4 Mode II critical curve Shear strength

300 200 100 0 0

0.2

0.4

0.6

0.8

1

Distance (mm) Fig. 6. Finite element curves – from mesh No. 1 (coarse) to mesh No. 4 (refined) – for a 744 MPa applied gross section stress versus the critical shear stress curve.

and both stress components rF.E. and sF.E. occuring at the stress concentration. To account for the interaction between both these stresses, the equivalent stress introduced in Eq. (8) is used. For n = 1, the mode I equivalent finite element stress may be written as:

req

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi GIc ¼ ðMaxðrF:E: ; 0ÞÞ2 þ gksF:E: k2 with g ¼ GIIc

4.2. Validation using a unidirectional curved beam specimen with cut plies The applicability of the method for a mixed-mode problem is illustrated by considering a unidirectional glass fibre/epoxy resin curved beam in four-point bending, see Fig. 7. This kind of geometry with and without cut plies across the width has been studied extensively by Wisnom and co-workers, see [30–34] for more details. This test produces theoretically a pure bending moment over the inner span. The material properties are given in Table 1. The  4 =022 Þ from the lowermost to the uppermost ply lay-up is ð06 =0 with a bar indicating the selected plies cut. Linear elastic finite element simulations were performed using a plane stress model of the half geometry. Five different meshes were built with Abaqus/CAE using fully integrated quadratic elements of type CPS8 to assess the influence of the refinement when applying the HSC method. A 58 N load per unit width on the upper roller has been applied. This value is equivalent to a 1740 Nmm bending moment per unit width which is lower than the experimental 1802 Nmm bending moment measured at failure [30,32]. To account for the expected high stress concentration, a fine mesh was used in the vicinity of the cut plies. In each of the five meshes, respectively 4, 8, 16, 32 and 160 elements per ply of size 31.25, 15.62, 7.81, 3.91 and 0.78 lm in the radial direction and approximately 27.5, 13.74, 6.87, 3.44 and 0.34 lm in the circumferential direction were used. In the region far from the cut, a coarse mesh was used to reduce the number of degrees of freedom in the whole of the numerical system to be solved. Results indicated that an extreme mesh refinement is required to get a smooth stress distribution in the region above the material

4 cut plies

R8

R10

R6 R10

60

ð13Þ

It predicts inter-laminar failure once the contribution of the finite element normal and shear stress components generate a mode I equivalent stress which is higher than or equal to the mode I critical stress of a crack on the point of propagating within the material considered. Therefore, applying the HSC methodology, the plot of the mode I equivalent finite element curve gives one of the two cases highlighted in Fig. 4.

4 cut plies

R10

sffiffiffiffiffi Ez Ex

R10

120 Fig. 7. Curved beam specimen dimensions (mm) and location of the upper and lower interfaces subjected to delamination.

1161

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

strength line. They showed two high localised stress concentration phenomena, one at the upper interface subjected to delamination and the other at lower one, the former being more important than the later. Applying the HSC method using GIc = 0.25 N/mm and GIIc = 1.08 N/mm [32,33], the result at the lower interface indicates that failure would not be expected at an applied bending moment per unit width of 1740 Nmm, see Fig. 8. The finer finite element curve from mesh No. 5 is just under the critical one, indicating that the specimen is very close to failure. The very high stresses that would still be present at slightly lower loads may therefore be ignored. It is worth noting that the use of any of the coarser meshes – Nos. 1–4 – would predict a failed specimen because their finite element curves all cross the critical curve. It was also found that the coarser the mesh the higher the level of conservativeness. This means that a conservative assessment may be obtained even with a relatively coarse mesh. 5. Method for high localised stresses between dissimilar material 5.1. Principle There are often high localised stresses at the free edge of a multi material interface for general lay-ups modelled on a ply-level basis [35] or in finite element models using blocks of dissimilar materials [4]. This is due for example to the mismatch of Poisson’s ratio and coefficient of thermal expansion. Since the considerations (c) and (d) made in Section 2.1 no longer hold because the assumed inter-laminar crack is no longer within a single material, it is necessary to reconsider the definition of the critical curve, see Eq. (3). Many authors have addressed the theoretical problem of an interface crack between two bonded dissimilar homogeneous isotropic [36–40], orthotropic and anisotropic materials [38,41–46]. In all of these developments, inspired by Williams’ pioneering work [47], they have found the presence of an oscillating singularity of the stress field of the form rie confined very near to the crack 600 Mesh No.1

Stress (MPa)

500

Mesh No.2 Mesh No.3

400

Mesh No.4

300

Mesh No.5 Critical curve

200 100

tip. In addition, they have found the stress intensity factor to be a complex number and to be determined using an arbitrary characteristic length parameter. Because of the mathematical complexity of the equations, a reasonable way of applying the HSC method to a bi-material interface for engineering purposes is to consider only a mixed-mode loading condition whose level of conservatism is the highest. This is achieved by selecting the material, matA or matB, giving the lowest mode I critical curve i.e. by plotting:

1

n

rc ¼ pffiffiffiffiffiffiffiffiffi min K matA ; K matB Ic Ic 2pr

o

ð14Þ

and by selecting the material giving the highest equivalent stress curve for mixed-mode i.e. by plotting:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðMaxðrF:E: ; 0ÞÞ2 þ gksF:E: k2 where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! GmatA EmatA GmatB EmatB Ic Ic z z ; matB g ¼ Max matA min matB min matA GIIc GIIc ðEx Þ ðEx Þ

req ¼

ð15Þ

The material strength line to be plotted is the one from the lower material transverse tensile strength between matA and matB. For an off-axis direction different from any axis of orthotropy, the same approach presented in Section 2.5 holds. 5.2. Validation using a (±45/0/90)S quasi-isotropic plate The case chosen for demonstrating the applicability of the HSC method at a bi-material interface is a (±45/0/90)S quasi-isotropic plate loaded in pure tension whose lay-up and dimensions are given in Fig. 9. It is made of HTA/922, a carbon/epoxy prepreg, whose typical material properties are given in Table 2. Its length L is such that the thickness and the width are negligible in comparison. Each ply is nominally 0.125 mm thick. Under such conditions, it is possible to model a slice of the plate assuming generalised plane strain, see Fig. 9, accounting for the half thickness because of the symmetry. Provided that the ±45° plies are modelled by using an equivalent orthotropic material, so-called here Equiv45, a plane of symmetry can also be introduced across the width. Failure occurred by delamination between the 0° and the 90° plies once the experimentally measured overall longitudinal deformation, exx, reached 0.75% [48]. Four different finite element meshes of the slice were built with Abaqus/CAE by using fully integrated elements of type C3D20 to investigate its influence when applying the proposed HSC method to assess the criticality of the high stresses at the 0°/90° bi-material interface. One element per ply was used, but the free edge was

0 0

0.02

0.04

0.06

0.08

0.1

Distance (mm)

Table 2 Thermo-elastic properties of HTA/922.

Fig. 8. Equivalent stresses distributions at the lower interface susceptible to delamination – from mesh No. 1 (coarse) to mesh No. 5 (finer) at the upper interface – versus the mode I critical curve for a 1740 N mm/mm per unit width applied bending moment.

E1 = 140.77 GPa E2 = 8.85 GPa E3 = 8.85 GPa

Plane of symmetry

m23 = 0.43 m13 = 0.28 m12 = 0.28

G23 = 3.10 GPa G13 = 4.59 GPa G12 = 4.59 GPa

Plane of symmetry

20 mm

+ 45 - 45

Equiv45 0 1 mm

z

90 0.1 mm

x Fig. 9. Illustration of the considered generalised plane strain geometric approach.

a1 = -0.5  106 °C1 a2 = 37.5  106 °C1 a3 = 37.5  106 °C1

1162

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163

180 160 140 120 100 80 60 40 20 0 -20

b

Mesh No.1 Mesh No.2 Mesh No.3 Mesh No.4

60 50

TauYZ (MPa)

SigZZ (MPa)

a

Transverse tensile strength (65 MPa)

40 30

Mesh No.1 Mesh No.2 Mesh No.3 Mesh No.4

20 10 0

0

0.2

0.4

0.6

0.8

1

0

0.2

Distance (mm)

0.4

0.6

0.8

1

Distance (mm)

Fig. 10. Stress distributions from mesh No. 1 (coarse) to mesh No. 4 (finer) at the 0°/90° interface starting from the free edge. (a) Through-thickness stress distribution (SigZZ). (b) Through thickness shear stress distribution (SigYZ).

more refined with respectively for each mesh 3, 4, 7 and 8 elements per ply because of the high stress gradients expected in this zone due to the free edge effect [5,7]. The element size across the width was for each mesh respectively 83, 40, 20 and 10 lm near to the edge. Symmetric boundary conditions were applied at the bottom and on the centre of the slice model. The nodes on the left hand side plane are constrained in the x-direction whilst MPCs were applied on the opposite right hand side plane to model the generalised plane strain condition. After analysing the model by the finite element method, it was found that the most critical stress tensor components were both the through-thickness stress rzz and the inter-laminar shear stress syz in the neighbourhood of the free edge at the 0°/90° interface for a uniform 402 MPa applied tensile stress rxx corresponding to the 0.75% strain measured at failure, see Fig. 10. To apply the HSC method, the critical curve has to be plotted. The typical fracture toughness values to be used are GIc = 0.15 N/ mm and GIIc = 0.51 N/mm for the material considered. By computing the critical stress intensity factors KIc for both the 0° and the 90° plies it is found that the 0° ply gives the lowest critical curve, see Fig. 11. By combining the picked-up stresses into equivalent stresses at each point of the selected inter-laminar line from the free edge, it is found that the material properties from the 0° ply give the highest finite element curves. This was as expected since the ratio aII/aI = Ez/Ex, see Eq. (11), is equal to 1 which is higher than the one for an orthotropic material where Ex > Ez. In contrast to the maximum stress at the free edge it was found that a converged solution is obtained with mesh refinement if one considers the intersection point between the critical curve and the material strength line. This is achieved by the level of refinement shown in mesh 4. This finite element curve is found to be 8.3% below this intersection point (distance: 0.05 mm) which can still be considered a reasonably accurate result, especially considering the engineering nature of the HSC method. It is perhaps also important to

recall that the fracture toughness values used were measured from a unidirectional HTA/922 specimen whereas here they have been applied to a bi-material interface. It has been shown [26] that fracture toughness values may indeed vary from unidirectional results at bi-material interfaces. As with previous examples shown, it can be seen that the level of conservatism increases with a less refined mesh for the assessment of the high free edge stresses at a bimaterial interface. 6. Conclusion A methodology allowing the assessment of the criticality of high localised stress concentrations in linear elastic finite element analyses has been defined. Based on linear fracture mechanics theory for orthotropic materials, it aims to help stress engineers design complex composite structural components by identifying when high inter-laminar stresses in a finite element analysis can be safely ignored. Applications to some finite element test cases have proved the approach to give very good results. Mode II was validated by using a unidirectional specimen with cut central plies loaded in tension and mixed-mode was validated by using a unidirectional curved beam specimen with cut plies under a four-point bending load. The extension to the case of a high localised stress occurring at a bi-material interface gave satisfactory results validated by using a (±45/0/90)S quasi-isotropic plate in tension. The next step is to present an extension of the HSC method to fatigue loading. Acknowledgement The authors would like to acknowledge Rolls-Royce Plc for their support of this research.

Stress (MPa)

References 200 180 160 140 120 100 80 60 40 20 0

Mesh No.1 Mesh No.2 Mesh No.3 Mesh No.4 Critical curve

0

0.05

0.1

0.15

0.2

Distance (mm) Fig. 11. Comparison between the critical curves – from mesh No. 1 (coarse) to mesh No. 4 (finer) and the mode I equivalent finite element curve for a 402 MPa applied tensile stress.

[1] Botting AD, Vizzini AJ, Lee SW. Effect of ply-drop configuration on delamination strength of tapered composite structures. AIAA J 1996;34(8):1650–6. [2] Her S-C. Stress analysis of ply drop-off in composite structures. Compos Struct 2002;57:235–44. [3] Raju IS, O’Brien TKO. Fracture mechanics concepts, stress fields, strain energy release rates, delamination initiation and growth criteria. In: Sridharan S, editor. Delamination behaviour of composites. Woodhead Publishing; 2008. p. 3–27. [4] Jansson NE, Lutz A, Wolfahrt M, Sjunnesson A. Testing and analysis of a highly loaded composite flange. In: Proceedings of ICCM-13 conference, Beijing; June 25–29 2001 [CD-ROM]. [5] Kant T, Swaminathan K. Estimation of transverse/interlaminar stresses in laminated composites – a selective review and survey of current developments. Comput Struct 2000;49(1):65–75. [6] Soutis C. Analytical and numerical modelling: interlaminar stresses and freeedge effects. In: Matthews FL et al., editors. Finite element modelling of composite materials and structures. Woodhead Publishing; 2004. p. 123–38.

F. Hélénon et al. / Composites: Part A 41 (2010) 1156–1163 [7] Mittelstedt C, Becker W. Interlaminar stress concentrations in layered structures: Part I – a selective literature survey on the free-edge effect since 1967. J Compos Mater 2004;38(12):1037–62. [8] Pagano NJ, Pipes RB. The influence of stacking sequence on laminate strength. J Compos Mater 1971;5(1):50–7. [9] Herakovich CT. On the relationship between engineering properties and delamination of composites materials. J Compos Mater 1981;15:336–48. [10] Whitney JM, Nuismer RJ. Stress fracture criteria for laminated composites containing stress concentrations. J Compos Mater 1974;8:253–65. [11] Kim RY, Soni SR. Experimental and analytical studies on the onset of delamination in laminated composites. J Compos Mater 1984;18(1):70–80. [12] Kim RY, Soni SR. Failure of composite laminates due to combined interlaminar normal and shear stresses. In: Kawat K, Umekawa S, Kobayashi E, editors. Composites’86: recent advances in Japan and the United States. Proceedings of the 3rd Japan–US conference on composites materials, Tokyo; 23–25 June 1986. p. 341–50. [13] Brewer JC, Lagace PA. Quadratic stress criterion for initiation of delamination. J Compos Mater 1988;22:1141–55. [14] Sun CT, Zhou SG. Failure of quasi-isotropic composite laminates with freeedges. J Reinf Plast Compos 1988;7(6):515–57. [15] Lecuyer F, Engrand D. A methodology for the identification of a criterion for delamination initiation. In: Proceedings of JNC-8 conference, Palaiseau; November 16–18 1992. p. 751–62. [16] Tsai SW, Wu EM. A general theory of strength for anisotropic materials. J Compos Mater 1971;5:58–80. [17] Whitcomb JD. Analysis of instability-related growth of a through width delamination. NASA Technical Memorandum No. 86301; September 1984. [18] Wu EM, Reuter RC. Crack extension in fibreglass reinforced plastics, T&AM. Report No. 275. University of Illinois; 1965. [19] Paris PC, Sih GC. Stress analysis of cracks, fracture toughness testing and its application. ASTM Special Technical Publication No. 381. Sixty-Seventh annual meeting ASTM, Chicago, USA; June 21–26 1964. p. 30–83. [20] Andersons TL. Fracture mechanics: fundamentals and application. 3rd ed.; 2005. 640 p. [21] Shi GC, Paris PC, Irwin GR. On crack in rectilinearly anisotropic bodies. Intervent J Fract 1965;1(3):189–203. [22] Tada H, Paris PC, Irwin GR. The stress analysis of cracks handbook. 3rd ed. New-York, London: ASME Press, Professional Engineering Publishing; 2000. 677 p. [23] Hoenig A. Near-tip behaviour of a crack in a plane anisotropic elastic body. Eng Fract Mech 1982;16(3):393–403. [24] Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 1957;24:361–4. [25] Reeder JR. A bilinear failure criterion for mixed-mode delamination. In: Camponeschi Jr ET, editor. Composite materials: testing and design (7th vol.), ASTM Special Technical Publication No. 1206. Philadelphia: American Society for Testing and Materials; 1993. p. 303–22. [26] Andersons J, König M. Dependence of fracture toughness of composite laminates on interface ply orientations and delamination growth direction. Compos Sci Technol 2004;64(13–14):2139–52.

1163

[27] Jones RM. Mechanics of composite materials. New-York: Hemisphere; 1999. [28] Lekhnitskii SG. Theory of elasticity of an anisotropic elastic body; 1963. [29] Cui W, Wisnom MR, Jones M. An experimental and analytical study of delamination of unidirectional specimens with cut central plies. J Reinf Plast Compos 1994;13(8):722–39. [30] Wisnom MR. Delamination in tapered unidirectional glass fibre–epoxy under static tension loading. In: AIAA 32nd structures, structural dynamics and materials conference, Baltimore (MD), USA; April 8–10 1991. p. 1162–72. [31] Wisnom MR, Jones MI. Delamination due to interaction between curvature induced interlaminar tension and stresses at terminating plies. Comput Struct 1995;32(1):615–20. [32] Kaczmarek K, Wisnom MR, Jones MI. Edge delamination in curved (O4/±456)s glass-fibre/epoxy beams loaded in bending. Compos Sci Technol 1998;58(1):155–61. [33] Petrossian Z, Wisnom MR. Prediction of delamination initiation and growth from discontinuous plies using interface elements. Compos Part A 1998;29(5):503–15. [34] Wisnom MR. 3-D finite element analysis of curved beams in bending. J Compos Mater 1996;30(11):1178–90. [35] Lessard LB, Schmidt AS, Shokrieh MM. Three-dimensional stress analysis of free-edge effects in a simple composite cross-ply laminate. Int J Solids Struct 1996;33(15):2243–59. [36] Erdogan F. Stress distribution in bonded dissimilar materials with cracks. J Appl Mech 1965;32:403–10. [37] England AH. A crack between dissimilar media. J Appl Mech 1965;32(2):400–2. [38] Gotoh M. Some problems of bonded anisotropic plates with cracks along the bond. Int J Fract 1967;3:253–65. [39] Rice JR, Sih GS. Plane problems of cracks in dissimilar media. ASME J Appl Mech 1965;32:418–23. [40] Shih CF. Cracks on bimaterial interfaces: elasticity and plasticity aspects. Mater Sci Eng 1991;A143:77–90. [41] Clements DL. A crack between dissimilar anisotropic media. Int J Eng Sci 1971;9:257–65. [42] Willis JR. Fracture mechanics of interfacial cracks. J Mech Phys Solids 1971;19:353–68. [43] Bassani JR, Qu J. Finite crack on biomaterial and bicrystal interface. J Mech Phys Solids 1971;37:435–53. [44] Ting TCT. Interface cracks in anisotropic bimaterials. J Mech Phys Solids 1990;38:505–13. [45] Suo Z. Singularities, interfaces and cracks in dissimilar anisotropic media. Proc Royal Soc London 1990;A427:331–58. [46] Qu J, Bassani JL. Interfacial fracture mechanics for anisotropic bimaterials. ASME J Appl Mech 1993;60:422–31. [47] Williams ML. The stress around a fault of crack in dissimilar media. Bull Seismolog Soc Am 1959;49(2):199–204. [48] Rezaifard A. Characterisation of damage development and residual properties in F922/HTA 5131 quasi-isotropic laminates. PREDICT II project. Sowerby Research Centre, Report No. JS13484; May 1996. 27 p.