Mesoscale modelling of polymer-mediated adhesion: application to tack tests

Abstract

We develop a generic computational methodology to understand and predict adhesion between polymers and solid substrates. The motion of coarse-grained polymer segments is tracked via a hybrid particle-field mesoscopic simulation method (BD/kMC) combining Brownian dynamics (BD) and kinetic Monte Carlo (kMC) for the entanglement dynamics as described by the slip-spring model. The method addresses entangled polymer films capped between solid surfaces under both quiescent and nonequilibrium conditions. The latter entail imposing constant rate extension along the aperiodic (normal) direction, while keeping the lateral dimensions constant. Experimentally relevant length scales and elongation rates can be addressed thanks to the coarse-graining inherent in the approach. These simulations are representative of “tack” tests, employed routinely for assessing the performance of soft adhesive materials. The performance of each interface is characterized by the stress–strain curves, yield stress, and toughness. The failure mechanism is determined upon analyzing the evolution of the stress–strain curve and the morphology of the fractured interfaces. The simulations are conducted over a broad parameter space by varying the rate of elongation, the rate constants for attachment/detachment of polymer segments to/from the surface, and the activation length. The latter describes the coupling with the pulling forces exerted on the particles at the interface by the rest of the polymer. Setting the activation length to zero is suitable for describing strong adhesives or highly compressible materials (foams). Under these conditions, toughness is maximized and increases significantly with elongation rate, sometimes leading to chain fracture. With increasing activation length the toughness of the interface decreases and detachment becomes more efficient at higher elongation rates since the increased stress accelerates the detachment process. In all cases considered here, toughness increases monotonically with adhesion. Furthermore, the yield stress increases consistently with increasing elongation rate due to the inability of the polymer to relax the imposed stress.

---

Design, System, Application

Designing adhesives with tailor-made properties for specific applications constitutes a complex challenge involving chemistry, physical chemistry, and mechanics. We develop a mesoscopic computational strategy for quantifying the adhesion of long polymer chains on solid surfaces, with relaxation times spanning the millisecond to second time regime. The short-time dynamics of the coarse-grained moieties evolves with Brownian dynamics. The entanglement dynamics of the polymer chains is described with the Slip-Spring model. Computational specimens are subjected to constant rate elongation experiments, mimicking the setup of experimental tack tests. The performance of the interfaces is quantified by the evolution of stress-strain curves, toughness, and the morphology of the fractured surface. Toughness depends on a complex interplay of elongation rate, attachment/detachment rate constants, and the coupling of the pulling force with the rate constants. High elongation rates lead to increased yield stress, enhanced toughness, and pronounced fibrillation/cavitation, especially for strongly adhering polymers. In addition, the polymer chains detach more efficiently (at lower strains) with increasing elongation rate, especially as the coupling between the detachment rate and the pulling force improves. Our computational strategy is versatile, generic and can be applied to a range of industrial polymers and solid surfaces.

1. Introduction

The optimal design of adhesives and adhesive joints for specific materials is a very challenging problem involving chemistry, physical chemistry, and mechanics.^1–6 Apart from chemical bonding and physical interactions at the adhesive/substrate interface, the topography of the substrate surface and the elastic, viscoelastic, and ultimate properties of adhesive and substrate play a significant role. Tackiness is the characteristic property of pressure-sensitive adhesives (PSA) enabling them to adhere to surfaces upon simple contact, without requiring heat or chemical reactions.⁵ In other words, low contact pressure and short contact time are sufficient for bond formation with sufficient adhesive strength, referred to as ‘tack’ (ref. 4 cited in ref. 7). The higher the work of detachment per unit area the higher the tackiness of the PSA.⁷

Tack experiments are employed routinely for the characterization of the adhesive properties and toughness of polymer/solid interfaces.^5,8–11 A weight, with a flat or spherical tip, is brought into contact with the adhesive film and kept in contact for a given time. Subsequently, the weight is lifted at a steady rate, thus subjecting the interface to tensile deformation. By normalizing⁹ the recorded force and displacement throughout the course of the experiment one can obtain the tensile stress (σ) as a function of the engineering strain (ε). Accompanying the tack experiments with real-time in situ optical observations allows determining the failure mechanism and the form of the fractured interface.⁸ The mode of failure depends on interplay between the cohesive interaction of the polymer and the adhesive strength of the polymer/solid interface. Good adhesion promotes cohesive failure and rough end surfaces with bridging fibrils (fibrillation). Weak adhesion, on the other hand, favors adhesive failure; i.e., the polymer detaches completely from the solid surface. The pulling rate (or debonding rate)⁵ affects significantly the fracture mechanism as well; in general, fast (slow) pulling rates promote adhesive (cohesive) failure.^5,9

Modeling adhesion has been dominated by continuum fracture mechanics analyses and numerical simulations (e.g., by the finite element method), using criteria based on the local stress state to describe crack propagation.^12,13 Atomistic molecular dynamics simulations have been invoked to understand crack propagation and cohesive failure in two-dimensional materials such as graphene,¹⁴ brittle solids such as NaCl,¹⁵ a silica crystal (α-cristobalite) and a nanoporous, highly crosslinked polymer (saccharose-based heat-activated carbon).¹⁶ The latter study is particularly thorough in that it compares loading curves from molecular dynamics against the corresponding continuum mechanics predictions. It clearly shows that linear fracture mechanics works well in the case of the silica crystal, a brittle solid, but fails for the porous polymer, where the process zone is an order of magnitude larger. The introduction of a cohesive zone to describe plastic deformation past the crack tip in the continuum analysis reconciles the continuum description with the atomistic molecular dynamics.¹⁶

Failure at bi-material interfaces has been simulated with atomistic molecular dynamics for various systems. An interface between asphalt and crystalline silica was subjected to tension normal to the interface by Xu and Wang¹⁷ under elongation rates on the order of 10⁹ s⁻¹; it was shown that fast/slow pulling rates promote adhesive/cohesive failure in accordance with experimental tack tests.^5,9 An effort to simulate pulling, peeling, and sliding in polyimide/silica glass systems using steered molecular dynamics and model systems of finite (nanoscopic) size was presented by Min et al.¹⁸ The molecular dynamics simulations of Gersappe and Robbins¹⁹ demonstrated that failure occurs in the region with the lowest initial yield stress and cannot be solely predicted by equilibrium interfacial free energies. As a result, failure can arise either at the polymer/solid interface (adhesive) or within the polymer film (cohesive).

In a very recent set of atomistic nonequilibrium molecular dynamics simulations, Kallivokas et al.²⁰ predicted the conditions of local failure for graphene/epoxy interfaces subjected to shear stress over a range of temperatures, using periodic computational specimens. They found that failure times exhibited an extended Boltzmann–Arrhenius–Zhurkov dependence on stress and temperature. They also determined how the activation energy and activation volume governing the kinetics of failure changed upon introduction of defects or polar functional groups on the graphite surface. A similar study was conducted on carbon nanotube/epoxy interfaces.²¹

In general, atomistic molecular dynamics simulations of adhesion are limited by the high elongation rates and the small model system sizes they employ. When polymeric adhesives are present, atomistic simulations are typically confined to short chain lengths. Results are strongly rate dependent and cannot capture the extended plastic deformation zones associated with crack propagation in real polymeric materials.

Mesoscopic simulations offer the opportunity to address larger length and time scales. An important early work was that of Stevens.²² He undertook molecular dynamics simulations of interfacial fracture between highly cross-linked polymer networks, intended to mimic epoxy resins, and a solid surface, using a Kremer–Grest bead and spring model for the chains. Some of the wall particles were assumed to bond to polymer beads by breakable bonds; they were arranged in regular patches on the surface. The surface fraction occupied by these patches was varied systematically and its impact on the stress–strain curve and failure mode was examined. More recently, mesoscopic molecular dynamics of cohesive-adhesive failure of melts with Kremer–Grest chains slightly above the entanglement molar mass sandwiched between two “chemically heterogeneous” surfaces have been presented by Baggioli et al.²³ The model surfaces consisted of strongly and weakly adsorbing sites in 1 : 1 proportion, but were characterized by different degrees of segregation of these sites. When strongly and weakly adsorbing sites were randomly interspersed on the surface, the surface behaved as a homogeneous one and adhesion was strong. As the distribution of strongly (and weakly) adsorbing sites became more “patchy”, characterized by larger length scales, adhesion weakened. In addition, Baljon et al.²⁴ demonstrated that adhesive energy increases with prolonged contact duration between the polymeric material and the surface due to the redistribution of the polymer segments across the surface.

A set of mesoscopic entanglement network-based simulations explored adhesion between polypropylene and polyamide 6 in the presence of polypropylene chains terminally grafted to the polyamide surface.²⁵ It was determined that there is an optimal areal density of grafted chains, for which the fracture energy is maximized and which can play a role when barrier properties or penetration issues are regarded for light elements.²⁶

The effect of varying the crosslinking density has been investigated by Jin et al.²⁷ and Solar et al.¹¹ In the absence of crosslinks the interface fails cohesively, whereas increasing the crosslinking degree enhances the tendency for adhesive failure. The performance of the interface was maximized for intermediate crosslinking degrees; at high degrees of crosslinking, even though the plateau stress increased, the sample became more brittle and failed at lower strains. The shape of the stress–strain curve at small strains was not affected by the crosslink density.

Subjecting highly entangled polymer chains to shear^28,29 and extensional^30,31 flows can result in bimodal chain end-to-end distance distributions under specific conditions. When specimens are subjected to extensional flow, they can undergo a flow-induced phase separation, resulting in the formation of distinct domains: one dominated by coiled macromolecules and another by significantly stretched, chain-like structures.^30,31

A recent work by some of the authors³² proposed a multiscale simulation strategy for investigating and predicting interfacial structure and dynamics at entangled molten polymer/solid interfaces over a wide range of time and length scales, using polyethylene/graphite as an example. At the mesoscopic level, Brownian dynamics was used to follow the motion of beads in three-dimensional space while, in parallel, kinetic Monte Carlo simulation tracked slippage across entanglements, generation and loss of entanglements at chain ends, and adsorption/desorption events of coarse-grained segments on the solid substrate. The multi-chain slip-spring model was employed to represent entanglements,^33–45 imposing constraints on the lateral motion of chains relative to their contour. The slip-spring model, a tube theory-based model, has demonstrated remarkable success in reproducing the experimental viscoelastic properties of polymers.^{29,34,46–52} The mesoscopic polyethylene/graphite model, extending the engine for mesoscopic simulations of polymer networks (EMSiPoN)^29,34,53 to interfaces, was parameterized entirely on the basis of atomistic molecular dynamics simulations of shorter-chain melts at the same interface.^32,54 The residence time of chains in the adsorbed state and their relaxation time under equilibrium conditions were determined as functions of molar mass, in good agreement with experiments.

The performance of the mesoscopic model has been tested extensively over a wide range of system geometries (e.g., bulk conditions,^29,34 free films⁵⁵ and polymer-solid³² interfaces), chain architectures,^51–53 and flow conditions.^29,52 The variation of the viscoelastic properties with increasing chain molecular weight conforms with scaling theories and data from atomistic simulations and experiments.^29,56–58 Under high shear rate deformations, the model reproduces the startup stress overshoot observed in both simulations^59,60 and experiments.^61,62 Furthermore, it exhibits shear viscosities that scale with Weissenberg number similarly to those observed in slip-link simulations⁶³ and experiments.^62,64 The transition of chain end-to-end distance distributions from unimodal to broad bimodal,²⁹ indicative of chain extension followed by rapid retraction, aligns with observations from detailed atomistic simulations.²⁸

Single-chain slip-link^65–67 and slip-spring^33,43 models describe accurately the long-time dynamics and relaxation phenomena in highly entangled polymers. Despite being computationally efficient,⁶⁷ single-chain models exhibit limitations, such as neglecting the spatial correlations between chains, the spatial inhomogeneities at the interfaces, and other effects.^40,43 On the contrary, multi-chain slip-link⁴⁵ and slip-spring^{29,34,38–40,63} models incorporate three-dimensional correlations. Being more detailed, multi-chain models can be connected more readily to atomistic calculations, and applied straightforwardly to a variety of situations, including free polymer surfaces,⁵⁵ branched polymers,⁶⁸ nanocomposites,³² etc. According to Masubuchi and Uneyama,⁶⁹ the chain conformations, diffusion and relaxation modulus are compatible among multichain slip-link and slip-spring models; there are, however, significant discrepancies in inter-chain cross-correlations within the relaxation modulus. A key advantage of slip-spring models over slip-link models is their foundation in a well-defined free-energy description, which facilitates the incorporation of interchain interactions.⁴³

In this work we develop a mesoscopic simulation framework to model the constant-rate elongation of polymer films confined between adhesive solid surfaces. Throughout the simulation the lateral dimensions of the system are kept constant. The aforementioned setup is similar to the configuration employed in experimental tack tests.⁵ A step-by-step procedure is provided for parameterizing the mesoscopic model in a bottom-up or top-down fashion based on observables from atomistic simulations or experiments, respectively. For demonstration purposes, we parameterize the mesoscopic model using a minimal set of data representative of a noncrystallizable fluoropolymer. The response of the stress–strain curves and the corresponding failure mode are explored across a broad parameter space in terms of varying the attachment/detachment rate constants, the elongation rate and the coupling of the interfacial kinetics with the pulling force exerted by the polymer segments within the polymer film.

The article is structured as follows: methodological details are provided in section 2. In particular, section 2.1 develops the mesoscopic model for adhesive surfaces, section 2.2 provides a bottom-up parameterization framework for mapping the atomistic observables to the mesoscopic model, and section 2.3 illustrates the equilibration protocol preceding the tack tests. Section 3 discusses the findings from the tack test simulations. Section 3.1 illustrates the setup of tack simulations. Section 3.2 discusses results from the model with constant-rate interfacial kinetics. Section 3.3 introduces the effect of coupling the attachment/detachment rates with the pulling forces exerted by the polymer within the film. Finally, section 4 concludes the manuscript.

2. Methods

2.1 Mesoscopic model

Each bead describes a polymer subchain of several Kuhn segments. Neighboring beads in the same chain are connected with entropic springs (hereafter referred to as strands) modelling the entropic elasticity of the chains, as shown in Fig. 1. Different chains can be connected with permanent (crosslinks)⁵³ or temporary springs; the latter will be referred to as “slip-springs” (SS) and mimic the effect of chain entanglement.

Fig. 1 Schematic illustration of the mesoscopic interface.

The beads can be attached to/detached from a solid surface subject to the transition rates:

with a/f denoting the attached/free state, E_i→j a barrier energy between states i and j and v an oscillation frequency. l_i→j is a characteristic length scale (activation length) tuning the coupling of the transition rate with the normal force F_z toward the bulk phase, exerted on the bead by the rest of the polymer. F_z is averaged over a short time interval (details in section 3.3) to reduce high frequency fluctuations. is the transition rate in the absence of normal force coupling. The transitions take place within an interfacial domain of thickness L_d. The attached beads are fixed in relation to the nearest solid surface for as long as they reside in the attached state.

The thermodynamics is described by a Helmholtz free energy functional of the form:^29,32

A({r_ij}, {ρ(r)}, T) = A_b({r_ij}, T) + A_SS({r_ij}, T) + A^ex_nb({ρ(r)}, T) + A_wall{r_ij}. (3)

A_b in eqn (3) denotes the contribution of the strands,
anddenotes the contribution of the n_ss slip-springs (SS). The free energy of the strands/slip-springs is the following:rendering eqn (5) continuous and differentiable at x = x_{tol,strand/SS}, with

being Kröger's approximation of the inverse Langevin function,^29,70 λ_strand/SS the slope of the force at small extensions, and x = r_i/r_{max,strand/SS} the extension ratio, where r_{max,strand/SS} and r_i are the maximum extension and length of the strand/SS, respectively. A breaking criterion was incorporated to simulate chain scission, triggered when a strand's length exceeded its contour length, L_{contour,strand}; see section 2.2.3.

A^ex_nb is the contribution from the nonbonded interactions,

It is expressed in terms of an excess free energy density f computed from the local density distribution ρ(r). The excess free energy density f can be retrieved from an equation of state (EOS). Here we will invoke the Sanchez–Lacombe EOS:⁷¹cast in terms of the reduced quantities [small rho, Greek, tilde] = ρ/ρ*, P̃ = P/P*, T̃ = T/T*, and r_SL = P*M_bead/(ρ*RT*), with ρ*, P* and T* being the characteristic mass density, pressure and temperature of the equation of state, and M_bead the molar mass of the bead. The rightmost term in eqn (6) multiplied by P* corresponds to the ideal gas contribution; P*T̃[small rho, Greek, tilde]/r_SL = n_beadk_BT, with

n_bead = ρ/m_bead,

being the bead number density, m_bead = M_bead/N_Avo, and N_Avo Avogadro's number. Although only the excess part contributes to conservative forces, the ideal part must be included in calculating the total stress in position Langevin dynamics where there are no velocities.⁷²

The spatial discretization of the domain is realized via a kernel-based discretization scheme,⁷² where the local bead density (n_bead,i = ρ_i/m_bead,i) is assigned to individual particles through a weighting kernel with finite support, r_c:⁷²

withbeing Lucy's function.⁷³

Finally, A_wall is the contribution from the polymer-wall interactions:

A_wall = N_aw_a + N_fw_f (9)

where N_a and N_f is the number of adsorbed and free beads within an interfacial region, and w_a, w_f are the free energies per bead in each state. The latter can be related as follows:³²

w_a − w_f = −k_BT ln(K_ads) (10)

where

K_ads = k_f→a/k_a→f (11)

is the equilibrium constant for attachment (adsorption), and k_f→a and k_a→f are the rate constants for attachment and detachment, respectively.

The concentration of the attached (C_a) and detached (C_f = 1 − C_a) beads at the interface is dictated by the microscopic reversibility condition, C_ak_a→f = C_fk_f→a, or equivalently, by the equilibrium constant for attachment,

Although w_a and w_f do not affect the forces (and consequently, the dynamics), they contribute indirectly to the adhesion tension of the interface. Supposing that the solid surface is incompressible, the adhesion tension can be defined as the difference between the surface free energy of the solid substrate (γ_sv) and the interfacial energy (γ_sf):⁷⁴

γ_adh ≡ (γ_sv − γ_sf). (13)

The adhesion tension of the mesoscopic interface can be expressed using a modified^32,75 Kirkwood–Buff relation of the form:with S = 2L_xL_y being the total area of the interface of a capped polymer film,⁵⁴ V the system volume, and 〈τ_T〉 and 〈τ_N〉 the average tangential and normal stress, respectively. By fixing γ_adh to a target value (e.g., from experiments or atomistic simulations⁵⁴) and estimating 〈τ_T〉, 〈τ_N〉, N_a and N_f from an equilibrium simulation, one can determine the free energies per bead in each state upon solving a system of two equations (eqn (10) and (14)) with two unknowns (w_a and w_f). For more details the reader is referred to section 5.5.2 in ref. 32.

The mesoscopic particles execute Langevin dynamics in the high friction limit, subject to the stochastic equation of motion:^29,32

where F_i = −∇_riA is the total conservative force exerted upon each bead. ξ_i is a stochastic force with zero mean obeying the fluctuation–dissipation theorem, and ζ_bead,i is the friction factor per bead.

The dynamics of the slip-springs are dictated by a microscopically reversible kinetic Monte Carlo scheme (hopping scheme)^29,34 for non-constant number of slip-springs (ensemble is grand-canonical with respect to the slip-springs) subject to an activity z_activ = e^μ/k_BT, where μ is the corresponding SS chemical potential. The hopping probability of a slip-spring end to an adjacent bead along the chain on which it lies is calculated according to the following equation:³⁴

with v_hop being a frequency factor, A_SS the instantaneous free energy of the SS, and Δt_kMC the elapsed time between hopping attempts. New slip-springs can be created at the chain ends subject to the probability

P_form,i = v_hopz_activn_cands,iΔt_kMC (17)

where n_cands,i is the number of candidate beads within an attempt radius (r_attempt) that are capable of forming slip-springs with terminal bead i. For additional details, the interested reader is referred to previous publications by the authors.^29,32,34

2.2 Bottom-up parameterization framework

We have developed a parameterization framework for mapping the mesoscopic model to observables from atomistic simulations and/or experimental measurements. As a case study, we will parameterize the mesoscopic model based on the minimal data set illustrated in Tables 1–3.

Table 1 Properties of the effective monomer

Mmon (Da)	50
lc–c (Å)	1.54
θc–c–c (°)	112

---

Table 2 Structural and dynamic properties of unentangled polymers at 300 K and 1 bar

Nmon/chain	Mchain (kDa)	Rete^2 [Å]	Nmon/chainD (10^12 m^2 s^−1)
50	2.5	619	24.863
100	5	1395	45.972
150	7.5	2113	44.264
200	10	2843	44.798

---

Table 3 Volumetric properties

Mchain (kDa)	P (bar)	T (K)	ρ (kg m^−3)
10	1	300	1500
10	100	300	1545
10	200	300	1557
7.5	1	300	1499
7.5	1	400	1425
7.5	1	500	1328

---

In particular, Table 1 displays the molar mass of the effective monomer M_mon (average molar mass per backbone carbon bond), the bond length l_c–c between successive backbone carbon atoms and the bond angle θ_c–c–c formed by triplets of backbone carbon atoms. Table 2 depicts indicative structural (mean squared end-to-end distance) and dynamical properties (diffusion coefficient) of the polymer across the Rouse regime. Table 3 illustrates the volumetric properties of the reference polymer as a function of molar mass, temperature and pressure. The dataset of the reference sample is indicative of a noncrystallizable fluoropolymer and was obtained via atomistic simulations. Table 4 presents key quantities extracted from the dataset in Tables 1–3.

Table 4 Properties of the reference system

Quantity	Value	Equation	Description
C∞	6	23	Flory's characteristic ratio
bk (Å)	11.15	24	Kuhn length
Me (kDa)	22.56	20	Molar mass between entanglements
atube (Å)	79.33	30	Tube diameter
ζmon (10^−12 kg s^−1)	92.13	29	Monomeric friction coefficient
ρ300 K (kg m^−3)	1504.1	6	Mass density (limit Mchain → ∞)
κex,300 K (MPa^−1)	0.00585	22	Excess isothermal compressibility

---

The mesoscopic model can be parameterized via a six-step procedure which has as follows:

1. Coarse graining

2. Volumetric properties

3. Conformational statistics

4. Bead dynamics

5. Slip-spring kinetics

6. Kernel optimization

For the convenience of the reader, the parameters of the mesoscopic model are illustrated in Table 5. The following sections provide a detailed description of each step.

Table 5 Parameters of the mesoscopic model

Scheme	Quantity	Value	Description
Coarse-graining	Nmon/bead	160	Monomers per bead
	Mbead (kDa)	8	Bead molar mass
Strand	λstrand (kcal mol^−1 Å^−2)	0.0031	Strand stiffness
	rmax,strand (Å)	204.27	Strand maximum extension
	xtol,strand	0.99	Strand extension tolerance
Slip-spring	λss (kcal mol^−1 Å^−2)	0.0004	Slip-spring stiffness
	rmax,SS (Å)	119	Slip-spring maximum extension
	xtol,SS	0.99	Slip-spring extension tolerance
Kinetic Monte Carlo	rattempt (Å)	117.8	Slip-spring formation radius
	vhop (s^−1)	65.21	Hop pre-exponential factor
	zactiv	0.002235365	Slip-spring activity
EOS	ρ* (kg m^−3)	1559.83	EoS characteristic density
	P* (MPa)	20.3	EoS characteristic pressure
	T* (K)	759.86	EoS characteristic temperature
Kernel	rc (Å)	50	Kernel support
	CW	0.953	Kernel weighting factor
	Cζ	0.900	Kernel friction factor
Brownian dynamics	ζbead (10^−9 kg s^−1)	14.74	Bead friction coefficient
	Δtcrit (ns)	58	Maximum time step

---

2.2.1 Coarse-graining. The first stage deals with setting the coarse graining degree of the mesoscopic model, that is, the molar mass ascribed to each bead, M_bead. Aggressive coarse graining allows for considerable simulation speedup. Care must be exercised, however, in order to avoid coarse-graining away length scales which might be relevant for the mesoscopic phenomena under study. For example, it is advisable that the bead molar mass be lower than the molar mass between entanglements:

M_bead < M_e (18)

so as to respect the picture of the tube model.^58,76 Indeed, the investigations by Masubuchi and Uneyama⁵⁸ have shown that, in most cases, the dynamics of Rouse chains interconnected by virtual springs are largely independent on the level of coarse-graining, provided appropriate rescaling is applied. Several strategies have been devised for obtaining M_e.⁷⁷ Here we estimate it from the packing length,as follows:

M_e = N_An_t²ρp³ (20)

where n_t has a universal value of 20.6 (±8%) across a range of polymers.⁷⁷

The packing length of the longest chain considered in Table 2 is p = 3.89 Å at 300 K, resulting in M_e = 22.56 kDa. The bead mass was set to M_bead = 8 kDa per bead, roughly three times smaller than M_e. This corresponds to N_mon/bead = M_bead/M_mon = 160 monomers per bead.

2.2.2 Volumetric properties. The characteristic coefficients of the equation of state⁷¹ can be determined through volumetric measurements (i.e., density as a function of temperature, pressure and chain length) either from experiment, or from atomistic simulations. Here, we opted to calculate them based on the reference volumetric measurements in Table 3. The characteristic EoS parameters were derived numerically in order to match the density from eqn (6) with the density from the reference sample over the whole (T, P, M_chain) range considered here. In particular, the optimization involves minimization of the following objective function:where ρ_EoS is determined by solving eqn (6). The optimization was performed with the Nelder–Mead method.^78,79 Fig. 2a illustrates comparisons between the reference densities and the evaluations with the equation of state (eqn (6)) using the optimized set of parameters. Fig. 2b shows the evolution of the merit function (eqn (21)) during the course of the optimization. The inverse excess compressibility of the polymer equals

κ_ex⁻¹ = n_bead(∂P_ex/∂n_bead)|_T (22)

with P_ex = P − n_beadk_BT being the excess pressure. At 300 K, ρ ≃ 1500 kg m⁻³ and eqn (22) yields κ_ex = 0.00585 MPa⁻¹.

Fig. 2 (a) Density versus pressure, temperature and chain molar mass from the reference sample (solid, Table 3), and the equation of state (diagonal pattern) with the optimized EoS parameters listed in Table 5. (b) Evolution of the merit function (eqn (21)) during the course of the optimization starting with the initial set of parameters, (ρ*,P*,T)_init = (1692.77 kg m⁻³, 43.6 MPa, 769.81 K).

2.2.3 Conformational statistics. Subsequently we parameterize the conformational free energy terms of the mesoscopic model based on conformational statistics of the reference system. The stiffness of a polymer chain can be quantified in terms of Flory's characteristic ratio, C_∞:⁸⁰with 〈R_ete²〉 being the mean squared end-to-end distance of a polymer chain consisting of N_mon/chain skeletal bonds.

A useful quantity for the mapping atomistic segments to the mesoscopic level, is the Kuhn length (b_k):⁸¹

where γ = sin(θ_c–c–c/2). The number of Kuhn segments per mesoscopic bead can be retrieved as follows:⁸¹The mean squared strand length equals:

〈R_{ete, strand}²〉 = N_k/beadb_k² (26)

and the contour length:

L_{contour,strand} = N_k/beadb_k (27)

The bonded terms of the free energy were parameterized based on the conformational statistics of Table 2. For clarity, the variation of the end-to-end vector with molar mass of the reference system is illustrated in Fig. 3. Flory's characteristic ratio was estimated as C_∞ = 6 for the largest chain considered here (eqn (23)). The corresponding Kuhn length equals b_k = 11.15 Å (eqn (24)) and the number of Kuhn segments per bead is N_k/bead = 18.33 (eqn (25)). Modelling finitely extensible strands, as in eqn (5), necessitates defining both their maximum extension (r_max,strand) and stiffness (λ_strand). The former was set equal to the contour length of the strand, r_max,strand = L_{contour,strand} = 204.27 Å (eqn (27)). The maximum extension tolerance was set to x_tol,strand = 0.99. The stiffness was optimized to λ_strand = 0.0031 kJ mol⁻¹ Å⁻²) upon matching the second moment of the strand length distribution with eqn (26), 〈r_strand²〉 ∼ 〈R_ete,strand²〉, using Brent's optimization algorithm.⁷⁹ As shown in Fig. 3, the mesoscopic model faithfully reproduces the conformations of the reference dataset and the corresponding random walk model, eqn (26).

Fig. 3 Mean square end-to-end vector of the reference polymer from Table 2 (circles), mesoscopic model (squares), and the random flight model, eqn (26) (dashed line). The mesoscopic simulations were conducted using the parameters in Table 5.

2.2.4 Bead dynamics. The bead friction coefficient (ζ_bead,i in eqn (15)) can be determined from the monomeric friction coefficient (ζ_mon) according to eqn (28).By mapping the atomistically computed dynamics onto the Rouse model (unentangled melts, M_chain < M_e), one can relate ζ_mon with the self-diffusion coefficient D according to eqn (29).⁸²

Fig. 4 illustrates the self-diffusion coefficient scaled by the number of monomers for each chain in Table 2. The friction coefficient was estimated via eqn (29) across the plateau region in Fig. 4 to ζ_mon = 92.13 × 10⁻¹² kg s⁻¹. As shown in Fig. 4, the mesoscopic model reproduces the short-time dynamics from the reference system and the corresponding Rouse model (eqn (29)). It is noted that the model exhibits limitations in accurately capturing short-time dynamics, particularly those associated with segmental motions of the polymer chains. This shortcoming is evident in the limited accuracy of predictions for G′ and G′′ at high frequencies.³⁴

Fig. 4 Self-diffusion coefficient scaled by the number of monomers per chain of the reference polymer from Table 2 (circles), the mesoscopic model (squares) and the Rouse model, eqn (29) (dashed line). The mesoscopic simulations were conducted using the parameters in Table 5 without slip-springs (since M_chain < M_e).

2.2.5 Slip-spring kinetics. The final stage deals with restoring the viscoelastic properties of the polymer. The slip-spring model requires determination of the: tube diameter (a_tube), slip-spring free energy (A_SS,i), slip-spring formation radius (r_attempt), number of slip-springs per bead (N_SS/bead), pre-exponential factor (v_hop) and slip-spring activity (z_activ).

The tube diameter can be estimated from M_e based on eqn (30):⁸³

Resulting in a_tube = 79.33 Å. The maximum extensibility of the slip-springs (eqn (5)) was set to r_max,SS = 1.5a_tube = 119.0 Å and the stiffness was optimized with Brent's optimization algorithm⁷⁹ to λ_SS = 0.0004 kJ mol⁻¹ Å⁻² by matching the second moment of the length distribution with the squared tube diameter (i.e., 〈r_SS²〉 ∼ a_tube²). The maximum extension tolerance was set to x_tol,SS = 0.99. Subsequently, the hopping attempt radius was set to r_attempt = x_tol,SSr_max,SS.

The number of slip-springs per bead can be estimated based on the number of entangled strands per chain:

as,with λ_ES being a factor in the order of 0.3–0.5;^29,84 here it was set to 0.5.

The preexponential factor v_hop is a chain-length-independent quantity dictating the hopping frequency. An upper bound can be obtained by the following equation:³⁴

In a previous investigation,²⁹ v_hop was optimized for reproducing viscoelastic properties based on quiescent (stress relaxation function) and shear-flow (shear viscosity) simulations. The optimal value was lower than the prediction from eqn (33) by ∼6.5 times. Here, we will employ the same scaling and set it to v_hop ~ v⁺_hop/6.5 = 65.21 s⁻¹. Nevertheless, v_hop can be fine-tuned for matching viscoelastic properties, if available.

Simulations in the grand-canonical ensemble require inputting the slip-spring activity which can be estimated as follows:²⁹

with 〈N_create〉 being the average number of slip-spring creation or destruction events by kMC within time Δt_kMC from a single simulation carried out under bulk conditions with constant number of slip-springs (canonical ensemble with respect to slip-springs). 〈n_cands〉 = 4πr_attempt³n_bead/3 is the average number of slip-springs that can be formed from a chain end. Note that the number of slip-springs per bead scales roughly proportionally with the activity, which increases slightly with chain length.²⁹ For a mesoscopic chain with 15 beads (2400 effective monomers), the activity was estimated to z_activ = 0.002235. For additional details the reader is referred to Appendix B of ref. 29.

2.2.6 Kernel optimization. The local density is computed on the level of individual beads with the discretization scheme illustrated in eqn (7). A finite support of r_c = 50 Å was used. In doing so, the local density is shaped by the contribution of approximately ∼ 60 beads. We introduce the following reference units of length, σ_r = [m_bead/ρ_{300 K}]^1/3 = 2.067 × 10⁻⁹ m; mass, m_r = m_bead = 1.328 × 10⁻²³ kg; energy, ε_r = k_BT = 4.14 × 10⁻²¹ J; and time, = 1.17063 × 10⁻¹⁰ s. According to ref. 72, an upper bound for the (reduced) time step can be obtained from the empirical formula, Δt̃_crit ∼ (1/6)r̃_c^2.4[small zeta, Greek, tilde][small kappa, Greek, tilde]_{ex,T=300 K} = 495; or, in real units, Δt_crit ∼ 58 ns. The compensating factors for the weighting kernel (w′ → C_Ww, eqn (8)) and friction ( eqn (28)) were determined from a bulk simulation with 800 beads at 300 K as C_W = ρ_CW=1/ρ_EoS = 0.953 and C_ζ = D_Cζ=1/D_Einstein = 0.900, where D_Einstein = k_BT/(N_mon/beadζ_bead). For additional details the reader is referred to section 3.2 in ref. 72.

Even though successful parameterization of the kernel scheme can restore the desired volumetric properties, the scheme exhibits limitations in accurately capturing structural features at short length scales. As demonstrated in ref. 72, the radial distribution function softens at small interparticle distances, and this effect intensifies with increasing coarsening. The range over which errors in representing structural features are introduced is commensurate with the characteristic length scale of the density discretization scheme.

2.3 Equilibration protocol

An initial configuration with 400 chains, each one with molecular weight 120 kDa (15 beads per chain) was generated with the EMSiPoN code^29,34,85 at T = 300 K and ρ = 1500 kg m⁻³ subject to the constraint that the polymer be confined between the solid surfaces along the z-axis. The connectivity of the strands and slip-springs was described with the fixed image convention.⁸⁶ The generation was realized with a Monte Carlo scheme which conducts reverse sampling from the probability distribution dictated by the free energy of the strands (eqn (5)); for details see section S2.3 in ref. 87.

A common criterion for minimizing finite-size effects is a box dimension-to-chain radius of gyration ratio of approximately four. In our simulations, the box dimension along each direction was set to ∼37.6 nm; roughly five times larger than the chain radius of gyration (∼7.5 nm); test simulations with lateral dimensions of 3.5 and 7 times the chain's R_g produced very similar results.

Fig. 5 illustrates the evolution of box size, normal pressure and number of slip-springs during the stages of the equilibration protocol. Initially, a sample of noninteracting chains (deactivated nonbonded interactions) was simulated in the NVT ensemble in order to efficiently relax their conformations. Next, the nonbonded interactions were activated and the sample was simulated in the NVT ensemble and subsequently in the NL_xL_yσ_zzT ensemble to relax the normal stress. The simulation continued with the grand canonical slip-spring scheme activated subject to the activity listed in Table 5 until the number of slip-springs fluctuated about a plateau value. A hopping frequency 100 times higher than the standard value in Table 5 was used in this step to accelerate slip-spring dispersion. Note that the final slip-spring distribution is a thermodynamic property, independent of the hopping frequency (see Appendix B in ref. 29). Also note that the equilibrium number of slip-springs of the capped polymer is lower than that of a bulk counterpart of the system with the same temperature and number of chains (dashed line in Fig. 5c), because the number of entanglements is lower at the polymer/solid interfaces.^32,36,37

Fig. 5 Evolution of (a) L_z, (b) P_zz and (c) N_SS during the course of the relaxation simulation protocol. The first phase has a duration of 0.003 s in the NVT ensemble without nonbonded interactions (up to the red dashed line). In the second phase, the polymer is simulated in the presence of nonbonded interactions in the NVT ensemble for 0.00025 s (up to the blue dashed line). In the third phase, the system is simulated in the NL_xL_yP_zzT ensemble for 0.001 s (up to the green dashed line). Finally, the system is simulated with the grand canonical slip-spring scheme for 0.03 s.

After thoroughly equilibrating the sample, the surface attachment/detachment scheme was activated subject to the prescribed transition rates (k_a→f, k_f→a) and activation lengths (l_a→f, l_f→a) until the populations of attached/detached beads converged. The resulting concentration depends strictly on K_ads = k_f→a/k_a→f and not on the individual values of k_f→a and k_a→f. This is demonstrated in Fig. 6a which depicts the evolution of the concentration of adsorbed beads (C_a) for various values of K_ads = k_f→a/k_a→f and k_a→f. At long times, C_a plateaus to the concentration dictated by K_ads in eqn (12), regardless of the magnitude of k_a→f. Decreasing k_a→f results in the same concentration, albeit equilibrium is achieved in longer times. According to Fig. 6b, the total bead population close to the interfacial region increases slightly with K_ads, because the effective polymer/surface interactions become more attractive. The increase in density near the solid surface depends on the strength of the attractive interactions from the solid surface relative to the compressibility of the polymer.^88,89

Fig. 6 Evolution of the (a) concentration of adsorbed beads (C_a) and (b) population of adsorbed and free beads within an interfacial domain with thickness r_ads = 19.42 Å, following the relaxation protocol described in Fig. 5. Several cases are examined here where (K_ads, k_a→f/s⁻¹) = (1/4, 10⁶), (1, 10⁶), (4, 10⁶), (1, 10⁵) and (1, 10⁴) corresponding to green, light blue, red, blue and dark blue colours, respectively. The dashed lines in (a) illustrate the prediction of the concentration from eqn (12); C_a ∼ K_ads/(1 + K_ads).

The mesoscopic representations were visualized with the open source Visual Molecular Dynamics (VMD, v1.9.4)⁹⁰ package and the Photopea Photo Editor (last accessed 15/12/2024).⁹¹ The graphs were prepared using the Veusz package (v.3.6.2).⁹²

3. Results

3.1 Principles of tack tests

Tack tests were simulated by elongating the simulation domain along the z-direction at a constant elongation rate according to eqn (34):

L_z(t) = L_z,0(1 + [small epsi, Greek, dot above]_zzt). (34)

The system is aperiodic along the normal (z) direction. The periodic lateral dimensions (L_x, L_y) are fixed, under the assumption that the solid surface is incompressible in comparison with the polymer. The total duration of each elongation experiment is ε_max/[small epsi, Greek, dot above]_zz, with ε_max being the maximum engineering strain.

Unless otherwise stated, the time step was set to Δt = min(Δt_crit/2, k_a→f⁻¹, k_f→a⁻¹). By default, the maximum value was set to half the critical time step (Δt_crit = 59 ns) predicted by the analysis in section 2.2.6.⁷² However, in cases with fast interfacial kinetics, Δt was decreased to the lowest value of the inverse detachment or attachment rate constant to accurately capture attachment and detachment events.

Fig. 7a displays the evolution of the normal stress during the course of the tack test. Fig. 7b illustrates snapshots at strain increments of 0.2. During the initial phase, the stress increases at a rate commensurate with the elastic modulus of the material. This continues until the stress reaches a maximum value, at which point the material can no longer sustain higher loads. Thereafter, stress decreases progressively and cavities expand toward the center of the film until the eventual failure of the interface, either cohesive or adhesive.

Fig. 7 Illustration of an indicative tack test simulation. (a) Stress–strain curve of the tack test. The thin blue curve indicates the instantaneous stress. The thick line depicts the smoothed stress–strain curve with a Savitzky–Golay filter^79,93 over a window Δε = 0.012. The dotted lines indicate the point where the stress becomes maximum. (b) Snapshots of the simulated trajectory at strain increments of 0.2. The system comprises 400 chains, each with a molecular weight of 120 kDa. The periodic x and y box dimensions remain constant. The aperiodic z-dimension of the box is elongated at a steady rate [small epsi, Greek, dot above]_zz = 10³ s⁻¹. Beads within an interfacial region of thickness ∼19.42 Å near each surface (yellow dashed rectangles) can reversibly attach and detach with rates k_a→f = k_f→a = 10⁴ s⁻¹ (K_ads = 1). Here, the activation lengths are set to zero (no force coupling). Attached/detached beads are displayed with red/gray colour. The simulations are conducted with the parameters listed in Table 5.

3.2 Constant-rate interfacial kinetics

A series of tack tests have been realized, spanning a broad parameter range in terms of varying elongation rate [small epsi, Greek, dot above]_zz, detachment rate constant k_a→f, and equilibrium constant for adsorption K_ads, as shown in Table 6. The elongation rate spans three orders of magnitude and ranges from 10² s⁻¹ to 10⁴ s⁻¹. The desorption rate constant spans four orders of magnitude, ranging between 10² s⁻¹ to 10⁵ s⁻¹. For reference, the desorption rate constant of the weakly adsorbing PE/graphite interface³² is ∼1.73 × 10⁷ s⁻¹. The adsorption equilibrium constant has been set to 1/4, 1 and 4, values indicative of weak, moderate and good adhesion, respectively. The thickness of the interfacial region was set to the radius of gyration corresponding to the mesoscopic beads, . The activation lengths in eqn (1) and (2) will be set to zero for the time being (l_f→a = l_a→f = 0); therefore, attachment/detachment rates remain constant throughout the simulation; and in eqn (1) and (2). This assumption is relevant in situations where strong bonds are formed between the polymer and the solid surface; e.g., detachment of strong adhesives. Given that the energy barriers are significantly high, the effect of the activation length term in eqn (1) and (2) becomes minimal. The maximum engineering strain will be set to ε_max = 1 which is large enough for the cavitation/fibrillation processes to kick in and the damage of the polymer/solid interface to be considered irreversible. This contraint on the maximum strain will be relaxed in the next section. For each combination of parameters we conducted three individual experiments starting with a different configuration. A total of 108 simulated tack tests have been conducted serially, taking a total of ∼20 days on a Intel(R) Xeon(R) Silver 4314 CPU@2.40GHz processor with 16 physical cores.

Table 6 Parameters of tack tests

Quantity	Value	Description
[small epsi, Greek, dot above]zz (s^–1)	10^2–10^4	Elongation rate
ka→f (s^–1)	10^2–10^5	Detachment rate
Kads = kf→a/ka→f	0.25, 1, 4	Equilibrium adsorption constant
rads = Rg,bead (Å)	19.42	Interfacial thickness
lf→a = la→f	0	Activation lengths

---

Fig. 8 illustrates the stress–strain curves from the tack tests in dependence of the elongation rate [small epsi, Greek, dot above]_zz (increasing from top-to-bottom), detachment rate k_a→f (increasing from left-to-right) and equilibrium constant for adsorption K_ads (green < blue < red). The corresponding snapshots for K_ads = 1/4, 1 and 4 are illustrated in Fig. 9, 10 and 11, respectively.

Fig. 8 Stress–strain curves in dependence of the detachment rate k_a→f (increasing from left to right), elongation rate [small epsi, Greek, dot above]_zz (increasing from top to bottom) and adsorption constant K_ads = k_f→a/k_a→f = 1/4 (green), 1 (blue) and 4 (red). Scenarios where [small epsi, Greek, dot above]_zz is smaller than (larger than or equal to) k_a→f are indicated with green (red) shades. The stress has been smoothed with a Savitzky–Golay filter^79,93 over a window Δε = 0.012.

Fig. 9 Snapshots of the final mesoscopic configurations at various detachment rates (increasing from left to right) and elongation rates (increasing from top to bottom) for K_ads = 1/4. Green/red shaded plots indicate total/partial detachment until ε_max = 1.

Fig. 10 Same as Fig. 9 but for K_ads = 1.

Fig. 11 Same as Fig. 9, for K_ads = 4.

In situations where the strain-rate is slower than the detachment rate constant ([small epsi, Greek, dot above]_zz < k_a→f) the interface has adequate time to respond to the deformation and relax the imposed tensile stress. Such conditions promote fast detachment in general, albeit the effect is sensitive to the equilibrium adsorption constant. To elaborate, for the smallest equilibrium adsorption constant considered here (K_ads = 1/4), we find a consistent behavior where failure is strictly adhesive across this regime; e.g., panels with green shade in Fig. 8 and 9. For intermediate values (K_ads = 1) the polymer detaches fully at ε = 1, with two exceptions; ([small epsi, Greek, dot above]_zz, k_a→f) = (10², 10³) and (10³, 10⁴) s⁻¹ in Fig. 10. However, for K_ads = 4, no instances of complete polymer detachment from the surface were observed at ε = 1 within the parameter range investigated here; see mesoscopic configurations in Fig. 11.

Scenarios where the elongation rate is faster than the detachment rate constant, [small epsi, Greek, dot above]_zz ≥ k_a→f (red-shaded plots in Fig. 8), promote partial detachment, in general, because mesoscopic segments do not have adequate time to respond to the deformation. For K_ads = 1, the adhesion between the polymer and the surface is already quite strong; therefore, further increasing K_ads does not significantly affect the system's response. Interestingly, when the timescale separation between deformation and detachment rate is large ([small epsi, Greek, dot above]_zz ≫ k_a→f) the system's response becomes insensitive to the latter. For example, note that the two leftmost panels in the last row of Fig. 8 are practically identical; the same applies in Fig. 9–11 as well. This occurs because the detachment rate becomes negligible in relation to the elongation rate and thus only the initial concentration of adsorbed segments (dictated by K_ads) affects the system's response.

According to our simulations, the general trend is that strongly adsorbing surfaces (high K_ads) result in slow detachment (cohesive failure is promoted), whereas weakly adsorbing ones (low K_ads) in fast detachment (promotion of adhesive failure). According to ref. 23 and 27, increasing the crosslink density (which in turn, improves the cohesion of the polymer network) enhances the tendency for adhesive failure, and vice versa. It appears that the failure mechanism is a strong function of the ratio between strength of the polymer–polymer (cohesion) and surface-polymer (adhesion) interactions. Supposing that the adhesion of the polymer to the solid is much stronger than the cohesive interactions, the sample is expected to fail cohesively. On the other hand, a weakly adsorbing surface is expected to fail adhesively.

According to experimental tack tests,^5,9 the roughness and fibrillation degree at the fractured material surface depend strongly on the imposed elongation rate. Low deformation rates result in sparse distributions of thick fibers and large cavities. High deformation rates, on the other hand, result in dense distributions of thin fibers.

The mesoscopic model yields a similar behavior:

• At slow elongation rates we observe thick and sparse fibers forming bridges between the edges of the simulation box. This response is illustrated clearly in the top row of Fig. 9–11 for [small epsi, Greek, dot above]_zz = 10² s⁻¹ with red shades (partial detachment).

• At moderate rates ([small epsi, Greek, dot above]_zz = 10³ s⁻¹), we observe denser distributions of thinner fibers, as shown in the second row of Fig. 9–11.

• At the highest rates considered here ([small epsi, Greek, dot above]_zz = 10⁴ s⁻¹), the rate of elongation of the box boundaries (and the velocity of the adsorbed beads) becomes commensurate with the relaxation rate on the strand level (∼ζ_bead/m_bead). This rapid elongation induces a shock, where strands lack sufficient time to respond, leading to the formation of a low density region between the adsorbed layer and the bulk polymer. This region is populated by highly ordered strands and is indicative of extreme fibrillation conditions.

Overall, we observe a trend toward smaller cavities/denser distributions of thinner fibers as the elongation rate increases and a progressive orientation of the polymer chains^5,94 along the pulling direction. Similar fibril-like nanomorphologies and cavitation tendencies have been reported in relevant simulation studies.^11,22,23,27

Fig. 12 illustrates the yield tensile stress (σ_y) at the end of the reversible portion (also see Fig. 7a) of the stress–strain curve, averaged over three individual tack tests performed for each scenario considered in Table 6. The yield stress is an increasing function of the elongation rate [small epsi, Greek, dot above]_zz, e.g., inspect Fig. 12 from top to bottom. This happens because the polymer cannot relax the imposed stress efficiently under fast elongation rates. This response conforms with reports from experimental tack tests^5,9 and relevant atomistic^17,18 and mesoscopic²³ simulations. In addition, the maximum stress reported from experimental tack tests⁵ is in the order of 0.1–1 MPa; in good quantitative agreement with our findings. K_ads indicates higher adsorption rates and enhanced concentration of adsorbed beads (see eqn (12)), therefore increased adhesion. In situations where [small epsi, Greek, dot above]_zz ≥ k_a→f (red background in Fig. 12) σ_y increases from K_ads = 1/4 to 1. Further increasing K_ads does not increase stress, indicating that above a threshold value the concentration of the adsorbed beads at the interface is high enough and cavitation is realized slightly further towards the bulk phase.

Fig. 12 Yield stress σ_y as a function of the detachment rate k_a→f, elongation rate [small epsi, Greek, dot above]_zz and adsorption constant K_ads = k_f→a/k_a→f = 1/4 (green), 1 (blue) and 4 (red). Scenarios where [small epsi, Greek, dot above]_zz is smaller (larger or equal) than k_a→f are indicated with green (red) background. The error bars correspond to the standard deviation from three individual experiments.

In cases where [small epsi, Greek, dot above]_zz < k_a→f, σ_y increases monotonically in most cases (panels with green background in Fig. 12). In cases with pronounced timescale separation k_a→f ≪ [small epsi, Greek, dot above]_zz, the maximum stress is practically insensitive to detachment rate. For instance, the panels ([small epsi, Greek, dot above]_zz/s⁻¹, k_a→f/s⁻¹) = (10⁴, 10²), (10⁴, 10³), (10⁴, 10⁴) are identical with each other in Fig. 12.

The evolution of stress during a tack test is directly related to the strength of the cohesive interactions within the polymer. If adhesion to the surface is significant, the larger the elastic modulus of the polymer, the higher the resistance to deformation (i.e., the induced stress).^5,9 This effect is demonstrated in Fig. 13 which compares the stress–strain curves for samples with increasing bulk modulus, or conversely, decreasing compressibility. The red curve shows the case where we consider the default parameters in Table 5, whereas the blue and green curves illustrate cases where the compressibility has been decreased by 10 and 100 times, respectively. It is evident that decreasing compressibility results in steeper stress increase at low strains and higher yield stress. In addition, the samples become more brittle with decreasing compressibility, as the location of yield stress moves toward zero.

Fig. 13 Stress–strain curve from a tack experiment with [small epsi, Greek, dot above]_zz = k_a→f = 10³ s⁻¹ and K_ads = 1. Red curve illustrates a case using the default parameters in Table 1. Blue/green curve illustrates the impact of decreasing compressibility by a factor of 10/100 upon multiplying the characteristic pressure P* by 10/100, respectively.

3.3 Effect of pulling force

The present section relaxes the assumption incorporated in the previous section and considers the effect of coupling the normal forces exerted on the beads to the interfacial kinetics. The coupling is tuned by the activation lengths in eqn (1) and (2). The higher the activation length, the higher the coupling with the pulling forces. As a result, beads near the interface that experience pulling forces toward the bulk region are more likely to detach, and less likely to attach to the interface. On the contrary, beads that experience forces toward the solid surface will attach more frequently and detach at a slower pace. According to Transition State Theory, the activation length is commensurate with the distance between the coordinate of the current state (attached/detached) and the coordinate at the free energy barrier beyond which a transition in state occurs (detached/attached). A recent model for the fracture of polymer chains⁹⁵ demonstrated that, even though this approach correctly estimates the order of magnitude of the activation length, it can underestimate it by a factor 38–60%.

To simplify our analysis we will consider the extreme scenarios shown Table 6. Specifically, we will address weakly ( = 1/4) and strongly ( = 4) adsorbing interfaces exhibiting either slow ( = 10² s⁻¹) or fast ( = 10⁵ s⁻¹) kinetics. Note that denotes the normal force-independent factor of the transition rate, which increases with decreasing transition barriers E_i→j, as indicated in eqn (1) and (2). Several values for the activation length will be considered here, with a maximum value set to R_g,bead = 19.42 Å. We will consider the faster elongation rates in Table 6 ([small epsi, Greek, dot above]_zz = 10³ s⁻¹ and 10⁴ s⁻¹) in order to access very large deformations. The time step is set to 29 ns for the case [small epsi, Greek, dot above]_zz = 10³ s⁻¹, and 2.9 ns for the case [small epsi, Greek, dot above]_zz = 10⁴ s⁻¹ and the normal forces are averaged for 2000 time steps in each case. In doing so, the normal forces are averaged within the same number of steps and “strain window”, Δε = 2000[small epsi, Greek, dot above]_zzdt = 0.058. For each scenario, two simulations were performed, each initialized with a different initial condition. A total of 32 simulated tack tests have been conducted serially, taking a total of ∼10 days on a Intel(R) Xeon(R) Silver 4314 CPU@2.40GHz processor with 16 physical cores.

Fig. 14a illustrates the stress–strain curves of weakly (K_ads = 1/4) and strongly (K_ads = 4) adsorbing interfaces as a function of the activation length at various elongation rates. The detachment rate has been set to = 10² s⁻¹ indicative of slow interfacial kinetics (high transition barriers).

Fig. 14 (a) Stress–strain curves and (b) toughness as a function of the activation length l_a→f, elongation rate [small epsi, Greek, dot above]_zz = 10³ s⁻¹ (solid line, squares), 10⁴ s⁻¹ (dashed line, circles) and equilibrium adsorption constant = 1/4 (green) and 4 (red). In all cases = 10² s⁻¹. The legends in (a) denote the failure mode (adhesive or cohesive) and the percentage of raptured strands for each case.

In situations where l_a→f = 0 (leftmost panels of Fig. 14a), there is delayed failure at higher elongation strains and significantly increased stress. Interestingly, for the fastest elongation rate, ∼1% of chains fractured during the course of the simulation, indicating a tendency for partial cohesive failure. The stress is maximized in this case and reaches ca. 4 MPa. This behavior is expected, since the detachment rate is two orders of magnitude lower than the elongation rate; therefore, the segments cannot respond efficiently to the rapidly imposed deformation. This is also illustrated by Fig. 15 which illustrates the evolution of a system with K_ads = 4 and l_a→f = 0 for [small epsi, Greek, dot above]_zz = 10³ s⁻¹ and [small epsi, Greek, dot above]_zz = 10⁴ s⁻¹ (see panels (a) and (b), respectively).

Fig. 15 Evolution of the configurations for [l_a→f/Å, [small epsi, Greek, dot above]_zz/s⁻¹] equal to (a) [0, 10³], (b) [0, 10⁴], (c) [9.71, 10³] and (d) [9.71, 10⁴]. In all cases, = 4 and = 10² s⁻¹.

With increasing activation length the polymer detaches from the solid surface at much lower strains. The effect is very pronounced in the weakly adsorbing interface where K_ads = 1/4. In addition, the trend observed in the previous case (l_a→f = 0) reverses and detachment is realized at lower strains with increasing elongation rate. E.g., compare panels (c) and (d) of Fig. 15. This response is in accordance with theoretical¹⁷ and experimental^5,9 observations.

The reversible portion of the stress–strain curve preceding the yield strain (ca. 0.06) does not depend on the value of the activation length. That is to say, the low strain portion in Fig. 14a is practically the same from left to right. Finally, the maximum value of the stress is an increasing function of the elongation rate for all cases considered in Fig. 14a.

Fig. 16a depicts the same quantities as in Fig. 14a but for interfaces with fast detachment kinetics ( = 10⁵ s⁻¹), indicative of low transition barriers. Because detachment rate is much faster than deformation rate, failure occurs prematurely at significantly lower elongation strains. With increasing activation length detachment is facilitated, though the enhancement is less pronounced than in cases with slow interfacial kinetics. Failure is always adhesive for this detachment rate.

Fig. 16 As in Fig. 14 for = 10⁵ s⁻¹. In all cases, the samples exhibit cohesive failure, and no chain rupture events are observed.

The toughness of the interfaces can be quantified by the integral of the stress–strain curve:

up to the failure strain ε_f. Evaluations of the toughness for the systems considered in Fig. 14a and 16a are illustrated in Fig. 14b and 16b, respectively. Toughness increases consistently as the adsorption becomes stronger (i.e., as K_ads increases). Indeed, for a given elongation rate, systems with K_ads = 4 (red curves) exhibit higher toughness than those with K_ads = 1/4 (green curves). The activation length plays a crucial role in the detachment mechanism and the resulting toughness of the interface. When the activation length is zero, toughness increases considerably and becomes a strongly increasing function of elongation rate. Conversely, toughness decreases with increasing activation length because the polymer detaches at progressively lower elongation strains in most cases. The strain-rate dependence of toughness decreases at high activation lengths because, even though stress becomes higher with increasing strain, the film detaches at lower strains.

4. Concluding remarks

The present work develops a generic mesoscopic model, capable of addressing high molar mass amorphous polymers above their glass temperature capped between two solid surfaces, under quiescent and nonequilibrium conditions. The latter are imposed by subjecting the samples to constant rate elongation across the aperiodic (normal) direction, as in the so-called tack tests. A detailed six-stage parameterization framework is established for deriving essential parameters of the mesoscopic model, including the optimization of the equation of state coefficients, the mesoscopic force fields, and chain dynamics in the Rouse and entangled regimes. The framework is generic and applicable to polymers of arbitrary architecture and chemical constitution. The parameterization can be performed either in a bottom-up or a top-down fashion, based on input from atomistic simulations or experimental data, respectively.

As a case study, we apply the parameterization framework based on a minimal data set from an indicative noncrystallizable fluoropolymer. The agreement of the volumetric properties, conformational statistics and dynamics between the mesoscopic model and the atomistic reference data set is excellent.

We present results from simulated tack experiments conducted over a wide parameter space, including: the rate of elongation [small epsi, Greek, dot above]_zz, detachment rate constant k_a→f, activation lengths of the attachment/detachment transitions (l_a→f, l_f→a), and the adsorption equilibrium constant K_ads = k_f→a/k_a→f. The resulting failure mechanism depends on a complex interplay between these parameters. The attachment/detachment rate constants dictate the frequency at which the polymer chains attach/detach from the surface. The magnitudes of the activation lengths tune the coupling of the attachment/detachment rates with the normal forces exerted on the mesoscopic beads at the interface. The adsorption equilibrium constant is a thermodynamic property and determines the adhesive strength of the interface. The performance of each interface was characterized by stress–strain curves from which yield strain and toughness were extracted.

We begin the recapitulation of our findings with situations where the activation length is negligible in relation to the transition energy barriers or, alternatively, the normal forces generated by the polymer at the interface are weak. This scheme is suitable for modeling the detachment of strong adhesives and also highly compressible materials (such as foams) that generate weak normal forces.

When the detachment rate is faster than the elongation rate (k_a→f > [small epsi, Greek, dot above]_zz), the interface can respond efficiently and relax the imposed stress. Weakly adsorbing interfaces (K_ads = 1/4) promote cavitation very close to the surface, enabling easy detachment of the polymer at short times. The corresponding toughness becomes minimal in these cases. With increasing K_ads, the toughness increases slightly.

In situations where the elongation rate is faster than the detachment rate ([small epsi, Greek, dot above]_zz > k_a→f), pronounced fibrillation tendencies are observed. Interestingly, when the timescale separation is very large ([small epsi, Greek, dot above]_zz ≫ k_a→f), the magnitude of the detachment rate becomes irrelevant and the outcome depends strictly on the initial concentration of attached beads, dictated by the equilibrium adsorption constant. The toughness of the interface is maximized under these conditions and becomes a strongly increasing function of elongation rate.

We now turn our attention to situations with finite activation lengths. Note that the effect of force coupling becomes more pronounced as the transition barriers decrease (low E_i→j in eqn (1) and (2)) and also when modelling incompressible materials. Contrary to the previous case, the polymer detaches more efficiently (at lower strains) with increasing elongation rate, because the attached beads are pulled more strongly toward the polymer film with increasing strain. Indeed, when the normal force coupling is very strong (l_a→f ≫ E_a→f/〈F_z〉) toughness decreases considerably, since the polymer detaches from the surface very readily, even for the smallest [small epsi, Greek, dot above]_zz/k_a→f ratio considered here. This response conforms with the picture from experimental tack tests.^5,9

We note a consistent trend where toughness is a strictly increasing function of K_ads. In addition, regardless of the activation length, the yield stress increases monotonically with elongation rate due to the inability of the interface to efficiently relax the imposed stress. At high elongation rates the polymer becomes stiff, leading to the formation of dense distributions of thinner fibers. At low elongation rates, on the other hand, the polymer can relax the imposed stress efficiently and the behavior becomes more ductile, leading to the formation of large cavities—in other words thick and sparse fibers, forming bridges between the edges of the simulation box. Furthermore, the yield stress exhibits an inverse relationship with compressibility. Low compressibility is indicative of stiffer samples that resist deformation and yield at increased stress under load.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

A. P. Sgouros thanks Syensqo S.A., Bollate, Italy, for financial support. The authors would like to congratulate Prof. J. J. de Pablo on his birthday and wish him many fruitful years in science to come. S. K. and S. C. are grateful to Prof. Theodorou and Dr. Sgouros for their productive stay at NTUA in January 2023.

References

1. A. J. Kinloch, Adhesion and adhesives: science and technology, Springer Science& Business Media, 2012.
2. R. A. Pethrick, Proc. Inst. Mech. Eng., Part L, 2015, 229, 349–379.
3. H. R. Gualberto, F. do Carmo Amorim and H. R. M. Costa, J. Braz. Soc. Mech. Sci. Eng., 2021, 43, 389.
4. D. Son, V. Liimatainen and M. Sitti, Small, 2021, 17, 2102867.
5. C. Creton and P. Fabre, in Adhesion Science and Engineering, ed. D. A. Dillard, A. V. Pocius and M. Chaudhury, Elsevier Science B.V., 2002, pp. 535–575.
6. C. Creton and M. Ciccotti, Rep. Prog. Phys., 2016, 79, 46601.
7. A. Zosel, Colloid Polym. Sci., 1985, 263, 541–553.
8. A. Zosel, J. Adhes., 1989, 30, 135–149.
9. H. Lakrout, P. Sergot and C. Creton, J. Adhes., 1999, 69, 307–359.
10. A. Chiche, P. Pareige and C. Creton, C. R. Acad. Sci., Ser. IV:Phys., Astrophys., 2000, 1, 1197–1204.
11. M. Solar, Z. Qin and M. J. Buehler, J. Mater. Res., 2014, 29, 1077–1085.
12. T.-P. Fries and M. Baydoun, Int. J. Numer. Methods Eng., 2012, 89, 1527–1558.
13. P.-A. Guidault, O. Allix, L. Champaney and C. Cornuault, Comput. Meth. Appl. Mech. Eng., 2008, 197, 381–399.
14. Z. Zhang, X. Wang and J. D. Lee, J. Appl. Phys., 2014, 115, 114314.
15. S.-H. Cheng and C. T. Sun, J. Nanomechanics Micromech., 2014, 4, 1–8.
16. L. Brochard, G. Hantal, H. Laubie, F. J. Ulm and R. J. M. Pellenq, Int. J. Fract., 2015, 194, 149–167.
17. G. Xu and H. Wang, Constr. Build. Mater., 2016, 121, 246–254.
18. K. Min, A. R. Rammohan, H. S. Lee, J. Shin, S. H. Lee, S. Goyal, H. Park, J. C. Mauro, R. Stewart, V. Botu, H. Kim and E. Cho, Sci. Rep., 2017, 7, 1–11.
19. D. Gersappe and M. O. Robbins, Europhys. Lett., 1999, 48, 150–155.
20. S. V. Kallivokas, A. P. Sgouros and D. N. Theodorou, Phys. Rev. E, 2020, 102, 30501.
21. I. Delasoudas, S. V. Kallivokas and V. Kostopoulos, Polymer, 2024, 296, 126830.
22. M. J. Stevens, Macromolecules, 2001, 34, 2710–2718.
23. A. Baggioli, M. Casalegno, A. David, M. Pasquini and G. Raos, Macromolecules, 2021, 54, 195–202.
24. A. R. C. Baljon, J. Vorselaars and T. J. Depuy, Macromolecules, 2004, 37, 5800–5806.
25. A. F. Terzis, D. N. Theodorou and A. Stroeks, Macromolecules, 2002, 35, 508–521.
26. E. Voyiatzis and A. Stroeks, J. Phys. Chem. B, 2022, 126, 6102–6111.
27. K. Jin, D. López Barreiro, F. J. Martin-Martinez, Z. Qin, M. Hamm, C. W. Paul and M. J. Buehler, Polymer, 2018, 154, 164–171.
28. M. H. Nafar Sefiddashti, B. J. Edwards and B. Khomami, J. Rheol., 2015, 59, 119–153.
29. A. P. Sgouros, G. Megariotis and D. N. Theodorou, Macromolecules, 2017, 50, 4524–4541.
30. M. H. Nafar Sefiddashti, B. J. Edwards, B. Khomami and E. S. G. Shaqfeh, J. Chem. Phys., 2021, 154, 204907.
31. M. H. Nafar Sefiddashti, B. J. Edwards and B. Khomami, Phys. Rev. Lett., 2018, 121, 247802.
32. A. P. Sgouros, G. G. Vogiatzis, G. Megariotis, C. Tzoumanekas and D. N. Theodorou, Macromolecules, 2019, 52, 7503–7523.
33. A. E. Likhtman, Macromolecules, 2005, 38, 6128–6139.
34. G. G. Vogiatzis, G. Megariotis and D. N. Theodorou, Macromolecules, 2017, 50, 3004–3029.
35. T. Uneyama and Y. Masubuchi, J. Chem. Phys., 2011, 135, 184904.
36. A. Ramírez-Hernández, B. L. Peters, M. Andreev, J. D. Schieber and J. J. de Pablo, J. Chem. Phys., 2015, 143, 243147.
37. J. Kirk, Z. Wang and P. Ilg, J. Chem. Phys., 2019, 150, 094906.
38. M. Langeloth, Y. Masubuchi, M. C. Böhm and F. Müller-plathe, J. Chem. Phys., 2013, 138, 104907.
39. T. Uneyama and Y. Masubuchi, J. Chem. Phys., 2012, 137, 154902.
40. V. C. Chappa, D. C. Morse, A. Zippelius and M. Müller, Phys. Rev. Lett., 2012, 109, 148302.
41. A. Ramírez-Hernández, B. L. Peters, L. Schneider, M. Andreev, J. D. Schieber, M. Müller and J. J. de Pablo, J. Chem. Phys., 2017, 146, 014903.
42. A. Ramírez-Hernández, F. A. Detcheverry, B. L. Peters, V. C. Chappa, K. S. Schweizer, M. Müller and J. J. de Pablo, Macromolecules, 2013, 46, 6287–6299.
43. Y. Masubuchi, Annu. Rev. Chem. Biomol. Eng., 2014, 5, 11–33.
44. Y. Masubuchi, J. Chem. Phys., 2015, 143, 224905.
45. Y. Masubuchi, J. I. Takimoto, K. Koyama, G. Ianniruberto, G. Marrucci and F. Greco, J. Chem. Phys., 2001, 115, 4387–4394.
46. S. F. Edwards, Proc. Phys. Soc., 1967, 92, 9–16.
47. P. G. de Gennes, J. Chem. Phys., 1971, 55, 572.
48. M. Doi, J. Polym. Sci., Polym. Phys. Ed., 1983, 21, 667–684.
49. S. T. Milner and T. C. B. McLeish, Phys. Rev. Lett., 1998, 81, 725–728.
50. A. E. Likhtman and T. C. B. McLeish, Macromolecules, 2002, 35, 6332–6343.
51. N. Evangelou, G. Megariotis, A. P. Sgouros, G. G. Vogiatzis, N. A. Romanos and D. N. Theodorou, AIP Conf. Proc., 2021, 2343, 130008.
52. A. P. Philippas, A. P. Sgouros, G. Megariotis and D. N. Theodorou, AIP Conf. Proc., 2021, 2343, 130003.
53. G. Megariotis, G. G. Vogiatzis, A. P. Sgouros and D. N. Theodorou, Polymers, 2018, 10, 1156.
54. A. P. Sgouros, G. G. Vogiatzis, G. Kritikos, A. Boziki, A. Nikolakopoulou, D. Liveris and D. N. Theodorou, Macromolecules, 2017, 50, 8827–8844.
55. A. P. Sgouros, A. T. Lakkas, G. Megariotis and D. N. Theodorou, Macromolecules, 2018, 51, 9798–9815.
56. V. A. Harmandaris, V. G. Mavrantzas, D. N. Theodorou, M. Kröger, J. Ramírez, H. C. Öttinger and D. Vlassopoulos, Macromolecules, 2003, 36, 1376–1387.
57. M. Ansari, A. Alabbas, S. G. Hatzikiriakos and E. Mitsoulis, Int. Polym. Process., 2010, 25, 287–296.
58. Y. Masubuchi, Y. Doi and T. Uneyama, J. Phys. Chem. B, 2022, 126, 2930–2941.
59. P. Ilg and M. Kröger, J. Rheol., 2011, 55, 69.
60. P. S. Stephanou, I. C. Tsimouri and V. G. Mavrantzas, Macromolecules, 2016, 49, 3161–3173.
61. T. Schweizer, J. van Meerveld and H. C. Öttinger, J. Rheol., 2004, 48, 1345.
62. S.-Q. Wang, Nonlinear polymer rheology: macroscopic phenomenology and molecular foundation, John Wiley & Sons, 2018.
63. A. Ramírez Hernández, M. Müller and J. J. de Pablo, Soft Matter, 2013, 9, 2030–2036.
64. D. Nichetti and I. Manas-Zloczower, J. Rheol., 1998, 42, 951–969.
65. M. Doi and S. F. Edwards, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 1802–1817.
66. C. C. Hua and J. D. Schieber, J. Chem. Phys., 1998, 109, 10018–10027.
67. J. D. Schieber and M. Andreev, Annu. Rev. Chem. Biomol. Eng., 2014, 5, 367–381.
68. Y. Masubuchi, Macromolecules, 2018, 51, 10184–10193.
69. Y. Masubuchi and T. Uneyama, Soft Matter, 2018, 14, 5986–5994.
70. M. Kröger, J. Nonnewton. Fluid Mech., 2015, 223, 77–87.
71. I. C. I. C. Sanchez and R. H. R. H. Lacombe, Macromolecules, 1978, 11, 1145–1156.
72. A. P. Sgouros and D. N. Theodorou, J. Phys. Chem. B, 2024, 128, 6907–6921.
73. L. B. Lucy, Astron. J., 1977, 82, 1013–1024.
74. K. F. Mansfield and D. N. Theodorou, Macromolecules, 1991, 24, 4295–4309.
75. J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1949, 17, 338–343.
76. J. T. Padding and W. J. Briels, J. Chem. Phys., 2001, 114, 8685–8693.
77. L. J. Fetters, D. J. Lohse and R. H. Colby, Phys. Prop. Polym. Handb., 2006, 445–452.
78. J. A. Nelder and R. Mead, Comput. J., 1965, 7, 308–313.
79. P. Virtanen, R. Gommers, T. E. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser and J. Bright, Nat. Methods, 2020, 17, 261–272.
80. D. N. Theodorou, in Computer Simulations of Surfaces and Interfaces, ed. B. Dünweg, D. P. Landau and A. I. Milchev, 2003, NATO Science Series, pp. 329–419.
81. D. N. Theodorou, in Polymers at Surfaces and Interfaces, ed. B. Dünweg, D. P. Landau and A. I. Milchev, 2003, NATO Science Series, pp. 329–419.
82. V. A. Harmandaris, V. G. Mavrantzas and D. N. Theodorou, Macromolecules, 1998, 31, 7934–7943.
83. Y. H. Lin, Macromolecules, 1987, 20, 3080–3083.
84. Y. Masubuchi, G. Ianniruberto, F. Greco and G. Marrucci, J. Chem. Phys., 2003, 119, 6925–6930.
85. G. Megariotis, G. G. Vogiatzis, L. Schneider, M. Müller and D. N. Theodorou, J. Phys.: Conf. Ser., 2016, 738, 012063.
86. A. P. Sgouros and D. N. Theodorou, Computation, 2023, 11, 106.
87. A. P. Sgouros, S. Knippenberg, M. Guillaume and D. N. Theodorou, Soft Matter, 2021, 17, 10873–10890.
88. A. P. Sgouros, C. J. Revelas, A. T. Lakkas and D. N. Theodorou, J. Phys. Chem. B, 2022, 126, 7454–7474.
89. C. J. Revelas, A. P. Sgouros, A. T. Lakkas and D. N. Theodorou, Computation, 2021, 9, 57.
90. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
91. I. Kuckir, https://www.photopea.com/.
92. J. Sanders, Veusz-3.6.2, https://veusz.github.io/.
93. A. Savitzky and M. J. E. Golay, Anal. Chem., 1964, 36, 1627–1639.
94. R. J. Good and R. K. Gupta, J. Adhes., 1988, 26, 13–36.
95. D. N. Theodorou, J. Chem. Phys., 2024, 161, 184905.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4me00199k

---