Numerical Study on Surface Roughness Measurement Based on Nonlinear Ultrasonics in Through-Transmission and Pulse-Echo Modes

Ultrasonic is one of the well-known methods for surface roughness measurement, but small roughness will only lead to a subtle variation of transmission or reflection. To explore sensitive techniques for surfaces with small roughness, nonlinear ultrasonic measurement in through-transmission and pulse-echo modes was proposed and studied based on an effective unit-cell finite element (FE) model. Higher harmonic generation in solids was realized by applying the Murnaghan hyperelastic material model. This FE model was verified by comparing the absolute value of the nonlinearity parameter with the analytical solution. Then, random surfaces with different roughness values ranging from 0 μm to 200 μm were repeatedly generated and studied in the two modes. The through-transmission mode is very suitable to measure the surfaces with roughness as small as 3% of the wavelength. The pulse-echo mode is sensitive and effective to measure the surface roughness ranging from 0.78% to 5.47% of the wavelength. This study offers a potential nondestructive testing and monitoring method for the interfaces or inner surfaces of the in-service structures.


Introduction
In practical engineering problems, the structure surface is always far from ideally flat and smooth and varies from specimen to specimen. Such an imperfect surface might affect the mechanical performance, including wear resistance, fatigue strength, and corrosion resistance of the components. Although various surface treatments such as polishing, lapping, and milling could be applied to reduce the surface roughness, it may not be easy to prepare perfect surfaces to the required level due to the high need for manpower, machinery, and economic cost. Moreover, the surface state of in-service structure probably changes with time. Therefore, as one of the key indicators of surface texture, it is of importance to evaluate surface roughness during the manufacturing and the service process.
Typical roughness measurement methods include stylus sensing, optical interferometry, and scanning tunneling microscopy. These profiling methods can provide detailed information of the surface topology. But abrasion and scratches might be introduced to the surface during stylus sensing and optical methods are applicable primarily to the roughness less than 1 µm due to the limitation of light wavelength. For the industrial components manufactured by common processes like milling, shot peening, and chemical deposition, the surface roughness ranges from about 0.5 µm to 50 µm. Moreover, most current methods are applicable only to accessible top surfaces. In a practical engineering application, the information of interfaces [1,2] or inner surfaces [3] are sometimes also demanded. To balance the resolution and penetration capability, ultrasonic method is the will demonstrate the feasibility of surface roughness measurement based on nonlinear ultrasonics. Besides, the study on ultrasonic wave scattering from rough surfaces is also beneficial for the enhancement of ultrasonic thickness measurements [3] and the detection of rough cracks or delamination in real cases [29][30][31].
In this study, we focus on the investigation of the feasibility of surface roughness measurement based on the nonlinear ultrasonic in both through-transmission and pulseecho modes. First, the basic theory for ultrasonic higher harmonic generation based on the plane wave theory and ultrasonic scattering from rough surfaces will be briefly given in Sections 2 and 3. An efficient FE model for nonlinear ultrasonic testing will be described in Section 4. Based on this verified FE model, the results will be discussed for surface roughness measurement in both through-transmission and pulse-echo modes. Finally, conclusions will be drawn in Section 6.

Higher Harmonic Generation in Solid Media
Numerous works had introduced the derivation and solution of the nonlinear ultrasonic wave equation. Following Green's approach and notations [32], the differential equation of wave motion in a form with separate linear and nonlinear terms is written as Equation (1): where ρ is the density; u is the particle displacement in the x-direction; E 1 = E and E 2 = βE 1 give the expressions of elastic constants. It should be noted that E 1 is expressed in terms of second-order elastic constants only, while E 2 is expressed in terms of both second-and third-order elastic constants. Using the perturbation technique, the displacement solution of Equation (1) can be expressed as u = u 0 + u 1 , where u 0 is the linear response; u 1 is the perturbation term; and u 1 u 0 . The trial solution, u 0 = Asin(kx − ωt) is used, where A is the incident displacement amplitude, and k and ω are the wavenumber and the angular frequency of the incident wave, respectively.
When a semi-infinite material is taken into consideration, the approximate solution of Equation (1) involving the second harmonic is given as Equation (2) [32]: where E = ρc 2 is the Young's modulus and c is the wave speed. Therefore, when a continuous plane wave is incident into the solid, the ultrasonic nonlinearity parameter β for the lossless material can be obtained by measuring the fundamental A 1 and second harmonic displacement A 2 after a propagation distance x. The resulted expression is thus derived as Equation (3): Equation (3) is widely used in nonlinear ultrasonic measurement to calculate the absolute ultrasonic nonlinearity. Once the amplitudes of fundamental A 1 and second harmonic A 2 are obtained, the nonlinearity parameter β can be calculated as a material property indicator.
On the other hand, various analytical expressions derived from the hyperelastic theory have been developed to account for the elastic contribution to the ultrasonic nonlinearity, which is related to the third-order elastic constants and is sensitive to the material microstructural variations. The well-known hyperelastic constitutive models, including the Neo-Hookean, Mooney-Rivlin, Yeoh, Ogden, and Murnaghan models [33,34], can describe the behaviors of nonlinear materials. Among them, Murnaghan's model is a popular and classical hyperelastic model used to study the wave propagation in a quadratic nonlinear  [33] as Equation (4): where λ and µ are the Lame constants; l and m are the Murnaghan third-order elastic constants. With these elastic constants, the ultrasonic nonlinearity parameter can be calculated explicitly for numerical model verification.

Ultrasonic Scattering from Rough Surfaces
Generally, rough surfaces can be categorized into simplified regular surfaces, such as sinusoidal, triangular, rectangular, and stochastically dominated random surfaces. Different roughness parameters have been defined to characterize the surfaces. The peak R p , valley R v , average R a , and total R t are defined to characterize the surface h(x), as shown in Figure 1. Moreover, commonly used root-mean-square (RMS) roughness R q are defined within a sampling length L s as Equation (5): = − 3( + 2 ) + 2( + 2 ) + 2 (4) where and are the Lame constants; and are the Murnaghan third-order elastic constants. With these elastic constants, the ultrasonic nonlinearity parameter can be calculated explicitly for numerical model verification.

Ultrasonic Scattering from Rough Surfaces
Generally, rough surfaces can be categorized into simplified regular surfaces, such as sinusoidal, triangular, rectangular, and stochastically dominated random surfaces. Different roughness parameters have been defined to characterize the surfaces. The peak , valley , average , and total are defined to characterize the surface ℎ( ), as shown in Figure 1. Moreover, commonly used root-mean-square (RMS) roughness are defined within a sampling length as Equation (5): The general statistical characteristic of the rough surface ℎ( ) can be expressed using the probability density function (PDF). For classic random cases, the central limit theorem implies that the PDF can be described as the Gaussian distribution, as expressed in Equation (6): The PDF only deal with the amplitudes variation, so another function should be applied to represent the lateral variation. An autocorrelation function is defined for a Gaussian surface ℎ( , ) as Equation (7): where ⃗ refers to the separation distance between two points; ⃗ indicates a radius vector from the origin. and refer to the correlation lengths in the -and -directions, defined as the distances where the correlation function falls to 1/e of its maximum value. For isotropic rough surfaces, = = . By controlling the RMS value and the correlation length , commonly used moving average method can be applied to generate random surfaces [13]. Examples of 2D random surfaces were generated with a The general statistical characteristic of the rough surface h(x) can be expressed using the probability density function (PDF). For classic random cases, the central limit theorem implies that the PDF can be described as the Gaussian distribution, as expressed in Equation (6): The PDF only deal with the amplitudes variation, so another function should be applied to represent the lateral variation. An autocorrelation function is defined for a Gaussian surface h(x, y) as Equation (7): where → r 0 refers to the separation distance between two points; → r indicates a radius vector from the origin. L cx and L cy refer to the correlation lengths in the xand y-directions, defined as the distances where the correlation function falls to 1/e of its maximum value. For isotropic rough surfaces, L cx = L cy = L c . By controlling the RMS value R q and the correlation length L c , commonly used moving average method can be applied to generate random surfaces [13]. Examples of 2D random surfaces were generated with a Gaussian height density function and an exponential autocorrelation function, as shown in Figure 2. Clearly, the surface profile changes significantly with R q value. What should be mentioned is that the generated surface profiles with the same parameters are different from each process due to the random nature. Therefore, numerous studies should be conducted for the same rough surface to obtain a general rule with respect to the statistical parameters R q and L c . Gaussian height density function and an exponential autocorrelation function, as shown in Figure 2. Clearly, the surface profile changes significantly with value. What should be mentioned is that the generated surface profiles with the same parameters are different from each process due to the random nature. Therefore, numerous studies should be conducted for the same rough surface to obtain a general rule with respect to the statistical parameters and . In ultrasonic nondestructive testing and monitoring, the rough surface may render the ultrasonic signals unpredictable due to the complicated wave scattering. When the slope of the surface is sufficiently small, i.e., ≪ , The wave scattering from random rough surface can be simplified into phase changes as PSA theory to deal with the problems of reflection and transmission as reported by Nagy and Rose [6]. For a normally incident longitudinal wave, the roughness modified reflection coefficient and transmission coefficient can be expressed as Equations (8) and (9) [6]: where and refer to the reflection and transmission coefficients from a smooth surface; and indicate the wavenumber in the incident and transmitted media. When nonlinear effect is considered, scaling factors In ultrasonic nondestructive testing and monitoring, the rough surface may render the ultrasonic signals unpredictable due to the complicated wave scattering. When the slope of the surface is sufficiently small, i.e., R q L c , The wave scattering from random rough surface can be simplified into phase changes as PSA theory to deal with the problems of reflection and transmission as reported by Nagy and Rose [6]. For a normally incident longitudinal wave, the roughness modified reflection coefficient and transmission coefficient can be expressed as Equations (8) and (9) [6]: where R 0 and T 0 refer to the reflection and transmission coefficients from a smooth surface; k and k T indicate the wavenumber in the incident and transmitted media. When nonlinear effect is considered, scaling factors exp −8R 2 q k 2 and exp −2R 2 q (k T − k) 2 should be applied to the amplitude of the second harmonic in the pulse-echo and through transmission modes. As the nonlinear parameter is calculated as β = 8 , a scaling factor s should be considered for nonlinearity parameter calculation with rough surfaces. In the through- transmission mode with direct contact measurement, the factor s T can be expressed as . Therefore, the roughness-induced variation in conventional ultrasonic testing is much larger than that in nonlinear ultrasonic testing. Figure 3 shows the comparison of the roughness-induced change of the fundamental and second harmonic, and the nonlinearity parameter in the through-transmission mode. As expected, the second harmonic and nonlinearity are much affected by the surface roughness. Therefore, applying nonlinear ultrasonic testing has high potential to conduct surface roughness measurement with small RMS and to bridge the gap between conventional ultrasonic method and the optical methods.
Materials 2021, 14, 4855 6 of 16 transmission modes. As the nonlinear parameter is calculated as = , a scaling factor should be considered for nonlinearity parameter calculation with rough surfaces.
In the through-transmission mode with direct contact measurement, the factor can be expressed as = − 2 . Therefore, the roughness-induced variation in conventional ultrasonic testing is much larger than that in nonlinear ultrasonic testing. Figure  3 shows the comparison of the roughness-induced change of the fundamental and second harmonic, and the nonlinearity parameter in the through-transmission mode. As expected, the second harmonic and nonlinearity are much affected by the surface roughness. Therefore, applying nonlinear ultrasonic testing has high potential to conduct surface roughness measurement with small RMS and to bridge the gap between conventional ultrasonic method and the optical methods.

Unit-Cell Model Setup
The FE methods are sufficiently powerful to calculate problems of linear and nonlinear elastic wave propagation with complicated boundaries [18,35]. In this section, the rough surfaces with Gaussian random profiles were taken into consideration to investigate the variation on roughness. Firstly, Murnaghan material model from solid mechanic module in a commercial FE software Comsol Multiphysics v5.4 (Burlington, MA, USA) was adopted to simulate the nonlinear wave propagation [33]. The specimen is Al-1200 and its material properties are shown in Table 1 [36].

Density Lame Parameters Murnaghan Third-Order Elastic Constants
To avoid the wave scattering from the irrelevant boundaries, a large-size model with non-reflecting boundaries is usually applied [12,29,35]. However, the non-reflecting boundary may also bring in unpredictable numerical errors, making the received signals difficult to analyze. Unit-cell models, calculating only a represented domain, have been proposed and proven to be effective and efficient to simulate the wave propagation in infinite periodic boundaries [28,35]. As the correlation length in random surface describes periodicity of the surface, unit-cell model was applied in this study to facilitate a series of

Unit-Cell Model Setup
The FE methods are sufficiently powerful to calculate problems of linear and nonlinear elastic wave propagation with complicated boundaries [18,35]. In this section, the rough surfaces with Gaussian random profiles were taken into consideration to investigate the variation on roughness. Firstly, Murnaghan material model from solid mechanic module in a commercial FE software Comsol Multiphysics v5.4 (Burlington, MA, USA) was adopted to simulate the nonlinear wave propagation [33]. The specimen is Al-1200 and its material properties are shown in Table 1 [36]. Table 1. Material properties of Al-1200 [36].

Density Lame Parameters Murnaghan Third-Order Elastic Constants
To avoid the wave scattering from the irrelevant boundaries, a large-size model with non-reflecting boundaries is usually applied [12,29,35]. However, the non-reflecting boundary may also bring in unpredictable numerical errors, making the received signals difficult to analyze. Unit-cell models, calculating only a represented domain, have been proposed and proven to be effective and efficient to simulate the wave propagation in infinite periodic boundaries [28,35]. As the correlation length in random surface describes periodicity of the surface, unit-cell model was applied in this study to facilitate a series of simulations for different rough surfaces. Moreover, although the 3D model is able to include the phenomena as much as possible, the 2D model with 1D surfaces has been proven to be effective to investigate the interaction between acoustic waves and rough surfaces [28][29][30].
These results have been well explained by the Kirchhoff approximation and validated by experimental results in previous literatures. In this study, a lot of simulations (50 realizations for each roughness) should be conducted for the rough surfaces. To complete all the simulation with proper accuracy at an acceptable computation cost, only 2D model is applied. Therefore, the schematic diagram of 2D FE model was shown in Figure 4a. A 10-cycle Gaussian-modulated sinusoidal wave with a central frequency of f c = 5 MHz was applied as the incident wave to suppress the sidelobes. The incident waveform with a normalized amplitude is shown in Figure 4b. In order to generate an obvious nonlinear effect, the incident pressure of relatively high amplitude should be applied, and the peak amplitude of the pressure was set as 10 MPa. Symmetric boundaries were set for upper and lower boundaries. A symmetry condition is free in the plane and fixed in the out-of-plane direction, which is defined by the equation u·n = 0. Besides, low-reflecting boundary was set at the left side to prevent the wave reflection. simulations for different rough surfaces. Moreover, although the 3D model is able to include the phenomena as much as possible, the 2D model with 1D surfaces has been proven to be effective to investigate the interaction between acoustic waves and rough surfaces [28][29][30]. These results have been well explained by the Kirchhoff approximation and validated by experimental results in previous literatures. In this study, a lot of simulations (50 realizations for each roughness) should be conducted for the rough surfaces. To complete all the simulation with proper accuracy at an acceptable computation cost, only 2D model is applied. Therefore, the schematic diagram of 2D FE model was shown in Figure 4a. A 10-cycle Gaussian-modulated sinusoidal wave with a central frequency of = 5 MHz was applied as the incident wave to suppress the sidelobes. The incident waveform with a normalized amplitude is shown in Figure 4b. In order to generate an obvious nonlinear effect, the incident pressure of relatively high amplitude should be applied, and the peak amplitude of the pressure was set as 10 MPa. Symmetric boundaries were set for upper and lower boundaries. A symmetry condition is free in the plane and fixed in the out-ofplane direction, which is defined by the equation ⋅ = 0. Besides, low-reflecting boundary was set at the left side to prevent the wave reflection.
(a) (b) As the FE methods numerically solve the partial differential equations by discretizing the problem domain into small elements, the careful choice of mesh size is critical to obtain accurate FE simulation results, especially for the nonlinear wave propagation problem. In order to investigate the variation as a function of the mesh, different sizes of quadratic elements were set to mesh the domain and get the solutions. Element number per wavelength = / was chosen as 1-30 to implement parametric study. Both the variations of and were investigated for nonlinear ultrasonic wave propagation. The amplitudes of and converge when the mesh is refined as = 20 [18]. As the longitudinal wavelength at 5 MHz is = 1.28 mm, the element size was chosen as = 0.06 mm in this study. For the irregular domain near the rough surface, triangular element should be applied, and a further mesh refinement is required in order to ensure the convergence of the model [18]. The element size was set as 0.003 mm around the surface boundary, as shown in Figure 4a. Moreover, a generalized-α method was adopted as the time-dependent solver in this FE model to solve the transient wave propagation problem. To ensure the calculation convergence, the time increment should satisfy the Courant-Friedrichs-Lewy (CFL) condition, which is defined as is the mesh size, and is the element order. For quadratic elements, a value of 0.2 for has been demonstrated to be accurate for the final results [34]. Considering the finer mesh around the rough surface, the time increment ∆ was chosen to be 1 ns in our model, which is also the sampling interval in the following data acquisition.
To balance the computation amount and the accuracy, the width of model was set as 6 mm (≈5 ). The correlation length was set as = 0.1 mm (≈ /12) to ensure that the generated periodic surfaces have the same statistical characteristics. The RMS height As the FE methods numerically solve the partial differential equations by discretizing the problem domain into small elements, the careful choice of mesh size is critical to obtain accurate FE simulation results, especially for the nonlinear wave propagation problem. In order to investigate the variation as a function of the mesh, different sizes ∆x of quadratic elements were set to mesh the domain and get the solutions. Element number per wavelength N x = λ w /∆x was chosen as 1-30 to implement parametric study. Both the variations of A 1 and A 2 were investigated for nonlinear ultrasonic wave propagation. The amplitudes of A 1 and A 2 converge when the mesh is refined as N x = 20 [18]. As the longitudinal wavelength at 5 MHz is λ w = 1.28 mm, the element size was chosen as ∆x = 0.06 mm in this study. For the irregular domain near the rough surface, triangular element should be applied, and a further mesh refinement is required in order to ensure the convergence of the model [18]. The element size was set as 0.003 mm around the surface boundary, as shown in Figure 4a. Moreover, a generalized-α method was adopted as the time-dependent solver in this FE model to solve the transient wave propagation problem. To ensure the calculation convergence, the time increment should satisfy the Courant-Friedrichs-Lewy (CFL) condition, which is defined as CFL m = c∆t/∆x, where ∆t is the time step, ∆x is the mesh size, and m is the element order. For quadratic elements, a value of 0.2 for CFL 2 has been demonstrated to be accurate for the final results [34]. Considering the finer mesh around the rough surface, the time increment ∆t was chosen to be 1 ns in our model, which is also the sampling interval in the following data acquisition.
To balance the computation amount and the accuracy, the width of model was set as 6 mm (≈5 λ). The correlation length was set as L c = 0.1 mm ( ≈ λ w /12) to ensure that the generated periodic surfaces have the same statistical characteristics. The RMS height R q varying from 0 to 200 µm (≈ λ w /6) were applied to investigate the influence of roughness on nonlinearity parameters. The random rough surfaces were generated by moving average method as described in Section 3. The total time was set as 9 µs to record both transmission and reflected signals. The typical results of simulated wave propagation at different times were shown in Figure 5. Besides the major energy of longitudinal wave reflecting from the rough surface, a small energy of shear wave can also be seen due to mode conversion at the surface. varying from 0 to 200 μm (≈ /6) were applied to investigate the influence of roughn on nonlinearity parameters. The random rough surfaces were generated by moving av age method as described in Section 3. The total time was set as 9 μs to record both tra mission and reflected signals. The typical results of simulated wave propagation at dif ent times were shown in Figure 5. Besides the major energy of longitudinal wave refle ing from the rough surface, a small energy of shear wave can also be seen due to mo conversion at the surface.

Model Verification
To demonstrate the feasibility of this model for nonlinear ultrasonics, the nonline ity parameter should be compared with the analytical results. Therefore, a smooth a flat surface was studied at first. When = 0, the horizontal displacement at the differ positions were recorded, as shown in Figure 6a. Both transmitted and reflected sign were obtained and two components are separated when the receiving position is aw from the surface. Otherwise, the transmitted and reflected waves will be overlapped, m ing it difficult to separate the corresponding higher harmonics. Therefore, only 0~14 mm were applied for later analysis. The results of fast Fourier transform (FFT the transmitted waves are shown in Figure 6b. Higher harmonics, especially the seco harmonic, are clearly generated in the spectrum. With the increase of propagation d tance, the amplitude of the second harmonic increases, while the fundamental h monic remains unchanged, as expected from Equation (2).

Model Verification
To demonstrate the feasibility of this model for nonlinear ultrasonics, the nonlinearity parameter should be compared with the analytical results. Therefore, a smooth and flat surface was studied at first. When R q = 0, the horizontal displacement at the different positions were recorded, as shown in Figure 6a. Both transmitted and reflected signals were obtained and two components are separated when the receiving position is away from the surface. Otherwise, the transmitted and reflected waves will be overlapped, making it difficult to separate the corresponding higher harmonics. Therefore, only x = 0 ∼ 14 mm were applied for later analysis. The results of fast Fourier transform (FFT) of the transmitted waves are shown in Figure 6b. Higher harmonics, especially the second harmonic, are clearly generated in the spectrum. With the increase of propagation distance, the amplitude of the second harmonic A 2 increases, while the fundamental harmonic A 1 remains unchanged, as expected from Equation (2). Both the absolute amplitude of second harmonic for reflected and transmitted waves were calculated for different propagation distance. The results are shown in Figure 7. For the transmitted wave, the amplitude of the second harmonic increases linearly with the propagation distance, which is the same as Equation (3). On the other hand, the amplitude of the reflected wave decreases with the propagation distance with the same change rate. When the reflected wave propagates back to the source position, the second harmonic decreases to zero. This similar phenomenon has been reported in [23,24]. To evaluate the Both the absolute amplitude of second harmonic for reflected and transmitted waves were calculated for different propagation distance. The results are shown in Figure 7. For the transmitted wave, the amplitude of the second harmonic increases linearly with the propagation distance, which is the same as Equation (3). On the other hand, the amplitude of the reflected wave decreases with the propagation distance with the same change rate. When the reflected wave propagates back to the source position, the second harmonic decreases to zero. This similar phenomenon has been reported in [23,24]. To evaluate the accuracy of this model, the calculated nonlinearity parameter β is compared with the analytical result using Equation (3). The analytical value with Murnaghan constants is β = 16.84. The calculated value using Equation (4) is β = 17.48. Although there is little difference between both values (relative error is 3.8%) due to the errors from numerical calculation and the input signal with finite cycles, the nonlinearity parameters calculated from FE model and analytical solution are well agreed. Therefore, this FE model is accurate enough for modeling higher harmonic generation in a nonlinear elastic material and will be applied for rough surface evaluation. Both the absolute amplitude of second harmonic for reflected and transmitted waves were calculated for different propagation distance. The results are shown in Figure 7. For the transmitted wave, the amplitude of the second harmonic increases linearly with the propagation distance, which is the same as Equation (3). On the other hand, the amplitude of the reflected wave decreases with the propagation distance with the same change rate. When the reflected wave propagates back to the source position, the second harmonic decreases to zero. This similar phenomenon has been reported in [23,24]. To evaluate the accuracy of this model, the calculated nonlinearity parameter is compared with the analytical result using Equation (3). The analytical value with Murnaghan constants is = 16.84. The calculated value using Equation (4) is = 17.48. Although there is little difference between both values (relative error is 3.8%) due to the errors from numerical calculation and the input signal with finite cycles, the nonlinearity parameters calculated from FE model and analytical solution are well agreed. Therefore, this FE model is accurate enough for modeling higher harmonic generation in a nonlinear elastic material and will be applied for rough surface evaluation.

Through-Transmission Mode
To simulate the nonlinear ultrasonic testing in the through-transmission mode, a perfect coupled water film was added at the right boundary to mimic the couplant. To reduce

Through-Transmission Mode
To simulate the nonlinear ultrasonic testing in the through-transmission mode, a perfect coupled water film was added at the right boundary to mimic the couplant. To reduce the computation amount, the length of the solid material was reduced to 5 mm. The total time was set as 4 µs to record the transmission signals. Different rough surfaces were generated with RMS roughness R q ranging from 0 to 200 µm. A line average in the water layer was recorded as the transmitted signal and three typical signals for the same roughness R q = 50 µm as shown in Figure 8a. Compared to the smooth surface, a ringdown effect occurs for rough surfaces, making the received waveforms more complicated. Moreover, there are remarkable differences in the amplitude among the three signals due to the mutual interference between components of the signal transmitted through different parts of the rough surface. Figure 8b shows the corresponding spectra. Both fundamental and second harmonics can be clearly identified, while the amplitudes slightly vary among each realization. It demonstrates that rough surface will bring in some noises both in the linear and nonlinear ultrasonic measurement. effect occurs for rough surfaces, making the received waveforms more complicated. Moreover, there are remarkable differences in the amplitude among the three signals due to the mutual interference between components of the signal transmitted through different parts of the rough surface. Figure 8b shows the corresponding spectra. Both fundamental and second harmonics can be clearly identified, while the amplitudes slightly vary among each realization. It demonstrates that rough surface will bring in some noises both in the linear and nonlinear ultrasonic measurement. the number is few, the coherent amplitudes are not stable for both fundamental and second harmonics. When the number increases to more than 30, the coherent amplitudes of the fundamental wave are close to convergence, with a relatively small variation less than 5%. However, the coherent amplitudes of the second harmonic converge at a larger number around 50. The reason lies in that the second harmonic contains the information of the interaction between waves with smaller wavelength and rough surfaces. Therefore, to ensure the accuracy of the coherent wave for both fundamental and second harmonics, 50 In order to investigate the influence of surface profiles on the amplitude of transmitted wave, different realizations were conducted for the rough surface with the same roughness R q = 200 µm. The normalized amplitudes of fundamental and second harmonics are shown in Figure 9a,b. Due to the random process of the rough surface realization, both amplitudes vary from each realization. For normalized fundamental amplitude, the arithmetic mean of the 50 realizations is 0.4730 and the standard deviation is 0.1867. For the second harmonic, the mean value and standard deviation are 0.3508 and 0.1866, respectively. Compared to the fundamental harmonic, the amplitudes of the second harmonic decrease much more. The fluctuations of both amplitudes are very large; therefore, sufficient realizations of rough surfaces for the same roughness are required to obtain the coherent signal, calculated as the average of multiple simulations. To fully remove the out-of-phase components from the scattered signals, a proper number of simulations should be determined to balance the accuracy and computation amount. The influence of the increasing number of realizations on the amplitude of coherent signal was investigated. Figure 9c,d show the variations of coherent fundamental and second harmonic wave for the roughness R q = 200 µm, which is the maximum value in this study. When the number is few, the coherent amplitudes are not stable for both fundamental and second harmonics. When the number increases to more than 30, the coherent amplitudes of the fundamental wave are close to convergence, with a relatively small variation less than 5%. However, the coherent amplitudes of the second harmonic converge at a larger number around 50. The reason lies in that the second harmonic contains the information of the interaction between waves with smaller wavelength and rough surfaces. Therefore, to ensure the accuracy of the coherent wave for both fundamental and second harmonics, 50 simulations will be conducted repeatedly for each surface with the same roughness R q ranging from 0 µm to 200 µm. The normalized amplitude variations as a function of the ratio of roughness to the incident wavelength are shown in Figure 10. When the roughness increases, both fundamental and second harmonic amplitudes vary for each realization. The average of 50 simulations was calculated as the coherent amplitude for each roughness to compare with the analytical solutions from the PSA theory. For the normalized fundamental amplitude, i.e., transmission coefficient in this case, the coherent amplitudes from FE simulation generally agree with the PSA results when the roughness is small. However, an obvious deviation appears when the roughness is larger than 0.07 and increases further with the increasing roughness value . It means that PSA fails to predict the transmission coefficient when the roughness is large, which is consistent with the result in the previous literature [2]. Therefore, as the roughness increases, Kirchhoff based analytical solution becomes increasingly inaccurate and a fully numerical approach is required. The normalized amplitude variations as a function of the ratio of roughness to the incident wavelength are shown in Figure 10. When the roughness increases, both fundamental and second harmonic amplitudes vary for each realization. The average of 50 simulations was calculated as the coherent amplitude for each roughness to compare with the analytical solutions from the PSA theory. For the normalized fundamental amplitude, i.e., transmission coefficient in this case, the coherent amplitudes from FE simulation generally agree with the PSA results when the roughness is small. However, an obvious deviation appears when the roughness is larger than 0.07 λ w and increases further with the increasing roughness value R q . It means that PSA fails to predict the transmission coefficient when the roughness is large, which is consistent with the result in the previous literature [2]. Therefore, as the roughness increases, Kirchhoff based analytical solution becomes increasingly inaccurate and a fully numerical approach is required. For the second harmonic amplitude, as shown in Figure 10b, the FE results also show a generally decreasing trend with increasing roughness and starts to deviate from the PSA result at a smaller roughness value of around 0.03 . Moreover, the fluctuation of nonlinearity parameter is much larger than that of transmission coefficient due to the existence of the second harmonic, which is half the wavelength of the incident wave. Therefore, surface roughness will induce much more noises in the nonlinear ultrasonic measurement, which is the reason why careful surface treatment must be applied to the specimen prior to the measurement. Meanwhile, the second harmonic changes with the roughness at a higher rate than transmission coefficient, as predicted in Figure 3. When a roughness is as small as 20 μm corresponding to 1.56% , only a small change of 0.52 dB will be introduced to the transmission coefficient. Such subtle variation brings in difficulty to accurately measure this small roughness by conventional ultrasonic method. For the same roughness, the amplitude of second harmonic decreases to −1.70 dB. As the numerical results agree well with the PSA result when the roughness is smaller than 0.03 , PSA solution can be applied as the calibration curve for uncertainty analysis. Thus, when the second harmonic is applied for the measurement of nominal roughness for = 20 μm, the standard uncertainty is calculated as = ( ) ∑ ( , ) = 0.45 μm with = 50 observations. It indicates that the second harmonic generation is much more sensitive to the surface roughness. Therefore, nonlinear ultrasonic is very suitable to characterize the surface roughness, especially for the surfaces with small RMS.

Pulse-Echo Mode
The nonlinear ultrasonic testing in the pulse-echo mode was also built and the model is the same as described in Section 4.1. A line average near the excitation position was recorded and three typical signals for the same roughness = 50 μm as shown in Figure 11a. Compared to the results from the smooth surface, a ring-down effect is also observed for rough surfaces. Similarly, there are obvious differences in the amplitude among the three realizations due to mutual interference from different parts of the rough surface. Both fundamental and second harmonics of the reflected wave vary from each realization, as shown in the corresponding spectra in Figure 11b. As illustrated in Figure 7, the second harmonic almost decreases to zero when the receiving position is close to the excitation. However, the second harmonics are clearly observed when it reflected from the rough surface. It means that the second harmonic can also be received For the second harmonic amplitude, as shown in Figure 10b, the FE results also show a generally decreasing trend with increasing roughness and starts to deviate from the PSA result at a smaller roughness value of around 0.03 λ w . Moreover, the fluctuation of nonlinearity parameter is much larger than that of transmission coefficient due to the existence of the second harmonic, which is half the wavelength of the incident wave. Therefore, surface roughness will induce much more noises in the nonlinear ultrasonic measurement, which is the reason why careful surface treatment must be applied to the specimen prior to the measurement. Meanwhile, the second harmonic changes with the roughness at a higher rate than transmission coefficient, as predicted in Figure 3. When a roughness R q is as small as 20 µm corresponding to 1.56% λ w , only a small change of 0.52 dB will be introduced to the transmission coefficient. Such subtle variation brings in difficulty to accurately measure this small roughness by conventional ultrasonic method. For the same roughness, the amplitude of second harmonic decreases to −1.70 dB. As the numerical results agree well with the PSA result when the roughness is smaller than 0.03 λ w , PSA solution can be applied as the calibration curve for uncertainty analysis. Thus, when the second harmonic is applied for the measurement of nominal roughness for R q = 20 µm, the standard uncertainty is calculated as s Rq = 1 n(n−1) ∑ i=n i=1 R q,i − R q 2 = 0.45 µm with n = 50 observations. It indicates that the second harmonic generation is much more sensitive to the surface roughness. Therefore, nonlinear ultrasonic is very suitable to characterize the surface roughness, especially for the surfaces with small RMS.

Pulse-Echo Mode
The nonlinear ultrasonic testing in the pulse-echo mode was also built and the model is the same as described in Section 4.1. A line average near the excitation position was recorded and three typical signals for the same roughness R q = 50 µm as shown in Figure 11a. Compared to the results from the smooth surface, a ring-down effect is also observed for rough surfaces. Similarly, there are obvious differences in the amplitude among the three realizations due to mutual interference from different parts of the rough surface. Both fundamental and second harmonics of the reflected wave vary from each realization, as shown in the corresponding spectra in Figure 11b. As illustrated in Figure 7, the second harmonic almost decreases to zero when the receiving position is close to the excitation. However, the second harmonics are clearly observed when it reflected from the rough surface. It means that the second harmonic can also be received with planar transducer rather than focused transducer or dual-element transducer as reported in the previous literatures [24][25][26][27]. As the accumulated second harmonic is zero for smooth surface, these remaining second harmonics only result from the roughness. Therefore, there is potential in measuring the surface roughness in the pulse-echo mode based on second harmonic generation. with planar transducer rather than focused transducer or dual-element transducer as reported in the previous literatures [24][25][26][27]. As the accumulated second harmonic is zero for smooth surface, these remaining second harmonics only result from the roughness. Therefore, there is potential in measuring the surface roughness in the pulse-echo mode based on second harmonic generation.
(a) (b) Similarly, rough surfaces with roughness ranging from 0 to 200 μm were investigated in the pulse-echo mode and 50 simulations were conducted repeatedly for each roughness. The normalized amplitude variations as a function of the ratio of roughness to the incident wavelength are shown in Figure 12. For the normalized fundamental amplitude, i.e., reflection coefficient in the pulse-echo mode, the coherent amplitudes from FE simulation generally agree with the PSA results when the roughness are smaller than 0.007 . However, an obvious deviation increases with the increasing roughness value . Also, PSA fails to predict the reflection coefficient when the roughness is large. On the other hand, the second harmonic amplitude is close to zero when the roughness is small, which is like the reflection from smooth surface as shown in Figure 12b. When the roughness increases to 10 μm corresponding to 0.78% , the second harmonic amplitudes begin to increase quickly with the roughness. When the roughness is too large, the second harmonic amplitudes of reflected waves keep almost constant. When = 70 μm, the averaged amplitudes increase to 17.33 dB. In contrast, the reflection coefficient decreases to −1.85 dB when = 70 μm. It indicates that the second harmonic generation in the pulse-echo mode is very sensitive to measure the surface roughness within the range of = 0.78~5.47% . Similarly, rough surfaces with roughness R q ranging from 0 to 200 µm were investigated in the pulse-echo mode and 50 simulations were conducted repeatedly for each roughness. The normalized amplitude variations as a function of the ratio of roughness to the incident wavelength are shown in Figure 12. For the normalized fundamental amplitude, i.e., reflection coefficient in the pulse-echo mode, the coherent amplitudes from FE simulation generally agree with the PSA results when the roughness are smaller than 0.007 λ w . However, an obvious deviation increases with the increasing roughness value R q . Also, PSA fails to predict the reflection coefficient when the roughness is large. On the other hand, the second harmonic amplitude is close to zero when the roughness is small, which is like the reflection from smooth surface as shown in Figure 12b. When the roughness R q increases to 10 µm corresponding to 0.78% λ w , the second harmonic amplitudes begin to increase quickly with the roughness. When the roughness is too large, the second harmonic amplitudes of reflected waves keep almost constant. When R q = 70 µm, the averaged amplitudes increase to 17.33 dB. In contrast, the reflection coefficient decreases to −1.85 dB when R q = 70 µm. It indicates that the second harmonic generation in the pulse-echo mode is very sensitive to measure the surface roughness within the range of R q = 0.78 ∼ 5.47% λ w . Therefore, second harmonic generation can be applied to measure the surface roughness, not only in the through-transmission mode but also in the pulse-echo mode. Compared to the transmission coefficient or reflection coefficient, the second harmonic generation is more remarkably influenced by the rough surface. Therefore, surface roughness with relatively small RMS can be measured by nonlinear ultrasonics with a higher sensitivity. In the through-transmission mode, the surface roughness within the range of R q = 0 ∼ 3% λ w can be measured. The effective range of surface roughness is R q = 1 ∼ 5% λ w in the pulseecho mode. Therefore, second harmonic generation can be applied to measure the surface roughness, not only in the through-transmission mode but also in the pulse-echo mode. Compared to the transmission coefficient or reflection coefficient, the second harmonic generation is more remarkably influenced by the rough surface. Therefore, surface roughness with relatively small RMS can be measured by nonlinear ultrasonics with a higher sensitivity. In the through-transmission mode, the surface roughness within the range of = 0~3% can be measured. The effective range of surface roughness is = 1~5% in the pulse-echo mode.

Conclusions
In this study, a numerical study has been carried out to investigate the surface roughness measurement based on nonlinear ultrasonic method. The Murnaghan hyperelastic material model was adopted in a unit-cell FE model to simulate higher harmonic generation in solids. The absolute value of the nonlinearity parameter was recovered, agreeing well with the analytical solution. Therefore, this FE model provides an effective way to analyze the higher harmonic generation in solids. Based on the verified FE model, the ultrasonic wave scattering with rough surfaces in both through-transmission and pulse-echo modes were investigated. Rough surfaces with Gaussian random profiles were generated by the moving average method and the RMS height varying from 0 to 200 μm (≈ /6) were studied. 50 simulations were conducted repeatedly for each roughness to remove the out-of-phase components.
In the through-transmission mode, the transmission coefficient from FE simulation generally agrees with the PSA prediction when the roughness is smaller than 0.07 , while the second harmonic deviates at a smaller roughness value of around 0.03 . The higher decreasing rate of the second harmonic as a function of the roughness demonstrated that the nonlinear ultrasonic testing is very suitable to measure the surface roughness, especially for the surfaces with RMS as small as 0.03 . In the pulse-echo mode, the second harmonic is clearly observed when it is reflected from the rough surface with RMS larger than 0.007 and its amplitude increases quickly with the roughness. As the accumulated second harmonic is zero reflected from smooth surface, the second harmonic generation in the pulse-echo mode is very sensitive to measure the surface roughness within the range of = 0.78~5.47% .
The numerical model introduced in this paper offers an efficient way to study the interaction of nonlinear ultrasonic wave with rough surfaces. The results can provide the correction terms for the absolute nonlinear ultrasonic measurement. The feasibility of the surface roughness with small RMS in the pulse-echo mode offers a potential early-stage monitoring method

Conclusions
In this study, a numerical study has been carried out to investigate the surface roughness measurement based on nonlinear ultrasonic method. The Murnaghan hyperelastic material model was adopted in a unit-cell FE model to simulate higher harmonic generation in solids. The absolute value of the nonlinearity parameter was recovered, agreeing well with the analytical solution. Therefore, this FE model provides an effective way to analyze the higher harmonic generation in solids. Based on the verified FE model, the ultrasonic wave scattering with rough surfaces in both through-transmission and pulse-echo modes were investigated. Rough surfaces with Gaussian random profiles were generated by the moving average method and the RMS height R q varying from 0 to 200 µm (≈ λ w /6) were studied. 50 simulations were conducted repeatedly for each roughness to remove the out-of-phase components.
In the through-transmission mode, the transmission coefficient from FE simulation generally agrees with the PSA prediction when the roughness is smaller than 0.07 λ w , while the second harmonic deviates at a smaller roughness value of around 0.03 λ w . The higher decreasing rate of the second harmonic as a function of the roughness demonstrated that the nonlinear ultrasonic testing is very suitable to measure the surface roughness, especially for the surfaces with RMS as small as 0.03 λ w . In the pulse-echo mode, the second harmonic is clearly observed when it is reflected from the rough surface with RMS larger than 0.007 λ w and its amplitude increases quickly with the roughness. As the accumulated second harmonic is zero reflected from smooth surface, the second harmonic generation in the pulse-echo mode is very sensitive to measure the surface roughness within the range of R q = 0.78 ∼ 5.47% λ w .
The numerical model introduced in this paper offers an efficient way to study the interaction of nonlinear ultrasonic wave with rough surfaces. The results can provide the correction terms for the absolute nonlinear ultrasonic measurement. The feasibility of the surface roughness with small RMS in the pulse-echo mode offers a potential early-stage monitoring method for the interfaces or inner surfaces of the in-service structures. In the future, work will be extended on a more realistic 3D model with 2D surface, and the experimental validation based on the dual-frequency ultrasonic transducer.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.