On the Onset of Cracks in Arteries

We present a theoretical approach to study the onset of failure localization into cracks in arterial wall. The arterial wall is a soft composite comprising hydrated ground matrix of proteoglycans reinforced by spatially dispersed elastin and collagen fibers. As any material, the arterial tissue cannot accumulate and dissipate strain energy beyond a critical value. This critical value is enforced in the constitutive theory via energy limiters. The limiters automatically bound reachable stresses and allow examining the mathematical condition of strong ellipticity. Loss of the strong ellipticity physically means inability of material to propagate superimposed waves. The waves cannot propagate because material failure localizes into cracks perpendicular to a possible wave direction. Thus, not only the onset of a crack can be analyzed but also its direction. We use the recently developed constitutive theories of the arterial wall including 8 and 16 structure tensors to account for the fiber dispersion. We enhance these theories with energy limiters. We examine the loss of strong ellipticity in uniaxial tension and pure shear in circumferential and axial directions of the arterial wall. We find that the vanishing longitudinal wave speed predicts the appearance of cracks in the direction perpendicular to tension. We also find that the vanishing transverse wave speed predicts the appearance of cracks in the the direction inclined (non-perpendicular) to tension. The latter result is counter-intuitive yet it is supported by recent experimental observations.


Introduction
In this paper, we show how to predict the onset of cracks in arterial tissue. Our approach requires: (a) development of constitutive models of arteries including a failure description and (b) analysis of the mathematical condition of strong ellipticity for the deformation under consideration.
Constitutive models of the arterial wall have a long lasting history. Arteries are anisotropic and heterogeneous [1] and they undergo large deformations at relatively low stresses. These features require the use of nonlinear elasticity models. Fung with collaborators pioneered these models by proposing an exponential form for strain energy density function in terms of the Green strain tensor [2,3]. They assumed that characteristic material directions coincided with the radial, circumferential and axial directions of artery. Various Fung-type phenomenological theories were developed since then [4][5][6][7][8]. Later on, these models were modified by introducing the so-called structure tensors [9][10][11][12] which provided frameinvariant formulations. Alternatively, more physically appealing structural models were proposed [13][14][15][16][17], which accounted for the angular dispersion of collagen fibers in the strain energy density. Unfortunately, the analytical formulation of the fiber dispersion needed computationally intensive techniques for angular integration on a unit sphere. To compromise between approaches based on structure tensors and analytical dispersion functions the so-called generalized structure tensors were introduced [18,19]. The latter approach was computationally attractive yet it did not allow for an easy exclusion of compressed fibers. In view of the said, a handier approach was developed recently in [20] that used 16 and 8 specially chosen structure tensors to describe the fiber dispersion. This latter approach allowed retaining the physical appeal of the model and, at the same time, easily excluding compressed fibers.
Incorporation of the failure description in the model of the arterial wall is not self-evident. The standard approach of continuum damage mechanics (CDM) is based on the introduction of internal damage variables to decrease the material stiffness [21][22][23]. Evolution of the damage variable is governed by an equation, which describes the gradual degradation of material properties. However, it is not easy to give a physical interpretation of the damage parameter and as a result of that, for example, the use of CDM for the analysis of the condition of strong ellipticity becomes nontrivial. A simpler method to describe material failure, which does not require internal variables was proposed in [24][25][26], which was based on the introduction of the energy limiters in the strain energy functions. This approach is readily suitable for the analysis of the violation of the strong ellipticity condition [27].
We use the recently developed constitutive theories of the arterial wall including 8 and 16 structure tensors to account for the fiber dispersion [20]. We enhance these theories with energy limiters. We examine the loss of strong ellipticity in uniaxial tension in circumferential and axial directions of the arterial wall. We also study the effect of the incompressibility constraint on the analysis of strong ellipticity in uniaxial tension, pure shear and equibiaxial tension. We find that the enforcement of the incompressibility constraint can significantly affect the crack direction.

Constitutive Theory
According to the continuum mechanics approach, the space occupied by a body is continuously filled with the so-called material particles representing on average the assemblies of real physical particles. Material particle placed at x in the initial reference configuration 0 moves to position yðxÞ in the current configuration . The deformation in the neighborhood of the generic material particle is described by the deformation gradient: F ¼ Grady ¼ @y=@x: We assume that the mechanical response of arterial wall is hyperelastic. Arterial wall consists of collagen fibers dispersed in isotropic ground matrix. The strain energy function W of the intact artery wall comprises two terms: is the neo-Hookean strain energy for the isotropic ground matrix and is the energy of the dispersed fibers with the integration on a unit sphere.
Here c and K are the shear and bulk modulus respectively; I 1 ¼ trC and I 3 ¼ det C, where C ¼ F T F is the right Cauchy-Green strain tensor; Ä is a solid angle; is the angular density of the fiber distribution, which is normalized as R dÄ ¼ 4; and w is the strain energy of a single fiber per unit reference volume.
In the case where the incompressibility condition is imposed, we assume I 3 ¼ 1 and the second term in (2) vanishes. We designate a unit vector in the direction of a generic material fiber in the initial configuration as where 0 È 2 and 0 Â .
Note that È is the angle in the tangent plane marked by the local unit vectors e 1 and e 2 denoting the circumferential and axial directions of the arterial wall respectively. È is measured from the circumferential direction e 1 . Â is the angle in the plane perpendicular to this tangent plane, which is measured from the radial direction e 3 in this plane.
The form for the angular distribution function of fibers is chosen following [20,28], based on the experiments by [29], where ip ðÈÞ and op ðÂÞ are the in-plane and out-of-plane dispersion functions defined as follows where material constants are a ¼ 2:54, b ¼ 19:44 and ¼ 47:99 is the angle between the circumferential direction e 1 and the mean fiber direction. The strain energy of single fiber per unit reference volume is chosen as an exponential function of the form [9], where k 1 and k 2 are material parameters; I 4 ¼ C : a a > 1; and triangular brackets denote Macaulay brackets to exclude the fiber response in compression Note that a a is called structure tensor.
Unfortunately, the strain density function of the dispersed fibers (3) can rarely be calculated analytically and a numerical integration procedure should be used. Such numerical integration is called the cubature formula and it can be written in the following form [20] f ¼ where ðiÞ is the weight coefficient and w ðiÞ ¼ wða ðiÞ a ðiÞ Þ; a ðiÞ ¼ aðÈ ðiÞ ; Â ðiÞ Þ; where integration points ðÈ ðiÞ ; Â ðiÞ Þ are taken on the unit sphere. Note that a finite number of structure tensors a ðiÞ a ðiÞ account for anisotropy induced by dispersed fibers.
In introducing the strain energy we should note that the number of physical particles in a representative volume of real materials is limited. Consequently, the bond energy of these particles is limited. Thus, the strain energy density on the macroscopic scale should also be limited. Obviously, a limited strain energy automatically provides a limited stress, which a material can sustain. The latter condition of the limited or bounded stress is a qualitative indicator of the material failure.
The limiter can be introduced in the strain energy in the following form [26], for example, where and Àðs; xÞ ¼ R 1 x t sÀ1 e Àt dt is the upper incomplete gamma function. Here, f is the failure energy; e ðFÞ is the elastic energy; 1 is identity tensor; is the energy limiter (average bond energy); and m is a material parameter.
We emphasize that the present formulation is valid for the analysis of the onset of failure. If the failure localization and propagation are of interest, then a regularized formulation should be used as in [30,31], for example.

Strong Ellipticity Condition
The motion and deformation of the body is described by equations of momenta balance and the constitutive law where 0 is the referential mass density; P is the first Piola-Kirchhoff stress tensor; ðDivPÞ i ¼ @P ij =@x j ; is the strain energy; Å is the Lagrange multiplier enforcing the incompressibility condition det F ¼ 1 if necessary.
We emphasize that braces f…g are used for the terms that should be counted for incompressible material models only.
Designating increments with tildes it is possible to derive incremental equations of momenta balance [25] and the incremental constitutive law P ¼ ð@ 2 =@F@FÞ :F þ fÅF ÀTFT F ÀT ÀÅF ÀT g; fF : F ÀT ¼ 0g: (16) Alternatively, these incremental equations can be reformulated in the Eulerian form where the current configuration is referential and σ ¼ A :L þ fÅL T ÀÅ1g; fL : 1 ¼ 0g; PF T is the Cauchy stress tensor and ðdivσÞ i ¼ @σ ij =@y j ;σ ¼ I À1=2 3P F T is the incremental Cauchy stress;L ¼FF À1 is the incremental velocity gradient; A is the fourth order elasticity tensor with Cartesian components For the strain energy defined by (11), we further calculate and Substitution of (21) in (19) yields F js F lm @ 2 W @F is @F km À mW mÀ1 Àm @W @F km @W @F is exp½ÀW m Àm : We look for a plane wave solution of the incremental initial-boundary-value problem y ¼ rgðs Á y À vtÞ; fÅ ¼ Çg 0 ðs Á y À vtÞ; g where r and s are the unit vectors in the directions of wave polarization and wave propagation respectively; v is the wave speed; g 0 denotes the differential of g with respect to the argument of the function.
Substituting forÅ andL ¼ gradỹ ¼ @ỹ=@y from (23) to (18) 1 , we get the incremental stressσ. Then, substituting forỹ from (23) andσ to the linear momentum balance (17) 1 , we get where ÃðsÞ is the acoustic tensor with Cartesian components Taking scalar product of (24) with r, we obtain for the wave speed where According to the assumptions of the previous section the intact strain energy function W ðI 1 ; I 3 ; I ðrÞ 4 Þ depends on the first and third principal invariants as well as a number of (pseudo) invariants I ðrÞ 4 ¼ C : a ðrÞ a ðrÞ > 1 (no sum over r). Assuming that mixed derivatives of the intact strain energy @ 2 W =@I 1 @I ðrÞ 4 and @ 2 W =@I 3 @I ðrÞ 4 vanish, we can derive factors f 3 and f 4 analytically as follows and where . The positive wave speed corresponds to the mathematical condition of the strong ellipticity of the incremental initial boundary-value-problem. Zero wave speed mathematically means violation of the strong ellipticity condition and, physically, it means inability of the material to propagate a wave in direction s. The latter notion can also be interpreted as the onset of a crack perpendicular to s.
In the present work, we will consider longitudinal wave (P-wave) and transverse wave (S-wave) in the plane of the arterial sheet for the vanishing wave speed.
Denoting the unit tangent vector in the circumferential and axial directions of the arterial wall by e 1 and e 2 respectively, we can write for the P-wave and s ¼ cos e 1 þ sin e 2 ; r ¼ À sin e 1 þ cos e 2 (32) for the S-wave, where is unknown angle in the tangent plane.

Specialization of Material
We use the material models experimentally calibrated in [20]. These models describe the intact material behavior and we enhance them with a failure description via energy limiters.
The first model includes 16 structure tensors with out of plane fiber dispersion and its parameters are given in Tabs.1 and 2.
Uniaxial tension in circumferential and axial directions of the arterial wall for this model are shown in Fig. 1. The second model includes 8 structure tensors without out of plane fiber dispersion and its parameters are given in Tabs. 3

and 4.
Uniaxial tension in circumferential and axial directions of the arterial wall for this model are shown in Fig. 2.

Results
In this section, we report the results of the analysis of the loss of strong ellipticity for the vanishing wave speed in the following three deformation modes: uniaxial tension in circumferential and axial directions; pure shear in circumferential and axial directions; equibiaxial tension.
We consider the "stiff" displacement-controlled loading. For a plane sheet of arterial wall, the unknown out-of-plane principal stretch 3 is determined in terms of the known in-plane stretches 1 and 2 assuming zero out-of-plane stresses. Three deformation modes in circumferential direction are described as follows: 2 ¼ À0:5 1 in uniaxial tension; 2 ¼ 0 1 ¼ 1 in pure shear; and 2 ¼ 1 1 ¼ 1 in equibiaxial tension. Three deformation modes in axial direction are described as follows: 1 ¼ À0:5 2 in uniaxial tension; 1 ¼ 0 2 ¼ 1 in pure shear; and 1 ¼ 1 2 ¼ 2 in equibiaxial tension. From the condition of the vanishing wave speed: I 1=2 3 v 2 ¼ f 1 f 2 ¼ 0; we find the critical stretches that determine the loss of strong ellipticity. We consider longitudinal P-waves and transverse S-waves in slightly compressible and incompressible (for S-waves only) material models. The results obtained for the slightly compressible and incompressible materials are numerically very close.

Discussion
In the present work, we studied a possible onset of cracks in the arterial wall. For that purpose, we developed two constitutive models with 16 and 8 structure tensors to account for anisotropy and failure of the wall. The intact material behavior was calibrated based on the experimental data for adventitia and the energy limiters were introduced to describe failure.
These models were used in analysis of the loss of strong ellipticity in uniaxial tension and pure shear in circumferential and axial directions of the artery and in biaxial tension. Directions of possible cracks were obtained from the condition of the vanishing speed of the superimposed longitudinal (P-) and transverse (S-) waves. The main findings can be summarized as follows: 1. The vanishing P-wave speed predicts the appearance of cracks in the direction perpendicular to tension in uniaxial tension and pure shear. Remarkably, such prediction would be suppressed by the imposition of  the incompressibility constraint, which is often used in calculations. This constraint sorts out the longitudinal wave and, consequently, we lose the possibility to predict cracks associated with it. The incompressibility constraint becomes a Trojan Horse for the analysis of the onset of cracks. Similar results for purely isotropic soft material have been reported in [32]. 2. The vanishing S-wave speed predicts the appearance of cracks in the the direction inclined (nonperpendicular) to tension in uniaxial tension and pure shear. Such cracks are counter-intuitive at first glance. However, the the recent experimental work [33] reported an observation of cracks in the direction of tension in a silicone elastomer. The authors of the work attributed these "sideways" cracks to "microstructural anisotropy (in a nominally isotropic elastomer)". 3. Equibiaxial stretching can lead to the appearance of cracks in any direction despite the anisotropy of material. The inclined cracks oriented along the bundles of collagen fibers have been found in experiments [34].
We should note that the presented methodology can possibly give new insights in experiments with cracking. The character and direction of cracks can provide information about material anisotropythe inverse problem. This circle of questions is beyond the scope of the present work, however.
Finally, we emphasize that the proposed methods are valid for the analysis of the onset of cracks only. Tracking the crack development would require regularized formulations (e.g., [30,31]), which were not considered in this work.