Viscoplastic cohesive zone modelling of interfacial fibrillation in soft tissue

In the present study we propose two phenomenological CZM formulations to represent fibrillation during dissection of soft tissue. We firstly consider elastic fibrillation. 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 fibrillation 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 observation 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 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 loading 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 dissipation (representing 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


General Background
Cohesive zone models are commonly used to model the delamination of an interface [1][2][3][4][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 models (CZMs) have been applied to model interfacial fibrillation in a variety of engineered materials including polymers, elastomers, paper, and plastics [11][12][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 displacement problems, such as fibrillation [10,14].This model was implemented by Vossen et al. as 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 orientations, with debonding of a fibril based on the mechanics of peeling of an elastic film from a substrate [15].Viscoplastic CZM models of crazing 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 cohesive zone formulations have previously been used to model fracture in biological tissues such as arteries [22][23][24], and cells [25], but viscoplastic toughening mechanisms have not been previously considered.

Examination of fibrillation during mode II delamination of arteries uncovered in recent experimental study (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 uniform and that significant fibrillation occurs, both during initiation and propagation (see Figure 1).While a detailed description of the novel mode II soft tissue fracture 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 along 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 radial 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 have undergone partial pull-out, such that removal of elastic strain energy does not result in significant fibril shortening and closure of crack surfaces.Finally, 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 (GII) 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.05N/mm captures only the slow growth rate at the start of the experiment, whereas GII=0.0005N/mm captures the fast crack growth rate in regime (ii), but not the slow crack growth in regimes (i) and (iii) (Figure 1(f)).

Structure and aims of the study
In the current study we aim to develop new CZM approaches 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).Following 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 arteries.Furthermore, analysis of our experimental fracture surfaces suggests that fibrils 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.

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 ii.
Analysis of the finite element implementation of the model to simulate previously reported mode II fracture experiments [21].
Mixed mode separation-specified loading paths: Proportional and non-proportional loading paths are considered in the present 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 nonproportional loading; an initial tangential separation (Δ  = Δ   ) followed by a subsequent normal separation to failure (Δ  = ∞), as shown in the illustration in Figure 2(b).For the remainder of the present study this path is referred to as non-proportional loading.Model 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 interface using user-defined interfacial constitutive behaviour 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 exponential constitutive laws typically used for arterial tissue [27][28][29]..The isochoric strain energy density function is given as: 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 C 0 and C 1 continuity the following relations must be enforced: 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).

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 to 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 ultimate 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 fibrils 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 stiffness (  ), 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 nonuniform crack propagation rates observed experimentally in our artery fracture experiments.

EF-CZM formulation
For a given interface displacement vector  with normal and tangential (shear) components Δ  and Δ  , respectively, the corresponding displacement magnitude is given as 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: 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 separation 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 in the    fibril damage/rupture region (red).The interface strength is defined as a function of the mode as follows: 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 where  ̂() =    () +    () is the mode-dependent fracture energy prior to fibrillation, given as: We may specify the mode-dependence of  ̂() using the following function: 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: 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: or we may consider the mode dependence of the fibril fracture energy including only the fibril damage energy We specify the mode dependence of the fibril fracture energy  ̂  () using the following function: where  t  is the mode II fibril fracture energy,  n  is the mode I fibril fracture energy, and the parameter Ω   sets the nonlinearity of the transition from mode II to mode I.
Finally, we complete the description of the cohesive zone formulation by decomposing   into the normal and tangential components,   and   , respectively, such that 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).

EF-CZM parametric investigation of mixed mode behaviour
Proportional loading path response: In Figure 4 we explore the 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 addition to the total energy, the individual energy contribution of each region (previously outlined in Figure 3) are also shown.Monotonic increases/decreases as a function of 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 contributions, 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 fibrils 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 regimes, while providing physically appropriate mode dependent traction-separation relationships and fracture energies.

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 value 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  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.

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 shown 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 performed to explore the influence of the fibril 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-uniformity 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 increase 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 crackpropagations rates are not sensitive to the fibril fracture energy.Figure 8(C) shows that 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 shortens 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, in the following section we propose a visco-plastic cohesive zone formulation in which fibrillation is dependent on the rate of interface separation.

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 observations from FitzGibbon and McGarry [21]: (i) fibrils do not merely deform elastically; unloading of the visible fibrils 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 fibrils near the crack-tip that have not ruptured and are attached to opposing fracture surfaces, are observed to have permanently "pulledout" 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 fibrillation initiation.We hypothesise that a visco-plasticity mechanism can capture this phenomenon, whereby high rates of interface separation will delay the onset of plastic yielding (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 fibrils will undergo early yield and significant plasticity (pullout) before the ultimate fibril strength is reached.Essentially, a higher level of fibril plasticity (pull-out) is hypothesised to 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 separationrate 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 yielding and hardening is governed by both the interface separation and separation-rate.When the ultimate stress (  ) is reached an exponential damage phase initiates, resulting in ultimate interface rupture.

VP-CZM formulation
We propose a simple 1D elasto-visco-plastic formulation to describe elasto-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]: 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): The axial fibril stress is obtained from the elastic effective fibril strain, such that Strain-rate dependent visco-plasticity is implemented through the following relationship: where   ̇ is the effective fibril axial plastic strain rate and   is the fibril stress.The 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: where    ∞ is the saturation yield stress,    0 is the initial yield stress, and   is a hardening parameter.We assume a criticalstress 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: 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 fibrillation model.The total stress (   ) is expressed as the sum of the incremental increase in stress (Δ  ) and the stress at the previous increment (  −1 ): The incremental increase in stress is given as: where Δ is the incremental 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: 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) Using the rate tangent solution as an initial guess, we can then iterate using a Newton-Raphson scheme using the following function: where  denotes the iteration number, and  − 1 is the previous iteration. ′ () is given as: and the updated increment of plastic strain is then given as: This formulation for axial elasto-visco-plastic fibril deformation is generalised within the framework of a CZM such that: Mode-dependence of fibril visco-plasticity is specified through the following expressions for    ,    , and   * as a function of mode angle,  : 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: 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 where  ̂() is the mode-dependent fibril fracture energy The mode dependence of the fibril fracture energy  ̂() is described using the following expression: 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).

Thermodynamic consideration in relation to a VP-CZM
From Gurtin (1979), the first and second principles of thermodynamics are expressed as: where   ,   and (  )  are the surface densities of internal energy, entropy, and internal entropy, respectively.⟦Δ⟧ represents the interface separation vector at the interface, T s is the stress vector defined over the discontinuity 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: Using equations ( 35) and ( 36) one can obtain the following expression for the free surface energy: This leads to the expression (originally derived by Cazes et al. [31]) for an increment of dissipated energy, given as: The total interface separation vector ⟦Δ⟧ is divided into elastic and plastic components.
The free surface energy can trivially be defined as: Introducing equation ( 41) into equation ( 39), an expanded expression for an increment of dissipated energy is obtained, such that: 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: Where the dissipated energy due to the plastic interface separation ⟦Δ⟧  is given as: and the dissipated energy due to the elastic interface separation ⟦Δ⟧  is given as: 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 expression to estimate the thermodynamic legitimacy of the proposed model:

VP-CZM parametric investigation of mixed-mode behaviour
The VP-CZM is shown for a range of axial fibril strain rates () 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 the 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 levels of plastic deformation (fibril pull-out) and dissipation occur during low rates of interface separation.

Non-proportional loading path response:
The mixed-mode behaviour of the VP-CZM is next explored by considering nonproportional loading paths subject to a range of strain rates for each path.The non-proportional loading paths consist of an initial mode II displacement to a prescribed value (Δ     ⁄ = 0 → 10) followed by complete normal separation (Δ  → ∞).
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 cases if significant plastic deformation occurs during the initial mode II separation, the computed level 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 entails 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 separation 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 separation 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.

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 energy is also modeindependent.The stored elastic energy is also shown in Figure 12(B).Tangential and normal components of the tractionseparation relationships are shown in Figure 12(C,D).Corresponding relationships are shown in Figure 12(E-H) for the cases 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 ).The calculated fracture energy is show in the case of: (B)   =   ,   =   and (F)   =   ,   =   .The total fracture energy is shown as well as the individual energy contributions of the elastic recoverable energy (  ) and the plastic irrecoverable (dissipated) energy (  ).

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 for 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 tissue 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-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 consistent behaviour with (A-B) in mixed mode scenarios.

Simulation of artery crack propagation experiment using VP-CZM
In Figure 15 we assess the ability of the VP-CZM to predict non-uniform crack propagation rates in arteries.The model 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 16 shows the results of a finite element cohesive zone analysis of the experiments 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 followed 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 undergo no plastic deformation as the strain rate results in rapid fracture.The 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 phases mirror the trends observed in the experiment: a prolonged region of slow growth 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).

Concluding 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 initial 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 propagation 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 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 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 loading 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 dissipation (representing fibril pull-out) at low interface separation rates.
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 develops 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 instantaneous incremental dissipation and compliance with the Clausius-Duhem inequality is demonstrated for a wide range of proportional and non-proportional mixed mode loading 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 propagation [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][37][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 constructs [19,20,39,40].The framework can also be used to predict fibre bridging and mode-dependant crack kinking cracking in bone [41][42][43][44][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].

Figure 1 .
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 circumferential strain.Three characteristic regimes observed and are approximately indicated on the curves; (blue) initial crack growth is slow due to fibrillation; (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 show detailed views of interfacial fibrillation.Experimental image of the ruptured aorta specimen 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 growth predictions of the exponential CZM for two different mode II fracture energy (GII)

Figure 2 .
Figure 2. (a) illustration 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 configureations are shown./  = . is the point of crack initiation and /  = . shows the progression of the mode II crack.

Figure 3 .
Figure 3. (a) Schematic of the EF-CZM undergoing a separation.  is the fibrillation stiffness,    is the fibrillation 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 parameters are shown as a function of mode angle.Interface strength, elastic limit, separation at fibrillation initiation, maximum fibril extension, and fracture energy are all shown.

Figure 4 .
Figure 4. Response of EF-CZM to proportional loading paths with a constant mode angle .(A)Traction-separation 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   =       =    (     )

2 . 2 .
the following cases: (A-B)   =   ,  0 =  0 ; (D-E) 5  =   ,  0 =   0 ( As shown in Figure5(A), the occurrence of damage or fibrillation during the initial mode II separation depends on the prescribed value of Δ  /  .Traction-separation during the subsequent normal separation is shown in Figure5(B).If 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 fracture energy as a function of the applied value of Δ   /  is shown in Figure5(C).While the non-linear influence of fibril formation and fibril stretching is computed, the total energy is mode independent, as prescribed, and monotonic increases/decreases in the normal and tangential contributions are observed.Corresponding non-proportional loading path responses are shown in Figure5(D-F) for the case of 5  =   and  0 =   0 Mode mixity of the total fracture energy is specified using equation (9).As shown in Figure5(D,E), the occurrence of fibrillation during the initial mode II separation results in reduced elastic fibril deformation in the subsequent normal separation, with early onset of fibril damage and rupture.The total energy of separation is dominated by the mode II contribution during the initial phase of the loading path, as shown in Figure5(F).

Figure 5 :
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 results 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   =      =    (     )

Figure 7 .
Figure 7. Exploration of loading-unloading paths in the normal direction (A) and tangential direction (C).Positive instantaneous incremental dissipation for the normal loading path (B) and for the tangential loading path (D).

Figure 8 .
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 fibril extension (   ) on predicted crack lengths.

Figure 10 .
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.

Figure 12 .
Figure 12.Traction separation response for the proportional loading paths subject to a constant separation rate (̇    ⁄ = .) for

Figure 14 .
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

Figure 15 (
I) shows 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 (Figure15(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 (Figure15(III), resulting in a slowing of the rate of crack growth.

Figure 16 .
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.