Contact modelling of highly heterogeneous friction material for braking applications

. Friction brakes are increasingly undergoing considerable development to improve their durability, efficiency, maintenance costs and environmental impact. Nevertheless, to achieve this, it is necessary to understand the different mechanisms involved in contact that are multi-scale and multi-physical in nature. On the multi-scale aspect, it is well known experimentally that heterogeneities have a pre-weighted role on performance without being able to explain it. Thus, modelling seems to be a good way to better understand the influence of these heterogeneities, provided that we have a multi-scale method to consider them. The objective of this article is therefore to propose a methodology for simulating contact in the presence of heterogeneous materials. The strategy consists in enriching the contact rigidity in terms of behaviour by a method of numerical homogenization. The significant advance of this article lies in the consideration of contact within the technique of numerical homogenization of a heterogeneous material. The strategy is then validated by comparing the mechanical fields between the proposed method and an explicitly meshed case. One of the main contributions of this work is the reduction in computing time compared to the traditional FEM method.


Introduction
Increasingly in the development of materials, manufacturers tend to exploit heterogeneous materials because they must respect a certain number of performance criteria. For example, in the field of transport, friction materials are made up of about twenty components in order to obtain a stable friction coefficient at any temperature level, good durability, controlled wear and suitable thermal behaviour. A major problem in the field of transport, particularly for braking applications, is the generation of particle emission and noise pollution. The latter is a major ecological issue. In order to improve the development of these materials, it is then essential to understand what is happening within a contact interface and in particular the role of heterogeneities. The difficulties lie in the distribution of stresses within the contact, which depends on both the loads applied to the brake system and heterogeneities. This localization of stresses within the interface leads to the formation of wear, tear and third body [1]; strongly coupled to the presence of heterogeneities. Current modelling of brake systems generally considers contact materials to be homogeneous and the interface to be continuous, which is insufficient. From an experimental * e-mail: essoarfa@yahoo.fr point of view, there are two types of heterogeneities within the contact interface: materials (components) and surface heterogeneities (roughness, flatness, etc.).
On surface heterogeneities, a contact model between rough surfaces was first proposed by Greenwood and Williamson [2] and developed by authors such as Nayak [3], Bush et al. [4], Persson [5], Ciavarella et al. [6], etc. These authors focus on surface scale and neglect the system. Such models are often based on Hertz theory for elastic solids [7]. Recently, Waddad [8] proposed an analytical solution to take into account the interactions between the different roughnesses. In this paper, a multi-scale method incorporating these surface defects into a large-scale problem was proposed. This strategy was used to highlight the size effect of roughness on vibration instabilities in a complete brake system [9].
On material heterogeneities, Magnier et al. [10] showed that the presence of porosities changes the natural frequencies of the system. Nevertheless, the technique used by [10] consists in explicitly meshing the porosities, making computation time-consuming and therefore limited to a characteristic length of large porosity. Thus, it is essential to work with homogenization methods that take into account the different components. Homogenization methods are common for heterogeneous materials [11,12] with either numerical or analytical techniques. However, these methods are not suitable when there is contact because once the actual mechanical properties have been calculated and the macroscopic computation is performed, the relocation of macroscopic mechanical fields to microscopic ones is not in conformity with explicit results. As mentioned above, the location of the load, which also depends on heterogeneities, is an important parameter which should be considered in the search for good performance (wear, squeal, etc.). Some studies, such as those developed by Mbodj et al. [13], Peillex et al. [14], studied the mechanical properties of heterogeneous material using a homogenization technique. The majority of these models did not take into account the location of heterogeneities because of no relocation of stresses. A problem would arise if we want to consider tribological phenomena, such as wear and tear, on the scale of heterogeneities. Indeed since wear depends on contact pressure distribution [15], the location of contact heterogeneities is important for computing loss of wear volume. From a local point of view, some authors like Leroux et al. [16] presented a method for solving contact problems when one of the contact bodies contains multiple heterogeneous inclusions. These methods are fast but focus on a given scale. From a numerical point of view, Temizer and P. Wriggers [17] proposed a multi-scale finite element method to integrate the heterogeneity of materials in a contact problem. The above developments assume the material in contact as homogeneous and most of them do not take into account the wear mechanisms.
In this article, using a multi-scale strategy, a model that takes into account a heterogeneous material and integrating this aspect into a simulation of a complete braking system is proposed. The strategy is based on a numerical homogenization technique considering contact under sliding conditions in the determination of effective properties. Then, via a simple example, the multi-scale methodology is validated. The results are compared with those of a simulation in which the heterogeneities are explicitly meshed. One of the main contributions of the strategy proposed in this work is the reduction of computation time in comparison to the traditional FEM method.

Multi-scale approach
This section is devoted for a multiscale methodology leading to contact modelling of highly heterogeneous friction material. The multi-scale strategy shown in Figure 1 is based on 4 steps: -The discretization of the heterogeneous material into a set of patches The heterogeneous material, like the one shown in Figure 1, is discretized into a set of heterogeneous cubical patches. Each patch constitutes our Representative Volume Element(RVE).

-Homogenization with contact
In this work, an interface of the heterogeneous material is in contact with another material, so one of the faces of the patches is supposed to be in contact with a flat disc in order to consider the non-linear effects of contact. The computation of the effective properties of each patch considering these non-linearities requires a new homogenization method which will later be called KUBC-Contact -Macroscale embedding The effective property, of each heterogeneous patch, obtained by the mean of the KUBC-Contact method is enriched at the macroscale level in terms of behaviour and a computation is performed to determine the macroscopic mechanical fields. -Re-localization After the macroscopic computation, the relocation of the mechanicals fields is done using macro contact forces and displacements at the heterogeneous micro patch level.

Contact homogenization: KUBC-CONTACT
Our enrichment method is based on the numerical homogenization used in the work of Lejeunes et al. [18] consisting in modelling RVE and extracting its effective properties. There are many boundary conditions (KUBC, SUBC, PBC etc.). Nevertheless, these boundary conditions have very different apparent properties, which raise the question of the representativeness of the samples, particularly when an interface of this sample is in contact with another material. Especially if we take the conditions at the known limits, the relocation of the stresses does not allow us to take into account the engagement effect of the contact, the presence of the friction coefficient, etc. Thus, we propose to take the KUBC boundary conditions on 5 sides of the RVE and to introduce on the 6th side a boundary condition with contact. In other words, contact constraints will be imposed on the contact side of the patch. This method will later be called KUBC-Contact. We denote V the domain, ∂V the free boundary, ∂V c the constraint boundary and | V | the volume of the patch. The patch is presented in Figure 2. The equation governing the KUBC-Contact method is presented in equations (1) and (2): where x is a point of the volume of the patch,¯ d is the imposed mean strain, and U c is determined by the contact conditions. The subscript 1 and 3 represent the transverse direction and the subscript 2 represents the normal direction. Once the exact micro fields in the patch, σ ij and ij , under the applied boundary conditions, is obtained, the averaged stresses and strains are computed over the patch from (3) and (4) In order to determine the effective modulus, equations (3) and (4) are used. The traction (in the three directions) and shear (in the three dimensions) boundary conditions are applied in the bulk modulus K and shear modulus G computation (isotropic linear elasticity). The bulk modulus K is described in equation (5): Once shear stresses and strains are obtained, G is computed from (6): Using the results of equations (5) and (6), Young modulus (Ē) and Poisson ratio (ν) can be computed from equations (7) and (8) respectively: These properties will be reinjected into the macroscopic calculation. Once the latter has been performed, it is possible to relocate from macro to micro in order to have a state of local mechanical fields. This step will be performed in the following. Macro strain (E XX , E Y Y , E ZZ ) are applied as boundary conditions of the patch for bulk computation. The displacements (U X , U Y , U Z ) field of the heterogeneous patch is computed from equation (9): with X Y Z a vector of nodal coordinates.  Macro shear strain (E XY , E XZ , E Y Z ) are applied as boundary conditions of the heterogeneous patch. The displacements (U X , U Y , U Z ) field are computed from equation (10):  Table 1 and are in a good agreement with the classical KUBC method.

Description of the model
In this subsection, the 3D FEM model is inspired by the experimental pin-on disc system developed in the PhD thesis of Duboc [19]. The numerical model, shown in Figure 3, is composed of a steel disc, a friction pin, two steel pin-housing and a steel thin plate. The outer and inner radiuses of the disc are 107.5 mm and 12.5 mm respectively. A prescribed displacement of 0.2 mm is applied at the extremities of the thin plate which enforces contact between the pin and the disc at a time period of t = 1s. The lateral displacements of the plates extremities are fixed. Moreover, the disc is rotating with a constant velocity of 0.1 rad/s. The contact surface of the pin is flat with 20*20 mm 2 of dimensions. In this work, the friction pin contains one inclusion near the inner side of the leading edge of the contact surface (see Fig. 4 ). The spherical inclusion, of radius 0.35 mm at 0.15 mm from the contact surface, is made of steel. Tie elements are used between inclusion and the pin, the pin and the plate, and the pin-housing and the plate in order to ensure continuity between all these components. Material properties of all the components are reported in Table 2. Figure 4 shows a heterogeneous patch, containing the spherical inclusion placed at the leading edge of the friction pad, which will be used in part of the multi-scale strategy proposed in this paper.
All simulations are performed by the commercial Finite element code ABAQUS with 20 CPU cores. Two types of models are proposed here : -The first, called later "explicit method", is based on the conventional method where the inclusion is explicitly meshed. This will be our reference solution in this work. Contact is frictional and the penalty method is used to solve contact constraints. Friction coefficient is set to 0.3. The complete model is meshed with 155447 hexahedrons and 175859 nodes. -The second model is based on the multi-scale strategy presented in the previous section. The inclusion is not physically present in the complete model but is replaced by the effective module obtained by using the KUBC-Contact method on the heterogeneous patch. The spherical inclusion is centered within the square of 1 mm 3 volume. The effective behaviour of the square part is obtained using the KUBC-Contact method (see Fig. 5). After the macroscopic calculation, a relocation of the mechanical fields is performed. In this relocation step, macro contact forces are reapplied at the microscopic scale.

Results with the explicit method
The results of the explicit model illustrated in Section 4.1 are presented here. Contact between the disc and the pin is our main interest here. Therefore in this work, the contact pressure distribution over the contact surface of the pin is presented. Contact pressure distribution, for heterogeneous result compared to the perfect case in which the pin is supposed to be homogeneous, is shown for time period t = 1s in Figure 6. Contact pressure distribution in the presence of the inclusion is represented in Figure 6a. Facing each other on the Figure 6b is presented a case where the patch inclusion has a property of 3GPa, thus making the friction material completely homogeneous. Figure 6b shows a conventional pressure distribution where the maximum value is reached on the leading edge of the pin and decreases in the sliding direction. At the back of the pin, disbonding is observed. The introduction of an inclusion disrupts the pressure distribution where it is logically concentrated on the area where the spherical inclusion is placed. Logically, there is a spherical pressure distribution. This is due to the form of inclusion. Figure 7 shows a zoom of the contact pressure distribution in the heterogeneous patch area.

KUBC-contact versus explicit method
Three results are presented in this section, considering the multiscale strategy we presented in Figure 1: -The homogenization result based on the KUBC-Contact method Here, the overall property of the heterogeneous material is computed using the KUBC-Contact method presented in Section 2. The computation of the Young modulus and the Poisson ratio of the effective material require the computation of the bulk modulus and shear modulus. Figure 8 shows the KUBC-Contact model configuration.

-The macroscale embedded result
Here, the effective modulus obtained in Table 1 will be used to enrich the macro patch according to Figure 5. The same model configuration as described in Figure 3 is used. The heterogeneous patch containing the spherical inclusion is replaced by a homogeneous   patch of the same volume. This means that inclusion in the model will not be physically present and that this model will be faster than the explicit model. The complete model is meshed with 134971 hexahedrons and 154092 nodes. The results are presented in Figure 9.
Contact pressure distribution is conventional, the maximum contact pressure is reached at the leading edge of the contact surface and decreases with the sliding direction. However, the maximum contact pressure in the patch area is higher due to its higher rigidity in comparison to the rest of the material. Contact pressure evolution, in the patch area, of the enriched model does not have the same evolution as the explicit model (see Fig. 7) because the inclusion is not physically present, but the average contact pressure in the patch area of the explicit and enriched models is identical, see Figure 10.

-The re-localization result
This result is compared to the explicit contact pressure distribution shown in Figure 7. As seen in Figure 10, contact pressure distribution between explicit and integrated results is different. In all homogenization based work's, this relocation step is almost forgotten, but as we see the contact pressure distribution is not well estimated. Although the average pressure is the same, a relocation strategy, based on the KUBC-Contact model and boundary conditions previously explained, is proposed. The macroscopically determined contact pressure is reinjected at the patch scale where heterogeneity is explicitly represented. The  relocation boundary conditions and its contact pressure are shown in Figure 11. Contact pressure distribution between explicit strategy and relocation strategy presented here is shown in Figure 12. Contact pressure distribution obtained via our relocation strategy is in good agreement with the explicit pressure distribution. Figure 12c shows the difference between the two pressure fields. The most significant errors are located on the edges of the patch. At the point where the heterogeneity is located, the pressure difference between the two configurations is almost the same. The proposed strategy therefore seems relevant to address contact problems in the presence of heterogeneities.

Parametric study of the friction coefficient
In this subsection, in order to show the robustness of the multi-scale method proposed in this paper, a parametric study on the friction coefficient is performed. The above results are established for a friction coefficient set at 0.3. Here, we vary it from 0.3 to 1.0 and compare the results between the explicit strategy and our multi-scale strategy presented in Figure 1. Results are shown in Figure 13. For each friction coefficient, results are presented on the same scale. Maximum pressure distribution is achieved on the leading edge of the pin. As the coefficient of friction increases, the maximum pressure also increases due to the greater tilting of the pad. This can be seen in particular by the increasing extent of the area where the pressure is zero.
Between explicit and multi-scale methods, it appears that there is a good agreement of the results, particularly on the central part of the patch, whatever the friction coefficient. The most significant errors (in the order of a few %) are found on the edge of the patch and especially on the leading edge (Fig. 13). This parametric study then makes it possible to consolidate the strategy implemented.

Conclusions
In this work, a strategy leading to contact modelling of heterogeneous materials was presented. Firstly, an explicit reference model using a method where inclusion  is explicitly meshed is presented for the validation of our strategy. The inclusion is placed on the leading edge of the contact surface at 0.15 mm from the contact surface to highlight its effect. The result is quite predictive, because the inclusion is rigid compared to the rest of the pin material. A localization of the contact pressure is observed in the contact area of the patch containing the inclusion.
Secondly, in the context of the multi-scale strategy presented in Figure 1, a new contact-based homogenization method is implemented on the basis of the uniform kinematic boundary conditions method. KUBC-Contact method, which is inspired by the KUBC method to which a contact face is added, has the advantage of taking into account all localizations induced by the presence of heterogeneities in contact conditions. The effective elastic properties obtained are in a good agreement with the conventional KUBC method. The addition of the boundary condition in contact allows a good relocation of the mechanical fields.
Finally, after the macroscopic computation, the relocation, using contact pressure at each node of the contact surface of the macro patch, shows a good agreement with an explicit model even varying friction coefficient. The strategy proposed here can be parallelized. Then, it is possible to generalize for highly heterogeneous materials and to save a lot of computation time.
The simulation is conducted using 20 CPU cores. The simulation time of the explicit and embedded results is 22 min and 11 min respectively. Despite the heterogeneous model configuration that we choose in this work is quite simple, 50 percent reduction of computation time is obtained. For highly heterogeneous material, computation time will be reduced more.
Another advantage of the relocation strategy is a better estimation of the contact pressure evolution, thus allowing precision in tribological phenomena consideration, such as wear and tear. Since wear depends on the distribution of contact pressure [15], the location of heterogeneities due to contact forces is important in computing the loss of wear volume. Our strategy allows not only to consider wear modeling, but also to introduce other mechanisms such as decohesion, cracking, etc.