# Effects of Delamination on Guided Waves in a Symmetric Laminated Composite Beam

Authors: Seungwan Kim; Usik Lee
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

For successful structural health monitoring and structural integrity evaluation of a laminated composite structure, it is important to study the effects of delamination on the propagations of the guided waves in a delaminated composite beam by using an accurate and computationally efficient method. Thus, we developed a “frequency-domain” spectral element model for the symmetric composite beams. First-order-shear-deformation-theory (FSDT) based Timoshenko beam theory and Mindlin-Herrmann rod theory are adopted for the flexural (bending) waves and axial (extensional) waves, respectively. A spectral element model is derived from the governing equations of motion by using the variation method in the frequency domain. After validating the accuracy of the proposed spectral element model, the model is used to investigate the effects of delamination on the propagation of guided waves in examples of composite beams.


---

## Body

## 1. Introduction

Due to the many inherent advantages of composite materials over metallic materials such as high performance, stiffness, strength, and low weight, composite materials have been increasingly applied in the fields of aerospace, mechanical engineering, civil engineering, and other technologies over the past decades [1]. However, composite materials are still susceptible to delamination and other modes of damage. Therefore, there is a growing need for structural health monitoring (SHM) of composite structures.Recently, guided wave-based SHM techniques have received considerable attention due to the guided wave characteristics of long-range propagation and sensitivity to various modes of damage [2]. The smaller wavelength of guided waves is more effective for detecting or monitoring incipient small damages; thus, high frequency guided waves are essential for most SHM techniques. Piezoelectric transducers have been used to excite or sense high frequency guided waves (on the order of MHz).Knowledge of the effects of damage on guided wave propagation is very useful for successful structural health monitoring, and wave propagation analysis has been extensively reported in the literature. In practice, ultrasonic guided waves in damaged composite structures are very complicated due to anisotropic material properties, multilayered construction, and multiple reflections and mode conversions at material and geometric discontinuities. Thus it is essential to develop an efficient, low-cost, and accurate wave propagation analysis method to investigate ultrasonic guided waves in damaged composite structures.As described in the literature, a number of computational methods have been used for the analysis of guided wave propagation. These include the finite element method (FEM) [3–6], strip element method [7], boundary element method [8], local interaction simulation approach (LISA) [9], and “time-domain” spectral element method [10–12]. Among these approaches, the FEM, which is normally used in the time-domain, is certainly one of the most popular methods for wave propagation analysis due to its particular advantages. To successfully capture high frequency waves, the finite element mesh size must be smaller than the wavelengths involved, which can be very small at high frequencies. This can result in an extremely large finite element model. Hence, the FEM becomes increasingly inaccurate and inefficient as the frequency increases. This problem can be resolved by adopting the “frequency-domain” spectral element method (SEM), wherein frequency-dependent interpolation functions are used to generate an exact dynamic stiffness matrix (called the “spectral element matrix”) by which high frequency waves can be successfully captured [13, 14]. The SEM provides solutions that are often referred to as “exact” due to the use of an exact dynamic stiffness matrix. It also has the advantage of reducing computational effort due to a significant reduction in the number of degrees of freedom (DOFs) by representing a uniform structural member as a single spectral element model, regardless of its dimension. It is worth mentioning that, though the “time-domain” SEM initiated by Patera [11] is also often called by the same name, “SEM,” in the literature, it is completely different from the “frequency-domain” SEM because it belongs to a p-version FEM where the higher-order polynomials such as orthogonal Chebyshev polynomials are used as the shape functions.In the literature, some researchers have used the “frequency-domain” SEM to study the propagations of guided waves in a cracked rod [15], cracked beam [16], delaminated beam [17], debonded damaged reinforced concrete beam [18], cracked composite beam [19], and delaminated composite beam [20, 21]. For composite beams, Palacz et al. [21] used a spectral element model based on the first-order shear deformation theory (FSDT), wherein the axial displacement and the lateral contraction in the thickness direction (i.e., the thickness contraction) were not considered. Nag et al. [20] considered axial-flexural-shear coupling in their FSDT-based spectral element model but ignored the effects of thickness contraction. The thickness contraction is known to be important for axial (or extensional) waves at high frequencies via the axial-thickness contraction coupling, which can be represented by the Mindlin-Herrmann rod theory [22]. Thus, it is very important to fully take into account the effects of axial-thickness contraction coupling in the spectral element model in order to successfully investigate the effects of delamination on axial and flexural (bending) waves in a delaminated composite beam at high frequencies.Thus, the purposes of this study are (1) to derive governing equations of motion for a symmetric composite beam by adopting both Timoshenko beam theory and Mindlin-Herrmann rod theory, (2) to develop a “frequency-domain” spectral element model for a finite composite beam element, and (3) to investigate the effects of delamination on guided waves in some examples of composite beams.


## 2. Derivation of the Equations of Motion

Consider a symmetric composite beam as shown in Figure1, where spatial coordinates x and z represent the axis and thickness directions of the composite beam, respectively. The composite beam has plane motion in the x-z plane. Axial and flexural displacements are represented by u(x,z,t) and w(x,z,t), respectively. Based on FSDT-based Timoshenko beam theory [1], these displacement fields can be represented by
(1)w(x,z,t)=w0(x,t)+zψ(x,t),u(x,z,t)=u0(x,t)-zθ(x,t),
where θ(x,t) is the rotation of the cross-section of the composite beam and ψ(x,t) is the thickness contraction. The subscript “0” denotes a quantity at the midplane of the composite beam. The equations of motion for a composite beam element of length l can be derived from Hamilton’s principle [1]:
(2)∫0t(δT-δV+δW)dt=0,
where T, V, and δW are the kinetic energy, the potential energy, and the virtual work done by external forces, respectively. By using the displacement fields given by (1) and the constitutive relations for the composite laminates provided in Appendix A, the kinetic and strain energies can be given by
(3)T=12∫0l(Iθ˙2ρAw˙02+2ρRw˙0ψ˙+ρIψ˙2+ρAu˙02-2ρRu˙0θ˙+ρIθ˙2)dx,V=12∫0l(bA11au0′2+κ1bA55bw0′2-2κ1bA55bw0′θ+κ1bA55bθ2+bD11bθ′2+κ2bD55aψ′2+bA33aψ2+2bA13au0′ψ)dx,
where
(2)(Aija,Bija,Dija)=∑k=1N∫zk-1zkC-ij*(k)(1,z,z2)dz,(Aijb,Bijb,Dijb)=∑k=1N∫zk-1zkQij*(k)(1,z,z2)dz,(ρA,ρR,ρI)=b∑k=1Nρ(k)(zk-zk-1,12(zk2-zk-12),13(zk3-zk-13)),
where b is the width, ρb(k) is the mass density of the kth laminate of the composite beam, and material properties C-ij* and Qij* are defined in Appendix A. Finally the virtual work δW can be represented as
(5)δW=N1(t)δu0(0)+N2(t)δu0(l)+M1(t)δθ(0)+M2(t)δθ(l)+Q1(t)δw0(0)+Q2(t)δw0(l)+∫0l[fu(x,t)δu0+fw(x,t)δw0+fθ(x,t)δθ]dx,
where Nj(t), Mj(t), and Qj(t) (j=1,2) are, respectively, the resultant axial forces, bending moments, and the transverse shear forces applied at the boundaries x=0 and l. fu(x,t), fw(x,t), and fθ(x,t) are externally applied distributed forces.Laminated composite beam: (a) geometry and (b) deformations.
(a)
(b)By substituting (3) and (5) into (2) and applying integration by parts, the governing equations of motion and associated natural boundary conditions can be derived simultaneously. The equations of motion are obtained in the following form:
(6)D2w′′+D1w′+D0w+μw¨+p(x,t)=0,
where
(7)w(x,t)={u0(x,t),w0(x,t),θ(x,t),ψ(x,t)}T,p(x,t)={fu(x,t),fw(x,t),fθ(x,t),0}T,
and D0,D1,D2, and μ are four-by-four matrices defined in Appendix B. The natural boundary conditions are also obtained in the form
(8)D2w′(0,t)+D1Uw(0,t)=-f1(t),D2w′(l,t)+D1Uw(l,t)=f2(t),
where D1U is the four-by-four matrix defined in Appendix B and
(9)f1(t)={N1(t),Q1(t),M1(t),0}T,f2(t)={N2(t),Q2(t),M2(t),0}T.


## 3. Spectral Element Formulation and Analysis

### 3.1.  Spectral Element Model for a Composite Beam Element

To formulate a spectral element model for a composite beam element of lengthl, the equations of motion and the associated boundary conditions must first be transformed in the frequency domain. To that end, the time-domain displacement fields and forcing terms are all represented in spectral forms by applying the discrete Fourier transform theory [23] as follows:
(10){w(x,t),p(x,t),f1(x,t),f2(x,t)}=1N∑n=0N-1{w-n(x;ωn),p-n(x;ωn),f-1n(x;ωn),f-2n(x;ωn)}eiωnt,
where i=-1 is the imaginary unit and the quantities with overbars are the spectral (Fourier) components of the corresponding time-domain quantities. N/2 is the number of spectral components up to the Nyquist frequency to be considered in the FFT-based spectral analysis. By using (10), the time-domain equations of motion given by (6) can be transformed in the frequency domain as follows:
(11)D2w-′′+D1w-′+D0w--ω2μw-+p-(x;ω)=0.
Similarly the boundary conditions can be also transformed into the frequency domain as
(12)D2w-′(0;ω)+D1Uw-(0;ω)=-f-1(ω),D2w-′(l;ω)+D1Uw-(l;ω)=f-2(ω).
In (11) and (12) and also in the following derivations, the subscript n is omitted for the sake of brevity.Consider the homogeneous frequency-domain equations of motion that can be reduced from (11) by setting p-(x;ω)=0 as follows:
(13)D2w-′′+D1w-′+D0w--ω2μw-=0.
The general solutions of (13) are assumed in the form
(14)w-(x;ω)=rae-ikx,
where a is a constant and k is the wavenumber. The vector r is defined by
(15)r={1,rw,rθ,rψ}T.
Substituting (14) into (13) yields the following eigenvalue problem:
(16)(-k2D2-ikD1+D0-ω2μ)ra=0.
The wavenumber-frequency dispersion equation can be derived from (16) in the form
(17)k8+c1k6+c2k4+c3k2+c4=0.
By using eight wavenumbers kj (j=1,2,3,…,8) computed from (17), the general solutions of (13) can be written as
(18)w-(x;ω)=∑j=18rjaje-ikjx=RE(x;ω)a,
where
(19)a={a1,a2,a3,…,a8}T,(20)R=[r1,r2,r3,…,r8],(21)E(x;ω)=diagonal[e-ikjx](j=1,2,3,…,8).
The jth eigenvectors rj={1,rwj,rθj,rψj}T in (20) can be readily computed from (16) by specifying that k=kj.For a finite composite beam element of lengthl as shown in Figure 2, the nodal degrees of freedom in the frequency domain (called spectral nodal DOFs) are defined as
(22)d-(ω)={d-1(ω)d-2(ω)}={w-(0;ω)w-(l;ω)},
where
(23)d-1(ω)={u-01,w-01,θ-1,ψ-1}T,(24)d-2(ω)={u-02,w-02,θ-2,ψ-2}T.Finite composite beam element: (a) spectral nodal DOFs and (b) spectral nodal forces and moments.
(a)
(b)By applying (18) to (22), the spectral nodal DOFs vector d-(ω) can be represented as
(25)d-(ω)=[R(ω)E(0;ω)R(ω)E(l;ω)]a=H(ω)a.
Using (25), the constant vector a can be removed from (18) to write the general solutions in terms of d-(ω) as follows:
(26)w-(x;ω)=N(x;ω)d-(ω),
where
(27)N(x;ω)=R(ω)E(x;ω)H(ω)-1.To formulate the spectral element model by using the variational method, the weak form of (6) is obtained from its weighted-integral statement, by integrating by parts and then applying the boundary conditions, as follows:
(28)∫0l(δw-′TD2w-′+δw-′TD1Uw-+δw-TD1UTw-′-δw-TD0w-+δw-Tω2μw-)dx=∫0lδw-Tp-dx+δd-Tf-b,
where
(29)p-(x;ω)={f-u(x;ω),f-w(x;ω),f-θ(x;ω),0}T,f-b(ω)={N-1(ω),Q-1(ω),M-1(ω),0,N-2(ω),Q-2(ω),M-2(ω),0}T.
By substituting (26) into (28), the spectral element equation can be obtained in the following form:
(30)S(ω)d-(ω)=f-(ω),
where
(31)f-(ω)=f-b(ω)+∫0lN(x;ω)Tp-(x;ω)dx
and S(ω) is the spectral element matrix defined in Appendix C.


### 3.2. Spectral Element Model for a Delaminated Composite Beam Element

Figure3 shows a composite beam element that contains a delamination inside. Following the approaches introduced in [20, 21], a spectral element model for the delaminated composite beam element can be obtained by assembling spectral element equations for the four segments shown in Figure 3. The spectral element equations for four segments can be written in partitioned form as
(32)[S11(q)S12(q)S21(q)S22(q)]{d-1(q)d-2(q)}={f-1(q)f-2(q)}(q=1,2,U,L),
where
(33)d-1(q)={u-01(q),w-01(q),θ-1(q),ψ-1(q)}T,d-2(q)={u-02(q),w-02(q),θ-2(q),ψ-2(q)}T
and the superscripts q=1,2,U, and L represent the left, right, upper, and lower segments in Figure 3, respectively. The kinematic connectivity conditions at the vertical interfaces between four segments can be represented as
(34){d-1(U)d-2(U)}=[TU00TU]{d-2(1)d-1(2)},{d-1(L)d-2(L)}=[TL00TL]{d-2(1)d-1(2)},
where
(35)TU=[10-hU0010hU00100001],TL=[10hL0010-hL00100001].
By using (34), the spectral element equations for four segments represented by (32) can be assembled in the following form:
(36)[S11(1)S12(1)00S21(1)S22(1)+S11(U+L)S12(U+L)00S21(U+L)S22(U+L)+S11(2)S12(2)00S21(2)S22(2)]{d-1(1)d-2(1)d-1(2)d-2(2)}={f-1(1)f-2(1)+f-1(U+L)f-1(2)+f-2(U+L)f-2(2)},
where
(37)Sij(U+L)=TUTSij(U)TU+TLTSij(L)TLf-i(U+L)=TUTf-i(U)+TLTf-i(L)(i,j=1,2).Figure 3
Finite composite beam element with a delamination.It is very often to use piezoelectric (PZT) transducers to excite or measure the high frequency dynamic responses of a structure. To that end, the PZT transducers are usually bonded on a surface of the structure. The spectral element model for the composite beam element whose upper surface is bonded with a PZT transducer is available from [24].


### 3.3. Assembly and Spectral Element Analysis

The assembly of the spectral elements for a composite beam structure is quite straightforward because it can be readily accomplished with the same approach that is commonly used in a standard FEM. The spectral element analysis methods for both frequency-domain solutions (i.e., frequency response functions and natural frequencies) and time-domain solutions (i.e., time histories of dynamic responses and waves) have been well developed in the literature (see [13, 14]).


## 3.1.  Spectral Element Model for a Composite Beam Element

To formulate a spectral element model for a composite beam element of lengthl, the equations of motion and the associated boundary conditions must first be transformed in the frequency domain. To that end, the time-domain displacement fields and forcing terms are all represented in spectral forms by applying the discrete Fourier transform theory [23] as follows:
(10){w(x,t),p(x,t),f1(x,t),f2(x,t)}=1N∑n=0N-1{w-n(x;ωn),p-n(x;ωn),f-1n(x;ωn),f-2n(x;ωn)}eiωnt,
where i=-1 is the imaginary unit and the quantities with overbars are the spectral (Fourier) components of the corresponding time-domain quantities. N/2 is the number of spectral components up to the Nyquist frequency to be considered in the FFT-based spectral analysis. By using (10), the time-domain equations of motion given by (6) can be transformed in the frequency domain as follows:
(11)D2w-′′+D1w-′+D0w--ω2μw-+p-(x;ω)=0.
Similarly the boundary conditions can be also transformed into the frequency domain as
(12)D2w-′(0;ω)+D1Uw-(0;ω)=-f-1(ω),D2w-′(l;ω)+D1Uw-(l;ω)=f-2(ω).
In (11) and (12) and also in the following derivations, the subscript n is omitted for the sake of brevity.Consider the homogeneous frequency-domain equations of motion that can be reduced from (11) by setting p-(x;ω)=0 as follows:
(13)D2w-′′+D1w-′+D0w--ω2μw-=0.
The general solutions of (13) are assumed in the form
(14)w-(x;ω)=rae-ikx,
where a is a constant and k is the wavenumber. The vector r is defined by
(15)r={1,rw,rθ,rψ}T.
Substituting (14) into (13) yields the following eigenvalue problem:
(16)(-k2D2-ikD1+D0-ω2μ)ra=0.
The wavenumber-frequency dispersion equation can be derived from (16) in the form
(17)k8+c1k6+c2k4+c3k2+c4=0.
By using eight wavenumbers kj (j=1,2,3,…,8) computed from (17), the general solutions of (13) can be written as
(18)w-(x;ω)=∑j=18rjaje-ikjx=RE(x;ω)a,
where
(19)a={a1,a2,a3,…,a8}T,(20)R=[r1,r2,r3,…,r8],(21)E(x;ω)=diagonal[e-ikjx](j=1,2,3,…,8).
The jth eigenvectors rj={1,rwj,rθj,rψj}T in (20) can be readily computed from (16) by specifying that k=kj.For a finite composite beam element of lengthl as shown in Figure 2, the nodal degrees of freedom in the frequency domain (called spectral nodal DOFs) are defined as
(22)d-(ω)={d-1(ω)d-2(ω)}={w-(0;ω)w-(l;ω)},
where
(23)d-1(ω)={u-01,w-01,θ-1,ψ-1}T,(24)d-2(ω)={u-02,w-02,θ-2,ψ-2}T.Finite composite beam element: (a) spectral nodal DOFs and (b) spectral nodal forces and moments.
(a)
(b)By applying (18) to (22), the spectral nodal DOFs vector d-(ω) can be represented as
(25)d-(ω)=[R(ω)E(0;ω)R(ω)E(l;ω)]a=H(ω)a.
Using (25), the constant vector a can be removed from (18) to write the general solutions in terms of d-(ω) as follows:
(26)w-(x;ω)=N(x;ω)d-(ω),
where
(27)N(x;ω)=R(ω)E(x;ω)H(ω)-1.To formulate the spectral element model by using the variational method, the weak form of (6) is obtained from its weighted-integral statement, by integrating by parts and then applying the boundary conditions, as follows:
(28)∫0l(δw-′TD2w-′+δw-′TD1Uw-+δw-TD1UTw-′-δw-TD0w-+δw-Tω2μw-)dx=∫0lδw-Tp-dx+δd-Tf-b,
where
(29)p-(x;ω)={f-u(x;ω),f-w(x;ω),f-θ(x;ω),0}T,f-b(ω)={N-1(ω),Q-1(ω),M-1(ω),0,N-2(ω),Q-2(ω),M-2(ω),0}T.
By substituting (26) into (28), the spectral element equation can be obtained in the following form:
(30)S(ω)d-(ω)=f-(ω),
where
(31)f-(ω)=f-b(ω)+∫0lN(x;ω)Tp-(x;ω)dx
and S(ω) is the spectral element matrix defined in Appendix C.


## 3.2. Spectral Element Model for a Delaminated Composite Beam Element

Figure3 shows a composite beam element that contains a delamination inside. Following the approaches introduced in [20, 21], a spectral element model for the delaminated composite beam element can be obtained by assembling spectral element equations for the four segments shown in Figure 3. The spectral element equations for four segments can be written in partitioned form as
(32)[S11(q)S12(q)S21(q)S22(q)]{d-1(q)d-2(q)}={f-1(q)f-2(q)}(q=1,2,U,L),
where
(33)d-1(q)={u-01(q),w-01(q),θ-1(q),ψ-1(q)}T,d-2(q)={u-02(q),w-02(q),θ-2(q),ψ-2(q)}T
and the superscripts q=1,2,U, and L represent the left, right, upper, and lower segments in Figure 3, respectively. The kinematic connectivity conditions at the vertical interfaces between four segments can be represented as
(34){d-1(U)d-2(U)}=[TU00TU]{d-2(1)d-1(2)},{d-1(L)d-2(L)}=[TL00TL]{d-2(1)d-1(2)},
where
(35)TU=[10-hU0010hU00100001],TL=[10hL0010-hL00100001].
By using (34), the spectral element equations for four segments represented by (32) can be assembled in the following form:
(36)[S11(1)S12(1)00S21(1)S22(1)+S11(U+L)S12(U+L)00S21(U+L)S22(U+L)+S11(2)S12(2)00S21(2)S22(2)]{d-1(1)d-2(1)d-1(2)d-2(2)}={f-1(1)f-2(1)+f-1(U+L)f-1(2)+f-2(U+L)f-2(2)},
where
(37)Sij(U+L)=TUTSij(U)TU+TLTSij(L)TLf-i(U+L)=TUTf-i(U)+TLTf-i(L)(i,j=1,2).Figure 3
Finite composite beam element with a delamination.It is very often to use piezoelectric (PZT) transducers to excite or measure the high frequency dynamic responses of a structure. To that end, the PZT transducers are usually bonded on a surface of the structure. The spectral element model for the composite beam element whose upper surface is bonded with a PZT transducer is available from [24].


## 3.3. Assembly and Spectral Element Analysis

The assembly of the spectral elements for a composite beam structure is quite straightforward because it can be readily accomplished with the same approach that is commonly used in a standard FEM. The spectral element analysis methods for both frequency-domain solutions (i.e., frequency response functions and natural frequencies) and time-domain solutions (i.e., time histories of dynamic responses and waves) have been well developed in the literature (see [13, 14]).


## 4. Numerical Results and Discussion

The composite materials (AS/3501-6) for the examples of composite beams considered in this study have the following material properties:E1=144.48 GPa, E2=9.632 GPa, E3=9.412 GPa, G12=4.128 GPa, G13=7.457 GPa, G23=6.516 GPa, ν12=0.3, ν13=0.3, ν23=0.49, and ρc=1389.2 kg/m3 [25].


### 4.1. Validation of Spectral Element Model

To evaluate the accuracy of the present spectral element model, a cantilevered cross-ply symmetric composite beam[0/90]s without delamination was considered as the first example problem. The examples of composite beam has thickness of 5 mm and length of 600 mm. The natural frequencies and frequency response functions (FRFs) obtained using the present spectral element model and the ANSYS commercial finite element (FE) analysis package [26] are compared in Table 1 and Figure 4, respectively. A one-element model was used to obtain the SEM results. The number of finite elements used in the ANSYS FE analysis was increased step by step to investigate the convergence of the ANSYS results. For the ANSYS FE analysis, eight-node quadratic plane-stress elements were used. Table 1 and Figure 4 show that the ANSYS results converged to the SEM results as the number of finite elements (i.e., the number of DOFs) was increased. Small deviations are probably due to differences in the kinematic assumptions and constitutive relations adopted in the present SEM and ANSYS formulations.Table 1
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the composite beam without delamination.
Mode number
1
3
5
10
SEM (4)
21.498
375.43
1204.3
4155.9
ANSYS (140)
21.502
379.98
1353.6
9313.0
ANSYS (560)
21.499
375.60
1206.2
4193.6
ANSYS (3360)
21.499
375.56
1205.5
4168.8
Note: (  ) is the number of degrees of freedom used in the analysis.Figure 4
Comparison of the frequency response functions (FRFs) obtained by SEM and ANSYS for the composite beam shown in Figure5: n = the number of finite elements used in the analysis.As the second example problem, a cantilevered cross-ply symmetric composite beam[0/90]s with a delamination was considered. It was assumed that the length (L) of the composite beam is 600 mm. The delamination is located in the middle (x=300 mm) of the composite beam, as shown in Figure 5, and it has length lw=5 mm. The natural frequencies obtained using the present spectral element model and ANSYS are compared in Table 2. It is obvious from Table 2 that the ANSYS results converged to the SEM results as the number of finite elements was increased. These observations made from two example problems certainly validate the high accuracy of the present spectral element model.Table 2
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the delaminated composite beam shown in Figure5 (lw = 5 mm, r=0.25).
Mode number
1
3
5
10
SEM (12)
21.482
375.43
1204.3
4155.3
ANSYS (174)
21.504
385.52
1421.9
10490
ANSYS (622)
21.499
375.58
1206.0
4190.6
ANSYS (3422)
21.499
375.55
1205.3
4166.9
Note: (  ) is the number of degrees of freedom used in the analysis.Figure 5
Cantilevered composite beam with a delamination (unit: mm).


### 4.2. Effects of the Delamination Length and Location

To investigate the effects of the delamination length and location on the dynamic characteristics and waves of a composite beam, we reconsidered the cantilevered cross-ply symmetric composite beam[0/90]s shown in Figure 5. The composite beam has the length (L) of 1400 mm and a delamination of length lw is located in the middle (x=700 mm) of the composite beam. The location of the delamination in the thickness direction is represented by r (the ratio of the thickness of the upper part of the delamination to the composite beam thickness). It was assumed that examples of composite beams considered in this paper are all initially at rest.Table4 shows the effects of delamination length on the natural frequencies of the composite beam. Four delamination lengths (lw=0, 5, 50, and 100 mm) were considered and it was assumed that they were located at r=0.25. Table 4 certainly shows that the natural frequencies tend to decrease as the delamination length gets increased.To investigate the effects of the delamination length on the guided waves whenr=0.25, the flexural (bending) and axial (extensional) waves that propagate forward and backward through the upper and lower parts of the delamination were observed at x=500 mm for four values of delamination length (lw). The guided waves were excited by applying an input point force in the form of a Morlet wavelet signal with a 100 kHz center frequency and having five cycles at the free end (x=1400 mm) of the composite beam. The results are compared in Figure 6. The center frequency of the excitation force was chosen to be 100 kHz. Accordingly, wavelengths of the flexural wave w0(x,t) and the axial wave u0(x,t) are 18.6 mm and 74.4 mm, respectively. The axial displacement and flexural displacement are not coupled for a symmetric composite beam. Only the flexural waves are initially excited by the transverse loading. This is why only flexural wave packets were observed at t=412 μs when no delamination existed (i.e., lw=0). On the other hand, Figure 6 shows that axial wave packets can be observed when delamination exists (i.e., lw≠0). Axial wave packets are, in fact, generated by delamination-induced axial-bending coupling. From Figure 6, we note the followings. (1) When the delamination length is larger than the wavelength of a flexural wave, a part of the flexural wave that propagates across the delamination is repeatedly reflected backward and forward within the delamination zone. This seems to result in tail waves that appear after the flexural wave packet at t=412 μs. On the other hand, when the delamination length is smaller than the wavelength of a flexural wave, the aforementioned wave reflections seem to be negligible. (2) The axial wave packets that appear before the flexural wave packets at t=412 μs are generated in part by the mode conversion that may occur via delamination-induced axial-bending coupling when a flexural wave arrives at the delamination and in part by multiple wave reflections within the delamination zone.Flexural and axial waves observed atx=0.5 m versus the delamination length lw.
(a)
(b)Details of the flexural waves that propagate through the upper and lower parts of a delamination located atr=0.25 are shown in Figure 7 for two cases: (1) when the delamination length is 5 mm and (2) when the delamination length is 50 mm. When the delamination lengthis smaller than the wavelength of the flexural wave (i.e., 18.6 mm), wave reflections within the delamination zone are so weak that they are not visible, as shown in Figure 7(a). On the other hand, Figure 7(b) shows that there are strong multiple wave reflections within the delamination zone when the delamination length is larger than the wavelength of the flexural wave. Comparing the wave reflections in the upper part and lower part of the delamination shown in Figure 7(b), we note that the wave reflections seem to be stronger in the thinner (upper) part than in the thicker (lower) part of the delamination.Flexural waves propagating through the upper and lower parts of a delamination that is located atr=0.25.
(a)
Whenlw=5 mm
(b)
Whenlw=50 mmTo investigate the effects of the delamination location in the thickness direction on the guided waves, we assumed that 5 mm long delaminations are located atr=0.25, 0.5, and 0.75. We note that r is defined as the ratio of the thickness of the upper part of the delamination to the total thickness of the whole composite beam. For example, r=0.5 implies that a delamination is located at the midplane of the composite beam. The guided waves to be observed at x=0.5 m were simulated for three different delamination locations and the results are compared in Figure 8. The large wave packets at t=412 μs are the flexural waves that start from the excitation point at the free-end tip and propagate directly to the observation point at x=0.5 m. Figure 8 shows that the flexural wave packets are very similar in both amplitude and phase, regardless of the delamination location. The wave packets observed at t=350 μs are axial waves that are converted from flexural waves via the delamination-induced axial-bending coupling. Axial wave packets can be observed (see Figure 8) when the delamination is located at r=0.25 and 0.75 but not when the delamination is located at r=0.5. This is because the delamination-induced axial-bending coupling diminishes when r=0.5 due to the material and geometric symmetry of the composite beam under consideration. From Figure 8, we also note that axial wave packets when the delamination is located at r=0.25 and at r=0.75 are out of phase with each other. This is because when a symmetric composite beam is subjected to an axial force, the delamination at r=0.25 tends to deform the composite beam downward, while a delamination at r=0.75 tends to deform the composite beam upward, and vice versa.Figure 8
Effects of delamination location (r) in the thickness direction on the guided waves in the symmetric composite beam shown in Figure 5.


### 4.3. Effects of the Lay-Up of Composite Laminates

Figure9(a) compares the flexural waves observed when a delamination of length lw=50 mm is located at r=0.25 to those observed when no delamination exists in a symmetric composite beam as shown in Figure 5 for four different lay-ups of composite laminates: (1) [0/90]s, (2) [30/90]s, (3) [45/90]s, and (4) [45/90]s. Table 3 shows that, regardless of which part of the delamination a flexural wave propagates through, the group velocity of the flexural wave tends to decrease as the fly orientation angle α is increased from 0° to 60°. Thus, Figure 9(a) shows that the first arrival of the flexural waves is delayed longer at the larger value of α. Both the wavenumber dispersion curves and the group velocity dispersion curves shown in Figures 9(b) and 9(c), respectively, demonstrate that the flexural waves at the larger value of α are more dispersive in the frequency range of excitation. As a result, Figure 9(a) shows that the widths of flexural wave packets at the larger value of α are larger than those at the smaller value of α. Compared with the flexural waves observed when no delamination exists in the symmetric composite beam (blue solid lines), the flexural waves observed when a delamination exists (red dotted lines) are found to have tail waves that appear after the main flexural wave packets. The tail waves are generated via the multiple wave reflections within the delaminated zone, as discussed previously. Table 3 shows that the difference between the group velocity of a flexural wave that propagates through the upper part of the delaminated zone and through the lower part of the delaminated zone becomes larger at a larger value of α. Faster waves that travel through the upper part and slower waves that travel through the lower part of the delaminated zone combine to form the flexural waves observed at x=500 mm. This is why more complex tail waves appear at larger value of α.Table 3
Group velocities (m/s) of flexural waves versus the fly orientation angleα for the lay-up [α/90]s.
α (degree)
Part without delamination
Part with delamination
Upper part
Lower part
0
2178
2131
2194
30
2087
1531
1934
45
1968
1382
1845
60
1912
1326
1811Table 4
Natural frequencies (Hz) obtained by SEM for the delaminated composite beam shown in Figure5 (lw = 5 mm, r=0.25).
Mode
Delamination length (mm)
0
5
50
100
1
3.9496
3.9483
3.9366
3.9230
3
69.244
69.244
69.198
68.887
5
223.90
223.90
223.42
220.29
10
989.25
987.91
977.62
962.29Effects of lay-up of composite laminates on flexural waves: (a) flexural waves, (b) wavenumber dispersion curves, and (c) group velocity dispersion curves.
(a)
(b)
(c)


## 4.1. Validation of Spectral Element Model

To evaluate the accuracy of the present spectral element model, a cantilevered cross-ply symmetric composite beam[0/90]s without delamination was considered as the first example problem. The examples of composite beam has thickness of 5 mm and length of 600 mm. The natural frequencies and frequency response functions (FRFs) obtained using the present spectral element model and the ANSYS commercial finite element (FE) analysis package [26] are compared in Table 1 and Figure 4, respectively. A one-element model was used to obtain the SEM results. The number of finite elements used in the ANSYS FE analysis was increased step by step to investigate the convergence of the ANSYS results. For the ANSYS FE analysis, eight-node quadratic plane-stress elements were used. Table 1 and Figure 4 show that the ANSYS results converged to the SEM results as the number of finite elements (i.e., the number of DOFs) was increased. Small deviations are probably due to differences in the kinematic assumptions and constitutive relations adopted in the present SEM and ANSYS formulations.Table 1
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the composite beam without delamination.
Mode number
1
3
5
10
SEM (4)
21.498
375.43
1204.3
4155.9
ANSYS (140)
21.502
379.98
1353.6
9313.0
ANSYS (560)
21.499
375.60
1206.2
4193.6
ANSYS (3360)
21.499
375.56
1205.5
4168.8
Note: (  ) is the number of degrees of freedom used in the analysis.Figure 4
Comparison of the frequency response functions (FRFs) obtained by SEM and ANSYS for the composite beam shown in Figure5: n = the number of finite elements used in the analysis.As the second example problem, a cantilevered cross-ply symmetric composite beam[0/90]s with a delamination was considered. It was assumed that the length (L) of the composite beam is 600 mm. The delamination is located in the middle (x=300 mm) of the composite beam, as shown in Figure 5, and it has length lw=5 mm. The natural frequencies obtained using the present spectral element model and ANSYS are compared in Table 2. It is obvious from Table 2 that the ANSYS results converged to the SEM results as the number of finite elements was increased. These observations made from two example problems certainly validate the high accuracy of the present spectral element model.Table 2
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the delaminated composite beam shown in Figure5 (lw = 5 mm, r=0.25).
Mode number
1
3
5
10
SEM (12)
21.482
375.43
1204.3
4155.3
ANSYS (174)
21.504
385.52
1421.9
10490
ANSYS (622)
21.499
375.58
1206.0
4190.6
ANSYS (3422)
21.499
375.55
1205.3
4166.9
Note: (  ) is the number of degrees of freedom used in the analysis.Figure 5
Cantilevered composite beam with a delamination (unit: mm).


## 4.2. Effects of the Delamination Length and Location

To investigate the effects of the delamination length and location on the dynamic characteristics and waves of a composite beam, we reconsidered the cantilevered cross-ply symmetric composite beam[0/90]s shown in Figure 5. The composite beam has the length (L) of 1400 mm and a delamination of length lw is located in the middle (x=700 mm) of the composite beam. The location of the delamination in the thickness direction is represented by r (the ratio of the thickness of the upper part of the delamination to the composite beam thickness). It was assumed that examples of composite beams considered in this paper are all initially at rest.Table4 shows the effects of delamination length on the natural frequencies of the composite beam. Four delamination lengths (lw=0, 5, 50, and 100 mm) were considered and it was assumed that they were located at r=0.25. Table 4 certainly shows that the natural frequencies tend to decrease as the delamination length gets increased.To investigate the effects of the delamination length on the guided waves whenr=0.25, the flexural (bending) and axial (extensional) waves that propagate forward and backward through the upper and lower parts of the delamination were observed at x=500 mm for four values of delamination length (lw). The guided waves were excited by applying an input point force in the form of a Morlet wavelet signal with a 100 kHz center frequency and having five cycles at the free end (x=1400 mm) of the composite beam. The results are compared in Figure 6. The center frequency of the excitation force was chosen to be 100 kHz. Accordingly, wavelengths of the flexural wave w0(x,t) and the axial wave u0(x,t) are 18.6 mm and 74.4 mm, respectively. The axial displacement and flexural displacement are not coupled for a symmetric composite beam. Only the flexural waves are initially excited by the transverse loading. This is why only flexural wave packets were observed at t=412 μs when no delamination existed (i.e., lw=0). On the other hand, Figure 6 shows that axial wave packets can be observed when delamination exists (i.e., lw≠0). Axial wave packets are, in fact, generated by delamination-induced axial-bending coupling. From Figure 6, we note the followings. (1) When the delamination length is larger than the wavelength of a flexural wave, a part of the flexural wave that propagates across the delamination is repeatedly reflected backward and forward within the delamination zone. This seems to result in tail waves that appear after the flexural wave packet at t=412 μs. On the other hand, when the delamination length is smaller than the wavelength of a flexural wave, the aforementioned wave reflections seem to be negligible. (2) The axial wave packets that appear before the flexural wave packets at t=412 μs are generated in part by the mode conversion that may occur via delamination-induced axial-bending coupling when a flexural wave arrives at the delamination and in part by multiple wave reflections within the delamination zone.Flexural and axial waves observed atx=0.5 m versus the delamination length lw.
(a)
(b)Details of the flexural waves that propagate through the upper and lower parts of a delamination located atr=0.25 are shown in Figure 7 for two cases: (1) when the delamination length is 5 mm and (2) when the delamination length is 50 mm. When the delamination lengthis smaller than the wavelength of the flexural wave (i.e., 18.6 mm), wave reflections within the delamination zone are so weak that they are not visible, as shown in Figure 7(a). On the other hand, Figure 7(b) shows that there are strong multiple wave reflections within the delamination zone when the delamination length is larger than the wavelength of the flexural wave. Comparing the wave reflections in the upper part and lower part of the delamination shown in Figure 7(b), we note that the wave reflections seem to be stronger in the thinner (upper) part than in the thicker (lower) part of the delamination.Flexural waves propagating through the upper and lower parts of a delamination that is located atr=0.25.
(a)
Whenlw=5 mm
(b)
Whenlw=50 mmTo investigate the effects of the delamination location in the thickness direction on the guided waves, we assumed that 5 mm long delaminations are located atr=0.25, 0.5, and 0.75. We note that r is defined as the ratio of the thickness of the upper part of the delamination to the total thickness of the whole composite beam. For example, r=0.5 implies that a delamination is located at the midplane of the composite beam. The guided waves to be observed at x=0.5 m were simulated for three different delamination locations and the results are compared in Figure 8. The large wave packets at t=412 μs are the flexural waves that start from the excitation point at the free-end tip and propagate directly to the observation point at x=0.5 m. Figure 8 shows that the flexural wave packets are very similar in both amplitude and phase, regardless of the delamination location. The wave packets observed at t=350 μs are axial waves that are converted from flexural waves via the delamination-induced axial-bending coupling. Axial wave packets can be observed (see Figure 8) when the delamination is located at r=0.25 and 0.75 but not when the delamination is located at r=0.5. This is because the delamination-induced axial-bending coupling diminishes when r=0.5 due to the material and geometric symmetry of the composite beam under consideration. From Figure 8, we also note that axial wave packets when the delamination is located at r=0.25 and at r=0.75 are out of phase with each other. This is because when a symmetric composite beam is subjected to an axial force, the delamination at r=0.25 tends to deform the composite beam downward, while a delamination at r=0.75 tends to deform the composite beam upward, and vice versa.Figure 8
Effects of delamination location (r) in the thickness direction on the guided waves in the symmetric composite beam shown in Figure 5.


## 4.3. Effects of the Lay-Up of Composite Laminates

Figure9(a) compares the flexural waves observed when a delamination of length lw=50 mm is located at r=0.25 to those observed when no delamination exists in a symmetric composite beam as shown in Figure 5 for four different lay-ups of composite laminates: (1) [0/90]s, (2) [30/90]s, (3) [45/90]s, and (4) [45/90]s. Table 3 shows that, regardless of which part of the delamination a flexural wave propagates through, the group velocity of the flexural wave tends to decrease as the fly orientation angle α is increased from 0° to 60°. Thus, Figure 9(a) shows that the first arrival of the flexural waves is delayed longer at the larger value of α. Both the wavenumber dispersion curves and the group velocity dispersion curves shown in Figures 9(b) and 9(c), respectively, demonstrate that the flexural waves at the larger value of α are more dispersive in the frequency range of excitation. As a result, Figure 9(a) shows that the widths of flexural wave packets at the larger value of α are larger than those at the smaller value of α. Compared with the flexural waves observed when no delamination exists in the symmetric composite beam (blue solid lines), the flexural waves observed when a delamination exists (red dotted lines) are found to have tail waves that appear after the main flexural wave packets. The tail waves are generated via the multiple wave reflections within the delaminated zone, as discussed previously. Table 3 shows that the difference between the group velocity of a flexural wave that propagates through the upper part of the delaminated zone and through the lower part of the delaminated zone becomes larger at a larger value of α. Faster waves that travel through the upper part and slower waves that travel through the lower part of the delaminated zone combine to form the flexural waves observed at x=500 mm. This is why more complex tail waves appear at larger value of α.Table 3
Group velocities (m/s) of flexural waves versus the fly orientation angleα for the lay-up [α/90]s.
α (degree)
Part without delamination
Part with delamination
Upper part
Lower part
0
2178
2131
2194
30
2087
1531
1934
45
1968
1382
1845
60
1912
1326
1811Table 4
Natural frequencies (Hz) obtained by SEM for the delaminated composite beam shown in Figure5 (lw = 5 mm, r=0.25).
Mode
Delamination length (mm)
0
5
50
100
1
3.9496
3.9483
3.9366
3.9230
3
69.244
69.244
69.198
68.887
5
223.90
223.90
223.42
220.29
10
989.25
987.91
977.62
962.29Effects of lay-up of composite laminates on flexural waves: (a) flexural waves, (b) wavenumber dispersion curves, and (c) group velocity dispersion curves.
(a)
(b)
(c)


## 5. Conclusions

In this study, a frequency-domain spectral element model was developed for symmetric composite beams. The governing equations of motion were derived by adopting Timoshenko beam theory and Mindlin-Herrmann rod theory for flexural waves and axial waves, respectively. By using a spectral element model formulated using the variation method, the following phenomena were numerically investigated.(1)
Delamination-induced axial-bending coupling can generate axial waves even in a symmetric composite beam. The axial waves in a delaminated symmetric composite beam are generated by mode conversion via the delamination-induced axial-bending coupling.(2)
When the delamination length is larger than the wavelength of a flexural wave that is propagating through a delamination zone, the multiple reflections of the flexural wave within the delamination zone were found to be more significant.(3)
The wave reflections are in general more severe in the thinner (upper) part of the delamination zone than in the thicker (lower) part.(4)
When a transverse loading is applied to a symmetric composite beam, no matter where the delamination is located in the thickness direction, the flexural waves are very similar in both amplitude and phase, which is not true for the axial waves.(5)
For a composite beam with a lay-up of[α/90]s, the flexural waves become more dispersive when the fly orientation angle α is increased, with the increase in their group velocities. This makes wave propagation in a delaminated composite beam more complicated.



---

*Source: 956043-2014-01-12.xml*