A Multi-Scale Investigation to Predict the Dynamic Instabilities Induced by Frictional Contact

Abstract

We propose a new variational formulation to model and predict friction-induced vibrations. The multi-scale computational framework exploits the results of (i) the roughness measurements and (ii) the micro-scale contact simulations, using the boundary element method, to enrich the contact zone of the macroscopic finite element model of rubbing systems with nominally flat contact boundaries. The resulting finite elements at the contact interface of the macroscopic model include (i) a modified normal gap and (ii) a micro-scale description of the contact law (i.e., pressure gap) derived by solving the frictionless contact problem on a rough surface indenting a rigid half-plane. The method is applied to a disc brake system to show its robustness in comparison with classical deterministic formulations. With respect to the traditional complex eigenvalues analysis, the proposed multi-scale approach shows that the inclusion of roughness significantly improves the results at low frequencies. In this panorama, any improvement of dynamic instabilities predictions should be based on an uncertainty analysis incorporating roughness combined with other parameters such as friction coefficient and shear moduli of the pads, rather than on roughness itself.

1. Introduction

The frictional contact problem plays a major role in engineering applications, whether in studying the fretting fatigue in dovetail blade roots in aeronautics [1], valve systems for nuclear cooling [2] or braking systems in the automotive industry [3,4,5]. The mechanical behavior of frictional systems mentioned above is micro-scale dependent. Their performance is conditioned by numerous microscopic phenomena, such as roughness, progressive damage of asperities as well as resulting debris [6,7]. For instance, it was shown in [8] that the macroscopic behavior of the braking system is strongly influenced by the microscopic contact properties, which involves the roughness of the sliding interfaces. Hence, it is important to review and to overcome the assumption of perfectly flat surfaces by taking into account the statistical characteristics of rough surfaces.

In general, the modeling of rough surfaces requires the knowledge of two functions: (i) the height distribution function (HDF) [9] and (ii) the spatial function, also known by the autocorrelation function (ACF) or its Fourier transform, which is called the power spectral density (PSD) function [10]. The first one characterizes the roughness through the average height parameters [11,12], while the second describes the spacial arrangement and the variation of the asperities, which often depicts a fractal behavior [13,14]. Most rough surfaces have a self-affine property, which means that the profile remains similar under different magnification [13,15]. The latter can be decoded by a PSD having a power-law shape as a function of the wavenumber (or wavelength) with a potential plateau characterized by a longer wavelength, namely the long-distance roll-off wavevector [16,17,18].

With regard to the mechanical contact models, the first contribution goes back a century with the widely used Hertz theory [19,20]. This theory assumes perfect smooth contacting surfaces, and hence, it is not valid for rough surfaces since it failed to give a reasonable estimation for the real contact area and for the pressure-gap behavior. The strong deviation from the perfect contact assumption led researchers to develop more realistic models to take into account the roughness. The first asperity-based models, pioneered by research work [9,21], are based on a statistical method and the random process theory introduced for the first time in [10,22]. The concept of their approach is to transform the problem of two contacting rough surfaces to a deformable rough one in contact with a rigid plane. Then, the roughness is modeled by a number of isolated spherical asperities with the same radius, and their height is described randomly following a given probability density function. The asperity-based models provide a good approximate solution for the evolution of the contact area especially in the case of a small load [23]. However, their main drawback is the omission of interactions between the asperities, which can lead to an overestimation of the contact pressure. With the development of advanced analytical models [24,25,26,27], interactions between asperities can be taken into account. The coalescence of asperities is another effect that has to be considered for a realistic representation. In fact, the contact spots do not grow independently in reality. They merge and, therefore, contact patches develop. Afferante et al. [28] proposed a suitable approach to address the coalescence of two contacting spots. The concept of their approach is that the asperities with overlapping contacts spots are eliminated and corrected with a single contact patch, namely the equivalent asperity. It should be noted that the contact area of the new asperity includes the area of the suppressed contact spots. In parallel to the aforementioned analytical models, Person proposed an ingenious theory [13,14,29] to tackle the contact problem for rough surfaces. The fundamental concept of the fractal approach is to solve the contact problem at different scale by introducing the evolution of contact pressure probability density. The author mentioned in [14] that the latter satisfies a diffusion-like equation with an appropriate boundary conditions. Persson’s approach has been criticized on the pretext that the proof of the diffusion-like equation is not rigorous, since it is derived assuming full contact and used with the boundary condition on its solution to model the partial contact. In general, it provides a solution that is quantitatively inconsistent with the available solutions, which leads to somewhat smaller contact areas and, thus, greater pressures than those found in the literature [30,31].

The development of computing facilities has enabled the emergence of accurate numerical models that, unlike previous analytical models, are free of assumptions. A couple of numerical techniques can be distinguished: (i) the finite element method (FEM) [2,32], (ii) the boundary element method (BEM) [33] and (iii) Green’s function molecular dynamics (GFMD) [34,35]. The first class of methods is the most used and the most accurate technique. It is based on the optimization of the variational formulation of the contact problem [36,37]. In general, it is applied to a representative volume element (RVE) of the rough surface and requires intensive computing, since the element size should be small enough to capture the asperity waviness and the roughness [2]. The last two classes are more attractive for modeling rough contacts. This is essentially due to the fact that only the rough surface needs to be discretized, and not the bulk, as required by FEM, which allows increasing the mesh density and, hence, perform more accurate studies. Despite these advances, BEM and GFMD are based on the fundamental elasticity theory. The generalization of these approaches to take into account non-linearities is sometimes possible but is not an easy task. For a detailed review, see [32,38].

In the industrial context, FEM is mainly used to model and predict the behavior of frictional systems. In most applications, the contact interfaces are assumed to be flat and have dimensions generally much larger than the microscopic scale of the roughness. Therefore, explicit integration of the roughness in these interfaces is not possible, since it generates extremely dense meshes which, in turn, will increase the computation time. In consequence, authors in literature abandon the classical path of explicitly modeling surface roughness by proposing multi-scale embedded strategies to enhance the computational cost of the direct FEM. For instance, the authors in [39,40,41] suggested integrating the micro-mechanical contact laws in the sliding interfaces. The fundamental concept behind this approach is to derive and then integrate the constitutive contact equations using the homogenization process, which describes the microscopic behavior within the contact zone. Bonari et al. [42] also developed a novel multi-scale formulation that takes advantage of micro-scale simulations to enrich the contact interfaces of macroscopic FE models. The developed method allows taking into account any desirable topography to model efficiently and within a reasonable CPU time the frictional system. In the same spirit, Paggi et Reinoso [43] introduced a new idea that takes into account the shape of the roughness to correct the normal separation in macroscopic contact interfaces. Tison et al. [3,44] investigated the use of rough contact to model the dynamic behavior of a disc brake system. They used a multi-scale framework based on the random field theory [45,46]. The integration of the latter within the complex eigenvalues analysis (CEA) provides a remarkable prediction, which allowed for better correlation with the experimental results.

Regarding the dynamic instabilities that represent the core of this paper, it is well accepted that the frictional contact influences drastically the global dynamic behavior of a rubbing system [47,48]. Taking into account the contribution of the frictional contact to the system’s overall stiffness matrix, dynamic instabilities are predicted using two methods: (i) transient analysis and (ii) CEA. The latter is the most widely used in the industry, as it offers a very good compromise between computation time and accuracy. Indeed, CEA approach uses the mode coupling theory to evaluate the stability of the rubbing system. Thus, it will be unstable if and only if at least one predicted vibratory mode has a strictly positive real part [49]. Since the frictional contact problem is closely linked to dynamic instabilities, the contribution of the roughness should be taken into account and investigated in detail for a potential enhancement of the deterministic CEA. As pointed above, the classic FEM is time consuming since the mesh should be fine enough to capture all the physics of the rough surface. A novel strategy has therefore to be found to fulfill two main criteria: (i) an accurate prediction of the dynamic instabilities and (ii) a reasonable computational cost. In this light, the main contribution of this paper is to suggest a simple alternative approach, highly flexible, that allows accurately and reliably predicting the dynamic instabilities of an industrial braking system. The proposed multi-scale finite element formulation, presented in detail in Section 2, abandons the traditional path of explicitly introducing roughness by means of the asperities in the contacting interfaces. In fact, it considers that the contact interfaces of a large-scale FE model are perfectly flat and smooth. But each element (contact patch) is enriched by adding two contributions, namely (i) the pressure-gap contact law and (ii) the contact element activation. The first one is derived using BEM on a generated rough surface and is assigned to each patch. The use of micro-scale contact law, for each patch, will modify the penalty method and will introduce the roughness in an implicit manner. The second enhancement is introduced within the normal gap function. Indeed, a critical gap is added to the classic formulation. It is derived from HDF resulting from the topography measurement. The embedded strategy is then applied to a disc brake system in Section 3 to investigate the prediction of dynamic instabilities.

2. Governing Equation

The first section introduces the basic governing equations for the non-linear contact problem. After recalling the general framework of the boundary value problems with constraints, the variational formulation of the weak problem will be presented for the case of frictionless contact problems of two deformable solids. This is the first step of a complex eigenvalues analysis (CEA) where the equilibrium position, $u_e$, is determined. Finally, we will present the principle of our approach, which consists of integrating roughness in macroscopic FE models in order to predict dynamic instabilities of a frictional system in a more realistic way.

2.1. General Framework

In the framework of continuum mechanics, a solid is modeled by an open domain $\mathrm{\Omega }$, which is assumed to be bounded and regular. The boundary of $\mathrm{\Omega }$ is divided in two parts: (i) the first one denoted $\partial {\mathrm{\Omega }}_u$ where the Dirichlet boundary conditions are applied and (ii) the second one denoted $\partial {\mathrm{\Omega }}_f$ where the Neumann boundary conditions are applied. So, without any loss of generality, the formulation of the problem when two deformable bodies ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ (see Figure 1), with $\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2$, come in contact at a single contact zone $\partial {\mathrm{\Omega }}_c$ is given by a system of elliptical partial differential equations that can be categorized as follows:

• static balance of momentum equation,

  $$\left\{\begin{array}{c}\begin{array}{cccc}\nabla ·\underline{\underline{\mathit{\sigma }}}& =& 0\hfill & \mathrm{in}\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2\hfill \\ \underline{\underline{\mathit{\sigma }}}.\underline{n}& =& {\underline{\sigma }}_0\hfill & \mathrm{at}\partial {\mathrm{\Omega }}_f\hfill \end{array},\hfill \end{array}\right.$$

  (1)
• compatibility equations,

  $$\left\{\begin{array}{c}\begin{array}{cccc}\underline{\underline{\mathit{\epsilon }}}\left(\underline{u}\right)& =& \frac{1}{2}(\nabla ·\underline{u}+\nabla ·{\underline{u}}^T)\hfill & \mathrm{in}\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2\hfill \\ \underline{u}& =& {\underline{u}}_0\hfill & \mathrm{at}\partial {\mathrm{\Omega }}_u\hfill \end{array},\hfill \end{array}\right.$$

  (2)
• constitutive relation,

  $$\underline{\underline{\mathit{\sigma }}}={}^4\underline{\underline{\mathbb{C}}}:\underline{\underline{\mathit{\epsilon }}}\left(\underline{u}\right)\mathrm{in}\mathrm{\Omega }={\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2,$$

  (3)
  where $\underline{\underline{\mathit{\sigma }}}$ is the Cauchy stress tensor and $\underline{\underline{\mathit{\epsilon }}}$ is the so-called small strain tensor, which is related to the displacement field $\underline{u}$. ${\underline{\sigma }}_0$ and ${\underline{u}}_0$ are the stress and the displacement field imposed on the boundaries $\partial {\mathrm{\Omega }}_f$ and $\partial {\mathrm{\Omega }}_u$, respectively. ${}^4\underline{\underline{\mathbb{C}}}$ is the Hook fourth-order tensor for the considered domain.

In the case of a frictionless contact between two bodies, a relevant complimentary conditions should be formulated. The latter defines the geometrical and the mechanical state of the two contacting surfaces. The complimentary conditions are called Hertz–Signorini–Moreau conditions and can be written as

$$g\ge 0,{\sigma }_n\le 0,g.{\sigma }_n=0,{\underline{\sigma }}_t=0\mathrm{at}\partial {\mathrm{\Omega }}_c={\mathrm{\Gamma }}_c^1\cup {\mathrm{\Gamma }}_c^2,$$

(4)

where ${\sigma }_n$ and ${\underline{\sigma }}_t$ refer to the normal and the tangential contact stresses, respectively. g refers to the gap function between the slave and master surfaces. Hence, the Hertz–Signorini–Moreau conditions can be read as the non-penetration–non-adhesion conditions. In other words, if the two bodies ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ are in contact, then the gap function g is equal to 0 and the normal contact stress ${\sigma }_n$ is below 0; otherwise, $g>0$ and ${\sigma }_n=0$.

2.2. Variation Formulation of the Frictionless Contact Problem

The statement of the weak formulation for the contact problem is obtained by multiplying the strong form in Equation (1) by a test function $\underline{v}$ and integrating over the domain $\mathrm{\Omega }$, and we write

$${\int }_{\mathrm{\Omega }}\nabla ·\underline{\underline{\mathit{\sigma }}}·\underline{v}d\mathrm{\Omega }=0\forall \underline{v}\in \mathbb{V}.$$

(5)

By performing integration and applying the divergence theorem, Equation (5) becomes

$${\int }_{\partial \mathrm{\Omega }}\underline{n}·\underline{\underline{\mathit{\sigma }}}·\underline{v}d\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}\underline{\underline{\mathit{\sigma }}}:\nabla \underline{v}d\mathrm{\Omega }=0,$$

(6)

where: is the contraction operation of two second-order tensors and $\underline{n}$ denotes the outward normal at the slave surface ${\mathrm{\Gamma }}_c^1$. We have $\underline{n}=-\underline{\nu }$ where $\underline{\nu }$ is the outward normal at master surface ${\mathrm{\Gamma }}_c^2$ (see Figure 1).

Note that by introducing the traction vector ${\underline{\sigma }}_0$, the contact boundary conditions and the virtual displacement $\delta \underline{u}$ in the test function $\underline{v}$, the first term in Equation (6) may be rewritten as the sum of two main parts: (i) the frictionless contact contribution and (ii) the Neumann boundary condition contribution. Symbolically, we write

$${\int }_{\partial {\mathrm{\Omega }}_c}\underline{n}·\underline{\underline{\mathit{\sigma }}}·\delta (\underline{\rho }-\underline{r})d\mathrm{\Gamma }+{\int }_{\partial {\mathrm{\Omega }}_f}{\underline{\sigma }}_0·\delta \underline{u}d\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}\underline{\underline{\mathit{\sigma }}}:\delta \nabla \underline{u}d\mathrm{\Omega }=0,$$

(7)

where $\underline{\rho }$ and $\underline{r}$ are the displacement vectors of master and slave points, respectively. Furthermore, the quantity $\underline{r}-\underline{\rho }$ is the gap vector describing the position of the slave point $\underline{r}$ and his projection into the master surface $\underline{\rho }$.

In the case of a normal gap, $g_n\underline{n}=\underline{r}-\underline{\rho }$, the variational form in Equation (7) may be written as follows:

$${\int }_{\partial {\mathrm{\Gamma }}_c^1}{\sigma }_n\delta g_{n}d\mathrm{\Gamma }-{\int }_{\partial {\mathrm{\Omega }}_f}{\underline{\sigma }}_0·\delta \underline{u}d\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}\underline{\underline{\mathit{\sigma }}}:\delta \nabla \underline{u}d\mathrm{\Omega }=0.$$

(8)

It should be noted that the mathematical constraints of the strong and weak forms are not the same. For the strong form in Equation (1), the Cauchy tensor is required to be smooth enough; i.e., $\underline{\underline{\sigma }}\in {\mathcal{C}}^1\left(\mathrm{\Omega }\right)$, which is not the case for the weak form in Equation (6). The order of differentiation in the weak formulation is lower than the strong one. In fact, the order of differentiability in the weak integral form is distributed between the Cauchy stress tensor and the test function.

Finally, the balance of virtual work, called also the weak form, for the frictionless contact problem can be formulated as follows:

Find $\underline{u}\in \mathbb{U}=\left\{\underline{u}\in {\mathbb{H}}^1\left(\mathrm{\Omega }\right)\mid \underline{u}={\underline{u}}_0\mathrm{on}\partial {\mathrm{\Omega }}_u\right\}$ such as

$$\left\{\begin{array}{c}\begin{array}{c}{\displaystyle {\int }_{{\mathrm{\Gamma }}_c^1}{\sigma }_n\delta g_{n}d\mathrm{\Gamma }-{\int }_{\partial {\mathrm{\Omega }}_f}{\underline{\sigma }}_0·\delta \underline{u}d\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}\underline{\underline{\mathit{\sigma }}}:\delta \nabla \underline{u}d\mathrm{\Omega }=0}\\ \mathbb{V}=\left\{\delta \underline{u}\in {\mathbb{H}}^1\left(\mathrm{\Omega }\right)\mid \delta \underline{u}=0\mathrm{on}\partial {\mathrm{\Omega }}_u\right\}\\ \mathbb{C}=\left\{\delta \underline{u}\in \mathbb{V}\mid (\underline{r}+\delta \underline{r}-\underline{\rho }+\delta \underline{\rho })·\underline{n}\ge -g_{n0}\right\}\end{array},\hfill \end{array}\right.$$

(9)

where ${\mathbb{H}}^1\left(\mathrm{\Omega }\right)$ denotes the Hilbert space of the first order, $\delta \underline{r}$ and $\delta \underline{\rho }$ refer to $\delta \underline{u}$ but in the contact interface. And, finally, $g_{n0}$ represents the initial gap between contacting surfaces.

2.3. Overview of the Embedded Computational Strategy to Include the Roughness on Macro-Scale Model

Under the framework of the finite element method, the formulation in Equation (9) will be solved using the most popular scheme: the penalty method. The goal as mentioned above is to find the equilibrium position due to the frictionless contact (see Section 1 in [4] for more detail).

To approximate the Hertz–Signorini–Moreau condition, the normal contact stress arising at the contact interface will be expressed as a function of the gap using the following approximation:

$${\sigma }_n\left(g\right)={\epsilon }_n(<-g>)=\left\{\begin{array}{cc}0\hfill & g>0,\mathrm{there}\mathrm{is}\mathrm{no}\mathrm{contact}\hfill \\ {\epsilon }_n(-g)\hfill & g\le 0,\mathrm{there}\mathrm{is}\mathrm{contact}\hfill \end{array}\right.,$$

(10)

where ${\epsilon }_n$ is a non-positive, continuous and strictly decreasing function defining the evolution of the contact stresses, ${\sigma }_n$, as a function of the gap, g, between the master and the slave surfaces. In the case of the linear penalty method, the function ${\epsilon }_n$ will simply denote the classic contact stiffness (i.e., penalty coefficient). Finally, $<·>$ refers to $max\{·,0\}$.

The penalty method can be seen as an approximation of the contact constraints. It leads to a small penetration between the slave and master surfaces (see the case when $g<0$ in Equation (10)). It should be noted that the definition in Equation (4) will be met if the contact stress is higher for a small penetration. In other words, the penalty method does not restrict the penetration between the contacting surfaces but it resists it. If it is deeper (i.e., $g<0$), the value of ${\epsilon }_n(-g)$ will be higher and, hence, the real contact stress (i.e., reaction) will appear.

Under the last assumption, the virtual work in Equation (9) can be divided into the classic solid mechanics part, $\delta W_s$, and the contribution from the contact problem, namely $\delta W_c$. The weak formulation of the system is obtained using the penalty method. We have

$$\delta W_c+\delta W_s=0,$$

(11)

where

$$\delta W_c={\int }_{{\mathrm{\Gamma }}_c^1}{\epsilon }_n(<-g_n>)\delta g_{n}d\mathrm{\Gamma },$$

(12)

and

$$\delta W_s={\int }_{\mathrm{\Omega }}\underline{\underline{\mathit{\sigma }}}:\delta \nabla \underline{u}d\mathrm{\Omega }-{\int }_{\partial {\mathrm{\Omega }}_f}{\underline{\sigma }}_0·\delta \underline{u}d\mathrm{\Gamma }.$$

(13)

In this paper, the the potential energy of the contact interaction will be formulated as follows:

$${\mathrm{\Pi }}_c=\frac{1}{2}{\int }_{{\mathrm{\Gamma }}_c^1}{\epsilon }_n{(<-(g_n-g_n^{\mathrm{c}\mathrm{r}})>)}^{2}d\mathrm{\Gamma },$$

(14)

where $g_n^{\mathrm{c}\mathrm{r}}$ is called the critic gap. It is a value that will be defined for each contact element and from which the latter will be activated. The first enrichment introduced in Equation (14) will allow for a non-homogeneous activation of the contact status for all nodes of the contact interface.

Considering Node-To-Surface (N2S) discretization, the slave surface ${\mathrm{\Gamma }}_c^1$ can be presented by nodes denoted by ${\underline{r}}^i$ and the master one, ${\mathrm{\Gamma }}_c^2$, can be discretized as set of segments denoted by ${\mathrm{\Gamma }}_{c,j}^2$. Let us assume that we have N slave points and M master segments. One can observe from Figure 2 that for all slave nodes ${\underline{r}}^i$, one or more master segments ${\mathrm{\Gamma }}_{c,j}^2$ can be determined using the normal projection. Hence, the i-th contact element can be defined by a combination of one slave node ${\underline{r}}^i$ with its correspondent master segment ${\mathrm{\Gamma }}_{c,j}^2$. By using the last definition and the fact that we have $n_c$ contact finite elements, the discretized form of Equation (14) can be rewritten as follows:

$$\left\{\begin{array}{c}{\mathrm{\Pi }}_c\approx \frac{1}{2}\sum _{i=1}^{n_c}{\int }_{{\mathrm{\Gamma }}_{c,i}^2}{\epsilon }_n^i{(<-(g_n^i-g_n^{cr,i})>)}^{2}d\mathrm{\Gamma }\hfill \\ {\epsilon }_n^i\ne {\epsilon }_n^j\forall i\ne j\hfill \end{array}\right.,$$

(15)

where ${\epsilon }_n^i$ denotes the contact law for the i-th contact element. It pairs up the normal contact stress with the gap distance. Note that the contact law from one element to another is not the same. Each contact element will have its own law constituting the second enrichment.

The idea behind the proposed strategy is to integrate the heterogeneity observed in the asperity scale into macroscopic models. In general, the present methodology incorporates three main steps: (i) The first one consists in characterizing the rough surface (the objective of Section 3.1). (ii) The second one aims to solve, at the roughness scale, the frictionless contact problem. This step is important as it allows deriving the pressure-gap’s contact law by taking into account the roughness. Finally, (iii) the third one consists in introducing the height distribution as well as the local behavior of the contact stiffness to enrich the gap, equivalently the contact detection, and the penalty function. The overall procedure is summarized in Figure 3.

In the first step (Figure 3a), the rough surface is characterized for example by using a confocal laser scanning microscopy (CLSM) technique in order to measure the height z of each point. Since the roughness is random [10], the height of the rough surface is considered, in this work, as a random variable (RV) with an independent Cartesian coordinates. The measured realization of the RV will be used to estimate the HDF denoted by $P\left(z\right)$. The latter will be introduced in the third step into the macroscopic FE model through the gap function. Practically, the statistics of the rough surface will be transmitted to the FE model using the new form of the normal gap in Equation (14) or (15). Hence, we write

$$P\left(z\right)\equiv P(\Delta ),$$

(16)

where $\Delta =g_n-g_n^{\mathrm{c}\mathrm{r}}$ is a RV. The realization ${\Delta }_i^j$ of $\Delta$ will be assigned for each contact element i.

In a second step (Figure 3b), several rough surfaces will be generated numerically based on the statistics and the spectral analysis of the measured surface. The objective is to apply the BEM in order to solve microscopically the contact problem of the rough surfaces. The purpose of this micro-scale simulation is to compute the micro-mechanical behavior of the contact interface involving the roughness. The latter is modeled by the pressure-gap law for each generated rough surface and injected into the weak form through the function ${\epsilon }_n$ of Equation (10). It is agreed that this requires intensive computation on several rough surfaces, especially if advanced BEM solvers are used to take into account interfacial or material non-linearities such as adhesion or plasticity. Here, the crucial goal is to derive a set of local contact laws, by taking into account the roughness, for each contact element (patch).

Finally, in the last step (Figure 3c), the contact laws as well as the enriched gap function will be assigned to each patch. The integration procedure is the same for both parameters. Indeed, assuming that the contact interface is smooth and flat, the sampled results resulting from the BEM simulations and the measured surface characterization will be embedded within the contact elements. Following the logic of the embedded approach, the roughness will not be represented explicitly in the FE model. It will be included implicitly through the new formulation of the contact gap and the modified penalty function. Next, the enriched FE model can be solved using CEA to compute the dynamic instabilities. In general, the procedure falls into a fully parallelizable loop. After the surface characterization, steps 2 and 3 (Figure 3b,c) can be repeated as many times as needed in order to obtain a representative statistical prediction.

3. Application to a Disc Brake System

In this section, the previously proposed computational strategy is implemented to predict the dynamic behavior of an automotive disc brake system. First, the pad surface is measured by a CLSM. The objective is to characterize the pad surface topography by means of (i) the height distribution function (HDF) and (ii) the power spectral density (PSD) function. Second, a batch of artificial surfaces will be randomly generated based on the previously measured functions (HDF and PSD). The aim is to apply BEM, developed in [50], to predict the contact laws, i.e., the evolution of the load as a function of the gap. These elements will be integrated, toward the end, into a macroscopic FE model of a disc brake to reproduce the heterogeneity observed in the contact interfaces in order to predict the dynamic instabilities induced by friction.

3.1. Pad Surface Characterization

The pad’s surface is measured experimentally through a CLSM. Figure 4 shows an example of the measured topography over a scan area of $1.8\times 1.4$ mm${}^2$ with a resolution of $256\times 256$ pixels. The measured topography (see Figure 4b) appears irregular and depicts characteristics of randomness. Hence, the random process theory, widely applied to analyze roughness, can be used to model accurately the probability function of the surface heights also known by the HDF. The latter holds only the out-of-plan information. To complete the characterization, the PSD is used to describe the spatial arrangement in the plane.

Figure 5a illustrates the computed 2D-PSD of the pad’s surface using the developments in our previous research work [18,38]. It shows clearly that the pad’s roughness is highly isotropic. Such a result allows computing the radial PSD (averaged) in order to study the effect of resolution and scan length. Moreover, it can be seen, in Figure 5b, that the behavior of the measured PSD is not sensitive to resolution or scan length. It conserves its shape despite the change of the measurement parameters. With regard to the PSD behavior, Figure 5b depicts roughly two linear regions in the log–log plot. It starts from the lower frequency, which is inversely correlated to the scan length $q_L=\frac{2\pi }{L}$, up to a high measured frequency, related to the short-distance cut-off wavevector, defined by $q_s=\frac{2\pi }{\Delta }$ where $\Delta$ refers to the sampling length. In fact, the radial PSD has a shape similar to the bi-fractal surfaces defined in [51]. The difference lies in the definition of the slopes of the linear regions. In fact, bi-fractal surfaces have a theoretical PSD form that can be defined as follows:

$$\mathrm{PSD}\left(q\right)=\left\{\begin{array}{cc}{C}_0{\left(\frac{q}{q_L}\right)}^{c_1}\hfill & \mathrm{if}{q}_L<q<q_c\hfill \\ C_1{\left(\frac{q}{q_c}\right)}^{c_2}\hfill & \mathrm{if}{q}_c\le q<q_s\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.,$$

(17)

where the exponents $c_i$, $i\in \{1,2\}$ (the slope of the linear regions in log–log plot) are related to the Hurst exponent $H_i$ as follows: $c_i=-2(H_i+1)$ with $0<H_i<1$. In our case, the two regions have a slope of −2 and −4.5, which correspondents to Hurst exponents of $H_1=0$ and $H_2=1.25$, respectively. The first region defines a borderline case of fractality, while the second is completely outside the fractal or the self-affine framework. However, the model defined in Equation (17) will be used to characterize the PSD of the pad’s surface. The parameters $c_1$, $c_2$, $q_L$, $q_c$, $C_0$ and $C_1$ are computed using a suitable optimization scheme on the radially averaged PSD exactly as defined in Equation (22) of [18].

Figure 6 complements the information given by the PSD. It presents the HDF of the measured pad’s surface in a normalized manner and compares it with the centered and reduced normal law. As a first observation, one can see that the surface is approximately Gaussian. A notable difference is observed for the higher heights which makes the HDF asymmetric. This behavior was expected since automotive disc brake pads are usually compacted and ground to correct height. From Figure 6, the roughness of the pad can be defined by means of several parameters. For example, the mean (or median) plan $\overline{z}=0.0256$ mm, the root mean square (rms) roughness $\sigma =0.0331$ mm, the skewness $s=-$0.2410 and the kurtosis $\kappa =3.6549$.

Note that in the following, the fitted PSD and HDF will be taken into account to numerically generate rough surfaces having the same spatial arrangement and roughness as that measured experimentally. The numerical procedure to generate artificial rough surface is inspired by Wu’s algorithm [52]. The latter is adapted to take into account PSDs with shapes similar to what has been defined in Equation (17).

3.2. Contact Law of Pad/Disc Interfaces

In the previous section, the surface of the pad has been fully described by two quantities: (i) PSD and (ii) HDF. These will act as an input data to artificially generate a batch of random rough surfaces. The artificial surfaces are intended to mimic the roughness behavior of the automotive disc brake pads so that they can be used to compute contact laws, i.e., load-separation curves. An example of an artificial rough surface is presented in Figure 7a. It has the same spectral and statistical specifications as the measured pad’s surface, namely (i) PSD properties (compare Figure 7b with Figure 5a) and (ii) the roughness parameters (the mean plan and the rms roughness).

After generating batch of artificial random surfaces, the next step is to compute the evolution of contact stresses against contact kinematics for each generated pad’s topography. At this stage, many models can be used. The study of the latter has been addressed in our previous research [38], where authors tried to explain the philosophy behind each theory with a comparative study. The concern here is to apply directly the result of [38] in order to solve frictionless rough contact and thus, obtain the evolution of contact stresses against the surface separation, commonly known as the gap. More precisely, BEM solver, introduced by Bemporad and Paggi [50], will be used to solve the contact problem between the artificial rough surfaces and the rigid surface of the disc. The choice of using BEM is motivated by the fact that it offers the best trade-off between FEM and semi-analytical models. Indeed, it is free of any kind of assumption unlike semi-analytical models, and it allows discretizing only the rough surface without the bulk, which can save considerably the CPU cost.

As explained above, the curve load-separation for each generated rough surface is obtained using BEM code developed in [50]. BEM simulation is performed on a generated rough surface of $2.825\times 2.825$ mm${}^2$ with a resolution of 128 pixels in both directions. This numerical setup is chosen because (i) BEM solvers show strong convergence on meshes containing at least $128\times 128$ elements [38] and (ii) the single finite element (patch) of the pad’s surface has a averaged area equal to $2.825\times 2.825$ mm${}^2$ (see the element j defined in Figure 3c). Since the pad’s material involves an orthotropic behavior, only Young’s modulus, $E_z=1983.3$ MPa, in the most stressed direction, in this case the z-direction, and the Poisson’s ratio ${\nu }_{xz}=0.4054$ are taken into account as inputs for BEM simulations. It should be noted that these micro-scale simulations assume a homogeneous and linearly elastic rough surface. Therefore, geometric (i.e., large deformations), material (i.e., plasticity) and contact interface non-linearities (i.e., adhesion and friction) are not considered in this computational framework.

Figure 8 shows BEM results for 100 generated rough surfaces. It depicts the evolution of the normalized contact pressure, $\frac{P}{E^{\ast }}$ (where $E^{\ast }$ is the composite modulus) versus the separation. The latter is defined as a distance between the imposed displacement of the rigid surface (i.e., the surface of the disc) and the median plan, $\overline{z}$, of the rough surface. Each simulation, performed on a generated rough surface of $2.825\times 2.825$ mm${}^2$, demonstrates a non-linear mechanical behavior of the considered set of asperities. For a small displacement, the separation between the two surfaces is large, which implies that only few asperities will be in contact, thus generating a low contact pressure. As long as the imposed displacement increases (synonymous to a decrease of the separation between the rigid surface and the midplane), the real contact area evolves, leading to an increase in contact pressure and, hence, the contact stiffness. This behavior, i.e., pressure-gap law, will be integrated into each patch (i.e., finite element) of the pad’s surface of the macroscopic FE model. Following this logic and according to the strategy defined in Figure 3, the macroscopic FE model of a disc brake system will have flat and smooth contact interfaces including the finite element patches. It should be noted that each patch will be driven by a contact law (BEM solution on single generated rough surface), instead of the classic penalty coefficient. Moreover, the gap of each patch is enriched by a threshold according to Equation (16). The new gap formulation will benefit from the HDF of Figure 6 to add a contribution (i.e., threshold) modeling the randomness of the measured height z. This enrichment is intended to ensure the activation of the contact elements in a non-uniform way just like what happens in surfaces with low-scale asperities.

3.3. Application to Dynamic Instabilities Prediction

The multi-scale approach is implemented to predict dynamic instabilities of a disc brake system model composed from three main structural subsystems (see Figure 3c): two pads and a rotating disc. The pad subsystem is broken down into three elements, namely the lining material, the backplate (lining support) and the shim that tends to reduce vibrations. The braking action is represented by a 5 bar pressure applied to the pads on a circular surface. This pressure replaces the action of the hydraulic system which tends to move the piston to squeeze pads against the rotating disc. The current FE model of the braking system has been used to conduct several numerical tests in our previous work [4,5].

The applied multi-scale analysis is divided into four major steps that are part of parallelized loop. After the pad’s topography characterization, the loop is started by:

1.

Artificial rough surface generation: The purpose of this step is to use the characterized HDF and PSD to generate a batch of artificial rough surfaces similar to the measured topography.

2.

Micro-scale contact simulations: Here, BEM solver introduced in [38,50] is used to solve the elastic contact problem considering artificial surface asperities and a rigid half-plane. The resolution is performed on the whole batch of the generated rough surfaces. Solving the micro-scale contact problem minimizes the convex quadratic program (QP) (see Equation (50) in [38]). The obtained results are the contact load (or the contact pressure) and the separation between the two contacting bodies.

3.

Enrichment of the contact element: The objective of the third step is to assign each micro-scale contact law (obtained in the second step) to each flat patch (i.e., contact element). Moreover, the gap between slave and master nodes of the macro-scale FE model (disc brake system FE model) is modified by adding a threshold to each gap in order to activate the contact elements in a non-uniform manner. As mentioned above, the added threshold depends closely on the measured HDF.

4.

Complex eigenvalues analysis: At this stage, the traditional CEA is performed. The beginning of the last step starts with a quasi-static analysis. Its goal is to solve, progressively as the load increases, the frictional contact using the enriched contact elements. At the end of this step, the well-known complex modal analysis is performed to compute both complex eigenvalues and eigenvectors and identify the unstable modes (those with negative damping). For more details, see Section 2.1 in [4].

In the studied example, 300 iterations are used to compute the dynamic instabilities of the braking system.

The result of the multi-scale approach is presented in Figure 9 (blue marker) where the negative damping, a ratio between the real and the imaginary part of the complex eigenvalues, is plotted as a function of the frequency. This figure shows only dynamic instabilities, which means unstable modes with a positive real part (i.e., negative damping below zero). In order to compare the robustness and limitations of the multi-scale approach, the result of FEM considering perfect contact with a pad/disc friction coefficient of $0.5$ (red square marker) as well the results of a stochastic model based on a FAST-FEM solver (black marker) developed by Maaboudallah et al. [4] are added in the same figure. It is known [3,44] that stochastic modeling leads to better correlation with experimental data. At first sight, it can be seen that the introduced contact elements enrichment affects the dynamic instability predictions. In particular, the unstable modes predicted by FAST-FEM in the low frequency range are reproduced by the multi-scale approach. However, the latter fails to reproduce the instabilities predicted by FAST-FE or even by standard finite elements with perfect contact at high frequencies. Despite these observations, some unstable family modes around 13 and 15 kHz are roughly highlighted by the multi-scale approach. This unconventional behavior tends to demonstrate that the roughness affects mainly low frequencies. In other words, the addition of the roughness contribution, under the framework mentioned above, can be seen as a necessary but not sufficient condition. Indeed, a robust prediction of dynamic instabilities requires the integration of the roughness but also the randomness of the most sensitive variables like the friction coefficient of pad/disc interface and the normal Young’s modulus of the pad.

4. Conclusions

In this paper, we have presented a multi-scale computational method to predict friction-induced vibrations. The method uses micro-mechanic contact simulations on characterized rough surface to enrich the macro-scale finite element of a rubbing system. The idea behind the proposed approach lies in four essential points:

1.

Roughness characterization using the power spectral density function and the height distribution function;

2.

Micro-scale contact simulations on the characterization roughness using the boundary element method;

3.

Enrichment of the contact finite element using (i) the micro-scale contact laws (obtained from BEM) and (ii) the modified gap functions (obtained from HDF);

4.

Performing the complex eigenvalue analysis on the updated stiffness matrix.

The multi-scale method was applied to a scale model of a braking system. The results show that the proposed approach is more accurate than the classical one using a perfect contact law, based on a comparison with a stochastic model accounting for the uncertainties on several parameters such as friction. In addition, it is found that roughness has an effect on the predictions of dynamic instabilities. Indeed, taking into account the roughness brings out new unstable modes in the low frequencies. However, the effect of the integrated roughness is not clear in the high frequencies, since many stochastically predicted modes were not reproduced by the multi-scale strategy. Results of the investigation demonstrate that roughness is only one element among others that must be taken into account for a robust and accurate prediction of the dynamic instabilities. An uncertainty analysis must be conducted by varying the most important parameters including roughness to correctly model and predict the dynamic instabilities of frictional systems.