A robust anisotropic hyperelastic formulation for the modelling of soft tissue
David R. Nolan, Artur L. Gower, Michel Destrade, Ray W. Ogden, Pat McGarry
AA robust anisotropic hyperelastic formulationfor the modelling of soft tissue
David R. Nolan a , Artur L. Gower b ,Michel Destrade b , Ray W. Ogden c , Pat McGarry aa Biomedical Engineering,National University of Ireland, Galway, Galway, Ireland b School of Mathematics, Statistics and Applied Mathematics,National University of Ireland, Galway, Galway, IrelandSchool of Mathematics and Statistics,University of Glasgow, Glasgow, Scotland
Abstract
The Holzapfel–Gasser–Ogden (HGO) model for anisotropic hyper-elastic behaviour of collagen fibre reinforced materials was initiallydeveloped to describe the elastic properties of arterial tissue, but isnow used extensively for modelling a variety of soft biological tissues.Such materials can be regarded as incompressible, and when the in-compressibility condition is adopted the strain energy Ψ of the HGOmodel is a function of one isotropic and two anisotropic deformationinvariants. A compressible form (HGO-C model) is widely used in fi-nite element simulations whereby the isotropic part of Ψ is decoupledinto volumetric and isochoric parts and the anisotropic part of Ψ isexpressed in terms of isochoric invariants. Here, by using three simpledeformations (pure dilatation, pure shear and uniaxial stretch), wedemonstrate that the compressible HGO-C formulation does not cor-rectly model compressible anisotropic material behaviour, because theanisotropic component of the model is insensitive to volumetric defor-mation due to the use of isochoric anisotropic invariants. In order to a r X i v : . [ c ond - m a t . s o f t ] S e p orrectly model compressible anisotropic behaviour we present a mod-ified anisotropic (MA) model, whereby the full anisotropic invariantsare used, so that a volumetric anisotropic contribution is represented.The MA model correctly predicts an anisotropic response to hydro-static tensile loading, whereby a sphere deforms into an ellipsoid. Italso computes the correct anisotropic stress state for pure shear anduniaxial deformation. To look at more practical applications, we de-veloped a finite element user-defined material subroutine for the sim-ulation of stent deployment in a slightly compressible artery. Signifi-cantly higher stress triaxiality and arterial compliance are computedwhen the full anisotropic invariants are used (MA model) instead ofthe isochoric form (HGO-C model). Keywords
Anisotropic, Hyperelastic, Incompressibility, Finite element,Artery, Stent
Nomenclature I – identity tensorΨ – Helmholtz free-energy (strain-energy) functionΨ vol – volumetric contribution to the free energyΨ aniso – anisotropic contribution to the free energyΨ iso – isotropic contribution to the isochoric free energyΨ aniso – anisotropic contribution to the isochoric free energy σ – Cauchy stress σ (cid:48) – deviatoric Cauchy stress q – von Mises equivalent stress σ hyd – hydrostatic (pressure) stress F – deformation gradient J – determinant of the deformation gradient; local ratio of volume change C – right Cauchy–Green tensor I – first invariant of C I , – anisotropic invariants describing the deformation of reinforcing fibres F – isochoric portion of the deformation gradient C – isochoric portion of the right Cauchy–Green deformation tensor I – first invariant of C I , – isochoric anisotropic invariants a i , i = 4 , a i , i = 4 , Fa i ) κ – isotropic bulk modulus µ – isotropic shear modulus k i , i = 1 , ν – isotropic Poisson’s ratioBold uppercase symbols represent second order tensors, bold lowercase sym-bols represent vectors and un-bold symbols represent scalars. The anisotropic hyperelastic constitutive model proposed by Holzapfel et al.(2000) (henceforth referred to as the HGO model) is used extensively tomodel collagen fibre-reinforced biological materials, even more so now that ithas been implemented in several commercial and open-source Finite Element(FE) codes for the simulation of soft tissue elasticity.The constitutive equation builds upon previously published transverselyisotropic constitutive models (e.g. Weiss et al. (1996)) and reflects the struc-tural components of a typical biological soft tissue, and hence its strain-energy density consists of two mechanically equivalent terms accounting forthe anisotropic contributions of the reinforcing fibre families, in additionto a term representing the isotropic contribution of the ground matrix inwhich the fibres are embedded. Also, it assumes that the collagen fibres donot support compression, and hence they provide a mechanical contributiononly when in tension (this may be taken care of by pre-multiplying eachanisotropic term with a Heaviside, or “switching”, function).For the original incompressible HGO model the strain energy Ψ is ex-pressed as a function of one isochoric isotropic deformation invariant (denoted I ) and two isochoric anisotropic invariants (denoted I and I ). A Lagrangemultiplier is used to enforce incompressibility (Holzapfel et al., 2000). Onceagain it should be stressed that the original HGO model is intended only forthe simulation of incompressible materials.A modification of the original HGO model commonly implemented in fi-nite element codes entails the replacement of the Lagrange multiplier penaltyterm with an isotropic hydrostatic stress term that depends on a user speci-fied bulk modulus. This modification allows for the relaxation of the incom-3ressibility condition and we therefore refer to this modified formulation asthe HGO-C (compressible) model for the remainder of this study.The HGO-C model has been widely used for the finite element simulationof many anisotropic soft tissues. For example, varying degrees of compress-ibility have been reported for cartilage in the literature (e.g. Guilak et al.(1995); Smith (2001)). It has been modelled as a compressible material usingthe HGO-C model (e.g. Pe˜na et al. (2007) used a Poisson’s ratio, ν = 0 . ν = 0 . ν = 0 . ν = 0 . − .
43 and Iannaccone et al. (2014) ν = 0 . ν = 0 . ν = 0 . • The first objective of this study is to demonstrate that the HGO-Cformulation does not correctly model an anisotropic compressible hy-perelastic material.Recently, Vergori et al. (2013) showed that under hydrostatic tension, asphere consisting of a slightly compressible HGO-C material expands into alarger sphere instead of deforming into an ellipsoid. It was suggested that thiseffectively isotropic response is due to the isochoric anisotropic invariants I i being used in the switching function instead of the full invariants I i , i = 4 , • The second objective of the study is to implement a modification of theHGO-C model so that correct anisotropic behaviour of compressiblematerials is achieved.This modified anisotropic (MA) model uses the full form of the anisotropicinvariants and through a range of case studies we show this leads to thecorrect computation of stress in contrast to the widely used HGO-C model.The paper is structured as follows. In Section 2 we demonstrate and high-light the underlying cause of the insensitivity of the anisotropic componentof the HGO-C model to volumetric deformation in compressible materials.We demonstrate that the modification of the model to include the full formof the anisotropic invariants corrects this deficiency. In Section 3 we showhow the HGO-C model yields unexpected and unphysical results for pure in-plane shear and likewise in 4 for simple uniaxial stretching, in contrast to themodified model. We devote Section 5 to two Finite Element biomechanicscase studies, namely pressure expansion of an artery and stent deployment inan artery, and illustrate the significant differences in computed results for theHGO-C model and the modified model. Finally, we provide some concludingremarks and discussion points in Section 6.
The original HGO model is intended for incompressible materials. Howevera variation of the HGO model whereby a bulk modulus is used instead of apenalty term has been implemented in a number of FE codes. Several authorshave used this formulation t model compressible anisotropic materials buusing a relatively low value of bulk modulus. An important objective of thispaper is to highlight that this HGO-C formulation does not correctly modelcompressible anisotropic material behaviour.The kinematics of deformation are described locally in terms of the defor-mation gradient tensor, denoted F , relative to some reference configuration.5he right Cauchy–Green tensor is defined by C = F T F , where T indicatesthe transpose of a second-order tensor.Hyperelastic constitutive models used for rubber-like materials often splitthe local deformation into volume-changing (volumetric) and volume-preserving(isochoric, or deviatoric) parts. Accordingly the deformation gradient F isdecomposed multiplicatively as follows: F = (cid:16) J I (cid:17) F , (1)where J is the determinant of F . The term in the brackets represents thevolumetric portion of the deformation gradient and F is its isochoric portion,such that det( F ) = 1 at all times.Suppose that the material consists of an isotropic matrix material withinwhich are embedded two families of fibres characterized by two preferreddirections in the reference configuration defined in terms of two unit vectors a i , i = 4 ,
6. With C , J and a i are defined the invariants I = tr( C ) , I = a · ( Ca ) , I = a · ( Ca ) , (2) I = J − / I , I = J − / I , I = J − / I , (3)where I i ( i = 1 , ,
6) are the isochoric counterparts of I i . The HGO modelproposed by Holzapfel et al. (2000) for collagen reinforced soft tissues ad-ditively splits the strain energy Ψ into volumetric, isochoric isotropic andisochoric anisotropic terms,Ψ ( C , a , a ) = Ψ vol ( J ) + Ψ iso (cid:0) C (cid:1) + Ψ aniso (cid:0) C , a , a (cid:1) , (4)where Ψ iso and Ψ aniso are the isochoric isotropic and isochoric anisotropicfree-energy contributions, respectively, and C = J − / C is the isochoric rightCauchy–Green deformation tensor.In numerical implementations of the model (Abaqus, 2010; ADINA, 2005;Gasser et al., 2002), the volumetric and isochoric isotropic terms are repre-sented by the slightly compressible neo-Hookean hyperelastic free energyΨ vol ( J ) = 12 κ ( J − , Ψ iso ( C ) = 12 µ ( I − , (5)where κ and µ are the bulk and shear moduli, respectively, of the softisotropic matrix. Of course one may write (5) in terms of the full invariantsalso, using the results from (3). 6he isochoric anisotropic free-energy term is prescribed asΨ aniso (cid:0) C , a , a (cid:1) = k k (cid:88) i =4 , { exp[ k (cid:0) I i − (cid:1) ] − } , (6)where k and k are positive material constants which can be determinedfrom experiments.For a general hyperelastic material with free energy Ψ the Cauchy stressis given by σ = 1 J F ∂ Ψ ∂ F . (7)For the Cauchy stress derived from Ψ above, we have the decomposition σ = σ vol + σ iso + σ aniso , where σ vol = κ ( J − I , σ iso = µ J − (cid:0) B − I I (cid:1) , (8)with B = F F T , and σ aniso = 2 k J − (cid:88) i =4 , (cid:0) I i − (cid:1) exp[ k (cid:0) I i − (cid:1) ]( a i ⊗ a i − I i I ) , (9)where a i = Fa i . This slightly compressible implementation is referred to asthe HGO-C model henceforth.The original incompressible HGO model by Holzapfel et al. (2000) speci-fied that for arteries the constitutive formulation should be implemented forincompressible materials. In that limit, κ → ∞ , ( J − → p , and the volumetric stress assumes the form, σ vol = − p I . Indeed theoriginal incompressible HGO model can equally be expressed in terms of thefull invariants I and I (with J →
1) (e.g., Holzapfel et al. (2004)).However, in the case of the HGO-C implementation, if κ is not fixednumerically at a large enough value, then slight compressibility is intro-duced into the model. The key point of this paper is that the isochoricanisotropic term Ψ aniso defined in (6) does not provide a full representationof the anisotropic contributions to the stress tensor for slightly compress-ible materials. In Section 2.4 we introduce a simple modification of theanisotropic term to account for material compressibility.7 .2 Pure dilatational deformation First we consider the case of the HGO-C material subjected to a pure dilata-tion with stretch λ = J / , so that F = λ I , C = λ I , J = λ . (10)We expect that an anisotropic material requires an anisotropic stress state tomaintain the pure dilatation. However, calculation of the invariants I i and I i yields I i = a i · ( Ca i ) = λ , I i = J − / I i = 1 , i = 4 , , (11)so that while I i is indeed the square of the fibre stretch and changes withthe magnitude of the dilatation, its isochoric counterpart I i is always unity.Referring to (9), it is clear that the entire anisotropic contribution to thestress (7) disappears (i.e. σ aniso ≡ ), and the remaining active terms arethe isotropic ones. Thus, under pure dilatation, the HGO-C model computesan entirely isotropic state of stress. Now we investigate the reverse question: what is the response of the HGO-Cmaterial to a hydrostatic stress , σ = σ I , (12)where σ > σ < C , the squared principal stretches, λ , λ , λ say, to be distinct. Hence, if the material is slightly compressible,then a sphere should deform into an ellipsoid (Vergori et al., 2013) and acube should deform into a hexahedron with non-parallel faces (N´ı Annaidhet al., 2013b).However, in the HGO-C model the Ψ aniso contribution is switched on onlywhen I i (not I i ) is greater than unity. Vergori et al. (2013) showed that infact I i is always less than or equal to one in compression and in expansionunder hydrostatic stress, so that the HGO-C response is isotropic, contrary tophysical expectations. Then we may ask if removal of the switching functioncircumvents this problem so that anisotropic response is obtained.8ith the fibres taken to be mechanically equivalent and aligned with a = (cos Θ , sin Θ ,
0) and a = (cos Θ , − sin Θ ,
0) in the reference config-uration, we have, by symmetry, I = I and I = I and Ψ = Ψ , wherethe subscripts 4 and 6 on Ψ signify partial differentiation with respect to I and I , respectively. Similarly, in the following the subscript 1 indicatesdifferentiation with respect to I . For this special case, Vergori et al. (2013)showed that the stretches arising from the application of a hydrostatic stressare λ = J / (cid:20) Ψ (Ψ + 2Ψ sin Θ)(Ψ + 2Ψ cos Θ) (cid:21) ,λ = J / (cid:20) Ψ (Ψ + 2Ψ cos Θ)(Ψ + 2Ψ sin Θ) (cid:21) ,λ = J / (cid:34) Ψ + 2Ψ Ψ + Ψ sin (cid:35) . (13)Explicitly,Ψ = ∂ Ψ ∂I = 12 µ , Ψ = ∂ Ψ ∂I = k ( I −
1) exp[ k (cid:0) I − (cid:1) ] . (14)Looking at (13), we see that there is a solution to the hydrostatic stressproblem where the stretches are unequal, so that a sphere deforms into anellipsoid. However, there is also another solution: that for which I ≡
1, inwhich case, Ψ ≡ λ = λ = λ = J / by(13). Thus, a sphere then deforms into another sphere.Of those (at least) two possible paths, FE solvers converge upon theisotropic solution. One possible explanation for this may be that the ini-tial computational steps calculate strains in the small-strain regime. In thatregime, Vergori et al. (2013) showed that all materials with a decoupled vol-umetric/isochoric free-energy behave in an isotropic manner when subjectto a hydrostatic stress. Hence the first computational step brings the de-formation on the isotropic path, and I = 1 then, and subsequently. InSection 2.2 and Section 2.3 we have thus demonstrated that the use of anisochoric form of the anisotropic strain energy Ψ aniso from the HGO modelin the HGO-C model cannot yield a correct response to pure dilatation orapplied hydrostatic stress. 9 .4 Modified Anisotropic Model for Compressible Ma-terials In order to achieve correct anisotropic behaviour for compressible materi-als we introduce a modification to the anisotropic term of the HGO model,whereby the anisotropic strain energy is a function of the ‘total’ right Cauchy–Green deformation tensor C , rather than its isochoric part C , so thatΨ ( J, C , a , a ) = Ψ vol ( J ) + Ψ iso ( J, C ) + Ψ aniso ( C , a , a ) , (15)where the expressions for strain energy density terms Ψ vol and Ψ iso arethe same as those in (5), andΨ aniso ( C , a , a ) = k k (cid:88) i =4 , { exp[ k ( I i − ] − } . (16)This modification to the HGO-C model is referred to as the modified anisotropic (MA) model hereafter. Combining (5), (15) and (16), the Cauchy stressfor the MA model is determined using (7) and the decomposition σ = σ vol + σ iso + σ aniso resulting in the expression: σ = κ ( J − I + µ J − / (cid:0) B − I I (cid:1) + 2 k (cid:88) i =4 , ( I i −
1) exp[ k ( I i − ] a i ⊗ a i . (17)where a i = Fa i , i = 4 ,
6. Now it is easy to check that in the cases ofa pure dilatation and of a hydrostatic stress, the MA model behaves in ananisotropic manner, because the term I i − (cid:54) = 0 and hence Ψ aniso (cid:54) = 0 and σ aniso (cid:54) = . This resolves the issues identified above for the HGO-C model.We have developed a user-defined material model (UMAT) Fortran sub-routine to implement the MA formulation for the Abaqus/Standard FE soft-ware. The FE implicit solver requires that both the Cauchy stress and theconsistent tangent matrix (material Jacobian) are returned by the subrou-tine. Appendix A gives the details of the consistent tangent matrix.We have used the above subroutine to repeat the simulations of expansionof a sphere under hydrostatic tension of Vergori et al. (2013), this time usingthe MA formulation. Again two families of fibres are assumed, lying in the(1 ,
2) plane and symmetric about the -axis (the sphere and axes are shown10igure 1: A) Schematic of an undeformed sphere highlighting three radiion orthogonal axes, 1-2-3, centred at the sphere origin. Two families of fi-bres are contained in the (1 ,
2) plane and symmetric about the 1-axis. B) Computed (deformed/undeformed) ratios ( r/r ) of the orthogonal radii forboth MA and HGO-C models versus the ratio σ hyd /σ maxhyd . Note that the de-formation computed for the HGO-C model incorrectly remains spherical. C) Deformed ellipsoidal shape computed for the MA model; contours illustratethe inhomogeneous distribution of stress triaxiality ( σ hyd /q ) throughout thedeformed body.in Figure 1A). The displacements of points on the surface of the sphere at theends of three mutually orthogonal radii with increasing applied hydrostatictension are shown in Figure 1B. Clearly the sphere deforms into an ellipsoidwith a major axis oriented in the 3-direction and a minor axis oriented inthe 1-direction, confirming the simulation of orthotropic material behaviour.The distribution of stress triaxiality in the deformed ellipsoid, measured by σ hyd /q , is shown in Figure 1C, where σ hyd ≡ tr( σ ) / q ≡ (cid:112) / σ (cid:48) : σ (cid:48) is the von Mises equivalent stress, σ (cid:48) being thedeviatoric Cauchy stress tensor. Clearly an inhomogeneous stress state iscomputed in the deformed body.The results shown in Figure 1 contrast sharply with the equivalent simu-lations using the HGO-C model (Vergori et al., 2013) superimposed in Figure1B for comparison. In that case a similar fibre-reinforced sphere is shown todeform into a larger sphere with a homogeneous stress distribution, indicativeof isotropic material behaviour. 11 Analysis of Pure Shear
A pure dilatation and a hydrostatic stress each represent a highly idealizedsituation, unlikely to occur by themselves in soft tissue in vivo . This sec-tion highlights the unphysical behaviour can also emerge for common modesof deformation if the anisotropic terms are based exclusively on the iso-choric invariants. Considering once again the general case of a compressibleanisotropic material, we analyse the response of the HGO-C and MA modelsto pure in-plane shear. Regarding the out-of-plane boundary conditions, wefirst consider the case of plane strain (Section 3.1). Even though this defor-mation is entirely isochoric the HGO-C model yields incorrect results. Wethen consider the case of plane stress (Section 3.2), and again demonstratethat the HGO-C model yields incorrect results. By contrast, we show thatthe MA model computes a correct stress state for all levels of compressibilityand specified deformations. In the following calculations we assume a shearmodulus, µ = 0 .
05 MPa and anisotropic material constants k = 1 MPa and k = 100. With restriction to the (1 ,
2) plane we now consider the plane strain deforma-tion known as pure shear , maintained by the application of a suitable Cauchystress. In particular, we take the deformation gradient for this deformationto have components F = (cid:112) F + 1 F F (cid:112) F + 1 00 0 1 , (18)where F is a measure of the strain magnitude. Figure 2A depicts the defor-mation of the (1 ,
2) square cross section of a unit cube, which deforms intoa parallelogram symmetric about a diagonal of the square. The deformationcorresponds to a stretch λ = (cid:112) F + 1 + F along the leading diagonal witha transverse stretch λ − = (cid:112) F + 1 − F . We can think of the deformationarising from displacement components applied to the vertices of the square,as indicated in Figure 2A. Two families of fibres, with reference unit vectors a and a are assumed to lie in the (1 ,
2) plane, as illustrated in Figure 2A,oriented with angles ± θ to the 1 axis. We perform some calculations for arange of fibre orientations for each of the HGO-C and MA models.12igure 2: A) Schematic illustrating the kinematics of the pure shear de-formation of the (1 ,
2) section of a unit cube. Note the rotated coordinatesystem (1 (cid:48) , (cid:48) ), orientated at 45 ◦ to the (1 ,
2) axes, used to specify the vertexdisplacement components u (cid:48) and u (cid:48) . Note also the vectors a i , i = 4 , θ . Resultsare displayed for a range of fibre orientations with θ from ± ◦ to ± ◦ withrespect to the (1 ,
2) coordinate system. B) Computed stress ratio σ /σ versus F for the HGO-C model, illustrating significant negative (compres-sive) stresses in the out-of-plane direction. C) Computed stress ratio versus F for the MA model, illustrating very small negative (compressive) stressesin the out-of-plane direction (an order of magnitude lower than for the HGO-C model).First we note that although, for this specific case, the free energies ofthe HGO-C and the MA models coincide (because J = 1 and hence I = I ), the corresponding stress tensors are very different. This is due to the“deviatoric” form of the anisotropic stress contribution that emerges for theHGO-C model, as in the final term of (9), compared with the final term of(17). It gives rise to a significant negative (compressive) out-of-plane stresscomponent σ which is comparable in magnitude to σ , as shown in Figure2B. Such a negative stress is anomalous in the sense that for large κ theresult for the incompressible limit should be recovered, but it is not. Indeed,if we start with the incompressible model we obtain σ = µ − p , which is13igure 3: Dimensionless stress components σ ij /k versus F for the case ofa single family of fibres orientated at θ = 30 ◦ . A) MA model; B) HGO-Cmodel.independent of σ . However, as (18) represents a kinematically prescribedisochoric deformation, the volumetric stress in the HGO-C model goes tozero and does not act as the required Lagrange multiplier.By contrast, the out-of-plane compressive normal stress component σ computed for the MA model is at least an order of magnitude lower than thein-plane shear stress component σ (Figure 2C), and is close to zero for mostfibre orientations. This is consistent with the incompressible case because,since p is arbitrary it may be chosen to be µ so that σ = 0. This is whatmight be expected physically, given that the fibres and the deformations areconfined to the (1 ,
2) plane.Because of the deviatoric component of the stress tensor emerging fromthe HGO-C model, the trace of the Cauchy stress is always zero when J = 1as equations (8) and (9) will confirm. By contrast, the trace of the Cauchystress is not zero for the MA model. Hence the in-plane stress components aresignificantly different from those for the HGO-C model, as shown in Figures3A and 3B, respectively, for the case of a single fibre family with θ = 30 ◦ .14 .2 Plane stress pure shear The kinematically prescribed isochoric deformation in Section 3.1 is volumeconserving and makes the Ψ vol terms equal to zero. We modify the out-of-plane boundary condition to enforce a plane stress ( σ = 0) simulation Thisallows a compressible material to deform of out-of-plane.A plane stress pure shear deformation is given as F = (cid:112) F + 1 F F (cid:112) F + 1 00 0 F , (19)where the out of plane stretch component F in general is not equal to 1,so that the deformation is not in general isochoric. If the bulk modulus κ is very large compared with the initial shear modulus µ , then it acts as aLagrange multiplier to enforce incompressibility, such that F = 1 (at leastapproximately). If the magnitude of the bulk modulus is reduced, then thematerial becomes slightly compressible and F (cid:54) = 1. Here we investigate thesensitivity of the stress computed for the HGO-C and MA models to themagnitude of the bulk modulus κ .First, we consider the almost incompressible case where the ratio of bulkto shear modulus is κ /µ = 2 × for the isotropic neo-Hookean componentof the model, equivalent to a Poisson ratio of ν = 0 . J = F = 1 . κ /µ = 50( ν = 0 . σ ij /k versus F for the case of a single family of fibres orien-tated at θ = 30 ◦ . A) Computed stresses for both the HGO-C and MA modelswith a large bulk modulus κ /µ = 2 × (equivalent to a Poisson ratio of0.49999975). B ) Computed stresses for the HGO-C model with κ /µ = 50(equivalent to a Poisson ratio of 0.490). Note that the stresses computed forthe HGO-C model are an order of magnitude lower in the slightly compress-ible small bulk modulus case than in the almost incompressible large bulkmodulus case.modulus and, consequently, incompressibility must be enforced by choosing avery large magnitude for the bulk modulus in order to avoid the computationof erroneous stress states.By contrast, the MA model computes identical stress states for κ /µ =2 × and κ /µ = 50 (Figure 4A in both cases). This response highlightsthe robustness of the MA model, which computes correct results for all levelsof material compressibility (including the incompressible limit). We now consider a confined uniaxial stretch, as illustrated in Figure 5A,where a stretch is imposed in the 2-direction ( λ = λ >
1) and no lateraldeformation is permitted to occur in the 1- and 3-directions ( λ = λ = 1).16uch a simple deformation may have biomechanical relevance as, for example,in a blood vessel undergoing large circumferential strain, but little or no axialor radial strain.Figure 5: A) Schematic of confined uniaxial stretch ( λ = λ > , λ = λ = 1), showing the fibre family reference directional vector a in the (1 , σ /σ is computed basedon a model with a single fibre family and plotted as a function of λ . Re-sults are displayed for a range of fibre orientations θ from 0 ◦ to 90 ◦ . B) Computed results for the HGO-C model, illustrating negative (compressive)lateral stresses. C) Computed results for the MA model, all lateral stressesbeing positive (tensile).We derive analytically the stress components for the HGO-C and MAmodels using the formulas of Section 2. We assume there is a single familyof parallel fibres aligned with the reference unit vector a in the (1 ,
2) planeand with orientation θ relative to the 1-axis ranging from 0 ◦ to 90 ◦ . Wetake µ = 0 .
05 MPa, κ = 1 MPa for the slightly compressible neo-Hookeanisotropic matrix, and material constants k = 1 MPa and k = 100 for thefibre parameters.The ratio of the lateral to axial Cauchy stress components, σ /σ , isplotted as a function of applied stretch λ for the HGO-C model (Figure5B) and the MA model (Figure 5C). Results for the HGO-C model exhibitnegative (compressive) stresses in the lateral direction for certain fibre ori-entations. This auxetic effect suggests that the material would expand in17he lateral direction in the absence of the lateral constraint and is contraryto expectations, particularly for fibre orientations closer to the axial direc-tion. In fact, here the computed lateral compressive force is most pronouncedwhen the fibre is aligned in the direction of stretch ( θ = 90 ◦ ), where a trans-versely isotropic response, with exclusively tensile lateral stresses, shouldbe expected. For all fibres orientated within about 45 ◦ of the direction ofstretch, the lateral stress changes from tensile to compressive as the appliedstretch increases. By contrast to the HGO-C model, the MA model yieldsexclusively tensile lateral stresses for all fibre orientations (Figure 5C). Following from the idealized, analytical deformations considered above, wenow highlight the practical significance of the errors computed by using theHGO-C model for slightly compressible tissue. We consider, in turn, twoFinite Element case studies using Abaqus (2010) to implement the HGO-Cand MA models with user-defined material subroutines (see Appendix A).
First we simulate the deformation of an artery under a lumen pressure ( LP ).A schematic of a quarter artery is shown in Figure 6A. The vessel has aninternal radius r i of 0.6 mm and an external radius r e of 0.9 mm. The lengthof the artery in the z -direction is 0.3 mm with both ends constrained in the z -direction.We model the wall as a homogeneous material with two families of fibreslying locally in the ( θ, z ) plane, where ( r, θ, z ) are cylindrical polar coordi-nates. The fibre families are symmetric with respect to the circumferentialdirection and oriented at ± ◦ measured from the circumferential direction.For the fibres, the material constants are k = 1 MPa and k = 2, and forthe neo-Hookean matrix, they are µ = 0 .
03 MPa, κ = 1 MPa, resulting in aslightly incompressible material (corresponding to a Poisson ratio of 0.485).A mesh sensitivity study confirms a converged solution for a model using atotal of 1,044 eight-noded full-integration hexahedral elements.The (dimensionless) changes in the internal and external radii ∆ r/r asfunctions of increasing dimensionless lumen pressure LP/LP max are plotted in18igure 6B. They reveal that the HGO-C model predicts a far more compliantartery than the MA model.Notable differences in the arterial wall stress state arise between the HGO-C and MA models. Figures 6C, D and E present the von Mises stress, pressurestress and triaxiality, respectively, in the arterial wall. The magnitude andgradient through the wall thickness of both the von Mises stress and pressurestress differ significantly between the HGO-C and MA models. This contrastis further highlighted by the differing distributions of triaxiality for bothmodels, confirming a fundamental difference in the multi-axial stress statecomputed for the two models.
The final case study examines the deployment of a stainless steel stent ina straight artery. Nowadays most medical device regulatory bodies insiston computational analysis of stents (FDA , 2010) as part of their approvalprocess. Here we demonstrate that the correct implementation of the con-stitutive model for a slightly compressible arterial wall is critical for thecomputational assessment of stent performance.We use a generic closed-cell stent geometry (Conway et al., 2012) with anundeformed radius of 0.575 mm. It is made of biomedical grade stainless steelalloy 316L with Young’s modulus of 200 GPa and Poisson’s ratio 0.3 in theelastic domain. We model plasticity using isotropic hardening J − plasticitywith a yield stress of 264 MPa and ultimate tensile strength of 584 MPa ata plastic log strain of 0.274 (McGarry et al., 2007). We mesh the stentgeometry with 22,104 reduced integration hexahedral elements. We modela balloon using membrane elements, with frictionless contact between themembrane elements and the internal surface of the stent. Finally, we simulatethe balloon deployment by imposing radial displacement boundary conditionson the membrane elements.For the artery, we take a single layer with two families of fibres symmet-rically disposed in the ( θ, z ) plane. The fibres are oriented at ± ◦ to thecircumferential direction and material constants and vessel dimensions arethe same as those used in Section 5.1. Here the FE mesh consists of 78,100full integration hexahedral elements; a high mesh density is required due tothe complex contact between the stent and the artery during deployment.“Radial stiffness”, the net radial force required to open a stent, is a com-monly cited measure of stent performance (FDA , 2010). Figure 7 presents19igure 6: A) Schematic illustrating the geometry, lines of symmetry andboundary conditions for modelling the inflation of an artery under a lumenpressure LP . B) Prediction of the internal ( r i ) and external ( r e ) radialstrain ∆ r/r = ( r − r ) /r in the artery under a normalized lumen pressure LP/LP max for the HGO-C and MA models. Panels C) , D) and E) arecontour plots illustrating the von Mises ( q ), pressure ( − σ hyd ) and triaxiality( σ hyd /q ) stresses, respectively, in the artery wall for the HGO-C and MAmodels. 20igure 7: Plot of the dimensionless radial force ( F − F ) /F required todeploy a stent in an artery with increasing stent radial expansion. Radialforce is normalized by the radial force at the point immediately before contactwith the artery ( F ). The radial expansion is normalized using the initialundeformed internal radius ( r i ) and the final fully deployed internal radius( r f ). Note that the HGO-C model predicts a more compliant artery than theMA model.plots of the predicted net radial force as a function of radial expansion forthe HGO-C and MA models. The predicted radial force required to expandthe stent to the final diameter is significantly lower for the HGO-C modelthan for the MA model. This result correlates with the previous finding inSection 5.1 that the HGO-C model underestimates the arterial compliance,with significant implications for design and assessment of stents.Figure 8 illustrates the notable differences that appear in the artery stressstate between the HGO-C and MA models. Again, higher values of von Misesstress (Figure 8A) and pressure stress (Figure 8B) are computed for the MAmodel. Both the triaxiality (Figure 8C) and the ratio of axial to circumferen-tial stress (the stress ratio in the plane of the fibres) (Figure 8D) confirm thatthe nature of the computed multi-axial stress state is significantly different21igure 8: Contour plots illustrating differences in the stresses computed forthe HGO-C and MA models after stent deployment. A) von Mises stress q , B) pressure stress - σ hyd , C) triaxiality, D) ratio of axial stress to thecircumferential stress σ zz /σ θθ .between the MA and HGO-C models.A detailed examination of the stress state through the thickness (radialdirection) of the artery wall is presented in Figure 9. A comparison betweenHGO-C and MA simulations in terms of the ratios of the Cauchy stress com-ponents emphasizes further the fundamentally different stresses throughoutthe entire artery wall thickness. It is not merely that the MA model calcu-lates a different magnitude of stress, rather the multi-axiality of the stressstate has been altered. 22igure 9: Stress measures computed through the arterial wall from the inter-nal ( r i ) to external radius ( r e ) at full deployment of the stent for the HGO-Cand MA models. A) Triaxiality ratio σ hyd /q of the pressure stress to vonMises stress. B) Ratio σ zz /σ θθ of the axial to circumferential stress. C) Ratio σ rr /σ zz of the radial to axial stress. D) Ratio σ rr /σ θθ of the radial tocircumferential stress. The original HGO model (Holzapfel et al. (2000)) is intended for modelling ofincompressible anisotropic materials. A compressible form (HGO-C model)is widely used whereby the anisotropic part of Ψ is expressed in terms ofisochoric invariants. Here we demonstrate that this formulation does not cor-rectly model compressible anisotropic material behaviour. The anisotropic23omponent of the model is insensitive to volumetric deformation due to theuse of isochoric anisotropic invariants. This explains the anomolous finiteelement simulations reported in Vergori et al. (2013), whereby a slightly com-pressible HGO-C sphere was observed to deform into a larger sphere undertensile hydrostatic loading instead of the ellipsoid which would be expectedfor an anisotropic material. In order to achieve correct anisotropic compress-ible hyperelastic material behaviour we present and implement a modified(MA) model whereby the anisotropic part of the strain energy density is afunction of the total form of the anisotropic invariants, so that a volumet-ric anisotropic contribution is represented. This modified model correctlypredicts that a sphere will deform into an ellipsoid under tensile hydrostaticloading.In the case of (plane strain) pure shear, a kinematically enforced iso-choric deformation, we have shown that a correct stress state is computedfor the MA model, whereas the HGO-C model yields incorrect results. Cor-rect results are obtained for the HGO-C model only when incompressibilityis effectively enforced via the use of a large bulk modulus, which acts as aLagrange multiplier in the volumetric contribution to the isotropic terms (inthis case HGO-C model is effectively the same as the original incompressibleHGO model). In the case of a nearly incompressible material (with Poisson’sratio = 0.490, for example) we have shown that the in-plane stress compo-nents computed by the HGO-C model are reduced by an order of magnitude.Bulk modulus sensitivity has been pointed out for isotropic models by Gentet al. (2007) and Destrade et al. (2012), and for the HGO-C model by N´ı An-naidh et al. (2013b). Here, we have demonstrated that a ratio of bulk to shearmodulus of κ /µ = 2 × (equivalent to a Poisson’s ratio of 0.49999975)is required to compute correct results for the HGO-C model. By contrast,the MA model is highly robust with correct results being computed for alllevels of material compressibility during kinematically prescribed isochoricdeformations.From the view-point of general finite element implementation, red therequirement of perfect incompressibility (as in the case of a HGO material)can introduce numerical problems requiring the use of selective reduced inte-gration and mixed finite elements to avoid mesh locking and hybrid elementsto avoid ill-conditioned stiffness matrices. Furthermore, due to the complexcontact conditions in the simulation of balloon angioplasty (both between theballoon and the stent, and between the stent and the artery), explicit FiniteElement solution schemes are generally required. However, Abaqus/Explicit24or example has no mechanism for imposing an incompressibility constraintand assumes by default that κ /µ = 20 ( ν = 0 . κ /µ > ν = 0 . I i is appropriate.The current paper demonstrates the importance of a volumetric anisotropiccontribution for compressible materials, highlighting the extensive range ofnon-physical behaviour that may emerge in the simulation of nearly incom-pressible materials if the HGO-C model is used instead of the MA model.Examples including the Finite Element analysis of artery inflation due toincreasing lumen pressure and stent deployment. Assuming nearly incom-pressible behaviour ( ν = 0 . Acknowledgements
DRN and ALG wish to acknowledge the receipt of their PhD scholarshipsfrom the Irish Research Council. MD and RWO wish to thank the RoyalSociety for awarding them an International Joint Project. This research was25lso funded under Science Foundation Ireland project SFI-12/IP/1723.
ReferencesReferences
Abaqus/Standard users manual, Ver. 6.10 (2010) Dassault Syst`emes SimuliaCorporation, Pawtucket.ADINA theory and modeling guide (2005). ADINA R&D, Inc., Watertown.Conway C, Sharif F, McGarry JP, McHugh P (2012) A computational test-bed to assess coronary stent implantation mechanics using a population-specific approach. Cardiovasc Eng Technol 3:374-387Cardoso L, Kelly-Arnold A, Maldonado N, Laudier D, Weinbaum S (2014)Effect of tissue properties, shape and orientation of microcalcifications onvulnerable cap stability using different hyperelastic constitutive models. JBiomech 47:870-877Destrade M, Gilchrist MD, Motherway J, Murphy JG (2012) Slight com-pressibility and sensitivity to changes in Poisson’s ratio. Int J Num MethEng 90:403-411Dowling EP, Ronan W, McGarry JP Computational investigation of in situchondrocyte deformation and actin cytoskeleton remodelling under physi-ological loading Acta Biomater 9:5943-5955Famaey N, Sommer G, Vander Sloten J, Holzapfel GA (2012) Arterial clamp-ing: finite element simulation and in vivo validation. Journal mech behavbiomed 12:107-118.FDA (2010) Guidance for Industry and FDA Staff - Non-Clinical EngineeringTests and Recommended Labeling for Intravascular Stents and AssociatedDelivery Systems.Federico S (2010) Volumetric-distortional decomposition of deformation andelasticity tensor. Math Mech Solids 15:672-69026arc´ıa A, Pe˜na E, Mart´ınez MA (2012) Influence of geometrical parameterson radial force during self-expanding stent deployment. Application fora variable radial stiffness stent. J Mech Behav Biomed Mat 10:166-175.doi:http://dx.doi.org/10.1016/j.jmbbm.2012.02.006Gasser, T.C., Holzapfel G.A. (2002) A rate-independent elastoplastic con-stitutive model for biological fiber-reinforced composites at finite strains:continuum basis, algorithmic formulation and finite element implementa-tion. Comput Mech 29:340-360. doi: 10.1007/s00466-002-0347-6Suh JB, Gent AN, Kelly SG (2007) Shear of rubber tube springs. Int J Non-Linear Mech 42:1116-1126Guilak F, Ratcliffe A, Mow VC (1995) Chondrocyte deformation and localtissue strain in articular cartilage: a confocal microscopy study. J OrthopRes, 13: 410-421.Helfenstein J, Jabareen M, Mazza E, Govindjee S (2010) On non-physicalresponse in models for fiber-reinforced hyperelastic materials. Int J SolidsStruct 47:2056-2061Holzapfel GA, Gasser TC, Ogden RW (2000) A new constitutive frameworkfor arterial wall mechanics and a comparative study of material models. JElast 61:1-48Holzapfel GA, Gasser TC, Ogden RW (2004) Comparison of a Multi-LayerStructural Model for Arterial Walls With a Fung-Type Model, and Issuesof Material Stability. J Biomech Eng 126:264-275Huang R, Becker AA, Jones IA (2012) Modelling cell wall growth usinga fibre-reinforced hyperelasticviscoplastic constitutive law. J Mech PhysSolids 60:750-783Iannaccone F, Debusschere N, De Bock S, De Beule M, Van Loo D, VermassenF, Segers P, Verhegghe B (2014) The Influence of Vascular Anatomy onCarotid Artery Stenting: A Parametric Study for Damage Assessment JBiomech 47:890-898Kiousis DE, Wulff AR, Holzapfel GA (2009) Experimental studies andnumerical analysis of the inflation and interaction of vascular ballooncatheter-stent systems. Annals Biomed Eng 37:315-33027aquer G, Laurent M, Brandejsky V, Pretterklieber ML, Zysset PK (2014)Finite Element Based Nonlinear Normalization of Human Lumbar Inter-vertebral Disc Stiffness to Account for Its Morphology. Journal of Biomedeng 136McGarry J, O’Donnell B, McHugh P, O’Cearbhaill E, McMeeking R (2007)Computational examination of the effect of material inhomogeneity on thenecking of stent struts under tensile loading. J Appl Mech 74:978-989N´ı Annaidh A, Bruy`ere K, Destrade M, Gilchrist MD, Maurini C, Ott´enio M,Saccomandi G (2013a). Automated estimation of collagen fibre dispersionin the dermis and its contribution to the anisotropic behaviour of skin,Annals Biomed Eng 40:1666-1678N´ı Annaidh A, Destrade M, Gilchrist MD, Murphy JG (2013b) Deficienciesin numerical models of anisotropic nonlinearly elastic materials, BiomechModel Mechanobiol 12:781-791Pe˜na E, Del Palomar AP, Calvo B, Mart´ınez M, Doblar´e M (2007) Compu-tational modelling of diarthrodial joints. Physiological, pathological andpos-surgery simulations. Arch Comp Methods Eng 14:47-91P´erez del Palomar A, Doblar´e M (2006) On the numerical simulation ofthe mechanical behaviour of articular cartilage. Int J Num Methods Eng67:1244-1271Pierce DM, Trobin W, Raya JG, Trattnig S, Bischof H, Glaser C, HolzapfelGA (2010) DT-MRI based computation of collagen fiber deformation inhuman articular cartilage: a feasibility study. Annals Biomed Eng 38:2447-2463Roy D, Kauffmann C, Delorme S, Lerouge S, Cloutier G, Soulez G(2012) A literature review of the numerical analysis of abdominal aorticaneurysms treated with endovascular stent grafts. Comp Math MethodsMed 2012:820389Sansour C (2008) On the physical assumptions underlying the volumetric-isochoric split and the case of anisotropy. Eur. J. Mech. A/Solids 27:28-3928mith HE, Mosher TJ, Dardzinski BJ, Collins BG, Collins CM, Yang QX,Schmithorst VJ, Smith, M. B. (2001). Spatial variation in cartilage T2 ofthe knee. J Magn Reson Im, 14: 50-55.Sun W, Chaikof EL, Levenston ME. (2008) Numerical approximation of tan-gent moduli for finite element implementations of nonlinear hyperelasticmaterial models. J Biomech eng 130:061003Thury A, van Langenhove G, Carlier SG, Albertal M, Kozuma K, RegarE, Sianos G, Wentzel JJ, Krams R, Slager CJ (2002) High shear stressafter successful balloon angioplasty is associated with restenosis and targetlesion revascularization. Am Heart J 144:136-143Vergori L, Destrade M, McGarry P, Ogden RW (2013) On anisotropic elas-ticity and questions concerning its Finite Element implementation. CompMech 52 :1185-1197Weiss JA, Maker BN, Govindjee S (1996) Finite element implementation ofincompressible, transversely isotropic hyperelasticity. Computer Methodsin Applied Mechanics and Engineering, 135(1-2):107-128, 1996.Wentzel JJ, Gijsen FJ, Stergiopulos N, Serruys PW, Slager CJ, KramsR (2003) Shear stress, vascular remodeling and neointimal formation. JBiomech 36:681-688Xenos M, Rambhia SH, Alemu Y, Einav S, Labropoulos N, Tassiopoulos A,Ricotta JJ, Bluestein D (2010) Patient-based abdominal aortic aneurysmrupture risk prediction with fluid structure interaction modeling. AnnalsBiomed Eng 38:3323-3337
A Consistent Tangent Matrix
To write a UMAT, we need provide the Consistent Tangent Matrix (CTM)of the chosen model. When expressed in terms of Cauchy stress the CTMgiven in Abaqus (2010) may be written as C ijkl = σ ij δ kl + 12 (cid:18) ∂σ ij ∂F kα F lα + ∂σ ij ∂F lα F kα (cid:19) , (20)29hich has both the i ↔ j and k ↔ l minor symmetries.The CTM may estimated using either numerical techniques or an analyt-ical solution. Here we first describe a numerical technique for estimation ofthe CTM. We then present the analytical solution for the MA and HGO-CCTM. Numerical Approximation of the CTM
The CTM may be approximated numerically (Sun et al. (2008)), and a shortoverview is presented here. This numerical approximation is based on alinearised incremental form of the Jaumann rate of the Kirchhoff stress:∆ τ − ∆ W τ − τ ∆ W T = C : ∆ D , (21)where τ is the Kirchhoff stress, ∆ τ is the Kirchhoff stress rate, ∆ D the rate-of-deformation tensor and ∆ W the spin tensor are the symmetric and anti-symmetric parts of the spatial velocity gradient ∆ L (where ∆ L = ∆ FF − ),and C is the CTM.To obtain an approximation for each of components of the CTM, a smallperturbation is applied to (21) through ∆ D . This is achieved by perturb-ing the deformation gradient six times, once for each of the independentcomponents of ∆ D , using∆ F ( ij ) = (cid:15) e i ⊗ e j F + e j ⊗ e i F ) , (22)where (cid:15) is a perturbation parameter, e i is the basis vector in the spatialdescription, ( ij ) denotes the independent component being perturbed.The ‘total’ perturbed deformation gradient is given by ˆF ( ij ) = ∆ F ( ij ) + F . The Kirchhoff stress is then calculated using this perturbed deformationgradient ( τ ( ˆF ij )). The CTM is approximated using C ( ij ) ≈ J (cid:15) ( τ ( ˆF ( ij ) ) − τ ( F )) , (23)where J is the determinant of the deformation gradient. Each perturbationof (23) will produce six independent components. This is performed six timesfor each independent ( ij ), giving the required 6 × nalytical solutions for the MA and HGO-C CTM Here we present an analytical solution for the CTM for the MA and HGOmodels. For convenience we give the volumetric, isotropic and anisotropiccontibutions separately.For the MA model the stress is given by equations (8) and (17). We cancalculate C ijkl from( σ vol ) ij δ kl + ∂ ( σ vol ) ij ∂F kα F lα = κ (2 J − δ ij δ kl , (24)( σ iso ) ij δ kl + ∂ ( σ iso ) ij ∂F kα F lα = µ J − (cid:0) B jl δ ik + B il δ jk − B ij δ kl − B kl δ ij + I δ ij δ kl (cid:1) , (25)( σ aniso ) ij δ kl + ∂ ( σ aniso ) ij ∂F kα F lα = 2 k J − (cid:88) n =4 , ( I n −
1) exp[ k ( I n − ] ( a nj a nl δ ik + a ni a nl δ jk )+ 4 k J − (cid:88) n =4 , [2( I n − k + 1] exp[ k ( I n − ] a ni a nj a nk a nl , (26)where we have used a ni , n = 4 , , i = 1 , ,
3, is the i th component of a n = Fa n .For the HGO-C model the stress is given by equations (8) and (9). Onceagain the isotropic contributions to C ijkl are given by equations (24) and (25).The anisotropic contribution to C ijkl for the HGO-C model is given as:( σ aniso ) ij δ kl + ∂ ( σ aniso ) ij ∂F kα F lα = 4 k J − (cid:88) n =4 , [1 + 2 k (cid:0) I n − (cid:1) ] exp[ k (cid:0) I n − (cid:1) ] × (cid:0) a ni a nj − I n δ ij (cid:1) (cid:0) a nk a nl − I n δ kl (cid:1) + 2 k J − (cid:88) n =4 , ( I n −
1) exp[ k (cid:0) I n − (cid:1) ] ( δ ik a nj a nl + δ jk a ni a nl − δ kl a ni a nj − δ ij a nk a nl + I n δ ij δ kl (cid:1) , (27)where a ni is the i th component of a n = Fa nn