*2.3. D Numerical Model*

/௫ୀ

୦ ଶ୫னబ

*ξ*=

A one-dimensional numerical model of the experimental set-up was implemented (Figure 2) consisting of the guide bar Ω<sup>1</sup> of the experimental set-up and the tested aluminium sample Ω2, modelled by unidimensional deformable bodies.

**Figure 2.** Diagram (**a**) and numerical model (**b**) of the set-up for contact stiffness experiments.

Two contacting interfaces are considered: the first between the guide bar Ω<sup>1</sup> and the tested sample Ω2, and the second between the tested sample and the frame, considered as infinitely rigid in the model (tribometer disc). The model parameters are provided in Table 2. Note that, in reality, the cross-sections of the guide bar Ω<sup>1</sup> and aluminium sample Ω<sup>2</sup> are different. Numerically, an equivalent 1D-model with the same cross-section S<sup>2</sup> is considered, corresponding to the contact surface S2.

∂ ଶu x<sup>ଶ</sup> + ℎ

∂ ଶ ∂t<sup>ଶ</sup> − <sup>ଶ</sup>

 ଶ

∂u ∂t =

*ξ*


**Table 2.** Geometry and material parameters of Ω<sup>1</sup> and Ω<sup>2</sup> in the numerical model.

For both interfaces, a nonlinear contact law is introduced. This law includes a nonlinear stiffness function of the contact pressure.

Once a force is applied to the top of the system, longitudinal waves propagate along the vertical x-axis and excite the longitudinal modes of the system, including the mass-spring mode, where the mass is the guide bar and the spring is due to the stiffness in series of the sample and the two contact interfaces.

The numerical modelling assumes the following hypotheses:


$$
\sigma\_{/x} = 0 \,\,=\,\frac{F}{S\_2}
$$

where σ/*<sup>x</sup>* <sup>=</sup> <sup>0</sup> is the stress at the upper side of the guide bar.

a. Governing equations

In the absence of elastic wave, the system is assumed to be free at the upper part and subject only to its own weight. In a first static step, the equilibrium position of the system can be calculated.

For the longitudinal waves propagating along the *x*-direction (longitudinal), in both the guide bar and the aluminum sample, the equation of motion is the following

$$\frac{\partial^2 \mathbf{u}}{\partial \mathbf{t}^2} - c^2 \frac{\partial^2 \mathbf{u}}{\partial \mathbf{x}^2} + h \frac{\partial \mathbf{u}}{\partial \mathbf{t}} = \text{ g} \tag{1}$$

where *u* (*x*,t) is the displacement in the *x*-direction at time *t*, from the non-deformed configuration, *c* is the wave velocity, *g* is the gravity acceleration, *h* is the viscosity factor and *t* denotes the time.

The viscosity factor *h* is constant in the model and was determined from the damping factor ξ = <sup>h</sup> 2mω<sup>0</sup> , where ω<sup>0</sup> is the angular frequency of interest. The damping factor value ξ = 0.035 was identified experimentally by logarithm decrement at the frequency of interest, which is the frequency of the corresponding mass-spring mode. In fact, the main mode of vibration, the focus of the analysis, is the response of the guide mass to the stiffness obtained by the bulk stiffness of the sample in series with the stiffness of the two contact interfaces (Figure 2).

The model is discretized in both time and space using the Euler scheme [31].

The unidimensional wave Equation (1) has been discretized using a second order centered finite difference (FD) scheme, explicit in time, as follows

$$\begin{array}{ll} \frac{u\_k^{m+1} - 2u\_k^m + u\_k^{m-1}}{\delta t^2} & -c^2 \frac{u\_{k+1}^m - 2u\_k^m + u\_{k-1}^m}{\delta x^2} + h \frac{u\_k^{m+1} - u\_k^m}{\delta t} \\ & = g \ \ \forall k \varepsilon [2, n\_1 - 1] \cup \ [n\_1 + 2, n\_1 + n\_2 - 1] \end{array} \tag{2}$$

where the subscript *k* indicates the node number and *m* indicates the time step number, with δ*x* and δ*t* corresponding, respectively, to the space and time discretisation steps. *n*<sup>1</sup> and *n*<sup>2</sup> are, respectively, the total number of nodes contained in solids Ω<sup>1</sup> and Ω2. This second-order general scheme is conditionally stable under the classic Courant–Friedrichs–Lewy (CFL) condition:

$$\frac{c\delta t}{\delta x} < 1\tag{3}$$

Equation (2) can be re-written to obtain the displacement at the next time step *m* + 1 for each internal node, as follows

$$u\_k^{m+1} = u\_k^m \left(\frac{2(1-\alpha\_1) + \alpha\_2}{1+\alpha\_2}\right) + \frac{\alpha\_1}{1+\alpha\_2} \left(u\_{k+1}^m + u\_{k-1}^m\right) + \frac{1}{1+\alpha\_2} \left(gdt^2 - u\_k^{m-1}\right) \tag{4}$$

$$\forall k \epsilon \left[2, n\_1 - 1\right] \cup \left[n\_1 + 2, n\_1 + n\_2 - 1\right]$$

where α<sup>1</sup> = *<sup>c</sup>* <sup>2</sup>δ*t* 2 δ*x* <sup>2</sup> and α<sup>2</sup> = *h*. *dt*

At contact interfaces (at x = *xC*<sup>1</sup> between the two solids Ω<sup>1</sup> and Ω<sup>2</sup> and x = *xC*<sup>2</sup> between Ω<sup>2</sup> and the frame), the contact stresses σ(x = *xC*1,t) and σ(x = *xC*2,t) has the following expression in compression

$$
\sigma(\mathbf{x}\_{\text{Ci}}, t) = K(\sigma(\mathbf{x}\_{\text{Ci}}, t)) \Delta u\_{\text{Ci}} \tag{5}
$$

where K denotes the nonlinear stiffness of the contact interface depending on the contact stress and ∆*uCi* the gap distance at interfaces. i indicates either the contact between the two solids Ω<sup>1</sup> and Ω<sup>2</sup> (i = 1) or the contact between Ω<sup>2</sup> and the frame (i = 2).

In order to evaluate properly the solution at both contact interfaces (σ(*xCi*, *t*), *K*(σ(*xCi*, *t*)) and ∆*uCi*), it is necessary to solve Equation (5), and thus apply an implicite scheme (finite difference), using the relation between the contact pressure and the gap distance extracted from the implemented nonlinear relation between the stiffness and the pressure, presented in what follows.
