Damage Identification and Assessment Using Lamb Wave Propagation Parameters and Material Damping in FRP Composite Laminates

A methodology for identify damage in the fiber reinforced polymer (FRP) composite has been proposed in this article. The Lamb wave dispersion theory was used to find the fitted peak frequency and loss less finite element model was used to find the modal frequencies in the composite laminates. The change in modal parameters with respect to undamaged and damaged specimen has been considered for the structural diagnosis. The combined finite element and Lamb wave method has been used to obtain damping parameters. The damping capacity was calculated at higher frequency and smaller amplitudes by using hybrid method. The Lamb waves were generated using ultrasonic pulse generator setup. The proposed method was implemented on FRP laminates (CFRP and GFRP) and the results were compared with bandwidth method.


Introduction
Composite materials with advanced properties like, high specific strength and fatigue resistance are being used for many components of aircraft structures in recent era. However, they have high chances of failure when they are subjected to low velocity-impact, which could lead to barely visible impact damage (BVID). BVID are considered to be internal defects which can lead to catastrophic accidents in service. Damage is defined as the changes introduced into a system that leads to affect adversely to its current or future performance. Damage assessment techniques in the structural dynamics have been divided into three categories such as linear, non-linear and transient vibrational measurements. The change in modal parameters such as natural frequencies and material damping can be considered as the prevalent damage detection methods in structural assessment procedure. The existing damage in a structure, leads to the reduction in stiffness and consequently decreasing of the natural frequencies of the system. According to Doebling et al. [1] damage detection by using vibration measurement in elastomers has been first reported by Lifshitz and Rotem [2]. They used the changes in the dynamic moduli and change in the natural frequencies to detect damage.
Many researchers identified crack depth and location from the dependency of the first two structural Eigen frequencies and presented contour graph. The superposed contour of the frequencies variations between the undamaged and damaged structures is used to identification the damage. The intersection point of superposed contour allows identifying both the crack depth and location [3].
According to Gillich and Praisach [4] the change in damping and the friction between crack surfaces lead to dissipative effects. The advantages of using changes in damping is that the cracks allows changes in natural frequencies due to uncertainties and cause important changes in the damping factor allowing damage detection. Kyriazoglou and Guild [5] predicted damping parameters of GFRP and CFRP laminates by using finite element model. Most of the damping related calculations and experiments were carried out by Berthelot and Sefrani [6][7][8] for various composites using Ritz Method. They have performed damping analysis of composite plate and structures by using this method [9][10][11].
Chen and Gibson [12] have studied damping mechanisms in composites which involves a variety of energy dissipation mechanisms. The vibrational parameters such as frequency and amplitude in fiber-reinforced polymers that depend on energy dissipation mechanisms are studied with nondestructive evaluation.
Many researchers have used nondestructive evaluation (NDE) techniques to characterize the fiber-reinforced composites [13,14]. Energy dissipation theory has been used for measuring damping capacity based on vibration damping method.
Damping measurement techniques often deal with natural frequency or resonant frequency of a system. The experimental setup to find the vibration is categorized as free vibration (or free decay) and forced vibration. Free-free beam technique and the piezoelectric ultrasonic composite oscillator technique (PUCOT) are forced vibration techniques. Dynamic mechanical analysis (DMA) uses these techniques to find damping characteristics of the material. However, the instrument is relatively expensive and it cannot be operated at higher frequency and low amplitudes where more information from the tested materials can be obtained [15].
Guan and Gibson [16] have developed micromechanical models for damping in woven fabric-reinforced polymer matrix composites. Where as many other researchers has published results for continuous FRP composites that show damping characteristics of composite material come from microplastic or viscoelastic phenomena associated with the matrix and slippage at the interface between the matrix and the reinforcement [17,18].
The article presents the methodology to find viscous damping, which is the dominant mechanism in FRP composites vibrating at small amplitudes. A hybrid method employing combined finite element and frequency response has been developed to measure the damping properties of composite material.

Ultrasonic pulse generator
In this work the specimens are carbon fiber/epoxy (CFRP) and glass fiber/epoxy (GFRP) laminates. The laminate contains woven fiber with 12 plays with average epoxy layer thickness of 0.2 mm. Figure 1 shows the specimen with the dimensions 250 Â 50 Â 2 mm in damaged state where a cut (30 Â 2 Â 2 mm) at a distance 120 mm away from free end has been introduced. Figure 2 shows an experimental setup to develop Lamb waves using ultrasonic equipment NDT™ EPOCH 4PLUS. The specimen is supported as cantilever; pitch-catch radio frequency (RF) test method was used. Dual-element transducers (DIC-0408) where one element transmits and the other element receives the burst of acoustic waves generated in the specimen. 1 kHz to 4 MHz frequency range transducers was used and they were placed 80 mm apart from each other. Honey glycerine couplant made by Panametrics has been used to couple the two sensors on to the test specimen. The couplant helps to transmit a normal incident shear wave to propagate across the test piece between the transducer tips. Scanview plus™ software was used to acquire and process the data obtained from the test specimens [19,20].

Generation of optimal Lamb wave
Acoustic impedance of the material plays important role in deciding the selection of transducer frequency, low impedance require lower frequency transducers. The materials such as carbon fiber or glass fiber have low impedance thus low frequency transducers will be used to generate Lamb waves. Whereas metal skin layers have high impedance therefore higher frequencies transducers are used for thinner and metallic layers.  The Panametrics-NDT™ EPOCH 4 PLUS is used to generate acoustic waves and the device is equipped with four channels. Optimal Lamb wave propagation parameters were arrived through calibration of the device using editable parameters as shown in Table 1.
The parameters have been varied one by one and arrived to a conclusion of optimal driving frequency for different specimens. Figure 3 shows the optimal driving frequency of different materials calibrated through ultrasonic pulse generator test setup. With the help fitted peak value and percentage amplitude at constant gain 55 (db) Lamb waves generation frequency is identified for different materials. Figure 4 shows the histogram representation of % amplitude of the waveform with a bin range of 0-820 kHz of pulser [21,22]. The most effective range of frequencies to generate Lamb waves is 140-420 kHz. The frequency range to generate Lamb waves for GFRP is 170-190 kHz where as it is 260-280 kHz for CFRP and for Aluminum (Al) it is 360-380 kHz.

Material properties
The Young's modulus for the specimens is calculated from Eq. 1, the calibrated experimental setup is used to generate Acoustic Emission (AE) on to the test specimen. AE velocities travelling in the material are determining from the instrument EPOCH 4 PLUS [22,23].

Pulser
where V L is longitudinal velocity of AE and V T is transverse velocity of AE traveling in the material respectively, ρ is material density and υ 12 is Poisson's ratio. The material properties of specimens arrived from the experimental setup is presented in Table 2.

Methodology
In this work a hybrid method has been proposed for identify change in damping capacity of a material using combined finite element and Lamb wave method. The process diagram of the hybrid method for the dynamic mechanical analysis is shown in Figure 5.
The group velocity (c g n ), modal frequency (f n ) are determined from Lamb wave model and loss less finite element model respectively.

Lamb wave model for laminated composite plate
Lamb waves can be of two groups, symmetric and anti-symmetric, these waves propagate independently of the other and boundary conditions of the wave equation are being satisfied by both of them for this problem. Actuating frequency relating the velocity of Lamb wave propagation has been derived in the following section. Dispersion curves of Lamb wave in a particular material, which plot the  phase and group velocities versus the excitation frequency given by Dalton et al. [24]. The anti-symmetric Lamb wave solution formulated as seen in Eq. 2: individual laminate stress-strain relationship is given by where σ is normal stress, τ represent shear stress, ε is normal strain and ϒ represent the shear strain. Reduced stiffness components Q ij are defined in terms of the engineering constants as where E 1 Young's moduli in the longitudinal and E 2 Young's moduli in the transverse directions. The major and minor Poisson's ratios represented by ν 12 and ν 21 respectively. The relation between Poisson's ratios in Eq. (4) is given by: A 11 and A 22 , are in-plane stiffnesses of plate and these are obtained by integrating the Q ij across the thickness of the plate [21]. These stiffness values are given as: where plate thickness is represented by "h" and "k" represents each individual lamina. The transformed stiffness coefficients Q 0 ij are defined as where m = cos(θ) and n = sin(θ), the angle θ is taken positive for counterclockwise rotation and it is considered from the primed (laminate) axes to the unprimed (individual lamina) axes. Q 0 ij for the 0°and 90°laminas are given by The extensional plate mode velocity is related to the in-plane stiffness of a composite [14]. A 11 and A 22 are the stiffnesses propagating in the 0°and 90°d irections respectively. The relation between extensional plate mode velocity and stiffness is given by: for 90°direction c l ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi The inplane stiffnesses A 11 and A 22 are calculated using Eqs. (4)-(8) by substituting engineering stiffnesses of the composite. The extensional plate mode velocities are substituted into Eq. (2) and it is solved numerically for phase velocity in Mathematica TH . Phase velocity (c phase ) is the dependent variable being solved for the independent variable being iteratively supplied is the frequency-thickness product, where ω is the driving frequency in radians. Group velocity dispersion curve, which are derived from the phase velocity curve using Eq. (11): where f is the frequency in Hz.

Finite element model for free vibration of a laminated composite plate
Natural frequency is the phenomenon that occurs with oscillatory motion at certain frequencies known as characteristic values, and it follows well defined deformation pattern known as mode shapes or characteristic modes. The study free vibration is important in finding the dynamic response of elastic structures. It is assumed that the external force vector P ! to be zero and the harmonic displacement as: (12) and the free vibration is given by: where Q where [k] is stiffness matrix and [M] is mass matrix, which are derived through finite element formulation. Figure 6 shows the plate bending formulation where x, y, and z describes the global coordinate of the plate whereas u, v, and w are the displacements, h represents plate thickness. The xy plane is parallel to the midsurface plane prior to deflection. The displacements in the plate at any point is expressed as The plane displacement u and v vary through the plate thickness as well as with in the xy-plane while the transverse displacement w remains constant through the plate thickness. In order to develop shape functions two different interpolations are used one interpolation within the xy-plane and the other in the z-axis. For the xy-plane interpolation, shape function N i (x,y) are used where subscript i varies depending on the number of nodes on the xy-plane. Shape function H j (z) is used for interpolation along the z-axis, where subscript j varies depending on the number of nodes along the plate thickness. Since two inplane displacement are functions of x, y, and z, both shape functions are used while the shape functions N i (x,y) was used for transverse displacement. The mapping of ξ, η-plane onto xy-plane and ζ-axis to zaxis, was done using isoparametric element and the three displacements are expressed as w ¼ where N 1 represents the number of nodes in xy-plane (ξ, η-plane) and N 2 represents the number of nodes in z-axis (ζ-axis). The first subscript for u and v denotes the node numbering in terms of xy-plane (ξ, η-plane) and the second subscript indicates the node numbering in terms of z-axis (ζ-axis). Four-node quadrilateral shape function is considered for the xy-plane (ξ, η-plane) interpolation, i.e., N 1 = 4 and N 2 = 2 that is linear shape function which is considered for the z-axis (ζ-axis) interpolation. Nodal displacement u i1 and v i1 are displacement on the bottom surface of the plate element and u i2 and v i2 are displacement on the top surface. As seen in Eqs. (18)- (20), there is no rotational degree of freedom for the present plate bending element were as both bending strain energy and transverse shear strain energy are included.
The relation between bending strains and transverse shear strain with respect to displacements is given by: where ε b f g and ε s f g are the bending strain and transverse shear strain respectively. The normal strain along the plate thickness ε z is not considered.
Substitution of Eqs. (18)- (20), into the Eqs. (25) and (26), with N 1 = 4 and N 2 = 2 gives: where The constitutive equation is For the bending components where where Eq. (33) is for the plane stress condition for the plate bending theory and for a FRP composite, is given by In which Here, the longitudinal direction is represented with 1 and transverse direction is represented with 2 for the FRP composite. Further E i and G ij are the elastic modulus and shear modulus respectively, whereas υ ij is Poisson's ratio for strain in the jdirection. Five independent material properties will be considered for Eqs. (37)-(42) because of the reciprocal relation The stiffness matrix of the element k e ð Þ h i is expressed as: where Ω represents the plate domain.
Similarly the mass matrix is given by where A is area of the element, t is the element thickness and ρ density of material. Natural frequencies are arrived for composite plates from the lossless finite element formulation. The Mathematica TH software has been used to compute the Eigen values using the inputs taken from the experimental data discussed in the previous sections. The analytical model developed was correlated ANSYS model. Block Lancozs method was used to carry out modal analysis in ANSYS. The analysis was done for 30 subsets and shell-190 has been used as meshing element. Figures 7  and 8 shows the first and twentieth mode of natural frequency of the undamaged specimen, i.e., GFRP and CFRP respectively and similarly Figures 9 and 10 shows for damaged specimen.

Bandwidth method
The damping parameters in FRP composites are based on the energy dissipation mechanism. Vibrational parameters such as frequency and amplitude are used to determine the dynamic characteristics of a system. The best practice to study vibrational parameters is with nondestructive evaluation. Damping characteristics of a system can be determined by the maximum response, i.e., the response at the resonance frequency as indicated by the maximum value of R v . Figure 11 illustrates the Bandwidth method of damping measurement where, damping in a system is indicated by the sharpness or width of the response curve in the vicinity of a resonance frequency ω r , designating the width as a frequency increment (i.e., Δω=Δω 2 -ω 1 ) measured at the "half-power point" (i.e., at a value (R= ffiffi ffi 2 p )) and the damping ratio ζ can be estimated by using band width in the relation given by   for i th mode damping ratio is given by The equation of motion of a system with viscous damping, when the excitation is a force F ¼ F o sin ωt applied to the system, is given by Eq. (48) represents the forced vibration of a damped system and the resulting motion occurs at the forcing frequency ω. The damping coefficient c is greater than zero, leads to change in the phase between the force and resulting motion. The phase change is termed as phase angle δ which is a function of the frequency ratio ω/ω r and for several values of the fraction of critical damping ζ, given by [25].

Results and discussion
The phenomenon of change in modal parameters has been used to identify the damage in the specimens. The damage in the specimen is identified by change in damping capacity with respect to undamaged specimen. The first order Lamb wave equation is used to determine the storage modulus. Lamb wave propagating is quite complex to understand, i.e., an increase in modulus slightly speeds the wave velocity. An increase in the density would have the opposite effect slowing wave velocity, as it appears in all the same terms as the modulus but on the reciprocal side of the divisor.
The AE velocities of the specimens were arrived experimentally using ultrasonic pulse generator test setup and the engineering constants, Young's modulus and poison's ratio were calculated. The engineering constants are substituted in Lamb wave model discussed in previous sections for finding dispersion characteristics shown in Figure 12. The same material properties are used for finite element model to determine natural frequencies.
The group velocity (c g n ) at natural frequency (f n ) and thickness (h) is substituted in Eq. (50) to determine the phase shift and thus finding material damping capacity Tan δ ð ). Dynamic mechanical analysis can be carried out using the same procedure by getting the E o value from group velocity dispersion at iteratively supplied frequencies.
δ ¼ 2πf n h=c g n (50) Damping is the term used in vibration and noise analysis to describe any mechanism whereby mechanical energy in the system is dissipated. The damping properties of so-called damping materials, such as elastomeric materials, are usually temperature and frequency dependent, so the experimental determination of damping material properties requires a long and repeating process.
In dynamic mechanical analysis damping measurements is done in temperature sweep mode whereas in this work frequency sweep mode is used. In the present work damping measurements were carried out using combined finite element and Lamb wave method and the results were compared with bandwidth method.
The modal analysis was carried out using developed finite element model and it was correlated with ANSYS. The waveform from the instrument is processed through virtual controlling software and the continuous waveform is subjected to fast Fourier transform (FFT) which yield a single peak from the calibrated optimal driving frequency, however for a few finite cycles, the FFT appears as a Gaussian curve. The response curve of the undamaged and damaged specimens being tested for damping capacity using bandwidth method is shown in Figure 13.
The Lamb wave dispersion curves have been obtained from the iterative supply of the frequency using Mathematica TH code. The group velocity dispersion curve of the specimens used in this research is shown in Figure 12. The group velocities and the natural frequencies obtained from modal analysis are used to determine damping capacity at various mode of interest. Table 3 shows the damping capacities of the undamaged specimen in comparison at critical modes similarly for damaged specimen it has been reported in Table 4. It is observed that the natural frequencies of the damaged specimen fell down and the damping capacities have increased slightly with respect to undamaged specimens. Figure 14 shows the damping capacities and dynamic  storage modulus of the tested specimens with respect to their natural frequencies.
The material GFRP and CFRP exhibits similar damping property to a certain range of frequency, and in between 2 and 8 kHz GFRP has better damping property among the two and at higher range of frequencies CFRP is found to be good in damping characteristics.

Conclusions
Dynamic mechanical analysis is a technique used to study and characterize damping behavior of materials. It is most useful for studying the viscoelastic behavior of polymers. The tests were conducted on polymer composites CFRP and GFRP laminates in their undamaged and damaged state. A hybrid method has been explored in this work and the materials have been characterized for damping parameters at their mode frequencies. The change in the modal parameters (i.e., natural frequencies and damping capacity) can be used to identify and assesses the health of the structures. It is very advantageous method to obtain damping characteristics of the materials at higher frequency and at relatively low amplitudes.