- Email: [email protected]
A viscoelastic approach to the modelling of hydrocephalus
Applied Mathematics and Computation 163 (2005) 1097–1107
www.elsevier.com/locate/amc
A viscoelastic approach to the modelling of hydrocephalus S. Sivaloganathan a
a,*
, M. Stastna a, G. Tenti a, J.M. Drake
b
Department of Applied Mathematics, University of Waterloo, Waterloo, Ont., Canada N2L 3G1 b Division of Neurosurgery, Hospital for Sick Children, University of Toronto, Toronto, Ont., Canada M5G 1X8
Abstract In earlier work [Appl. Math. Comput. 94 (1998) 243, Can. Appl. Math. Quart. 7 (1999) 111], we modelled the brain as a porous elastic medium containing an incompressible fluid and solved the simplified problem analytically using consolidation theory. In this paper we apply linear viscoelasticity theory to the modelling of the brain in the same simplified geometry. The aim in both approaches is to be able to eventually accurately predict the ventricular wall configuration in shunted hydrocephalus which would be of practical use to the neurosurgeon in locating the optimal position for shunt implantation in the treatment and management of hydrocephalus. We use the elastic– viscoelastic analogy (EVA) to obtain analytic solutions for the case when the brain is considered to be composed of a viscoelastic solid (the so-called ‘‘standard solid’’). We finally conclude with suggestions of how experimental data could be used to obtain realistic values for the mechanical parameters that appear in the constitutive relation for the viscoelastic solid. 2004 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail address: [email protected] (S. Sivaloganathan).
0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.06.037
1098
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1. Introduction The cerebral system is a complex continuum where a solid phase (brain tissue) interacts with fluid phases (essentially blood and cerebrospinal fluid) in different types of vessels and in the ventricles. A true mathematical description of the physics of such an interaction poses enormous challenges. In part due to this and the inherent complexity of the problem, early researchers replaced the brain with one or more units, or compartments, whose behaviour is represented by a single value of the pressure and by values of the flux exchanged with adjacent compartments. All these values may be time dependent but do not vary spatially. The resistance to the flow from one compartment to an adjacent one is lumped at the boundary between the two compartments. The genesis of this approach can be traced back over two centuries to Munro [7], since then numerous compartment models have appeared in the literature (e.g. see [4,6,9]). These have been recently reviewed by Sivaloganathan et al. [10] where the full mathematical solutions of the most relevant with a critical assessment of their validity was given. It was further demonstrated in [10] that all of these models can be derived as special cases of a general differential equation for the CSF (cerebrospinal fluid) pressure P(t). The weakness of these models, however, is their inability to describe the stress, strain distribution within the brain tissue, as well as their inability to say anything about the ventricular wall configuration in a hydrocephalic brain. Motivated in part by this, a new and more sophisticated generation of mathematical models has emerged over the past decade (see [8]). The starting point, for this approach, is the assumption that the brain parenchyma can be schematized as a sponge-like material. This approach relies heavily on the theory of consolidation introduced by Biot [1], and from a physical standpoint falls squarely into the category of flow through porous media. It is a continuum theory and leads to a set of PDEs which are usually solved by standard numerical techniques. However, Nagashima et al encountered serious problems in their numerical solution of the PDEs. In Tenti et al. [12] and Stastna et al. [11] we addressed these questions and elucidated the underpinnings of consolidation theory as applied to the brain. By solving the model analytically in a simplified geometry, we were able to link the numerical problems to the value of the Poisson ratio m used and provided a new and widely accepted estimate. In this paper we take yet a different approach and treat the brain parenchyma as a viscoelastic solid. We start again with a simplified picture of the brain, by assuming that the brain parenchyma is an impermeable, viscoelastic solid. In order to simplify the analysis as much as possible, we assume that the brain can be represented as a thick-walled cylindrical tube the inside of which is kept at a constant pressure pi (the ventricular CSF pressure) while the outside is kept at a lower constant pressure p0, which w.l.o.g. can be assumed to be zero. Although the geometry may appear rather simplistic, we are assured by neuro-
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1099
surgeons that in the case of hydrocephalus the shape of the lateral ventricular space filled with CSF is very approximately cylindrical. A natural question that arises, is what would be an appropriate viscoelastic model to use to represent the brain parenchyma. The basic Maxwell and Kelvin-Voigt models appear to be inappropriate; the former because it fails to capture the decreasing strain rate under constant stress characteristic of primary creep and the latter since it fails to represent time independent, instantaneous strain, or recovery of strain upon loading or unloading. Further neither of the models captures the very steep initial strain rates characteristic of many viscoelastic materials. In this paper, we propose the ‘‘standard solid’’ (see [2,3]) as an appropriate model. This is schematically represented in Fig. 1. Using Fig. 1, it is straightforward to derive the following constitutive equation for the standard solid: 1 1 g or g oe ¼eþ : ð1Þ þ rþ E1 E2 E1 E2 ot E2 ot Eq. (1) can be written in the general form: P r ¼ Qe; where P and Q are differential operators given by P¼
n X i¼0
Q¼
n X i¼0
pi
oi ; oti
oi qi i : ot
ð2Þ
Generally, viscoelastic stress analysis problems are more involved than elasticity problems due to the inclusion of the time variable in the DEs in addition to the space variables. However, in many problems (including the present one) where the types of BCs and temperature remain constant in time, the time variable can be removed by using Laplace transforms. We thus have an equivalent elastic problem. When a solution has been found in terms of the Laplace transform variable s, the inverse Laplace transform gives the required solution in the
Fig. 1. Schematic representation of standard solid.
1100
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
time variable t for the time dependent behaviour in the viscoelastic problem. This is, in essence, the elastic-viscoelastic analogy (EVA) or the ‘‘correspondence principle’’. For the theoretical development and mathematical justification of the EVA, the reader is referred to Gurtin and Sternberg [5].
2. Quasi-static approximation Inertia terms in the equations of motion are neglected, based on the assumption that the dependent variables vary slowly with respect to time so that the rate of change of momentum can be neglected. Thus we have the equilibrium equation orij þ F j ¼ 0; oxi
ð3Þ
where rij = rij(x, t) is the stress and F the body force. Due to the fact that stress relaxation and heat conduction waves are damped out very quickly in viscoelastic materials, disturbances (from boundary forces or changes in displacements) are transformed into static deformation and heat. Hence one assumes that the viscoelastic body is in equilibrium at every moment. This is the case when changes in b.c.s or changes due to stress relaxation take place slowly enough so that there is no significant change during the decay of the disturbance. An attractive feature in using the quasi-static approximation is that many corresponding elastic problems have been solved––the solution of the elastic problem can then be used to find the required viscoelastic solution. The boundary conditions can be defined by specifying both tractions and surface displacements. We divide the boundary C into two portions, CT where tractions are specified and Cu = CnCT where the displacements are specified, i.e.: T i ðx; tÞ
¼ rij ðx; tÞnj
ui ðx; tÞ
¼ fi ðx; tÞ on
on
CT ;
Cu :
However, the boundary C must be independent of t. It should be noted that if C varies with t, transform methods can still be used; however, Wiener–Hopf techniques will then be required and the great simplicity of the EVA is lost.
3. Mathematical formulation The assumption of a cylindrical geometry naturally suggests the use of cylindrical polar co-ordinates (r, h, z). The cylinder is assumed to be tethered, there are thus no deformations produced as a result of contact forces on its ends. As
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1101
a consequence, the absence of deformations in the z-direction implies the strain distribution will be planar (i.e. the same in any horizontal cross-section). The cylindrical annulus is assumed to be a homogeneous isotropic solid. The displacement and distribution of stresses and strains for this problem have been well known since the time of Lame´ (1852) (see [5]) and these results can be used to solve the problem, where the material is viscoelastic, through the ‘‘correspondence principle’’. Since the elastic problem is very important for the analysis of the viscoelastic problem, we review it in detail before applying the EVA. Apart from the assumption of negligible inertial effects, we neglect for the time being the external body forces, so that Eq. (3) becomes orij ¼ 0: oxi The Cauchy strain tensor is given by 1 oui ouj eij ¼ þ 2 oxj oxi
ð4Þ
and the stress tensor rij in terms of the strain is given by rij ¼ 2Geij þ kekk dij
ð5Þ
substituting from Eq. (4) into (5), we obtain ! o2 ui o2 uj o2 uk G þ ¼ 0; þ k ox2j oxi oxj oxk oxi )
G
o2 ui o2 uj þ ðG þ kÞ ¼ 0: ox2j oxj oxi
ð6Þ
In vector form this becomes Gr2 u þ ðG þ kÞrðr uÞ ¼ 0
ð7Þ
which can be written equivalently in the form: ð2G þ kÞrðr uÞ Gr ðr uÞ ¼ 0:
ð8Þ
Since the cylinder is tethered (i.e. planar stress distribution), we assume (by symmetry) that the displacements are radial uðr; hÞ ¼ uðrÞ^er it is easy to show that Eq. (8) reduces to d 1 d ðruÞ ¼ 0: ð2G þ kÞ dr r dr
ð9Þ
ð10Þ
1102
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
Integrating twice with respect to r, we obtain c2 uðrÞ ¼ c1 r þ ; r where c1, c2 are arbitrary constants using the boundary conditions
ð11Þ
rrr ðaÞ ¼ pi ;
ð12Þ
rrr ðbÞ ¼ 0:
ð13Þ
It is straightforward to show that a2 p i 1þm c1 ¼ 2 ð1 2mÞ; E b a2 a2 b2 1þm p c2 ¼ 2 i E b a2 and therefore Eq. (11) can be written as a2 p 1 þ m b2 þ ð1 2mÞr : uðrÞ ¼ 2 i 2 E r b a
ð14Þ
Using Eqs. (4) and (9), it is clear that erh = ehr = 0 and that a2 p i 1þm b2 err ¼ 2 ð1 2mÞ 2 E r b a2 and
a2 p i 1þm b2 ehh ¼ 2 ð1 2mÞ þ 2 E r b a2
and using (5) in polar co-ordinates a2 p b2 rrr ¼ 2 i 2 1 2 ; r b a 2 a pi b2 1þ 2 : rhh ¼ 2 r b a2
ð15Þ
We now proceed to solve the viscoelastic problem (posed earlier) by using the Elastic-Viscoelastic analogy. We assume that the deviatoric effects are modelled by the ‘‘standard’’ solid while the dilatational effects are that of an isotropic elastic solid. Thus the coefficients of the time operators in Eq. (2) are given by 1 1 g P 1 : p0 ¼ þ ; ; p1 ¼ E E2 E1 E 2 ð16Þ g Q1 : q0 ¼ 1; q1 ¼ E2 for the standard solid, and by
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
P 2 : p 0
¼ 1;
Q2 : q 0
¼ 3K
1103
ð17Þ
for the isotropic elastic solid. The stress–strain relation of an isotropic elastic material can be expressed by two equations (separating the deviatoric sij, eij from dilatational rii, eii stresses and strains) as follows: sij ¼ 2Geij ;
ð18Þ
rii ¼ 3Keii ;
ð19Þ
where G is the shear modulus and K the bulk modulus. By considering the analogy between Eqs. (16), (17) and (18), (19) and also employing the superposition principle for both linearly elastic and linearly viscoelastic coefficient solids, we can establish the following relationships among the viscoelastic coefficients and the time operators: 1 þ m P1 ¼ ; E Q1 1 2m ¼
3Q1 P 2 : Q1 P 2 þ 2P 1 Q2
ð20Þ
ð21Þ
Using (20) and (21) in Eq. (14), u becomes a function of time, and we take the Laplace transform of (14) to obtain a2 p i P 1 b 2 3Q1 P 2 r þ uðr; sÞ ¼ 2 ; ð22Þ b a2 Q1 r Q1 P 2 þ 2P 1 Q2 where u, pi , P 1 , Q1 , P 2 , Q2 are the Laplace transforms of u, pi, P1, Q1, P2 and Q2, respectively. For the standard solid: 1 1 g P1 ¼ þ s; þ E 1 E2 E1 E 2 g Q1 ¼ 1 þ s E2 and for linear elastic solid P2 Q2
¼ ¼
1; 3K
substituting into Eq. (22), we obtain a2 1=g 1 b2 þ uðr; sÞ ¼ 2 pi ðsÞ s þ E2 =g E1 r b a2 ðE1 þ E2 Þ=g þ s r : þ3 1=gðE1 E2 þ 6KðE1 þ E2 ÞÞ þ ð6K þ E1 Þs
ð23Þ
1104
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
The viscoelastic solution is now found by taking the inverse Laplace transform of Eq. (23) and using the convolution theorem: 2Z t a2 b 1 E2 =gs 1 e uðr; tÞ ¼ 2 p ðt sÞ þ dðsÞ ds i g E1 b a2 r 0 Z t 3 a a=bs þ r pi ðt sÞ c þ dðsÞ ds ; ð24Þ e b 0 b where E1 E2 þ 6KE1 þ 6KE2 ; g E1 þ E2 b ¼ E1 þ 6K and c ¼ : g a¼
If the ventricular pressure is kept constant i.e. pi(t) = pi = const, then (24) becomes 2 a2 p i 1 1 b 1 c 1 E2 =gt a=bt þ 3r þ ð1 e uðr;tÞ ¼ 2 þ ð1 e Þ Þ b a b r b a2 E 1 E 2 ð25Þ and using the fact that err ¼
ou or
and
ehh ¼
u r
the viscoelastic strains are given by 2 E a2 p i 1 c 1 1 1 b bat g2 t ð1 e Þ err ðr;tÞ ¼ 2 3 þ þ ð1 e Þ 2 2 b a b E1 E2 r b a ð26Þ and a2 p ehh ðr;tÞ ¼ 2 i 2 b a
2 E2 1 1 b 1 c 1 bat t g þ ð1 e Þ 2 þ 3 þ 1 e : E1 E2 b a b r ð27Þ
From Eq. (15) note that no material constants appear in the expressions for rrr and rhh and thus the radial and tangential stresses in both the elastic and viscoelastic cases will be identical, i.e. a2 p i b2 rrr ¼ 2 1 2 ; r b a2 ð28Þ a2 p i b2 1þ 2 : rhh ¼ 2 r b a2
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1105
Once the strains have been calculated, it is a simple matter to compute the creep compliances. From (26), the radial creep compliance function J(t) is given by err/pi and we examine the limiting cases err ðr; 0Þ a2 3ð1 2mÞ 1 b2 þ J ð0 Þ ¼ ¼ 2 ; pi b a2 2E þ E1 ð1 2mÞ E1 r2 where E is the dilatational spring constant K¼
E 3ð1 2mÞ
also a2 c a 1 b2 2 ; b b2 g r2 b2 a2 a2 c 1 1 b2 3 þ J ð1Þ ¼ 2 a E1 E2 r2 b a2
J_ ð0þ Þ ¼
and J_ ð1Þ ¼ 0; where J_ ðtÞ signifies differentiation with respect to t. Typical qualitative behaviour of the creep function J(t) is shown in Fig. 2. From which it is evident that the point A is given by J(0+), the point B by the intercept of the t ! 1 asymptote with the vertical axes while J(0+) = tan /2 and J(1) = tan /1. Thus it appears that the four parameters of the standard solid E1, E2, g and k can be determined from the measurements of 0A, 0B, /1, and /2.
Fig. 2. Qualitative behaviour of creep function.
1106
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
4. Conclusion A simple model of the brain parenchyma, treated as an impermeable, viscoelastic solid, has been presented. This model treats the brain/ventricular cavity as a thick walled cylinder which is tethered at the ends and filled with CSF at a fixed pressure. Since this model yields explicit, analytical solutions, it can be used to provide benchmark solutions for numerical procedures capable of handling more complex, realistic geometries. A method of determining the mechanical parameters (for the standard solid), using it is creep compliance curve, has been discussed. However, initial time measurements are notoriously difficult to record accurately, due partly to inherent instrument error and partly to the fact that viscoelastic materials in general exhibit rapid initial creep rates (see [5]). Although this poses problems in the calculation of the intersection point A in Fig. 2, it can be obtained by interpolation techniques. Another approach to determining the parameters is to use dynamic-load experiments (i.e. oscillatory loading of the inner tube wall). An attractive feature of this type of experiment is the time scale involved: whereas creep and relaxation experiments may accurately provide data from 10 s to about 3 · 108 s the data from dynamic experiments extends over a time range from 103 s down to about 10 8 s (which is essentially the ‘‘initial time’’). From these, we may hope to compute explicitly the complex modulus and complex compliance from which the storage and loss modulus and compliance, and the loss angles may be extracted, and hence the mechanical parameters.
Acknowledgements This work was completed with support from the National Science and Engineering Research Council of Canada through a Collaborative Health Research Grant.
References [1] M.A. Biot, General theory of three dimensional consolidation, J. Appl. Phys. 12 (1941) 155– 164. [2] R.M. Christensen, Theory of Viscoelasticity: An Introduction, Academic Press, New York, 1982. [3] J.M. Golden, G.A.C. Graham, Boundary Value Problems in Linear Viscoelasticity, Springer Verlag, 1988. [4] J.E. Guinane, An equivalent circuit analysis of cerebrospinal fluid hydrodynamics, Am. J. Physiol. 223 (1972) 425–430. [5] M.E. Gurtin, E. Sternberg, On the linear theory of viscoelasticity, Arch. Ration. Mech. Anal. 11 (1962) 291–356.
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1107
[6] A. Marmarou, K. Shulman, R.M. Rosende, A non-linear analysis of the cerebrospinal fluid system and intracranial pressure dynamics, J. Neurosur. 49 (1978) 332–344. [7] J. Munro, Observations of the Structure and Functions of the Nervous System, Scotland, W. Greach et al., 1783. [8] T. Nagashima, N. Tanaki, S. Matsumoto, B. Horwitz, Y. Seguchi, Biomechanics of hydrocephalus: a new theoretical model, Neurosurgery 21 (1987) 898–904. [9] R.H. Simon, R.A.W. Lehman, S.J. OÕConnor, A comparison of pressure–volume models in hydrocephalus, Neurosurgery 15 (1984) 694. [10] S. Sivaloganathan, G. Tenti, J.M. Drake, Mathematical pressure–volume models of the cerebrospinal fluid, Appl. Math. Comput. 94 (1998) 243–266. [11] M. Stastna, G. Tenti, S. Sivaloganathan, J.M. Drake, Brain biomechanics: consolidation theory of hydrocephalus. Variable permeability and transient effects, Can. Appl. Math. Quart. 7 (1999) 111–124. [12] G. Tenti, S. Sivaloganathan, J.M. Drake, Brain biomechanics: steady-state consolidation theory of hydrocephalus, Can. Appl. Math. Quart. 7 (1999) 93–110.
www.elsevier.com/locate/amc
A viscoelastic approach to the modelling of hydrocephalus S. Sivaloganathan a
a,*
, M. Stastna a, G. Tenti a, J.M. Drake
b
Department of Applied Mathematics, University of Waterloo, Waterloo, Ont., Canada N2L 3G1 b Division of Neurosurgery, Hospital for Sick Children, University of Toronto, Toronto, Ont., Canada M5G 1X8
Abstract In earlier work [Appl. Math. Comput. 94 (1998) 243, Can. Appl. Math. Quart. 7 (1999) 111], we modelled the brain as a porous elastic medium containing an incompressible fluid and solved the simplified problem analytically using consolidation theory. In this paper we apply linear viscoelasticity theory to the modelling of the brain in the same simplified geometry. The aim in both approaches is to be able to eventually accurately predict the ventricular wall configuration in shunted hydrocephalus which would be of practical use to the neurosurgeon in locating the optimal position for shunt implantation in the treatment and management of hydrocephalus. We use the elastic– viscoelastic analogy (EVA) to obtain analytic solutions for the case when the brain is considered to be composed of a viscoelastic solid (the so-called ‘‘standard solid’’). We finally conclude with suggestions of how experimental data could be used to obtain realistic values for the mechanical parameters that appear in the constitutive relation for the viscoelastic solid. 2004 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail address: [email protected] (S. Sivaloganathan).
0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.06.037
1098
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1. Introduction The cerebral system is a complex continuum where a solid phase (brain tissue) interacts with fluid phases (essentially blood and cerebrospinal fluid) in different types of vessels and in the ventricles. A true mathematical description of the physics of such an interaction poses enormous challenges. In part due to this and the inherent complexity of the problem, early researchers replaced the brain with one or more units, or compartments, whose behaviour is represented by a single value of the pressure and by values of the flux exchanged with adjacent compartments. All these values may be time dependent but do not vary spatially. The resistance to the flow from one compartment to an adjacent one is lumped at the boundary between the two compartments. The genesis of this approach can be traced back over two centuries to Munro [7], since then numerous compartment models have appeared in the literature (e.g. see [4,6,9]). These have been recently reviewed by Sivaloganathan et al. [10] where the full mathematical solutions of the most relevant with a critical assessment of their validity was given. It was further demonstrated in [10] that all of these models can be derived as special cases of a general differential equation for the CSF (cerebrospinal fluid) pressure P(t). The weakness of these models, however, is their inability to describe the stress, strain distribution within the brain tissue, as well as their inability to say anything about the ventricular wall configuration in a hydrocephalic brain. Motivated in part by this, a new and more sophisticated generation of mathematical models has emerged over the past decade (see [8]). The starting point, for this approach, is the assumption that the brain parenchyma can be schematized as a sponge-like material. This approach relies heavily on the theory of consolidation introduced by Biot [1], and from a physical standpoint falls squarely into the category of flow through porous media. It is a continuum theory and leads to a set of PDEs which are usually solved by standard numerical techniques. However, Nagashima et al encountered serious problems in their numerical solution of the PDEs. In Tenti et al. [12] and Stastna et al. [11] we addressed these questions and elucidated the underpinnings of consolidation theory as applied to the brain. By solving the model analytically in a simplified geometry, we were able to link the numerical problems to the value of the Poisson ratio m used and provided a new and widely accepted estimate. In this paper we take yet a different approach and treat the brain parenchyma as a viscoelastic solid. We start again with a simplified picture of the brain, by assuming that the brain parenchyma is an impermeable, viscoelastic solid. In order to simplify the analysis as much as possible, we assume that the brain can be represented as a thick-walled cylindrical tube the inside of which is kept at a constant pressure pi (the ventricular CSF pressure) while the outside is kept at a lower constant pressure p0, which w.l.o.g. can be assumed to be zero. Although the geometry may appear rather simplistic, we are assured by neuro-
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1099
surgeons that in the case of hydrocephalus the shape of the lateral ventricular space filled with CSF is very approximately cylindrical. A natural question that arises, is what would be an appropriate viscoelastic model to use to represent the brain parenchyma. The basic Maxwell and Kelvin-Voigt models appear to be inappropriate; the former because it fails to capture the decreasing strain rate under constant stress characteristic of primary creep and the latter since it fails to represent time independent, instantaneous strain, or recovery of strain upon loading or unloading. Further neither of the models captures the very steep initial strain rates characteristic of many viscoelastic materials. In this paper, we propose the ‘‘standard solid’’ (see [2,3]) as an appropriate model. This is schematically represented in Fig. 1. Using Fig. 1, it is straightforward to derive the following constitutive equation for the standard solid: 1 1 g or g oe ¼eþ : ð1Þ þ rþ E1 E2 E1 E2 ot E2 ot Eq. (1) can be written in the general form: P r ¼ Qe; where P and Q are differential operators given by P¼
n X i¼0
Q¼
n X i¼0
pi
oi ; oti
oi qi i : ot
ð2Þ
Generally, viscoelastic stress analysis problems are more involved than elasticity problems due to the inclusion of the time variable in the DEs in addition to the space variables. However, in many problems (including the present one) where the types of BCs and temperature remain constant in time, the time variable can be removed by using Laplace transforms. We thus have an equivalent elastic problem. When a solution has been found in terms of the Laplace transform variable s, the inverse Laplace transform gives the required solution in the
Fig. 1. Schematic representation of standard solid.
1100
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
time variable t for the time dependent behaviour in the viscoelastic problem. This is, in essence, the elastic-viscoelastic analogy (EVA) or the ‘‘correspondence principle’’. For the theoretical development and mathematical justification of the EVA, the reader is referred to Gurtin and Sternberg [5].
2. Quasi-static approximation Inertia terms in the equations of motion are neglected, based on the assumption that the dependent variables vary slowly with respect to time so that the rate of change of momentum can be neglected. Thus we have the equilibrium equation orij þ F j ¼ 0; oxi
ð3Þ
where rij = rij(x, t) is the stress and F the body force. Due to the fact that stress relaxation and heat conduction waves are damped out very quickly in viscoelastic materials, disturbances (from boundary forces or changes in displacements) are transformed into static deformation and heat. Hence one assumes that the viscoelastic body is in equilibrium at every moment. This is the case when changes in b.c.s or changes due to stress relaxation take place slowly enough so that there is no significant change during the decay of the disturbance. An attractive feature in using the quasi-static approximation is that many corresponding elastic problems have been solved––the solution of the elastic problem can then be used to find the required viscoelastic solution. The boundary conditions can be defined by specifying both tractions and surface displacements. We divide the boundary C into two portions, CT where tractions are specified and Cu = CnCT where the displacements are specified, i.e.: T i ðx; tÞ
¼ rij ðx; tÞnj
ui ðx; tÞ
¼ fi ðx; tÞ on
on
CT ;
Cu :
However, the boundary C must be independent of t. It should be noted that if C varies with t, transform methods can still be used; however, Wiener–Hopf techniques will then be required and the great simplicity of the EVA is lost.
3. Mathematical formulation The assumption of a cylindrical geometry naturally suggests the use of cylindrical polar co-ordinates (r, h, z). The cylinder is assumed to be tethered, there are thus no deformations produced as a result of contact forces on its ends. As
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1101
a consequence, the absence of deformations in the z-direction implies the strain distribution will be planar (i.e. the same in any horizontal cross-section). The cylindrical annulus is assumed to be a homogeneous isotropic solid. The displacement and distribution of stresses and strains for this problem have been well known since the time of Lame´ (1852) (see [5]) and these results can be used to solve the problem, where the material is viscoelastic, through the ‘‘correspondence principle’’. Since the elastic problem is very important for the analysis of the viscoelastic problem, we review it in detail before applying the EVA. Apart from the assumption of negligible inertial effects, we neglect for the time being the external body forces, so that Eq. (3) becomes orij ¼ 0: oxi The Cauchy strain tensor is given by 1 oui ouj eij ¼ þ 2 oxj oxi
ð4Þ
and the stress tensor rij in terms of the strain is given by rij ¼ 2Geij þ kekk dij
ð5Þ
substituting from Eq. (4) into (5), we obtain ! o2 ui o2 uj o2 uk G þ ¼ 0; þ k ox2j oxi oxj oxk oxi )
G
o2 ui o2 uj þ ðG þ kÞ ¼ 0: ox2j oxj oxi
ð6Þ
In vector form this becomes Gr2 u þ ðG þ kÞrðr uÞ ¼ 0
ð7Þ
which can be written equivalently in the form: ð2G þ kÞrðr uÞ Gr ðr uÞ ¼ 0:
ð8Þ
Since the cylinder is tethered (i.e. planar stress distribution), we assume (by symmetry) that the displacements are radial uðr; hÞ ¼ uðrÞ^er it is easy to show that Eq. (8) reduces to d 1 d ðruÞ ¼ 0: ð2G þ kÞ dr r dr
ð9Þ
ð10Þ
1102
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
Integrating twice with respect to r, we obtain c2 uðrÞ ¼ c1 r þ ; r where c1, c2 are arbitrary constants using the boundary conditions
ð11Þ
rrr ðaÞ ¼ pi ;
ð12Þ
rrr ðbÞ ¼ 0:
ð13Þ
It is straightforward to show that a2 p i 1þm c1 ¼ 2 ð1 2mÞ; E b a2 a2 b2 1þm p c2 ¼ 2 i E b a2 and therefore Eq. (11) can be written as a2 p 1 þ m b2 þ ð1 2mÞr : uðrÞ ¼ 2 i 2 E r b a
ð14Þ
Using Eqs. (4) and (9), it is clear that erh = ehr = 0 and that a2 p i 1þm b2 err ¼ 2 ð1 2mÞ 2 E r b a2 and
a2 p i 1þm b2 ehh ¼ 2 ð1 2mÞ þ 2 E r b a2
and using (5) in polar co-ordinates a2 p b2 rrr ¼ 2 i 2 1 2 ; r b a 2 a pi b2 1þ 2 : rhh ¼ 2 r b a2
ð15Þ
We now proceed to solve the viscoelastic problem (posed earlier) by using the Elastic-Viscoelastic analogy. We assume that the deviatoric effects are modelled by the ‘‘standard’’ solid while the dilatational effects are that of an isotropic elastic solid. Thus the coefficients of the time operators in Eq. (2) are given by 1 1 g P 1 : p0 ¼ þ ; ; p1 ¼ E E2 E1 E 2 ð16Þ g Q1 : q0 ¼ 1; q1 ¼ E2 for the standard solid, and by
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
P 2 : p 0
¼ 1;
Q2 : q 0
¼ 3K
1103
ð17Þ
for the isotropic elastic solid. The stress–strain relation of an isotropic elastic material can be expressed by two equations (separating the deviatoric sij, eij from dilatational rii, eii stresses and strains) as follows: sij ¼ 2Geij ;
ð18Þ
rii ¼ 3Keii ;
ð19Þ
where G is the shear modulus and K the bulk modulus. By considering the analogy between Eqs. (16), (17) and (18), (19) and also employing the superposition principle for both linearly elastic and linearly viscoelastic coefficient solids, we can establish the following relationships among the viscoelastic coefficients and the time operators: 1 þ m P1 ¼ ; E Q1 1 2m ¼
3Q1 P 2 : Q1 P 2 þ 2P 1 Q2
ð20Þ
ð21Þ
Using (20) and (21) in Eq. (14), u becomes a function of time, and we take the Laplace transform of (14) to obtain a2 p i P 1 b 2 3Q1 P 2 r þ uðr; sÞ ¼ 2 ; ð22Þ b a2 Q1 r Q1 P 2 þ 2P 1 Q2 where u, pi , P 1 , Q1 , P 2 , Q2 are the Laplace transforms of u, pi, P1, Q1, P2 and Q2, respectively. For the standard solid: 1 1 g P1 ¼ þ s; þ E 1 E2 E1 E 2 g Q1 ¼ 1 þ s E2 and for linear elastic solid P2 Q2
¼ ¼
1; 3K
substituting into Eq. (22), we obtain a2 1=g 1 b2 þ uðr; sÞ ¼ 2 pi ðsÞ s þ E2 =g E1 r b a2 ðE1 þ E2 Þ=g þ s r : þ3 1=gðE1 E2 þ 6KðE1 þ E2 ÞÞ þ ð6K þ E1 Þs
ð23Þ
1104
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
The viscoelastic solution is now found by taking the inverse Laplace transform of Eq. (23) and using the convolution theorem: 2Z t a2 b 1 E2 =gs 1 e uðr; tÞ ¼ 2 p ðt sÞ þ dðsÞ ds i g E1 b a2 r 0 Z t 3 a a=bs þ r pi ðt sÞ c þ dðsÞ ds ; ð24Þ e b 0 b where E1 E2 þ 6KE1 þ 6KE2 ; g E1 þ E2 b ¼ E1 þ 6K and c ¼ : g a¼
If the ventricular pressure is kept constant i.e. pi(t) = pi = const, then (24) becomes 2 a2 p i 1 1 b 1 c 1 E2 =gt a=bt þ 3r þ ð1 e uðr;tÞ ¼ 2 þ ð1 e Þ Þ b a b r b a2 E 1 E 2 ð25Þ and using the fact that err ¼
ou or
and
ehh ¼
u r
the viscoelastic strains are given by 2 E a2 p i 1 c 1 1 1 b bat g2 t ð1 e Þ err ðr;tÞ ¼ 2 3 þ þ ð1 e Þ 2 2 b a b E1 E2 r b a ð26Þ and a2 p ehh ðr;tÞ ¼ 2 i 2 b a
2 E2 1 1 b 1 c 1 bat t g þ ð1 e Þ 2 þ 3 þ 1 e : E1 E2 b a b r ð27Þ
From Eq. (15) note that no material constants appear in the expressions for rrr and rhh and thus the radial and tangential stresses in both the elastic and viscoelastic cases will be identical, i.e. a2 p i b2 rrr ¼ 2 1 2 ; r b a2 ð28Þ a2 p i b2 1þ 2 : rhh ¼ 2 r b a2
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1105
Once the strains have been calculated, it is a simple matter to compute the creep compliances. From (26), the radial creep compliance function J(t) is given by err/pi and we examine the limiting cases err ðr; 0Þ a2 3ð1 2mÞ 1 b2 þ J ð0 Þ ¼ ¼ 2 ; pi b a2 2E þ E1 ð1 2mÞ E1 r2 where E is the dilatational spring constant K¼
E 3ð1 2mÞ
also a2 c a 1 b2 2 ; b b2 g r2 b2 a2 a2 c 1 1 b2 3 þ J ð1Þ ¼ 2 a E1 E2 r2 b a2
J_ ð0þ Þ ¼
and J_ ð1Þ ¼ 0; where J_ ðtÞ signifies differentiation with respect to t. Typical qualitative behaviour of the creep function J(t) is shown in Fig. 2. From which it is evident that the point A is given by J(0+), the point B by the intercept of the t ! 1 asymptote with the vertical axes while J(0+) = tan /2 and J(1) = tan /1. Thus it appears that the four parameters of the standard solid E1, E2, g and k can be determined from the measurements of 0A, 0B, /1, and /2.
Fig. 2. Qualitative behaviour of creep function.
1106
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
4. Conclusion A simple model of the brain parenchyma, treated as an impermeable, viscoelastic solid, has been presented. This model treats the brain/ventricular cavity as a thick walled cylinder which is tethered at the ends and filled with CSF at a fixed pressure. Since this model yields explicit, analytical solutions, it can be used to provide benchmark solutions for numerical procedures capable of handling more complex, realistic geometries. A method of determining the mechanical parameters (for the standard solid), using it is creep compliance curve, has been discussed. However, initial time measurements are notoriously difficult to record accurately, due partly to inherent instrument error and partly to the fact that viscoelastic materials in general exhibit rapid initial creep rates (see [5]). Although this poses problems in the calculation of the intersection point A in Fig. 2, it can be obtained by interpolation techniques. Another approach to determining the parameters is to use dynamic-load experiments (i.e. oscillatory loading of the inner tube wall). An attractive feature of this type of experiment is the time scale involved: whereas creep and relaxation experiments may accurately provide data from 10 s to about 3 · 108 s the data from dynamic experiments extends over a time range from 103 s down to about 10 8 s (which is essentially the ‘‘initial time’’). From these, we may hope to compute explicitly the complex modulus and complex compliance from which the storage and loss modulus and compliance, and the loss angles may be extracted, and hence the mechanical parameters.
Acknowledgements This work was completed with support from the National Science and Engineering Research Council of Canada through a Collaborative Health Research Grant.
References [1] M.A. Biot, General theory of three dimensional consolidation, J. Appl. Phys. 12 (1941) 155– 164. [2] R.M. Christensen, Theory of Viscoelasticity: An Introduction, Academic Press, New York, 1982. [3] J.M. Golden, G.A.C. Graham, Boundary Value Problems in Linear Viscoelasticity, Springer Verlag, 1988. [4] J.E. Guinane, An equivalent circuit analysis of cerebrospinal fluid hydrodynamics, Am. J. Physiol. 223 (1972) 425–430. [5] M.E. Gurtin, E. Sternberg, On the linear theory of viscoelasticity, Arch. Ration. Mech. Anal. 11 (1962) 291–356.
S. Sivaloganathan et al. / Appl. Math. Comput. 163 (2005) 1097–1107
1107
[6] A. Marmarou, K. Shulman, R.M. Rosende, A non-linear analysis of the cerebrospinal fluid system and intracranial pressure dynamics, J. Neurosur. 49 (1978) 332–344. [7] J. Munro, Observations of the Structure and Functions of the Nervous System, Scotland, W. Greach et al., 1783. [8] T. Nagashima, N. Tanaki, S. Matsumoto, B. Horwitz, Y. Seguchi, Biomechanics of hydrocephalus: a new theoretical model, Neurosurgery 21 (1987) 898–904. [9] R.H. Simon, R.A.W. Lehman, S.J. OÕConnor, A comparison of pressure–volume models in hydrocephalus, Neurosurgery 15 (1984) 694. [10] S. Sivaloganathan, G. Tenti, J.M. Drake, Mathematical pressure–volume models of the cerebrospinal fluid, Appl. Math. Comput. 94 (1998) 243–266. [11] M. Stastna, G. Tenti, S. Sivaloganathan, J.M. Drake, Brain biomechanics: consolidation theory of hydrocephalus. Variable permeability and transient effects, Can. Appl. Math. Quart. 7 (1999) 111–124. [12] G. Tenti, S. Sivaloganathan, J.M. Drake, Brain biomechanics: steady-state consolidation theory of hydrocephalus, Can. Appl. Math. Quart. 7 (1999) 93–110.