Rate-, state-, and pressure-dependent friction model based on the elastoplastic theory

Adhesion is one of essences with respect to rubber friction because the magnitude of the friction force is closely related to the magnitude of adhesion on a real contact area. However, the real contact area during sliding depends on the state and history of the contact surface. Therefore, the friction force occasionally exhibits rate-, state-, and pressure dependency. In this study, to rationally describe friction and simulate boundary value problems, a rate-, state-, and pressure-dependent friction model based on the elastoplastic theory was formulated. First, the evolution law for the friction coefficient was prescribed. Next, a nonlinear sliding surface (frictional criterion) was adopted, and several other evolution laws for internal state variables were prescribed. Subsequently, the typical response characteristics of the proposed friction model were demonstrated, and its validity was verified by comparing the obtained results with those of experiments conducted considering the contact surface between a rough rubber hemisphere and smooth acrylic plate.


Introduction
The interactions involving the frictional contact phenomena between solids constitute a fundamental topic in engineering and science. The prediction and control of friction at the maximum, minimum, and optimum levels are required for the design and maintenance of advanced equipment and structures. Under this background, considerable research efforts to understand and control friction have been made in various fields such as tribology, solid dynamics, and geophysics, in a wide range of scales ranging from the molecular level to that of a continental plate. Pioneering contributions such as those of Vinci L da, Amontons G D, and Coulomb C A indicate that frictional sliding has been focused on for a long time, however, several unexplained aspects remain even today. A likely reason for this is that the sliding friction phenomenon depends on the rate and state (or history) of the contact surface. In particular, when two solid bodies in contact slide slowly past each other, an intermittent vibration phenomenon might occur. Such rate-dependent frictional behaviors are referred to as stick-slip motions [1,2], and they can impair the stability of machines and structures. Meanwhile, stick-slip motions are also strongly related to the occurrence mechanism of earthquakes [3].
Since the 1700s, a scalar quantity called the friction coefficient has been used as an engineering indicator for the sliding friction phenomenon. Despite being an empirical parameter, the friction coefficient is used in a wide range of fields. However, the friction coefficient is not a parameter in a strict sense, and it is well known that the rate-and state-dependency of the friction coefficient is quite complex [4−9]. For instance, when sliding between bodies commences, a high friction coefficient is first observed, termed the static friction. Next, the friction coefficient decreases and approaches its lowest stationary value, which is called the kinetic friction. Thereafter, if the sliding stops for a while and restarts, the friction coefficient recovers, and a behavior similar to that pertaining to the initial sliding is reproduced. In addition, in general, the friction coefficient exhibits velocity-weakening or -strengthening behaviors. Therefore, as indicated by Oden and Martins [2], in the analysis of the rate-dependent frictional sliding, including stick-slip motion, instead of a friction model that is simple but deviates from experimental facts, a model that can appropriately describe the relationship between the friction force and sliding displacement is required, even if it is somewhat complex to use.
The "state variable approach," termed the Dieterich-Ruina friction (DRF) law [10−13], is an approach to model the complex friction phenomena described above. In this approach, the rate-and state-dependency are directly introduced into the friction model, and the internal state variables are introduced to express the dependency of the friction coefficient on the relative velocity and contact history. The DRF law was originally developed to describe the seismic dynamic process, and it has been used in the context of frictional contact of rock, as well as that of a wide range of materials such as steel, glass, plastic, and wood [14−16].
Different from the state variable approach, another modeling method, termed the rate form approach, which describes the sliding friction behavior in the same manner as in the elastoplastic constitutive model, considering the relationship between the frictional stress rate and sliding velocity, was proposed in the 1980s. This approach, until recently, has been employed mainly in the domain of computational solid mechanics [17−19]. Since the implementation method of the frictional contact behavior as a constraint condition into the finite element method was established for such a friction model, it has generally been used in the analysis of frictional contact boundary value problems.
Based on the approach of the elastoplastic theory, a rate-and state-dependent friction model, called the rate-dependent subloading-friction model (SF model) has been proposed [20−22]. According to this model, it is possible to express a smooth transition from the stick state to the sliding state. Further, by prescribing an evolution rule for isotropic hardening or softening of the sliding surface (frictional criterion), it is also possible to reasonably express the mutual transition of static-kinetic friction and velocity-weakening of the frictional resistance. Moreover, the anisotropic friction behavior can also be described using the SF model by introducing an orthotropic frictional criterion and rotational hardening [23].
Both the state variable approach and the rate form approach are based on the adhesive friction theory related to the real contact area, and they can describe the fundamental rate-and state-dependent behavior. However, in the case of a soft material such as rubber, because the ratio of the real contact area A r to the apparent contact area Aa is relatively large, the proportional relationship of  r A W (W: normal load) is no longer guaranteed in the high-contact pressure regime [24−29]. As the normal load increases, the distance between the real contact points decreases, and the mutual interference between the contact points slows down the increase in the real contact area, resulting in  a r A A [26,27,29]. Therefore, the friction coefficient decreases with increasing normal load (contact pressure). Thus, a friction model capable of describing these parameters and relationships is required, because the friction of soft materials exhibits not only rate-and state-dependency but also pressure dependency.
In this study, we formulate a rate-, state-and pressure-dependent friction model based on the rate form approach, in which a pressure-dependent frictional criterion can easily be introduced. Further, we demonstrate the typical response characteristics of the proposed model for various frictional sliding behaviors, including stick-slip motion. Moreover, the validity of the proposed model is demonstrated by comparing the test results pertaining to sliding between a rough rubber plate and smooth acrylic plate to the analytical results.
In this paper, superscripts e ( ) and p ( ) respectively denote elastic and plastic components, and the subscripts n ( ) and t ( ) respectively denote the normal and tangential components.
where,  min and  max are the minimum and maximum values of  .  is a parameter prescribing the characteristic length, and  is a parameter denoting a delay time. It is understood that the friction coefficient is a state variable that depends on the sliding velocity, time, and its own state. The first term in Eq. (3) contributes to the reduction of the static friction to the kinetic friction as softening behavior owing to the plastic sliding. The second term in Eq. (3) contributes to the recovery of the friction from kinetic friction to static friction as hardening behavior due to the creep deformation of asperities. These terms lead to a competition between the deterioration and formation of adhesions, inevitably leading towards the velocityweakening of the friction coefficient. On the other hand, kinetic friction laws having a similar physical and mathematical nature as Eq. (3) were proposed by Ostermeyer [30].
Meanwhile, based on the adhesive friction theory, the tangential stress vector can be defined as where a A and r A are the apparent and real contact areas, where respectively.  is the shear strength of the adhesive part. If we define the ratio of the real contact area to the apparent contact area as  r r a / S A A , the Eq. (4) can be rewritten as where r S is regarded as a function of the normal stress n f , the irreversible (plastic) sliding displacement p u , and the contact time t . The consistency condition can be obtained as the material-time derivative of Eq. (5), as follows Regarding the second and third terms on the righthand side of the Eq. (6), we adopt the same evolution law of Eq. (3), as follows where, min S and max S are the minimum and maximum values of r S . The Appendix contains the physical motivation of the derivations for Eqs. (3), (7) and (8).
On the other hand, in general, the shear strength  is a function of the sliding velocity [31], as follows Thus, we can consider that the shear strength does not vary at a certain moment of constant sliding velocity, i.e., at   v 0.

Formulation of friction model
In this section, we describe the formulation of the rate-, state-, and pressure-dependent SF model based on the concept of the evolution law described in Section 2.

Decomposition of sliding velocity
The sliding velocity v between contact surfaces is additively decomposed into the normal component n v and the tangential component t v , as follows where I is the identity tensor. The symbol  denotes the tensor product. Based on the elastoplastic theory, v is assumed to be additively decomposed into the elastic-sliding velocity e v and the plastic-sliding velocity p v , as follows [17][18][19][20][21][22][23]32] First, consider that the elastic part can be defined using the hypo-elasticity property, as follows , ( ) f C v C n n I n n (13) where e C is the contact elastic tensor, and  n and  t respectively denote the contact elastic moduli in the normal and tangential directions to the apparent contact surface. These moduli are, at times, referred to as the penalty coefficients, because the sliding velocity is equivalent to the gap velocity in finite element analyses based on the penalty method [17−19]. The variable f is the contact stress vector applied to a unit area of the apparent contact surface, and  ( ) denotes the co-rotational rate with objectivity.

Normal-sliding and subloading-sliding surfaces
Let Eq. (5) be assumed as one pertaining to an isotropic sliding surface considering the pressure dependency.
Here, based on the concept of SF model within the framework of the unconventional elastoplastic theory, it is assumed that the interior of the sliding surface is not a purely elastic domain, and plastic sliding also occurs by the change in the contact stress inside the sliding surface [20][21][22][23]32]. Thus, we call the Eq. (5) as the normal-sliding surface.
Next, we introduce the subloading-sliding surface, which always passes through the contact stress f, and maintains a shape similar to that of the normal-sliding surface in the tangential stress plane ( t1 f , t2 f ), as shown in Fig. 1(a). The subloading-sliding surface can then be described as where   (0 1) R R is called the normal-sliding ratio. R plays the role of a three-dimensional measure of the degree of approach to the normal-sliding state (R = 1) corresponding to a gross sliding state. Referring to Eq. (6), the material-time derivative of Eq. (14) can be written as  where the first term on the right side of Eq. (15) corresponds to a variation in the subloading-sliding surface resulting from microscopic sliding, the second term corresponds to the pressure dependency, the third term corresponds to the sliding-weakening (deterioration of the real contact area), and the fourth term corresponds to the recovery of the contact state (increase in the real contact area, due to creep). As described in Section 2, the velocity-weakening of the frictional resistance is naturally described by the competition between the third and fourth terms. Note that the velocitystrengthening can be described by separately considering the rate dependency of the shear strength  , as in Eq. (9).

Evolution laws of internal state variables
According to the progress of plastic sliding inside the normal-sliding surface, the normal-sliding ratio R gradually approaches 1. In other words, when the tangential stress increases under constant normal stress, the tangential stress increases almost elastically when the tangential stress is zero, and thereafter, it gradually increases to approach the normal-sliding surface but does not increase further after reaching the normalsliding surface [20][21][22][23]32]. Thus, the evolution law of the normal-sliding ratio can be defined as follows where the function U is the monotonically decreasing function fulfilling the following condition ( Fig. 1 Regarding the pressure dependency, we adopt the following equation: Furthermore, for the third and fourth terms on the right-hand side of Eq. (15), we adopt Eqs. (7) and (8), respectively, which have already been derived previously.

Relationships of contact stress rate and sliding velocity
We assume the following sliding-flow rule for the plastic-sliding velocity where   (> 0) is the magnitude of the plastic-sliding velocity, and it is often termed as the plastic positive proportionality factor or the plastic multiplier [32]. Note that, in this study, the normal plastic-sliding velocity is not taken into account, i.e.,  p n v 0. By substituting Eq. (20) into Eq. (19), the magnitude of the plastic-sliding velocity   can be obtained as Meanwhile, to derive the relationship between the contact stress rate and the sliding velocity, the plastic positive proportionality   must be expressed as a function of the sliding velocity v . The substitution of Eqs. (12) and (13) into Eq. (19) yields which, using Eqs. (11) and (20), can be further rewritten as The magnitude of the plastic-sliding velocity, expressed in terms of the sliding velocity, and denoted by   instead of   , can be obtained as follows: Thus, the sliding-flow rule in Eq. (20) can be written as Figure 2 shows the schematic diagram of relations between the elastic stress increment e df and the plastic relaxation stress increment p df during the hardening and softening processes. When the point "a" is a start point, the line "a-b-c" is a stress cycle during loading and unloading. Here, where du and p du are the sliding increment and the plastic where sliding increment, respectively. Based on Eqs. (12) and (13), the contact stress increment df is given by the sum of the elastic stress increment and the plastic relaxation stress increment. Thus, the formulation of elastoplastic constitutive equation holds Further, Eq. (27) is equivalent to the following equation, which is also given by Eqs. (12) and (13) Consequently, the contact stress rate can be obtained by substituting Eq. (25) into Eq. (28) and considering the relations of Eqs. (10) and (11), as follows Here, if we ignore the rate dependency of the real contact area, Eq. (29) degenerates to the following equation The loading criterion is given by where the judgement of whether the contact stress reaches the sliding surface is not required because the plastic-sliding velocity is induced continuously as the contact stress approaches the normal-sliding surface [20][21][22][23]32]. Figure 3 shows the schematic for a typical variation   of the tangential stress with sliding displacement during a sliding-holding-sliding state. Figure 3 also shows the relevant variations of the normal-sliding and subloading-sliding surface in the ( r S , n f ) plane. When we input a constant tangential sliding velocity under a constant normal stress, first, preliminary microscopic sliding occurs. In this (a) sub-sliding state (R < 1), the contact stress vector lies on the subloadingsliding surface. Further, the subloading-sliding surface expands and the normal-sliding surface shrinks owing to the occurrence of the plastic sliding. Hence, the system demonstrates a smooth transition from the stick to the sliding state. Subsequently, the contact stress reaches (b) the maximum friction state after going through the normal-sliding state (R = 1). Next, the contact stress demonstrates (c) the sliding-weakening (softening) state, while the normal-sliding and subloading-sliding surfaces overlap and shrink together towards (d) the minimum friction state. After reaching the kinetic friction value, the tangential stress is unloaded to zero. During the cessation of sliding ((e) holding state), the normal-sliding surface expands with the elapsed time while R = 0. Hence, if a sliding velocity is provided again, larger recovery of the static friction occurs for a longer holding time.

Concrete functions
As described above, we adopted Eqs. (7) and (8) to prescribe the evolution law of the real contact area. In the following section, we explain other concrete functions for evolution laws, which were applied to the numerical analysis performed in this study. We adopted the following function for the evolution law of the normal-sliding ratio: where r is a parameter, and its dimension is the reciprocal of length. The use of Eq. (34) leads to notable advantages in the numerical simulation [20−23, 32]: 1) The realistic smooth transition from the stick to the sliding state, i.e., the preliminary microscopic sliding displacement, is rationally represented; 2) The accumulation of sliding is described for the cyclic loading process, even for a low amplitude of the tangential contact stress t f ; 3) The contact stress is automatically attracted to the normal sliding-surface in the frictional loading process (i.e., if  1, R then   0 R ) (Fig. 1). This aspect leads to high efficiency, and robustness in the numerical analysis conducted in finite calculation steps. In addition, the return mapping method for integrating Eq. (29) in quickly and accurately is applicable [32].
The transition from the preliminary microscopic sliding to gross sliding is caused by the deformation of microscopic asperities and the non-uniform occurrence of slip on the apparent unit contact area (from partial occurrence to total occurrence). Note that the parameter r is related to the characteristic length of preliminary microscopic sliding. Based on the tangential contact theory and the method of dimensionality reduction (MDR) [33], it was demonstrated that the length of preliminary sliding is purely of contact mechanical nature and can be calculated exactly if the elastic properties and the shape of the contacting bodies are known, -without requiring any fitting parameters [34,35]. Therefore, it is thought that the value of parameter r could be determined by referring to a series of works [34,35].
Regarding the pressure dependency, we adopt the following simple function for the ratio of the real contact area to the apparent contact area S r : where b is parameter, and its dimension is the reciprocal of stress. Then,  r S in Eq. (18) can be obtained as Meanwhile, referring to the results by Lorenz et al. [31], we adopt the following concrete function for the rate dependency of the shear strength  defined in Eq. (9).
where  0 is the value during quasistatic sliding. c and d are parameters. In the case that the influence of velocity strengthening is significant, the Eqs. (29) and (37) are simultaneously used in the numerical simulations.

Numerical analysis
This section describes the basic response of the rate-, state-, and pressure-dependent SF model, obtained by performing numerical experiments for the straight sliding phenomenon without a rigid-body rotation, i.e.,    f f. By adopting a two-dimensional coordinate system, we have For all calculations, we adopted the following model parameters: Smooth transition from the static friction to kinetic friction was observed, as shown in Fig. 4. Further, the velocity-weakening and velocity-strengthening of friction in the lower and higher sliding velocity regimes could be described by the coupling of Eqs. (29) and (37). Figure 5 shows the effect of the preliminary microscopic sliding prior to macroscopic sliding on the accumulation of sliding; here, the parameters and normal stress were the same as those used to obtain the results shown in Fig. 4. We realized the cyclic pulsating tangential stress via sliding velocity control. First, we applied frictional loading by using the constant sliding velocity (  t 0.1 mm/s v ); next, unloading was performed when the prescribed tangential stress was attained. Subsequently, reloading was performed when the value of the tangential stress reached zero. The above series of sliding velocity control was repeated. We set three levels of amplitude of the tangential stress, specifically,  t m a x / ( ) f S = 0.6, 0.7, and 0.8 (Eq. (14)). Figure 5 suggests that the accumulation of the sliding displacement under the cyclic frictional loading below the normal-sliding surface was represented appropriately by the proposed model, even though it could not be predicted by using the conventional friction model (   r ) [20,21,32]. Furthermore, the calculations indicate that the sliding accumulates according to the amplitude of tangential stress. After several cyclic loadings, because the normal-sliding   surface shrinks owing to the occurrence of plastic sliding, gross sliding eventually occurs. Note that the above cyclic behavior includes the influence of the rate dependency, as shown in Fig. 4.
Here, practically, the preliminary microscopic sliding up to the gross sliding depends on the normal load. That is, the length of microscopic sliding increases with normal load [34]. The proposed friction model can also consider this effect. Figure 6 shows the numerical results of the pressure dependency. In the calculation, we changed the normal stress for the sliding displacement, as shown in Fig. 6(a), while the sliding velocity remained constant at

Pressure dependency
In this case, the parameters were the same as those used to obtain the results shown in Fig. 4; however, to clearly observe the tendency of the pressure dependence, the rate dependence was ignored (i.e., c = 0,    ,    ). Further, we set three different values for the parameter b defined in Eq. (35). As can be seen from Fig. 6(b), even when the normal stress changed during sliding, the variation of the tangential stress could be analyzed, accompanying the change in the normal stress. In addition, because the nonlinear sliding surface defined using Eqs. (14) and (35) was adopted, the increase in the tangential stress reduced with the linear increase in the normal stress between points B and C. Here, the degree of nonlinearity could be controlled using parameter b.
As indicated by the dashed lines in Fig. 6(a), we conducted analyses in which the degree of increase in the normal stress between points B and C was changed. The relationships between the traction ratio t n / f f in the steady state (between points C and D) and normal stress are shown in Fig. 6(c). For this analysis, plots were evaluated considering the average values of contact stresses during sliding displacements of 0.3-0.4 mm under steady state conditions. As can be seen from Fig. 6(c), the traction ratio (friction coefficient) decreased with increase in the normal stress, reflecting the gradual reduction and saturation of the increase in the real contact area; this indicates that it was possible to represent the typical pressure dependence of the friction coefficient [29]. Figure 7 shows the response of the pressure dependency when considering the rate dependency. In the calculation, we changed the normal stress for the sliding displacement, as shown in Fig. 6(a), while several levels of the sliding velocity were set. In this case, the parameters were the same as those used to obtain the results shown in Fig. 3. As shown in Fig. 7(a), when the rate dependency was considered, slidingweakening was noted to occur between points A and B; however, thereafter, the tangential stress increased owing to increase in the normal stress between points B and C. However, the degree of increase in the tangential stress was minor compared to that when the rate dependency was ignored, because a competition with the effect of sliding-weakening existed. Subsequently, further sliding-weakening occurred in the steady state (points C-D). Figures 7(b) and 7(c) respectively indicate the relationship between the traction ratio and normal stress in the steady state for cases in which velocity-weakening and velocitystrengthening occurred. As can be verified from the results shown in the figures, both the rate dependency and pressure dependency of the frictional resistance could be expressed using only the combination of Eqs. (29) and (37).

Application to stick-slip motion
The stick-slip motion is a typical rate-and statedependent behavior in the frictional system. To | https://mc03.manuscriptcentral.com/friction demonstrate the applicability of the proposed model for the analysis of stick-slip motion, we performed a numerical analysis using a one-degree-of-freedom system [4,21], as shown in Fig. 8.
When a driver applies the prescribed constant driving velocity V at the spring-end, the slider slides, and the corresponding spring elongation is defined by the equation . Thus, the equation of motion can be described as where M is the mass of the slider, K is the spring stiffness, A a is the apparent contact area, and a and u are respectively the acceleration and displacement of the slider, which represent relative values with respect to the fixed base. For simplicity, we ignored the damping effect. The tangential stress t f was estimated using Eqs. (29) and (37). In this study, we solved Eq. (39) based on the Newmark- method [21,22].
For the calculation, we set the parameters of the dynamic system and the shear strength as follows, considering the values presented in previous studies [4,21,22,29]:  (37), and the other parameters for the friction model were set as the same as those used to obtain the results shown in Fig. 4. Figure 9 shows the variations of the spring force with the elapsed time. It is indicated that even if a harmonic vibration is not assumed as an input value, the typical stick-slip motion under a constant driving velocity can be reproduced. Further, in the case in which the velocity-weakening is dominant, corresponding Fig. 8 One-degree-of-freedom spring-mass system. Reproduced with permissions from Refs. [4,21]. Copyrights Nature, 1994, and Elsevier, 2010. to lower values of c, clear stick-slip fluctuation can be observed. The amplitude and period of the stick-slip motion reduce with increase in the value of c. However, if the value of c is extremely large, the rate dependency shifts to the velocity-strengthening regime, and stick-slip no longer occurs. Furthermore, as shown in Figs. 6 and 7, the effects of not only the rate-and state-dependency, but also the pressure dependency on the stick-slip motion can be examined by using the proposed model. By considering the pressure dependency, the difference in the behavior from that observed when using a linear sliding surface such as a Coulomb type surface, in a higher pressure regime, might be analyzed.
The above results demonstrate that the proposed friction model can describe various frictional behaviors via the unified form, and it can be applied to boundary value problems. Note that the form of Eq. (29) can be directly applied to the finite element method as a constraint condition due to the frictional contact [22,23].

Comparison of numerical analysis results with friction test results
This section discusses the validation of the rate-, stateand pressure-dependent SF model for rubber friction by comparing the obtained results with results of experiments performed under a prescribed sliding velocity and normal load. In particular, we focused on the pressure dependency of rubber friction. Figure 10 shows a schematic of the experimental apparatus used in this study. The apparatus employs the contact between a rough rubber hemisphere made of cross-linked polydimethylsiloxane (PDMS) and a smooth plate made of polymethylmethacrylate (PMMA). The PDMS hemisphere is fixed to a rigid base via a three-directional dynamometer and Z-directional motorized stage. Meanwhile, the PMMA plate is fixed to the rigid base via a motorized X-directional stage. Cross-linked PDMS (Dow Corning's SYLPOT 184) was used to create the rubber specimens. First, a mixture comprising a base prepolymer and cross-linker agent with a compounding ratio of 10:1 was poured into a hemispherical steel mold. Next, the mold was heated at a temperature of 120 °C for 30 min to cure the PDMS mixture. Subsequently, the mold was allowed to cool naturally at room temperature, and after a sufficient duration, the PDMS hemisphere was removed from the mold. Finally, the convex surface of the PDMS hemisphere was polished using sandpaper to introduce surface roughness. Two types of PDMS hemispheres with different surface roughnesses were prepared; i.e., for specimens I and II, the Ra values were 2.11 and 4.95 μm, respectively.

Outline of friction test
The test specimens were washed using ethanol before rubbing. Subsequently, the PDMS hemisphere was pressed to the PMMA plate using the Z-directional stage. After allowing sufficient resting to ignore the increase in the real contact area with the contact time [14,35], the PMMA plate was driven using the X-directional stage at the prescribed speed. During the sliding motion, temporal changes in the normal load W and tangential load Fx acting on the contact surface were recorded using the three-directional dynamometer at a sampling rate of 10 kHz.

Comparison
In the comparison of the experimental and analytical results, the stiffness of the experimental system, including that of the rubber, was ignored. We compared the coupling response of Eqs. (29) and (37) with the tangential load Fx, which is regarded as the friction force. In addition, the apparent contact area was assumed to be constant during sliding, and it was derived from Hertz's contact theory. Using the JKR test, Young's modulus of the PDMS hemisphere was estimated to be 1.2 MPa [36,37], and Poisson's ratio was assumed to be 0.5 [38]. Note that the distribution of the normal stress on the contact surface was also ignored, and we adopted average values of contact pressure according to the normal load W.
The model parameters were appropriately fitted to reproduce all experimental results by one parameter set. As an exception, to reflect the influence of the surface roughness, the parameter b defined in Eq. (35) representing the pressure dependency was set to be 7 and 11 MPa −1 for specimens I and II, respectively. For all calculations, the following model parameters were adopted: It should be noted that velocity-weakening does not occur when using this specific parameter set, as it was not observed in the conducted experiment. Figure 11 shows the variation of the friction force with the elapsed time under four levels of normal load, wherein the sliding velocity was 0.1 mm/s. As an example, the result for specimen I was shown. The figure confirms that the experimental results for various levels of normal load were simulated well by the proposed model. However, the increasing behavior of the friction force at the beginning of sliding was not in agreement because the stiffness and freedom of rubber were not considered in the numerical simulation. Here, the crucial importance of considering the contact stiffness to obtain an agreement between theoretical and experimental results was already noted by Teidelt [34]. Therefore, by implementing the friction model in an equation of motion or a finite element analysis model, the influences of the contact stiffness and freedom of system, Mindlin slip, and the distribution of normal stress can be analyzed.
Next, we verified the applicability of the proposed model to evaluate both the rate dependency and pressure dependency. Figure 12 shows the pressure dependency of the traction ratio in the steady state under three levels of sliding velocity. The abscissa of the graphs denotes the prescribed normal load W. It is seen that the traction ratio obtained using the experiment decreases with increase in the normal load and exhibits typical pressure dependency. In addition, the degree of pressure dependency varies with the surface roughness, reflecting the dullness and saturation of the increase in the real contact area. Furthermore, the traction ratio increases with increase in the sliding velocity and demonstrates velocity strengthening. As shown in the figures, and by comparing the experimental results under the considered conditions with the analytical results, the quantitative predictability of the proposed model was verified.

Conclusions
In this study, to establish a friction model that could easily be implemented into various numerical simulation methods, such as the finite element method, as a constraint condition, we formulated a rate-, state-and pressure-dependent SF model based on the elastoplastic theory.
First, we phenomenologically derived the evolution law of friction coefficient (real contact area) based on the consistency condition. Then, we formulated the friction model by introducing a nonlinear sliding surface. Subsequently, we demonstrated the typical response characteristics of the proposed model for various frictional sliding behaviors, including stick-slip motion. Furthermore, the validity of the proposed model was demonstrating by comparing the analytical results with results for a rubbing test between a rough rubber hemisphere and a smooth acrylic plate.
A limitation of the study is that we focused only on the frictional sliding caused by adhesion of the real contact area. It is widely known that the hysteresis friction owing to the viscoelastic property of rubber is also another critical aspect affecting the friction force. The effectiveness of the proposed model for the case in which both the adhesion friction and hysteresis friction are generated must thus be examined in the future. In such a case, the hysteresis friction could be separately simulated using different frameworks, such as the finite element method and the method of dimensionality reduction [33].