# Stability of Membrane Elastodynamics with Applications to Cylindrical Aneurysms

Authors: A. Samuelson; P. Seshaiyer
Journal: Journal of Applied Mathematics
 (2011)

---

## 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.


---

## Body

## 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 arelateral 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 wordmembrane 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 Figure1. 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].Membrane cross-sections.
(a)(b)Let the stretch ratio of the aneurysmal wall be denotedλ=r/R where r and R are the deformed and undeformed radii of the aneurysm. Also let h and H 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 t and radial position r, and axial velocity is constant, that is, u=u(r,t) and w≡w0, where u is the radial velocity, v is the circumferential velocity, w is the axial velocity, and the velocity vector is v=(u,v,w). 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 ∇·v=0, which in cylindrical coordinates (r,θ,z) is(2.1)1ξ∂∂ξ(ξu)+1ξ∂v∂θ+∂w∂z=0.Assuming that the axial velocity is independent ofz and the circumferential velocity is zero, we then have v=0 and ∂w/∂z=0. Equation (2.1) becomes(2.2)1ξ∂∂ξ(ξu)=0,which yields(2.3)u=g(t)ξ,
where g(t) 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)∂∂t(r-R)=R∂λ∂t=g(t)r.Solving forg(t), we have(2.5)g(t)=R2λ∂λ∂t.Combining (2.3) and (2.5) gives(2.6)u(ξ,t)=g(t)ξ=R2λξ∂λ∂t.The radial velocityu satisfies the Navier-Stokes radial equation given by(2.7)∂u∂t+u∂u∂ξ=-1ρf∂p∂ξ+ν(∂2u∂ξ2+1ξ∂u∂ξ-uξ2),
where p is the pressure and ρf and ν are the density and the kinematic viscosity, respectively, of the surrounding fluid. Here the fluid pressure P=P(ξ,t), is a function of ξ and t alone. Therefore, the only nontrivial Navier-Stokes equation is the radial equation. Using (2.3), one can compute the derivatives of u to show that the coefficient of ν(2.8)∂2u∂ξ2+1ξ∂u∂ξ-uξ2=0.Equation (2.7) then simplifies to(2.9)∂u∂t+u∂u∂ξ=-1ρf∂p∂ξ.We now substitute (2.6) into (2.9) and integrate both sides with respect to ξ for ξ∈[r,N], which yields an equation for Po, the pressure exerted on the outer surface of the membrane by the cerebrospinal fluid:(2.10)Po-pN=ρf∫rN[R2ξ(∂λ∂t)2+R2λξ∂2λ∂t2-R4λ2ξ3(∂λ∂t)2]dξ=ρf(R2ln⁡(Nr)(∂λ∂t)2+R2λln⁡(Nr)∂2λ∂t2+R4λ22N2(∂λ∂t)2-R4λ22r2(∂λ∂t)2)=ρf(R2λln⁡(NRλ)∂2λ∂t2+(R2ln⁡(NRλ)+R4λ22N2-R22)(∂λ∂t)2).Here the pressurepN represents the pressure far enough from the membrane that u(N,t)≪u(r,t) for all t.Next we use (2.10) for Po to compute the radial stress boundary condition on the outer wall. We use the well-known constitutive relation for a Newtonian fluid [1](2.11)t=-PI+2μD,
where t is the Cauchy stress tensor that is expressed in terms of the pressure at the outer wall boundary ro, D=1/2[∇v+(∇v)T], the stretching tensor, and μ, the viscosity of the CSF. For a detailed description of tensors, see [14]. In physical components, (2.11) becomes(2.12)tξξ(ro)=-Po+2μ∂u∂ξ(r,t).Substituting (2.6) in (2.11) yields(2.13)tξξ(ro)=-Po-2μR2λr2∂λ∂t,so that(2.14)tξξ(ro)=-Po-2μλ∂λ∂t.


## 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 ∇·v=0, which in cylindrical coordinates (r,θ,z) is(2.1)1ξ∂∂ξ(ξu)+1ξ∂v∂θ+∂w∂z=0.Assuming that the axial velocity is independent ofz and the circumferential velocity is zero, we then have v=0 and ∂w/∂z=0. Equation (2.1) becomes(2.2)1ξ∂∂ξ(ξu)=0,which yields(2.3)u=g(t)ξ,
where g(t) 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)∂∂t(r-R)=R∂λ∂t=g(t)r.Solving forg(t), we have(2.5)g(t)=R2λ∂λ∂t.Combining (2.3) and (2.5) gives(2.6)u(ξ,t)=g(t)ξ=R2λξ∂λ∂t.The radial velocityu satisfies the Navier-Stokes radial equation given by(2.7)∂u∂t+u∂u∂ξ=-1ρf∂p∂ξ+ν(∂2u∂ξ2+1ξ∂u∂ξ-uξ2),
where p is the pressure and ρf and ν are the density and the kinematic viscosity, respectively, of the surrounding fluid. Here the fluid pressure P=P(ξ,t), is a function of ξ and t alone. Therefore, the only nontrivial Navier-Stokes equation is the radial equation. Using (2.3), one can compute the derivatives of u to show that the coefficient of ν(2.8)∂2u∂ξ2+1ξ∂u∂ξ-uξ2=0.Equation (2.7) then simplifies to(2.9)∂u∂t+u∂u∂ξ=-1ρf∂p∂ξ.We now substitute (2.6) into (2.9) and integrate both sides with respect to ξ for ξ∈[r,N], which yields an equation for Po, the pressure exerted on the outer surface of the membrane by the cerebrospinal fluid:(2.10)Po-pN=ρf∫rN[R2ξ(∂λ∂t)2+R2λξ∂2λ∂t2-R4λ2ξ3(∂λ∂t)2]dξ=ρf(R2ln⁡(Nr)(∂λ∂t)2+R2λln⁡(Nr)∂2λ∂t2+R4λ22N2(∂λ∂t)2-R4λ22r2(∂λ∂t)2)=ρf(R2λln⁡(NRλ)∂2λ∂t2+(R2ln⁡(NRλ)+R4λ22N2-R22)(∂λ∂t)2).Here the pressurepN represents the pressure far enough from the membrane that u(N,t)≪u(r,t) for all t.Next we use (2.10) for Po to compute the radial stress boundary condition on the outer wall. We use the well-known constitutive relation for a Newtonian fluid [1](2.11)t=-PI+2μD,
where t is the Cauchy stress tensor that is expressed in terms of the pressure at the outer wall boundary ro, D=1/2[∇v+(∇v)T], the stretching tensor, and μ, the viscosity of the CSF. For a detailed description of tensors, see [14]. In physical components, (2.11) becomes(2.12)tξξ(ro)=-Po+2μ∂u∂ξ(r,t).Substituting (2.6) in (2.11) yields(2.13)tξξ(ro)=-Po-2μR2λr2∂λ∂t,so that(2.14)tξξ(ro)=-Po-2μλ∂λ∂t.


## 3. Model for the Blood

To represent theblood pressure on the inner wall Pi(t) we use a (finite) Fourier series expansion, with mean arterial pressure Pm, and circular frequency ω;(3.1)Pi(t)=Pm+∑n=110[Ancos⁡(nωt)+Bnsin⁡(nωt)].
Here An, Bn, Pm, and ω are experimentally determined constants which are well-established in the literature [15, 16].


## 4. Model for the Wall

We now use theequation for the motion of a membrane(4.1)κ1T1+κ2T2+qn=-ρhRλ̈,
where(4.2)qn=-Pi-tξξ(ro)=tξξ(ri)-tξξ(ro),ρ is the density of the membrane, and Ti are the principal stress resultants for i=1,2. Since the principal curvatures of a cylindrical membrane are κ1=0, κ2=1/r where r is the radius of the cylinder, and since h=H/λ2 for an incompressible material [1], we have(4.3)ρhRλ̈=-Tr-qn=-TRλ-tξξ(ri)+tξξ(ro).
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 (4.4)ρHRλ2∂2λ∂t2=-ρf[R2λln⁡(NRλ)∂2λ∂t2+(R2ln⁡(NRλ)+R4λ22N2-R22)(∂λ∂t)2]-pN+Pi-2μλ∂λ∂t-T(λ)Rλ,
which can be rewritten in the form(4.5)(ρHRλ2+ρfR2λln⁡(NRλ))∂2λ∂t2+ρf[R2ln⁡(NRλ)+R4λ22N2-12R2](∂λ∂t)2+2μλ(∂λ∂t)+T(λ)Rλ=Pi-pN.


### 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(R,Θ) to the deformed coordinates (r,θ) is given by(4.6)F=[∂r∂R1R∂r∂Θr∂θ∂RrR∂θ∂Θ]=[λ00λ].
Since we employ a membrane approximation, one can postulate the existence of a two-dimensional strain-energy function w. The constitutive relation for the membrane can be expressed as [1](4.7)Tαβ=1JFαaFβb∂w∂Eab,α,β,a,b=1,2,
where T is the two-dimensional Cauchy stress resultant tensor and J=detF. We consider a Fung-type strain-energy function w of the form [1, 17](4.8)w=c(eQ-1),Q=c1E112+c2E222+2c3E11E22,
where c and ci are material parameters and the Green strain is E=1/2(C-I) with the right Cauchy-Green Strain given by C=FTF. Substituting Γ≡c1+c3, (4.8), and (4.6) into (4.7) yields(4.9)T(λ)=cΓeQ(λ2-1),Q=12Γ(λ2-1)2.
Equation (4.5) can be nondimensionalized using the following scales for length (Ls), time (Ts), and mass (Ms):(4.10)Ls=R,Ts=ρR2Hc,Ms=ρR2H.
Using these scales reduces (4.5) to(4.11)(1x2+bx(β-ln⁡x))ẍ+(b(β-ln⁡x)+bax22-12b)ẋ2+2mẋx+f(x)x=F(τ),
where superimposed dots imply differentiation with respect to nondimensionalized time τ and(4.12)x≡λ,b≡ρfRρH,β≡ln⁡(NR)m≡μρcH,(4.13)f≡Tc,F≡(Rc)(Pi-pN),τ≡tcρR2H,a≡R2N2.


### 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 L=ḞF-1, the 2D stretching tensor (4.14)D=12(L+LT)=[1λdλdt001λdλdt],
and the viscohyperelastic response is given by(4.15)T=ht=1JF⋅∂w∂E⋅FT+2μhD=cf(λ)+2μH1λ3dλdt,
where T is the Cauchy stress resultant tensor, t is the Cauchy stress, f is the same as in the hyperelastic model, μ is the “solid viscosity” of the aneurysm, and h=H/λ2 since the material is incompressible. This yields the following nondimensionalized governing equation:(4.16)(1x2+bx(β-ln⁡x))ẍ+(b(β-ln⁡x)+bax22-12b)ẋ2+2mfẋx+f(x)x+2mdẋx4=F(τ),
where μf is the viscosity of the CSF, d=H/A, mf=μf/ρch, m=μ/ρch, and the other parameters are the same as in the hyperelastic case.


## 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(R,Θ) to the deformed coordinates (r,θ) is given by(4.6)F=[∂r∂R1R∂r∂Θr∂θ∂RrR∂θ∂Θ]=[λ00λ].
Since we employ a membrane approximation, one can postulate the existence of a two-dimensional strain-energy function w. The constitutive relation for the membrane can be expressed as [1](4.7)Tαβ=1JFαaFβb∂w∂Eab,α,β,a,b=1,2,
where T is the two-dimensional Cauchy stress resultant tensor and J=detF. We consider a Fung-type strain-energy function w of the form [1, 17](4.8)w=c(eQ-1),Q=c1E112+c2E222+2c3E11E22,
where c and ci are material parameters and the Green strain is E=1/2(C-I) with the right Cauchy-Green Strain given by C=FTF. Substituting Γ≡c1+c3, (4.8), and (4.6) into (4.7) yields(4.9)T(λ)=cΓeQ(λ2-1),Q=12Γ(λ2-1)2.
Equation (4.5) can be nondimensionalized using the following scales for length (Ls), time (Ts), and mass (Ms):(4.10)Ls=R,Ts=ρR2Hc,Ms=ρR2H.
Using these scales reduces (4.5) to(4.11)(1x2+bx(β-ln⁡x))ẍ+(b(β-ln⁡x)+bax22-12b)ẋ2+2mẋx+f(x)x=F(τ),
where superimposed dots imply differentiation with respect to nondimensionalized time τ and(4.12)x≡λ,b≡ρfRρH,β≡ln⁡(NR)m≡μρcH,(4.13)f≡Tc,F≡(Rc)(Pi-pN),τ≡tcρR2H,a≡R2N2.


## 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 L=ḞF-1, the 2D stretching tensor (4.14)D=12(L+LT)=[1λdλdt001λdλdt],
and the viscohyperelastic response is given by(4.15)T=ht=1JF⋅∂w∂E⋅FT+2μhD=cf(λ)+2μH1λ3dλdt,
where T is the Cauchy stress resultant tensor, t is the Cauchy stress, f is the same as in the hyperelastic model, μ is the “solid viscosity” of the aneurysm, and h=H/λ2 since the material is incompressible. This yields the following nondimensionalized governing equation:(4.16)(1x2+bx(β-ln⁡x))ẍ+(b(β-ln⁡x)+bax22-12b)ẋ2+2mfẋx+f(x)x+2mdẋx4=F(τ),
where μf is the viscosity of the CSF, d=H/A, mf=μf/ρch, m=μ/ρch, 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 y0≡x and y1≡dx/dτ. Then the final system of equations becomes
(5.1)dy0dτ=y1,dy1dτ=F(τ)+0.5by12-2m(y1/y0)-f(y0)/y0-bβy12+bln⁡y0y12-bay02y12/2y0-2+by0(β-ln⁡(y0)).
We consider an autonomous problem, that is, F(τ)≡F0 for some constant F0. Then we have(5.2)f(x)=ΓeQ(x2-1),Q=12Γ(x2-1)2.
The derivative can be computed as(5.3)∂f∂y0|(α,0)=2ΓαeQ[Γ(α2-1)2+1].
We now expand the first-order equations in Taylor series around the equilibrium point (α,0) to get (5.4)ẏ0(y0,y1)=G(α,0)+∂G∂y0|(α,0)(y0-α)+∂G∂y1|(α,0)(y1-0)+HOT,ẏ1(y0,y1)=H(α,0)+∂H∂y0|(α,0)(y0-α)+∂H∂y1|(α,0)(y1-0)+HOT,
where G and H represent the functions in the system (5.1) and HOT denotes higher-order terms, which we will neglect. Since G(α,0)=H(α,0)=0 the final linearized system of equations can be written as(5.5){ẏ0ẏ1}=[01F0-∂f/∂y0|(α,0)1/α+bα2(β-ln⁡α)-2m1/α+bα2(β-ln⁡α)]{y0-αy1},
where F0 is the constant pressure at the equilibrium (α,0); hence F0=f(α)/α. We denote the Jacobian matrix J. Let τr=trJ and let δ=detJ. Now τr<0 if and only if the CSF is viscous (i.e., m>0). Next, we will prove stability results for inviscid and viscous CSF.Lemma 5.1.
Letℍ={(x,y)∈ℝ2:x≥1}. Also, let b>0 and suppose 1/α+bα2(β-lnα)≠0. Then the nondimensionalized system (4.11) and (4.16) are C2(ℍ).Proof .
Sinceb>0, x≥1, and 1/α+bα2(β-lnα)≠0forallx∈ℍ, the first and second derivatives of (4.11) exist and are continuous on ℍ.Lemma 5.2.
Supposem>0. Then the point (α,0) is a hyperbolic equilibrium point.Proof.
Since the trace of the linearizationτr<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.
Letm=0 and suppose (α,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose N≥Rα, so that β≥lnα. Then for the system (5.1) one has the following.
(1) If(ΓeQ(2αΓ(α2-1)2+α+α-1))/(1/α+bα2(β-lnα))>0, the nonlinear system has a center at (α,0).
(2) If(ΓeQ(2αΓ(α2-1)2+α+α-1))/(1/α+bα2(β-lnα))<0, the nonlinear system has a saddle point at (α,0).Proof.
Sincem=0,  τr=0. Thus the eigenvalues
(5.6)Λ=τr±τr2-4δ2=±-δ=±-ΓeQ(2αΓ(α2-1)2+α+α-1)1/α+bα2(β-ln⁡α).
(1) SinceΓeQ(2αΓ(α2-1)2+α+α-1)/(1/α+bα2(β-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 (τ,y1)→(-τ,-y1), it is symmetric with respect to the x-axis. Now, the nonlinear system is C1(ℍ) by Lemma 5.1. Therefore, the point (α,0) must also be a center for the nonlinear system.
(2) SinceΓeQ(2αΓ(α2-1)2+α+α-1)/(1/α+bα2(β-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 C1(ℍ) 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,m>0. Also let Γ>0 and suppose (α,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose N≥Rα, 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(5.7)m2=[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable proper node.(2)
If(5.8)m2<[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable focus.(3)
If(5.9)m2>[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable node.Proof.
The eigenvalues of the associated linearized system (see (5.5)) are
(5.10)Λ=τr±τr2-4δ2=-2m/A±(-2m/A)2-4(ΓeQ(2αΓ(α2-1)2+α+α-1)/A)2,
where 𝔄 denotes (α-1+bβα2-bα2lnα). First, note that since m,b>0,β≥lnα,τr is negative for all three cases.
(1)
If(5.11)m2=[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then the linear system has the repeated real eigenvalue
(5.12)Λ=-mα-1+bβα2-bα2ln⁡α<0,
and the linearized system has an asymptotically stable proper node at (α,0). Now, the nonlinear system is C2(ℍ) 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(5.13)m2<[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
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 C2(ℍ) 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(5.14)m2>[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then the linear system has a pair of real eigenvalues. Since Γ>0,
(5.15)m>m2-[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α]
and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (α,0). Now, the nonlinear system is C2(ℍ) 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,m>0. Also suppose F(τ)<f(x)/x and N≥Rα, so that β≥lnα. Then the equilibrium point (α,0) is asymptotically stable.Proof.
Using a result of Freitas's [18],
(5.16)L(x,ẋ)=2∫0xc(r)B(r)dr+B(x)(ẋ)2,
with
(5.17)B(x)=e2∫0xb(r)dr,
is a Lyapunov function for systems of the form
(5.18)ẍ+[a(x)+b(x)ẋ]ẋ+c(x)=0,
where a,b,c satisfy the conditions
(1)a,b,c:ℝ→ℝ and are continuous,
(2)a(x)>0 except possibly at points where c(x)=0,
(3)xc(x)>0 outside a bounded interval.
Applying this to the nondimensionalized governing equation (4.11) with
(5.19)a(x)=2m/x1/x2+bx(β-ln⁡x),b(x)=b(β-ln⁡x)+bax2/2-(1/2)b1/x2+bx(β-ln⁡x),c(x)=f(x)/x1/x2+bx(β-ln⁡x)-F(τ)1/x2+bx(β-ln⁡x),
yields L(α,0)=0 and L̇(x,ẋ)<0 on an open region around the equilibrium and the result follows.Corollary 5.6.
Suppose the CSF is viscous, that is,m>0. Also suppose F(τ)≡0 and β≥0. Then the equilibrium point (1,0) is asymptotically stable.Proof.
Sincef(x)/x>0 on {(x,y)∈ℝ2:x>0&(x,y)≠(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(5.20)dy0dτ=y1dy1dτ=F(τ)+0.5by12-2mf(y1/y0)-f(y0)/y0-bβy12+bln⁡y0y12y0-2+by0(β-ln⁡(y0))-(bay02y12/2)+2md(y1/y04)y0-2+by0(β-ln⁡(y0)).
Linearization of this system of equations yields(5.21){ẏ0ẏ1}=[01F0-∂f/∂y0|(α,0)α-1+bα2(β-bln⁡α)-2mfα3-2mdα2+bβα5-bα5ln⁡α]{y0-αy1}.
Here the Jacobian matrix J has trace τr=trJ=(-2mfα3-2md)/(α2+bβα5-bα5lnα) and determinant δ=detJ=(-F0+(∂f/∂y0)|(α,0))/(α-1+bα2(β-blnα)).Here we must haveF0=f(α)/α at the equilibrium (α,0) and using f(x)=ΓeQ(x2-1), we have (5.22)δ=ΓeQ(2αΓ(α2-1)2+α+α-1)α-1+bα2(β-bln⁡α)>0,
and we can analyze the stability of the linearized system.We now have an immediate analogue of Lemma5.2.Lemma 5.7.
Supposem>0 or mf>0 and β≥lnα. Then the point (α,0) is a hyperbolic equilibrium point.Proof.
Since the trace of the linearizationτr<0, and δ≠0, all of the eigenvalues of the linearized system (5.21) have nonzero real part.Also, note that the traceτr<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,m,mf>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(5.23)mf2α6+m2d2+2mfmdα3=(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable proper node.(2)
If(5.24)mf2α6+m2d2+2mfmdα3<(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable focus.(3)
If(5.25)mf2α6+m2d2+2mfmdα3>(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable node.Proof.
The eigenvalues of the associated linearized system (see (5.21)) are
(5.26)Λ=τr±τr2-4δ2=-mfα3-mdB±((-2mfα3-2md)/B)2-4((-F0+∂f/∂y0|(α,0))/(α-1+bα2(β-bln⁡α)))2,
where 𝔅 denotes (α2-b(βα5-α5lnα)). First, note that since m,mf>0,τr is negative for all three cases.
(1)
If(5.27)mf2α6+m2d2+2mfmdα3=(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
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 C2(ℍ) 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(5.28)mf2α6+m2d2+2mfmdα3<(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
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 C2(ℍ) 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(5.29)mf2α6+m2d2+2mfmdα3>(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),then the linear system has a pair of real eigenvalues. Since
(5.30)Γ>0,mfα3+md>mf2α6+m2d2+2mfmdα3-(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1))
and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (α,0). Now, the nonlinear system is C2(ℍ) 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.


## 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 y0≡x and y1≡dx/dτ. Then the final system of equations becomes
(5.1)dy0dτ=y1,dy1dτ=F(τ)+0.5by12-2m(y1/y0)-f(y0)/y0-bβy12+bln⁡y0y12-bay02y12/2y0-2+by0(β-ln⁡(y0)).
We consider an autonomous problem, that is, F(τ)≡F0 for some constant F0. Then we have(5.2)f(x)=ΓeQ(x2-1),Q=12Γ(x2-1)2.
The derivative can be computed as(5.3)∂f∂y0|(α,0)=2ΓαeQ[Γ(α2-1)2+1].
We now expand the first-order equations in Taylor series around the equilibrium point (α,0) to get (5.4)ẏ0(y0,y1)=G(α,0)+∂G∂y0|(α,0)(y0-α)+∂G∂y1|(α,0)(y1-0)+HOT,ẏ1(y0,y1)=H(α,0)+∂H∂y0|(α,0)(y0-α)+∂H∂y1|(α,0)(y1-0)+HOT,
where G and H represent the functions in the system (5.1) and HOT denotes higher-order terms, which we will neglect. Since G(α,0)=H(α,0)=0 the final linearized system of equations can be written as(5.5){ẏ0ẏ1}=[01F0-∂f/∂y0|(α,0)1/α+bα2(β-ln⁡α)-2m1/α+bα2(β-ln⁡α)]{y0-αy1},
where F0 is the constant pressure at the equilibrium (α,0); hence F0=f(α)/α. We denote the Jacobian matrix J. Let τr=trJ and let δ=detJ. Now τr<0 if and only if the CSF is viscous (i.e., m>0). Next, we will prove stability results for inviscid and viscous CSF.Lemma 5.1.
Letℍ={(x,y)∈ℝ2:x≥1}. Also, let b>0 and suppose 1/α+bα2(β-lnα)≠0. Then the nondimensionalized system (4.11) and (4.16) are C2(ℍ).Proof .
Sinceb>0, x≥1, and 1/α+bα2(β-lnα)≠0forallx∈ℍ, the first and second derivatives of (4.11) exist and are continuous on ℍ.Lemma 5.2.
Supposem>0. Then the point (α,0) is a hyperbolic equilibrium point.Proof.
Since the trace of the linearizationτr<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.
Letm=0 and suppose (α,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose N≥Rα, so that β≥lnα. Then for the system (5.1) one has the following.
(1) If(ΓeQ(2αΓ(α2-1)2+α+α-1))/(1/α+bα2(β-lnα))>0, the nonlinear system has a center at (α,0).
(2) If(ΓeQ(2αΓ(α2-1)2+α+α-1))/(1/α+bα2(β-lnα))<0, the nonlinear system has a saddle point at (α,0).Proof.
Sincem=0,  τr=0. Thus the eigenvalues
(5.6)Λ=τr±τr2-4δ2=±-δ=±-ΓeQ(2αΓ(α2-1)2+α+α-1)1/α+bα2(β-ln⁡α).
(1) SinceΓeQ(2αΓ(α2-1)2+α+α-1)/(1/α+bα2(β-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 (τ,y1)→(-τ,-y1), it is symmetric with respect to the x-axis. Now, the nonlinear system is C1(ℍ) by Lemma 5.1. Therefore, the point (α,0) must also be a center for the nonlinear system.
(2) SinceΓeQ(2αΓ(α2-1)2+α+α-1)/(1/α+bα2(β-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 C1(ℍ) 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,m>0. Also let Γ>0 and suppose (α,0) is the equilibrium point in the region ℍ for the nondimensionalized system (5.1). Also suppose N≥Rα, 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(5.7)m2=[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable proper node.(2)
If(5.8)m2<[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable focus.(3)
If(5.9)m2>[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then (α,0) is an asymptotically stable node.Proof.
The eigenvalues of the associated linearized system (see (5.5)) are
(5.10)Λ=τr±τr2-4δ2=-2m/A±(-2m/A)2-4(ΓeQ(2αΓ(α2-1)2+α+α-1)/A)2,
where 𝔄 denotes (α-1+bβα2-bα2lnα). First, note that since m,b>0,β≥lnα,τr is negative for all three cases.
(1)
If(5.11)m2=[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then the linear system has the repeated real eigenvalue
(5.12)Λ=-mα-1+bβα2-bα2ln⁡α<0,
and the linearized system has an asymptotically stable proper node at (α,0). Now, the nonlinear system is C2(ℍ) 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(5.13)m2<[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
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 C2(ℍ) 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(5.14)m2>[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
then the linear system has a pair of real eigenvalues. Since Γ>0,
(5.15)m>m2-[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α]
and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (α,0). Now, the nonlinear system is C2(ℍ) 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,m>0. Also suppose F(τ)<f(x)/x and N≥Rα, so that β≥lnα. Then the equilibrium point (α,0) is asymptotically stable.Proof.
Using a result of Freitas's [18],
(5.16)L(x,ẋ)=2∫0xc(r)B(r)dr+B(x)(ẋ)2,
with
(5.17)B(x)=e2∫0xb(r)dr,
is a Lyapunov function for systems of the form
(5.18)ẍ+[a(x)+b(x)ẋ]ẋ+c(x)=0,
where a,b,c satisfy the conditions
(1)a,b,c:ℝ→ℝ and are continuous,
(2)a(x)>0 except possibly at points where c(x)=0,
(3)xc(x)>0 outside a bounded interval.
Applying this to the nondimensionalized governing equation (4.11) with
(5.19)a(x)=2m/x1/x2+bx(β-ln⁡x),b(x)=b(β-ln⁡x)+bax2/2-(1/2)b1/x2+bx(β-ln⁡x),c(x)=f(x)/x1/x2+bx(β-ln⁡x)-F(τ)1/x2+bx(β-ln⁡x),
yields L(α,0)=0 and L̇(x,ẋ)<0 on an open region around the equilibrium and the result follows.Corollary 5.6.
Suppose the CSF is viscous, that is,m>0. Also suppose F(τ)≡0 and β≥0. Then the equilibrium point (1,0) is asymptotically stable.Proof.
Sincef(x)/x>0 on {(x,y)∈ℝ2:x>0&(x,y)≠(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(5.20)dy0dτ=y1dy1dτ=F(τ)+0.5by12-2mf(y1/y0)-f(y0)/y0-bβy12+bln⁡y0y12y0-2+by0(β-ln⁡(y0))-(bay02y12/2)+2md(y1/y04)y0-2+by0(β-ln⁡(y0)).
Linearization of this system of equations yields(5.21){ẏ0ẏ1}=[01F0-∂f/∂y0|(α,0)α-1+bα2(β-bln⁡α)-2mfα3-2mdα2+bβα5-bα5ln⁡α]{y0-αy1}.
Here the Jacobian matrix J has trace τr=trJ=(-2mfα3-2md)/(α2+bβα5-bα5lnα) and determinant δ=detJ=(-F0+(∂f/∂y0)|(α,0))/(α-1+bα2(β-blnα)).Here we must haveF0=f(α)/α at the equilibrium (α,0) and using f(x)=ΓeQ(x2-1), we have (5.22)δ=ΓeQ(2αΓ(α2-1)2+α+α-1)α-1+bα2(β-bln⁡α)>0,
and we can analyze the stability of the linearized system.We now have an immediate analogue of Lemma5.2.Lemma 5.7.
Supposem>0 or mf>0 and β≥lnα. Then the point (α,0) is a hyperbolic equilibrium point.Proof.
Since the trace of the linearizationτr<0, and δ≠0, all of the eigenvalues of the linearized system (5.21) have nonzero real part.Also, note that the traceτr<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,m,mf>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(5.23)mf2α6+m2d2+2mfmdα3=(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable proper node.(2)
If(5.24)mf2α6+m2d2+2mfmdα3<(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable focus.(3)
If(5.25)mf2α6+m2d2+2mfmdα3>(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
then (α,0) is an asymptotically stable node.Proof.
The eigenvalues of the associated linearized system (see (5.21)) are
(5.26)Λ=τr±τr2-4δ2=-mfα3-mdB±((-2mfα3-2md)/B)2-4((-F0+∂f/∂y0|(α,0))/(α-1+bα2(β-bln⁡α)))2,
where 𝔅 denotes (α2-b(βα5-α5lnα)). First, note that since m,mf>0,τr is negative for all three cases.
(1)
If(5.27)mf2α6+m2d2+2mfmdα3=(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
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 C2(ℍ) 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(5.28)mf2α6+m2d2+2mfmdα3<(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),
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 C2(ℍ) 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(5.29)mf2α6+m2d2+2mfmdα3>(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1)),then the linear system has a pair of real eigenvalues. Since
(5.30)Γ>0,mfα3+md>mf2α6+m2d2+2mfmdα3-(α2+bβα5-bα5ln⁡α)2α-1+bα2(β-bln⁡α)(ΓeQ(2αΓ(α2-1)2+α+α-1))
and the eigenvalues are both negative. Thus the linearized system has an asymptotically stable node at (α,0). Now, the nonlinear system is C2(ℍ) 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 c and Γ taken from [19]: ρ=1050 kg/m3, R=3×10-3 m, H=27.8×10-6 m, ρf=1000 kg/m3, pN=399.9375 N/m2,c=10.18 N/m, Γ=23.6, μf=1.26×10-4 Ns/m2, N=105m, and the value for the “solid viscosity” μ=331.827 Ns/m2 The mean blood pressure Pm=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.Figure2 demonstrates that greatly increasing the distance N 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 N to vary, as long as R≪N. The physical meaning of the parameter N 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 N 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 N.Table 1
Eigenvalues for models with varyingN.
N= 105 mEigenvaluesHyperelastic model-0.000001±0.189731iViscoelastic model-0.002453±0.189715iN= 103 mEigenvaluesHyperelastic model-0.000001±0.221537iViscoelastic model-0.003344±0.221512iN= 1 mEigenvaluesHyperelastic model-0.000001±0.351127iViscoelastic model-0.008398±0.351026iN= 0.1 mEigenvaluesHyperelastic model-0.000002±0.732381iViscoelastic model-0.036509±0.731470iStretch ratio with varyingN.
(a)
Viscoelastic model,N=0.1m(b)
Viscoelastic model,N=1m(c)
Viscoelastic model,N=103m(d)
Viscoelastic model,N=105mFor the values ofN chosen in the following numerical experiments, the nondimensional parameter a is very small. This causes the nonlinear term involving y02y12 to have only a small effect on the solution. For example, for N=105m, as used to generate the solution plots, we have a=9×10-16, which makes the contribution due to this term nearly negligible.We considered an initial disturbance atτ=0 by allowing (y0,y1) to vary between (1.01,0) and (1.02,0). Using these values, and N=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.Phase portraits for inviscid CSF, viscous CSF, and viscoelastic model.
(a)
Inviscid CSF(b)
Viscous CSF(c)
Viscous CSF,mf=10(d)
Viscoelastic modelThe numerically generated phase portraits for the inviscid and viscoelastic cases (Figures3(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-6, which is very close to zero. Thus the solution converges to the equilibrium very slowly relative to the other cases. Using the value mf=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 y0≈0.0240839, y1=0, which is not physically relevant since y0<1. The eigenvalues of the Jacobian matrix can be computed as Λ1≈-0.255257, Λ2≈0.255259, indicating that the linearized system indeed has a saddle at the equilibrium.Table2 shows the different possibilities for the eigenvalues of the Jacobian matrix for the hyperelastic model for various values of mf with N=105 m. The theoretical value given is computed from the condition
(6.1)m2=[ΓeQ(2αΓ(α2-1)2+α+α-1)][α-1+bβα2-bα2ln⁡α],
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 mf computed from the numerics beyond which the equilibrium becomes a node is close to the theoretical value, lying between mf=373 and mf=373.666.Table 2
Eigenvalues for hyperelastic model with varyingmf.
mf=372.3979903Eigenvalue 1Eigenvalue 2Theoretical value−.1897424481−.1897424481mf=375Eigenvalue 1Eigenvalue 2Numerical model−0.174348−0.206469mf=373.666Eigenvalue 1Eigenvalue 2Numerical model−0.189063−0.190400mf=373Eigenvalue 1Eigenvalue 2Numerical model−0.189393 + 0.011303i−0.189393 − 0.011303iPlots of the stretch ratioλ are shown in Figures 5, 6, and 7 for various values of m, R, and ω. The actual equilibrium point is at approximately (y0,y1)=(1.0531,0); however, we examine stability by considering a perturbed system with initial condition (y0,y1)=(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.Stretch ratio for inviscid CSF, viscous CSF, and viscoelastic model.
(a)
Inviscid CSF(b)
Viscous CSF(c)
Viscoelastic model(d)
Viscoelastic model,mf=10Stretch ratio for viscoelastic model with varyingm.
(a)
m = 1(b)
m = 10(c)
m = 100(d)
m=4mfStretch ratio for viscoelastic model with varying radii.
(a)
R=0.003m(b)
R=0.005m(c)
R=0.03m(d)
R=0.05mStretch ratio with varyingω.
(a)
ω=4πHz(b)
ω=10πHz(c)
ω=100Hz(d)
ω=400HzFigure5 shows the effect of varying the nondimensional parameter m on the stretch ratio λ. Larger values of m cause the solution to approach the equilibrium state more rapidly. The value m=4mf was used in [11] to illustrate the effect of a viscoelastic wall on the stability of the solution. However, the much larger values of m 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 radiusR in the viscoelastic model is seen in Figure 6. Larger values of R 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.Figure7 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 (y0,y1)=(1.029,0); however, we follow [11] and examine stability by considering a perturbed system with initial condition (y0,y1)=(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.
CylinderEigenvaluesHyperelastic model-0.000001±0.189731iViscoelastic model-0.002453±0.189715iSphereEigenvaluesHyperelastic model-0.000005±0.987053iViscoelastic model-0.094307±0.982538iStretch ratio for cylindrical versus spherical geometry.
(a)
Viscous CSF, cylinder(b)
Viscous CSF, sphere(c)
Viscoelastic model, cylinder(d)
Viscoelastic model, sphere


## 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.



---

*Source: 906475-2011-09-25.xml*