Lamb waves are a growing method for the Non -Destructive Testing and Evaluation (NDT &E) and Structural Health Monitoring (SHM) of aerospace vehicles. These guided waves can propagate over large distances and have a strong tendency to interact with damage. Whilst several me thods exist for the modelling of Lamb wave propagation, this paper is the first to introduce a first principles numerical model that can efficiently and accurately predict the behaviour o f Lamb waves. The numerical model is easier to understand and impleme nt compared with analytical solutions and significantly faster than discretised numerical methods. The numerical model is presented in detail for an isotropic and homogenous plate , along with validation against the industry accepted, WaveForm Revealer 3 (W FR3) software. The results show a mean correlation across all assessed parameters of 90.4% and 96.6% for the symmetric and antisymmetric modes, respectively. Further discussion is provided on future developments to the model, including on the topic of high temperature effects, anisotropic materials, and edge reflections. 2 Introduction Non -Destructive Testing (NDT) has been around for over 100 year s, with o ne of the very first applications being in the railroad industry [1]. The 1970s significantly changed the field of NDT with the emergence of the discipline of fracture mechanics [1]. This enabled the prediction of whether a crack of a given size would fail under a given load if the relevant material property or fracture toughness was known. Laws to predict the rate of growth of cracks under cyclic loading were also developed. With these laws, it became possible to use structures with known defects, if the sizes of the defects were known [2]. That is, design incorporated a new approach, and they became “damage tolerant designs” [3]. Components with known defects could be used as long the defects did not grow past a critical limit, which could result in a catastrophic failure. This presented a new challenge for NDT, ‘simple’ detection of flaws was no longer enough , quantitative information was needed about the size of flaws to enable pr edictions of the component’s residual life. That is, Non-Destructive Evaluation (NDE) [1, 4]. Structural Health Monitoring (SHM) is an emerging area which brings together NDE with smart technologies and materials [5]. The drive for SHM in the aerospace industry lies with the need to improve safety , for both aerospace manufacturers and operators. Maintenance activities can account for as much as 20% of the direct operating costs , and huma n error in maintenance contributes to 15% of all aircraft accidents [5]. Beyond that, reliable SHM would allow for less conservative design philosophies to be used, leading to significant weight savings on future aircraft whilst improving the safety of passengers. The rapid advances that have been made in recent decades in sensor technology, smart materials , and innovative structural concepts are leading to the practical use of SHM. One such area is the use of ultrasonic guided w aves [6, 7], a direct evolution of traditional ultrasoni c NDT. Simulation of ultrasonic guided wave based NDT and NDE for SHM is onerous [6]. While the use of Finite Element Analysis (FEA) has become more common [8], there are limitations and issues with such a method . Specifically, FEA can be computationally intensive [9]. Reducing computational effort is essential for airframe digital twins [10, 11], especially if an intelligent digital twin is used onboard, in flight [12]. In this work, a simple numerical model for guided wave propagation in 2D plate structures is presented. These structures made specifically of alum inium alloys are representative of those used to skin aircraft [13, 14]. The Numerical Model The original analysis of guided waves in thin plates was published by Lamb in 1917 [15] and derives itself from an interrogation of the elastodynamic Navier -Cauchy equations for isotropic materials [16]. 𝜇𝑣𝑖,𝑗𝑗+(𝜆+𝜇)𝑣𝑗,𝑗𝑖+𝜌𝑓𝑖=𝜌𝑣𝑖̈ (1) The geom etry is simple, the thin plate is a solid medium which is bounded on the top and bottom by parallel planes which are a distance 2ℎ apart , as depicted in Fig. 1 [17]. In the horizontal distance within the plate, the dimension s are unbounded, and it is also assumed that there is no medium outside the upper and lower bounds of the solid plate. In a bulk solid medium, mechanical waves can propagate as longitudinal (pressure) waves, or as lateral (shear) waves. In a plate, the clo se nature of bounding planes results in a mixing of these types 3 of waves, and interestingly results in the production of two unique propagation modes. These are the symmetric wave which is an even function of the displacement of the medium in 𝑧, where 𝑢(𝑧)=𝑢(−𝑧), and the antisymmetric wave which is an odd function of the displacement of the medium i n 𝑧, where −𝑢(𝑧)=𝑢(−𝑧) [17]. The longitudinal component of the displacement 𝑢 and a transverse component 𝑤 have the form, 𝑢=𝑖(𝑘𝐴cos(𝑝𝑧)+𝑞𝐵cos(𝑞𝑧))𝑒𝑖𝑘(𝑥−𝑐𝑡) (1.1) 𝑤=(−𝑝𝐴sin(𝑝𝑧)−𝑘𝐵sin(𝑞𝑧))𝑒𝑖𝑘(𝑥−𝑐𝑡)(1.2) for symmetric modes, and 𝑢=𝑖(𝑘𝐶sin(𝑝𝑧)−𝑞𝐷sin(𝑞𝑧))𝑒𝑖𝑘(𝑥−𝑐𝑡)(2.1) 𝑤=(𝑝𝐶cos(𝑝𝑧)−𝑘𝐷cos(𝑞𝑧))𝑒𝑖𝑘(𝑥−𝑐𝑡)(2.2) for antisymmetric modes. While A, B, C, and D are arbitrary constants, the coefficients 𝑝 and 𝑞 represent the dispersive wavenu mbers, and are given by, 𝑝2=𝜔2 𝑐𝐿2−𝑘2(3.1) 𝑞2=𝜔2 𝑐𝑇2−𝑘2(3.2) where 𝜔 is the circular frequency of the wave, 𝑘 is the wavenumber determined by the phase velocity of the wave at 𝜔, 𝑐𝐿 is the longitudinal wave velocity of the bulk material, and 𝑐𝑇 is the transverse wave velocity of the bulk material. 𝑐𝐿=√𝜆+2𝜇 𝜌(4.1) 𝑐𝑇=√𝜇 𝜌(4.2) The variables 𝜆 and 𝜇 represent t he Lamé constants as given below [18]. 𝜆=𝐸𝜈 (1+𝜈)(1−2𝜈)(5.1) 𝜇=𝐸 2(1+𝜈)(5.2) The boundary conditions are that 𝑡𝑥𝑧=𝑠𝑧𝑧=0 on 𝑧=±ℎ [19]. The result of imposing these boundary conditions is that (1 .1) to ( 2.2) only have nontrivial solution s if, tan(𝑞ℎ) tan(𝑝ℎ)=−4𝑘2𝑝𝑞 (𝑞2−𝑘2)2(6.1) tan(𝑝ℎ) tan(𝑞ℎ)=−4𝑘2𝑝𝑞 (𝑞2−𝑘2)2(6.2) for the symmetric and antisymmetric modes, respectively [19]. 4 Figure 1. The general setup and co -ordinate system for a guided Lamb wave. The important result is that Lamb waves are dispersive, given the velocity depends on the product of ℎ and 𝜔, or the thickness frequency product. That i s, any arbitrary wave, which based on Fourier’s theorem can be described as the sum of many sinusoidal waves, will experience dispersion, as well as generating two wave modes that travel independently of each other. The group velocity given as 𝑐𝑔=𝑑𝜔 𝑑𝑘 can be related to the phase velocity through the equation (7). 𝑐𝑔=𝑐𝑝2[𝑐𝑝−(𝑓ℎ)𝑑𝑐𝑝 𝑑(𝑓ℎ)]−1 (7) Where 𝑓ℎ represents the frequency -thickness product of the case. The frequency -thickness product is often used as the independent variable as it allows for calculation of the phase velocity independent of the thickness of the plate. Model ling the propagation of an ultrasonic guided wave in a plate structure can be achieved using the previous stated fact. Any signal is a sum of sinusoidal waves. Each sinusoida l wave will have a velocity dictated by the dispersive relationship, which is a function of the material properties and the plate thickness. Consider a typical Gaussian modulated sinusoidal signal , 𝑓(𝑡)=(𝐴cos(𝜔𝑡))𝑒𝑡2 2𝑐2⁄(8) The Fourier transform of this signal will be a Gaussian distribution of frequencies about the central frequency, 𝜔 [20]. Rather than defining analytically what the Fourier transform of ( 8) is, a discrete numerical approximation of the function can be defined with a set sampling rate. The fast Fouri er transform (FFT) of this discrete signal will then provide the phase and amplitude information for all the sinusoids that sum together to form this signal. In the frequency domain, the spacing of these constituent signals ( 𝑑𝑓) is derived from the sampl ing rate (1/𝑑𝑡) and the total number of samples ( 𝑛), such that, Δ𝑓=1 𝑛Δ𝑡(9) Knowing each frequency, the dispersion for each can be determined. There are two options, root finding algorithms can be used for ( 6.1) and ( 6.2) to determine the valid wave sp eed at each frequency, or a lookup table for the given material with a given thickness can be utilised. 5 The second provides a more computationally efficient means of determining the dispersion but is less general purpose . For each signal, a time component can be determine based on the distance between the transmitter and receiver, and the wave speed at that frequency. Importantly, the relative phase information from the FFT is preserved and used as the initial phase offset of the sinusoid before adding the traveling wave component. For the received signal, there is no inverse problem that needs to be solved. That is, each of the signals at the receiver location do not have their amplitude and phase information taken to for the FFT of the received signal, whi ch is the inverted (iFFT) to give the received time domain signal. Rather, the direct sum of the array of sinusoids is the time domain signal. Effectively, the computationally intensive solution to the inverse problem is solved with the novel application o f first principles physics. The overall procedure is illustrated in Fig. 2 . Figure 2. The process, 1) transmitted signal (Tx), 2) taking the fast Fourier transform (FFT(Tx)), 3) giving an array of amplitudes and phases applied to a sum of sinusoids, 4) each with a different frequency for which the wave speed can be looked up, 5) which can then be used to travel the signal to the rece iver at a specified distance away, and 6) finally, all of the travelled signals can be summed up to give the received signal (Rx). Other Modelling Approaches Several numerical methods already exist for the modelling of Lamb waves in plate -like materials su ch as the Finite Element Method (FEM), Spectral Element Method (SEM), Semi- Analytical Finite Element method (SAFE). These numerical modelling techniques are often 6 time intensive in both setup and operation, requiring extensive resources. However, the use of such methods permits analysis of complex structures such as composite materials, stepped geometr ies, assemblies, and damage interaction. Finite Element Analysis (FEA) The finite element approach for the computation of Lamb wave behaviour is the most vers atile and resource intensive. Using FEA, the entire structural model is spatially discretised, and a series of partial differential equations are solved using adequate boundary conditions for the time response of the model. FEA is a powerful tool for the e xploration of Lamb wave responses to various geometries [21], materials [22], damages [23-25] and sensor -material interactions [26]. Analysis of Lamb waves using FEA is also one of the most accessible methods due to the availability of commercial software. However, the high resource demand has led to the development of various other numerical methods. Spectral Element Analysis (SEA) An alternative method to FEA i s that of SEA developed principally for work with fluid flows [27]. SEA has a high degree of similarity to FEA, but offers large computational savings due to the use of Lagrange interpolation functions [28]. The SEA method has been used in similar applications to FEA such as material modelling [22], damage [29] and sensor -material interactions [30]. Ultimately, SEA offers significant benefits over FEA but is otherwise less accessible du e to the lack of commercial software available. Semi -Analytical Finite Element Analysis (SAFE) The semi -analytical method for the simulation of Lamb wave propagation in elastic wave guides has been studied intensively. The SAFE method offers computational savings over other, pure, finite element methods [31]. The SAFE method was first formulated by Aalami [32] in 1973. Using SAFE, the waveguide is discretised across the cross -section whilst an exponential function is used in the propagation direction [31]. The SAFE method has been used for the efficient calculation of the dispersive properties of a Lamb wave in anisotropic laminates and has been shown to have exact solutions for homogenous materials [33]. The SAFE method is an effective tool for calculating the propagation of Lamb waves in irregular and non -planar geometries [34-36]. Recent work with SAFE has involved the modelling of more complex materials and interactions such as viscoelastic effects in composite mediums [37], boundary reflection interactio ns [32], reflection and transmission behaviour for damage and discontinuities [38], pre -stresse d materials [39] and leaky Lamb wave behaviour [40]. The SAFE method requires integration of analytical equations with typically finite element methods and is generally developed for individual cases . WaveForm Revealer 3 As an alternative, WaveForm Revealer 3 (WFR3) developed by the University of South Carolina is an analytical method for the modelling of Lamb wave propagation by the convolution of a user input function with a structural transfer function in the frequency domain [41]. The structural transfer function is derived by Giurgiutiu [18, 41] and given below. 𝐺(𝑢,𝜔)=𝑆(𝜔)𝑒𝑖𝑘𝑆𝑢+𝐴(𝜔)𝑒𝑖𝑘𝐴𝑢(10) 7 Where the superscript S and A and functions 𝑆(𝜔) and 𝐴(𝜔) correspond to the symmetric and anti-symmetric components, respectively. These are given below. 𝑆(𝜔)=−𝑖𝑎𝜏0 𝜇sin(𝑘𝑆𝑎)𝑁𝑆(𝑘𝑆) 𝐷𝑆′(𝑘𝑆)(11.1) 𝐴(𝜔)=−𝑖𝑎𝜏0 𝜇sin(𝑘𝐴𝑎)𝑁𝐴(𝑘𝐴) 𝐷𝐴′(𝑘𝐴)(11.2) The analytical model for WFR3 is programmed in MATLAB and ope rated through a graphical user interface (GUI). WFR3 is also capable of modelling Lamb wave damage interactions such as transmission, reflection, and conversion behaviour. WFR3 has been validated against FEM and experimental testing such that it is now fav oured by industry and academia [41]. Moreover, the analytical method offered significant computational time savings over typical numerical techniques. As a result of the low -resource requirement of WFR3, it has been used in later sections to verify the proposed first principles numerical model. Validation The first step in deploying the numerical model is to determine the d ispersion relationships for the given material. Firstly, a range of phase velocities must be selected in which the solution to the Rayleigh -Lamb equation is predicted to exist. These selected phase velocities are used to calculate the predicted wave number s for the inputted frequency range which are then used to determine the 𝑝 and 𝑞 constants as given in equations (3.1) and (3.2). To begin solving the Rayleigh -Lamb equation, it is better to rearrange it into a more useful form suitable for solving with a numerical method. The rearranged forms of equation (6.1) and (6.2) are given below in subscript notation. tan(𝑞𝑛ℎ) tan(𝑝𝑛ℎ)+4𝑘𝑛2𝑞𝑛ℎ (𝑘𝑛2−𝑞𝑛2) 2=𝑓𝑛𝑆(12.1) tan(𝑞𝑛ℎ) tan(𝑝𝑛ℎ)+(𝑘𝑛2−𝑞𝑛2) 2 4𝑘𝑛2𝑞𝑛ℎ=𝑓𝑛𝐴(12.2) The exact solution for the wavenumber and hence, phase velocity occurs when 𝑓𝑛=0, representing a root of the rearranged Rayleigh -Lamb equation. A variety of numerical solvers can be applied at this point such as the bisection, secant, or Newton -Raphson m ethod. For most solvers, an interval must be determined in which the root is known to lay. Hence, the range [n, n+1] in which a root is found is determined for the case that 𝑓𝑛+1⋅𝑓𝑛<0. Discovery of this range restricts the search field, rendering the a bility to apply the desirable solver. For the case of higher order wave modes, it must be acknowledged that more than one root may exist within the search fiel d, requiring more rigorous methods. The search resolution for the root can be coarse as the Rayle igh-Lamb equation has high smoothness and hence, a cubic spline may be fitted to the roots to improve the efficiency of the process. Once the phase velocities have been found for the given frequency range, the group velocities may be determined by applying a suitable finite difference approximation to equation (7). 8 Figure 3. Flow chart depicting the first principles numerical method. Following determination of the dispersive rela tionship, the numerical model may be applied as shown in Fig. 2 and Fig. 3 . To validate the proposed method and evaluate it against common practices, a parametric study was performed for a range o f variables. The numerical model was compared with the industry accepted and peer -reviewed software WFR3, developed by the University of South Carolina. The study was conducted using a Gaussian windowed sinusoid centred at a peak frequency of 200 kHz with a fractional bandwidth of 30% as the input signal , shown in Fig. 4. Figure 4. The time domain input si gnal (left) and its frequency spectrum (right). In performing the parametric study, the symmetric and antisymmetric wave modes were compared separately. The choice of sensor has a dramatic impact upon the amplitude ratio between the two modes, most notably ; fibre Bragg gratings are much more sensitive to the symmetric mode as they measure in -plane displacement whilst piezoelectric transducers are generally more sensitive to the antisymmetric mode due to their measurement of out -of-plane displacement [42]. Hence, the amplitude relationship between the two modes would great ly differ dependent upon the sensor choice. Moreover, it has been shown that piezoelectric transducers can be carefully selected such that they may be tuned for the detection of selected wave modes [43]. Ultimately, separation of the two modes allows for normalisation and for a sensor -independent comparison. The relationship between the two wave modes is commutative and henc e, separation will not affect the results. User Input Signal FFT 𝑢𝑛 𝑈𝑘 arg (𝑈𝑘) mod (𝑈𝑘) Phase Velocity Information 𝑣𝑚=෍𝑨⋅sin (𝒌𝐿−2𝜋𝑭𝑡𝑚+𝝓)𝑀 𝑚=1𝑒−𝜇𝐿 𝒌=2𝜋𝑭𝒄𝒑⁄ 𝝓 𝑨 9 Data from the WFR3 software is exported with a sampling rate of 10 Mega Samples per Second (MSPS) . Hence, as a matter of consistency, the sampling frequency of the proposed model has been set to match this. Whils t the sampling frequency from the WFR3 software does not affect the overall results, the sampling frequency in the numerical model is used for the constitutive frequency selection. As such, varying the sampling frequency of the numerical model could have a minor effect on the results of the study. Nonetheless, as the sampling frequency is increased, it produces results that are more like the analytical solution. At a sampling frequency of 10 MSPS, it is expected that the numerical model results should offer a high degree of similarity to the analytical solution. A comparison of three sampling frequencies is shown below for the numerical model output in Fig. 5 . Figure 5. A comparison of the propagated wave using the numerical model with differe nt sampling frequencies. At a value of 1 MSPS it is evident that the wave packet is not smooth . At both sampling frequencies of 10 MSPS and 100 MSPS, the results are near identical, and the wave packets can hardly be distinguished from one another . The eff ect of frequency selection between 1 MSPS and 10/100 MSPS is evident as there is a phase delay between the signals. Ultimately, a sampling frequency between 1 0 MSPS and 10 0 MSPS should provide signal outputs that are akin to the analytical solution and any aliasing effects are negligible. To evaluate the similarity between WFR3 and the numerical model, two methods have been employed. Firstly, the cross -correlation function has been used in the time domain. The cross - correlation function is a method of determining the lag between location of highest correlation of two signals. A small lag indicates that the signals are very close to their highest point of correlation. If the lag value across all the parametric trials is relatively constant, this is a possible indicator that there is a constant phase delay between the two methods but that they otherwise produce very similar results. The cross -correlation is a useful tool in det ermining this phase delay but is not a definite method for determining similarity. The method for calculating cross -correlation is given below for two, time dependent signals 𝑥1(𝑡) and 𝑥2(𝑡) where signal two is a copy of signal one separated by a time de lay (lag) 𝜏12 and attenuation 𝐴. Both signals are assumed to have independent, Gaussian noise components given by 𝑛1(𝑡) and 𝑛2(𝑡). For the parametric study both signals have no noise and hence, 𝑛1(𝑡)=𝑛2(𝑡)=0. 10 𝑥1=𝑠1(𝑡)+𝑛1(𝑡) (13.1) 𝑥2=𝐴𝑠1(𝑡+𝜏12)+𝑛2(𝑡) (13.2) A perfect cross -correlation between two signals can only be obtained over an infinite period and hence, any calculation possesses a degree of truncation . The method for determining the cross -correlation function for a pair of discrete signals is given below where 𝑥𝑛− is the complex conjugate of 𝑥𝑛. 𝜓𝑗=lim 𝑁 →∞෍𝑥1𝑘+𝑗⋅𝑥2𝑘−𝑁 𝑘=−𝑁(14) The cross -correlation can be more efficiently calculated in the frequency domain where the lag can be determined as follows where 𝑋1(𝑓) is the Fourier transform of signal one and 𝑋2(𝑓)− is the complex conjugate of the Fourier transform of signal two. 𝜏12=maxiFFT(𝑋1(𝑓)⋅𝑋2(𝑓)−) (15) The second measure of similarity used is Pearson’s linear correlation. The correlation coefficient is a useful metric in determining a linear relationship between sets of data. As the time information, given as the independent variable, will be constant between data sets, the correlation coeffic ient may be used to compare the amplitude response of the two signals. However, if the signals are not properly aligned due to a possible phase delay, this metric would not prove useful. Hence, the cross -correlation is initially performed and then used to align the signals to the point of their maximum correlation before the correlation coefficient is calculated to determine the magnitude of correlation at this location. In this way, the combination of the two metrics should act as a robust way for determin ing similarity in both amplitude and phase . The Pearson’s correlation coefficient is given below for two signals 𝑥1 and 𝑥2 with means of 𝑥1̅̅̅ and 𝑥2̅̅̅. 𝑟=∑(𝑥1𝑖−𝑥1̅̅̅)(𝑥2𝑖−𝑥2̅̅̅)𝑁 𝑖=1 √∑(𝑥1𝑖−𝑥1̅̅̅)2∑(𝑥2𝑖−𝑥2̅̅̅)2𝑁 𝑖=1𝑁 𝑖=1(16) A comparison of the dispersive relationships calculated from the numerical model and obtained from WFR3 are shown in Fig. 6 . The resultant curves are for aluminium a lloy 2024 - T3. Figure 6. A comparison of dispersion properties of aluminium from the numerical model and WFR3. 11 The Pearson’s linear correlation coefficient has been used to compare the two dispersion relationships . As can be seen , both wave modes are perfectly correlated with a correlation value of one. Hence, the dispersion information from both the proposed technique in the numerical model as well as the WFR3 software can be considered identical. The remainder of the numerical model is validated against WFR3 over a range of material properties and variables as noted in Table 1. These properties are assessed over a range that is suitable for most materials . Symmetric and antisymmetric modes are compared independently as previously discussed. Table 1. The range of variables assessed in the numerical model validation. Lower Limit Spacing Upper Limit Young’s Modulus ( 𝐸) 50 GPa 25 GPa 300 GPa Density ( 𝜌) 500 kg/m3 500 kg/m3 5000 kg/m3 Plate Thickness ( ℎ) 1 mm 1 mm 8 mm Poisson’s Ratio ( 𝜈) 0.05 0.05 0.45 Distance ( 𝐷) 50 mm 50 mm 500 mm The results of the parametric study are shown below in Fig. 7 with phase error, depicted as a time delay on the left and amplitude error, depicted as the Pearson’s correlation coefficient shown on the right. Both sets of data show that the numerical model is comparatively very similar to the results obtained from WF R3. There are notably few outliers. WFR3 struggles to provide a solution as material properties approach non -physical values. Most notably, a Poisson’s ratio of 0.05 produced a bad result. However, the numerical model continues to perform well in this regi me and yields qualitatively accurate results. Nonetheless, the results indicate that the numerical model is capable of accurately calculating the propagation of Lamb waves. Figure 7. Time delay comparison between the numerical model and WFR3 (left) and the Pearson’s correlation coefficient of the signals (right). Antisymmetric modes are shown as solid lines and symmetric as dashed lines. The following sections detail several proposed improvements to the basic version of the num erical model, including how higher levels of fidelity may be incorporated. 12 Modelling Extensions High Temperature Effects Whilst numerous advancements have been made in the fields of NDE and SHM, there remains several challenges, including operation in extr eme environments. High temperature sensors that can withstand such conditions are difficult and expensive to produce . Modelling of high temperature Lamb waves allows for the optimisation of sensing systems to overcome some of these challenges. Aerospace fi elds that would see the biggest improvement from research in this area includes the health monitoring of propulsive systems and hypersonic vehicles. Propulsion systems have recently been shown as the most common cause of fatigue failure in aircraft [44]. Moreover, w hilst hypersonic flight has seen deployment of experimental vehicles (~TRL 5), mass production has been delayed, partly owing to the low reliability of some systems [45]. The extreme environments in which these vehicles operate introduce large uncertainties and sensitiviti es of which the use of SHM can aid in overcoming. Whilst a suitable body of literature exists that explores the thermal effects on Lamb wave propagation, little work has been completed at temperatures of interest (> 1000°C). Within this higher temperature range, it is expected that additional phenomena may be encountered that have not been considered in previous work. Materials for use in this temperature range including nickel -based superalloys, titanium aluminides, and Ceramic Matrix Composites (CMCs) ex hibit phenomena such as viscoelastic behaviour and creep that has not been characterized or modelled [46]. Moreover, accumulation of creep damage results in higher void content and cracking that can further influence acoustic behaviour and further engenders reflection -transmission and mode conversion [46]. First attempts can be made to model high temperature Lamb waves by solving a temperature dependent Rayleigh -Lamb equation. Here, 𝜆𝑇 and 𝜇𝑇 are the first and second temperature dependent Lamé’s constants, respectively. The temperature change Δ𝑇 is taken as the difference between the current temperature 𝑇 and that of a reference temperature 𝑇0. The half - thickness of the plate ℎ with respect to the thickness at the reference temperature ℎ0 is then determined by considering the linear coefficient of thermal expansion 𝛼. ℎ(Δ𝑇)=ℎ0(1+𝛼Δ𝑇) (17) Likewise, the density of the material with respect to the density at the reference temperature can be determined. 𝜌(Δ𝑇)=𝜌0 (1+𝛼Δ𝑇)3(18) The variation in Young’s modulus and Poisson’s ratio were empirically derived by Augereau, et al. [47] for aluminium 6061 -T6 over a temperature range of 20°C to 220°C as given below. Material properties at high temperatures is arguably one of the biggest limitations in this area. 𝐸=77.59−27.03⋅𝑇 103−13.78⋅𝑇2 106(19.1) 𝜈=0.317+54.79⋅𝑇 106+6.5⋅𝑇2 109(19.2) Taking the derivative of 𝑝 and 𝑞 with respect to temperature provides the following relationships. 13 𝑝𝜕𝑝 𝜕𝑇=𝑘𝜕𝑘 𝜕𝑇−𝑘𝐿𝜕𝑘𝐿 𝜕𝑇(20.1) 𝑞𝜕𝑞 𝜕𝑇=𝑘𝜕𝑘 𝜕𝑇−𝑘𝑇𝜕𝑘𝑇 𝜕𝑇(20.2) Where 𝑘𝐿 and 𝑘𝑇 are the longitudinal and transverse wave numbers, respectively. Finally, the group velocity 𝐶𝐺 may be calculated as follows from the phase velocity 𝐶𝑃. 𝜕𝐶𝐺 𝜕𝑇=𝜕𝐶𝑃 𝜕𝑇(1−𝑘 𝐶𝑃𝜕𝐶𝑃 𝜕𝑘)+𝑘(𝜕2𝐶𝑃 𝜕𝑘𝜕𝑇) (21) A plot of the dispersive phase velocities for aluminium 6061 -T6 as well as the corresponding temperature sensitivities are shown in Fig. 8 . It is clear that an increase in temperature results in a decrease in the phase velocity of both wave modes but is more pronounced for the symmetric mode. Figure 8. A comparison of dispersive phase velocities for aluminium (left) and the c orresponding sensitivity (right). The change in the time -of-flight (TOF) of an acoustic signal as well as the corresponding increase in frequency is shown in Fig. 9 . An example of a 200 kHz wave in a 3 mm plate that is transmitted over 250 mm with an increase from the reference temperature of 20 °C would result in a change in the TOF of the symmetric mode by 0.134 µs and the antisymmetric mode by 0.239 µs. Figure 9. Change in the signal time -of-flight (TOF) (left) and corresponding increase in sampling frequency (right). 14 Anisotropic Materials A further extension of the numerical model may include consideration of anisotropic materials. Such materials are commonplace in aerospace settings, primarily in the form of fibre reinforced polymers (FRPs) and cold rolled aluminium [48]. Adjustments can be made to the original model by including these effects in the form of directionally dependent material properties. The model can easily incorporate this modification by calculating the dispersion relationships for each direction of interest. Whilst such a m ethod would account for the variation in wave speed, it would not capture the complex lamb wave and shear horizontal wave coupling effects [49]. To model the anisotropic material properties in FRPs, classical laminate theory (CLT) is used to determine the directionally dependent material properties. We begin with the reduced stiffness coefficients of the material along the ply axis , knowing that 𝜈12 𝐸1=𝜈21 𝐸2 [50]. 𝑄11(0)=𝐸1 1−𝜈12𝜈21(22.1) 𝑄22(0)=𝐸2 1−𝜈12𝜈21(22.2) 𝑄12(0)=𝜈21𝐸1 1−𝜈12𝜈21(22.3) 𝑄66(0)=𝐺12 (22.4) These coefficients must then be transformed along the laminate axis, where 𝑐=cos(𝜃),𝑠= sin (𝜃). 𝑸(𝜃)=𝑻𝑸(0), (23) 𝑻= [ 𝑐42𝑐2𝑠2𝑠44𝑐2𝑠2 𝑐2𝑠2𝑐4+𝑠4𝑐2𝑠2−4𝑐2𝑠2 𝑠42𝑐2𝑠2𝑐44𝑐2𝑠2 𝑐3𝑠−𝑐𝑠(𝑐2−𝑠2)−𝑐𝑠3−2𝑐𝑠(𝑐2−𝑠2) 𝑐𝑠3𝑐𝑠(𝑐2−𝑠2)−𝑐3𝑠2𝑐𝑠(𝑐2−𝑠2) 𝑐2𝑠2−2𝑐2𝑠2𝑐2𝑠2(𝑐2−𝑠2)2] ,𝑸(0)= 𝑄11(0) 𝑄12(0) 𝑄22(0) 𝑄66(0)(24) Next, we calculate the in -plane stiffness matrix terms for the stress resultants . 𝐴𝑖𝑗∗=1 𝑡෍𝑄𝑖𝑗(𝜃𝑘)(ℎ𝑘−ℎ𝑘−1)𝑛 𝑘=1(25) Finally, the stiffness properties of the material along the laminate axis may be obtained. 𝐸𝑥=𝐴𝑥𝑥∗−𝐴∗𝑥𝑦2 𝐴𝑦𝑦∗(26.1) 𝐸𝑥=𝐴𝑦𝑦∗−𝐴∗𝑥𝑦2 𝐴𝑥𝑥∗(26.2) 𝜈𝑥𝑦=𝐴𝑥𝑦∗ 𝐴𝑦𝑦∗(26.3) 𝜈𝑦𝑥=𝐴𝑥𝑦∗ 𝐴𝑥𝑥∗(26.4) 𝐺𝑥𝑦=𝐴𝑠𝑠∗(26.5) Let us consider an orthotropic laminate with the following layup [0, 90, 0, 90] 2S, and material properties given in Table 2 below. 15 Table 2. Material properties o f the composite for the given example. Young’s Modulus 𝐸1 137.4 GPa Young’s Modulus 𝐸2 11.71 GPa Shear Modulus 𝐺12 5.510 GPa Poisson’s Ratio 𝜈12 0.250 Ply Thickness 𝑡𝑝𝑙𝑦 0.200 mm The resulting CLT analysis provides polar plots of the material properties as shown below in Fig. 10 . Figure 10. Polar plots for modulus (left) and Poisson’s ratio (right) dependent upon the material direction of orientation. The resulting phase velocities for b oth the symmetric and antisymmetric modes at 200 kHz are plotted below in Fig. 11 . As evident, the change in material properties has a marked effect upon the acousti c response of the panel and should be considered in any form of analysis. Figure 11. Polar plot of the phase velocities of both Lamb wave modes. 16 Edge Reflections Edge reflections of lamb waves has not been implemented but a proposed addition to the original numerical model is a raytracing method that would predict the path of the Lamb wave through a discrete number of angles , like that sho wn in Fig. 12 . Boundary conditions can be implemented by considering wheth er the plate is fixed or free at an edge, resulting in an inverted or non -inverted reflection, respectively . Likewise, a relationship for reflected energy could be used to predict the amplitude of reflected wav es as well as models for predicting mode conve rsion . Edge reflections are an important phenomenon, as the colocation of waves would result in constructive and destructive interference, affecting the ability to analyse and localise the emission. Figure 12. The general princi ples of ray tracing edge reflections. Previous approaches including Monte Carlo simulations are computationally expensive, and hence, the use of ray tracing algorithms that can accurately and efficiently analyse the path of the wave from source to sink are desirable. Heinze, et al. [51] developed a ray tracing algorithm for anisotropic composite materials based on Fermat’s principle . Zhang, et al. [52] proposed a novel version of the linear travel time interpolation method which was also used to predict and visualise geometric defects in plates. Balvantín, et al. [53] also developed a ray tracing method based on Fermat’s principle by discretising the plate into a grid of constant material properties and applying Snell’s law. A thickness reduction defect was experimentally and numerically studied which showed the ability to reconstruct the tomography of a plate based upon the wave’s time of flight. Conclusion Ultimately , this work has presented a novel first principles numerical model for predicting the behaviour of Lamb wave propagation in 2D plate structures. The model is easier to understand and implement compared to analytical solutions and faste r than traditional numerical methods. Validation against WFR3 shows excellent correlation. The presented numerical model opens opportunities for improved NDT&E and SHM of aerospace vehicles. The study also identifies potential avenues for future research and improvements to the code , including the modelling of high temperature effects, anisotropic materials, and edge reflections. Overall, this research paper has contributed to the field of Lamb wave propagation 17 modelling and provides valuable insights for r esearchers and engineers in the aerospace industry. 18