Wave propagation analysis in adhesively bonded composite joints using the wavelet spectral finite element method

Composite Structures 122 (2015) 271–283

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Wave propagation analysis in adhesively bonded composite joints using the wavelet spectral finite element method Dulip Samaratunga a,⇑, Ratneshwar Jha a, S. Gopalakrishnan b a b

Raspet Flight Research Laboratory, Mississippi State University, 114 Airport Road, Starkville, MS 39759, USA Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India

a r t i c l e

i n f o

Article history: Available online 3 December 2014 Keywords: Wavelet spectral finite element Single lap joint Bonded beams Wave propagation

a b s t r a c t A wavelet spectral finite element (WSFE) model is developed for studying transient dynamics and wave propagation in adhesively bonded composite joints. The adherands are formulated as shear deformable beams using the first order shear deformation theory (FSDT) to obtain accurate results for high frequency wave propagation. Equations of motion governing wave motion in the bonded beams are derived using Hamilton’s principle. The adhesive layer is modeled as a line of continuously distributed tension/compression and shear springs. Daubechies compactly supported wavelet scaling functions are used to transform the governing partial differential equations from time domain to frequency domain. The dynamic stiffness matrix is derived under the spectral finite element framework relating the nodal forces and displacements in the transformed frequency domain. Time domain results for wave propagation in a lap joint are validated with conventional finite element simulations using AbaqusÒ. Frequency domain spectrum and dispersion relation results are presented and discussed. The developed WSFE model yields efficient and accurate analysis of wave propagation in adhesively-bonded composite joints. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Adhesive bonding of composite structural components has gained considerable attention during the last few decades due to several advantages over mechanical fastening. Major advantages of adhesive bonding include higher fatigue resistance and longer fatigue life, light weight, ability to join thin and dissimilar components, good sealing, low manufacturing cost, and good vibration and damping properties. In fact, adhesive bonding has already become a common practice for joining carbon fiber reinforced polymer (CFRP) in secondary aircraft structures due to these advantages [1–4]. However, adhesive bonding is yet to be adopted for application in primary structural joints. Current aircraft certification requirements mandate proof that adhesively bonded joints will maintain their integrity at critical design loads. Presently these requirements are complied using thousands of mechanical fasteners along with adhesives. The two modern day landmark largely composite commercial aircraft, Boeing 787 Dreamliner and Airbus A350 XWB, are no exception [5]. Addition of mechanical fasteners obviously hinders the realization of full cost and weight savings of

⇑ Corresponding author. Tel.: +1 662 325 5004; fax: +1 662 325 3864. E-mail addresses: [email protected] (D. Samaratunga), jha@raspet. msstate.edu (R. Jha), [email protected] (S. Gopalakrishnan). http://dx.doi.org/10.1016/j.compstruct.2014.11.053 0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.

composites as a structural material. Alternatively, primary structures with bonded joints can be certified by establishing repeatable and reliable nondestructive evaluation (NDE) or structural health monitoring (SHM) procedure to ensure joint integrity. Ultrasonic wave based methods are considered among most suitable techniques for NDE/SHM of composite structures. Adams and Cawley [6] reviewed ultrasonic and other NDE methods investigated until late 1980s for evaluating composites and bonded joint integrity. Several of the early methods were able to detect delaminations in composites and disbonds in adhesive joints, but poor joint adhesion was very difficult to detect. Later investigations focused on adhesive joint quality inspections using ultrasonic testing methods including acoustic emission, bulk and guided (longitudinal, shear and Lamb) waves, nonlinear solitary waves, etc. [7–14]. Crom and Castaings [9] presented a one-dimensional semi-analytical finite element model for predicting dispersion curves and mode shapes of shear-horizontal (SH) guided waves propagating along multi-layered plates made of anisotropic and viscoelastic materials (composite and adhesive layer). Hanneman et al. [12–13] derived the transfer function (TF) of a multilayer medium based on the exact solution of 1-D wave equation in order to evaluate the TF sensitivity with adhesive layer thickness and modulus and verified the TF through experiments. Saito and Tani [15] investigated the natural frequencies and loss factors of the

272

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

coupled longitudinal and flexural vibrations of parallel and identical cantilevers lap-joined using a viscoelastic bonding material. A complete set of equations of motion and boundary conditions governing the dynamics of the system based on the Euler–Bernoulli beam theory was derived. In addition, finite element modeling of bonded lap joints for transient response prediction is reported by a number of authors [16–18]. Structural health monitoring methods aim to perform nondestructive damage diagnosis through actuators/sensors integrated with the structure. SHM is an inverse problem where the measured data are used to predict the condition of the structure. Therefore, to be able to differentiate between damage and structural features, prior information is required about the structure in its undamaged state. This is typically in the form of a baseline data measured from the so called ‘‘healthy state’’ to be used as a reference for comparison with the test case. Alternatively, a validated physics-based model for wave propagation combined with experimental measurements can be used for complete characterization (presence, location, and severity) of damages. The modeling of wave propagation in composites presents complexities beyond that for isotropic structures. Analytical solutions for wave propagation are not available for most practical structures due to complex nature of governing differential equations and boundary/initial conditions. The finite element method (FEM) is the most popular numerical technique for modeling wave propagation phenomena. However, for accurate predictions using FEM, typically 10–20 elements should span a wavelength [19,20], which results in very large system size and enormous computational cost for wave propagation analysis at high frequencies. In addition, solving inverse problems (as required for SHM) is very difficult using FEM. Spectral finite element (SFE), which follows FEM modeling procedure in the transformed frequency domain, is highly suitable for wave propagation analysis [21–25]. SFE models are many orders smaller than FEM and highly suitable for efficient SHM. Frequency domain formulation of SFE enables direct relationship between output and input through system transfer function (frequency response function). SFE has very high computational efficiency since nodal displacements are related to nodal tractions through frequency-wave number dependent stiffness matrix. In SFE mass distribution is captured exactly and the accurate elemental dynamic stiffness matrix is derived. Subsequently, in the absence of any discontinuities, one element is sufficient to model a beam or plate structure of any length. SFE was initially developed using fast Fourier transform as the integral transformation technique. Much of the early work in this field performed by Doyle and associates is reported in Doyle [21]. Later, Gopalakrishnan and associates [22] extensively investigated Fourier Spectral Finite Element (FSFE) method for beams and plates with anisotropic and inhomogeneous material properties. Although the FSFE method is very efficient for wave motion analysis and suitable for solving inverse problems, it is not capable of modeling short waveguides. Further, for 2-D problems, FSFE are essentially semi-infinite, that is, they are bounded only in one direction and the effect of lateral boundaries cannot be captured due to the global basis functions of the Fourier series approximation used in the spatial dimension. In addition to that, the required assumption of periodicity in time approximation results in ‘‘wraparound’’ problem for smaller time window distorting the response. The wavelet spectral finite element (WSFE) developed by Gopalakrishnan and Mitra [23] overcomes these shortcomings by using Daubechies compactly supported wavelet scaling functions as the basis functions. The WSFE method is highly efficient and suitable for wave propagation analysis in short waveguides as well as finite plate structures. The present authors developed a new 2-D WSFE model (published in this Journal, [24]) based on the first order shear deformation theory (FSDT) which yields accurate

results for wave motion at high frequencies (unlike the previous classical laminate theory based model, [23]). The present authors also modeled a 2-D laminate with transverse crack wherein separate displacement fields were assumed on either side of the crack (penetrating part-through the thickness) and the crack was represented as a continuous line spring [25]. To the authors’ knowledge, no model has been developed so far to simulate adhesively bonded composite joints using the spectral finite element method. In this work, we report a newly developed wavelet spectral finite element model for studying transient dynamics and wave propagation across adhesively bonded composite joints. The model is highly efficient computationally and can be utilized in solving inverse problems such as damage detection and force reconstruction. Analytical equations of motion governing the bonded double beam system are based on the first order shear deformation theory and derived using the dynamic version of the principle of virtual work (Hamilton’s principle). The adhesive layer is modeled as a line of continuously distributed tension/compression and shear springs. Daubechies compactly supported wavelet scaling functions are used to transform the governing partial differential equations (PDEs) into ordinary differential equations (ODEs). The ODEs are solved exactly for wavenumbers and the dynamics stiffness matrix is derived relating the nodal forces and displacements in the transformed domain. Frequency domain spectrum and dispersion relation results are presented and discussed, including the influence of bondline degradation on wave modes. The bonded double beam WSFE model is then used to investigate wave propagation in adhesively bonded single lap joints. Time domain wave propagation results are compared with conventional finite element simulations using AbaqusÒ. Concluding remarks are given in the final section. 2. Mathematical formulation for bonded single lap joint A schematic of the single lap joint (SLJ) studied in this work is shown in Fig. 1. The two beam adherands with sections lying on the X-Y plane are bonded using adhesive along the overlap zone (shown in gray in Fig. 1). Mathematical formulation starts with the derivation of equations of motion (EOM) for each adherand.

Fig. 1. Bonded single lap joint.

Fig. 2. Double beam system with stress resultants and interaction forces.

273

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

The equations of motion of the adherands inside the overlap zone (the area between the dashed lines I and II in Fig. 1) are first derived using Hamilton’s principle. The adherand portion outside the overlap area can be treated as a beam free of external interactions and EOM are reduced from the formulation for overlapped region by relieving adhesive effect. 2.1. Equations for wave motion in adhesively bonded double beam system The free body diagram of the double beam system (representing the overlap zone of the single lap joint in Fig. 1) is shown in Fig. 2. Normal, shear and bending stress resultants are denoted by Ni, Qi and Mi, respectively, for each of the beams numbered by i (=1, 2). Normal and shear interaction forces between the adherands and the adhesive layer are given by p and s measured in force per unit length. Transverse shear in the adhesive layer is assumed to be negligible since its thickness is much smaller than the adherands. The cross sections of the adherands are considered to be rectangular with area A = 2hib. In a beam, the displacement fields for axial and transverse motion based on FSDT [26] are given by,

uðx; y; tÞ ¼ uc ðx; tÞ  y/ðx; tÞ

v ðx; y; tÞ ¼ v c ðx; tÞ

ð1Þ

8 8 9 9 > Z h> < Nxx > < rxx > = = yrxx dy; M xx ¼ b > > > > h : : ; rxy ; Qx

Here L and b are the length and width of the beam, respectively. Integration by parts is used to relieve variations from any differentiation and thus Eq. (4) can be written as

dU ¼ 

dU ¼

rij deij dV

ð3Þ

V

where d is the variational operator [27] and rij are stress terms corresponding to strain terms in Eq. (2). Substituting strains of Eq. (2) into Eq. (3), we get

dU ¼

Z 0

¼

Z

0

L

Z

L

A

  @N xx @M xx @Q x duc  d/ þ dv c þ Q x d/ dx @x @x @x

þ ½Nxx duc L0  ½Mxx d/L0 þ ½Q x dv c L0

ð6Þ

In the absence of body forces, virtual work due to external tractions, ^t; causing virtual displacement, du, over the portion S2 of total boundary S can be written as

dV ¼ 

Z

^t dudS

ð7Þ

S2

The external forces on the top adherand are s and p, and the virtual work due to external forces can be written as

dV ¼ 

Z

L

½sðduc þ d/hÞ  pdv c dx

ð8Þ

The virtual kinetic energy of the body is defined as

dK ¼

Z

q

V

@u @du dV  @t @t

ð9Þ

where q is the mass density. Upon substitution of the displacement field given by Eq. (1) into Eq. (9), the virtual kinetic energy, dK, can be written as

dK ¼

Z 0

L

Z

h

h

    @uc @/ @duc @d/ @ v c @dv c þ dydx y y @t @t @t @t @t @t

qb

ð10Þ This expression also contains terms with derivatives of variations. Integration by parts is used to obtain the final form of virtual kinetic energy.

# Z L" 2 @ uc @ 2 uc @2/ @2/ @2v c I0 2 duc  I1 2 d/  I1 2 duc þ I2 2 d/ þ I0 2 dv c dx @t @t @t @t @t 0  L  L  L  L @uc @uc @/ @/ duc  I1 d/  I1 duc þ I2 d/ þ I0 @t @t @t @t 0 0 0 0  L @v c þ I0 dv c @t 0

dK ¼ 

ð11Þ where the inertial terms are defined as

8 9 8 9 > Z h> < I0 > = <1> = I1 ¼ qb y dy; > h > : > ; : 2> ; y I2

ð12Þ

The inertial term I1 would be zero when the mid plane is taken as the reference. According to Hamilton’s principle the relationship of kinetic and total potential energies for a conservative elastic body can be written as

Z

½rxx dexx þ rxy dcxy dAdx

 @duc @d/ @dv c Nxx  Mxx þ Qx  Q x d/ dx @x @x @x

L

ð2Þ

Total potential energy is the sum of the strain energy and potential energy due to internal and external forces. For a body whose volume is denoted by V, the virtual strain energy can be written as

Z

Z 0

exx ¼ eyy

ð5Þ

0

where u and v are axial and transverse displacements, respectively, at any material point on the beam. The terms uc and vc are the beam axial and transverse displacements along the reference (mid) plane. Anticlockwise rotation of the beam cross-section about Z-axis is given by /. Note that the subscript of the terms denoting adherand number (appearing in Fig. 2) has been omitted for simplicity. The derivation of EOM for each of the adherands of the above double beam system starts by expressing the strain in terms of the displacement field. Once the strains are determined, then virtual kinetic and total potential energies are calculated. According to the dynamic version of the principle of virtual work (also called Hamilton’s principle) the variation of line integral of the kinetic and total potential energy difference, also known as the Lagrangian, between two arbitrary time instants is equal to zero for a conservative body [27]. Thus the Euler–Lagrange EOM for beam adherands, including the effect of adhesive layer, can be derived. Even though one can use the general nonlinear strain–displacement relations, the present work restricts the development to small strains and displacements for simplicity since our main focus is WSFE modeling. The linear strains associated with the above displacement field (Eq. (1)) are

@u @uc @/ @ v @u @ v c ¼ y ; cxy ¼ þ ¼ / @x @x @x @x @y @x ¼ ezz ¼ cxz ¼ cyz ¼ 0

where the force resultants are defined as

t2

d

½K  ðU þ VÞdt ¼ 0

ð13Þ

t1

ð4Þ

where t1 and t2 are two arbitrary time instants. By substituting Eqs. (6), (8) and (11) into Eq. (13) and rearranging the terms, we get

274

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

Z t2 (Z L "

!

! @Nxx @ 2 uc @2/ @Q x @2v c þ s  I0 2 þ I1 2 duc þ  p  I 0 2 dv c @x @x @t @t @t 0 t1 ! #  @M xx @ 2 uc @2/ þ  þ Q x þhs þ I1 2 I2 2 d/ dx @x @t @t ) ½N xx duc L0 þ ½M xx d/L0  ½Q x dv c L0 dt

ð14Þ

The terms obtained in volume V but evaluated at t1 and t2 were set to zero because the virtual displacements are zero there (by assumption). The Euler–Lagrange equations for the top adherand based on FSDT are obtained by setting the coefficients of duc ; dv c and d/ in Eq. (14) to zero:

du1c :

@N xx1 @x

dv 1c :

@Q x1 @x

2 I1 @@t/21

þs

2 I0 @@tu21c

þ

p

2 I0 @ @tv21c

¼0

¼0 ð15Þ

ð21Þ The boundary conditions can be written as

Mxx2

M xx1 ¼ M xx1

Np X ðkÞ  A11 ¼ b Q 11 ykþ1  yk ;

D11 ð16Þ

Similar procedure is followed to obtain EOM and the boundary conditions for the bottom adherand. The final form of the EOM for the bottom adherand is 2

2

du2c :

@N xx2 @x

 s  I0 @@tu22c þ I1 @@t/22 ¼ 0

dv 2c :

@Q x2 @x

þ p  I0 @ @tv22c ¼ 0

2

ð17Þ 2

2

d/2 :  @M@xxx2 þ Q x2 þ h2 s þ I1 @@tu22c  I2 @@t/22 ¼ 0 and the natural boundary conditions for the bottom adherand are written as

Nxx2 ¼ Nxx2 ;

Q x2 ¼ Q x2 ;

M xx2 ¼ M xx2

ð18Þ

Subscripts denoting top and bottom adherands are included in Eqs. (15)–(18) for clarity. Since the development so far did not utilize the constitutive equations, Eqs. (15) and (17) can be tailored to address nonlinear elastic body problems by incorporating nonlinear strain terms in Eq. (2). However, present work is focused on linear aspect of the problem, as noted earlier. Next, force and moment resultants given in Eq. (5) are related to the strains of the laminate through constitutive equations. To this end, it is assumed that each layer is orthotropic with respect to its material symmetry lines and obeys the Hook’s law. Then the strain expressions given in Eq. (2) are used to express the above EOM and boundary conditions (Eqs. (15)–(18)) in displacement terms. The EOM of the top adherand expressed in displacement terms can be written as 2

2

2

2

A11 @@xu21c  B11 @@x/21 þ s  I0 @@tu21c þ I1 @@t/21 ¼ 0 2  2 jA55 @@xv21c  @/@x1  p  I0 @ @tv21c ¼ 0  v 1c 2 2 2 2 B11 @@xu21c þ D11 @@x/21 þ A55 @@x  /1 þ h1 s þ I1 @@tu21c  I2 @@t/21 ¼ 0 ð19Þ The term j is the shear correction factor which is introduced to mitigate the constant transverse strain state assumed by the FSDT [26]. The boundary conditions can be written as

  @u1c @/ @ v 1c  B11 1 ; Q x1 ¼ jA55  /1 ; @x @x @x @u1c @/1 ¼ B11 þ D11 @x @x

Nxx1 ¼ A11 M xx1

ð22Þ

The axial, axial-bending, and bending stiffness terms appearing in the formulation are defined as

B11 ¼

k¼1

Eq. (14) contains the boundary condition terms given by

Q x1 ¼ Q x1 ;

  @u2c @/ @ v 2c  B11 2 ; Q x2 ¼ jA55  /2 ; @x @x @x @u2c @/2 ¼ B11 þ D11 @x @x

Nxx2 ¼ A11

d/1 :  @M@xxx1 þ Q x1 þ h1 s þ I1 @@tu21c  I2 @@t/21 ¼ 0

Nxx1 ¼ Nxx1 ;

2

2

 v 2c 2 2 2 2  /2 þ h2 s þ I1 @@tu22c  I2 @@t/22 ¼ 0 B11 @@xu22c þ D11 @@x/22 þ A55 @@x

2

2

2

2

A11 @@xu22c  B11 @@x/22  s  I0 @@tu22c þ I1 @@t/22 ¼ 0 2  2 jA55 @@xv22c  @/@x2 þ p  I0 @@tv22c ¼ 0

Np 1 X ðkÞ  b Q 11 y2kþ1  y2k ; 2 k¼1

Np 1 X ðkÞ  ¼ b Q 11 y3kþ1  y3k 3 k¼1

ð23Þ

ðkÞ

where Q 11 represents the stiffness of the kth ply of the laminate [26] and Np is the total number of plies. The two sets of EOM given by Eqs. (19) and (21) are coupled due to the presence of normal and shear forces (p and s) at the adhesive and adherand interfaces. By setting these adhesive-adherand interaction forces to zero in one set of the equations, EOM for a simple beam (that is, in the absence of distributed external forces) can be recovered. These forces can be derived through constitutive relations of the adhesive layer. The constitutive relations for the adhesive layer are established using two parameter elastic foundation approach, where the adhesive layer is assumed to be composed of continuously distributed shear and tension/compression springs. Similar approach has been taken by Mortensen [28] for solving static problems and by Saito and Tani [15] for studying coupled longitudinal and flexural vibrations of adhesively bonded joints. From Fig. 2, the interaction forces p and s are found to be

s ¼ ks ðu2c  h2 /2  u1c  h1 /1 Þ; p ¼ kt ðv 1c  v 2c Þ;

kt ¼

ks ¼

Ga ; ta

Ea ta

ð24Þ

where Ea and Ga are Young’s and shear moduli, respectively. The terms ks and kt are spring constants into X and Y directions. This completes the derivation of partial differential equations for the bonded double beam system shown in Fig. 2. It should be mentioned that the procedure outlined above can be used to formulate EOMs for other bonded structures (with larger number of members), such as a double lap joint. As already mentioned, the governing differential Eqs. (19) and (21) represent a system of coupled linear partial differential equations (PDEs), which are difficult to solve exactly in the time domain for all boundary conditions. Here we employ WSFE method to solve the governing PDEs in the frequency domain. The solution for given loading and boundary conditions is obtained in the frequency domain and inverse wavelet transform is performed to obtain time domain results. The next section describes the implementation of WSFE method for solving EOMs presented above. 2.2. Wavelet approximation of governing equations

ð20Þ

Similarly, the EOM for the bottom adherand in displacement terms can be written as

The WSFE formulation begins with transformation of the field variables (displacements) into the frequency domain using the wavelet transform. Daubechies compactly supported wavelet scaling functions [29] are used for approximation in time, which

275

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

treatment of finite length data and uses polynomial to extrapolate coefficients at boundaries either from interior coefficients or boundary values. The method is particularly suitable for approximation in time for the ease of initial value imposition. After treating the boundaries using wavelet extrapolation technique, the ODEs in Eq. (27) can be written in matrix form as

A11 Fig. 3. Nodal representation of the spectral element for the bonded double beam system.

reduces the PDEs to ODEs in the spatial dimension. Compactly supported scaling functions have only a finite number of filter coefficients with non-zero values, which enables easy handling of finite geometries and imposition of boundary conditions. Mitra and Gopalakrishnan [30] provide a complete prescription of the use of Daubechies compactly supported wavelets for solving 1-D wave equations. The steps involved in the WSFE method are briefly described here; readers interested in further details may refer to [23,30]. Let the time–space displacement variable u(x, t) be discretized at n points in the time window (0, tf) and s = 0, 1,..., n1 be the sampling points, then t ¼ Dt s where Dt is the time interval between two sampling points. The function u(x, t) can be approximated at an arbitrary scale as

uðx; tÞ ¼ uðx; sÞ ¼

X uk ðxÞ uðs  kÞ;

k2Z

ð25Þ

k

where uk ðxÞ are the approximation coefficients at a certain spatial dimension and uðsÞ are scaling functions associated with Daubechies wavelets. The other translational and rotational displacements v(x,y,t) and /ðx; y; tÞ are approximated similarly. By substituting these approximations into the first equation of Eq. (19) we get

A11

X d2 uk 2

uðs  kÞ  B11

X @ 2 /k @x2

uðs  kÞ þ

X sk uðs  kÞ

k dx k k I0 X I1 X 00 00 uk u ðs  kÞ þ 2 /k u ðs  kÞ ¼ 0  2 Dt k Dt k

2

d uj 2

dx

 B11

ð28Þ

2

where ½C1  is the first order connection coefficient matrix. Next, the coupled ODEs in Eq. (28) are decoupled using eigenvalue analysis of C1 given by

C1 ¼ UPU1

ð29Þ

where P is the diagonal matrix with eigenvalues and U is the eigenpffiffiffiffiffiffiffi vectors matrix of C1 . By letting the eigenvalues be icj ði ¼ 1Þ, the final decoupled form of the reduced ODEs given in Eq. (29) can be written as

A11

2 ^j d u 2

dx

 B11

2^ d / j 2

dx

þ ^sj þ I0 c2j uj  I1 c2j /j ¼ 0

ð30Þ

^ j is defined as u ^ j ¼ U1 uj . It should be mentioned here that where u the sampling rate Dt should be less than a certain value to avoid spurious dispersion in simulations using WSFE. In Mitra and Gopalakrishnan [33], a numerical study has been conducted from which the required Dt can be determined depending on the order N of the Daubechies scaling function and frequency content of the input load. The governing PDEs for each adherand are transformed following similar steps as in Eqs. (25)–(30). For further steps, the subscript j (of Eq. (30)) is dropped for simplified notations and all of the following equations are valid for j = 0, 1, ..., n  1. Also the subscripts denoting the adherand number is reintroduced for easy reference. Following these, the transformed governing ODEs for the top adherand are written as 2

2^

2

ð26Þ

jþN2 @ 2 /j 1 X 2 þ sj  2 X ðI0 uk  I1 /k Þ ¼ 0; 2 @x Dt k¼jNþ2 jk

j ¼ 0; 1; 2; . . . ; n  1

2

dx

( 2 ) @ /j 2 2  B11 þ fsj g  I0 ½C1  uj þ I1 ½C1  /j ¼ 0 @x2

^ ^1 ¼ 0 ^ 1c þ I1 c2 / A11 ddxu21c  B11 ddx/21 þ ^s þ I0 c2 u

Here the subscripts denoting the adherand have been dropped again for simplified notation although the terms represent mid plane displacements of the top adherand. Taking inner product on both sides of Eq. (26) with translations of the scaling function (uðs  jÞ for j = 1, 2, . . ., n-1) and using the orthogonal property of Daubechies scaling function results in the cancelation of all the terms except when j = k and yields n simultaneous equations. Thus Eq. (26) reduces to

A11

( 2 ) d uj

^

jA55 ðddxv^21c  ddx/1 Þ  p^ þ I0 c2 v^ 1c ¼ 0 2^

2^

B11 ddxu21c þ D11 ddx/21 þ A55



dv^ 1c dx

 ^1 ¼ 0 ^ 1c þ I2 c2 /  /1 þ h1^s  I1 c2 u ð31Þ

Similarly the governing PDEs for the bottom adherand take the form 2^ 2^ ^2 ¼ 0 ^ 2c þ I1 c2 / A11 ddxu22c  B11 ddx/22  ^s þ I0 c2 u 2  ^ jA55 ddxv^22c  ddx/2 þ p^ þ I0 c2 v^ 2c ¼ 0   2^ 2^ v^ 2c ^2 ¼ 0 ^ 2c þ I2 c2 / B11 ddxu22c þ D11 ddx/22 þ A55 ddx  /2 þ h2^s  I1 c2 u

ð32Þ ð27Þ

where N is the order of the Daubechies wavelet and X2jk are the connection coefficients. For compactly supported wavelets, X2jk are nonzero only in the interval k = j  N + 2 to k = j + N  2. A complete description for the evaluation of connection coefficients of different orders of derivative is given in [31]. It can be observed in the ODEs of Eq. (27) that certain coefficients uj near the vicinity of the boundaries (j = 0 and j = n  1) lie outside the time window [0, tf] defined by j = 0, 1, ..., n  1. Therefore it is necessary to treat the boundaries appropriately. Among several methods available [23], we use wavelet extrapolation technique [32] for solving the boundary value problem for aperiodic boundaries. This approach allows

and the transformed boundary conditions for the two adherands are,

^1 ^ 1c d/ b ¼ A du N  B11 ; xx1 11 dx dx

  dv^ 1c ^ b Q x1 ¼ jA55  /1 ; dx

^1 b ^2 ^ 1c ^ 2c du d/ du d/ c M xx1 ¼ B11 þ D11 N xx2 ¼ A11  B11 ; dx dx dx dx   ^2 ^ 2c dv^ 2c ^ du d/ b Q x2 ¼ jA55  /2 ; c M xx2 ¼ B11 þ D11 dx dx dx

ð33Þ

Once the ODEs are derived for the beam adherands the next step is to compute wavenumbers and wave amplitudes which are used in deriving shape functions of the spectral element.

276

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

The nodal forces given by force boundary conditions in Eq. (33) and unknown coefficients can be related as

^e g ¼ ½T 2 fag fF Fig. 4. Nodal representation of single lap joint.

2.3. Wavenumber computation and spectral element formulation Since the ODEs in Eqs. (31)–(33) are analogous to that in FSFE [34], the WSFE formulation follows similar procedure. Only one spectral element is needed to represent the bonded region of the lap joint (Fig. 2) when there is no discontinuity (such as a damaged region). There are four nodes associated with the bonded double beam element (denoted by dark circles in Fig. 3). Each node has ^ lm ^ lm ; v ^ lm and / 3 degrees of freedom which are labeled as u ðl; m ¼ 1; 2Þ where l and m stand for adherand and node numbers, respectively. The ODEs given by Eqs. (31) and (32) are solved for ^ lm in the frequency domain and their time domain ^ lm ; v ^ lm and / u solutions ulm ðx; tÞ; v lm ðx; tÞ and /lm ðx; tÞ are obtained using inverse wavelet transform. The exact interpolating functions for an element of length L, obtained by solving Eqs. (31) and (32), respectively are

fu1c ; v 1c ; /1 ; u2c ; v 2c ; /2 gT ¼ ½R½Hfag

ð34Þ

where ½H is a diagonal matrix with the diagonal terms given by ½eik1 x ; eik1 ðLxÞ ; eik2 x ; eik2 ðLxÞ ; eik3 x ; eik3 ðLxÞ ; eik4 x ; eik4 ðLxÞ ; eik5 x ; eik5 ðLxÞ ; eik6 x ; eik6 ðLxÞ , and ½R is a 6 x 12 matrix composed of amplitude ratios corresponding to wave numbers k1 , k2 ;    ; k6 . The wavenumbers are obtained by substituting Eq. (34) into Eqs. (31) and (32) which gives the characteristic equation

½WðkÞ½Rfag ¼ f0g

ð35Þ

where W(k), known as wave matrix, is of size 6  6 with non-zero terms given as 2

2

Wð1;1Þ ¼ A11 k þ I0 c2  ks ; Wð1;3Þ ¼ B11 k  I1 c2  h1 ks ; Wð1;4Þ ¼ ks ; Wð1;6Þ ¼ h2 ks 2

Wð2;2Þ ¼ A55 k þ I0 c2  kt ; Wð2;3Þ ¼ iA55 k; Wð2;5Þ ¼ kt 2

Wð3;1Þ ¼ B11 k  I1 c2  h1 ks ; Wð3;2Þ ¼ iA55 k; 2

2

Wð3;3Þ ¼ D11 k  A55 þ I2 c2  h1 ks Wð3;4Þ ¼ h1 ks ; Wð3;6Þ ¼ h1 h2 ks ; 2

Wð4;1Þ ¼ ks ; Wð4;3Þ ¼ h1 ks ; Wð4;4Þ ¼ A11 k þ I0 c2  ks ; 2

Wð4;6Þ ¼ B11 k  I1 c2 þ h2 ks 2

Wð5;2Þ ¼ kt ; Wð5;5Þ ¼ A55 k þ I0 c2  kt ; Wð5;6Þ ¼ iA55 k 2

Wð6;1Þ ¼ h2 ks ; Wð6;3Þ ¼ h1 h2 ks ; Wð6;4Þ ¼ B11 k  I1 c2 þ h2 ks 2

2

Wð6;5Þ ¼ iA55 k; Wð6;6Þ ¼ D11 k  A55 þ I2 c

2  h2 ks

xx22

xx2

A relationship between nodal forces and the displacements can be obtained using Eqs. (36) and (37) as

^e g ¼ ½T 2 ½T 1 1 fu ^ e g ¼ ½Kfu ^ eg fF

ð38Þ

where [K] is the dynamic stiffness matrix of the spectral element. ^e g are The order of the matrix is 12 x 12. When the nodal forces fF ^ e g can be obtained from the known, the nodal displacements fu above relationship. Then the unknown constants fag can be determined using Eq. (36). The transformed displacements at any point along the adherands can then be computed using Eq. (34). With the above step, WSFE formulation for a bonded double beam system is completed. As mentioned earlier, the formulation can be extended to any number of bonded members to study transient dynamics efficiently. In the present work, we study a single lap joint which has attracted attention due to increased use of composites in structural applications. 2.4. Single lap joint: dynamic stiffness matrix derivation A single lap joint consists of two adherands bonded using adhesive as shown in Fig. 1. The adherands extending beyond the bonded regions are modeled as beams where the governing PDEs can be reduced from Eqs. (19) or (21) with no adhesive effect (that is, terms s = p = 0). The force boundary conditions remain the same as in Eq. (20). The elemental dynamic stiffness matrices of individual beams and bonded double beam system are assembled at the appropriate nodes in order to generate global structural stiffness matrix. The assembling procedure is similar to conventional finite element method except that assembly is performed in the transformed frequency domain. Examples of assembling elemental stiffness matrices in the frequency domain to form global stiffness matrix are presented in [21,35]. A nodal representation of the SLJ (in Fig. 1) is shown in Fig. 4. Here only three spectral elements are used to represent the lap joint with two beam elements (nodes 1–2 and 5–6) and a bonded double beam element (nodes 2–3 to 4–5) developed above. The elements are assembled at nodes 2 and 5. Nodal forces and displacements for the SLJ can be written in global coordinate system as

^ eg g fF^eg g ¼ ½Kg fu

ð39Þ

where ½Kg 1818 is the global structural dynamic stiffness matrix defined as

By solving the characteristic equation Eq. (35) the wavenumbers are determined. Wave amplitudes corresponding to each wave number are determined using singular value decomposition and stored in [R]. The term {a} includes 12 unknown constants grouped as {a1, a2, . . ., a12}T which are determined from transformed nodal ^ e } = {u ^ 11 ; v ^ 11 ; quantities. Let’s take nodal displacements {u ^ 11 ; u ^ 21 ; u ^ 12 ; u ^ 22 }T where u ^ 21 ; v ^ 21 ; / ^ 12 ; v ^ 12 ; / ^ 22 ; v ^ 22 ; / ^ 11 ¼ u1c ð0Þ; / v^ 11 ¼ v 1c ð0Þ; /^ 11 ¼ /^ c ð0Þ; u^21 ¼ u2c ð0Þ; v^ 21 ¼ v 2c ð0Þ; /^ 21 ¼ /^ c ð0Þ; u^12 ¼ ^ 12 ¼ / ^ c ðLÞ; u ^ 22 ¼ / ^ c ðLÞ. ^ 12 ¼ v 1c ðLÞ; / ^ 22 ¼ u2c ðLÞ; v ^ 22 ¼ v 2c ðLÞ; / u1c ðLÞ; v From the relationship of the displacements and unknown constants in Eq. (34), the nodal displacements can be written as

^ e g ¼ ½T 1 fag fu

ð37Þ

b b b b c b c b c ^e g ¼ f N where fF xx11 ; Q x11 ; M xx11 ; N xx21 ; Q x21 ; M xx21 ; N xx12 ; Q x12 ; M xx12 ; b b c T b c b ^ ^ N and N xx22 ; Q x22 ; M xx22 g xx11 ¼ N xx1 ð0Þ; Q x11 ¼ Q x1 ð0Þ; M xx11 ¼ b b c b ^ x2 ð0Þ; M xx21 ¼ M ^ xx1 ð0Þ; N xx21 ¼ N ^ xx2 ð0Þ; Q x21 ¼ Q ^ xx2 ð0Þ; N M xx12 ¼ b b c b ^ ^ ^ ^ ^ Nxx1 ðLÞ; Q x12 ¼ Q x1 ðLÞ; M xx12 ¼ Mxx1 ðLÞ; N xx22 ¼ N xx2 ðLÞ; Q x22 ¼ Q x2 ðLÞ; c ^ ðLÞ. M ¼M

ð36Þ

2

K11 b1 6 21 6 Kb1 6 6 6 0 ½Kg  ¼ 6 6 0 6 6 4 0 0

K12 b1

0

0

0

11 K22 b1 þ K

K12

K13

K14

K21

K22

K23

K24

31

32

33

K34

K

K

K

K41

K42

K43

K44 þ K11 b2

0

0

0

K21 b2

0

3

7 0 7 7 7 0 7 7 0 7 7 12 7 Kb2 5

ð40Þ

K22 b2

where each sub matrix has a size of 3  3. Here [Kb1] and [Kb2] are the elemental stiffness matrices of the individual beam elements and [K] is the stiffness matrix of the bonded double beam section of the SLJ. Once the global quantities are determined using

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

Eq. (39) the local displacements can be obtained using respective nodal quantities and elemental shape functions. 3. Results and discussion The WSFE model developed above is used for frequency and time domain wave propagation behavior analysis of composite single lap joints. Frequency domain spectrum and dispersion relations are investigated first. Thereafter, time domain wave propagation results are compared with conventional finite element simulations performed using commercial FE software AbaqusÒ [36]. Example lap joints are made up of AS4/3501-6 graphite-epoxy composite with layup sequence as specified in each case. Each of the beam adherands has a depth of h1 = h2 = 0.01 m and width of b = 0.02 m. Lamina properties are as follows: E1 = 144.48 GPa, E2 = E3 = 9.63 GPa, G23 = G13 = G12 = 4.128 GPa, m23 = 0.3, m13 = m12 = 0.02, and q = 1389 kgm3. Epoxy based structural paste adhesive Hysol EA 9394 with Young’s modulus and Poison’s ratio of 4.24 GPa and 0.45, respectively, is taken as the bondline material. The bondline thickness ta = 0.005 m unless specified otherwise. Shear correction factor j appearing in the formulation is taken as 5/6 in all the examples. The order of Daubechies scaling function used is N = 22, if not mentioned specifically. 3.1. Wave propagation: spectrum and dispersion relations In order to understand the wave propagation behavior of bonded joints it is important to get an insight of frequency dependent wave parameters, spectrum and dispersion relations, which are readily available in the WSFE method. WSFE can be used to obtain the frequency dependent wave characteristics up to a fraction of Nyquist frequency where the fraction is dependent on the order of the Daubechies wavelet used. In order to be able to use WSFE in frequency domain analysis the boundaries arising from finite data lengths need to be treated as periodic. For further details

277

of the use of WSFE for studying frequency dependent wave characteristics, readers may refer to [23,33]. For studying wave propagation in bonded single lap joint, it is necessary to treat the overlapped region and the rest separately. Wave motion of the adherands outside the overlap region is governed by Eq. (15) without the bondline effects. Mahapatra and Gopalakrishnan [34] discuss the wave propagation behavior of such standalone beams (also referred to as single beams elsewhere in the menuscript) in detail. Therefore we focus here on the overlapped region or bonded double beam region of the joint. There are six modes associated with bonded double beam system where the characteristic equation is given by Eq. (35). The wavenumbers are computed by setting the determinant of the wave matrix ðjWðkÞjÞ to zero. The resulting 12th order characteristic equation is very complex and cumbersome to solve analytically. Here the companion matrix method is used to solve the characteristic equation where the roots are determined as eigenvalues of the companion matrix. Once the eigenvalues are obtained they can be used to obtain eigenvectors using singular value decomposition [22]. Some a priori tests can be performed on the characteristic equation in order to get an idea about the presence of propagating and evanescent modes as well as cutoff frequencies. By setting the frequency ðcÞ to zero and solving the characteristic equation for wavenumbers (k), it can be observed that for the uncoupled case (Bij = 0) there are three imaginary wavenumbers corresponding to three evanescent modes. The cutoff frequencies (where the evanescent modes turn into propagating modes) can be evaluated using the characteristic equation by setting k = 0 and solving for frequency. Fig. 5 shows the spectrum relations of bonded double beam with layup sequence of [0]10 for each adherand. This relationship is obtained by solving the characteristic equation with a sampling time of Dt ¼ 0:1 ls: Real and imaginary wavenumbers are plotted along positive and negative Y-axis, respectively. There are three modes with zero wavenumber values at zero frequency, (c ¼ 0), meaning none of these three will propagate at that frequency. These fundamental modes resemble axial, flexural and shear

Fig. 5. Spectrum relations for bonded double beam with layup sequence of [0]10 (zoomed image on the right).

278

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

Fig. 6. Dispersion relations for bonded double beam with layup sequence of [0]10 (a) phase velocity cp ¼ c=k (b) group velocity cg ¼ dc=dk.

12

Group velocity (km/s)

10

Axial

baseline 10% reduction 20% reduction

Shear

8 6

higher order modes

4

Flexural 2 0 0

100

200 Frequency (kHz)

300

400

Fig. 7. Group velocity dispersion comparison of bonded double beam with degraded bondline properties. [0]10 layup sequence is used for each adherand.

modes similar to the ones present in standalone beams [34]. The axial and flexural wavenumbers become real values as the frequency increases while the shear mode does so after a certain cutoff frequency. In a physical sense, the axial and flexural waves propagate at any frequency above zero, while the shear mode awaits propagation until the excitation frequency reaches the corresponding cutoff frequency. The three imaginary modes which start with complex wavenumbers at zero frequency ðc ¼ 0Þ and remain complex until certain cutoff frequencies to become real represent higher order modes present in the double beam system. As mentioned already, these cutoff frequencies can also be determined by substituting k = 0 in the characteristic equation and solving for c. For this example, the cutoff frequency for shear mode is

40.7 kHz and that for the three higher order modes are 86.8, 175.8 and 220.2 kHz. The modes with complex wavenumbers are called evanescent modes which tend to attenuate rapidly. These higher order modes present in the bonded double beam do not propagate as long as the excitation signal has a frequency below cutoff frequencies. Fig. 6 shows phase and group velocities for the bonded double beam, with layup sequence of [0]10 for both top and bottom adherands. Fundamental axial and flexural waves propagate at speeds depending on the excitation frequency and the fundamental shear mode appears in the response after the cutoff frequency mentioned above (40.7 kHz). Mechanical properties of the adhesive layer affect spectrum and dispersion relations. The spring constants, kt and ks , related to adhesive layer moduli and thickness (Eq. (24)) are reduced by a certain percentage and the resulting group velocity diagram is shown in Fig. 7. It is observed that the fundamental modes remain intact irrespective of the slight variation in adhesive properties or the bondline thickness. However, higher order modes show significant decrease in cutoff frequencies, with the exception of the first higher order mode. 3.2. Wave propagation: time domain analysis The WSFE model for single lap joint is validated with conventional FE simulation results. Composite material AS4/3501-6 with mechanical properties as given above is used with a ply layup sequence of [0]10. For this example, the adherands 1 and 2 (Fig. 1) have length of 0.35 m with an overlap length of 0.1 m. Thickness of each adherand is taken as 0.01 m and the bondline thickness is 0.005 m. Fixed boundary conditions are applied at the left end of adherand 1.

Fig. 8. Input impulse load in (a) time domain (b) frequency domain.

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

279

Fig. 9. (a) Transverse (b) axial response of the single lap joint for transverse impulse input load applied and measured at the free end.

Fig. 10. (a) Transverse response for single lap joint measured at far end of the bonded region due to transversely applied impulse input load at the free end.

Impulse load with a frequency content of 44 kHz (Fig. 8) is applied at the free end of adherand 2. The single lap joint is modeled in WSFE using two beam elements and one double beam element developed in this work. AbaqusÒ model for the SLJ is meshed with 5700, 2D plane stress elements with 8 elements along the thickness of each adherand. Adhesive bondline consists of 100 elements

(one element along the thickness) with isotropic properties (Hysol EA 9394) as specified above. FE simulations are performed using AbaqusÒ Explicit dynamic solver. Fig. 9 shows the transverse and axial response comparisons for the impulse input load applied into transverse direction. Here we can see that WSFE and FE (AbaqusÒ) predictions are in very good agreement. Transverse response at the excitation location (Fig. 9 (a)) shows peak 1 as the excitation input and peak 2 is the reflection from the overlap boundary (interface II in Fig. 1). Once the transmitted waves enter the overlapped area, another boundary (interface I in Fig. 1) causes partial reflection shown as peak 3. Due to the complicated nature of the boundaries in the overlapped area partial reflections and transmissions are observed repeatedly. The largest response appearing after peak 3 is the reflection from the fixed boundary located at far left end of the SLJ. Although the individual adherands consist of symmetric balanced layup, wave propagation across the lap joint exhibits axialtransverse coupling. It means that both axial and flexural waves are generated irrespective of the input excitation direction. Fig. 9 (b) shows the axial response measured at the free end of adherand 2 (Fig. 1) for a transversely excited impulse input. Since the adherand laminates are symmetric and balanced, no coupling exists as long as the wave is propagated through standalone beam. Once the flexural wave reaches the overlapped region boundary (interface II in Fig. 1), it undergoes mode conversion resulting in axial wave propagation. This mode converted and reflected wave (peak 1 in Fig. 9 (b)) reaches the free end first, followed by wave reflections from boundary at the far end of the

Fig. 11. (a) Axial (b) transverse response for single lap joint for axial impulse input load applied and measured at the free end.

280

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

Fig. 12. Axial response of the lap joint measured at far end of the bonded region due to axially applied impulse input load at the free end.

[24,25] compared simulations times of AbaqusÒ with standard solver and observed over two orders of magnitude reduction in computation time when compared with WSFE. Fig. 11 shows the axial and transverse response comparisons for an axial impulse input load. For this example the lengths of adherands 1 and 2 are extended to 1.5 m and 1.0 m, respectively, with 0.5 m bondline length in order to clearly visualize high velocity axial wave reflections from each boundary. Again, it is observed that WSFE predictions are in very good agreement with FE results. We have already observed through group velocity dispersion diagrams that the axial waves are non-dispersive. Fig. 11(a) shows that reflected waveforms preserve their shape irrespective of the propagation lengths, which is a major feature of non-dispersive signals in time domain. In the axial response shown, peak 1 is input load and peak 2 is the fixed boundary reflection. The smaller amplitude peaks between those two are created by boundaries in the overlapped region. Similar to the observations made in the previous example (Fig. 9), mode converted transverse waves are reflected as shown in Fig. 11 (b). Due to the longer wave guide and considerably smaller group velocity the mode converted transverse response reaches the free end at a much later time. The axial velocity response of the top adherand (adherand 1 in Fig. 1) measured at the far end of the overlapped region (interface I in Fig. 1) for the same axial input at free end is shown in Fig. 12. Here too the WSFE prediction is in very good agreement with FE (AbaqusÒ) results. 3.3. Coupled wave propagation in bonded double beam system

Fig. 13. Semi-infinite bonded double beam for studying wave propagation in different modes.

overlapped region and the fixed boundary. The transverse velocity response of the top adherand (adherand 1 in Fig. 1) measured at the far end of the overlapped region (interface I in Fig. 1) for the same input is shown in Fig. 10. Here also the WSFE prediction is in very good agreement with FE (AbaqusÒ) results. All the transverse response results (Figs. 9 and 10) are obtained from a single simulation performed in a computer with IntelÒ Core™ i7-3770 (3.40 GHz) processor. Computational times taken to complete analysis for WSFE and AbaqusÒ (explicit solver) are recorded as 3.81 s and 74 s, respectively. It can be noticed that there is a significant reduction in computation time for WSFE when compared with conventional FE (explicit) analysis. Authors’ previous studies

Examination of Eq. (35) indicates that wave propagation across the bonded double beam system is coupled in nature, which leads to propagation of multimodal waves for a given input load. Fig. 13 shows the semi-infinite double beam considered in this study in order to understand the presence of these modes. Here the input tone burst with 250 kHz central frequency is applied at the free end of the beam (nodes n1, n2 in Fig. 13). This frequency is selected since all the modes present in the double beam system are propagating modes at this frequency. The beam is modeled using semiinfinite or throw-off elements available in spectral finite element analysis. Hence, there are no boundary reflections present in the response which makes the signal analysis and interpretation less complicated. In reality, this beam configuration corresponds to a very long cantilever beam whose response at the said locations does not contain any fixed boundary reflections inside the time window of interest. Studying wave propagation using semi-infinite or throw-off elements is fairly straight forward using the spectral finite element method as reported previously [21,22]. Similar

( Δd = 2m )

( Δd = 2m )

( Δd = 4m )

( Δd = 4m )

Fig. 14. Axial and transverse velocity responses measured at 2 m and 4 m away from the free end for axial input load (250 kHz toneburst).

281

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

( Δd = 2m )

( Δd = 2m )

( Δd = 4m )

( Δd = 4m )

Fig. 15. Axial and transverse velocity responses measured at 2 m and 4 m away from the free end for transverse input load (250 kHz toneburst).

approach is used here in order to model the semi-infinite beam without repeating the formulation. Use of this throw-off element at the structural boundary amounts to absorption of all the energy from the structure and hence maximum damping artificially to the response. In the examples, the order of Daubechies scaling function used for time approximation is set to be N = 36 and all other parameters are same as above examples. The modes that may appear in the time domain responses is dependent on the nature of loading conditions as well. Therefore the toneburst input load is first applied at both nodes (n1, n2 in Fig. 13) into axial and transverse (X and Y) directions. The responses are recorded at a distance ðDd ¼Þ 2 and 4 m away from the excitation location (free end) along the beam length. Fig. 14 shows the axial and transverse responses recorded at the said locations for an axial input. Here the axially excited toneburst travels non-dispersively and no coupling exists (no transverse response is observed). However, when the beam is excited transversely the axial-flexural coupling comes into effect as shown in Fig. 15. Here the amplitudes of wave packets appeared in axial response is much smaller than that of flexural wave. Therefore the Y-axis scale is zoomed in by one order of magnitude compared to the flexural response scale. Presence of adhesive layer leads to this coupling

behavior which causes wave propagation in multiple modes. In the axial response the first wave packet is an axial mode and the second wave packet is a shear mode. According to the group velocity dispersion curves for the bonded double beam system (Fig. 6) the axial and shear modes have higher group velocities (which are close to each other in magnitude) compared to the flexural mode. The axial and flexural modes are non-despersive while the shear mode displays dispersive behavior at this excitation frequency. Therefore the shear mode is dispersed continuously as it propagates along the beam while the other two modes retain their wave shape. As mentioned previously, the coupling behavior is also dependent on the input loading conditions. To investigate this effect the input tone burst with 250 kHz central frequency is applied to one of the beams only and the resulting wave propagation on the top beam is studied. Excitation given to one adherand resulting in wave propagation in the other adherand simulates single lap joints, which is very important from practical applications point of view. Figs. 16 and 17 show the axial and transverse responses for the input load applied to node n1. The presence of axial-flexural modes is observed even when the excitation is given axially into only one of the two beams. This

( Δd = 2m )

( Δd = 2m )

( Δd = 4m )

( Δd = 4m )

Fig. 16. Axial and transverse velocity responses measured at 2 m and 4 m away from the free end for axial input load (250 kHz toneburst) at n1.

282

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283

( Δd = 2m )

( Δd = 2m )

( Δd = 4m )

( Δd = 4m )

Fig. 17. Axial and transverse velocity responses measured at 2 m and 4 m away from the free end for transverse input load (250 kHz toneburst) at n1.

observation is different from Fig. 14 where both adherands were excited. In the axial response (Fig. 16), the non-dispersive axial mode is followed by a dispersive shear mode. Due to the presence of coupling behavior the two modes appear in the transverse response too. In addition, the axial excitation causes flexural wave propagation which appears as the third wave packet in the transverse response. Here also the amplitudes of transverse waves are much smaller than the axial waves. Hence the Y-axis is enlarged for better visualization. When the top adherand is excited transversely using the same toneburst input, the response is obtained as shown in Fig. 17. Similar to previous observations, coupling effect causes the input wave packet to propagate in multiple modes. Axial response is similar to previous observation (Fig. 15) wherein both adherands were excited. However, there is an additional dispersive wave packet in the transverse response which appears following the non-dispersive flexural mode. This mode is one of the higher order flexural modes appearing in the dispersion curves (Fig. 6) which shows dispersive behavior at the given excitation frequency. It is necessary to mention that the nature of wave propagation along adhesively bonded beams is highly dependent on the adhesive properties which give rise to the coupling terms in the wave matrix (Eq. (35)). Therefore depending on the adhesive properties we may or may not observe certain modes in the time domain responses. Hence it is important to perform frequency domain analysis to understand and resolve complicated wave behavior present in bonded beam systems. 4. Conclusions A new wavelet spectral finite element (WSFE) model was developed for studying transient dynamics and wave propagation across adhesively bonded double beams. The beam formulation was based on first order shear deformation theory which yields accurate results at relatively high frequencies, unlike previous models based on the classical laminate theory. Equations of motion for the bonded double beams were derived using Hamilton’s principle. The adhesive layer was modeled as a line of continuously distributed tension/compression and shear springs. Daubechies compactly supported scaling functions were used to transform the governing partial differential equations into ordinary differential equations. The ODEs were solved using spectral finite element technique where the dynamic stiffness matrix relates the nodal forces and displacements in the transformed domain. Frequency

domain spectrum and dispersion relation results were presented for gaining insight into the wave propagation phenomena. It was observed in the frequency domain simulations that the fundamental axial, flexural and shear modes resembled that of standalone beams. Higher order modes were imaginary at low frequencies and hence did not participate in propagation at low frequencies. However, after a certain cutoff frequency these higher order modes started propagating, thereby complicating the wave propagation behavior. Time domain wave propagation results were in very good agreement with conventional finite element simulations using AbaqusÒ. Time domain wave propagation analysis of single lap joint examples showed that mode conversion occurred once the waves reached the overlapped area. Therefore it was concluded that coupled axial - flexural wave propagation was present in the joints irrespective of the layup sequence of the adherands. A numerical example was also presented to resolve complex coupled wave propagation in bonded beams. Results showed that the adhesive layer caused coupling terms resulting in multi-modal wave propagation. References [1] Wahab AMM. Fatigue in adhesively bonded joints: a review. ISRN Mater Sci 2012;2012:1–25. [2] Katnam KB, Da Silva LFM, Young TM. Bonded repair of composite aircraft structures: A review of scientific challenges and opportunities. Prog Aerosp Sci 2013;61:26–42. [3] Pantelakis S, Tserpes KI. Adhesive bonding of composite aircraft structures: Challenges and recent developments. Sci China Phys Mech Astron 2014;57(1):2–11. [4] Baker AA. In: Rose LRF, Jones R, editors. Advances in the Bonded Composite Repair of Metallic Aircraft Structure, vol. 1. Oxford, UK: Elsevier; 2003. [5] G. Gardiner, Certification of bonded composite primary structures. HighPerformance Composites 2014: pp 50 – 57, http://www.compositesworld. com/articles/certification-of-bonded-composite-primary-structures. [6] Adams RD, Cawley P. A review of defect types and nondestructive testing techniques for composites and bonded joints. NDT Int 1988;21(4):208–22. [7] Gozdecki C, Smardzewski J. Detection of failures of adhesively bonded joints using the acoustic emission method. Holzforschung 2005;59(2):219–29. [8] Di Scalea FL, Rizzo P, Marzani A. Propagation of ultrasonic guided waves in lapshear adhesive joints: case of incident A0 lamb wave. J Acoust Soc Am 2003;115(1):146–56. [9] Le Crom B, Castaings M. Shear horizontal guided wave modes to infer the shear stiffness of adhesive bond layers. J Acoust Soc Am 2010;127(4):2220–30. [10] Matt H, Bartoli I, di Scalea FL. Ultrasonic guided wave monitoring of composite wing skin-to-spar bonded joints in aerospace structures. J Acoust Soc Am 2005;118(4):2240–52. [11] Lowe MJS, Challis RE, Chan CW. The transmission of Lamb waves across adhesively bonded lap joints. J Acoust Soc Am 2000;107(3):1333–45.

D. Samaratunga et al. / Composite Structures 122 (2015) 271–283 [12] Hanneman SE, Kinra VK. A new technique for ultrasonic nondestructive evaluation of adhesive joints: Part I. Theory. Exp Mech 1992;32(4): 323–31. [13] Hanneman SE, Kinra VK, Zhu C. A new technique for ultrasonic nondestructive evaluation of adhesive joints: part II experiment. Exp Mech 1992;32(4): 332–9. [14] Ni X, Rizzo P. Highly nonlinear solitary waves for the inspection of adhesive joints. Exp mech 2012;52(9):1493–501. [15] Saito H, Tani H. Vibrations of bonded beams with a single lap adhesive joint. J Sound Vib 1984;92(2):299–309. [16] Higuchi I, Sawa T, Suga H. Three-dimensional finite element analysis of singlelap adhesive joints under impact loads. J Adhes Sci Technol 2002;16(12): 1585–601. [17] Sawa T, Higuchi I, Suga H. Three-dimensional finite element stress analysis of single-lap adhesive joints of dissimilar adherands subjected to impact tensile loads. J Adhes Sci Technol 2003;17(16):2157–74. [18] Asgharifar M, Kong F, Carlson B, Kovacevic R. Dynamic analysis of adhesively bonded joint under solid projectile impact. Int J Adhes Adhes 2014;50: 17–31. [19] Alleyne D, Cawley P. A two-dimensional Fourier transform method for the measurement of propagating multimode signals. J Acoust Soc Am 1991;89(3): 1159–68. [20] Moser F, Jacobs LJ, Qu J. Modeling elastic wave propagation in waveguides with the finite element method. Ndt E Int 1999;32(4):225–34. [21] Doyle JF. Wave Propagation in Structures: Spectral Analysis Using Fast Discrete Fourier Transforms. 2nd ed. New York: Springer-Verlag; 1997. [22] Gopalakrishnan S, Roy Mahapatra D, Chakraborty A. Spectral Finite Element Method. London (United Kingdom): Springer-Verlag; 2007. [23] Gopalakrishnan S, Mitra M. Wavelet Methods For Dynamical Problems. Boca Raton (FL): CRC Press; 2010. [24] Samaratunga D, Jha R, Gopalakrishnan S. Wavelet spectral finite element for wave propagation in shear deformable laminated composite plates. Compos Struct 2014;108:341–53.

283

[25] Samaratunga D, Jha R, Gopalakrishnan S. Wave propagation analysis in laminated composite plates with transverse cracks using wavelet spectral finite element method. Finite Elem Anal Des 2014;89:19–32. [26] Reddy JN. Mechanics of Laminated Composite Plates and Shells. Boca Raton: CRC Press; 2004. [27] Reddy JN. Energy Principles and Variational Methods in Applied Mechanics. John Wiley & Sons; 2002. [28] F. Mortensen, Development of Tools for Engineering Analysis and Design of High-Performance FRP-Composite Structural Elements. Aalborg: Aalborg Universitetsforlag, 1998. (Special Report: Institute of Mechanical Engineering; No. 37). [29] Daubechies I. Ten Lectures on Wavelets, vol. 61. Philadelphia: Society for Industrial and Applied Mathematics; 1992. [30] Mitra M, Gopalakrishnan S. Spectrally formulated wavelet finite element for wave propagation and impact force identification in connected 1-D waveguides. Int J Solids Struct 2005;42:4695–721. [31] Beylkin G. On the representation of operators in bases of compactly supported wavelets. SIAM J Numer Anal 1992;6:1716–40. [32] Williams JR, Amaratunga K. A discrete wavelet transform without edge effects using wavelet extrapolation. J Fourier Anal Appl 1997;3:435–49. [33] Mitra M, Gopalakrishnan S. Extraction of wave characteristics from wavelet based spectral finite element formulation. Mech Syst Signal Process 2006;20: 2046–79. [34] Roy Mahapatra D, Gopalakrishnan S [A spectral finite element model for analysis of axial-flexural-shear coupled wave propagation in laminated composite beams]. Compos. Struct 2003;59(1):67–88. [35] Samaratunga D, Jha R, Gopalakrishnan S [SHM of composite skin-stiffener structures using wavelet spectral finite element method]. In: Bakis Charles E, editor. Proceedings of the American Society for Composites: Twenty-eighth Technical Conference on Composite Materials. College Park (PA): Destech Pubns Inc; 2013. p. 9–11. [36] Abaqus Users’ Manual, Version 6.12-1, Dassault Systémes Simulia Corp., Providence, RI.