A Fully Coupled Normal and Tangential Contact Model to Investigate the Effect of Surface Roughness on the Partial Slip of Dissimilar Elastic Materials

By adopting a global search method for all shear tractions in the contacting area and using an alternative convergence criterion regarding load balance in the lateral direction, the newly developed algorithm provides stable solutions to partial-slip problems of elastically dissimilar materials. The model is validated via the comparison between the simulation and literature results for a sphere-on-flat problem under fully coupled conditions. It is then employed to investigate the influence of surface roughness parameters including the root mean square (RMS) roughness and RMS slope on partial-slip solutions under coupled conditions. Since the gross sliding condition is modified under coupling effects, the relationship between the ratio of the stick area to the contacting area (stick ratio) and tangential load in the coupled case, unlike that in the uncoupled one, becomes non-linear for all tested rough surfaces. Under low or medium tangential loads, the surface with a higher RMS gradient or a lower RMS roughness experiences more stick regions within the contacting area. This trend then becomes irregular at higher tangential loads.


Introduction
Frictional contact has always been an area of concern within the framework of contact mechanics since the phenomenon of the relative motion being impeded by friction is commonly encountered in practice when two contacting surfaces are subjected to tangential loads. Common examples of such a scenario are the rubber friction in automotive vehicles [1] and fretting wear in orthopaedic implants [2]. When the tangential load is not high enough to induce gross sliding, the interaction zone between surfaces is separated into localised stick and slip regions. The solution to this so-called 'partialslip' problem depends on the properties of both materials in contact. When they are identical as assumed in the wellknown Cattaneo-Mindlin solution [3 4] or where the two materials are incompressible, the normal and tangential contact problems are decoupled. In this case, the normal pressure causes no relative tangential displacement and correspondingly the shear traction does not impose any normal displacement. However, when it comes to the contact of two dissimilar and compressible materials, there exists mutual interaction between the normal pressures and shear tractions. As the corresponding parts of the two surfaces experience different tangential displacements arising from the mismatch of their elastic responses, a normal load could solely cause a certain tangential displacement, which is resisted by friction and consequently results in shear tractions and micro-slip within the contacting areas.
The development of a closed-form analytical solution to such a problem can be troublesome because of the interrelation of the normal and tangential contact problems. A method known as the Goodman approximation [5] was frequently applied in early stick-slip analysis, which neglects the effect of shear tractions on normal pressure but takes that of pressure on tangential tractions into account. The normal contact problem is uncoupled from the tangential one in this way. Due to this limitation, the solution based on such an approximation is only valid for plane (cylindrical) contact 1 3 98 Page 2 of 24 problems in which the material Poisson's ratios are large, according to Hills' later analysis [6]. Spence [7] adopted an iterative scheme to analyse the normal indentation contact of dissimilar materials in detail, where the indenter was limited to be a rigid and flat-ended punch. The solution to a spherical indentation problem was later developed by Spence [8] through a transformation technique. Hills [9] derived a similar coupled solution to elastic cylindrical contact problems and analysed the relevant fretting contact problem under oscillating tangential loading conditions.
Considering the simplified contact geometries in analytical solutions along with their mathematical complexity during applications, recently great efforts have been made to solve the coupled partial-slip problem numerically, for which the boundary element method (BEM) has been widely applied. Based on the Gauss-Seidel method, Willner [10] developed a fully coupled model for the point contact of dissimilar elastic materials. A similar coupled model was developed by Chen and Wang [11] using the conjugate gradient method (CGM). Their model was later extended to analyse multi-layered fretting contact problems with the addition of twisting moments by Wang et al. [12][13][14]. They also investigated the effects of inhomogeneous materials [15] and elastoplasticity [16] on the stick-slip behaviour. Different CGM-based algorithms were proposed by Gallego [17] and Spinu [18,19] to investigate the fretting contact of distinct elastic materials. The former was then further developed by Leroux and Nelias [20] to analyse the stick-slip between a rigid sphere and an inhomogeneous cylinder.
The studies described above each investigated the partial slip of the corresponding smooth elastic surfaces only. Nevertheless, it has long been recognised that surface roughness plays an essential role in contact problems of real engineering surfaces. In terms of uncoupled rough surface contact, fretting wear analysis of sinusoidal rough surfaces was initiated by Ciavarella [21], in which a closed-form analytical solution to the shear traction distribution was proposed. Kasarekar et al. [22] investigated how fretting wear evolves in a rough circular stick-slip contact, where the effect of surface roughness was highlighted in their wear calculation. A similar uncoupled stick-slip contact analysis was conducted by Dini and Hills [23], who focused on the relationship between the roughness and frictional energy dissipation in a rough Hertzian-type contact problem. The tangential stiffness and energy loss in a rough fretting contact were investigated by Eriten et al. [24] using a physics-based frictional contact model, for which the effectiveness of provided solutions is limited to small loads. A study concerning the fretting wear of rough point contacts was conducted by Lehtovaara and Lonnqvist [25], who reported that the ratio of the tangential force to the static friction controls the separation of stick and slip regions under uncoupled partialslip conditions. A similar finding of the linear relationship between the ratio of the stick area to the contacting area and the ratio of the applied tangential load to the static frictional load was reported in the studies of Paggi et al. [26] and Wang et al. [27]. Based on their respective BEM models, the effect of the Hurst exponent of rough surfaces on the transition of the contact state (from full stick to full slip) was analysed by Paggi et al. [26], whilst Wang et al. [27] conducted similar research focusing on the role of the root mean square (RMS) roughness, RMS slope, skewness and kurtosis. Effects of the RMS roughness on the preliminary displacement (the displacement at which full slip starts) in the frictional contact of self-affine fractal surfaces have also been analysed by Grzemba et al. [28], where the preliminary slip of surfaces with a longer cut-off wavelength was found to have an order of magnitude of the RMS roughness multiplied with the coefficient of friction.
The reported research on the coupled partial-slip contact of rough surfaces is relatively limited compared to the uncoupled approaches as discussed. Chen and Wang [29] simulated the partial slip of a rough sphere against a flat half-space, where the rough surfaces were generated using a digital filtration algorithm. The relationship between the surface roughness and the static friction coefficient was analysed using their coupled stick-slip algorithm. Similar to Chen and Wang [29], Bazrafshan et al. [30] included surface adhesion in the partial-slip contact problem and reported the effects of adhesion and surface roughness on the contact solution. It was reported in their study that the surface adhesion increases the pre-sliding distance and the static friction, whilst the surface roughness only affects the former factor. To date, the contact of two dissimilar elastic rough surfaces remains an under-exploited field of study. Few studies shed light on the role of surface roughness in the separation of stick and slip regions under coupling circumstances. This is the problem we intend to address in more detail in the current paper.
Inspired by the coupled partial-slip algorithm of Chen and Wang [11], a three-level iteration model is developed to solve non-conformal contact problems of two dissimilar elastic bodies with arbitrary geometries and surface roughness. Within our framework, the shear traction in the stick region must be less than the local Coulomb friction. CGM is used as the relaxation method in the model considering its guaranteed convergence for quadratic optimization problems. A global adjustment method for all shear tractions within the contacting area is adopted during the design of the CGM-based algorithm. This novel approach to updating contact tractions helps to avoid oscillating results at the interface between stick and slip regions. The Discretised Convolution Fast Fourier Transform (DC-FFT) is used in the algorithm to accelerate the relevant convolution operation. The combination of these two computational techniques (DC-FFT and CGM) has been widely applied in many numerical modelling of contact problems in the past few decades [31]. To avoid convergence issues arising from the lack of enforcement of load balance in the tangential direction, which could be encountered in traditional models [11,12,17,32], a new alternative criterion was proposed here to improve the stability of the algorithm. The model is well validated by simulating a sphere-on-flat problem under coupled partial-slip conditions, to which the relevant solutions have been reported previously by Chen and Wang [11]. It is then extended to simulate the coupled partial slip of rough surfaces. By comparing the modelling results with our former uncoupled study [27], the influence of surface roughness on the contact solutions is investigated when coupling effects are included.

Methods
On the basis of the coupled relationship between surface displacements and contact tractions, the following fully deterministic partial-slip contact model has been developed, which refers to Chen and Wang's [11] CGM algorithm for the point contact of dissimilar materials.

Problem Formulation
A typical partial-slip contact problem of a smooth sphere (Body 1) against a flat elastic half-space (Body 2) is illustrated in Fig. 1a, in which the material properties of the two bodies ( E 1 , E 2 , 1 , 2 ) are dissimilar. The x-and y-axes are positioned laterally on the surface, whilst the z-axis is perpendicular to those and normal to the surfaces. The sphere behaves as an indenter to compress the half-space under a normal force (W). Meanwhile, the tangential loads ( F x , F y ), which are parallel to the x-and y-axes, respectively, are applied in the plane of contact aiming to generate a relative movement between the two bodies. The surface interaction resulting from the input load brings about contact tractions (q x , q y , p) and surface displacements ( u x , u y , u z ). The relations between these tractions and displacements at any point of the surface are summarised by Wang and Zhu [31] as follows: in which the operation '*' denotes the continuous convolution and G ij (i, j = x, y, z) denotes the Green function. Physically, the Green function G ij characterises the displacement on the i-axis caused by the unit traction on the j-axis and mathematically, it is affected by the material properties of the two contacting bodies and the position of the surface point being considered. The continuous formulas to determine all the relevant Green functions are given by Wang and Zhu [31] and shown in Appendix 1.
For the sake of simplicity, contact is assumed to be purely elastic. Besides, Coulomb's friction law is applied and the coefficient of friction f is assumed to be constant in the following formulation and modelling work. According to this friction law, every point within the contacting area is either in a stick state, in which no relative motion occurs for the corresponding parts of the surface and the shear traction is of magnitude less than ' f p ' or in a slip state, in which there exists relative tangential motion and the shear traction opposes the instantaneous direction of such slip motion Fig. 1 a Geometry of the contacting bodies: W, F x , F y are the input normal load and tangential loads in x-and y-directions, respectively, and E and ν are the elastic modulus and Poisson's ratio of the material, respectively, with subscripts denoting the independent bodies, and b displacement condition on the x-z plane: h i , h denote the gaps between the unloaded surfaces and deformed surfaces, respectively, s x is the slip distance in the x-direction, δ is the rigid body displacement, and u is the surface deformation with subscripts denoting the direction of relevant vectors with magnitude ' f p '. It should be noted that partial slip is intrinsically an incremental process, to which the solution depends on the loading history. The separation of the stick and slip regions is expected to be identified based on the change in magnitude and direction of tangential loads. In this way, the state determined for any point within the contacting area on a loading curve relies on the previous solution despite a purely elastic analysis being conducted. In such context, one of the prerequisites for the effectiveness of the partial-slip solution provided in the manuscript is that the load (normal or tangential) is applied at a sufficiently slow speed such that the stick-slip contact problem can be treated as a quasi-static process in which inertial effects are negligible. Hence, the system experiences a series of equilibrium states and the contact problem can be reduced to an equivalent static problem. In this case, the linear theory of elasticity can be applied and the final state of the surface points being assessed can be determined by applying the load in one step. Such an assumption is valid when the loading is monotonic and proportional. As the conducted stick-slip analysis is not load path dependent, the problem of the dissipative friction and its irreversibility, which require an incremental application of the oscillating tangential load that occurs in fretting contacts, is not addressed in the paper. It is noted that although the assumptions including a purely elastic contact and a constant coefficient of friction are commonly implemented in numerical studies [10-12, 17-19, 33], they may not be suitable when addressing more realistic contact problems such as the junction growth derived from the elastic-plastic contact of surfaces [34,35].
To solve the contact problem based on those assumptions described above using BEM, the potential contact area between the two bodies is meshed into equally spaced rectangular elements with the size of 2a × 2b and number of N 1 and N 2 in x-and y-directions, respectively. The control points for the established grid are the centroids of the elements. To assume the element size is small enough that the contact traction and surface deformation are constant inside any element, the fully coupled relation between the tractions and displacements can now be expressed in the following matrix form: where C ij ( i, j = x, y, z ) is an influence matrix and has a similar physical meaning as that of the Green function but in a discretised form. Via the full solution of the Boussinesq-Cerruti integral equations [36,37], the closed-form solutions of these influence coefficient matrices are provided by Ghanbarzadeh et al. [38] and given in Appendix 2. (2) For the discretised partial-slip contact problems under coupled conditions, a sequence of equations and inequalities summarised by Johnson [39] must be satisfied to obtain valid contact solutions:

Static force equilibrium:
Sum of contact tractions must be equal to the specified force in the corresponding direction: where Δ(Δ = 4ab) is the area of each mesh element, I c is the domain of contacting area and i and j here are the indices of the element (i,j) ( 1 ≤ i ≤ N 1 and 1 ≤ j ≤ N 2 ).

Geometrical condition of surface displacement:
In the normal direction, the following equation should be satisfied: where I p is the whole computational domain and h i , h and z are the initial surface gap before loading, the gap between contacting surfaces under loading and the rigid body indentation, respectively, as illustrated in Fig. 1b.
In lateral directions, the following equation should be satisfied: where s x and s y denote the slip distances in x-and y-directions, respectively, whilst x and y denote the tangential displacements in x-and y-directions, respectively, as shown in Fig. 1b. 3. The following contact conditions must be met for the two surfaces: In the normal direction, the Kuhn-Tucker complementary conditions should be satisfied: where I p − I c is the non-contacting region.
The boundary conditions (Eq. 6) indicate that problems of surface adhesion (p < 0) and penetrable surfaces (h < 0) are not accounted for in the current study.
In tangential directions, the norm of shear tractions in the stick zone is limited to always be less than the local friction (product of the coefficient of friction f and nodal pressure) and there should exist no slip distance as specified by Coulomb's friction law. As to shear tractions in the slip zone, the norm should be equal to the local friction. In addition, the directions of the shear stress and slip distance should be opposite. Any point within the contacting domain I c must be in one of the two states (stick or slip) and such a condition is expressed by the following equation: where I s is the stick domain and I c − I s is the slip domain.
It is noted that although Coulomb's friction law is applied in our model as well as Chen and Wang's model, the actual formulas of the complementary conditions used in the two algorithms are different. The norm of shear tractions in the stick zone is allowed to be equal to the local friction as long as the nodal slip distance vanishes in Chen and Wang's model [11]. A discussion on the different complementary conditions and the potential conflicts caused by the equal sign was given in our former uncoupled partial-slip study [27]. To facilitate the design of the later algorithm, we have maintained the same complementary condition applied in our former research.

Algorithm Description
To determine the contact traction distribution and surface deflection in the contacting area, of which the location is not known a priori, iterative methods are frequently applied during the development of contact models. Considering that the mathematical form of the established linear system deriving from Eqs. 4 and 5 tends to be a numerical optimization problem, where the surface gap h and slip distance s are the system residual to be minimised, the linear search method CGM is adopted to develop the partial-slip contact model. An overview of the developed 3-level iteration algorithm based on the mutual adjustment between normal pressure and shear traction is shown in Fig. 2. The normal pressures and shear tractions are determined in two separate algorithms (normal contact solver and tangential contact solver). As illustrated in Fig. 2, the search for contact solutions starts with the determination of the normal pressure, where the shear traction vanishes during the normal contact analysis in the first iteration of the outmost loop. Hence, the contact problem being addressed here becomes decoupled temporarily. The CGM-based algorithm employed to search for the pressure solution is that given by the classic model of Polonsky and Keer [40]. Similar to the Goodman approximation [5], the obtained normal pressure is then used as the input to determine the local friction ' f p ' at each node and to distinguish the existence of slip from that of stick following Coulomb's friction law. The search of shear tractions within the stick and slip regions is executed by the two-level iteration tangential contact algorithm, and it is detailed as follows: Fig. 2 Overview of the algorithm for the coupled partial-slip solver 1. Determine the components of surface displacements in the x-and y-directions arising from the normal pressure: Estimate the rigid body displacement in the x-and y-directions based on Johnson's analytical solution [39]: where 0 is the tangential rigid body translation when the uncoupled contact is under gross sliding conditions Initialise the shear traction q , the search direction d and the square Euclidean norm of the system residual R: It is noted that the initialisation for nodal shear traction is necessary for each iteration of the outmost loop to obtain pressure convergence. The coupling effect changes the contact area and the range of the shear tractions covered must correspond to this change.

Determine the complete surface displacements in the
x-and y-directions, the slip distance for each node and the square Euclidean norm: where the slip distance in the stick region is the system residual and its square Euclidean norm is to be minimised (approximated to zero) by the CGM algorithm. This is expressed as follows: If R < 1 is the inner loop of the tangential algorithm stops and then proceeds to step 7. 4. Determine the search directions ( d x and d y ) as well as search step ( ) for the shear tractions and modify the shear traction using these two variables. Two different methods are commonly used to search for valid shear tractions. The first approach is to only adjust shear tractions in the stick region and maintain those in the slip region (search direction for the nodes (8) u xz = IC xz * p, u yz = IC yz * p, in slip region is 0), such as the algorithms of Spinu [18] and Chen and Wang [11]. The second one is to conduct a global adjustment for all the shear tractions in the contacting region, for example, see the algorithms of Gallego et al. [17] and Wang et al. [12]. It was reported by Wang et al. [12] and verified by the authors of this manuscript that the latter method provides a more stable solution. A smoother transition from the stick to slip regions was achieved in this way. Inspired by the method used in Wang's layered contact model [12], the global search direction and step for the shear tractions are calculated as follows: Shear tractions within the contacting area are adjusted in the following way:

Check the complementary conditions: If
√ q x (i, j) 2 + q y (i, j) 2 > f p(i, j) , the node (i, j) is moved from the stick domain to the slip domain. The nodal shear traction is updated to enforce the norm of the shear traction to be equal to the local friction: 6. Check the convergence criterion for the shear traction: Or else, move back to step 3. 7. Determine the sum of the shear tractions and check the convergence criterion for the loads in the tangential direction: , the iteration stops, or else, move to the following step. 8. Update the rigid body displacements in the x-and y-direction as follows and move back to step 2: It is noted that in some contact problems the convergence for the outer loop regarding the load may be difficult to achieve as the load balance is not enforced in this tangential contact algorithm. The unguaranteed convergence speed of this loop relies heavily on how the updated shear tractions and rigid body displacements affect each other. An alternative to the convergence criterion for the load balance is to compare the shear tractions after the updated rigid body displacement ( q 1 x , q 1 y ) and those before the relevant adjustment ( q 2 x , q 2 y ), expressed as follows: By employing Eq. 20 as the substitute convergence criterion, the convergence speed for the outer loop of the tangential contact solver is improved and the load balance is satisfied.
As indicated in Fig. 2, pressure convergence is checked as follows after the determination of the shear tractions: In the first iteration of the outmost loop, the old pressure p old is the pressure initialised during the normal contact analysis, in which the whole simulation domain is assumed to be in contact with a uniform pressure of W N 1 * N 2 * Δ . And the updated pressure p new is the output of the normal contact analysis. Hence, the pressure convergence can hardly be satisfied. When it comes to the second iteration of the loop, the variable p old becomes the p new in the last iteration. Shear tractions obtained from the tangential contact analysis in the last iteration are then used to determine the resulting component of the normal displacement, a new contact pressure profile p new is subsequently determined in the normal contact analysis with the updated surface displacement. Afterwards, a new tangential contact analysis is conducted with the updated pressurep new . This whole process keeps running until the convergence criterion of the outmost loop (Eq. 21) is satisfied.

Results and Discussion
The non-conformal contact model based on the algorithm described above was validated by comparing the output solutions with the available literature in the case of a smooth contact problem (a sphere against a flat halfsurface) under coupled partial-slip conditions. It was then extended to multi-asperity contact problems (a nominally flat surface against a flat half-space) to analyse the effects of surface roughness on the separation of stick and slip regions.

Model Validation
To validate the developed model, the partial-slip contact between a rigid smooth sphere and a flat elastic half-space (carbon steel) was simulated here, to which the solution is compared to that derived from Chen and Wang's model [11] using the identical material parameters as shown in Table 1.
The whole computational domain ( 3a 0 × 3a 0 ) is discretised by a mesh system with 256 × 256 uniformly distributed elements. The coefficient of friction f is set to be the value of the Dundurs constant , which is a parameter characterising the level of dissimilarity of two contacting materials as given by [11]: It is of note that the ratio of the frictional coefficient to the Dundurs constant ( f ∕ ) plays an essential role in partial-slip analysis. However, in the following contact simulation, this parameter is maintained to be unity as the study of its effects is beyond the scope of the current paper.
To simply the contact problem, only a tangential load in the x-direction was applied to the two contacting bodies. The simulation results presented are non-dimensionalised by the Hertzian solutions to the sphere indentation problem, where the contact tractions and nodal coordinates are normalised by the peak pressure (p 0 ) and contacting radius (a 0 ) , respectively. The evolution of stick and slip regions with increasing tangential load is visualised in Fig. 3. The distribution of the stick region and slip region under no tangential load is similar to that of an uncoupled case, on which occasion the central stick region is surrounded by a slip annulus. Nevertheless, the stick region starts to recede and it tends to shift to the rear part of the contacting area with increasing tangential load. This transition is peculiar, in that the decreasing stick region keeps moving in the direction opposite to the input load when the load is relatively low and up until it reaches the edge of the contacting area when the load is approximately 0.6 f W . Once the input tangential load exceeds this value ( 0.6 f W ), the stick region shrinks dramatically in the opposing direction until gross sliding is reached.
Such separation patterns of the stick and slip regions can be expected based on the corresponding distribution of normal pressure shown in Fig. 4a. The area in which comparatively low pressure exists tends to slip as it is easy for the shear traction to reach the limit of local friction following Coulomb's law of friction. The distribution of stick and slip regions is still of the uncoupled type when only normal load is applied as the pressure distribution in this case is symmetric about x-and y-axes). However, the increase of the specified tangential load breaks down this symmetry such that it reduces the pressure at the contact edge in the force-leading direction but enhances that in the opposite direction. Hence, shrinking stick regions tend to appear at the contacting edge in the friction trailing direction before vanishing.
Regarding the distribution of shear tractions in the x-direction shown in Fig. 4b, apart from the same loss of axial symmetry, there exist both positive and negative shear tractions under smaller tangential loads ( F x ≤ 0.6 f W ) as shown in Fig. 4b. Globally positive shear tractions are obtained when the load is higher than 0.6 f W.
The dynamic response of the stick ratio ( ing tangential load is shown in Fig. 5. The load required for a vanishing stick area is found to be lower ( 0.9463 f W ) when compared with the uncoupled case. This mainly results from the addition of shear traction in the y-direction when identifying the stick and slip region. The gross slip obtained in advance (before f W ) and the existing micro-slip under normal load alone illustrated in Fig. 5 are two essential contact phenomena arising from the coupling between shear traction and normal pressure. It is noted that the difference between uncoupled and coupled solutions to smooth sphere problems depends on the level of the dissimilarity of the two contacting materials (Dundurs constant). As labelled in Fig. 5, the non-linear relationship between the stick ratio and ratio of the tangential load to static friction is following an analytical solution under uncoupled circumstances. Assuming the same coefficient of friction, a material combination with a higher Dundurs constant will exhibit a more significant difference between uncoupled and coupled solutions. Good agreement can be found between the contact solutions including the pressure, shear traction and stick ratio derived from the present model and those from Chen and Wang [11], as indicated in Figs. 4a and b and 5, respectively.

Rough Surface Study
The current study on the effects of surface roughness on the separation of stick and slip regions under coupled circumstances is an important extension work to our former study which concerned uncoupled situations [27]. Amongst the known roughness parameters, RMS roughness R q is firstly investigated considering that it is one of the factors frequently employed to quantify the surface roughness and there exist several studies reporting the essential role of RMS roughness in determining the contact area [41,42], surface adhesion [43,44], contact friction [29,30,45] and material wear [46]. Besides, as reported in our former uncoupled partial-slip study [27], the RMS gradient (g) was found to influence the contact region and stick region, where explicit mathematical formulas were given. Analytical solutions linking the real area of contact to the RMS  gradient were also proposed by Pastewka and Robbins [47] and Muser [48], respectively. Thus, the effect of the RMS slope on the coupled partial-slip solutions is investigated here as well.
To facilitate the analysis of additional coupling effects in nominally flat-flat contact problems, the rough surfaces used in the following simulations are identical to those that were used previously. The roughness parameters of the twelve surfaces analysed in this study are given in Table 2. The relevant surface profiles are shown in Appendix 3. The reader could refer to Ref. [27] to review the detail on the fractal method used to generate them.
Note that in our former uncoupled surface study [27] the nominally flat rough surfaces are assumed to be incompressible to decouple the normal and tangential contact problems. The coefficient of friction used in the former uncoupled case is different from that used in the current study, as a coefficient of friction ( f = 0.2857 ) related to the Dundurs constant is intended to be used here to facilitate the later comparison between smooth and rough contact problems. Therefore, instead of comparing the magnitudes of simulation results, the trend difference of those contact responses with respect to surface roughness parameters is highlighted when comparing solutions to contact problems of rough surfaces under coupled and uncoupled conditions.

Contact of Rough Surfaces Under Constant Loads
To investigate the role played by RMS roughness, partialslip problems of rough surfaces with the same RMS slope but different values of RMS roughness (rough surfaces 1-6 given in Table 2) were simulated. The material properties and load input are given in Table 3.
As shown in Fig. 6a, the contacting areas tend to remain almost unchanged with varying RMS roughness, whilst the relationships between the stick zone and RMS roughness are different in the coupled and uncoupled cases. Compared with the uncoupled solution, the stick region in the coupled case decreases with the increasing RMS roughness significantly, which leads to the decline of the stick ratio ( with the increasing RMS roughness as shown in Fig. 6b. It must be highlighted that those surface parameters including the RMS slope, RMS roughness and cut-off frequency ( q 0 ) are symbiotic. When the RMS slopes of the rough surfaces 1-6 were adjusted to be identical, the cut-off frequency q 0 of each surface was changed simultaneously. Therefore, it is the synergy of the RMS roughness and the cut-off frequency that really determines the contact area. A similar conclusion was reached in the recent theoretical work regarding rough surface contact by Violano et al. [49] and Ciavarella [44].  Although the interaction between shear traction and normal pressure is considered in the current coupled partial-slip simulation, the coupling effects on the contacting area tend to be trivial, which can also be observed in the former study on the smooth sphere contact (Figs. 3 and 4). This explains the insignificant variation of contacting area with the increasing RMS roughness under coupled conditions. On the occasion when the coupling effect is absent, the variation of the contacting region in stick is insensitive to the change of the surface roughness as shown in Fig. 6b, where the stick ratio fluctuates around 0.5. Different from the response of the contacting area, the separation of stick and stick regions is more vulnerable to the coupling effects since the contact condition of gross sliding and the distribution pattern as well as the magnitude of shear tractions are now altered significantly.
Regarding the study on the effect of RMS slope (g) on coupled partial-slip solutions, the contacts of these surfaces (rough surfaces 7-12 as given in Table 2) with the identical RMS roughness but different values of RMS slope are simulated under the inputs given in Table 3. The relationship between the ratio of the contacting or stick area to the computational domain and the RMS slope is observed to be an exponential decay as indicated in Fig. 7a. The trend lines using exponential decay were found to perfectly fit the plotted data following the theory of least-squares regression, where r 1 and r 2 denote the ratio of the contacting nodes and stick nodes to the whole computational domain, respectively, and R 2 is the coefficient of determination characterising how close the data are to the fitted regression line globally. It is of note that trend lines, which were generated using quadratic fit and perfectly fit the uncoupled data, were also reported in our former uncoupled rough surface study concerning the RMS slope. The reader could refer to Ref. [27] to review the relevant equations. As mentioned above, concerning the explicit mathematical relationship between the contacting area and RMS slope, a well-validated formula for uncoupled rough sphere indentation problems was proposed by Muser [48] based on the parameter-free equations developed by Pastewka and Robbins [47]. The increasing RMS slope of rough surfaces employed here is achieved by increasing the cut-off wavevector q 0 , whilst keeping the RMS roughness constant, which again explains that the change of contacting area results from the aforementioned synergy of the RMS roughness and cut-off frequency. As to the stick ratio, it increases significantly with the growing RMS slope in the coupled case compared with the relatively insensitive response in the uncoupled case as shown in Fig. 7b.
These results suggest that rough surfaces with different RMS roughness or RMS slopes can have different responses to the coupling condition, leading to different stick ratios under the same loading circumstances. To visualise such behaviours, partial-slip solutions of rough surfaces 1, 5, 7, and 10 are presented here to discuss in detail (note that surfaces 1 and 5 have the same RMS slope but different RMS roughness, whilst surfaces 7 and 10 have the same RMS roughness but different RMS slopes).
As shown in Fig. 8, lateral sizes of asperities of the surfaces 5 and 10 are significantly smaller compared with that of surfaces 1 and 7. As the highest asperities of rough surfaces always come into contact first when the normal load is imposed, relatively large contacting islands are observed for surfaces 1 and 7 ( Fig. 9a and c, respectively). On the contrary, the height distributions of surfaces 5 and 10 lead to smaller and more uniformly distributed contacting spots as illustrated in Fig. 9b and d, respectively.  Under the same tangential load, the separations of stick and slip regions within the four different contacting areas are shown in Fig. 10. Although uneven pressure profiles are usually observed for rough contact problems, surfaces 1 and 7 exhibit more pronounced non-uniform pressure distributions compared with those of surfaces 5 and 10. Two regions are selected on each surface to zoom in and to have a detailed analysis. For surface 1, the pressure in region 1 is significantly higher than that in region 2 shown in Fig. 11a and b, respectively. Since the slip tends to occur in the area where a lower local pressure profile is experienced, region 1 is still in a partial-slip state (stick ratio is 0.1728 in region 1), whilst region 2 tends to almost reach gross sliding (stick ratio is 0.0131 in region 2) as indicated in Fig. 11c and d, respectively. On the other hand, for surface 5 which exhibits the asperities that are comparatively smaller and distributed more intensively as shown in Fig. 8b, the discrepancy between the pressures of the two zoomed-in regions is less significant as shown in Fig. 12a and b. This results in similar partial-slip states experienced by the two regions in surface 5, where stick ratios in regions 1 and 2 are 0.1784 and 0.1868 as labelled in Fig. 12c and d, respectively. Regarding surfaces 7 and 10 that show the same RMS roughness but different RMS slopes, for similar reasons, the same contact phenomena are observed for the selected regions as demonstrated in Figs. 13 and 14, respectively.
The root cause of these results is the introduction of the coupling between shear tractions and pressure, which changes the contact condition of gross sliding to make the surface easier to reach the sliding state with a lower tangential load. When comparing the contact solutions of two surfaces with the same RMS slope or RMS roughness under the current medium loading condition ( F x = 0.5 f W ), a more unevenly distributed pressure profile is often obtained for the surface with the higher RMS roughness or lower RMS slope. This leads to early slip for the contacting island where lower pressure is experienced and consequently brings about a lower stick ratio. It is noted that gross sliding should be more difficult to achieve for such a surface, in that higher loads are required to be in the slip state for the areas exhibiting higher pressure. Hence, a different response of the stick ratio to the surface roughness parameters shall be expected under a higher loading condition.

Contact of Rough Surfaces under Increasing Loads
A linear relationship between the stick ratio and ratio of input tangential load to static friction was reported in our former uncoupled rough surface study [27]. To verify if such a relation still holds under coupled conditions, the evolution of stick and slip regions in rough surface contact with increasing tangential load was simulated. The rough surface 1 presented in Fig. 8a and a higher normal load of 0.0025N were used to facilitate the observation of the dynamic transition from partial slip to gross slip within the contacting area. The separation of the stick (dark) and slip (grey) zones under no tangential load is given in Fig. 15a and an area was chosen and zoomed in to analyse. The stick zone within this area is found to keep decreasing and shifting to the rear part of the contacting region with increasing tangential load as shown in Fig. 15b-g, where such a change is relatively slow and steady when the load is low ( F x ≤ 0.7 f W ) compared with the significant shrink of the stick zone at a higher load ( F x > 0.7 f W ). The evolution of stick ratio with the increasing tangential load is shown in Fig. 15h, where the approximated linear relationship reported in our former uncoupled study [27] is apparently lost due to the coupling effects. Compared with the solution to the coupled smooth sphere contact problem investigated in Sect. 3.1, the rough surface experiences a lower portion of contacting area in stick when the normal load is applied alone but it requires a higher tangential load to reach the gross sliding state as labelled in Fig. 15h. To analyse the dynamic transition further, the partialslip contacts of different rough surfaces (Rough surfaces 1, 3, 5, 9, 10, 11 given in Table 2) were simulated under the same normal load of 0.0003 N but increasing tangential load. The relevant contact solutions are given in Fig. 16. The input normal load is found to affect the separation of the stick and slip regions since the contact of rough surface 1 ( R q = 1.0 μm) provides dissimilar solutions under normal loads of 0.0025 N (blue line in Fig. 15h) and 0.0003 N (orange line in Fig. 16). This is different from what was reported in our former uncoupled rough surface study [27]. Our interpretation of this is that a higher normal load brings about higher normal pressures due to the load balance, which in turn leads to more significant coupling effects on the shear traction (e.g. the magnitude of the surface displacement u xz increases). The tangential loads required to achieve gross sliding in all the simulated rough surface contact problems were found to always be higher than that in of an equivalent smooth sphere contact problem ( F x = 0.9463 f W).
Regarding the trend difference for surfaces with different roughness parameters under the increasing tangential load, one with the higher RMS roughness tends to always have less proportion of contacting area in stick when the load is low or medium, as illustrated in Fig. 16. This is consistent with the results obtained from the previous rough surface study under constant loading conditions (Sect. 3.2.1).
However, when the input tangential load is high (approximately 0.76-0.9725 as shown in zoomed view 1 in Fig. 16), the surface with the highest RMS roughness ( R q = 1.0 μm μm) exhibits an intermediate value of stick ratio. The trend then becomes similar to that under small and medium loading conditions once the load exceeds this range. The solutions presented regarding RMS slope are irregular as well when it comes to a high loading condition. As shown in Fig. 16, the surface with a higher RMS slope always experiences more stick ratio in the early period of the dynamic loading process. However, the surface with the lowest RMS slope ( g = 0.6455 ) turns out to require the highest tangential load to induce gross sliding. On the other hand, the surface with the highest RMS slope ( g = 1.7250 ) only requires an intermediate value of tangential load for global slip as labelled in the zoomed view 2 in Fig. 16. Such inconsistent trends can be related to the more non-uniformly distributed pressure experienced by some rough surfaces On one hand, when the tangential load is low or medium, it is easier for the area where pressure is relatively low to be in the slip state and this leads to a smaller stick ratio. On the other hand, it is more difficult for the entire surface to reach gross sliding since the area exhibiting higher pressure now requires a higher load to slip.
Although the trend between the surface roughness parameters (RMS roughness and RMS slope) and the stick ratio turns out to depend on the loading condition, these results tend to prove that the effects of the coupling between normal pressure and shear traction on the separation of the stick and slip regions depend on the roughness parameters of surfaces.

Concluding Remarks
By adjusting the shear tractions within the contacting region globally and adopting an alternative convergence criterion regarding the load balance in the tangential direction, the partial-slip contact of two elastically dissimilar materials was successfully reproduced by the newly developed model. The identification of stick and slip regions within the contacting area could be accomplished in a short time. (Usually, 120 s for smooth sphere contact problems with 256 × 256 nodes and 300 s for rough surface contact problems with 512 × 512 nodes using a desktop PC with four cores).
The numerical results of rough stick-slip contacts under coupled partial-slip conditions lead to the following conclusions: • Regardless of the roughness parameters, the tangential load needed for a nominally flat surface to reach gross sliding is higher than that needed for a smooth sphere. • With an identical RMS slope, the contacting region is insensitive to the change of RMS roughness under constant loads (low or medium). However, a higher RMS roughness leads to the decrease of stick region. • The contact area or the stick area is found to decay exponentially with the RMS slope under constant loads (low or medium); equations are proposed to determine the proportion of stick and contacting area to the whole simulation domain. • The linear response of the stick ratio to the input tangential load that occurs in the rough surface contact under uncoupled partial-slip conditions is now lost due to the additional interaction between shear tractions and normal pressures. Such coupling effects on the separation of stick and slip regions are found to vary with the surface roughness parameters. During the transition from partial slip to gross sliding, there are multiple observations. Under low and medium tangential loads, the surface with a higher RMS gradient or a lower RMS roughness exhibits more ratio of contacting area in stick as mentioned above. However, the trend becomes irregular when the load keeps increasing to the final loading point needed to achieve the gross sliding. .