## Abstract

Biological adhesive contacts are usually of hierarchical structures, such as the clustering of hundreds of sub-micrometre spatulae on keratinous hairs of gecko feet, or the clustering of molecular bonds into focal contacts in cell adhesion. When separating these interfaces, releasable adhesion can be accomplished by asymmetric alignment of the lowest scale discrete bonds (such as the inclined spatula that leads to different peeling force when loading in different directions) or by elastic anisotropy. However, only two-dimensional contact has been analysed for the latter method (Chen & Gao 2007 *J. Mech. Phys. Solids* **55**, 1001–1015 (doi:10.1016/j.jmps.2006.10.008)). Important questions such as the three-dimensional contact morphology, the maximum to minimum pull-off force ratio and the tunability of releasable adhesion cannot be answered. In this work, we developed a three-dimensional cohesive interface model with fictitious viscosity that is capable of simulating the de-adhesion instability and the peripheral morphology before and after the onset of instability. The two-dimensional prediction is found to significantly overestimate the maximum to minimum pull-off force ratio. Based on an interface fracture mechanics analysis, we conclude that (i) the maximum and minimum pull-off forces correspond to the largest and smallest contact stiffness, i.e. ‘stiff-adhere and compliant-release’, (ii) the fracture toughness is sensitive to the crack morphology and the initial contact shape can be designed to attain a significantly higher maximum-to-minimum pull-off force ratio than a circular contact, and (iii) since the adhesion is accomplished by clustering of discrete bonds or called bridged crack in terms of fracture mechanics terminology, the above conclusions can only be achieved when the bridging zone is significantly smaller than the contact size. This adhesion-fracture analogy study leads to mechanistic predictions that can be readily used to design biomimetics and releasable adhesives.

## 1. Introduction

Two clean and atomically smooth surfaces, when being brought together, will have strong adhesive force because of intermolecular interactions such as van der Waals force. Macroscopic surfaces do not retain this property because surface roughness results in a tiny fraction of surface area in true contact, and thus the overall adhesion force is negligibly small. Biological systems have evolved into intriguing surface structures to maximize the adhesion/de-adhesion forces by hierarchical structures. For instance, some animals and insects have keratinous hairs, called setae, on their feet, which then have hundreds of small hairs, called spatulae, of sub-micrometre diameters [1,2]. The small size in each hair enables a nearly uniform stress distribution during de-adhesion, thus avoiding crack-like behaviour that usually corresponds to low pull-off force. Also the long hairs can easily bend to conform to the surrounding rough surface and thus to increase the fraction of true contact area. For another example, cell adhesion is accomplished by micrometre-sized focal contacts, each of which contains the clustering of a larger number of molecular receptor–ligand bonds. The adhesion and de-adhesion of these focal contacts lead to elastic deformation of the surrounding material, which is termed mechanosensitivity [3–6]. There are also many other examples in which the adhesive contacts comprise discrete ‘bonds' in biological systems such as placental villi in embryo implantation, and in artificial systems such as micromachined pillar arrays or microfludics with subsurface microstructures [7], and in geological systems such as multi-asperity contact [8].

A central problem in adhesion and de-adhesion of biological surface structures is how these long-chain discrete bonds (being spatulae or receptor–ligand bonds) are broken. Would these bonds be broken simultaneously, or the debonding starts somewhere and propagates into a certain spatial pattern? How do these discrete bonds feel the neighbouring bonds and the faraway applied load? We note that the de-adhesion process from a finite contact area is equivalent to an interface fracture problem, where the crack is external and can be of various shapes. In the limit of Griffith crack, the stress near the contact edge diverges to infinity, and the fracture occurs when the applied energy release rate (that is, the energy drop with a unit length increase of the crack front) reaches the interface fracture energy [9]. This scenario applies to breaking brittle materials such as glass and ceramics. When the interface is made of long-chain molecules, the low interface stiffness may significantly reduce the stress singularity at the contact edge or crack front, with the limit case being that the entire interface has a uniform stress that equals the interface theoretical strength [10–12]. This scenario applies to tearing a Velcro patch or pulling off gecko toes that comprise a huge number of hairs. The transition from the Griffith limit to the uniform stress limit is governed by the ratio of the crack bridging zone size to the contact radius [10–12], which is also called the stress concentration index [13]. Therefore, how the long-chain bonds break spatially is closely connected to the above transition, which has important consequences in designing surface structures to achieve optimized adhesion force.

Another central problem in bio-adhesive contact is the orientation-dependent adhesion strength, so that releasable adhesion can be easily achieved when loading in various directions. In the hierarchical structure, this is accomplished by the alignment of discrete bonds (such as seta–spatula structure in gecko toes) in inclined directions, so that the pull-off force varies significantly with respect to the pulling direction [2]. This is like peeling a Velcro or Scotch tape at different attack angles from the substrate. Besides this structural point of view, the orientation dependence of de-adhesion force can be obtained by the elastic anisotropy of the substrate [14,15]. In the context of fracture energy, the energy release rate is inversely proportional to the moduli in the plane spanned by the crack growth direction and the crack plane normal. This line of thought has been used in natural systems such as the attachment pad of cicada and grasshopper which have anisotropic microstructure and thus elastic anisotropy. A plane-strain elastic analysis of circular adhesive contact with a substrate of transverse isotropy found that the maximum pull-off force can be achieved if pulling in the direction normal to the plane of transverse isotropy [2,14]. However, there are two severe limitations of these works. First, their solutions are based on the so-called Johnson–Kendall–Roberts adhesive contact, which corresponds to the Griffith crack limit [16]. Since many bio-adhesive contacts will be close to the uniform stress limit, it is unclear whether the degree of orientation dependence of pull-off forces in the above calculations still remains large. Second, both elastic anisotropy and actual bio-adhesive contacts are of three-dimensional nature. The pull-off force is closely connected to the evolution of crack morphology which in turn depends on elastic anisotropy and initial contact zone shape. For example, the long-chain molecular bonds may break non-uniformly in space so that an initial circular contact may evolve into an elliptical-shaped crack. The synergy of all these factors can certainly lead to a releasable adhesion that is very different from predictions in [14].

In this work, we develop a cohesive interface model that describes the interface fracture process of the long-chain molecules in bio-adhesive contacts. It should be noted that the Griffith limit corresponds to a sudden release of applied energy and a snap-back instability will occur, while no instability occurs for the uniform stress limit. To examine the onset of instability and simulate the post-stability evolution, we have adopted a fictitious viscosity approach in the numerical implementation of the cohesive interface model into a nonlinear finite-element framework [12], as will be discussed in §2. A parametric study has been conducted to study (i) the transition from Griffith to uniform stress limits as governed by the ratio of crack bridging zone to contact radius, (ii) the degree of elastic anisotropy on the ratio of maximum-to-minimum pull-off forces, and (iii) the evolution of the crack front or peripheral morphology with respect to the crack bridging characteristics and the degree of elastic anisotropy. Representative examples include the twofold symmetric surface of Cu (110) in §3 and the transverse isotropic material in §4. Limitations of the two-dimensional results in [14] have been identified. Asymptotic limits of interface fracture and contact stiffness analysis discussed in §5, so as to develop design rules that can be used to tune the orientation dependence of pull-off forces. Conclusions are given in §6.

## 2. Problem definition and interface model

As shown by the boundary value problem in figure 1, we aim to calculate the pull-off force that separates a finite contact (circular, elliptical or arbitrary shape) that is connected to an elastically anisotropic substrate by discrete bonds. Rather than treating these bonds individually, we assume there is a uniform distribution of these bonds and the interface constitutive law can be described by an interface traction–separation relationship in figure 2. The fracture mechanics analysis dictates that the de-adhesion force depends on the fracture energy and the fracture strength of the interface, but not on the detailed shape of the traction–separation law [12,17]. The use of bilinear model in figure 2 applies well for long range forces such as van der Waals interaction, or nonlinear springs such as spatulae or micromachined pillar arrays. As pointed out by Freund & Lin [18], the bio-adhesive contacts in cells are realized by mobile receptor–ligand bonds on the interface, but not by any kind of chemical signalling process. Therefore, cell adhesion is primarily governed by thermodynamic and mechanical driving forces. The receptor–ligand bonding energy is about 1–10 times *k*_{B}*T* with *k*_{B} being the Boltzmann constant and *T* being absolute temperature. Therefore, one immediate consequence of thermodynamic fluctuation is the stochastic behaviour of the receptor–ligand molecular bonds [19,20]. The lifetime of an individual bond is of the order of 100 s, while the de-adhesion process is much faster. Consequently, a mean-field stochastic elastic model of dissociation and association of receptor-ligand bond [13] will be adopted. Details of this model will not be pursued here, since essential features can also be well captured in the traction–separation law in figure 2. It should be noted that the same constitutive parameters are used across the entire interface, which implies that we do not consider a spatial variation of the bond density. In reality, the bond density can be higher at the contact boundary because of entropic effects. Such an observation can be easily incorporated in our model by assigning different constitutive parameters across the interface, which however is not pursued in this work for the sake of clarity [21,22]. We also note that the adhesive interaction can also be modelled in a body force formulation [23], which will not be pursued in this work.

The bilinear model in figure 2 has three governing parameters: the interface strength *σ*_{max}, the characteristic length scale *δ*_{max} and the decay length *δ _{c}*. The total separation can be projected into the normal and two tangential directions in local coordinates. When Δ

*< 0, the two contacting surfaces are not penetrable, so only tangential separations contribute to the total separation. Following [24], we write the relationship as 2.1and 2.2where the total separation is 2.3In the above model,*

_{n}*k*is taken as orders of magnitude higher than

_{n}*k*. When implementing the above relationship into a nonlinear finite-element framework, we need to use the material tangent, given by

_{c}*J*=

_{kl}*∂T*/

_{k}*∂*Δ

*, where subscripts*

_{l}*k*and

*l*run as (

*n*,

*t*

_{1},

*t*

_{2}). A fictitious viscosity,

*ζ*, is introduced in order to successfully simulate the snap-back instability that is intrinsic to Griffith-type cracks [12,24] and adhesive contacts at the Griffith limit [16].

In the geometric set-up in figure 1, we differentiate the applied force or displacement conditions, which are not necessarily in the same direction because of the elastic anisotropy. Note that our finite-element method is displacement based. We consider a circular or elliptical surface contact on a half space, and the boundary extends 20 times the contact radius. The normalized fictitious viscosity is optimized to be ; a smaller value than this will lead to smaller time steps and longer calculation time, and a larger value will lead to excessive stress due to the viscous force term in equations (2.1) and (2.2). Particularly, this choice of normalized fictitious viscosity only introduces negligible effect on the calculated pull-off or adhesion force [12]. The shape of the interface traction–separation law is fixed, e.g. *δ*_{max}/*δ _{c}* = 0.5, and the interface fracture energy, , is kept constant in our simulations so that the variations of

*σ*

_{max}and

*δ*

_{max}will be inversely proportional to each other. A total number of about 50 000 three-dimensional eight-node elements are used. The cohesive interface model is implemented into a commercial finite-element software, ABAQUS, by the user defined element subroutine.

Our numerical simulations will be specifically interested in the following questions. First, the pull-off force depends on the loading direction because of the elastic anisotropy. We will explore representative anisotropic materials including cubic crystal and transversely isotropic material. One question that naturally arises is in which loading directions we can find the pull-off force extrema. Second, during de-adhesion, the contact periphery will evolve into a shape dictated by the elastic anisotropy, which can be monitored by finding the contours of Δ = *δ*_{max}. It is particularly interesting to find the relationship between the evolution of peripheral morphology and the pull-off forces.

## 3. Anisotropic elastic substrate with cubic symmetry

Representative load–displacement curves are plotted in figure 3 for Cu (110) surface with respect to various loading directions and interface properties. The elastic constants are given table 1 as Material (A), where the Voigt contracted form, *c _{ij}*, is used. The twofold surface of this material exhibits a significant elastic anisotropy. Note that the elastic constants of typical biological materials are of the order of kPa–GPa, as opposed to approximately 100 GPa in metals and ceramics. However, the purpose of our study is merely to illustrate the role of the degree of elastic anisotropy, and the absolute value of elastic constants will be embedded in the normalization factor.

As the first observation, after extensive numerical tests, we found that an applied tangential displacement along the [110] direction gives the minimum pull-off force, and an applied normal displacement in the direction corresponds to the maximum pull-off force. To explain this observation, we need to refer to the analytical solutions of the contact problem [16]. Regardless of the indenter shape, a circular contact under shear loading is equivalent to an external circular crack under shear force, with the interface shear stress distribution given by
3.1where *r* is the radial coordinate, . At the contact edge, the stress singularity is the same as that in a crack. The mode II (in-plane shear) and mode III (anti-plane shear) stress intensity factors are
3.2where . The non-adhesive normal contact depends on the indenter shape, but a normal pull-off force applied on an adhesive contact is of the same form as the flat-ended punch contact problem. We have the normal stress distribution and the mode I (in-plane tension) stress intensity factor
3.3The above solutions apply to any arbitrarily anisotropic solid as long as the contact shape is circular. The de-adhesion process is thus equivalent to a fracture problem with mixed modes as in equations (3.1)–(3.3).

In the linear fracture mechanics (i.e. the Griffith limit of the crack characteristics), fracture occurs when the applied energy release rate, *G*, reaches the interface fracture energy, *Γ*. Although the stress intensity factors are purely governed by the applied load, crack size and crack shape, the energy release rate depends additionally on the elastic constants. For instance, the isotropic elasticity gives
3.4where *E* and *ν* are Young's modulus and Poisson's ratio, respectively. Substituting equations (3.2) and (3.3) into (3.4), one can easily prove that for an isotropic elastic solid subjected to a constant force with varying loading directions, the tangential load will give the highest *G* and the normal load gives the lowest *G*. Therefore, de-adhesion is the most difficult in normal direction, and the easiest in tangential direction. Our simulations for isotropic elasticity (not shown here for brevity) have confirmed this conclusion, in which the pull-off force depends on *β* and *ν*. For an anisotropic elastic solid, the relationship between *G* and *K* is quite complicated, which is out of the scope of this work. Considering the limit cases of normal and tangential loading conditions, an elementary tensor transformation shows that the tensile modulus in normal direction (which enters in the relationship of *G* and *K*_{I}) is higher than other directions, and the shear modulus in (*x*_{2}, *x*_{3}) plane is higher than that in (*x*_{1}, *x*_{3}) plane (which enter in the relationship of *G* and *K*_{II}, *K*_{III}). It is thus anticipated that tangential loading in *x*_{1}-direction gives the largest *G* that corresponds to the lowest pull-off force *F*_{min}, and the normal loading gives the lowest *G* that corresponds to the largest pull-off force *F*_{max}. However, the contact deformation fields are not simply directional or shear, and a more detailed analysis based on contact stiffness will be given in §5.

The second important observation is the role played by the crack bridging characteristics. The stress singularity in the linear elastic fracture mechanics will be limited by the interface strength. The relative significance of the surrounding elastic stiffness and the interface stiffness defines a length scale, *c*_{11}*δ*_{max}/*σ*_{max}, which determines the size of the crack bridging zone. If , i.e. the small-scale bridging (SSB) behavior, the crack approaches the Griffith crack. The large-scale bridging (LSB) behaviour, when *c*_{11}*δ*_{max}/*σ*_{max} approaches or is larger than the crack size, approaches the uniform stress limit. The releasable adhesion, as characterized by the ratio of *F*_{max}*/F*_{min}, varies significantly with respect to the ratio of crack bridging zone size to contact radius, *c*_{11}*δ*_{max}/*σ*_{max}*a*. As shown in figure 4*a*, the SSB limit at corresponds to the Griffith crack. From equations (3.1)–(3.4), one can derive the following dimensionless relationship:
3.5where the dimensionless pre-factor, *γ*, depends on the elastic anisotropy and loading direction. In the LSB limit, the pull-off force simply approaches the product of uniform stress and contact area
3.6which is independent of elastic anisotropy and loading direction. Obviously, the degree of releasable adhesion, *F*_{max}/*F*_{min}, only deviates from unity when the controlling parameter *c*_{11}*δ*_{max}/*σ*_{max}*a* is far away from the LSB limit, as shown in figure 4*b*.

The transition from SSB to LSB correlates closely to the crack shape, as the third important observations from our numerical simulations. As shown in figure 5, the peripheral morphology, or the crack front shape, is illustrated by the contours of crack opening displacement, Δ/*δ*_{max}. The left column is close to the LSB asymptote, so that the crack bridging zone is large and the crack shape exhibits a strong dependence on loading direction and elastic anisotropy. That is, the tangential loading in figure 5*a* leads to an elliptical shape with major axis in the compliant direction in shear, the inclined loading in figure 5*b* leads to a rotated elliptical shape (but the major axis does not necessarily align with the loading direction), and the normal loading in figure 5*c* gives rise to an elliptical shape with major axis in the compliant direction in tension. The right column is close to the SSB asymptote, so that the crack bridging zone is narrow but the degree of eccentricity and rotation of the elliptical crack morphology is the same as that in the left column. In all these contour plots, the limits are set as Δ*/*Δ*δ*_{max} ∈ [0.9, 1.1] for better visualization.

## 4. Anisotropic elastic substrate with transverse isotropy

Many biological materials are transversely isotropic [5,6]. Following the notations in Bower [25] as shown in figure 6*a*, the elastic constants include the directional modulus in *p* and *t* directions, *E _{p}* and

*E*, in-plane and out-of-plane Poisson's ratios,

_{t}*ν*and

_{p}*ν*, and shear modulus

_{tp}*μ*. The degree of elastic anisotropy can be characterized primarily by the following ratio [14]: 4.1

_{t}A representative set of parameters is given in table 1 as Material (B) [6]. Keeping all other elastic constants the same but varying *E _{t}*, we get a fictitious Material (C) with a very large degree of elastic anisotropy,

*L*

_{22}/

*L*

_{11}= 440.

To ensure a two-dimensional contact problem between a cylinder and the transversely isotropic material, the cylinder axis lies along the *x*_{1}-direction in figure 6*a*, but the plane of transverse isotropy can rotate about *x*_{1}. Chen & Gao [14] solved the corresponding adhesive contact problem at the Griffith crack limit and determined the ratio of *F*_{max}*/F*_{min} when loading in various directions in the (*x*_{2}, *x*_{3}) plane. It is found that *F*_{max}*/F*_{min} is about 4 when *L*_{22}*/L*_{11} = 100 and about 12 when *L*_{22}*/L*_{11} = 1000, and the scaling relation can be approximated by *F*_{max}*/F*_{min} ∼ ln(*L*_{22}*/L*_{11}). However, in three-dimensional contacts, because of the finite size in the *x*_{1}-direction and the possibility of loading in arbitrary direction not in the (*x*_{2}, *x*_{3}) plane, their conclusions need to be revisited.

For a circular contact loaded on the plane, our extensive simulations find that the maximum pull-off force corresponds to tangential loading in the *t*-direction and the minimum one to tangential loading in the *p*-direction in figure 6*a*. As shown in figure 7, the dependence of *F*_{max}*/F*_{min} on that describes the cracking bridging behaviour is exactly the same as that for Cu (110) surface. For the Material (C), the transition from LSB to SSB occurs at a very low ratio of , suggesting that the degree of elastic anisotropy actually plays a minor role in a wide range of cohesive zone size.

The ratio of *F*_{max}*/F*_{min} from a three-dimensional contact analysis is much smaller than the prediction from [14]. We now examine if the two-dimensional limit solution can be approached if the initial contact shape is elongated as opposed to circular. This line of thought is motivated by the crack morphology study in figure 8. In calculations of figures 7 and 8, *L*_{22}/*L*_{11} = 440 so that . Consequently, when loaded in *t*-direction (here a tangential load in *x*_{2}-direction), the crack will initiate at *θ* = 0° and 180°, and fracture will be extremely difficult at *θ* = ±90° because the extremely large *E _{t}* leads to a very small energy release rate. For tangential loading along

*x*

_{1}-direction, the shear modulus anisotropy leads to a lens-shaped crack morphology in figure 8

*a*. In order to maximize (or minimize) the pull-off force, we can ensure a large fraction of the contact edge lying in the stiff (or compliant) direction. For instance, the right column plot in figure 8

*a*can be changed to an elliptical contact (with the same contact area and the major axis lying in the compliant direction), so that the uncracked fraction will be maximized at the peak pull-off force. This leads to a much larger pull-off force in the right column plot in figure 8

*b*than that in the right column plot in figure 8

*a*. Similarly, the left column plot in figure 8

*b*increases the fraction of the contact edge with low shear modulus, which leads to much lower pull-off force than the left column plot in figure 8

*a*. It should be noted that the contour plots in figure 8 have different bounds of Δ/

*δ*

_{max}for better illustration of the crack fronts.

Our results find that the elongated shape with *a*_{major}/*a*_{minor} = 4 increases the ratio of *F*_{max}*/F*_{min} by 60%, as shown in figure 7*b*. However, the limit values at are still far less than the two-dimensional results in [14].

## 5. Contact stiffness governs the directions of *F*_{max} and *F*_{min}

In the above sections, the effects of elastic anisotropy are discussed in terms of the directional modulus and shear modulus, which will enter the energy release rate representation (e.g. equation (3.4)). The detailed relationship between *G* and *K* is complicated for anisotropic solids and varies significantly with respect to the crack orientations [9]. The contact problem by itself also leads to complicated stress fields that combine both tensile/compressive and shear components. All the above simulation results and the contact-fracture analogy suggest the validity of the ‘stiff-adhere and compliant-release’ concept [2], in which the orientation-dependent adhesion maximizes in the stiffest direction—for strong bonding, and minimizes in the most compliant direction—for easy de-adhesion. The contact stiffness analysis has been given in [26] by
5.1where is the contact area, *G*_{eff,x1} and *G*_{eff,x2} are effective shear moduli, and *E*_{eff} is the effective normal stiffness. These effective moduli depend on the elastic constants and the orientations of the three axes and can be calculated from a line integral along the contact periphery [26]. Referring to figure 1 where the direction of the applied displacement is described by *α _{U}* and

*β*, we have 5.2Therefore, finding the pull-off force extrema is equivalent to finding the extrema of stretching stiffness, d

_{U}*F/*d

*U*, as shown in figure 9.

For isotropic material, and , and the stretching stiffness is given in the polar plot in figure 9*a*. The normal contact gives the highest contact stiffness so that the pull-off force maximizes, while the minimum pull-off force corresponds to the tangential loading. For Cu (110) surface in figure 9*b*, the stiffest direction is normal contact and the most compliant one is shear in [110] direction. For the two transversely isotropic materials, the stiffest direction is always shear in *t*-direction and the most compliant direction is shear in *p*-direction. Note that these contact stiffness calculations in figure 9 are all for circular contacts. One can optimize the shape of the contact so as to obtain a pull-off range [*F*_{min}, *F*_{max}] that is much wider than that of the circular contact.

## 6. Conclusion

Interface fracture mechanics and cohesive interface simulations have been conducted to analyse the non-uniform breaking of molecular bonds in bio-adhesive contacts. Motivated by a recent two-dimensional analytical solution of pulling off circular contacts from a transverse isotropic material, we systematically investigate the roles played by three-dimensional contact shape, the degree of elastic anisotropy and the crack bridging characteristics.

— Since the adhesion is accomplished by clustering of discrete bonds or so-called bridged crack, the de-adhesion process is controlled by the transition from SSB behaviour (i.e. the Griffith crack limit) to LSB behaviour (i.e. the uniform stress limit), as the dimensionless parameter

*c*_{11}*δ*_{max}/*σ*_{max}*a*increases.— The maximum and minimum pull-off forces correspond to the largest and smallest contact stiffness, i.e. ‘stiff-adhere and compliant-release’, which confirms the previous two-dimensional investigations in Yao & Gao [2].

— The pull-off extrema are closely related to the crack morphology. That is, the compliant direction will fracture more easily, and the resulting crack morphology corresponds to, but inversely, the contact stiffness contours.

— The two-dimensional analysis in Chen & Gao [14] significantly overpredicts the pull-off force extrema, due to the restricted contact shape into an infinitely long strip and due to the mere consideration of the Griffith limit.

— These predictions shed interesting light on the design of biomimetics and releasable adhesives. The maximum (or minimum) pull-off force can be optimized by searching the largest (or lowest) contact stiffness, and by arranging the initial contact shape so that a large fraction of the contact periphery lies in the stiffest (or most compliant) directions.

## Funding statement

Y.L. acknowledges the support from Tianjin Bureau of Public Health (project 2012KZ025), as well as a visiting scholar position provided by the University of Tennessee during the course of this study. Y.G. is grateful for the financial support from NSF CMMI 0900027 and CMMI 1300223.

- Received September 17, 2014.
- Accepted October 20, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.