6th Australasian Congress on Applied Mechanics, ACAM 6 12–15 December 2010, Perth, Australia

Wave propagation in an elastic waveguide: fluid-structure interactions in a spinal disease Novak S. J. Elliott1 *, Anthony D. Lucey1 , Duncan A. Lockerby2 and Andrew R. Brodbelt3 1 Fluid Dynamics Research Group, Curtin University of Technology, Perth, Australia 2 Fluid Dynamics Research Centre, School of Engineering, University of Warwick, Coventry, U.K. 3 The Walton Centre for Neuroradiology and Neurosurgery NHS Trust, Liverpool, U.K.

* Corresponding author. Email: [email protected] Abstract: Syringomyelia is a disease in which fluid-filled cavities, called syrinxes, form in the spinal cord (SC). The progressive expansion of syrinxes over many years compresses the surrounding nerve fibres and blood vessels, which is associated with neurological damage. In the present work we aim to elucidate the mechanics underlying syrinx formation and expansion by investigating the wave-propagation characteristics of the spinal system in healthy and diseased configurations. We use the standard biomechanical analogue consisting of cylindrical, axisymmetric solid and fluid layers. Specifically, the SC is represented as an elastic cylinder, which becomes an annulus containing inviscid fluid when a syrinx is included, and this is surrounded by inviscid fluid representing the cerebrospinal fluid (CSF) occupying the subarachnoid space, bound by a rigid dura. The model is formulated as a system of Helmholtz equations which describe axisymmetric harmonic motion of the cylindrical layers. These equations are discretised using Chebyshev polynomials and then solved as a generalised eigenvalue problem. This linear algebra approach gives explicit access to wave properties like traditional root-finding methods but without the need for a long wave assumption, and is also more computationally efficient than finite element/volume methods used in other spinal models. Our results reproduce the wave speeds of other syringomyelia models and the dispersion diagrams are qualitatively similar to other acoustic models with like topologies. This demonstrates the applicability of the numerical method to the biological problem. Additionally we are able to recover the associated displacement and stress modes from the eigenvectors. This investigation serves as a framework for studying cylindrical waveguides in biological systems. Keywords: Syringomyelia, Helmholtz, Chebyshev, Eigenvalue

1

Introduction

The SC lies within the fluid-filled passageway formed through the vertebrae and is held in place by ligaments and supporting elastic structures. This assembly of soft tissues and CSF is a wave-bearing system that may be excited by pressure sources such as the cardiac cycle at low frequencies and coughs/sneezes at higher frequencies. Previous investigators have studied the wave-propagation characteristics of the spinal anatomy with an analogue model based on a system of axisymmetric cylindrical fluid and solid layers. The mathematical modelling strategies can be grouped into those that used analytical methods [2, 7, 8, 10, 13], such as asymptotic analysis or root finding, and those that used computational methods [3–6, 14], such as finite difference/volume/element. Models of the former type had semi-infinite geometries (1-d, 2-d) and were non-dispersive due to a simplifying long wavelength assumption. In contrast, models of the latter type necessarily had finite geometries (1-d, 2-d) and the full frequency spectrum was included. Syringomyelia is a disease that progresses over months to years. In contrast, mechanical disturbances of the cerebrospinal system are largely driven by percussive and cardiac pressure sources operating over fractions of a second up to several seconds, respectively. This difference in timescales suggests that the role of biomechanical phenomena in syrinx formation and expansion may be a gradual accumulation of very small effects, rather than the effect of an isolated critical event. On this basis we seek

Dura (rigid)

Dura (rigid)

SSS (fluid)

SSS (fluid) Syrinx (fluid)

SC (solid)

SC (solid)

(a) Healthy

(b) Diseased

Figure 1: The two spinal canal configurations studied here.

to study the wave propagation characteristics of the spinal canal over the whole frequency spectrum. It is also desirable to obtain this information explicitly and with minimal computational effort. To achieve these aims a model is formulated as a system of Helmholtz equations describing axisymmetric harmonic motion of the cylindrical layers. These equations are discretised using Chebyshev polynomials and then solved as a generalised eigenvalue problem. This linear algebra approach was devised [1] with industrial problems in mind, and has been applied to geophysical borehole analysis [11]; i.e. nondestructive testing of very stiff materials (rock) in the megahertz frequency spectrum. Here we demonstrate the applicability of this numerical method to very soft materials in the sub-kilohertz range more typical of biological systems.

2

Method

We consider two cross-sectional configurations of the spinal canal as shown in Fig. 1. The ‘healthy’ model consists of an elastic cylinder representing the SC tissue, surrounded by an annulus of inviscid fluid, representing the CSF in the spinal subarachnoid space (SSS), which in turn has a rigid boundary representing the dura mater. The ‘diseased’ model includes a central cylinder of inviscid fluid representing a syrinx.

2.1

Governing equations

We present the equations for the propagation of time-harmonic waves through a cylindrically layered structure, restricting our attention to the axisymmetric case where displacement is in the r-z plane and wave travel is in the axial (z) direction; i.e., u(r, z, t) = U (r)ei(kz z−ωt) ,

(1)

with kz being the axial wavenumber and ω the angular frequency. We consider two layer types in order of increasing complexity: (i) inviscid fluid, (ii) elastic solid. The Helmholtz equations for wave propagation are derived from displacement potentials.

2.1.1

Fluid

The linearised Navier-Stokes equation for an inviscid fluid is −∇p = ρf u ¨

(2)

and the pressure may be expressed in terms of the fluid displacement, p = −Kf ∇ · u,

(3)

where ρf and Kf are the fluid density and bulk modulus, respectively. The fluid velocity may be expressed in terms of a scalar potential, u˙ = ∇ϕ, ˜ (4) which satisfies (2) if ∇2 ϕ˜ =

1 ¨ ϕ, ˜ c2f

(5)

Kf ; ρf

(6)

where the speed of sound in the fluid is given by s cf =

∇2 is the scalar Laplace operator. The potential ϕ˜ takes the same separable form as the displacement in (1), i(kz z−ωt) ˜ ϕ(r, ˜ z, t) = Φ(r)e , (7) which reduces the wave equation (5) into a Helmholtz ordinary differential equation, 

 1 d ω2 d2 + + 2 Φ = kz2 Φ, dr2 r dr cf | {z }

(8)

Lf

˜ where Φ = Φ/(−iω). Substituting (7) and (8) into (4), and noting ∂/∂t = −iω for our harmonic system (7), we arrive at the respective radial and axial displacement components, dΦ ur = , dr

 u ˆz = −

d2 1 d ω2 + + dr2 r dr c2f

 Φ,

(9a,b)

where u ˆz = ikz uz is introduced to avoid functions of complex numbers. Similarly the pressure may be obtained by substituting (3), (4), (6) and (7) into (2), p = ω 2 ρf Φ.

2.1.2

(10)

Solid

The linear theory of elasticity gives the momentum equation as ∇ · σ = ρ¨ u,

(11)

where the stress tensor for an isotropic solid is σ = λ(∇ · u)I + µ(u∇ + ∇u),

(12)

and λ and µ are Lamé’s first parameter and the shear modulus (or Lamé’s second parameter), respectively. For spinal tissues it proves more convenient to work from the Young’s modulus E and the Poisson’s ratio ν; i.e., Eν E λ= , µ= . (13a,b) (1 + ν)(1 − 2ν) 2(1 + ν) The displacement may be decomposed into a scalar potential and a vector potential, u = ∇ϕ + ∇ × ψ,

∇ · ψ = 0,

(14a,b)

which satisfies (11) if these potentials obey the respective wave equations, 1 ϕ, ¨ c2d

∇2 ϕ = where

s cd =

λ + 2µ , ρ

∇2 ψ =

1 ¨ ψ, c2s r

cs =

µ , ρ

(15a,b)

(16a,b)

correspond to the speeds of dilatational (d) and shear (s) waves in an infinite elastic medium∗ ; ∇2 is the vector Laplace operator. For axisymmetric motion ψ = {0, ψθ , 0} and thus we denote ψ = ψθ . By expressing the pair of scalar potentials in harmonic form, ϕ = Φ(r)ei(kz z−ωt) ,

ψ = Ψ(r)ei(kz z−ωt) ,

(17a,b)

the wave equations (15a,b) become a pair of Helmholtz equations  1 d d2 ω2 + + Φ dr2 r dr c2d | {z }

= kz2 Φ,

(18a)

 1 d 1 ω2 d Ψ + − + dr2 r dr r2 c2s {z } |

= kz2 Ψ.

(18b)



Ld



2

Ls

Substituting (17a,b) into (14a) gives the radial and axial displacements, ur u ˆz

dΦ ˆ − Ψ, dr  2    d 1 d ω2 d 1 ˆ = − + + 2 Φ+ + Ψ, dr2 r dr cd dr r =

(19a) (19b)

likewise substituting (14a) and (17a,b) into (12) gives the normal and shear stress, σrr

=

σ ˆrz

=

  ˆ d2 ω2 dΨ , 2µ 2 − λ 2 Φ − 2µ dr cd dr  3   2  d 1 d2 1 d ω2 d 2 d 2 ω2 ˆ d −2µ + − 2 + 2 − + 2 Ψ; Φ+µ 2 2 + dr3 r dr2 r dr cd dr dr r dr r2 cs

(20a) (20b)

ˆ = ikz Ψ and σ Ψ ˆrz = ikz σrz are introduced for mathematical convenience.

2.1.3

Boundary conditions

At the interfaces between adjacent layers we require continuity of displacement and stress. Of relevance here is the interface between a solid layer and a fluid layer thus ur |s = ur |f ,

σrr |s = −p|f ,

σrz |s = 0.

(21a,b,c)

Likewise we require boundary conditions at the outer surface of our cylindrically layered structure. For an elastic solid cord in vacuo the stress-free condition is σrr = σrz = 0, ∗ Strictly

(22a,b)

speaking dilatational waves are best described as irrotational and shear waves as equivolumal, although they are variously referred to in different literatures as longitudinal/compressional/P waves and transverse/distortional/S waves, respectively.

while for the rigid boundary of a fluid-filled subarachnoid space, representing the dura, the boundary condition is simply ur = 0. (23)

2.2

Numerical implementation

We follow the method developed by Adamou and Craster [1] and Karpfinger et al. [11], which is briefly described here. An nth -order derivative operating over a finite interval, in our case a radial thickness, (n) can be approximated at N interpolation points using Chebyshev polynomials; i.e., DN ≈ dn /drn [16]. Using these differentiation matrices we can form discrete approximations to the differential operators Lf , Ls1 and Ls2 in (8), (18a) and (18b), respectively. E.g., the fluid operator becomes (2)

(1)

Lf = DN + diag(1/RN )DN + diag(ω 2 /c2f ) ≈ Lf ,

(24a,b)

where RN ≈ r is the vector of radial interpolation points. For the two-layered case of an elastic cylinder (SC) surrounded by an annulus of fluid (SSS) the set of Helmholtz equations may be expressed as LΘ = kz2 Θ, where 

Ls1 L= 0 0

0 Ls2 0

 0 0 , Lf

(25)    x(s1)  Θ= x  (s2)  x(f)

(26a,b)

and xs1 ≈ Φs , xs2 ≈ Ψs and xf ≈ Φf are the discrete approximations of the displacement potential amplitudes. Following a similar approach to that outlined above the matrices for the displacement (T ) and stress (S) differential operators may be deduced. To incorporate the boundary conditions into (25) the rows of L corresponding to the interfaces and outer surface are substituted for differences of the ˜ By zeroing out the same rows of an identity respective rows of T and S, giving the modified matrix L. matrix of the same dimension, M , the boundary value problem can be written as ˜ = kz2 M Θ. LΘ

(27)

In other words, a generalised eigenvalue problem where the phase velocities are given by the eigenvalues, ω , (28) c= Re(kz ) and the displacement and stress mode shapes by the eigenvectors. We use standard MATLAB routines to compute the solution.

2.3

Physiological parameters

The values of the physiological parameters are summarised in Table 1. The fluid, representing the CSF in the subarachnoid space and also fluid in a syrinx, has been assigned the density and bulk modulus of water in accordance with physical measurements [12]. The solid corresponds to the combined SC parenchyma and pial membrane surrounding it. Tissue is largely water by mass hence the density, the value of Young’s modulus is taken from Ozawa et al. [15], and we follow Bertram et al. [5] in choosing a Poisson’s ratio for an almost perfectly incompressible material. The cross-sectional dimensions of the SC and dura are taken from Elliott et al. [10] while the radius of the syrinx is simply chosen to be half of that of the SC.

Table 1: Parameter values. Medium Fluid Solid

(Geometry)

Parameter ρf Kf ρ E ν rsyrinx R rdura

Value 1000 kg/m3 2.3 GPa 1000 kg/m3 15 kPa 0.49 2 mm 4 mm 7.5 mm

Description Density Bulk modulus Density Young’s modulus Poisson’s ratio Radius of syrinx Radius of spinal cord Radius of dura mater

c / c0

1.5

Experiment Numerical

1

0.5 0

0.5

1

1.5

2

R/Λ

Figure 2: Validation of the numerical code against published experiment data [9] in computing the dispersion diagram of axisymmetric waves in a steel cylinder; R is the radius, Λ is the wavelength, c is the phase speed and c0 is the bar velocity.

3

Results and Discussion

To validate our code we compute the dispersion relation for the benchmark case of a steel cylinder (radius R) and compare it to published experimental results [9]. As can be seen in Fig. 2 the first three modes agree well. At long wavelengths (R/Λ  p 1) only the first mode is propagatory and the phase speed approaches the ‘bar velocity’, c → c0 = E/ρ. At shorter wavelengths there are an infinite number of propagatory modes, all of which approach the speed of Raleigh surface waves at R/Λ  1. The equivalent plot for a cylinder having the material properties of SC is shown in Fig. 3(a); the dashed p line is E/ρ. Using this plot as a reference we introduce the dispersion diagrams for the healthy and diseased models in Fig. 3(b) and (c), respectively. The pairs of dashed lines in each of these plots correspond to the long-wave solutions derived analytically by Cirovic [8]. In Fig. 3(b) the lower dashed line denotes the speed of shear waves in the SC while the upper dashed line denotes the bar speed scaled by the ratio of cross-sectional areas of the SSS to the SC. The shear wave does not appear in the eigenvalue solution as this only applies to an elastic solid annulus (SC) with a vanishingly small fluid core (syrinx), as opposed to the topologically different case of there being no fluid core. In Fig. 3(c) the lower and upper dashed lines denote the modified Young’s mode and Lamb’s mode of the SC, respectively. Figures 4 and 5 show the displacement and stress mode shapes for the first mode in the healthy and diseased models, respectively. The wavelengths are chosen to correspond to the length of the SC (R/Λ = 0.025), the radius of the SC (R/Λ = 1) and a value in between (R/Λ = 0.13). We make the following observations: • At long wavelengths axial displacement dominates over radial (uz  ur ), normal stress dominates over shear (σrr  σrz ) but these become of the same order of magnitude as the wavelenth approaches the radius of the SC.

5

c (m/s)

4

3

2

1

0

0

0.2

0.4

0.6

0.8

1

0.2

0.4

R/Λ

5

5

4

4

c (m/s)

c (m/s)

(a) SC in vacuo

3

3

2

2

1

1

0

0

0.2

0.4

0.6

R/Λ

(b) Healthy

0.8

1

0

0

0.6

0.8

1

R/Λ

(c) Diseased

Figure 3: Dispersion plots for spinal canal models; R is the radius of the spinal cord, Λ is wavelength, c is phase speed.

• Fluid in the syrinx and SSS moves in the opposite axial direction to the SC (uz ) indicating a sloshing motion, with the syrinx fluid moving faster than the SSS fluid relative to the SC. Williams [17] proposed that such a sloshing motion of syrinx fluid may be a mechanism for syrinx elongation. The SSS fluid slows down when a syrinx is introduced due to the inertial contribution of the syrinx. • The peak normal stress occurs at the syrinx wall. If this exceeds the yield stress of the cord material then the accumulative effect of repeated disturbances, such as pressure impulses generated from coughs/sneezes or the cardiac cycle, may be to expand the syrinx. Similarly, the normal stress gradients in the SC are much higher for the diseased case — these stress gradients would promote the degradation of the tissue through material fatigue that might permit a syrinx to expand.

4

Conclusion

Our results reproduce the wave speeds of other syringomyelia models and the dispersion diagrams are qualitatively similar to other acoustic models with like topologies. This demonstrates the applicability of the numerical method to the biological problem. This investigation serves as a framework for studying cylindrical waveguides in biological systems.

3000

0.8

2500

0.6

2000

stress / max(σrz)

displacement / max(uz)

1

0.4 0.2 0

1500 1000 500

ur

−0.2

σrr

0

σrz

uz −0.4

0

0.5

1

1.5

−500

2

0

0.5

r/R

1

1.5

2

r/R

(a) Displacement mode, R/Λ = 0.025

(b) Stress mode, R/Λ = 0.025

1

600

stress / max(σrz)

displacement / max(uz)

500

0.5

0

400 300 200 100

ur

σrr

0

σrz

uz −0.5

0

0.5

1

1.5

−100

2

0

0.5

r/R

1

1.5

2

r/R

(c) Displacement mode, R/Λ = 0.13

(d) Stress mode, R/Λ = 0.13 1.4

1.5

stress / max(σrz)

displacement / max(uz)

1.2 1

0.5

0

1 0.8 0.6 0.4

−0.5

ur

σrr

0.2

σrz

uz −1

0 0

0.5

1

1.5

r/R

(e) Displacement mode, R/Λ = 1

2

0

0.5

1

1.5

2

r/R

(f) Stress mode, R/Λ = 1

Figure 4: Mode shapes for mode 1 of the healthy model; R is the radius of the spinal cord, Λ is wavelength; SC: 0 < r/R < 1, SSS: 1 < r/R < 1.875.

4

x 10 9

1

8 7

0.6

stress / max(σrz)

displacement / max(uz)

0.8

0.4 0.2 0 −0.2 −0.4 ur

−0.6 −0.8

uz

6 5 4 3 2 1

σrr

0

σrz

−1 0

0.5

1

1.5

2

0

0.5

1

1.5

2

r/R

r/R

(a) Displacement mode, R/Λ = 0.025

(b) Stress mode, R/Λ = 0.025

1

1200

0.8

stress / max(σrz)

displacement / max(uz)

1000

0.6 0.4 0.2 0 −0.2

800 600 400 200

−0.4 ur

−0.6 −0.8

σrr

0

σrz

uz 0

0.5

1

1.5

−200

2

0

0.5

r/R

(c) Displacement mode, R/Λ = 0.13

1.5

2

(d) Stress mode, R/Λ = 0.13 10

0.6 0.4

8

0.2

stress / max(σrz)

displacement / max(uz)

1

r/R

0 −0.2 −0.4

6

4

2

−0.6

−1

σrr

0

ur

−0.8

σrz

uz

−2 0

0.5

1

1.5

r/R

(e) Displacement mode, R/Λ = 1

2

0

0.5

1

1.5

2

r/R

(f) Stress mode, R/Λ = 1

Figure 5: Mode shapes for mode 1 of the diseased model; R is the radius of the spinal cord, Λ is wavelength; Syrinx: 0 < r/R < 0.5, SC: 0.5 < r/R < 1, SSS: 1 < r/R < 1.875.

Acknowledgements We acknowledge the support of both the Australian Research Council through the project DP0559408 and the WA State Centre of Excellence in eMedicine.

References [1] A. T. I. Adamou and R. V. Craster. Spectral methods for modelling guided waves in elastic media. Journal of the Acoustical Society of America, 116(3):1524–1535, 2004. [2] K. Berkouk, P. W. Carpenter, and A. D. Lucey. Pressure wave propagation in fluid-filled co-axial elastic tubes part 1: basic theory. Journal of Biomechanical Engineering, 125(6):852–856, 2003. [3] C. D. Bertram. A numerical investigation of waves propagating in the spinal cord and subarachnoid space in the presence of a syrinx. Journal of Fluids and Structures, 25:1189–1205, 2009. [4] C. D. Bertram. Evaluation by fluid/structure-interaction spinal-cord simulation of the effects of subarachnoid-space stenosis on an adjacent syrinx. Journal of Biomechanical Engineering, 132(6):061009–1–15, 2010. [5] C. D. Bertram, A. R. Brodbelt, and M. A. Stoodley. The origins of syringomyelia: numerical models of fluid/structure interactions in the spinal cord. Journal of Biomechanical Engineering, 127(7):1099–1109, 2005. [6] C. D. Bertram, L. E. Bilston, and M. A. Stoodley. Tensile radial stress in the spinal cord related to arachnoiditis or tethering: a numerical model. Medical & Biological Engineering & Computing, 46:701–707, 2008. [7] P.W. Carpenter, K. Berkouk, and A.D. Lucey. Pressure wave propagation in fluid-filled co-axial elastic tubes part 2: mechanisms for the pathogenesis of syringomyelia. Journal of Biomechanical Engineering, 125(6):857–863, 2003. [8] S. Cirovic. A coaxial tube model of the cerebrospinal fluid pulse propagation in the spinal column. Journal of Biomechanical Engineering, 131(2):021008, 2009. [9] R. M. Davies. A critical study of the hopkinson pressure bar. Philosophical Transactions of the Royal Society of London. Series A, 240(821):375–457, 1948. [10] N. S. J. Elliott, D. A. Lockerby, and A. R. Brodbelt. The pathogenesis of syringomyelia: A re-evaluation of the elastic-jump hypothesis. Journal of Biomechanical Engineering, 131(4):044503–1–6, 2009. [11] F. Karpfinger, B. Gurevich, and A. Bakulin. Modeling of wave dispersion along cylindrical structures using the spectral method. Journal of the Acoustical Society of America, 124(2):859–865, 2008. [12] J. A. Kiernan. Barr’s the Human Nervous System, an Anatomical Viewpoint. Lippincott-Raven, Philadelphia, 7th edition, 1998. [13] P. Lockey, G. Poots, and B. Williams. Theoretical aspects of the attenuation of pressure pulses within cerebrospinal fluid pathways. Medical & Biological Engineering & Computing, 13:861–869, 1975. [14] F. Loth, M. A. Yardimci, and N. Alperin. Hydrodynamic model of cerebrospinal fluid motion within the spinal cavity. Journal of Biomechanical Engineering, 123:71–79, 2001. [15] H. Ozawa, T. Matsumoto, T. Ohashi, M. Sato, and S. Kokuban. Mechanical properties and function of the spinal pia mater. Journal of Neurosurgery (Spine 1), 1:122–127, 2004. [16] J. A. C. Weideman and S. C. Reddy. A MATLAB differentiation matrix suite. ACM Transactions on Mathematical Software, 26(4):465–519, 2000. [17] B. Williams. On the pathogenesis of syringomyelia: a review. Journal of the Royal Society of Medicine, 73(11):798, 1980.