Stability of Membrane Elastodynamics with Applications to Cylindrical Aneurysms

Abstract

The enlargement and rupture of intracranial and abdominal aortic aneurysms constitutes a major medical problem. It has been suggested that enlargement and rupture are due to mechanical instabilities of the associated complex fluid-solid interaction in the lesions. In this paper, we examine a coupled fluid-structure mathematical model for a cylindrical geometry representing an idealized aneurysm using both analytical and numerical techniques. A stability analysis for this subclass of aneurysms is presented. It is shown that this subclass of aneurysms is dynamically stable both with and without a viscoelastic contribution to the arterial wall.

1. Introduction

An arterial aneurysm is a focal dilatation or lesion of the vascular wall. It is thought that two to five percent of the general population of the Western world harbors ruptured or unruptured saccular aneurysms. Such aneurysms may remain dormant for years or even decades. These lesions are more common in women than in men, and they tend to manifest themselves during the fifth to seventh decades of life. The mean age of rupture is around 52 years old [1].

Rupture of intracranial aneurysms is the leading cause of spontaneous subarachnoid hemorrhage (SAH), which is fatal in 35–50% of cases. Many of the survivors suffer functional deficits [1]. There are approximately 27,000 patients per year in the US with reported ruptured aneurysms, or about 10 in every 100,000 persons per year [2].

The three general types of aneurysms are lateral aneurysms, fusiform aneurysms, and saccular aneurysms [2]. A saccular aneurysm is a pouch-like sac of blood attached to an artery or a branch of a blood vessel by a neck or stem. A lateral aneurysm is a bulge on one wall of the blood vessel while a fusiform aneurysm bulges along all walls of the aneurysm [2]. Fusiform aneurysms usually occur in the basilar artery and become symptomatic by pressing on adjacent tissue. Saccular aneurysms usually occur at the apex of a bifurcation in or near the circle of Willis and often remain asymptomatic until they rupture [1]. Aneurysms can also be classified by size. Small aneurysms are less than 11 millimeters in diameter, larger aneurysms are 11–25 millimeters in diameter, and giant aneurysms are greater than 25 millimeters in diameter [2]. Intracranial aneurysms are usually saccular and are often modeled using spherical geometry [3]. However, abdominal aortic aneurysms are usually fusiform and should be modeled using nonspherical geometry. This motivated us to consider cylindrical geometry, which is the focus of this paper.

Much like intracranial aneurysms, abdominal aortic aneurysms are focal dilatations of the arterial wall. They are most problematic in the elderly population, and the vast majority of them are of the fusiform variety [1]. They also tend to be much larger than their intracranial counterparts. Abdominal aortic aneurysms affect 6–9% of people aged 65 or older in the industrialized world. The rate of rupture has been estimated as 3–9% per year, with mortality rates between 65% and 85% [4].

While all aneurysms have the potential to rupture, most intracranial aneurysms show no symptoms unless they become very large or rupture. Many aneurysms go unnoticed until they rupture or are detected by brain imaging, possibly obtained for another condition [2]. In a meta-analysis of a number of studies, the overall rupture risks were 1.2% for a followup within 5 years of detection, 0.6% for a followup 5 to 10 years after detection, and 1.3% for a followup of 10 years or more after detection [5]. The decision whether to treat the aneurysm is currently made based primarily on the size of the lesion. If detected before rupture, there are several treatment options, including clipping the aneurysm to prevent further blood flow into the lesion and endovascular embolization to fill the aneurysm with coils or latex balloons [2]. For abdominal aortic aneurysms, surgery to remove the lesion is also possible. However, this is a high-risk procedure with a 5% mortality rate [6]. Since these treatment options all have associated risks, the need for better information about the causes of aneurysmal growth and rupture to guide treatment decisions is evident.

There are several major hypotheses which have been proposed for why aneurysms rupture. The presence of limit point instabilities, that is, mathematical bifurcations in the aneurysm's quasistatic response to increased distension pressure, was an early hypothesis for a cause of rupture. However, the studies which suggested such instabilities, such as Austin et al. [7] and Akkas [8], either used elastomeric membranes or neo-Hookean constitutive relations, where one should expect such instabilities. Since collagenous soft tissues tend to exhibit an exponential behavior, rather than a rubber-like response, these results need to be revisited [9]. Another rupture hypothesis is the possibility of dynamic instabilities, such as turbulence or resonance within the aneurysm. However, most older studies of the elastodynamics were based on classical shell theory and are thus linearized models [9]. Moreover, the geometries considered in the few nonlinear models for aneurysms were spherical. Analysis of nonlinear elastodynamics in aneurysms along with fluid-structure coupling was first performed by Shah and Humphrey [10] for a spherical geometry with homogeneous, isotropic properties. A viscoelastic wall response was considered in a later paper by David and Humphrey [11], also for a spherical geometry with homogeneous, isotropic properties.

In this paper, we develop hyperelastic and viscoelastic membrane models which incorporate fluid-structure interactions for a cylindrical geometry. These models are for aneurysms with homogeneous, isotropic properties undergoing radial inflation. We also perform stability analysis on the governing equations derived for the models. We then use numerical techniques to further examine and validate the analytical results, as well as comparing the numerical results with the known results for a pulsating sphere.

2. Model and Governing Equations

Membranes often imply a thin layer of tissue covering a surface or dividing a space. Some examples include cell membranes and walls covering major organs in the body. In mechanics, the word membrane refers to a thin flexible structure which offers negligible resistance to bending. As a result, the in-plane stresses are assumed to be constant throughout the thickness of the membrane. The mathematical theory and governing equations for membrane deformations can be derived from the general theory of plates and shells [12, 13].

For simplicity, consider Figure 1. It illustrates cross-sections of a cylindrical blood vessel consisting of three components, namely, the arterial wall, cerebrospinal fluid (CSF), and blood flow. Next we will describe the mathematical models for each of these components and derive a radial inflation model for a cylindrical aneurysm analogous to a spherical inflation model derived by Shah and Humphrey [10].

(a)

(b)

Figure 1 

Membrane cross-sections.

Let the stretch ratio of the aneurysmal wall be denoted 𝜆=𝑟/𝑅 where 𝑟 and 𝑅 are the deformed and undeformed radii of the aneurysm. Also let ℎ and 𝐻 be the deformed and undeformed thicknesses of the aneurysm wall. We consider radial motion of the inner and outer fluids. Thus, we assume radial velocity is a function of only time 𝑡 and radial position 𝑟, and axial velocity is constant, that is, 𝑢=𝑢(𝑟,𝑡) and 𝑤≡𝑤₀, where 𝑢 is the radial velocity, 𝑣 is the circumferential velocity, 𝑤 is the axial velocity, and the velocity vector is 𝐯=(𝑢,𝑣,𝑤). In this work, we will also assume the CSF is a viscous, incompressible, Newtonian fluid.

2.1. Model for the CSF

Next, we describe the governing equations, derived from the conservation of mass and conservation of linear momentum. Assuming a flow in the radial direction 𝜉 in the cylindrical domain, the conservation of mass leads to the continuity equation ∇⋅𝐯=0, which in cylindrical coordinates (𝑟,𝜃,𝑧) is1𝜉𝜕()+1𝜕𝜉𝜉𝑢𝜉𝜕𝑣+𝜕𝜃𝜕𝑤𝜕𝑧=0. (2.1)

Assuming that the axial velocity is independent of 𝑧 and the circumferential velocity is zero, we then have 𝑣=0 and 𝜕𝑤/𝜕𝑧=0. Equation (2.1) becomes1𝜉𝜕()𝜕𝜉𝜉𝑢=0, (2.2)

which yields𝑔(𝑡)𝑢=𝜉, (2.3) where 𝑔(𝑡) is a function only depending on time. Since any particle on the membrane wall will have the same radial velocity as an adjacent fluid particle in the cerebrospinal fluid, we match the velocities on the outer boundary of the aneurysm to yield𝜕()𝜕𝑡𝑟−𝑅=𝑅𝜕𝜆=𝑔(𝑡)𝜕𝑡𝑟. (2.4)

Solving for 𝑔(𝑡), we have𝑔(𝑡)=𝑅²𝜆𝜕𝜆.𝜕𝑡 (2.5)

Combining (2.3) and (2.5) gives𝑢()=𝑔(𝑡)𝜉,𝑡𝜉=𝑅²𝜆𝜉𝜕𝜆.𝜕𝑡 (2.6)

The radial velocity 𝑢 satisfies the Navier-Stokes radial equation given by𝜕𝑢𝜕𝑡+𝑢𝜕𝑢1𝜕𝜉=−𝜌_𝑓𝜕𝑝𝜕𝜕𝜉+𝜈²𝑢𝜕𝜉²+1𝜉𝜕𝑢−𝑢𝜕𝜉𝜉², (2.7) where 𝑝 is the pressure and 𝜌_𝑓 and 𝜈 are the density and the kinematic viscosity, respectively, of the surrounding fluid. Here the fluid pressure 𝑃=𝑃(𝜉,𝑡), is a function of 𝜉 and 𝑡 alone. Therefore, the only nontrivial Navier-Stokes equation is the radial equation. Using (2.3), one can compute the derivatives of 𝑢 to show that the coefficient of 𝜈𝜕²𝑢𝜕𝜉²+1𝜉𝜕𝑢−𝑢𝜕𝜉𝜉²=0. (2.8)

Equation (2.7) then simplifies to𝜕𝑢𝜕𝑡+𝑢𝜕𝑢1𝜕𝜉=−𝜌_𝑓𝜕𝑝.𝜕𝜉 (2.9)

We now substitute (2.6) into (2.9) and integrate both sides with respect to 𝜉 for 𝜉∈[𝑟,𝑁], which yields an equation for 𝑃_𝑜, the pressure exerted on the outer surface of the membrane by the cerebrospinal fluid:𝑃_𝑜−𝑝_𝑁=𝜌_𝑓^𝑁_𝑟𝑅²𝜉𝜕𝜆𝜕𝑡²+𝑅²𝜆𝜉𝜕²𝜆𝜕𝑡²−𝑅⁴𝜆²𝜉³𝜕𝜆𝜕𝑡²𝑑𝜉=𝜌_𝑓𝑅²𝑁ln𝑟𝜕𝜆𝜕𝑡²+𝑅²𝑁𝜆ln𝑟𝜕²𝜆𝜕𝑡²+𝑅⁴𝜆²2𝑁²𝜕𝜆𝜕𝑡²−𝑅⁴𝜆²2𝑟²𝜕𝜆𝜕𝑡²=𝜌_𝑓𝑅²𝑁𝜆ln𝜕𝑅𝜆²𝜆𝜕𝑡²+𝑅²𝑁ln+𝑅𝑅𝜆⁴𝜆²2𝑁²−𝑅²2𝜕𝜆𝜕𝑡². (2.10)

Here the pressure 𝑝_𝑁 represents the pressure far enough from the membrane that 𝑢(𝑁,𝑡)≪𝑢(𝑟,𝑡) for all 𝑡.

Next we use (2.10) for 𝑃_𝑜 to compute the radial stress boundary condition on the outer wall. We use the well-known constitutive relation for a Newtonian fluid [1]𝐭=−𝑃𝐈+2𝜇𝐃, (2.11) where t is the Cauchy stress tensor that is expressed in terms of the pressure at the outer wall boundary 𝑟_𝑜, 𝐃=1/2[∇𝐯+(∇𝐯)^𝑇], the stretching tensor, and 𝜇, the viscosity of the CSF. For a detailed description of tensors, see [14]. In physical components, (2.11) becomes𝑡_𝜉𝜉𝑟_𝑜=−𝑃_𝑜+2𝜇𝜕𝑢().𝜕𝜉𝑟,𝑡 (2.12)

Substituting (2.6) in (2.11) yields𝑡_𝜉𝜉𝑟_𝑜=−𝑃_𝑜𝑅−2𝜇²𝜆𝑟²𝜕𝜆,𝜕𝑡 (2.13)

so that𝑡_𝜉𝜉𝑟_𝑜=−𝑃_𝑜−2𝜇𝜆𝜕𝜆.𝜕𝑡 (2.14)

3. Model for the Blood

To represent the blood pressure on the inner wall 𝑃_𝑖(𝑡) we use a (finite) Fourier series expansion, with mean arterial pressure 𝑃_𝑚, and circular frequency 𝜔;𝑃_𝑖(𝑡)=𝑃_𝑚+¹⁰_𝑛=1𝐴_𝑛cos(𝑛𝜔𝑡)+𝐵_𝑛sin(𝑛𝜔𝑡). (3.1) Here 𝐴_𝑛, 𝐵_𝑛, 𝑃_𝑚, and 𝜔 are experimentally determined constants which are well-established in the literature [15, 16].

4. Model for the Wall

We now use the equation for the motion of a membrane𝜅₁𝑇₁+𝜅₂𝑇₂+𝑞_𝑛=−𝜌ℎ𝑅𝜆̈, (4.1) where𝑞_𝑛=−𝑃_𝑖−𝑡_𝜉𝜉𝑟_𝑜=𝑡_𝜉𝜉𝑟_𝑖−𝑡_𝜉𝜉𝑟_𝑜, (4.2)

𝜌 is the density of the membrane, and 𝑇_𝑖 are the principal stress resultants for 𝑖=1,2. Since the principal curvatures of a cylindrical membrane are 𝜅₁=0, 𝜅₂=1/𝑟 where 𝑟 is the radius of the cylinder, and since ℎ=𝐻/𝜆² for an incompressible material [1], we have𝜌ℎ𝑅𝜆̈=−𝑇𝑟−𝑞_𝑛𝑇=−𝑅𝜆−𝑡_𝜉𝜉𝑟_𝑖+𝑡_𝜉𝜉𝑟_𝑜. (4.3) Using (2.14) and substituting (4.3) into (4.1) yields the equation of motion for the membrane model of cylindrical inflation of an aneurysm as 𝜌𝐻𝑅𝜆²𝜕²𝜆𝜕𝑡²=−𝜌_𝑓𝑅²𝑁𝜆ln𝜕𝑅𝜆²𝜆𝜕𝑡²+𝑅²𝑁ln+𝑅𝑅𝜆⁴𝜆²2𝑁²−𝑅²2𝜕𝜆𝜕𝑡²−𝑝_𝑁+𝑃_𝑖−2𝜇𝜆𝜕𝜆−𝑇(𝜆)𝜕𝑡,𝑅𝜆 (4.4) which can be rewritten in the form𝜌𝐻𝑅𝜆²+𝜌_𝑓𝑅²𝑁𝜆ln𝜕𝑅𝜆²𝜆𝜕𝑡²+𝜌_𝑓𝑅²𝑁ln+𝑅𝑅𝜆⁴𝜆²2𝑁²−12𝑅²𝜕𝜆𝜕𝑡²+2𝜇𝜆𝜕𝜆+𝑇(𝜆)𝜕𝑡𝑅𝜆=𝑃_𝑖−𝑝_𝑁. (4.5)

4.1. Hyperelasticity

In order to analyze this model, we need to consider a specific constitutive relation for the membrane. The two-dimensional deformation gradient which maps the undeformed coordinates (𝑅,Θ) to the deformed coordinates (𝑟,𝜃) is given by𝐅=𝜕𝑟1𝜕𝑅𝑅𝜕𝑟𝑟𝜕Θ𝜕𝜃𝑟𝜕𝑅𝑅𝜕𝜃=𝜕Θ𝜆00𝜆. (4.6) Since we employ a membrane approximation, one can postulate the existence of a two-dimensional strain-energy function 𝑤. The constitutive relation for the membrane can be expressed as [1]𝑇_𝛼𝛽=1𝐽𝐹_𝛼𝑎𝐹_𝛽𝑏𝜕𝑤𝜕𝐸_𝑎𝑏,𝛼,𝛽,𝑎,𝑏=1,2, (4.7) where 𝐓 is the two-dimensional Cauchy stress resultant tensor and 𝐽=det𝐅. We consider a Fung-type strain-energy function 𝑤 of the form [1, 17]𝑤=𝑐𝑒^𝑄−1,𝑄=𝑐₁𝐸²₁₁+𝑐₂𝐸²₂₂+2𝑐₃𝐸₁₁𝐸₂₂, (4.8) where 𝑐 and 𝑐_𝑖 are material parameters and the Green strain is 𝐄=1/2(𝐂−𝐈) with the right Cauchy-Green Strain given by 𝐂=𝐅^𝑇𝐅. Substituting Γ≡𝑐₁+𝑐₃, (4.8), and (4.6) into (4.7) yields𝑇(𝜆)=𝑐Γ𝑒^𝑄𝜆²1−1,𝑄=2Γ𝜆²−1². (4.9) Equation (4.5) can be nondimensionalized using the following scales for length (𝐿_𝑠), time (𝑇_𝑠), and mass (𝑀_𝑠):𝐿_𝑠=𝑅,𝑇_𝑠=𝜌𝑅²𝐻𝑐,𝑀_𝑠=𝜌𝑅²𝐻. (4.10) Using these scales reduces (4.5) to1𝑥²+𝑏𝑥𝛽−ln𝑥𝑥̈+𝑏𝛽−ln𝑥+𝑏𝑎𝑥²2−12𝑏𝑥̇²𝑥̇+2𝑚𝑥+𝑓(𝑥)𝑥=𝐹(𝜏), (4.11) where superimposed dots imply differentiation with respect to nondimensionalized time 𝜏 and𝜌𝑥≡𝜆,𝑏≡^𝑓𝑅𝑁𝜌𝐻,𝛽≡ln𝑅𝜇𝑚≡,𝑇𝜌𝑐𝐻𝑓≡𝑐𝑅,𝐹≡𝑐𝑃_𝑖−𝑝_𝑁𝑡√,𝜏≡𝑐𝜌𝑅²𝐻𝑅,𝑎≡²𝑁². (4.12)(4.13)

4.2. Viscoelasticity

In this section we consider the effects of a short-term viscoelastic wall response on the stability of the aneurysm. We will describe the mathematical models for these components and derive a radial inflation model for a cylindrical aneurysm analogous to a viscoelastic spherical inflation model derived by David and Humphrey [11]. Using the velocity gradient 𝐋=𝐅̇𝐅⁻¹, the 2D stretching tensor 1𝐃=2𝐋+𝐋^𝑇=1𝜆𝑑𝜆001𝑑𝑡𝜆𝑑𝜆𝑑𝑡, (4.14) and the viscohyperelastic response is given by1𝐓=ℎ𝐭=𝐽𝐅⋅𝜕𝑤𝜕𝐄⋅𝐅^𝑇(𝜆)1+2𝜇ℎ𝐃=𝑐𝑓+2𝜇𝐻𝜆³𝑑𝜆,𝑑𝑡 (4.15) where 𝐓 is the Cauchy stress resultant tensor, 𝐭 is the Cauchy stress, 𝑓 is the same as in the hyperelastic model, 𝜇 is the “solid viscosity” of the aneurysm, and ℎ=𝐻/𝜆² since the material is incompressible. This yields the following nondimensionalized governing equation:1𝑥²+𝑏𝑥𝛽−ln𝑥𝑥̈+𝑏𝛽−ln𝑥+𝑏𝑎𝑥²2−12𝑏𝑥̇²+2𝑚_𝑓𝑥̇𝑥+𝑓(𝑥)𝑥𝑥̇+2𝑚𝑑𝑥⁴=𝐹(𝜏), (4.16) where 𝜇_𝑓 is the viscosity of the CSF, 𝑑=𝐻/𝐴, 𝑚_𝑓=𝜇_𝑓/𝜌𝑐ℎ, 𝑚=𝜇/𝜌𝑐ℎ, and the other parameters are the same as in the hyperelastic case.

5. Stability Analysis

5.1. Hyperelastic Model

In order to analyze (4.11), we transform it into a system of first-order equations using a change in variable. Let 𝑦₀≡𝑥 and 𝑦₁≡𝑑𝑥/𝑑𝜏. Then the final system of equations becomes 𝑑𝑦₀𝑑𝜏=𝑦₁,𝑑𝑦₁=𝐹(𝜏)𝑑𝜏+0.5𝑏𝑦²₁−2𝑚𝑦₁/𝑦₀−𝑓𝑦₀/𝑦₀−𝑏𝛽𝑦²₁+𝑏ln𝑦₀𝑦²₁−𝑏𝑎𝑦²₀𝑦²₁/2𝑦₀⁻²+𝑏𝑦₀𝛽−ln𝑦₀. (5.1) We consider an autonomous problem, that is, 𝐹(𝜏)≡𝐹₀ for some constant 𝐹₀. Then we have𝑓(𝑥)=Γ𝑒^𝑄𝑥²1−1,𝑄=2Γ𝑥²−1². (5.2) The derivative can be computed as𝜕𝑓𝜕𝑦_0(𝛼,0)=2Γ𝛼𝑒^𝑄Γ𝛼²−1²+1. (5.3) We now expand the first-order equations in Taylor series around the equilibrium point (𝛼,0) to get 𝑦̇₀𝑦₀,𝑦₁=𝐺(𝛼,0)+𝜕𝐺𝜕𝑦_0(𝛼,0)𝑦₀−𝛼+𝜕𝐺𝜕𝑦_1(𝛼,0)𝑦₁𝑦̇−0+HOT,₁𝑦₀,𝑦₁=𝐻(𝛼,0)+𝜕𝐻𝜕𝑦_0(𝛼,0)𝑦₀−𝛼+𝜕𝐻𝜕𝑦_1(𝛼,0)𝑦₁−0+HOT, (5.4) where 𝐺 and 𝐻 represent the functions in the system (5.1) and HOT denotes higher-order terms, which we will neglect. Since 𝐺(𝛼,0)=𝐻(𝛼,0)=0 the final linearized system of equations can be written as𝑦̇₀𝑦̇₁=𝐹01₀−𝜕𝑓/𝜕𝑦₀|_(𝛼,0)1/𝛼+𝑏𝛼²𝛽−ln𝛼−2𝑚1/𝛼+𝑏𝛼²𝑦𝛽−ln𝛼⁰𝑦−𝛼₁, (5.5) where 𝐹₀ is the constant pressure at the equilibrium (𝛼,0); hence 𝐹₀=𝑓(𝛼)/𝛼. We denote the Jacobian matrix 𝐉. Let 𝜏_𝑟=tr𝐉 and let 𝛿=det𝐉. Now 𝜏_𝑟<0 if and only if the CSF is viscous (i.e., 𝑚>0). Next, we will prove stability results for inviscid and viscous CSF.

Lemma 5.1. Let ℍ={(𝑥,𝑦)∈ℝ²∶𝑥≥1}. Also, let 𝑏>0 and suppose 1/𝛼+𝑏𝛼²(𝛽−ln𝛼)≠0. Then the nondimensionalized system (4.11) and (4.16) are 𝐶²(ℍ).

Proof . Since 𝑏>0, 𝑥≥1, and 1/𝛼+𝑏𝛼²(𝛽−ln𝛼)≠0forall𝑥∈ℍ, the first and second derivatives of (4.11) exist and are continuous on ℍ.

Lemma 5.2. Suppose 𝑚>0. Then the point (𝛼,0) is a hyperbolic equilibrium point.

Proof. Since the trace of the linearization 𝜏_𝑟<0, and 𝛿≠0, all of the eigenvalues of the linearized system (5.5) have nonzero real part.

For the inviscid case we have the following theorem.

Theorem 5.3. Let 𝑚=0 and suppose (𝛼,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose 𝑁≥𝑅𝛼, so that 𝛽≥ln𝛼. Then for the system (5.1) one has the following.
(1) If (Γ𝑒^𝑄(2𝛼Γ(𝛼²−1)²+𝛼+𝛼⁻¹))/(1/𝛼+𝑏𝛼²(𝛽−ln𝛼))>0, the nonlinear system has a center at (𝛼,0).
(2) If (Γ𝑒^𝑄(2𝛼Γ(𝛼²−1)²+𝛼+𝛼⁻¹))/(1/𝛼+𝑏𝛼²(𝛽−ln𝛼))<0, the nonlinear system has a saddle point at (𝛼,0).

Proof. Since 𝑚=0, 𝜏_𝑟=0. Thus the eigenvalues 𝜏Λ=^𝑟±𝜏²_𝑟−4𝛿2=±√−𝛿=±−Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹1/𝛼+𝑏𝛼².𝛽−ln𝛼 (5.6)
(1) Since Γ𝑒^𝑄(2𝛼Γ(𝛼²−1)²+𝛼+𝛼⁻¹)/(1/𝛼+𝑏𝛼²(𝛽−ln𝛼))>0, we have 𝛿>0, and the linearized system has a pair of purely imaginary complex conjugate eigenvalues and (𝛼,0) is a center. Since the nonlinear system is invariant under the transformation (𝜏,𝑦₁)→(−𝜏,−𝑦₁), it is symmetric with respect to the 𝑥-axis. Now, the nonlinear system is 𝐶¹(ℍ) by Lemma 5.1. Therefore, the point (𝛼,0) must also be a center for the nonlinear system.
(2) Since Γ𝑒^𝑄(2𝛼Γ(𝛼²−1)²+𝛼+𝛼⁻¹)/(1/𝛼+𝑏𝛼²(𝛽−ln𝛼))<0, we have 𝛿<0, and the system has two real eigenvalues of opposite sign. This shows that (𝛼,0) is a saddle point for the linearized system. Now, the nonlinear system is 𝐶¹(ℍ) by Lemma 5.1. Also, (𝛼,0) is a hyperbolic equilibrium point since the eigenvalues of its linearization both have nonzero real part. Therefore, the point (𝛼,0) must also be a saddle point for the nonlinear system.

Also, for a viscous CSF we have the following theorem.

Theorem 5.4. Let the CSF be viscous, that is, 𝑚>0. Also let Γ>0 and suppose (𝛼,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose 𝑁≥𝑅𝛼, so that 𝛽≥ln𝛼 and that 2Γ𝛼−(1/𝛼)>0. Then for the nondimensionalized system (5.1), there are three possible cases for the stability of the equilibrium point (𝛼,0).
(1)If 𝑚²=Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.7) then (𝛼,0) is an asymptotically stable proper node. (2)If 𝑚²<Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.8) then (𝛼,0) is an asymptotically stable focus. (3)If 𝑚²>Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.9) then (𝛼,0) is an asymptotically stable node.

Proof. The eigenvalues of the associated linearized system (see (5.5)) are 𝜏Λ=^𝑟±𝜏²_𝑟−4𝛿2=−2𝑚/𝔄±(−2𝑚/𝔄)²−4Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹/𝔄2, (5.10) where 𝔄 denotes (𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼). First, note that since 𝑚,𝑏>0,𝛽≥ln𝛼,𝜏_𝑟 is negative for all three cases.
(1)If 𝑚²=Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.11) then the linear system has the repeated real eigenvalue 𝑚Λ=−𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼<0, (5.12) and the linearized system has an asymptotically stable proper node at (𝛼,0). Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1. Also, (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.2. Therefore the point (𝛼,0) must be an asymptotically stable proper node for the nonlinear system.(2)If 𝑚²<Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.13) then the linear system has a pair of complex conjugate eigenvalues with nonzero real part. Thus (𝛼,0) is an asymptotically stable focus for the linear system. Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1. Also, (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.2. Therefore the point (𝛼,0) must be an asymptotically stable focus for the nonlinear system.(3)If 𝑚²>Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (5.14) then the linear system has a pair of real eigenvalues. Since Γ>0, 𝑚>𝑚²−Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼 (5.15) and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (𝛼,0). Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1. Also, (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.2. Therefore the point (𝛼,0) must be an asymptotically stable node for the nonlinear system.

Theorem 5.5. Suppose the CSF is viscous, that is, 𝑚>0. Also suppose F(𝜏)<𝑓(𝑥)/𝑥 and 𝑁≥𝑅𝛼, so that 𝛽≥ln𝛼. Then the equilibrium point (𝛼,0) is asymptotically stable.

Proof. Using a result of Freitas's [18], 𝐿(𝑥,𝑥̇)=2^𝑥₀𝑐(𝑟)𝐵(𝑟)𝑑𝑟+𝐵(𝑥)(𝑥̇)², (5.16) with 𝐵(𝑥)=𝑒^{2∫𝑥0𝑏(𝑟)𝑑𝑟}, (5.17) is a Lyapunov function for systems of the form 𝑥̈+[𝑎(𝑥)+𝑏(𝑥)𝑥̇]𝑥̇+𝑐(𝑥)=0, (5.18) where 𝑎,𝑏,𝑐 satisfy the conditions
(1) 𝑎,𝑏,𝑐∶ℝ→ℝ and are continuous,
(2) 𝑎(𝑥)>0 except possibly at points where 𝑐(𝑥)=0,
(3) 𝑥𝑐(𝑥)>0 outside a bounded interval.
Applying this to the nondimensionalized governing equation (4.11) with 𝑎(𝑥)=2𝑚/𝑥1/𝑥²,𝑏(𝑥)=𝑏+𝑏𝑥𝛽−ln𝑥𝛽−ln𝑥+𝑏𝑎𝑥²/2−(1/2)𝑏1/𝑥²,𝑐(𝑥)=𝑓(𝑥)+𝑏𝑥𝛽−ln𝑥/𝑥1/𝑥²−𝐹(𝜏)+𝑏𝑥𝛽−ln𝑥1/𝑥²,+𝑏𝑥𝛽−ln𝑥 (5.19) yields 𝐿(𝛼,0)=0 and 𝐿̇(𝑥,𝑥̇)<0 on an open region around the equilibrium and the result follows.

Corollary 5.6. Suppose the CSF is viscous, that is, 𝑚>0. Also suppose 𝐹(𝜏)≡0 and 𝛽≥0. Then the equilibrium point (1,0) is asymptotically stable.

Proof. Since 𝑓(𝑥)/𝑥>0 on {(𝑥,𝑦)∈ℝ²∶𝑥>0&(𝑥,𝑦)≠(1,0)}, applying Theorem 5.5 to the nondimensionalized governing equation (4.11) yields the desired result.

5.2. Viscoelastic Model

For the viscoelastic case, we can also rewrite the governing equation as a system of first-order ODEs, yielding𝑑𝑦₀𝑑𝜏=𝑦₁𝑑𝑦₁=𝐹(𝜏)𝑑𝜏+0.5𝑏𝑦²₁−2𝑚_𝑓𝑦₁/𝑦₀−𝑓𝑦₀/𝑦₀−𝑏𝛽𝑦²₁+𝑏ln𝑦₀𝑦²₁𝑦₀⁻²+𝑏𝑦₀𝛽−ln𝑦₀−𝑏𝑎𝑦²₀𝑦²₁/2+2𝑚𝑑𝑦₁/𝑦⁴₀𝑦₀⁻²+𝑏𝑦₀𝛽−ln𝑦₀. (5.20) Linearization of this system of equations yields𝑦̇₀𝑦̇₁=𝐹01₀−𝜕𝑓/𝜕𝑦_0(𝛼,0)𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼−2𝑚_𝑓𝛼³−2𝑚𝑑𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵𝑦ln𝛼⁰𝑦−𝛼₁. (5.21) Here the Jacobian matrix 𝐉 has trace 𝜏_𝑟=tr𝐉=(−2𝑚_𝑓𝛼³−2𝑚𝑑)/(𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼) and determinant 𝛿=det𝐉=(−𝐹₀+(𝜕𝑓/𝜕𝑦₀)|_(𝛼,0))/(𝛼⁻¹+𝑏𝛼²(𝛽−𝑏ln𝛼)).

Here we must have 𝐹₀=𝑓(𝛼)/𝛼 at the equilibrium (𝛼,0) and using 𝑓(𝑥)=Γ𝑒^𝑄(𝑥²−1), we have 𝛿=Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼>0, (5.22) and we can analyze the stability of the linearized system.

We now have an immediate analogue of Lemma 5.2.

Lemma 5.7. Suppose 𝑚>0 or 𝑚_𝑓>0 and 𝛽≥ln𝛼. Then the point (𝛼,0) is a hyperbolic equilibrium point.

Proof. Since the trace of the linearization 𝜏_𝑟<0, and 𝛿≠0, all of the eigenvalues of the linearized system (5.21) have nonzero real part.

Also, note that the trace 𝜏_𝑟<0 if either the CSF or the membrane is viscous. For a viscous CSF and viscoelastic wall we have the following theorem.

Theorem 5.8. Let the CSF be viscous and the wall be viscoelastic, that is, 𝑚,𝑚_𝑓>0. Also let Γ>0 and suppose (𝛼,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.20). Finally, suppose 𝛽≥ln𝛼. Then for the nondimensionalized system (5.20), there are three possible cases for the stability of the equilibrium point (𝛼,0).
(1)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³=𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.23) then (𝛼,0) is an asymptotically stable proper node. (2)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³<𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.24) then (𝛼,0) is an asymptotically stable focus. (3)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³>𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.25) then (𝛼,0) is an asymptotically stable node.

Proof. The eigenvalues of the associated linearized system (see (5.21)) are 𝜏Λ=^𝑟±𝜏²_𝑟−4𝛿2=−𝑚_𝑓𝛼³−𝑚𝑑𝔅±−2𝑚_𝑓𝛼³−2𝑚𝑑/𝔅²−4−𝐹₀+𝜕𝑓/𝜕𝑦₀|_(𝛼,0)/𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼2, (5.26) where 𝔅 denotes (𝛼²−𝑏(𝛽𝛼⁵−𝛼⁵ln𝛼)). First, note that since 𝑚,𝑚_𝑓>0,𝜏_𝑟 is negative for all three cases.
(1)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³=𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.27) then the linear system has repeated negative real eigenvalues and the linearized system has an asymptotically stable proper node at (𝛼,0). Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1 and (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.7. Therefore the point (𝛼,0) must also be an asymptotically stable proper node for the nonlinear system.(2)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³<𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.28) then the linear system has a pair of complex conjugate eigenvalues with nonzero real part. Thus (𝛼,0) is an asymptotically stable focus for the linear system. Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1 and (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.7. Therefore the point (𝛼,0) must also be an asymptotically stable focus for the nonlinear system.(3)If 𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³>𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹, (5.29)then the linear system has a pair of real eigenvalues. Since Γ>0,𝑚_𝑓𝛼³>+𝑚𝑑𝑚²_𝑓𝛼⁶+𝑚²𝑑²+2𝑚_𝑓𝑚𝑑𝛼³−𝛼²+𝑏𝛽𝛼⁵−𝑏𝛼⁵ln𝛼²𝛼⁻¹+𝑏𝛼²𝛽−𝑏ln𝛼Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹ (5.30) and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (𝛼,0). Now, the nonlinear system is 𝐶²(ℍ) by Lemma 5.1 and (𝛼,0) is a hyperbolic equilibrium point by Lemma 5.7. Therefore the point (𝛼,0) must also be an asymptotically stable node for the nonlinear system.

6. Numerical Findings

Next we solve the systems (5.1) and (5.20) numerically for the following parameter values, taken from [1, 10] with the newer values for 𝑐 and Γ taken from [19]: 𝜌=1050 kg/m³, 𝑅=3×10⁻³ m, 𝐻=27.8×10⁻⁶ m, 𝜌_𝑓=1000 kg/m³, 𝑝_𝑁=399.9375 N/m²,𝑐=10.18 N/m, Γ=23.6, 𝜇_𝑓=1.26×10⁻⁴ Ns/m², 𝑁=10⁵m, and the value for the “solid viscosity” 𝜇=331.827 Ns/m² The mean blood pressure 𝑃_𝑚=90 mmHg, the Fourier coefficients for the blood pressure were taken from [11], and all values were nondimensionalized. The systems (5.1) and (5.20) are solved numerically using initial conditions such that 𝜆=1 and 𝜆̇=0. We then considered the effects of small perturbations of the initial conditions on the behavior of the two components of the solution, that is, the stretch ratio 𝜆 and the stretch rate 𝜆̇. We were primarily interested in the behavior of the stretch ratio over time for different parameter values, but also looked at the behavior of the system in two-dimensional phase space.

Figure 2 demonstrates that greatly increasing the distance 𝑁 to the outer boundary of the exterior fluid region leads to a decrease in the rate at which the solution approaches the equilibrium. This is the case for both the hyperelastic model and the viscoelastic model. However, the qualitative behavior remains the same, as the solution approaches an asymptotically stable equilibrium point. As can be seen in Table 1, the eigenvalues of the linearizations demonstrate that the equilibrium is an asymptotically stable focus in each case. This provides justification for allowing the value of 𝑁 to vary, as long as 𝑅≪𝑁. The physical meaning of the parameter 𝑁 is a distance far from the membrane, at which the pressure due to the surrounding material is much less than the pressures exerted on the membrane by the blood and the CSF on the membrane boundaries. Note that 𝑁 was taken to be an infinite quantity in spherical models, for example, see [3, 10, 11]. This justifies the use of a large value for 𝑁.

Table 1 

Eigenvalues for models with varying 𝑁.

(a) Viscoelastic model, 𝑁=0.1m

(b) Viscoelastic model, 𝑁=1m

(c) Viscoelastic model, 𝑁=103m

(d) Viscoelastic model, 𝑁=105m

Figure 2 

Stretch ratio with varying 𝑁.

For the values of 𝑁 chosen in the following numerical experiments, the nondimensional parameter 𝑎 is very small. This causes the nonlinear term involving 𝑦²₀𝑦²₁ to have only a small effect on the solution. For example, for 𝑁=10⁵m, as used to generate the solution plots, we have 𝑎=9×10⁻¹⁶, which makes the contribution due to this term nearly negligible.

We considered an initial disturbance at 𝜏=0 by allowing (𝑦₀,𝑦₁) to vary between (1.01,0) and (1.02,0). Using these values, and 𝑁=1m, appropriate phase portraits were generated to validate the analytical results proved in the previous section. The governing equations were solved numerically using Runge-Kutta methods available via the X-Windows Phase Plane (XPP) software package. Figure 3 shows the phase portraits for the inviscid and viscous cases for the CSF, as well as for the viscoelastic model. Figure 3 was generated using 𝜔 = 2𝜋 Hz, as were the other figures unless otherwise indicated.

(a) Inviscid CSF

(b) Viscous CSF

(c) Viscous CSF, m𝑓=10

(d) Viscoelastic model

Figure 3 

Phase portraits for inviscid CSF, viscous CSF, and viscoelastic model.

The numerically generated phase portraits for the inviscid and viscoelastic cases (Figures 3(a) and 3(d), resp.) appear to show a center and a stable focus, consistent with the analytical results in the stability analysis section. The figure showing the phase portrait for the hyperelastic model with viscous CSF (see Figure 3(b)) seems to show a center; however, this is due to the fact that the real parts of the eigenvalues for this case are approximately −1×10⁻⁶, which is very close to zero. Thus the solution converges to the equilibrium very slowly relative to the other cases. Using the value 𝑚_𝑓=10, we were able to generate a phase portrait (see Figure 3(c)) which better represents the theoretical dynamical behavior of the hyperelastic model.

We also investigated the behavior for Γ=−10 to try and locate a saddle point for the inviscid case, as should be expected from the analytical results. The equilibrium point is now located at 𝑦₀≈0.0240839, 𝑦₁=0, which is not physically relevant since 𝑦₀<1. The eigenvalues of the Jacobian matrix can be computed as Λ₁≈−0.255257, Λ₂≈0.255259, indicating that the linearized system indeed has a saddle at the equilibrium.

Table 2 shows the different possibilities for the eigenvalues of the Jacobian matrix for the hyperelastic model for various values of 𝑚_𝑓 with 𝑁=10⁵ m. The theoretical value given is computed from the condition 𝑚²=Γ𝑒^𝑄2𝛼Γ𝛼²−1²+𝛼+𝛼⁻¹𝛼⁻¹+𝑏𝛽𝛼²−𝑏𝛼²ln𝛼, (6.1) which gives the theoretical cutoff point for the equilibrium to be a node instead of a focus. As demonstrated in the table, the value of 𝑚_𝑓 computed from the numerics beyond which the equilibrium becomes a node is close to the theoretical value, lying between 𝑚_𝑓=373 and 𝑚_𝑓=373.666.

Table 2 

Eigenvalues for hyperelastic model with varying 𝑚_𝑓.

Plots of the stretch ratio 𝜆 are shown in Figures 5, 6, and 7 for various values of 𝑚, 𝑅, and 𝜔. The actual equilibrium point is at approximately (𝑦₀,𝑦₁)=(1.0531,0); however, we examine stability by considering a perturbed system with initial condition (𝑦₀,𝑦₁)=(1.05,0). Figure 4(a) shows the neutrally stable behavior of the solution for the inviscid case while Figures 4(b)–4(d) show the decay of the solution to an equilibrium state for the viscous and viscoelastic cases. This supports the conclusion that the equilibrium is asymptotically stable if either viscosity is nonzero.

(a) Inviscid CSF

(b) Viscous CSF

(c) Viscoelastic model

(d) Viscoelastic model, m𝑓=10

Figure 4 

Stretch ratio for inviscid CSF, viscous CSF, and viscoelastic model.

(a) m = 1

(b) m = 10

(c) m = 100

(d) 𝑚=4𝑚𝑓

Figure 5 

Stretch ratio for viscoelastic model with varying 𝑚.

(a) 𝑅=0.003m

(b) 𝑅=0.005m

(c) 𝑅=0.03m

(d) 𝑅=0.05m

Figure 6 

Stretch ratio for viscoelastic model with varying radii.

(a) 𝜔=4𝜋Hz

(b) 𝜔=10𝜋Hz

(c) 𝜔=100Hz

(d) 𝜔=400Hz

Figure 7 

Stretch ratio with varying 𝜔.

Figure 5 shows the effect of varying the nondimensional parameter 𝑚 on the stretch ratio 𝜆. Larger values of 𝑚 cause the solution to approach the equilibrium state more rapidly. The value 𝑚=4𝑚_𝑓 was used in [11] to illustrate the effect of a viscoelastic wall on the stability of the solution. However, the much larger values of 𝑚 consistent with the wall viscosity used in [3] and seen in Figure 4(c) lead to much faster decay to the equilibrium state. This larger viscosity value is likely more in line with the biological structure of the arterial wall.

The effect of varying the initial radius 𝑅 in the viscoelastic model is seen in Figure 6. Larger values of 𝑅 lead to slower decay to the equilibrium state. The data used in [10] is for saccular aneurysms, which are represented by a spherical geometry. For a cylindrical geometry representing a fusiform aneurysm, it may be more appropriate to consider a larger initial radius. Figures 6(c)-6(d) show the effects of a more realistic aneurysm size.

Figure 7 demonstrates that greatly increasing the circular frequency 𝜔 leads to a slight increase in the rate at which the solution approaches the equilibrium. However, it seems that the frequency is not a major factor in the rate of decay; changing the viscosities has a much larger effect.

The hyperelastic and viscoelastic aneurysm models derived in this paper are analogous to the models derived for a spherical geometry in [10, 11], respectively. The actual equilibrium point for the spherical case is at approximately (𝑦₀,𝑦₁)=(1.029,0); however, we follow [11] and examine stability by considering a perturbed system with initial condition (𝑦₀,𝑦₁)=(1.031,0). The comparison of these models for the two geometries is in Figure 8. In both cases, the magnitude of the solution and the size of the oscillations are larger in the cylindrical case than in the spherical case. The eigenvalues of the linearization for a cylinder and for a sphere are given in Table 3. The eigenvalues are complex and have negative real parts for each model, supporting the conclusion that the equilibrium is an asymptotically stable focus for each case.

Table 3 

Eigenvalues for spherical and cylindrical membrane models.

(a) Viscous CSF, cylinder

(b) Viscous CSF, sphere

(c) Viscoelastic model, cylinder

(d) Viscoelastic model, sphere

Figure 8 

Stretch ratio for cylindrical versus spherical geometry.

7. Discussion and Conclusions

Shah and Humphrey [10] were the first to consider the nonlinear elastic behavior of saccular aneurysms over finite strains as well as the effects of a viscous CSF on the coupled fluid-structure nature of such lesions for a spherical geometry. In this paper, we have extended results from [10, 11] for a spherical geometry to a cylindrical one. This may allow the present results to be applied to fusiform as well as saccular aneurysms while including the short-term effects of viscoelasticity. This is especially important in view of the prevalence of fusiform aneurysms in the abdominal aorta. The simplifications of the present model allowed us to investigate the elastodynamics of a subclass of aneurysms. However, aneurysms come in a large variety of shapes and sizes, and the need for more general and realistic models is evident.

We have also performed stability analysis of the hyperelastic and viscoelastic models using both numerical experiments and analytical techniques. The present results support the earlier conclusion [10, 11] that aneurysms are stabilized both by a viscous CSF and by viscoelasticity in the wall while under the dynamical loads of the cardiac cycle. This seems to indicate that the pulsatile blood pressure is not a major cause of aneurysm enlargement and rupture. The numerical results also demonstrate that the oscillations and wall stretch are of greater magnitude in the cylindrical case than they were for a sphere. This indicates that a cylindrically shaped (possibly fusiform) aneurysm may expand farther than a spherical saccular aneurysm with similar material properties. However, both subclasses of aneurysms are stable under the loads of the cardiac cycle.

Future research could involve extending these results using the full 3D elasticity theory in place of the membrane approach as was done for a spherical geometry in [3]. The need for better estimation of material parameters is also evident, especially in view of the differences in size and shape between abdominal aortic and intracranial aneurysms. For example, the Fourier coefficients for blood pressure used in these models were based on a study of anesthetized patients, thus, the values lead to a lower blood pressure than would be expected for a fully conscious subject. There is also the need for consideration of much longer time scales, as aneurysm growth and rupture usually takes place over a period of years rather than seconds as considered in this work. The material parameters, blood pressure, and longer time scales will be considered in a forthcoming paper.

Acknowledgment

This work is supported in part by National Institutes for Health grant R21 EB008804.