On some model equations for pulsatile flow in viscoelastic vessels

TL;DR

This paper derives and analyzes asymptotic PDE models for pulsatile flow in viscoelastic vessels, incorporating viscous effects, and explores wave propagation numerically.

Contribution

It introduces new PDE models for pulsatile flow in viscoelastic vessels, including viscous effects, and investigates wave behavior numerically.

Findings

Viscoelasticity affects wave propagation characteristics.

Viscous stresses influence pulse attenuation and speed.

Numerical simulations reveal complex wave interactions.

Abstract

Considered here is the derivation of partial differential equations arising in pulsatile flow in pipes with viscoelastic walls. The equations are asymptotic models describing the propagation of long-crested pulses in pipes with cylindrical symmetry. Additional effects due to viscous stresses in bio-fluids are also taken into account. The effects of viscoelasticity of the vessels on the propagation of solitary and periodic waves in a vessel of constant radius are being explored numerically.

Equations

$$u_t+u{u}_x+v{u}_r+\frac{1}{\rho }{p}_x+\kappa u^w=0,$$

$$v_t+u{v}_x+v{v}_r+\frac{1}{\rho }{p}_r=0,$$

$$u_x+v_r+\frac{1}{r}v=0,$$

$$v(x,r,t)=0,\text{mbox}forr=0.$$

$$v(x,r,t)={\eta }_t(x,t)+(r_0(x)+\eta (x,t){)}_{x}u(x,r,t),\text{mbox}forr=r^w(x,t),$$

$${\rho }^{w}h{\eta }_{tt}(x,t)=p^w(x,t)-\frac{E_{\sigma }h}{r_0^2(x)}(\eta (x,t)+\gamma {\eta }_t(x,t)),$$

$${\varphi }_t+\frac{1}{2}{\varphi }_x^2+\frac{1}{2}{\varphi }_r^2+\frac{1}{\rho }p+\kappa \varphi =0,\text{mbox}forr=r^w.$$

$$r{\varphi }_{xx}+(r{\varphi }_r{)}_r=0,\text{mbox}0<r<r^w,$$

$${\varphi }_r=0,\text{mbox}forr=0,$$

$${\varphi }_r={\eta }_t+(r_0(x)+\eta {)}_x{\varphi }_x,\text{mbox}forr=r^w(x,t).$$

$${\eta }^{\ast }=\frac{\eta }{a},x^{\ast }=\frac{x}{\lambda },r^{\ast }=\frac{r}{R},t^{\ast }=\frac{t}{T},{\varphi }^{\ast }=\frac{1}{\lambda \epsilon \tilde{c}}\varphi ,p^{\ast }=\frac{1}{\epsilon \rho {\tilde{c}}^2}p,$$

$$\epsilon =\frac{a}{R},\delta =\frac{R}{\lambda }.$$

$${\varphi }_{t^{\ast }}^{\ast }+\frac{\epsilon }{2}{{\varphi }^{\ast }}_{x^{\ast }}^2+\frac{\epsilon }{{\delta }^2}\frac{1}{2}{{\varphi }^{\ast }}_{r^{\ast }}^2+p^{\ast }+\epsilon {\kappa }^{\ast }{\varphi }^{\ast }=0,\text{mbox}for{r}^{\ast }={r^{\ast }}^w,$$

$${\delta }^2{r}^{\ast }{\varphi }_{x^{\ast }{x}^{\ast }}^{\ast }+(r^{\ast }{\varphi }_{r^{\ast }}^{\ast }{)}_{r^{\ast }}=0,\text{mbox}0<r^{\ast }<{r^{\ast }}^w,$$

$${\varphi }_{r^{\ast }}^{\ast }=0,\text{mbox}for{r}^{\ast }=0,$$

$${\varphi }_{r^{\ast }}^{\ast }={\delta }^2{\eta }_{t^{\ast }}^{\ast }+{\delta }^2(r_0^{\ast }(x^{\ast })+\epsilon {\eta }^{\ast }{)}_{x^{\ast }}{\varphi }_{x^{\ast }}^{\ast },\text{mbox}for{r}^{\ast }={r^{\ast }}^w,$$

$$p^{\ast }={\alpha }^{\ast }{\delta }^2{\eta }_{t^{\ast }{t}^{\ast }}^{\ast }+{\beta }^{\ast }({\eta }^{\ast }+{\delta }^2{\gamma }^{\ast }{\eta }_{t^{\ast }}^{\ast }),\text{mbox}for{r}^{\ast }={r^{\ast }}^w,$$

$$\varphi (x,r,t)=\sum _{m=0}^{\infty }{r}^m{\varphi }_m(x,t).$$

$${\delta }^2{\partial }_x^2{\varphi }_{2m}+(2m+2{)}^2{\varphi }_{2m+2}=0,{\varphi }_{2m+1}=0,$$

$${\varphi }_2=-\frac{{\delta }^2}{4}{\partial }_x^2{\varphi }_0,$$

$${\varphi }_{2m+2}=\frac{{\delta }^4}{(2m+2{)}^2(2m{)}^2}{\partial }_x^4{\varphi }_{2m-2}=O({\delta }^4),$$

$${\varphi }_{2m}=(-1{)}^m\frac{{\delta }^{2m}}{2^{2m}}{\partial }_x^{2m}{\varphi }_0,$$

$$\varphi (x,r,t)=\sum _{m=0}^{\infty }{r}^{2m}(-1{)}^m\frac{{\delta }^{2m}}{2^{2m}}{\partial }_x^{2m}{\varphi }_0.$$

$$\varphi (x,r,t)={\varphi }_0(x,t)-{\delta }^2\frac{r^2}{4}{\partial }_{xx}{\varphi }_0(x,t)+O({\delta }^4).$$

$${\eta }_t+r_x^w{{\varphi }_0}_x+\frac{r^w}{2}{{\varphi }_0}_{xx}-{\delta }^2\frac{r_0^2{r_0}_x}{4}{{\varphi }_0}_{xxx}-{\delta }^2\frac{r_0^3}{16}{{\varphi }_0}_{xxxx}=O({\delta }^4,\epsilon {\delta }^2).$$

$${\varphi }_t+\frac{\epsilon }{2}{\varphi }_x^2+\frac{\epsilon }{{\delta }^2}\frac{1}{2}{\varphi }_r^2+\epsilon \kappa \varphi +\alpha {\delta }^2{\eta }_{tt}+\beta (\eta +\gamma {\eta }_t)=0.$$

$${{\varphi }_0}_t-{\delta }^2\frac{r_0^2}{4}{{\varphi }_0}_{xxt}+\frac{\epsilon }{2}{{\varphi }_0}_x^2+\epsilon \kappa {\varphi }_0+\alpha {\delta }^2{\eta }_{tt}+\beta (\eta +{\delta }^2\gamma {\eta }_t)=O({\delta }^4,\epsilon {\delta }^2).$$

$${\eta }_t+r_x^{w}w+\frac{r^w}{2}{w}_x-{\delta }^2\frac{r_0^2{r_0}_x}{4}{w}_{xx}-{\delta }^2\frac{r_0^3}{16}{w}_{xxx}=O({\delta }^4,\epsilon {\delta }^2),$$

$$w_t+(\beta \eta {)}_x+\epsilon w{w}_x+\epsilon \kappa w-{\delta }^2\frac{r_0{r_0}_x}{2}{w}_{xt}-{\delta }^2\frac{r_0^2}{4}{w}_{xxt}+\alpha {\delta }^2{\eta }_{xtt}+{\delta }^2\gamma (\beta {\eta }_t{)}_x=O({\delta }^4,\epsilon {\delta }^2).$$

$$w_t=-[\beta \eta {]}_x+O(\epsilon ,{\delta }^2),{\eta }_t=-r_x^{w}w-\frac{r^w}{2}{w}_x+O(\epsilon ,{\delta }^2).$$

$$[1-{\delta }^2\alpha {r_0}_{xx}]w_t$$

$$+\epsilon \kappa w-{\delta }^2\gamma [\beta ({r_0}_{x}w+\frac{r_0}{2}{w}_x){]}_x=O({\delta }^4,\epsilon {\delta }^2).$$

$$w(x,t)=u^{\theta }(x,t)+{\delta }^2{\theta }^2\frac{{r_0}^2}{4}{u}_{xx}^{\theta }(x,t)+O(\epsilon {\delta }^2,{\delta }^4).$$

$${\eta }_t+r_x^w{u}^{\theta }+\frac{r^w}{2}{u}_x^{\theta }+{\delta }^2\frac{(2{\theta }^2-1)r_0{r_0}_x}{4}{u}_{xx}^{\theta }+{\delta }^2\frac{(2{\theta }^2-1)r_0^3}{16}{u}_{xxx}^{\theta }=O(\epsilon {\delta }^2,{\delta }^4).$$

Full text

1. Introduction

The description of fluid flows in pipes with viscoelastic wall material is motivated mainly by the studies of hemodynamics [16]. A cardiac cycle consists of the systolic phase where the heart ventricles contract and pump blood to the arteries and the diastolic phase where the heart ventricles are relaxed and the heart fills with blood again. During the systolic phase the large arteries are deformed and store elastic energy that is released during the diastolic phase. This property of the vessels is usually referred to as the compliance of the vessels. Modelling the viscoelastic properties of the vessels appears to have significant difficulties of mathematical and numerical nature [16, 8, 27, 34, 23, 30]. The mathematical modelling of such flows suggests the use of the equations of continuum mechanics for incompressible fluid flow known as the Navier-Stokes equations.

The Navier-Stokes equations in three dimensions are too complicated to be used in practical situations and for this reason several simplified mathematical models have been derived [37, 26, 33, 14, 3, 13]. The main simplifications that have been made are based on the assumption of the axial (or cylindrical) symmetry of the vessels. This assumption and using approximations of the averaged velocity of the 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], and bidirectional models [4, 22, 7]. Bidirectional models can approximate accurately reflections of pulses occurred in the presence of non-uniformities in the vessels wall as opposed to unidirectional models such as the KdV equation ([38, 10]) where it is assumed that the pulses propagate mainly in one direction. In this paper we focus on the derivation of bidirectional equations. Further simplifications were made by assuming that the velocity of the fluid is very large. This assumption led to lumped-parameter (or zero-dimensional) models [2, 19] and the references therein. The later models are used in practice to predict the flow and the pressure of the blood in operational situations [34].

For practical reasons, the inclusion of the dissipative effects in the flow can be done by assuming a laminar flow and small viscosity. For example assuming a parabolic profile for the horizontal velocity of the fluid it has been shown that the Navier-Stokes equations can be reduced to a modified system which is very similar to the Euler equations [30], (we also refer to the Poiseuille solution for the 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}})^{2}-r^{2})u^{w}(x,t)$ these equations can be written in cylindrical coordinates in the form:

[TABLE]

where $p=p(x,r,t)$ is the pressure of the fluid, $\rho$ is the constant density of the fluid and $\kappa$ is the viscous frequency parameter (also known as the Rayleigh damping coefficient) with dimensions $[s^{-1}]$.

A sketch of the physical domain for this problem is presented in Fig. 1, where the distance of vessel’s wall from the centre of the vessel in a cross section is denoted by $r^{w}(x,t)$ and depends on $x$ and $t$ while the radius of the vessel at rest is given by the function $r_{0}(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 $\eta(x,t)$ then the vessel wall radius can be written as $r^{w}(x,t)=r_{0}(x)+\eta(x,t)$.

The governing equations (1.1)–(1.3) combined with initial and boundary conditions form a closed system. A compatibility condition is also applied at the centre of the vessel (due to cylindrical symmetry). Specifically, we assume that

[TABLE]

In general we assume for consistency purposes that $v(x,r,t)=O(r^{1+\tau})$, with $\tau\geq 0$ as $r\rightarrow 0$. The impermeability of the vessel wall can be described by the equation:

[TABLE]

and expresses that the fluid velocity equals the wall speed $v=r^{w}_{t}$. The system is accompanied also by a second boundary condition, which is the Newton second law applied on the vessel wall:

[TABLE]

where $\rho^{w}$ is the wall density, $p^{w}$ is the transmural pressure, $h$ is the thickness of the vessel wall, $E_{\sigma}=E/(1-\sigma^{2})$ where $E$ is the Young modulus of elasticity with $\sigma$ denoting the Poisson ratio of the viscoelastic wall. In this study we assume that $E$ is a constant and in general we simplify the notation by denoting $E_{\sigma}$ with $E$. The last term in (1.6) models the viscous nature of the vessel wall and can be derived by using a simple Kelvin-Voigt model (spring-dashpot model). In this setting $\gamma=n_{w}/E_{v}$ where $n_{w}$ is the dashpot coefficient of viscosity and $E_{v}$ is the Young modulus of the viscous part of Kelvin body. In practical situations the parameter $\gamma$ is very small and usually can be taken to be of order $O(10^{-4})$. It is noted that because the flow is pressure-driven the effect of gravity is neglected. For more information about the derivation of the Euler equations and the boundary conditions we refer to [40, 8].

Although the dispersion of the flow can be ignored from the majority of mathematical models derived from the Navier-Stokes equations or from the Euler equations resulting into very simple systems of conservation laws, the need for more accurate description of the waves and their reflections suggests the inclusion of this fundamental property. A first attempt towards the derivation of bi-directional weakly-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 terms can 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 effects of 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). The derivation is based on formal expansions of the velocity potential as in [12]. The new systems generalise the previously derived Boussinesq systems of [20] as they coincide with them when the viscoelastic property of the vessel wall is ignored. The new mathematical models are of significant importance since they include all the necessary ingredients for the accurate description of regular fluid flows in pipes with viscoelastic walls.

Since the dissipation caused by the fluid viscosity and the dissipation caused by the viscoelasticity of the vessel wall are different in their nature there is a question on whether the different dissipative terms have also different effects on the propagation of pulses. The answer to this question is explored computationally by studying the effects 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 new system of Boussinesq type for the description of the velocity and the deviation of the vessel wall for fluid flow in a viscoelastic vessel. This system is further improved in Section 3 by computing the fluid velocity at different levels of radius. In Section 4 further simplifications lead to unidirectional equations that depend only on the deviation of the vessel wall, while the velocity of the fluid can be computed explicitly using a simple asymptotic formula. Section 5 demonstrates the dissipation effects on the propagation of solitary and periodic waves. We close this paper with some conclusions and perspectives.

2. Derivation of the new mathematical models

Here we proceed with the derivation of the new equations. The derivation is based on the assumption that the flow is irrotational and therefore we assume the existence of a smooth velocity potential $\phi(x,r,t)$ such that $(u,v)^{T}=\nabla\phi$, i.e. we assume that $u=\phi_{x}$ and $v=\phi_{r}$. Then, as in [36] the velocity potential can be chosen appropriately such that the equations (1.1)–(1.2) can be integrated into the generalised Cauchy-Lagrange integral:

[TABLE]

The mass conservation (continuity) equation is then reduced to the elliptic equation

[TABLE]

and boundary conditions for the velocity are written as

[TABLE]

and

[TABLE]

In order to make simplifications to the previous equations we consider the following non-dimensional (scaled) variables:

[TABLE]

where $a$ is a typical deviation of the vessel wall from its rest position, $\lambda$ a typical wavelength of a pulse, $R$ is a vessel’s typical radius, $T=\lambda/\tilde{c}$ the characteristic time scale, while $\tilde{c}=\sqrt{Eh/2\rho R}$ is the Moens-Korteweg characteristic speed, [16]. It is noted that the external pressure is considered zero and is neglected. The parameters $\varepsilon$ and $\delta$ characterise the nonlinearity and the dispersion of the system:

[TABLE]

Usually, $\varepsilon$ and $\delta$ are very small. Specifically, we assume that $\varepsilon\ll 1$, $\delta^{2}\ll 1$, while the Stokes-Ursell number is of order 1: $\varepsilon/\delta^{2}=O(1)$. The system of equations (2.1)–(2.5) along with the boundary condition (1.6) is then written in dimensionless variables in the form:

[TABLE]

where $\kappa^{\ast}=\kappa\lambda/\tilde{c}\varepsilon$, $\alpha^{\ast}=\rho^{w}h/\rho R$, $\beta^{\ast}=\beta^{\ast}(x)=2R^{2}/r_{0}^{2}(x)$ and $\gamma^{\ast}=\gamma/\delta^{2}T$. For the sake of simplicity in the notation, we drop the asterisk from the new variables in the following derivations for the non-dimensional variables except if it is stated otherwise.

Following standard asymptotic techniques, cf. [5], we consider a formal expansion of the velocity potential [18]:

[TABLE]

Demanding $\phi$ to satisfy the equation (2.8) leads to the following recurrence relation

[TABLE]

for $m=0,1,2,\cdots$, where $\partial_{x}^{j}$ denotes the $j$-th order derivative with respect to $x$. A direct application of the last relation is

[TABLE]

and

[TABLE]

for $m=1,2,\cdots$. The last relation ensures that the terms $\phi_{m}$ of the velocity potential expansion for $m\geq 4$ are negligible. More general, we observe that

[TABLE]

for $m=1,2,\cdots$, and therefore

[TABLE]

A 2nd order asymptotic approximation of the velocity potential is

[TABLE]

Using the previous observations on the expansion of the velocity potential we observe that (2.10) can be approximated by the relation

[TABLE]

Since the momentum balance laws were reduced to the Cauchy-Lagrange integral equation (2.7) for $r=r^{w}$ we can eliminate the pressure using (2.11) and obtain

[TABLE]

Substituting (2.18) into (2.20) we obtain the approximate momentum equation

[TABLE]

Denoting the horizontal velocity at the centre of the vessel $u(x,0,t)={\phi_{0}}_{x}(x,t)$ by $w(x,t)$ we rewrite the equations (2.19)–(2.21) in the following form:

[TABLE]

Although the system (2.22)–(2.23) is a valid approximation of the Euler equations for the prescribed asymptotic regime, it is not of much practical use due to the temporal derivatives of the wall deviations in the momentum equation (2.23). For this reason we proceed with further simplifications using low-order approximation for the velocity $w$ and the deviation of the vessel wall.

From the equations (2.22)–(2.23) we observe that

[TABLE]

Substituting these low-order approximations into (2.23) we obtain the simplified momentum equation

[TABLE]

The system (2.22)–(2) can be the base to other more amenable Boussinesq systems along the lines of [5]. In the next section we derive a simplified Boussinesq systems with favourable properties in analogy to the classical Boussinesq system for 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) based on 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)-\delta^{2}\frac{r^{2}}{4}\;w_{xx}(x,t)+O(\delta^{4})$. From this we observe that considering the fluid velocity at any radius $r=\theta r^{w}$ for $0\leq\theta\leq 1$ we have $u^{\theta}(x,t)\doteq u(x,\theta r^{w},t)=w(x,t)-\delta^{2}\theta^{2}\frac{{r_{0}}^{2}}{4}\;w(x,t)+O(\varepsilon\delta^{2},\delta^{4})$ and thus

[TABLE]

Substitution of (3.1) into (2.22) and (2) leads to the more general Boussinesq system:

[TABLE]

[TABLE]

If we take $\gamma=0$ in the system (3.2)–(3.3) then we recover the Boussinesq system of KdV-BBM type of [20].

The linear dispersion properties of the system (3.2)–(3.3) depend on the choice of the parameter $\theta$. Taking $\theta^{2}=1/2$ the resulting system has the simplest form since the third order derivative term in (3.2) is canceled, and moreover, its linear dispersion relation is very close to the linear dispersion relation of the Euler equations [20]. Specifically, after taking $\theta^{2}=1/2$ and rearranging terms, the system (3.2)–(3.3) is simplified to the classical Boussinesq system for viscoelastic vessels:

[TABLE]

where we drop the $\theta$ in this notation by taking $u=u^{\theta}$ for $\theta=1/2$. This system is very similar to the Peregrine system of water wave theory and we usually call it 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 dimensional variables form

[TABLE]

where

[TABLE]

while all the variables $x$, $t$, $r_{0}(x)$, $\eta(x,t)$ and $u(x,t)$ are in dimensional form.

We underline that the new dissipative terms $\kappa u$ and $-\gamma[\bar{\beta}({r_{0}}_{x}u+\frac{r_{0}}{2}\;u_{x})]_{x}$ are totally different in nature and in mathematical properties. Their coefficients appeared in 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 and can be specified by the wall material properties.

System (3.8) – (3.9) extends the classical Boussinesq system derived and analysed in [20] in the case of vessels with elastic walls. Using low order corrections to the dispersive terms one can extend the whole class of the Boussinesq systems derived in [20]. The extended systems will differ only by the additional viscoelastic term and 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 assuming that the undisturbed radius of the vessel is constant. The resulting system takes the form

[TABLE]

where $\bar{\alpha}$ and $\bar{\beta}$ are constants.

The dispersion relation $\omega\ =\ \omega\,(k)$ of the derived system can be easily computed. It is given by the following quadratic algebraic equation:

[TABLE]

Here $k$ is the wavenumber and $\omega\ =\ \omega\,(k)$ is the corresponding wave frequency. The dispersion relation (3.10) can be solved explicitly for $\omega$ using the standard formulas for the roots of a quadratic equation. It will contain two branches whose expressions we omit here for brevity.

3.1. Solitary wave solutions

In the case where all the dissipative terms are neglected (i.e. when $\kappa=\gamma=0$) the system possesses classical solitary wave solutions that satisfy the following speed-amplitude relation [20]

[TABLE]

where $\zeta=1-r_{0}^{2}/(a+r_{0})^{2}$, and $s$ denotes the speed and $a$ the amplitude of the solitary wave. The graph of the speed-amplitude relationship for values of the parameters $\bar{\alpha}$ and $\bar{\beta}$ for a typical blood vessels of different radius $r_{0}$ is presented in Figure 2. We observe that there is as the amplitude of the solitary waves grows their speed tend to approach an upper bound. The upper bound for the solitary wave speed can be computed by taking the limit $a\to\infty$ in (3.11), and which is $\sqrt{3\bar{b}r_{0}}$.

For the numerical generation of these solitary waves we refer to [20]. Other more general systems similar to the systems derived in [20] can also be derived but since they only differ on the viscoelastic term we omit their derivation here and we refer to [20] for more information.

3.2. Symmetries and conservation laws

In the presence of any dissipative term (i.e. $\kappa^{\,2}\ +\ \gamma^{\,2}\ \neq\ 0$), the system (3.8), (3.9) possesses only two point symmetries: translations in time and in space. If both dissipative terms are neglected (i.e. $\kappa^{\,2}\ +\ \gamma^{\,2}\ =\ 0$), then we gain an additional point symmetry transformation whose infinitesimal generator is given by

[TABLE]

The corresponding symmetry transformation can be readily computed:

[TABLE]

where $c$ is an arbitrary constant.

For arbitrary $\kappa$ and $\gamma$ system (3.8), (3.9) can be written in an elegant conservative form:

[TABLE]

[TABLE]

For instance, these conserved quantities can be used in numerical simulations to check method accuracy.

4. Unidirectional models

In this section we present unidirectional models, namely the BBM and KdV equations by considering waves that propagate mainly in one direction. In order to derive such models we consider the following dimensionless variables:

[TABLE]

where here $c_{0}=\frac{a}{r_{0}}\sqrt{\frac{2Eh}{\rho r_{0}}}$ is a modified Moens-Korteweg characteristic speed and $T=2\frac{a\lambda}{r_{0}c_{0}}\,$. The system (3.8)–(3.9) in dimensionless variables then is written:

[TABLE]

where here $\kappa^{\ast}=2\frac{\kappa}{\sigma}$, $\gamma^{\ast}=\frac{1}{2}\frac{\sigma}{\varepsilon\delta^{2}}\gamma$, $\sigma=\frac{c_{0}}{\lambda}$ and $\varepsilon=\frac{a}{r_{0}}$, $\delta=\frac{r_{0}}{\lambda}$. Considering a flow mainly towards one direction, we can use the low-order approximation for unidirectional wave propagation [36]

[TABLE]

where $F$ and $G$ are unknown functions of $x$, $t$. Substitution of (4.4) into (4.2)–(4.3) gives the equations:

[TABLE]

Choosing appropriate functions $F$ and $G$

[TABLE]

equations (4.5) and (4.6) coincide up to the order $O(\varepsilon,\delta^{2})$ to a single equation for $\eta^{\ast}$, namely, the dimensionless BBM equation:

[TABLE]

In dimensional variables (4.8) takes the form

[TABLE]

where here $\tilde{c}=\sqrt{\frac{Eh}{2\rho r_{0}}}$ is the standard Moens-Korteweg characteristic speed. The dispersion relation $\omega\ =\ \omega\,(k)$ can be easily computed:

[TABLE]

In the absence of any form of dissipation, it is known that the BBM equation possesses classical solitary waves propagating with speed $c_{s}$ given by the formula, [20],

[TABLE]

with $a=\tilde{c}$, $b=5\tilde{c}/2r_{0}$, and $c=\tilde{c}(4\bar{\alpha}+r_{0})r_{0}/16$.

Observing that $\eta^{\ast}_{x^{\ast}}=-\eta^{\ast}_{t^{\ast}}+O(\varepsilon,\delta^{2})$ from (4.8) and modifying accordingly the dispersive term of the BBM equation we obtain the analogous KdV equation:

[TABLE]

which in dimensional form becomes

[TABLE]

The dispersion relation $\omega\ =\ \omega\,(k)$ of the derived KdV equation (4.13) can be easily computed:

[TABLE]

The imaginary part $\Im\omega\,(k)$ comes with the negative sign, which indicates that we have 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 translations symmetries when $\kappa^{\,2}\ +\ \gamma^{\,2}\ \neq\ 0\,$. However, when we neglect completely the dissipation (i.e. $\kappa^{\,2}\ +\ \gamma^{\,2}\ =\ 0$), another symmetry appears with the following infinitesimal generator:

[TABLE]

The corresponding symmetry transformation can be readily computed:

[TABLE]

where $c$ is again an arbitrary constant. The BBM equation (4.9) admits also the following conservative form:

[TABLE]

The point symmetry group of the KdV equation (4.13) has one extra symmetry transformation even in the dissipative case (i.e. $\kappa^{\,2}\ +\ \gamma^{\,2}\ \neq\ 0$). It is given by the following infinitesimal generator:

[TABLE]

The corresponding symmetry transformation can be readily obtained:

[TABLE]

with $c$ an arbitrary constant. If we neglect the dissipative effects in KdV equation (4.13) (i.e. $\kappa^{\,2}\ +\ \gamma^{\,2}\ =\ 0$), there are two additional point symmetry transformations (we always keep time and space translations) given by infinitesimal generators:

[TABLE]

The corresponding transformations can be readily computed:

[TABLE]

where $c$ is an arbitrary constant.

The KdV equation (4.13) (with dissipative terms) admits the following conservative form:

[TABLE]

In the next Section we compare the BBM equation with the Boussinesq system and we study the dissipation effects due to the viscosity of the fluid and the viscoelastic walls.

5. The effect of viscoelasticity

In this section we study the effects of the viscoelasticity and the dissipation due to the viscous nature of a fluid in the propagation of solitary and periodic waves and compare the BBM equation with the classical Boussinesq system.

5.1. Dissipative effects on solitary waves

We consider the system (3.8)–(3.9) for a vessel with undisturbed radius $r=0.01~{}m$, wall thickness $h=0.0003~{}m$, Young modulus $E=4.1\times 10^{5}~{}kg/m\cdot s^{2}$, wall density $\rho^{w}=1000~{}kg/m^{3}$, fluid density $\rho=1060~{}kg/m^{3}$ and length $0.4~{}m$ for the propagation of a solitary wave with amplitude $a=0.0035~{}m$ obtained using the Petviashvili method for the classical Boussinesq system described in [20]. For the solitary waves of the BBM equation we used the exact formula (4.11). We discretised both models using the standard pseudo-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 from being realistic, the parameters are chosen to resemble a large blood vessel and it serves 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 vessel wall on the propagation of the solitary wave we consider three cases: (i) a vessel with elastic walls but with no viscosity and an inviscid fluid ($\kappa=0,~{}\gamma=0$); (ii) a vessel with elastic walls but with no viscosity ($\gamma=0$) and a viscous fluid ($\kappa=1~{}s^{-1}$); and (iii) a vessel with viscoelastic walls ($\gamma=10^{-4}~{}s$) and an inviscid fluid ($\kappa=0$). We integrate the system numerically until the maximum time $t=0.08~{}s$. The results obtained at time $t=0.08~{}s$ are presented in Figure 3.

The amplitude of the solitary wave was reduced approximately by $6\%$ when $\kappa=1$ and $\gamma=0$ and by $17\%$ when $\kappa=0$ and $\gamma=10^{-4}$ for both models. We also observe that although the initial solitary waves have slightly different shape and speed, the shape of the solutions at $t=0.08~{}s$ is very similar in both cases.

In order to complete the study of the importance of the new dissipative terms we also performed an experiment combining the two dissipative terms with $\kappa=1$ and $\gamma=10^{-4}$. 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 should be included in future studies. Moreover, we observe that small amplitude waves propagate in the opposite direction with respect to the direction of the propagation of the solitary wave. For the accurate description of waves propagating in different directions it is required the model to be able to describe two-way propagation of the waves, and the new models are capable of doing that. On the contrary unidirectional models such as the KdV or BBM equations although they require the knowledge of only one quantity (namely the initial condition for the initial disturbance of the wall or the pressure), they have certain disadvantages, especially if it is required to consider reflections due to branching or due to other forms of obstacles.

Although the propagation of a solitary wave in a dissipative environment is interesting for the investigation of the effects of the viscosity and viscoelasticity, it is even interesting to examine the interaction of two solitary waves in the same environment. Even if the interaction of two solitary waves is governed by the nonlinear properties of the model, the linear dissipative terms can affect the interaction quite dramatically.

In the next experiment we consider the interaction of two solitary waves with speeds amplitudes $0.0035~{}m$ and $0.0005~{}m$. All the parameters of the equations are the same as before whereas $\kappa=1$ and $\gamma=10^{-4}$. For this we take the spatial interval $[-0.5,0.5]$ where the solitary waves are located at $x=-0.1$ and $x=0.1$ respectively. Although the solitary waves of the BBM equation and the classical Boussinesq system have similar shapes they differ slightly in their speeds. For this reason certain discrepancies related to the location of the solitary waves after the interaction can be observed. Figure 4 (a) depicts the initial location of the two solitary waves. In Figure 4 (b) the solution at $t=0.7~{}sec$ is presented in the presence of dissipation, and Figure 4 (c) the evolution of the same initial conditions without the presence of any dissipation is presented at the same time $t$. The dissipation not only affect the amplitudes of the solitary waves but also the length of the interaction. Because the solitary waves of the BBM equation propagate faster than the corresponding solitary waves of the classical Boussinesq system we observe that at the same time $t=0.7~{}sec$ the interaction of the two solitary waves of the BBM equation is a more advanced stage of the separation phase compared to the analogous results of the classical Boussinesq system. Figure 5 shows a phase diagram of the two solitary waves during the interaction for both systems.

5.2. Dissipative effects on periodic waves

In order to further assess the effects of the dissipation due to the viscoelastic walls one can examine the linear dispersion relations of the equations at hand. All the dispersion relations (3.10), (4.10) and (4.14) have a nonzero imaginary part due to the dissipative terms. Thus, plane waves of the form $\eta(x,t)=A_{0}e^{i(kx-\omega(k)t)}$ will have a decaying amplitude of the form $A_{0}e^{\Im\omega t}$. Of course the new dissipative equations are all non-linear. In fact in this section we show that the dissipation rate estimate we obtain from the linear dispersion relation of the models predict very accurately the actual dissipation for the nonlinear equations. For the purposes of this experiment we consider $\kappa=1$ and $\gamma=10^{-4}$. In order to study the combined effects of the nonlinearities with the dissipation we consider a simple periodic wave $\eta_{0}(x)=A_{0}\sin(kx)$ with $A_{0}=10^{-3}~{}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 cases, namely, the classical Boussinesq and the BBM equations. For the classical Boussinesq system an approximate velocity profile was chosen to simulate one-way propagation of the periodic wave. Because this formula is not exact we observe discrepancies between the computed and the theoretical amplitude of the classical Boussinesq system. The predicted amplitude from the linear theory apparently is almost identical to the numerical amplitude of the nonlinear BBM equation. This fact shows that dissipative effects are linear processes and thus can be included in other derived 0D linearized models. For the damping rate of dissipative KdV equations we refer to [9]. Figure 6 (b) presents the profile of the $\eta$-solution at $t=1~{}sec$. We observe that although the initial condition for the $u$-component of the classical Boussinesq solution is not exact, the profiles are very similar while the effect of dissipation is very strong. It is noted that the pulses have complete one period up to $t=1~{}sec$. It is worth to mention that the wave profile will break in later time generating dispersive shock waves as it is well-known [36].

6. Conclusions

In this paper we derived new weakly nonlinear and weakly dispersive asymptotic equations that describe the irrotational and dissipative flow of a fluid in pipes with viscoelastic walls. We also derived unidirectional equations of BBM and KdV type when the undisturbed radius is constant along the pipe. In order to study the dissipative effects due to fluid viscosity and the viscoelastic walls, we considered solitary and periodic waves propagating in a vessel of constant undisturbed radius and with parameters that resemble a large blood vessel. We observed that the dissipation caused by the viscoelastic wall is equally important compared to the dissipation caused by the viscosity of the fluid or more important, and therefore should not be neglected. It is also observed that the dissipative effects can be described very accurately by linear approximations. The new asymptotic models have the potential to contribute in the derivation of new lumped parameter models that can be used in operational situations where measurements of the pressure and flow of the fluid are required.

Acknowledgments

The work of D. Mitsotakis was supported by the Marsden Fund administered by the Royal Society of New Zealand. The authors would like to thank sincerely the anonymous referees who helped us to improve this paper with their comments and suggestions.