Linear and Nonlinear Normal Interface Stiffness in Dry Rough Surface Contact Measured Using Longitudinal Ultrasonic Waves

: When two rough surfaces are loaded together contact occurs at asperity peaks. An interface of solid contact regions and air gaps is formed that is less stiff than the bulk material. The stiffness of a structure thus depends on the interface conditions; this is particularly critical when high stiffness is required, for example in precision systems such as machine tool spindles. The rough surface interface can be modelled as a distributed spring. For small deformation, the spring can be assumed to be linear; whilst for large deformations the spring gets stiffer as the amount of solid contact increases. One method to measure the spring stiffness, both the linear and nonlinear aspect, is by the reﬂection of ultrasound. An ultrasonic wave causes a perturbation of the contact and the reﬂection depends on the stiffness of the interface. In most conventional applications, the ultrasonic wave is low power, deformation is small and entirely elastic, and the linear stiffness is measured. However, if a high-powered ultrasonic wave is used, this changes the geometry of the contact and induces nonlinear response. In previous studies through transmission methods were used to measure the nonlinear interfacial stiffness. This approach is inconvenient for the study of machine elements where only one side of the interface is accessible. In this study a reﬂection method is undertaken, and the results are compared to existing experimental work with through transmission. The variation of both linear and nonlinear interfacial stiffnesses was measured as the nominal contact pressure was increased. In both cases interfacial stiffness was expressed as nonlinear differential equations and solved to deduce the contact pressure-relative surface approach relationships. The relationships derived from linear and nonlinear measurements were similar, indicating the validity of the presented methods. Author Contributions: Conceptualization, S.T. and R.S.D.-J.; methodology, S.T. and R.S.D.-J.; soft-ware, S.T.; validation, S.T.; formal analysis, S.T. and R.S.D.-J.; investigation, S.T. and R.S.D.-J.; data curation, S.T.; writing—original draft preparation, S.T.; writing—review and editing, R.S.D.-J.; visual-ization, S.T.; supervision, R.S.D.-J.


Introduction
When two solid surfaces come into contact, it is the asperity peaks that touch, and the real area of contact is substantially less than the nominal contact area [1]. Applying a pressure to the surfaces causes a small approach of the mean lines of their roughness. The interfacial stiffness (per unit area) is then defined as the rate of change of contact pressure with approach of the mean lines of the roughness of the contacting surfaces [2]. If the real contact area is low, then a low nominal contact pressure is all that is required to deflect the asperities. Increasing the nominal contact pressure brings more asperities into contact and the interface becomes stiffer; hence interface stiffness is a nonlinear variable. Interfacial stiffness is important in tribology because it affects machine element deflection, wear, and friction, since stiffer interfaces tend to less deflection. Machining accuracy, for example, depends on the stiffness of the joints in the machine tool assembly [2].
Several studies have focused on the interfacial stiffness of rough contact from both theoretical and experimental perspectives. An early example from Mindlin [3] provided a theoretical model for elastic contact of a smooth curved surface. Later, Greenwood and Williamson (GW model) [4] presented a statistical model for the relation between the contact area of the rough elastic mating surfaces and applied pressure. The model is limited by two assumptions: firstly asperities are treated as a distribution of hemispherical capped peaks, and secondly those peaks are treated as independent such that the deformation of one peak does not affect its neighbors. The relationship between the contact area and applied pressure also presented by models such as the Whitehouse-Archard-Onion model (WAO model) [5,6] and the Buh-Gibson-Thomas model (BGT model) [7]. These models have been widely used to predict contact stiffness. Campana et al. [8] used Mindlin's theory to measure the stiffness ratio of stationary contact.
Digital image correlation (DIC) [9][10][11] has been used to measure stiffness in joints. In this method, a series of images prior to applying load and during loading are taken from the interface and surrounding area. The comparison of the images shows the relative displacement. A plot of the applied load measured with a load cell against the relative surface displacement is created; the gradient gives the contact interface. Although it has advantages such as a full field capability and high resolution, the main limitation is that a surface near the contacting surfaces must be accessible for the camera to capture images. In many items of engineering machinery optical access is limited and this limits application of DIC outside a laboratory setting. Another method is to make direct measurements of the displacement of the interface under applied load. Burdekin et al. [12] used relative displacement transducers between the fixed frame and the stack of rings to measure direct displacement of the interface. Chikate et al. [13] also directly measured displacement of a smooth surface using strain gauges and spring-loaded pins. For example, Handzel-Powierza et al. [14] used a tenso-metric bridge to measure the applied load and an indicator sensor to measure the surface approach. They compared their results to theoretical predictions from the Greenwood-Williamson model in the range of elastic deformation for quasi-isotropic surfaces. However, this technique is restricted to low contact stiffness (1 GPa/µm) [15].
It is also possible to measure contact stiffness from a structure's resonance frequency. This contact frequency method uses a relationship between the contact stiffness and the natural frequency of the mating surfaces [16]. One of the disadvantages of this method is the need for vibration response equations of the experimental tools [17].
The approach used here, also a vibrational method, is through the reflection or transmission of ultrasound [11,[18][19][20][21][22][23]. This has the advantage being non-invasive, nondestructive, and can measure interface inside a machine where desired surfaces are inaccessible. Reflection and transmission coefficients are defined as the proportion of an incident wave reflected from or transmitted through the interface, respectively. The reflection coefficient depends on the amount of solid contact. The acoustic impedance of air is significantly small, compared to solid materials, and so the air gap at the interface causes more incident energy to reflect. Increasing the nominal contact pressure leads to an increase in the number of the asperities in contact, therefore less incident energy is reflected and more energy transmitted [18]. Kendal and Tabor [24] and Tattersall [19] used a spring model, where the interface is treated as a distributed stiffness, to demonstrate that the reflection coefficient of an imperfect contact (incomplete contact) is dependent on the interfacial stiffness.
Krolikowski et al. [15,28] presented a spring-damper model in a parallel configuration to address the effect of ultrasound attenuation. However, the results were still far beyond those predicted by the GW model. Rokhlin et al. [30] modeled the interface as a layer of equivalent materials such as viscoelastic materials. Their study is accurate in cases where the Poisson's ratio of the layer model is known. The principle has been modified by Du. et al. [31], Xiao et al. [32] and Sun et al. [17].
These studies have been based around the linear response of an interface to an ultrasonic wave. The pressures generated by the wave and the resulting deflections are normally very small and so the process is elastic and reversible. However, in recent years there have been a number of studies based around higher power ultrasound that causes nonlinear behaviour. The nonlinear behaviour may originate from the material itself, as the stress-strain relationship of the medium exceeds Hooke's law (nonlinear elasticity) [33,34]. Alternatively, the interface may also be a source of nonlinearity. If the wave is of sufficient power, it can cause opening and closing of the air gaps and this can lead to the generation of higher harmonics in the reflected/transmitted waves. This is known as Contact Acoustic Nonlinearity CAN [35][36][37]. In Section 2, a description of the higher harmonics generation (mechanics of CAN) is described.
Richardson [38] studied one-dimensional nonlinear wave propagation through identical contacting materials by considering the concept of CAN. He defined an analytical relation between the incident wave and the relative surface approach of the interface. Biwa et al. [39] used Richardson's findings to propose a nonlinear nominal contact pressurerelative surface approach relationship using a polynomial Taylor series to define the linear and nonlinear interfacial stiffnesses using the reflection coefficient. Nonlinear interfacial stiffness is important when a high-power ultrasound is employed to define the interface. The wave itself causes a distortion of the interface and hence a linear model can no longer describe the proper interface response. In previous studies only through transmission methods had been used to measure the nonlinear interfacial stiffness. The aim of this paper is to evaluate linear and nonlinear interfacial stiffnesses measured with a reflection method to demonstrate that they are both based on the same fundamental relationship between contact pressure and relative surface approach. The result is compared to published experimental work using through transmission.
Both the linear and nonlinear interfacial stiffness have been deduced using a reflection coefficient R and a second order nonlinear parameter for reflected ultrasound from the interface γ. A fourth-order polynomial expression for the linear and nonlinear interfacial stiffness variation with nominal contact pressure based on experimental data is used. Since the accuracy of the nonlinear stiffness measured with ultrasound has not been studied, a first order nonlinear homogenous differential equation of the linear interfacial stiffness and a second order nonlinear homogenous differential equation of the nonlinear interfacial stiffness have been used to derive the relationship between nominal contact pressure and approach of the surface. The study shows the nominal contact pressure-relative surface approach measured from the linear and nonlinear interfacial stiffness both give similar results. This indicates the equations proposed by [39] for the linear and nonlinear interfacial stiffness can accurately measure the interfacial stiffness.

Governing Equation
Consider two linear elastic media in contact subjected to a one-dimensional elastic longitudinal ultrasonic wave ( Figure 1).
The direction of ultrasound propagation is through the thickness of the media which is defined as the x-axis. The origin of the x-axis is the point of first contact of the surfaces. The roughness mean lines above and below the interface are located at X − and X +, respectively ( Figure 1a). The separation (mean gap) between the two surfaces under zero load is then given by (X + − X − ). When there is no elastic wave the surfaces are in static equilibrium under a nominal contact pressure p 0 causing an equilibrium separation of h 0 (Figure 1b,c). When an elastic wave passes through the interface the surface separation, h(t) is cyclically reduced and then increased about the equilibrium value h 0 . The relative surface approach Y(t) is defined as the difference between the separation variation with time h(t) and equilibrium separation h 0 , (Y(t) = h(t) − h 0 ). In the absence of an elastic wave, the relative surface approach Y(t) is zero. The sign of relative surface approach Y(t) is opposite to the direction of loading/unloading by the wave. During the loading process by the elastic wave, the relative surface approach Y(t) is negative, whereas during the unloading process it is positive. In order to clearly illustrate the relative approach of the surfaces (assuming the lower body is fixed and the upper moving), the upper body is considered to be rigid and flat whilst the lower body is elastic and rough (Figure 1c). The direction of ultrasound propagation is through the thickness of the media which is defined as the -axis. The origin of the -axis is the point of first contact of the surfaces. The roughness mean lines above and below the interface are located at and , respectively ( Figure 1a). The separation (mean gap) between the two surfaces under zero load is then given by ( − ). When there is no elastic wave the surfaces are in static equilibrium under a nominal contact pressure causing an equilibrium separation of ℎ (Figure 1b,c). When an elastic wave passes through the interface the surface separation, ℎ( ) is cyclically reduced and then increased about the equilibrium value ℎ . The relative surface approach ( ) is defined as the difference between the separation variation with time ℎ( ) and equilibrium separation ℎ , ( ( ) = ℎ( ) − ℎ ). In the absence of an elastic wave, the relative surface approach ( ) is zero. The sign of relative surface approach ( ) is opposite to the direction of loading/unloading by the wave. During the loading process by the elastic wave, the relative surface approach ( ) is negative, whereas during the unloading process it is positive. In order to clearly illustrate the relative approach of the surfaces (assuming the lower body is fixed and the upper moving), the upper body is considered to be rigid and flat whilst the lower body is elastic and rough ( Figure 1c). The elastic wave pressure superimposed on the nominal contact pressure generates a relative surface approach ( ) which varies with time (shown in Figure 1d). Increasing nominal contact pressure decreases separation of the surfaces. The equations described in this section follow the method presented in [38,39]. The equation of the interface subjected to a one-dimensional longitudinal ultrasonic wave is given by: The stress-displacement relation is given by: The elastic wave pressure superimposed on the nominal contact pressure p 0 generates a relative surface approach Y(t) which varies with time (shown in Figure 1d). Increasing nominal contact pressure p 0 decreases separation of the surfaces.
The equations described in this section follow the method presented in [38,39]. The equation of the interface subjected to a one-dimensional longitudinal ultrasonic wave is given by: The stress-displacement relation is given by: where ρ is the density of the media, σ(x, t) is the normal stress generated by the ultrasonic wave and p 0 in the media, u(x, t) is the displacement of the ultrasonic wave and p 0 is the nominal contact pressure. The general solution of Equation (1) results in the displacement of each side of the interface: where u(x, t) is the displacement of the interface, c is the speed of sound in the media which is E/ρ, t is the time of flight of ultrasound, f (x − ct) is the incident wave, G(x + ct) and H(x − ct) are the reflected and transmitted wave, respectively. The displace-ment through the upper and lower media is as x < X − and x > X + , respectively. The boundary condition at the interface is defined by: where h(t) is the rough surface mean line separation (surface separation) variation with time and h 0 is the separation under static equilibrium p 0 in the absence of the elastic wave (as shown in Figure 1). A negative sign pressure term indicates compression at the interface. Substituting Equation (4) into Equation (2) gives: The relative surface approach Y(t) of the mating surfaces is defined as the difference in the wave displacement u(x, t) of both sides of the interface ( Figure 1d): Equation (6) shows that the relative surface approach Y(t) can be used instead of the surface separation h(t). The translational motion of the interface is defined by: In order to simplify the calculation, the displacement of the interface is represented only in terms of incident wave [38]: Differentiation of Equations (6) and (8) with respect to time gives: where dots denote the differentiation with respect to time.
Equation (10) is expressed in terms of incident, reflected, and transmitted ultrasonic waves; it is necessary to reduce this in terms of only the incident wave and nominal contact pressure. To do this Equation (3) is differentiated with respect to x: Rearranging Equation (11) to make the reflected and transmitted waves the subject gives: Substituting Equation (5) into Equation (12) gives: Appl. Sci. 2021, 11, 5720 6 of 20 Substituting Equation (13) into Equation (10) gives: Equation (14) is a first order differential equation describing the relative surface approach in terms of the incident wave and the nominal contact pressure p(h 0 + Y). Solving this differential equation results in an equation for the relative separation of the interface. Equation (14) can be written in terms of surface separation h(t): The reflected ultrasonic wave is now defined in terms of relative surface approach. Equations (5) and (6) give the relationship between relative surface approach Y(t) and translational motion of the contact interface X: Substituting Equation (3) (for x < X − ) and Equation (7) into Equation (15) gives: It can be seen from Equation (17) that the reflected ultrasonic wave G(x + ct) is expressed only in terms of the relative surface approach. The negative sign of the reflected wave indicates propagation in the opposite direction to the incident wave through the upper body.

Determination of Contact Acoustic Nonlinearity (CAN)
An interface consists of asperity contacts and air gaps. In the regions of the interface with the air gap, there is no contact between the surfaces (Figure 2a). The approach of the two surfaces depends on the magnitude of the applied pressure [40]. The applied pressure is the sum of the externally applied pressure and the pressure applied by the incident wave. A longitudinal incident cosinusoidal wave with centre angular frequency and amplitude strikes the interface at time . A longitudinal incident cosinusoidal wave with centre angular frequency ω and amplitude A strikes the interface at time t 0 .
It is assumed at time t 0 = 0 that the gap is closed and remains closed until just before time t 1 . For the sake of derivation, the gap is assumed to be fully closed; however, in reality even under high pressure there will remain air gaps and the contact does not fully close. Therefore, the gap separation h and relative surface separation Y are zero, subsequently, the initial surface separation h 0 is zero. When the gap is open at time t 1 , the pressure at the interface is zero, so the boundary condition in Equation (4) is: Substituting Equation (20) into Equation (19) and integrating with respect to time gives time t 1 : Equation (21) is different from those presented in [38]. Equation (21) is a nonlinear equation and in the present work it is solved by the Newton-Raphson method.
The gap remains open until time t 2 . Integration of Equation (19) with respect to time and substituting Equation (20) in the interval [t 1 , t 2 ] gives opening interval: Hence at time t 2 the gap is just about to be closed; the right-hand side of Equation (21) is zero. Time t 2 can be found using the Newton-Raphson method: The gap remains closed until time t 3 . The opening-closing gap transition is repeated periodically: where n = 0, 1, 2, . . . and T is the period of the incident wave. A single longitudinal incident ultrasonic wave cycle consists of two parts causing tension and compression (Figure 2b). Figure 3 shows surface separation at the interface of aluminium blocks subjected to a longitudinal wave with centre frequency 2 MHz and amplitude 20 nm. The gap is initially assumed to be closed at time t 0 = 0. It is seen from Figure 3a that, when the compressional part of wave reaches the interface (corresponding to zero to point a on the incident wave), the applied pressure is sufficiently large to keep the gap closed in the interval [t 0 , t 1 ] (also see Figure 2c). As the wave passes (corresponding to part of the wave between point a and b on the incident wave), the applied pressure reduces so the surfaces of the gap are open at time t 1 and be pulled apart away from the mean position (also see Figure 2d). After that, the applied pressure increases (corresponding to part of the wave between point b and c on the incident wave) and presses the surfaces of the gap closer together. Figure 3a that, when the compressional part of wave reaches the interface (corresponding to zero to point a on the incident wave), the applied pressure is sufficiently large to keep the gap closed in the interval [ , ] (also see Figure 2c). As the wave passes (corresponding to part of the wave between point a and b on the incident wave), the applied pressure reduces so the surfaces of the gap are open at time and be pulled apart away from the mean position (also see Figure 2d). After that, the applied pressure increases (corresponding to part of the wave between point b and c on the incident wave) and presses the surfaces of the gap closer together.  In the time interval [t 2 , t 3 ] the gap is closed as the applied pressure increases (corresponding to part of the wave between point c and d on the incident wave). The openingclosing gap transition is repeated periodically. Figure 3b shows solid-solid contact and solid-air contact. Time domain signals presented in Figure 3b,c were converted to frequency domain using Fast Fourier Transform (FFT). It is seen that the surface of a solid-air contact freely moves, therefore no higher harmonic is generated. However, higher harmonics are generated (both even and odd) in solid-solid contact (Figure 3d). Figure 3c shows the effect of nominal contact pressure on the opening and closing the gap in the interface. It is seen that the higher nominal contact pressure closes the gap for longer period and the surfaces are separated less.

Theoretical Analysis
Evaluation of the general equations presented in Section 2.1 requires further theoretical analysis to drive the equations of interfacial stiffness.

Linear and Nonlinear Interfacial Stiffness
The stiffness of an interface can be modelled as a series of springs created by the individual asperity contacts [20,24]. The series of the springs is then equivalent to a single distributed spring with an equivalent stiffness K (expressed per unit area) [2]: Appl. Sci. 2021, 11, 5720 9 of 20 For practical rough surface contacts, the spring is nonlinear since the relation between the applied pressure and surface separation is not directly proportional. As the contact pressure is increased, more asperity contact occurs and the interface becomes stiffer. If the deflection is small (and less than the elastic limit of the spring material [42]) the nonlinearity can be neglected. For large deflections, the pressure-separation relationship and the resulting nonlinear stiffness must be considered. One of the approaches to deal with the nonlinear contact pressure-relative surface approach of the interface (boundary or contact nonlinearity) is to express it as a series of polynomial terms in a Taylor series around h = h 0 [39]: where K 1 and K 2 are the linear and second order nonlinear interfacial stiffness and p 0 is the nominal contact pressure at h 0 . The linear K 1 and nonlinear K 2 terms of the interfacial stiffness are then defined by:

Incident and Reflected Elastic Waves
Differentiation of Equation (18) with respect to argument of the function at the interface Substituting Equations (26) and (29) into (14) gives: which is a first order nonlinear differential equation. An approximate solution, such as the homotopy perturbation method (HPM) [43], is required to derive the particular solution of the relative approach of the contacting surfaces. To do this, the solution of the differential equation is defined as the sum of the linearized equations of the differential equation which is known as a perturbation series [43]: where Y 0 , Y 1 , Y 2 and Y 3 are the first, second, third and fourth terms of the perturbation series and α is an imbedding parameter which is assumed to be 1 to derive the particular solution of the differential equation. It should be noted that α is a small parameter in the range of [0, 1] [43]. As α approaches 1, Equation (31) gives an approximate solution to the nonlinear differential Equation (30) [43]. In order to solve the nonlinear differential equation (Equation (30)), the differentiation of Equation (31) with respect to time is substituted into Equation (30). A series of differential equations in terms of only Y 0 , Y 1 , Y 2 and Y 3 are derived. The differential equations containing Y 0 and Y 2 are in transient state (homogeneous differential equations). The initial surface approach is zero, subsequently, Y 0 = Y 2 = 0. Therefore, the transient behaviour of the differential equation was ignored and only the steady state behaviour was considered.
Substituting Equations (32) and (33) into Equation (31) gives the relative surface approach Y(t): where δ 1 = tan −1 (ρcω/2K 1 ) and δ 2 = tan −1 (K 1 /ρcω). Substituting Equation (34) into Equation (17) results in an analytical equation for the reflected wave: The terms in this equation with a frequency ω and 2ω represent the fundamental and second order harmonics, respectively. So, the coefficient of terms containing ω and 2ω are the amplitude of the fundamental frequency A 1 and second harmonic A 2 , respectively. The reflection coefficient R is defined as the ratio of the amplitude of the reflected pulse from the interface to the amplitude of the incident wave [2]. The reflection coefficient R at imperfect contact is dependent on the interfacial stiffness per unit area of the interface (where the wavelength of elastic wave is large compared to the air gap). The reflection coefficient of the fundamental frequency is given by [19]: where A 1 is the amplitude of the fundamental frequency of reflected pulse from the contact interface. In practice, A 1 can be found from the experimental pulses from the amplitude of Fast Fourier Transform (FFT) of the fundamental frequency. Rearranging Equation (36) gives the linear interfacial stiffness K 1 in terms of reflection coefficient: The second order nonlinear parameter for the reflected ultrasound from the interface γ (derived from Equation (35)) is defined by [39]: where A 2 is the amplitude of the second harmonic of the reflected wave. In the analytical approach (Equation (38)), A 2 is the coefficient of the 2ωt term. Again, experimentally, A 2 can be obtained from the amplitude of the second order harmonics obtained using a Fast Fourier Transform (FFT). It should be noted that the amplitude of the third harmonic is small compared to the amplitude of the second harmonic, therefore only the second harmonic is considered here. K 2 is the second order nonlinear interfacial stiffness. The second order nonlinear parameter for reflected ultrasound γ can be measured by the experiment with the values of A 1 and A 2 . In order to derive the nonlinear stiffness per unit area of the interface K 2 , Equation (38) is rearranged: The nonlinear nominal contact pressure relationship (Equation (26)) can then be defined at any nominal contact pressure p 0 with the corresponding linear and nonlinear interfacial stiffness, K 1 and K 2 from Equations (37) and (39).

Experimental Approach
In this experiment the reflection coefficient from an interface subjected to normal applied load was measured. Two transducers were placed on the upper specimen to act in pitch-catch mode. Reflections were recorded as loads were applied to the specimen pair. Figure 4 illustrates the experimental apparatus designed for capturing reflected pulses from the interface. Two cylindrical aluminium 6082 T6 blocks, of diameter 39 mm and thickness 32 mm, were loaded together creating a dry contact interface. The load was applied using with a hydraulic material testing machine and measured using a load cell. A spacer was used to protect the sensor location as shown in the figure.  The two specimen contact faces were polished with silicon carbide papers to P2500 grit size and lapped. The arithmetic mean deviation , root mean square and mean roughness depth were measured according to ISO 4287 with an optical roughness measurement device (Alicona SL) as shown in Table 1. Whilst the mean contact pressure is well below the yield stress of the materials, the asperity contacts are at much higher pressure [18] and so experience plastic flow. The plastic flow is greatest at the first loading and then the contact points conform or build up residual stress so subsequent loading is elastic. To ensure all measurements were recorded for an elastic contact, prior to the capturing, 10 loading/unloading cycles at a nominal contact pressure up to 5.5 MPa (which is beyond the nominal contact pressure decided to compress the contacting surfaces (5 MPa)) were applied. Further deformation of the asperity contacts was then expected to be fully elastic.

Instrumentation
Two piezoelectric longitudinal transducers (18 mm diameter) with centre frequencies of 2 MHz and 5 MHz were attached on the upper surface of the top aluminium blocks using a couplant gel. The transducer generated the incident wave and caught the signal reflected from the interface. A 2 MHz longitidinal toneburst with 11 cycles at different peak-to-peak excitation voltages of 90 V, 140 V or 210 V was generated and amplified using a high-power amplifier (RITEC RAM-5000). The power generated at these voltages was sufficiently large to generate higher harmonics in reflected pulses from the interface. In this experiment, higher harmonics were observed when the incident wave had a voltage greater than 60 V peak-to-peak. The reflected pulses from the contact in- The two specimen contact faces were polished with silicon carbide papers to P2500 grit size and lapped. The arithmetic mean deviation R a , root mean square R q and mean roughness depth R z were measured according to ISO 4287 with an optical roughness measurement device (Alicona SL) as shown in Table 1. Whilst the mean contact pressure is well below the yield stress of the materials, the asperity contacts are at much higher pressure [18] and so experience plastic flow. The plastic flow is greatest at the first loading and then the contact points conform or build up residual stress so subsequent loading is elastic. To ensure all measurements were recorded for an elastic contact, prior to the capturing, 10 loading/unloading cycles at a nominal contact pressure up to 5.5 MPa (which is beyond the nominal contact pressure decided to compress the contacting surfaces (5 MPa)) were applied. Further deformation of the asperity contacts was then expected to be fully elastic.

Instrumentation
Two piezoelectric longitudinal transducers (18 mm diameter) with centre frequencies of 2 MHz and 5 MHz were attached on the upper surface of the top aluminium blocks using a couplant gel. The transducer generated the incident wave and caught the signal reflected from the interface. A 2 MHz longitidinal toneburst with 11 cycles at different peak-to-peak excitation voltages of 90 V, 140 V or 210 V was generated and amplified using a high-power amplifier (RITEC RAM-5000). The power generated at these voltages was sufficiently large to generate higher harmonics in reflected pulses from the interface. In this experiment, higher harmonics were observed when the incident wave had a voltage greater than 60 V peak-to-peak. The reflected pulses from the contact interface were captured with the longitudinal transducer with centre frequency 5 MHz and stored using a digital oscilloscope. A reference signal was captured from a solid-air contact to eliminate the inherent effect of the materials and consider only the influence of the interface. A PC running LabView was used to trigger the pulsing and receive the digitized reflection data. Figure 5a shows the first and second reflected pulses from the interface of the solid-air contact (reference). This solid-air reflection acts as a reference. Since almost all the sound energy will be reflected from a solid-air interface, this is equivalent to the incident signal. Figure 5d shows the same pulses from a solid-solid contact at nominal contact pressure 4 MPa. In the present analysis, only the first reflected pulse was used as this contains most information from the interface.

Signal Processing
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 21 pressure 4 MPa. In the present analysis, only the first reflected pulse was used as this contains most information from the interface. This first peak was extracted and a Hanning window function and a zero-pad function were applied to increase the resolution of the signals (Figure 5b,e). It should be noted that the type of window function as well as the length of zero-pad function are dependent on the quality of captured signals. In this study, other window functions, such as Hamming window function and Blackman window function were employed, but showed insignificant difference in the results, so the Hanning window function was adopted throughout.
The time domain signals cannot precisely (in most cases this is impossible) indicate the existence of a super-harmonic (higher harmonics), subharmonic, or a mix of both. Therefore, a Fast Fourier Transform was employed to convert the time domain reflected signal to the frequency domain as shown in Figure 5c,f. These peaks are clearly visible corresponding to the fundamental frequency, as well as both, and even and odd (second) harmonics. Since Figure 5c shows solid-air contact, the source of the higher order harmonics in this case cannot be due to the interfacial contact. This nonlinearity might be generated by the bulk material, bonding nonlinearity (transducers to the specimen), the transducers themselves, or in the amplifier circuitry [29,33,44]. A test to see whether the interface (in This first peak was extracted and a Hanning window function and a zero-pad function were applied to increase the resolution of the signals (Figure 5b,e). It should be noted that the type of window function as well as the length of zero-pad function are dependent on the quality of captured signals. In this study, other window functions, such as Hamming window function and Blackman window function were employed, but showed insignificant difference in the results, so the Hanning window function was adopted throughout.
The time domain signals cannot precisely (in most cases this is impossible) indicate the existence of a super-harmonic (higher harmonics), subharmonic, or a mix of both. Therefore, a Fast Fourier Transform was employed to convert the time domain reflected signal to the frequency domain as shown in Figure 5c,f. These peaks are clearly visible corresponding to the fundamental frequency, as well as both, and even and odd (second) harmonics.
Since Figure 5c shows solid-air contact, the source of the higher order harmonics in this case cannot be due to the interfacial contact. This nonlinearity might be generated by the bulk material, bonding nonlinearity (transducers to the specimen), the transducers themselves, or in the amplifier circuitry [29,33,44]. A test to see whether the interface (in the solid-solid case) is a source of harmonics is to increase the nominal contact pressure; this causes a reduction in the amplitude of higher harmonics indicating that the interface is also the source of nonlinearity generated in the reflected signals. To eliminate all sources of nonlinearity apart from those generated by the contact [45,46], Equation (38) is reformulated as: where A 1 re f and A 2 re f are amplitudes of fundamental frequency and second order harmonic from a solid-air contact, respectively. In this way all nonlinearities not associated with the interface are subtracted out. The amplitude of the fundamental frequency of the reflected pulse in solid-solid contact (Figure 5f) was divided by the amplitude of the corresponding harmonic from the reference solid-air contact (Figure 5c) to derive the reflection coefficient R (Equation (36)). For the second order nonlinear parameter γ, the amplitude of the fundamental frequency was divided by the square of the amplitude of the second harmonic (Equation (40)).

Voltage to Meter Conversion
Reflected and transmitted signals received by the transducer from the interface are in units of voltage, therefore, the second order nonlinear parameter for the reflected ultrasound from the interface γ is in V −1 . To compute the nonlinear stiffness K 2 , the nonlinear parameter for the reflected ultrasound amplitude from the interface γ must be in units of meters. A method for determining the amplitude (in meters) from the piezo voltage was required. A Laser Doppler vibrometer (Polytec with laser sensor head OFV 354 and the vibrometer controller OFV 2500) was employed, as shown in Figure 6. A laser beam was incident on the free surface of the aluminium block with a transducer connected to the back face. While the transducer generated longitudinal waves in the aluminium block at different excitation voltages from 70 V to 350 V at center frequency 2 MHz, the laser beam was reflected from the free surface and captured by the laser sensor head. The captured laser beam returns velocity information. A numerical integral with respect to time results in the displacement of the surface in meters. The same signal processing procedure as explained in Section 3.3 was used to measure the displacement of the surface at different piezo voltages and hence incident amplitudes. Figure 7 presents the linear variation of amplitude of the reflected signal from the interface in voltage unit against the displacement of solid-air contact in meter. The slope of the curve δ (V/m) gives voltage to meter conversion [47]: block at different excitation voltages from 70 V to 350 V at center frequency 2 MHz, the laser beam was reflected from the free surface and captured by the laser sensor head. The captured laser beam returns velocity information. A numerical integral with respect to time results in the displacement of the surface in meters. The same signal processing procedure as explained in Section 3.3 was used to measure the displacement of the surface at different piezo voltages and hence incident amplitudes.   Figure 8a shows the reflection coefficient for incident waves of excitations 90 V, 140 V and 210 V as a function of the nominal contact pressure, determined from the amplitude of the first order reflection from the solid interface divided by that from the solidair interface according to Equation (36). They were measured after 10 loading/unloading cycles to remove any asperity plasticity effects. It can be seen that the reflection coefficient decreases with increasing the nominal contact pressure due to increase in the real contact area.  Figure 8a shows the reflection coefficient R for incident waves of excitations 90 V, 140 V and 210 V as a function of the nominal contact pressure, determined from the amplitude of the first order reflection from the solid interface divided by that from the solid-air interface according to Equation (36). They were measured after 10 loading/unloading cycles to remove any asperity plasticity effects. It can be seen that the reflection coefficient decreases with increasing the nominal contact pressure due to increase in the real contact area.

Results and Discussion
The data of Figure 8a was used to determine the first order linear interfacial stiffness K 1 using Equation (37) (Figure 8b). Higher pressures cause a stiffer interface. As expected, the reflection coefficient and linear stiffness determined at any voltage is identical. The amplitude of the ultrasonic wave does not change the linear stiffness component.
140 V and 210 V as a function of the nominal contact pressure, determined from the amplitude of the first order reflection from the solid interface divided by that from the solidair interface according to Equation (36). They were measured after 10 loading/unloading cycles to remove any asperity plasticity effects. It can be seen that the reflection coefficient decreases with increasing the nominal contact pressure due to increase in the real contact area. The data of Figure 8a was used to determine the first order linear interfacial stiffness using Equation (37) (Figure 8b). Higher pressures cause a stiffer interface. As expected, the reflection coefficient and linear stiffness determined at any voltage is identical. The amplitude of the ultrasonic wave does not change the linear stiffness component.
A best fit curve to the linear interfacial stiffness shown in Figure 8b is a fourth-order polynomial expression: A best fit curve to the linear interfacial stiffness shown in Figure 8b is a fourth-order polynomial expression [22]: where a n is a constant coefficient. This coefficient will depend on factors such as surface topography and elastoplastic deformation of the asperities in contact.
The nominal contact pressure p 0 can be expressed in terms of the relative surface approach Y of the contacting surfaces. To do this, the linear interfacial stiffness K 1 (Equation (42)) was substituted into Equation (27) creating a first order nonlinear homogeneous differential equation: dp dY The initial condition indicates that in the absence of elastic wave (ultrasound), there is zero relative surface approach. The homotopy perturbation method (HPM) was employed to solve Equation (43) with respect to relative surface approach: where b n is a constant coefficient. This coefficient, as with a n , will also depend on the surface asperity and material properties. A numerical differential equation (ode45 function in MATLAB [48]) of the data in Figure 8b according to Equation (43) was used to determine the nominal contact pressure-relative surface approach relationship, as shown in Figure 8c. The shape of the nonlinear nominal contact pressure-relative surface approach curve demonstrates the stiffening behaviour of the interfacial spring under increasing contact pressure. In other words, increasing the nominal contact pressure leads to a reduction in the gap, more solid contact, and hence a stiffer interface. In the absence of an elastic wave, the relative surface approach was zero (Y = 0) under nominal contact pressure p 0 = 0 MPa (nominal contact pressure generated by the weight of upper body is negligible). The negative sign of the relative surface approach Y indicates the reduction of the interface separation (greater deflection in the asperities). The gradient of the curve increases with increase in nominal contact pressure demonstrating stiffer interfacial contact. Figure 9a shows the second order nonlinear parameter γ at excitations 90 V, 140 V, and 210 V as a function of the nominal contact pressure. This nonlinear parameter is determined from the amplitude of the second harmonic divided by the square of the first harmonic (Equation (40)). To evaluate the second order nonlinear stiffness K 2 , the experimental results of the second order nonlinear parameter γ were substituted into Equation (39) (Figure 9b). It is seen from Figure 9b that the second order nonlinear stiffness increases as the nominal contact pressure increases up to 2.5 MPa. However, beyond this value the nonlinearity decreases, so the nonlinear stiffness decreases. The comparison shows that is independent of the amplitude of the incident wave. A best fit curve to the second order nonlinear interfacial stiffness shown in Figure 9b is a fourth-order polynomial expression: where is a constant coefficient. This coefficient depends on some factors such as surface topography and elastoplastic deformation of the asperities in contact. However, further studies are required. As before, the nominal contact pressure is expressed in terms of the relative surface approach of the contacting surfaces. To do this, the nonlinear interfacial stiffness (Equation (45)) was substituted into Equation (28) creating a second order nonlinear homogeneous differential equation: The initial conditions indicate that in the absence of an elastic wave (ultrasound), there is zero relative surface approach and linear interfacial stiffness. The homotopy perturbation method (HPM) was again employed to solve Equation (46) with respect to relative surface approach. A numerical differential equation (ode45 function in MATLAB [48]) of the data in Figure 9b according to Equation (46) was used to determine the nominal contact pressure-relative surface approach relationship, as shown in Figure 9c. The result was the same as Equation (44).
The comparison between the pressure-relative surface approach derived from linear and nonlinear interfacial stiffness is presented in Figure 10. Figure 10a shows the pressurerelative surface approach derived from the present experimental data. Figure 10b is included for comparison, and shows the pressure-relative surface derived from the experimental data from [22]. It should be noted that the experiment presented in [22] was through transmission while the current paper considered reflected signal. Both approaches yield a similar relationship. It is seen from Figure 9b that the second order nonlinear stiffness increases as the nominal contact pressure increases up to 2.5 MPa. However, beyond this value the nonlinearity decreases, so the nonlinear stiffness decreases. The comparison shows that K 2 is independent of the amplitude of the incident wave. A best fit curve to the second order nonlinear interfacial stiffness K 2 shown in Figure 9b is a fourth-order polynomial expression [22]: where c n is a constant coefficient. This coefficient depends on some factors such as surface topography and elastoplastic deformation of the asperities in contact. However, further studies are required. As before, the nominal contact pressure p 0 is expressed in terms of the relative surface approach Y of the contacting surfaces. To do this, the nonlinear interfacial stiffness K 2 (Equation (45)) was substituted into Equation (28) creating a second order nonlinear homogeneous differential equation: The initial conditions indicate that in the absence of an elastic wave (ultrasound), there is zero relative surface approach and linear interfacial stiffness. The homotopy perturbation method (HPM) was again employed to solve Equation (46) with respect to relative surface approach. A numerical differential equation (ode45 function in MATLAB [48]) of the data in Figure 9b according to Equation (46) was used to determine the nominal contact pressure-relative surface approach relationship, as shown in Figure 9c. The result was the same as Equation (44).
The comparison between the pressure-relative surface approach derived from linear and nonlinear interfacial stiffness is presented in Figure 10. Figure 10a shows the pressure-relative surface approach derived from the present experimental data. Figure 10b is included for comparison, and shows the pressure-relative surface derived from the experimental data from [22]. It should be noted that the experiment presented in [22] was through transmission while the current paper considered reflected signal. Both approaches yield a similar relationship. It is seen that the pressure-approach relationship derived from the linear stiffness and second order nonlinear stiffness are very similar for all excitation voltage cases. As expected, the pressure-approach response should not depend on whether it was determined from either linear or nonlinear stiffness measured data. This gives confidence in the robustness of the measurement method and the analysis approach. The interfacial spring can be modelled with a stiffness described by a Taylor expansion with first and second order components (Equation (26)). Further, an experimental measurement of either one of these components, using ultrasonic reflection, is sufficient to deduce the other and hence characterise the spring stiffness.

Conclusions
In this work, the stiffness of a rough surface contact has been expressed in terms of a linear and second order nonlinear stiffness. Both these interfacial stiffness components were determined from the refection of a high-power ultrasonic wave. The high-power wave causes a non-linear response of the interface and generates higher order harmonics. The first order linear stiffness was deduced from the reflection of the fundamental frequency wave; whilst the second order nonlinear stiffness was deduced from the second harmonic frequency reflection. A fourth-order polynomial expression for the linear and nonlinear interfacial stiffness as a function of nominal contact pressure were defined. A first order nonlinear homogeneous differential equation of linear stiffness and a second order nonlinear homogeneous differential equation of the nonlinear stiffness with respect to relative surface approach was used to determine the relationship between the nominal contact pressure and relative surface approach. The results derived from both linear and nonlinear interfacial stiffness were similar. This indicates both linear and nonlinear interfacial stiffness subjected to the same nominal contact pressure can be measured by ultrasound. Further, the Taylor series used to represent the contact pressure-relative surface approach relationship can accurately represent the interface behaviour.
Author Contributions: Conceptualization, Taghizadeh, S. and Dwyer-Joyce, R.S; methodology, Taghizadeh, S. and Dwyer-Joyce, R.S; software, Taghizadeh, S; validation, Taghizadeh, S; formal analysis, Taghizadeh, S and Dwyer-Joyce, R.S; investigation, Taghizadeh, S. and Dwyer-Joyce, R.S; It is seen that the pressure-approach relationship derived from the linear stiffness K 1 and second order nonlinear stiffness K 2 are very similar for all excitation voltage cases. As expected, the pressure-approach response should not depend on whether it was determined from either linear or nonlinear stiffness measured data. This gives confidence in the robustness of the measurement method and the analysis approach. The interfacial spring can be modelled with a stiffness described by a Taylor expansion with first and second order components (Equation (26)). Further, an experimental measurement of either one of these components, using ultrasonic reflection, is sufficient to deduce the other and hence characterise the spring stiffness.

Conclusions
In this work, the stiffness of a rough surface contact has been expressed in terms of a linear and second order nonlinear stiffness. Both these interfacial stiffness components were determined from the refection of a high-power ultrasonic wave. The high-power wave causes a non-linear response of the interface and generates higher order harmonics. The first order linear stiffness was deduced from the reflection of the fundamental frequency wave; whilst the second order nonlinear stiffness was deduced from the second harmonic frequency reflection. A fourth-order polynomial expression for the linear and nonlinear interfacial stiffness as a function of nominal contact pressure were defined. A first order nonlinear homogeneous differential equation of linear stiffness and a second order nonlinear homogeneous differential equation of the nonlinear stiffness with respect to relative surface approach was used to determine the relationship between the nominal contact pressure and relative surface approach. The results derived from both linear and nonlinear interfacial stiffness were similar. This indicates both linear and nonlinear interfacial stiffness subjected to the same nominal contact pressure can be measured by ultrasound. Further, the Taylor series used to represent the contact pressure-relative surface approach relationship can accurately represent the interface behaviour.