Attenuation of Rayleigh waves due to three-dimensional surface roughness: A comprehensive numerical evaluation

Rayleigh waves exist on the surfaces of elastic solids, at traction-free boundaries. Surfaces of propagation are well known to attenuate the Rayleigh waves, due to the presence of roughness. Since all solids unavoidably possess a certain degree of roughness,¹ it is important to understand how this roughness influences the propagation of Rayleigh waves, especially when attenuation is used as a nondestructive evaluation measure.^2–5 Additionally, good understanding of this phenomenon can lay the groundwork for utilizing it as an efficient method for characterising the degree of roughness of a material, when it is unknown.

There exist multiple theoretical studies that derive expressions predicting the attenuation coefficient of a Rayleigh wave traveling over a rough surface. Some of the most well-known results come from the work in Refs. 6 and 7. In their work, the authors derived expressions relating the attenuation coefficient of Rayleigh waves with the frequency (f) and the root mean square (RMS) height (δ) and correlation length (Λ) of the surface. The latter two variables are metrics characterising the height of the peaks/troughs of the rough surface and the spacing between them, respectively. The authors show that in the long wavelength limit, known as the Rayleigh regime, the attenuation coefficient is proportional to $δ 2 Λ 2 f 5$⁠. The results from these authors were obtained using small perturbation methods; in Ref. 8, the same power relationships were retrieved using the Rayleigh–Born approximation.

Most of the literature in rough surface scattering focuses on the study of the Rayleigh (long wavelength) regime, which was discussed in the previous paragraph. As the frequency increases, the findings in Ref. 7 and, later, Ref. 9 show that in the stochastic regime (medium wavelength), the attenuation coefficient is proportional to $δ 2 Λ − 1 f 2$⁠. Finally, in the geometric regime (short wavelength), we rely on the findings of our previous work in Ref. 10 combined with the dimensional analysis in Ref. 11 and the results in Refs. 12 and 13, where it was shown that the attenuation coefficient is proportional to $Λ − 1$ and independent of δ and f.

However, experimental validation of these relationships has proven unfeasible due to a combination of the difficulty in reproducing the required rough surfaces in the real world and the high attenuation values that they cause, which makes the waves very difficult to detect in experimental investigations. Additionally, the theoretical models require the use of various approximations and mathematical techniques to predict the attenuation coefficient in each scattering regime, further increasing the complexity in the predictions. In our previous work,¹⁰ we were able to verify the relevant theory in two-dimensional (2D) space and expand it to roughness combinations beyond its limits of validity. The verification was completed using results from finite element (FE) modeling. FE modeling allowed us to complete the verification across all three scattering regimes under a single, unified approach, as we were able to use the same model regardless of the scattering regime and the roughness parameters. Additionally, physical assumptions beyond the fundamental discretisation of the underlying displacement field, which can be addressed with suitable mesh refinement, are not made during FE modeling—this eliminates any additional assumptions that were made in past theoretical work.

Completing a three-dimensional (3D) study previously was considered not to have been possible due to the computational burden of solving multiple 3D models to determine roughness parameters suitable for our study, as, to our knowledge, there did not exist another FE study that could be used as a reference. Now, having identified particular regions of parameter space that are accessible by our models and that satisfy the theory's requirements, and tuning our models, we extend to three dimensions. The results here both strengthen our findings in Ref. 10 and allow us to gain further insight with regard to the theory in Refs. 7 and 9 in 3D space, eliminating the need for using analogies between the 2D and 3D theories. Additionally, we demonstrate the power relationships between the attenuation coefficient and the surface's statistical parameters and frequency in the geometric regime, which, to our knowledge, is relatively unexplored in terms of scattering from rough surfaces; in our 2D work, we used an analogy with grain scattering and dimensional analysis to derive the expected power relationships; however, here we are able to investigate this relationship and the analogy directly. Also, following the demonstration in all scattering regimes similarly to Ref. 10, we extend our results to regions beyond the postulated limits of validity of the existing theory, showing the capabilities of FE modeling to produce results for a plethora of roughness cases without the relatively low roughness restriction (⁠ $δ / Λ < 0.3$⁠) that the theory requires. This advancement in understanding lays the groundwork for either future use of FE modeling for predicting (and correcting for) the attenuation coefficient caused by a given roughness profile or an inverted approach, where the attenuation coefficient is used to quantify the surface roughness of a component, without the relatively low roughness restriction imposed by the theory. It is also worth noting that the theory in Refs. 7 and 9 also predicts a frequency shift caused by the roughness; however, our work here focuses solely on the behavior of the attenuation coefficient. Dispersion and diffraction effects can be also formulated in the general framework of a “global theory” based on dimensional analysis and similitude, similarly to our 2D work in Ref. 10. From the practical point of view of ultrasonic nondestructive evaluation (NDE), those effects are less significant than the attenuation effect studied in this paper. However, there cannot be any doubt that the numerical technique presented in this paper is suitable to study those effects as well, although sufficiently accurate assessment of weaker effects naturally requires higher numerical accuracy (i.e., higher meshing refinement, smaller time steps, and hence more computational power etc.).

This paper is organised as follows: Sec. II presents a summary of the theory used in our work, and Sec. III explains the details behind generating the FE models. Then Sec. IV presents the asymptotic power relationships obtained from solving the FE models, between the attenuation coefficient and δ, Λ, and f, comparing them with theoretical predictions, and finally, Sec. V summarises the work.

In this section, we present a brief summary of the theory used in this work; first, we discuss the relevant theory for defining and generating statistically rough surfaces, for later use in our FE models, and subsequently we present the theory and resulting power relationships for describing the behavior of the attenuation coefficient in the three scattering regimes.

A surface can be thought of as rough when its local height over a reference plane varies with space. There are two statistical parameters with which the roughness of a surface is often mathematically quantified—the RMS height δ and correlation length Λ. The former parameter can be thought of as a metric of the height/depth of the peaks and troughs of the surface, while the latter is a measure of the spacing between them.

As shown in Fig. 1, this method generates a rough surface with smooth peaks and troughs, and the windowing reduces the amplitude of those to zero, at the edges of the rough zone, to ensure a smooth joining with the smooth parts of the FE domain.

In the stochastic regime (medium wavelength), analytical expressions for the behavior of the attenuation coefficient have been derived in Ref. 9; the authors have shown that as the frequency increases, the $α R ∝ δ 2 Λ 2 f 5$ proportionality becomes $α S ∝ δ 2 Λ − 1 f 2$⁠, where the superscript S denotes the stochastic regime. This power relationship is also inferred in Ref. 7, as the ω₂ function can be seen to vary linearly with $ω Λ$⁠, which, when substituted in Eq. (8), yields the same power relationship.

It is worth noting that the power relationships between the attenuation coefficient and the statistical parameters/frequency are expected to be identical between two and three dimensions in the stochastic and geometric regimes. This phenomenon has been discussed in Ref. 11 and demonstrated in our work in Ref. 10, where we have used 3D theory to predict the 2D power relationships, which we were subsequently able to verify. Therefore, in the geometric regime, it is expected that $α G ∝ Λ − 1$⁠, where the subscript G denotes the geometric regime.

A summary of the expected power relationships between α and δ, Λ and f, which our work attempts to demonstrate, is shown in Table I.

We used the 2D model in our previous work in Ref. 10 as a basis for constructing the 3D model for this study. In the 2D model, a rough surface was inserted at the bottom boundary of a rectangular FE domain, with monitor nodes on both sides of it, and a source line, composed of multiple source nodes, further to the left. Source nodes with pre-determined displacements were used to excite waves in a model, and monitor nodes measure displacements that are stored by the FE solver and can be used during post-processing. For the initial 3D model with flat surfaces, prior to the insertion of the rough area, we used an “extruded” version of the 2D model, i.e., the rough surface and source lines were converted into rough and source areas, respectively, while the monitor nodes were converted into monitor lines. A planar schematic of this configuration is shown in Fig. 2. Similarly to our 2D work, the source nodes were chosen to be outside the rough area, to ensure that a clean Rayleigh wave would arrive at the rough area; although placing the source area within the rough region would be more realistic, it would prevent the consistent creation of a clean Rayleigh wave. It should also be emphasised that our goal is not to perfectly mimic a real-world experiment, but instead to establish the underlying principles of wave interaction with surface roughness. All of these features were defined on the top surface of a rectangular block, while all other surfaces and the volume of the domain were left unaffected.

Obtaining statistically and ergodically meaningful results requires averaging the attenuation coefficient values over multiple simulations, with identical statistical roughness parameter values but unique rough surface profiles. This can create significant computational burden, especially for 3D simulations. As a result, the dimensions of the FE domain were selected dynamically based on the frequency (hence, wavelength) of the simulation and Λ, instead of being fixed, to create a model with the minimum possible volume while ensuring adequate spacing for the source, monitor, and rough and absorbing layer regions. The dimensions of the FE domain are also included in the schematic of Fig. 2.

As shown in Fig. 2, the dynamic model sizing was implemented by defining the size of the domain to be a function of the wavelength and the correlation length; a script would accept the frequency and correlation length as an input and then calculate the suitable domain size, as per the dimensions shown in Fig. 2, which was subsequently used to generate a mesh of that size. Namely, absorbing layers with a thickness equal to $3 λ R$ were defined around the FE domain, to absorb any unwanted reflections,^23,24 where λ_R is the Rayleigh wavelength. Additionally, the theory in Refs. 7 and 9 assumes a space of infinite extent in both the x and y directions, and therefore, the addition of absorbing layers around the model more closely models this scenario. A square source area, with side length equal to $4 λ R$ was defined at a distance of one λ_R from the right edge of the left absorbing layer; the size of this gap was measured from the left edge of the source square, and not its centre, as demonstrated in Fig. 2. Then the monitor line before the rough surface was placed at a distance of $3 λ R$ from the right edge of the source area. The length of the rough surface itself was set to be at least equal to $50 Λ$⁠, as advised by Ref. 1. Another monitor line was placed just after the rough surface, and then, after a gap of length equal to λ_R, a $3 λ R$ thick absorbing layer was defined. For the width of the domain, a λ_R gap was defined above and below the source area, followed again by $3 λ R$ thick absorbing layers. The thickness of the domain (z direction) was set to be at least $3 λ R$⁠, so that the boundary opposite the source area (i.e., the top boundary) would not interfere with the Rayleigh wave.

By following this approach, where all dimensions of the domain were not absolute but rather were determined by the particular frequency and correlation length that were used each time, it was possible to create the minimum possible model size for each particular roughness/frequency case, reducing the computational expense of our simulations. A typical, moderately sized model was of the order of $5.5 × 10 7$ degrees of freedom and required approximately 6–8 min to solve using three Nvidia (Santa Clara, CA) RTX 2080 Ti cards, with 11 GB of memory each.

Steel was used for all simulations (Young's modulus, E = 200 GPa, density, ρ = 7800 kg/m³, and Poisson's ratio, $ν = 0.29$⁠) with longitudinal wave speed C_L = 5797 m/s and shear wave speed C_T = 3153 m/s. The Rayleigh wave speed, calculated from the approximation in Refs. 19 and 25, was found to be 2915 m/s. Here it is worth noting that the material choice is not expected to affect the asymptotic power relationships, as the theory suggests that they are independent of it.

When this map contains values corresponding to a 2D projection of a 3D rough surface, it creates a mesh that is locally rough, within the bounds shown in Fig. 2. An illustration of this process is shown in Fig. 4.

In Fig. 4(a), the lowest z position of an undeformed mesh is plotted at each (x, y) location; as the mesh is undeformed, and in the shape of a cuboid, all points have the same z value. Figure 4(b) shows the matrix that contains the values determining the desired final z position—those warping factor values are given as a ratio between the final and the initial position as per Eq. (10). A value of this ratio greater than one indicates extrusion, while a value smaller than one indicates compression.

Also, given that we required our model to only be rough within a predefined area, as shown in Fig. 2, we populated the warping factor matrix in Fig. 4(b) with ones everywhere besides the desired rough region. Then Fig. 4(c) shows the lowest z position of a deformed mesh at each (x, y) location, where the deformation has occurred as instructed by the warping factor map in Fig. 4(b); the resulting z value at each point in Fig. 4(c) is simply the result of multiplying the respective value in Fig. 4(a) with the respective warping factor in Fig. 4(b). In the FE domain, the mesh in each (x, y) location is warped by displacing the nodes from their initial position, according to the value of the warping factor at each location, distorting the mesh to achieve the desired z value at that location.

This number of realisations was deemed sufficient for the purpose of our study, in all scattering regimes, as a convergence study showed that the value of the attenuation coefficient had converged prior to the 30 realisation limit. An example of a convergence plot, illustrating the evolution of the ensemble average attenuation coefficient, at three roughness scenarios (each corresponding to one of the three scattering regimes) is shown in Fig. 5.

As shown in Fig. 5, the attenuation coefficient appears to have converged within 30 realisations; this relatively fast convergence of ensemble averaging is due to the fact that the length and width of the rough area was large enough for fairly efficient spatial averaging even without significant ensemble averaging. Therefore, this number was judged to be sufficient for the purposes of our study.

As discussed in Sec. II, in the Rayleigh regime, we expect $α R ∝ δ 2 Λ 2 f 5$⁠. For uniformity with our work in Ref. 10, we plot the FE results in their normalised form where possible, i.e., we plot α_n vs δ_n, and α_n vs Λ_n, where $α n = α λ R , δ n = δ / λ R , Λ n = Λ / λ R$⁠, and λ_R is calculated at the central frequency of the simulation. This does not affect the power relationship between α and either δ or Λ, as all simulations in the respective figure are completed at the same frequency, and therefore, all values are normalised by the same constant. For asymptotic analysis relating to f, we are plotting against absolute values; this is because the frequency is not constant in this case, and therefore, the wavelength varies.

The results relating to the Rayleigh regime are shown in Fig. 6, with the results relating to δ, Λ, and f shown in Figs. 6(a)–6(c), respectively. For investigating the $δ 2$ relationship, Λ was fixed to 400 μm, f was set to 1 MHz (⁠ $λ R = 2.92$ mm), and δ was varied from 100 to 250 μm. This normalisation value (⁠ $λ R = 2.92$ mm) was used in all subsequent analysis where the frequency was kept constant, as all such simulations used a Hann-windowed signal with a central frequency of 1 MHz. The gradient of the best fit line, on a log-log scale, in Fig. 6(a) is 2.09, which is close to the theoretical value of 2. For the $Λ 2$ relationship, δ was fixed to 100 μm, f was set again to 1 MHz, and Λ was varied from 100 to 350 μm. The gradient of the best fit line from this set of simulations, as shown in Fig. 6(b), was found to be 1.75, which is again close to the theoretical value. Finally, for the f⁵ relationship, the statistical parameters of the random rough surfaces were set to $δ = 100$ and $Λ = 300$ μm, while f was varied from 0.7 to 1.2 MHz (λ_R varying from 4.16 to 2.43 mm). As shown in Fig. 6(c), a gradient of 5.00 was retrieved from the best fit line.

It appears that the FE results match the theory closely, for all three power relationships. In this regime, as discussed in Sec. II, $ω 2 ∝ ( ω Λ ) 4$⁠, which, when substituted in Eq. (8), yields the $δ 2 Λ 2 f 5$ proportionality; physically, in the low frequency limit of this regime, the wavelength becomes so large, compared with the features of the rough surface, that the rough surface appears flat to the Rayleigh wave, and the attenuation coefficient tends to zero.

Here, it is important to add that a beamspread correction was applied to these results. The theory in Ref. 7 assumes a plane wave for the derivations. We have found that a square source with a side length of $4 λ R$ was sufficiently close to a plane wave source; however, in the Rayleigh regime, where the attenuation values are small, the small beamspread effect was no longer negligible. To account for this, a model with no roughness was run for each frequency and x_r case, and the beamspread was converted into an equivalent attenuation coefficient, which was then subtracted from the value obtained from each rough model, to give the true α value. There is no roughness-induced attenuation in a flat model; therefore, any loss in Rayleigh wave amplitude is solely due to beamspread. The principle of similitude, as discussed in this context in Ref. 10, implies that the rough case must suffer the same loss in amplitude due to beamspread, in addition to the attenuation caused by the roughness. For the frequencies presented here (0.7–1.2 MHz), the beamspread correction values were in the range from 0.91 to 2.96  $m − 1 .$

In the stochastic regime, the theory predicts that $α S ∝ δ 2 Λ − 1 f 2$⁠. The FE results relating to the stochastic regime are shown in Fig. 7, with the results relating to δ, Λ, and f shown in Figs. 7(a)–7(c), respectively. For the investigation of the $δ 2$ relationship, Λ was fixed to 800 μm, f was set to 1 MHz, and δ was varied from 100 to 400 μm. For the $Λ − 1$ relationship, δ was fixed to 200 μm, f was set to 1 MHz, and Λ was varied from 800 to 1600 μm. Finally, for the f² relationship, δ and Λ were fixed to 200 and 800 μm, respectively, while f was varied from 0.8 to 1.8 MHz (λ_R varying from 3.64 to 1.62 mm).

The gradients of the best fit lines were found to be 1.86, –1.06, and 2.15 for δ, Λ, and f, respectively, all of which agree well with the values proposed by the theory. This shows that our model was able to retrieve all the power relationships in the stochastic regime too.

In the geometric regime, the theory predicts that $α G ∝ Λ − 1$⁠, with the attenuation coefficient being independent of δ and f. The results relating to the geometric regime are shown in Fig. 8. For the investigation of this power relationship, δ and f were fixed at 600 μm and 1 MHz, respectively, and Λ was varied from 2000 to 3000 μm.

As shown in Fig. 8(b), the inverse proportionality between the attenuation coefficient and the correlation length is demonstrated well by our FE results, as the gradient of the best fit line is −0.93. This shows that our model is able to encapsulate the physical phenomena of the geometric scattering well; in the geometric regime, the wavelength has become so short that the only important metric in the scattering is the number of peaks/troughs, and not their size, as the wave can traverse their height without getting scattered from them. A long correlation length implies fewer obstacles, and therefore, the attenuation coefficient reduces with an increase in correlation length. This is also reflected in the power relationships relating to δ [Fig. 8(a)] and f [Fig. 8(c)], which have both significantly diminished, from their value of 2 in the stochastic region to values closer to zero.

A summary of our results, compared with the theoretical asymptotic power relationships is shown in Table II.

As shown in Table II, in the Rayleigh regime, the powers of $δ$⁠, $Λ$⁠, and f were 4.5%, −12.5%, and 0.0% away from the theoretical value, respectively. For the stochastic regime, the same differences for the same variables were equal to 7.0%, 6.0%, and 13.0%, while for the power of Λ in the geometric regime, the difference was found to be 7.0%. Overall, we were able to achieve good agreement between our FE results and the theoretical asymptotic power relationships in all scattering regimes. Also, we were able to show that the power relationships hold true and preserve their theoretical value even in cases where the statistical parameters fall outside the region of validity of the theory, i.e., for cases when $δ / Λ > 0.3$⁠.⁷ 

In this work, we have completed a study of the asymptotic behavior of the attenuation coefficient of a Rayleigh wave with respect to the RMS height and correlation length of a statistically rough surface and the wave's frequency. Results were obtained in all three scattering regimes (Rayleigh, stochastic, and geometric), using high-fidelity FE modeling. By utilising the findings in our previous work,¹⁰ we were now able to complete this investigation in 3D space.

We found good agreement between our results and the values proposed by the theory in all scattering regimes. As FE modeling makes no further assumptions apart from the discretisation of the domain into finite areas/volumes and the approximation of derivatives as quotients, it provides robust further insight into the theory, which has been derived under different assumptions. Additionally, we were able to provide further confidence to the $Λ − 1$ relationship between the attenuation coefficient and the correlation length, which was an unknown variable in our previous work. Our results also extended to roughness scenarios that lie outside the regions of validity of the small perturbation theory, showing that the power relationships hold true even in those.

By utilising FE modeling, we were able to obtain these results under a unified approach, as we used the same FE modeling techniques regardless of the scattering regime. This contrasts with previous theoretical work, which requires different mathematical tools and approaches, tailored to the scattering regime being investigated. Also, FE modeling removes the necessity of a priori knowledge of the ω₂ function proposed in Ref. 7, which is complicated to calculate, unique to each roughness/material case, and potentially impractical to use in a realistic environment. Additionally, 3D modeling removes the need for using analogies between 2D and 3D theory and allows for direct investigation of real-world scenarios. Finally, this agreement between FE modeling and the existing theory lays the foundation for potentially using solely FE results for attenuation coefficient predictions for roughness cases where the theory is limited and experimental measurement is difficult.

In terms of NDE applications, if the statistical parameters characterising a surface are known, FE modeling can be used to predict the corresponding attenuation coefficient and, hence, apply corrections to experimental measurements; such a capability has been demonstrated here, as it is possible to replace the random rough surfaces used in our study with rough surfaces obtained experimentally. This can be useful for cases where attenuation is used as an NDE metric. For instance, it is possible to use attenuation for fatigue state characterisation;^3,26,27 however, the presence of roughness would result in higher attenuation values than what can solely be attributed to fatigue and, thus, erroneous fatigue life predictions. Additionally, an “inverted” approach could be used, where an experimentally measured attenuation coefficient is resolved into the corresponding severity of surface roughness. Such an assessment is important for parts made with additive manufacturing, as their layer-by-layer creation results in inherent surface roughness,^28–31 which affects their performance.^32–34 Results for both examples are attainable with the existing theory; however, the FE method presented here removes the low roughness restriction imposed by it and allows for the prediction of the attenuation coefficient in a unified approach, regardless of the magnitude of the scattering regime and the statistical parameters characterising the rough surface. Finally, numerical modeling can quite readily handle cases where the roughness exists on curved components (e.g., shot-peened cylindrical shafts), for which currently there are no theoretical models predicting the attenuation coefficient.

Georgios Sarris was funded by the UK Research Centre in Nondestructive Evaluation, iCASE No. 17000191, with contributions from Rolls-Royce Holdings plc and Jacobs Engineering Group Inc. The authors have no conflicts of interest to disclose. The data are available on request from the authors.