## 2 Materials and Methods
### Virtual Material-Based Acoustic Model
#### 2.1.2 Equivalent Elastic Modulus  
Based on the Hertz elastic contact theory and the normal stress-strain relation of an asperity, the elasticity modulus of a single asperity can be expressed as  
\[e=\frac{2E}{3\pi}\sqrt{\frac{\beta}{w}} \tag{14}\]  
where \(w\) is the deflection of the asperity, which is affected by the asperity interactions (see Figure 2).  
The elasticity modulus of virtual material can be expressed as  
\[E_{0}=N\int_{h_{n}-d_{n}}^{\infty}e\phi_{n}(z_{n})\mathrm{d}z_{n} \tag{15}\]  
Since the asperity interactions can increase the surface separation and reduce the asperity deflection, the calculation of \(E_{0}\) should consider the influence of asperity interactions.  
1.3 P\({}_{\mathrm{m}}\)-K\({}_{\mathrm{n}}\)-R Relationship Considering Asperity Interactions  
Considering the influence of asperity interactions, the asperity deformation due to the contact load can be expressed as  
\[w=z+d_{d}-h+\left(u_{m}-u_{g}\right) \tag{16}\]  
where \(u_{m}\) is the displacement of the mean of asperity heights. According to the asperity interactions model proposed by [22], the deflection of the asperity can be expressed as  
\[w=z+d_{d}-h+1.3\sqrt{\beta^{1/2}w^{3/2}\frac{p_{m}}{E}} \tag{17}\]  
where \(p_{m}\) is the global mean contact stress. All the length parameters and variables in Equation (17) are normalized by dividing \(\sigma\). The resulting dimensionless equation is given by  
\[w_{n}=z_{n}+d_{n}-h_{n}+1.3\sqrt{\beta^{1/2}w^{3/2}\frac{p_{m}}{E}} \tag{18}\]Combining Equations (5), (14)-(17), the elasticity modulus can be expressed as  
\[E_{0}=\frac{2 NE}{3\pi}\sqrt{\frac{\beta}{\sigma}}\int_{h_{n}-d_{n}}^{\infty}w_{n }{}^{-1/2}\phi_{n}(z_{n})\mathrm{d}z_{n} \tag{19}\]  
With the ZMC model [23], the dimensionless contact load \(W_{tn}\) can be expressed as  
\[W_{tn}=\frac{p_{m}}{E}=\frac{4n\sigma\beta}{3\pi}\sqrt{\frac{\sigma}{\beta}} \int_{h_{n}-d_{n}}^{\infty}w_{n}{}^{3/2}\varnothing_{n}(z_{n})\mathrm{d}z_{n} \tag{20}\]  
Since the contact load \(W_{tn}\) and surface separation \(h_{n}\) are coupled together, the Newton-Raphson iterative formula and five-point Legendre-Gaussian quadrature are used to solve the above governing equations. As a result, the relationship of \(P_{m}\)-\(K_{n}\)-\(R\) can be obtained.  
The following describes the computational procedure in this work.  
1. Input parameters of the problem (\(m_{0}\), \(m_{2}\), \(m_{4}\), \(R\), \(E\), \(G\), \(S_{y}\), \(A_{n}\)).
2. For a given load \(W_{tn}\), the value of \(w_{n}\) as a function of \(h_{n}\) and \(z_{n}\) can be solved from Equation (18) using the Newton-Raphson iterative formula.
3. By adjusting the value of \(h_{n}\), the calculated load of Equation (20) using the five-point Legendre-Gaussian quadrature can be made equal to the given load, and the relation between \(W_{tn}\) and \(h_{n}\) can also be determined.
4. By selecting a different load, the corresponding mean separation \(h_{n}\) can be obtained following steps 2 and 3. Therefore, the relations between \(P_{m}\) and \(h_{n}\) and between \(P_{m}\) and \(R\) can be obtained.
5. Substituting \(h_{n}\) into Equation (19), \(E_{0}\) can be calculated using the five-point Legendre-Gaussian quadrature. \(G_{0}\) can be obtained by integrating Equation (12).
6. Finally, \(R\), \(E_{0}\) and \(G_{0}\) are substituted into Equations (4) and (6) to calculate the contact stiffness \(K_{n}\) under the given contact condition, and the relation of \(P_{m}\)-\(K_{n}\)-\(R\) can be determined.