nents of composites, visual and tactile perception, thermal conductivity and even electrical conductivity in the case of some new smart materials and sensors. However, the quantification and modelling of the contact and tribological properties of soft polymers remains a significant challenge, necessitating a trial-and-error approach and a substantial number of measurements for the optimisation of product properties. Since an analytical model that would take into account all the special properties of soft polymers would be so complicated as to prevent its application to practical problems, it is desirable to find an applicable approximate model that adequately addresses the main issues of contacts with soft polymers. As a first step towards such an approximate model for soft polymer contacts, it is logi- Science and Research 30 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 Introduction The real contact area between two materials is typically orders of magnitude smaller than the apparent contact area, due to the ubiquitous roughness of all surfaces. However, it is the real contact area that determines the tribological phenomena occurring at the interface [1]. The present study considers the phenomenon of large deformation contacts, which are of particular practical relevance in the context of soft polymers. Here, the term ‘soft polymers’ shall refer to polymeric materials that are used above their glass transition temperature. In this so-called rubbery state, the Young’s modulus can be orders of magnitude lower than in the hard, glassy state below the glass transition temperature. Due to the small forces required for deformation, typical application forces can easily lead to large deformations. The present study does not focus on soft polymers such as silicones or rubber, for which the deformation behaviour can be adequately described by elastic or hyperelastic models and the contact properties have been extensively examined by other researchers [2, 3, 4]. The study instead focuses on soft polymers with non-ideal elastic deformation behaviour, including polyethylene (PE), plasticised polyvinyl chloride (p-PVC) and polyurethanes (PUR). Applied in foils, paints, varnishes and coatings on a wide range of substrates, these materials determine the surface properties of a vast variety of everyday products, including contact properties in contacts to other surfaces. These contact properties, in turn, influence a number of crucial product characteristics, including blocking and sliding behaviour, grip, the generation of disturbing noises (such as squeaking and creaking), durability and susceptibility to wear, adhesion between different compo- Insight into large deformation contacts of soft polymers with molecular dynamics simulations Susanne Fritz* submitted: 20.09.2024 accepted: 9.05.2025 (peer review) Presented at GfT-Conference 2024 The size of the real contact area between mating materials represents a crucial quantity in the context of all tribological problems. Consequently, it is a topic that has been the subject of intensive investigation. Nevertheless, the subject of contact formation with soft polymers remains largely uninvestigated. Soft polymers are frequently employed as foils or coatings, thereby influencing the surface characteristics of a vast variety of everyday products. Employed at temperatures above their glass transition point, soft polymers with low values of the Young’s modulus, permit significant deformation during contact formation. To assess the impact of such large deformations on the size of the real contact area, atomistic molecular dynamics (MD) simulations of nanoindentation tests have been conducted. Depending on the material properties and loading conditions, the formation of the contact area was observed across the entire range of indentation depths. Keywords contact, contact mechanics, polymers, nanoindentation, MD simulations, Hertz-model Abstract * Dr. rer. nat. Susanne Fritz Orcid-ID: https: / / orcid.org/ 0000-0002-5006-8826 Department Surfaces, FILK Freiberg Institute gGmbH Meißner Ring 1-5, 09599 Freiberg, Germany cal to begin with the Hertz model, given that the majority of models for the real contact area in contact mechanics are based on the Hertz model. Heinrich Hertz [5] derived an analytical solution to the simple contact problem between a rigid sphere and an elastic half-space. A numerous number of other authors proceeded to expand his theory to include the effects of adhesion [6], friction [7, 8], viscoelasticity [9, 10] or roughness [11, 12], for example. However, since Hertz built his equations upon a certain set of fundamental assumptions, all the other theories are also only valid under these assumptions, which are not always met by contacts with soft polymers. These assumptions are in particular: • homogeneous, isotropic solids that can be considered as a continuum, • linear elastic behaviour, and • deformations that are so small, that the shape of the surface in the contact region can be reasonably well approximated by a second-degree polynomial. In contrast, the behaviour of soft polymers is typically characterised by a pronounced viscoelastic-viscoplastic response and can even exhibit viscous flow to some extent. At the length scale of small surface asperities (micro-roughness), a polymeric material made of large chain molecules and maybe containing crystallites cannot simply be assumed to behave as a homogenous continuum and usually deforms anisotropically. But the main difference surely are the large deformations, that can easily occur with soft polymers, especially at small asperities. The question of high deformation Hertz contacts has already been addressed by other authors [13, 14, 15, 16, 17, 18, 19], with results ranging from the good applicability of the Hertz model to the introduction of correction terms and the development of new empirical equations. Nevertheless, as the actual contact area is not readily accessible through direct experimentation, all the referenced papers have employed analytical or finite element analysis (FEA). This involves defining a mathematical equation to describe the deformation behaviour of the material, along with the appropriate material properties, as the input for the model. This approach is effective when the deformation model accurately represents the material’s behaviour, but it can also introduce unpredictable errors when the deformation behaviour is not straightforward to describe mathematically, as is the case with soft polymers. For this reason, molecular dynamics (MD) simulations [20, 21] were used in the present study instead of FEA. In MD simulations, the material is described at the atomic level, and the behaviour of the material is calculated based on the known interactions between the atoms. As a consequence, only small time and length scales are accessible, but no models or material properties are needed to describe, for example, the deformation behaviour, and nanoscopic properties that are difficult to measure can be accessed. The use of MD simulations for the analysis of contact problems has already been demonstrated in other studies [22, 23, 24, 25, 26, 27]. The present study should be regarded as an empirical investigation utilising MD simulations as an alternative to experiments, assuming that the simulated material behaviour is identical to the real behaviour of soft polymers. The simulations were designed based on the ideal Hertz contact between a non-deformable spherical indenter and a surface of a soft polymer. Nanoindentation simulations were conducted under various conditions, with the calculated properties subsequently compared to the equations of the Hertz model. The objective was to evaluate the applicability of the Hertz equations for approximate predictions of the real contact area for soft polymers where the basic assumptions are not met. Methods MD simulations were conducted using the GROMACS software package [28, 29, 30], which was run on NVIDIA Quadro P2000 graphics cards. The leap-frog algorithm [31] was employed for the integration of the equation of motion with a time step of 3 fs (or 1 fs for equilibration runs). The application of three-dimensional periodic boundary conditions enabled the simulation of condensed matter at the nanoscale. Nonbonded interactions were cut off at a distance of 1 nm with the Verlet algorithm for neighbour searching [32]. Temperature and pressure control were realised by the velocity-rescaling thermostat [33], and the anisotropic Berendsen barostat [34], respectively. Interaction energies were calculated using the GROMOS 53a6 force field [35, 36, 37]. Polyethylene (PE) was selected as a representative of soft polymers for a number of reasons, which were beneficial for the investigations: Its glass transition point lies considerably below ambient conditions, it is chemically relatively simple, and, as a semi-crystalline polymer, it offers the possibility of considering a range of deformation properties without affecting the chemical properties of the material. Starting from linear, unbranched polyethylene chains comprising 2000 carbon atoms, bulk simulation boxes of dense PE with edge lengths of at least 50 nm were generated via a subsequent heating and cooling process, whereby the different polymer chains underwent coiling, mixing and entanglement. By varying the cooling rate between 100 and 0.04 K/ ns, the degree of crystallisation could be adjusted between 0 and 50 %. A higher degree of crystallisation of the PE is associated with a higher density, a higher Young’s modulus and a lower Poisson’s Science and Research 31 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 distances in the xy direction. These atoms were not subject to mutual interaction but were constrained to maintain their respective distances. With regard to the interaction with the PE surface, the indenter atoms were treated as PE atoms. For the nanoindentation setup the indenter was positioned above the PE surface, as illustrated in Figure 2b. The bottom atoms of the PE layer model were fixed in space to prevent the model as a whole from moving through the application of external forces. Distance-controlled indentation simulations were conducted using the slow-growth functionality of GROMACS [39] to move the indenter towards the surface by a minimal defined distance every time step, determining the indentation velocity υ, which was varied in 4 stages (0.5, 1, 10, and 20 nm/ ns). The indentation process was simulated up to the full penetration of the spherical tip (p = R), but to a maximum of p = 15 nm to ensure, that the deformation was not affected by the limited size of the PE surface layer. For each parameter variant (c, R, v), a minimum of three determinations were conducted with the tip position varying along the surface to ensure statistical reliability. The normal forces were calculated by summing all the Science and Research 32 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 ratio (Figure 1). A more detailed description of the box generation method can be found in [38]. For the nanoindentation simulations PE with degrees of crystallization c of 0, 12, 30, 40 and 50 % was used. In order to form a 2D surface model and still utilise the faster 3D periodic routines, a vacuum slice with a thickness of 150 nm was introduced by enlarging the simulation box in the z-direction without modifying the atomic positions. The vacuum slice is of sufficient size to prevent interactions between the periodic PE slices, ensuring that the z-direction periodicity does not affect the surface. The simulation boxes were equilibrated for 3 ns to allow recombination on the freshly cleaved surfaces. Due to the presence of stable crystallites within the surface plane, the obtained surface models exhibited a certain nanoscale roughness (S a ≈ 1 nm). For the nanoindentation simulations, non-deformable indenters with spherical tips were produced. The initial basis geometry was a cone with an opening angle of 25°, the apex of which was rounded with a curvature radius R (3, 8, 10, 12, 14, and 18 nm). The indenter was composed of atoms positioned according to a mathematical calculation and arranged in a regular pattern with fixed Figure 2: a) Illustration of the symbols used to characterize the Hertz contact, b) setup of the nanoindentation simulations, illustrated by a detail of a cross-section through the simulation box of a polyethylene surface with c = 0.5 and an indentation radius of R = 8 nm, and c) exemplary snapshot of tip atoms in contact with PE atoms projected onto the xy-plane, together with the according contact circle of area A and radius a Figure 1: Relation between the simulated values of crystallinity c, cooling time τ c , density p, Young’s modulus E and Poisson ration υ for the semi-crystalline PE model polymer forces acting on the PE atoms in the z-direction. In order to quantify the size of the real contact area, the number of tip atoms that met a specific distance criterion (0.65 nm) to the PE atoms was counted for each time step. As a result of the uniform distribution of the tip atoms, geometric considerations could be employed to calculate the penetration depth of the tip p, and both area A and radius a of the contact circle parallel to the xy plane from the determined contact number (Figure 2c). Results and discussion In order to ensure the greatest possible comparability between the Hertz model and the simulations, the present investigation focuses on the short-term deformation behaviour. Given the applied velocities, it can be reasonably assumed that slow processes, such as material deformation caused by viscoelasticity or adhesion, can be disregarded during the indentation part and are only of relevance for the retraction part of the nanoindentation, which shall not be discussed here. In the absence of viscoelasticity and adhesion, the dimensions of the contact area can be assumed to be independent of previous states, and only dependent on the current position of the indenter. This assumption was tested through a series of simulations, which were repeated identically, apart from the adhesive force between the indenter and the surface (realised by the variation of the interaction potential between the tip and the PE atoms by ±50 % relative to the original values). The results showed no effects on the formation of the real contact area or the acting forces during the indentation part of the simulation (while significant effects were observed during the retraction part, which will be discussed elsewhere). Thereby, In the following, the term ‘data set’ is used to refer to the corresponding values of p, d, A and a, resulting from an indentation process with specific values of c, R, and v after time t in determination x. The analysis of the present study is based on a total of 23.400 of such data sets. In comparison to experiments, simulations facilitate a straightforward quantification of the real contact area. Nevertheless, the indentation depth is less clearly defined on the nanoscale for a nanoscopically rough surface. In the present study, two distinct definitions were employed for the indentation depth. In the first definition, d c represents the distance in the z-direction between the indenter positions at the specified time step and the initial atomic contact time step. This definition allows for a unique determination of the indentation depth within the simulation, is reasonably comparable to the Hertz model and will be used in the first part of the presented analysis. Nevertheless, since the atomic contact is not easily accessible within experiments, experimental setups typically employ an increasing force to determine the initial position of the surface. This definition was also applied in the presented simulations, resulting in d f , which was then used in the second part of the comparison with the Hertz model. Figure 3 shows a scatter plot of corresponding values of d c and d f for all data sets. Despite the inherent uncertainty in d f , stemming from the conjunction of high velocities and the statistical noise associated with the calculated normal force, d c and d f exhibit a constant shift of approximately 1 nm. This shift is observed to be independent of the indenter radius and the crystallinity of the PE. The shift between d c and d f can be attributed to adhesion forces. The attractive interactions between the atoms of the indenter and the PE result in a negative normal force at the initial atomic contact. It is only when repulsive forces due to material displacement exceed the attractive forces that the normal force increases, which is interpreted as reaching the surface in terms of d f . In comparison to harder materials, the necessary penetration depth to overcome the attractive forces is significantly greater for soft polymers, due to the lower modulus of deformation. This is the reason for the discrepancy between d c and d f , which should have minimal impact at the macro level but significant effects at the nanoscale. According to the Hertz model [7], the indentation depth d of a spherical indenter with radius R results in a displacement u z of the elastic surface within the contact area, that depends on the distance r in xy-direction to the con- Science and Research 33 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 Figure 3: Relationship between the differently defined values for the indentation depth for all datasets; the values are shifted by a constant of ≈ 1 nm independently of indenter radius R or crystallinity c (6) In order to facilitate a comparison between these theoretical considerations and the conducted simulations, the proportionality a 2 = Rd, as defined in Equation (2), was examined for all data sets, with d c representing the indentation depth. At the nanoscale, the material does not exhibit the characteristics of a homogeneous, isotropic continuum. Instead, it deforms anisotropically due to the semi-crystalline structure, which gives rise to considerable fluctuations in the calculated quantities. Nevertheless, when considering the mean behaviour across all data sets, it can be seen from Figure 4a that a 2 = Rd seems to be a very good approximation for the whole data range. This means, the approximation is not exclusive to small deformations, as given by the Hertz model, but can also be reasonably well applied to large deformations extending up to the indenter radius, despite the violation of all crucial assumptions. This is an unexpected result, given that for large deformations, the linear relationship between a 2 and p is no longer valid and must be replaced by the non-linear relation described in Equation (6) (Figure 4b). Concurrently, however, an increasing deviation from the relationship p = d/ 2 (Equation (5)) is also observed with increasing indentation depth (Figure 4c). This appears to compensate for the non-linear relationship between a 2 and p, resulting in a seemingly linear relationship between a 2 and d. Thus, it can be seen that the violation of the small deformation condition results in significant discrepancies from the Equations (1) and (5), yet it nonetheless produces a seemingly identical relationship to Equation (2). As previously stated, this is not a consequence of material deformation caused by adhesion or viscoelasticity. Given that Equation (2) offers a satisfactory approximation of the simulated contacts, even in the context of significant deformations, Equation (4) was similarly considered by examining the proportionalities F ~ a 3 , F ~ 1⁄R, and F ~ E * . and . In conducting this analysis, only those simulations with a velocity of 20 nm/ ns were considered, as this velocity was also applied for the cal- 3 4 5 67 8 4 9 Science and Research 34 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 tact point at r = 0 (Figure 2a). Equation (1) is derived from the premise that a spherical surface can be approximated with a high degree of accuracy by a seconddegree polynomial in the vicinity of the contact point. This implies that the approximation is only applicable to small deformations, where the displacement around the contact point is the primary concern. (1) By identifying a pressure distribution that causes the displacement described in Equation (1) and integrating it over the contact area, the Equations (2) and (3) for the contact radius a and the normal force F, respectively, can be derived [7]. (2) (3) Equation (2) and (3) can be combined to give Equation (4), which relates the normal force with the contact radius. (4) By employing Equation (1) to ascertain the displacement u z at the boundary of the contact area r = a with the aid of Equation (2), one arrives at Equation (5), which asserts that the penetration depth p is precisely half of the indentation depth d. (5) Obviously, Equations (1) to (5) are only valid for values of d that are small compared to R. Considering the simple case of full penetration (p = R), it can be seen from Equation (2) that the approximation already causes an error of 50 % for the real contact area (a 2 = 2R 2 compared to the exact geometrical solution of a 2 = R 2 ). For geometric reasons, the precise relation between the contact radius a and the penetration depth p for sphere segment is given by Equation (6). : $ 5 8 ; 4 67 9 3 4 5 7 9 < 5 > ? @ A 8 B 4 C 7 D 9 < 5 > ? @ EA 8 B 4 F 3 D 7 5 > ? @ G 7 3 D 9 5 6 9 Figure 4: Observed relations between the contact radius a, the indentation depth compared to first atomic contact d c and the penetration depth p for all simulated data sets culation of the bulk deformation properties in Figure 1c, which are markedly velocity-dependent. As the example simulation results in Figure 5a show, there is a perfect correlation between F and a 3 not only for small deformations as proposed by Equation (4), but over the entire range of indentation depths up to the indenter radius. According to Equation (4), the slope of the regression line m should be proportional to the inverse of the indenter radius 1/ R. This hypothesis was tested and confirmed in Figure 5b for simulations with the same degree of crystallisation but different indenter radii. According to Equation ( 4 ), the slope of the regression line should now be solely dependent on the deformation characteristics of the surface E * = E ⁄ (1 - υ 2 ). E * was calculated in this way for all the considered PE surfaces with different degree of crystallisation. Despite the excellent reflection of the relationships expressed in Equation (4) in the simulation results, Figure 5c demonstrates a significant discrepancy between the calculated E *values from indentation simulations and the E * -values determined from bulk simulations (see Figure 1c). While a similar tendency is evident, a systematic deviation is also present. The discrepancy between the deformation properties determined from indentation and bulk simulations was also confirmed qualitatively by experiments, which will not be illustrated in detail here. Creep experiments were conducted on a variety of different samples of PE, p-PVC, PUR, and silicones, using spherical indenters with radii between 0.8 and 20 mm and constant loads ranging from 20 to 80 mN. In order to obtain the Young’s modulus, the measured creep curves were fitted with the equation provided by [40], which accounts for the viscoelasticity of the material. Despite the fact that the fits resampled the measured curves with great precision and the small deformation condition was fulfilled, the obtained Young’s moduli exhibited a notable and systematic discrepancy from the values obtained through standardised tensile tests. This finding aligns with the simulation results. The findings presented thus far have utilised the indentation depth defined by the initial atomic contact, d c . For d f values, which are more appropriate for experimental setups but exhibit a constant shift in comparison to d c (see Figure 3), there is a notable discrepancy between the analytical equations and the simulation results. Ne- Science and Research 35 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 Figure 5: Evaluation of the relationship in equation (4) for the conducted simulations by testing the proportionality between a) F and a 3 , b) F and 1/ R, and c) the resulting E * from indentation simulations and the calculated E * from bulk simulations Figure 6: Utilizing the definition of the indentation depth based on the force, d f , to get approximate values of the real contact area: a) approximation of d c via equation (7), b) and c) approximation of p and A via equation (8) Acknowledgement The research project “Determination of the real contact area of soft polymers” (49VF190053) was financially supported by German Federal Ministry for Economic Affairs and Climate Action (BMWK). References [1] Jacobs T D B and Martini A 2017 Measuring and Understanding Contact Area at the Nanoscale: A Review Appl. Mech. Rev. 69 060802. [2] Persson B N J, Albohr O, Tartaglino U, Volokitin A I and Tosatti E 2004 On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion J. Phys. Condens. Matt 17 R1. [3] Persson B N J 2006 Contact mechanics for randomly rough surfaces Surf. Sci. Rep. 61 201-227. [4] Klüppel M and Heinrich G 2000 Rubber Friction on selfaffine Road Tracks Rubber Chem. Technol. 73(4) 578- 606. [5] Hertz H 1881 Über die Berührung fester elastischer Körper J. F. reine u. angew. Math. 92 159-171. [6] Johnson K L, Kendall K and Roberts A D 1971 Surface Energy and the Contact of Elastic Solids Proc. R. Soc. Lond. A Math. Phys. Sci. 324(1558) 301-313. [7] Popov V L 2017 Contact Mechanics and Friction Springer Heidelberg Dordrecht London New York 2 nd edn.. [8] Johnson K L 1955 Surface interaction between elastically loaded bodies under tangential forces Proc. Roy. Soc. A320 531-548. [9] Banks H T, Hu S and Kenz Z R 2011 A Brief Review of Elasticity and Viscoelasticity for Solids Adv. Appl. Math. Mech. 3 1-51. [10] Chen W W, Wang Q J, Huan Z and Luo X 2011 Semi- Analytical Viscoelastic Contact Modeling of Polymer- Based Materials J. Tribol. 133 041404. [11] Greenwood J and Williamson J 1966 Contact of Nominally Flat Surfaces Proc. R. Soc. A 295(1442) 300-319. [12] Bush A W and Gibson R D 1975 The Elastic Contact of a Rough Surface Wear 35(1) 87-111. [13] Shull K R 2002 Contact mechanics and adhesion of soft solids Mater. Sci. Eng. R36 1-45. [14] Wu C-E, Lin K-H and, Juang J-Y 2016 Hertzian load-displacement relation holds for spherical indentation on soft elastic solids undergoing large deformations Tribol. Int. 97 71-76. [15] Dintwa E, Tijskens E and Ramon H 2008 On the accuracy of the Hertz model to describe the normal contact of soft elastic spheres Granul. Matter 10 209-221. [16] Lin Y-Y and Chen H-Y 2006 Effect of Large Deformation and Material Nonlinearity on the JKR (Johnson-Kendall- Roberts) Test of Soft Elastic Materials J. Polym. Sci. B44(19) 2912-2922. [17] Feng Z Q, Peyraut F and Labed N 2003 Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies Int. J. Eng. Sci. 41(19) 2213-2225. Science and Research 36 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010 vertheless, there are two distinct approaches that can be employed to obtain approximate values for the real contact area based on d f . Since these approximations are based on empirical evidence, it is yet to be determined whether they can be applied to other systems and, in particular, to different length scales. The most evident possibility is to calculate d c from d f via Equation (7), that can be derived from Figure 3a. This produces a satisfactory agreement between the simulation results and Equation (2), as illustrated by the correlation in Figure 6a. (7) Nevertheless, a somewhat stronger agreement between the A-values calculated from the simulations and approximated from d f (Figure 6c) can be achieved by employing the observed empirical relationship p ≈ 0.83 d f (Figure 6b) in conjunction with the precise analytical relationship between p and a 2 (Equation (8)). (8) Conclusion The empirical simulations demonstrate that, under the specified conditions of investigation (nanometre level, PE, high velocities and short-term deformation behaviour), despite significant deviations from the fundamental assumptions (i.e., small deformations and ideal elastic behaviour), the Hertz Equations (2) and (4) provide a reasonable approximation for the real contact area. Nevertheless, it should be noted that the Young’s modulus in Equation (4) is not identical to the usual Young’s modulus that can be obtained from experimental tensile tests or bulk simulations. Furthermore, the indentation depth d is defined from the position of the first atomic contact, which is not easily accessible in experiments. Further approximations have to be utilized to calculate the real contact area from more customary values for the indentation depth d f , defined from increasing forces. Based on findings of these initial investigations, future work will address remaining questions, including: a) How can a suitable Young’s modulus be determined for predictions? b) Are the findings also transferable to other soft polymers and conditions (velocity, length scales)? c) To what extent do alternative models based on the Hertz model (including, but not limited to, time dependency, adhesion, rough surfaces, friction, and so forth) apply to soft polymers? H 5 I3 4 5 I E 67 8 4 F J IKALMM7 N 8 OLMP N 4 Q W J N X A9 RSTU999 J OLV? N 9 [18] Dimitri R, De Lorenzis L, Scott M A, Wriggers P, Taylor R L and Zavarise G 2014 Isogeometric large deformation frictionless contact using T-splines Comput. Methods Appl. Mech. Eng. 269 394-414. [19] Lin D C, Shreiber D I, Dimitriadis E K and Horkay F 2009 Spherical indentation of soft matter beyond the Hertzian regime: numerical and experimental validation of hyperelastic models Biomech. Model. Mechanobiol. 8 345-358. [20] Frenkel D and Smit B 2023 Understanding molecular simulation: from algorithms to applications Academic Press London California 3 rd edn.. [21] Zhang J, Wang Z, Yan Y and Sun T 2016 Concise Review: Recent Advances in Molecular Dynamics Simulation of Nanomachining of Metals Curr. Nanosci. 12(6) 653-665. [22] Yang C and Persson B N J 2008 Contact mechanics: contact area and interfacial separation from small contact to full contact J. Phys.: Condens. Matter 20 215214. [23] Solhjoo S and Vakis A I 2015 Definition and Detection of Contact in Atomistic Simulations Comput. Mater. Sci. 109 172-182. [24] Vishnubhotla S B, Chen R, Khanal S, Hu X, Martini A and Jacobs T D B 2019 Matching Atomistic Simulations and In Situ Experiments to Investigate the Mechanics of Nanoscale Contacts Tribol. Lett. 67 97. [25] Li S, Li Q, Carpick R W, Gumbsch P, Liu X Z, Ding X, Sun J and Li J 2016 The Evolving Quality of Frictional Contact with Graphene Nature 539(7630) 541-545. [26] Spijker P, Anciaux G and Molinari J F 2012 The Effect of Loading on Surface Roughness at the Atomistic Level Comput. Mech. 50(3) 273-283. [27] Luan B and Robbins M O 2005 The Breakdown of Continuum Models for Mechanical Contacts Nature 435(7044) 929-932. [28] Abraham M J, Murtola T, Schulz R, Páll S, Smith J C, Hess B and Lindahl E 2015 GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers SoftwareX 1-2 19- 25. [29] Berendsen H J C, van der Spoel D and van Drunen R 1995 GROMACS: A message-passing parallel molecular dynamics implementation Comput. Phys. Commun. 91(1- 3) 43-56. [30] Gromacs. Fast. Flexible. Free 2024 https: / / www.gromacs.org/ . [31] Hockney R W, Goel S P and Eastwood J 1974 Quiet High Resolution Computer Models of a Plasma J. Comput. Phys. 14(2) 148-158. [32] Páll S and Hess B 2013 A flexible algorithm for calculation pair interactions on SIMD architectures Comput. Phys. Commun. 184(12) 2641-2650. [33] Bussi G, Donadio D and Parrinello M 2007 Canonical sampling through velocity rescaling J. Chem. Phys. 126(1) 014101. [34] Berendsen H J C, Postma J P M, van Gunsteren W F, Dinola A and Haak J R 1984 Molecular dynamics with coupling to an external bath J. Chemi. Phys. 81(8) 3684- 3690. [35] Daura X, Mark A E and van Gunsteren W F J 1998 Parametrization of Aliphatic CHn United Atoms of GRO- MOS96 Force Field J. Comput. Chem. 19(5) 535-547. [36] Oostenbrink C, Villa A, Mark A E and van Gunsteren W F 2004 A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6 J. Comp. Chem. 25(13) 1656-1676. [37] Schuler L D, Daura X and van Gunsteren W F J 2001 An Improved GROMOS96 Force Field for Aliphatic Hydrocarbons in the Condensed Phase J. Comput. Chem. 22(11) 1205-1218. [38] Fritz S 2021 Considering semi-crystallinity in molecular simulations of mechanical polymer properties - using nanoindentation of polyethylene as an example CMMS21(1) 35-50. [39] Gromacs 2024.2 Manual (https: / / doi.org/ 10.5281/ zenodo.11148638). [40] Cheng L, Xia X, Scriven L E and Gerberich W W 2005 Spherical-tip indentation of viscoelastic material Mech. Mater. 37 213-226. Science and Research 37 Tribologie + Schmierungstechnik · volume 72 · issue 2/ 2025 DOI 10.24053/ TuS-2025-0010