Nonlinear model of ice surface softening during friction

The ice surface softening during friction is shown as a result of spontaneous appearance of shear strain caused by external supercritical heating. This transformation is described by the Kelvin-Voigt equation for viscoelastic medium, by the relaxation equations of Landau-Khalatnikov-type and for heat conductivity. The study reveals that the above-named equations formally coincide with the synergetic Lorenz system, where the order parameter is reduced to shear strain, stress acts as the conjugate field, and temperature plays the role of a control parameter. Using the adiabatic approximation, the stationary values of these quantities are derived. The examination of dependence of the relaxed shear modulus on strain explains the ice surface softening according to the first-order transition mechanism. The critical heating rate is proportional to the relaxed value of the ice shear modulus and inversely proportional to its typical value.


Introduction
Ice and snow friction is of great importance in everyday life, sport, nature and industry [1,2]. The kinetics of ice friction is determined by such processes as adhesion, surface and pressure melting, frictional heating, creep, and fracture [3,4]. According to [5], the question whether temperature or yield stress of ice plays the crucial role at friction is still a subject of some debate. The study [6] for the first time concludes that the reason for reduced friction is a water film produced on the ice surface due to frictional heating. Several investigators, in particular [7][8][9], have extensively developed this idea, since understanding of liquid film formation conditions is necessary for practical applications. Due to [10], the premelting layer is formed with fluctuating domains of liquid water and solid ice that resembles the defect structure. The studies of [11] focused on the thermal conductivity effect that reduces with an increasing sliding velocity.
Consider some aspects of the theoretical models with quantitative estimations. It is generally recognized that ice surface melting during friction does not always take place due to perfect mechanical heating. The results of [7,12] explained this feature by thermodynamic arguments. Within the framework of approach [13], the dependencies of friction coefficient for mild steel, Perspex and copper on ice were described regarding the sliding velocity, melting and ambient temperatures, and thermal conductivity of slider. The paper provided the evidence that an enhanced softening of ice above −2 • C resulted in significant wear. The study [14] extended the Evans et al. theory [13] to the case of hydrodynamic lubrication. Two domains are revealed with different dependencies of friction coefficient on velocity. In the first domain, there is no melting of ice while in the second domain, friction is governed by the water film at the contact area appearing due to frictional heating. The idea of a prevailing thermal control mechanism during ice friction is confirmed in [12]. This survey considers the roles of softening in the wear of rubbing materials, hydrodynamic friction, and squeezing-out of a lubricant film. The investigations [15,16] expound the theory for melt lubrication including the case of squeezing-out. This approach, taking into account both hydrodynamic friction and surface roughness, permitted to obtain the expressions for the layer thickness, coefficient of friction and wear. Colbeck [17] considered that kinetic friction on snow is governed by three components: dry, lubricated, and capillary friction. He constructed several dependencies: 1) the dependence of water film thickness for perfectly insulated and aluminum sliders on temperature; 2) coefficient of the lubricated friction versus velocity at various temperatures; 3) the dependence of water film thickness on the distance along the lubricated area for plastic and aluminum sliders at various temperatures; 4) total friction versus the length along an aluminum slider at various temperatures. These calculations show that friction force influencing the slider changes slightly over a wide range of velocities and temperatures. This theory is restricted to the suppositions that friction is independent of load, exceptional attention is paid to water friction, approximate estimation of the heat current to the slider, etc.
Recently, there appeared atomistic simulations of ice friction including [18][19][20]. Using molecular dynamics method [18,19] and Ginzburg-Landau free energy (Hamiltonian) for the case of first-order phase transition [19], it is shown that premelted ice surface film consists of some molecular layers and its thickness grows with load and temperature. This leads to an increased lubrication and to lowered friction due to the weakening of hydrogen bonds between ice molecules. Since the ascent of sliding velocity results in the frictional heating and, eventually, in the growth of temperature and a softened film thickness, the former is maintained constant due to the introduction of a thermostat in simulations [18]. Thus, the calculated friction force vs velocity dependence increases linearly owing to a viscous component of stress arising in the liquid-like film of ice surface during shear.
The relaxation of the shear component of stress proceeds in course of time [21] τ = η/G, (1) where η is the dynamical shear viscosity, G is the shear modulus. The principal assumption of [22] is as follows: while the kinetic effect of liquid is freezing, the viscosity η becomes infinite at a finite shear modulus G. However, the situation is opposite to the usual second-order phase transition, where an infinite increase of the τ at a critical point is also observed. Actually, changing from viscoelastic liquid into a general case, expression (1) assumes the form τ = χ/γ, where χ is the generalized susceptibility, γ is the kinetic coefficient [in equation (1) these quantities are G −1 , η −1 , respectively] [23,24]. At the phase transition, an infinite increase in susceptibility χ is observed while the kinetic coefficient γ has no singularity. This is equivalent to the shear modulus G tending to zero at a finite viscosity η in equation (1). This situation meets the viscoelastic transition [25][26][27][28][29].
The underlying assumption of our approach is that ice softening during friction is ensured by selforganization of both stress σ and strain ε shear components, on the one hand, and the temperature T , on the other [30]. The relationship between components σ and ε is well-known, with the Kelvin-Voigt model describing its simplest case [31]. The temperature effect is caused by critical increase in the shear modulus G(T ) with a decrease in the temperature: G = 0 in the water, and G 0 in the ice. The governing equations are derived in section 2 considering the above mentioned circumstances. Section 3 presents the study of realization of the conditions of ice surface softening according to the mechanism of continuous second-order transition. The interaction of the mentioned factors results in the onset of the steady state at supercritical value of thermal energy imposed in the surface layer, where the shear strain can take anomalously large values. The lubrication ice friction regime is discussed here, i.e., the model applicable to dry ice friction when the temperature is too low for ice to melt. Section 4 is devoted to the description of ice surface softening by a scheme of discontinuous first-order transition that is observed experimentally in [4]. The exposition in this part is basically different from [30] since the dependence of a relaxed shear modulus on the strain is of another form (31). Therefore, the resultant figures are quantitatively different from those in [30]. Thus, the description of self-organization of adatoms on the semiconductor surface and the ice surface softening within the framework of a similar approach has many differences.

Basic equations
The configuration for which we present a solution, consists of two rubbing planes of ice or planes of ice and of other material (e.g., solid, rubber and so on), separated by a lubricating softened ice layer. It is widely accepted that the relaxation of the shear component of a strain tensor ε in the ice surface layer is determined by the Kelvin-Voigt equation for viscoelastic medium [31] ε = −ε/τ ε + σ/η ε , (2) where τ ε is the Debye relaxation time and η ε is the effective shear viscosity coefficient. The second term on the right-hand side describes the flow of a viscous liquid due to action of the corresponding shear component of the stress σ. In stationary state,ε = 0, equation (2) is reduced to the Hooke-type relationship is the relaxed value of shear modulus (ω is circular frequency of a periodic external effect).
Within the framework of the phenomenological Landau theory [23,25], phase transition is governed by free energy F that is expanded into power series over σ playing the role of an order parameter in study [26]: where G(T ) ≡ G(ω)| ω→∞ is the non-relaxed shear modulus that depends on the temperature, σε implies the external field ε effect, and A is the positive unharmonicity constant. The equilibrium value of σ is determined by the equality where F 0 is free energy at ε = 0. The relaxation transition to equilibrium state is described by the Landau-Khalatnikov-type equation [24,26,32,33 Here, η is a kinetic coefficient, which has the meaning of the shear viscosity. If σ is close to its equi- Hence, the relaxation equation (5) presumes the linear form Here, the first term on the right-hand side describes the relaxation during time τ σ ≡ η/G(T ). In a steady stateσ = 0, the kinetic equation (6) has the form of the Hooke's law Substituting ε/τ σ for ∂ε/∂t in equation (6) reduces it to a Maxwell-type equation for a viscoelastic matter [21].
Note that the effective viscosity η ε ≡ τ ε G ε and a relaxed modulus G ε ≡ η ε /τ ε do not coincide with the real viscosity η and non-relaxed modulus G(T ), respectively. This is caused by a different physical meaning of the Landau-Khalatnikov-type (6) and the Kelvin-Voigt (2) equations [21,26,31]. The values G ε , η, η ε very weakly depend on the ice surface layer temperature T , while the shear modulus G(T ) vanishes when the temperature decreases to T c [22,[34][35][36]. Further, the temperature dependencies are used for the approximation: In order to present the self-organization process [26,[35][36][37][38][39][40], the kinetic equation for temperature T is needed for completing the equations system (2) and (6), which contains the order parameter ε, conjugate field σ, and control parameter T . Employing the approach [30], based on the elasticity theory relationships in [21], § 32, the following equation can be derived: where c p is the heat capacity, κ is the heat conductivity. The last term on the right-hand side stands for dissipative heating of a viscous liquid, flowing under the effect of the stress σ, that can be neglected in 33002-3 the case under consideration. On the other hand, the one-mode approximation (κ/l 2 )(τ T Q − T ) ≈ κ∇ 2 T can be used with acceptable accuracy in equation (9) [21,26,35,[41][42][43]. Thus, we consider the thermal effect of friction surfaces whose value is not reduced to the Onsager component and is fixed by external conditions (l is the scale of heat conductivity, i.e., the distance into which heat penetrates ice, τ T ≡ l 2 c p /κ is the time of heat conductivity): Here, Q 0 is a heat flow from the surrounding solids to the surface layer. The square contribution of the stress is implied to be included in the rubbing surfaces temperature T e = τ T Q. The obvious account of this term leads to a significant complication of the subsequent analysis, though it results only in renormalization of the quantities. Therefore, for our further consideration, the component T e = τ T Q in equation (10) is presumed to be constant. It is noteworthy that during derivation of equation (10) we accepted the equilibrium value of the temperature of ice surface layer T 00 to be equal to zero. Evidently, contrary to the ice surface being heated initially to the temperature T 00 0, the term T 00 /τ T should enter equation (11). This term describes the relaxation of the current temperature of ice surface layer to its equilibrium value T 00 in the absence of the heat flow Q from the background solids.
It is convenient to introduce the following measure units: for the variables σ, ε, T , respectively. Then, the basic equations (2), (6), and (10) are reduced to a form applicable to any viscoelastic medium [30]: where the constant is introduced. The equations (13) -(15) have a form similar to the Lorenz scheme [37] which allows us to denote the thermodynamic phase and kinetic transitions [26,35,36,[38][39][40].

Continuous second-order transition
In general, equations (13) - (15) have no analytical solution, therefore, the adiabatic approximation is used for this purpose: This approach suggests that in the course of evolution, stress σ(t ) and temperature T (t ) follow the variation of strain ε(t ). The minimal relaxation time of strain τ ε is defined by time of reorientations of the water molecules at the freezing point of fresh water 2 × 10 −5 s and τ ε increases by several orders of magnitude at confinement of premelted ice layer between the rubbing surfaces [44,45]. The microscopic Debye time is estimated by relation τ σ ≈ a/c∼10 −12 s, where a ∼ 1 nm is the lattice constant or intermolecular distance, and c ∼ 10 3 m/s is the sound velocity. Therefore, the first of inequalities (17) is valid.
The second condition (17) can be written in the form l L, (18) where the maximal value of the characteristic length of heat conductivity
Then, equaling the left-hand sides of equations (14) and (15) to zero, we can express stress σ and temperature T in terms of strain ε: According to equation (21), at the important interval of the parameter T e = τ T Q > 1 values, the temperature T decreases monotonously with strain growth ε from the value T e at ε = 0 to (T e + 1)/2 at ε = ε m ≡ 1/g . Obviously, the negative feedback of the stress and strain on the temperature in equation (15) leads to this descent. Such an activity is connected with implementation of the Le Chatelier principle for this problem. Indeed, the positive feedback of strain and temperature on stress in equation (14) is the reason for ice melting. Consequently, the self-organization should occur more intensively with the increase in temperature. However, in accordance with equation (15), the system demonstrates such a pattern that the consequence of transition, i.e., ascent of strain, is a drop of temperature, whose growth is the reason for ε increasing. The stress vs strain dependence (20) has the linear Hooke's section at ε ε m with the effective shear modulus G ef ≡ g (T e − 1). At ε = ε m , the function σ(ε) rises to a maximum and at ε > ε m it reduces, which has no physical meaning. Thus, the constant ε m ≡ 1/g corresponds to maximal strain. The growth of a typical value of the modulus G 0 decreases the maximal strain ε m and increases the effective modulus G ef whose value is proportional to the background ice temperature T e .
Substitution of equation (20) into equation (13) produces the Landau-Khalatnikov-type equation [24,32,33,46] τ εε = −∂V /∂ε. (22) Here, the synergetic potential has the form that is reduced after expansion of logarithm over ε 2 to free energy used in [25] for a description of viscoelastic transition of unstructured condensed matter. At stationary stateε = 0, the potential (23) acquires a minimum. When the temperature T e is lower than the critical value this minimum corresponds to ε = 0, i.e., the ice surface is not softened. In the opposite situation T e > T c0 , the steady shear strain acquires a nonzero value increasing with T e growth in accordance with the root law. This causes the ice softening. Consideration of strain proportionally to the thickness of premelting ice layer reveals a qualitative agreement of equation (25) with the results of molecular dynamics simulations and statistical field theory [18,19]. Besides, in line with [18], we reckon that the friction reduces with the temperature growth because hydrogen couplings fail. Using equations (20) and (21) we get the stationary values of stress and temperature: It is noteworthy that, on the one hand, the steady temperature T 0 coincides with the critical value (24) and, on the other hand, its value differs from the temperature T e . Since T c0 is the minimal temperature

33002-5
at which the ice softening proceeds, this statement represents the effect of a negative feedback of stress σ and strain ε on temperature T [see the last term on the right-hand side of equation (15)]. Self-organization occurs in the limit because the sample temperature drops so much. At a stationary state, the non-relaxed shear modulus is equal to the relaxed one Thus, the surface softening is represented in the model by increasing ε and T e because the ice modulus is a fixed quantity.
Two cases can be distinguished by the parameter g = G 0 /G ε . At g 1, that meets the large value of the modulus G 0 , equations (24)-(26) assume the forms corresponding to the "ice (brittle)" limit. The reverse case g 1 (small modulus G 0 ) corresponds to the "strongly viscous liquid"

Discontinuous first-order transition due to deformational defect of modulus
By using the Kelvin-Voigt equation (2), we suggest the validity of the idealized Genki model for the stress σ vs strain ε dependence. It implies the realization of the Hooke's law σ = G ε ε at ε < ε m and the constant σ m = G ε ε m at ε ε m [σ m , ε m are the maximal stress and strain, σ > σ m results in a viscous flow with the deformation rateε = (σ − σ m ) /η ε ]. Actually, the σ(ε) dependence possesses two regions: the first one, Hookean, has a large tilt corresponding to the relaxed shear modulus G ε , then follows a gentler sloping section of the plastic deformation whose slope is fixed by the hardening factor Θ < G ε . Indeed, the described case means the relaxed shear modulus in equation (2) depending on the strain value. For example, the simplest approximation is used in [47,48] for the case of ultrathin lubricant film melting which describes the above mentioned transition from the elastic deformation mode to the plastic one (β is the positive constant). It takes place at a characteristic value of the strain ε p . It is noteworthy that a relationship of equation (30) type has been proposed, for the first time, by Haken [37] in order to represent the rigid mode of laser radiation. This equation was used for denoting the first-order phase transition [26,35,36,[38][39][40]. Experimental dependencies of shear force on displacement for friction of ice on ice demonstrate the similar peculiarities [4] but, as a rule, the plastic section is horizontal. To describe such a behavior, it is necessary to surmise that Θ = 0 and β = 1 in (30): Moreover, the interpretation of ice surface premelting as a plastic act due to first-order transition agrees with the results of atomistic and statistical field approaches, developed in studies [18,19].
Using the adiabatic approximation (17) for the Lorenz equations (13) - (15), where G ε is replaced by the dependence G ε (ε) (31), the Landau-Khalatnikov equation (22) is derived. The synergetic potential differs from (23) by the last term containing the constant α≡ε p /ε s . At a minor value of temperature T e , the dependence (32) is monotonously increasing with a minimum at ε = 0 corresponding to stationary state of ice (curve 1 in figure 1). As shown in figure 1, plateau (curve 2) appearing at  for T e > T 0 c is transformed into a minimum at the strain ε 0 0, and maximum at ε m that separates the minima at the values ε = 0 and ε = ε 0 (curve 3). With the further growth of the temperature T e , the "ordered" phase minimum, corresponding to the ice softened structure ε = ε 0 , becomes deeper, and the height of the interphase barrier diminishes, disappearing at the critical value T c0 = 1 + g −1 (24 where the upper sign indicates the stable softened ice structure and the lower sign denotes the unstable one. At T e T c0 , the shape of the V vs ε dependence is similar to the one at the absence of the deformational modulus defect (see curve 4 in figure 1). With a decrease of the rubbing surfaces temperature T e the interphase barrier disappears at its critical value T 0 c . At the same time, the potential minimum, corresponding to the ice softened structure ε = ε 0 , vanishes too and this phase transforms into solid ice meeting the minimum at ε = 0. It is noteworthy that temperature T 0 c represents the low limit of the range of realization of first-order phase transition T 0 c < T e < T c0 . At the point T 0 c , stable ε = ε 0 and unstable ε = ε m solutions of stationary states equation ∂V (ε)/∂ε = 0 are equal ε 0 = ε m .
The potential barrier, typical of a synergetic first-order transition, manifests itself only at the deformational defect of the modulus. Since it always takes place [4,19,47,48], the studied ice softening represents the synergetic first-order transition. This occurrence is more complex than the thermodynamic phase transition. Obviously, in the latter case, the stationary value of the softened layer temperature T 0 coincides with a thermostat value T e . In this investigation, the T 0 is equal to the critical value T c0 for the

33002-7
A.V. Khomenko, K.P. Khomenko, V.V. Falko  synergetic second-order transition (see section 3). When the modulus defect is examined, the temperature corresponding to the minimum of the dependence (32), presents itself. In accordance with equations (34) and (35), the magnitude T 0 smoothly descends from the value at T e = T 0 c , to 1 at T e → ∞. As delineated in figure 3, the steady temperature T 0 grows linearly from 0 to T c0 , with T e being in the same interval and, after the drop at T e = T c0 , the value T 0 smoothly decreases. The maximal softened ice temperature (36) is lower than the minimal temperature of the background ice T 0 c > 1 (33). As depicted in figure 3, at T e > T 0 c , the steady-state temperature T 0 of the ice surface layer is less than the T e .

Summary
This consideration shows that the ice surface softening during friction is conditioned by the selforganization of the strain and stress shear components, on the one hand, and by the layer temperature, on the other hand. At this, strain ε is the order parameter, stress σ plays the role of the conjugate field, and temperature T acts as the control parameter. The positive feedback of T and ε on σ [see equation (14)] leads to self-organization. The temperature dependence of shear modulus in equations (6) and (8) has a crucial role. The assumption about shear modulus vs strain dependence allows us to acquire the relationships for temperatures of absolute instability of the softened ice layer T 0 c [equation (33)] and stability limit of the solid ice T c0 [equation (24)]. The real temperature of transition, lying in the (T 0 c , T c0 ) range, can be extracted from the equality V (0) = V (ε 0 ) of potentials in various phases. The analysis of equation (24) demonstrates that the softening begins earlier in the systems with a large typical G 0 and minor relaxed G ε values of shear modulus. The kinetic Landau-Khalatnikov equation (22) with the synergetic potential (32) describes this first-order transition. The freezing of softened ice at τ ε = ∞ can occur (ε → 0) even in a nonstationary state ∂V /∂ε 0.
The obtained expressions for temperatures (24) and (33) are of an engineering character because they can be used to predict the friction reduction or increase as well as to remove the negative effect of the interrupted mode of ice friction being the main reason for destruction of the rubbing surfaces. At this temperature interval, the stick-slip mode of friction is possible, characterized by transitions between two dynamic states during the stationary sliding. The latter is observed due to the presence of rapidly fluctuating (in space and time) domains of ice and liquid-like ice [10,27,28,42]. In the coming work, we are aiming to theoretically define the parameters of the rubbing surfaces at which the ice surface layer is in a liquid-like state and the friction between the surfaces decreases or increases. It is planned to combine the contact mechanics, thermodynamic and nonlinear models in order to define the behavior of the system at intermediate velocity interval, where the frictional heating leads to the lowering of friction by either thermal softening of a thin surface layer, or by the formation of nonuniform thin surface film [2]. By joint use of these approaches, it is possible to construe why the ice friction sharply diminishes with ascent of shear rate before reaching the velocity at which a thin homogeneous water layer appears on the ice surface at the expense of frictional heating.