<figure><img src="image_1.png" /></figure>

# Numerical simulations of crack propagation tests in adhesive bonded joints

**Abstract**

Mainly due to their low weight, low cost and ease of assembly, the adhesive bonds have emerged as a promising technology. However, the lack of adequate tools of design and control remain an obstacle to the use of the adhesives. In this work a cohesive interface model formulated within the framework of damage mechanics is applied for the simulation of decohesion during crack propagation tests. Considering the mechanical tests of aluminium/epoxy specimens, comparisons between experimental and numerical results are presented.

**Keywords**

bonded joints; cohesive-zone model; crack propagation tests

Silvio de Barros<sup>a,*</sup>,
Laurent Champaney<sup>b</sup> and
Nunziante Valoroso<sup>c</sup>

Received 09 Apr 2012;
In revised form 24 Apr 2012

* Author email: silvio.debarros@gmail.com

# 1 INTRODUCTION

The widespread use of adhesive joints is indicative of the advantages of bonding techniques over welding and riveting ones. Contrary to holes, rivets, clamps and screws that have a tendency to cause stress concentration areas, adhesives can distribute the load over the entire bonded area. However, the use of adhesive bonding in aircraft structures and other safety critical applications has been limited due to the lack of adequate tools of design and control. The development of numerical tools of design is necessary to increase the utilization of bonded joints in the industry. Interface damage models have been extensively used for the non-linear incremental analysis of debonding in the last years [1, 3, 4]. This damage models use some parametres that can be identified from mechanical tests.

The paper investigates the validation of a damage interface model for the simulation of crack propagations test. Comparisons between experimental results in aluminium/epoxy specimens and numerical simulation are presented.

<footnote><sup>a</sup> Department of Mechanical Engineering. CE-FET/RJ, Rio de Janeiro Av. Maracan, 229 - CEP 20271-110 Brazil</footnote>

<footnote><sup>b</sup> LMT Cachan 61, av du Pt Wilson - 94235 Cachan Cedex France</footnote>

<footnote><sup>c</sup> Dipartimento per le Tecnologie, Università di Napoli Parthenope, Centro Direzionale Isola C4, 80143 Napoli, Italy</footnote>
## 2 INTERFACE MODEL

In the cohesive-zone approach the description of a state of damage along an interface relies upon the definition of a traction-separation law incorporating the dependence of the surface tractions on the corresponding displacement discontinuities $[u] = u^+ - u^-$ and the damage criterion to be met for the cohesive process zone to grow and the crack advance. In the simplest one-dimensional case the damage onset and decohesion propagation conditions only involve the single-mode displacement or energy release rate component; on the contrary, when considering the mixed-mode case these conditions have to properly account for the interaction of the pure-mode contributions. In this last case the work of separation per unit fracture area does actually results from the interplay of the I and II pure-mode contributions, that are not independent in that they evolve together as a consequence of the interaction of the traction-displacement jump relationships in two directions. In what follows we shall briefly discuss the cohesive-zone model used in this work. A more exhaustive presentation of this model can be found in [10].

### 2.1 Pure-mode model

The adhesive joint here considered consists of two elastic bodies (adherends) joined by a plane adhesive layer whose thickness is assumed to be negligible compared to both that of the joined bodies and to its in-plane dimensions. These features enable the adhesive layer to be conveniently schematized as an interface, i.e. as a zero-thickness surface entity which ensures displacement and stress transfer between the adherends, see Fig. 1.

<figure><img src="image_2.png" /><figcaption>Figure 1 Interface schematization.</figcaption></figure>

Assuming that the displacement jump $[u] = u^+ - u^-$ at the interface in one direction is small in the usual sense, the elastic damage model for the interface can be derived based on a stored energy function defined as:

$$ \psi ([u], D) = \frac{1}{2} (1 - D) k^+ \langle [u] \rangle_+^2 + \frac{1}{2} k^- \langle [u] \rangle_-^2 \quad (1) $$

where $D \in [0, 1]$ denotes a scalar damage variable in the usual sense, the symbols $\langle \cdot \rangle_+$ and $\langle \cdot \rangle_-$ stand for the positive and negative part of the argument $\langle \cdot \rangle$, defined as $\langle x \rangle_\pm = 1/2 (x \pm |x|)$, and $k^+$ and $k^-$ are the undamaged interface stiffnesses in tension and compression, respectively, the latter representing a penalty stiffness accounting for the impenetrability constraint.

The associated interface traction in the direction of the jump is then the following:
$$t = \frac{\partial \psi}{\partial [\mathbf{u}]} = (1 - D) k^{+} \langle [u] \rangle_{+} + k^{-} \langle [u] \rangle_{-} \quad (2)$$

The damage driving force is classicaly defined by:

$$Y = -\frac{\partial \psi}{\partial D} = \frac{1}{2} k^{+} \langle [u] \rangle_{+}^{2} \quad (3)$$

The damage evolution is subjected to the classical loading/unloading conditions:

$$f(Y) \le 0 \quad \dot{D} \ge 0 \quad f(Y) \dot{D} = 0 \quad (4)$$

$$f(Y) = Y - Y^* \quad \dot{D} \in [0, 1] \quad (5)$$

where the damage threshold $Y^*$ is defined by:

$$\left\{ \begin{array}{ll} Y^* = & G_o \quad \text{if } D = 0 \\ Y^* = & G_o + (Y_f - G_o) [-\log(1-D)]^N \quad \text{if } D \in ]0, 1[ \\ Y^* = & \max_{\tau \in [0, T]} Y(\tau) \quad \text{if } D = 1 \end{array} \right. \quad (6)$$

The energy dissipated in the decohesion process is:

$$\int_{0}^{1} Y^{*}(D) dD = G_{o} + (Y_{f} - G_{o}) \Gamma (N + 1) = G_{c} \quad (7)$$

where $\Gamma$ is the Gamma function [5], defined by:

$$\Gamma (N + 1) = \int_{0}^{+\infty} x^{N} e^{-x} dx = N \cdot \Gamma (N) \quad (8)$$

The traction-separation relationship for this model is depicted in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Figure 2 Traction-separation relationships.</figcaption></figure>
## 2.2 Mixed-mode model

As opposite to the single-mode case, where the criteria used for determining damage onset and propagation up to complete failure only involve one single component of the energy release rate, when considering mixed-mode conditions the total energy released during decohesion ($G_T$) results from the interplay of the I and II pure-mode contributions that evolve together as a consequence of the interactions between the traction-displacement jump in the normal and tangential directions.

$$G_T = G_I + G_{II} \qquad (9)$$

The stored energy function takes the following form:

$$\psi([\mathbf{u}], D) = \frac{1}{2}(1 - D) \left[ k_n^+ \langle [u_n] \rangle_+^2 + k_s [u_s]^2 \right] + \frac{1}{2} k_n^- \langle [u_n] \rangle_-^2 \qquad (10)$$

where $k_n^+$ and $k_n^-$ are the interface stiffness normal component in tension and compression, $k_s$ is the interface stiffness sliding (tangential) component, and $[u_n]$ and $[u_s]$ denote the normal and sliding components of the displacement jump vector $[u]$, i.e. $[u_n] = [u] \cdot n$, $[u_s] = [u] \cdot n$, $n$ and $s$ being the outward unit normal and the unit tangent vector to the interface, see also Fig. 1.

The constitutive equations for the interface traction vector $t$ and the damage driving force are obtained in the usual way as:

$$\left\{ \begin{array}{l} \mathbf{t} = \frac{\partial \psi}{\partial [\mathbf{u}]} = (1 - D) [k_n^+ \langle [u_n] \rangle_+ \mathbf{n} + k_s [u_s] \mathbf{s}] + k_n^- \langle [u_n] \rangle_- \mathbf{n} \\ Y_m = -\frac{\partial \psi}{\partial D} = Y_I + Y_{II} \end{array} \right. \qquad (11)$$

The energy release rate for the two modes are:

$$\left\{ \begin{array}{l} Y_I = \frac{1}{2} k_n^+ \langle [u_n] \rangle_+^2 \\ Y_{II} = \frac{1}{2} k_s [u_s]^2 \end{array} \right. \qquad (12)$$

In the above equation and in the remainder of the paper the subscript $m$ is appended to the mixed-mode variables in order to emphasize the difference with the analogous unsuffixed variables, that refer to the single-mode case.

Based on the above relationships, the equivalent mixed-mode energy release rate $Y_m$ can be expressed as

$$Y_m = \frac{1}{2} k_n^+ \delta^2 \qquad (13)$$

where $\delta$ is an equivalent opening displacement given by

$$\delta = (\langle [u_n] \rangle_+^2 + \alpha^2 [u_s]^2)^{1/2} \qquad (14)$$
being

$$ \alpha = \sqrt{\frac{k_s}{k_n^+}} \qquad (15) $$

a mixed-mode parameter $\beta$ can be defined as

$$ \beta = \alpha \tan(\varphi) \qquad (16) $$

where $\varphi$ is the loading angle

$$ \varphi = \arctan \left[ \frac{[u_s]}{\langle [u_n] \rangle_+} \right] \in [0, +\pi/2] \qquad (17) $$

The expressions of the pure-mode contributions (12) follow as

$$ \left\{ \begin{array}{l} Y_I = \frac{1}{1 + \beta^2} Y_m \\ Y_{II} = \frac{\beta^2}{1 + \beta^2} Y_m \end{array} \right. \qquad (18) $$

The cohesive relationship can thus be reformulated as

$$ t_{\delta} = (1 - D) k_{n}^{+} \delta \qquad (19) $$

where $t_{\delta}$ is an equivalent scalar traction

$$ t_{\delta} = \left( t_{n}^{2} + \frac{1}{\alpha^{2}} t_{s}^{2} \right)^{1/2} \qquad (20) $$

being the normal and sliding components of the traction vector

$$ \left\{ \begin{array}{l} t_n = \langle \mathbf{t} \cdot \mathbf{n} \rangle_+ = \frac{1}{(1 + \beta^2)^{1/2}} t_{\delta} \\ t_s = \mathbf{t} \cdot \mathbf{s} = \frac{\alpha \beta}{(1 + \beta^2)^{1/2}} t_{\delta} \end{array} \right. \qquad (21) $$

Having identified the damage-driving force as the mixed-mode energy release rate (13), one can specify the evolution equations as

$$ \dot{D} = \dot{\gamma} \frac{\partial \phi_m}{\partial Y_m} \qquad (22) $$

along with the KKT conditions:

$$ \phi_m \le 0; \quad \dot{\gamma} \ge 0; \quad \dot{\gamma} \phi_m = 0 \qquad (23) $$

for the damage mode:
$$ \phi_m = Y_m - Y_m^* \le 0 \qquad (24) $$

where, analogous to the one-dimensional case, $Y_m^*$ denotes the mixed-mode instantaneous critical energy release rate, whose evolution is governed by a monotonically increasing positive function $F_m$.

As opposite to the one-dimensional situation, where the damage onset is determined by comparing the energy release rate with the initial pure-mode threshold $G_o$, under mixed-mode loading damage can occur before any single-mode component attains its initial allowable value. Accordingly, the definition of the critical damage-driving force for a regularized mixed-mode model, that at least formally can be given as in the single-mode case, i.e. as

$$ \left\{ \begin{array}{l} Y_m^* = Y_{mo} \quad \Rightarrow \quad D = 0 \\ Y_m^* = F_m(D) \quad \Rightarrow \quad D \in ]0, 1[ \\ Y_m^* = \max_{\tau \le t} Y_m(\tau) \quad \Rightarrow \quad D = 1 \end{array} \right. \qquad (25) $$

The damage onset is obtained according to the following criteria:

$$ \left(\frac{Y_I}{G_{oI}}\right)^{\alpha_1} + \left(\frac{Y_{II}}{G_{oII}}\right)^{\alpha_2} = 1 \qquad (26) $$

where $G_{oI}$ and $G_{oII}$ are the initial pure-mode damage thresholds while $\alpha_1$ and $\alpha_2$ are model parameters to be chosen in accordance with experimental data, that are assumed to be both strictly positive and non-necessarily integer.

For $\alpha_1 = \alpha_2$ the initial mixed-mode threshold $Y_{mo}$ can be computed:

$$ Y_{mo} = \frac{(1 + \beta^2) G_{oI} G_{oII}}{[(G_{oII})^{\alpha_1} + (\beta^2 G_{oI})^{\alpha_1}]^{1/\alpha_1}} \qquad (27) $$

For the delamination propagation, the well-known ellipse criterion is assumed [9].

$$ \left(\frac{G_I}{G_{cI}}\right)^{\beta_1} + \left(\frac{G_{II}}{G_{cII}}\right)^{\beta_2} = 1 \qquad (28) $$

where the exponents $\beta_1$ and $\beta_2$ are strictly positive reals while the mode I and mode II released energies are given by

$$ G_i = \int_0^{+\infty} Y_i \dot{D} dt; \quad i \in \{I, II\} \qquad (29) $$

For the particular case of $\beta_1 = \beta_2$ the propagation of decohesion takes place for:

$$ G_{Tc} = \frac{(1 + \beta^2) G_{cI} G_{cII}}{[(G_{cII})^{\beta_1} + (\beta^2 G_{cI})^{\beta_1}]^{1/\beta_1}} \qquad (30) $$

where $G_{Tc}$ is computed as the total work of separation:

$$ G_{Tc} = \int_{0}^{+\infty} Y_{m}^{*} \dot{D} dt \qquad (31) $$
whose expression depends upon that of the function $F_m$ defining the critical damage-driving force in the range $D \in ]0, 1[$.

In particular, taking for $F_m$ one of the expressions used in the one-dimensional case, i.e. :

$$F_m(D) = Y_{mo} + (Y_{mf} - Y_{mo})[-\log(1 - D)]^N \quad (32)$$

According to the damage evolution law, one has the expression of the parameter $Y_{mf}$ as:

$$Y_{mf} = Y_{mo} + \frac{1}{\Gamma(N + 1)} [G_{Tc} - Y_{mo}] \quad (33)$$

where $\Gamma$ is the Gamma function [5].

One can see that the interface model takes into account the modification of the mixed mode ratio during the loading path. Figure 3 presents the behaviour of the model for mixed mode.

This model has been implemented in the Finite Element Code CAST3M, where it can be used for simulation of damage evolution in adhesively bonded joints.

<figure><img src="image_4.png" /><figcaption>Figure 3 Mixed-mode traction-separation relationships.</figcaption></figure>

## 3 MECHANICAL TESTS

The parameters of the interface model are the undamaged stiffnesses ($k_n$ and $k_s$), the activation energies for each pure mode ($G_{oI}$ and $G_{oII}$), the critical energies ($G_{cI}$ and $G_{cII}$) and the exponent of the ellipse criterion for activation and propagation ($\alpha$). The exponent $\alpha$ is classicaly set to 2.

The stiffnesses of a thin layer of adhesive can not be derived from the elastic properties of the adhesive itself. They can not be identified from mechanical tests on adhesively bonded assemblies as they have a small influence on the global response of the assembly [6]. $k_n$ and $k_s$ are identified from acoustical tests not presented in this paper [11].
The activation energies $G_{oi}$ and the critical energies $G_{ci}$ can be identified straight from classical crack propagation test results.

The tests depend on the load application mode used to propagate the crack. The double-cantilever beam (DCB) and the end-notched flexure (ENF) are pure mode I and pure mode II tests, respectively. We can also have mixed-mode tests like the mixed-mode flexure (MMF). These test are presented schematically in Fig. 4.

<figure><img src="image_5.png" /><figcaption>Figure 4 Crack propagation tests</figcaption></figure>

Due to their boundary conditions simpler than mode I tests, ENF and MMF tests were performed in this work.

The samples tested consist of two 3mm thick and 20mm wide aluminum plates bonded with a layer of 0.5mm of epoxy. They were tested using a traction/compression machine (MTS 816). The tests had been perfomed via displacement control with a three-point bending fixture with span L = 120mm. The elastic properties of the aluminum are E = 75000MPa and ν = 0.3. Figure 5 shows two results of ENF tests for two different initial crack lengths a.

Figure 6 shows some results of this mixed-mode tests. On the different curves, inclinations at the beginning of each curve correspond to different initial crack lengths a. The results show that the structure compliance depends on the length of the initial crack as expected.

<figure><img src="image_6.png" /><figcaption>Figure 5 ENF results</figcaption></figure>
MMF

<figure><img src="image_7.png" /><figcaption>Figure 6 MMF results</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7 ENF test</figcaption></figure>

## 4 NUMERICAL RESULTS

Numerical simulations and tests results must be compared to evaluate the parameters of the mechanical interface model. Owing to the softening behaviour of the local traction-separation relationship, the structural response will unavoidably suffer from some mesh-dependence. This problem has been reported by many analysts and the general recommendation is that of using a sufficiently refined mesh around the decohesion front in order to allow values of the peak stress of the local traction-separation relationship that are high enough to correctly predict the decohesion [1]. The finite-element (FE) mesh used to simulate an ENF test showed in Fig. 7 is composed of 528 quadratic elements with eight nodes (3 elements in the thickness of each plate) and 74 quadratic elements of interface.

Figures 8 and 9 give the comparison between experimental curve and simulation result after identification of the damage parameters in mode II $G_{oII}$ (and) $G_{cII}$. Just after the start of the crack propagation, effects of dynamic propagation not represented in this quasi-static model do not allow the correct representation of the structure answer. Besides, the computations have
been carried out by using a varying step size and a local-control-based arc-length algorithm in order to make convergence easier by taking into account the snap-back problem [2].

To simulate a MMF test, the FE mesh showed in Fig. 7 is used. It is composed by 504 quadratic elements with eight nodes (3 elements in the thickness of each plate) and 77 quadratic elements of interface.

<figure><img src="image_9.png" /><figcaption>Figure 8 ENF numerical and test results</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9 ENF numerical and test results</figcaption></figure>
Figures 11, 12 and 13 give the comparison between a experimental curve and the result of the numerical model after identification of the rest of the damage parameters ($G_{oI}$ and $G_{cI}$). Just after the start of the crack propagation, effects of great displacements of the lower plate not represented in this model do not allow the correct representation of the structure answer.

At least, the elastic characteristics of the bonded interface that were identified in the acoustic tests and its damage characteristics that were identified in mechanical tests in mode II (ENF) and mixed-mode (MMF) are:

$$ 
\begin{array}{l@{\quad}l}
\alpha_1 = \alpha_2 = 2 & \beta_1 = \beta_2 = 2 \\
k_n = 810 \, \text{N/mm}^3 & k_s = 760 \, \text{N/mm}^3 \\
G_{cI} = 0.02 \, \text{N/mm}^3 & G_{oI} = 0.4 \times G_{cI} \\
G_{cII} = 0.09 \, \text{N/mm}^3 & G_{oII} = 0.4 \times G_{cII}
\end{array}
$$

<figure><img src="image_11.png" /><figcaption>Figure 10 MMF test</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11 MMF numerical and test results</figcaption></figure>
MMF

<figure><img src="image_13.png" /><figcaption>Figure 12 MMF numerical and test results</figcaption></figure>

MMF

<figure><img src="image_14.png" /><figcaption>Figure 13 MMF numerical and test results</figcaption></figure>
5 CONCLUSION

In the present work, a damage interface model was validated as a possible tool applied for the
simulation of crack propagation during adhesive tests.

A cohesive zone model formulated within the framework of damage mechanics was briefly
presented. The model has been implemented in a Finite Element Code and numerical simula-
tions have been carried out for some examples referring to both single-mode and mixed-mode
solicitations.

ENF and MMF crack propagation tests have been performed in aluminium/epoxy samples
with different initial crack lengths. The results show that the structure compliance depends
on the length of the initial crack as expected.

Finally, comparasions with numerical results have shown a satisfactory agreement, once
the crack propagation point could be foreseen with a quite good accuracy.

References

[1] G. Alfano and M.A. Crisfield. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. *Int. J. Num. Meth. Eng.*, 50:1701–1736, 2001. http://onlinelibrary.wiley.com/doi/10.1002/nme.93/abstract.

[2] G. Alfano and M.A. Crisfield. Solution strategies for the delamination analysis based on a combination of local-control arclength and line searches. *Int J Numer Meth Engng*, 58:999–1048, 2003. http://onlinelibrary.wiley.com/doi/10.1002/nme.806/abstract.

[3] O. Allix and A. Corigliano. Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens. *Int. J. Fracture*, 77:111–140, 1996. http://www.springerlink.com/content/v41x6g53324v2507/.

[4] O. Allix, D. Lévêque, and L. Perret. Identification and forecast of delamination in composite laminates by an interlaminar interface model. *Composite Science and Technology*, 185:671–678, 1998. http://www.sciencedirect.com/science/article/pii/S0266353897001449.

[5] G.E. Andrews, R. Askey, and R. Roy. *Special functions*. Cambridge University Press, Cambridge, 1999. http://www.amazon.com/exec/obidos/ASIN/0521623219/ref=nosim/weisstein-20.

[6] S. De Barros and L. Champaney. *Crack propagation tests : analytical and numerical approaches*. Solid Mechanics in Brazil 2009, 2009. http://www.abcm.org.br/symposiumSeries/SSSMVol2/ssSolidMechanics.shtml.

[7] J.L. Chaboche, F. Feyel, and Y. Monerie. Interface debonding models : a viscous regularization with a limited rate dependency. *Int J Solids Struct*, 38(18):3127–3160, 2001. doi:10.1016/S0020-7683(00)00053-6.

[8] Y. Mi, M.A. Crisfield, G.A.O. Davies, and H.B. Hellweg. Progressive delamination using interface elements. *J. Composite Materials*, 32(14):1246–1272, 1998. http://jcm.sagepub.com/content/32/14/1246 abstract.

[9] J.R. Reeder. *An evaluation of mixed-mode delamination failure criteria*. Technical Memorandum 104210, NASA, 1992. http://portal.acm.org/citation.cfm?id=888243.

[10] N. Valoroso and L. Champaney. A damage-mechanicsbased approach for modelling decohesion in adhesively bonded assemblies. *Engineering Fracture Mechanics*, 73:2774–2801, 2008. http://www.sciencedirect.com/science/article/pii/S0013794406001871.

[11] V. Vlasie, S. De Barros, M. Rousseau, and Champaney. Ultrasonic rheological model of cohesive and adhesive zones in aluminum joints : validation by mechanical tests. *Arch Appl Mech*, 75:220–234, 2006. http://www.springerlink.com/content/p0278447866v2240/.
<figure><img src="image_15.png" /></figure>