Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

49, 1, pp. 17-29, Warsaw 2011

BEM AND SHEAR LAG METHOD FOR INTERFACE

PROBLEM OF BI-MATERIAL STRUCTURE UNDER

STATIC LOADING

Varbinka Valeva

Jordanka Ivanova

Institute of Mechanics, Sofia, Bulgaria

e-mail: valeva@imbm.bas.bg; ivanova@imbm.bas.bg

Barbara Gambin

Institute of Fundamental Technological Research, Warsaw, Poland

e-mail: bgambin@ippt.gov.pl

The behaviour of the interface of a pre-cracked bi-material ceramic-metal
structure under static axial loading is an object of interest in the present
paper. To solve the problem for interface delamination of the structure
and to determine the debond length along the interface, a 2D BEM co-
de was created and applied. The interface plate is assumed as a very thin
plate comparing with the others two. The parametric (geometric and ela-
stic) analysis of the debond length and interface shear stress is done. First,
the obtained numerical results are compared with analytical ones from 1D
Shear lag analysis of the considered structure.The respective comparison is
illustrated in figures and shows a good agreement. The comparison betwe-
en the calculated using 2D BEM code elastic-brittle debond lengths with
Song’s experimental data for the bi-material structure Zinc/Steel as well
as with respective results from FEM simulation shows good coincidence.

Key words:BEM, Shear lag analysis, bi-material structure, debond length

1. Introduction

The boundary element method (BEM) have been demonstrated to be a
viable alternative to the FEM for many engineering problems, due to its fe-
atures of boundary-only discretization and high accuracy (Mukherjee, 1982;
Cruse, 1988; Banerjee, 1994). The high accuracy and efficiency of the BEM
for stress analysis, especially in fracture mechanics (Cruse, 1988), is very well



18 V. Valeva et al.

recognised of its semi-analytical nature. The meshing for the BEM is also
muchmore efficient than those in other domain-basedmethods, especially for
problems with changing boundaries such as crack propagation problems. Re-
cently, it was shown in Liu (1998), both analytically and numerically, that the
conventional boundary integral equation can be successfully applied to thin
structures, such as layered structures, thin films and coatings. It was shown in
Luo et al. (1998) that very accurate numerical solutions can be obtained for
thin structureswith the thickness-to-length ratio in themicro- and even nano-
scales, using the newly developed BEM approach, without seeking refinement
of the BEMmesh as the thickness decreases.

The interface strength, toughness and stiffness are important factors af-
fecting the mechanical response of multi-material layered structures. A weak
interface induces loss of structure stiffness and strength. On the other hand, a
brittle and strong interfacemay induce excessive cracking of bonded elements.
Interfacial fracture of layered composite materials under mechanical loading
was analysed innumerouspapers (see, for exampleHutchinsonandSuo, 1992).
In most papers, the analysis of thin layer cracking combined with progressive
delamination is based on assumptions of the linear fracture mechanics. Such
behaviour is treated as a mixed mode crack propagation with critical condi-
tions expressed in terms of stress intensity factors (Chiang, 1991; Lemaitre et
al., 1996; Zhang, 2000). The crack singularities at bi-material interfaces were
analysed by Hw and Hutchinson (1989), Sternitzke et al. (1996). The present
literature review is not complete as the extensive research is progressing on
multilayer and graded layer systems.

SinceCox (1952) proposeda simple one-dimensional equation for analysing
the stress transfer between a fibre and a matrix, the Shear lag approximate
analysis has become a tool for stress analysis in composite materials as well
as in layered structures. The main idea of the Shear lag analysis is such an
assumptionwhich involves a simplification of the in-plane shear stress τxy and
decouples the 2D problem into two 1D ones. Hedgepeth (1961) was the first
whoapplied the Shear lagmodel to unidirectional composites. In theShear lag
models, the hypothesis that the load is transferred from broken fibres to the
adjacent ones by thematrix shear force is stated.Hence, thematrix shear force
is independent of the transverse displacements. In Ivanova et al. (2006), Niko-
lova et al. (2007, 2009), Nikolova (2008), the Shear lag approachwas applied to
a bi-material layered structurewith a pre-cracked first thin layer. Different lo-
adingswere considered: static, thermal and combined thermo-mechanical ones.
The elastic-brittle, sleep and cohesive behaviour of the interface was assumed
and the respective lengths of delamination were found.



BEM and Shear lag method... 19

The present paper is aimed at the behaviour of the interface of bi-material
ceramic-metal plates under static axial load. The interface plate is assumed
as a very thin plate comparing with the second one and is subject only to
shear stress. To validate the application domain of the Shear lag analysis, the
problem for delamination of the interface of the biomaterial structure, a BEM
code has been created and used. The numerical model of the structure is con-
sidered in a 2Dplane-strain state. Delamination starts at the assumed restrict
condition for the value of shear stress of the interface. The obtained numerical
results are compared with analytical ones from 1D Shear lag analysis, which
can give a clear picture of the application of 1DShear lag analysis. The second
comparison between the calculated using 2DBEMcode elastic-brittle debond
length with Song’s experimental data for bi-material structure Zinc/Steel as
well as with respective results fromFEM simulation (Song et al., 2006) shows
good agreement.

2. Shear lag analysis

Consider two elastic plates A and B with finite lengths 2L, thicknesses
2hA, 2hB, bonded by an interface I and tensionally loaded with a strain ε0,
and the zero thickness interface undergoing pure shear (Fig.1). Themodified
Shear lag model will be applied (Ivanova et al., 2006), taking into account
plasticity and damage of the interface. In the shear lag model, negligence of
the bending effects results in qualitative values of the stress-strain behaviour.
Themain purpose of this study is to give simple analytical solutions, helping
the design of graded materials.

Fig. 1.

Theposedproblemconsists in the following.A crack normal to the interfa-
ce in the layer A has reached the interface by propagation inmode I.We will



20 V. Valeva et al.

study the conditions for the crack to reinitiate in the layer B.We assume that
the conditions of debonding of the interface (mode II) occur before the condi-
tions for crack reinitiation (mode I) in the layer B. The debonding length l
will be defined as the length of the interface for which the shear stresses at
the interface I reach their critical value for the interface material.

The interface is supposed to be with a negligible thickness and the shear
modulus G, shear stress τI(x,t), where the superscript I for stresses and di-
splacement belongs to the interface. It is assumed that the layers and interface
are modeled as isotropic elastic materials.

The axial stresses and strains are uniform over the cross section of each
plate, working only on tension-pressure, while the interface works on shear.
The bending is neglected.

The origin of the Cartesian coordinate system is located at the artificial
crack tip. Due to the symmetry, a half of the structure will be considered.
The stress-strain behaviour of the structure is determined by σi(x), τ

I(x),
εi(x), uj(x) (i = A,B), (j = A,B,I), where by the subscript I we denote
the interface displacement uI(x). The superscript (e) for σi(x), εi(x), uj(x)
denotes that the elastic law is assumed to describe the interface behaviour.

The following ordinary equilibrium differential equations hold

dσeA
dx
=
τI

2hA

dσeB
dx
=−

τI

2hB
(2.1)

together with respective boundary conditions and constitutive equations.

In this case, the following boundary conditions and constitutive equations
for the interface are taken

σeA(0)= 0 ε
e
A(L)= ε0 ε

e
B(L)= ε0

ueB(0)= 0 u
e
I(0)=u

e
A(0)

σeA(x)=EAε
e
A(x) σ

e
B(x)=EBε

e
B(x)

τI(x)=GweI(x) w
e
I(x)=

ueA(x)−u
e
B(x)

hA+hB
=
ueI(x)

hA+hB

(2.2)

Introduce now non-dimensional parameters, as follows

(hA+hB)x=x (hA+hB)u
e
i =u

e
i EBε0σ

e
i =σ

e
i

EBε0τ
I = τI EBε0G=G ε0ε

e
i = ε

e
i

ξ=
hA
hB

η=
EA
EB

i=A,B,I

(2.3)



BEM and Shear lag method... 21

then (2.1), (2.2) become

dσeA
dx
= τI
1+ ξ

2ξ

dσeB
dx
=−τI

1+ ξ

2
(2.4)

and

ueI(0)=u
e
A(0) u

e
B(0)= 0 ε

e
A(L)= ε

e
B(L)= 1

σeA(0)= 0 σ
e
A(x)= ηε

e
A(x) σ

e
B(x)= ε

e
B(x)

τI(x)=GueI(x)

(2.5)

Equations (2.4) result in
d2ueI
dx2
=λ
2
ueI (2.6)

where

λ
2
=
G(1+ ξ)(1+ ξη)

2ξη

and (2.4) becomes

d2ueA
dx2
=
λ
2

1+ ξη
ueI

d2ueB
dx2
=−

λ
2

1+ ξη
ξηueI

The general solution to (2.6) is

ueI(x)=A1 sinh(λx)+A2cosh(λx) (2.7)

To find the stresses, strains and respective displacements in the plates, we
use equations (2.4), together with boundary conditions (2.5). In addition, the
substitution ueI(x)=u

e
A(x)−u

e
B(x) has to bemade.

We obtain the following expressions for the interfacial displacement and
shear stress in dimensionless parameters

ueI(x)=
1+ ξη

λ

cosh[λ(L−x)]
sinh(λL)

τI(x)=G
1+ ξη

λ

cosh[λ(L−x)]
sinh(λL)

(2.8)

The debond length le, which gives themagnitude of brittle cracking along the
interface can be calculated from (2.8) on the assumption that the shear stress
reaches its critical failure value τI = τcr and ueI(le)=u

cr = τcr/G. Then

τcr =G
1+ξη

λ

cosh[λ(L− le)]
sinh(λL)

(2.9)



22 V. Valeva et al.

Using the substitution exp[λ(L− le)] = y and (2.9), the following equation
for y is obtained

y2−2Ay+1=0 A=
λτcr sinh(λL)

G(1+ ξη)
(2.10)

Then two roots of (2.10) are available: y1,2 = A±
√
A2−1. Now using the

substitution exp[λ(L− le)] = y, we obtain

(le)1,2 =L−
1

λ
ln[A±

√

A2−1]

It is necessary that

A2−1=
[λτcr sinh(λL)

G(1+ ξη)

]2
−1> 0

This requirement poses some condition for the value of τcr.
Thenwe have to choose the length of the debonding zone froma condition

that this length must beminimum (Ivanova et al., 2006), i.e.

le =L−
1

λ
ln[A+

√

A2−1] (2.11)

3. BEM formulation

The following known boundary integral equations for two-dimensional elasti-
city problems can be applied in eachmaterial domain (index notation is used,
where repeated subscripts imply summation) (Mukherjee, 1982)

Cij(P0)u
(β)
j (P0)=

∫

Γ

[U
(β)
ij (P,P0)t

(β)
j (P)−T

(β)
ij (P,P0)u

(β)
j (P)] dΓ(P) (3.1)

in which u
(β)
j and t

(β)
j are the displacement and traction fields, respectively;

U
(β)
ij (P,P0) and T

(β)
ij (P,P0) the displacement and traction kernels (Kelvin’s

solution), respectively; P is the field point and P0 – the source point, and Γ –
the boundary of the single domain. Cij(P0) is a constant coefficient matrix
depending on the smoothness of the boundary Γ at the source point P0. The
superscript β on the variables in Eq. (3.1) signifies the dependence of these
variables on the individual domains β=A,B,I.



BEM and Shear lag method... 23

The two kernel functions U
(β)
ij (P,P0) and T

(β)
ij (P,P0) in boundary integral

equation (3.1) for plain-strain problems are given as follows

U
(β)
ij (P,P0)=

1

8πµ(β)(1−ν(β))

[

(3−4ν)δij ln
(1

r

)

+r,ir,j
]

(3.2)

T
(β)
ij (P,P0)=−

1

4πr(1−ν(β))
·

·
{

r,n
[(

1−2ν(β)
)

δij +2r,i r,j
]

+
(

1−2ν(β)
)

(r,jni−r,inj)
}

where µ(β) is the shear modulus and ν(β) Poisson’s ratio for three diffe-
rent domains, respectively; r is the distance from the source point P0 to
the field point P; ni is the i-th directional cosine of the outward normal n;
(·),i= ∂(·)/∂xi with xi being the coordinate of the field point P; and δij is
the Kronecker delta.
In Eq. (3.1), the integral containing the U

(β)
ij (P,P0) kernel is weakly sin-

gular, while the one containing T
(β)
ij (P,P0) is strongly singular and must be

interpreted in the Cauchy principal value sense. However, when the structure
becomes thin in shape, such as the interphase, both integrals are difficult to
dealwithwhen the sourcepoint is on theone sideand the integration is carried
out on the other side of the thin structure. These types of integrals are called
nearly singular integrals since the distance r is very small in this case but
is still not zero. Recently, several techniques, including singularity subtrac-
tions, analytical integration, and nonlinear coordinate transformations have
been developed to calculate the nearly singular integrals (Luo et al., 1998).
The combination of these techniques is found to be extremely effective and
efficient in computing the nearly singular integrals in two-dimensional boun-
dary integral equations, no matter how close the source point to the element
of integration is.
The discretization of BIE (3.1) using boundary elements follows the stan-

dard BEM procedures except for the nearly-singular integrals. For multi-
domain (material) problems, the resulting BEM equations for each material
domain are coupled together by the interface conditions (continuity of both
displacements and equilibriumof both tractions) and then solved to obtain the
displacement and traction vectors at each node on the boundary and interfa-
ces. In the BEM approach used in the present paper, for solving the nearly
singular integrals, subdivision of the element of integration and an adaptive
integration scheme are proposed. For the integration of a logarithmic function,
amodifiedGaussQuadrature (Gauss-Laguerre) is used. For the casewhen the
collocation point is located at the element or is very close to the element, the



24 V. Valeva et al.

integration of the nonsingular part (where the shape function value is zero at
the collocation point) and the singular part (for weakly singular behaviour of
the kernel) is performed separately. For the corners, discontinuous elements
are used.

4. Numerical example

The first numerical example is a comparison between the BEMand Shear lag
model results.

On the basis of the obtained analytical formulae for the assumed interface
shear stress law, the stress behaviour (especially the debond length on the
interface) of the two-plate structure with different mechanical and geometric
properties under tension ε0 will be studied. Using BEM, the bending of the
structure is avoidedby the imposedboundarycondition uB(x,−(2hB+t))= 0,
where t is the thickness of the interface I.

The following geometric andmechanical properties (Table 1) are used:

2L=24mm, 2hB =6
′,mm

{

2hA =2mm, (ξ=1/3)

2hA =1mm, (ξ=1/6)

τcr =18MPa, a=1mm, t=0.1mm, ε0 ∈ [0.001,0.008]

Table 1.Characteristics of the materials (Sternitzke et al., 1996)

Material
E ν
[GPa] [–]

Layer A
C84 [Al2O3/Al composites,

285 0.28
(C84=84vol%Al2O3+16vol%Al)]

Layer A Alumina [DEGUSSIT Al 23, Friatec.] 380 0.24

Layer B 100Cr6 [AISI 52100] 210 0.29

Interphase Polyacrylate thermoplastic glue 2.5 0.50

In Fig.2, a comparison between 1D Shear lag and BEM 2D interface de-
bond length predictions is shown. The values of debond length le versus the
applied load ε0 are obtained for two different ratios η of the elastic moduli
and for two different values of the thickness ratio ξ.

It can be seen that geometric characteristics much more influence the de-
bond length than thematerial characteristics. Considering the numerical and



BEM and Shear lag method... 25

analytical results (BEM, Shear lag), the bigger is the thickness ratio ξ, the
smaller is the value of applied load εcr0 needed for full delamination (degra-
dation) of the interface. The critical load εcr0 calculated using Shear lag at
different thickness ratios ξ is much smaller comparing with εcr0 obtained by
BEM. This difference for εcr0 can be explained with the presence of a normal
crack, which strongly reflects on the stress-strain behaviour (BEM) of the first
plate. On the other hand, the very thin first pre-cracked plate plays a signifi-
cant role in full degradation of the bi-material structure, allowing for a bigger
critical load.

Fig. 2. Comparison and parametric analysis (geometric and elastic properties)
between the BEM and Shear lag model for the debond length

The relative error of comparison between the BEM and Shear lag values
of the debond length

r=

√

√

√

√

1

M

M
∑

i=1

y
theory
i

yBEMi
·100%

is:



26 V. Valeva et al.

C84/100Cr6

ξ=1/3 η=1.36 r=10.75%

ξ=1/6 η=1.36 r=13.01%

Alumina/100Cr6

ξ=1/3 η=1.81 r=9.22%

ξ=1/6 η=1.81 r=11.64%

The investigations showthat thebigger is the load, thebigger is the relative
error. The decreasing of the thickness of the first plate ξ at a constant elastic
ratio η also leads to increment of the relative error. The increasing of the
value of the elastic ratio η at a constant value of the geometric ratio ξ leads
to decrement of the relative error. It is a consequence of negligible thickness
of the interface as well as the fact that the approximate analytical Shear lag
model is 1D.

InFig.3, thenumericalBEMresults for stressesof thebi-material structure
for ξ = 1/3, 1/6 and η = 1.81 are shown (as an example). The loading is
ε0 =0.0025.

Fig. 3. Plots of stresses σxx(x,y), σxy(x,y) for the pre-cracked bimaterial structure
for different thickness ratios ξ=hA/hB and the elastic moduli ratio η=1.81

The secondnumerical example is the comparisonbetween the experimental
data and FEM results for progressive elastic-brittle interface debond lengths
of zinc coatings on a steel substrate (Song et al., 2006) with our BEM results
for respective values of debond lengths. The calculations are performed for the
following geometrical andmechanical properties of the zinc coating (layer A)
and steel substrate (layer B): 2L = 100µm, 2hA = 10µm, 2hB = 50µm,
EA =70GPa, νA =0.3, EB =200GPa, νA =0.27. The delamination along
the zinc coating/steel substrate interface is simulated by deleting the elements



BEM and Shear lag method... 27

for which the shear stress is larger than the critical value representing the
interface strength – this is the criterion for crack initiation and propagation
along the zinc coating/steel interface with increasing tensile load. Figure 4
shows the interface debond length le as a function of the applied strain ε0,
comparing the experimental and calculated by FEM and BEM data for the
given interface shear strength of 180MPa.TheBEMresults for debond length
are in very good coincidence with the experimental data and FEM results by
Song et al. (2006).

Fig. 4. Measured average interface debond length le as a function of the applied
strain ε0 versus calculated interface debond length for the zinc coating shear

strength τcr =180MPa

5. Conclusions

In the paper, a comparison between the approximate Shear lag 1D method
and 2D BEM for interface delamination of the bi-material structure under
static load is done.The relative error between analytical and numerical results
confirms validity of the Shear lag approach. The obtained predictions can
be applied to a pre-cracked by an indentor bi-material structures undergoing
static tension for different mechanical behaviour andmaterials of plates.

Acknowledgement

The authors acknowledge the support from the bi-lateral project of BAS/PAS

”New composite materials, homogenization and macroscopic behaviour of structural

elements” 2009-2011.



28 V. Valeva et al.

References

1. Banerjee P.K., 1994, The Boundary Element Methods in Engineering, 2nd
ed., McGraw-Hill, NewYork

2. Chiang R., 1991, On the stress intensity factors of crack near an interface
between twomedia, Int. J. Fracture, 47, R55-R58

3. Cox L.H., 1952, The elasticity and strength of paper and other fibrous mate-
rials,Brit. J. Appl. Phys., 72, 3

4. Cruse T.A., 1988, Boundary Element Analysis in Computational Fracture
Mechanics, Kluwer Academic Publishers, Dordrecht, The Netherlands

5. He M.Y., Hutchinson J.W., 1989, Crack deflection at an interface between
dissimilar elastic materials, Int. J. Solids Structures, 25, 9, 1053-1067

6. Hedgepeth J.M., 1961,StressConcentrations in Filamentary Structures,NA-
SATND-882

7. Hutchinson J.W., Suo Z., 1992,Mixedmode cracking in layeredmaterials,
Adv. Appl. Mech., 29, 62-191

8. Ivanova J., Valeva V., Mroz Z., 2006, Mechanical modelling of the de-
lamination of bi-material plate structure, Journal of Theoretical and Applied
Mechanics, Sofia, 36, 4, 39-54

9. Lemaitre J., Desmorat R., Vidonne M.P., Zhang P., 1996, Reinitiation
of a crack reaching an interface, Int. J. Fracture, 80, 257-276

10. Liu Y.J., 1998, Analysis of shell-like structures by the boundary element me-
thod based on 3-D elasticity: formulation and verification, Int. J. Numer. Me-
thods Engrg., 41, 541-558

11. Luo J.F., Liu Y.J., Berger E.J., 1998, Analysis of two-dimensional thin
structures (from micro- to nano-scales) using the boundary element method,
Comput. Mechanics, 24, 404-412

12. Mukherjee S., 1982,Boundary Element Methods in Creep and Fracture, Ap-
plied Science Publisher, NewYork

13. Nikolova G., 2008, Thermo-Mechanical Behaviour of Thin Graded Layered
Structures, PhDThesis

14. Nikolova G., Ivanova J., 2009, Cracked biomaterial plates under thermo-
mechanical loading,Proceedings Fractography of AdvancedCeramics III, Trans.
Tech. Publications Ltd. In Key Material Series, 409, 406-413

15. NikolovaG., IvanovaJ., ValevaV.,MrozZ., 2007,Mechanical and ther-
mal loading of two-plate structure, Comptes Rendus De L’Academie Bulgare
Des Sciences, 60, 7, pp. 735



BEM and Shear lag method... 29

16. Song G.M., Sloof W.G., Pei Y.T., De Hosson J.Th.M., 2006, Interface
behaviour of zinc coating on steel: Experiments and finite element calculations,
Surface and Coating Technology, 201, 4311-4316

17. Sternitzke M., Knechtel M., Hoffman M., Broszeit E., Rödel J.,
1996,Wear properties of alumina/aluminum composites with interpenetrating
networks, J. Am. Ceram. Soc., 79, 1, 121-128

18. Zhang S., 2000, Thermal stress intensities at an interface crack between two
elastic layers, Int. J. Fract., 106, 277-290

Metoda Elementów Brzegowych i metoda ”shear lag” w zagadnieniu

delaminacji dwuwarstwowej struktury pod wpływem statycznego

obciążenia

Streszczenie

W pracy badano zachowanie się na granicy pomiędzy warstwami metalu i cera-
miki pod wpływem statycznego obciążenia przyłożonegow kierunku równoległymdo
połączenia pasm, w przypadku istnienia początkowego nacięcia w jednej z warstw
prostopadłego do powierzchni połączenia. W celu rozwiązania zagadnienia delami-
nacji wzdłuż powierzchni łączącej oba materiały i wyznaczenia długości odspojenia
został stworzony i zastosowany kod Metody Elementów Brzegowych – zagadnienie
2-wymiarowe.Warstwa łączącadwamateriały zostałapotraktowana jakobardzo cien-
ka płyta, w porównaniu do grubości obuwarstwmateriałowych. Przeprowadzonopa-
rametryczną (geometrycznaą i sprężystą) analizę długości odspojenia (delaminacji)
i naprężenia stycznego. Otrzymane rezultaty numeryczne porównano z analityczny-
mi rozwiązaniami 1-wymiarowej analizy tzw.metodą ”shear lag”.Otrzymanewyniki,
zilustrowane na rysunkach, wykazują wzajemną zgodność. Pokazano również, że wy-
niki uzyskane przy użyciu Metody Elementów Brzegowych są zgodne z wynikami
eksperymentu Song’a przeprowadzonego dla dwuwarstwowego elementu kompozyto-
wego pomiędzy warstwą cynku i stali, a także z wynikami uzyskanymi w Metodzie
Elementów Skończonych.

Manuscript received October 21, 2009; accepted for print March 15, 2010