Numerical study of Rayleigh wave interaction with wedge geometry

– This article aims to study the interaction of Rayleigh or surface wave with a varying angled wedge using numerical simulations. This work uses numerical tools to understand this complex problem and ﬁ lls some existing gaps such as the in ﬂ uence of frequency and geometry of the wedge (curved vs. sharp transition). Quantitative analysis was carried out by calculating the transmission and re ﬂ ection coef ﬁ cients, and qualitative analysis used displacement vector plots to study the scattering and mode conversion phenomenon. The results suggest a strong dependence of transmission and re ﬂ ection coef ﬁ cients on the frequency and geometry of transition, which has not been reported before in the literature.


Introduction
Surface or Rayleigh wave interaction with an elastic wedge is one of the interesting and unsolved classical problems in the geophysics and acoustics area [1][2][3][4][5][6]. Rayleigh surface waves can propagate only along the stress-free boundary of a half-space, with their energy confined to the sub-surface and decaying exponentially along the outof-plane direction [7][8][9][10][11][12]. Rayleigh wave interaction with discontinuities can produce mode conversions which will result in the transmission, reflection and scattering of the incident energy. This interaction has been of significant interest to several communities including geophysics, acoustics and nondestructive evaluation [3,7,[12][13][14][15]. Rayleigh wave propagation and interaction with a wedge like discontinuity can be simplified as two ends of an elastic boundary that meet at a non-orthogonal angle while satisfying stressfree boundary conditions. This approach has been used in geophysics to simplify the complex topography caused by mountains and other landmarks [3,6,13]. These studies provide a gateway to understanding the Rayleigh wave scattering on a seismic scale by analysis via analytical or numerical methods. Several foundational studies by researchers such as Mal [1,3], Knopoff [1][2][3][4], Gangi [4,5], Gautesen [16][17][18][19][20] and more [14,[21][22][23] have been crucial in understanding Rayleigh wave interaction with discontinuities, and especially elastic wedges. An incident Rayleigh wave mode converts into bulk modes upon interacting with the discontinuity. Both the transmission and reflection coefficients of the Rayleigh wave interaction with a step-change and at a corner has been reported by Mal [1,3]. These types of wedges can be similar to a 90°or a 270°wedge. The transmission and reflection coefficients for Rayleigh wave interaction with a varying angled wedge has also been observed experimentally by Mal [3], Knopoff [3,4] and Gangi [4,5]. One of the challenges with an experimental analysis is that it can only be used for validation of the transmitted and reflected wave modes and the limitations of the experimental setups do not allow for quantitative characterization of the scattering phenomenon.
Some of the first elastic wedge models were developed by Gautesen [16][17][18][19][20] using a set of Green's function integrals to capture the transmission (T c ) and reflection coefficients (R c ) and their respective phase for wedge angles ranging from 63°to 350°. Gautesen also carried out this analysis at different Poisson's ratios. They observed that the T c and R c exhibit fluctuations at angles h < 180°. Beyond wedge angles of h > 180°, Tc seems to reduce in value, and Rc experiences a small increase. Mal and Knopoff [1][2][3] also used a set of Green's function integrals derived from Huygens' principle to capture the displacements for the reflected and transmitted waves. This approach seems to reduce the error between the analytical and experimental results. Other researchers such as Fuji et al. [21,22] have proposed wedge models using sets of complex integrals and that the coefficients experienced high fluctuations from wedge angles of 36°-180°. Budaev [14] has proposed Sommerfeld's integrals to capture the transmitted, reflected and scattered energy. The experimental coefficients at small wedge angles less than or equal to 90°exhibited considerable error, thus they only considered wedge angles greater than 95°. One of the many challenges in using analytical models is being able to calculate the coefficients effectively when the wedge angle becomes very small (i.e. h < 90°), such that it would produce many irrelevant reflected waves from the surrounding boundaries [1,3,14,16,17]. Another challenge with analytical models is that the scattering phenomena may not be fully captured. Rayleigh wave interaction with discontinuities can result in a wide range of mode conversions and different incident waves on the half-space boundary. This can be more complicated for the case of an acute wedge, where multiple overlapping internal reflections can cause mode conversions and local constructive and destructive interference. These issues can be overcome by using numerical solvers such as finite element method (FEM) [12] or finite-difference time-domain (FDTD). There is an abundant literature using FEM and FDTD tools for acoustics and elastodynamics problems [24][25][26][27][28][29]. The advantage of using numerical methods is the ability to accurately capture complex boundaries and the interaction of elastic waves with these boundaries. While analytical methods are very fast, numerical solutions can capture full field interactions, which allows us to accurately capture the scattering and mode conversions. More recently, these numerical tools have become very powerful and capable of solving very large problems in acoustics and elastodynamics using tools like graphic processor unit acceleration and parallel computing. Therefore, revisiting the interaction of Rayleigh waves with elastic wedge using powerful numerical solvers can be beneficial to our understanding of the classical wedge problem.
While the previous research in literature have studied the wedge problem with different material properties, they have not studied the effect of frequency. The previous work by Gautesen [16,17] was setup up in such a way that the resulting reflection and transmission coefficients were invariant of the excitation frequency. This is typically expected for the standard reflection and transmission coefficients of both bulk and Rayleigh waves in half-spaces [30]. However, changing the excitation frequency will change the wavelength of the incident Rayleigh wave, which may change the interaction phenomenon primarily due to the amount of out-of-plane energy which interacts with the discontinuity as shown in Figure 1. This might have little effect on wedge angles h > 180°shown in Figure 1b, since the 2nd stress-free edge of the wedge is beyond the incident energy. But for wedge angles h < 180°shown in Figure 1a, the incident wavelength may have a larger effect, since the two stress-free edges wedge converge as shown in Figure 1a [1,3,14,16,17]. Therefore, we can hypothesize that the wavelength of the incident Rayleigh wave will have an effect on the interactions especially for wedge angles h < 180°. A change in wavelength may also affect the phase of the transmitted wave, since as the wedge angle reduces to <90°, superposition and multiple interactions of waves can affect the phase of the transmitted waves. It is also interesting to note that the primary source of reflection, transmission and scattering of an incident Rayleigh wave is the discontinuity (D1 as shown in Figs. 1a and 1b). Most existing literature has used a sharp or discrete discontinuity to transition between the two half-spaces, as shown in Figure 1. Rayleigh wave interaction with a discrete discontinuity will result in bulk mode conversions and re-conversion into Rayleigh waves transmit and propagate along the next half-space. If this discrete transition is replaced by a smooth or curve like transition, then the Rayleigh wave interaction can be different. Several researchers including Gangi [5], Wong [13] and Sánchez-Sesma et al. [31] have explored the various mode conversions of Rayleigh waves when interacting with a changing semi-elliptical geometry. Experimentally, the observations from Gangi [5] was that at smaller wedge angles, the mode conversion efficiency from Rayleigh to bulk are higher. Previous work on Rayleigh wave propagation with curved geometry shows dispersion, i.e. frequency dependent velocity depending on a frequency-radius of curvature criteria [6,[31][32][33]. One can hypothesize that a similar effect can occur at the discontinuity, if the sharp discontinuity is replaced by a smooth curve like geometry.
The objective of this work is to revisit Gautesen's models [16][17][18][19][20] using numerical simulations to validate the previous results and gain a deeper understanding by exploring the effects of frequency and curvature. This allows us to quantitatively compare the present results against Gautesen's results using transmission, reflection and scatter coefficients and the phase of the transmitted and reflected waves. First, the hypothesis of frequency and wavelength dependence on the Rayleigh wave interaction was tested using the numerical models. Next, the influence of curvature on the interaction was tested by developing curved wedge models and comparing the T c , R c and scattered power (P s ) values between the curved and sharp transition models for different wedge angles. The reasons for possible discrepancy between present results and earlier results were analyzed both qualitatively using vector plots and quantitatively using coefficients. Finally, using the numerical data, a functional relationship between transmission coefficient, frequency and wedge angle was obtained for a narrow range of wedge angles, which has not been shown before. When the wedge angle is between 330°and 350°, it can be considered as a horizontal semi-infinite crack which has been previously reported by Chakrapani [34]. Each model was sufficiently large to ensure reflections from the boundaries do not interfere with the incident and scattered wave packets within the prescribed time window. Where the scattered waves are the mode converted bulk waves that propagate through the half-space after interaction at D1. The models were meshed with 8 noded quadrilateral elements with quadratic shape functions to satisfy the convergence criteria of 10 elements per wavelength (k) [34][35][36], and the total solve time and time step was varied to satisfy the Nyquist criteria. When the excitation frequency was changed, the model was re-meshed to re-satisfy the convergence and Nyquist criteria. The equation of motion for the system is given by:

Numerical model
where M is the mass matrix, C is the damping matrix and K is the stiffness matrix. ü, _ u and u are the nodal acceleration, velocity and displacement vectors, respectively, and F(t) is the load vector. Using the standard steel properties, the elemental stiffness matrix was obtained: m = 0.3, E = 210 GPa. The mass matrix is a function of the density: q = 7850 kg/m 3 , and the damping matrix was not used in the present case, so it was set to 0. The objective of the 2D plane strain approximation is to focus only on understanding the interaction and mode conversion phenomenon. A transient analysis was carried out with an initial input of 1 MHz, 7 cycle tone-burst [34,37] with a total simulation time of 50 ls and 50 ns time step. Newmark's time integration scheme [38] was used for time stepping, and the direct sparse solver was used to solve the series of equations. Once the model was solved at each node; the temporal profile of the displacement components was extracted to measure the coefficients. The load vector was prescribed with the timedependent excitation source, which was an out-of-plane displacement (U y ) prescribed on the horizontal boundary at a fixed distance away from D1, such that a discrete reflected wave packet without any overlap could be extracted. The position of the receiver was selected based on the shape of the transmitted wave packet, i.e. a complete wave packet without any additional overlapping waves. The y-displacement is shown in Figures 2a and  2b where, A i , A r and A t are the incident, reflected and transmitted amplitudes. To keep it consistent with earlier work [16,17], the transmission and reflection coefficients were calculated by: T c = A t /A i , and R c = A r /A i . This has been used by several researchers to study Rayleigh wave interaction with discontinuities [39,40]. Since the incident, transmitted and reflected modes are all Rayleigh waves, this definition will be valid. The coefficients can also be calculated using the definition of power flow [30], and we had carried out a validation to ensure that both definitions produce the same value of the coefficients.

Effect of frequency
The results of the numerical simulations carried out at different frequencies: 0.5 MHz (k % 6 mm), 1 MHz (k % 3 mm) and 2 MHz (k % 1.5 mm) are shown in Figures 3a-3c. Gautesen had previously presented the coefficients as a function of wedge angle from 63°to 350° [16,17]. The present results have been overlaid with Gautesen's It can be observed that the coefficients also show a dependence on frequency. The T c at 2 MHz agrees well with Gautesen's results for a wider range of wedge angles rather than just h > 180°. The R c and P s at 2MHz also match with Gautesen's results. The largest difference can be observed at 0.5 MHz, which does not seem to match with Gautesen's results. These results confirms the hypothesis that the coefficients have a frequency dependence.

Effect on phase
In the previous work by Gautesen [16], the effect of wedge angle on the phase of the transmitted and reflected waves was also carried out. For a quantitative comparison of the phase term at different wedge angles, it is important to extract the phase at a constant separation between the source and the transmitter location for all wedge angles. However, since the transmitter location was chosen to minimize overlap or superposition, any quantitative difference in phase term will be lost. Therefore, this work uses a modified approach of using the time-of-flight (TOF) data to indirectly determine phase differences. First, the numerical TOF was calculated by performing a Hilbert transform [41] on the time-amplitude data shown in Figures 2a and 2b by a red dotted line. From the Hilbert data, the TOF of the transmitted wave can be calculated by taking the time corresponding to the maximum amplitude for the incident (t 1 ) and transmitted (t 2 ) positions shown in Figures 2a and 2b using DT num = t 2t 1 . The numerical TOF of the reflected wave mode can be calculated using DT num = t 3t 1 , where t 3 is the time corresponding to the maximum amplitude A r of the reflected wave shown in Figure 2a. The theoretical TOF of the transmitted wave can be determined by calculating the distance of the incident Rayleigh mode x 1 to the discontinuity D1, then to the transmitted mode x 2 shown in   Figure 4 using DT th = (x 1 + x 2 )/V r , where V r is the Rayleigh wave velocity [7,42] for a given material. The difference in the theoretical and numerical TOF: DT = DT num À DT th , will be proportional to the change in the phase, which can be calculated using: where DT is the difference between numerical and theoretical TOF as a function of wedge angle (h) and f is the excitation frequency. Figures 5a and 5b show the phase of the reflected and transmitted wave mode as a function of wedge angle for 0.5 MHz, 1 MHz, and 2 MHz. For the transmitted wave mode, the phase shown in Figure 5a shows relatively little to no change from a value of Àp/2 to p/2 across all frequencies except at h = 60°for 0.5 MHz. The reflected wave mode shown in Figure 5b is observed to show a similar effect as the transmitted wave mode from a wedge angle of 60°< h < 150°and 240°< h < 350°across all excitation frequencies. For wedge angles between 160°< h < 230°, there is a significant phase change of À3p/2 at a wedge angle of h = 170°at 0.5 MHz; and a phase change of 2p at a wedge angle of h = 190°at 1 MHz. This shows that phase of the transmitted and reflected waves are also influenced by the wedge angle.

Effect of curvature
To evaluate the effect of the geometry of the discontinuity, a curved geometry as the diffraction source (D2) was modeled as shown in Figure 6. The same material properties described in Section 2 were used and, 0.5 MHz was chosen as the excitation frequency based on its large Rayleigh wavelength described in Section 3.1.
The radius of curvature (R) or the fillet radius was created using a Computer Aided Design (CAD) program, NX. First, a prescribed distance of 1 Â k from the discontinuity D1 was set along the two connecting boundaries shown in Figure 6. Next, the fillet radius was determined as a result of the two end points of the fillet satisfying the tangential criterion at each line segment denoted by a light blue dashed line in Figure 6. Finally, the initial value for the fillet radius R was scaled by an arbitrary factor (F) to get a smooth transition instead of an abrupt transition. The dashed lines shows where the sharp discontinuity D1 is in relation to the fillet D2 based on F Â k is also shown in Figure 6.
A total of four different wedge angles were modeled and tested with a factor F = 2.5 and the resulting radius of curvature values R are shown in Table 1.
Comparing the T c values in Figure 7a, it can be observed that the sharp discontinuities have the lowest T c compared to smoother discontinuity with the exception of 80°. For angles h > 180°, the difference between sharp and smooth geometry is smaller. The wedge with a D1 diffraction source has a higher R c and P s values as shown in Figures 7b and 7c compared to the curved geometry.   Finally, comparing the TOF values, with the exception of 80°, the curve geometry shows no change in TOF as a function of the angle. However, the TOF of the sharp geometry, changes upto 5% as a function of wedge angle.

Superposition of waves
One of the reasons for the difference between Gautesen's results and the present results could be due to the superposition of the waves. When extracting the U y displacement, the position of the receiver node was chosen such that a single Rayleigh wave packet without superposition or interference can be extracted. The discrete wave packets will look similar to the wave packets shown in Figures 2a and 2b. The local interference between multiple waves after the discontinuity at different time steps is shown in Figures 8a and  8b. The vector plots were obtained by extracting the displacement vectors at a given time instance for the entire half-space. As observed in Figure 8a after interaction the compression (P), shear (SV), and Rayleigh (R) wave modes have not fully separated into clear packets, which is referred to as "undeveloped". In Figure 8b, the transmitted Rayleigh (R t ), reflected Rayleigh (R r ), and shear (SV) wave modes have clearly separated, which is referred to as "fullydeveloped".
By following this criteria of the transmitted wave being undeveloped vs. fully-developed, Table 2 compares the difference in T c between the two conditions and Gautesen's results for 110°wedge angle. It can be observed that the T c of the undeveloped wave approximately matches with Gautesen's results, compared to the fully developed wave. This suggests, that the "undeveloped" wave could possibly be the reason behind the discrepancy for the T c , R c values. However, this is valid only at 1 MHz or lower, and at 2 MHz, the T c values seem to match well with Gautensen's work. The increase in wavelength at lower frequency could possibly result in higher spatial superposition. But, while extracting displacements, only the fully-developed i.e. discrete wave packets were used for the analysis.

Effect of frequency
From the results shown in Section 3.1, it can be observed that there is a frequency dependence on the coefficients at certain wedge angles with the largest excitation frequency being 2 MHz. For h > 180°, there seems to be no dependence on the frequency, and the T c values seem follow closely as observed in Figure 3a. For 100°< h < 180°, the 0.5 MHz and 1 MHz results match closely, and the 2MHz results matches with Gautesen's results. However, looking at the R c results for 2 MHz in Figure 3b, at a wedge angle of 60°the value seems to exceed a value of 1 which violates the conservation of power. This can possibly be due to constructive interference or overlapping of modes, which can violate the conservation principles. While the transmitted waves modes were meticulously chosen such that discrete wave packets without overlap were extracted, the reflected modes could possible have an overlap due to the relative shorter distance between the discontinuity and the reflected mode. Interestingly, this violation of the conservation principles can be found in Gautesen's work at wedge angles of 60°and 70°as well. Understanding this will require a combined numerical-analytical analysis which can decompose any overlapping modes, which is beyond the scope of this work.
The values for T c between 100°< h < 180°, have been plotted as a function of the excitation frequency as shown in Figure 9. The resulting data as shown in Figure 9 was then fitted with different functions (such as power law, polynomial etc.) to determine a functional dependence between T c , f, and h. For the T c versus h relationship, a power law fit had the best coefficient of determination (R 2 ), and a quadratic dependence was observed; T c / h 2 . To determine a functional relationship between T c and frequency, the data in Figure 9 was fitted with an exponential function, which was once again found to have the highest R 2 value. The frequency dependence was determined to be: T c / e f/3 . For 60°< h < 90°the T c follows a weak parabolic dependence with frequency, but does not show any clear dependence. This could possibly be due to the strong interaction of the Rayleigh wave between the two half-space surfaces.

Phase term
As the wedge angle decreases, the number of reflected wave modes increases [1,3,14,16,17] which can cause multi-mode interactions between the transmitted Rayleigh and the re-reflected wave modes. At 1 MHz, the change in phase for a varying angled wedge of the transmitted wave mode shown in Figure 5a, displays no significant change when compared across all excitation frequencies except when h = 60°. When comparing 1 MHz to 0.5 MHz and   2 MHz, there is a phase difference of roughly À3p for h = 60°. This could possibly be attributed to the constructive and destructive interference of the many re-reflected wave modes that occur at shallow wedge angles. Spatially, these constructive and destructive interferences could result in delayed or advanced formation of a coherent Rayleigh wave, since this work uses the peak amplitude to determine the TOF. A similar trend can be observed for the reflected wave modes in Figure 5b, where there is no significant difference in phase across all excitation frequencies except at 170°< h < 190°where the phase varies from À3p/2 to 2p. Compared to earlier work by Gautesen [16], the trend of the phase plots do not seem to agree. This could be due to the use of phase of the transmission and reflection coefficient rather than the current definition, which is expected to produce a different result.

Geometry of the discontinuity
The influence of the type of discontinuity on the Rayleigh wave interaction can be observed in Figures 7a and 7c. The sharp transition can act as a single diffraction source resulting in higher scattering. This is evident in the increase of R c and P s values for the sharp geometry vs. curved geometry. The R c and P s values follow an inverse trend compared to T c , which is once again consistent with the hypothesis of higher scattering from sharp discontinuity. However, at 80°and 300°, this observation is not true, with the curved geometry having a higher P s value. To further analyze this qualitatively, displacement vector plots were extracted for different wedge angles as shown in Figure 10a-10h. The mode conversion of an incident Rayleigh wave into bulk modes (L, S), and the diffraction of the bulk modes is well captured in the vector plots.
It is interesting to note that in all the sharp discontinuity (D1) plots, the directivity of the bulk modes, specifically the shear wave window is more narrow than the curved case. The reflected bulk and Rayleigh modes are also higher in amplitude in the case of D1. Interestingly, this phenomenon seems to be also dependant on the angle. For example, in the case of h < 90°and h > 270°, the diffraction of the bulk waves seems greater for the curved geometry than the sharp geometry. It can be hypothesized that since the wave is propagating from left to right towards the discontinuity, in the vector plot at certain angles, the mode conversion of incident Rayleigh waves into bulk modes will be based on the curvature. Mode conversion of Rayleigh waves into bulk modes has been explored extensively by several authors [34,43,44]. The geometry of the wedge can be approximated to a semi-infinite crack when h ? 360°. Previous studies show that an incident Rayleigh wave will diffract into bulk modes with a strong shear window. The curved discontinuity can be approximated into a series of step discontinuities. This will result in each point along the curvature acting like a diffraction source, unlike the sharp discontinuity, where only D1 acts as a diffraction source. As seen in Figures 10a, 10b and 10g, 10h, when the pointing vector of the scattered mode is opposite to the incident mode, i.e. backward propagating waves [34], (not be confused with reflected waves), the curved discontinuity seems to scatter/diffract bulk modes along the entire curve. This is evident in the shear window being much broader for the curved geometry as opposed to the sharp geometry. On the contrary, when the pointing vector is parallel with the incident angle as shown in Figures 10c, 10d and 10e, 10f, the overall scatter from the discontinuity is less, and the shear window is smaller as well. This qualitatively agrees with our hypothesis of the effect of the curved geometry. Unfortunately, it is challenging to do a comprehensive quantitative analysis because establishing a relationship between curvature radius, wedge angle and Tc can be tedious due to the large number of combinations that are possible. The objective of this study is simply to show the effect of curvature on the mode conversion and scattering for wedge like geometry. An analytical model that can capture these effects would be a more efficient solution compared to numerical analysis, which are more suited to capturing the full displacement field. The scatter coefficient can be evaluated from the full scatter field. However, capturing the mode conversion into bulk modes and the power coefficients as a function of angle as performed earlier by Chakrapani [34], is not possible due to the number of spatially, and temporarily overlapping modes. The superposition of the modes makes it impossible to calculate the power flux of the scattered modes.
The current work demonstrates the need for better analytical models which can account for various mode interaction effects. Future work focusing on developing improved and robust analytical models can be very useful especially when exploring a wide range of parameters. For example, different values for F for curved geometries, or the effect of the frequency on F. Both of which have not been explored in the current work due to the limitations on the problem size and number of combinations of frequency, wedge angle and fillet radius, which can be computationally expensive. Therefore, an analytical model which can accurately capture the effects is more desirable. Understanding the mode conversion phenomenon can be useful in a wide range of fields including geophysics, acoustics and nondestructive evaluation.

Summary
This article revisits the classical problem of surface/ Rayleigh wave interaction with wedge like geometry using numerical methods. The objective of the article was to study the transmission reflection and scatter of Rayleigh waves by wedge like discontinuity. Two hypothesis were tested: (a) the effect of frequency on the interaction, and (b) effect of discontinuity type, i.e. curved vs. sharp discontinuities. Quantitative analysis was carried out by extracting the T c , R c , and P s , and the apparent phase of the transmitted and reflected waves for different wedge angles and frequency. When compared to previously reported results, a frequency dependence was observable at certain wedge angles. Further, extracting functional relationships for T c showed that: T c / e f/3 h 2 , where f is the excitation frequency and h is the wedge angle. Furthermore, the phase of the transmitted and reflected waves also seems to have a frequency dependence. But the phase vs. wedge angle does not show any clear relationship. Finally, introducing a curved transition; D2 was shown to affect the interaction phenomenon. Interestingly, the effect including higher scattering/diffraction was more pronounced at certain wedge angles (h < 90°and h > 270°). This was shown qualitatively using displacement vector plots, which showed a change in the size of the scattered shear window due to the curved geometry.