In the present study we propose two phenomenological CZM formulations to represent fibrillation during dissection of soft tissue. We firstly consider elastic fibrillati on. In this formulation fibrillation initiates during an initial phase of interface damage/softening. Fibrils then deform elastically until an ultimate strength is reached, followed by fibril damage/ rupture. Simulations reveal that such elastic fibrillatio n does not accurately predict the highly non -uniform crack propagation rates measured experimentally for arterial tissue. We then propose a phenomenological visco -plastic cohesive zone model (VP- CZM) for fibrillation. This approach is motivated by the obse rvation that fibrils undergo partial pull -out (in addition to elastic stretching) during dissection of arteries. Simulations reveal that the VP -CZM provides a reasonable prediction of non -uniform crack propagation rates in arteries. Significant plasticity (representing fibril pull -out) during fracture initiation results in an initial slow phase of crack growth. The transition to fast crack growth is facilitated by the prediction of fibril rupture pr ior to extensive plasticity. The current implementation of fibril visco -plastic behaviour in a cohesive zone framework is limited to a phenomenological representation of axial fibril deformation and axial fibril pull -out. However, we provide extensive parametric exploration of the model behaviour under a range of mixed -mode loading paths, representing fibril reorientation in addition to fibril extension. We demonstrate that the formulation produces consistent behaviour for all loading paths and tha t positive incremental instantaneous dissipation is computed in all cases. The VP -CZM represents an advance on standard elastic -damage cohesive zone formulations. While the VP -CZM effectively produces such standard CZM behaviour at high interface separation rates, it predicts plastic deformation and energy dissipation (re presenting fibril pull -out) at low interface separation rates. The authors are not aware of a previous implementation of fibril visco -plastic behaviour in a CZM framework. Keywords : Cohesive zone, fibrillation, collagen, soft tissue, viscoplasticity 1 Introd uction 1.1 General Background Cohesive zone models are commonly used to model the delamination of an interface [1–5]. CZMs can be formulated as continuous functions [3,6,7] , or they may be expressed as piecewise functions [8,9] . Interfacial fibrillation is a mechanism whereby delamination takes place by initiation, growth, extension, and ultimately fracture of fibrils between fracture surfaces [10]. Previously, standard cohesive zone mode ls (CZMs) have been applied to model interfacial fibrillation in a variety of engineered materials including polymers, elastomers, paper, and plastics [11–13]. Van den Bosch et al., adapting earlier work, proposed the use of a single relation between the opening displacement and the traction to account for large displ acement problems, such as fibrillation [10,14] . This model was implemented by Vossen et al. a s part of a multiscale modelling approach to describe the fibrillation micromechanics in a delaminating copper -rubber interface [11]. In this approach idealised hyperelastic fibrils are explicitly modelled as part of the multiscale modelling approach. The intrinsic adhesion between the fibrils and the copper is described by the previously developed standard cohesive zone model [10]. Another study examining fibrillar metal -elastomer interfaces explored the role of discrete fibrils within the fracture process zone [13]. Several discrete hyperelastic elastomer fibrils were explicitly modelled without the use of a cohesive zone model. An anisotropic CZM for modelling fibrillar and crazing interfaces in cellulose nanopaper was proposed by Paggi & Reinoso [12]. This model considers a statistical distribution of fibril o rientations, with debonding of a fibril based on the mechanics of peeling of an elastic film from a substrate [15]. Viscoplastic CZM models of c razing in polymers are developed in the studies of Tijssens et al. [16] and Saad -Gouider et al. [17]. In addition to crazing and fibrillation in polymers , interfacial fibrillation has also been experimentally observed during fracture of biological materials such as hair [18], skin [19] and other collagenous tissues such as arterial tissue [20,21] . Standard c ohesive zone formulations have previously been used to model fracture in biological tissues such as arteries [22–24], and cells [25], but viscoplastic toughening mechanism s have not been previously considered . 2 1.2 Examination of fibrillation during mode II delamination of arteries uncovered in recent experimental s tudy (McGarry and FitzGibbon (2020 ) The series of mode II fracture tests on arterial tissue presented in FitzGibbon and McGarry [21] demonstrates the rate of crack propagation is non unifo rm and that significant fibrillation occurs, both during initiation and propagation (see Figure 1). While a detailed description of the novel mode II soft tissue fractu re methodology is beyond the scope of the current paper, the key evidence of fibrillation and nonuniform crack propagation rate is outlined in Figure 1. A schematic of the experiment is shown in Figure 1(a). Briefly, an excised ring of aorta is placed on two circular loading bars. A radial notch is inserted to the side of the loading bars where high shear stress occurs. Displacement of loading bars (characterised by an effective circumferential strain Δ𝐶/𝐶0) to a critical level results in kinking of the crack with mode II initiation and propagation in the circumferential direction al ong a circumferential -axial plane. A J -integral analysis demonstrates that crack -kinking and mode II fracture occurs because collagen fibre alignment in front of the crack tip provides a significant toughening mechanism against mode I crack growth in the r adial direction , through the thickness of the sample. Measured mode II crack propagation rates are shown as a function of the applied circumferential strain in Figure 1(b). Three characteristic regimes are observed; (i) initial slow crack growth with fibrillation (blue); (ii) fast crack growth over a long distance (red) without significant fibrillation ; (iii) slow crack growth with significant fibrillation (green). Significant fibre formation is observed during the test (Figure 1(c)) and observed in the region of the crack tip following specimen removal from the test rig ( Figure 1(c)). Fibrils are not simply elastically stretched between opposing crack surfaces; rather, they appeared to ha ve undergone partial pull -out, such that removal of elastic strain energy does not result in significant fibril shortening and closure of crack surfaces. Fi nally, a standard cohesive zone model (CZM) that implements critical traction damage initiation, and exponential softening/damage evolution was used to simulate these mode II fracture tests (Figure 1(d-e)). The mode II fracture strength (𝜏𝑚𝑎𝑥) primarily governs crack initiation, whereas the mode II fracture energy (G II) primarily governs the crack propagation rate. However, the use of a standard CZM is shown to compute an approximately uniform crack growth rate , providing only approximate predictions of experimental measurements without capturing the three regimes of crack growth rate described above ( GII=0.05 N/mm captures only the slow growth rate at the start of the experiment, whereas GII=0.0005 N/mm captures the fast crack growth rate in r egime (ii), but not the slow crack growth in regimes (i) and (iii) (Figure 1(f)). 1.3 Structure and aims of the study In the current study we aim to develop new CZM approa ches to explore the relationship between fibrillation and the experimentally observed non -uniform crack propagation rates in arteries . In Section 2 we propose a new elastic fibrillation cohesive zone model ( EF-CZM). F ollowing initial damage evolution, this formulation incorporates elastic stretching and fracture of fibrils. This approach results in improved predictions compared to a standard exponential damage CZM, but it does not accurately capture experimentally measured non -uniform propagation rates in a rteries. Furthermore, analysis of our experimental fracture surfaces suggests that fibril s undergo pull -out, in addition to elastic stretching, during propagation. Based on this observation we propose a novel visco -plastic CZM (VP-CZM) formulation to represent rate dependent fibril pull -out and associated dissipation during crack propagation. We demonstrate that the VP -CZM provides an improved prediction of non-uniform crack propagation during mode II dissection of arterial tissue. 3 Figure 1. (a) J -Integral analysis of the novel mode II fracture experiment illustrating significant mode I toughening and fibre alignment behind the crack tip. (b) Measured crack length is shown as a function of the applied circumf erential strain. Three characteristic regimes observed and are approximately indicated on the curves; (blue) initial crack growth is slow due to fib rillation; (red) crack growth velocity increases due to the rupture of fibres concentrating stresses at the crack tip; (green) crack growth slows due to fibrillation. (c) Experimental image showing fibrillation occurring, a large fibril is highlighted (red box). Images s how detailed views of interfacial fibrillation. Experimental image of the ruptured aorta spec imen showing significant permanent pull -out of fibrils; (d) Traction -separation response of the exponential CZM used previously; (e) image of the model in its undeformed configuration and two example contour plots showing crack propagation; (f) Crack growt h predictions of the exponential CZM for two different mode II fracture energy (G II) 1.4 Analysis procedures and finite element implementation In this study we develop two new CZM formulations, to analyse these CZMs we use two assessment techniques : i. Analysis of the model in response to mixed mode proportional and non -proportional specified loading paths in which the interface separation history is fully prescribed 4 ii. Analysis of the f inite element implementation of the model to simulate previously reported mode II fracture experiments [21]. Mixed mode s eparation -specified loading paths : Proportional and non -proportional loading paths are considered in the pre sent study. We describe a proportional loading path as a separation in which the mode angle (𝜑) remains constant with increasing magnitude of interface separation (Δ𝑚), this is illustrated in Figure 2(a). We describe a non -proportional loading path as a separation in which 𝜑 may change throughout the separation history, in the present study we consider one form of non- proportional loading; an initial tangential separation (Δ𝑡=Δ𝑡𝑚𝑎𝑥) followed by a subsequent normal separation to failure (Δ𝑛= ∞), as shown in the illustration in Figure 2(b). For the re mainder of the present study this path is referred to as non -proportional loading. M odel performance is assessed in terms of tangential and normal coupling, mixed mode fracture energy , and instantaneous incremental energy dissipation [26]. Finite element implementation: A mesh of ~2350 generalised plane strain elements is used to simulate the mode II fracture experiment (as shown in Figure 2(d)). Convergence of the predicted crack length was achieved for the mesh used. A cohesive zone model is implemented at the medial interf ace using u ser-defined interfacial constitutive behavio ur subroutine (UINTER) . The use of the recently developed constitutive model from Fereidoonnezhad et al. to model the aorta allows for the control of unphysical auxetic behaviour observed in exponentia l constitutive laws typically used for arterial tissue [27–29].. The isochoric strain energy density function is given as : 𝜓𝑖𝑠𝑜(𝜆̅1,𝜆̅2,𝜆̅3)=∑𝜓̅(𝜆̅𝑖)3 𝑖=1, 𝜓̅(𝜆̅𝑖)= { 𝐸1(𝜆̅𝑖−ln𝜆̅𝑖−1) |𝜆̅𝑖−1|≤𝐷1 𝑝(𝜆̅𝑖2 2−2𝜆̅𝑖+ln𝜆̅𝑖)+𝑞(𝜆̅𝑖−ln𝜆̅𝑖)+𝑟ln𝜆̅𝑖+𝜓01𝐷1<|𝜆̅𝑖−1|<𝐷2 𝐸2(𝜆̅𝑖−(1+𝐷2)ln𝜆̅𝑖)+(𝑝𝐷22+𝑞𝐷2+𝑟)ln𝜆̅𝑖+𝜓02|𝜆̅𝑖−1|≥𝐷2 (1) where 𝐷1,𝐷2,𝐸1,and 𝐸2 are material parameters, 𝜆̅𝑖 (𝑖=1,2,3) are the isochoric principal stretches, 𝐽 =𝜆1𝜆2𝜆3 is the jacobian , and 𝜓01 and 𝜓02 are two constants which ensure the continuity of strain energy. Moreover 𝑝,𝑞, and 𝑟 are not independent parameters; in order to maintain C0 and C1 continuity the f ollowing relations must be enforced: 𝑝=𝐸1−𝐸2 2(𝐷1−𝐷2), 𝑞=𝐸1−2𝐷1𝑝, 𝑟=(𝐸1−𝑞)𝐷1−𝑝𝐷12 (2) For further details of the constitutive law used see Fereidoonnezhad et al. [28]. We assume a slight material compressibility (𝜈𝑒𝑓𝑓=0.445) based on the compressibility experiments of Nolan et al. [30]. The ability of the calibrated constitutive law to predict the non -linear anisotropic behaviour of aortic tissue has been previously demonstrated by FitzGibbon and McGarry [21] and is shown here in Figure 2(c). 5 Figure 2. (a) illust ration of proportional loading paths examined in the present study; (b) illustration of non -proportional loading paths examined in the present study; (c) best fit of the material constitutive behaviour to arterial tissue; (d) illustration of the mesh used with the radial notch (highlighted red), CZM (highlighted blue) and loading bars shown. The undeformed and deformed configureation s are shown. 𝚫𝑪/𝑪𝟎=𝟎.𝟔𝟖 is the point of crack initiation and 𝚫𝑪/𝑪𝟎=𝟎.𝟖𝟖 shows the progression of the mode II crack. 2 Elastic fibrillation cohesive zone model (EF -CZM) In an effort to phenomenologically capture the fibril stretching between opposing crack flanks, we propose a new elastic fibrillation CZM (EF -CZM). As illustrated in Figure 3, this framework entails an initial intrinsic elastic region (shown in blue) , followed by a standard exponential softening/damage region (yellow) . Fibril formation is assumed to occur as part of the damage process, with fibrillation initiating when traction softens to a critical value ( 𝑇𝑓𝑖𝑛𝑖𝑡). Fibrils are assumed t o extend elastically (green) until they reach an ultimate fibril stress ( 𝑇𝑓𝑢𝑙𝑡), beyond which, fibril damage and softening initiates, ultimately leading to complete fibril rupture (red). Critical tractions associated with fibrillation initiation and u ltimate rupture can be specified as a function of the mode mixity. As an example, Figure 3(b) illustrates a case of mode dependent fibrillation whereby fibril s extend to a higher traction under mode II loading than under Mode I loading. A full calibration of the model requires the specification of the following parameters for both mode I and mode II: interface strength (𝑇𝑚 𝑚𝑎𝑥), interface stiffness ( 𝐾𝑚), fibril sti ffness (𝐾𝑓), fibril strength (T𝑚𝑓𝑢𝑙𝑡), interface fracture energy prior to fibrillation (𝐺̂𝑚), and fibrillation fracture energy (𝐺𝑚𝑑). The model also requires the specification of the mode dependence of the aforementioned parameters. In Section 2.1 we present full details of the EF-CZM formulation. In Section 2.2 we present rigorous parametric investigation of the model behaviour under mixed mode conditions. In Section 2.3 we explore the capability of the EF -CZM to predict the non - uniform crack propagation rates observed experimentally in our artery fracture experiments. 6 Figure 3. (a) Schematic of the EF -CZM undergoing a separation. 𝑲𝒇 is the fibrillation stiffness, 𝑻𝒎𝒇 is the fi brillation strength, and 𝜹𝒇𝒆𝒍 is the elastic limit of separation of the fibrillation regime. (b) 3D traction -separation schematic of the EF -CZM showing a mode II response (red), mixed mode response (black), and a mode I response (blue). (c) CZM paramet ers are shown as a function of mode angle. Interface strength, elastic limit, separation at fibrillation initiation, maximum fibril extension, and fracture energ y are all shown. 2.1 EF-CZM formulation For a given interface displacement vector 𝚫 with normal a nd tangential (shear) components Δ𝑛 and Δ𝑡, respectively , the corresponding displacement magnitude is given as Δ𝑚=(Δ𝑛2+Δ𝑡2)1/2 and the mode angle is given as 𝜑=tan−1(Δ𝑡Δ𝑛⁄). The magnitude of the interface traction is expressed as a function of Δ𝑚 and 𝜑 by the following formulation : 7 T𝑚(Δ𝑚,𝜑)= { 𝐾𝑚Δ𝑚 ,Δ𝑚<T𝑚𝑚𝑎𝑥(𝜑)/𝐾𝑚 T𝑚𝑚𝑎𝑥(𝜑)(exp(−𝛥𝑚𝑚𝑎𝑥−T𝑚𝑚𝑎𝑥(𝜑)/𝐾𝑚 𝛿𝑚𝑑∗(𝜑))Δ𝑚 Δ𝑚𝑚𝑎𝑥) ,T𝑚𝑓𝑖𝑛𝑖𝑡(𝜑)/𝐾𝑓>Δ𝑚≥T𝑚𝑚𝑎𝑥(𝜑)/𝐾𝑚 𝐾𝑓Δ𝑚 ,T𝑚𝑓𝑖𝑛𝑖𝑡(𝜑)/𝐾𝑓<Δ𝑚≤T𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓𝛿𝑓𝑒𝑙 T𝑚𝑓𝑢𝑙𝑡(𝜑)(exp(−𝛥𝑚𝑚𝑎𝑥−T𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑚 𝛿𝑚𝑓∗(𝜑))Δ𝑚 Δ𝑚𝑚𝑎𝑥) ,T𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓≤Δ𝑚 (3) 𝐾𝑚 is the intrinsic stiffness of an intact interface. In the current study we assume that 𝐾𝑚 is mode -independent . T𝑚𝑚𝑎𝑥(𝜑) is the specified mode -dependent interface strength; the corresponding displacement (𝛿𝑚𝑒𝑙(𝜑)=T𝑚𝑚𝑎𝑥(𝜑)/𝐾𝑚) represents the elastic limit of interface se paration at the point of damage initiation . The mode -dependent parameter 𝛿𝑚𝑑∗(𝜑) governs the rate of softening in the initial 𝐺𝑚𝑏(𝜑) damage region ( yellow ). Note that damage depends on the maximum separation throughout the displacement history, 𝛥𝑚𝑚𝑎𝑥, such that a reduction in interface separation during the damage region will directly result in elastic unloading. T𝑚𝑓𝑖𝑛𝑖𝑡(𝜑) is the specified mode -dependent traction at which fibrillation initiates following initial damage and softening. 𝐾𝑓 is the stiffness of a fibril, which is assumed to be mode -independent in the current implementation. T𝑚𝑓𝑢𝑙𝑡(𝜑) is the mode -dependent ultimate strength of a fibril, at which fibril damage initiates. By extension, the elastic limit of fibril extension is also mode -dependent ( 𝛿𝑚𝑓𝑢𝑙𝑡(𝜑)=T𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓). Finally , the mode -dependent parameter 𝛿𝑚𝑓∗(𝜑) governs the rate of softening i n the 𝐺𝑚𝑑 fibril damage/rupture region ( red). The interface strength is defined as a function of the mode as follows: T𝑚𝑚𝑎𝑥(𝜑)=𝜏𝑚𝑎𝑥− ( 𝜏𝑚𝑎𝑥−𝜎𝑚𝑎𝑥 1−exp(−𝜋 2 Ω𝑇) ) (1−exp(−𝜑 Ω𝑇)) (4) where 𝜏𝑚𝑎𝑥 is the mode II interface strength, Ω𝑇 is a CZM mode mixity parameter , and 𝜎𝑚𝑎𝑥 is the mode I interface strength. The mode mixity of the initial interface damage parameter 𝛿𝑚∗(𝜑) is obtained from 𝛿𝑚∗(𝜑)=𝐺̂𝑚(𝜑) T𝑚𝑚𝑎𝑥(𝜑)−T𝑚𝑚𝑎𝑥(𝜑) 2𝐾𝑚 (5) where 𝐺̂𝑚(𝜑)=𝐺𝑚𝑎(𝜑)+𝐺𝑚𝑏(𝜑) is the mode -dependent fracture energy prior to fibrillation , given as: 𝐺̂𝑚(𝜑) =1 2𝐾𝑚𝛿𝑚𝑒𝑙(𝜑)2+∫ T𝑚𝑚𝑎𝑥(𝜑)exp(−𝛥𝑚𝑚𝑎𝑥−𝛿𝑚𝑒𝑙 𝛿𝑚𝑑∗)𝑑Δ𝑚T𝑚𝑓𝑖𝑛𝑖𝑡(𝜑)/𝐾𝑓 T𝑚𝑚𝑎𝑥(𝜑)/𝐾𝑚 (6) We may specify the mode -dependence of 𝐺̂𝑚(𝜑) using the following function: 8 𝐺̂𝑚(𝜑)=𝐺𝑡0− ( 𝐺𝑡0−𝐺𝑛0 1−exp(−𝜋 2 Ω𝐺) ) (1−exp(−𝜑 Ω𝐺)) (7) where 𝐺𝑡0 is the mode II fracture energy prior to fibril formation and 𝐺𝑛0 is the mode I fracture energy prior to fibril formation , and the parameter Ω𝐺 sets the non -linearity of the transition from mode II to mode I . We specify the mode -dependence of 𝑇𝑚𝑓𝑢𝑙𝑡(𝜑) using the following function: T𝑚𝑓𝑢𝑙𝑡(𝜑)=𝜏𝑓𝑢𝑙𝑡− ( 𝜏𝑓𝑢𝑙𝑡−𝜎𝑓𝑢𝑙𝑡 1−exp(−𝜋 2 Ω𝑇𝑓) ) (1−exp(−𝜑 Ω𝑇𝑓)) (8) where 𝜏𝑓𝑢𝑙𝑡 is the mode II fibril strength , 𝜎𝑓𝑢𝑙𝑡 is the mode I fibril strength , and the parameter Ω𝑇𝑓 sets the non -linearity of the transition from mode II to mode I . We may consider the mode dependence of the fibril fracture energy 𝐺̂𝑚𝑓(𝜑) where we define the fracture energy as the stored elastic energy of the fibril and the fibril damage energy: 𝐺̂𝑚𝑓(𝜑)=∫ 𝐾𝑓𝛥𝑚𝑑𝑇𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓 𝑇𝑚𝑓𝑖𝑛𝑖𝑡(𝜑)/𝐾𝑓𝛥𝑚+∫ 𝐾𝑓𝛿𝑓𝑒𝑙(𝑒𝑥𝑝(−𝛥𝑚𝑚𝑎𝑥−𝑇𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓 𝛿𝑓∗))𝑑∞ 𝑇𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓𝛥𝑚 (9) or we may consider the mode dependence of the fibril fracture energy including only the fibril damage energy 𝐺̂𝑚𝑓(𝜑)=∫ 𝐾𝑓𝛿𝑓𝑒𝑙(𝑒𝑥𝑝(−𝛥𝑚𝑚𝑎𝑥−𝑇𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓 𝛿𝑓∗))𝑑∞ 𝑇𝑚𝑓𝑢𝑙𝑡(𝜑)/𝐾𝑓𝛥𝑚 (10) We specify the mode depe ndence of the fibril fracture energy 𝐺̂𝑚𝑓(𝜑) using the following function: 𝐺̂𝑚𝑓(𝜑)=𝐺t𝑑− ( 𝐺t𝑑−𝐺𝑛𝑑 1−exp(−𝜋 2 Ω𝐺𝑓) ) (1−exp(−𝜑 Ω𝐺𝑓)) (11) where 𝐺t𝑑 is the mode II fibril fracture energy, 𝐺n𝑑 is the mode I fibril fracture energy, and the parameter Ω𝐺𝑓 sets the non - linearity of the transition from mode II to mode I . Finally, we complete the description of the cohesive zone formulation by decompos ing 𝑇𝑚 into the normal and tangential components , 𝑇𝑛 and 𝑇𝑡, respectively, such that 9 𝑇𝑛={𝐾𝑛𝑜𝑐Δ𝑛,Δ𝑛<0 𝑇𝑚sin(𝜑),Δ𝑛≥0 (12) 𝑇𝑡=𝑇𝑚cos (𝜑) (13) where 𝐾𝑛𝑜𝑐 is the overclosure penalty stiffness. The EF -CZM is implemented in the finite element solver Abaqus through a User Defined Interface Subroutine (UINTER ). The consistent tangent matrix is derived analytically from (3). 2.2 EF-CZM parametric investigation of mixed mode behaviour Proportional loading path response : In Figure 4 we explore t he mixed mode coupling of the CZM under proportional loading paths. Figure 4(A) shows the traction -separation response of the EF -CZM for the case of 𝜎𝑚𝑎𝑥=𝜏𝑚𝑎𝑥 and 𝐺𝑛0=𝐺𝑡0, such that a mode -independent traction -separation relationship is obtained . Figure 4(B) shows the fracture energy as a function of 𝜑. In additio n to the total energ y, the individual energy contribution of each region (previously outlined in Figure 3) are also shown . Monotonic increases/decreases as a function o f mode angle are obtained for all energy contributions. In Figure 4(C-F) we consider cases in which the degree of fibrillation is mode dependant. Figure 4(C) we consider the case in which mode dependence is chosen so that self -similar scaling of the traction -separation relationship is obtained as a function of mode -angle (i.e. the shape of each curve is similar, with highest tractions and fibrillation for pure mode II separations). Such self -similar scaling is achieved by prescribing mode -dependence of fibril damage energy through equation (10). Fibrillation occurs at all mode -angles, which higher fibril tractions and stretching with increasing 𝜑. The resultant mode -dependence of total fracture energy, including normal and tangential contribution s, is shown in Figure 4(D). For the case shown in Figure 4(E-F) a different form of fibrillation mode dependence is shown. In this case, the mode -dependence of the total fibril energy is specified through equation (9). In this case, significant fibrillation is obtained only for mode angles close to mode II. For lower values of 𝜑, close to mode I, from the point of fibril initiation the fibrils immediately enter the damage regime, with the result that stretching of fibri ls at a constant stiffness is not computed. Again, the resultant mode -dependence of total fracture energy, including normal and tangential contributions, is shown in Figure 4(F). The three cases presented in Figure 4 merely serve to demonstrate the flexibility of the EF -CZM to represent a range of possible mode -dependent fibrillation r egimes, while providing physically appropriate mode dependent traction -separation relationships and fracture energies. 10 Figure 4. Response of EF -CZM to proportional loading paths with a constant mode angle 𝝋. (A)Traction -separa tion relationship for the case of 𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏=𝑮𝒕. (B) Computed fracture energies corresponding to (A). Individual energy contributions of each region of the traction -separation curve (as illustrated in Figure 3) are also shown; (C) Traction -separation relationship for the case of 𝟓𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎(𝝈𝒎𝒂𝒙 𝝉𝒎𝒂𝒙)𝟐 with mode dependence of fibril damage energy prescribed through equation (8); (D) Computed fracture energies corresponding to (C); (E) Traction -separation relationship for the case of 𝟓𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎= 𝑮𝒕𝟎(𝝈𝒎𝒂𝒙 𝝉𝒎𝒂𝒙)𝟐 with mode dependence of fibril total energy prescribed through equation (7); (F) Computed fracture energies corresponding to (E). Non-proportional loading path response : In this section the EF -CZM traction -separation response is examined in the case of a mode II separation to a specified finite val ue of Δ𝑡/𝛿𝑡, followed subsequently by complete normal separation , Δ𝑛=0→∞. Computed traction -separation relationships for a range of such non -proportional mixed -mode paths are shown in Figure 5 for the following cases: (A-B) 𝜎𝑚𝑎𝑥=𝜏𝑚𝑎𝑥,𝐺𝑛0=𝐺𝑡0 ; (D-E) 5𝜎𝑚𝑎𝑥=𝜏𝑚𝑎𝑥,𝐺𝑛0=𝐺𝑡0(𝜎𝑚𝑎𝑥 𝜏𝑚𝑎𝑥)2 . As shown in Figure 5(A), the occurrence of damage or fibrillation during the initial mode II separation depends on the prescribed value of Δ𝑡/𝛿𝑡. Traction - 11 separation during the subsequent normal separation is s hown in Figure 5(B). I f a fibril is formed during initial mode II separation, elastic fibril stretching continues immediately during subsequent normal separation. On the other hand, if damage and fibril initiation does not occur during the initial mode II separation, then damage is computed during the initial phase of the subsequent normal separation, followed by fibril initiation and fibril stretching. The computed fract ure energy as a function of the applied value of Δ𝑡𝑚𝑎𝑥/𝛿𝑡 is shown in Figure 5(C). While the non -linear influence of fibril formation and fibril stretching is comp uted, the total energy is mode independent, as prescribed, and monotonic increases/decreases in the no rmal and tangential contributions are observed. Corresponding non -proportional loading path responses are shown in Figure 5(D-F) for the case of 5𝜎𝑚𝑎𝑥=𝜏𝑚𝑎𝑥 and 𝐺𝑛0=𝐺𝑡0(𝜎𝑚𝑎𝑥 𝜏𝑚𝑎𝑥)2 . Mode mixity of the total fracture energy is specified using equation (9). As shown in Figure 5(D,E), the occurrence of fibrillation during the initial mode II separation results in reduced elastic fi bril deformation in the subsequent normal separation, with early onset of fibril damage and rupture. The total energy of separatio n is dominated by the mode II contribution during the initial phase of the loading path, as shown in Figure 5(F). Figure 5: Traction -separation response of the EF -CZM to non -proportional loading paths in which an initial mode II separation up to a specified finite value of 𝚫𝒕𝒎𝒂𝒙/𝜹𝒕𝒆𝒍 is followed by a subsequent normal separation to ultimate rupture. (A -C) shows res ults for the case of 𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎 in terms of, (A) tangential tractions, (B) normal tractions and (C) fracture energies. (D -F) shows results for the case of 𝟓𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎(𝝈𝒎𝒂𝒙 𝝉𝒎𝒂𝒙)𝟐 in terms of, (D) tangential tractions, (E) normal tractions and (F) fracture energies. Compliance of the EF -CZM with the requirement that instantaneous incremental dissipation is positive throughout a loading path is investigated in Figure 6. Instantaneous incremental dissipation is computed as 𝑑𝜙𝑖=0.5(𝑇𝑑Δ−Δ𝑑𝑇), in accordance with Cazes et al. (2009) . Analyses are presented for all proportional and non -proportional loading paths considered in Figure 4 and Figure 5, respectively . For all parameter sets, analyses reveal that 𝑑𝜙𝑖≥0 throughout all loading paths. 12 Figure 6. Instantaneous incremental dissipation ( 𝒅𝝓𝒊=𝟎.𝟓(𝑻𝒅𝚫−𝚫𝒅𝑻)) duing mixed mode loading paths. (A) Proportional loading paths for 𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎; (B) Proportional loading paths for 𝟓𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎(𝝈𝒎𝒂𝒙 𝝉𝒎𝒂𝒙)𝟐 ; (C) Non-proportional loading paths for 𝝈𝒎𝒂𝒙=𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎; (D) Non -proportional loading paths for 𝟓𝝈𝒎𝒂𝒙= 𝝉𝒎𝒂𝒙 𝐚𝐧𝐝 𝐚𝐧𝐝 𝑮𝒏𝟎=𝑮𝒕𝟎(𝝈𝒎𝒂𝒙 𝝉𝒎𝒂𝒙)𝟐 . Irreversible damage : The response of the EF -CZM to load -unload paths is explored in Figure 7. The response of the EF -CZM to mode I loading -unloading is shown in Figure 7(A). Two unloading paths are considered: one during the initial damage phase, and one during the final damage phase when the fibril s ultimate strength has been exceeded. In both cases the model imposes linear elastic unloading with the reduced modulus determined by the degree of damage . Negative values of Δ𝑛 are prevented by an overclosure penalty. As shown in Figure 7(B), positive instantaneous incremental dissipation is computed throughout all loading -unloading paths. The slightly more complex response of the model to mode II loading -unloading paths is shown in Figure 7(C). Tangential separation is increased such that the mode II strength is exceeded, and initial damage occurs. The direction of interface separation is then reversed and elastic unloading occurs at a reduced damage modulus , passing through the undeformed position (Δ𝑡=0), and continuing to deform elastically (Δ𝑡<0) until the previously computed damage stress magnitude is again reached. At this point, damage and softening resumes until the point of fibr il initiation. Elastic f ibril stretching then occurs until the maximum value of tangential traction is reached and fibril damage initiates. The direction of separation is again reversed and elastic fibril unloa ding-reloading occurs at a further reduced modulus, again until the reduced ultimate strength is reached and damage, and eventual fibril rupture, is computed. Positive instantaneous incremental dissipation is computed throughout this mode II loading -unload ing path ( Figure 7(D)). 13 Figure 7. Exploration of loading -unloading paths in the normal direction (A) and tangential direction (C). Positiv e instantaneous incremental dissipation for the normal loading path (B) and for the tangential loading path (D). 2.3 Simulation of artery crack propagation experiment using EF -CZM We next explore the ability of the EF -CZM to provide improved predictions (compared to a standard exponential damage CZM (recently implemented by FitzGibbon and McGarry (2020) and show n in Figure 1(D-F)). We focus specifically on the capability of the EF -CZM to predict the pattern of non -uniform crack propagation rate shown in Figure 1(F), i.e. initial slow crack propagation followed by fast crack propagation over a long distance, followed ultimately by slow crack propagation . A parametric investigation is performe d to explore the influence of the fibr il strength , the fibril fracture energy, and the maximum fibril extension on computed crack propagation rates . In all simulations presented in Figure 8(A-C) it can be noted that the EF - CZM introduces a strong non -uniform ity of crack propagation rate, with an initial slow crack rate followed by rapid crack growth, followed ultimately by slow crack growth. However, one significant limitation observed is that the formulation does not predict fast crack growth over a sufficiently long distance, and an alteration of model parameters does not facilitate an extension of this region. As shown in Figure 8(A), an increas e in the ultimate fibrillation strength merely delays the initiation of the fast crack growth phase but does not extend it over a longer distance. Figure 8(B) shows that the predicted crack - propagations rates are not sensitive to the fibril fracture energy . Figure 8(C) shows th at computed propagation rates are insensitive to the prescribed value of maximum fibril extension for values of 𝛿𝑓𝑒𝑙/𝛿𝑚𝑒𝑙≤75. Increasing the value of maximum fibril extension merely delays the onset of fast crack growth and sho rtens the distance over which it occurs. In summary, while the EF -CZM introduces highly non -uniform crack propagation, with slow crack growth predicted both before and after a phase of fast crack growth, the formulation is not capable of predicting a sufficiently long period of fast crack growth, compared to experimental observations. We suggest that this may be related to non -uniform fibrillation observed experimentally, with significantly lower fibrillation observed during fast crack growth. Therefore, i n the following section we propose a visco -plastic cohesive zone formulation in which fibrillation is dependent on the rate of interface separation. 14 Figure 8. Elastic fibrillation EF -CZM predicted crack length vs strain curves. (A) Effect of ultimate fibril strength (𝝉𝒇𝒖𝒍𝒕) on the predicted crack lengths. (B) Effect of fibril fracture energy (𝑮𝒕𝒅) on predicted crack lengths. (C) Effect of maximum fibri l extension (𝜹𝒇𝒆𝒍) on predicted crack lengths. 3 Visco -plastic cohesive zone model (VP -CZM) As demonstrated in the previous section, the EF -CZM does not accurately capture the non -uniform crack propagation rates reported for mode II artery dissection in the experimental study [21]. Specifically, while the EF -CZM correctly predicts initial and ultimate slow crack growth, it does not predict fast crack growth over a sufficiently long distance. We next propose a visco -plastic mode for fibrillation during crack growth in arteries based on the following experimental observation s from FitzGibbon and McGarry [21]: (i) fibril s do not merely deform elastically; unloading of the visible fibril s does not lead to fibril contraction and removal of the mode II separation between the fracture surfaces . Instead fibrils are observed to have undergone pull -out. Even fibril s near the crack -tip that have not ruptured and are attached to opposing fracture surfaces , are observed to have permanently “pulled - out” of the bulk material and do not exhibit any stored elastic strain energy ( Figure 1(h)). We hypothesise that a plasticity mechanism should be used to capture the dissipative process of fibril pull -out and associated permanent fibril deformation. (ii) Higher levels of fibrillation are observed during the initial and final phases of slow crack growth. This suggests a link between the rate of interface separation and fibril lation initiation. We hypothesise that a visco -plasticity mechanism can capture this phenomenon, whereby high rates of interface separation will delay the onset of plastic yiel ding (fibril pull -out) and result in higher fibril tractions so that the fibril ultimate strength is reached without significant plasticity (pull -out). Contrastingly, this hypothesis suggests that at slow interface separation rates fibril s will undergo ear ly yield and significant plasticity (pull - out) before the ultimate fibril strength is reached. Essentially, a higher level of fibril plasticity (pull -out) is hypothesised to 15 provide a significant dissipative mechanism, providing a link between non -uniform crack growth rates and non -uniform levels of fibrillation. In Figure 9 we present a schematic of the proposed visco -plastic cohesive zone model (VP -CZM) . A high interface separation - rate is shown in grey and a low interface separation -rate is shown in black. The initial regime of elastic deformation entails a linear increase in traction with increasing separation, characterised by the interface intrinsic elastic stiffness 𝑘𝑚. Plastic yield ing and hardening is governed by both the interface separation and separat ion-rate. When the ultimate stress (𝜎𝑢𝑡𝑠) is reached an exponential damage phase initiates, resulting in ultimate interface rupture. Figure 9. Schematic of the visco -plastic cohesive zone (VP -CZM) formulation 3.1 VP-CZM formulation We propose a simple 1D elasto -visco -plastic formulation to describe elast o-plastic axial deformation of a fibril . For a given interface displacement vector 𝚫, the fibril axial deformation is given as the separation magnitude, Δ𝑚. To describe permanent fibril extension due to pull -out we perform an additive decomposition of Δ𝑚 into an elastic recoverable separation , Δ𝑚𝑒𝑙, and an irrecoverable plastic (pull -out) separation , Δ𝑚𝑝 (based on the approach of Cazes et al. [31]: Δ𝑚 𝛿𝑚𝑒𝑙=Δ𝑚𝑒𝑙 𝛿𝑚𝑒𝑙+Δ𝑚𝑝 𝛿𝑚𝑒𝑙 (14) where 𝛿𝑚𝑒𝑙 is a characteristic length. For convenience, and in keeping with standard plasticity notation , in the remainder of this section we refer to these non-dimensional fibril separations as effective fibril strains (elastic and plastic): 16 Δ𝑚 𝛿𝑚𝑒𝑙=𝜀=𝜀𝑒𝑙+𝜀𝑝 (15) The axial fibril stress is obtained from the elastic effective fibril strain, such that 𝜎𝑓=𝑘𝑚(𝜀−𝜀𝑝) (16) Strain-rate dependent visco -plasticity is implemented through the following relationship: 𝜀̇𝑝=𝜀̇0(𝜎𝑓 𝜎𝑓ℎ(𝜀𝑝))𝑚 (17) where 𝜀𝑝̇ is the effective fibril axial plastic strain rate and 𝜎𝑓 is the axial fibril stress. T he parameter 𝜀̇0 is a reference strain rate, 𝑚 is the strain -rate hardening exponent, and 𝜎𝑓ℎ(𝜀𝑝) is the corresponding fibril yield stress (as a function of effective fibril axial plastic strain) at the reference strain -rate 𝜀̇0. The form of 𝜎𝑓ℎ(𝜀𝑝) is specified through the following strain hardening law: 𝜎𝑓ℎ(𝜀𝑝)=𝜎𝑓𝑦0+(𝜎𝑓𝑦∞−𝜎𝑓𝑦0)(1−𝑒𝑥𝑝(−𝜀𝑝 𝜀𝑛) ) (18) where 𝜎𝑓𝑦∞ is the saturation yield stress, 𝜎𝑓𝑦0 is the initial yield stress, and 𝜀𝑛 is a hardening parameter . We assume a critical - stress failure initiation criterion, whereby fibril rupture initiates when the fibril stress reaches a critical value; 𝜎𝑓=𝜎𝑓𝑢𝑡𝑠 . Fibril rupture is described using a standard exponential damage formulation: 𝜎𝑓=𝜎𝑓𝑢𝑡𝑠exp(−𝜀𝑚𝑎𝑥−𝜀𝑢𝑡𝑠 𝜀𝑓∗) (19) where 𝜀𝑢𝑡𝑠 is the effective fibril axial strain at which 𝜎𝑓𝑢𝑡𝑠 is reached, 𝜀𝑚𝑎𝑥 is the maximum strain reached, and the parameter 𝜀𝑓∗ governs the rate of damage/softening . Numerical implementation: Here we outline the numerical implementation of the visco-plastic fi brillation model. The total stress (𝜎𝑓𝑖) is expressed as the sum of the incremental increase in stress (Δ𝜎𝑓) and the stress at the previous increment (𝜎𝑓𝑖−1): 𝜎𝑓𝑖=𝜎𝑓𝑖−1+Δ𝜎𝑓 (20) The incremental increase in stress is given as: Δ𝜎𝑓=𝑘𝑚(Δ𝜀−Δ𝜀𝑝) (21) where Δ𝜀 is the increment al total strain and Δ𝜀𝑝 is the incremental change in plastic strain. The change in plastic strain ( Δ𝜀𝑝) is calculated through a Newton -Raphson scheme. The initial value/ guess for Δ𝜀𝑝 in the scheme is estimated using a rate tangent method. Briefly, the plastic strain increment is given as: 17 d𝜀𝑝=𝜕𝜀̇𝑝 𝜕𝜎𝑓Δ𝜀𝑘𝑚𝜃Δ𝑡𝑖𝑚𝑒+Δ𝑡𝑖𝑚𝑒𝜀̇𝑝 𝜕𝜀̇𝑝 𝜕𝜎𝑓𝑘𝑚𝜃Δ𝑡𝑖𝑚𝑒−𝜕𝜀̇𝑝 𝜕𝜎𝑓ℎ(𝜀𝑝)d𝜎̇𝑓ℎ(𝜀𝑝)𝜃Δ𝑡𝑖𝑚𝑒+1 (22) Where Δ𝑡𝑖𝑚𝑒 is the current time increment , and 𝜃 determines the rate tangent solution method (a value of 𝜃=0.5 is used (Crank –Nicolson method )). The following expressions can be derived from equations ( 16)-(18) 𝜕𝜀̇𝑝 𝜕𝜎𝑓=𝜀̇0 𝑚 𝜎𝑓(𝑚−1) 𝜎𝑓ℎ(𝜀𝑝)−𝑚 (23) 𝜕𝜀̇𝑝 𝜕𝜎𝑓ℎ=−𝜀̇0 𝑚 (𝜎𝑓𝜎𝑓ℎ(𝜀𝑝)−1)𝑚 𝜎𝑓ℎ(𝜀𝑝)−1 (24) d𝜎̇𝑓ℎ(𝜀𝑝)=(𝜎𝑓∞−𝜎𝑓𝑦 𝜀𝑛)exp(−𝜀𝑝 𝜀𝑛) (25) Using the rate tangent solution as an initial guess, we can then iterate using a Newton -Raphson scheme using the following function: 𝑓(𝑥)=Δ𝜀𝑝−(1−𝜃)Δ𝑡𝑖𝑚𝑒𝜀̇−𝜃Δ𝑡𝑖𝑚𝑒𝜀̇0(𝜎𝑓𝑖−1+Δ𝜎𝑓 𝜎𝑓ℎ(𝜀𝑝)𝑖−1+Δ𝜎𝑓)𝑚 (26) where 𝑖 denotes the iteration number , and 𝑖−1 is the previous iteration . 𝑓′(𝑥) is given as : 𝑓′(𝑥)=1−𝑚𝜃Δ𝑡𝑖𝑚𝑒𝜀̇0(𝜎𝑓𝑖−1+Δ𝜎𝑓𝑖−1 𝜎𝑓ℎ(𝜀𝑝)𝑖−1)𝑚−1 ∗ ( 𝜎𝑓ℎ(𝜀𝑝)𝑖−1(−𝑘𝑚)−(𝜎𝑓𝑖−1+Δ𝜎𝑓)(𝜎𝑓∞−𝜎𝑓𝑦) (−𝜀𝑛)exp(−𝜀𝑖−1𝑝+Δ𝑡𝑖𝑚𝑒Δ𝜀 𝜀𝑛) 𝜎𝑓ℎ(𝜀𝑝)𝑖−12 ) (27) and the updated increment of plastic strain is then given as: Δ𝜀𝑝𝑖+1=Δ𝜀𝑝−𝑓(𝑥) 𝑓′(𝑥) (28) This formulation for axial elasto -visco -plastic fibril deformation is generalised within the framework of a CZM such that : 18 𝑇𝑛=(𝜎𝑓)sin(𝜑); 𝑇𝑡=(𝜎𝑓)cos (𝜑) (29) Mode -dependence of fibril visco -plasticity is specified through the following expressions for 𝜎𝑓𝑦,𝜎𝑓𝑢𝑡𝑠, and 𝜀𝑓∗ as a function of mode angle, 𝜑 : 𝜎𝑓𝑦(𝜑)=𝜏𝑦− ( 𝜏𝑦−𝜎𝑦 1−exp(−𝜋 2 Ωy) ) (1−exp(−𝜑 Ωy)) (30) where 𝜏𝑦 and 𝜎𝑦 are the values of initial yield stress under mode II and mode I loading , respectively , and the parameter Ω𝑦 sets the non -linearity of the transition from mode II to mode I . Similarly, the mode -dependence of the ultimate strength of fibrils is given as: 𝜎𝑓𝑢𝑡𝑠(𝜑)=𝜏𝑢𝑡𝑠− ( 𝜏𝑢𝑡𝑠−𝜎𝑢𝑡𝑠 1−exp(−𝜋 2 Ωuts) ) (1−exp(−𝜑 Ωuts)) (31) where 𝜏𝑢𝑡𝑠 and 𝜎𝑢𝑡𝑠 are the values of initial yield stress under mode II and mode I loading, respectively , and the parameter Ω𝑢𝑡𝑠 sets the non -linearity of the transition from mode II to mode I . The description of the mode mixity of 𝜀𝑓∗(𝜑) is obtained from 𝜀𝑓∗(𝜑)=𝐺̂𝑚(𝜑) 𝜎𝑓𝑢𝑡𝑠(𝜑)−𝜎𝑓𝑢𝑡𝑠(𝜑) 2𝐾𝑚 (32) where 𝐺̂𝑚(𝜑) is the mode -dependent fibril fracture energy 𝐺̂𝑚(𝜑)=∫𝜎𝑓𝑢𝑡𝑠exp(−𝜀𝑚𝑎𝑥−𝜀𝑢𝑡𝑠 𝜀𝑓∗)∞ 𝜀𝑢𝑡𝑠d𝜀 (33) The mode dependence of the fibril fracture energy 𝐺̂𝑚(𝜑) is described using the following expression : 𝐺̂𝑚(𝜑)=𝐺𝑡𝑑− ( 𝐺𝑡𝑑−𝐺𝑛𝑑 1−exp(−𝜋 2 ΩGf) ) (1−exp(−𝜑 ΩGf)) (34) where 𝐺𝑡𝑑 is the fibril fracture energy in the tangential direction, 𝐺𝑛𝑑 is the fibril fracture energy in the normal direction, and the parameter Ω𝐺𝑓 sets the non -linearity of the transition from mode II to mode I . The VP -CZM is implemented in the finite element solver Abaqus through a User Defined Interface Subroutine (UINTER ). The consistent tangent matrix is derived estimated from (22). 19 3.2 Thermodynamic consideration in relation to a VP -CZM From Gurtin (1979) , the first and second principles of thermodynamics are expressed as: d𝑒𝑆 d𝑡=Tsd⟦Δ⟧ d𝑡−⟦𝑞⟧𝑛 (35) d𝑠𝑠 d𝑡=−⟦𝑞⟧𝑛 H+d(𝑠𝑖)𝑠 d𝑡 (36) where 𝑒𝑠,𝑠𝑠 and (𝑠𝑖)𝑠 are the surface densities of internal energy, entropy, and internal entropy, respectively. ⟦Δ⟧ represents the interface separation vector at the interface, Ts is the stress vector defined over the di scontinuity surface, ⟦𝑞⟧ is the heat flow jump, 𝑡 is the time and H is the absolute temperature . The surface energy dissipated by the cohesive zone is then: d𝜙𝑠=Hd(𝑠𝑖)𝑠 (37) Using equations (35) and (36) one can obtain the following expression for the free surface energy: 𝜓𝑠=𝑒𝑠−H𝑠𝑠 (38) This leads to the expression (originally derived by Cazes et al. [31]) for an increment of dissipated energy, given as: 𝑑𝜙𝑠=Tsd⟦Δ⟧−d𝜓𝑠−𝑠𝑠dH (39) The total interface separation vector ⟦Δ⟧ is divided into elastic and plastic components. ⟦Δ⟧=⟦Δ⟧𝑒𝑙+⟦Δ⟧𝑝 (40) The free surface energy can trivially be defined as: 𝜓𝑠=1 2Tsd⟦Δ⟧𝑒𝑙 (41) Introducing equation (41) into equation (39), an expanded expression for an increment of dissipated energy is obtained, such that: d𝜙𝑠=Tsd⟦Δ⟧𝑝+1 2(Tsd⟦Δ⟧𝑒𝑙−⟦Δ⟧𝑒𝑙dTs) (42) This can be rewritten in terms of the contribution of the evolution of the plastic interface separation ⟦Δ⟧𝑝 and of the elastic interface separation ⟦Δ⟧𝑒 as follows: d𝜙𝑠=d𝜙𝑠𝑝+d𝜙𝑠𝑒𝑙 (43) Where the dissipated energy due to the plastic interface separation ⟦Δ⟧𝑝 is given as: 20 d𝜙𝑠𝑝=Tsd⟦Δ⟧𝑝 (44) and the dissipated energy due to the elastic interface separation ⟦Δ⟧𝑒𝑙 is given as: d𝜙𝑠𝑒𝑙=1 2(Tsd⟦Δ⟧𝑒𝑙−⟦Δ⟧𝑒𝑙dTs) (45) The Clausius -Duhem inequality states that, in a thermodynamically allowable material, the dissipation is always positive. This requires the VP-CZM to exhibit instantaneous positive incremental dissipation over the entirety of a loading path. This leads to the following expr ession to estimate the thermodynamic legitimacy of the proposed model: d𝜙𝑠= Tsd⟦Δ⟧𝑝+1 2(Tsd⟦Δ⟧𝑒𝑙−⟦Δ⟧𝑒𝑙dTs)≥0 (46) 3.3 VP-CZM parametric investigation of mixed -mode behaviour The VP -CZM is shown for a range of axial fibril strain rate s (𝜀̇) in Figure 10. For high strain rates the VP-CZM behaviour is similar to a typical exponential -type model with negligible amounts of plasticity prior to reaching t he ultimate strength. In contrast, for low strain rates the model predicts lower initial yielding, low hardening rates, and, consequently, high levels of plastic deformation prior to reaching the ultimate strength. This predicts that significantly higher l evels of plastic deformation (fibril pull -out) and dissipation occur during low rates of interface separation. Figure 10. Traction as a function of separation for a range of strain rates. 𝜺̇𝟎 is a VP -CZM parameter which controls the strain rate threshold. Fast strain rates exhibit negligible amounts of fibril pull -out, slow strain rates result in significant levels of fibril pull -out as indicated by the black curve. Non-proportional loading p ath response : The mixed-mode behaviour of the VP -CZM is next explored by considering non - proportional loading paths subject to a range of strain rates for each path. The non -proportional loading paths consist of a n initial mode II displacement to a prescri bed value (Δ𝑡𝛿𝑡𝑒𝑙⁄=0→10) followed by complete normal separation (Δ𝑛→∞). 21 We consider three cases in terms of imposed non -uniform strain rates: 𝜀̇𝑛=𝜀̇𝑡; 𝜀̇𝑛>𝜀̇𝑡; 𝜀̇𝑡>𝜀̇𝑛. In Figure 11 we consider the case of mode -independent material parameters and we present computed traction separation relationships for the three aforementioned non -uniform strain rates ((A, D, G) Tangential components; (B, E, H) normal components; (C, F, I) magnitudes (effective fibril axial traction -separation relationship)). In all case s if significant plastic deformation occurs during the initial mode II separation, the computed leve l of plasticity (fibril pull -out) during subsequent normal separation is less pronounced. In terms of the effective axial fibril behaviour (traction -separation magnitude curves), transition from the mode II phase to the subsequent normal separation phase e ntails a temporary reduction in the strain -rate magnitude. This results in a reduction in the effective yield stress (in accordance with equation (17)) as normal separat ion increases , 𝜀̇→𝜀̇𝑛, and the final hardening curve prior to rupture depend on the value of 𝜀̇𝑛. For example, in Figure 11(C) the hardening curve for the normal sepa ration is lower than that for the initial mode II separation because 𝜀̇𝑡>𝜀̇𝑛. In the case of the two highest prescribed values of initial mode II separation, fibril rupture initiates prior to the transition to normal separation. Figure 11. Normalised traction separation curves for the VP -CZM undergoing non -proportional loading paths in the case of 𝜺̇𝒏= 𝜺̇𝒕=𝜺̇𝟎 (A-C), 𝜺̇𝒏=𝟑𝜺̇𝒕=𝜺̇𝟎 (D-F), 𝟑𝜺̇𝒏=𝜺̇𝒕=𝜺̇𝟎 (G-I). Other VP -CZM parameters are as follows: 𝒌𝒎=𝟏𝟎𝑴𝑷𝒂,𝝈𝒚=𝝉𝒚,𝑮𝒏𝟎= 𝑮𝒕𝟎,𝝈𝒖𝒕𝒔=𝝉𝒖𝒕𝒔,𝝈𝒖𝒕𝒔=𝟏.𝟐𝟓𝝈𝒚. Proportional loading path response : The VP -CZM is now explored under proportional loading paths subject to a constant separation rate (𝜀̇𝑚𝜀̇0⁄=2.5) for each path. Figure 12(A) confirms that the effective fibril axial traction -separation relationship is mode -independent when 𝜎𝑦=𝜏𝑦,𝐺𝑛0=𝐺𝑡0. Figure 12(B) shows the computed fracture energ y is also mode - independent . The stored elastic energy is also shown in Figure 12(B). Tangential and normal components of the traction - separation relationships are shown in Figure 12(C,D) . Corresponding relationships are shown in Figure 12(E-H) for the case s where the values of mode II yield stress and ultimate strength are higher than the respective mode I values. The VP -CZM formulation produces consistent yield, hardening, and rupture behaviour for all mode angles. Furthermore, monotonic 22 increases/decreases in the tangential/normal components of fracture energy are computed as a function of mode angle ( Figure 12(F)). Normal and tangential traction components also exhibit consistent behaviour ( Figure 12(G, H)). 23 24 Figure 12. Traction separation response for the proportional loading paths subject to a constant separation rate (𝜺̇𝒎𝜺𝒓𝟎𝒑⁄=𝟐.𝟓) for the VP-CZM with : (A-D): 𝝈𝒚=𝝉𝒚,𝑮𝒏𝟎=𝑮𝒕𝟎, (E-H): 𝟓𝝈𝒚=𝝉𝒚,𝟓𝑮𝒏𝟎=𝑮𝒕𝟎. In all cases 𝑻𝒎𝒖𝒕𝒔=𝟏.𝟐(𝑻𝒎𝒚𝟎 ). The calculated fracture energy is show in the case of: (B) 𝝈𝒚=𝝉𝒚,𝑮𝒏𝟎=𝑮𝒕𝟎 and (F) 𝟓𝝈𝒚=𝝉𝒚,𝟓𝑮𝒏𝟎=𝑮𝒕𝟎. The total fracture energy is shown as well as the indiv idual energy contributions of the elastic recoverable energy (𝑮𝒆𝒍) and the plastic irrecoverable (dissipated) energy (𝑮𝒑). The instantaneous incremental dissipation for the visco -plastic CZM is explored below in Figure 13 in the case of proportional loading (A, C, E) and non -proportional loading (D, D, F) in the case of 𝜎𝑦=𝜏𝑦,𝐺𝑛=𝐺𝑡 (A, B), 5𝜎𝑦=𝜏𝑦,5𝐺𝑛0=𝐺𝑡0 (C, D), and 𝜎𝑦=5𝜏𝑦,𝐺𝑛0=5𝐺𝑡0 (E, F ). Instantaneous incremental dissipation is positive throughout all loading scenarios explored; this satisfies the Clausium -Duhem inequality (equation (46)). 25 Figure 13. Demonstration of positive instantaneous incremental dissipation throughout loading history for the VP -CZM undergoing (A,C,E) Proportional loading and; (B,D,F) Non -proportional loading when; (A -B) 𝝈𝒚=𝝉𝒚,𝑮𝒏𝟎=𝑮𝒕𝟎 and; (C -D) 𝟓𝝈𝒚=𝝉𝒚,𝟓𝑮𝒏𝟎= 𝑮𝒕𝟎, and; (E -F) 𝝈𝒚=𝟓𝝉𝒚,𝑮𝒏𝟎=𝟓𝑮𝒕𝟎. Irreversible damage and plastic deformation : The response of the VP -CZM to load -unload paths is explored in Figure 14. Figure 14(A) shows a mode I load -unload path (the path follows the red arrows, blue arrows, then green arrows). An illustration of the fibrillation process is shown f or a mode I load -unload scenario (i -viii). Collagen fibrils are initially embedded in the tissue matrix (i), upon stretching the fibrils stiffen and yield (ii -iii) with increasing interface separation. Fibrils then begin to plastically pull out of the tiss ue matrix (iv). Upon unloading, fibrils recover elastically (v) with significant permanent fibril pull-out (vi). Upon reloading, fibrils stretch elastically until yield is reached and begin to pull -out permanently until the point of rupture initiation (vii ). Upon rupture, fibrils undergo total fibril pull -out or fibril damage (viii). The tangential traction - 26 separation response to the same loading regime in mode II is shown in Figure 14(B) with positive instantaneous incremental dissipation demonstrated throughout (Figure 14(C)). Figure 14(D) shows traction -separation curves for a range of mixed mode proportional loading for the load -unload regime outlined in (A). Figure 14(E) shows the tangential traction -separation response of mixed -mode proportional loading for the load -unload regime outlined in (B). Figure 14(D-E) demonstrate c onsistent behaviour with (A -B) in mixed mode scenarios. Figure 14. (A) Demonstration of permanent plastic deformation and fibril pull -out for a mode I loading -unloading scenario; (B) Demonstration of permanent plastic deformation and fibril pull -out for a mode II loading -unloading scenario following the arrows 27 in order of (red, blue, green); (C) demonstration of positive instantaneous incremental energy dissipation; (D) Normal traction - separation curves in mixed mode loading unloading scenarios following the path outlined in (A); (E) Tangential traction -separation curves in mixed mode loading unloading scenarios following the path outlined in (B). 3.4 Simulation of artery crack propagation experiment using VP -CZM In Figure 15 we assess th e ability of the VP -CZM to predict non -uniform crack propagation rates in arteries. The m odel predictions are superimposed upon experimental results from FitzGibbon and McGarry [21] (presented in terms of measured crack length as a function of applied curcumferential strain in the artery ). Three values of 𝜀̇0 are presented (𝜀̇0=0.1,𝜀̇0= 0.01,𝜀̇0=0.001). Similar predictions are obtained in all cases. An initial phase of slow crack growth is computed, followed by a phase of fast crack growth over a long distance, followed by a phase of slower crack growth. Figure 15(I) show s that significant plastic deformation (fibril pull -out) results in initially show crack growth. During the fast crack growth phase, as the ultimate fibril strength being reached without significant plasticity (Figure 15(II)) the effective VP-CZM behaviour in this region resembles a standard elastic -exponential damage formulation. Finally, after significant crack growth, plastic deformation re -emerges ( Figure 15(III), resulting in a slowing of the rate of crack growth. Figure 15: Predicted crack lengths superimposed over the experimentally measured crack lengths as a function of the circumferential strain for three values of 𝜺̇𝟎 (𝜺̇𝟎=𝟎.𝟏,𝜺̇𝟎=𝟎.𝟎𝟏,𝜺̇𝟎=𝟎.𝟎𝟎𝟏). Figure 16 shows the results of a finite element cohesive zone analysis of the experiment s of FitzGibbon and McGarry [21]. Traction as a function of separation is shown for three phases along the cohesive interface throughout the simulation. Extensive plastic deformation is observed at phase I where the strain rate is lowest . Node A undergoes significant plastic deformation due to the slow strain rate, it is foll owed by Node B which undergoes the most plastic deformation, followed by nodes C and D which undergo little plastic deformation and no plastic deformation, respectively due to the increasing strain rate. Phase II shows three nodes (E -G), all of which under go no plastic deformation as the strain rate results i n rapid fracture. T he final phase (III) shows three nodes (H -J) which increase from no plastic deformation to a moderate amount of plastic deformation as the strain rate begins to slow. These three phas es mirror the trends observed in the experiment : a prolonged region of slow growth 28 due to fibrillation (I), followed by a prolonged region of rapid crack growth (II), followed by a region of slow crack growth (III). This demonstrates the capability of the VP-CZM approach to generate a wide range of effective interface behaviours, resulting in non -uniform crack growth and a non -uniform prediction of plasticity (fibril pull -out). Figure 16. Normalised traction as a function of separation for three crack growth regimes selected along the cohesive zone interface . (I) Initial phase of crack growth: significant plasticity followed by rapid crack growth; (II) Rapid crack growth phase; (III) Final phase of increasing plasticity. 4 Concl uding remarks In the present study we propose two phenomenological CZM formulations to represent fibrillation during dissection of soft tissue. We first consider pure elastic fibrillation (EF -CZM). In this formulation fibrillation initiates during an initi al phase of interface damage/softening. Fibrils then deform elastically until an ultimate strength is reached, followed by fibril rupture . Simulations reveal that such elastic fibrillation does not accurately predict the highly non -uniform crack propagatio n rates measured experimentally for arterial tissue. We then propose a phenomenological visco -plastic cohesive zone formulation for fibrillation (VP -CZM). This approach is motivated by the experimental observation that fibrils undergo partial pull -out (in addition to elastic stretching) during dissection of arteries. Simulations reveal that this visco -plastic cohesive zone formulation provides a reasonable prediction of non -uniform crack propagation rates in arteries. Significant plasticity (representing fi bril pull-out) during fracture initiation results in an initial slow phase of crack growth. The transition to fast crack growth is facilitated by the prediction of fibril rupture prior to extensive plasticity. The current implementation of fibril visco -plastic behaviour in a cohesive zone framework is limited to a phenomenological representation of axial fibril deformation and axial fibril pull -out. However, we provide extensive parametric exploration of the model behaviour under a range of mixed -mode loadi ng paths, representing fibril reorientation in addition to fibril extension. We demonstrate that the formulation produces consistent behaviour for all loading paths and that positive incremental instantaneous dissipation is computed in all cases. The VP-CZM represents an advance on standard elastic -damage cohesive zone formulations. While the VP -CZM effectively produces such standard CZM behaviour at high interface separation rates, it predicts plastic deformation and energy dissipati on (representing fibril pull-out) at low interface separation rates. 29 Previous CZMs have incorporated viscoplasticity to model crazing in polymers [16,17] . The role of bridging fibres has also previously been examined in materials such as cellulose nanopaper [33]. The current study develop s the first visco -plastic CZM for fibril pull out during soft tissue fracture. A detailed analysis of visco -plastic mixed mode fibrillation for a range of mode I and mode II yield stresses and ultimate stresses. Positive i nstantaneous incremental dissipation and compliance with the Clausius -Duhem inequality is demonstrated for a wide range of proportional and non -proportional mixed mode loadi ng paths . Future models of artery dissection should incorporate rate-dependence of the bulk arterial tissue, in addition to investigating the role of active contractility of smooth muscle cells, fibroblasts, and endothelial cells on crack initiation and pr opagation [34,35] . Residual strain due to elastin should also be considered. Residual strain has been shown to be heterogeneously distributed axially along the length of the aorta, with additional variations between layers of the artery wall [36–38]. In addition to the demonstrated application of artery dissection, the VP -CZM framework could be used to simulate fibre alignment and fibrillation as a toughening mechanism in a range of soft tissues and tissue engineered collagen construc ts [19,20,39,40] . The framework can also be used to predict fibre bridging and mode -dependant crack kinking cracking in bone [41–45]. Furthermore, the consistent load -unload behaviour of the proposed visco -plastic CZM could provide improved predictions of fatigue and fretting failure of metals [46,47] . 30 5