Ansys Ogden Model Tissue

Abiomedical Engineering,

bSchool of Mathematics, Statistics and Applied Mathematics,

Abstract

The Holzapfel–Gasser–Ogden (HGO) model for anisotropic hyper- elastic behaviour of collagen fibre reinforced materials was initially developed to describe the elastic properties of arterial tissue, but is now 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 HGO model is a function of one isotropic and two anisotropic deformation invariants. A compressible form (HGO-C model) is widely used in fi- nite element simulations whereby the isotropic part of Ψ is decoupled into volumetric and isochoric parts and the anisotropic part of Ψ is expressed in terms of isochoric invariants. Here, by using three simple deformations (pure dilatation, pure shear and uniaxial stretch), we demonstrate that the compressible HGO-C formulation does not cor- rectly model compressible anisotropic material behaviour, because the anisotropic component of the model is insensitive to volumetric defor- mation due to the use of isochoric anisotropic invariants. In order to

Arxiv:2009.08754V1 [Cond-Mat.Soft] 18 Sep 2020

correctly model compressible anisotropic behaviour we present a mod- ified anisotropic (MA) model, whereby the full anisotropic invariants are 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. It also computes the correct anisotropic stress state for pure shear and uniaxial 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 computed when the full anisotropic invariants are used (MA model) instead of the isochoric form (HGO-C model).

Keywords Anisotropic, Hyperelastic, Incompressibility, Finite element,

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

F – Deformation Gradient

J – determinant of the deformation gradient; local ratio of volume change

I1 – First Invariant Of C

I4,6 – 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

I4,6 – Isochoric Anisotropic Invariants

a0i, i = 4, 6 – unit vector aligned with a reinforcing fibre in the reference

Configuration

ai, i = 4, 6 – updated (deformed) fibre direction (= Fa0i)

Ν – Isotropic Poisson’S Ratio

Bold uppercase symbols represent second order tensors, bold lowercase sym- bols represent vectors and un-bold symbols represent scalars.

Introduction

The anisotropic hyperelastic constitutive model proposed by Holzapfel et al. (2000) (henceforth referred to as the HGO model) is used extensively to model collagen fibre-reinforced biological materials, even more so now that it has 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 transversely isotropic 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 for the anisotropic contributions of the reinforcing fibre families, in addition to a term representing the isotropic contribution of the ground matrix in which the fibres are embedded. Also, it assumes that the collagen fibres do not support compression, and hence they provide a mechanical contribution only when in tension (this may be taken care of by pre-multiplying each anisotropic 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 I1) and two isochoric anisotropic invariants (denoted I4 and I6). A Lagrange multiplier is used to enforce incompressibility (Holzapfel et al., 2000). Once again it should be stressed that the original HGO model is intended only for the 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 penalty term with an isotropic hydrostatic stress term that depends on a user speci- fied bulk modulus. This modification allows for the relaxation of the incom-

3

pressibility condition and we therefore refer to this modified formulation as the HGO-C (compressible) model for the remainder of this study. The HGO-C model has been widely used for the finite element simulation of 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 using the HGO-C model (e.g. Pe˜na et al. (2007) used a Poisson’s ratio, ν = 0.1 and P´erez del Palomar et al. (2006) used ν = 0.1 and ν = 0.4). To date, material compressibility of arterial tissue has not been firmly established.

Incompressibility was assumed by the authors of the original HGO model and in subsequent studies (e.g. Kiousis et al. (2009)). However many studies model arteries as compressible or slightly compressible (e.g. Cardoso et al.

(2014) ν = 0.33 −0.43 and Iannaccone et al. (2014) ν = 0.475). In addition to arterial tissue the nucleus pulposus of an inter-vertebral disc has been modelled as a compressible anisotropic material using the HGO-C model (e.g. (Maquer et al., 2014) ν = 0.475). Furthermore the HGO-C formulation has been used to simulate growth of anisotropic biological materials, where volume change is an intrinsic part of a bio-mechanical process (e.g. (Huang et al., 2012) ν = 0.3). However, the enforcement of perfect incompressibility may not be readily achieved in numerical models. As an example, the fi- nite element solver Abaqus/Explicit assigns a default Poisson’s ratio of 0.475 to ”incompressible” materials in order to achieve a stable solution (Abaqus, 2010) and in this case the HGO-C model must be used (e.g. Conway et al. (2012); Famaey et al. (2012)). Despite the widespread use of the HGO- C model, its ability to correctly simulate anisotropic compressible material behaviour has not previously been established.

• The first objective of this study is to demonstrate that the HGO-C formulation does not correctly model an anisotropic compressible hy- perelastic material.

Recently, Vergori et al. (2013) showed that under hydrostatic tension, a sphere consisting of a slightly compressible HGO-C material expands into a larger sphere instead of deforming into an ellipsoid. It was suggested that this effectively isotropic response is due to the isochoric anisotropic invariants Ii being used in the switching function instead of the full invariants Ii, i = 4, 6.

However, in the current paper we show that the problem emerges fundamen- tally because there is no dilatational contribution to the anisotropic terms of

4

Ψ. In fact, modifying only the “switching criterion” for fibre lengthening is not sufficient to fully redress the problem. • The second objective of the study is to implement a modification of the HGO-C model so that correct anisotropic behaviour of compressible materials is achieved.

This modified anisotropic (MA) model uses the full form of the anisotropic invariants and through a range of case studies we show this leads to the correct 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 component of the HGO-C model to volumetric deformation in compressible materials.

We demonstrate that the modification of the model to include the full form of the anisotropic invariants corrects this deficiency. In Section 3 we show how 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 the modified model. We devote Section 5 to two Finite Element biomechanics case studies, namely pressure expansion of an artery and stent deployment in an artery, and illustrate the significant differences in computed results for the HGO-C model and the modified model. Finally, we provide some concluding remarks and discussion points in Section 6.

Hgo-C Model For Compressible Materials

The original HGO model is intended for incompressible materials. However a variation of the HGO model whereby a bulk modulus is used instead of a penalty term has been implemented in a number of FE codes. Several authors have used this formulation t model compressible anisotropic materials bu using a relatively low value of bulk modulus. An important objective of this paper is to highlight that this HGO-C formulation does not correctly model compressible 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.

5

The right Cauchy–Green tensor is defined by C = FTF, where T indicates the transpose of a second-order tensor. Hyperelastic constitutive models used for rubber-like materials often split the local deformation into volume-changing (volumetric) and volume-preserving (isochoric, or deviatoric) parts. Accordingly the deformation gradient F is

(1)

where J is the determinant of F. The term in the brackets represents the volumetric 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 within which are embedded two families of fibres characterized by two preferred directions in the reference configuration defined in terms of two unit vectors a0i, i = 4, 6. With C, J and a0i are defined the invariants

(3)

where Ii (i = 1, 4, 6) are the isochoric counterparts of Ii. The HGO model proposed by Holzapfel et al. (2000) for collagen reinforced soft tissues ad- ditively splits the strain energy Ψ into volumetric, isochoric isotropic and

(4)

where Ψiso and Ψaniso are the isochoric isotropic and isochoric anisotropic free-energy contributions, respectively, and C = J−2/3C is the isochoric right Cauchy–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

(5)

where κ0 and µ0 are the bulk and shear moduli, respectively, of the soft isotropic matrix. Of course one may write (5) in terms of the full invariants also, using the results from (3).

6

The isochoric anisotropic free-energy term is prescribed as

(6)

where k1 and k2 are positive material constants which can be determined from experiments. For a general hyperelastic material with free energy Ψ the Cauchy stress

(7)

For the Cauchy stress derived from Ψ above, we have the decomposition

(9)

where ai = Fa0i. This slightly compressible implementation is referred to as the 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 for incompressible materials. In that limit, κ0 →∞, (J −1) →0 while the product of these two quantities becomes an indeterminate Lagrange multi- plier, p, and the volumetric stress assumes the form, σvol = −pI. Indeed the original incompressible HGO model can equally be expressed in terms of the full invariants I4 and I6 (with J →1) (e.g., Holzapfel et al. (2004)).

However, in the case of the HGO-C implementation, if κ0 is not fixed numerically at a large enough value, then slight compressibility is intro- duced into the model.

The Key Point Of This Paper Is That The Isochoric

anisotropic term Ψaniso defined in (6) does not provide a full representation of the anisotropic contributions to the stress tensor for slightly compress- ible materials.

In Section 2.4 we introduce a simple modification of the anisotropic term to account for material compressibility.

Pure Dilatational Deformation

First we consider the case of the HGO-C material subjected to a pure dilata-

C = Λ2I,

J = λ3.

(10)

We expect that an anisotropic material requires an anisotropic stress state to maintain the pure dilatation. However, calculation of the invariants Ii and

(11)

so that while Ii is indeed the square of the fibre stretch and changes with the magnitude of the dilatation, its isochoric counterpart Ii is always unity. Referring to (9), it is clear that the entire anisotropic contribution to the stress (7) disappears (i.e. σaniso ≡0), and the remaining active terms are the isotropic ones. Thus, under pure dilatation, the HGO-C model computes an entirely isotropic state of stress.

Applied Hydrostatic Stress

Now we investigate the reverse question: what is the response of the HGO-C

(12)

where σ > 0 under tension and σ < 0 under pressure? In an anisotropic material, we expect the eigenvalues of C, the squared principal stretches,

2, Λ2

3 say, to be distinct. Hence, if the material is slightly compressible, then a sphere should deform into an ellipsoid (Vergori et al., 2013) and a cube should deform into a hexahedron with non-parallel faces (N´ı Annaidh et al., 2013b).

However, in the HGO-C model the Ψaniso contribution is switched on only when Ii (not Ii) is greater than unity. Vergori et al. (2013) showed that in fact Ii is always less than or equal to one in compression and in expansion under hydrostatic stress, so that the HGO-C response is isotropic, contrary to physical expectations. Then we may ask if removal of the switching function circumvents this problem so that anisotropic response is obtained.

8

With the fibres taken to be mechanically equivalent and aligned with a04 = (cos Θ, sin Θ, 0) and a06 = (cos Θ, −sin Θ, 0) in the reference config- uration, we have, by symmetry, I6 = I4 and I6 = I4 and Ψ6 = Ψ4, where the subscripts 4 and 6 on Ψ signify partial differentiation with respect to I4 and I6, respectively. Similarly, in the following the subscript 1 indicates differentiation with respect to I1. For this special case, Vergori et al. (2013) showed that the stretches arising from the application of a hydrostatic stress

6

.

I4 −1

2].

(14)

Looking at (13), we see that there is a solution to the hydrostatic stress problem where the stretches are unequal, so that a sphere deforms into an ellipsoid. However, there is also another solution: that for which I4 ≡1, in which case, Ψ4 ≡0 by the above equation, and then λ1 = λ2 = λ3 = J1/3 by (13). Thus, a sphere then deforms into another sphere.

Of those (at least) two possible paths, FE solvers converge upon the isotropic solution. One possible explanation for this may be that the ini- tial computational steps calculate strains in the small-strain regime. In that regime, Vergori et al. (2013) showed that all materials with a decoupled vol- umetric/isochoric free-energy behave in an isotropic manner when subject to a hydrostatic stress. Hence the first computational step brings the de- formation on the isotropic path, and I4 = 1 then, and subsequently.

In

Section 2.2 and Section 2.3 we have thus demonstrated that the use of an isochoric form of the anisotropic strain energy Ψaniso from the HGO model in the HGO-C model cannot yield a correct response to pure dilatation or applied hydrostatic stress.

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, a04, a06) = Ψvol (J) + Ψiso(J, C) + Ψaniso(C, a04, a06),

(15)

where the expressions for strain energy density terms Ψvol and Ψiso are

I=4,6

{exp[k2 (Ii −1)2] −1}.

(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 Stress

for the MA model is determined using (7) and the decomposition σ =

I=4,6

(Ii −1) exp[k2(Ii −1)2]ai⊗ai.

(17)

where ai = Fa0i, i = 4, 6. Now it is easy to check that in the cases of a pure dilatation and of a hydrostatic stress, the MA model behaves in an anisotropic manner, because the term Ii −1̸ = 0 and hence Ψaniso̸ = 0 and σaniso̸ = 0. 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 the consistent 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 expansion of a sphere under hydrostatic tension of Vergori et al. (2013), this time using the MA formulation. Again two families of fibres are assumed, lying in the (1, 2) plane and symmetric about the 1-axis (the sphere and axes are shown

Figure 1:

A) Schematic of an undeformed sphere highlighting three radii on 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/r0) of the orthogonal radii for both MA and HGO-C models versus the ratio σhyd/σmax

Hyd . Note That The De-

formation computed for the HGO-C model incorrectly remains spherical. C) Deformed ellipsoidal shape computed for the MA model; contours illustrate the inhomogeneous distribution of stress triaxiality (σhyd/q) throughout the deformed body.

in Figure 1A). The displacements of points on the surface of the sphere at the ends of three mutually orthogonal radii with increasing applied hydrostatic tension are shown in Figure 1B. Clearly the sphere deforms into an ellipsoid with a major axis oriented in the 3-direction and a minor axis oriented in the 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(σ)/3 is the hydrostatic

P

3/2 σ′ : σ′ is the von Mises equivalent stress, σ′ being the deviatoric Cauchy stress tensor. Clearly an inhomogeneous stress state is computed 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 Figure 1B for comparison. In that case a similar fibre-reinforced sphere is shown to deform into a larger sphere with a homogeneous stress distribution, indicative of isotropic material behaviour.

Analysis Of Pure Shear

A pure dilatation and a hydrostatic stress each represent a highly idealized situation, unlikely to occur by themselves in soft tissue in vivo. This sec- tion highlights the unphysical behaviour can also emerge for common modes of deformation if the anisotropic terms are based exclusively on the iso- choric invariants. Considering once again the general case of a compressible anisotropic material, we analyse the response of the HGO-C and MA models to pure in-plane shear. Regarding the out-of-plane boundary conditions, we first consider the case of plane strain (Section 3.1). Even though this defor- mation is entirely isochoric the HGO-C model yields incorrect results. We then consider the case of plane stress (Section 3.2), and again demonstrate that the HGO-C model yields incorrect results. By contrast, we show that the MA model computes a correct stress state for all levels of compressibility and specified deformations. In the following calculations we assume a shear modulus, µ0 = 0.05 MPa and anisotropic material constants k1 = 1 MPa and k2 = 100.

Plane Strain Pure Shear

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 Cauchy stress. In particular, we take the deformation gradient for this deformation

(18)

where F12 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 into a parallelogram symmetric about a diagonal of the square. The deformation

12 + 1 −F12. We Can Think Of The Deformation

arising from displacement components applied to the vertices of the square, as indicated in Figure 2A. Two families of fibres, with reference unit vectors a04 and a06 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 a range of fibre orientations for each of the HGO-C and MA models.

Figure 2:

A) Schematic illustrating the kinematics of the pure shear de- formation of the (1, 2) section of a unit cube. Note the rotated coordinate system (1′, 2′), orientated at 45◦to the (1, 2) axes, used to specify the vertex displacement components u1′ and u2′. Note also the vectors a0i, i = 4, 6, indicating the directions of the two families of fibres, with angle θ. Results are displayed for a range of fibre orientations with θ from ±45◦to ±90◦with respect to the (1, 2) coordinate system. B) Computed stress ratio σ33/σ12 versus F12 for the HGO-C model, illustrating significant negative (compres- sive) stresses in the out-of-plane direction. C) Computed stress ratio versus F12 for the MA model, illustrating very small negative (compressive) stresses in 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 of the HGO-C and the MA models coincide (because J = 1 and hence I4 = I4), the corresponding stress tensors are very different. This is due to the “deviatoric” form of the anisotropic stress contribution that emerges for the HGO-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 stress component σ33 which is comparable in magnitude to σ12, as shown in Figure 2B. Such a negative stress is anomalous in the sense that for large κ0 the result for the incompressible limit should be recovered, but it is not. Indeed, if we start with the incompressible model we obtain σ33 = µ0 −p, which is

13

Figure 3: Dimensionless stress components σij/k1 versus F12 for the case of a single family of fibres orientated at θ = 30◦. A) MA model; B) HGO-C model.

independent of σ12. However, as (18) represents a kinematically prescribed isochoric deformation, the volumetric stress in the HGO-C model goes to zero and does not act as the required Lagrange multiplier.

By contrast, the out-of-plane compressive normal stress component σ33 computed for the MA model is at least an order of magnitude lower than the in-plane shear stress component σ12 (Figure 2C), and is close to zero for most fibre orientations. This is consistent with the incompressible case because, since p is arbitrary it may be chosen to be µ0 so that σ33 = 0. This is what might be expected physically, given that the fibres and the deformations are confined to the (1, 2) plane.

Because of the deviatoric component of the stress tensor emerging from the HGO-C model, the trace of the Cauchy stress is always zero when J = 1 as equations (8) and (9) will confirm. By contrast, the trace of the Cauchy stress is not zero for the MA model. Hence the in-plane stress components are significantly different from those for the HGO-C model, as shown in Figures 3A and 3B, respectively, for the case of a single fibre family with θ = 30◦.

Plane Stress Pure Shear

The kinematically prescribed isochoric deformation in Section 3.1 is volume conserving and makes the Ψvol terms equal to zero. We modify the out-of- plane boundary condition to enforce a plane stress (σ33 = 0) simulation This allows a compressible material to deform of out-of-plane.

(19)

where the out of plane stretch component F33 in general is not equal to 1, so that the deformation is not in general isochoric. If the bulk modulus κ0 is very large compared with the initial shear modulus µ0, then it acts as a Lagrange multiplier to enforce incompressibility, such that F33 = 1 (at least approximately). If the magnitude of the bulk modulus is reduced, then the material becomes slightly compressible and F33̸ = 1. Here we investigate the sensitivity of the stress computed for the HGO-C and MA models to the magnitude of the bulk modulus κ0.

First, we consider the almost incompressible case where the ratio of bulk to shear modulus is κ0/µ0 = 2×106 for the isotropic neo-Hookean component of the model, equivalent to a Poisson ratio of ν = 0.49999975. The stress components are shown in Figure 4A. An important point to note is that in this case the deformation is effectively isochoric, because we find J = F33 = 1.00006, and yet the HGO-C model predicts an entirely different stress state from that for the kinematically constrained isochoric deformation of the previous section shown in Figure 3B. This is because the volumetric term of the free energy now contributes to the trace of the stress tensor, and therefore the high magnitude of bulk modulus effectively acts as a Lagrange multiplier to enforce incompressibility.

Indeed For These Conditions The Hgo-C And

MA models behave identically to the original HGO model. However, unlike the HGO-C model, the MA model computes identical stress components for both the kinematically constrained isochoric deformation (18) and for the Lagrange multiplier enforced volume preserving deformation (19).

If the incompressibility constraint is slightly relaxed, so that κ0/µ0 = 50 (ν = 0.490) the HGO-C model computes a very different stress state, as shown in Figure 4B, with stress components being reduced by an order of magnitude. Thus the HGO-C model is very sensitive to changes in the bulk

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.