## 3 Contact Model
### Determination of Stiffness, \(k_{n}\) and \(k_{t}\)  
The global matrix \([(K_{GLOB})_{FF}]\) is a function of the parameters \(k_{n}\) and \(k_{t}\). To compare the natural frequencies with experimental outcomes, they need to be measured during the complete adhesion of the contact regions. Specifically, by examining the concerning curves (see Figures 9 and 11), the bandwidth in which the frequency has a rather constant trend was selected, to indicate the complete contact adhesion. The frequencies of FEM model (for each axial traction force) were obtained from the classical linear modal analysis. Considering the global matrices \([M]=[(M_{GLOB})_{FF}]\) and \([K]=[(K_{GLOB})_{FF}]\), the governing equations of motion are:  
\[[M]\{\tilde{q}(t)\}+[K_{i}(k_{n},\ k_{t})]\{q(t)\}=0\ \ \ \forall F_{ax}\ i-th \tag{16}\]  
Figure 17: Assembly of contact elements into global stiffness matrix. With the following linear harmonic solution:  
\[\{q(t)\}=\{q\}\ e^{i\omega t}\ \Rightarrow\ \left(-\omega^{2}M+K_{i}\right)\{q\}=0 \tag{17}\]  
Excluding the trivial solution \(\{\mathsf{q}\}=0\), other possible solutions (eigenvectors) are obtained for the eigenvalues of the generalized problem \(\omega^{2}\), for which the matrix \(\left(-\omega^{2}M+K\right)\) is not invertible. This condition is verified when:  
\[det\left(-\omega^{2}M+K_{i}(k_{n},\ k_{t})\right)=0\quad\forall F_{i} \tag{18}\]  
from which \(\omega_{n,i}^{2}\) are the roots of the characteristic polynomial and the eigenvalues:  
\[\sum_{j=0}^{Ndef}a_{j}(k_{n},\ k_{t})\Big{(}\omega_{n,\ i}^{2}\Big{)}^{j}=0\ \ \ \forall F_{i} \tag{19}\]  
The shape of each depends on the coefficients \(a_{j}\), which are nonlinear functions of the parameters \(k_{n}\) and \(k_{t}\). Thus, the frequencies depend non-linearly on these parameters, which required an iterative local or global method to find the solution \((k_{n},k_{t})\). By making a comparison between the measured frequencies and those from the FEM model, for the first and second flexural modes (1B and 2B):  
\[ \begin{array}{l}\left(f_{1B}\right)_{FEM}(k_{n},\ k_{t})-\left(f_{1B}\right) _{Measured}=0\\ \left(f_{2B}\right)_{FEM}(k_{n},\ k_{t})-\left(f_{2B}\right)_{Measured}=0\end{array}  \tag{20}\]  
which can be written in the following vector form:  
\[\{\Delta f\}=\left\{ \begin{array}{c}\left(\Delta f\right)_{1B}(k_{n},\ k_{t })\\ \left(\Delta f\right)_{2B}(k_{n},\ k_{t})\end{array} \right\}=\left\{ \begin{array} []{c}0\\ 0\end{array} \right\} \tag{21}\]  
The function \(\{\Delta f\}\) does not vanish at any point, i.e., there are no \((k_{n},k_{t})\) such as to satisfy both Equation (20); however, it can be minimized. Therefore, the system of equations becomes:  
\[\{|\Delta f|\}<\delta\{\Delta f\} \tag{22}\]  
This system was solved with a global search within the domain \((k_{n},k_{t})\) for satisfying the following condition:  
\[||\{\Delta f\}||<\delta||\{\Delta f\}|| \tag{23}\]  
where \(||\{\Delta f\}||\) is the Euclidean norm of the vector \(\{\Delta f\}\) and \(\delta||\{\Delta f\}||\) is the error, calculated with the theory of error propagation and \(\sigma\) defines the stress. Following the procedure explained above, the parameters \(k_{n}\) and \(k_{t}\) are obtained for i-th axial force and reported in Table 1.  
Eventually the frequencies obtained from the experiment are compared with those from the presented numerical approach and the analytical method in Figure 18. The