Modeling Wear and Friction Regimes on Ceramic Materials with Positive and Negative Thermal Expansion

Computer modeling of rubbing between two surfaces with microasperities capable of expanding or contracting under conditions of frictional heating (i.e., possessing either positive and negative coefficient of thermal expansion (CTE)) allowed for the identification of wear-and-friction regimes on model ceramic materials. Assuming that no adhesion was involved in the interaction between asperities, two wear regimes—i.e., wear-free and continuous wear—have been revealed in both materials as dependent on the applied normal stress level and sliding velocity. The effect of the normal load on wear rate was similar for both positive and negative thermal expansion ceramics. Sliding velocity has a qualitatively different effect on the wear of materials with either positive or negative thermal expansion. The results indicated that the feasibility of reconstructing wear maps was common for both positive and negative CTE ceramics in terms of dimensionless mechanical and thermophysical characteristics.


Introduction
Ceramic materials possess high thermal stability, wear resistance and are thus widely used for fabricating components intended for high-temperature loading conditions [1][2][3][4][5] such as sealants, bearings, gaskets, etc., [6]. Ceramics are, as a rule, brittle or quasi-brittle materials of low-surface energy whose shear strength is greatly dependent on contact pressure. In addition, ceramics are commonly characterized by poor thermal conductivity. Such a combination of characteristics must be taken into account when considering their behavior in sliding contact.
First, their ability to relieve high local stresses by nucleation and propagation of crystalline lattice defects is strongly limited. Therefore, no adaptation mechanisms connected with generation of larger scale load-carrying interfacial elements (e.g., prows, wedges, and flakes) [7,8] that would completely separate the interacting surfaces and form the actual relief of the contact zone (a set of spots of local contact of surfaces with interfacial elements) are feasible [9,10]. Correspondingly, both friction and wear of ceramic materials are determined mainly by conditions created in direct contacting between tribologically mated surfaces of the rubbing bodies. Secondly, poor thermal conductivity results in intensive frictional heating of ceramics in sliding with flash temperatures as high as 1000-1500 K [11][12][13]. In its turn, frictional heating is accompanied by thermal expansion and additional increasing of the contact stress on the real contact areas. Thermal stresses arise also at the grain or phase boundaries below the worn surfaces.
Furthermore, it should be noted that problems related to thermal stresses occurring at the interior boundaries have been widely discussed in the corresponding literature for many years, while thermal expansion of the surface asperities in their contacting in sliding

Materials and Methods
Here, we numerically studied the tangential contact of two bodies made of materials having the same mechanical and thermophysical characteristics. These materials were represented by PTE and NTE model ceramics with thermophysical and mechanical properties close to those of ZrO 2 and scandium tungstate Sc 2 W 3 O 12 , respectively (Table 1). For ZrO 2 , there is a large amount of experimental data on mechanical properties, though this is not the case for the scandium tungstate. In fact, there are no data on its strength and thermal conductivity. The compressive strength shown in Table 1 is an estimation based on a comparison of the experimental hardness values of ZrO 2 and Sc 2 W 3 O 12 . The value of the thermal conductivity was taken from considering the characteristic values of k for various tungstate's (they vary in the range from 1 to 2) [51]. Data on ceramic material strength, as noted in the literature, refer mainly to macroscopic samples. At the same time, it is known that not only the absolute values of strength but also their ratios (in particular, the value of b) depend on the scale of the samples. Local (meso-or microscale) values of the parameter b can be significantly lower than those of macroscopic samples, since the latter may contain large internal defects. Therefore, Table 1 shows the range of the values of the parameter b from 5 (typical macroscopic strength ratio for many ceramics) to 1 (reasonable lower limit). Since internal defects have the greatest  (Table 1 shows the corresponding range of values).
The choice of ZrO 2 and Sc 2 W 3 O 12 for the numerical study and comparative analysis is substantiated by the following reasons. First, they are characterized by the stability of their structure and properties over a wide temperature range from the room temperature to 1300 • C. Second, they have rather high values of elastic moduli and strength, which allow using them in highly loaded friction units. These ceramic materials also have close values of specific heat capacity and thermal conductivity, which make them especially promising as components of Sc 2 W 3 O 12 /YSZ composites intended for various high-temperature applications [52]. Finally, the CTE values of these ceramics are close in absolute value but have the opposite sign (Sc 2 W 3 O 12 is NTE material, while ZrO 2 is a traditional PTE material).
Theoretical investigations into wear on PTE and NTE ceramic materials have been performed utilizing a numerical thermomechanical model that can be related to a group of bearing-area (asperity-based) contact models [48,49]. Two bodies were contacting with each other by means of asperities on their surfaces ( Figure 1).  1 shows the range of the values of the parameter b from 5 (typical macroscopic strength ratio for many ceramics) to 1 (reasonable lower limit). Since internal defects have the greatest effect on the strength of the material under tensile conditions, the variation of the parameter b was performed by varying the tensile strength σtens (Table 1 shows the corresponding range of values). The choice of ZrO2 and Sc2W3O12 for the numerical study and comparative analysis is substantiated by the following reasons. First, they are characterized by the stability of their structure and properties over a wide temperature range from the room temperature to 1300 °C. Second, they have rather high values of elastic moduli and strength, which allow using them in highly loaded friction units. These ceramic materials also have close values of specific heat capacity and thermal conductivity, which make them especially promising as components of Sc2W3O12/YSZ composites intended for various high-temperature applications [52]. Finally, the CTE values of these ceramics are close in absolute value but have the opposite sign (Sc2W3O12 is NTE material, while ZrO2 is a traditional PTE material).
Theoretical investigations into wear on PTE and NTE ceramic materials have been performed utilizing a numerical thermomechanical model that can be related to a group of bearing-area (asperity-based) contact models [48,49]. Two bodies were contacting with each other by means of asperities on their surfaces ( Figure 1). Each surface was characterized by the following roughness parameters: where m is the number of surface asperities in the model, i Y is the position of i-th asperity tip along the Y-axis, is the mean position of asperity tip along the Y-axis.
The initial surfaces were characterized by the same initial values of the statistical roughness. Particular initial values of asperity height ( i Y 0 ) were set out using a stochastic algorithm, which provided prescribed initial values of Ra and Rmax. The initial values of the asperity effective contact area ( i S 0 ) were set out randomly within the ranges limited Each surface was characterized by the following roughness parameters: where m is the number of surface asperities in the model, Y i is the position of i-th asperity tip along the Y-axis, Y mean = 1/m∑ m i=1 Y i is the mean position of asperity tip along the Y-axis.
The initial surfaces were characterized by the same initial values of the statistical roughness. Particular initial values of asperity height (Y i 0 ) were set out using a stochastic algorithm, which provided prescribed initial values of R a and R max . The initial values of the asperity effective contact area (S i 0 ) were set out randomly within the ranges limited by prescribed minimum and maximum values [S min 0 , S max 0 ]. The limits S min 0 and S max 0 were the same for both surfaces.
All asperities had the same width ∆x, i.e., the horizontal step of the roughness was kept constant along the X-axes of both surfaces. Each of the asperities was capable of experiencing deformation independently of the other ones and was characterized by the following geometrical parameters: asperity center position along the X-axis; position of the asperity tip along the Y-axis; and effective contact area S. The latter has the meaning of the real contact area as formed between the asperity and plane rigid punch. Note that the instant value of S can change during sliding of the punch along the asperity, and its variation depends on a number of factors, including form-factor of asperity, smaller scale roughness, applied pressure, etc. That is why S would be rather an effective integral parameter to characterize a time average value of the real contact area.
The keystones of the asperity-based thermomechanical model are as follows: (1) Adhesion between the surfaces was not considered because ceramic materials have low surface energy. An assumption was made that interaction between the asperities was limited to only opposite asperities, which contacted with zero gap between their tips. The force of interaction between two opposite asperities included two components: normal interaction force F i N along the Y-axis (resistance to compression) and tangential interaction force F i τ along the X-axis. The latter is a complex one as it consists of two components: the shear resistance force F i sh limited above dry friction force and the asperity bending resistance force F i bend . All the interaction forces are formulated in the approximation of linear elastic materials.
(2) Dry friction of opposite asperities generates heat on the real contact areas so that the asperities experience thermal expansion (including the change in the positions of their tips along the Y-axis). Local heating is accompanied by establishing subsurface thermal gradients. The hear redistribution is computed by seeking a solution to the corresponding heat conductance equation.
(3) Asperities are assumed to be destructible with the fracture criterion in the form of the two-parameter Drucker-Prager criterion [57,58], which is widely used to determine the limiting state of brittle and quasi-brittle materials including ceramics. Fracture of asperity results in removal of a material fragment from the asperity. The height and volume of the asperity fragment removed is determined by overlapping between the asperities.
(4) Wear debris are carried away from the contact zone during relative tangential displacement of the bodies by a distance ∆x (surface roughness step) and would have no effect on further course of surface interaction.
Further details on the asperity-based thermomechanical model are provided in Appendix A. Contact interaction of the surfaces is provided by loading them with external compression force F load applied to the upper body along the Y-axis ( Figure 1). The lower boundary of the lower body is fixed with respect to both axes (X and Y). A periodic boundary condition along the X-axis was applied to the bodies.
Sliding friction was modeled by discretely moving the upper body over the lower one along the X-axis at a constant sliding speed V sl (Figure 1). The discreteness of the horizontal movement of the upper body implies that each asperity of the upper body is above some asperity of the lower body during the time interval ∆t contact = ∆x/V sl . During the time ∆t contact , the asperity-asperity interaction force and the heat release power are calculated. Additionally, the heat conduction equation is solved. A change in the temperature of asperities leads to a change in their heights. Thus, during ∆t contact , the ensemble of contacts is in a state of dynamic equilibrium. At the end of the time interval ∆t contact , the upper body shifts by ∆x along the X-axis so that the whole system of asperity contacts is fully modified.
Vertical movement of the upper body is calculated from its inertia and the static force equilibrium. The static force equilibrium condition means that the magnitude of the applied force F load is equal to the lower body reaction force F repulse : where F i N is the normal interaction force in the pair of i-th opposite asperities. The vertical position of the upper body according to the static force equilibrium condition F Load = F N is limiting one, achieved in the absence of tangential movement. During the tangential movement of the upper body, its real vertical position "lags behind" the limiting (equilibrium) value due to inertial effect caused by the finite height and mass of the body. To take into account the inertial effect, the velocity of the vertical displacement of the upper body is calculated as follows: V y = (Y limit − Y cur )/t inertia , where Y cur is the current Y-coordinate of the upper surface of the upper body, Y limit is its equilibrium value corresponding to the condition F Load = F N , and t inertia is the time parameter characterizing the inertia of the upper body. The value of t inertia was assumed to be equal to twice the travel time of the longitudinal elastic wave from the lower (contact) surface to the upper surface of the upper body: t inertia = 2L upper /V P , where L upper is the height of the upper body, V P is the velocity of the longitudinal elastic wave in the material of the upper body. The value of Y limit is determined at each time step to calculate the current value of V y .
The linear dimensions of each body were 100∆x along the X-axis (the contact surface of each body contains 100 asperities) and 200∆x along the Y-axis. The statistical parameters of the initial roughness of the contact surfaces were also the same: ∆x = 25 µm, R max ≈ 1.2 µm, R a ≈ 0.3 µm. The same value of the local COF (coefficient of dry friction of a pair of interacting asperities) µ = 0.1 was used for PTE and NTE ceramics.
The effect of the loading parameters on friction steadiness, wear rate, and COF was studied. The following parameters were varied:

•
The normal pressure σ Load = F Load /S ext , where S ext = m(∆x) 2 is the area of the outer surface of the bodies.

•
The sliding velocity of the upper body is V sl .
As soon as the wear rate is determined, among others, by the strength characteristics, the tensile strength σ tens was considered as a parameter varying in the ranges as shown in Table 1.
The stability of sliding friction was analyzed using the Y mean (t) and COF(t) dependencies for different values of applied normal pressure σ Load and sliding velocity V sl . Here, Y mean is the mean coordinate of the bottom body asperity tips, and COF is the coefficient of were F i τ is the tangential interaction force in the pair of i-th opposite asperities.

Results
The study showed that thermally induced change in the linear dimensions (and particularly the height) of asperities significantly affected the process of friction and wear, and this effect was qualitatively different for either PTE or NTE ceramics. Linear dimensions of asperities change as a result of frictional heat release (see Equation (A11) in Appendix A) and following heat removal according to heat flux redistribution (Equation (A12) in Appendix A). Not only does the heat sink from the surface into the contacting bodies but it also redistributes temperature between the asperities; furthermore, an equalization of their temperatures also occurs. The heat flux balance includes components such as (i) heating real contact areas, (ii) heat removal, and (iii) temperature equalization. The balance of these thermal processes determines the mean change of the asperity height ∆Y mean therm = 1/m∑ m i=1 ∆Y i therm (∆Y i therm is the thermally induced asperity height change), as well as the homogeneity of the ∆Y i therm distribution on the worn surface. The average heat release power is controlled by the values of the applied pressure σ Load and sliding velocity V sl , while the homogeneity of the worn surface heating is determined by a specific number of local surface contacts (the number of contacts per unit length of the surface) and distribution of contact pressures on the surface. These characteristics depend on both thermal expansion of the asperities and their wear (local failure when the condition (A15) is met). As shown by [50], in which the authors examine classical PTE materials, thermal expansion contributes to enhancing wear by reducing the specific number of contacts and thus increasing the nonuniformity of local contact pressure distribution. For NTE materials, thermal contraction of asperities, on the contrary, suppresses wear due to asperity height leveling, increasing the real contact area and more even distribution of contact pressure. However, the contribution of wear to the worn surface leveling depends on the material characteristics, primarily on the strength properties of the material (in particular, on the parameter b = σ comp /σ tens ).
Therefore, loading conditions such as σ Load and V sl , as well as the material properties (including the strength parameter b), may be responsible for the established regime of friction and wear in friction pairs composed of PTE or NTE materials. The numerical study carried out in this work was aimed at revealing the contribution of each of these parameters to the frictional behavior of brittle PTE and NTE materials (including the possibility of implementing a wear-free friction mode). Figure 2 shows typical wear vs. time dependencies for the PTE model material at different applied pressures (σ Load ) and the same sliding velocity (V sl ). The wear is characterized by a change in the mean position of the surface asperity tips along the Y-axis

Friction and Wear of PTE Material
Here, as well as below, we will present the simulation results only for one of the bodies; for the lower body, it can be observed that the wear mode qualitatively changes with increasing the pressure.
wear due to asperity height leveling, increasing the real contact area and more even distribution of contact pressure. However, the contribution of wear to the worn surface leveling depends on the material characteristics, primarily on the strength properties of the material (in particular, on the parameter Therefore, loading conditions such as σLoad and Vsl, as well as the material properties (including the strength parameter b), may be responsible for the established regime of friction and wear in friction pairs composed of PTE or NTE materials. The numerical study carried out in this work was aimed at revealing the contribution of each of these parameters to the frictional behavior of brittle PTE and NTE materials (including the possibility of implementing a wear-free friction mode). Figure 2 shows typical wear vs. time dependencies for the PTE model material at different applied pressures (σLoad) and the same sliding velocity (Vsl). The wear is characterized by a change in the mean position of the surface asperity tips along the Y-axis

Friction and Wear of PTE Material
is the initial value of Ymean). Here, as well as below, we will present the simulation results only for one of the bodies; for the lower body, it can be observed that the wear mode qualitatively changes with increasing the pressure. At some "small" σLoad magnitudes, i.e., below a certain threshold value increases when contact pressure increases. However, in all cases, it stays less than the roughness Ra of the initial surface. When crit Load  is exceeded, wear becomes continuous, and its rate in the steady state of friction increases non-linearly with the applied pressure. Note that the transition from the running-in stage to the steady wear regime also occurs during the first hundredths of a second.
Wear is accompanied by smoothing the initial surface roughness (Figure 3) as shown by reduction of the Ra values in sliding. It can be observed that the steady-state value of At some "small" σ Load magnitudes, i.e., below a certain threshold value σ crit Load (σ crit Load ≈ 30 MPa for the PTE material under consideration), wear occurs only at the initial stage of sliding, for less than 0.1 s (Figure 2b). After such a running-in stage, the surfaces can interact without wear if no adhesion occurs between them. The amount of wear |∆Y mean | increases when contact pressure increases. However, in all cases, it stays less than the roughness R a of the initial surface. When σ crit Load is exceeded, wear becomes continuous, and its rate in the steady state of friction increases non-linearly with the applied pressure. Note that the transition from the running-in stage to the steady wear regime also occurs during the first hundredths of a second.
Wear is accompanied by smoothing the initial surface roughness (Figure 3) as shown by reduction of the R a values in sliding. It can be observed that the steady-state value of R a consistently decreases with increasing the pressure. At the same time, at low pressures σ Load when the "wear-free" friction mode is established, the R a reduction is not more than 30-35%. With σ Load > σ crit Load , the R a value dramatically decreases to a certain lower limit (R a ≈ 10 nm). Thus, the steady wear regime is characterized by the nanoscale roughness of the worn surfaces. In this regard, two facts should be noted. First, the rate of roughness reduction at a constant sliding speed depends on the degree of excessing the critical pressure. As the ratio σ Load /σ crit Load increases, the characteristic time needed for the R a to reach a steady level decreases from 10s of seconds to hundredths of a second (the speed of surface smoothing increases by three orders of magnitude). Secondly, the steady nanoscale value of R a is the same for all values of σ Load /σ crit Load . Its specific value does not depend on the sliding velocity (this will be shown below) and is determined by the material parameters of the rubbing bodies. sure. As the ratio σLoad/ crit Load  increases, the characteristic time needed for the Ra to reach a steady level decreases from 10s of seconds to hundredths of a second (the speed of surface smoothing increases by three orders of magnitude). Secondly, the steady nanoscale value of Ra is the same for all values of σLoad/ crit Load  . Its specific value does not depend on the sliding velocity (this will be shown below) and is determined by the material parameters of the rubbing bodies.  Wear occurs by failure of the highest asperities, which are therefore most highly loaded. In this case, the local wear volume depends on value of the asperity overlapping, i.e., on the local contact pressure. Local contact pressures, in turn, are determined not only by the applied pressure σLoad but also by the number of contact spots. Smoothing the surface during the running-in stage leads to a multiple increase in the number of contact spots and, accordingly, to reducing the local contact stress and stress leveling over the asperity ensemble. Figure 4 shows the examples of changes occurring to the asperities and local contact pressures during sliding at various values of the applied pressure σLoad.
Only the highest asperities experience wear at low values of the applied load σLoad < crit Load  . The corresponding increase in the number of contact spots ensures a decrease in the maximum local pressures to a level at which no further local asperity failure occurs. Such a situation can be characterized as wear-free sliding even if the proportion of contact spots to the total number of asperities is rather small (less than 20%).
When exceeding the critical value of the applied load (σLoad > crit Load  ), the wear of the highest asperities can no longer provide the necessary reduction in local contact pressures. In other words, there is always a certain number of contact spots at which local stresses exceed their critical value (condition (A15) is fulfilled) and local failure occurs. It is important to note that the wear-free sliding does not occur even when the surface roughness is reduced to the nanoscale value and the percentage of contact spots is more than 50%. Such an unobvious effect is a consequence of frictional heating and thermal expansion of the asperities. Indeed, in the process of relative sliding, the highest asperities of the lower surface occasionally contact those of the upper surfaces and local stresses exceed their critical value with ensuing local failure. At the same time, other asperities may interact without fracture and experience heating with the corresponding thermally induced increasing of their heights. Later, these asperities can be also destroyed when they meet the Wear occurs by failure of the highest asperities, which are therefore most highly loaded. In this case, the local wear volume depends on value of the asperity overlapping, i.e., on the local contact pressure. Local contact pressures, in turn, are determined not only by the applied pressure σ Load but also by the number of contact spots. Smoothing the surface during the running-in stage leads to a multiple increase in the number of contact spots and, accordingly, to reducing the local contact stress and stress leveling over the asperity ensemble. Figure 4 shows the examples of changes occurring to the asperities and local contact pressures during sliding at various values of the applied pressure σ Load .  In Figure 5, it is shown how such integral characteristics as mean contact pressure Pmean and number of contacts Ncont normalized to the number of surface asperities change with the sliding time. It can be observed that as the applied pressure σLoad increases, the dynamics of variation of the average contact pressure undergoes qualitative changes. At small loads σLoad < crit Load  (wear-free sliding), the average contact pressure reaches a Only the highest asperities experience wear at low values of the applied load σ Load < σ crit Load . The corresponding increase in the number of contact spots ensures a decrease in the maximum local pressures to a level at which no further local asperity failure occurs. Such a situation can be characterized as wear-free sliding even if the proportion of contact spots to the total number of asperities is rather small (less than 20%).
When exceeding the critical value of the applied load (σ Load > σ crit Load ), the wear of the highest asperities can no longer provide the necessary reduction in local contact pressures. In other words, there is always a certain number of contact spots at which local stresses exceed their critical value (condition (A15) is fulfilled) and local failure occurs. It is important to note that the wear-free sliding does not occur even when the surface roughness is reduced to the nanoscale value and the percentage of contact spots is more than 50%. Such an unobvious effect is a consequence of frictional heating and thermal expansion of the asperities. Indeed, in the process of relative sliding, the highest asperities of the lower surface occasionally contact those of the upper surfaces and local stresses exceed their critical value with ensuing local failure. At the same time, other asperities may interact without fracture and experience heating with the corresponding thermally induced increasing of their heights. Later, these asperities can be also destroyed when they meet the highest asperities of the opposite surface. The balance of processes of "thermal growth" of asperities and their wear determines the steady value of roughness and wear rate.
In Figure 5, it is shown how such integral characteristics as mean contact pressure P mean and number of contacts N cont normalized to the number of surface asperities change with the sliding time. It can be observed that as the applied pressure σ Load increases, the dynamics of variation of the average contact pressure undergoes qualitative changes. At small loads σ Load < σ crit Load (wear-free sliding), the average contact pressure reaches a steady level during the running-in stage within the first tenths of a second. In this case, the dependence P mean (σ Load ) demonstrates a non-linear growth. When critical load (σ Load > σ crit Load ) is exceeded, the value of P mean gradually decreases and tends to a certain lower limit, which is due to gradual wear and surface smoothing (Figures 2, 3 and 4b). The characteristic time required for P mean to reach a steady state is 1-2 orders of magnitude longer than that in the wear-free regime as well as it is decreased with the pressure. In this case, the steady-state value of P mean is practically independent of σ Load . It is important to note that as σ Load increases, the specific number of contact spots is increased too. When N cont approaches 100% level, the dynamics of P mean undergoes a further qualitative change. Like with the wear-free friction, the value of P mean reaches its steady value already during the first tenths of a second; however, the oscillation amplitude is an order of magnitude higher than that of with the wear-free sliding. Further increase in σ Load does not lead to any qualitative changes in the dynamics, yet the steady-state value of P mean increases in proportion to σ Load because almost all asperities interact with each other. It follows from Figure 5 that two wear regimes can be distinguished at applied pressures σLoad > crit Load  , such as mild wear and catastrophic wear. At moderately high loads, when the steady-state specific number of contacts is no more than 80% or 90%, the failure criterion is met only for a small number of asperities. Accordingly, the wear rate is low It follows from Figure 5 that two wear regimes can be distinguished at applied pressures σ Load > σ crit Load , such as mild wear and catastrophic wear. At moderately high loads, when the steady-state specific number of contacts is no more than 80% or 90%, the failure criterion is met only for a small number of asperities. Accordingly, the wear rate is low and slowly increased with σ Load . When a certain threshold level σ Load is reached, the number of contacts approaches 100%, and, at the same time, the number of simultaneously worn asperities increases multiply, which determines the multiple increase in the total wear rate. It should be noted that, in contrast to the distinguished (sharp) boundary between the wear-free and continuous wear friction regimes (σ crit Load ), the boundary between mild wear and catastrophic wear is diffuse within the interval ∆σ Load ≈ 10 MPa.
The temperature-induced change in the linear dimensions of asperities is determined by the rate of their frictional heating, which depends not only on applied pressure, but also on the sliding velocity. Indeed, increasing V sl will cause the enhanced heat release on the contact spot thus accelerate its thermal expansion. The simulation results show that in sliding between two PTE materials, this effect leads to a threshold dependence of the wear rate on V sl . Figure 6 shows an example of wear curves obtained at σ Load = 10 MPa and various V sl . It follows from Figure 5 that two wear regimes can be distinguished at applied pressures σLoad > crit Load  , such as mild wear and catastrophic wear. At moderately high loads, when the steady-state specific number of contacts is no more than 80% or 90%, the failure criterion is met only for a small number of asperities. Accordingly, the wear rate is low and slowly increased with σLoad. When a certain threshold level σLoad is reached, the number of contacts approaches 100%, and, at the same time, the number of simultaneously worn asperities increases multiply, which determines the multiple increase in the total wear rate. It should be noted that, in contrast to the distinguished (sharp) boundary between the wear-free and continuous wear friction regimes ( crit Load  ), the boundary between mild wear and catastrophic wear is diffuse within the interval ΔσLoad ≈ 10 MPa.
The temperature-induced change in the linear dimensions of asperities is determined by the rate of their frictional heating, which depends not only on applied pressure, but also on the sliding velocity. Indeed, increasing Vsl will cause the enhanced heat release on the contact spot thus accelerate its thermal expansion. The simulation results show that in sliding between two PTE materials, this effect leads to a threshold dependence of the wear rate on Vsl. Figure 6 shows an example of wear curves obtained at σLoad = 10 MPa and various Vsl. At small Vsl, below its critical value, at the initial stage of friction (during the first tenths of a second), only the highest asperities are subjected to wear and the surface roughness decreases by several tens of percent ( Figure 7). This ensures that the pressures in the local contact spots are equalized and reduced to a level where no further wear may occur At small V sl , below its critical value, at the initial stage of friction (during the first tenths of a second), only the highest asperities are subjected to wear and the surface roughness decreases by several tens of percent ( Figure 7). This ensures that the pressures in the local contact spots are equalized and reduced to a level where no further wear may occur (Figure 4a,c and Figure 8). In other words, at V sl < V crit sl , the wear-free regime may be established (in the given example V crit sl ≈ 1 m/s). When the threshold value of the sliding speed is exceeded, the wear of the contacting surfaces becomes continuous, and the wear rate increases non-linearly with the V sl . At the same time, similarly to the results described above for various σ Load , the roughness value decreases to a certain minimum nanoscale value (Figure 7). The quantitative value of this steady level does not depend on the sliding velocity (as well as on σ Load ), that is, it is a material characteristic. Similarly to what happens in the case of varying σ Load , the characteristic time for R a to reach its minimum level decreases with the sliding velocity. exceeded, the wear of the contacting surfaces becomes continuous, and the wear rate increases non-linearly with the Vsl. At the same time, similarly to the results described above for various σLoad, the roughness value decreases to a certain minimum nanoscale value ( Figure 7). The quantitative value of this steady level does not depend on the sliding velocity (as well as on σLoad), that is, it is a material characteristic. Similarly to what happens in the case of varying σLoad, the characteristic time for Ra to reach its minimum level decreases with the sliding velocity. Despite the apparent similarity of the effects that are determined either by applied pressure σLoad or sliding velocity Vsl, the mechanisms that govern the transition from the wear-free to the steady wear regime are different. Thus, the existence of threshold type dependence of the wear rate on the applied load σLoad is associated with achieving a certain threshold level of local pressures, at which contact spots regularly appear where local stresses exceed their critical value and wear of the corresponding asperities occurs. Therefore, the characteristic crit Load  is determined by the strength of the material. Thermal expansion has a significant effect on the crit Load  , i.e it reduces it as will be shown in the next section by the example of an NTE material, nevertheless it is not the main reason for the transition from wear-free to steady wear. Thermally induced change in the linear dimensions of asperities is the reason that determines the sensitivity of the sliding regime to the sliding velocity. In particular, as can be seen from Figure 8, both the number of contact exceeded, the wear of the contacting surfaces becomes continuous, and the wear rate increases non-linearly with the Vsl. At the same time, similarly to the results described above for various σLoad, the roughness value decreases to a certain minimum nanoscale value ( Figure 7). The quantitative value of this steady level does not depend on the sliding velocity (as well as on σLoad), that is, it is a material characteristic. Similarly to what happens in the case of varying σLoad, the characteristic time for Ra to reach its minimum level decreases with the sliding velocity. Despite the apparent similarity of the effects that are determined either by applied pressure σLoad or sliding velocity Vsl, the mechanisms that govern the transition from the wear-free to the steady wear regime are different. Thus, the existence of threshold type dependence of the wear rate on the applied load σLoad is associated with achieving a certain threshold level of local pressures, at which contact spots regularly appear where local stresses exceed their critical value and wear of the corresponding asperities occurs. Therefore, the characteristic crit Load  is determined by the strength of the material. Thermal expansion has a significant effect on the crit Load  , i.e it reduces it as will be shown in the next section by the example of an NTE material, nevertheless it is not the main reason for the transition from wear-free to steady wear. Thermally induced change in the linear dimensions of asperities is the reason that determines the sensitivity of the sliding regime to the sliding velocity. In particular, as can be seen from Figure 8, both the number of contact Despite the apparent similarity of the effects that are determined either by applied pressure σ Load or sliding velocity V sl , the mechanisms that govern the transition from the wear-free to the steady wear regime are different. Thus, the existence of threshold type dependence of the wear rate on the applied load σ Load is associated with achieving a certain threshold level of local pressures, at which contact spots regularly appear where local stresses exceed their critical value and wear of the corresponding asperities occurs. Therefore, the characteristic σ crit Load is determined by the strength of the material. Thermal expansion has a significant effect on the σ crit Load , i.e., it reduces it as will be shown in the next section by the example of an NTE material, nevertheless it is not the main reason for the transition from wear-free to steady wear. Thermally induced change in the linear dimensions of asperities is the reason that determines the sensitivity of the sliding regime to the sliding velocity. In particular, as can be seen from Figure 8, both the number of contact spots and the mean contact pressure change only slightly when passing over V crit sl . At V sl > V crit sl , the asperities grow faster compared to their critical growth value and the contact spots regularly appear, wherein local stresses reach their critical value and asperities fail. As the sliding velocity increases, the number of such contact spots increases too, and the wear rate increases despite that the total number of contact spots grows and average contact pressure reduces. In other words, if the V sl > V crit sl and the sliding velocity increases, both the inhomogeneity of the contact pressure distribution and the number of "critically" loaded contact spots increase too ( Figure 9). In this range of V sl values, the wear of the most highly loaded asperities and corresponding smoothing of surfaces do not compensate for the effect of their thermal expansion.
As the sliding velocity increases, the number of such contact spots increases too, and the wear rate increases despite that the total number of contact spots grows and average contact pressure reduces. In other words, if the Vsl > crit sl V and the sliding velocity increases, both the inhomogeneity of the contact pressure distribution and the number of "critically" loaded contact spots increase too (Figure 9). In this range of Vsl values, the wear of the most highly loaded asperities and corresponding smoothing of surfaces do not compensate for the effect of their thermal expansion. When critical value of either applied pressure or sliding velocity is exceeded, a qualitative change in the wear regime is also manifested in the sliding characteristics ( Figure  10) when dynamics of the coefficient of friction (COF) is dependent on the σLoad. At low pressures, when wear-free wear is realized, the COF value reaches a steady value already during the first seconds or even tenths of a second, after which it remains constant (the amplitude of COF dynamic oscillations is insignificant). The steady value of COF slightly increases with σLoad (ΔCOF ≈ 0.02 in the interval 0.3 crit ). When the threshold value of the applied pressure (σLoad > crit Load  ) is exceeded, there is an abrupt increase in the average COF value (by 0.04-0.1 depending on σLoad). In this case, the COF value is unstable. Irregular oscillations in the friction force occur that are determined by wear and tear of asperities as well as due to continuous changes of the worn surface topology. As in the wear-free mode, the average COF increases non-linearly with σLoad (the COF growth rate decreases as the applied load increases). When critical value of either applied pressure or sliding velocity is exceeded, a qualitative change in the wear regime is also manifested in the sliding characteristics ( Figure 10) when dynamics of the coefficient of friction (COF) is dependent on the σ Load . At low pressures, when wear-free wear is realized, the COF value reaches a steady value already during the first seconds or even tenths of a second, after which it remains constant (the amplitude of COF dynamic oscillations is insignificant). The steady value of COF slightly increases with σ Load (∆COF ≈ 0.02 in the interval 0.3σ crit Load < σ Load < σ crit Load ). When the threshold value of the applied pressure (σ Load > σ crit Load ) is exceeded, there is an abrupt increase in the average COF value (by 0.04-0.1 depending on σ Load ). In this case, the COF value is unstable. Irregular oscillations in the friction force occur that are determined by wear and tear of asperities as well as due to continuous changes of the worn surface topology. As in the wear-free mode, the average COF increases non-linearly with σ Load (the COF growth rate decreases as the applied load increases).
As the sliding velocity increases, the number of such contact spots increases too, and the wear rate increases despite that the total number of contact spots grows and average contact pressure reduces. In other words, if the Vsl > crit sl V and the sliding velocity increases, both the inhomogeneity of the contact pressure distribution and the number of "critically" loaded contact spots increase too (Figure 9). In this range of Vsl values, the wear of the most highly loaded asperities and corresponding smoothing of surfaces do not compensate for the effect of their thermal expansion. When critical value of either applied pressure or sliding velocity is exceeded, a qualitative change in the wear regime is also manifested in the sliding characteristics ( Figure  10) when dynamics of the coefficient of friction (COF) is dependent on the σLoad. At low pressures, when wear-free wear is realized, the COF value reaches a steady value already during the first seconds or even tenths of a second, after which it remains constant (the amplitude of COF dynamic oscillations is insignificant). The steady value of COF slightly increases with σLoad (ΔCOF ≈ 0.02 in the interval 0.3 crit ). When the threshold value of the applied pressure (σLoad > crit Load  ) is exceeded, there is an abrupt increase in the average COF value (by 0.04-0.1 depending on σLoad). In this case, the COF value is unstable. Irregular oscillations in the friction force occur that are determined by wear and tear of asperities as well as due to continuous changes of the worn surface topology. As in the wear-free mode, the average COF increases non-linearly with σLoad (the COF growth rate decreases as the applied load increases).  Figure 11 shows the COF dynamics obtained at various sliding velocities V sl . Similarly to the dependence on σ Load considered above, COF is stable in the range of V sl values corresponding to the wear-free regime, and its steady-state value is almost independent of the sliding velocity (the increase is less than 0.02). When the critical sliding velocity is exceeded, sliding becomes unstable, and the amplitude of dynamic oscillations increases multiply with the sliding velocity. In particular, at V sl = 2 m/s, the scatter of dynamic oscillations reaches 0.1, which is five times as much than for the cases shown in Figure 10.
Alike the wear-free mode, the average COF value increases with the V sl , although less significantly than in the above-discussed case for variation of the applied pressure σ Load . Figure 11 shows the COF dynamics obtained at various sliding velocities Vsl. Similarly to the dependence on σLoad considered above, COF is stable in the range of Vsl values corresponding to the wear-free regime, and its steady-state value is almost independent of the sliding velocity (the increase is less than 0.02). When the critical sliding velocity is exceeded, sliding becomes unstable, and the amplitude of dynamic oscillations increases multiply with the sliding velocity. In particular, at Vsl = 2 m/s, the scatter of dynamic oscillations reaches 0.1, which is five times as much than for the cases shown in Figure 10. Alike the wear-free mode, the average COF value increases with the Vsl, although less significantly than in the above-discussed case for variation of the applied pressure σLoad. Numerically derived effects of normal load and sliding velocity are in qualitative agreement with the experimental data obtained on traditional ceramic materials (oxides, carbides, and nitrides). It is a traditional approach to delineate two wear regimes, such as mild wear (low wear rate) and severe wear. The mild wear regime occurs usually at low contact pressure and accompanied by smoothing the worn surfaces [5,59]. The COF values are also at a minimum but may grow with pressure [60,61]. Severe wear regime is established at higher contact pressures close in the order to the tensile strength of the materials [5,59]. Not only the wear rate but also the size of wear debris increases with the load. Wear mechanism is changed when reaching the contact pressure close to the tensile strength so that whole grains may be pulled out of the surface, and both wear rate and surface roughness increase drastically [62]. Experimental results show that severe wear regime is typical during initial (running-in) stage sliding at both low and high contact pressures. The sliding velocity has its influence on both wear rate and COF when they increase with it [61]. The transition from mild-to-severe wear with increasing pressure or sliding velocity occurs rather sharply in the narrow ranges of these parameters [59]. Characteristic limiting values of sliding velocity differ from ceramics to ceramics, as well as for different geometry of experiments but stay within the 10 −2 -10 0 m/s range.
The results of modeling are consistent with those obtained from experiments that are not just qualitative. In particular, COF values in the range 0.3-0.4 were obtained on ZrO2 using a pin-on-disk scheme and both sliding speed and pressures close to those used in this work [63].
As the above-discussed results suggest, the asperity-based model used in this work not only was adequate for describing the wear and friction on ceramics but also allowed for the obtainment of reasonable quantitative characteristics except for the very high contact pressure values when some dramatic changes of wear mechanism may occur that Numerically derived effects of normal load and sliding velocity are in qualitative agreement with the experimental data obtained on traditional ceramic materials (oxides, carbides, and nitrides). It is a traditional approach to delineate two wear regimes, such as mild wear (low wear rate) and severe wear. The mild wear regime occurs usually at low contact pressure and accompanied by smoothing the worn surfaces [5,59]. The COF values are also at a minimum but may grow with pressure [60,61]. Severe wear regime is established at higher contact pressures close in the order to the tensile strength of the materials [5,59]. Not only the wear rate but also the size of wear debris increases with the load. Wear mechanism is changed when reaching the contact pressure close to the tensile strength so that whole grains may be pulled out of the surface, and both wear rate and surface roughness increase drastically [62]. Experimental results show that severe wear regime is typical during initial (running-in) stage sliding at both low and high contact pressures. The sliding velocity has its influence on both wear rate and COF when they increase with it [61]. The transition from mild-to-severe wear with increasing pressure or sliding velocity occurs rather sharply in the narrow ranges of these parameters [59]. Characteristic limiting values of sliding velocity differ from ceramics to ceramics, as well as for different geometry of experiments but stay within the 10 −2 -10 0 m/s range.
The results of modeling are consistent with those obtained from experiments that are not just qualitative. In particular, COF values in the range 0.3-0.4 were obtained on ZrO 2 using a pin-on-disk scheme and both sliding speed and pressures close to those used in this work [63].
As the above-discussed results suggest, the asperity-based model used in this work not only was adequate for describing the wear and friction on ceramics but also allowed for the obtainment of reasonable quantitative characteristics except for the very high contact pressure values when some dramatic changes of wear mechanism may occur that have never been considered in the numerical model. Moreover, such a model validation improves the reliability of results obtained in the following section on the NTE materials. It is of high importance especially if there is no experimental data on wear and friction of the NTE materials.

Friction and Wear of NTE Material
It was shown above that thermal expansion of the surface asperity had a significant effect on both wear rate and COF value of the PTE materials. In particular, thermal expansion leads to a narrowing of the ranges of pressures and sliding velocities at which a wear-free friction regime can be observed. The results obtained allow suggesting that the use of NTE materials can contribute to the expansion of these intervals. To verify this assumption, an NTE model material was considered, whose physical and mechanical properties corresponded to those of scandium tungstate (ScW 3 O 12 ). It is important to note that the elastic and strength characteristics, heat capacity and thermal conductivity of this NTE material are comparable with the corresponding characteristics of the PTE model material considered above (ZrO 2 ), and the CTE values of these two materials are opposite in sign and close in modulus. This allows analyzing the effect of CTE in various sliding conditions using real, rather than abstract model materials [50]. Moreover, the same characteristic intervals of applied pressure σ Load and sliding velocity V sl were considered for this model NTE material. Figure 12 shows wear curves for model NTE material at different σ Load values and V sl kept constant. It can be seen that alike the PTE material (Figure 2), at the end of the short initial running-in stage (wear of the highest initial asperities), a steady sliding regime was established. Two wear modes can also be distinguished: wear-free and continuous wear regimes. However, for NTE materials, the value of the critical applied stress σ crit Load is almost by a factor of three greater than that obtained for the above considered PTE pair (80 MPa and 30 MPa, respectively), despite the fact that the strength of the NTE material is only about 35% lower than the strength of the PTE material. Obviously, this effect is associated with the thermal contraction of the asperities.

Friction and Wear of NTE Material
It was shown above that thermal expansion of the surface asperity had a significant effect on both wear rate and COF value of the PTE materials. In particular, thermal expansion leads to a narrowing of the ranges of pressures and sliding velocities at which a wearfree friction regime can be observed. The results obtained allow suggesting that the use of NTE materials can contribute to the expansion of these intervals. To verify this assumption, an NTE model material was considered, whose physical and mechanical properties corresponded to those of scandium tungstate (ScW3O12). It is important to note that the elastic and strength characteristics, heat capacity and thermal conductivity of this NTE material are comparable with the corresponding characteristics of the PTE model material considered above (ZrO2), and the CTE values of these two materials are opposite in sign and close in modulus. This allows analyzing the effect of CTE in various sliding conditions using real, rather than abstract model materials [50]. Moreover, the same characteristic intervals of applied pressure σLoad and sliding velocity Vsl were considered for this model NTE material. Figure 12 shows wear curves for model NTE material at different σLoad values and Vsl kept constant. It can be seen that alike the PTE material (Figure 2), at the end of the short initial running-in stage (wear of the highest initial asperities), a steady sliding regime was established. Two wear modes can also be distinguished: wear-free and continuous wear regimes. However, for NTE materials, the value of the critical applied stress crit Load  is almost by a factor of three greater than that obtained for the above considered PTE pair (80 MPa and 30 MPa, respectively), despite the fact that the strength of the NTE material is only about 35% lower than the strength of the PTE material. Obviously, this effect is associated with the thermal contraction of the asperities. One can see from Figure 13 that for σ Load < σ crit Load (σ crit Load ≈ 80 MPa), the steady value of R a decreases by an order of magnitude with σ Load even becoming the nanoscale and reaches its minimum value when σ Load →σ crit Load . The steady value of R a increases with transition to the continuous wear regime (σ Load > σ crit Load ), albeit remaining within the nanoscale. At the same time, in the case of the PTE material, there are two characteristic minima of the steady R a values (Figure 3). The first minimum refers to the wear-free sliding (σ Load →σ crit Load ≈ 30 MPa). The minimum steady R a value (σ Load →σ crit Load ) is approximately 40% lower than the initial R a value. With transition to the continuous wear regime (σ Load >σ crit Load ), the steady value of R a (10 1 nm) is the same for all σ Load values and is an order of magnitude lower than the initial value of R a . Thus, thermally induced contraction of surface asperities is a factor contributing to extreme smoothing of surfaces even in the wear-free sliding. Figure 14 shows the dynamics of mean contact pressure and number of contact spots at different values of σLoad. As σLoad increases, the steady-state values of both Pmean and Ncont increase monotonically (Ncont reaches 100% at σLoad ≈ 50 MPa). Comparison with similar dependences for PTE material ( Figure 5) shows several qualitative differences. First, for the NTE material, the dynamics of the mean contact pressure and the number of contact spots is the same for both wearfree and continuous wear regimes. After a short initial running-in stage, both Pmean and Ncont reach their steady-state values and subsequently undergo only minor oscillations around a constant average level. At the same time, on the PTE material, crossing the boundary between the two modes leads to a strong instability of the Pmean (significant irregular fluctuations and a gradual overall decrease to a minimum level that is the same for different σLoad > crit ). Such a twofold difference cannot be explained only by the lower stiffness of the NTE material (its Young's modulus is only 25% less than the Young's modulus of the PTE material). The above-described differences are determined by the effective smoothing of the surfaces due to thermal contraction of the surface of the NTE material. Thus, thermally induced contraction of surface asperities is a factor contributing to extreme smoothing of surfaces even in the wear-free sliding. Figure 14 shows the dynamics of mean contact pressure and number of contact spots at different values of σ Load . As σ Load increases, the steady-state values of both P mean and N cont increase monotonically (N cont reaches 100% at σ Load ≈ 50 MPa). Comparison with similar dependences for PTE material ( Figure 5) shows several qualitative differences. First, for the NTE material, the dynamics of the mean contact pressure and the number of contact spots is the same for both wear-free and continuous wear regimes. After a short initial running-in stage, both P mean and N cont reach their steady-state values and subsequently undergo only minor oscillations around a constant average level. At the same time, on the PTE material, crossing the boundary between the two modes leads to a strong instability of the P mean (significant irregular fluctuations and a gradual overall decrease to a minimum level that is the same for different σ Load > σ crit Load ). In addition, full surface contact (N cont ≈ 100%) is already achieved in the wear-free mode on the NTE material, long before σ crit Load is reached (at σ Load ≈ 50 MPa ≈ 0.6σ crit Load ), while in the case of PTE material, full surface contact occurs already in the continuous wear mode (at σ Load ≈ 100 MPa ≈ 3.3σ crit Load ). Such a twofold difference cannot be explained only by the lower stiffness of the NTE material (its Young's modulus is only 25% less than the Young's modulus of the PTE material). The above-described differences are determined by the effective smoothing of the surfaces due to thermal contraction of the surface of the NTE material. the initial Ra value. With transition to the continuous wear regime (σLoad> crit Load  ), the steady value of Ra (10 1 nm) is the same for all σLoad values and is an order of magnitude lower than the initial value of Ra. Thus, thermally induced contraction of surface asperities is a factor contributing to extreme smoothing of surfaces even in the wear-free sliding. Figure 14 shows the dynamics of mean contact pressure and number of contact spots at different values of σLoad. As σLoad increases, the steady-state values of both Pmean and Ncont increase monotonically (Ncont reaches 100% at σLoad ≈ 50 MPa). Comparison with similar dependences for PTE material ( Figure 5) shows several qualitative differences. First, for the NTE material, the dynamics of the mean contact pressure and the number of contact spots is the same for both wearfree and continuous wear regimes. After a short initial running-in stage, both Pmean and Ncont reach their steady-state values and subsequently undergo only minor oscillations around a constant average level. At the same time, on the PTE material, crossing the boundary between the two modes leads to a strong instability of the Pmean (significant irregular fluctuations and a gradual overall decrease to a minimum level that is the same for different σLoad > crit ). Such a twofold difference cannot be explained only by the lower stiffness of the NTE material (its Young's modulus is only 25% less than the Young's modulus of the PTE material). The above-described differences are determined by the effective smoothing of the surfaces due to thermal contraction of the surface of the NTE material.  Figure 15 shows examples of the NTE material surface profiles obtained in wear-free ( Figure 15a) and continuous (Figure 15b) wear regimes. Not only the highest asperities are subjected to wear at this running-in stage and at "small" applied pressures (σ Load < σ crit Load ) but all of them are reduced in height (Figure 15a,c). Figure 15 shows examples of the NTE material surface profiles obtained in wear-free ( Figure 15a) and continuous (Figure 15b) wear regimes. Not only the highest asperities are subjected to wear at this running-in stage and at "small" applied pressures (σLoad < crit Load  ) but all of them are reduced in height (Figure 15a,c).
More efficient smoothing is achieved when extra contact spots are generated, which reduce the contact pressure compared to those obtained on the PTE material. This explains, in particular, high stability of Pmean and Ncont during sliding, smoother profile of the worn surfaces, and multiple increase in the value of the critical applied pressure during the transition from positive (PTE) to negative (NTE) sign of thermal expansion. These clear advantages of NTE material may, however, become a negative factor at high pressures and in continuous wear regime (σLoad > crit Load  ). Indeed, at σLoad crit Load  , a much greater real contact area is involved in sliding if local contact pressures were distributed almost uniformly. Accordingly, the wear rate in the NTE friction pair becomes high immediately after overcoming the crit Load  value. This qualitatively distinguishes NTE friction pairs from those made of PTE materials. In the latter case, the wear rate at a slight excess of crit Load  is low and gradually increases with increasing σLoad. Moreover, for the same large values of σLoad (corresponding to continuous wear in both NTE and PTE friction pairs, when in both cases the specific number of contacts is about 100%), the wear rate in the NTE friction pair is many times higher than that of the PTE friction pair. This is also a consequence of the different sign of the CTE. Thermal expansion of surface asperities in PTE material helps reducing the number of contacts (and hence the number of contacts under wear at the same time) and thus prevents an increase in the wear rate. In contrast, More efficient smoothing is achieved when extra contact spots are generated, which reduce the contact pressure compared to those obtained on the PTE material. This explains, in particular, high stability of P mean and N cont during sliding, smoother profile of the worn surfaces, and multiple increase in the value of the critical applied pressure during the transition from positive (PTE) to negative (NTE) sign of thermal expansion.
These clear advantages of NTE material may, however, become a negative factor at high pressures and in continuous wear regime (σ Load > σ crit Load ). Indeed, at σ Load →σ crit Load , a much greater real contact area is involved in sliding if local contact pressures were distributed almost uniformly. Accordingly, the wear rate in the NTE friction pair becomes high immediately after overcoming the σ crit Load value. This qualitatively distinguishes NTE friction pairs from those made of PTE materials. In the latter case, the wear rate at a slight excess of σ crit Load is low and gradually increases with increasing σ Load . Moreover, for the same large values of σ Load (corresponding to continuous wear in both NTE and PTE friction pairs, when in both cases the specific number of contacts is about 100%), the wear rate in the NTE friction pair is many times higher than that of the PTE friction pair. This is also a consequence of the different sign of the CTE. Thermal expansion of surface asperities in PTE material helps reducing the number of contacts (and hence the number of contacts under wear at the same time) and thus prevents an increase in the wear rate. In contrast, the thermal contraction of the surface asperities of the NTE of the material tends to increase the number of contacts and, accordingly, stimulate increasing the wear rate (due to the overall high level of stress, Figure 16). Thus, NTE materials are "genetically predisposed" to high wear when a critical applied load is exceeded. When the critical load is exceeded, wear mode switches to catastrophic one (while in PTE materials this leads to mild wear, which gradually turns into catastrophic wear with increasing σ Load ). the thermal contraction of the surface asperities of the NTE of the material tends to increase the number of contacts and, accordingly, stimulate increasing the wear rate (due to the overall high level of stress, Figure 16). Thus, NTE materials are "genetically predisposed" to high wear when a critical applied load is exceeded. When the critical load is exceeded, wear mode switches to catastrophic one (while in PTE materials this leads to mild wear, which gradually turns into catastrophic wear with increasing σLoad). The NTE materials revealed one more specificity in comparison with the classical PTE materials, i.e., the absence of a critical sliding velocity. This may be a consequence of the fact that the velocity dependence of wear is determined by enhanced heat release on the contact spots with increasing Vsl and corresponding height contraction. As an example, Figure 17 shows surface roughness vs. time curves for different sliding velocities and the same σLoad < crit Load  value. An increase in Vsl leads to the accelerated crossover of the friction pair to the stationary wear-free regime and, at the same time, does not change the roughness parameters as well as the amount of wear at the running-in stage. The better steadiness of sliding on the NTE materials both in wear-free and catastrophic wear modes is due to thermally induced smoothing of the contact surfaces, which is also manifested in the high stability of COF ( Figure 18). The value of COF increases nonlinearly with increasing σLoad, while the amplitude of oscillations near the mean value decreases. Crossing the crit Load  threshold is not accompanied by the visible changes in COF dynamics, despite the catastrophic wear. The NTE materials revealed one more specificity in comparison with the classical PTE materials, i.e., the absence of a critical sliding velocity. This may be a consequence of the fact that the velocity dependence of wear is determined by enhanced heat release on the contact spots with increasing V sl and corresponding height contraction. As an example, Figure 17 shows surface roughness vs. time curves for different sliding velocities and the same σ Load < σ crit Load value. An increase in V sl leads to the accelerated crossover of the friction pair to the stationary wear-free regime and, at the same time, does not change the roughness parameters as well as the amount of wear at the running-in stage.
the thermal contraction of the surface asperities of the NTE of the material tends to increase the number of contacts and, accordingly, stimulate increasing the wear rate (due to the overall high level of stress, Figure 16). Thus, NTE materials are "genetically predisposed" to high wear when a critical applied load is exceeded. When the critical load is exceeded, wear mode switches to catastrophic one (while in PTE materials this leads to mild wear, which gradually turns into catastrophic wear with increasing σLoad). The NTE materials revealed one more specificity in comparison with the classical PTE materials, i.e., the absence of a critical sliding velocity. This may be a consequence of the fact that the velocity dependence of wear is determined by enhanced heat release on the contact spots with increasing Vsl and corresponding height contraction. As an example, Figure 17 shows surface roughness vs. time curves for different sliding velocities and the same σLoad < crit Load  value. An increase in Vsl leads to the accelerated crossover of the friction pair to the stationary wear-free regime and, at the same time, does not change the roughness parameters as well as the amount of wear at the running-in stage. The better steadiness of sliding on the NTE materials both in wear-free and catastrophic wear modes is due to thermally induced smoothing of the contact surfaces, which is also manifested in the high stability of COF ( Figure 18). The value of COF increases nonlinearly with increasing σLoad, while the amplitude of oscillations near the mean value decreases. Crossing the crit Load  threshold is not accompanied by the visible changes in COF dynamics, despite the catastrophic wear. The better steadiness of sliding on the NTE materials both in wear-free and catastrophic wear modes is due to thermally induced smoothing of the contact surfaces, which is also manifested in the high stability of COF ( Figure 18). The value of COF increases non-linearly with increasing σ Load , while the amplitude of oscillations near the mean value decreases. Crossing the σ crit Load threshold is not accompanied by the visible changes in COF dynamics, despite the catastrophic wear.

Wear-Free Regime in Sliding on PTE and NTE Materials
The results of modeling showed that for any finite value of applied pressure σLoad there is a finite value of critical sliding velocity crit sl V , which divides between two wear modes (wear-free and continuous wear). The reverse statement is also true when for any finite Vsl there is a finite critical pressure

Wear-Free Regime in Sliding on PTE and NTE Materials
The results of modeling showed that for any finite value of applied pressure σ Load there is a finite value of critical sliding velocity V crit sl , which divides between two wear modes (wear-free and continuous wear). The reverse statement is also true when for any finite V sl there is a finite critical pressure σ crit Load . In other words, critical parameters σ crit Load and V crit sl are not independent from each other since both contact pressure (it is directly related to σ Load ) and sliding velocity V sl determine frictional heating rate. Intuitively, it can be suggested that σ Load and V sl should have equal effects on the wear-free-to-continuous-wear transition and that V crit sl (σ crit Load ) should be a linear dependence. In fact, the numerically obtained V crit sl (σ crit Load ) dependence for the PTE friction pair is non-linear ( Figure 19) and can be acceptably approximated using an exponential function V sl = A + e −Bσ Load , where A ≈ 0.15 m/s, B ≈ 0.0675 MPa −1 for the PTE material discussed. The existence of such a non-linearity is a direct indication that the applied pressure and sliding velocity have different contributions to wear-free-to-continuous-wear transition. Whereas sliding velocity determines mainly the power of heating, the applied pressure σ Load determines both the total number of real contact spots as well as the fraction of these spots where local shear stress achieves the local value of the dry friction force F dr . The local values of F dr are different for different contact spots, yet their average value is proportional to σ Load . Such a complex and non-linear effect of the applied normal load on the heating of particular asperities and increasing their heights results in the non-linearity of the relationship V crit sl (σ crit Load ).

Wear-Free Regime in Sliding on PTE and NTE Materials
The results of modeling showed that for any finite value of applied pressure σLoad there is a finite value of critical sliding velocity crit sl V , which divides between two wear modes (wear-free and continuous wear). The reverse statement is also true when for any finite Vsl there is a finite critical pressure   For low σ Load , the contact spot percentage with respect to total number of asperities is considerably less than 50%. This fact offers great potential for redistribution of contact stresses and hence provides relatively high values of V crit sl . Increasing the σ Load results in corresponding increasing the contact spot number that limits the efficiency of stress redistribution caused by thermal expansion of the asperities. That is why V crit sl reduction rate is decreased as a function of σ Load .
At σ Load →40 MPa, there occurs some percolation transition when the contact spot percentage becomes equal to 50%. At larger σ Load , the major part of surface asperities are in contact with the opposite ones. Accordingly, this increases the percentage of highly loaded and frictionally heated spots. Therefore, high σ Load values may provide more uniform distribution of contact stresses and more homogeneous heating of the worn surfaces. Moreover, high contact pressures and corresponding high magnitudes of F dr determine intensive heating of asperities. That is why the value of V crit sl is minimal and changes negligibly at high values of σ Load .
At σ Load ≈ 100 MPa, all asperities become contact spots. A slight additional increase in σ Load leads to local stresses in some fraction of the contact spots are always higher the material strength (the criterion of local fracture is satisfied), and the wear becomes continuous for all V sl . The limiting value of applied normal stress σ max Load is 105-110 MPa for the considered PTE material. The parameter σ max Load has the meaning of the maximum applied normal stress to achieve wear-free regime of sliding friction. Figure 19 is a wear map where different wear regimes may be distinguished. Experimental wear maps were also obtained on particular ceramic materials, as well as generalized dimensionless ones [59] when boundaries between mild and severe wear regimes were shown for several traditional PTE ceramics (including ZrO 2 ) using a coordinate axis such as the mean Hertzian pressure and sliding velocity, and where a boundary curve looks analogous to that obtained in this work for the model PTE ceramics (Figure 19). A qualitative agreement was achieved between mean Hertzian pressure and σ crit Load for ZrO 2 at the same V sl .
Another example of wear maps reconstructed in terms of dimensionless parameters such as mechanical severity of contact and thermal severity of contact was reported [6,59]. These parameters allowed for the consideration of both material characteristics and loading parameters. The above-identified examples may testify on the applicability of the wear maps for different ceramics with indication to their insufficient versatility and potential using other dimensionless variables relating to mechanical and thermophysical characteristics like those obtained using the subsurface temperature gradient established under steady friction conditions [50].
Since the above-discussed PTE material has low thermal conductivity, it is interesting how intensive the worn surface heating under limiting loading conditions would be, at least as determined according to boundary curves in Figure 19. Figure 20 shows the steady averaged surface temperatures, T steady , for boundary points (V crit sl , σ crit Load ) in Figure 19. Generally, it should be noted that T steady gradually increases with the applied pressure. This step-by-step growth is determined by changing the contacting conditions and by increasing the percentage of the contact spots. The existence of the first plateau at σ Load ≈ 20/40 MPa corresponds to a "percolation transition" of contact spot percentage N cont over the 50% threshold (30-60% for the applied pressure interval). The second one at σ Load > 70 MPa appears when N cont approaches to 100%. In this connection, it is important to note that temperature interval of stable functioning of ZrO 2 is limited from above by 1400/1700 K. One can see in Figure 20 that steady worn surface temperatures for the boundary points (V crit sl , σ crit Load ) do not exceed this limit so that this PTE material may be capable of demonstrating the wear-free sliding, as determined by the boundary in Figure 19. The high-temperature region at the boundary line between two wear regimes begs a question about the effect of temperature-dependent material characteristics on both profile and position of this boundary in the parametric plane (σLoad, Vsl). Some considerations may be derived even now from analyzing the dimensionless wear map parameters.
Let us take in to account the mechanical parameter [50]: . In this situation, temperature dependencies of heat capacity and heat conductivity of ZrO2 will grow and fall, respectively [68,69]. The CTE of ZrO2 is stable in the wide temperature range and because of these opposite temperature trends, their common influence on both Mcrit and crit Load  will be negligible. Thus, the main factor will be temperature dependence of the ratio Two regions may be delineated in the parametric plane (σLoad, Vsl), such as closed and open regions corresponding to wear-free and continuous wear, respectively. (Figure 19). Crossing the σLoad< max Load  boundary means transitioning to the continuous low wear (mild) regime. Increasing the distance from this boundary results in accelerated (non-linear) increasing the wear rate up to achieving the catastrophic wear (Figures 2 and 6). At the same time, crossing the boundary at σLoad = max Load  results in a sharp increase of the wear rate and crossover from a wear-free to a catastrophic regime.
For comparison, Figure 19 identifies an analogous boundary for the NTE material under study. As noted above, thermal height reduction of asperities in heating may serve for inducing the wear-free regime with the sliding velocity. Therefore, the wear regime boundary for NTE materials is in the form of a vertical line in the parametrical plane (σLoad, Vsl). In fact,  The high-temperature region at the boundary line between two wear regimes begs a question about the effect of temperature-dependent material characteristics on both profile and position of this boundary in the parametric plane (σ Load , V sl ). Some considerations may be derived even now from analyzing the dimensionless wear map parameters.
Let us take in to account the mechanical parameter [50]: M = σ Load E/σ 2 comp . The wearfree regime occurs at M < M crit = σ crit Load E/σ 2 comp . The magnitude of M crit is independent of E and σ comp , yet it depends on the dimensionless thermophysical parameter. For a great number of ceramics, including ZrO 2 , the temperature rises to~1000 K provides reduction in both Young's modulus and strength [64][65][66][67]. By definition of M crit , drawing into the temperature dependence of mechanical characteristics results in reducing σ crit Load . Moreover, since different contact spots have different temperatures the reduction in σ crit Load should be accompanied by blurring the boundary line between the wear regimes. In addition, both M crit and σ crit Load depend on thermophysical characteristics of the material and, therefore, their temperature dependencies will influence σ crit Load . In this situation, temperature dependencies of heat capacity and heat conductivity of ZrO 2 will grow and fall, respectively [68,69]. The CTE of ZrO 2 is stable in the wide temperature range and because of these opposite temperature trends, their common influence on both M crit and σ crit Load will be negligible. Thus, the main factor will be temperature dependence of the ratio E/σ 2 comp . Two regions may be delineated in the parametric plane (σ Load , V sl ), such as closed and open regions corresponding to wear-free and continuous wear, respectively. (Figure 19). Crossing the σ Load <σ max Load boundary means transitioning to the continuous low wear (mild) regime. Increasing the distance from this boundary results in accelerated (non-linear) increasing the wear rate up to achieving the catastrophic wear (Figures 2 and 6). At the same time, crossing the boundary at σ Load = σ max Load results in a sharp increase of the wear rate and crossover from a wear-free to a catastrophic regime.
For comparison, Figure 19 identifies an analogous boundary for the NTE material under study. As noted above, thermal height reduction of asperities in heating may serve for inducing the wear-free regime with the sliding velocity. Therefore, the wear regime boundary for NTE materials is in the form of a vertical line in the parametrical plane (σ Load , V sl ). In fact, σ crit Load = const = σ max Load for the NTE ceramic materials. Since the values of maximum critical pressure (σ max Load ) are close for the considered PTE and NTE materials, it seems appropriate to compare the corresponding values of the above considered dimensionless mechanical parameter M crit [50]. For the model, NTE material with mechanical characteristics is close to those of scandium tungstate, i.e., M crit ≈ 5.5. The physical meaning of the limiting σ Load values (σ max Load ) for both NTE and PTE materials is the same; it is relevant to estimate the limiting value of M crit for the PTE material used in this study: M crit = σ max Load E/σ c 2 ≈ 7.1. The specified values of the dimensionless mechanical parameter for PTE and NTE materials with comparable mechanical and thermophysical characteristics differ by less than 30%. Therefore, for a wide range of PTE and NTE ceramic materials, it may be assumed that M crit = 1/10.
Ceramics are pressure-sensitive materials whose shear strength depends on the local pressure and, therefore, may be characterized by a dimensionless parameter: b = σ comp /σ tens . It was in Section 2 that this parameter depends on the size of samples, i.e., that this ratio is reduced with scaling down the samples. The above-discussed results including those in Figure 19 have been obtained at b = 5, the value typical for the macroscopic samples. However, the surface asperities are the microscale objects and it becomes important analyzing the effect of the b-parameter magnitude on the wear regimes occurring on both PTE and NTE sliding pairs. Figure 21 shows a wear map reconstructed for three b values with b = 1 being the lowest limiting one not achievable on the ceramic materials. The b = σ comp /σ tens parameter was varied by varying the material's tensile strength while keeping its compressive strength constant. This approach allowed accounting for the fact that tensile strength is more sensitive to the internal defects compared to that of the compressive strength. In the case of PTE friction pair it can be observed that despite a five-fold differences in b, the wear-free/continuous wear region boundaries are close to each other as viewed in in the parametric (V crit sl , σ crit Load ) plane. Significant differences are observed only at low contact pressures (V crit sl at b = 1 is by a factor of 1.5-2 higher than that of at b = 5). At the higher pressures (σLoad > 30/40 MPa), all curves get closer and almost coincide with each other. Another feature is that σ max Load is reduced when the b-parameter increases at high applied pressures which is not the case with wear maps obtained at low pressures. This may be due to the enhancing effect of averaged contact spot stress on local fractures, i.e., increasing its contribution to the criterion of local fracture (A14) with b-parameter growth. An analogous trend is for NTE material, i.e., when σ max Load is reduced with the increase in b approximately by the same value as in the case of PTE material. Ceramics are pressure-sensitive materials whose shear strength depends on the local pressure and, therefore, may be characterized by a dimensionless parameter: b = σcomp/σtens. It was in Section 2 that this parameter depends on the size of samples, i.e., that this ratio is reduced with scaling down the samples. The above-discussed results including those in Figure 19 have been obtained at b = 5, the value typical for the macroscopic samples. However, the surface asperities are the microscale objects and it becomes important analyzing the effect of the b-parameter magnitude on the wear regimes occurring on both PTE and NTE sliding pairs. Figure 21 shows a wear map reconstructed for three b values with b = 1 being the lowest limiting one not achievable on the ceramic materials. The b = σcomp/σtens parameter was varied by varying the material's tensile strength while keeping its compressive strength constant. This approach allowed accounting for the fact that tensile strength is more sensitive to the internal defects compared to that of the compressive strength. In the case of PTE friction pair it can be observed that despite a five-fold differences in b, the wear-free/continuous wear region boundaries are close to each other as viewed in in the parametric ( crit  In summary, it can be concluded that the limiting (maximal) values of loading parameters σLoad and Vsl, at which the wear-free wear regime on PTE sliding pairs is still established, are determined mainly by the material's compression strength (σcomp). The sensitivity of the material's shear strength to the mean stress magnitude expressed in the form of b-parameter is of the first order of vanishing compared to that of σcomp, which controls In summary, it can be concluded that the limiting (maximal) values of loading parameters σ Load and V sl , at which the wear-free wear regime on PTE sliding pairs is still established, are determined mainly by the material's compression strength (σ comp ). The sensitivity of the material's shear strength to the mean stress magnitude expressed in the form of b-parameter is of the first order of vanishing compared to that of σ comp , which controls the maximum of σ Load and where wear-free regime is feasible. The latter is true for both PTE and NTE materials.

Conclusions
The results of this numerical study allow us to stress that the key role of CTE is to expose loading condition intervals where ceramic materials can function without catastrophic wear. For material with positive CTE values, the maximal value of sliding velocity, which still provides the wear-free regime of friction, exponentially depends on the applied pressure. In addition, there is a maximum critical applied pressure whose overcoming leads to establishing the continuous wear regime even at low sliding speeds. Changing the CTE sign to negative leads to the qualitative transformation of the boundary between two wear regimes. There is no wear of NTE ceramics at every sliding velocity value, especially if the applied pressure does not exceed some critical value. The value of critical applied pressure for an NTE material is comparable with that of PTE ceramics, provided that the mechanical and thermophysical characteristics of these ceramics are close in their value.
It should be noted that the maximal sliding velocity corresponding to the wear-free regime on the PTE ceramic materials depends not only on the applied pressure but also on the surface temperature. The latter determines the particular values of the material's mechanical characteristics because of their dependence on temperature and, more importantly, the stability of structures and phases. For the considered PTE material (ZrO 2 ), all of the numerically derived values of maximal sliding velocity lie within the temperature range of structure stability, although at high applied pressures they approach its upper boundary.
Given that the NTE material (Sc 2 W 3 O 12 ) is characterized by the temperature stability boundary close to that of ZrO 2 , the actual conditions for establishing the wear-free regime are determined not only by the critical applied pressure but also by the critical sliding velocity. The latter should be found from a condition that the steady worn surface temperature approaches the phase transformation temperature. The "phase stability-defined" dependence V crit sl (σ Load ) takes form the same as in the case of PTE material (Figure 19), although absolute values of V crit sl are by a factor of 1.5 lower because of lower heat conductivity and slow cooling.
The presented numerical results in combination with the results of earlier studies [6,50,59] indicate the feasibility of reconstructing wear maps for both PTE and NTE materials in terms of dimensionless mechanical and thermophysical characteristics. Analyzing these wear maps allow for the establishment of wear regime transitions and key parameters responsible for occurring wear regime transitions. Practical use of such maps should include analysis of possible temperature-induced structural and phase transformations, i.e., stability of the material under preset temperature and loading conditions. NTE materials are used mainly as NTE components in the composite materials in combinations with the PTE ones. The wear reduction effect of negative CTE makes the NTE materials extremely attractive for applications in both low-and heavy-loaded tribological devices-e.g., as components of composite friction (brake) materials with PTE and NTE components. Study of the behavior of such composites during friction requires using full thermomechanical numerical models of surface layers that consider the features of heterogeneous stress, strain distributions, and interface stresses. Further efforts will be focused on studying wear behavior of such composites. Nevertheless, the asperity-based model can be used, albeit with some caution, to estimate the mild wear regime conditions in composites with nanosized PTE and NTE components (provided that the surface roughness significantly exceeds the characteristic particle sizes). The effective properties of the composite should be used in this case.
Another important problem in tribological applications of the NTE materials and NTE/PTE composites is their low heat conductivity while efficient heat removal in highly loaded pairs is necessary for their structural stability. Looking for effective solutions to this problem will be another task of our future investigations.

Data Availability Statement:
The authors declare that all data and methods used in the research are presented in sufficient detail in this paper so that other researchers can replicate the work including all raw data.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Mechanical interaction of the surfaces is the sum of contact interactions of opposite asperities ( Figure 1). The total force of interaction is calculated as follows: Here, the subscript "N" and "τ" denote normal and tangential interaction, respectively. I is the number of asperity, m is the number of surface asperities (m is the same for both surfaces), F i N and F i τ are absolute values of the forces of normal and tangential interaction of i-th asperity of the lower body with i-th asperity of the upper body, and → n x and → n y are unit vectors along X-axis and Y-axis, respectively.
The material of the interacting bodies is assumed to be brittle, which is true with the majority of low porosity and pre-free ceramics, and is modeled as a linearly elastic one. The adhesion between the surfaces is assumed negligible and is not taken into account. Surface asperity dimensions are assumed to be much smaller than the height of the body. That is why the interaction force ratio of opposite asperities are based on classic analytical expression, which establish a relationship between a force homogeneously acting on a circle area of the free surface of an elastic half space and deflection at the central point of this circle area.
In the absence of adhesion-i.e., the normal interaction of asperities-is described with only their reaction to compression loading. Assuming that asperities are small enough, the asperity resistance to compression from the opposite one is determined from a classical linear relationship between the normal compression force (F) homogeneously acting on a circle area S of the free surface of an elastic plane's half space and normal displacement (∆H) at the central point of this circle [70,71]: where p 0 is a homogeneously distributed pressure, R is the radius of a circle area, and E and ν are the Young's modulus and the Poisson's ratio of the material of a half space, respectively. This expression can be formally rewritten in a simpler form similar to resistance force of a linearly elastic spring to compression: where K = πRE/2 1 − ν 2 is the stiffness of a virtual (effective) spring. Normal contacting between two opposite asperities is described as the asperity resistance to uniaxial compression from the force exerted by the opposite asperity. According to Newton's third law, the normal force for i-th pair of opposite asperities may be put as follows: Here, the superscripts "b" and "u" denote contacting asperities belonging to the bottom and the upper surfaces, respectively. ∆H b,i and ∆H u,i are normal displacements (indentations). K b,i N and K u,i N are stiffness values for the bottom and the upper bodies. ∆H i is the total value of indentation. Another interpretation of ∆H i is the overlap of two asperities if considering them (and the half-spaces) as rigid bodies. Using this interpretation, the ∆H i magnitude can be calculated at each moment of time ( Figure 1): ∆H i = Y b,i − Y u,i , where Y b,i and Y u,i are positions of the tops of i-th opposite asperities along the Y-axis. Asperities do not interact (F i N = F i τ = 0) unless they are overlapped (∆H i ≤ 0). The known value of the total overlap of opposite asperities ∆H i is used in (A4) for calculating the current values of both ∆H b,i and ∆H u,i , as well as the interaction force: The values of asperity normal stiffness are defined similar to (A2), as well as (A3) with consideration of the conditions as follows. Firstly, each of the asperities is characterized by its own value of the effective contact area. Secondly, the geometry of this contact area is not specified. Therefore, both K b,i N and K u,i N are computed in accordance with the following equations: where E b(u) and ν b(u) are the Young's modulus and Poisson's ratio of the material of bottom or upper body, respectively, and S i 0 is the minimum from S b,i 0 and S u,i 0 (the effective contact areas of the opposite asperities).
Tangential interaction of the contacting asperities is modeled by two components according to equation as follows: The first tangential force component F i sh is considered as approximation of shear resistance between two half-infinite bodies contacting with dry friction by some area S i 0 on their plane surfaces [71,72]. The equations to calculate F i sh are analogous to Equation (A4). They are formulated in hypoelastic (incremental) form: where K b,i τ = 0.5π S i 0 G b and K u,i τ = 0.5π S i 0 G u are the tangential stiffnesses of the bottom and the top asperities, respectively. G b(u) = E b(u) /2 1 + ν b(u) is the shear modulus of the material of the bottom or top body. V sl is the relative tangential velocity of contacting surfaces (sliding velocity). ∆t is the time increment (time step of the numerical scheme). ∆l b,i and ∆l u,i are the contributions of asperities into the total tangential displacement increment V sl ∆t.
Similar to normal interaction, the known value of the total tangential displacement increment V sl ∆t is used in (A8) for calculating the current values of both ∆l b,i and ∆l u,i , as well as an increment of the shear resistance force: The current value of F i sh is limited above by the current local value of the dry friction force. If F i sh > F i dr = µF i N (µ is the local dry friction coefficient), then set F i sh = F i dr . The second tangential force component F i bend serves for effective taking into account interaction between lateral surfaces of the overlapped parts of opposite asperities in a layer of thickness ∆H i during sliding. This type of interaction is a complex one since it combines localized shear of the overlapped asperity parts, as well as their bending and indentation into the bodies during relative sliding of asperities. Such a complex interaction is effectively described using the following equations expressed in hypoelastic (incremental) form: The combination ∆H i (V sl ∆t/∆x) in (A10) effectively describes an increase in the asperity volume involved into localized deformation with the relative sliding of the opposite asperities. By analogy with shear resistance force, the value of F i bend is bounded above: It is well-known that dry friction work dissipates mainly into heat so that the local heat release power on the contact spot formed by opposite i-th asperities can be obtained from equation: Redistribution of the heat fluxes between two contacting bodies is modeled by numerical solving classical 2D heat conductance equation: where a = k/(cρ) is the thermal diffusivity, c is the heat capacity, k is the heat conductivity, ρ is the material's density, and Q is the local heat source. This equation is solved for each body using the finite difference method on a square 2D grid with the grid size numerically equal to the asperity period ∆x. Boundary conditions are as follows: a heat flux (Neumann) condition −k ∂T ∂y = 0.5q (∆x) 2 is applied to the contact surface while Dirichlet condition (constant temperature T = T 0 ) is applied to the external surfaces of both bodies. The time integration step ∆t was by a factor of 1000 less than the critical one determined from the Courant's condition.
A change in temperature entails a change in the linear dimensions of the surface layers, including a change in the level (and topography) of the surface. Since internal layers of the material are under mechanically constrained condition, it was assumed that the local change in the surface level (the position of the top of asperity along the Y-axis) is determined mainly by the temperature change in the linear size of the first (surface) cell: where T i and T i 0 are the current and initial values of the temperature of the i-th asperity, Y i T and Y i T 0 are the current position of asperity top along the Y-axis (at temperature T i ) and position of this asperity in a cooled state (at temperature T i 0 ), α is the thermal expansion coefficient (TEC). ∆Y i therm = Y i T − Y i T 0 is the thermally induced asperity height change. The effect of temperature on the effective contact area S was not considered (S i = S i 0 ). If the stresses in the contact spot reach a critical value, then asperity experiences fracture. The asperity fracture by the shear mechanism is considered because stressed asperities are under the "compression + shear" loading condition. Ceramic materials are characterized by a strong dependence of the shear strength on the magnitude of the applied pressure. Therefore, the two-parameter Drucker-Prager criterion was used as a criterion for local failure [57,58]: where b = σ comp /σ tens , σ comp , and σ tens are the asperity material uniaxial compressive and tensile strengths, respectively. σ mean and σ eq are mean stress and equivalent stress.
Assuming the simplified description of the stress state of asperity (uniaxial compression plus tangential interaction of shear type), formulation of σ mean and σ eq is reduced to the following form: σ mean = −σ N /3 and σ eq = 1 √ 2 2σ 2 N + 6σ 2 τ , where σ N = F N /S 0 , and σ τ = F τ /S 0 . Note that the Drucker-Prager criterion is extensively used to determine the limiting state of brittle and quasi-brittle materials (in particular, ceramics).
The asperity fracture condition is determined by the following inequality: The criterion (A14) and condition (A15) are calculated for each asperity being in contact with opposite asperity (σ N > 0) because the mechanical characteristics (σ comp and σ tens ) of the opposite asperities may generally differ. The values of σ mean and σ eq are the same for both asperities in the contact.
Fracture of asperity means that some material fragment (wear debris) is separated from the asperity. This leads to decreasing the position of the top of asperity along the Y-axis (wear of asperity) by some value ∆Y i w . The corresponding value of wear volume can be evaluated as ∆W i = ∆Y i w S i 0 . It is assumed that the effective asperity contact area (S b,i 0 for the i-th asperity of the lower body or S u,i 0 for the i-th asperity of the upper body) changes in wear. Since it is hardly possible to predict the asperity contact area change, the new value of asperity contact area was set randomly within the range limited by minimal (S min 0 ) and maximal (S max 0 ) values after each act of wear. The value of ∆Y i w , which determines the elementary volume of wear, generally depends on many factors and can vary over a wide range. It is known from the results of numerous experimental and theoretical studies that the wear of a single asperity can occur through various mechanisms [36,38,41,42,73], as determined by the ductility/brittleness ratio, the magnitude of the applied normal stress, the degree of lateral overlap of asperities, the features of the geometry of asperities, etc.
Since the asperities are considered in a simplified manner (without taking into account the features of their individual geometry), it is assumed that if condition (A15) is met for some asperity, the broken (separated) part of this asperity has height ∆H b,i (for the asperity of the lower body) or ∆H u,i (for the asperity of the upper body). In other words, the asperity wear is localized in the asperity pair overlap (∆H i ), and the height of wear debris ∆Y i w is equal to the contribution of this asperity to ∆H i . Note that in the case of simultaneous fulfillment of conditions (A15) for both contacting asperities, the Y-coordinates of both tops change by ∆H b,i and ∆H u,i , respectively.