Description of sand–metal friction behavior based on subloading-friction model

A subloading-friction model is formulated to describe the smooth transient variation from static friction to kinetic friction, the recovery to static friction after the sliding velocity decreases, and the accumulation of sliding displacement under the cyclic loading of contact stress. In the past relevant studies, however, the model formulation used for simulations is limited to the hypoelastic-based plasticity framework, and the validation of the model is limited to simulations of the test data for metal-to-metal friction. In this study, the formulation of the subloading-friction model based on a hyperelastic-based plasticity framework is adopted. In the fields of civil, geotechnical, agricultural engineering, and terramechanics, the interaction between soils and metals is critical, as reflected in construction and agricultural machinery, foundation piles, and retaining walls. The validity of the model for describing the friction between various sands and metals is verified by simulations of the experimental data under monotonic and cyclic loadings.


Introduction
It is well acknowledged that (1) when a body at rest begins to slide, a high friction coefficient appears first, which is known as static friction; (2) subsequently, the friction coefficient decreases toward its lowest stationary value, which is known as kinetic friction; (3) when the sliding suspends and then restarts, the friction coefficient approaches the static friction coefficient, and a behavior similar to the initial sliding is observed [1][2][3][4][5][6][7][8][9][10]. The magnitude of the difference between static friction and kinetic friction can reach up to several tens of percent, depending on the materials as well as the sliding velocity and stress levels [11]. Hence, a mathematical model that considers these characteristics is required in engineering analyses involving frictional boundaries.
The recovery of the static friction coefficient has been investigated theoretically and modeled via mathematical expressions directly including the elapsed time after sliding ceases [2,3,[6][7][8]10]. However, the evaluation of the duration of sliding suspension is often accompanied by the ambiguity in judging when sliding ceases and restarts, particularly when the sliding velocity varies gradually within a low velocity range. Therefore, the inclusion of the duration time of sliding suspension in the constitutive modeling of the friction law may result in the loss of objectivity [12]. Meanwhile, the variation in frictional properties on the contact surfaces, including the above-mentioned transitional behavior, should be expressed in terms of the contact traction, sliding velocity, and internal variables representing the state and/or history of a sliding process without including the suspension duration.
In the 1970s and early 1980s, the "rate-and-state friction" (RSF) approach was proposed for modeling static and kinetic friction. Dieterich [13][14][15] proposed a function for the friction coefficient based on the time of the stationary contact and the sliding velocity, as well as discussed the mechanics of stick-slip. Ruina and Rice [16,17] proposed a friction law that depends on the state and slip history by introducing internal state variables. Meanwhile, Popova and Popov [18] presented a historical overview of studies pertaining to friction laws proposed in early works by Coulomb and Amontons, as well as subsequent developments leading to the RSF approach. However, the RSF model was formulated in terms of scalar-valued quantities based on only one-dimensional considerations. Therefore, the application of RSF models is limited to one-dimensional linear sliding.
For general three-dimensional engineering analyses, a friction model applicable to multidirectional frictional sliding on a contact surface is required. In this regard, constitutive relations to describe multi-dimensional frictional sliding have been formulated in the framework of elastoplasticity (e.g., Refs. [19,20]). The key idea of the friction models based on the elastoplastic approach is twofold: (1) the decomposition of the sliding displacement into an elastic (reversible or stick/adhesive) component and a plastic (irreversible slip) component to formulate plastic internal variables that are irrelevant to elastic sliding but change with the plastic sliding history; (2) the expression of a yield function in terms of the contact traction and internal variables to define the loading criterion of whether plastic sliding proceeds in multi-dimensional sliding processes. Computational treatment of friction constitutive laws within the finite element formulation has been established [21,22] and is widely used in engineering analyses involving frictional contact.
Within the framework of elastoplasticity theory, the aforementioned fundamental friction behaviors, i.e., the decrease in the friction coefficient from the static to the kinetic friction coefficient via plastic sliding and the recovery of the friction coefficient by the decrease in the sliding velocity, have been formulated in the subloading-friction model [23][24][25] based on the concept of a subloading surface [26]. In this concept, a subloading surface with a shape and a direction similar to those of the sliding-yield surface is introduced, and the plastic (irreversible) sliding velocity is assumed to be induced gradually as the subloading surface approaches the sliding-yield surface. Consequently, the subloading-friction model is endowed with the capability to describe the gradual evolution of plastic sliding as the contact traction approaches the slidingyield limit and hence the accumulation of sliding displacement under the cyclic loading of the contact traction inside the sliding-yield surface. The validity of the subloading-friction model for describing actual frictional sliding behavior was verified by comparison with the test data [23][24][25]. However, the validation of the model in these past studies was limited to the sliding behavior of metal interfaces. Recently, Ozaki et al. [27] demonstrated the applicability of the subloading-friction model to frictional sliding between a rough rubber hemisphere and a smooth acrylic plate.
In the field of civil and geotechnical engineering, the interaction between soils and metals is critical, as reflected in the driving installation process of the pile foundations and their subsidence/loosening under traffic loads, earthquakes, etc. [28][29][30][31]; retaining walls; construction machineries; and the effect of friction between a soil specimen and a testing apparatus in laboratory tests for geomaterials [32][33][34]. In the field of terramechanics, frictional sliding between soils and metals is relevant to the interactions between wheels/blades and soil terrains [35][36][37][38][39][40][41][42][43][44][45]. Extensive experimental research has been performed to precisely measure the friction property of soil-steel interfaces under monotonic or cyclic loading [46][47][48][49][50][51][52][53], and to elucidate the factors affecting the interface friction property, such as the testing methods, materials, sliding velocity, interface roughness, particle geometry, gradation, and grain crushing [54][55][56][57]. Recently, advanced experimental investigations on the frictional properties of the geological materials have been conducted, e.g., micromechanical experiments of inter-granule loading [58], sliding friction tests of the interface between shale rock and dry quartz sand particles [59], and granular friction tests under linearly reciprocating sliding [60]. It is noteworthy that a mathematical model is required for the successful analysis and prediction of problems in these areas. Hence, the abovementioned features of the subloadingfriction model are advantageous for describing the multidirectional and cyclic frictional sliding behavior at the interfaces between soils and metals, which are often encountered in the fields of civil and geotechnical engineering, construction and agricultural machinery, and terramechanics. Currently, however, the validation and application of the subloadingfriction model are limited to frictional sliding between manufacturing materials and have not been done for soil-metal friction.
In this study, the formulation of the subloadingfriction model based on the hyperelastic-based plasticity framework was adopted. The elastic-plastic additive decomposition of the sliding displacement vector is introduced in Section 2. Section 3 describes the hyperelastic equation used to define the relation between the contact traction vector and the elastic (reversible) sliding displacement vector. In the previous subloading-friction models [23][24][25]27], a hypoelasticbased plastic sliding formulation was adopted, in which the rate of elastic sliding displacement is associated with the co-rotational (objective) rate of the contact traction vector via a rate or incremental equation of the hypoelastic law [61,62]. However, the hypoelastic-based framework is limited to the description of frictional sliding involving infinitesimal elastic sliding. Furthermore, it requires a cumbersome time-integration scheme for the rate equation of the hypoelastic law to calculate the sliding displacement vs. the contact traction relationship; hence, meticulous treatment is necessary in the numerical time-integration scheme to guarantee the objectivity of the constitutive relation [63,64]. More importantly, it is well acknowledged that the hypoelastic-based framework has several serious drawbacks in terms of physical and theoretical aspects, namely, energy dissipation and its accumulation within a purely elastic range during cyclic loadings [65,66], and the arbitrariness or non-uniqueness regarding the selection of a co-rotational rate [67]. By contrast, the hyperelastic equation defines the direct relation between the contact traction and elastic sliding displacement, intrinsically assuring the objectivity property as the fundamental requirement for constitutive laws. Its computer implementation is simple and straightforward as it does not require a cumbersome time-integration scheme. Consequently, the hyperelastic-based formulation presented herein offers significant advantages in terms of both physical and numerical aspects.
In this study, we assumed an isotropic friction-yield criterion with tangential associative sliding evolution law, particularly to validate the predictive capability of the subloading-friction model for the evolution and accumulation of sliding displacement under monotonic/cyclic loadings. This is further described in Section 4. The validity of the subloading-friction model for the description of the frictional sliding behavior between sands and metals was verified by simulations of various test data under monotonic and cyclic sliding displacements, and the details are presented in Section 5.

Sliding displacement and contact traction
We consider the sliding displacement of the counter (slave) body relative to the main (master) body, denoted by u , and assume that it is decomposed orthogonally into the normal component n u and tangential component t u to the contact surface in the additive form, as Eq. (1): n n n n n I u u u u u u n u n (2) where I denotes the second-order identify tensor, n denotes the unit outward normal vector of the surface of the main body, and the normal component is expressed as where the minus sign for n u is introduced to yield a positive number when the counter body approaches the main body.
We assume the additive decomposition of the sliding displacement vector u into the elastic sliding displacement e u and the plastic (irreversible) sliding displacement p u in the following form: This additive decomposition exactly holds even for a finite (large) sliding displacement, which is not limited to an infinitesimal range. Subsequently, the elastic and plastic components are further decomposed into normal ( e and The elastic sliding displacement vector e u is directly associated with the current contact traction vector f via the hyperelastic equation described in Section 3. The elastic-plastic additive decomposition of the sliding displacement vector in Eq. (4) holds exactly even in a large sliding displacement, which is in contrast to the fact that the deformation gradient tensor, defined by the ratio of the current (deformed) infinitesimal line-element vector to the initial one, must be decomposed into elastic and plastic components in the multiplicative form [26]. Therefore, an exact mathematical description of the finite frictional sliding behavior involving large nonlinear sliding displacements and/or rotations of the soil-metal interface can be accomplished using the hyperelastic-based plastic sliding formulation presented herein.
The contact traction vector f for the main body is additively decomposed into the normal traction vector n f and tangential traction vector t f , as Eq. (9): (11) where f t means the unit direction vector of t f . The minus sign for n f is introduced to yield a positive number when compressive traction is exerted on the main body by the counter body.
The contact traction vectors f, n f , and t f can be calculated using the Cauchy stress tensor σ applied to the contact bodies using Cauchy's fundamental theorem [26], as Eq. (12):

Hyperelastic sliding displacement
We consider the hyperelastic model that defines the relationship between the elastic sliding displacement vector e u and the contact traction vector f via the elastic sliding displacement energy function  e ( ) u as Eq. (13): As a specific functional form of  e ( ) u , we adopt the quadratic form as Eq. (14): (14) where E denotes the second-order symmetric tensor of the elastic contact tangent stiffness modulus (  T E E ). In the present formulation, we assume a constant elastic modulus. In this case, substituting Eq. (14) into Eq. (13) We assume isotropy for the mechanical properties of the contact surface. In this case, the frictional property is independent of the sliding direction on the contact surface. We further introduce the normalized rectangular coordinate system,  , fixed to the contact surface. Under these settings, the elastic contact tangent stiffness modulus tensor E , with its inverse tensor, is expressed as Eq. (16): where n  and t  denote the normal and tangential contact elastic moduli, respectively. Applying Eq. (16) to Eq. (15) yields e e t t n n e t n t n 1 1 The contact traction and the elastic sliding displacement are described in terms of their components with respect to the rectangular coordinate system as Eq. (18) The substitution of Eqs. (16) and (18) and its inverse relation is where       . It is noteworthy that the hyperelastic equation provides a one-to-one correspondence between the elastic sliding displacement e u and the contact traction f, from which one can calculate f directly by substituting e u . Hence, neither the rate equation of the hypoelastic law in terms of a co-rotational rate of f nor its time integration is required.

Elastoplastic sliding velocity
The subloading-friction model is formulated based on the concept of a subloading surface. The model is composed of the sliding-yield criterion, the hardening/ softening rule of the friction coefficient, the evolution rule for the plastic sliding velocity, and a variable associated with the sliding-subloading surface.

Sliding normal-yield surface and slidingsubloading surface
Let us consider the sliding-yield surface defined by the following friction-yield function ( ) f f with isotropic hardening/softening (Fig. 1): which specifies the criterion of friction yielding, with the friction coefficient  being the function describing the friction hardening/softening, i.e., the variation in the size of the sliding-yield surface. The friction-yield function ( ) f f in terms of the contact traction f for the Coulomb friction law is expressed as The variation in  depending on the state of contact traction and/or the history of the sliding process is Fig. 1 Coulomb-type sliding normal-yield and sliding-subloading surfaces.
www.Springer.com/journal/40544 | Friction defined by the friction hardening/softening rule, which is described in Section 4.2. Next, we introduce the subloading surface concept: The plastic sliding velocity evolves as the magnitude of contact traction approaches the sliding-yield limit. Based on this notion, we introduce a sliding-subloading surface that has a shape and orientation similar to those of the sliding-yield surface and passes through the current contact traction. Hereinafter, the sliding-yield surface is renamed as the sliding normal-yield surface to distinguish it from the sliding-subloading surface. Furthermore, we introduce a scalar parameter, named the sliding normal-yield ratio, denoted by r (   0 1 r ), which defines the ratio of the size of the slidingsubloading surface to the size of the sliding normalyield surface. It is vital to the subloading-friction model as it designates the approaching degree of contact traction to the sliding-normal yield surface. Subsequently, the sliding-subloading surface is represented as At this stage, we formulate the isotropic hardening/ softening rule to define the friction coefficient  and the evolution rule of the sliding normal-yield ratio r, which will be presented in Sections 4.2 and 4.3, respectively.

Friction hardening/softening rule for friction coefficient
The following can be inferred based on the experimental findings: 1) The friction coefficient  transiently gains the maximum value of static friction, and then decreases to the minimum stationary value of kinetic friction. Subsequently, we assume that plastic sliding causes a decrease in the friction coefficient, namely, plastic softening, which corresponds to the transition from static friction to kinetic friction.
2) The friction coefficient gradually recovers its initial value at static friction during sliding suspension. Subsequently, we assume that the duration causes an increase in the friction coefficient, i.e., viscoplastic-like hardening.
Considering the above, we postulate the following isotropic hardening/softening function describing the variation in the friction coefficient  [23,27]: (24) or in the differential form: are the material constants representing the maximum and minimum values of  corresponding to static friction and kinetic friction, respectively. The constant  specifies the rate of decrease in  during the plastic sliding process. The constant  specifies the rate of recovery of  for the suspension duration. The first and second terms in Eq. (24) and their differential counterparts in Eq. (25) are associated with the destruction and reconstruction, respectively, of the engagements of surface asperities. The variations in the sliding coefficient, based on Eq. (25), are shown in Fig. 2.

Evolution rule of sliding normal-yield ratio
Based on the abovementioned hypothesis regarding elastoplastic sliding, the evolution of the sliding normal-yield ratio r must fulfill the following conditions for p  0  u : , 0 : quasi-elastic sliding state 0, 0 1 : sliding sub-yield state 0, 1 : sliding normal-yield state 0, 1 : sliding over normal-yield state Therefore, we apply the following evolution rule: (27) with ( ) U r being the monotonically decreasing function Subsequently, we employ the following specific form of ( ) U r : where  u is the evolution coefficient. The sliding normal-yield ratio r evolves along with the plastic sliding velocity, obeying the evolution rule shown in Eq. (27). Equation (29) implies that the increase in r with the progress of plastic sliding displacement becomes less significant as the value of  u decreases. Furthermore, it is noteworthy that the plastic sliding displacement required for a certain increase in r increases with r because of cot((π/2) ) r , and thus ( ) U r is the decreasing function of r.
Meanwhile, in the case of purely elastic sliding, the sliding normal-yield ratio r is calculated by substituting the contact traction at a current state obeying the elastic sliding constitutive relation, along with the current value of  , into the equation of the slidingsubloading surface shown in Eq. (22), and then solving it for r.
According to Eq. (27) and using Eqs. (28) 2 and (28) 3 , contact traction is automatically attracted onto the sliding normal-yield surface during plastic sliding. Meanwhile, contact traction is pulled back onto the sliding normal-yield surface when the contact traction exceeds the sliding normal-yield surface because   0 r for  1 r , based on Eqs. (27) and (28) 4 . The characteristics of Eq. (28), as illustrated in Fig. 3, are advantageous in numerical calculations, as they effectively reduce the numerical errors in the incremental step analysis.

Plastic sliding velocity
The partial derivative of the sliding-yield function is expressed as Subsequently, based on Eq. (11), I n n f f f I n n t f f f (32) Substituting Eqs. (24) and (27) (33) Fig. 3 Contact traction controlling feature in the evolution rule of r. Contact traction is automatically attracted to sliding normal-yield surface during plastic sliding.
www.Springer.com/journal/40544 | Friction Next, we assume that the direction of the plastic sliding velocity is tangential to the contact surface, and further assume that the plastic sliding velocity evolves in the direction of outward normal to the curve generated by the intersection of the sliding-subloading surface and the constant normal traction plane defined by n constant,  f which results in the tangential associated flow rule (Fig. 1), as Eq. (34): n n n f f (35) and where   and t n are the magnitude and direction of the plastic sliding velocity, respectively.
Substituting Eq. (34) into Eq. (33) yields which are associated with the plastic and the creep sliding velocities, respectively.

Isotropic sliding-yield surface
The traction function for the isotropic sliding-yield surface is expressed as (44) for which the following partial derivatives hold with the reference to Eq. (32): 1 (52)   t f t n n n t t f n n n ( ) 1  u u f I n n n t n t n (55) where  is the Macaulay's bracket, i.e.,   s for  0 s and  0 s for  0 s , where s is a scalar argument. In the coordinate system with base 1 2 ( , , ) e e n , by considering: we have

Application to description of soil-metal friction
The applicability of the subloading-friction model for describing the soil-metal friction behavior is examined in this section. All the test data were one-dimensional, i.e., linear sliding, and they included the data from monotonic and cyclic sliding tests. Seven material constants, i.e., s  , k  ,  ,  , u  , n  , and t  , and the initial value of the friction coefficient, denoted by  0 , were included in the present subloading-friction model. The initial friction coefficient was set as the static friction coefficient, i.e., 0 s    , in all simulations.
The material parameters used in the present model were determined as follows: 1) The normal and tangential contact elastic moduli n  and t ,  respectively, may be sufficiently large because they are relevant to the surface asperities of the contact interface. However, n  and t α do not significantly affect the contact stress vs. sliding displacement relationship provided that their values are sufficiently large. Therefore, a sufficiently large value can be set provided that it allows analysis to be performed stably.
2) The static and kinetic friction coefficients  s and  k (    k s ) were determined, respectively, by the peak and bottom values of the ratio of the tangential contact stress to the normal contact stress. Therefore, these values can be determined directly from the experimental data.
3) The material constants  and  were determined to reproduce the decrease in the friction coefficient due to plastic sliding and the recovery rate of the friction coefficient due to the pause of sliding. They interact mutually in the contact stress vs. sliding displacement relationship. If the experimental data of sufficiently fast and slow slidings are provided, and then the values of  and  can be determined independently.
4) The material constant  u is determined to reproduce the smooth transition from the elastic to plastic state.
These material constants are independent of each other in terms of their physical meanings; therefore, they can be determined definitely, although the model simulation of the contact stress vs. sliding displacement relationship involves their combination.

Monotonic sliding behavior
First, we examine the simulation of the monotonic sliding behavior by comparing it with the test data measured by Uesugi and Kishida [47] using a direct shear apparatus. Dry Toyoura sand with  50 D 0.18 mm (average diameter) was used as the soil specimen. Low-carbon structural steel (American Society for Testing and Materials (ASTM) Specification for Structural Steel A36) plates were used for the metal specimens with three levels of roughness, i.e.,  max 0.6 R μm (smooth),  max 6.8 R μm (medium), and  max 9.8 R μm (rough). Here, max R designates the relative height between the highest peak and the lowest valley along a surface profile over the gauge length [48]. A gauge length  0.2 L mm was used to measure the surface roughness, considering the average diameter of the Toyoura sand (  200 μm). The normal contact stress was 78.4 kPa. The model simulation is shown in Fig. 4, and its material parameters are listed in Table 1.
Next, the simulation of the test result for the monotonic sliding behavior measured by Yoshimi and Kishida [68] using a ring torsional shear apparatus is shown in Fig. 5. Tone River sand with a mean grain size of 0.27 mm (the relative density  r 60% D ) was used as the soil specimen. The metal specimen was the same steel presented in Fig. 4. The surface roughness index of steel was max 3.3 3.6 R   μm (the gauge length  0.25 L mm). The normal contact stress was 105 kPa. The material parameters used are listed in Table 2. Fig. 4 Comparison with the test data by Uesugi and Kishida [47] in monotonic sliding between sand and steel in direct shear test.  [68] in monotonic sliding between sand and steel in ring torsional shear test. Table 2 Material parameters used in the simulation of the monotonic sliding in ring torsional shear test [68]

Cyclic sliding behavior
The simulation of the cyclic sliding behavior was analyzed by comparing the result with the test data obtained by Uesugi et al. [50]. Toyoura sand with   The surface roughness index of the steel was  max 28 R μm (the gauge length  0.2 L mm) for two-way sliding. The normal contact stress was 98 kPa. The model simulations for the pulsating and two-way cyclic sliding are shown in Figs. 6 and 7, respectively. The material parameters used in these analyses are listed in Tables 3 and 4.
Precise simulations were performed for all the test data obtained from the monotonic and cyclic (pulsating and two-way) sliding behaviors up to the finite sliding displacements between sands and metals using the subloading-friction model. Fig. 6 Comparison with the test data by Uesugi et al. [50] in one-way cyclic (pulsating) sliding between sand and steel.

Fig. 7
Comparison with the test data by Uesugi et al. [50] in two-way cyclic siding between sand and steel. Table 3 Material parameters used in the simulation of the oneway cyclic (pulsating) sliding test [50] shown in Fig Table 4 Material parameters used in the simulation of the twoway cyclic sliding test [50] shown in Fig. 7.

Conclusions
In this study, we verified the applicability of the subloading-friction model for describing the frictional sliding behavior on soil-metal interface by comparing its results with the experimental data of monotonic and cyclic sliding behaviors under various conditions. The main features of the proposed model are summarized as follows: 1) The irreversible sliding caused by the change in contact traction below the sliding yield traction was successfully described via the subloading surface concept, i.e., the irreversible sliding rate develops gradually as the contact traction approaches the sliding-yield surface.
2) The smooth transition from static friction to kinetic friction was described.
3) The recovery of the friction coefficient after the sliding suspended, or the sliding velocity decreased to an extremely slow velocity. 4) Not only the monotonic sliding but also the cyclic sliding behavior can be described.
5) The finite sliding behaviors involving large