*2.2. MPJR Formulation*

For the solution of the contact problem, the MPJR interface finite element that is exposed in [6,8] is employed. It consists in a 4-nodes, zero-thickness element mutuated by the Cohesive Zone Model (CZM) and used in the context of nonlinear fracture mechanics. The framework is applied to the problem of a rigid body with a complex boundary making contact with another deformable body characterised by a smooth interface. The core of the approach is re-casting the original geometry of

the problem into a simpler one, consisting only in the deformable bulk and a single layer of interface elements disposed at its boundary, where contact is supposed to take place, as in Figure 1. The actual shape of the indenter is stored nodal-wise in each interface finite element employed for the interface discretization, and it is used to correct the normal gap function that is computed from the flat–flat configuration, in order to account for the exact geometry. This requires a preliminary step, which consists in mapping the indenter profile elevation in correspondence to the right node of the boundary. If the profile has analytic expression, this can be done right at the finite element level exploiting the global coordinates, otherwise the elevation field can be stored in a proper history variable and every entry associated with the correct node. It has to be remarked that, in spite of the present formulation being 2*D*, the proposed framework can be extended to 3*D* problems, provided that, for example, a 8-nodes interface element is used to discretize a surface, instead of a profile, equipped with a suitable friction law.

**Figure 1.** Profile discretization and equivalent interface definition. The interface element Γ is defined with the lower two nodes that belong to the deformable bulk, and the others placed at a given offset normal to the lower boundary. An abscissa *s* can be defined along the boundary to map the indenter's elevation field, which is stored inside the element and it is used to correct the normal gap.

Figure 2 shows the kinematics of the element. A vector of unknown nodal displacements **u** = [*u*1, *v*1, . . . , *u*4, *v*4] *T* is introduced for the evaluation of the tangential and normal gaps, collected in the vector **g** = [*gx*, *gz*] *T* , which reads:

$$\mathbf{g} = \mathbf{QNLu} \tag{1}$$

where **L** is a linear operator for computing the relative displacements across the interface, **N** is the shape functions matrix, and **Q** is a rotation matrix for transforming displacements from the global to the local reference frame of the element defined by the unit vectors *n* and *t*. The original geometry can be restored with a suitable correction of the normal gap, in the form:

$$\mathbf{g}^\* = \begin{bmatrix} \mathbf{g}\_x \\ \mathbf{g}\_z + h(\tilde{\xi}\_\prime t) \end{bmatrix} \,\prime \tag{2}$$

where *h*(*ξ*, *t*) maps the profile's shape and position in time. With respect to the formulation that is presented in [8], here the profile shape has been made time-dependent, in order to also account for finite sliding of the rigid indenter. For example, in the case of a flat interface, the result of the indenter sliding with a given constant velocity *v*<sup>0</sup> can be achieved by setting *h*(*ξ*, *t*) = *h*(*ξ* − *v*0*t*).

**Figure 2.** Four-nodes, zero-thickness eMbedded Profiles for Joint Roughness (MPJR) interface finite element.

*Lubricants* **2020**, *8*, 107

The current value of *t* is stored at the interface element level while using a time history variable, and it is updated every time step. A standard penalty approach is used in order to enforce the normal contact constraint, leading to

$$p\_z = \begin{cases} \arg\_z^\* \text{ if } \mathcal{g}\_z^\* < 0, \\ 0, \text{ if } \mathcal{g}\_z^\* \ge 0, \end{cases} \tag{3}$$

where *α* is the penalty parameter.

To deal with coupled frictional problems, a regularized Coulomb friction law [9] is used to set the interface constitutive equation in the tangential direction:

$$q\_{\mathbf{x}} = f p\_{\mathbf{z}} \tanh\left(\frac{\mathbf{g}\_{\mathbf{x}}}{\mathfrak{E}}\right),\tag{4}$$

where *q<sup>x</sup>* is the tangential traction and *g*˙*<sup>x</sup>* is the sliding velocity, as given by the difference between the velocity of the indenter and the horizontal velocity of the corresponding node. Finally, *ε*˙ is a parameter governing the slope of the regularised friction law.

The contribution of a single interface finite element to the variational formulation of the bulk material is expressed by its variation in terms of density of energy content, integrated over the domain Γ<sup>e</sup> that denotes the element itself:

$$
\delta \Pi\_{\mathbf{e}} = \int\_{\Gamma\_{\mathbf{e}}} \delta \mathbf{g}^\*(\mathbf{v})^T \mathbf{p}(\mathbf{u}, \dot{\mathbf{u}}) \, \mathrm{d}\Gamma\_{\mathbf{e}}.\tag{5}
$$

where **v** is the virtual displacement field, and the vector **p** collects the normal and tangential tractions. As a final step, the variation can be expanded and the integral set to zero, leaving the expression of the nonlinear residual vector, which reads:

$$\mathbf{R}\_{\mathbf{e}}(\mathbf{u}, \dot{\mathbf{u}}) = \int\_{\Gamma\_{\mathrm{e}}} \mathbf{L}^{T} \mathbf{N}^{T} \mathbf{Q}^{T} \mathbf{p}(\mathbf{u}, \dot{\mathbf{u}}) \, \mathrm{d}\Gamma\_{\mathrm{e}} = \mathbf{0}. \tag{6}$$

Because of the nonlinearity of **R**e, the Newton–Raphson iterative method has been applied, together with a backward Euler method for time integration.

#### *2.3. Rheological Model*

Three different Prony series models with a number of arms increasing from one to three are examined in order to assess the effect of viscoelasticity modelling on the overall contact mechanical response. The general equation for the shear relaxation modulus reads:

$$\frac{G(t)}{G^{\infty}} = \mu\_0 + \sum\_{n=1}^{3} \mu\_n \exp\left(-\frac{t}{\tau\_n}\right),\tag{7}$$

where *G* <sup>∞</sup> is the instantaneous shear modulus (evaluated at *t* = 0), *µ<sup>n</sup>* are the relaxation coefficients, and *τ<sup>n</sup>* are the corresponding relaxation times. Equation (7) has been tuned to fit the experimental values of EVA [4]. The model parameters for 1, 2 and 3 arms are collected in Table 1.


**Table 1.** Rheological parameters for Ethylene Vynil Acetate (EVA), where *n* is the number of Prony series' arms.

The identification of the above parameters has been carried out through a regression over the experimental data that were acquired in the time range *t* = 10[−1,...,+1] s. The following approach has been pursued in order to attain a high degree of accuracy. Firstly, trial relaxation times have been set and a preliminary linear regression has been performed involving *G* <sup>∞</sup> and *µ<sup>i</sup>* only. The objective function to be minimised reads:

$$\Pi(\mathbf{x}) = \sum\_{k=1}^{N} \left(\mathbf{g}\_k \cdot \mathbf{x} - \mathbf{G}\_k\right)^2,\tag{8}$$

where, for the three arms model, **g***<sup>k</sup>* = h 1, *e* (−*tk*/*τ*<sup>1</sup> ) , . . . ,*e* (−*tk*/*τ*3) i , *G<sup>k</sup>* is the value of the objective function at the sampling point and *N* is the number of samplings. The global minimiser **x** ∗ = arg min*<sup>x</sup>* Π(**x**) is evaluated and the values of the constants *µ<sup>i</sup>* and *G* <sup>∞</sup> are obtained according to:

$$\mathbf{G}^{\infty} \begin{bmatrix} \mu\_0 \\ \mu\_1 \\ \mu\_2 \\ \mu\_3 \end{bmatrix} = \mathbf{x}^\* \,\prime \,\tag{9}$$

together with the condition ∑*<sup>i</sup> µ<sup>i</sup>* = 1, related to the shear modulus at *t* = 0. The obtained coefficients, together with their respective relaxation times, have been used in order to define a vector of guess values **x**<sup>0</sup> for a second nonlinear regression, in which the relaxation times were also included in the optimisation vector **x**. The problem has been solved iteratively, updating the starting vector **x**<sup>0</sup> every cycle using the results that were obtained in the previous. Convergence is achieved within 5 iterations, when considering a relative error that is given by (**x** <sup>∗</sup> <sup>−</sup> **<sup>x</sup>**0)/**x**<sup>0</sup> and a tolerance *<sup>ε</sup>* <sup>=</sup> <sup>10</sup>−<sup>15</sup> . This procedure has also been repeated in the same way for the 1 and 2 arms models.

Once the parameters are identified, the Young's relaxation modulus *E*(*t*) can be obtained from *G*(*t*), and the behaviour of the three models can be investigated in time and frequency domains. The analysis in the frequency domain can be performed by defining a complex modulus *E*ˆ(*ω*), obtained via a Fourier transform of *E*(*t*), which can be expressed as:

$$\frac{\triangle(\omega)}{E^{\infty}} = \mu\_0 + \sum\_{i=1}^{n} \mu\_i \frac{\tau\_i^2 \omega^2}{1 + \tau\_i^2 \omega^2} + \imath \sum\_{i=1}^{n} \mu\_i \frac{\tau\_i \omega}{1 + \tau\_i^2 \omega^2}. \tag{10}$$

In the expression above, *ı* denotes the imaginary unit and the index *k* defines the number of arms being considered. It can be easily noticed that, for the single arm model, the maximum viscoelastic effect manifests in correspondence to the critical excitation frequency *ω*<sup>⋆</sup> = √*µ*0/*τ*1.

Figure 3a shows the plot of *E*(*t*). Figure 3b,c show the values of the loss modulus and the storage modulus, which were obtained as the imaginary part <sup>ℑ</sup>*E*ˆ(*ω*) and the real part <sup>ℜ</sup>*E*ˆ(*ω*) of the complex modulus *E*ˆ(*ω*), respectively. Finally, Figure 3d shows the loss tangent, given as the ratio of the imaginary part over the real part. As a comparison, the same quantities are also plotted for the relaxation modulus obtained for a model that is based on fractional calculus, which reads:

$$E\_{\mathbf{f}}(t) = \frac{E\_{\mathbf{f},\mathbf{a}}t^{-\alpha}}{\Gamma(1-\alpha)}.\tag{11}$$

In Equation (11), *E*f,*<sup>α</sup>* = 814.7 Pa s*<sup>α</sup>* and *α* = 0.226 have been chosen in order to fit the experimental data in [4], being <sup>Γ</sup>(·) the gamma function.

The simulation of the power-law viscoelastic response seen in the experiments, which is well approximated by the fractional calculus model, is progressively improved by increasing the number of terms in the Prony series representation. It has to be remarked that, since the Fourier transform of a power-law is a power-law itself, both loss and storage modulus in the frequency domain are represented, on a logarithmic scale, as straight lines.

Their trend can be satisfactory modelled with Prony series only for a narrow band of the whole spectrum, based on the relaxation time(s) employed. Therefore, the relaxation times entering Prony series have to be regarded as design parameters, to be chosen based on the loading history experienced by the viscoelastic material, rather than material parameters. With the values that were chosen here, an accurate estimation of the material response can be expected, at most, over two orders of magnitude, centred on a frequency of 1 Hz.

**Figure 3.** Relaxation modulus in time and frequency domain.
