- Email: [email protected]
Cohesive and continuum damage models applied to fracture characterization of bonded joints
ARTICLE IN PRESS
International Journal of Mechanical Sciences 48 (2006) 493–503 www.elsevier.com/locate/ijmecsci
Cohesive and continuum damage models applied to fracture characterization of bonded joints M.F.S.F. de Moura, J.A.G. Chousal Faculdade de Engenharia da Universidade do Porto, R. Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal Received 20 April 2005; received in revised form 27 November 2005; accepted 17 December 2005 Available online 7 February 2006
Abstract In this work, two different methods for simulating damage propagation are presented and applied to fracture characterization of bonded joints in pure modes I and II. The cohesive damage model is based on a special developed interface finite element including a linear softening damage process. In the continuum damage model the softening process is performed by including a characteristic length associated with a given Gauss point. The models were applied to the simulation of ‘‘double cantilever beam’’ (DCB) and ‘‘end notched flexure’’ (ENF) tests used to obtain the critical strain release rates in mode I and II of bonded joints. In mode I it was observed, under certain conditions, a good agreement between the results obtained by the two models with the reference value of critical strain energy release rate in mode I ðG Ic Þ, which is an inputted parameter. However, in mode II some discrepancies on the obtained GIIc values were observed between the two models. These inaccuracies can be explained by the simplifying assumptions inherent to the cohesive model. Better results were achieved considering the crack equivalent concept. r 2005 Elsevier Ltd. All rights reserved. Keywords: Damage models; Bonded joints; Fracture characterization
1. Introduction Adhesive joints are widely used in structural applications. However, it must be recognized that in the majority of the cases the applications are restricted to secondary non-critical structures. As a result, several different theories and models are being developed in order to increase the mechanical and structural performance of these joints and, obviously, the confidence of the designers. The joint strength prediction is usually based on the strength of materials or on the fracture mechanics models. Those based on the strength of materials generally use a stress or a strain failure criterion [1–3]. The fracture mechanics models are usually based on energetic criteria [4–10], namely on the relationship between energy release Abbreviations: FPZ, fracture process zone; DCB, double cantilever beam; ENF, end notched flexure; CBT, corrected beam theory; CCM, compliance calibration method; DBT, direct beam theory Corresponding author. Tel.: +351 225081727; fax: +351 225081584. E-mail address: [email protected] (M.F.S.F. de Moura). 0020-7403/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2005.12.008
rates with their critical values. Pradhan et al. [7,8] modelled the initiation of the debonding at the interfaces between the adherends and the adhesive. The authors considered pairs of nodes at these critical regions, which were released in function of the strain release rate. They concluded that this model can be used in the prediction of crack initiation and growth leading to failure in adhesive joints. Bogdanovich et al. [9] simulated the progressive failure analysis of double lap adhesive composite joints using the strain energy release rates to predict different possibilities of damage propagation: cohesive, adhesive or interlaminar crack propagation. The locus of initial failure and failure mode are determined by analysing the computed strains. Following, a progressive failure analysis is performed from the initial crack, which can be located at the interface adhesive–adherend, inside the adhesive or at interfaces between layers of the composite adherend. Gonc- alves et al. [10] used interface finite elements including a mixed-mode damage model to simulate progressive debonding in single lap adhesive joints. The damage model combines aspects of strength-based analysis and fracture mechanics and
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
494
Nomenclature a aeq a0 Da B C C D di D d di do;i du;i dr E E u;i F
crack length equivalent crack length initial crack length correction to the real crack length specimen width stiffness matrix for isotropic material compliance diagonal matrix containing the damage parameter damage parameter for mode i loading (i ¼ I,II) correction to account for crack tip rotation and deflection in mode I load-point displacement measured in the fracture test current relative displacement onset relative displacement of the softening process ultimate relative displacement of the softening process vector of relative displacements Young’s modulus diagonal matrix containing stiffness parameters ultimate strain on i direction correction factor accounting for large displacements for the moment arm
includes a linear softening process, which allows a smoothly decrease of stresses as the relative displacements grow. The area under the softening curve is equated to the critical energy release rate. The authors concluded that the model is accurate in predicting the joints mechanical behaviour. Li et al. [11] applied a three parameter cohesive-zone model to the study of mode-I fracture of adhesive joints made from a polymer–matrix composite. Cohesive-zone parameters were obtained by comparing numerical results to experimental observations. It was shown that there is a distinction between the characteristic strength of the interface associated with the toughness, and the intrinsic cohesive strength of the interface. Different cohesive-zone models were used for the interface and for the composite. It was concluded that cohesive-zone models accurately predict the joint strengths and deformations, and also the transition between failure of the composite and failure of the interface, which cannot be predicted by conventional fracture mechanics. The accuracy of the energetic criteria is highly dependent on the critical strain release rates of the adhesive joints. Several authors studied different specimens and methodologies to characterize the fracture mechanics of bonded joints. De et al. [12] developed an alternate specimen for measuring the fracture toughness of adhesively bonded joints. The specimen allows obtaining fracture toughness data over the entire range of mode mixity from pure mode
Fe fv G Gi G ic h I L lc n N N n P r su;i U u uL Vi
vector of external loads correction factor accounting for transverse shear effects shear modulus strain energy release rate for mode i loading critical strain energy release rate for mode i loading specimen half-thickness identity matrix length of DCB specimen and half-length of the ENF specimen characteristic length Poisson’s ratio shape function matrix correction factor accounting for large displacements for the compliance number of points at the FPZ load applied in the fracture test vector of interface stresses strength on i direction nodal displacement vector displacement vector at the integration points displacement vector at the integration points in local coordinates unit vectors of the local directions for mode i loading
I to pure mode II using a single geometry. The specimen is calibrated using the finite element method and mode mixity solutions assuming linear elastic fracture mechanics conditions. Yang et al. [13] performed a numerical study of the elastic–plastic mode-II fracture of adhesive joints. A traction-separation law was used to simulate the mode-II interfacial fracture of adhesively bonded end notch flexure (ENF) specimens with extensive plastic deformation accompanying failure. The authors found that the numerical predictions for the loads and deformation were in excellent agreement with the corresponding experimental results. The mode-II fracture was found to be much higher than mode-I and it was concluded that fracture is highly mode-dependent. A major contribution to the high value of GIIc is the higher capability of the adhesive layer to deform plastically in response to shear. Similar results were also observed by Chai [14] in his systematic studies of the fracture of elastic adhesive joints. These studies showed that, except for extremely thin adhesive layers ðo10 mmÞ, the typical value of GIIc is about an order of magnitude higher than the corresponding GIc value. The objective of this work is to evaluate the adequacy of the double cantilever beam (DCB) and ENF tests to obtain the critical strain energy release rates of bonded joints. A cohesive and a continuum damage models were developed and applied to the simulation of the experimental tests in order to understand the details of the
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
mechanical behaviour of the used specimens. Some conclusions were drawn about the satisfactoriness of these tests in fracture characterization of adhesive joints. 2. Damage models Damage can be considered as the progressive deterioration which occurs in materials before failure. It can be constituted by micro-cracking, voids nucleation and growth, and several inelastic processes that weaken the material. The analysis of cumulative damage is fundamental in life prediction of components and structures under loading. Kachanov [15] and Rabotnov [16] introduced the concept of simulating material damage by introducing damage variables into the constitutive equations. With this approach, the exact description of each individual micro-crack evolution is replaced by an averaging analysis describing the material macroscopic behaviour. In this work a cohesive and a continuum damage models are presented. In both cases, the material failure law can be defined as a relationship between normal (mode I) or shear (mode II) stresses on the crack plane and the corresponding normal and shear deformation variables (relative displacements in the cohesive model and strains in the continuum model). The objective is to evaluate the adequacy of each one in order to capture the phenomena related with the fracture characterization tests in pure modes of adhesive joints.
The cohesive damage model is based on using interface finite elements including a linear softening damage process. A user element subroutine UEL was developed to formulate a six-node interface finite element for two-dimensional problems (see Fig. 1) with the main ABAQUSs code. The stresses at the interface finite element are calculated from (1)
r ¼ Edr
being dr the vector of relative displacements between homologous points belonging to the top and bottom surfaces of the interface finite element in local coordinates. These are defined as associated to the failure modes (II, I) (2)
dr ¼ uL;top uL;bot or dr ¼
The matrix E is a diagonal matrix containing the stiffness parameters introduced by the user. They have to be carefully chosen to obtain a good performance of the model. As small values induce physically unacceptable interpenetrations and large values produce numerical problems, the optimum interface stiffnesses are the largest values that do not produce numerical problems. The interface stiffnesses were taken equal to 107 N=mm3 , which was shown to be a value satisfying the two referred conditions [17,18]. The displacement field expressed in the global system ðx; yÞ and associated to the element faces (top and bottom) can be obtained from utop;bot ¼ Ntop;bot Utop;bot ,
where U is the nodal displacement vector and N is the shape function matrix. The shape quadratic functions of the interface element are similar to those used in the eight node isoparametric plane solid element. The relation between displacements in local and global coordinates can be written as uL ¼ hT u,
dII dI
( ¼
uII uI
)
( L;top
(5)
being h ¼ ½VII ; VI ,
(6)
where VII ; VI are unit vectors of the local directions. VII can be written as T qx=qII; qy=qII; 0 VII ¼ (7) qx=qII; qy=qII; 0
VI ¼ k VII .
dr ¼ hT Ntop Utop hT Nbot Ubot
uII uI
. L;bot
I II
Fig. 1. Interface element.
(3)
(9)
or dr ¼ BUe ,
(10)
where B ¼ hT ½Ntop ; Nbot and U ¼
)
(8)
Substituting Eqs. (4) and (5) in Eq. (2) it can be established a relation between the local relative displacements in function of the nodal displacement vector
e
)
(4)
and VI is normal to the plane defined by VII and k ¼ ð0; 0; 1ÞT
2.1. Cohesive damage model
(
495
(
Utop Ubot
(11)
) .
(12)
The interface stresses can now be obtained from Eq. (1). The stiffness matrix can be obtained from the functional of the potential energy Z 1 e dT r dA UeT Fe , (13) pðU Þ ¼ 2 A r where Fe is the vector of external loads. The state of equilibrium is achieved when the displacement field Ue
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
496
minimizes the potential energy which leads to Z BT EB dA Ue , Fe ¼
(14)
A
where the integral between parentheses corresponds to the stiffness matrix which, with Fe , represent the contributions of the interface element to the global problem. The interface element incorporates a pure mode damage model to simulate damage initiation and growth. The damage model combines aspects of strength-based analysis and fracture mechanics. It is based on a softening process between stresses and interfacial relative displacements. To avoid the singularity at the crack tip and its effects, a gradual degradation is considered instead of a sudden one which would result in mesh-dependency. It is assumed that failure occurs gradually as energy is dissipated in a cohesive zone behind the crack tip. This is equivalent to the consideration of a ‘‘fracture process zone’’ (FPZ), defined as the region in which the material undergoes softening damage by different ways, e.g., micro-cracking and inelastic processes. This is achieved by including a damage parameter which moves smoothly from zero to unity as the material fails. For pure mode (i ¼ I or II) loading, the linear softening process starts when the interfacial stress reaches the respective strength su;i (see Fig. 2). The softening relationship can be written as r ¼ ðI DÞEdr ,
(15)
where I is the identity matrix and D is a diagonal matrix containing on the i-position the damage parameter which represents the damage accumulated at the interface di ¼
du;i ðdi do;i Þ , di ðdu;i do;i Þ
(16)
where di is the current relative displacement and do;i , du;i the onset and ultimate relative displacements of the softening process. This damage parameter is zero for the undamaged state and reaches one when the material is fully damaged. In pure modes, the other component of stress is not important and can be abruptly cancelled equating to one the respective element of the diagonal matrix D. Due to i u,i
B
A
0
C o,i
u,i
i
Fig. 2. Softening stresses/relative displacements relationship for the cohesive damage model.
the irreversibility of the damage process, the damage parameter d i for a material point at a given increment j is determined using the following equation: d i ¼ maxðd ji ; d j1 i Þ,
(17)
which avoids an healing phenomenon. If the load reverses during the softening process, the unloading and reloading procedures are represented by the line O–A in Fig. 2, and the released energy is given by Gi ¼ 12du;i ðsu;i si Þ.
(18)
The maximum relative displacement du;i , for which complete failure occurs, is obtained by equating the area under the softening curve (triangle OBC of Fig. 2) to the respective critical fracture energy Gic ¼ 12su;i du;i .
(19)
Due to the softening process is usual that convergence difficulties arise during the non-linear numerical analysis which is performed by using the Newton–Raphson method. In fact, although loading displacement was applied gradually considering small increments to ensure a smooth propagation process, a lot of iterations are needed for each increment. To overcome this problem the tangent stiffness matrix is considered. This matrix is defined from the tangent modulus matrix that can be obtained from the variation of Eq. (15) qr ¼ q½ðI DÞEdr þ ðI DÞE qdr
(20)
or qr ¼
dr;o E qdr , dr;u dr;o
(21)
where dr;o and dr;u represent the onset and ultimate relative displacements vectors. The tangent modulus matrix ðEtan Þ is the operator that relates qr with qdr and is used in Eq. (14) instead of E to define the tangent stiffness matrix. 2.2. Continuum damage model Finite element simulations using continuum damage models are known to be susceptible to mesh sensitivity. This means that refining the finite element discretization does not allow the numerical solution to converge to a physically meaningful solution of the problem. To overcome this problem, the continuum isotropic pure mode damage model presented here is based on a similar softening law used in the cohesive damage model and it was implemented via USDFLD user subroutine of ABAQUS Standards software. The basic assumption is to establish a correlation between the increase of damage and the decrease in stiffness at a given material point of a plane solid element. The adopted solution is based on the consideration of a damage parameter in a particular direction related to a given length associated with a Gauss point. The stiffness degradation is related to the critical strain energy release rate. However, in this case
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
i
497
tion between crack surfaces is avoided considering the initial modulus. Quadratic solid plane elements with eight nodes and four Gaussian points were considered in the simulations. The damage parameter is updated in the user subroutine during the incremental non-linear process affecting the referred elastic properties.
u,i
A
3. Models validation 0
o,i
u,i
i
Fig. 3. Softening stresses/strains relationship for the continuum damage model.
a characteristic length l c must be introduced to transform the crack opening (mode I) or shear (mode II) displacements into an equivalent strain (see Fig. 3) G ic ¼
1 2su;i u;i l c .
(22)
This parameter was considered to be equal to the length of influence of a Gauss point in the given direction and physically can be regarded as the dimension at which the material acts homogeneously. A specific algorithm was implemented to evaluate the characteristic length corresponding to each Gauss point. Basically it consists in computing the dimension of the element from the nodal coordinates and taking into account the number of Gauss points in the referred direction. Considering isotropic damage evolution, the stress– strain relation can be written considering an equation similar to Eq. (15) r ¼ ð1 d i ÞCe,
(23)
being d i the scalar damage parameter which grows from zero (undamaged initial state) to one (complete loss of integrity), and C the stiffness matrix of the undamaged material. The damage parameter is calculated by an expression similar to Eq. (16) but considering strains instead of relative displacements (see Fig. 3). According to the model assumptions the material stiffness remains isotropic during its progressive degradation and the Poisson’s ratio is not affected by damage. However, there is an implicit shear coupling. In fact, due to damage growth, the Young’s modulus decreases and, consequently, the shear modulus of the material also decreases due to the relationship G¼
E . 2ð1 þ nÞ
(24)
The softening law is similar to the one used in the cohesive model and the healing phenomenon is also avoided using the statement of Eq. (17). The properties are smoothly reduced due to the energy released at the FPZ. When contact occurs during crack propagation, the interpenetra-
In the present work a beam theory and fracture mechanics analyses are applied to obtain analytical P2d curves for the DCB and ENF specimens in order to validate the numerical models. The transverse shear deformations effects were neglected. 3.1. DCB specimen Considering the simple beam theory the relationship between the applied load and deflection at the end of the specimen, before damage starts to grow, can be written as 8a30 P, (25) EBh3 being B the width and h the half height of the specimen. During damage propagation ða4a0 Þ the Irwin–Kyes expression
d¼
P2 dC , (26) 2B da where C is the compliance ðC ¼ d=PÞ, is used to obtain the crack length a in function of G Ic 1=2 B Eh3 G Ic a¼ . (27) P 12 Gic ¼
Substituting Eq. (27) in Eq. (25) the P2d relationship during propagation is B2 EhG Ic 3=2 . (28) d¼ 3 EP2 3.2. ENF specimen Following the beam theory the displacement at the loading point, before damage initiation, can be written as 3a30 þ 2L3 P. (29) 8EBh3 During crack growth and for crack length smaller than L ðaoLÞ the displacement at the mid-span is obtained from Eq. (29) considering a instead of a0 . Combining Eqs. (29) and (26) the crack length can be obtained from 1=2 4B Eh3 G IIc a¼ , (30) 3P d¼
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
498
4.1. Double cantilever beam test
which, after substitution in Eq. (29) leads to P 32B3 3 3=2 3 d¼ L þ ðEh G IIc Þ . 9P3 4EBh3
(31)
For crack length larger than L ða4LÞ the beam theory originates d¼
P 3 3 3 ð2L aÞ L 8 EBh3
(32)
and the critical strain energy release rate is G IIc ¼
9P2 ð2L aÞ2 . 16B2 Eh3
(33)
The DCB test (see Fig. 4) is used to obtain the mode I critical strain release rate ðG Ic Þ. In the cohesive damage model, interface elements with null thickness were placed between the adherends along the specimen length to simulate the adhesive (see Fig. 5). In the continuum damage model the adhesive was modelled by two-dimensional solid elements with 0.25 mm thickness (see Fig 6). During load history the values of force, applied displacement and crack length ðP2d2aÞ were registered in order to calculate the critical strain energy release rate by using the compliance calibration method (CCM) using Eq. (26). Cubic polynomials ðC ¼ C 3 a3 þ C 2 a2 þ C 1 a þ C 0 Þ were used to fit the C ¼ f ðaÞ curves. The corrected beam theory (CBT) was also used to predict the R-curve
Substituting ð2L aÞ of Eq. (33) in Eq. (32) gives P 8B3 3 3=2 3 d¼ L 3 ðEh G IIc Þ . 9P EBh3
GIc ¼ (34)
3Pd , 2Bða þ DÞ
(35)
P,
4. Results The presented damage models were applied to the simulation of fracture characterization tests of bonded joints. One of the main differences between the two models is related with the thickness of the adhesive. In the cohesive model the adhesive thickness is neglected, i.e. the interface elements connecting the adherends have null thickness. On the contrary, the continuum model considers the real adhesive thickness (a typical value of 0.25 mm was considered). The analyses were performed considering plane strain conditions and non-linear geometrical behaviour. A loading displacement was applied progressively considering small increments to ensure a smooth propagation process. The used meshes have 2 280 solid plane elements in each arm. The specimens’ geometry and materials properties are presented in Tables 1 and 2 and Figs. 4 and 9.
t
h h a0
L Fig. 4. The DCB test.
Fig. 5. Mesh of the DCB specimen used in the cohesive damage model.
Fig. 6. Mesh of the DCB specimen used in the continuum damage model.
Table 1 Specimens geometry
30 100 50 1.5 0.25 30 10
25 20 P (N)
DCB - L (mm) ENF-L (mm) h (mm) t (mm) a0 (mm) B (mm)
15 Analytical model Cohesive damage model Continuum damage model (E=400 MPa) Continuum damage model (E=1000 MPa) Continuum damage model (E=4000 MPa) Continuum damage model (2 el., E=400 MPa)
10 Table 2 Materials properties
5
Adherends
Adhesive
0
E(MPa) n
E (MPa) n
69 000
0.33 4000
su (MPa) su (MPa) GIc (N/mm) GIIc (N/mm)
0.3 20
40
0.3
0.6
0
1
2
(mm)
3
Fig. 7. P2d curves for the DCB test.
4
5
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
where D represents a correction for crack tip rotation and deflection and is obtained from a linear regression analysis of C 1=3 versus a data [19]. Figs. 7 and 8 show the P2d and G ic ¼ f ðaÞ relations for the cohesive and continuum models. Good agreement between the analytical, cohesive and continuum ðE ¼ 400 MPaÞ damage models was observed in the post-peak region of the P2d curves, which corresponds to damage propagation. The higher peak load and stiffness value obtained by the analytical model can be explained by the softening onset that occurs in the earlier stages of the numerical P2d curves. In fact, it was observed that several points initiate the softening process a lot before the maximum load is attained, which obviously affects the stiffness. These effects are not accounted for in the analytical approach. In the post-peak region, it can also be observed some differences between the P2d curves of the cohesive and the continuum model. The agreement between the two models increases if we consider adhesives of lower modulus (E ¼ 1000 and 400 MPa) in the continuum damage model. A similar tendency can be observed for the R-curves considering different data reduction schemes. Table 3 presents the results obtained for both damage models. It can be concluded that cohesive model presents good accuracy relatively to the reference value, which is a material property introduced in the simulations (see Table 2). In the continuum model the difference relatively to the reference value increases with the Young’s modulus of the adhesive. In fact for E ¼ 4000 MPa the Young’s modulus of the adhesive becomes
499
non-negligible relatively to the flexural modulus of the adherends. It must be noted that this behaviour cannot be simulated in the cohesive model as it neglects the adhesive thickness. Similar results were obtained by Blackman et al. [20] on the analysis of the tapered DCB adhesive joint specimen. On the other hand, non-symmetrical damage propagation is simulated. In fact, the crack growth occurs at the Gauss points, which are located out of the mid-plane of the specimen. It must be noted that this non-symmetrical damage propagation occurs frequently in real tests. The consequence is the presence at the FPZ of non-negligible normal and shear stresses parallel to the crack inducing local mixed-mode loading and influencing the load–displacement relationship. Once again it must be emphasized that these effects are not captured when a cohesive damage model is used because adhesive thickness is neglected. To evaluate the influence of mesh refinement, an analysis considering two elements along the adhesive thickness direction was carried out. As observed in Fig. 7, the obtained P2d curve agrees quite well with the single element one, showing the model insensitivity to mesh refinement. 4.2. End notch flexure test The ENF test (see Fig. 9) is used to obtain the mode II critical strain release rate ðG IIc Þ. At the initial crack contact
P, a0
0.35
GIc(N/mm)
0.3
h t h
0.25 GIc (reference value) GIc (Cohesive model-CCM) GIc (Cohesive model-CBT) GIc (CCM, E=4000 MPa) GIc (CBT, E=4000 MPa) GIc (CCM, E=1000 MPa) GIc (CBT, E=1000 MPa) GIc (CCM, E=400 MPa) GIc (CBT, E=400 MPa)
0.2 0.15 0.1 0.05
L
L
Fig. 9. The ENF test.
0 30
32
34
36 a (mm)
38
40
42
Fig. 10. Detail of the ENF specimen in the cohesive damage model.
Fig. 8. The R-curves for DCB test.
Table 3 Results obtained for the DCB specimen Reference value G Ic ¼ 0:3 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
E ¼ 400 MPa
CCM
CBT
CCM
CBT
CCM
CBT
CCM
CBT
0.306 1.86 0.83
0.302 0.7 0.04
0.229 10.2 0.16
0.27 10.1 0.13
0.287 4.19 0.24
0.288 3.95 0.126
0.299 0.41 0.17
0.299 0.44 0.14
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
500
conditions were imposed to prevent interpenetration of the cracked parts. Details of the deformed shapes in both models are shown in Figs. 10 and 11. Fig. 12 presents the P2d curves obtained from the analytical and numerical approaches. In the post-peak region excellent agreement between the cohesive, continuum ðE ¼ 4000 MPaÞ and analytical models. The initial stiffness in the continuum damage model is greater than in the cohesive damage model. Considering a lower modulus adhesive ðE ¼ 1000 MPaÞ a slightly decrease on the initial stiffness was observed. The differences relatively to the cohesive and analytical models can be explained by the higher specimen thickness (plus 0.25 mm) considered in the continuum damage model. As observed in the DCB case, the mesh refinement along adhesive thickness (two elements) has a small influence on the P2d curve (see Fig. 12). Three different data reduction schemes were used to obtain the R-curves in mode II:
The CCM was also applied considering Eq. (26) and similar procedures described in the DCB test
Fig. 11. Details of the ENF specimen in the continuum damage model.
P (N)
200
G IIc ¼
9Pda2 . 2Bð2L3 þ 3a3 Þ
(36)
The CBT [21,22] 9Pa2 F f v, (37) 16B2 Eh3 N where F and N are corrections to account for large displacements for the moment arm and for the compliance, respectively, and f v comprise the effect of transverse shear G IIc ¼
From Fig. 13 and Table 4 we can conclude that in the cohesive model the best agreement was obtained by the CCM, and the worst by the CBT. However, in the continuum model the beam theories present better agreement with the reference value than the CCM (see Fig. 14 and Table 4). As in the DCB test non-symmetrical damage propagation was simulated. Consequently there exist at the FPZ non-negligible normal stresses parallel and perpendicular to the crack inducing local mixed-mode loading influencing the compliance, which explain the poor agreement of the CCM. The decrease of the resin modulus increases the G IIc value obtained from the different data reduction schemes (see Table 4). However, contrary to the DCB test the agreement does not increase. In fact, the lower modulus induces a more pronounced local deformation at the crack tip leading to more important normal
Analytical model Cohesive damage model Continuum damage model (E=4000MPa) Continuum damage model (E=1000MPa) Continuum damage model (2 el., E=1000 MPa)
0.7 0.6 GIIc (N/mm)
250
The direct beam theory (DBT) based on elementary beam theory [21]
150 100
0.5 0.4
GIIc (ref) GIIc (CBT)
0.3
GIIc (CCM)
0.2
50
GIIc (DBT)
0.1 0 0.0
1.0
2.0
3.0
4.0
0 32
34
36
(mm) Fig. 12. P2d curves for the ENF test.
38 a (mm)
40
42
44
Fig. 13. R-curves in mode II obtained in the cohesive model.
Table 4 Results obtained for the ENF specimen Reference value GIIc ¼ 0:6 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
CCM
CBT
DBT
CCM
CBT
DBT
CCM
CBT
DBT
0.61 1.67 0.28
0.524 12.66 0.44
0.583 2.8 0.8
0.713 18.88 0.88
0.573 4.47 0.62
0.592 1.34 1.66
0.798 33.04 0.8
0.632 5.41 0.91
0.661 10.11 2.16
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
stresses perpendicular to crack direction. The increase of mixed-mode local loading originates higher values of G IIc due to the pure mode damage propagation model used which, in these cases, is a conservative approach. In fact, the presence of a non-negligible local mixed-mode loading diminishes the shear component (mode II) at failure, which is the only one used in the pure damage model. Another effect that must be understood is the presence of the FPZ, which corresponds to a non-negligible region of damaged material. The size of this region influences the material compliance but is not accounted for in the beam theories for mode II. Numerically, the size of the FPZ corresponds to the group of points located after the crack tip undergoing the softening process and can be easily accessed. It was verified that the damage growth is selfsimilar for the propagation length considered, which means that the size of the FPZ is constant in this interval. To verify the influence of the FPZ on the G IIc using the beam theories, the energy released at each point included in this region is calculated. This energy corresponds to the triangular area (OBA) of Fig. 2. For a given point ðmÞ belonging to the FPZ it can be written as G m;II ¼ GIIc
du;II sm;II . 2
GIIðFPZÞ 1 ¼n ðs1;II þ s2;II þ þ sn;II Þ. su;II G IIc
(40)
This ratio represents the number of points completely damaged that must be added to the crack length to account for the energy released in the FPZ. This allows the definition of a correction Da to the real crack length or, in other words, the introduction of a concept of an equivalent crack, aeq ¼ a þ Da.
(41)
Fig. 15 shows the R-curves obtained using aeq in the cohesive damage model. In this case the FPZ measured 4.5 mm which led to Da ¼ 1:75 mm. Comparing with Fig. 13, it can be seen that the agreement with the reference value increased relatively to a non-corrected crack length (see also Tables 4 and 5). For the continuum damage model no remarkable differences were achieved in considering the concept of aeq (see Fig. 16 and Table 5). This reinforces the idea that the predominant effect is the presence of the non-negligible normal stresses at the crack tip. 5. Conclusions
du;II ðs1;II þ s2;II þ þ sn;II Þ 2
(39)
A two-dimensional numerical study was performed on the fracture tests characterization of adhesive joints. Two different models considering pure mode damage criterion were used to simulate crack propagation in the DCB 0.7 0.6 GIc (N/mm)
GIIc (N/mm)
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 32
or
(38)
Considering n points at the FPZ, the total energy released in this region is G IIðFPZÞ ¼ nGIIc
501
GIIc (ref) GIIc (CBT; E=4000 MPa) GIIc (DBT; E=4000 MPa) GIIc (CCM; E=4000 MPa) GIIc (CBT; E=1000 MPa) GIIc (DBT; E=1000 MPa) GIIc (CCM; E=1000 MPa)
0.5 0.4 GIIc (ref)
0.3
GIIc (CBT-a corr)
0.2
GIIc (DBT-a corr)
0.1 0 32
34
36
38
40
42
34
36
38 a (mm)
44
40
42
44
a (mm) Fig. 14. R-curves in mode II obtained in the continuum model.
Fig. 15. R-curves in mode II obtained in the cohesive model considering aeq .
Table 5 Results obtained for the ENF specimen considering the concept of aeq Reference value GIIc ¼ 0:6 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
CBT
DBT
CBT
DBT
CBT
DBT
0.591 1.49 0.12
0.606 0.94 0.08
0.614 2.39 0.55
0.609 1.48 1.21
0.643 7.09 0.85
0.666 10.93 2.03
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
502
the introduction of a concept of an equivalent crack. This approach resulted in better agreements between the beam theories and the reference value in the cohesive model. For the continuum model the effect was less prominent, which reinforces the conclusion that the mixed-mode local loading is the predominant effect.
0.7
GIIc (N/mm)
0.6 0.5 GIIc (ref) GIIc (CBT; E=4000 MPa) GIIc (DBT; E=4000 MPa) GIIc (CBT; E=1000 MPa) GIIc (DBT; E=1000 MPa)
0.4 0.3 0.2 0.1
Acknowledgements
0 32
34
36
38
40
42
44
a (mm) Fig. 16. R-curves in mode II obtained in the continuum model considering aeq .
The authors thank the Portuguese Foundation for Science and Technology (FCT, research project POCTI/ EME/45573/2002) for supporting the work here presented. References
and ENF tests. The cohesive damage model based on interface finite elements neglects the adhesive thickness. The continuum damage model considers the real thickness of the adhesive, which was simulated by solid plane elements. A bilinear softening model was considered in both cases. In the continuum model a characteristic length l c was introduced to transform the crack displacement into an equivalent strain. The advantage of the continuum over to the cohesive damage model is that it allows detecting the presence of other non-negligible stress components in the FPZ under non-symmetrical propagation. These are the normal and shear stresses acting in the crack direction in mode I and the normal stresses parallel and perpendicular to the crack plane in mode II. In the DCB test a good accuracy was obtained between the cohesive model and the reference value. Considering the continuum model it was verified that adhesive Young’s modulus is an important parameter. The agreement increases as the adhesive modulus decreases, independently of the data reduction scheme used. For higher values, the Young’s modulus of the adhesive becomes non-negligible relatively to the flexural modulus of the adherends. The conclusion is that the experimental fracture characterization in mode I using the DCB test can originate nonnegligible errors for lower ratios between the modulus of the adherends and that of the adhesive. Consequently, the fracture characterization of high modulus adhesives should be performed considering higher modulus adherends to increase the measurements accuracy. On the other hand, as the Gauss points are located out of the mid-plane of the specimen, the continuum model promotes non-symmetrical damage propagation typical of real experimental tests. The consequence is a mixed-mode loading at the crack tip that must be accounted for in damage propagation. These important phenomena are obviously not observed in the cohesive damage model. The same aspect was observed to be even more important in the ENF test. This reinforces the need of developing a mixed-mode continuum damage propagation model. Another feature studied was the effect of the FPZ size on the beam theories to obtain G IIc . The ratio between the energy released at the FPZ and the critical strain energy release rate in mode II allowed to
[1] Harris JA, Adams RD. Strength prediction of bonded single lap joints by non-linear finite element methods. International Journal of Adhesion and Adhesives 1948;4:65–78. [2] Crocombe AD, Bigwood DA, Richardson G. Analyzing structural adhesive joints for failure. International Journal of Adhesion and Adhesives 1990;10:167–78. [3] Czarnocki P, Piekarski K. Non-linear numerical stress analysis of a symmetric adhesive bonded lap shear joint. International Journal of Adhesion and Adhesives 1986;3:157–60. [4] Hamaush SA, Ahmad SH. Fracture energy release rate of adhesive joints. International Journal of Adhesion and Adhesives 1989; 9:171–8. [5] Anderson GP, Brinton SH, Ninow KJ, DeVries KL. A fracture mechanics approach to predicting bond strength. In: Advances in adhesively bonded joints. New York: ASME; 1988. p. 93–101. [6] Fernlund G, Papini M, McCammond D, Spelt JK. Fracture load predictions for adhesive joints. Composites Science and Technology 1994;51:587–600. [7] Pradhan SC, Iyengar NGR, Kishore NN. Parametric study of interfacial debonding in adhesively bonded composite joints. Composite Structures 1994;29:119–25. [8] Pradhan SC, Iyengar NGR, Kishore NN. Finite element analysis of crack growth in adhesively bonded joints. International Journal of Adhesion and Adhesives 1995;15:33–41. [9] Bogdanovich AE, Yushanov SP. 3-D Progressive failure analysis of bonded composite joints, In: AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics and materials conference, vol. 2. New York: AIAA; 1998, p. 1616–26. [10] Gonc- alves JPM, de Moura MFSF, Magalha˜es AG, de Castro PMST. Application of interface finite elements to threedimensional progressive failure analysis of adhesive joints. Fatigue and Fracture of Engineering Materials and Structures 2003;26: 479–86. [11] Li S, Thouless MD, Waas AM, Schroeder JA, Zavattieri PD. Use of mode-I cohesive-zone models to describe the fracture of an adhesively-bonded polymer-matrix composite. Composites Science and Technology 2005;65:281–93. [12] De D, Narasimhan R. Analysis of an interface fracture specimen for adhesively bonded joints. International Journal of Fracture 1998;92: L35–40. [13] Yang QD, Thouless MD, Ward SM. Elastic–plastic mode-II fracture of adhesive joints. International Journal of Solids and Structures 2001;38:3251–62. [14] Chai H. Micromechanics of shear deformations in cracked bonded joints. International Journal of Fracture 1992;58:223–39. [15] Kachanov LM. Time of the rupture process under creep conditions. Izvestiya Akademii Navk Tekhnicheskikh USSR 1958;8:26–31. [16] Rabotnov YN. Creep problems in structural members. Amsterdam, The Netherlands: North-Holland; 1969.
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503 [17] de Moura MFSF, Gonc- alves JPM, Marques AT, de Castro PMST. Modelling compression failure after low velocity impact on laminated composites using interface elements. Journal of Composite Materials 1997;31:1462–79. [18] Gonc- alves JPM, de Moura MFSF, de Castro PMST, Marques AT. Interface element including point-to-surface constraints for threedimensional problems with damage propagation. Engineering Computations: International Journal for Computer-Aided Engineering and Software 2000;17:28–47. [19] ISO/DIS 15024. Fibre-reinforced plastic composites: determination of mode I interlaminar fracture toughness, GIc, for unidirectionally reinforced materials, 2001
503
[20] Blackman BRK, Hadavinia H, Kinloch AJ, Paraschi M, Williams JG. The calculation of adhesive fracture energies in mode I: revisiting the tapered double cantilever beam (TDBC) test. Engineering Fracture Mechanics 2003;70(2):233–48. [21] Davies P, Sims GD, Blackman BRK, Brunner AJ, Kageyama K, Hojo M, Tanaka K, Murri G, Rousseau C, Giseke B, Martin RH. Comparision of test configurations for determination of mode II interlaminar fracture toughness results from international collaborative test programme. Plastics, Rubber and Composites 1999;28:432–7. [22] Carlsson LA, Gillespie JW, Pipes RB. On the analysis and design of the end notched flexure (ENF) specimen for mode II testing. Journal of Composite Materials 1986;20:594–604.
International Journal of Mechanical Sciences 48 (2006) 493–503 www.elsevier.com/locate/ijmecsci
Cohesive and continuum damage models applied to fracture characterization of bonded joints M.F.S.F. de Moura, J.A.G. Chousal Faculdade de Engenharia da Universidade do Porto, R. Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal Received 20 April 2005; received in revised form 27 November 2005; accepted 17 December 2005 Available online 7 February 2006
Abstract In this work, two different methods for simulating damage propagation are presented and applied to fracture characterization of bonded joints in pure modes I and II. The cohesive damage model is based on a special developed interface finite element including a linear softening damage process. In the continuum damage model the softening process is performed by including a characteristic length associated with a given Gauss point. The models were applied to the simulation of ‘‘double cantilever beam’’ (DCB) and ‘‘end notched flexure’’ (ENF) tests used to obtain the critical strain release rates in mode I and II of bonded joints. In mode I it was observed, under certain conditions, a good agreement between the results obtained by the two models with the reference value of critical strain energy release rate in mode I ðG Ic Þ, which is an inputted parameter. However, in mode II some discrepancies on the obtained GIIc values were observed between the two models. These inaccuracies can be explained by the simplifying assumptions inherent to the cohesive model. Better results were achieved considering the crack equivalent concept. r 2005 Elsevier Ltd. All rights reserved. Keywords: Damage models; Bonded joints; Fracture characterization
1. Introduction Adhesive joints are widely used in structural applications. However, it must be recognized that in the majority of the cases the applications are restricted to secondary non-critical structures. As a result, several different theories and models are being developed in order to increase the mechanical and structural performance of these joints and, obviously, the confidence of the designers. The joint strength prediction is usually based on the strength of materials or on the fracture mechanics models. Those based on the strength of materials generally use a stress or a strain failure criterion [1–3]. The fracture mechanics models are usually based on energetic criteria [4–10], namely on the relationship between energy release Abbreviations: FPZ, fracture process zone; DCB, double cantilever beam; ENF, end notched flexure; CBT, corrected beam theory; CCM, compliance calibration method; DBT, direct beam theory Corresponding author. Tel.: +351 225081727; fax: +351 225081584. E-mail address: [email protected] (M.F.S.F. de Moura). 0020-7403/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2005.12.008
rates with their critical values. Pradhan et al. [7,8] modelled the initiation of the debonding at the interfaces between the adherends and the adhesive. The authors considered pairs of nodes at these critical regions, which were released in function of the strain release rate. They concluded that this model can be used in the prediction of crack initiation and growth leading to failure in adhesive joints. Bogdanovich et al. [9] simulated the progressive failure analysis of double lap adhesive composite joints using the strain energy release rates to predict different possibilities of damage propagation: cohesive, adhesive or interlaminar crack propagation. The locus of initial failure and failure mode are determined by analysing the computed strains. Following, a progressive failure analysis is performed from the initial crack, which can be located at the interface adhesive–adherend, inside the adhesive or at interfaces between layers of the composite adherend. Gonc- alves et al. [10] used interface finite elements including a mixed-mode damage model to simulate progressive debonding in single lap adhesive joints. The damage model combines aspects of strength-based analysis and fracture mechanics and
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
494
Nomenclature a aeq a0 Da B C C D di D d di do;i du;i dr E E u;i F
crack length equivalent crack length initial crack length correction to the real crack length specimen width stiffness matrix for isotropic material compliance diagonal matrix containing the damage parameter damage parameter for mode i loading (i ¼ I,II) correction to account for crack tip rotation and deflection in mode I load-point displacement measured in the fracture test current relative displacement onset relative displacement of the softening process ultimate relative displacement of the softening process vector of relative displacements Young’s modulus diagonal matrix containing stiffness parameters ultimate strain on i direction correction factor accounting for large displacements for the moment arm
includes a linear softening process, which allows a smoothly decrease of stresses as the relative displacements grow. The area under the softening curve is equated to the critical energy release rate. The authors concluded that the model is accurate in predicting the joints mechanical behaviour. Li et al. [11] applied a three parameter cohesive-zone model to the study of mode-I fracture of adhesive joints made from a polymer–matrix composite. Cohesive-zone parameters were obtained by comparing numerical results to experimental observations. It was shown that there is a distinction between the characteristic strength of the interface associated with the toughness, and the intrinsic cohesive strength of the interface. Different cohesive-zone models were used for the interface and for the composite. It was concluded that cohesive-zone models accurately predict the joint strengths and deformations, and also the transition between failure of the composite and failure of the interface, which cannot be predicted by conventional fracture mechanics. The accuracy of the energetic criteria is highly dependent on the critical strain release rates of the adhesive joints. Several authors studied different specimens and methodologies to characterize the fracture mechanics of bonded joints. De et al. [12] developed an alternate specimen for measuring the fracture toughness of adhesively bonded joints. The specimen allows obtaining fracture toughness data over the entire range of mode mixity from pure mode
Fe fv G Gi G ic h I L lc n N N n P r su;i U u uL Vi
vector of external loads correction factor accounting for transverse shear effects shear modulus strain energy release rate for mode i loading critical strain energy release rate for mode i loading specimen half-thickness identity matrix length of DCB specimen and half-length of the ENF specimen characteristic length Poisson’s ratio shape function matrix correction factor accounting for large displacements for the compliance number of points at the FPZ load applied in the fracture test vector of interface stresses strength on i direction nodal displacement vector displacement vector at the integration points displacement vector at the integration points in local coordinates unit vectors of the local directions for mode i loading
I to pure mode II using a single geometry. The specimen is calibrated using the finite element method and mode mixity solutions assuming linear elastic fracture mechanics conditions. Yang et al. [13] performed a numerical study of the elastic–plastic mode-II fracture of adhesive joints. A traction-separation law was used to simulate the mode-II interfacial fracture of adhesively bonded end notch flexure (ENF) specimens with extensive plastic deformation accompanying failure. The authors found that the numerical predictions for the loads and deformation were in excellent agreement with the corresponding experimental results. The mode-II fracture was found to be much higher than mode-I and it was concluded that fracture is highly mode-dependent. A major contribution to the high value of GIIc is the higher capability of the adhesive layer to deform plastically in response to shear. Similar results were also observed by Chai [14] in his systematic studies of the fracture of elastic adhesive joints. These studies showed that, except for extremely thin adhesive layers ðo10 mmÞ, the typical value of GIIc is about an order of magnitude higher than the corresponding GIc value. The objective of this work is to evaluate the adequacy of the double cantilever beam (DCB) and ENF tests to obtain the critical strain energy release rates of bonded joints. A cohesive and a continuum damage models were developed and applied to the simulation of the experimental tests in order to understand the details of the
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
mechanical behaviour of the used specimens. Some conclusions were drawn about the satisfactoriness of these tests in fracture characterization of adhesive joints. 2. Damage models Damage can be considered as the progressive deterioration which occurs in materials before failure. It can be constituted by micro-cracking, voids nucleation and growth, and several inelastic processes that weaken the material. The analysis of cumulative damage is fundamental in life prediction of components and structures under loading. Kachanov [15] and Rabotnov [16] introduced the concept of simulating material damage by introducing damage variables into the constitutive equations. With this approach, the exact description of each individual micro-crack evolution is replaced by an averaging analysis describing the material macroscopic behaviour. In this work a cohesive and a continuum damage models are presented. In both cases, the material failure law can be defined as a relationship between normal (mode I) or shear (mode II) stresses on the crack plane and the corresponding normal and shear deformation variables (relative displacements in the cohesive model and strains in the continuum model). The objective is to evaluate the adequacy of each one in order to capture the phenomena related with the fracture characterization tests in pure modes of adhesive joints.
The cohesive damage model is based on using interface finite elements including a linear softening damage process. A user element subroutine UEL was developed to formulate a six-node interface finite element for two-dimensional problems (see Fig. 1) with the main ABAQUSs code. The stresses at the interface finite element are calculated from (1)
r ¼ Edr
being dr the vector of relative displacements between homologous points belonging to the top and bottom surfaces of the interface finite element in local coordinates. These are defined as associated to the failure modes (II, I) (2)
dr ¼ uL;top uL;bot or dr ¼
The matrix E is a diagonal matrix containing the stiffness parameters introduced by the user. They have to be carefully chosen to obtain a good performance of the model. As small values induce physically unacceptable interpenetrations and large values produce numerical problems, the optimum interface stiffnesses are the largest values that do not produce numerical problems. The interface stiffnesses were taken equal to 107 N=mm3 , which was shown to be a value satisfying the two referred conditions [17,18]. The displacement field expressed in the global system ðx; yÞ and associated to the element faces (top and bottom) can be obtained from utop;bot ¼ Ntop;bot Utop;bot ,
where U is the nodal displacement vector and N is the shape function matrix. The shape quadratic functions of the interface element are similar to those used in the eight node isoparametric plane solid element. The relation between displacements in local and global coordinates can be written as uL ¼ hT u,
dII dI
( ¼
uII uI
)
( L;top
(5)
being h ¼ ½VII ; VI ,
(6)
where VII ; VI are unit vectors of the local directions. VII can be written as T qx=qII; qy=qII; 0 VII ¼ (7) qx=qII; qy=qII; 0
VI ¼ k VII .
dr ¼ hT Ntop Utop hT Nbot Ubot
uII uI
. L;bot
I II
Fig. 1. Interface element.
(3)
(9)
or dr ¼ BUe ,
(10)
where B ¼ hT ½Ntop ; Nbot and U ¼
)
(8)
Substituting Eqs. (4) and (5) in Eq. (2) it can be established a relation between the local relative displacements in function of the nodal displacement vector
e
)
(4)
and VI is normal to the plane defined by VII and k ¼ ð0; 0; 1ÞT
2.1. Cohesive damage model
(
495
(
Utop Ubot
(11)
) .
(12)
The interface stresses can now be obtained from Eq. (1). The stiffness matrix can be obtained from the functional of the potential energy Z 1 e dT r dA UeT Fe , (13) pðU Þ ¼ 2 A r where Fe is the vector of external loads. The state of equilibrium is achieved when the displacement field Ue
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
496
minimizes the potential energy which leads to Z BT EB dA Ue , Fe ¼
(14)
A
where the integral between parentheses corresponds to the stiffness matrix which, with Fe , represent the contributions of the interface element to the global problem. The interface element incorporates a pure mode damage model to simulate damage initiation and growth. The damage model combines aspects of strength-based analysis and fracture mechanics. It is based on a softening process between stresses and interfacial relative displacements. To avoid the singularity at the crack tip and its effects, a gradual degradation is considered instead of a sudden one which would result in mesh-dependency. It is assumed that failure occurs gradually as energy is dissipated in a cohesive zone behind the crack tip. This is equivalent to the consideration of a ‘‘fracture process zone’’ (FPZ), defined as the region in which the material undergoes softening damage by different ways, e.g., micro-cracking and inelastic processes. This is achieved by including a damage parameter which moves smoothly from zero to unity as the material fails. For pure mode (i ¼ I or II) loading, the linear softening process starts when the interfacial stress reaches the respective strength su;i (see Fig. 2). The softening relationship can be written as r ¼ ðI DÞEdr ,
(15)
where I is the identity matrix and D is a diagonal matrix containing on the i-position the damage parameter which represents the damage accumulated at the interface di ¼
du;i ðdi do;i Þ , di ðdu;i do;i Þ
(16)
where di is the current relative displacement and do;i , du;i the onset and ultimate relative displacements of the softening process. This damage parameter is zero for the undamaged state and reaches one when the material is fully damaged. In pure modes, the other component of stress is not important and can be abruptly cancelled equating to one the respective element of the diagonal matrix D. Due to i u,i
B
A
0
C o,i
u,i
i
Fig. 2. Softening stresses/relative displacements relationship for the cohesive damage model.
the irreversibility of the damage process, the damage parameter d i for a material point at a given increment j is determined using the following equation: d i ¼ maxðd ji ; d j1 i Þ,
(17)
which avoids an healing phenomenon. If the load reverses during the softening process, the unloading and reloading procedures are represented by the line O–A in Fig. 2, and the released energy is given by Gi ¼ 12du;i ðsu;i si Þ.
(18)
The maximum relative displacement du;i , for which complete failure occurs, is obtained by equating the area under the softening curve (triangle OBC of Fig. 2) to the respective critical fracture energy Gic ¼ 12su;i du;i .
(19)
Due to the softening process is usual that convergence difficulties arise during the non-linear numerical analysis which is performed by using the Newton–Raphson method. In fact, although loading displacement was applied gradually considering small increments to ensure a smooth propagation process, a lot of iterations are needed for each increment. To overcome this problem the tangent stiffness matrix is considered. This matrix is defined from the tangent modulus matrix that can be obtained from the variation of Eq. (15) qr ¼ q½ðI DÞEdr þ ðI DÞE qdr
(20)
or qr ¼
dr;o E qdr , dr;u dr;o
(21)
where dr;o and dr;u represent the onset and ultimate relative displacements vectors. The tangent modulus matrix ðEtan Þ is the operator that relates qr with qdr and is used in Eq. (14) instead of E to define the tangent stiffness matrix. 2.2. Continuum damage model Finite element simulations using continuum damage models are known to be susceptible to mesh sensitivity. This means that refining the finite element discretization does not allow the numerical solution to converge to a physically meaningful solution of the problem. To overcome this problem, the continuum isotropic pure mode damage model presented here is based on a similar softening law used in the cohesive damage model and it was implemented via USDFLD user subroutine of ABAQUS Standards software. The basic assumption is to establish a correlation between the increase of damage and the decrease in stiffness at a given material point of a plane solid element. The adopted solution is based on the consideration of a damage parameter in a particular direction related to a given length associated with a Gauss point. The stiffness degradation is related to the critical strain energy release rate. However, in this case
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
i
497
tion between crack surfaces is avoided considering the initial modulus. Quadratic solid plane elements with eight nodes and four Gaussian points were considered in the simulations. The damage parameter is updated in the user subroutine during the incremental non-linear process affecting the referred elastic properties.
u,i
A
3. Models validation 0
o,i
u,i
i
Fig. 3. Softening stresses/strains relationship for the continuum damage model.
a characteristic length l c must be introduced to transform the crack opening (mode I) or shear (mode II) displacements into an equivalent strain (see Fig. 3) G ic ¼
1 2su;i u;i l c .
(22)
This parameter was considered to be equal to the length of influence of a Gauss point in the given direction and physically can be regarded as the dimension at which the material acts homogeneously. A specific algorithm was implemented to evaluate the characteristic length corresponding to each Gauss point. Basically it consists in computing the dimension of the element from the nodal coordinates and taking into account the number of Gauss points in the referred direction. Considering isotropic damage evolution, the stress– strain relation can be written considering an equation similar to Eq. (15) r ¼ ð1 d i ÞCe,
(23)
being d i the scalar damage parameter which grows from zero (undamaged initial state) to one (complete loss of integrity), and C the stiffness matrix of the undamaged material. The damage parameter is calculated by an expression similar to Eq. (16) but considering strains instead of relative displacements (see Fig. 3). According to the model assumptions the material stiffness remains isotropic during its progressive degradation and the Poisson’s ratio is not affected by damage. However, there is an implicit shear coupling. In fact, due to damage growth, the Young’s modulus decreases and, consequently, the shear modulus of the material also decreases due to the relationship G¼
E . 2ð1 þ nÞ
(24)
The softening law is similar to the one used in the cohesive model and the healing phenomenon is also avoided using the statement of Eq. (17). The properties are smoothly reduced due to the energy released at the FPZ. When contact occurs during crack propagation, the interpenetra-
In the present work a beam theory and fracture mechanics analyses are applied to obtain analytical P2d curves for the DCB and ENF specimens in order to validate the numerical models. The transverse shear deformations effects were neglected. 3.1. DCB specimen Considering the simple beam theory the relationship between the applied load and deflection at the end of the specimen, before damage starts to grow, can be written as 8a30 P, (25) EBh3 being B the width and h the half height of the specimen. During damage propagation ða4a0 Þ the Irwin–Kyes expression
d¼
P2 dC , (26) 2B da where C is the compliance ðC ¼ d=PÞ, is used to obtain the crack length a in function of G Ic 1=2 B Eh3 G Ic a¼ . (27) P 12 Gic ¼
Substituting Eq. (27) in Eq. (25) the P2d relationship during propagation is B2 EhG Ic 3=2 . (28) d¼ 3 EP2 3.2. ENF specimen Following the beam theory the displacement at the loading point, before damage initiation, can be written as 3a30 þ 2L3 P. (29) 8EBh3 During crack growth and for crack length smaller than L ðaoLÞ the displacement at the mid-span is obtained from Eq. (29) considering a instead of a0 . Combining Eqs. (29) and (26) the crack length can be obtained from 1=2 4B Eh3 G IIc a¼ , (30) 3P d¼
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
498
4.1. Double cantilever beam test
which, after substitution in Eq. (29) leads to P 32B3 3 3=2 3 d¼ L þ ðEh G IIc Þ . 9P3 4EBh3
(31)
For crack length larger than L ða4LÞ the beam theory originates d¼
P 3 3 3 ð2L aÞ L 8 EBh3
(32)
and the critical strain energy release rate is G IIc ¼
9P2 ð2L aÞ2 . 16B2 Eh3
(33)
The DCB test (see Fig. 4) is used to obtain the mode I critical strain release rate ðG Ic Þ. In the cohesive damage model, interface elements with null thickness were placed between the adherends along the specimen length to simulate the adhesive (see Fig. 5). In the continuum damage model the adhesive was modelled by two-dimensional solid elements with 0.25 mm thickness (see Fig 6). During load history the values of force, applied displacement and crack length ðP2d2aÞ were registered in order to calculate the critical strain energy release rate by using the compliance calibration method (CCM) using Eq. (26). Cubic polynomials ðC ¼ C 3 a3 þ C 2 a2 þ C 1 a þ C 0 Þ were used to fit the C ¼ f ðaÞ curves. The corrected beam theory (CBT) was also used to predict the R-curve
Substituting ð2L aÞ of Eq. (33) in Eq. (32) gives P 8B3 3 3=2 3 d¼ L 3 ðEh G IIc Þ . 9P EBh3
GIc ¼ (34)
3Pd , 2Bða þ DÞ
(35)
P,
4. Results The presented damage models were applied to the simulation of fracture characterization tests of bonded joints. One of the main differences between the two models is related with the thickness of the adhesive. In the cohesive model the adhesive thickness is neglected, i.e. the interface elements connecting the adherends have null thickness. On the contrary, the continuum model considers the real adhesive thickness (a typical value of 0.25 mm was considered). The analyses were performed considering plane strain conditions and non-linear geometrical behaviour. A loading displacement was applied progressively considering small increments to ensure a smooth propagation process. The used meshes have 2 280 solid plane elements in each arm. The specimens’ geometry and materials properties are presented in Tables 1 and 2 and Figs. 4 and 9.
t
h h a0
L Fig. 4. The DCB test.
Fig. 5. Mesh of the DCB specimen used in the cohesive damage model.
Fig. 6. Mesh of the DCB specimen used in the continuum damage model.
Table 1 Specimens geometry
30 100 50 1.5 0.25 30 10
25 20 P (N)
DCB - L (mm) ENF-L (mm) h (mm) t (mm) a0 (mm) B (mm)
15 Analytical model Cohesive damage model Continuum damage model (E=400 MPa) Continuum damage model (E=1000 MPa) Continuum damage model (E=4000 MPa) Continuum damage model (2 el., E=400 MPa)
10 Table 2 Materials properties
5
Adherends
Adhesive
0
E(MPa) n
E (MPa) n
69 000
0.33 4000
su (MPa) su (MPa) GIc (N/mm) GIIc (N/mm)
0.3 20
40
0.3
0.6
0
1
2
(mm)
3
Fig. 7. P2d curves for the DCB test.
4
5
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
where D represents a correction for crack tip rotation and deflection and is obtained from a linear regression analysis of C 1=3 versus a data [19]. Figs. 7 and 8 show the P2d and G ic ¼ f ðaÞ relations for the cohesive and continuum models. Good agreement between the analytical, cohesive and continuum ðE ¼ 400 MPaÞ damage models was observed in the post-peak region of the P2d curves, which corresponds to damage propagation. The higher peak load and stiffness value obtained by the analytical model can be explained by the softening onset that occurs in the earlier stages of the numerical P2d curves. In fact, it was observed that several points initiate the softening process a lot before the maximum load is attained, which obviously affects the stiffness. These effects are not accounted for in the analytical approach. In the post-peak region, it can also be observed some differences between the P2d curves of the cohesive and the continuum model. The agreement between the two models increases if we consider adhesives of lower modulus (E ¼ 1000 and 400 MPa) in the continuum damage model. A similar tendency can be observed for the R-curves considering different data reduction schemes. Table 3 presents the results obtained for both damage models. It can be concluded that cohesive model presents good accuracy relatively to the reference value, which is a material property introduced in the simulations (see Table 2). In the continuum model the difference relatively to the reference value increases with the Young’s modulus of the adhesive. In fact for E ¼ 4000 MPa the Young’s modulus of the adhesive becomes
499
non-negligible relatively to the flexural modulus of the adherends. It must be noted that this behaviour cannot be simulated in the cohesive model as it neglects the adhesive thickness. Similar results were obtained by Blackman et al. [20] on the analysis of the tapered DCB adhesive joint specimen. On the other hand, non-symmetrical damage propagation is simulated. In fact, the crack growth occurs at the Gauss points, which are located out of the mid-plane of the specimen. It must be noted that this non-symmetrical damage propagation occurs frequently in real tests. The consequence is the presence at the FPZ of non-negligible normal and shear stresses parallel to the crack inducing local mixed-mode loading and influencing the load–displacement relationship. Once again it must be emphasized that these effects are not captured when a cohesive damage model is used because adhesive thickness is neglected. To evaluate the influence of mesh refinement, an analysis considering two elements along the adhesive thickness direction was carried out. As observed in Fig. 7, the obtained P2d curve agrees quite well with the single element one, showing the model insensitivity to mesh refinement. 4.2. End notch flexure test The ENF test (see Fig. 9) is used to obtain the mode II critical strain release rate ðG IIc Þ. At the initial crack contact
P, a0
0.35
GIc(N/mm)
0.3
h t h
0.25 GIc (reference value) GIc (Cohesive model-CCM) GIc (Cohesive model-CBT) GIc (CCM, E=4000 MPa) GIc (CBT, E=4000 MPa) GIc (CCM, E=1000 MPa) GIc (CBT, E=1000 MPa) GIc (CCM, E=400 MPa) GIc (CBT, E=400 MPa)
0.2 0.15 0.1 0.05
L
L
Fig. 9. The ENF test.
0 30
32
34
36 a (mm)
38
40
42
Fig. 10. Detail of the ENF specimen in the cohesive damage model.
Fig. 8. The R-curves for DCB test.
Table 3 Results obtained for the DCB specimen Reference value G Ic ¼ 0:3 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
E ¼ 400 MPa
CCM
CBT
CCM
CBT
CCM
CBT
CCM
CBT
0.306 1.86 0.83
0.302 0.7 0.04
0.229 10.2 0.16
0.27 10.1 0.13
0.287 4.19 0.24
0.288 3.95 0.126
0.299 0.41 0.17
0.299 0.44 0.14
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
500
conditions were imposed to prevent interpenetration of the cracked parts. Details of the deformed shapes in both models are shown in Figs. 10 and 11. Fig. 12 presents the P2d curves obtained from the analytical and numerical approaches. In the post-peak region excellent agreement between the cohesive, continuum ðE ¼ 4000 MPaÞ and analytical models. The initial stiffness in the continuum damage model is greater than in the cohesive damage model. Considering a lower modulus adhesive ðE ¼ 1000 MPaÞ a slightly decrease on the initial stiffness was observed. The differences relatively to the cohesive and analytical models can be explained by the higher specimen thickness (plus 0.25 mm) considered in the continuum damage model. As observed in the DCB case, the mesh refinement along adhesive thickness (two elements) has a small influence on the P2d curve (see Fig. 12). Three different data reduction schemes were used to obtain the R-curves in mode II:
The CCM was also applied considering Eq. (26) and similar procedures described in the DCB test
Fig. 11. Details of the ENF specimen in the continuum damage model.
P (N)
200
G IIc ¼
9Pda2 . 2Bð2L3 þ 3a3 Þ
(36)
The CBT [21,22] 9Pa2 F f v, (37) 16B2 Eh3 N where F and N are corrections to account for large displacements for the moment arm and for the compliance, respectively, and f v comprise the effect of transverse shear G IIc ¼
From Fig. 13 and Table 4 we can conclude that in the cohesive model the best agreement was obtained by the CCM, and the worst by the CBT. However, in the continuum model the beam theories present better agreement with the reference value than the CCM (see Fig. 14 and Table 4). As in the DCB test non-symmetrical damage propagation was simulated. Consequently there exist at the FPZ non-negligible normal stresses parallel and perpendicular to the crack inducing local mixed-mode loading influencing the compliance, which explain the poor agreement of the CCM. The decrease of the resin modulus increases the G IIc value obtained from the different data reduction schemes (see Table 4). However, contrary to the DCB test the agreement does not increase. In fact, the lower modulus induces a more pronounced local deformation at the crack tip leading to more important normal
Analytical model Cohesive damage model Continuum damage model (E=4000MPa) Continuum damage model (E=1000MPa) Continuum damage model (2 el., E=1000 MPa)
0.7 0.6 GIIc (N/mm)
250
The direct beam theory (DBT) based on elementary beam theory [21]
150 100
0.5 0.4
GIIc (ref) GIIc (CBT)
0.3
GIIc (CCM)
0.2
50
GIIc (DBT)
0.1 0 0.0
1.0
2.0
3.0
4.0
0 32
34
36
(mm) Fig. 12. P2d curves for the ENF test.
38 a (mm)
40
42
44
Fig. 13. R-curves in mode II obtained in the cohesive model.
Table 4 Results obtained for the ENF specimen Reference value GIIc ¼ 0:6 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
CCM
CBT
DBT
CCM
CBT
DBT
CCM
CBT
DBT
0.61 1.67 0.28
0.524 12.66 0.44
0.583 2.8 0.8
0.713 18.88 0.88
0.573 4.47 0.62
0.592 1.34 1.66
0.798 33.04 0.8
0.632 5.41 0.91
0.661 10.11 2.16
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
stresses perpendicular to crack direction. The increase of mixed-mode local loading originates higher values of G IIc due to the pure mode damage propagation model used which, in these cases, is a conservative approach. In fact, the presence of a non-negligible local mixed-mode loading diminishes the shear component (mode II) at failure, which is the only one used in the pure damage model. Another effect that must be understood is the presence of the FPZ, which corresponds to a non-negligible region of damaged material. The size of this region influences the material compliance but is not accounted for in the beam theories for mode II. Numerically, the size of the FPZ corresponds to the group of points located after the crack tip undergoing the softening process and can be easily accessed. It was verified that the damage growth is selfsimilar for the propagation length considered, which means that the size of the FPZ is constant in this interval. To verify the influence of the FPZ on the G IIc using the beam theories, the energy released at each point included in this region is calculated. This energy corresponds to the triangular area (OBA) of Fig. 2. For a given point ðmÞ belonging to the FPZ it can be written as G m;II ¼ GIIc
du;II sm;II . 2
GIIðFPZÞ 1 ¼n ðs1;II þ s2;II þ þ sn;II Þ. su;II G IIc
(40)
This ratio represents the number of points completely damaged that must be added to the crack length to account for the energy released in the FPZ. This allows the definition of a correction Da to the real crack length or, in other words, the introduction of a concept of an equivalent crack, aeq ¼ a þ Da.
(41)
Fig. 15 shows the R-curves obtained using aeq in the cohesive damage model. In this case the FPZ measured 4.5 mm which led to Da ¼ 1:75 mm. Comparing with Fig. 13, it can be seen that the agreement with the reference value increased relatively to a non-corrected crack length (see also Tables 4 and 5). For the continuum damage model no remarkable differences were achieved in considering the concept of aeq (see Fig. 16 and Table 5). This reinforces the idea that the predominant effect is the presence of the non-negligible normal stresses at the crack tip. 5. Conclusions
du;II ðs1;II þ s2;II þ þ sn;II Þ 2
(39)
A two-dimensional numerical study was performed on the fracture tests characterization of adhesive joints. Two different models considering pure mode damage criterion were used to simulate crack propagation in the DCB 0.7 0.6 GIc (N/mm)
GIIc (N/mm)
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 32
or
(38)
Considering n points at the FPZ, the total energy released in this region is G IIðFPZÞ ¼ nGIIc
501
GIIc (ref) GIIc (CBT; E=4000 MPa) GIIc (DBT; E=4000 MPa) GIIc (CCM; E=4000 MPa) GIIc (CBT; E=1000 MPa) GIIc (DBT; E=1000 MPa) GIIc (CCM; E=1000 MPa)
0.5 0.4 GIIc (ref)
0.3
GIIc (CBT-a corr)
0.2
GIIc (DBT-a corr)
0.1 0 32
34
36
38
40
42
34
36
38 a (mm)
44
40
42
44
a (mm) Fig. 14. R-curves in mode II obtained in the continuum model.
Fig. 15. R-curves in mode II obtained in the cohesive model considering aeq .
Table 5 Results obtained for the ENF specimen considering the concept of aeq Reference value GIIc ¼ 0:6 N=mm
Cohesive model
Continuum model E ¼ 4000 MPa
GIIc (N/mm) Error (%) Std. dev. (%)
E ¼ 1000 MPa
CBT
DBT
CBT
DBT
CBT
DBT
0.591 1.49 0.12
0.606 0.94 0.08
0.614 2.39 0.55
0.609 1.48 1.21
0.643 7.09 0.85
0.666 10.93 2.03
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503
502
the introduction of a concept of an equivalent crack. This approach resulted in better agreements between the beam theories and the reference value in the cohesive model. For the continuum model the effect was less prominent, which reinforces the conclusion that the mixed-mode local loading is the predominant effect.
0.7
GIIc (N/mm)
0.6 0.5 GIIc (ref) GIIc (CBT; E=4000 MPa) GIIc (DBT; E=4000 MPa) GIIc (CBT; E=1000 MPa) GIIc (DBT; E=1000 MPa)
0.4 0.3 0.2 0.1
Acknowledgements
0 32
34
36
38
40
42
44
a (mm) Fig. 16. R-curves in mode II obtained in the continuum model considering aeq .
The authors thank the Portuguese Foundation for Science and Technology (FCT, research project POCTI/ EME/45573/2002) for supporting the work here presented. References
and ENF tests. The cohesive damage model based on interface finite elements neglects the adhesive thickness. The continuum damage model considers the real thickness of the adhesive, which was simulated by solid plane elements. A bilinear softening model was considered in both cases. In the continuum model a characteristic length l c was introduced to transform the crack displacement into an equivalent strain. The advantage of the continuum over to the cohesive damage model is that it allows detecting the presence of other non-negligible stress components in the FPZ under non-symmetrical propagation. These are the normal and shear stresses acting in the crack direction in mode I and the normal stresses parallel and perpendicular to the crack plane in mode II. In the DCB test a good accuracy was obtained between the cohesive model and the reference value. Considering the continuum model it was verified that adhesive Young’s modulus is an important parameter. The agreement increases as the adhesive modulus decreases, independently of the data reduction scheme used. For higher values, the Young’s modulus of the adhesive becomes non-negligible relatively to the flexural modulus of the adherends. The conclusion is that the experimental fracture characterization in mode I using the DCB test can originate nonnegligible errors for lower ratios between the modulus of the adherends and that of the adhesive. Consequently, the fracture characterization of high modulus adhesives should be performed considering higher modulus adherends to increase the measurements accuracy. On the other hand, as the Gauss points are located out of the mid-plane of the specimen, the continuum model promotes non-symmetrical damage propagation typical of real experimental tests. The consequence is a mixed-mode loading at the crack tip that must be accounted for in damage propagation. These important phenomena are obviously not observed in the cohesive damage model. The same aspect was observed to be even more important in the ENF test. This reinforces the need of developing a mixed-mode continuum damage propagation model. Another feature studied was the effect of the FPZ size on the beam theories to obtain G IIc . The ratio between the energy released at the FPZ and the critical strain energy release rate in mode II allowed to
[1] Harris JA, Adams RD. Strength prediction of bonded single lap joints by non-linear finite element methods. International Journal of Adhesion and Adhesives 1948;4:65–78. [2] Crocombe AD, Bigwood DA, Richardson G. Analyzing structural adhesive joints for failure. International Journal of Adhesion and Adhesives 1990;10:167–78. [3] Czarnocki P, Piekarski K. Non-linear numerical stress analysis of a symmetric adhesive bonded lap shear joint. International Journal of Adhesion and Adhesives 1986;3:157–60. [4] Hamaush SA, Ahmad SH. Fracture energy release rate of adhesive joints. International Journal of Adhesion and Adhesives 1989; 9:171–8. [5] Anderson GP, Brinton SH, Ninow KJ, DeVries KL. A fracture mechanics approach to predicting bond strength. In: Advances in adhesively bonded joints. New York: ASME; 1988. p. 93–101. [6] Fernlund G, Papini M, McCammond D, Spelt JK. Fracture load predictions for adhesive joints. Composites Science and Technology 1994;51:587–600. [7] Pradhan SC, Iyengar NGR, Kishore NN. Parametric study of interfacial debonding in adhesively bonded composite joints. Composite Structures 1994;29:119–25. [8] Pradhan SC, Iyengar NGR, Kishore NN. Finite element analysis of crack growth in adhesively bonded joints. International Journal of Adhesion and Adhesives 1995;15:33–41. [9] Bogdanovich AE, Yushanov SP. 3-D Progressive failure analysis of bonded composite joints, In: AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics and materials conference, vol. 2. New York: AIAA; 1998, p. 1616–26. [10] Gonc- alves JPM, de Moura MFSF, Magalha˜es AG, de Castro PMST. Application of interface finite elements to threedimensional progressive failure analysis of adhesive joints. Fatigue and Fracture of Engineering Materials and Structures 2003;26: 479–86. [11] Li S, Thouless MD, Waas AM, Schroeder JA, Zavattieri PD. Use of mode-I cohesive-zone models to describe the fracture of an adhesively-bonded polymer-matrix composite. Composites Science and Technology 2005;65:281–93. [12] De D, Narasimhan R. Analysis of an interface fracture specimen for adhesively bonded joints. International Journal of Fracture 1998;92: L35–40. [13] Yang QD, Thouless MD, Ward SM. Elastic–plastic mode-II fracture of adhesive joints. International Journal of Solids and Structures 2001;38:3251–62. [14] Chai H. Micromechanics of shear deformations in cracked bonded joints. International Journal of Fracture 1992;58:223–39. [15] Kachanov LM. Time of the rupture process under creep conditions. Izvestiya Akademii Navk Tekhnicheskikh USSR 1958;8:26–31. [16] Rabotnov YN. Creep problems in structural members. Amsterdam, The Netherlands: North-Holland; 1969.
ARTICLE IN PRESS M.F.S.F. de Moura, J.A.G. Chousal / International Journal of Mechanical Sciences 48 (2006) 493–503 [17] de Moura MFSF, Gonc- alves JPM, Marques AT, de Castro PMST. Modelling compression failure after low velocity impact on laminated composites using interface elements. Journal of Composite Materials 1997;31:1462–79. [18] Gonc- alves JPM, de Moura MFSF, de Castro PMST, Marques AT. Interface element including point-to-surface constraints for threedimensional problems with damage propagation. Engineering Computations: International Journal for Computer-Aided Engineering and Software 2000;17:28–47. [19] ISO/DIS 15024. Fibre-reinforced plastic composites: determination of mode I interlaminar fracture toughness, GIc, for unidirectionally reinforced materials, 2001
503
[20] Blackman BRK, Hadavinia H, Kinloch AJ, Paraschi M, Williams JG. The calculation of adhesive fracture energies in mode I: revisiting the tapered double cantilever beam (TDBC) test. Engineering Fracture Mechanics 2003;70(2):233–48. [21] Davies P, Sims GD, Blackman BRK, Brunner AJ, Kageyama K, Hojo M, Tanaka K, Murri G, Rousseau C, Giseke B, Martin RH. Comparision of test configurations for determination of mode II interlaminar fracture toughness results from international collaborative test programme. Plastics, Rubber and Composites 1999;28:432–7. [22] Carlsson LA, Gillespie JW, Pipes RB. On the analysis and design of the end notched flexure (ENF) specimen for mode II testing. Journal of Composite Materials 1986;20:594–604.