Numerical Investigation of Presliding in Viscoplastic Spherical Contacts

Presliding interfaces such as in underplatform dampers of aircraft engines [1] and bolted joints in assembled structures [2] provide nonlinear structural compliance, damping, and modal interactions [3]. Nonlinearities are often related to the coupled nature of frictional and material properties, surface geometry, and external loads [4,5]. Cattaneo and Mindlin (CM model) independently showed the coupling of frictional slip with shearing of preloaded elastic spheres and thus the first slip-induced nonlinearity in presliding contacts [6,7]. Mindlin later demonstrated that frictional-slip-induced energy dissipation under cyclic loading scaled nonlinearly with maximum tangential force in a loading cycle [8]. Those classical findings serve as limiting cases where material nonlinearities and dissimilarities, surface irregularities (e.g., roughness), and failure can be neglected in a presliding contact scenario. Influence of several of those neglected factors was revisited by experiments and modeling [9–18].

Load-dependent nonlinearities are very common in structural dynamics. For example, bolted joints introduce damping that varies with the bolt preload and amplitude of imposed vibrations [19–21]. As an example to coupled load and material-dependent nature of contacts, recent experiments by Wei et al. [22] are illustrative. In that study, the authors conducted friction and wear tests on ultra-fine-grained aluminum and observed grain refinement or growth depending on the stress levels imposed on the contact. Changes in grain sizes altered the hardness of the contact, which in turn influenced stresses and further changed the grain sizes. These wear (plasticity) induced changes continued until an equilibrium is reached. Yet, there are materials where a similar equilibrium could not be reached quickly under contact loading. For instance, Eriten et al. [23] conducted presliding tests on flat high density polyethylene (HDPE) and polyurethane (PU) elastomer samples by an in-situ scanning electron microscope (SEM) picoindenter (PI88 PicoIndenter, Bruker Inc.) and observed significant deviations from the classical findings of CM model in terms of load-dependent nonlinearities in tangential stiffness and damping capacity. In particular, PU samples remained mostly fully adhered to the spherical rigid probe and so tangential contact stiffness remained fairly constant at values comparable to CM model predictions. Consistently, energy loss per loading cycle scaled with square of the maximum tangential force, which clearly indicated that material damping was the dominant dissipation mechanism and contributions of frictional slip and other nonlinearities were negligibly small. In contrast, HDPE samples delivered strong nonlinearities in presliding response. Large amount of plastic deformations and pile-ups were observed for HDPE under heavy normal preloads. Those geometric alterations of the contact lead to nonlinear trends in tangential contact stiffness and energy dissipation with varying loads.

Tribological responses of HDPE and similar polymers under loading conditions other than presliding were reported in numerous past works. For instance, in Ref. [24] the researchers found that inserting multiwalled carbon nanotubes into HDPE matrix enhanced hardness, thermal stability, and increased wear performance. The friction and the onset of wear in HDPE were studied by Looijmans et al. [25]. They found that due to lower yield strength, the frictional resistance of HDPE was lower compared to that isotactic polypropylene and due to the lack of strain softening, and that HDPE prevented strain localization and eventually delayed wear during scratching. Briscoe et al. [26] showed that ductile ploughing was a dominant mode of deformation mechanism for polyethylene (PE). Several examples of tribological response of similar polymers can be found in the literature [27–30]. A critical aspect of HDPE’s tribological response is its nonlinear viscoplastic response [31].The reflection of viscoplasticity was clear in stress-dependent creep response of the HDPE samples under especially heavy loading conditions [23]. In limited number of cycles (30) and durations (200 s), the nonlinear rate-dependent HDPE response resulted in continued accumulation of wear and energy losses even under presliding contact conditions. Such formation of wear in HDPE-like polymers was reported by others as well [25,32,33]. Longer durations under a constant load or deformation cannot eliminate this effect as shear stresses during presliding perturbs the stress state and triggers transient viscoplastic response.

The aim of the current paper is to model coupled friction and material nonlinearities observed in the aforementioned in-situ SEM experiments on HDPE samples [23] and relate them to load-dependent nonlinearities in tangential stiffness and energy dissipation. In particular, a 3D finite element (FE) model of a rigid spherical indenter pressed and sheared on a half-space is developed. A two-layer viscoplasticity model is adopted for the constitutive response of the half-space. The indenter geometry and boundary conditions are set in accordance to the experiments. Then, the parameters of the constitutive model are calibrated against normal contact responses observed experimentally. With the tuned material properties, the model is shown to predict the tangential (shear) response in close agreement with the experiments. The model further reveals that nonlinear viscoplasticity dominates tangential stiffness and dissipation at high normal preloads. Accumulated plastic deformations are also found to correlate well with the experimental wear profiles. Our results confirm further that nonlinear material response contributes significantly to the load-dependent nonlinearities in presliding contacts. The paper outline is as follows: Sec. 2 sets up the problem with the details of the finite element model, material properties, and boundary conditions. Section 3 features a comparison of simulation results with experimental findings of normal and tangential contacts. Section 4 discusses the limitations of the study, nonlinearities observed, and their influence on structural dynamics. Finally, the conclusions are listed in Sec. 5.

A 3D FE model was developed in which a rigid spherical indenter 17 μm in radius was set in contact with a half-space of dimensions 100 × 80 × 80 μm³ along xyz axes. The half-space is representative of the flat HDPE samples with sufficiently long dimensions to ensure no boundary effects. Once the components of the model (indenter and HDPE) were positioned properly in assembly, boundary conditions and contact interaction properties were assigned as shown in Fig. 1. The bottom surface of HDPE was fixed in all directions and the surface outside the contact was set as free. The front face of the substrate with a normal in +z-axis was assigned a translational symmetry, i.e., U3 = UR1 = UR2 = 0, where U3, UR1, UR2 are the displacement in z and rotations along x and y axes, respectively. A constant preload, P in the normal direction was applied to the reference point on the indenter. Subsequently, a cyclic tangential displacement δ was applied to the same reference point in congruence to the experiments reported in Ref. [23]. The simulated outputs to be compared with the experiments are normal penetration w and tangential force Q. See Fig. 3 for examples of the imposed loading and simulated outputs.

Finite sliding, surface-to-surface interaction was used to define the contact at the interface of indenter and HDPE. Penalty-based normal and tangential behavior was assigned as the interaction properties. The friction coefficient, μ was set as 0.64, which corresponded to the gross sliding limit observed from the experiments at the lowest preload, 1 mN. The indenter in this model was considered as the master surface whereas the part of HDPE in contact with indenter was assigned as the slave. Spatial discretization employed biased seeding to create finer meshing in the vicinity of contact as seen in Fig. 1. The half-space (HDPE) was discretized into linear hexahedral (C3D8R) elements and the rigid indenter was discretized into a combination of linear quadrilateral (S4R) and linear triangular (S3) elements with almost 50 nodes at the contact diameter for the lowest normal preload case. The discretization was further verified against predictions of the Hertzian contact model for normal loading and Cattaneo–Mindlin model for uniaxial tangential loading with linear elastic material properties [34]. Contact pressure, normal load–displacement, and stress results obtained at the plane of symmetry were all within 3% of the Hertz solutions. For a very small strain rate (a/R = 0.05), the tangential load–displacement and associated stiffness obtained were within 5% of the Cattaneo–Mindlin solution.

In summary, there are seven parameters to be determined in the two-layer viscoplastic constitutive law. Those parameters, their units, and actual values used in this work are listed in Table 1. The total Young’s modulus, E, the Poisson’s ratio, ν, and the yield strength of HDPE were set as 1400 MPa, 0.45, and 23 MPa, respectively, in line with the reported values in the literature [23,41,42]. Using Eq. (10), the parameter f = 0.42. In that equation, ω₀ was taken from the normal displacements observed at the end of normal loading stage of 1 mN normal load cases and ω_∞ was taken at t = 200 s, i.e., after significant creep has happened (see Fig. 3(b)). For a fixed Young’s modulus, E, and parameter f, A = 8 e − 5 MPa⁻ⁿ s⁻¹ was determined by relating it to the time constant (τ = 1/EAf) observed from the experiments and by setting stress and time-hardening powers to n = 1 and m = 0 in Eq. (4), i.e., linear viscoelasticity. However, the long-term creep observed in measured normal displacements exhibited stress and time-dependent nonlinearities, which then required curve-fitting of those Norton–Hoff rate law power parameters. Therefore, a FE-based optimization was used to determine the coefficients n and m by setting their search boundaries as [0, 5], and [−1, 0], respectively while keeping other constitutive parameters constant as listed in Table 1. The optimization routine minimized the normalized root mean square error between the measured and FE-simulated normal displacements for all four preloads and delivered the best fit for n and m as listed in Table 1. Note that this fit is not unique especially since in this material model, total deformation consists of elastic, viscous, and plastic contributions (Eq. (8)). Hence, changes in elastoplastic parameters would have immediate influence on the viscous parameters of the Norton–Hoff rate law. Since the literature reports elastoplastic properties of HDPE more readily than viscous properties, we chose to curve-fit the latter. Besides, satisfactory correlations of FEM predictions with the measured contact stiffness and plastic deformations observed in SEM images support our choice of curve-fitting the viscous properties while fixing elastoplastic parameters of the constitutive law.

Finally, a nonlinear quasistatic analysis was conducted with the above determined viscoplastic parameters. The entire simulation in this study was conducted on commercially available software Abaqus 2018. A typical simulation takes around 5 h on a PC with a 3.60 GHz Intel(R) Core(TM) i7-7770 processor and 48 GB RAM.

Figure 4 shows the time evolution of the normal displacement predictions of the FEM at different preloads in comparison with the experimental responses reported in Ref. [23]. The design of experiment simulation is intended toward determining the load-dependent nonlinearity and plastic deformation effects in stiffness and dissipation of HDPE. Due to the limitations of experiments and to focus more on the other aspects of contact physics and mechanics, the study is designed in this way. At all preloads, normal displacements exhibited a sharp increase initially (mostly due to elasticity), which corresponded to the ramping of the normal force in Fig. 3(a). A nonlinear creep in normal displacements followed that sharp increase. Nonlinear creep period coincided with the constant normal force (hold) stage of loading (Fig. 3(b)). The rates of nonlinear creep observed in experiments were found to depend on both mean contact pressures and time. Hence, the time-hardening power law with nonzero stress and time-hardening powers, Eq. (4) was employed in modeling the viscoelasticity of HDPE samples. The FEM with the constitutive parameters listed in Table 1 delivered reasonable predictions for nonlinear creep at all preloads. At the lowest and highest preloads, the FE predictions overshot and undershot the experimental observations, respectively. The correlations at intermediate preloads (5 and 10 mN) were much better, and for all normal displacement evolutions, correlations with experiments were reasonable (R² ≥ 0.9). In summary, the FEM results simulated the elasticity-dominated sharp increase and nonlinear creep in normal displacements observed experimentally, verifying further the two-layer viscoplastic material model and parameters (E, ν, σ_y, A, n, m, f) listed in Table 1. Note that both FE simulations and experiments exhibited small oscillations during the hold stage of the normal loading. Those oscillations had 6 s period as the cyclic tangential displacements. Neither the FEM nor experiments could fully decouple normal and tangential loading responses. Elastic mismatch and imposed boundary conditions on the rigid probe were primary sources of coupling, which lead to small perturbations on the nonlinear creep response [13,43].

The contact radii predicted by the FEM at normal preloads 1, 5, 10, and 15 mN are compared to the experimental [23] and analytical (truncated) contact radii in Fig. 5. The truncated contact radii were provided to estimate the extent of plastic deformations and were calculated as $at=R2−(R−ω0)2$⁠, where R is the radius of the indenter, and ω₀ is the normal displacement at the end of ramping of the normal force (at t = 11 s). As expected, the contact radii increased with the normal preload. Both FE predictions and experimental values were in close agreement with the truncated contact radii. That close agreement was due to the large amount of plastic deformations that occurred in the vicinity of the contact [23]. Since the FEM predicted contact radii congruent with those observations, the yield strength σ_y extracted from the literature (Table 1) could be accepted with confidence. Note that mean contact pressures estimated from the experiments ranged from 24 MPa at 1 mN to 53 MPa at 15 mN. Therefore, elastoplastic response was expected at all preloads. The FEM predicted the corresponding mean contact pressures as 21 MPa and 61 MPa, within 20% of the experimental values. Further validation of the FEM and the selected constitutive parameters will be provided by the tangential response of the contact discussed in Sec. 3.2.

Representative hysteresis loops formed by plotting tangential force as a function of tangential displacement are presented in Fig. 6. In particular, the hysteresis loops obtained from the experiments and FE simulations are shown at two extreme preloads (1 and 15 mN). Five loops in each plot correspond to the output of five different maximum tangential displacement amplitudes. The fifth cycle at each maximum tangential displacement amplitude case is plotted as the representative steady-state response.

The loops at 1 mN preloads predominantly exhibited gross sliding, i.e., constant tangential force for most of the loading cycle [6,7]. At high normal preloads, presliding response dominated with seemingly linear tangential force and displacement relation. Besides, energy dissipation at high normal preloads decreased as can be inferred from narrower loops. Similar trends in energy dissipation were also observed at fixed normal force and decreasing tangential displacements. Lower loads and smaller tangential displacements result in the increased stiffness (slope) due to the increase of contact radius as the time approaches from high amplitude to low amplitude states (see Fig. 3). The asymmetric hysteresis loops obtained were previously attributed to uneven plastic deformations and pile-ups over the contact in Ref. [23]. These types of asymmetric hysteresis loops are also observed in the earlier works in which elastoplastic shakedown happens at the very early stages of fretting [17,44,45]. In general, Fig. 6 indicates that the FE-simulated and experimental hysteresis loops agree well at the lowest and highest preload cases. Further comparison of load-dependent tangential stiffness and energy dissipation extracted from similar loops at all preload cases studied follows next.

According to CM model [8], the linearized contact stiffness for fully adhered elastic spherical contacts is 8G*a, where G* is the equivalent shear modulus, a is the Hertzian contact radius. The tangential stiffness values obtained from Eq. (11) are normalized to that asymptotic stiffness and plotted with respect to maximum tangential forces normalized to the preloads used. At the lowest tangential forcing, normalized tangential stiffness values decrease from about 0.7 at 1 mN preload (Fig. 7(a)) to about 0.4 at 15 mN preload (Fig. 7(d)) gradually. This is because of the evolution of contact from elastic–plastic to fully plastic deformation states by an increase in preloads.

At 1 and 5 mN preloads, the tangential stiffness decreases with tangential forcing thanks to frictional slip. At higher preloads, frictional slip is suppressed. Instead, softening effect of viscoplasticity competes with increasing contact areas due to creep, junction growth, pile-ups and debris formation. Thus, the tangential stiffness at 10 and 15 mN preloads remains relatively constant with increasing tangential forcing. Mean contact pressures for those cases range from 50 to 61 MPa, both of which are significantly higher than the elastic shakedown limit of σ_y/2μ ≈ 16 MPa [46,47]. Therefore, plastic shearing over the contact surface sustains over the whole duration of the presliding tests under high preloads. This is in line with earlier observations on plastic shearing of elastoplastic contacts under heavy preloads [11,48]. At all preload cases, normalized tangential stiffness predicted by the FE simulations shows reasonable agreement with the experimental values. Hence, we will utilize FE-predicted dissipation breakdown among frictional slip, viscoelastic and plastic deformations to describe qualitatively the dominant mechanisms governing presliding response of HDPE. Besides, we will compare von Mises stresses from the FE model to SEM images of the contact area after 15 mN experiments. This comparison will further validate the FE model’s capability of accounting for plastic deformations that were observed experimentally.

The SEM images of the contact taken after the experiments exhibit accumulated damage in the form of pile-ups and radial and circumferential cracks (Fig. 9(a)). The stresses obtained via the FEM support the observed damage accumulation. In particular, Figs. 9(b)–9(d) show FE-simulated von Mises (SVon), maximum normal (S11), and hydrostatic (SHyd) stresses in the normal direction to the contact that is under 15 mN preload and at maximum tangential displacement (see Fig. 3). Extreme values of stresses are obtained in the range of −8.3 to 9.2 times the yield strength of the HDPE (SY ≡ σ_y given in Table 1), which clearly suggests that damage is inevitable. For instance, von Mises stress reaches 9.2σ_y close to the trailing edge of the contact (Fig. 9(b)). In the two-layer viscoplastic representation (Fig. 2), such high stresses would lead to significant plastic flow in the material. Plastic flow over cycles of tangential un/loading accumulates at the contact edge and forms the pile-ups as seen in the postmortem SEM image (Fig. 9(b)). Similarly, maximum tensile stresses reach about 4.4σ_y, and thus they are higher or comparable to the ultimate tensile strengths reported for HDPE with different densities [41,50]. Those stresses can explain the ring cracks shown in the SEM image. Besides, circumferential residual stresses during cyclic un/loading and upon retraction of the probe apply tension about several times the yield strength, which would result in radial cracks observed in Fig. 9(a). Similarly high tensile hydrostatic stresses also suggest the potential void-growth over the contact. Note that hydrostatic stresses were recently used as an indicator of wear in polymer tribology by others too [25]. In summary, accumulated damage observed in the SEM can be explained by the simulated stress field localized to the contact. This further validates the constitutive fitting. However, a more complete simulation would integrate damage initiation and propagation models and simulate actual crack nucleation and growth. Such a model is expected to deliver slightly different stresses and strains. Overall, this is not expected to change the relative importance of dissipation mechanisms on observed stiffness and damping nonlinearities as demonstrated by one of the authors by extended FEM of presliding spherical contacts [12].

Load-dependent nonlinearities observed in tangential contact stiffness (Fig. 7) and damping capacities (Fig. 8) can be qualitatively explained by the breakdown of dissipation into different mechanisms. Note that the classical CM model employs linear elastic material response and thus attributes all the softening and dissipation nonlinearities to frictional slip. In the constitutive and contact models we adopted and validated here against presliding tests on HDPE samples, frictional and viscoplastic material losses are expected to be the primary sources of dissipation [5,23]. Integration of dissipation rates over the volume of the simulated deformable body for a period of tangential loading (T = 6 s) delivers dissipation breakdown as $∫0T∫Vσi:ε˙idVdt$⁠, where i = Fr, Vi, Pl refer to the frictional, viscous, and plastic stresses and strain/slip rates, respectively. Figure 10 shows the energy dissipation breakdown in percentages as measured from the FE simulations at the fifth cycle of the lowest and highest tangential loads and two extreme preloads (1 and 15 mN). The breakdown for the intermediate cases follows a gradual transitioning between the extreme cases and thus is omitted for brevity. It is noticeable that friction losses are found to be the dominant contributors of energy dissipation only at low preloads (1 and 5 mN) and high tangential forcing (Fig. 10(b)). Those are the loading scenarios corresponding to significant softening observed in Figs. 7(a) and 7(b), and a sharp increase in damping capacity (Figs. 8(a) and 8(b)). This is in line with predominantly elastic frictional contacts under gross sliding. Material losses dominate all other loading scenarios. For instance, regardless of the normal preloads, viscous losses are the largest damping source at small tangential forcing (Figs. 10(a) and 10(c)). Higher tangential forcing at high preloads induces more plastic shearing and thus the viscous and plastic contributions in damping become comparable (Fig. 10(d)). This increase in the relative contribution of plasticity is also pronounced in the nonlinear trends observed in the damping capacity, e.g., Fig. 8(c). More plasticity leads to pile-ups and asymmetric alterations in geometry, which then results in asymmetric stiffness increase and thus reduction in total damping capacity. Nevertheless, domination of material losses at high preloads (10 and 15 mN) explains the relatively lower band of damping capacities shown in Figs. 8(c) and 8(d) (Ψ ranging from 0.1 to 1). In a similar fashion to the current study, one of the authors, Eriten and his collaborators had studied relative contributions of frictional slip, crack formation, and plasticity on the onset of sliding in elastoplastic spherical contacts [12]. They had also observed the increasing influence of plasticity on frictional sliding thresholds at high normal preloads. The formation of pile-ups and junction growth due to excessive plasticity was also known to alter contact geometry [9,51–55]. So, the current findings on HDPE samples confirm those earlier observations while adding the influence of nonlinear viscoplastic deformations on key contact parameters: tangential stiffness and damping capacities.

The aforementioned results explain governing mechanisms behind nonlinear load-dependent trends observed in presliding response of HDPE samples. Selection of HDPE was intentional due to two particular reasons as mentioned in our previous work with more details [23]: (i) machine compliance of in-situ SEM picoindenter prevented gathering presliding displacement data for the stiffer materials such as ceramics or metals and (ii) stress-dependent elasto-viscoplastic response of HDPE also reported elsewhere [56] was expected to deliver material losses and related nonlinearities in contact stiffness and damping in presliding. Being approximately two orders of magnitude more compliant than typical engineering ceramics and metals, HDPE helped us to resolve the machine compliance issue. However, HDPE being an insulator posed the imaging challenge via electron microscopy. To circumvent that challenge, we coated HDPE samples with 5-nm-thick silver and used low beam energies (5 kV) for SEM imaging [23]. Previous studies on fretting of layered elastic materials show the diminishing influence of coating layers for small layer thickness to contact radius ratios [57,58]. Note that the contact radii in the experiments and simulations discussed here range from about 4 to 9 μm (Fig. 5), and thus that ratio is smaller than 0.0013. Therefore, we omit the silver layer in the numerical modeling of the half-space.

Previous studies reporting stress-dependent elasto-viscoplasticity of glassy polymers and HDPE had employed uniaxial tension/compression, creep, and relaxation tests [56,59]. Here, we observe and model the load-dependent nonlinearities using a generic two-layer elasto-viscoplastic constitutive law (Fig. 2). That law needs seven parameters to describe elastic, viscous, and plastic responses. Table 1 lists one set of those parameters curve-fit to normal displacement responses obtained in the experiments. As mentioned in Sec. 2, elastoplastic properties E, ν, σ_y found correlate well with the values reported in the literature. Besides, close agreements between simulated values and experimental tangential contact stiffness (Fig. 7) and plastic deformations (Fig. 9) validate further the values used for the elastoplastic properties. Viscous parameters A, n, m, f of the Norton–Hoff rate law, however, are not verified independently in this work partially because the experiments were conducted at a single rate (cyclic loading period T = 6 s). Nevertheless, long creep history extracted from the experimental normal displacements at different preloads enabled us to identify the need for a stress and time-dependent creep law for the viscous part of the two-layer viscoplastic model. Besides, we estimated time constants and un/relaxed moduli from the experimental normal displacements at small preloads and resort to an approximate linear viscoelastic analyses to set two of those viscous parameters, A and f. That asymptotic analyses reduced the number of curve-fit parameters to two: stress and time-hardening powers, n and m, respectively. Note that curve-fit n = 1.05, and thus the creep rates were fairly close to linear viscoelasticity in terms of stress. However, significant reduction of creep rate occurred as a function of time in each experiment, explaining the curve-fit m = −0.85. Note that viscous parameters chosen deliver tangential stiffness and damping that are heavily influenced by the viscous material losses. An optimization and curve-fitting against both normal and tangential experimental responses could have delivered another set of viscous parameters, and thus the energy breakdown values reported in Fig. 10 could change. That change, however, is expected to be limited thanks to consistent and validated elastoplastic formulation chosen. Nevertheless, the main conclusion of the work remains unaltered; i.e., load-dependent nonlinearities in tangential stiffness and damping capacity are largely governed by constitutive effects rather than friction in viscoplastic spherical contacts.

Material-response dominated contact behavior is common to soft materials. For instance, one of the authors had reported the influence of poroviscoelastic response in rate-dependent work of adhesion and presliding response in articular cartilage [60,61]. Both work of adhesion and static friction thresholds were shown to increase orders of magnitude with un/loading rates, and thus softening and dissipation responses could cover fully adhered to fully lubricated contact scenarios for the same contacting material pair. Besides, changes in contact geometry such as pile-up during loading add more complication to the contact response as documented in Fig. 9, and elsewhere [55] especially under heavy preloads. Such constitutive and geometric nonlinearities and their temporal evolution continue to hinder predictive contact modeling and health monitoring under presliding conditions. Focusing on friction alone as adopted by many researchers including the author and his collaborators [16,62,63,] will not resolve those issues. More advanced models coupling interfacial and material nonlinearities with failure mechanics are promising but come with an expense of computational inefficiency [13,64] and thus cannot be adopted in structural dynamics simulations. Note that the load-dependent trends in stiffness (Fig. 7) and damping capacities (Fig. 8) exhibit clear deviations from the classical friction-dominated nonlinearities. This means that the dominant frequencies and damping of a dynamic system with presliding contacts would change with the vibration amplitudes and that change would vary depending on the dominant mechanism of nonlinearity. Therefore, combination of data-driven nonlinear model updating [65] and system identification tools [3] can provide real-time monitoring and insights into what mechanisms govern interface response. The current work of the authors focuses on identification and distinguishing the sources of such nonlinearities in composite soft-hard material systems.

The aim of the current paper is to model coupled friction and material nonlinearities in HDPE samples and to relate them to load-dependent nonlinearities in tangential stiffness and damping capacities as previously observed in in-situ SEM experiments. In particular, a 3D finite element model of a rigid spherical indenter pressed and sheared on a half-space is developed. A two-layer viscoplasticity model is adopted for the constitutive response of the HDPE samples. The indenter geometry and boundary conditions are set in accordance with the experiments. Then, the parameters of the constitutive model are calibrated against normal contact responses observed experimentally. With the tuned material properties, the model is shown to predict the tangential (shear) response in close agreement with the experiments. The model further reveals that nonlinear viscoplasticity dominates tangential stiffness and damping capacities at high normal preloads. Accumulated plastic deformations also correlate well with the wear profiles obtained by the SEM images. In summary, our results show that nonlinear material response contributes significantly to the load-dependent nonlinearities in viscoplastic presliding contacts, and those nonlinearities deviate considerably from the predictions of the classical friction-dominated contact models. Those findings further assert that accurate constitutive and failure laws are required for a contact model to predict apparent nonlinearities in stiffness and damping in systems involving viscoplastic materials.

This work is partially supported by the National Science Foundation CMMI-CAREER (Grant No. 1554146).

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.