Finite element evaluation of the effects of curvature in Lamb waves for composites structural health monitoring

Abstract

This work presents the effect of the high curvature to thickness ratio on the characteristics of Lamb Waves propagating over the skins of composite structures. It is accessed how the curvature on composite skins affects the group velocity of the symmetrical (S0) and asymmetrical (A0) wave modes. It is also accessed the gradient of curvature effect, when the wave propagates from a flat to a curved part on the skin. The results are intended to be used for the improvement of structural health monitoring of wings and wind turbine blades. This is accomplished through dynamic explicit linear finite element method simulations of plates with flat and curved parts made of Eglass-epoxy bidirectional laminate. As the skin structures are often designed to withstand torsional loads, the fibers are aligned with 45 degrees from the leading edge (curved region) throughout the analysis. Several skins with different curvature to thickness ratios are generated and simulated. Results and trends are presented and can be used to improve the algorithms for damage detection on those structures. It is shown that the group velocity of the incoming waves change considerably with the presence of curvature, for both main modes of vibration.

Keywords
TDFEM; Lamb waves; curvature; glass fiber; composites; GFRP; SHM

1 INTRODUCTION

The end of the 20^th century and the beginning of the 21^st has seen an exponential growth in composite structures usage. The aerospace and marine industries dominated the former era. But on the latter era, the wind energy industry increased a lot the demand for composite materials and is now the main user (Ellabban et al., 2014). On its roadmap that set the framework for the cost reduction and future increase in wind energy usage, the International Energy Agency stated that operation and maintenance costs (O&M) need to reduce drastically (IEA, 2013). According to the agency, due to potential catastrophic failure aftermaths and inherent operational hazards, the wind turbine blade skins need to be inspected periodically. That poses additional difficulties and costs, especially for offshore wind turbines, where the access for workers is extremely complex and depends on good weather. Because of those difficulties, among others, the IEA suggested the usage of real time predictive maintenance.

There are some ways to perform real time damage prediction on the wind turbine blade structure. Some rely on the deformation measurement at some points at the blade (Schubel et al., 2013). Some of them rely on the measurement of vibration characteristics of the blades, comparing the modal analysis of a baseline with the real vibration of the blade over time (Lading et al., 2002 and Fan & Qiao, 2011). A big part of the predictive maintenance in the field is done with acoustic emissions (AE) on the surface of the wind turbine blades. Some of the methods rely on capturing the waves in the failure moment, which might be generated by an impact, yield, crack, or other events (Sengupta et al., 2015). In this method, it is possible to perform the predictive maintenance only with receptor sensors, without the need of actuators. Other methods utilize source actuators and are capable of pinpointing the damages at any moment, not only at the moment they happen. Through an array of actuator and receptors (mostly piezoelectric - PZT), it is possible to generate surface elastic waves and detect changes on the signature, from a baseline set of data. These changes in signature can be used to identify surface damage, existence, location and size, performing, in real time, what is known as Structural Health Monitoring (SHM). The longitudinal and bending surface acoustic waves are called Lamb Waves, as the wave equations were deduced by Horace Lamb, for isotropic infinite flat plates (Lamb, 1917). The field of structural health monitoring that uses Lamb Waves is called LW-SHM.

As discovered by Lamb (1917) and observed on future experimental studies, the Lamb waves are dispersive by nature, which means the velocity of each wave component varies with its source frequency. Also, there are two main modes of vibration with its harmonics: the symmetric “S” mode, which is the surface acoustic wave that propagates radially and also vibrates radially (also called longitudinal wave) and the asymmetric “A” mode, which is the surface acoustic wave that propagates radially and vibrates perpendicularly to the plate (also called bending wave). Due to the inherent physics, generally the A wave is slower than the S wave, and for orthotropic materials (Daniel & Ishai, 1994) such as fiberglass, the wave velocity is very dependent on fiber direction for the S wave, whereas not so dependent of fiber direction in the case of A wave. The A wave has considerably bigger amplitudes for smaller travel distances (order of magnitude of 0.1 m travel distance) than the S wave, although the latter type has smaller attenuation over distance (Ono & Gallego, 2012). It has been shown that when a failure is present in the surface of a plate, the failure acts as a source, also known as scatter signal (Wang et al., 2004). Hence, knowing the group velocity of the wave, it is possible to analyze the arriving signal to receptor sensors, considering the signal originated from the source, traveled until the failure and traveled from the failure to the receptor. The total time of travel is called time of flight (TOF). Capturing the signal in real time and subtracting from a baseline, it is obtained a scatter signal, which, together with the TOF concept, can be used to obtain the position and size of damages, even for many damages happening at the same time (Qiu et al., 2013).

The framework done by Qiu et al. (2013) represents a relatively simple and effective way to perform LW-SHM, but it is necessary to have information regarding the wave. It is necessary to know the group velocity of the wave in any direction, as well as if there is a change in the group velocity due to geometry. The group velocity is the atribute that is normally used for TOF based LW-SHM, because it is the velocity of the modulation of the wave, which contains the wave packet energy. As pointed out by Sengupta et al. (2015), there is a lack of studies regarding the effect of geometry on the Lamb Wave behaviour. Most of the studies are done for flat plates, including the work of Qiu et al. (2013), they do not address the effect of high curvature in the Lamb waves and, more importantly, in the group velocities. This high curvatures are present in wings and wind turbine blades, in the leading edge region of the airfoil. Furthermore, no studies are made for the abroupt curvature change (gradient) near those leading edges, specifically for the LW-SHM on composite material structures. There are some articles that access the curvature changes in Lamb waves, though. Marston (1989) presents the mathematical treatment of the phase velocity of Lamb waves in a spherical shell, showing the dependence with curvature. Towfighi et al. (2002) present the mathematical treatment of Lamb waves traveling in the circunferential direction on hollow anisotropic cylinders. Nazeer et al. (2017) present a damage detection analysis in a plate with multiple bends (top hat stiffener) using the shear surface waves, also called Love waves. They assessed the sensitivity of the waves to find damages in isotropic studies for varying frequencies and had good agreement with finite element method simulations.

The importance of these works and the ones before them are paramount for the field of SHM, but there is still a lack of understanding of the characteristics of group velocity for Lamb waves in orthotropic bent plates, both inside the curved region and around it. The present work aimed at solving this issue, so that future TOF algorithms for LW-SHM can be used more generally, with more accuracy and efficiency, using the underlying physics of the wave. In order to achieve this goal, we created a simplified 3D model of a bent plate and used time domain finite element analysis (TDFEM) and analyzed the behavior of Lamb waves for many values of leading edge curvature.

2 METHODOLOGY

2.1 Choosing the most suited TDFEM simulation of waves

In order to perform the TDFEM simulation of the structure, first we had to choose between two types of integration methods. One is the explicit method and the other is the implicit method. The explicit method is the regular “time marching” method, were the solutions (in this case nodal DOF displacements) for all nodes on the next step are determined based on the previous steps conditions. In the implicit method, the solution for the next time step involves not only the solution of the previous time steps, but also unknown solutions. The latter one involves a matrix solution on time, whereas the former one doesn´t need the matrix solution (He et al., 2012). For the same time increment of the simulation and same element sizes, the explicit method is more computational efficient (at least one order of magnitude), because it avoids solving for the matrix calculation over time. However, if no detail regarding the wave propagation is needed and bigger time steps can be used (e.g.: quasi steady state solution) the explicit method starts having stability problems (Jiao & Jin, 2001) and the implicit method becomes more advantageous. Usually, explicit analysis is very useful for issues where fast loading happens, including explosions, impact or progressive failures (Rahnavard et al., 2018 and Bettinotti et al., 2017). Additionally, implicit dynamic analysis is applied to quasi-static problems such as cyclic pushover analysis, force-control, displacement-control, and thermal analyzes (Lee et al, 2017, Rahnavard et al., 2017).

For the case of the present work, it is needed to know the details of the transient part of the simulation and the wave front forming on the receptor points with adequate precision in space and time. Hence, the explicit method was used. When the finite element method with explicit analysis is used to simulate the behavior of elastic waves, the simulation parameters of time increment and element size must be carefully chosen. This is necessary to guarantee five aspects: i) numerical stability; ii) wave time accuracy; iii) curvature geometrical accuracy; iv) Mesh quality and v) time of flight (TOF) accuracy between analysis points. The last one is important for the simulations of the present work, where it was analyzed the difference in time of flight and, hence, wave speed between different points in the geometry, to be used in SHM algorithms. Each one of the five aspects must be attained without making it impractical, in terms of computational time, to simulate many curvatures. Figure 1 shows a chart with all five aspects and how element size and time increment are related to them.

Figure 1
Element size and time increment magnitude effects.

In figure 1, curvature accuracy is the ability to represent the curved part of the geometry. The smaller the element size, the most accurate is the curvature geometry modelled and the bigger it is the curvature accuracy. Wave time accuracy is the ability to represent the sine wave over time accurately. Too large time increments compared to the wave period decrease wave time accuracy. Time of flight accuracy is the ability to measure changes in wave arrival time between source and receptor. The smaller the time increment, the bigger this accuracy. For mesh quality, it is considered that there is only one element on the plate thickness. Hence, the element size in figure 1 is related only with the surface dimension of the hex element, since the thickness is fixed. Smaller values of element size increase the element aspect ratio, decreasing mesh quality. But given that the minimum mesh requirements are attained, the smaller the element size, the most accurate is the solution. Referring to figure 1, it can be concluded that the best option is to use relatively smaller values of time increment, but not necessarily decreasing element size to its minimum.

Regarding numerical stability, it is important to guarantee the CFL or Courant-Friederichs-Lewy condition, following the works of Courant et al. (1967) and others (Sanz-Serna & Spijker, 1986). This condition was devised for finite difference numerical solutions of hyperbolic partial differential equations, such as those governing elastic waves on solids. According to the authors’ knowledge, if regular time marching explicit analysis is done, the solution “locks” or find a singularity, whenever the wave travels more than one finite difference in any given time increment. This condition can be better understood visualizing figure 2, depicting one dimensional problems.

Figure 2
CFL condition for one-dimension hyperbolic equations.

It was found, by trial and error, that the same behavior can be extended to TDFEM of Lamb Waves and the following equation applies:

$\frac{u△t}{L}<Cmax.$ (1)

The value of “u” is the wave speed in the inspected direction, “$△t$” is the time increment of the simulation and “L” is the element size in the direction of wave travel, “Cmax” is the CFL criterion, which needs to be smaller than 1 for explicit analysis.

Another aspect of TDFEM of waves is that the boundary conditions are “felt” only when the wave reaches the boundary and returns to a given node. In the case of the present work, the scope of the analysis is the first incoming wave and the results of it slightly after it reached the history points. Hence, the type of boundary conditions doesn´t matter for the solution.

2.2 FEM simulation

In order to analyze the effects of curvature in the Lamb waves for curved structures, a baseline geometry was generated, as it is shown in figure 3.

Figure 3
Baseline geometry for simulation.

Numerical simulations were performed using commercial software Abaqus^® (Abaqus, 2014a). A Python-based script was written to do the numerical analysis, aiming to easily change the curvature radius and store desired output parameters (Abaqus, 2014b). For the baseline geometry shown in figure 3, for different external radii, “r”, the simulation was run and history values as well as field variables were recorded. The source point and the History points were located in the positions shown in figure 4. These history points were used to characterize the behavior of the wave speed before, along, and after the curved region, as well as in the plain section (between source and plain receptor).

Figure 4
Source point and history points.

The python script that we used can be seen in detail in the ^{Supplementary Material}. We chose the source excitation frequency of 50 kHz, the same used in the work of Qiu et al. (2013). This is the smallest ultrasound frequency that we saw been used in experimental studies of Lamb Waves in composite materials.

The material for the simulation is a 4-ply bidirectional laminate of glass fiber reinforced polymer (GFRP) with epoxy matrix. The material orientation is such as the fibers are +45 degrees and -45 degrees with the leading edge (+45 and -45 degrees rotation through the x axis), with symmetric laminate scheme. The whole simulation was done with the compatible units: N, mm, MPa, tons (1e3 kg). The lamina elastic properties were determined considering a fiber volumetric ratio of 0.41, which is normally obtained in our laboratory through vacuum bagging without autoclave. Table 1 presents the material properties of the lamina.

Table 1
Lamina properties for E-glass/epoxy and 41% fiber volume ratio.

Eglass lamina (41% vf)
Property	Value
E1 (MPa)	2.90E+04
E2 (MPa)	6.39E+03
G12 (MPa)	2.82E+03
ν12	0.28
ν21	0.06
Density (Tons/mm3)	1.51E-09
Thickness (mm)	0.12

In order to start accessing the total simulation time, as well as the time increment, we used the wave velocity equations determined by Wang & Yuan (2007). From the G12 and density properties, we determined the shear wave velocity. With the shear wave velocity, the plate thickness and excitation frequency, we used the charts of dimensionless velocities versus dimensionless frequencies from Wang & Yuan (2007) to estimate the phase and group velocities for the A0 and S0 modes. The faster wave had an order of magnitude of about 4 km/s and the slower wave had a value of 1 km/s using the numbers from this paper. It was just a first guess, since there was no data for the (+45/-45)_s laminate with waves on 45 degrees. With the estimate of velocity at hand, the size of the model needed to be determined, in order for us to access the total simulation time. Through trial and error, visualizing some field results for one leading edge curvature, we determined that the edge with 800 mm (1600 mm plus the curved region on lateral perimeter) or more could sufficiently separate the A0 and S0 waves in space. Also, a value of time increment of 140 ns was enough to get a good sine wave peak resolution (140 points per period) in order to differentiate the arrival of the wave before and after the bend. Hence, the total time of 1.155 ms was deemed enough to visualize the S0 and A0 wave on the history points. For the element size, it was used the curvature criterion h/L of 2.5%. Where h is the maximum distance from the element edge to the arc geometry and L is the element base edge dimension. The curvature criterion meant that the bigger the radius, the bigger the element base size. The element size was varied for many values of radii but limited in 2 mm for bigger values (10 mm radius or more). Using this curvature criterion, the estimated maximum wave velocity and the time increment, a minimum value of 3 mm for the leading edge radius could be simulated inside the Courant Fredrich Lewis condition, as explained on section 2.1.

For the geometry creation, the distance between the source point to the receptor at the beginning of the radius was the same as the distance from the receptor at the end of the radius to the curved receptor (most distant receptor).

The type of element used is the continuum shell SC8, with second order accuracy and the laminate was assembled using each lamina and the continuum shell option (Abaqus, 2014c).

A concentrated load of 0.1 N maximum amplitude was applied acting at the y direction in the source point. The amplitude variation followed a sin function with 50 kHz frequency and a Hann window function (1 + cosine) with three positive peaks, forming a wavelet. Figure 5 shows the load over time. Pinned-type boundary conditions were applied on the four points, as shown in figure 6.

Figure 5
Loading profile applied on the source point.

Figure 6
Loading and boundary point positions.

The total number of elements were around 800,000. Figure 7 shows the schematic of the discretized problem. For the output, we chose to display the 3 displacements for all time increments on history points (source and 4 receptors). For the field results, all node displacements were recorded just for a few key time increments as will be discussed in the results.

Figure 7
Schematic of the discretized problem.

After running the problem, a Matlab^® code was used to interpret the output text file (containing parameter values with time) and calculate the speed of the wave between each set of history points, both for the S₀ and the A₀ modes.

3 RESULTS

3.1 Field Results

First, we start accessing the simulation from the field output results. It was created a time table comprising more time resolution when the wave is close to the radius and less resolution far from the radius. This approach helped to avoid too big result files generation. The field result file for each radius amounts for 20 GB, following this framework. Without following this framework, it could go as much as 200 GB per radius. Figures 8 to 15 show the field results. In those figures, it was used a scale magnification factor of 8e5 for the total displacement and a saturation of the color bar at 8e-9 mm (except for figure 9, where it was used a saturation of 5e-6 mm). Those values were found through trial and error aiming for a better visualization of the results. The objective of the color bar saturation was to differentiate the S0 mode (symmetrical - lesser amplitude) from the A0 mode (asymmetrical - larger amplitude). The wave that is depicted in the colored part has a smaller amplitude and is faster, which are characteristics from S0 modes. The waves that appear with gray color and shadows (saturated color bar) have bigger amplitude and appear latter, which are characteristic from A0 waves.

Figure 8 shows the points used to calculate the S0 mode wavelengths at the 0 degrees, 45 degrees and 30 degrees directions. The values of wavelength obtained for the S0 mode were 56 mm (0 degree), 62.5 mm (30 degrees) and 67.1 mm (45 degrees, fiber direction).

Figure 8
Query points for obtaining S0 mode wavelengths in two directions (frontal view). 1258^th time increment.

Figure 9 shows the analogous points for the A0 mode. In the case of the A0 wave, the wave front is slightly deformed. The wave is faster on the direction of the lower right side of the picture (+45 degrees rotating from x axis). That happens because we used symmetrical laminates, which means that the most outward laminas are aligned in the same direction (+45 degrees) and, thus, the bending stiffness (and phase velocity) is slightly higher on that direction, because of higher inertia. This does not happen for the S0 mode, because the laminate longitudinal stiffness does not change with the lamina order. In real world applications, the bidirectional laminates present the fiber crimp and this effect for A0 will be averaged and nearly equal between +45 and -45 degrees. Hence, we performed the average of values of figure 9 to obtain the wavelengths. Those values from the averages for A0 were 9.6 mm (0 degree), 9.75 mm (30 degrees) and 10.2 mm (45 degrees). From the work of Wang & Yuan (2007), the ratio of S0 phase velocity to A0 phase velocity obtained for a [-45₆/45₆] carbon fiber laminate in the 30 degrees direction was 6.6. The phase velocity is calculated multiplying the frequency of excitation by the wavelength. Hence, for the same frequency, the phase velocity ratio is the same as the wavelength ratio. Through inspection of figures 8 and 9, for 30 degrees direction, we obtained the value of 6.41 for the same ratio. This value is within 3% of the value found by Wang & Yuan (2007). Even though the materials are different, comparing the ratios is an accurate way to check for consistency and accuracy of the simulation.

Figure 9
Query points for obtaining A0 mode wavelengths in two directions (frontal view). 972^nd time increment.

Another way to check the simulation for physical consistency is through the calculation of longitudinal S0 phase velocities. It is known that the phase velocity (Cp) of the fundamental longitudinal mode (acoustic wave) equals:

$Cp=\sqrt{\frac{Q}{\rho }},$ (2)

where “Q” is the longitudinal stiffness of the laminate on the direction investigated and $\rho$ is the laminate density. In order to find “Q”, it is necessary to calculate the stiffness matrix of the lamina and then calculate the laminate stiffness by rotating each individual lamina matrix forming the laminate stack. Following Daniel & Ishai (1994), we have the following elements for the stiffness matrix for each lamina:

$Q11=\frac{E_1}{1-{\vartheta }_{12}.{\vartheta }_{21}},Q22=\frac{E_2}{1-{\vartheta }_{12}.{\vartheta }_{21}},Q12=\frac{{\vartheta }_{21}{E}_1}{1-{\vartheta }_{12}.{\vartheta }_{21}}andQ66=G_{12}.$ (3)

The values of material properties used on equation 3 are present on table 1. For the bi-directional laminate of the present simulation, which is [-45/45/-45/45]_s, the -45 and 45 degrees laminate stiffness will be an average of the values of Q11 and Q22. The 0 and 90 degrees laminate stiffness will be a 45 degrees rotation for the lamina matrix. According to Daniel & Ishai (1994) the value of the rotated longitudinal stiffness is:

$Qxx=m^4.Q11+n^4.Q22+2.m^2.n^2.Q12+4.m^2.n^2.Q\mathrm{66,}$ (4)

where “m” is the cosine of the rotated angle (45 degrees for this case) and “n” is the sine of the rotated angle. The results for the computation of equation 4 and the values of theoretical phase velocities are found on table 2 below.

Table 2
Values for Phase velocity calculation versus fiber angle.

	Laminate angle (degrees)
Value	0	45
Q11 (Gpa)	29.5	-
Q22 (Gpa)	6.50	-
Q12 (Gpa)	1.77	-
Q66 (Gpa)	2.82	-
Qxx (Gpa)	12.7	18.0
r (kg/m3)	1.51E+03	1.51E+03
Cp calculated (m/s)	2.90E+03	3.45E+03
Wavelength S0 (mm)	56.0	67.1
Frequency (kHz)	50.0	50.0
Cp from simulation	2.80E+03	3.36E+03
Difference (%)	3.6%	2.9%

The differences between the simulation values and theoretical values of phase velocity for both directions are within 4%, which could be due to the discretization of the simulation and due to the graphical method error. But we find them sufficient to attest the simulation consistency and accuracy for engineering usage.

Figure 10 shows the first incoming wave (S0) right before reaching the bend at the leading edge for three values of radii. No apparent change in the wave front is detected by the solely analysis of the field graph. Figures 11 through 13 show the S0 wave reaching the leading edge for three values of radii. It is possible to see a distortion of the wave front at the leading edge, especially for the radius value of 6 mm and gradually less distortion as the radius increases. That is a qualitative evidence of wave speed change at the leading edge bend, causing refraction of the wave.

After the S0 wave passes through the leading edge bend, there is a hiatus until the A0 wave arrives. That portion is depicted in figure 14, for three values of radii. Here, some interesting patterns arise. Another lower amplitude mode of the wave can be seen, it generates near elliptical wave fronts, with the main axes aligned with the -45 and +45 degrees directions of the laminate. At 0 and 90 degrees direction, the two elliptical modes interfere with each other, creating different patterns on these directions. We hypothesize that these are each lamina wave front for the longitudinal (symmetrical) vibration. It is possible to observe, also, that two small wavelength waves appear from the leading edge in the format of fringes, as depicted in the right upper portion of figure 14. The leading edge acts as a source and those fringes propagate away from the leading edge bend on both senses. This is analogous to the Scatter waves that are generated by the presence of a material failure, delamination or hole in the structure. As the leading edge perimeter increases to numbers close to the S0 wavelength, those fringes begin to disappear, as can be seen on the right lower picture of figure 14.

Figure 10
Total displacement at the 472nd time increment (isometric view). Clockwise: 6 mm, 16 mm and 36 mm radii.

Figure 11
Total displacement at the 600th time increment (isometric view). Clockwise: 6 mm, 16 mm and 36 mm radii.

Figure 12
Total displacement at the 672nd time increment (45 degree rotated in z axis). Clockwise: 6 mm, 16 mm and 36 mm radii.

Figure 13
Total displacement at the 600th time increment (front view). 6 mm at the back, 16 mm at the top and 36 mm radius at the bottom.

Figure 14
Total displacement at the 1258th time increment (isometric view). Clockwise: 6 mm, 6 mm detailed, 16 mm and 36 mm radii.

The same fringe pattern at the leading edge appears for the A0 wave, as can be seen in figure 15. However, as the wavelengths for the A0 modes are smaller, the fringes disappear for a smaller value of radius. It is possible to see that no elliptical mode is evident and the overall aspect of the wave is “cleaner” making this type of mode easier to use for LW-SHM, as it is already widely used throughout the literature (Qiu et al., 2013).

Figure 15
Total displacement at the 2972nd time increment (isometric view). Clockwise: 6 mm, 6 mm detailed, 16 mm and 36 mm radii.

3.1.1 Impact of findings for LW-SHM

According to the previous findings, for one to use the S0 mode for damage detection on (+45/-45)_S laminates, it is important to have in mind the underlying physics involved. In order to consider the curved part of the structure (e. g.: Leading edge on wings and wind turbine blades), the change in velocity (distortion of the wave front) for the curved part must be considered in the algorithms for finding location and size. Also, one needs to consider the change in velocity versus direction of travel along the laminate. Furthermore, it is needed to choose a wavelength that is smaller than the leading-edge perimeter, in order to avoid the fringes that act as damage signals, even for the healthy structure. Depending on the direction of detection versus source position, three distinct groups of waves shall appear in a signal over time, the first group (envelope) is the S0 wave, the second group is the elliptical interference wave from each lamina and the third group is the A0 wave, for relatively small frequencies (50 kHz or less for sure). The first group needs to be considered for the analysis of LW-SHM using the S0 mode and the third group needs to be considered for the A0 mode usage. If the distances between the source and the detector are small (less than 5 S0 wavelengths) all the wave modes mix, including the elliptical mode. This fact might make it difficult to use the S0 mode for LW-SHM.

The LW-SHM using the A0 mode is simpler, due to the following facts: the wave velocity is almost unchanged with direction, so that the algorithms for finding position and size just need to use a constant value of wave velocity; the wavelength is lower than the S0 mode, giving more options for choosing the wave frequency for adapting to the leading-edge radius; there are no elliptical secondary modes nor the interference region found for the S0 mode.

3.2 History points results

The applied load versus time that was described in the Methodology section is depicted in figure 16. The red trace line represents the peak of the envelope of load, that happens on 4e-5 seconds after the beginning of the vibration. This value of time was subtracted from the arrival time (peak) of the envelope of waves on each receptor point, in order to obtain the time of flight (TOF).

Figure 16
Load applied at the source point.

Figure 17 below shows the displacements over time for the point right before the bend (H1) and right after the bend (H2) for a 12 mm bend radius. The figure focus on the displacement y for the point before the bend (red solid line) and the displacement x for the point after the bend (blue solid line). To get the envelopes of the displacements, the positive peaks of displacements were connected by a spline curve. The envelope splines of H2 and H1 displacements are shown in the yellow and purple solid lines respectively. The peaks of the envelopes are depicted in the dashed lines and can be clearly seen in the zoomed details. It is possible to see that time of peak for the envelope is not intuitive to determine without a proper envelope algorithm, as it happens between the peaks of the displacement wave per see.

Figure 17
Displacements of interest for the intermediary points for 12 mm radius.

Figure 18 shows the same displacements of interest, but for the plain and curved receptors. Due to the bigger distance from the source, it is possible to see more clearly the elliptical mode interference around 4e-4 seconds.

Figure 19 shows the compiled data of S0 group velocity before the radius versus leading edge radius. This is the average velocity of the envelope of the wave in the portion between the source and history point H1. It is possible to see that the bigger the radius, the smaller the group velocity, asymptotically. The group velocity in that region is slightly bigger than the plain section velocity, which is the velocity between the source and the plain receptor (blue dashed line). The red line represents the plotting of a least squares regression yielding the following empirical equation:

$Cg=\frac{1}{{\left(r-8.35\right)}^{1.29}}+\mathrm{3.00.}$ (5)

Cg is the group velocity in km/s and “r” is the radius, in mm. The red dashed lines are the 95% confidence limits for the obtained equation from the least squares algorithm.

The A0 group velocity before the radius is practically unchanged and was not plotted.

Figure 18
Displacements of interest for the plain and curved receptors for 12 mm radius.

Figure 20 shows the S0 group velocity at the bend versus the radius. This is the average velocity of the envelope of the wave between H1 and H2 receptors. It is possible to observe that the group velocity increases with radius asymptotically. Also, the velocities are smaller than the benchmark (plain S0 velocity). The red line is the determined by the least square regression equation below:

$Cg=0.154\frac{r-15.8}{\sqrt{1+{\left(0.211\left(r-15.8\right)\right)}^2}}+\mathrm{2.01.}$ (6)

The red dashed lines are the 95% confidence limits for the obtained equation.

Figure 19
S0 group velocity versus radii before the bend.

Figure 20
S0 group velocity versus radii at the bend.

Figure 21 is analogous to figure 20, but for the A0 group velocity. This velocity decreases with radius, also asymptotically, and is always bigger than the benchmark. The least squares equation is shown below:

$Cg=\frac{1}{{\left(r-10.4\right)}^{1.46}}+\mathrm{0.833.}$ (7)

Figure 21
A0 group velocity versus radii at the bend.

Figure 22 shows the S0 group velocity after the leading edge, between H2 and the curved receptor. The equation follows:

$Cg=\frac{1}{{\left(r-7.38\right)}^{0.757}}+\mathrm{2.78.}$ (8)

As it happened for the velocity before the leading edge, the A0 group velocity is also practically unchanged after the leading edge and it was not plotted.

Figure 22
S0 group velocity versus radii after the bend.

The reasons why the S0 group velocity changes before and after the leading edge are unknown and need to be accessed in future works. They are probably related to the scatter interference originating from the leading edge which can impact the envelope of the wave before and after the bend. Regarding the behavior at the leading edge, we hypothesize that the bending stiffness decreases with an increasing radius, whereas the longitudinal stiffness increases with an increasing radius. Those values of stiffness are directly related to the wave velocity and, hence, that would explain the change in velocity acting differently for the two modes. But that also needs to be accessed on future works.

4 CONCLUSION

The present work showed the changes in the group velocity of the main Lamb wave modes at and near the bend (leading edge) for a glass fiber/epoxy bi-directional laminate. It was assessed a laminate with a 90 degrees bend using time domain finite element analysis for many values of leading edge radii. The findings can impact the field of Structural Health Monitoring using Lamb waves (LW-SHM), because the algorithms for LW-SHM rely on the properties of the wave, specially the wave speed, through the use of time of flight (TOF) algorithms. LW-SHM is a field of engineering that has a huge potential for decreasing operational and maintenance costs in the fields of aerospace (mainly for wing skin and airframe damage prediction) and renewable energy (wind turbine blade skin damage detection), approaching the framework of predictive maintenance.

Evaluating the field results from the finite element simulations, some quantitative and qualitative conclusions could be drawn. The wavelengths of the two main modes S0 (extension-contraction wave) and A0 (bending wave) were determined and their ratio compared with the literature. Another vibration mode was observed, with an intermediary velocity between the A0 and S0 mode. This mode was called Elliptical mode and its interference can affect the analysis of the S0 and A0 waves for small distances between source and receptor (less than 5 S0 wavelengths). It was also found that for values of leading edge perimeter closer or smaller than the wavelength, a scatter signal arises at the leading edge, which can be detrimental to the LW-SHM. Hence, a good practice is to choose the wave with a wavelength smaller than the leading edge perimeter to perform the SHM. In this regard, the A0 mode is easier, because it inherently has a smaller wavelength than the S0. Observing the wave fronts, it was possible to see evidence of refraction at the leading edge, with a bigger distortion of the wave fronts happening for smaller radii.

Through the analysis of the history points versus time, it was found that the S0 mode velocity changes near the leading edge, it decreases velocity with the radius asymptotically, whereas the A0 mode velocity is practically unchanged near the leading edge. For the curved portion (leading edge), the S0 wave mode increases the velocity with radius, asymptotically, whereas the A0 wave mode has the inverted behavior, it decreases velocity with radius, asymptotically.

Future SHM algorithms that are to be used on actual structures (not just plane plates) need to account for those behaviors. The group velocity variations with curvature can be an important source of error for damage positioning and sizing.

5 ACKNOWLEDGEMENTS

The authors would like to acknowledge CNPQ and CAPES for the funding of the present work, as well as the Mecbio laboratory at the Departamento de Engenharia de Estruturas (DEES) of UFMG that contributed with the Software licensing used for all the simulations inside this work. The second author acknowledge support from the São Paulo State Research Foundation (FAPESP, Grant no. 2017/20994-1).

References

• Abaqus (Dassault systemes simulia corp.). (2014a). Abaqus 6.14 analysis user guide. Providence, RI, USA.
• Abaqus (Dassault Systemes Simulia Corp.). (2014b). Abaqus 6.14 scripting reference guide. Providence, Ri, USA.
• Abaqus (Dassault Systemes Simulia Corp.). (2014c). Abaqus 6.14 analysis user manual. Providence, Ri, USA.
• Bettinotti O., Allix O., Perego U., Oancea V., Malherbe B., (2017). Simulation of delamination under impact using a global-local method in explicit dynamics. Finite Elements in Analysis and Design, Volume 125, Pages 1-13
• Courant, R., Friedrichs, K., & Lewy, H. (1967). On the partial difference equations of mathematical physics. IBM journal.
• Daniel, I. M., & Ishai, O. (1994). Engineering mechanics of composite materials. New York: Oxford university press.
• Ellabban, O., Abu-Rub, H., & Blaabjerg, F. (2014). Renewable Energy Resouces: Current Status, future prospects and their enabling technology. Renewable and Sustainable Energy Reviews, pp. 748-764.
• Fan, W., & Qiao, P. (2011). Vibration-based damage identification methods: a review and comparative study. Structural Health Monitoring, 10(1), pp. 83-111.
• He, Q., Gan, H., & Jiao, D. (2012). Explicit time-domain finite-element method stabilized for an arbitrarily large time step. IEEE Transactions on Antennas and Propagation, 60(11), pp. 5240-5250.
• IEA. (2013). Technology Roadmap - Wind Energy. International Energy Agency.
• Jiao, D., & Jin, J. (2001). A general approach for the stability analysis of time-domain finite element methods. Antennas and Propagation Society International Symposium. 4, pp. 506-509. IEEE.
• Lading, L., McGugan, M., Sendrup, P., RheinLander, J., & Rusborg, J. (2002). Fundamentals for remote structural health monitoring of wind turbine blades - a preproject. Annex B. Sensors and non-destructive testing. Forskningscenter Risoe. Risoe-R, 1341 (EN).
• Lamb, H. (1917). On waves in an elastic plate. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 93, pp. 114-128. The Royal Society.
• Lee T. Y., Chung K. J., Chang H., (2017). A new implicit dynamic finite element analysis procedure with damping included. Engineering Structures, Volume 147, Pages 530-544.
• Marston, P. L. (1989). Phase velocity of Lamb waves on a spherical shell: Approximate dependence on curvature from kinematics. The journal of the Acoustical Society of America, 85, pp. 2663 - 2665.
• Nazeer, N., Ratassepp, M., & Fan, Z. (2017). Damage detection in bent plates using shear horizontal guided waves. Ultrasonics, 75, pp. 155-163.
• Ono, K., & Gallego, A. (2012). Attenuation of Lamb Waves in CFRP Plates. Journal of Acoustic emission, 30, pp. 109-123.
• Qiu, L., Menglong, L., Xinlin, Q., & Shenfang, Y. (2013). A quantitative multidamage monitoring method for large-escale complex composite. Structural Health Monitoring, 3, pp. 183 - 196.
• Rahnavard R., Hassanipour A., Suleiman M., Mokhtari A., (2017). Evaluation on eccentrically braced frame with single and double shear panels. Journal of Building Engineering 10 13-25.
• Rahnavard R., Fard F. F. Z., Hosseini A., Suleiman M., (2018). Nonlinear analysis on progressive collapse of tall steel composite buildings. Case Studies in Construction Materials 8 359-379.
• Sanz-Serna, J., & Spijker, M. (1986). Regions of stability, equivalence theorems and the Courant-Friedrichs-Lewy condition. Numerische Mathematik, 49.2, pp. 319-329.
• Schubel, P. J., Crossley, R., Boateng, E., & Hutchinson, J. (2013). Review of structural health and cure monitoring techniques for large wind turbine blades. Renewable energy, 51, pp. 113-123.
• Sengupta, S., Datta, A., & Topdar, P. (2015). Structural damage localisation by acoustic emission technique: A state of the art review. Latin American Journal of Solids and Structures, 12, pp. 1565-1582.
• Towfighi, S., Kundu, T., & Ehsani, M. (2002). Elastic Wave Propagation in circumferential Direction in anisotropic cylindrical curved plates. Journal of Applied Mechanics - ASME, 69, pp. 283 - 291.
• Wang, C. H., Rose, J. T., & Chang, F. K. (2004). A synthetic time-reversal imaging method for structural health monitoring. Smart Materials and Structures, 13(2), pp. 415 - 423.
• Wang, L., & Yuan, F. G. (2007). Group velocity and characteristic wave curves of Lamb waves in composites: Modeling and experiments. Composites science and technology, 67(7-8), pp. 1370-1384.

• Available online: August 30, 2018

Publication Dates

• Publication in this collection
  2018

History

• Received
  18 July 2018
• Reviewed
  21 July 2018
• Accepted
  23 Aug 2018