Topographic de-adhesion in the viscoelastic limit

The superiority of many natural surfaces at resisting soft, sticky biofoulants have inspired the integration of dynamic topography with mechanical instability to promote self-cleaning artificial surfaces. The physics behind this novel mechanism is currently limited to elastic biofoulants where surface energy, bending stiffness and topographical wavelength are key factors. However, the viscoelastic nature of many biofoulants causes a complex interplay between these factors with time-dependent characteristics such as material softening and loading rate. Here, we enrich the current elastic theory of topographic de-adhesion using analytical and finite-element models to elucidate the nonlinear, time-dependent interaction of three physical, dimensionless parameters: biofoulant’s stiffness reduction, the product of relaxation time and loading rate, and the critical strain for short-term elastic de-adhesion. Theoretical predictions, in good agreement with numerical simulations, provide insight into tuning these control parameters to optimize surface renewal via topographic de-adhesion in the viscoelastic regime.


Introduction
Many natural surfaces including airways, arteries and intestines have dynamically actuated geometries [1][2][3][4][5][6][7][8][9][10][11]. Specifically, pulse pressure in arteries is believed to drive the luminal arterial geometry between low and high curvature states, generating an actuating topography [8][9][10][11]. Figure 1 schematically demonstrates the change in arterial topography during the cardiac cycle, with a flat luminal surface at systole (figure 1b, top) and a wrinkled surface at diastole (figure 1b, bottom). The existence of these states has been experimentally validated ex vivo in fresh arterial segments in a previous work [12]. Such topographic variations may play an important role in keeping these surfaces clean from biofouling on macroscopic scales and potentiate nature's multi-scale anti-fouling strategies at complex interfaces. Inspired by these natural phenomena, a new de-adhesion mechanism, using dynamic surface topography, was discovered: topography-driven delamination (figure 1), where foulant deformation, driven by evolving surface curvature, generates an energy release mechanism by balancing elastic energy with adhesion strength [6][7][8][9][10]. Thus far, topography-driven delamination has only been explored in the idealized elastic foulants, where all stored elastic energy is available to drive surface renewal once a critical surface curvature, κ c , is reached [8].
Biofoulants, such as platelets, thrombus, biofilms and bacteria, are complex, dynamic materials. In order to integrate topographic anti-fouling strategies into biomaterials and medical devices with specific engineering design limits with given loads and achievable actuation strains, topography-driven delamination must be understood in the viscoelastic regime that more accurately represents the material properties of biological foulants [2,8,[13][14][15][16][17]. Viscoelasticity allows a continuous softening of material properties, a dependence on loading history and an intrinsic mechanism for dissipating elastic energy within a stressed material [18,19]. Since topographydriven delamination is an energy release mechanism where available foulant elastic energy drives interfacial fracture, the existence of a competing dissipation mechanism, intrinsic to the material, enriches the problem substantially, compared with the elastic limit.
A model system with a thin foulant on a bilayer composed of a thin elastic film on a soft, thick elastic substrate is adopted from the previous work [8]. The system is subjected to continuous loading on the two ends with the deformation profile e ¼ _ et, where t is the loading time. If the foulant behaves elastically, details of the delamination process have been presented in this previous work. In our Results section, we will provide detail calculations for the viscoelastic foulant. Figure 2 shows a striking difference between the detachment process of an elastic (or equivalently, a viscoelastic foulant layer with a large intrinsic relaxation time τ R relative to the intrinsic loading rate _ e) and a viscoelastic foulant layer with a short relaxation time scale (see also electronic supplementary material, video S1). In both cases, the foulant layers initially follow the wrinkled surface, which is the source of changing curvature κ = A/λ 2 , where A is wrinkle amplitude and λ is wrinkle wavelength. They ultimately delaminate; however, the one with fast intrinsic relaxation, _ et R ( 1, remains attached until higher critical curvatures (critical strains). These examples show that the presence of significant viscoelasticity, on the time scale of loading, stabilizes the foulant layer/substrate interface. Of note, the nominal compressive critical strain more than doubles from approximately 5% to 12% in these two cases.
For optimal design of self-cleaning surfaces in contact with viscoelastic biofoulants, the interaction of the timedependent characteristics of the viscoelastic foulant layer with the geometry set by the wrinkled topography and the surface energy needs to be understood, especially for dynamic topography-driven anti-fouling strategies in biomaterials and medical devices [6,7,9,10]. In this paper, we tackle the complex interaction between foulant layer's viscoelasticity, rate of topographic changes, mechanical instability and surface adhesion in the wrinkle-induced delamination model [8]. Using analytical methods based on energy minimization [8,18,[20][21][22][23][24][25][26][27][28] and finite-element analysis with cohesive zone modelling [8,25,[29][30][31][32][33][34][35][36][37], we reduce the parameter space that controls the delamination of a thin, viscoelastic foulant layer from a wrinkled surface to three physical dimensionless parameters: magnitude of foulant layer's relaxation, rate of material relaxation relative to loading rate and the critical strain for delamination in the short-term (instantaneous elastic) limit. The last parameter incorporates the surface energy, the bending stiffness of the foulant layer and the surface geometry via the wrinkle wavelength. Our analytical model for viscoelastic topographic deadhesion is able to collapse numerical data from a large range of the dimensionless parameter space. Our analysis provides insight into viscoelastic delamination for geometrically nonlinear interfaces that inevitably exist in biological systems and for artificial materials incorporating topography as self-cleaning strategies [6-10].

Scaling analysis
A viscoelastic foulant layer (adherent biofoulant) of thickness h is attached to the surface of an elastic bilayer designed to mimic the multi-layered structure of the arterial wall [8,12]. The resultant tri-layer system is subjected to compression as shown in figure 3a. Under increasing applied compression on the two ends and assuming plane strain in the orthogonal z-direction, the bilayer composed of a stiff film of thickness h f attached to a soft substrate of thickness h s ≫ h f , undergoes wrinkling with increasing wrinkle amplitudes as described previously [3,8,[38][39][40][41][42][43]. Specifically, the critical strain for wrinkle onset is set by the mismatch stiffness between the film and the substrate: e w ðE s =E f Þ 2=3 where E s and E f are the The SEM images show that the arterial wrinkle surface is covered with platelets and thromboses at an early stage. The schematic drawing of the artery provides a continuum model for their interaction that has been proposed in previous works [8][9][10] at the stage where biofoulants grow and expand over multiple wrinkle wavelengths. (b) Cross-sectional views of the postulated detachment process of biofoulants from the wrinkled topography at different levels of actuated pulse pressure through the cardiac cycle illustrated in black on the right. During systole (top), the high pressure distends the artery and flattens the luminal wrinkling. As the cardiac cycle progresses, the arterial pressure drops gradually until the diastolic pressure (bottom) is reached. As the pressure decreases, the artery contracts and the amplitude of the wrinkles grows. The cardiac cycle then repeats when the left ventricle again contracts pumping blood into the arterial system. Current attempts with experiment and computational modelling focus on elastic behaviours of the biofoulant and show a scaling dependence of the critical surface curvature on the surface energy, bending stiffness and topographic geometrical wavelength [8].
substrate and the film moduli, respectively, and E f ≫ E s . Surface topography is set by its curvature κ = A/λ 2 , which depends on wrinkle wavelength , where e is the applied nominal strain. This wrinkle pattern is not affected by the foulant layer because the foulant layer is significantly softer than the constituents of the bilayer [8]. To reduce the effect of the pre-wrinkled state on the subsequent delamination process, we study the regime where the wrinkle onset strain e w is small, thereby A l ffiffi ffi e p . As a critical curvature k c A c =l 2 ffiffiffiffi e c p =l is reached, the viscoelastic foulant layer starts to de-adhere from the bilayer surface. In the elastic limit for the foulant layer solved by Pocivavsek et al. [8], no time-dependent quantity enters the solution e c ¼ e de ¼ al 2 G=B, where G is the adhesion energy, B = Eh 3 / 12(1 − ν 2 ) is the bending stiffness of the elastic foulant layer (ν is the Poisson's ratio of the foulant layer) and a is a numerical pre-factor. In the current case of viscoelastic foulant layer, the intrinsic relaxation time of the material interacts with the time scale set by the rate of topographic loading to give rise to a time-dependent delamination condition [44]. Here, for the simplest viscoelastic foulant layer, the relaxation is described using a single term Prony series: where τ R is the relaxation time and E ∞ ≤ E 0 are the long-term and instantaneous elastic stiffnesses of the foulant layer, respectively. A mechanical analogue for this representation is a Maxwell model with two springs of stiffnesses E ∞ and E 0 − E ∞ and a dashpot to dissipate energy (see electronic supplementary material, appendix 1). With this decaying function of the modulus, a simple mathematical extension of the result of Pocivavsek et al. [8] for compression applied at a constant rate e ¼ _ et is to substitute E(t) into the scaling law derived in the elastic limit, which leads to the identity : ð2:1Þ In terms of time scales, this equation is written as bending stiffnesses of the foulant layer at the instantaneous and long-term elastic limits, respectively; and t de ¼ e de =_ e is the characteristic time for delamination in the instantaneous elastic limit.
Three physical dimensionless parameters emerge from these equations: E ∞ /E 0 (magnitude of the foulant layer's relaxation), _ et R (rate of material relaxation relative to loading rate, which is similar to the Weissenberg number used by rheologists to quantify viscoelastic effects [45]) and e de (critical delamination strain at the instantaneous elastic limit, which depends on wrinkle wavelength, adhesion strength and the instantaneous elastic modulus of the foulant layer as quantitatively described by the elastic model of Pocivavsek et al. [8]). The ratio between the last two parameters τ de /τ R is a control parameter for the system: when τ de /τ R ≪ 1 we expect from the balance in equation (2.2) that t c /τ R ≪ 1. It implies that the denominator is approximately equal to one such that t c ≈ τ de . Similarly, when τ de / τ R ≫ 1 the same balance gives t c /τ R ≫ 1 and the result t c % t de =ð1 À bÞ ¼ al 2 G=ð_ eB 1 Þ. Thus, delamination is controlled by the bending stiffness B 0 for τ de /τ R ≪ 1 and the opposite limit corresponds to delamination controlled by B ∞ .

Energy analysis
Though the above mathematical approach indicates the presence of important dimensionless parameters controlling this complex viscoelastic system and their nonlinear coupling, to obtain a physical mechanism driving delamination in this system, we employ an analysis based on energy minimization. In the case of a purely elastic foulant layer without delamination, the work performed by the wrinkled surface generates stored elastic energy in the foulant layer _ W ¼ _ U, where the total work injected into the system is _ W ¼ Ð V s_ e dV and s ¼ Ee. In the case of a viscoelastic foulant layer, part of the total work is stored and the other part is dissipated through viscous relaxation: Here again _ W ¼ Ð V s_ e dV but the stress and strain are divided in two parts: s ¼ E 1 ðe e þ e v Þ þ ðE 0 À E 1 Þe e and e ¼ e e þ e v , with e e being the elastic strain and e v being the viscous strain. Using a mechanical analogue of linear viscoelasticity (see electronic supplementary material, appendix 1), the stored elastic energy of the foulant layer is derived as royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220598 Note that for t ≪ τ R , no viscous dissipation has yet occurred as no viscous strain has had time to develop, thus D = 0, On the other hand, for t ≫ τ R , the material has exhausted all possible sources of viscous dissipation, thus the foulant layer's elastic strain e e ¼ 0 and Since the total strain e ¼ e e þ e v , we see in the two limits the available energy is simply that of the purely elastic case with the two moduli E = E 0 and E = E ∞ in the single term Prony series. In other words, a physical interpretation of the stored energy is the elastic energy stored in the two elastic springs of the Maxwell model [18]. At the instantaneous limit t ≪ τ R , U and W reduce to the bending energy for the elastic case with modulus E 0 . At the long-term limit t ≫ τ R , U and W have the same form of the elastic bending energy but with modulus E ∞ . However, in the intermediate regime between these two limits, the difference between U and W arises due to active viscous dissipation in this system (see detailed derivations in electronic supplementary material, appendix 2). It is in this active regime, where two time-dependent mechanisms interact, that the analysis becomes interesting and complex.
To drive fracture, a balance between the energy available in the system and the surface energy is required. For a conservative system, dW = dU and a balance between the stored energy and the surface energy d(U + U S ) = d(W + U S ) = 0, where U S = Gl and l is the crack length, is widely used to compute variations of the crack length to study fracture under displacement controlled conditions [8,[23][24][25][26]46]. However, for a dissipative system, due to the presence of the material's dissipation energy, two arguments are presented in the literature for the source of energy to predict fracture. One approach balances the total work W against the surface energy as a condition for crack propagation [18,27]: d(W + U S ) = 0. A second approach assumes that only the stored energy U is the amount of energy available for fracture [28]: d(U + U S ) = 0. Note that both approaches are equivalent for a conservative system. The correct energy is debatable, and more accurate experimental data and understanding of other possible sources of energies that might play a role in real systems such as dynamic effects are needed to resolve this conflict. Yet, these two approaches provide insights into two limiting cases for the energy available for driving fracture in a non-conservative, viscoelastic system. In particular, the first one considers the largest possible amount of energy while the second one considers the smallest possible amount of energy available to overcome surface energy and cause fracture to occur. Thus, in order to physically understand the emergence and coupling of the three dimensionless parameters and how they govern the delamination of the viscoelastic foulant layer from the wrinkled surface, we use both approaches in our analytical models and compare them with the results obtained from numerical simulations.
Specifically, from the imposed wrinkled topography, the height of the mid-plane of the foulant layer is h m (x,t) = A(t)sin(kx), where t ¼ e=_ e and k = 2π/λ is the wavenumber. The curvature is computed as k ¼ @ 2 x h m ðx, tÞ ¼ Àk 2 AðtÞ sinðkxÞ giving the bending strain γ = κy at a point located at 0 ≤ x ≤ L along the length L and at a distance of y to the neutral axis (the mid-plane) of the thin foulant layer. Thus, the strain in the foulant layer is computed as γ = −yk 2 A(t)sin(kx) (see electronic supplementary material, appendix 2). Neglecting the small compressive strain prior to buckling, the viscous strain in the foulant layer becomes g v ¼ ð1=t R Þ Ð t 0 gðtÞ e ÀðtÀtÞ=tR dt and the elastic strain is γ e = γ − γ v . Both approaches lead to nonlinear, time-dependent equations whose solutions give the critical time (equivalently critical strain) to trigger delamination Both functions P ± are such that 0 < P ± < 1 and 10 -4 10 -2 1 10 -2 10 -1 1 increasing total work stored energy  [8] is reproduced as shown by the horizontal, flat curves. In these cases, the critical strain is equal to e de and is independent of _ et R .
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220598 for 0 < β < 1 (see electronic supplementary material, appendix 2). Thus, the critical strain for delamination is always larger for the stored energy than the total work. This agrees with physical intuition, since in the case of the stored energy, only the elastic energy at any given time is available to drive fracture. However, in the case of total work, some of the dissipated energy may have gone into fracture. As compared with the empirical derivation in equation (2.2), the solutions for critical time t c (or equivalently strain e c ) obtained from the energy balance approach involve more complex nonlinear time-convolution functions P ± . Nevertheless, they again reveal the dependence of t c (or e c ) on three physical dimensionless quantities. This further confirms the important role of these dimensionless parameters in controlling the wrinkle-induced delamination process of the viscoelastic foulant layer. Note also that in the elastic limit E 0 = E ∞ , thus β = 0, the critical strain for delamination becomes the same as presented in the previous work that focused on elastic foulant where the numerical scaling pre-factor a is used in the scaling law (see equation (30) in electronic supplementary appendix 2, figure S5 in electronic supplementary material, appendix 3 and [8]).
The wrinkling strain e w is assumed to be small, therefore the energy and relaxation in the foulant layer prior to wrinkling can be neglected in the above analysis. The effect of e w , however, is taken as a shifting parameter to the solution obtained above as suggested in Pocivavsek et al. [8]: e c À e w . Shown in figure 3b are the critical strains for delamination obtained from solving equation (2.4) when the three physical dimensionless parameters are varied. While a decrease in E ∞ /E 0 or in _ et R leads to an increase in the critical strain, a decrease in e de reduces the critical strain. However, as the figure shows, the effects are highly nonlinear in the intermediate regime between the two elastic limits. In addition, smaller E ∞ /E 0 values correspond to larger differences between the solutions using total work W and the stored energy U in the energy balance approach. When there is no intrinsic material relaxation, (E ∞ /E 0 = 1), the solutions from both approaches coincide with the prior linearly elastic case [8]. Furthermore, in the two elastic limits _ et R ) 1 and _ et R ( 1, the two limiting elastic solutions e c À e w ¼ e de and e c À e w ¼ ðE 0 =E 1 Þe de are obtained. Analysing the behaviour of the analytical functions P ± (z) in equation (2.4) confirms the same limiting behaviour (see electronic supplementary material, appendix 2). It is important to note that e c À e w does not simply increase linearly with the dimensionless parameter e de . Instead, we observe a right shift (in the arrow direction shown in figure 3b) for the transition from the instantaneous response to the other regimes which can be attributed to the nonlinear effect of τ de /τ R .

Finite-element analysis
To further study the roles of the dimensionless control parameters and their influence on e c , as well as verify the trends observed from the analytical method, a detailed parametric study is performed using finite-element method (FEM) with cohesive zone model (CZM) implemented in Abaqus (Dassault Systèmes, MA) [47]. In this approach, traction separation laws are prescribed between the foulant layer and the film interface to study the delamination process [8,25,[29][30][31][32][33][34][35][36][37]. Two input parameters, the cohesive strength σ c and the fracture energy G, together with a damage law are necessary to describe the computational cohesive laws. In this study, bilinear softening laws are employed to improve the numerical implementation in the previous study which employed linear elastic brittle laws [8] (see electronic supplementary material, appendix 3). The contributions of the two CZM parameters σ c and G to the fracture process are also conveniently studied with this implementation. Solving equilibrium equations set by balancing the total energy, including contributions from both the foulant layer and the cohesive interface, allows the determination of e c . In Abaqus, this solution process can be based on an implicit or explicit solver [47]. The implicit solver offers the advantage of solving this quasi-static problem without introducing additional dynamic energy. However, the presence and coupling of surface instability, material softening, contact conditions and interfacial delamination require tuning of various solver parameters and the introductions of certain artificial energies, such as damping, to resolve convergence issues of this iterative solution scheme. Therefore, in this study, we employ Abaqus dynamic explicit solver in order to conduct a parametric study with minimum adjustment of solver parameters over a large space of material properties and varying loading rates. Furthermore, in the CZM method, the process zone length along the interface L pz ¼ EG=s 2 c , which is the length over which CZM elements enter the degradation part of the traction separation law, affects the delamination mechanism in the problem. It has been shown in the literature for several classical interfacial geometries that L pz has to be smaller than a characteristic length of the system for CZM to produce the same solution as the energy-based approach, such as L pz < h for double cantilever beam, where h is the beam thickness [32,33], or L pz < L 2 /(2h) for edge delamination, where L, h are the length and thickness of the layer [29]. The role of L pz was not considered in prior topography-driven delamination work [8]. Thus, here we conducted a detailed sensitivity study for the effects of CZM parameters on e c in the instantaneous and long-term limits for the viscoelastic foulant layer. At these limits, the viscoelastic foulant layer can be treated as an elastic foulant layer with modulus E p = E 0 = E ∞ and the scaling law becomes: e c ¼ al 2 G=B p , where a is a constant pre-factor [8]. Figure 4 plots the normalized ratio ðe c À e w Þ=ðl 2 G=B p Þ, which provides the parameter a, as a function of L pz for different foulant layer's modulus E p . For each foulant layer's modulus, simulations with different sets of CZM parameters σ c and G were performed.
As shown in figure 4, the presence of a process zone in CZM might influence the transition of different delamination mechanisms. We can compare the predictions from our analytical methods, which use only fracture toughness, with the FEM simulations only in the regime where the FEM data collapse to a flat line corresponding to a constant parameter for a. In this regime, fracture is dominated by energy, and the cohesive strength has negligible effect. Outside this regime, the strength might play a significant role and hence the scaling law e c ¼ al 2 G=B p cannot be used to interpret the FEM data. As the scope of this paper is on the wrinkle-induced delamination mechanism of a viscoleastic foulant layer, focus on determination of the transition between these delamination mechanisms in FEM simulations will be presented in a separate publication. A summary is presented in electronic supplementary material, appendix 3 to emphasize that the CZM analysis is conducted with careful attention to important numerical aspects including mesh royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220598 refinement, process zone length, strength and energy dominated regimes [25,29,[31][32][33][34][35][36][37]. Our simulation results for the delamination of a viscoelastic foulant layer from a wrinkled surface in the energy regime are plotted in figure 5a. When compared with figure 3b, the FE data confirm a similar dependence of the critical strain on the three dimensionless control parameters. A consistent right-shift of the transition from the instantaneous response to the other regimes as e de increases is also observed.

General solution
Our analytical method and simulations indicate that the delamination of a viscoelastic foulant layer from a wrinkled surface shows a complex dependence on three control parameters. We, therefore, investigate the design parameter space by deriving a general fit to collapse the simulation data for a wide range of these parameters. Equation ( which predicts the collapse of all the data into a region defined by the bounded functions 0 < P ± < 1. With the use of this general solution, all data for the analytical method can be collapsed to two general curves corresponding to either the use of W or U, as shown in figure 5b. FEM data are more noisy due to various numerical factors in the simulations such as the determination of the onset of delamination; nevertheless, they also nicely collapse into the small area bounded by the two analytical solutions. The FEM data also show that the amount of energy available to drive fracture in the simulations is bounded by the two limiting cases considered in the analytical approaches. The upper-bound and lowerbound curves are physically intuitive as the stored energy provides the least amount of possible energy while the total work provides the largest possible energy available to drive fracture. As discussed above, we conducted an FEM analysis using an explicit solver. We controlled the dynamic effect to be small to represent a quasi-static condition; however, certain amount of dynamic energy due to the solution process and dynamic propagation of the crack tip may still be present in these highly nonlinear simulations, contributing to a difference in the amount of available energy in the simulations as compared with either limits in the analytical approach. Though neither of them exactly overlaps the numerical data from analytical methods, they provide good insight into the controlling factors in this system, and how to use them to control the delamination process. The general collapse obtained in figure  5b further confirms a reasonable agreement between the theoretical and numerical predictions in how the wrinkle-induced delamination mechanism of the viscoelastic foulant layer is controlled by different physical parameters related to the viscoelastic properties and loading process. Equation (2.6) shows that the three physical parameters can further be combined into two dimensionless parameters, _ et R =e de and β, significantly reducing the design parameter space needed for optimizing this complex system.

Discussion
Topography-driven surface renewal is a powerful new mechanism to drive interfacial fracture over large surfaces at risk of fouling. The initial theory was limited to an elastic response regime for the fouling layer [8]. Nevertheless, the mechanism shows promising applications in medical device design, specifically anti-thrombotic vascular grafts [9,10]. The application of actuation into medical devices requires precise knowledge of target strains, which will inform the choice of graft materials, sources of actuation loading, and range of biofoulants that the mechanism will effectively remove from the surface. In this paper, we substantially enrich the existing theory on actuating topography to account for both dynamics in the surface (strain rate _ e) and the viscoelastic nature of the fouling layer. The results of this paper are directly applicable to the design of actuating biomedical devices targeting viscoelastic foulants such as thrombus, bacteria and biofilms [6,7,9,10].
In summary, our FE results show consistent agreements with analytical predictions of the nonlinear interplay between the three control parameters, which suggest several conditions to promote delamination in the viscoelastic regime. With the same rate of loading and adhesion energy, the critical strain for delamination onset decreases with increasing relaxation time of the foulant layer τ R . In other words, a foulant layer that relaxes slowly is easier to de-adhere. The limit of very fast relaxation, in which the thin foulant layer approaches the fluid-like limit, needs further investigation, as FE simulations in this regime indicate that another mechanism involving the critical strength, rather than only fracture toughness, may play a role here. With the same adhesion energy and same relaxation time, increasing loading rate _ e reduces the critical strain for the wrinkle-induced royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220598 delamination onset. Furthermore, with the same relaxation time and rate of loading, decreasing adhesion energy also decreases the critical strain for the wrinkle-induced delamination onset. In addition, comparisons are also made for the energy balance approach using the stored energy and the total energy, and a general collapse of simulations data to a bounded regime between these analytical predictions helps to reduce the design parameter space. Thus, this work provides a first attempt towards understanding the delamination of realistic viscoelastic biofoulants. Ultimately, topographic surface renewal is based on an energy release mechanism whereby external loading and deformation of an adherent foulant leads to accumulated energy in the foulant, which beyond a critical value is able to drive interfacial crack propagation. In the purely elastic case, as shown in the key result of our prior work (see equation (2) in Pocivavsek et al. [8]), the critical strain density needed to drive surface renewal is given by e c =l 2 1=' 2 ec , where e c is the critical nominal compressive strain in the bilayer substrate driving surface wrinkling and foulant deformation and ℓ ec = (B 0 /G) 1/2 is the elasto-capillary length scale. Because energy dissipation only occurs with interfacial fracture, the elastic case is independent of loading history; as such it is independent of actuation in the dynamic sense.
Unlike in the elastic case, the presence of viscoelasticity in the foulant introduces a second intrinsic mode of dissipation in addition to interfacial fracture. Furthermore, the strain state of a viscoelastic foulant layer is a superposition of states because of the intrinsic material softening triggered by loading. The coupling of the nonlinear, wrinkle-induced strain field with the nonlinear dependence of the strain state on loading leads to a highly nonlinear equation to determine the critical delamination strain for a viscoelastic foulant layer. We show that this critical strain is controlled by three physical quantities: the magnitude of foulant layer relaxation E ∞ /E 0 , the rate of material relaxation relative to loading rate _ et R and the critical strain to delaminate the foulant layer in the instantaneous elastic limit e de l 2 G=B 0 . A fourth dimensionless parameter is the strain to initiate wrinkling e w , which can be tuned by device construction and made negligibly small in the case of vascular graft design [3,4,9,10]. In the first part of this paper, we show that viscoelasticity, characterized as the decay of foulant stiffness from E 0 to E ∞ over a characteristic time scale τ R , breaks the critical strain for topographic de-adhesion into two limiting regimes: e de e c À e w ðE 0 =E 1 Þe de . The smallest strain is set by the elastic limit [8]. However, the critical strain increases proportionally to the degree of degradation in stiffness to an upper bound set by E 0 /E ∞ . This inequality can be directly obtained from our general solution by solving for e c in equation (2.6) or equivalently by using figure 5b to determine the value of where f ∈ [0, 1] is the range of the vertical axis in figure 5b. By moving f, the critical strain moves between the short-and long-time limits which are connected by a highly nonlinear transition region. Equation (3.1), and equivalently figure 5b, provide a unified approach to optimize strategies for actuated topographic surfaces designed to remove viscoelastic foulants. The horizontal axis in figure 5b is controlled by where τ is the applied actuation time. The first set of parentheses contains parameters under direct control of the device designer, and the second set of parentheses incorporates intrinsic properties of the foulant.
An effective anti-fouling design must enforce the system to work under short time conditions to lower the critical strain to the value e de , which is achieved when f = 0. Physically, this is equivalent to a state where intrinsic viscoelastic dissipation has not had time to take effect and all accumulated strain energy is available to drive interfacial fracture. Following the results in figure 5b and equation (2.6), a sufficient condition to be in this regime is to satisfy _ et R =e de * 10. Thus, a system designed to be strained to a value e along a time scale τ must satisfy the condition e=t * 10e de =t R . Equivalently, the system must be designed to have a strain rate larger than e de =t R . In terms of the elasto-capillary length scale whose value depends on the specific response of the foulant, we obtain the following  Figure 5. Finite-element results show how the three controlling dimensionless parameters affect the critical strain (a) and general fits to collapse numerical and analytical data (b). Four marker types (downward pointing triangle, square, diamond and circle) correspond to four loading rates _ e=_ e el ¼ 0:75, 1, 1:5, 2, respectively, used in FE simulations, where _ e el is the loading rate adopted from the FE simulations in the previous study for elastic foulant layer [8] (see electronic supplementary material, appendix 3). Open and closed markers are for E ∞ /E 0 = 0.5 and E ∞ /E 0 = 0.25, respectively. Increasing levels of grey colour of the markers for the numerical data correspond to increasing values of e de ; however, the colour code and line styles in (b) for the analytical results are consistent with figure 3b.
In general, if one moves f up from 0, then the critical strain given by equation (3.1) increases. This is important from a design perspective as it provides a design criterion to improve the de-adhesion capability of the system if the optimal surface renewal condition cannot be achieved. For instance, if the system is forced to operate at _ et R =e de * 5 then the largest value of f becomes 0.3 and the critical strain needed will be given by equation (3.1).
In topography-driven surface renewal for viscoelastic foulants, the fouling layer accumulates strain density (e=l 2 ) over an applied actuation time (τ). We will take the example of vascular grafts as an illustration, although similar arguments can be constructed for any specific application. The designer has full control over all three parameters independently: e is the nominal strain in the graft substrate and will be set by E s and the load available in the system (for example, pulse pressure in vascular grafts [9,10]), surface wavelength λ is tunable via graft bilayer construction [3,4] and lastly τ is either set by the system (heart rate in the case of vascular grafts [9,10]) or, if externally driven, by an actuating power source (such as in soft robotic on-demand fouling-release urinary catheters [6,7]). These three parameters must combine such that the time-integrated accumulated strain density is within the limits set by purely intrinsic properties of the foulant layer: G, E 0 , E ∞ , h and τ R .
Actuated topographic vascular grafts (topografts) in contact with whole blood were shown to have optimal surface renewal when actuated at 1 Hz with a wavelength of 80 μm [9,10]. At a wavelength of 250 μm, the surface renewal capability was shown to be reduced by 50%. When the wavelength was increased to 1000 μm, the topography had no added effect on surface renewal (see fig. 5 in [9]). Surface renewal using just the elastic limit only predicts whether topographic delamination occurs or not, based on if the applied strain is larger than the critical threshold e de . Thus, it can be used to interpret the observations for the two wavelengths 80 and 1000 μm. However, it is not sufficient for explaining the cross-over delamination capability observed at the intermediate wavelength 250 μm. Figure 5b and equation (3.2) can be used to correlate all these observations. The lefthand side of equation (3.2) contains factors under control of the designer. In the case of topografts: e 0:15, λ = 80, 250, 1000 μm, τ = 1 s, so that the surface wavelength was studied as a design control parameter [9,10]. The right-hand side of equation (3.2) contains parameters intrinsic to the foulant. The whole blood thrombus is assumed to be incompressible with E 0 = 10 kPa, ν = 0.5, h = 100 μm [9,48]. The fracture energy G can be calculated from published cohesive strength of whole blood σ c ∼ 2 kPa [49]. Using the analysis outlined in §2.3 with L pz ≈ h, G ¼ hs 2 c =E 0 ¼ 0:04 N m À1 or 40 mJ m −2 , equivalent to a rubber/water interface [50]. Whole blood thrombus viscoelasticity has been studied and noted to fit well to a second-order Prony expansion with t 1 R 1 s and t 2 R 30 s with weak dependence on strain rates [48]. The e de ¼ al 2 G=B 0 corresponding to the three wavelengths λ = 80, 250, 1000 μm are 0.0058, 0.056 and 0.9, where B 0 = E 0 h 3 /12(1 − ν 2 ) is the bending stiffness and the numerical pre-factor a is determined in figure 4 and electronic supplementary material, appendix 3. Combining these terms provide the values for the left-hand side ðe=l 2 Þ=t of equation (3.2) to be in the range of 2.3 × 10 7 (m 2 s) −1 for λ = 80 μm, 2.4 × 10 6 (m 2 s) −1 for λ = 250 μm and 1.5 × 10 5 (m 2 s) −1 for λ = 1000 μm. The ability of the actuating topography to function as an optimal surface release mechanism according to equation (3.2) can only occur if ðe=l 2 Þ=t is larger than 10aG/τ R B 0 . A range of the value for this right-hand side for the whole blood thrombus is from 9 × 10 6 (m 2 s) −1 for t 1 R 1 s to 3 × 10 5 (m 2 s) −1 for t 2 R 30 s. In addition, with e w ¼ 0 and E ∞ = 0 kPa, the value in the vertical axis of figure 5b will be 0.96, 0.62 and −5 correspondingly for the three wavelengths λ = 80, 250, 1000 μm. With these computed properties, we can interpret the trend of the above experimental data for thrombus delamination.
For the wavelength λ = 1000 μm, the applied strain e ¼ 0:15 , e de ¼ 0:9 is not enough to induce delamination. This is also consistent with the negative value −5 of the ratio on the vertical axis of the master curve in figure 5b. The applied strain e ¼ 0:15 is big enough to induce delamination for the wavelength λ = 80 μm. Furthermore, at this smallest wavelength, the value of ðe=l 2 Þ=t ¼ 2:3 Â 10 7 ðm 2 sÞ À1 is larger than the upper bound 9 × 10 6 (m 2 s) −1 of the righthand side for the two relaxation times, suggesting that optimal condition for delamination (equation 3.2) is always satisfied. This confirms the capability of these foulants to fail in the elastic limit when the smallest wavelength is used. The high value 0.96 of the ratio plotted in figure 5b further shows that the critical strain is far above the threshold to induce delamination. The interesting phenomenon at the intermediate wavelength λ = 250 μm is also explained with the incorporation of viscoelastic foulant behaviour. The applied strain e ¼ 0:15 is still larger than the elastic delamination threshold e de ¼ 0:056, corresponding to the value 0.62 on the vertical axis of figure 5b. Thus, surface renewal using just the elastic limit predicts that the delamination occurs as in the optimal wavelength λ = 80 μm. This is not consistent with the experimental observation where a significant reduction in delamination capability is observed. However, by considering the inequality (3.2), at this wavelength λ = 250 μm, the left-hand side 2.4 × 10 6 (m 2 s) −1 lies in the middle of the right-hand side range 9 × 10 6 (m 2 s) −1 for t 1 R 1 s to 3 × 10 5 (m 2 s) −1 for t 2 R 30 s. This suggests that whether optimal delamination is satisfied or not depends on the relaxation properties of the foulants. In addition, the intermediate value 0.62 of the ratio plotted on the vertical axis of figure 5b suggests that delamination may or may not occur depending on how the dissipation energy plays a role in the delamination process. Specifically, at this wavelength, the value _ et R =e de on the horizontal axis of figure 5b is 2.67 and 80 for the two relaxation times t 1 R and t 2 R , respectively. Thus, the corresponding point on the master curve for t 1 R lies above the critical point predicted using the total energy approach, but is below the critical point predicted using the stored energy approach. This indicates that viscous dissipation is likely to play a role in dampening the effect of topographic surface renewal at this wavelength, which can be why in practice the capability for delamination to happen is reduced at λ = 250 μm when compared with at the optimal wavelength λ = 80 μm. This analysis shows that going to even smaller wavelengths would most improve surface renewal given the nonlinear dependence on λ. Thus, our royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20220598 analysis allows a concrete, quickly applicable methodology to pick design parameters for given anti-fouling applications of topography-driven surface renewal. It shows that to optimize topographic de-adhesion for viscoelastic foulants the condition set forth in equation (3.2) should be satisfied. This helps guide design of anti-fouling surfaces targeted for biomedical applications such as anti-thrombotic grafts.
The energy balance and numerical approaches in this work and the previous publication in the elastic limit [8] provide a general framework to study the interaction between the external loading, the foulant deformation and the surface energy. Thus, it will be the basis for future studies that focus on incorporating more complexity in loading conditions, material behaviours and interfacial properties. Specifically, the relaxation behaviour of the foulant will depend on the frequency of the external loading cycle, such as the constant loading and unloading (equivalently contraction and expansion) of the arterial wall [17][18][19]. Hence, the energy balance must consider how the energy in the foulant is accumulated and relaxed during multiple excitations to determine the critical number of cycles to induce foulant delamination. In our model, we assume that when delamination occurs, it is an irreversible process and the detached foulant does not reattach to the arterial wall when the wall relaxes. As such, the new configuration with this new crack length has to be used in the next contracting and expanding cycle of the arterial wall. This assumption should be experimentally tested and our model can be further enriched if partial reattachment phenomena occurs. Another interesting phenomenon that can be included to extend our studies is the relaxation behaviour in the arterial wall during a cycle. This will lead to changes in the wrinkle topography where the wrinkle wavelength and amplitude are influenced by the viscoelastic properties of the wall [51,52], consequently affecting the scaling for the critical strain. Furthermore, adhesion in biological materials can be a rate-dependent process where the fracture toughness G depends on the loading rates [53]. This complexity, if being added to our energy balance model, will lead to the dependence of e de on the loading rate even in the elastic foulant limit. When the foulant is viscoelastic, we anticipate another dimensionless parameter relating how fast the foulant relaxes to how fast G changes with loading will also become a control parameter. For the numerical approach, a rate-dependent cohesive zone model with mixed-mode failure conditions can be developed and incorporated into our model to tackle more complex interfacial delamination processes. Additionally, the flat configuration in this work is applicable to study the limit when the length of the foulant, which covers multiple wrinkle wavelengths, is much smaller than the radius of the arterial wall. The curvature of the artery can be conveniently introduced into our model as a next step to study its effect on the delamination strain when this limit does not apply.
Data accessibility. Data for this study can be found in the main manuscript and the electronic supplementary material files and video [54].