Computational homogenization of soft matter friction: Isogeometric framework and elastic boundary layers

Tam metin

(1)INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2014; 100:953–981 Published online 17 September 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.4778. Computational homogenization of soft matter friction: Isogeometric framework and elastic boundary layers ˙I. Temizer*,† Department of Mechanical Engineering, Bilkent University, 06800 Ankara, Turkey. SUMMARY A computational contact homogenization framework is established for the modeling and simulation of soft matter friction. The main challenges toward the realization of the framework are (1) the establishment of a frictional contact algorithm that displays an optimal combination of accuracy, efficiency, and robustness and plays a central role in (2) the construction of a micromechanical contact test within which samples of arbitrary size may be embedded and which is not restricted to a single deformable body. The former challenge is addressed through the extension of mixed variational formulations of contact mechanics to a mortar-based isogeometric setting where the augmented Lagrangian approach serves as the constraint enforcement method. The latter challenge is addressed through the concept of periodic embedding, with which a periodically replicated C 1 -continuous interface topography is realized across which not only pending but also ensuing contact among simulation cells will be automatically captured. Two-dimensional and three-dimensional investigations with unilateral/bilateral periodic/random roughness on two elastic micromechanical samples demonstrate the overall framework and the nature of the macroscopic frictional response. Copyright © 2014 John Wiley & Sons, Ltd. Received 25 November 2013; Revised 7 June 2014; Accepted 4 August 2014 KEY WORDS:. contact homogenization; friction; isogeometric analysis; periodic embedding. 1. INTRODUCTION The characterization of the link between microscopic surface roughness and the macroscopic frictional interaction among contacting surfaces has been one of the major challenges in tribology. In the regime of continuum mechanics where the present emphasis lies, efforts have historically been concentrated on the contact of metals, and the approaches developed therein have constituted a basis for the modeling of friction among hard non-metallic surfaces as well [1]. While soft matter friction shares similarities with friction among hard materials [2], notably through adhesion-type contributions, a number of different challenges arise that are primarily associated with the material response, and the addressing of which requires an efficient computational framework that is governed by contact detection and treatment. The development of such a framework in the context of contact homogenization is the main goal of the present work. For metals, the shearing of adhesive junctions across the rough interface and the plowing of one of the paired surfaces by the asperities of the other are among the major mechanisms responsible for the macroscopically observed frictional response, which is also strongly influenced by microscopic wear and lubrication [3, 4]. Fundamental micromechanical models for friction are typically amenable to an analytical contact homogenization treatment based on simplified geometric constructions, such as a single rigid asperity indenting a deformable surface, combined with the modeling of the plowing and adhesive forces [5]. These asperity-level interaction models may be comple*Correspondence to: ˙I. Temizer, Department of Mechanical Engineering, Bilkent University, 06800 Ankara, Turkey. † E-mail: temizer@bilkent.edu.tr Copyright © 2014 John Wiley & Sons, Ltd..

(2) 954. ˙I. TEMIZER. mented by numerical approaches that have the potential to resolve the complex material response in the vicinity of a contact zone [6–8]. Numerical approaches also assist in the incorporation of statistical models of multiscale random roughness toward a deeper understanding of the interfacial topography influence on the homogenized response. Beginning with the influential work of [9], a natural starting point for such investigations is to initially consider normal contact alone—see [10, 11] for reviews of the multiscale modeling of roughness and contact and [10, 12–14] for representative recent studies. A successful resolution of the microscopic contact area during normal loading enables the subsequent analysis of the macroscopic tangential response and, hence, of friction—see [15] for an early example and [16, 17] for further studies. Because of the complex molecular structure of soft matter, the domain of modeling often lies in continuum mechanics, and in view of the highly nonlinear response (see below), the finite element method appears to be the most relevant numerical technique. For hard materials such as metals, on the other hand, a wider range of modeling and simulation choices are available for the multiscale analysis of contact in addition to the finite element method [18, 19], such as molecular dynamics [20, 21], the discrete element method [22, 23], and the boundary element method [24, 25], as well as semi-analytical methods [26, 27]. Driven by technological interest, over the past few years, significant experimental effort has also been focused toward a better understanding of sliding (or kinetic) friction that is associated with biological and synthetic soft matter. For instance, skin friction plays an important role in tactile perception [28–31] while elastomeric friction governs traction [32–35]. In comparison with hard materials, the major difference in the response of such soft matter is the intrinsic viscoelasticity of the material and the very large deformations that the material can sustain without permanent strains. The microscopic viscous dissipation resulting from the continuous and cyclic loading of the boundary layers in the vicinity of the contact zone by the asperities of the rough surface has the potential to significantly augment the macroscopically observed friction. This essentially leads to a characteristically non-Amontons and non-Coulomb type behavior where the macroscopic friction depends, respectively, on the contact pressure and the slip velocity. The early works by Schallamach [36] and Grosch [37] had already highlighted this relation for rubber, a widely employed elastomeric material, while subsequent work additionally noted the significance of adhesive contributions [38, 39], which also intrinsically depend on the surface topography [40, 41]. Although the viscous and adhesive contributions may not always provide a complete basis for soft matter friction [42], they highlight the tribological complexity as well as the central role that surface roughness plays in the homogenized interface response. In order to numerically capture this response accurately, the microscopic contact area must be resolved properly, which necessitates reaching a steady-state cyclic viscoelastic response by ensuring sufficiently long sliding contact on the microscale. Assuming linear viscoelasticity and a unilateral roughness that is assigned to a rigid surface, most of these challenges have been elegantly addressed for rubber friction across multiple length scales within a fractal setting [43], possibly with anisotropy [44], and the developed approaches make fairly accurate predictions [45, 46]. While limited in number, extensions of the setup to a materially and geometrically nonlinear finite element framework, but preserving unilateral roughness on a rigid surface, have also been pursued in two-dimensional [47] and three-dimensional settings [48] as well as in the case when the interface heterogeneities are mobile rigid third bodies rather than the asperities of a rough rigid surface [49]. The common assumption of a unilateral roughness assigned to a rigid surface employed in the aforementioned studies is often justified by the frictionless normal contact mechanics of two linearly elastic bodies with rough surfaces. The primary step in the analysis of this problem is the calculation of an equivalent roughness, assigned to a rigid surface, and an equivalent elastic modulus, assigned to an elastic body [50, 51]. However, the theoretical support for this reduction is invalid under frictional sliding [52] and, in particular, in the nonlinear setting that is required for the treatment of soft matter contact. Moreover, while the summarized approaches constitute a suitable first step toward the characterization of friction among soft materials in general, the necessity to consider two deformable surfaces in a kinematically and geometrically nonlinear setting places stringent requirements on the contact algorithms employed, demanding a higher degree of simultaneous accuracy, efficiency, and robustness [53, 54], which poses a challenge. As an additional challenge, long Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(3) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 955. dragging durations as well as the appropriate statistical sampling of roughness requires the construction of large interface topographies, which is not a major issue with unilateral roughness where the rough surface is not numerically discretized but can easily be prohibitive when both surfaces are deformable if, for instance, the stationary surface is chosen to be sufficiently large so that the sliding one never leaves its span. The construction of a computational contact homogenization framework that is capable of addressing these two central and interactive challenges for soft matter friction is the present goal. The remainder of this work is organized as follows. In Section 2, a mixed variational framework for mortar-based frictional contact is extended to the augmented Lagrangian (AL) setting in order to improve accuracy, efficiency, and robustness. The framework benefits from an isogeometric setting that allows using higher-order contact variable discretizations that emanate from the geometric description of the micromechanical samples. Representative test cases demonstrate the satisfactory performance of the contact algorithm in comparison with pure penalty or Uzawa-staggered forms. Section 3 introduces the concept of periodic embedding that enables the construction of a micromechanical contact test within which arbitrary samples sizes may be employed and arbitrarily long durations of microscopic sliding contact can be achieved. Through geometric periodicity constraints and cell-periodic closest-point projection, a periodically replicated C 1 -continuous interface topography will be realized across which not only pending but also ensuing contact among simulation cells will be automatically captured. The emphasis of the work is on the construction and the demonstration of the computational contact homogenization framework. Consequently, the numerical examples in Section 4 will concentrate on how the microscale frictional response is reflected to the macroscale after being altered by the interface topography rather than on the augmentation of the frictional response through the accumulation of viscous bulk or adhesive interface dissipation across the scales. For this purpose, a purely elastic bulk response and an Amontons–Coulomb type microscopic interface friction will be admitted. Here, following most studies on rubber friction, dynamic effects in soft matter friction [39, 55, 56] will be omitted. Finally, an outlook toward a more comprehensive computational tribology approach for soft matter contact is provided in Section 5. 2. MIXED VARIATIONAL FRAMEWORK FOR ISOGEOMETRIC CONTACT 2.1. Mortar-based contact with friction Two-dimensional and three-dimensional mortar-based isogeometric schemes for frictional contact have been presented in [57] and [58], respectively, where the penalty method was employed for enforcing the contact constraints. A mixed variational framework, based on the classical work [59], that delivers such mortar-based schemes while benefiting from isogeometric discretizations has been discussed in [60]. Therein, the AL method was employed for constraint enforcement in a staggered setting via Uzawa iterations [61]. Presently, this variational framework is extended to the original AL scheme [62, 63]. The presentation is compact, in particular omitting linearization, because this extension follows mostly from the tangential contribution to the mixed contact functional in the Uzawa setting and the framework in [63]. The procedure for normal contact is similar and, as a special case, also delivers the isogeometric contact scheme of [64]. For a compact presentation, standard computational contact mechanics terminology is employed [53, 54]. In particular, all contact integrals will be evaluated on the image oc of the potential contact interface mapped to the reference configuration of the slave body via Z dA: hi D (2.1) oc. This choice is accompanied by the decomposition of the contact (Piola) traction acting on the slave into normal and tangential parts such that the contact pressure is pN > 0 while pT ˛ (˛ D 1; 2) are the covariant components of the tangential traction. If the outward unit normal is and the contravariant basis vectors are a˛ at the closest-point projection y of a slave point x onto the master, p D pN pT ˛ a ˛ Copyright © 2014 John Wiley & Sons, Ltd.. (2.2) Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(4) ˙I. TEMIZER. 956. is the total traction acting on the slave. This kinetic construction is complemented by the kinematic quantities of normal (gN ) and tangential (gT˛ ) gap measures: gN D .x y/ and gT˛ D ˛;nC1 ˛;n numerically represents the change in the convected curvilinear coordinates ˛ of y from time c C ı GTc to the weak form of the step n to n C 1. With this setup, the contact contribution ı G c D ı GN linear momentum balance has the normal and tangential parts ˛ ˝ c ı GN D hıgN pN i; ı GTc D ıgT˛ pT ˛ : (2.3) An independent Lagrange multiplier contribution has a similar decomposition: ı Lc D ı LcN C ı LcT . In order to discretize the contact contributions, the approach pursued will particularly benefit from the mortar schemes advocated in [65, 66] and from the idea of choosing the integration surface coincident with the slave [67] instead of carrying out segmentation for accurate integration. The choice of the integration surface delivers implementation convenience for the class of problems considered but sacrifices machine precision accuracy for the flat interface patch test. However, this is not a significant shortcoming because Temizer [60] demonstrates that patch tests can be passed to very good accuracy even for inclined and curved interfaces. Sacrificing machine precision is acceptable because segmentation approaches applied to the types of non-flat NURBS discretizations encountered in this work also deliver high accuracy but not machine precision [68, 69]. For frictional mortar contact approaches employing alternative interface variable discretizations and constraint enforcement algorithms, the reader is referred to recent works [70–72] and reviews [73, 74]. See also [75, 76] for non-mortar isogeometric contact. 2.2. Normal contact functional A mixed kinematic variable N and a normal Lagrange multiplier N are introduced and, together with pN , assigned an interface discretization using the shape functions N I of the slave surface: X X X I I pN D N I pN ; N D N I N ; N D N I IN : (2.4) I. I. I. In terms of the normal contact variables, one may introduce the functional 1 ˝ 2˛ C hpN .N gN /i hN N i CN ŒgN ; pN ; N ; N D N N 2 where N is the normal penalty parameter. For notational compactness, it is useful to define ˛ ˛ ˝ ˝ ˝ ˛ g IN D N I gN ; ÎJ D N I N J ; w I D N I : Now, upon making use of the mortar-based discretization (2.4), CN takes the discrete form ´ ! μ X X X X 1 I I IJ J I IJ J I IJ J CN D ˆ N C pN ˆ N g N N ˆ N : N N 2 I. J. J. (2.5). (2.6). (2.7). J. Clearly, only the projected normal gap g IN is a known input to CN . For numerical purposes, although not required [60], it is advantageous to eliminate the coupling among the contact DOFs induced by the overlap matrix ÎJ through row-sum lumping, following similar ideas for the treatment of the mass matrix ([77], p.444), thereby obtaining ³ X² 1 I I I I I N N : CN D w I N C pN g IN IN w I N w I N (2.8) 2 I. c The variation of CN with respect to g IN delivers the contribution ı GN to ı G c while the remaining I I I variations with respect to pN , N and N induce, respectively, the contact equalities X ® ¯ I I I I N D 0: D g IN =w I ; pN D IN C N N ; ı LcN D ıIN w I N (2.9) I. ı LcN is the weak form contribution associated with the multiplier variations and indicates the weak enforcement of gN D 0 on oc , satisfied at convergence. It is emphasized that the mixed Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(5) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 957. kinematic variables only serve to establish a link between the actual interface geometry and this weak constraint enforcement. They do not appear in the system of equations: the only unknowns are the displacement and Lagrange multiplier DOFs. Overall, this approach benefits from the underlying isogeometric discretization that guarantees N I > 0, thereby allowing the application of the same contact algorithm independent of the discretization order [58]. Higher-order Lagrange discretizations do not display this property and hence require special mortar-based treatment [71, 78]. CN and, therefore, ı LcN are incomplete in their present forms because the discrete functional (2.8) I > 0. The following modification expliis only applicable to I belonging to the active set A: pN citly takes into account the active set and resets a Lagrange multiplier to zero upon contact loss while ensuring the continuity of the contributions to the finite element residual through the activeto-inactive transition, that is, as the contact state switches between active and inactive [63, 79]: ³ X² ³ X² 1 1 –1 I I I I I I I I I I I I I CN D w N : N N w N C pN w N g N N w N C 2 2 N N I 2A I …A (2.10) 2.3. Tangential contact functional The tangential contribution introduces the mixed kinematic variables T˛ and tangential Lagrange multipliers T ˛ in addition to pT ˛ . The relevant tangential contact functional is E ˝ ˛ ˛ ˝ 1 D ˇ (2.11) CT ŒgT˛ ; pT ˛ ; T˛ ; T ˛ D T T˛ a˛ˇ T C pT ˛ .T˛ gT˛ / T ˛ T˛ ; 2 where T is the tangential penalty parameter and a˛ˇ are the covariant metric components on the master. The tangential contact variables are subject to an interface discretization similar to that of their normal counterparts in (2.4). In addition to (2.6), the definitions ˝ I ˛ ˝ I ˛˛ ˝ I ˛ N a˛ˇ ˛;I IJ J I (2.12) g T D N gT ; a˛ˇ D N a˛ˇ N ; m˛ˇ D ; kvI k2 D v˛I m˛ˇ;I vˇI wI are introduced where m˛ˇ;I is associated with the inverse of mI˛ˇ and vI symbolically indicates a tangential vector with covariant components v˛I . The discrete form of CT is then ´ ! μ X X X X 1 ˇ;J T T˛;I IT ˛ CT D C pTI ˛ ÎJ TJ ˛ g ˛;I ÎJ T˛;J (2.13) aIJ ˛ˇ T T 2 I. J. J. J. while the final lumped form that is convenient for numerical purposes is ³ X² 1 ˇ;I ˛;I I I ˛;I I I I I I ˛;I T T w m˛ˇ T C pT ˛ w T ˛ g T : CT D T ˛ w T 2. (2.14). I. In addition to delivering the contribution ı GTc to ı G c , the variations of CT effectively induce T˛;I D I g ˛;I T =w and, thereby, . ° ± X ˛;I I c I pTI ˛ D IT ˛ C T mI˛ˇ g ˇ;I ; ı L D ı w g D 0: T T˛ T T (2.15) I. In the first iteration of step n C 1, IT ˛ inherits its value I;n T ˛ from step n. Let S A denote the set associated with stick. Above, I 2 S was implicitly assumed. To handle active-to-inactive or stick-to-slip transitions, the functional CT must be extended [63, 79]: ° ± P (1) Active-to-Inactive: The inactive set contribution is I …A 12 T –1 w I kIT k2 . If contact is reinitiated in a next iteration in step n C 1, IT ˛ must initially be reassigned I;n T˛. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(6) ˙I. TEMIZER. 958. (2) Stick-to-Slip: The Coulomb criterion does not have a unique numerical expression in the mortar setting because of different constructions of the mortar algorithms [60, 66]. Here, I I with I D kpIT k I 6 0 is employed ° k as the friction coefficient ± and D kpN (fixed P during variation). I 2AnS 12 T–1 w I kIT k2 2 I kpIT k C I I is added when I > 0, effectively assigning IT ˛ its slip value I s˛I with s˛I D pTI ˛ =kp IT k as the trial slip direction. To benefit from the standard tangential variable update algorithm [53, 54], the metric variation is omitted from ı CT , leading to the well-known loss of symmetry in linearization—see [80] for an improved update algorithm that preserves symmetryPin stick. Additionally, the active-to-inactive contribution to the weak form may be replaced by I …A ıIT ˛ ¹T–1 w I IT ˛ º and the stick-to-slip P contribution to the ı LcT portion by I 2AnS ıIT ˛ ¹T–1 w I .IT ˛ I s˛I /º. While such a modification would sacrifice the continuity of ı LcT through the transitions [63], the linearization of ı LcT D 0 is simplified with no numerical reliability degradation for the class of problems investigated, which will be justified by the numerical examples that follow.. 2.4. Representative test case I: rotating frictional contact In order to demonstrate the AL formulation outlined earlier and to highlight its advantages toward the realization of the micromechanical investigations where persistent frictional sliding contact conditions prevail, a representative test case is considered where an elastic body is compressed against another in frictional contact, and subsequently, an increasing twisting torque is applied with displacement control until sliding dominates at the interface. Instances from the simulation are provided in Figure 1. The problem setup, including the loading/geometry parameters and the material model, is borrowed from [60] and is not repeated. In this example, as well as in upcoming ones, the time step size is refined/coarsened according to a basic adaptivity criterion such that it is halved upon convergence failure and, whenever a refinement has been carried out, doubled after a preset number of consecutively successful steps. A number of numerical observations regarding this test case are summarized in Figure 2. An appropriate mesh is chosen from three resolutions. The coarse resolution is represented in Figure 1 while the medium/fine resolution corresponds to two/three times the coarse one in each parametric direction. The twisting torque evolution is sufficiently accurately captured at a medium mesh resolution, which is chosen as the default. When there is no friction, the twisting torque is zero, apart from negligibly small oscillations about zero, which can be further reduced by knot refinement or order elevation within an isogeometric setting [57, 58]. For small k, sliding conditions quickly dominate the interface response, and hence, the twisting torque saturates quickly. The saturation is delayed with increasing k. The default is chosen as k D 1. On the basis of these default choices, two sets of investigations are carried out. First, the AL scheme is degraded to a pure penalty form where N is varied at a fixed N =T ratio. Second, the AL scheme is degraded to a staggered form in which Uzawa iterations augment the Lagrange multipliers to a tolerance T OL where the error. Figure 1. Simulation instances from the rotating frictional contact problem of Section 2.4. The bodies are discretized (coarse resolution) with quadratic NURBS in each parametric direction. The red spheres represent the corresponding control points. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(7) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 959. Figure 2. (a and b) The mesh resolution and coefficient of friction effects are summarized for the AL scheme. (c and d) Pure penalty and Uzawa-staggered forms of the AL scheme are assessed.. corresponds to the 2-norm of the multiplier updates scaled by that of the multipliers from the last iteration. Well-known observations regarding these alternative approaches are briefly summarized: (1) Pure penalty form: Choosing an appropriate penalty number is not straightforward, with too small or too large values typically leading to convergence problems. Presently, the two smallest choices fail to converge even with time step adaptivity while the largest choice only leads to convergence after aggressive refinement such that the total number of calls to the linear system solver exceeds twice that of the AL scheme for the same setup. The difficulties in identifying appropriate penalty parameters lead to uncertainties in accuracy assessment. (2) Uzawa-staggered form: These uncertainties may be alleviated by augmentation with a sufficiently small Uzawa tolerance. Robustness is also enhanced because convergence may now be achieved for previously unsuccessful pure penalty setups. A large tolerance (T OL D 101 ) keeps the number of augmentations to a minimum but does not ensure an accurate result if the penalty parameters are not sufficiently large, while a small tolerance (T OL D 103 ) well approximates the AL solution but may lead to a large number of augmentations (> 10). In view of these observations, the AL scheme provides a framework within which, independent of the choices of the penalty parameters to a large extent, the solution is obtained with an optimal combination of accuracy, efficiency, and robustness. Overall, it provides satisfactory numerical reliability for the micromechanical contact investigations to be carried out in the upcoming sections. It is remarked that algorithmically consistent linearization is carried out such that asymptotically quadratic convergence rates are achieved. This is demonstrated in Appendix C, along with the number of Uzawa iterations needed for the staggered form. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(8) 960. ˙I. TEMIZER. 2.5. Representative test case II: sliding frictional contact In order to demonstrate the AL formulation further, the sliding contact of two thick-walled shells is considered as a second test case. Similar to the first test case, the setup is borrowed from an earlier work [81], so the details are not repeated. However, it is extended to include friction. In the rotating contact problem, the active-to-inactive transitions occur only during the compression phase. During rotation, significant stick-to-slip transitions take place. In the present problem, both types of transitions take place simultaneously. The problem is very relevant to the goals of this work because the geometry of contact is reminiscent of two highly deformable asperities in sliding contact. Simulation instances in Figure 3 summarize the test case for k D 1:00. At this value of k, Figure 4 first investigates the effect of mesh resolution by monitoring the normal and tangential forces measured. The fine resolution represents sufficiently converged results. At this resolution, k is varied to monitor its effect. By construction, the frictionless case gives rise to a contact problem. Figure 3. Simulation instances for the sliding frictional contact of Section 2.5 with k D 1:00. The bodies are discretized (fine resolution) with quadratic NURBS in each parametric direction.. Figure 4. The mesh resolution and coefficient of friction effects are summarized for the problem of Figure 3 by monitoring the normal (FN ) and tangential (FT ) forces measured. The problem is solved only with the AL scheme. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(9) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 961. Figure 5. Simulation instances for different values of k at D 0:5 from the analysis of Figure 4. For all k values, initial ( D 0) and final ( D 1) configurations are the same (see Figure 3).. that is symmetric about D 0:5. With increasing k, the tangential resistance increases such that at k D 1:00, a tangential force must continuously be applied all the way until separation. This increasing tangential resistance is also observed from the simulation instances at D 0:5 for different k values, as summarized in Figure 5. In all cases, asymptotically quadratic convergence rates are achieved (Appendix C).. 3. COMPUTATIONAL CONTACT HOMOGENIZATION 3.1. Micromechanical testing procedure Once a specific frictional interaction is admitted between the contacting asperities of rough surfaces, the overall topography of the contact interface will alter the reflection of this microscale behavior onto the macroscale. On the macroscale, one observes a homogenized response that can be extracted via a micromechanical test. The testing procedure, summarized in Figure 6, reflects the mechanics of contact on the macroscale both numerically and physically; the aspects of which are individually addressed below. 3.1.1. Numerical setup. Numerically, the classical designation of slave and master [53, 54] will be transferred to samples associated with the contacting surfaces (blue: slave; yellow: master). These samples represent the interaction of thin layers of material in the vicinity of the contact interface (boundary layers [82]). The material response outside the boundary layers is not expected to affect the macroscopic interface response [52, 82]. The raw geometry of the samples exclude surface roughness. In a first step, the dimensions of the samples in the plane of macroscopic contact are chosen to represent the statistical characteristics of roughness, corresponding to the period of roughness in the simplest scenario of periodicity. The height of the samples should likewise be chosen appropriately, leading to height-sensitive results if not chosen sufficiently large—see [82, 83] and Section 4.3.1. In a second step, the surface roughness is assigned to the potential contact surfaces such that the raw surface corresponds to the mean plane of roughness (see Figure 6 and Section 3.2.1). Two types of surface roughness will be considered: (i) periodic (sinusoidal) and. Figure 6. The micromechanical testing procedure is summarized: (left) raw geometry before the assignment of surface roughness and (right) the heterogeneous sample geometry on which the boundary conditions are imposed (blue: slave; yellow: master). Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(10) 962. ˙I. TEMIZER. Figure 7. Linear Lagrange (L1 ) and quadratic NURBS (N 2 ) discretizations are compared at different resolutions, demonstrating the ability of isogeometric discretizations to rapidly capture the underlying smooth surface topography with an acceptable quality even at low resolutions.. (ii) Gaussian random (random-field model [84, 85]). Surfaces are characterized with respect to their out-of-plane (spatial) and in-plane (spectral) properties [86]. Among the out-of-plane statistical properties, only the root-mean-square (RMS) roughness will be controlled, following the generation procedure in [87]. With respect to the in-plane ones, it is noted that the employed periodic surfaces will be anisotropic while the random ones can display both isotropy and anisotropy. Only isotropic random surfaces will be employed in this work. In this case, an important consequence is that the macroscopic frictional response will be isotropic provided the surfaces are not subjected to macroscopic in-plane deformation (Section 3.1.2). However, macroscopic anisotropy can also be observed in the case of statistical anisotropy (Section 4.4.1). The computational framework will benefit from three advantages associated with isogeometric NURBS-based discretizations [88]. As noted in Section 2.2, one advantage is the non-negativity of isogeometric basis functions, which enables the uniform applicability of the contact algorithm to higher-order discretizations while most mortar-based approaches are based on linear Lagrange polynomials. A second advantage is the additional underlying smoothness of isogeometric discretizations that will be further employed in Section 3.2.1. This leads to the third advantage, namely, that the topography of the rough surface can be represented with high accuracy even at coarse discretizations (Figure 7). While the overall computational framework is to a large extent suitable for employing NURBS discretizations of arbitrary order, numerical examples are based on quadratic and cubic NURBS with unit weights (i.e., quadratic/cubic B-splines [89]) and a uniform knot vector because such a structure is favorable for enforcing geometric periodicity (Section 3.2.1). 3.1.2. Physical parameters. Physically, the macroscopic contact conditions will be projected onto the samples by imposing them as boundary conditions (BCs) on their outer surfaces. Specifically, the master will be held fixed while the slave is subjected to the macroscopic contact pressure p N and slip velocity vT (Figure 6). Because a classical Amontons–Coulomb frictional response is admitted Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(11) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 963. and the boundary layers are elastic, there is no rate-dependent dissipation mechanism so that kvT k does not affect the homogenized response [90]. The rough surfaces are subject to frictional contact conditions (Section 2) while periodic BCs will be imposed on the remaining inner surfaces (Figure 6): (3.1) x C x D F c X C X ; pC C p D 0: Here, x/X denotes the current/reference configuration position vector, F c represents the surfacial (in-plane) deformation of the macroscopic contact interface [8, 82, 87], and pC= indicates values of the Piola traction on periodically linked points (Figures 6 and 10). Explicitly, F c D e˛ ˝ E ˛ C ˝ . (3.2). where is the unit normal to the master at the macroscopic point of contact (Figure 6) and E ˛ (˛ D 1; 2) are tangential unit vectors, constituting an orthonormal basis together with , while e ˛ D F c˛ˇ E ˇ. (3.3). are tangential (non-unit) vectors on the deformed configuration (Figure 10). F c alters the statistical properties of the interface in two ways: (i) if F c˛ˇ D ı˛ˇ , then the RMS roughness will decrease (increase), and the period will increase (decrease) if the stretch is greater (less) than unity and (ii) if F c D I C E 1 ˝ E 2 , as a specific example, then the surfaces are sheared such that a surface that initially responds isotropically will display a deformation-induced anisotropy. The former will affect the macroscopic frictional response only quantitatively while the latter will also cause a qualitative change. It is noted that, for numerical convenience, the macroscopic deformations of the slave and master surfaces are presently assumed to be the same although they can be significantly different in general. Figure 10 employs surfacial deformation with F c˛ˇ D ı˛ˇ C 0:25. 3.2. Periodic embedding The micromechanical test will be decomposed into three phases: (1) deformation phase where F c is imposed; (2) compression phase where p N is imposed and maintained throughout the simulation; and (3) dragging phase where macroscopic slip is imposed on the slave sample. Figure 6 already depicts the embedding of the original samples among their periodic images throughout the micromechanical test, the need for which is highlighted in a sample computation summarized in Figure 8. The projection of the slave rough surface coincides perfectly with the master at the beginning of the simulation but not after the compression phase starts. Moreover, the slave sample must continuously interact with the master surface for a sufficiently long dragging duration in order to ensure an accurate characterization of the macroscopic frictional response. For both periodic and random surfaces, one numerical approach to enable continuous interaction is to choose the master surface sufficiently large so that the slave never leaves its span. However, this is inconvenient, particularly when both surfaces are deformable or if the generation of a large random master surface is required. Moreover, the determination of an appropriate size for the master surface is not straightforward because one. Figure 8. Cell-periodic sliding contact conditions are summarized in a two-dimensional setting. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(12) ˙I. TEMIZER. 964. typically monitors the macroscopic frictional response until a suitable moving time average converges (Section 4.2). Periodic embedding, which requires more than just imposing periodic BCs on the deformation, circumvents these difficulties by facilitating continuous slave–master sample interaction for arbitrarily long durations while enabling the use of arbitrarily sized samples—see [91] for an early example with rigid periodic samples. In order to achieve periodic embedding, the C 1 -continuity of the surfaces is advantageous for smooth contact interactions as well as for the projection algorithm used in detecting contact. The projection of a slave point onto the periodic images of the master is realized by cell periodicity and benefits from the continuity across the images that are enforced by geometric periodicity, which are discussed next. 3.2.1. Geometric periodicity. The mapping of the rough surfaces onto the raw geometry (Figure 6) is realized by a least-squares fit. Considering the slave surface as an example and recalling the P discretization X D I N I X I (2.4), the minimization is carried out with respect to control point positions X I :. 1 2. Z oc. .X /2 dA ! min:. (3.4). Here, is the desired vertical position at X , which is an input from the underlying periodic or random surface generation algorithm. Although not pursued in this work, a more accurate fit can be obtained by optimizing additionally with respect to control point weights of the underlying NURBS description [92], assumed to be unity throughout this work. Once the surface mesh fitting is completed, the displacements of the raw surface control points with respect to the mean plane are applied to the remaining bulk control points with a weight that decreases with the distance from the mean plane (Figure 6). A direct application of this minimization delivers incorrect geometric periodicity because the procedure does not even guarantee C 0 -continuity. Such results are compared with the desired correct geometric periodicity for periodic and random surfaces in Figure 9. C 1 -continuous geometric periodicity is achieved on the basis of the setting depicted in Figure 10 by adding penalization terms to. Figure 9. Two-dimensional periodic and random surface topography generation examples are provided, without (incorrect) and with (correct) the imposition of geometric BCs that ensure C 1 -continuity across the boundaries among the periodic images of the central cell. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(13) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 965. Figure 10. The variables employed in the enforcement of geometry-periodic and cell-periodic sliding contact conditions are summarized on undeformed and deformed master surfaces. Here, F c˛ˇ D ı˛ˇ C 0:25.. the minimization problem (3.4) that are associated with continuity constraints. Considering a sample set of points along a given axis with position vectors A (or, B), the first constraint enforces C (3.5) C 0 -continuity: A A D 0: Only the out-of-plane components, with respect to the mean plane, are considered because the inplane periodicity is satisfied by the construction of the raw geometry (Figure 6). The C 0 -continuity of the resulting geometry is maintained throughout the deformation of the micromechanical samples by the periodicity condition (3.1). Similarly, employing the out-of-plane components only, the second constraint enforces ¯ ® C (3.6) D 0: C 1 -continuity: A AC o C A Ao Here, the subscript ./o refers to the control point that is adjacent to the one on the boundary with respect to the perpendicular parametric direction [89]. Along the vertical edges in Figure 10, as well as within a surface, horizontal C 1 -continuity is ensured by the B-spline construction. The constraint stated earlier additionally ensures C 1 -continuity across the edge, that is, across patches. A similar constraint in terms of the point set B ensures C 1 -continuity across the horizontal edges. Using the notation introduced in Figure 10, C 1 -continuity across the sample boundaries is maintained throughout the deformations by complementing (3.1) with the enforcement of the generalized condition C (3.7) x xC o C x xo D 0 which is an adaptation of the method employed in [93]—see [94, 95] for alternative schemes for shells and membranes. These methods employ non-periodic (or, open) knot vectors [89] that render the control points on the boundary interpolatory. Geometric periodicity can also be ensured by using periodic knot vectors with which a k-th order discretization automatically ensures C k1 -continuity not only within the cell mesh but also across the cell images [96]. The reference configurations of three-dimensional samples generated via geometric periodicity are provided in Figure 11 for periodic and random surface topographies. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(14) ˙I. TEMIZER. 966. Figure 11. Three-dimensional periodic and random surface topographies with geometric periodicity.. 3.2.2. Cell periodicity. As the slave sample slides over the master, the computation of the contact contributions to the weak form requires the determination of the closest-point projection of a slave point in a first step. Because it is not known beforehand on which image of the master the projection lies in, it is necessary to search an appropriate set of images and subsequently carry out coupling between the moving slave and the stationary central master cell. Such cell-periodic projection is carried out on the basis of the setting depicted in Figure 10, independent of whether or not there is macroscopic surfacial deformation. The starting point is the use of a global parametrization of the master surface. The central master cell is a NURBS patch with parametric convected coordinates ˛ 2 Œ0; 1. Periodic images of the master cell therefore generate a master surface such that the addition of each image along a direction increases (or decreases) the parametric span by unity. Hence, a master point y is a function of the parametric position D . 1 ; 2 / such that the closest-point projection minimizes the distance of the slave point to the master on the basis of the standard criterion f˛ ./ D a˛ .x y/ D 0. (3.8). where a˛ are covariant basis vectors on the master. It is convenient to determine the parametric coordinates ˛ for the origin of a master image cell within which the point lies, such that ˛ D 0 corresponds to the central cell’s origin. Using bc as the floor function, these are ˛ D b ˛ c:. (3.9). On the basis of this setting, the algorithmic procedure for locating the projection point nC1 with global minimum distance at time step n C 1 is as follows: (1) Last projection cell: The last projection cell’s origin from step n is recorded: ˛;n D b ˛;n c. (2) Genetic search: In particular, for random rough surfaces, the direct application of a Newton– Raphson algorithm to (3.8) is not desirable because it is unlikely to converge to the global minimum. Instead, a genetic search (Appendix A) that only evaluates the cost function f˛ is carried out in a preliminary step by sampling the parametric domain ˛ 2 Œ ˛;n 1; ˛;n C2, that is, the one-cell neighborhood of the last projection cell. The algorithm delivers an updated 1 initial guess ˛;nC 2 . Within the algorithm, the evaluation of all master quantities (such as a˛ or y) requires the parametric coordinate to lie within the knot span Œ0; 1 of the master NURBS patch. This is realized with a two-stage algorithmic shift: (a) Shift-to-span: The original coordinates are shifted to the central master cell span, with ˛ D ˛ ˛. which all quantities are evaluated using the central master cell mesh: shift (b) Shift-to-image: Among the master quantities, one obtains the projection point y shift on the central mesh, which must be shifted back to the original image cell in order to evaluate the residual (3.8): y D y shift C Lo ˛ e ˛ . Here, Lo denotes the undeformed in-plane dimension of the master, assumed equal in both directions without loss of generality, so that kLo e ˛ k is indicative of the new length of the edge along E ˛ . Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(15) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 967. Figure 12. Cell-periodic sliding contact conditions are summarized in a three-dimensional setting. The images of the slave cell, appearing in Figure 8, are not displayed for clarity. The slave sample does not slide directly upwards but rather at an angle of 30ı . 1. (3) Gradient search: Using the updated guess ˛;nC 2 as the starting point and making use of the shifting algorithm, the classical Newton–Raphson algorithm is applied to (3.8) where the continuity across the boundaries of the master cells is crucial, ensured by geometric periodicity. Within the algorithm, the iteration error is checked with respect to the update

(16) ˛ for ˛ , which converges to ˛;nC1 . Quantities evaluated at the projection point, required for the residual and the tangent associated with the contact algorithm in Section 2, which lead to ˛;nC1 coupling among the slave and master meshes, are then evaluated at shift . A sample three-dimensional cell-periodic sliding contact computation is shown in Figure 12. As in the two-dimensional version (Figure 8), the central slave cell is continuously shifted to the one-cell neighborhood of the central master cell for visual clarity (Appendix B).. 4. ELASTIC BOUNDARY LAYERS 4.1. Material model On the basis of the mortar framework of Section 2 for handling frictional contact and on the periodic embedding approach of Section 3 to model and detect the persistent sliding contact between the micromechanical test samples, the macroscopic contact response of rough boundary layers will be investigated next. The emphasis is on the ability of the computational framework to capture the influence of the contact interface topography on the macroscopic frictional response, rather than on demonstrating an ability to predict the tribologically complex contact response of soft materials (cf. [90]). As such, the constitutive response of the boundary layers is limited to finite elasticity so that the only source of microscopic dissipation is due to friction. Consequently, the macroscopic frictional response is expected to retain the Coulomb character of the microscopic friction although it will possibly display a non-Amontons response, depending on the details of the interface topography. The hyperelastic response is modeled according to the isotropic Ogden material model, with the specific parameters reported in [97]. The controlling parameters in this model are the linear elasticity variables of the two bodies, namely, the shear ( .I / ) and bulk ( .I / ) moduli. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(17) ˙I. TEMIZER. 968 4.2. Macroscopic frictional response. Periodic and isotropic random rough surfaces will be employed (Figure 11), both of which display an isotropic macroscopic frictional response in the absence of macroscopic surfacial deformation. In this case, as well as in two dimensions, it is convenient to characterize the macroscopic frictional response via an instantaneous macroscopic coefficient of friction kI D. kf T k fN. (4.1). where fN (constant for a given p N ) and f T are the normal and tangential (with respect to ) forces measured on the outer surface of the slave (Figure 6). Whenever there are significant oscillations in f T , the macroscopic coefficient of friction is identified as the time average of its instantaneous counterpart, starting at a time to to exclude the initial transition to macroscopic slip and carried out for a sufficiently long duration Tavg to ensure saturation to a limit: 1 kD Tavg. toZ CTavg. kI dt :. (4.2). to. In the case of anisotropy, without attempting a rigorous characterization of the anisotropic frictional response [98–100], k will represent the magnitude of frictional resistance with direction d k D cos.k /E 1 C sin.k /E 2 such that f T D kfN d k. (4.3). where f T is the time average of f T . It is remarked that k does not necessarily align with the macroscopic slip velocity vT D v T d v. (4.4). where d v D cos.v /E 1 C sin.v /E 2 . An isotropic macroscopic friction implies a constant k independent of v D k whereas k varying with v ¤ k implies anisotropy. Instead of generating textures that display macroscopic anisotropy, the surfaces will be subjected to macroscopic surfacial deformation, and the induced anisotropy will be monitored to highlight finite deformation effects. 4.3. Two-dimensional contact interfaces 4.3.1. Unilateral roughness with periodicity. In order to demonstrate the fundamental features of the macroscopic response emanating from rough elastic boundary layers in frictional contact, twodimensional test cases are considered where periodic roughness will first be assigned to the slave (upper) surface only. Initially, the material parameters are chosen as .I / D 5 MPa and .I / = .I / D 2. The undeformed width (Lo ) and height (Ho , with respect to the mean plane of roughness for the rough side) of the samples are equal and set to 10 m while the RMS roughness is set to o D 1 m. The microscopic coefficient of friction is a constant at ko D 0:1. Contact at p N D 1 MPa is initiated in five normal steps after which the upper sample is displaced to the right in 20 tangential steps through a distance Lo . The samples are initially coarsely discretized, with five quadratic NURBS elements in the horizontal (H) and three in the vertical (V) direction. Figure 13 summarizes a series of four investigations involving the variation of the initial choices for the simulation parameters. First, the basic parameters are varied. It is observed that a cubic discretization alleviates the strong oscillations observed with quadratic elements at coarse discretizations. To obtain more accurate results, the number of cubic elements along each direction is additionally tripled, and the number of normal/tangential time steps is doubled. Second, on the basis of these choices, the mismatch between the slave and master elastic properties is varied at a constant .I / = .I / D 2 ratio. Depending on whether the properties of the rough or the smooth surface is kept fixed, k either approaches the limit ko that represents contact with a smooth rigid surface (mismatch > 1) or it reflects the indentation of a smooth soft surface with a rigid rough one Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(18) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 969. Figure 13. Simulation parameter effects are summarized in a two-dimensional setting with unilateral roughness.

(19) to represents the time step size for 5/20 normal/tangential steps. Beyond slave/master mismatch assessment, .2/ = .1/ D 1=8 is employed, and beyond compressibility assessment, .I / = .I / D 12 is employed.. (mismatch < 1). Clearly, in the large deformation regime considered, it does matter whether the soft or the stiff surface is rough, unlike the equivalent roughness approach that is typically suitable in the small deformation regime of contact mechanics [50]. In subsequent investigations, the mismatch is set to .2/ = .1/ D 1=8 such that the rough surface is slightly stiffer but far from being rigid. So far, the observed effect of roughness on the macroscopic response is small. In the third case, because biological or synthetic soft materials are nearly incompressible, the .I / = .I / ratio, which represents the compressibility of the samples, is varied and a strong influence on k is observed. Although the maximum value considered is not close to the incompressibility limit, possible locking effects are nevertheless checked for. For this purpose, quadratic discretizations are used on the basis of the Q2P0 mixed formulation. The Q2P0 approach does not deliver optimal convergence but successfully avoids locking [101]. Higher-order isogeometric discretizations can be handled with similar mixed formulations [102]. In view of the satisfactory performance of the cubic discretization, .I / = .I / D 12 will be employed in all subsequent investigations without the need for a mixed approach. Finally, length scale effects are considered in the fourth case. If all dimensions of a sample are varied, the macroscopic frictional response remains invariant. This is the case that motivates a contact homogenization approach based on the classical scale separation assumption. On the other hand, if only the width of a sample is scaled, the spectral properties of the rough surface changes while the spatial ones (e.g., o ) remain fixed such that k strongly varies. It is highlighted that for most surfaces, there is no separation between the individual scales of roughness. However, what is essential to contact homogenization is that there is a significant gap between the largest roughness scale and the representative dimension associated with the macrostructure. It is often reasonable to Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(20) 970. ˙I. TEMIZER. Figure 14. Studies complementing Figure 13 are summarized. For the height variation, below H=Ho D 0:6, only the lower sample height is reduced. For the pressure variation, the mesh resolution is doubled to 30/18 (H/V) cubic elements below p N D 0:5 MPa to ensure converged results and further increased to 75/45 (H/V) cubic elements for p N < 0:05 MPa.. assume the existence of such a gap for practical applications, and such an assumption is implicit in a significant portion of works on multiscale contact. It is this assumption that is made use of in this work, not one regarding the scales of roughness. Figure 14 summarizes a complementary set of four investigations. First, the strong effect of scaling only the spatial surface properties is observed. In subsequent investigations, o D 1:5 m will be employed. Second, it is observed that the height of the samples will affect the macroscopic response unless they are chosen sufficiently large. Presently, Ho D 6 m is found sufficient. The microscopic friction is varied from the default 0.1 in a third case. For vanishing ko , the effect of roughness approaches a supremum, and the effect significantly decreases with increasing friction. On the other hand, the fourth case shows that the decreasing effect can be offset by increasing the macroscopic contact pressure. In all remaining investigations, p N D 1 MPa will be employed. However, with decreasing pressure, k rapidly approaches ko . This final observation suggests that the observed deviation of k from ko is essentially due to the presence of large deformations and would vanish in a contact regime where linear elasticity is applicable. 4.3.2. Randomness and bilateral roughness. The case with random surfaces shows similar trends. Only the effect of randomness will be emphasized. The simulation parameters chosen from the preceding discussion are employed, except for the RMS roughness, which is chosen as o D 0:75 m. Presently, the default width Lo D 10 m will control the sample size: larger random samples will be generated by choosing the width to be multiples of Lo . For any given size, different realizations of the samples are possible because of randomness. This randomness reflects onto the frictional response, as summarized in Figure 15. In particular, at small sample sizes, there is a significant scatter in the macroscopic frictional response although this effect is spurious and can be alleviated by choosing larger samples that more closely represent the expected value of k. This trend is Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(21) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 971. Figure 15. Unilateral randomness effects are summarized. Ten samples were tested per sample size. The averages are shown with the standard deviation, slightly offset along the x-axis for clarity. As in Figure 14, roughness p is assigned only to the upper sample. T represents the Cauchy stress (MPa), and the magnitude kT k D T T of which is shown simply as an average over each isogeometic element only to highlight randomness.. Figure 16. The macroscopic frictional response for bilateral roughness with periodicity (o.1/ D 1:5 m, o.2/ D 0:2 m). Increasing friction hinders relaxation, leading to reduced oscillations. The dashed lines indicate the instantaneous response (kI =ko ) obtained with an initial phase shift imposed on the upper surface.. reflected by the converging ensemble (arithmetic) average of individual sample responses from a given sample size. Periodic or random, the instantaneous macroscopic friction coefficient is a constant with unilateral roughness once macroscopic sliding is initiated for sufficiently fine meshes (Figure 13). This behavior is significantly altered in the presence of even a small roughness on the other surface (Figure 16). The RMS roughness o.1/ on the upper surface will be retained from earlier investigations, and the lower one is assigned o.2/ D o.1/ =7:5. Despite this small roughness, significant oscillations in the instantaneous macroscopic coefficient of friction is observed, one period of oscillation corresponding to the period of lateral transition of the upper body over the lower one. One possible approach to identifying a unique k in this case is to compute a time average (Section 4.2). Another approach is to simulate various realizations of periodicity with initial phase shifts. While not explicitly calculated, it is observed that the average of these different responses at a given time will capture the time average obtained from any of the curves. The attractive feature of this observation is that the computation of the time average on large samples is not parallelizable while averaging among phase shifts is amenable to parallelization. In all cases, the source of the strong oscillations is interface Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(22) 972. ˙I. TEMIZER. Figure 17. For random surfaces with bilateral roughness (o .1/ D 0:75 m, o.2/ D 0:1 m), time averaging, effect of different realizations for a given sample size and the effect of increasing sample size are summarized. Except for time averaging, the samples are dragged through a distance 1:5Lo , which is sufficient to traverse one period of oscillations. Larger samples display milder oscillations.. relaxation where the lower surface slides back to relieve the shear stresses. Consequently, increasing microscopic friction arrests this relaxation so that the magnitude of the oscillations with respect to ko decreases. In the case of random bilateral roughness (Figure 17), significant oscillations can also be observed, and there is still an inherent periodicity in the oscillation pattern because of periodic embedding. Moreover, the magnitude and the frequency of the oscillations additionally significantly depend on the particular realization for a given sample size. Some realizations induce excessive shear loading due to conforming peaks and valleys while others show a lower degree of conformation. With increasing sample size, on the other hand, the influence of a single conforming peak–valley combination on the instantaneous macroscopic response becomes more isolated and less frequent such that the overall response looks smoother. This trend is already observed for various realization choices from different sample sizes. Because the sample sizes tested are still limited in size, different choices of realizations may not reveal as clear a trend. However, a comparison of individual sample responses at the smallest and largest sample sizes tested clearly demonstrates that larger sample sizes display significantly less scatter in the instantaneous macroscopic response. This leads to the interesting observation that although random surfaces are computationally more expensive to characterize, it may be possible to identify a time-independent instantaneous macroscopic coefficient of friction for a random interface topography with bilateral roughness while such an identification is not possible for the case of periodicity. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(23) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 973. 4.4. Three-dimensional contact interfaces Two-dimensional investigations satisfactorily demonstrate the computational framework and the nature of the macroscopic frictional response with unilateral/bilateral periodic/random roughness on two elastic micromechanical samples. The three-dimensional investigations of this section concentrate on the remaining aspect associated with macroscopic anisotropy. 4.4.1. Deformation-induced anisotropy. Macroscopic anisotropy can be induced by the imposition of macroscopic surfacial deformation F c . For this purpose, it is sufficient to assign a unilateral periodic roughness where the final set of parameters chosen earlier in the two-dimensional investigations will be employed. The simulation instances in Figure 18 for this setting display the mesh resolution, geometry, and deformation of the samples for the case when F c D I. Two types of surfacial deformation are considered such macroscopic deformation is area P that the ˛ E ˝ E ˛ C ˝ with 1 2 D 1, preserving (det F c D 1): (i) axial, where F c D ˛ ˛ and (ii) shear, where F c D I C E 1 ˝ E 2 . In both cases, the direction v of macroscopic slip is varied in increments of 30 degrees from 0 to 330. The corresponding k is summarized in. Figure 18. Simulation instances from the investigations of Section 4.4.1 display the resolution, geometry, and deformation of the samples. Only the central slave/master mesh is shown.. Figure 19. The anisotropic macroscopic frictional response (k=ko ) under axial-type macroscopic surfacial deformation is summarized. The first case employs 1 D –1 2 D 1:2, whereas the second one employs 2 D –1 1 D 1:2. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(24) 974. ˙I. TEMIZER. Figure 20. The anisotropic macroscopic frictional response (k=ko ) under shear-type macroscopic surfacial deformation is summarized. The first case employs D 0:2, whereas the second one employs D 0:4.. the polar plots of figures 19 and 20 where the angular position indicates the direction k of the frictional resistance (Section 4.2). The macroscopic response of the undeformed samples significantly differs from the microscopic friction coefficient but is isotropic. Axial and shear deformations induce macroscopic anisotropy such that friction decreases along the direction of the principal stretch exceeding unity because of increasing wavelength of periodicity and increases along the other principal direction because of decreasing wavelength—(cf. Figure 13(d)). Non-matching slip and frictional resistance directions are also reflected by v ¤ k . Clearly, larger deformations induce stronger anisotropy. A qualitatively similar situation is encountered in volumetric homogenization. A periodic microstructure with isotropic phases displaying cubic symmetry is sta-. Figure 21. The macroscopic response (k=ko ) for the case of unilateral random roughness without macroscopic surfacial deformation for the samples of Figure 22. Varying the mismatch ratio scales the response, which is observed to display anisotropy because of the small sample size. Moreover, for different realizations at the same sample size, the direction and magnitude of anisotropy also changes. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(25) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 975. Figure 22. Sample random surfaces generated and simulation instances from computations where they have been used to model unilateral roughness. See Figure 21 for the macroscopic responses.. tistically anisotropic but displays macroscopic isotropy with respect to conduction ([85], p.320). However, deformation destroys this isotropy [103]. Presently, statistical anisotropy for surfaces may or may not lead to a macroscopically anisotropic response. Deformation acts as a parameter that continuously varies the degree of macroscopic anisotropy by altering the statistical properties of rough surfaces. 4.4.2. Size-induced anisotropy. Macroscopic anisotropy may also be observed because of the choice of the sample size in the case of randomness. In order to demonstrate the capability of the numerical framework to address such sample size effects, unilateral random roughness is considered. Again, the final set of parameters chosen earlier in the two-dimensional investigations will be employed. The macroscopic frictional response is monitored for different mismatch ratios and for different realizations (Figure 21). Because of the statistical properties of the generated random surfaces (Section 3.1.1), the corresponding macroscopic response is expected to be isotropic for large sample sizes. Consequently, the observed anisotropy should be attributed to the small sample size employed (Figure 22). Analysis of the macroscopic frictional response for larger sample sizes and with more realistic random roughness descriptions [10] is left for a future study. 5. CONCLUSION The goal of this work was the establishment of a computational contact homogenization framework for the modeling and simulation of soft matter friction. The crucial ingredients toward the realization of this goal were (i) the establishment of a frictional contact algorithm that displays an optimal combination of accuracy, efficiency, and robustness and plays a central role in (ii) the construction of a micromechanical contact test within which samples of any size may be embedded and which is not restricted to a single deformable body. The former was realized through the extension of mixed variational formulations of contact mechanics to a mortar-based isogeometric setting where the augmented Lagrangian approach serves as the constraint enforcement method. The reliable performance of the resulting algorithm was highlighted through its comparison with pure penalty and Uzawastaggered versions. The latter was realized with periodic embedding, which ensures C 1 -continuous geometric periodicity across computational cells in addition to capturing the contact interactions among distant cells through cell periodicity. The overall homogenization framework was demonstrated by characterizing the influence of microscopic roughness on the macroscopic friction coefficient for elastic boundary layers. Twodimensional and three-dimensional samples with unilateral and bilateral roughness were considered on the basis of the variation of material and interface properties, highlighting observations that Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(26) ˙I. TEMIZER. 976. are unique to finite deformations. As a particular case of interest, the numerical results indicated that it may be possible to assign a well-defined macroscopic friction coefficient to a random interface topography with bilateral roughness while this is not possible in the case of periodicity. This observation is in contrast with the mechanics of random materials where periodicity facilitates the identification of the homogenized response but is in convenient agreement with the engineering demand for a uniquely identifiable friction coefficient. Periodic embedding enables achieving arbitrarily long dragging durations for two rough surfaces in sliding contact independent of the micromechanical sample size. This is advantageous for the macroscopic thermomechanical characterization of micro-rough soft contact interfaces, which entails the homogenization analysis of both the friction coefficient and the thermal contact parameters associated with frictional heating. The presented framework also allows the incorporation of microscopic lubrication and adhesion mechanisms in a numerical efficient manner, thereby constituting a basis toward a comprehensive computational tribology characterization of such interfaces. APPENDIX A: GENETIC ALGORITHM FOR CLOSEST-POINT PROJECTION A simple genetic search algorithm [104] is applied in order to approach the global minimum before switching to a gradient search within the closest-point projection algorithm of Section 3.2.2: (1) Initialize gene pool: A pool of N random points with parametric positions i are generated in the one-cell neighborhood of the last projection cell: i˛ 2 Œ ˛;n 1; ˛;n C 2, i D 1; : : : ; N . (2) Fitness evaluation: The fitness of each point is evaluated using the projection distance as the cost function: Fi D kx y i k with y i as the physical position vector of point i. (3) Create new generation: A new generation of guesses are generated: (a) Selection: i are sorted according to their fitness such that i D 1 corresponds to the point with the minimum distance. Subsequently, the best gene is retained: new D 1. 1 (b) Crossover: Top n < N genes generate offspring, for example, new D ˛ i i 1 C .1 i ˛i / i where ˛i 2 Œ0; 1 is a random number for i D 2; : : : ; n. (c) Resupply: Remaining genes are replaced by randomly generating new for i D n C 1; N . i (4) Convergence criterion: the genetic search is carried out for G generations or until the best gene 1 does not significantly change among subsequent generations, whichever is met earlier. Alternatively, one can adapt contact search algorithms [105] to approach the global minimum. APPENDIX B: VISUALIZATION OF PERIODIC EMBEDDING While the slave mesh continuously moves away from the central master mesh and interacts with it only through the images according to the algorithm of Section 3.2.2, it is advantageous to visualize the slave mesh in the one-cell neighborhood of the central master. Denoting the dragging duration with t , the slave control point positions with x I and the macroscopic slip velocity with vT , the algorithm for this purpose sequentially checks ˛ D 1 to 2 (no sum over ˛; bc is the floor function): (1) Basis change: using m˛ˇ D e ˛ e ˇ and with m˛ˇ as the components of its inverse, let p e ˛ D m˛ˇ e ˇ . Then, v ˛T D m˛˛ vT e ˛ are the components of vT with respect to e ˛ =ke ˛ k. (2) Cell location: let ˛ D 2b¹v ˛T t = kLo e ˛ k C 1º=2c and ˛ D ¹v ˛T t = kLo e ˛ k C 1º ˛ 2 Œ0; 2/. (3) Preliminary shift: the control points are shifted back to the central cell: x Io D x I vT t . (4) Shift adjustment: the preliminary image is adjusted for each tangential direction to obtain the final visualization positions: x Ishift D x Io C Lo .˛ 1/e ˛ . The one-cell neighborhood images of the slave, and similarly those of the master, are generated from the shifted central mesh via x Iimage D x Ishift C Lo ˛ e ˛ where ˛ D ˙1. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.

(27) COMPUTATIONAL HOMOGENIZATION OF SOFT MATTER FRICTION. 977. Figure C.1. Number of Uzawa iterations required to achieve the results in Figure 2(d) for N D 102 . Closer lines indicate time step size refinement. The average number of Newton–Raphson iterations per Uzawa iteration is 5/4/3 for T OL D 101 =102 =103 .. Figure C.2. Newton–Raphson iterations required for the AL scheme at different time instances: (a) iterations for Figure 2(b) and (b) iterations for Figure 4(b).. APPENDIX C: ITERATIVE CONVERGENCE PERFORMANCE Figure C.1 summarizes the Uzawa iterations required to achieve the results of Figure 2(d), demonstrating the high cost of the Uzawa-staggered form of the AL scheme for small Uzawa tolerances. The Newton–Raphson iterations for the AL scheme in Figures 2(b) and 4(b) are summarized in Figure C.2. Convergence is achieved when the error k

(28) uk2 =kuold k2 decreases below a tolerance of 107 . Here,

(29) u denotes the update to the vector u of control point positions and uold is the value of u from the last iteration. For both examples, asymptotically quadratic convergence rates are observed. Note that the chosen time instances involve active-to-inactive and stick-to-slip transitions. Hence, there appear to be no shortcomings that result from the simplifications introduced into the contact algorithm in Section 2.3 to handle these transitions. ACKNOWLEDGEMENT. Support for this work was provided by the Scientific and Technological Research Council of Turkey (TÜB˙ITAK) under the Career Programme Grant No. 110M661. The author would like to thank one reviewer for the meticulous reading. Copyright © 2014 John Wiley & Sons, Ltd.. Int. J. Numer. Meth. Engng 2014; 100:953–981 DOI: 10.1002/nme.