Direct Calculation of the Group Velocity for Two-Dimensional Complex, Composite and Periodic Structures Using a Wave and Finite Element Scheme

Guided waves have immense potential for structural health monitoring applications in numerous industries including aerospace. It is necessary to evaluate guided wave dispersion characteristics, i.e., group velocity and phase velocity profiles, for using them effectively. For complex structures, the profiles can have highly irregular shapes. In this work, a direct method for calculating the group velocity profiles for complex, composite, and periodic structures using a wave and finite element scheme is presented. The group velocity calculation technique is easy to implement, highly computationally efficient, and works with the standard finite element formulation. The major contribution is summarised in the form of a comprehensive algorithm for calculating the group velocity profiles. The method is compared with the existing analytical and numerical methods for calculation of dispersion curves. Finally, an experimental study in a multilayered composite plate is conducted and the results are found to be in good agreement. The technique is suitable to be used in all guided wave application areas such as material characterisation, non-destructive testing, and structural health monitoring.


Introduction
Structural health monitoring (SHM) in the aerospace sector has the potential to reduce costs and increase reliability [1]. A number of nondestructive inspection techniques are already applied in the aerospace industry such as radiography [2], acoustic emission [3], eddy current [4], and ultrasonic inspection [5]. The use of advanced composite structures has also introduced an added layer of complexity where the failure modes are complex [6] and damage may occur between layers. Guided waves present a suitable SHM technique for thin structures [7][8][9]. They can travel large distances and interact with defects. Hence, they can be used for detection, localisation, and identification of damage [10,11].
Ultrasonic guided waves are guided by the traction free top and bottom surfaces of thin plate-like structures. They are generated due to the constructive and destructive interference and superposition of pressure waves and shear vertical waves undergoing multiple reflections from the boundaries [12]. Guided waves in plates can be classified into three broad categories based on their polarisation [13]: (a) symmetric (S) waves polarised perpendicular to the plane of the plate with the plate surfaces displacing out of phase, (b) anti-symmetric (A) waves polarised perpendicular to the plane of the plate with the plate surfaces displacing in phase, and (c) shear horizontal (SH) waves polarised in the plane of the plate. SH waves can also be symmetric or anti-symmetric with respect to the mid-plane of the plate. In theory, an infinite number of guided wave modes can propagate in a plate. Hence, the conventional notation used to represent them is with a subscript n, where n = 0, 1, 2, 3 . . . , i.e., S n and A n . For SH n , an even n denotes symmetric and an odd n denotes anti-symmetric waves. In multilayered structures, the elastic anisotropy results in quasi-Lamb waves and quasi-SH waves, which are mutually coupled, and the distinction is not as clear as in isotropic structures. It is still standard practice to adopt the nomenclature for guided waves in isotropic structures, which is followed here.
Guided wave dispersion curves provide important insight into the wave propagation in a structure. The group velocity curve shows the value of velocity versus the frequency. It is, however, important to calculate the change in the propagation direction of the group velocity versus the phase velocity for anisotropic materials. For isotropic plates with constant thickness, the group and phase velocities do not have angular dependence and the 2D profile is circular. This quickly becomes non-trivial in complex multilayered structures where the directions of wave propagation and energy propagation are different. The directional nature of composite structures creates another increase in complexity, where the 2D profile can take any shape depending on the composite layup. Nevertheless, identifying these characteristics is important for applications that use guided waves.
There are two main approaches for computing the dispersion properties: (a) analytical and (b) semi-analytical. The analytical methods consist of solving the governing differential equations by applying the appropriate boundary condition. They can further be divided into exact methods based on 3D elasticity theory and semi-exact methods based on approximate plate theories. The 3D elasticity method was used to calculate the phase velocity of composite lamina and laminate [14,15]. The exact methods are computationally intensive due to the transcendental equations used for multilayered composites. Therefore, approximate plate theories are employed to ensure the solutions are practical [16,17]. Recently, Biot's general energy expression was applied to derive the group velocity vector of elastodynamic guided waves [18] in arbitrary layered structures. The semi-analytical methods use finite element (FE) approaches to model part of the structure and impose wave solutions in the propagation directions. The most well-known is the semi-analytical finite element (SAFE) method, where FE is used to model the cross-section, and a complex exponential function is used to describe the displacement field in the direction of wave propagation [19]. However, the SAFE method cannot be used to model truly periodic structures [20] that consist of identical substructures joined to form a continuous structure. The scaled boundary finite element method (SBFEM) was also applied in a similar semianalytical framework for guided wave simulations [21] but only in isotropic structures. The wave and finite element (WFE) method has been developed over the last decade [22]. It combines FE and periodic structure theory to model wave propagation in periodic structures. It has been used in a number of research areas like structural identification [23], damage detection [24], multi-scale wave propagation [20], transient wave propagation [25], and wave steering in composites [26].
The main issue associated with the existing group velocity calculation methods is that they require custom formulations, which are difficult to implement, and none of them are able to handle truly periodic structures. Moreover, the group velocity from WFE was approximated by the slope of the wavenumber-frequency dispersion curve using finite differentiation. The principal novelty of the work exhibited in this manuscript is the calculation of the group velocity curves for two dimensional complex periodic and/or composite structures using a wave and finite element scheme without involving any finite differentiation. The approach is therefore free from computational errors, requirement for convergence checks, and additional computational cost implied by finite differentiation schemes. The use of periodic structures such as textile composites has become increasingly common in automotive and aerospace industries [27]. In order to employ guided-wavebased damage localization and identification, it is necessary to calculate the dispersion properties of such structures. The WFE scheme based calculation of dispersion curves has the following advantages: • It is capable of handling both periodic and layered composite structures.
• It is relatively easier to implement and does not require custom FE formulation. • It is computationally efficient and has minimal memory requirements.
The novelty of this work is the form of a comprehensive algorithm to calculate the dispersion properties of arbitrary plate-like structures. The paper is organised as follows: A brief introduction to the WFE method is presented in Section 2, followed by the energy and power flow considerations necessary for dispersion curve calculation. At the end of Section 2, the main algorithm is presented. The symbols used in this section are listed in Table 1. The method is compared with existing dispersion properties calculation techniques in Section 3, highlighting the main advantages of this approach. The comparison with experimental results is described in Section 4. This is followed by our conclusions in Section 6. Group velocity and phase velocity θ Propagation direction

Dispersion Curve Calculation
The methodology for obtaining the dispersion curves of 2D periodic structures is presented in this section. The WFE method is used to obtain wave propagation characteristics, wavenumbers, and wave modes of the periodic waveguide. This is followed by the necessary energy and power flow considerations for calculating the group and phase velocities of guided wave modes. Finally, the algorithm for calculating the 2D group velocity profile is provided.

Wave and Finite Element Scheme
The WFE method has been developed over the last three decades and has been applied to a wide variety of structural calculations. The reader is referred to classical literature regarding this methodology [28][29][30] in order to obtain more information. With regard to twodimensional periodic media as shown in Figure 1a, and a close-up example of such a structure shown in Figure 1b, the first step is to model a single periodic section using standard finite element software. There is no limitation with regard to the complexity of the cross-section to be modelled, making the method suitable for modelling layered media. The finite element model of the periodic section is used to obtain the stiffness (K), mass (M), and damping (C) matrices to set up the dynamic equilibrium as follows: where  (1) is commonly known as the dynamic stiffness matrix D. The dimension of Equation (1) can be reduced by applying the periodicity condition in the y-direction. For a free wave propagation scenario, where k y refers to the component of wavenumber in the y-direction, a transformation matrix T can be specified that reduces the nodal degrees of freedom vector as follows: The entries of matrix T are determined from the wave propagation constant λ y , which, according to Bloch's theorem, is given as λ y = e −ik y l y , where l y is the height of the periodic section. The complete matrix depends on the ordering of the q vector; an example can be found in [28]. This projects all degrees of freedom on the nodes lying on the x-axis of the waveguide by imposing the periodic boundary conditions. A similar dimensional reduction is performed to obtain the reduced force vector f red . Since there are no forces on the internal nodes, the internal nodal force vector f b = 0. This helps in simplifying Equation (1) as follows [31]: with the dynamic stiffness matrix being partitioned with respect to the boundary nodes of the periodic section. Similar to λ y , the free wave propagation in a waveguide of length l x has a characteristic propagation constant λ x = e −ik x l x , which relates the displacements of the nodal degrees of freedom on the left-and right-hand sides by q rb = λ x q lb and f rb = −λ x f lb . An eigenvalue problem for λ x is formulated by substituting the above into Equation (3), that is, where where T is known as the transfer matrix [22]. The propagation constants λ x are eigenvalues of T. They exist in pairs of [λ + x andλ − x ] for the positive and negative travelling waves, respectively. A number of formulations may be employed to mitigate this effect as described in [32]. The wavenumbers k x are directly derived from the computed propagation constants. The corresponding eigenvector φ φ φ x is also extracted and represents the mode shape of the considered wave type. The mode shape vector is internally partitioned to represent displacements and forces such that

Energy and Power Considerations
The total energy flow and time-averaged power associated with wave propagation in a structure provide important insight into the properties of wave modes such as group velocity and the fundamental principle of energy conservation. The time-averaged total energy of a particular wave mode is the sum of kinetic (E k ) and potential (E p ) energies [31]. The relations in terms of the wave basis are shown below: Here, (.) H denotes a Hermitian conjugate and φ φ φ q,x is the displacement mode shape for a particular wave mode. For an element length ∆l in the propagation direction, the total time-averaged energy per unit length is given as: The time-averaged power (P) associated with a propagating wave is obtained as follows [30]: The group velocity c g for a guided wave is defined as the speed at which energy propagates in a particular wave. It has a number of implications for setting up a guided wave simulation, such as the maximum element size in the wave propagation direction and the maximum time step increment allowed. Normally, there should be at least 10 elements per wavelength and the wave should not propagate more than one element distance in a single time step [25]. The knowledge of c g as a function of frequency for a wave mode helps with fixing these parameters a priori with sufficient degree of confidence. It is obtained from the ratio of power and energy as follows: where c g,x denotes the x-component of group velocity of a particular wave mode. Notably, group velocity is only defined for a propagating wave as evanescent modes do not transfer any energy.

Velocity Dispersion Curves
The main contribution of this paper is Algorithm 1, which is presented in this section, for obtaining group velocity dispersion curves for periodic structures without finite differentiation using a WFE scheme. The algorithm starts following the basic steps of WFE method, where a periodic section of the structure is modelled in FE software and the K, M, and C matrices are extracted. They are used to calculate the transfer matrix T and obtain its eigenvalues and eigenvectors.
Algorithm 1: Group and phase velocity calculation 1 Model periodic section in FE and extract K, M, and C 2 Specify angular frequency ω 3 Initialise arrays k x , k y , θ θ θ, andc g,x , c g,y 4 while θ < θ max do 5 Calculate the transformation matrix T and transfer matrix T if k y = 0 and (θ n − θ n−1 ) > ∆θ then 10 Reduce step size ∆k y 11 Go to line 5 12 else 13 Append k x , k y , θ to k x , k y , θ θ θ respectively 14 Update k y ← k y + ∆k y end 15 Calculate c g,x using Equations (6)-(9) and append to c g,x end 16 Append 0 to k x 17 foreach k x in k x do 18 Perform steps 5 to 14, but this time with φ φ φ y as unknowns 19 Calculate c g,y using Equations (6)-(9) and append to c g,y end 20 Calculate group velocity c g (θ) and phase velocity c p (θ) using Equation (10) In the while loop starting on line 4, the reduction is performed as shown in Equation (2). The variable k y is initially assigned zero for a point in the k x − k y plane lying on the positive x−axis (θ = 0°). The value of k y is updated in subsequent iterations of the while loop by a suitable step size ∆k y , such that the polar angle between two consecutive points is less than a desired angular resolution ∆θ. At each step, θ is calculated for a particular wave mode by the wavenumbers of that mode using the relation θ = tan −1 (k y /k x ). Note that sorting and identifying a single mode may be difficult especially if a number of propagating modes exist at the excitation frequency. Normally, a modal assurance criterion (MAC) [33] is used to identify modes that can be used here, or a more advanced technique such as the partial wave method presented in [34]. Since a single mode is used for calculating θ, the algorithm can only be used to extract one mode at a time. The same algorithm can easily be extended to filter multiple modes but that is not presented here. The mode shape vector φ φ φ x is then used to calculate the x−component of group velocity c g,x (θ). The loop runs until θ ≥ θ max , where θ max is slightly below 90°. This essentially means that we are looking for part of the dispersion curve that lies in the first quadrant (i.e., between 0°and 90°). The reason for stopping before reaching 90°is the asymptotic behaviour of the tangent function at that angle, which makes it numerically difficult to compute. The exact value at 90°is computed later in the algorithm to bypass this issue.
The part of the algorithm discussed so far is used for calculating the x−component of the group velocity. To compute the actual group velocity c g (θ) of the mode, c g,y (θ) must also be calculated. We find that the ground work for the y−component has already been completed. The first zero is appended to the array k x to compute k y corresponding to 90°. It is important to understand that k x is zero at this angle. Then, a for loop is executed, which iterates through each element of k x . This time, the reduction is performed such that q red = {q T lb q T lt q T l } T . The entries of the transformation matrix are modified based on the new q red and now depend on the propagation constant λ x , as all the degrees of freedom are being projected on the nodes lying on the y−axis of the waveguide. The wave mode shape φ φ φ y is the main unknown quantity that needs to be computed. After that, the calculation of c g,y (θ) follows straightforwardly from Equations (6)- (9). Thereafter, the group and phase velocity profiles in the first quadrant are given by: where k(θ) = (k x (θ)) 2 + (k y (θ)) 2 .
The next step is to obtain the velocity profile in the quadrant of 270°to 360°. The same algorithm is used to calculate the group velocity profile using the prior knowledge that k y ≤ 0 in this quadrant. Hence, the step size ∆k y becomes negative. The rest of the steps for calculating the velocity profile remain the same as for the first quadrant. This essentially leads to the complete 360°profile by realising that the remaining two quadrants are mirror images (i.e., the profile in quadrant 180°to 270°is a mirror image of the profile in quadrant 0°to 90°and the profile in quadrant 90°to 180°is a mirror image of the profile in quadrant 270°to 360°).

Comparison with the Literature
The characteristic dispersion curves are important quantities needed for any guided wave application. Therefore, numerous methods have been developed over the years to calculate them. They range from purely analytical methods to purely numerical methods. In this section, the WFE-based method is compared with the most popular existing techniques for the calculation of dispersion curves, and its main advantages are highlighted. All the comparisons were performed against results that have been already published. Notably, none of these methods can be applied to truly periodic structures. The periodic sections for WFE were modelled with standard eight-node hexahedral elements with three degrees of freedom per node. They are shown in Figure 2. The layup was modelled by assigning directional material properties to elements through the thickness.

3D Elasticity Theory
In this approach [13], the exact dispersion relations are obtained using 3D elasticity theory. The transcendental equations are numerically solved for an infinite number of possible modes. The formulation consists of a robust method to separate symmetric and anti-symmetric modes in symmetric laminates. The wave curves are then obtained by performing the finite difference on the exact solutions of two adjacent slowness surfaces. The results for a 1 mm thick laminate [+45/ − 45/0/90] s are compared with the WFE scheme in Figure 3, where excellent agreement can be observed. The composite material used in this study was AS4/3502 graphite/epoxy with material properties described in [13]. The material properties were assigned to each layer of the periodic section shown in Figure 2a. The 3D elasticity model is limited to symmetric laminates, whereas the WFE-based approach can handle arbitrary layered structures. Moreover, the solution of transcendental equations is computationally taxing, whereas the developed approach is highly efficient, requiring minimal computational resources with the solution obtained in a fraction of a second.

Semi-Analytical Finite Element (SAFE) Method
The SAFE method is similar to the WFE scheme, where the finite element method is used for discretising the cross-section of the waveguide [35]. The SAFE method uses a harmonic exponential term to describe the wave behaviour in the propagation direction. The FE discretisation occurs at the cross-section of the waveguide. For plane wave propagation in a plate, a 1D discretisation across the thickness of the plate is sufficient. Combining the harmonic assumption and FE discretisation, the particle displacements are expressed in terms of shape functions. Custom stiffness matrices are assembled using this representation, which are used to set up an eigenvalue problem for the frequencydependent wavenumber calculation [36]. WFE has a major advantage as it is set up using standard stiffness and mass matrices and no custom formulation is required. Moreover, SAFE cannot be used to model the dispersion characteristics of periodic structures such as the one shown in Figure 1b. The SAFE and WFE models are here compared for the 1 mm thick laminate [+45/ − 45/0/90] s presented in [35]. The composite laminate in this case is composed of unidirectional T800/924 graphite/epoxy laminae with the material properties described in [35]. The WFE model of the periodic section is shown in Figure 2a, which is again modelled with standard eight-node hexahedral elements. The comparison results are shown in Figure 4, where a similar trend is observed in both approaches with a slight difference in the values.

Standard Finite Element
The standard finite element (FE) method for extracting the dispersion curves consists of performing a series of modal analyses of a representative part of the structure [37,38]. The dispersion relations are obtained by simply changing the length of the representative portion in order to allow the wavelength to change. The mode shapes and corresponding natural frequencies can be obtained by solving eigenvalue problems. Then, the dispersion curves are obtained from the eigensolutions. This method can only be used to extract propagating modes, whereas the WFE methodology is more general as it can extract both the propagating and evanescent modes. The two techniques were compared for a 2.16 mm thick laminate [0/90/0] s . The FE model is shown in Figure 5, where a coarse mesh is shown for visual representation. The laminate is modelled with 2D plane strain elements with the plane strain in the z-direction. The element length is 0.09 mm in the propagation direction and four elements per layer through the thickness. The length of the representative part needs to be changed to obtain the result for different wavelengths. The comparison of dispersion curves is shown in Figure 6, where a good agreement can be observed. The WFE periodic section is shown in Figure 2b, which is modelled using standard eight-node hexahedral elements and the composite material properties for each layer is taken from [37]. The WFE periodic section has orders of magnitude fewer degrees of freedom compared with the standard FE approach. Moreover, the standard FE provides discrete output for each wavelength by changing the length of the representative part, whereas a single periodic section is used in WFE, which can be used to extract dispersion curves for a range of frequencies.

Comparison with Experimental Results
In this section, the group velocity computed from WFE is compared with experimentally determined group velocity values for a 2 mm thick CFRP laminate [+45/0/ − 45/90/ + 45/0/ − 45/90] s . The composite laminate is composed of unidirectional T800S carbon fibre/epoxy laminae with the material properties described in Table 2. The dimension of the plate is 200 by 200 mm. The excitation signal was generated from the TG5011 waveform generator and amplified by the Ciprian amplifier. The pulse generator settings were set to a single cycle sine burst centred at 100 kHz for the first test ( Figure 7) and 200 kHz for the second test. The experimental setup is shown in Figure 8. The output of the amplifier was split using a Tjunction. One output was directed toward the excitation piezoelectric transducer, while the other was attached to the Picoscope signal acquisition system to act as a trigger. The second transducer was also connected to the Picoscope, which received the signal with a sampling frequency of 10 MS/s. All the data were stored in a laptop with the Picoscope 6 program installed. The waves were excited from the centre of the plate and received on a circular boundary 50 mm from the central point of excitation, as shown in Figure 8a. The circular boundary had sensing points in 15°intervals. The group velocity for the fastest wave mode at this frequency (i.e., S 0 ) was calculated by the time of flight (ToF) method. The ToF of the S 0 wave was determined by measuring the difference between the time the signal was transmitted and received by the transducer. The ToF and and known distance were used to calculate the velocity of the wave at that given point. The tests were repeated 10 times and the average values were computed. The pointwise difference between simulation and experimental data was calculated and found to be less than 5% for all observation points. The results are shown in Figure 9. Notably, the experimental approach used in this research is the simplest method to calculate the group velocity and is limited to the fundamental S 0 mode. More advanced techniques are available that can robustly capture higher-order modes such as the one presented in [13]. The WFE model is also an idealisation of the physical specimen, which might have inherent manufacturing flaws, which is common for composite structures but is not represented in the FE model. Despite these challenges, the experimental and numerical results were found to be in good agreement.

Discussion
The proposed WFE-based group velocity calculation technique was described and compared with existing techniques. The model is computationally efficient and can be used for plate-like structures of arbitrary complexity. However, the model is limited by the periodicity assumption in the plane of the plate and hence cannot model variation in cross-section geometry. Moreover, the material properties must also remain consistent in the plane. The WFE scheme also depends on the use of a FE modeller and cannot exist as a standalone platform [40]. The computational efficiency of the model is also reduced in the case of a highly complex single periodic section, which might warrant a further model reduction strategy. In case of a large number of degrees of freedom in a single periodic section, further postprocessing is required to identify and separate propagating and evanescent modes.

Conclusions and Outlook
This work presented a new approach for the direct calculation of group velocity curves using the WFE scheme. The WFE scheme was used to exploit the periodicity of the structure to extract the group and phase velocity curves from a representative periodic section. A comprehensive algorithm for calculating the velocity profiles was also constructed. Since the technique is based on FE method, it contains all the uncertainty associated with numerical modelling. This includes the effect of mesh size, element formulation, and model description. In this regard, all the standard best practices of FE modelling should be followed. Moreover, the accuracy of the results will only be as accurate as the model. In case of complex composite structures that have inherent manufacturing flaws, the results will always have uncertainty due to the model definition of the periodic section. Nevertheless, all modelling approaches face similar challenges. The method was compared with existing dispersion curve calculation techniques by comparing the results for complex multi-layered structures. A number of advantages were highlighted with respect to the existing techniques, which are summarised here: • The technique is easy to implement using the standard FE method without the need for custom formulations. • It is capable of calculating dispersion properties for arbitrary layered plate-like structures. • The periodic section is modelled with very few elements, which makes it highly efficient and computationally cheap.
• The theory of the WFE method enables the technique to be used for truly periodic structures, such as textile composites, which cannot be modelled with any existing technique.
Considering the strengths of this approach, it can be used in all guided wave application areas such as material characterisation, non-destructive testing, SHM, etc. It is especially useful in the NDE of aerospace structures, which are fabricated with layered composites and textile composites. It provides a suitable alternative to existing techniques while having the potential to hande more complex structures such as textile composites.