Across all the tests, peeling begins along the outer edge of the circular adhered surface at a site of imperfect adhesion between the sample and the glass plate, evident from the symmetry breaking in Fig. 2. This pattern is in accordance with analytical models of the delamination of tubular adhesive joints,²⁷ in which imperfections on the outer edge of the joint lead to higher stresses and increased likelihood of shear failure. Interestingly, our experiments show that the level of defect sensitivity in the initiation and propagation of delamination can be tuned by applying a normal force, as seen by the increasing repeatability among tests with increasing normal force in Fig. 3.

Once delamination initiates, it propagates circumferentially, as seen by the shrinkage of the red-shaded regions in Fig. 2 with increasing rotation angle. The propagation occurs in a stick-slip manner, such that regions that have previously delaminated re-adhere and later delaminate again. For samples with no applied normal force, delamination propagates through the sample's entire radius from the onset of peeling (Fig. 2a). For samples with an applied normal force, the stick-slip cycles start at the edge and spiral inward as they move around the circumference (Fig. 2b and c). The propagation of delamination that we observe is in line with results from prior torsion studies: the inward-radial progression matches the findings of Chateauminois et al.,²⁴ and the circumferential stick-slip cycles are broadly consistent with the work of Chaudhury and Chung.²⁸ This stick-slip process is visible as alternating jumps and plateaus in the delaminated radius percentage in Fig. 2 and corresponds to fluctuations in the otherwise-monotonic normal force and torque graphs (Fig. 3; also see ESI† for more detail).

In all samples, we eventually observe the termination of delamination. The reason for the termination depends on the applied normal force. For samples with non-zero applied force, the delamination stops at some radial distance from the center, giving way to fracturing. For samples without applied compression, the entire sample delaminates by the end of the test and does not fracture. This behavior may be seen by contrasting the observations in Fig. 2: in the absence of normal force, the sample fully delaminates by 92 degrees of rotation, while the plots for samples with an initial normal force reach a plateau prior to delaminating through the whole surface.

Aside from the three observations detailed above, we identified experimental evidence of several other phenomena. In Fig. 3, the normal compression is shown to increase as the rotation angle increases. This finding is consistent with the well-documented ‘Poynting effect’,³⁶ which predicts that torsion elongates cylinders, so compression must be applied to preserve their height.^37,38 At the very beginning of rotation, the normal force sometimes decreases before it begins to increase – one possible explanation is that any initial misalignment can cause the cylinder to reconfigure at onset of loading. This effect induces an offset of the normal force throughout the loading, which becomes more pronounced for larger samples, as shown in the ESI.† In addition, we observed the emergence of wrinkles on the lateral surface of the cylinders once the sample was rotated beyond a critical angle, in accordance with the theory proposed by Ciarletta and Destrade.³⁹ While the fractures on the lateral surface have been discussed briefly above, more investigation is required to understand the particular fracturing patterns that emerge in these compression/torsion tests. More detail on the wrinkling and fracturing patterns we observed is provided in the ESI.†

Based on the experimental data, delamination of a circular surface under combined axial compression and torsion is a process in which an initial imperfection in the adhesion between the PDMS sample and the glass plate creates a delamination front that propagates circumferentially and radially inward. The sample often re-adheres to the substrate after delamination, creating successive stick-slip waves. Fluctuations in the trends of normal force and torque as a function of the rotation angle correspond to snap-through events where delamination progresses rapidly. The delamination process is highly sensitive to imperfections, though application of a normal force is shown to reduce the level of imperfection sensitivity and make the process considerably more repeatable. Moreover, by resisting delamination, applied compression limits the radial range of delamination and eventually drives the sample to fracture rather than delaminate the innermost part of the surface. However, if no normal force is initially applied, the sample delaminates across the whole interface and does not fracture. The physical mechanisms behind these key observations are discussed through a theoretical lens below.

To describe the experimental system illustrated in Fig. 1, we consider an elastic cylinder made of a hyperelastic, isotropic, homogeneous, and incompressible material. Its undeformed axial length is denoted by H, and it has an undeformed radius of R₀. Material points are labeled using the cylindrical coordinate system, as X = X(R, Θ, Z) in the undeformed state, such that10 ≤ R ≤ R₀, 0 ≤ Θ < 2π, −H/2 ≤ Z ≤ H/2,as illustrated in Fig. 4. Upon deformation, a mapping function x = χ(X) assigns the material points to their location in the deformed cylindrical coordinates x = x(r, θ, z).

The cylinder is initially pushed against a rigid substrate (at Z = −H/2), by application of a normal stress at the top surface (Z = H/2), which is then rotated quasistatically to increase the prescribed rotation angle α. Delamination is permitted only at the bottom surface and, with the assumption of axial symmetry, is restricted to the annular region R_d ≤ R ≤ R₀, where the innermost radius of the delaminated region, R_d = R_d(α), can vary throughout the loading. Understanding the propagation of this delamination front and its dependence on model parameters is the central goal of this formulation.

To simplify the mathematical derivation, a set of kinematic assumptions are made. First, we assume that horizontal planes remain horizontal throughout the deformation. Second, we impose axial symmetry such that the deformation field is independent of Θ. As a result of these two assumptions, the vertical displacement u_z of any point on a horizontal plane is independent of the radial and angular coordinates (i.e. u_z = u_z(Z)). We restrict our attention to small initial vertical deformations imposed by the applied normal force, which results in a uniform vertical stretch (λ_z = ∂z/∂Z) that is held constant throughout the torsion process. Lastly, we assume that the cylinder preserves its cylindrical shape throughout the deformation process, namely that radial displacements are independent of the vertical coordinate, i.e. the radial displacement function becomes u_r = u_r(R). This assumption neglects any barreling effects that may occur in the initial compression, essentially implying that the initial compression permits sliding along the substrate. Additionally, wrinkling of the free surface of the cylinder is neglected. In our experiments, such wrinkling may emerge at large rotations, as described in the previous section, but capturing the wrinkling process is beyond the scope of this theory. With these kinematic assumptions, we are left with displacement fields of the form2r = R + u_r(R), θ = Θ + u_θ(R, Z), z = Z + u_z(Z),where the displacement functions are defined as the difference between the deformed (lowercase) and undeformed (uppercase) coordinates. The deformation gradient F = ∂χ/∂X can thus be written as3and incompressibility implies4Integrating the above formula, and eliminating radial translation of the cylinder (i.e. imposing u_r(0) = 0), yields the relations5Substituting the above result into (3), we obtain the final form of the deformation gradient as6where the only remaining unknown field variable is the angular displacement function u_θ.

The problem at hand is concerned with a non-conservative process, whereby torsional deformation induces delamination and frictional sliding. All of these energetic contributions must be considered to determine the propagation of delamination R_d(α). Nonetheless, if R_d and α are separately prescribed, the deformation of the elastic body, u_θ(R, Z; R_d, α) and the corresponding stored elastic energy E_e = E_e(R_d, α), can be derived for a known set of boundary conditions (note the notation used in defining u_θ: the variables before the semicolon are spatial coordinates, while those after the semicolon are prescribed system parameters). After obtaining E_e = E_e(R_d, α), we will also consider the surface energy E_d = E_d(R_d), i.e. the energy invested in delamination, and we will show in the next subsection that the work done by friction can be represented as an effective potential E_f = E_f(R_d). By combining all of these energetic contributions, we will determine the location of the delamination front as the one that minimizes the total energy in the system.

First, we employ the neo-Hookean model to describe the elastic strain energy density of the material as7where μ is the shear modulus. Substituting (6) in (7) gives8The total elastic energy of the system, E_e, is thus given by integration of the strain energy density over the entire volume of the cylinder:9Next, the virtual work W done by the traction force on the top and bottom surfaces, can be written as10where t_θz is the generalized shear traction corresponding to displacement u_θ, but is defined positive in the opposite direction. Here, without loss of generality, we assume an arbitrary externally applied shear traction and account for pre-deformation of the adhered surface, u⁰_θ = u⁰_θ(R). Though initially u⁰_θ = 0; if a ‘stick-slip’ event occurs, bonding of the interface can be re-established imposing a nonzero u⁰_θ(R). This will be explained later in more detail in describing the boundary conditions and solution procedure.

An additional energy contribution E_d comes from the creation of new surface area as the bottom surface of the cylinder delaminates – this contribution can be written as11E_d = πΓ(R₀² − R_d²),where Γ is the interface toughness per unit area.

An equilibrium solution minimizes the potential energy U in the system,12U = E_e + E_d − W.For a prescribed rotation angle α, the potential energy is a functional of u_θ and R_d, i.e. U = U(u_θ, R_d; α).

Using the tools of calculus of variations, we consider perturbations of the stationary function u_θ(R,Z) in the form ũ_θ = u_θ + εδu_θ, where δu_θ(R, Z) is an arbitrary function that must vanish on the boundaries of the body that are subjected to a displacement constraint, and ε is a small constant. Next, we substitute ũ_θ into the potential energy and require that13where we have used the notation Ũ = U(ũ_θ, R_d; α). Additionally, a full derivative in the second equation accounts for the dependence of the displacement field on the location of the delamination front, which may emerge from the boundary condition at Z = −H/2, such that u_θ = u_θ(R, Z; R_d).

By substituting (9) and (10) and performing integration by parts, the first of the above two equations can be recast in the following form14which does not depend on the interfacial energy contribution (E_d). For the above equality to hold for an arbitrary variation δu_θ, each term must vanish separately. The first term gives rise to the governing equation15The second integral implies that on the free boundary of the cylinder, we have16The last integral corresponds to the top and bottom surfaces of the cylinder. At Z = H/2, to enforce uniform rotation about the axis of symmetry, the displacement is prescribed as17Hence, the variation must vanish (δu_θ = 0), and thus the corresponding integral term vanishes. At Z = −H/2, we split the integral to distinguish the different regions of the surface. In the adhered region, the displacement is prescribed:18and thus the corresponding integral term vanishes.

In the delaminated region, sliding is permitted and can be balanced by a shear traction. The variation need not vanish in this region and thus, to ensure that the integral term vanishes, we enforce19The shear traction that emerges from the frictional sliding between the substrate and the cylinder acts as an external force on the system and needs to be prescribed to complete the boundary value problem.

Inspired by the work of Bigoni et al.,^40,41 we now perform the differentiation in (13)₂, to obtain the configurational force associated with propagation of delamination. Such forces are often named after Eshelby.^42,43

It is instructive to note that the derivative in (13)₂ is directly related to the “energy release rate”, which is commonly considered in fracture mechanics. The resulting integral that we write next is thus a specialized version of the J-integral, for the present problem. An illuminating discussion on the analogies between these different formulations for a 1D delamination problem is found in ref. 44. Here, the derivation is complicated by the fact that we are considering a 2D field with combined traction and displacement boundary conditions, as detailed next.

First, it is instructive to examine the differentiability of the displacement function, u_θ(R, Z; R_d). In the bulk of the homogeneous cylinder, we can assume the displacement field is smooth, namely that u_θ ∈ C². However, the transition from a displacement boundary condition (18) to a traction boundary condition (19) permits discontinuity in the slope of u_θ at (R, Z) = (R_d, −H/2), while continuity of u_θ must be preserved, namely20Here, we have introduced the bracket notation, which will be used in this work as a short-hand to denote the jump of a function exclusively at (R, Z) = (R_d, −H/2).

Accordingly, we perform the differentiation in (13)₂ while using the Leibniz integral rule to account for the discontinuity at (R, Z) = (R_d, −H/2), to write‡‡

Note that this formula has been divided by 2π for compactness.

By substituting (8) in (*), performing integration by parts, and making use of the governing eqn (15), we have22By direct integration of the above terms, and making use of the boundary condition (16) and the relation for the boundary traction (19), we find that the first integral cancels and the remaining term reads23Now, by substituting (23) in (21) and noticing that according to (17), ∂u_θ/∂R_d = 0 at Z = H/2, the derivative in (13)₂ can be recast in the form24To determine the Eshelby-like force that emerges from this analysis, we use the principle of virtual work to write a modified potential energy in the form25Û = U − 2πR_d 0000000111100 0000001100010 0000011001010 0000011000100 0000001000000 0000011000000 0001100010000 0111000101000 0110001001000 0110000010000 0010000100000 0001111100000 ,where is an Eshelby-like force, and the second term on the right-hand-side represents the virtual work done by this configurational force in propagating the delamination front along the negative R direction. In equilibrium we require26and thus from (24) we have27As seen from the above result, the Eshelby-like force is determined from the singularity at the delamination front.⁴² Namely, since u_θ − u⁰_θ → 0 as R → R_d the shear traction must become unbounded (i.e. t_θz → ∞). From a numerical point of view this becomes intractable and is thus avoided in this work by directly performing the differentiation in (13)₂, in analogy to the use of the J-integral.⁴⁵

The above derivation considered an arbitrary applied shear traction, t_θz, on the top and bottom boundaries of the cylinder (Z = ±H/2). Here we specialize the formulation to the specific boundary conditions. As a first step, we notice that the virtual work on the displacement controlled boundaries can be neglected. On the remaining boundary, at the limit of quasistatic motion, we can simplify the representation of the kinetic friction by assuming that it is constant throughout the sliding, and proportional to the initially applied normal force N through the constant k, such that t_θz = kN. Accordingly, (19) simplifies to a boundary condition on the slope of the displacement field, in the form28and the virtual work done by friction equals the true work (10) which can be alternatively written as29Since this result is path-independent§§

Note that here we restrict our attention to loading, without considering unloading.

We would like to emphasize that although the present model is limited to quasistatic motion, kinetic friction must be accounted for here if unstable delamination ensues, and will allow us to capture the newly established equilibrium state that follows. This will be explained in more detail in describing the solution procedure.

The above formulation is derived assuming a prescribed rotation angle – α, and deformed height – h = λ_zH. Hence, the corresponding generalized forces, i.e. the torque – T, and the normal force – N (defined positive in compression as shown in Fig. 1, where the compressive force F is equivalent to the quantity N discussed in this section), do not naturally emerge from the formulation. Nonetheless, a direct method to obtain them at a given state (i.e. with given R_d), is by considering perturbations of the elastic energy with respect to the corresponding generalized coordinate, which translates to the partial derivatives31Though these derivatives can only be calculated a posteriori, they will be useful for comparison with the experimental observations.

The above formulation consists of a single second order, partial differential eqn (15) for u_θ = u_θ(R,Z), complemented by conditions on the different regions of the boundary (16)–(18) and (28). This completes our boundary value problem, but for a prescribed location of the delamination front R_d and for a given α. To determine R_d(α), the location of the delamination front as rotation progresses, we must invoke the second requirement in (13). This requirement implies that the system should choose the location of the delamination front to minimize the total potential energy of the system. Though the dependence on R_d is not explicit in U and thus the differentiation cannot be done analytically, identifying the location of the front can be achieved numerically by considering all different candidate values of R_d and finding the one that corresponds to the minimal potential energy. This can alternatively be written in the form32

The delamination process is not gradual. The build-up of elastic energy in the cylinder is necessary to trigger delamination, which is then abrupt, leading to finite sliding at the interface, before arriving at a new equilibrium state for which kinetic friction is balanced by the elastic forces. Then, once at rest, the adhesive bond is re-established, akin to a static friction, and loading can proceed without immediately inducing additional delamination. It is instructive to notice that this sequence of events is analogous to a one-dimensional stick-slip system, with the applied displacement in the one-dimensional system representing the prescribed angular rotation of the top surface; the deformation of the spring in the one-dimensional system representing the torsional deformation of the cylinder; and the resulting uniaxial motion in the one-dimensional system representing the motion of material particles distributed throughout the surface, along circular trajectories.

In each of the steps of the deformation and delamination process described above (i.e. as α increases), the boundary conditions need to be adjusted accordingly in the solution procedure:

Step I – Initial deformation: In this initial stage of the process, no delamination has occurred. Accordingly, the boundary value problem is solved with R_d = R₀ and u_θ⁰ = 0. Nonetheless, at every increment of applied rotation α, solutions with R_d ≤ R₀ are examined, employing the shear traction boundary condition (28). Then, the minimum energy requirement (32) is used to determine the location of the propagation front. If R_d = R₀, the solution procedure remains in Step I; otherwise, it transitions to Step II.

Step II – Snap-through: Once R_d = R₀ no longer provides an energetically favorable solution, a delamination front will propagate, also leading to a drop in the applied torque. Considering a quasistatic process, our model captures the new equilibrium state that the system will choose, which results from the minimum energy requirement (32), as calculated in the previous step. Once the propagation stops, the adhesive bond is re-established and the displacement at the bottom surface is constrained to a new location, namely u_θ⁰ is re-assigned to u_θ⁰(R) = u_θ(R, −H/2).

Step III – Secondary deformation: In this step torsional deformation proceeds with no additional delamination. Hence, the boundary value problem is solved with R_d = R₀ and u_θ⁰(R) from Step II. Similar to Step I, in this stage of the process, no additional delamination occurs. Nonetheless, at every increment of applied rotation α we examine solutions with R_d ≤ R₀ and employ the shear traction boundary condition (28). The minimum energy requirement (32) is again used to determine the location of the propagation front. If R_d = R₀ then the solution procedure remains in Step III, otherwise it transitions back to Step II.

All of the numerical derivations are conducted in Matlab. The governing equation in (15) is integrated using a finite difference scheme.