3D simulation of dry friction of metal-based composites

A numerical 3D model of friction at mesoscale of the contact patch between steel samples and composite steel samples with carbon nanoinclusions has been developed based on the movable cellular automaton method. The response function of automata simulating steel corresponds to an elastic-plastic body with linear hardening, and the function for simulating carbon corresponds to an elastic-plastic body with bilinear hardening. Von Mises criterion is used for breaking bonds between the automata, and the criterion based on plastic heat is used for coupling the automata. The values of basic parameters of the model, as well as the loading conditions, that allow the formation of the characteristic quasiliquid layer in the friction zone are determined. It was shown that the presence of carbon nanoinclusions in the contact region leads to decreasing of the coefficient of friction, lower values of stress, and higher shear strain in the friction zone.


Introduction
It is known that one of the main reasons for machine breakdown is the wear of contacting joints. However, the phenomenon of friction has not only a negative effect. Thus, a number of units of modern machines are based on the use of the phenomenon of friction: mechanical brakes, friction mechanisms, gears; as well as some technological processes: friction stir welding and friction stir processing. That is why the study of the phenomena occurring in friction is one of the priority problems of modern technical science [1][2][3][4].
Today it is known that the processes of friction and wear are a combination of successive transitions of material from one body to another. In the contact zone of two bodies, intense processes of deformation, mixing and damage accumulation occur. Carrying out experimental measurements directly in this zone in dynamics is too difficult to implement, therefore, to study the processes taking place in the contact zone based on the methods of numerical simulation is in great demand [3][4][5][6]. The ability to explicitly describe the accumulation and generation of damage, as well as mass mixing, makes the method of movable cellular automata (MCA) a powerful tool for studying the dynamics of contact interaction [7][8][9].
Most of MCA applications for sliding friction are made in 2D statement [10][11][12][13], the aim of the work is to develop a 3D numerical model of the sliding friction of steel samples at the scale of contact patch and study on its basis the influence of carbon nanoinclusions on the sliding friction.
2 Description of the model

Movable cellular automaton method
Movable cellular automaton method is a numerical method in solid mechanics based on particle approach. It is assumed that the material considered consists of a certain number of particles (called movable cellular automata) of a finite size that can move, rotate and change their state due to interaction with neighbours thereby simulating real deformation and physical processes [7][8][9]. The motion of the particles is governed by Newton-Euler equations are the location, rotation velocity, mass and moment of inertia tensor of i-th automaton, pair ij F is the interaction force of the pair of i-th and j-th automata,  i F is the volume-dependent force acting on i-th automaton and depending on the interaction of its neighbours with the remaining automata. The upper dot denotes a time derivative. In the second equation, , here q ij is the distance from the centre of i-th automaton to the point of its interaction with j-th automaton, is the unit vector directed from the centre of i-th automaton to the j-th one and r ij is the distance between automata centres, ij K is the torque caused by relative rotation of automata in the pair.
It has to be noted, that automata of the pair may represent the parts of different bodies (contacted pair) or one body (bonded pair, its interaction really is not a contact one). The automaton size is characterized by one parameter d i , but it does not mean that the shape of the automaton is a sphere in 3D or a circle in 2D. The real shape of the automaton is determined by its "contacts" with neighbours.
For isotropic material the volume-dependent part of the total force can be written as follows [13][14][15]: where P j is the pressure in the volume of the neighbouring automaton j, S ij is the area of interaction surface of automata i and j, and A is defined by the ratio of elastic properties of the material. This allows us to rewrite the total force acting on automaton i as a sum of normal n ij F and tangential τ ij F components: where n pair, ij F and τ pair, are the corresponding pair interaction forces depending respectively on the automata overlap h ij and relative tangential displacement shear ij l . To compute components of the average stress tensor in the bulk of automaton i we can use homogenization procedure described in [14,15]: where  and  denote the global coordinate axes X, Y, Z, V i is the automaton volume, Because each automaton of a pair may represent different material, the pair overlap is somehow distributed between both i-th and j-th automata: where  denotes an increment per integration time step t  . The strain distribution rule in the pair is defined by the expression for computing interaction forces. For central interaction it is similar to Hooke's law for diagonal components of stress tensor: where K and G are the bulk and shear moduli of i-th automaton material, the pressure P i may be computed using (4) at the previous time step or by a predictor-corrector scheme.
A pair of automata can be in one of two states: bound and unbound. Thus, in MCA, fracture and coupling of fragments (crack healing, microwelding, etc) are simulated by the corresponding switching of the pair state. Switching criteria depend on physical mechanisms of the material behaviour. Note, that knowing stress and strain tensor in the bulk of an automaton makes possible direct application of conventional fracture criteria written in tensor form. Unbonded (i.e. contacting) automata do not resist to tension, and their tangential force is limited by the dry friction model (for example, Amonton's law, Dieterich's model and so on). To model a new bond formation in the pair of unbonded automata we use special criteria based on two threshold values. First is the magnitude of the central compression strength, it should be equal to yield strength, i.e. pair should experience plastic deformation. The second is the magnitude of plastic work (heat) in the pair, which means that to form a new pair of bonded automata the forces applied to the pair have to perform the work equal to the energy of the new chemical bond.
A detailed description of the inter-automata forces and torques, as well as plastic flow and transitions from bonded state to unbonded one and vice versa, can be found elsewhere [13,15]. Here we just note, that to describe elastoplastic flow of the material, it was proposed to use the plastic flow theory (namely von Mises model) by adapting the algorithm of Wilkins in the MCA method [15,16]. To integrate the motion equations, an explicit velocity Verlet algorithm is used. In this case, the value of the time step is limited above by the time of sound propagation throughout the bulk of automaton.

A model for sliding friction
In [10] it was shown that MCA method allowed studying numerically so-called quasi-fluid nano-layer between rubbing bodies. That and other MCA simulation of friction in contact patch were conducted in 2D statement [11][12][13]. 3D simulations of contact problems have been mainly devoted to indentation and scratching [14,16,17], as well as friction stir welding [14]. Herein, we developed 3D MCA model that allows simulating friction between steel bodies and steel-based composite containing carbon nano-inclusions at the scale of contact patch.
First, we describe the model for pure steel bodies. We consider steel as an elastic-plastic material with bilinear hardening. Plastic behaviour was described using plastic flow theory with Von Mises criterion. The main material properties are listed in Table 1, the response function of the corresponding automata, which actually describes stress-strain behaviour of the material in uniaxial tension, is depicted in Fig. 1,a. Geometrically the model consists of two interacting bodies (Fig. 1,b). Both of them are divided into two parts (shown in different colours in Fig. 1,b), one part contains the rough surface contacting with the other body, and another larger part is the base of the body. The model including its geometry as well as the rough interface is created using special preprocessor of the in-house software MCA3D_Friction. Dimensions of the bases are 1.4×1.4×1.9 µm, and dimensions of the contacting parts are 1.4×1.4×24 µm; automaton size is equal to 7.5 nm.

4
To simulate sliding friction, the following loading conditions are used. In the initial state, the interacted bodies have full contact along the rough surface. So, at the first stage the bodies are moved along the vertical axis Z in the opposite direction to separate them, and then shifted along axis Y at half of the bodies' width. For approaching the bodies to each other the velocity V z along Z-axis was applied to the automatons of the top and bottom layers of the specimen. After touching the upper and lower bodies, the forces directed along the Z-axis (compressing pressure), as well as the horizontal velocity V Y along Y-axis, were applied to the automatons of the top and bottom layers of the specimen. The specified loading conditions at the top and bottom layers of the model are symmetrical, i.e. they have the same magnitude, but the actions are applied in opposite directions. To avoid too hard loading in the friction zone, a gradual increase in the applied loads is implemented. To imitate the large size of real material along the vertical direction, it is necessary to damp the elastic waves generating in the friction zone and propagating along this direction. To this end, a special viscous damping force along Z-axis is applied for the automata of both loading layers. The magnitude of the damping coefficient was estimated in order to dissipate longitudinal elastic wave after propagating the distance of the double height of the specimen. To imitate the environment of the selected simulation box along the axes X and Y, periodic boundary conditions are set on the lateral faces of the specimen along these axes correspondingly.
To get a quasi-fluid layer in the contact zone it is necessary to allow the unbonded automata to switch to bonded state [10,15]. The criterion for such a transition is chosen based on the magnitude of plastic heat generated in the pair of unbonded automata as a result of plastic deformation while interacting with all their neighbours. The critical value of this plastic heat is equal to half of the energy required for breaking the bond in a pair of bonded automata.
The second model is aimed to simulate friction between steel-based composite. For this purpose, we generate carbon inclusions in the interacting parts of the bodies. Their shape was a cylinder with a radius of 10 and a length of 70 nm; the volume fraction of the inclusions relative to the volume of the corresponding part is about 2 %. The physicalmechanical properties of the model carbon correspond to carbon nanofibres and are listed in Table 2; the response function of such automata is depicted in Fig. 2

Simulation results and discussion
Let first consider the simulation results for the friction of pure steel bodies. Figure 3 depicts the distributions of equivalent stress and strain in the cross-section (actually, in the seven layers closed to the coordinate plane ZY) of the simulated specimen at the end of the simulation in quasi-steady regime. As one can see, high stresses are localized in the contact region. In Fig.3,a the automaton pairs, which are unbonded, are connected by lines; these automata form a quasi-fluid layer and are located in the region of high stresses. High strains form a narrower layer compared to stresses (Fig. 3,b), and contain more unbonded pairs. a) b) Fig. 3. Fields of equivalent stress (a, points are centres of automata) and strain (b, automata shown as spheres) in the cross-section of the steel specimen at the final stage of the simulation.
The presented results prove that all parameters of the model are correct and allow simulating dry friction of steel bodies at the scale of the contact patch with reproducing quasi-fluid layer formation similar to 2D simulations [10]. Figure 4 shows the fields of equivalent stress and strain in the cross-section of the rubbing composite bodies when their contacting parts contain carbon nanoinclusions. One can see that the presence of soft inclusions leads to a wider area of higher stress as well as to strain localization not only in the friction zone but also in neighbouring inclusions. Finally, let consider the dependencies of the coefficient of friction on time in both steel and composite specimens (Fig. 5). Fig. 5,a shows time dependence of the dynamic value of the friction coefficient for steel bodies for a total time of the process of 600 ns (yellow line) as well as the smoothed curve for this dependence (green line). According to this plot, one can conclude that the system reaches the quasi-steady regime at about 350 ns. Plots in Fig. 5,b depicts smoothed dependencies of the friction coefficients versus time for both steel and composite bodies and allows us to compare resistance to the relative sliding for these materials. One can see that the presence of 2 % of carbon inclusions in the contact region of steel bodies results in 5 % lower resistance to sliding friction at reaching a steady-state regime.

4 Conclusion
Herein we presented an MCA model for 3D simulation of friction between bodies of steel as well as steel based composite with carbon nanoinclusions at the scale of the contact patch. The model includes initial roughness of the contacting surfaces, and special loading conditions allowing to simulate quasi-fluid layer formation between rubbing bodies at the nanoscale. Simulation results show that the presence of 2 % of carbon inclusions in the contact region of steel bodies leads to 5 % lower resistance to sliding friction at reaching a steady-state regime.