A new anisotropic soft tissue model for elimination of unphysical auxetic behaviour

Auxetic behaviour , the unphysical transverse expansion during uniaxial tension, is a common and undesirable feature of classical anisotropic hyperelastic constitutive models for soft tissue. In this study we uncover the underlying mechanism of such behaviour; high levels of in-plane compaction occurs due to increasing tension in strain-stiffening ﬁbres, leading to unphysical out-of-plane expansion. We demonstrate that auxetic behaviour is primarily inﬂuenced by the ratio of ﬁbre to matrix stiffness, and is accentuated by strain-stiffening ﬁbres in a constant stiffness matrix (e.g., the widely used neo-Hookean matrix with exponentially stiffening ﬁbres). We propose a new bilinear strain stiffening ﬁbre and matrix (BLFM) model which allows close control of the ﬁbre-matrix stiffness ratio, thereby robustly eliminating auxetic behaviour. We demonstrate that our model provides accurate prediction of experimentally observed out-of-plane compaction, in addition to stress-stretch anisotropy, for arterial tissue subjected to uniaxial tension testing.


Introduction
The hyperelastic strain energy density function used to model fibrous soft tissues is typically decomposed into an anisotropic part (representing the contribution of aligned families of fibres), and an isotropic part (representing the contribution of non-fibrous ground matrix).Common examples include models of arterial tissue (Gasser et al., 2006;Nolan and McGarry, 2016a;Fereidoonnezhad et al., 2016;Holzapfel et al., 2015), cartilage (Dowling et al., 2013;Brown et al., 2009;Pierce et al., 2013;Deneweth et al., 2015), and the annulus fibrosis (Guo et al., 2006;Peng et al., 2006).The most commonly imple-mented form for the strain energy density function entails exponentially strain-stiffening fibres in parallel with a near-linear neo-Hookean ground matrix, following from the work of Holzapfel, Gasser and Ogden (2000).While such an additive decomposition of soft tissue stress into a fibre contribution and a ground matrix contribution is attractive in terms of physical interpretation of the constitutive law, auxetic behaviour (i.e.orthogonal expansion during uniaxial tension) is a common and undesired feature of such models (Skacel andBursa, 2016, 2019;Volokh, 2017).To better illustrate such auxetic behaviour , we reproduce a simulation from the 'Abaqus benchmarks Guide' (input file adventitia_circ.inp) of uniaxial tension of an adventitial strip cut along circumferential direction (Abaqus, 2017).The material is modelled using the well-known HGO constitutive model (Gasser et al., 2006), with fibre families prescribed at ±49.98 • , as shown in Figure 1A.As shown in Figure 1 thickness in the out-of-plane 3-direction increases as the specimen is stretched.
While auxetic behaviour has been reported for certain classes of synthetic materials (e.g.foams and lattice materials (Lim, 2015)), such auxetic behaviour is not observed in experimental tests of artery tissue (Skacel and Bursa, 2016).To the best of the authors knowledge, auxetic behaviour has not been reported in any naturally occurring soft tissue or tissue engineered collagen constructs.Therefore, in the context of soft tissue material models, we refer to the computation of auxetic behaviour as unphysical.
Auxetic behaviour of soft tissue hyperlastic models has been the subject of three recent studies.Skacel and Bursa (2016) demonstrated the computation of auxetic behaviour (which they refer to as a negative Poisson's ratio) for the anisotropic hyperelastic soft tissue models of Holzapfel et al. (2000) and Gasser et al. (2006;2012).The authors incorrectly suggest that such unphysical auxetic response might be the result of model incompressibility.In a follow-on study the authors later demonstrate that introduction of volumetric compressibility into constitutive model in fact increases auxetic behaviour of the models (Skacel and Bursa (2019)).The recent study by Volokh (2017) considers an anisotropic hyperelastic model with exponential stiffening and fibre dispersion.The author suggested that a more accurate consideration of the fibre dispersion could be used to reduce the auxetic effect.Fibre dispersion is also taken into account in the models investigated by Skacel andBursa (Skacel andBursa, 2016, 2019).However, in Section 7 of the current study, we demonstrate that dispersion can merely reduce, but not eliminate, auxetic behaviour.Furthermore, fibre disperion is a physically motivated model feature that must be calibrated to reproduce observed material anisotropy.It is therefore not an independent parameter that can be used to control auxetic behaviour.
The three previous studies (Skacel andBursa, 2016, 2019;Volokh, 2017) demonstrate that incorrect computation of auxetic behaviour is a general problem for established anisotropic hyperelastic soft tissue models.However, these studies do not uncover the underlying mechanistic cause of such auxetic behaviour, but merely focus on secondary factors such as compressiblity and fibre dispersion.No study to date has uncovered the root cause of computed auxetic behaviour in anisotropic soft tissue models, and no study has proposed a robust solution for such incorrect behaviour.
The current paper (i) provides, for the first time, an explanation for the underlying mechanism of auxetic behaviour in anisotropic soft tissue models, and (ii) proposes a new constitutive law for anisotropic soft tissue that robustly controls and eliminates auxetic behaviour.The paper is structured as follows.In Section 2 we provide a mechanistic explanation for the computation of auxetic behaviour in fibre reinforced soft tissue models.In Section 3 we demonstrate the key role of the fibre-matrix stiffness ratio in auxetic behaviour.In Section 4 we demonstrate that strain stiffening fibres in a constant stiffness matrix promotes auxetic behaviour at high strains.In Section 5 we propose a novel bilinear strain stiffening fibre and matrix model (BLFM model) and we demonstrate that this new formulation allows robust control of out-of-plane deformation and elimination of auxetic behaviour, resulting in accurate prediction of experimentally measured artery anisotropic strain stiffening and out-of-plane contraction.

Why does auxetic behaviour occur?
In order to provide a clear explanation of the origin of auxetic behaviour for constitutive models of fibrous tissue, we consider the simple case of a cube that consists of a matrix material and two families of embedded collagen fibres.Following the typical approach for artery models, we assume that the fibres are confined to the (1-2) plane and are symmetric about the 1-direction (see Figure 1A).Following the established approach for the modelling of hyperelastic soft tissue, the strain energy density function is additively decomposed into an isotropic matrix contribution and an anisotropic fibre contribution: This leads to an additive decomposition of the Cauchy stress tensor: Therefore, when the cube is subjected to uniaxial tension in the 1-direction, the equilibrium equations in the 2-and 3-directions can be written as: where σ m 2 and σ m 3 are the 2-and 3-components of the matrix Cauchy stress tensor, respectively.σ f 2 is the 2-component of the fibre Cauchy stress tensor.Because we assume that fibres are confined to the (1-2) plane without any out-of-plane contribution, σ f 3 is equal to zero.Additionally, due to symmetry of the two fibre families about the loading direction, the 1-, 2-, and 3-direction are the principal directions.Typically the constitutive law for fibres assume significant strain stiffening (e.g.exponential stiffening of collagen fibres) whereas matrix constitutive laws typically assume much lower strain stiffening (e.g.neo-Hookean model without progressive stiffening).Therefore, when a high level of stretching is imposed on the cube in the 1-direction, the fibres exert significant tension in the 2-direction (in addition to a larger contribution in the 1-direction).Therefore, to achieve stress equilibrium, the near constant stiffness matrix must develop large negative strains in the 2-direction such that σ m 2 = −σ f 2 .As a result of extreme matrix compaction in the 2-direction, the matrix may undergo expansion in the 3-direction in order to satisfy stress equilibrium.Such expansion in a lateral direction during uniaxial stretching is referred to as auxetic behaviour.Such behaviour is not observed experimentally for arterial tissue (Skacel and Bursa, 2016), or, to our knowledge, any other soft tissue.Therefore, new modelling strategies to control and eliminate the computation of non-physical auxetic behaviour must be developed for anisotropic fibrous soft tissue.

The simple case of linear fibres in a linear matrix
We begin our investigation with a simple parametric analysis of the sensitivity of auxetic behaviour to (i) fibre to matrix stiffness ratio, (ii) fibres orientation and (iii) matrix compressibility.In this Section we initially limit our analysis to linear fibre and linear matrix behaviour in order to elucidate the factors underlying auxetic behaviour.Once again, a cube with two fibre families confined to the (1-2) plane is considered, with applied uniaxial extension in the 1-direction (Figure 2).Hyperelastic constitutive models for arterial tissue and other fibrous soft tissues are defined through the strain energy density function Ψ as follows: Ψ m (C) and Ψ f i (C, a 0i ) are the strain energy density contributions of the ground matrix and collagen fibres, respectively.Two fibre families (i = 4, 6) are represented in equation (4), as commonly assumed for arterial tissue.C = F T F is the right Cauchy-Green tensor, where F is the deformation gradient tensor.Typically for soft tissue (e.g.arteries) two fibre families are represented, and a 0i are two fibre directions in the reference configuration.Typically Ψ m (C) is decomposed into isochoric and volumetric matrix contributions, such that (5) C is the isochoric modified right Cauchy-Green tensor, while J = det F > 0 is the volumetric jacobian.The simple neo-Hookean formulation is generally used to model the ground matrix of soft tissues, such that: where K is the effective bulk modulus and µ is the effective shear modulus.Ī1 = trC is the first invariant of C.This formulation provides near-linear stress-strain behaviour.It is convenient to define E m , the effective matrix stiffness under conditions of uniaxial stress, such that E m = 9Kµ/(3K + µ).Note that the volumetric contribution in equation ( 6) is appropriate for soft tissues that do not undergo large deformations (e.g.arteries (Nolan and McGarry, 2016b)).However, for materials that exhibit large volume changes the reader is referred to Moerman et al. (2020) for more complex volumetric strain energy density functions.For the initial investigation in this section we consider a simplified constant stiffness fibre strain energy density function is given as which leads to a linear relationship between fibre Cauchy-stress and fibre stretch, where E 1f is the fibre stiffness and λ f i is the fibre stretch obtained as: For θ = 30 • auxetic behaviour is computed from the onset of stretching.This occurs due to a high component of fibre tension in the 2-direction.Stress equilibrium therefore results in high matrix compaction in the 2-direction, and consequently an auxetic expansion in the 3-direction.λ 3 is observed to peak at an applied stretch of λ 1 ≈ 2. Further stretching results in further rotation of the fibres towards the loading direction and a consequent slight reduction in λ 3 .A similar trend is observed for a fibre angle of θ = 45 • , with λ 3 peaking at a higher stretch of λ 1 ≈ 2.8.
For fibre angles greater than θ = 45 • , the fibres are in compression at the start of the deformation.Since fibres are not capable of bearing compressive loads (McEvoy et al., 2018;Holzapfel and Ogden, 2017;Holzapfel and Fereidoonnezhad, 2017), during the initial stages of stretching the material behaves like an isotropic matrix with no fibre contribution and consequently no auxetic behaviour .However, as the deformation increases the fibres rotate towards the stretching (1-) direction and develop a tensile 22-Cauchy stress component.Eventually auxetic behaviour begins to appear and progressively gets worse.
The effect of material compressibility on auxetic behaviour is also investigated ( Figure 3 C, F).It is observed that auxetic behaviour increases with increasing material compressibility.Figure 3 presents a limited parametric investigation to illustrate the mechanisms underlying auxetic behaviour.It should be noted that the full range of interactions between the governing parameters is not presented.For example, auxetic behaviour will be computed for θ = 10 • if a fibre to matrix stiffness ratio of 4. Auxetic behaviour is increased by strain-stiffening fibres in a linear matrix The isotropic neo-Hookean model with a linear fibre model has been used for the illustration of the auxetic behaviour and parameter studies in Section 2. However, experimental data of fibrous soft tissue exhibits anisotropic strain-stiffening behaviour .Typically this is captured using a strain-stiffening constitutive model for fibres in parallel with a neo-Hookean matrix.
In order to systematically analyse the role of fibre strain stiffening on auxetic be- The proposed bilinear stress-strain relation for the collagen fibres is given as: where, E as straight stiff elastic fibres.Moreover, p f , q f and r f are not independent parameters; in order to maintain C 0 and C 1 continuity the following relations must be enforced: The corresponding fibre strain energy density function is given as: where, ψ 01 and ψ 02 are two constants which ensure the continuity of strain energy.The BLF model is placed in parallel with a linear neo-Hookean matrix and subjected to uniaxial tension boundary conditions.The influence of the fibre strain stiffening in regimes (ii) and (iii) is investigated in Figure 5.In the low-stiffness regime no auxetic behaviour is observed.The onset of fibre strain stiffening in the transition regime leads to an increased rate of compaction in the 2-direction and the onset of auxetic behaviour (λ 3 > 1).The degree of auxetic deformation increases with increasing fibre strain-stiffening.However, even at low ratios of fibre stiffness to matrix stiffness (e.g. ) auxetic deformation occurs upon the onset of fibre strain stiffening.In order to accurately model the anisotropic strain-stiffening behaviour observed in experimental testing of arteries, a high fibre stiffness is required at high strains, so that auxetic behaviour becomes unavoidable if a linear neo-Hookean (constant stiffness) matrix is used.Moreover, increasing the transition strain D 1f leads to postponing the auxetic behaviour as it postpones the fibre strain-stiffening (Figure 5 D, E, F).However, it is important to note that the parameters explored in Figure 5 (i.e. the fibre-matrix stiffness ratio (E 2f /E m ) and the values of transition strain (e.g.D 1f or D 2f ) are calibrated to fit experimental stress-strain curves, and cannot be arbitrarily chosen to eliminate auxetic behaviour.

A strain-stiffening matrix with a strain-stiffening fibre reduces auxetic behaviour
As demonstrated in Figures 3 and 5, auxetic deformation is strongly influenced by the ratio of fibre to matrix stiffness.Auxetic deformation is accentuated by strain-stiffening fibres in a constant stiffness matrix.Therefore, in this section we propose a strain-stiffening isotropic matrix model with near bilinear behaviour (BLM model), placed in parallel with the BLF model described in Section 4. This approach allows precise control of the fibre-matrix stiffness ratio at small and large deformations, facilitating the simulation of anisotropic strain-stiffening behaviour without unphysical auxetic behaviour.
Assuming the isotropic strain energy function Ψ m ( λ1 , λ2 , λ3 ) for matrix, where λi , i = 1, 2, 3 are the modified principal stretches, the Cauchy stress components, σ i are determined from the following relationship (Holzapfel, 2000): Using the Valanis-Landel hypothesise the strain energy function Ψ m ( λ1 , λ2 , λ3 ) is presented as In order to simulate an isotropic matrix material that exhibits bi-linear type strain-stiffening we propose a isochoric strain energy density function such that: where the parameters E 1m , E 2m , D 1m and D 2m have the same meaning as the corresponding parameters in the fibre model (equation ( 8)) and p m , q m , and r m are obtained in a similar manner as equation ( 9) by using the corresponding matrix parameters.After a straightforward manipulation the strain energy function associated with equation ( 13) is obtained as: where, ψ 01m and ψ 02m are two constants which ensure the continuity of strain energy.As an alternative formulation to the piece-wise strain-stiffening constitutive law above, we propose a non-piece-wise equivalent formulation in Appendix A (based on error function strain stiffening).It is shown that identical material behaviour can be obtained, both in terms of stress and strain energy density.
This proposed BLM model is used in parallel with our BLF (equation ( 10)) to investigate the control/elimination of auxetic deformation under uniaxial tension.Henceforth, this combination of matrix and fibre models is called BLFM model for convenience.In Figure 6 we explore the influence of matrix strain-stiffening on auxetic behaviour.Matrix strain-stiffening is characterised by the ratio of the matrix tangent modulus at high strains (E 2m ) to that at low strains (E 1m ).Increasing E 2m /E 1m > 1 results in reduced compaction in the 2-direction at higher applied stretch, and consequently a significant reduction in auxetic behaviour.For the parameter sets considered in Figure 6, values of λ 3 > 1 are computed only for the case of a non-strain stiffening matrix (E 2m /E 1m = 1).
We also explore the influence of the matrix transition strain parameter D 1m on auxetic behaviour (with E 2m /E 1m = 3).If the value of D 1m is increased, the matrix strainstiffening is delayed, and may be preceded by significant fibre strain-stiffening.Auxetic behaviour occurs if D 1m is significantly higher than D 1f .However, if the fibre and the matrix undergo strain-stiffening at similar levels of deformation, auxetic behaviour may be eliminated.
To demonstrate the practical importance of our findings, we use our new fibre/matrix strain-stiffening formulation to simulate experimental test data for aorta tissue (FitzGibbon and McGarry, 2020).The results, as shown in Figure 7, reveal that the BLFM model is not only able to acccurately capture experimentally measured anisotropic stress-stretch data, but also correctly predicts the non-auxetic (λ 3 < 1) out-of-plane stretch with reasonable accuracy.In contrast, if a non-strain stiffening linear matrix (L-BLF) model is used, significant auxetic behaviour (λ 3 > 1) is incorrectly computed.

Auxetic behaviour using an exponentially stiffening fibre model
The strain energy density function for the classical exponential form for collagen fibre strain-stiffening (Holzapfel et al. (2000)), is given as: where k 1 and k 2 are material parameters.Note that equation ( 15) is the modified anisotropic (MA) form for slightly compressible materials ( (Nolan et al., 2014)) requiring the use of the full anisotropic invariants I 4,6 instead of the isochoric anisotropic invariants Ī4,6 .
This exponential fibre component is placed in parallel with a neo-Hookean matrix (equation ( 6)).This exponential-fibre formulation combined with neo-Hookean matrix has been widely used following implementation in leading commercial FE software programs such as Abaqus (Abaqus, 2017), FEAP (Taylor, 2013), Ansys (Ansys, 2017) and ADINA (Adina, 2017).It worth to note that although the contribution of the fibres in compression is usually excluded, ANSYS implementation does not include the tension-compression switch of the fibres.
The response to a unixial tension test of this exponential-fibre formulation combined with neo-Hookean matrix is demonstrated in Figure 8 using material parameters calibrated for human coronary arteries (Holzapfel, 2002).Strong strain stiffening and anisotropy is evident in Figure 8A.However, significant auxetic behaviour is evident (Figure 8B), with λ 3 > 1 when λ 1 > 1.15.Based on our proposed concept of using a non-linear matrix to control auxetic behaviour (Section 5), here we next explore the combination of the exponential fibre model (equation ( 15)) with the established Ogden strain-stiffening isotripic hyperelastic model for the matrix.The Ogden hyperelastic strain energy function is given as: where, λ1 , λ2 and λ3 are isochoric principal stretches and µ i and α i are material constants and we assumed N = 2 in this paper.Additionally, we retain the isotropic volumetric strain energy density function (first term on RHS of equation ( 6)) and enforce near incompressibility with K = 1.4 MPa.The material parameters of the 'Ogden-Exponential' model was chosen to reproduce the anisotropic strain-stiffening behaviour of the aforementioned 'neo-Hookean-Exponential' model (see Figure 8A).Two parameter sets have been found for 'Ogden-Exponential' model which fit the corresponding stress-stain curve of the 'neo-Hookean-Exponential' curve, as shown in Figure 8A.In parameter set 1 a low rate of matrix strain-stiffening is considered.In parameter set 2 a higher rate of matrix strain-stiffening is considered.As shown in Figure 8B, using the Ogden model in place of a neo-Hookean model for the matrix does not eliminate auxetic behaviour, but merely postpones it (in this specific case λ 3 > 1 when λ 1 > 1.4.)The stiffness ratio between exponentially stiffening fibres and the power-law stiffening matrix will changes rapidly with increasing applied strain, resulting in uncontrolled auxetic behaviour.In contrast, our proposed BLFM model facilitates control of stiffness ratios and, consequently, elimination of auxetic behaviour.

Does including fibre dispersion reduce auxetic behaviour?
In this section the influence of fibre dispersion on auxetic behaviour is investigated.
To this end, we have implemented a discrete fibre dispersion method around two mean fibre directions for the exponential (MA) model and the BLF model via the user-defined subroutine UMAT in Abaqus where the mechanical contribution of all fibres under compression are excluded.We consider that fibres can exist in m discrete directions in a 3D sphere at each integration point of a finite element model.A non-uniform segmentation of the unit sphere has been done such that higher density of segmentation around the mean fibre direction is obtained.A sensitivity analysis of the number of fibre directions m revealed that convergence is achieved for m = 50.
For each fibre family, the transversely isotropic and π-periodic von-Mises probability distribution function (PDF) ρ i , as proposed by Gasser et al., (2006), is used to describe fibre dispersion, whereby where a 0,i is a unit vector indicating one of the m discrete directions, a 0,mean is a unit vector indicating the mean fibre direction, b is a constant dispersion parameter, and erfi(x) denotes the imaginary error function.The introduction of fibre dispersion results in a fibre contribution in the 3-direction.Therefore, in equation ( 3), the 33-component of the fibre stress tensor σ f 3 will not be zero, and may therefore influence the computation of auxetic behaviour.
The results for the BLF model and exponential (MA) fibre model, both with a neo-Hookean matrix, are presented in Figure 9. Figure 9 suggests that the inclusion of fibre dispersion does decrease the auxetic deformation in both BLF and exponential MA models.As fibre dispersion is increased the level of auxetic behaviour decreases.However, even for highly dispersed fibres (b = 10) auxetic behaviour occurs in both the BLF and exponential MA model.Only in the case of a nearly isotropic fibre distribution (i.e. an isotropic material with nearly equal fibre contribution in all directions) auxetic behaviour is eliminated.

Summary and conclusions
In this study we uncover the underlying mechanism of auxetic behaviour in anisotropic hyperelastic soft tissue models that comprise of strain-stiffening fibres in parallel with an isotropic matrix.High levels of in-plane matrix compaction due to increasing tension in strain-stiffening fibres leads to unphysical out-of-plane expansion, i.e. auxetic behaviour.Previous papers have reported the secondary influence of compressibility (Skacel andBursa, 2016, 2019) and fibre dispersion (Volokh, 2017) on the degree of auxetic behaviour, but no previous study has uncovered this primary underlying mechanism.Furthermore, we propose a new bilinear strain-stiffening fibre and matrix (BLFM) model which facilitates close control of the fibre-matrix stiffness ratio, thereby robustly eliminating auxetic behaviour.We demonstrate that our model provides accurate prediction of experimentally observed out-of-plane compaction, in addition to stress-stretch anisotropy, for arterial tissue subjected to uniaxial tension testing.The key findings of our study are summarised as follows: • Auxetic deformation is strongly influenced by the ratio of fibre to matrix stiffness.
• Auxetic deformation is accentuated by strain-stiffening fibres in a constant stiffness matrix (e.g., the widely used neo-Hookean matrix with exponentially stiffening fibres).
• Our proposed bilinear strain-stiffening fibre-matrix constitutive model (BLFM model) facilitates robust control of the fibre to matrix stiffness ratio in the low and high deformation regime, consequently eliminating the computation of unphysical auxetic behaviour.
• Our BLFM model is shown to accurately capture experimentally reported anisotropic stress-stretch data for aortic tissue, in addition to correctly predicting the degree of (non-auxetic) compaction in the out-of-plane direction.In contrast, if a non-strain stiffening linear matrix model is used, significant auxetic behaviour is incorrectly computed.
Several studies demonstrate that removal of active cellular contractility, using cytochalasin-D, exposes non-linear strain-stiffening passive behaviour, e.g.(Reynolds et al., 2014;Reynolds and McGarry, 2015).Therefore, the concept of a strain stiffening matrix introduced in our model may be motivated by the presence of the non-linear passive contribution of endothelial cells, smooth muscle cells, and fibroblasts in arterial tissue, or the passive contribution of cardiac myocytes in cardiac tissue.Of course, the complete description of the non-linear anisotropic behaviour of such tissues requires the addition of our anisotropic fibre component.
The new modelling framework will provide improved predictions of artery deformation during physiological loading and stent deployment without incorrectly predicting auxetic behaviour under such large-strain deformation.Therefore, the new BLFM model can potentially improve current practices in the field of design and regulatory approval of medical devices.Future extensions of the model should consider viscoelasticity (Jacopo and Paola, 2020), viscoplasticity (Weickenmeier and Jabareen, 2014), damage (Fereidoonnezhad et al., 2016;Holzapfel and Fereidoonnezhad, 2017), and non-affine matrixfibre interaction.The modelling framework should also be considered for anisotropic contractility of smooth muscle cells and fibroblasts (McEvoy et al., 2019), in addition to passive non-linear behaviour of the cell cytoplasm (Reynolds and McGarry, 2015;Weafer et al., 2013).While we specifically calibrate and validate our novel modelling approach using experimental data for aortic tissue, we also explore the ability of our model to robustly eliminate computation unphysical auxetic behaviour over a wide range of applied deformations (including very high deformations that may be encountered for skin (Lanir and Fung, 1974) or urinary bladder samples (Farhat et al., 2008).Future work should precisely calibrate the model using experimental data for several types of soft tissues, including myocardium (McEvoy et al., 2018), cartilage (Dowling et al., 2013) and skin (Annaidh et al., 2012b,a;Limbert, 2019).

Figure 1 :
Figure 1: Results of the Abaqus Benchmark Guide example Abaqus (2017) for uniaxial tension of an adventitial strip cut along the circumferential direction.(A) Two fibre families oriented as θ = ±49.98• are considered.(B) Increase in the thickness of the strip in the case of fully aligned fibres (κ = 0) at ∆L L 0 = 0.33.(C) An increase in the thickness of the strip is computed for ∆L L 0 > 0.2 in the case of fully aligned fibres (κ = 0) and for 0.09 < ∆L L 0 < 0.3 in the case of dispersed fibres (κ = 0.226) .The reader is refer to Abaqus (2017) for details of material parameters and dimensions.

Figure 2 :
Figure 2: Schematic of the uniaxial tension of a fibre-reinforced cube with two fibre families.(a) undeformed cube (b) deformed cube without auxetic behaviour , and (c) deformed cube with auxetic behaviour .The fibres are distributed in (1-2) plane.

Figure 3 :
Figure 3: Influence of the fibre to matrix stiffness ratio (A, D), fibre angle (B, E) and compressibility (C, F) on the in-plane (first row) and out-of-plane stretch λ 3 (Second row) for the uniaxial extension of a cubic element of fibrous materials.A neo-Hookean matrix with a linear fibre model with the reference parameters value of θ = ±45 • , E m = 0.5, E f /E m = 10 and K/µ = 100have been used in all plots (blue curves).λ 1 is the stretch in the loading direction.
haviour we propose a new bilinear strain-stiffening fibre (BLF) constitutive law.Fibre deformation is considered to consist of three distinct regimes: (i) an initial compliant quasi-linear region to model the stretching of coiled collagen fibres; (ii) a non-linear transition region to model the uncoiling of fibres; (iii) a stiff quasi-linear region to model fully uncoiled fibres.
1f and E 2f are slopes of the linear regimes, D 1f and D 2f are the values of the nominal fibre strain at the end of the first linear regime and at the beginning of the second linear regime, a i , i = 4, 6 are the directions of the fibres in deformed configuration and σ stands for the Cauchy stress.It is noted that the material parameters of the BLF model in equation (8) have a clear physical meaning that can be related to a soft tissue experimental stress-stretch curve.During the initial regime the collagen fibres are tightly coiled.As the initial regime reaches D 1f and enters the transition regime the collagen fibres uncoil which leads to the rapid stiffening of the stress-stretch curve.At D 2f , the point where the second linear regime is entered, the collagen fibres are completely uncoiled and behave

Figure 4
Figure 4 contains a schematic illustration of the material parameters of the proposed BLF model and functionality of stress and strain energy.An alternative form of the BLF model that uses a single function, rather than the piece-wise function, is presented in Appendix A.

Figure 4 :Figure 5 :
Figure 4: A schematic illustration of a typical σ f -λ f (A) and Ψ f -λ f (B) curves for the proposed BLF model.The nature of the parameters E 1f , D 1f , E 2f and D 2f , and the uncoiling of the collagen fibres during each deformation regime are illustrated.D 1f is the nominal strain at the end of the first linear regime which is controlled by the initial slope E 1f .D 2f is the nominal strain at the beginning of the second linear regime which is controlled by the slope E 2f .The lowstiffness and the high-stiffness linear regime are connected by a transition regime between D 1f and D 2f .

Figure 6 :
Figure 6: The influence of matrix stiffening and transition strain on the stress (A, D), inplane transverse stretch (B, E), and auxetic behaviour (C, F).Reference parameters D 1f =

Figure 7 :
Figure 7: The capability of the BLFM model to reproduce the stress-stretch experimental data of the uniaxial tension of ascending aorta FitzGibbon and McGarry (2020) (A) and out-of-plane stretch in circumferential test (B).FE results of the uniaxial tension tests for BLFM model and L-BLF models are also shown (C)."L-BLF" stands for a linear matrix model in parallel with the BLF model.

Figure 9 :
Figure 9: Effects of fibre dispersion on the stress-stretch curve and auxetic behaviour, for bilinear (B,C) and exponential (D,E) fibre models with neo-Hookean model for the matrix.The probability density function, ρ, for the fibres dispersion around the mean fibre directions (θ = ±45 • ) is also shown in (A).