On some model equations for pulsatile flow in viscoelastic vessels
aa r X i v : . [ phy s i c s . f l u - dyn ] M a y ON SOME MODEL EQUATIONS FOR PULSATILE FLOW INVISCOELASTIC VESSELS
DIMITRIOS MITSOTAKIS, DENYS DUTYKH, QIAN LI, AND ELIJAH PEACH
Abstract.
Considered here is the derivation of partial differential equationsarising in pulsatile flow in pipes with viscoelastic walls. The equations areasymptotic models describing the propagation of long-crested pulses in pipeswith cylindrical symmetry. Additional effects due to viscous stresses in bio-fluids are also taken into account. The effects of viscoelasticity of the vesselson the propagation of solitary and periodic waves in a vessel of constant radiusare being explored numerically. Introduction
The description of fluid flows in pipes with viscoelastic wall material is motivatedmainly by the studies of hemodynamics [16]. A cardiac cycle consists of the systolicphase where the heart ventricles contract and pump blood to the arteries and thediastolic phase where the heart ventricles are relaxed and the heart fills with bloodagain. During the systolic phase the large arteries are deformed and store elasticenergy that is released during the diastolic phase. This property of the vesselsis usually referred to as the compliance of the vessels. Modelling the viscoelasticproperties of the vessels appears to have significant difficulties of mathematical andnumerical nature [16, 8, 27, 34, 23, 30]. The mathematical modelling of such flowssuggests the use of the equations of continuum mechanics for incompressible fluidflow known as the Navier-Stokes equations.The Navier-Stokes equations in three dimensions are too complicated to be usedin practical situations and for this reason several simplified mathematical modelshave been derived [37, 26, 33, 14, 3, 13]. The main simplifications that have beenmade are based on the assumption of the axial (or cylindrical) symmetry of thevessels. This assumption and using approximations of the averaged velocity ofthe fluid led to the derivation of simple one-dimensional (1+1D) models [17, 39,32, 31, 24, 1, 28, 35]. The models include unidirectional [13, 38, 21, 11, 6], andbidirectional models [4, 22, 7]. Bidirectional models can approximate accuratelyreflections of pulses occurred in the presence of non-uniformities in the vessels wallas opposed to unidirectional models such as the KdV equation ([38, 10]) where it isassumed that the pulses propagate mainly in one direction. In this paper we focuson the derivation of bidirectional equations. Further simplifications were madeby assuming that the velocity of the fluid is very large. This assumption led tolumped-parameter (or zero-dimensional) models [2, 19] and the references therein.
Date : May 30, 2019.2000
Mathematics Subject Classification.
Key words and phrases.
Boussinesq systems, blood flow, viscoelastic vessel.
Figure 1.
Sketch of the physical domain for a single vessel segmentwith elastic and impenetrable wall.
The later models are used in practice to predict the flow and the pressure of theblood in operational situations [34].For practical reasons, the inclusion of the dissipative effects in the flow can bedone by assuming a laminar flow and small viscosity. For example assuming aparabolic profile for the horizontal velocity of the fluid it has been shown thatthe Navier-Stokes equations can be reduced to a modified system which is verysimilar to the Euler equations [30], (we also refer to the Poiseuille solution forthe justification of this parabolic profile of the horizontal velocity). Specifically,denoting u = u ( x, r, t ), v = v ( x, r, t ) the horizontal and radial velocity respectively,and u w ( x, t ) = u ( x, r w , t ) the horizontal velocity of the fluid on the vessel wall(at radius r = r w ( x, t )), then assuming that u ( x, r, t ) is proportional to (( r w ) − r ) u w ( x, t ) these equations can be written in cylindrical coordinates in the form: u t + uu x + vu r + 1 ρ p x + κu w = 0 , (1.1) v t + uv x + vv r + 1 ρ p r = 0 , (1.2) u x + v r + 1 r v = 0 , (1.3)where p = p ( x, r, t ) is the pressure of the fluid, ρ is the constant density of the fluidand κ is the viscous frequency parameter (also known as the Rayleigh dampingcoefficient) with dimensions [ s − ].A sketch of the physical domain for this problem is presented in Fig. 1, wherethe distance of vessel’s wall from the centre of the vessel in a cross section is denotedby r w ( x, t ) and depends on x and t while the radius of the vessel at rest is givenby the function r ( x ). In general the deformation of the wall will be a function of x and t . If we denote the radial displacement of the wall by η ( x, t ) then the vesselwall radius can be written as r w ( x, t ) = r ( x ) + η ( x, t ).The governing equations (1.1)–(1.3) combined with initial and boundary condi-tions form a closed system. A compatibility condition is also applied at the centreof the vessel (due to cylindrical symmetry). Specifically, we assume that(1.4) v ( x, r, t ) = 0 , for r = 0 . In general we assume for consistency purposes that v ( x, r, t ) = O ( r τ ), with τ ≥ r →
0. The impermeability of the vessel wall can be described by the equation:(1.5) v ( x, r, t ) = η t ( x, t ) + ( r ( x ) + η ( x, t )) x u ( x, r, t ) , for r = r w ( x, t ) , and expresses that the fluid velocity equals the wall speed v = r wt . The system isaccompanied also by a second boundary condition, which is the Newton second law N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 3 applied on the vessel wall:(1.6) ρ w hη tt ( x, t ) = p w ( x, t ) − E σ hr ( x ) ( η ( x, t ) + γη t ( x, t )) , where ρ w is the wall density, p w is the transmural pressure, h is the thickness ofthe vessel wall, E σ = E/ (1 − σ ) where E is the Young modulus of elasticity with σ denoting the Poisson ratio of the viscoelastic wall. In this study we assumethat E is a constant and in general we simplify the notation by denoting E σ with E . The last term in (1.6) models the viscous nature of the vessel wall and canbe derived by using a simple Kelvin-Voigt model (spring-dashpot model). In thissetting γ = n w /E v where n w is the dashpot coefficient of viscosity and E v is theYoung modulus of the viscous part of Kelvin body. In practical situations theparameter γ is very small and usually can be taken to be of order O (10 − ). It isnoted that because the flow is pressure-driven the effect of gravity is neglected. Formore information about the derivation of the Euler equations and the boundaryconditions we refer to [40, 8].Although the dispersion of the flow can be ignored from the majority of math-ematical models derived from the Navier-Stokes equations or from the Euler equa-tions resulting into very simple systems of conservation laws, the need for moreaccurate description of the waves and their reflections suggests the inclusion ofthis fundamental property. A first attempt towards the derivation of bi-directionalweakly-nonlinear and weakly-dispersive system of equations was presented in [7].Using asymptotic techniques more general asymptotic models were derived in [20].The systems derived in [20] appeared to justify the non-dispersive models of [15, 29]with asymptotic reasoning. It was also shown that the inclusion of dispersive termscan describe more accurately the effects of the vessel wall variations within the flow.One basic ingredient that was ignored in both works [7, 20] is the viscosity effectsof the vessels by assuming simple elastic vessels.In this paper we extend the work [20] and derive some new asymptotic one-dimensional equations of Boussinesq type (weakly non-linear and weakly dispersive)that approximate the system (1.1)–(1.3) with boundary conditions (1.4)–(1.6). Thederivation is based on formal expansions of the velocity potential as in [12]. Thenew systems generalise the previously derived Boussinesq systems of [20] as theycoincide with them when the viscoelastic property of the vessel wall is ignored. Thenew mathematical models are of significant importance since they include all thenecessary ingredients for the accurate description of regular fluid flows in pipes withviscoelastic walls.Since the dissipation caused by the fluid viscosity and the dissipation caused bythe viscoelasticity of the vessel wall are different in their nature there is a question onwhether the different dissipative terms have also different effects on the propagationof pulses. The answer to this question is explored computationally by studying theeffects caused by the two dissipative terms on the propagation of a solitary wave.The paper is organised as follows: In Section 2 we present the derivation of a newsystem of Boussinesq type for the description of the velocity and the deviation ofthe vessel wall for fluid flow in a viscoelastic vessel. This system is further improvedin Section 3 by computing the fluid velocity at different levels of radius. In Section4 further simplifications lead to unidirectional equations that depend only on thedeviation of the vessel wall, while the velocity of the fluid can be computed explicitlyusing a simple asymptotic formula. Section 5 demonstrates the dissipation effects DIMITRIOS MITSOTAKIS, DENYS DUTYKH, QIAN LI, AND ELIJAH PEACH on the propagation of solitary and periodic waves. We close this paper with someconclusions and perspectives.2.
Derivation of the new mathematical models
Here we proceed with the derivation of the new equations. The derivation isbased on the assumption that the flow is irrotational and therefore we assume theexistence of a smooth velocity potential φ ( x, r, t ) such that ( u, v ) T = ∇ φ , i.e. weassume that u = φ x and v = φ r . Then, as in [36] the velocity potential can bechosen appropriately such that the equations (1.1)–(1.2) can be integrated into thegeneralised Cauchy-Lagrange integral:(2.1) φ t + 12 φ x + 12 φ r + 1 ρ p + κφ = 0 , for r = r w . The mass conservation (continuity) equation is then reduced to the elliptic equation(2.2) rφ xx + ( rφ r ) r = 0 , < r < r w , and boundary conditions for the velocity are written as(2.3) φ r = 0 , for r = 0 , and(2.4) φ r = η t + ( r ( x ) + η ) x φ x , for r = r w ( x, t ) . In order to make simplifications to the previous equations we consider the fol-lowing non-dimensional (scaled) variables:(2.5) η ∗ = ηa , x ∗ = xλ , r ∗ = rR , t ∗ = tT , φ ∗ = 1 λε ˜ c φ, p ∗ = 1 ερ ˜ c p , where a is a typical deviation of the vessel wall from its rest position, λ a typicalwavelength of a pulse, R is a vessel’s typical radius, T = λ/ ˜ c the characteristictime scale, while ˜ c = p Eh/ ρR is the Moens-Korteweg characteristic speed, [16].It is noted that the external pressure is considered zero and is neglected. Theparameters ε and δ characterise the nonlinearity and the dispersion of the system:(2.6) ε = aR , δ = Rλ .
Usually, ε and δ are very small. Specifically, we assume that ε ≪ δ ≪
1, whilethe Stokes-Ursell number is of order 1: ε/δ = O (1). The system of equations(2.1)–(2.5) along with the boundary condition (1.6) is then written in dimensionlessvariables in the form: φ ∗ t ∗ + ε φ ∗ x ∗ + εδ φ ∗ r ∗ + p ∗ + εκ ∗ φ ∗ = 0 , for r ∗ = r ∗ w , (2.7) δ r ∗ φ ∗ x ∗ x ∗ + ( r ∗ φ ∗ r ∗ ) r ∗ = 0 , < r ∗ < r ∗ w , (2.8) φ ∗ r ∗ = 0 , for r ∗ = 0 , (2.9) φ ∗ r ∗ = δ η ∗ t ∗ + δ ( r ∗ ( x ∗ ) + εη ∗ ) x ∗ φ ∗ x ∗ , for r ∗ = r ∗ w , (2.10) p ∗ = α ∗ δ η ∗ t ∗ t ∗ + β ∗ ( η ∗ + δ γ ∗ η ∗ t ∗ ) , for r ∗ = r ∗ w , (2.11)where κ ∗ = κλ/ ˜ cε , α ∗ = ρ w h/ρR , β ∗ = β ∗ ( x ) = 2 R /r ( x ) and γ ∗ = γ/δ T . Forthe sake of simplicity in the notation, we drop the asterisk from the new variablesin the following derivations for the non-dimensional variables except if it is statedotherwise. N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 5
Following standard asymptotic techniques, cf. [5], we consider a formal expansionof the velocity potential [18]:(2.12) φ ( x, r, t ) = ∞ X m =0 r m φ m ( x, t ) . Demanding φ to satisfy the equation (2.8) leads to the following recurrence relation(2.13) δ ∂ x φ m + (2 m + 2) φ m +2 = 0 , φ m +1 = 0 , for m = 0 , , , · · · , where ∂ jx denotes the j -th order derivative with respect to x . Adirect application of the last relation is(2.14) φ = − δ ∂ x φ , and(2.15) φ m +2 = δ (2 m + 2) (2 m ) ∂ x φ m − = O ( δ ) , for m = 1 , , · · · . The last relation ensures that the terms φ m of the velocitypotential expansion for m ≥ φ m = ( − m δ m m ∂ mx φ , for m = 1 , , · · · , and therefore(2.17) φ ( x, r, t ) = ∞ X m =0 r m ( − m δ m m ∂ mx φ . A 2nd order asymptotic approximation of the velocity potential is(2.18) φ ( x, r, t ) = φ ( x, t ) − δ r ∂ xx φ ( x, t ) + O ( δ ) . Using the previous observations on the expansion of the velocity potential weobserve that (2.10) can be approximated by the relation(2.19) η t + r wx φ x + r w φ xx − δ r r x φ xxx − δ r φ xxxx = O ( δ , εδ ) . Since the momentum balance laws were reduced to the Cauchy-Lagrange integralequation (2.7) for r = r w we can eliminate the pressure using (2.11) and obtain(2.20) φ t + ε φ x + εδ φ r + εκφ + αδ η tt + β ( η + γη t ) = 0 . Substituting (2.18) into (2.20) we obtain the approximate momentum equation(2.21) φ t − δ r φ xxt + ε φ x + εκφ + αδ η tt + β ( η + δ γη t ) = O ( δ , εδ ) . Denoting the horizontal velocity at the centre of the vessel u ( x, , t ) = φ x ( x, t )by w ( x, t ) we rewrite the equations (2.19)–(2.21) in the following form: η t + r wx w + r w w x − δ r r x w xx − δ r w xxx = O ( δ , εδ ) , (2.22) w t + ( βη ) x + εww x + εκw − δ r r x w xt − δ r w xxt + αδ η xtt + δ γ ( βη t ) x = O ( δ , εδ ) . (2.23) DIMITRIOS MITSOTAKIS, DENYS DUTYKH, QIAN LI, AND ELIJAH PEACH
Although the system (2.22)–(2.23) is a valid approximation of the Euler equationsfor the prescribed asymptotic regime, it is not of much practical use due to thetemporal derivatives of the wall deviations in the momentum equation (2.23). Forthis reason we proceed with further simplifications using low-order approximationfor the velocity w and the deviation of the vessel wall.From the equations (2.22)–(2.23) we observe that(2.24) w t = − [ βη ] x + O ( ε, δ ) , η t = − r wx w − r w w x + O ( ε, δ ) . Substituting these low-order approximations into (2.23) we obtain the simplifiedmomentum equation[1 − δ αr xx ] w t + [ βη ] x + δ (3 α + r ) r x βη ] xx + εww x − δ (2 α + r ) r w xxt ++ εκw − δ γ [ β ( r x w + r w x )] x = O ( δ , εδ ) . (2.25)The system (2.22)–(2.25) can be the base to other more amenable Boussinesqsystems along the lines of [5]. In the next section we derive a simplified Boussinesqsystems with favourable properties in analogy to the classical Boussinesq systemfor fluid flow in purely elastic vessels derived in [20].3. The classical Boussinesq system
In this section we proceed with some improvements on system (2.22)-(2.25) basedon the evaluation of the horizontal velocity at any level of radius r . Specifically,from (2.18) we have that u ( x, r, t ) = w ( x, t ) − δ r w xx ( x, t ) + O ( δ ). From this weobserve that considering the fluid velocity at any radius r = θr w for 0 ≤ θ ≤ u θ ( x, t ) . = u ( x, θr w , t ) = w ( x, t ) − δ θ r w ( x, t ) + O ( εδ , δ ) and thus(3.1) w ( x, t ) = u θ ( x, t ) + δ θ r u θxx ( x, t ) + O ( εδ , δ ) . Substitution of (3.1) into (2.22) and (2.25) leads to the more general Boussinesqsystem:(3.2) η t + r wx u θ + r w u θx + δ (2 θ − r r x u θxx + δ (2 θ − r u θxxx = O ( εδ , δ ) . [1 − δ αr xx ] u θt + [ βη ] x + δ (3 α + r ) r x βη ] xx + εu θ u θx − δ [2 α − ( θ − r ] r u θxxt ++ εκu θ − δ γ [ β ( r x u θ + r u θx )] x = O ( δ , εδ ) . (3.3)If we take γ = 0 in the system (3.2)–(3.3) then we recover the Boussinesq systemof KdV-BBM type of [20].The linear dispersion properties of the system (3.2)–(3.3) depend on the choiceof the parameter θ . Taking θ = 1 / θ = 1 / N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 7 (3.2)–(3.3) is simplified to the classical Boussinesq system for viscoelastic vessels: η t + 12 ( r + εη ) u x + ( r + εη ) x u = O ( εδ , δ ) , (3.4) [1 − δ αr xx ] u t + [ βη ] x + εuu x − δ (4 α + r ) r u xxt +(3.5) + δ (3 α + r ) r x βη ] xx + εκu − δ γ [ β ( r x u + r u x )] x = O ( εδ , δ ) , where we drop the θ in this notation by taking u = u θ for θ = 1 /
2. This systemis very similar to the Peregrine system of water wave theory and we usually callit the classical Boussinesq system [25]. The classical Boussinesq system (3.4)–(3.5)after discarding the negligible terms on the right side, can be written in dimensionalvariables form η t + 12 ( r + η ) u x + ( r + η ) x u = 0 , (3.6) [1 − ¯ αr xx ] u t + [ ¯ βη ] x + uu x − (4 ¯ α + r ) r u xxt +(3.7) + (3 ¯ α + r ) r x βη ] xx + κu − γ [ ¯ β ( r x u + r u x )] x = 0 , where ¯ α = ρ w hρ and ¯ β ( x ) = Ehρr ( x ) , while all the variables x , t , r ( x ), η ( x, t ) and u ( x, t ) are in dimensional form.We underline that the new dissipative terms κu and − γ [ ¯ β ( r x u + r u x )] x are to-tally different in nature and in mathematical properties. Their coefficients appearedin these formulas are also important for the derivation of Windkessel (0D) models[23] as they are the result of asymptotic simplifications to the Euler equations andcan be specified by the wall material properties.System (3.8) – (3.9) extends the classical Boussinesq system derived and analysedin [20] in the case of vessels with elastic walls. Using low order corrections to thedispersive terms one can extend the whole class of the Boussinesq systems derivedin [20]. The extended systems will differ only by the additional viscoelastic termand because of their complexity we do not proceed with their derivation here.Further simplifications can be achieved in the system (3.8)–(3.9) by assumingthat the undisturbed radius of the vessel is constant. The resulting system takesthe form η t + 12 r u x + 12 ηu x + η x u = 0 , (3.8) u t + ¯ βη x + uu x − (4 ¯ α + r ) r u xxt + κu − γ ¯ β r u xx = 0 , (3.9)where ¯ α and ¯ β are constants.The dispersion relation ω = ω ( k ) of the derived system can be easily computed.It is given by the following quadratic algebraic equation:(3.10) (cid:16) α r k + r k (cid:17) ω + i (cid:16) κ + γ ¯ β r k (cid:17) ω − ¯ β r k . Here k is the wavenumber and ω = ω ( k ) is the corresponding wave frequency.The dispersion relation (3.10) can be solved explicitly for ω using the standard DIMITRIOS MITSOTAKIS, DENYS DUTYKH, QIAN LI, AND ELIJAH PEACH formulas for the roots of a quadratic equation. It will contain two branches whoseexpressions we omit here for brevity.3.1.
Solitary wave solutions.
In the case where all the dissipative terms are ne-glected (i.e. when κ = γ = 0) the system possesses classical solitary wave solutionsthat satisfy the following speed-amplitude relation [20](3.11) s = s βr (2 − ζ ) − √ − ζζ (3 − ζ ) , where ζ = 1 − r / ( a + r ) , and s denotes the speed and a the amplitude of thesolitary wave. The graph of the speed-amplitude relationship for values of theparameters ¯ α and ¯ β for a typical blood vessels of different radius r is presentedin Figure 2. We observe that there is as the amplitude of the solitary waves growstheir speed tend to approach an upper bound. The upper bound for the solitarywave speed can be computed by taking the limit a → ∞ in (3.11), and which is p br . Figure 2.
Speed-amplitude relation (3.11) for solitary waves of theclassical Boussinesq system for various values of r . For the numerical generation of these solitary waves we refer to [20]. Other moregeneral systems similar to the systems derived in [20] can also be derived but sincethey only differ on the viscoelastic term we omit their derivation here and we referto [20] for more information.3.2.
Symmetries and conservation laws.
In the presence of any dissipativeterm (i.e. κ + γ = 0), the system (3.8), (3.9) possesses only two point sym-metries: translations in time and in space. If both dissipative terms are neglected(i.e. κ + γ = 0), then we gain an additional point symmetry transformationwhose infinitesimal generator is given by X = t ∂∂t − η + r ) ∂∂η − u ∂∂u . The corresponding symmetry transformation can be readily computed: x ← x ,t ← e c t ,η ← ( r + η ) e − c − r ,u ← e − c u , N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 9 where c is an arbitrary constant.For arbitrary κ and γ system (3.8), (3.9) can be written in an elegant conservativeform: (cid:0) η ( η + 2 r ) (cid:1) t + (cid:2) ( r + η ) u (cid:3) x = 0 , (cid:16)(cid:0) u − r ( r + 4 ¯ α ) u x x (cid:1) e κ t (cid:17) t + h (cid:0) u + 24 ¯ β η + r κ u x + 4 r ¯ α κ u x − r u t x − r α u t x − γ ¯ β u x (cid:1) e κ t i x = 0 . For instance, these conserved quantities can be used in numerical simulations tocheck method accuracy. 4.
Unidirectional models
In this section we present unidirectional models, namely the BBM and KdVequations by considering waves that propagate mainly in one direction. In order toderive such models we consider the following dimensionless variables:(4.1) η ∗ = ηa , x ∗ = xλ , r ∗ = rr , t ∗ = tT , u ∗ = uc , where here c = ar q Ehρr is a modified Moens-Korteweg characteristic speed and T = 2 aλr c . The system (3.8)–(3.9) in dimensionless variables then is written: η ∗ t ∗ + u ∗ x ∗ + εη ∗ u ∗ x ∗ + 2 εη ∗ x ∗ u ∗ = 0 , (4.2) u ∗ t ∗ + η ∗ x ∗ + u ∗ u ∗ x ∗ − δ α ∗ + 18 u ∗ x ∗ x ∗ t ∗ + εκ ∗ u ∗ − δ γ ∗ u xx = 0 , (4.3)where here κ ∗ = 2 κσ , γ ∗ = σεδ γ , σ = c λ and ε = ar , δ = r λ . Consideringa flow mainly towards one direction, we can use the low-order approximation forunidirectional wave propagation [36](4.4) u ∗ = η ∗ + εF + δ G + O ( ε , εδ , δ ) , where F and G are unknown functions of x , t . Substitution of (4.4) into (4.2)–(4.3)gives the equations: η ∗ t ∗ + η ∗ x ∗ + ε ( F x ∗ + 3 η ∗ η ∗ x ∗ ) + δ G x = O ( ε , εδ , δ ) , (4.5) η ∗ t ∗ + η ∗ x ∗ + ε ( F t ∗ + 2 η ∗ η ∗ x ∗ + k ∗ η ∗ ) + δ ( G t ∗ − α ∗ + 18 η ∗ x ∗ x ∗ t ∗ − γ ∗ η ∗ x ∗ x ∗ ) = O ( ε , εδ , δ ) . (4.6)Choosing appropriate functions F and G (4.7) F = − η ∗ + κ ∗ Z η ∗ dx ∗ , G = − α ∗ + 116 η ∗ x ∗ t ∗ − γ ∗ η ∗ x ∗ , equations (4.5) and (4.6) coincide up to the order O ( ε, δ ) to a single equation for η ∗ , namely, the dimensionless BBM equation:(4.8) η ∗ t ∗ + η ∗ x ∗ + ε η ∗ η ∗ x ∗ − δ α ∗ + 116 η ∗ x ∗ x ∗ t ∗ + ε κ ∗ η − δ γ ∗ η ∗ x ∗ x ∗ = 0 . In dimensional variables (4.8) takes the form(4.9) η t + ˜ cη x + 52 1 r ˜ cηη x − (4 ¯ α + r ) r η xxt + κ η − ˜ c γη xx = 0 , where here ˜ c = q Eh ρr is the standard Moens-Korteweg characteristic speed. Thedispersion relation ω = ω ( k ) can be easily computed:(4.10) ω = ˜ c k r (4 ¯ α + r )16 k − i2 κ + γ ˜ c k r (4 ¯ α + r )16 k . In the absence of any form of dissipation, it is known that the BBM equationpossesses classical solitary waves propagating with speed c s given by the formula,[20],(4.11) η ( x, t ) = 3 c s − ab sech (cid:18)r c s − a c s c ( x − c s t ) (cid:19) , with a = ˜ c , b = 5˜ c/ r , and c = ˜ c (4 ¯ α + r ) r / η ∗ x ∗ = − η ∗ t ∗ + O ( ε, δ ) from (4.8) and modifying accordingly thedispersive term of the BBM equation we obtain the analogous KdV equation:(4.12) η ∗ t ∗ + η ∗ x ∗ + ε η ∗ η ∗ x ∗ + δ α ∗ + 116 η ∗ x ∗ x ∗ x ∗ + ε κ ∗ η − δ γ ∗ η ∗ x ∗ x ∗ = 0 , which in dimensional form becomes(4.13) η t + ˜ cη x + 52 1 r ˜ cηη x + ˜ c (4 ¯ α + r ) r η xxx + κ η − ˜ c γη xx = 0 . The dispersion relation ω = ω ( k ) of the derived KdV equation (4.13) can be easilycomputed:(4.14) ω = ˜ c k − r ˜ c (4 ¯ α + r )16 k − i2 ( κ + γ ˜ c k ) . The imaginary part ℑ ω ( k ) comes with the negative sign, which indicates that wehave effectively introduced dissipation into the model.4.1. Symmetries and conservation laws.
Similar to the Boussinesq-type model(3.8), (3.9), the BBM equation (4.9) possesses only space and time translationssymmetries when κ + γ = 0 . However, when we neglect completely thedissipation (i.e. κ + γ = 0), another symmetry appears with the followinginfinitesimal generator: X = t ∂∂t − (cid:16) r + η (cid:17) ∂∂η . The corresponding symmetry transformation can be readily computed: x ← x ,t ← e c t ,η ← (cid:16) r + η (cid:17) e − c − r , N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 11 where c is again an arbitrary constant. The BBM equation (4.9) admits also thefollowing conservative form: (cid:16)(cid:0) η − r ( r + 4 ¯ α )48 η x x (cid:1) e κ t (cid:17) t + h r (cid:16)
52 ˜ c η + 2 ˜ c r η + r (cid:0) r κ + 4 κ ¯ α r −
48 ˜ c γ (cid:1) η x − r ( r + 4 ¯ α )12 η x t (cid:17) e κ t i x = 0 . The point symmetry group of the KdV equation (4.13) has one extra symmetrytransformation even in the dissipative case (i.e. κ + γ = 0). It is given by thefollowing infinitesimal generator: X = e − κ t ∂∂x − κ r c e − κ t ∂∂η . The corresponding symmetry transformation can be readily obtained: x ← x + c e − κ t ,t ← t ,η ← η − r κ c c e − κ t , with c an arbitrary constant. If we neglect the dissipative effects in KdV equation(4.13) (i.e. κ + γ = 0), there are two additional point symmetry transformations(we always keep time and space translations) given by infinitesimal generators: X = t ∂∂x + 2 r c ∂∂η , X = x ∂∂x + t ∂∂t − (cid:16) η + 4 r (cid:17) ∂∂η . The corresponding transformations can be readily computed: x ← x , x ← x e c ,t ← t e c , t ← t e c ,η ← (cid:16) η + 2 r (cid:17) e − c − r , η ← (cid:16) η + 2 r (cid:17) e − c − r , where c is an arbitrary constant.The KdV equation (4.13) (with dissipative terms) admits the following conser-vative form: (cid:16) η e κ t (cid:17) t + h ˜ c (cid:16) η + 54 r η − γ η x + r (4 ¯ α + r )16 η x x (cid:17) e κ t i x = 0 . In the next Section we compare the BBM equation with the Boussinesq sys-tem and we study the dissipation effects due to the viscosity of the fluid and theviscoelastic walls. 5.
The effect of viscoelasticity
In this section we study the effects of the viscoelasticity and the dissipation dueto the viscous nature of a fluid in the propagation of solitary and periodic wavesand compare the BBM equation with the classical Boussinesq system.
Dissipative effects on solitary waves.
We consider the system (3.8)–(3.9)for a vessel with undisturbed radius r = 0 . m , wall thickness h = 0 . m , Youngmodulus E = 4 . × kg/m · s , wall density ρ w = 1000 kg/m , fluid density ρ = 1060 kg/m and length 0 . m for the propagation of a solitary wave withamplitude a = 0 . m obtained using the Petviashvili method for the classicalBoussinesq system described in [20]. For the solitary waves of the BBM equationwe used the exact formula (4.11). We discretised both models using the standardpseudo-spectral method in space and the fourth-order, four-stage, classical Runge-Kutta scheme for the integration in time. Although the particular case is far frombeing realistic, the parameters are chosen to resemble a large blood vessel and itserves only the purposes of the study of the dissipative effects of the new equations.In order to study the effects of the dissipation and viscoelasticity of the vesselwall on the propagation of the solitary wave we consider three cases: (i) a vessel withelastic walls but with no viscosity and an inviscid fluid ( κ = 0 , γ = 0); (ii) a vesselwith elastic walls but with no viscosity ( γ = 0) and a viscous fluid ( κ = 1 s − ); and(iii) a vessel with viscoelastic walls ( γ = 10 − s ) and an inviscid fluid ( κ = 0). Weintegrate the system numerically until the maximum time t = 0 . s . The resultsobtained at time t = 0 . s are presented in Figure 3. -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2-0.100.10.20.30.4-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2-0.100.10.20.30.4 Figure 3.
The effects of viscoelastic walls and viscosity of the fluidin the propagation of a solitary wave. ( r = 0 . m , h = 0 . m , E = 4 . × kg/m · s , ρ w = 1000 kg/m , ρ = 1060 kg/m ). The amplitude of the solitary wave was reduced approximately by 6% when κ = 1 and γ = 0 and by 17% when κ = 0 and γ = 10 − for both models. We alsoobserve that although the initial solitary waves have slightly different shape andspeed, the shape of the solutions at t = 0 . s is very similar in both cases.In order to complete the study of the importance of the new dissipative termswe also performed an experiment combining the two dissipative terms with κ = 1 N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 13 and γ = 10 − . The amplitude of the solitary wave in this case was reduced by 20%.It seems that the effects from the viscoelastic walls are very important and shouldbe included in future studies. Moreover, we observe that small amplitude wavespropagate in the opposite direction with respect to the direction of the propagationof the solitary wave. For the accurate description of waves propagating in differentdirections it is required the model to be able to describe two-way propagation of thewaves, and the new models are capable of doing that. On the contrary unidirectionalmodels such as the KdV or BBM equations although they require the knowledgeof only one quantity (namely the initial condition for the initial disturbance of thewall or the pressure), they have certain disadvantages, especially if it is required toconsider reflections due to branching or due to other forms of obstacles.Although the propagation of a solitary wave in a dissipative environment is in-teresting for the investigation of the effects of the viscosity and viscoelasticity, it iseven interesting to examine the interaction of two solitary waves in the same envi-ronment. Even if the interaction of two solitary waves is governed by the nonlinearproperties of the model, the linear dissipative terms can affect the interaction quitedramatically. -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.4-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.050.1-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.4 Figure 4.
The effects of viscoelastic walls and viscosity of the fluidin the interaction of two solitary waves at t = 0 . sec . ( r = 0 . m , h = 0 . m , E = 4 . × kg/m · s , ρ w = 1000 kg/m , ρ =1060 kg/m ). In the next experiment we consider the interaction of two solitary waves withspeeds amplitudes 0 . m and 0 . m . All the parameters of the equations Figure 5.
The effects of viscoelastic walls and viscosity of the fluidin the interaction of two solitary waves at t = 0 . sec . ( r = 0 . m , h = 0 . m , E = 4 . × kg/m · s , ρ w = 1000 kg/m , ρ =1060 kg/m ). are the same as before whereas κ = 1 and γ = 10 − . For this we take the spatialinterval [ − . , .
5] where the solitary waves are located at x = − . x = 0 . t = 0 . sec is presented in thepresence of dissipation, and Figure 4 (c) the evolution of the same initial conditionswithout the presence of any dissipation is presented at the same time t . Thedissipation not only affect the amplitudes of the solitary waves but also the lengthof the interaction. Because the solitary waves of the BBM equation propagatefaster than the corresponding solitary waves of the classical Boussinesq system weobserve that at the same time t = 0 . sec the interaction of the two solitary wavesof the BBM equation is a more advanced stage of the separation phase comparedto the analogous results of the classical Boussinesq system. Figure 5 shows a phasediagram of the two solitary waves during the interaction for both systems.5.2. Dissipative effects on periodic waves.
In order to further assess the effectsof the dissipation due to the viscoelastic walls one can examine the linear dispersionrelations of the equations at hand. All the dispersion relations (3.10), (4.10) and
N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 15 (4.14) have a nonzero imaginary part due to the dissipative terms. Thus, planewaves of the form η ( x, t ) = A e i ( kx − ω ( k ) t ) will have a decaying amplitude of theform A e ℑ ωt . Of course the new dissipative equations are all non-linear. In fact inthis section we show that the dissipation rate estimate we obtain from the lineardispersion relation of the models predict very accurately the actual dissipation forthe nonlinear equations. For the purposes of this experiment we consider κ = 1and γ = 10 − . In order to study the combined effects of the nonlinearities withthe dissipation we consider a simple periodic wave η ( x ) = A sin( kx ) with A =10 − m and k = 2 for the same mathematical models and parameters as before.Figure 6 (a) shows the amplitude of the periodic wave as a function of time for both Figure 6.
The effects of viscoelastic walls and viscosity of the fluidin the propagation of a periodic wave at t = 1 sec . ( r = 0 . m , h =0 . m , E = 4 . × kg/m · s , ρ w = 1000 kg/m , ρ = 1060 kg/m ). cases, namely, the classical Boussinesq and the BBM equations. For the classicalBoussinesq system an approximate velocity profile was chosen to simulate one-waypropagation of the periodic wave. Because this formula is not exact we observediscrepancies between the computed and the theoretical amplitude of the classicalBoussinesq system. The predicted amplitude from the linear theory apparently isalmost identical to the numerical amplitude of the nonlinear BBM equation. Thisfact shows that dissipative effects are linear processes and thus can be includedin other derived 0D linearized models. For the damping rate of dissipative KdVequations we refer to [9]. Figure 6 (b) presents the profile of the η -solution at t = 1 sec . We observe that although the initial condition for the u -component ofthe classical Boussinesq solution is not exact, the profiles are very similar while theeffect of dissipation is very strong. It is noted that the pulses have complete oneperiod up to t = 1 sec . It is worth to mention that the wave profile will break inlater time generating dispersive shock waves as it is well-known [36]. Conclusions
In this paper we derived new weakly nonlinear and weakly dispersive asymptoticequations that describe the irrotational and dissipative flow of a fluid in pipes withviscoelastic walls. We also derived unidirectional equations of BBM and KdV typewhen the undisturbed radius is constant along the pipe. In order to study thedissipative effects due to fluid viscosity and the viscoelastic walls, we consideredsolitary and periodic waves propagating in a vessel of constant undisturbed radiusand with parameters that resemble a large blood vessel. We observed that thedissipation caused by the viscoelastic wall is equally important compared to thedissipation caused by the viscosity of the fluid or more important, and thereforeshould not be neglected. It is also observed that the dissipative effects can bedescribed very accurately by linear approximations. The new asymptotic modelshave the potential to contribute in the derivation of new lumped parameter modelsthat can be used in operational situations where measurements of the pressure andflow of the fluid are required.
Acknowledgments
The work of D. Mitsotakis was supported by the Marsden Fund administered bythe Royal Society of New Zealand. The authors would like to thank sincerely theanonymous referees who helped us to improve this paper with their comments andsuggestions.
References [1] J. Alastruey, A.W. Khir, K.S. Matthys, P. Segers, S.J. Sherwin, P.R. Verdonck, K.H. Parker,and J. Peir´o. Pulse wave propagation in a model human arterial network: Assessment of 1-Dvisco-elastic simulations against in vitro measurements.
J. Biomech. , 44(12):2250–2258, 2011.[2] J. Alastruey, K.H. Parker, J. Peir´o, and S.J. Sherwin. Lumped parameter outflow models for1-D blood flow simulations: effect on pulse waves and parameter estimation.
Comm. Comp.Phys. , 4:317–336, 2008.[3] F. Berntsson, M. Karlsson, V. Kozlov, and S. Nazarov. A one-dimensional model of viscousblood flow in an elastic vessel.
Appl. Math. Comp. , 274:125–132, 2016.[4] D. Bessems, M. Rutten, and F. van de Vosse. A wave propagation model of blood flow in largevessels using an approximate velocity profile function.
J. Fluid Mech. , 580:145–168, 2007.[5] J.L. Bona, M. Chen, and J.-C. Saut. Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I: Derivation and linear theory.
J. Non-linear Sci. , 12:283–318, 2002.[6] R. Cascaval. Variable coefficient KdV equations and waves in elastic tubes. In G.R. Gordstein,R. Nagel, and S. Romanelli, editors,
Evolution Equations , pages 57–70. Marcel Dekker, NewYork, 2003.[7] R. Cascaval. A Boussinesq model for pressure and flow velocity waves in arterial segments.
Mathematics and Computers in Simulation , 82(6):1047–1055, 2012.[8] K.B. Chandran, S.E. Rittgers, and A.P. Yoganathan.
Biofluid mechanics: the human circu-lation . CRC press, Boca Raton, 2012.[9] J.-P. Chehab and G. Sadaka. On damping rates of dissipative KdV equations.
Discrete Con-tinuous Dyn. Syst. Ser. B. , 6:1487–1506, 2013.[10] H. Demiray. Solitary wave in prestressed elastic tubes.
Bulletin of Mathematical Biology ,58:939–955, 1996.[11] H. Demiray. Waves in fluid-filled elastic tubes with a stenosis: Variable coefficients KdVequations.
J. Comp. Appl. Math. , 202:328–338, 2007.[12] D. Dutykh and F. Dias. Dissipative Boussinesq equations.
C. R. M´ecanique , 335:559–583,2007.
N SOME MODEL EQUATIONS FOR PULSATILE FLOW IN VISCOELASTIC VESSELS 17 [13] H.A. Erbay, S. Erbay, and S. Dost. Wave propagation in fluid filled nonlinear viscoelastictubes.
Acta Mechanica , 95(1):87–102, 1992.[14] C. A. Figueroa, I.E. Vignon-Clementel, K.E. Jansen, T.J.R. Hughes, and C.A. Taylor. Acoupled momentum method for modeling blood flow in three-dimensional deformable arteries.
Comput. Methods Appl. Mech. Engrg. , 195(41):5685–5706, 2006.[15] L. Formaggia, D. Lamponi, and A. Quarteroni. One-dimensional models for blood flow inarteries.
J.Eng. Math. , 47(3-4):251–276, 2003.[16] Y. Fung.
Biomechanics: Circulation . Springer-Verlag, New York, 1997.[17] T.J.R. Hughes and J. Lubliner. On the one-dimensional theory of blood flow in the largervessels.
Mathematical Biosciences , 18(1-2):161–170, 1973.[18] J.-L. Lagrange. M´emoire sur la th´eorie du mouvement des fluides.
Nouv. M´em. Acad. Berlin ,695, 1781.[19] V. Mili˘si´c and A. Quarteroni. Analysis of lumped parameter models for blood flow simulationsand their relation with 1D models.
ESAIM: Math. Model. Num. Anal. , 38:613–632, 2004.[20] D. Mitsotakis, D. Dutykh, and L. Qian. Asymptotic nonlinear and dispersive pulsatile flowin elastic vessels with cylindrical symmetry.
Computers and Mathematics with Applications ,75:4022–4047, 2018.[21] L.O. M¨uller and E.F. Toro. Well-balanced high-order solver for blood flow in networks ofvessels with variable properties.
Int. J. Numer. Meth. Biomed. Engng. , 29(12):1388–1411,2013.[22] J. P. Mynard and P. Nithiarasu. A 1D arterial blood flow model incorporating ventricu-lar pressure, aortic valve and regional coronary flow using the locally conservative Galerkin(LCG) method.
Commun. Numer. Meth. Engng , 24:367–417, 2008.[23] W. Nichols, M. O’Rourke, and C. Vlachopoulos.
McDonald’s blood flow in arteries: theoret-ical, experimental and clinical principles . CRC press, Boca Raton, 2011.[24] M.S. Olufsen, C.S. Peskin, W. Y. Kim, E.M. Pedersen, A. Nadim, and J. Larsen. Numericalsimulation and experimental validation of blood flow in arteries with structured-tree outflowconditions.
Annals of Biomedical Engineering , 28(11):1281–1299, 2000.[25] D. H. Peregrine. Long waves on a beach.
J. Fluid Mech. , 27:815–827, 1967.[26] K. Perktold and G. Rappitsch. Computer simulation of local blood flow and vessel mechanicsin a compliant carotid artery bifurcation model.
J. Biomechanics , 28(7):845–856, 1995.[27] A. Quarteroni and L. Formaggia. Mathematical modelling and numerical simulation of thecardiovascular system.
Handbook of numerical analysis , 12:3–127, 2004.[28] P. Reymond, Y. Bohraus, F. Perren, F. Lazeyras, and N. Stergiopulos. Validation of a patient-specific one-dimensional model of the systemic arterial tree.
American journal of physiology.Heart and circulatory physiology , 301(3):H1173–H1182, 2011.[29] S.J. Sherwin, V Franke, J Peir`o, and K Parker. One-dimensional modelling of a vascularnetwork in space-time variables.
J. Eng. Math. , 47(3-4):217–250, 2003.[30] N.P. Smith, A.J. Pullan, and P.J. Hunter. An anatomically based model of transient coronaryblood flow in the heart.
SIAM J. Appl. Math. , 62(3):990–1018, 2002.[31] N. Stergiopulos, D. F. Young, and T. R. Rogge. Computer simulation of arterial flow withapplications to arterial and aortic stenoses.
Journal of Biomechanics , 25(12):1477–1488, 1992.[32] R.J. Tait, T. Bryant Moodie, and J.B. Haddow. Wave propagation in a fluid-filled elastictube.
Acta Mechanica , 38(1):71–83, 1981.[33] C.A. Taylor, T.J.R. Hughes, and C.K. Zarins. Finite element modeling of blood flow inarteries.
Comp. Methods Appl. Mech. Engrg. , 158(1-2):155–196, 1998.[34] F.N. van de Vosse and N. Stergiopulos. Pulse wave propagation in the arterial tree.
Ann.Rev. Fluid Mech. , 43(1):467–499, 2011.[35] X. Wang, J.-M. Fullana, and P.-Y. Lagr´ee. Verification and comparison of four numericalschemes for a 1D viscoelastic blood flow model.
Computer Methods in Biomechanics andBiomedical Engineering , 18(15):1704–1725, 2015.[36] G. B. Whitham.
Linear and nonlinear waves , volume 42. John Wiley & Sons, New York,2011.[37] N. Xiao, J. Alastruey, and C. A. Figueroa. A systematic comparison between 1-D and3-D hemodynamics in compliant arterial models.
Int. J. Numer. Meth. Biomed. Engng. ,30(2):204–231, 2014.[38] S. Yomosa. Solitary waves in large blood vessels.
Journal of the Physical society of Japan ,56(2):506–520, 1987. [39] M. Zagzoule and J.P. Marc-Vergnes. A global mathematical model of the cerebral circulationin man.
J. Biomechanics , 19(12):1015–1022, 1986.[40] M Zamir.
The Physics of Pulsatile Flow , volume 55. Springer-Verlag, New York, 2000.
D. Mitsotakis:
Victoria University of Wellington, School of Mathematics and Sta-tistics, PO Box 600, Wellington 6140, New Zealand
E-mail address : [email protected] D. Dutykh:
Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAMA, 73000Chamb´ery, France
E-mail address : [email protected]
Q. Li:
Victoria University of Wellington, School of Mathematics and Statistics,PO Box 600, Wellington 6140, New Zealand
E-mail address : [email protected] E. Peach:
Victoria University of Wellington, School of Mathematics and Statistics,PO Box 600, Wellington 6140, New Zealand
E-mail address ::