Characterization and dynamical stability of fully nonlinear strain solitary waves in a fluid-filled hyperelastic membrane tube

Abstract

We first characterize strain solitary waves propagating in a fluid-filled membrane tube when the fluid is stationary prior to wave propagation and the tube is also subjected to a finite stretch. We consider the parameter regime where all traveling waves admitted by the linearized governing equations have nonzero speed. Solitary waves are viewed as waves of finite amplitude that bifurcate from the quiescent state of the system with the wave speed playing the role of the bifurcation parameter. Evolution of the bifurcation diagram with respect to the pre-stretch is clarified. We then study the stability of solitary waves for a representative case that is likely of most interest in applications, the case in which solitary waves exist with speed c lying in the interval $[0, c_1)$ where $c_1$ is the bifurcation value of c, and the wave amplitude is a decreasing function of speed. It is shown that there exists an intermediate value $c_0$ in the above interval such that solitary waves are spectrally stable if their speed is greater than $c_0$ and unstable otherwise.

1 Introduction

Understanding pulse wave propagation and reflection in distensible fluid-filled tubes has a range of applications such as the detection of the site of obstruction and diagnosis of the health status of arteries [1]. A complete study would need to take into account the fact that blood is non-Newtonian, and the arterial wall is dissipative, dispersive and nonlinear. However, very often only some of these features are considered depending on the particular emphasis of the study. For instance, most of the early literature is based on a linear, long wavelength (non-dispersive) theory, first put forward in [2] and refined in [3]. A linear theory taking into account the leading-order dispersive effect was derived in [4] and was used in [5] to explain the dispersive effect observed in [6]. Weakly nonlinear solitary waves based on various approximate models have been studied in [7,8,9,10,11,12,13].

The development of nonlinear elasticity theory makes it possible to describe the constitutive behavior of distensible tubes such as arteries in a more precise manner [14]. Under this theory, the dispersion relation can be derived for a general material model and a finite deformation prior to wave propagation may be taken into account. Figure 1 shows a typical dispersion curve for the Ogden material model when the tube is stretched by 20% prior to wave propagation and the inertia of fluid is neglected. The character of the two branches can be understood by noting that in the limit $k \rightarrow \infty$ the $\rho c^2$ associated with the upper and lower branches tend to, respectively, the extensional modulus and tension in the axial direction. In other words, waves associated with the upper branch behave like longitudinal waves in a rod, whereas waves associated with the lower branch are like transverse waves in a stretched string. The classical long wavelength theory corresponds to the lower branch near $k=0$. For instance, the Moens–Korteweg formula for the pulse wave speed would correspond to the speed at $k=0$, which we denote by $c_1$. It will be shown later that the lower branch is also the branch of most interest as far as solitary waves are concerned.

It is expected that the nonlinear evolution of long wavelength traveling waves is governed by the Korteweg–de Vries (KdV) equation and so there exist two families of solitary waves. Our current study is concerned with the characterization and stability of these traveling wave solutions. It is also observed that the lower branch is very sensitive to pre-stress such as the internal pressure and axial stretch, whereas the upper branch is not and neither does it vary significantly over the entire wave number regime. The speed $c_1$ decreases with increasing internal pressure. When $c_1$ vanishes, the associated solitary wave becomes a standing wave, corresponding to a static bulged configuration. Characterization of these bulging (aneurysm) solutions is the object of a series of recent studies. For instance, spectral stability of the aneurysm solutions in the absence of fluid inside the tube (pressure-controlled case) is given in [15]. A bifurcation parameter was the inflation pressure, and the authors found that the entire family of the aneurysm solutions is spectrally unstable (i. e., a localized perturbation would grow exponentially with time). The authors of [16] studied the stability of the whole branch of aneurysm solutions in the presence of a fluid with zero mean flow. It was found that the aneurysm solutions are still unstable, but the presence of fluid has a strong stabilizing effect. In [17] a stability analysis of the aneurysm solutions in the presence of a mean flow was undertaken and it was found that if the speed of the mean flow is large enough, then the aneurysm solutions may be spectrally stable. It was found in [18] that for membrane tubes with localized wall thinning there exist two families of bulging solitary waves, and the lower family with amplitudes increasing with the inflation pressure is spectrally stable. In the latter case, we can speak about the standing wave (not orbital) stability, because the problem has no translational invariance any more.

Fig. 1

Dispersion curve when the tube is stretched by 20% prior to wave propagation. Left: the two branches with positive c; right: blow-up of the upper branch showing that it is actually monotonically decreasing

In the present paper, we examine the problem of stability of solitary waves, propagating with a nonzero velocity in a fluid-filled axisymmetric membrane tube. In this case, the bifurcation parameter is the solitary wave speed.

The paper is organized as follows. In Sect. 2 we give the formulation of the problem. Section 3 is devoted to the description and computation of solitary wave solutions and their variation with respect to the pre-stretch. In Sect. 4 we perform the spectral stability analysis for solitary waves corresponding to a representative case. The paper is concluded in Sect. 5 with a summary and additional discussion.

2 Formulation of the problem

We model the tube as an incompressible, isotropic, hyperelastic, cylindrical membrane with constant undeformed radius R, thickness H and density $\rho$. The tube is assumed to be infinitely long, and end conditions are imposed at infinity. We use cylindrical polar coordinates; for the undeformed configuration, the longitudinal direction corresponds to the coordinate Z.

We assume that the axisymmetry remains throughout the entire deformation; the deformed configuration is expressed using cylindrical polar coordinates r, z, where $r=r(Z,t)$, $z=z(Z,t)$, and t denotes time.

Since the deformation is axially symmetric, the principal directions of stretch correspond to the lines of latitude (the 1-direction), the meridian (the 2-direction) and the normal to the deformed surface (the 3-direction), and the principal stretches are given by [14]

$$\begin{aligned} \lambda _1=\frac{r}{R},\quad \lambda _2=(r'^2+z'^2)^{\frac{1}{2}},\quad \lambda _3=\frac{h}{H}, \end{aligned}$$

(1)

where a prime represents differentiation with respect to Z and h denotes the deformed thickness.

The principal Cauchy stresses $\sigma _1$, $\sigma _2$, $\sigma _3$ in the deformed configuration for an incompressible material are given by

$$\begin{aligned} \sigma _i=\lambda _i \hat{W}_i-\hat{p},\qquad i=1,2,3\quad \mathrm {(no\,\, summation)}, \end{aligned}$$

where $\hat{W}=\hat{W}(\lambda _1,\lambda _2,\lambda _3)$ is the strain–energy function, $\hat{W}_i=\partial \hat{W}/\partial \lambda _i$ and $\hat{p}$ is a Lagrange multiplier associated with the constraint of incompressibility. Utilizing the incompressibility constraint $\lambda _1\lambda _2\lambda _3=1$ and the membrane assumption of no stress through the thickness direction $\sigma _3=0$, we find

$$\begin{aligned} \sigma _i=\lambda _i {W}_i, \qquad i=1,2, \end{aligned}$$

where ${W}(\lambda _1,\lambda _2)=\hat{W}(\lambda _1,\lambda _2,\lambda _1^{-1}\lambda _2^{-1})$ and ${W}_1=\partial {W}/\partial \lambda _1$ etc.

For our numerical calculations, we shall use the Ogden material model [19] for which $\hat{W}$ is given by

$$\begin{aligned} \hat{W}=\mu \sum \limits _{r=1}^3\mu _r(\lambda _1^{\alpha _r}+\lambda _2^{\alpha _r} +\lambda _3^{\alpha _r}-3)/\alpha _r, \end{aligned}$$

(2)

with $\alpha _1=1.3$, $\alpha _2=5.0$, $\alpha _3=-2.0$, $\mu _1=1.491$, $\mu _2=0.003$, $\mu _3=-0.023$, and $\mu =4.225\, \mathrm{kg}\,\mathrm{cm}^{-2}$. The parameter $\mu$ is the ground-state shear modulus.

We assume that the fluid filling the tube is inviscid and incompressible with constant density $\rho _\mathrm{f}$. It is initially at rest, and during wave propagation its axial velocity is denoted by $v_\mathrm{f}$ and the pressure in the fluid is P. The equations of motion for the tube can be derived using the linear momentum balance applied to an infinitesimal volume element of the tube wall [14] . We shall non-dimensionalize the governing equations using the following scales: R for Z, z and r, $\mu$ for the Cauchy stresses $\sigma _{1}, \,\sigma _{2}$, and the energy function W, $\mu H/R$ for P, $\sqrt{\mu /\rho }$ for $v_\mathrm{f}$, and $R\sqrt{\rho /\mu }$ for the time. Using the same notation for the scaled dimensionless variables, we have [14, 17]

$$\begin{aligned} \left[ \sigma _2 \frac{z'}{\lambda _2^2}\right] '-P rr'= \ddot{z}, \quad \left[ \sigma _2\frac{r'}{\lambda _2^2}\right] '- \frac{\sigma _1}{\lambda _1}+P r z'= \ddot{r}, \end{aligned}$$

(3)

where the dot denotes differentiation with respect to time.

Expressed in terms of the Lagrangian coordinate Z and the time t, the dimensionless equations of motion for the fluid inside the tube read [14]

$$\begin{aligned} \dot{r}z'-r'\dot{z}+v_\mathrm{f} r'+\frac{1}{2} rv_\mathrm{f}'=0,\quad b_\mathrm{f}\bigl [\dot{v}_\mathrm{f} z'-v_\mathrm{f}'\dot{z}+v_\mathrm{f} v_\mathrm{f}'\bigr ]+P'=0, \end{aligned}$$

(4)

where $b_\mathrm{f}$ is defined by

$$\begin{aligned} b_\mathrm{f}=\frac{\rho _\mathrm{f} R}{\rho H}. \end{aligned}$$

The governing equations (3) and (4) admit the uniform solution

$$\begin{aligned} r=r_\infty =\lambda _{1\infty },\quad z'= \lambda _{2\infty },\quad v_\mathrm{f}=0,\quad P=P_\infty \equiv \frac{W_{1}^\infty }{r_{\infty }\lambda _{2\infty }}, \end{aligned}$$

(5)

where $W_1=\partial W/\partial \lambda _1$ and this notion of using a subscript to signify differentiation with respect to the associated stretch will be followed throughout this paper. We also use the superscript $\infty$ to signify evaluation at $\lambda _1=r_{\infty }$, $\lambda _2=\lambda _{2\infty }$.

Before characterizing fully nonlinear solitary wave solutions, we first consider small-amplitude traveling waves superposed on the above uniform state. Assuming that the small-amplitude perturbations are proportional to

$$\begin{aligned} \mathrm{exp} \left( \mathrm{i} \frac{\lambda _{2\infty }}{r_\infty } k (Z-\hat{c} t ) \right) , \end{aligned}$$

then the scaled wavenumber k and wave speed $c=\hat{c} \lambda _{2\infty }$ satisfy the dispersion relation [17, 20]

$$\begin{aligned}&\left( k^2 m+2 \right) c^4 -\left( m \alpha _0 k^2+m \gamma _1 k^2 -m \beta _0+m \beta _1+2 \gamma _1\right) c^2 \nonumber \\&\quad + m \gamma _1(k^2 \alpha _0 +\beta _1 -\beta _0 )-m (\alpha _1-\beta _0)^2 =0, \end{aligned}$$

(6)

where

$$\begin{aligned} m= & {} 1/(b_\mathrm{f} r_{\infty }^2 \lambda _{2\infty }), \quad \alpha _0=\lambda _{2\infty } W_{2}^\infty ,\quad \alpha _1=r_{\infty }\lambda _{2\infty }W_{12}^\infty ,\\ \beta _0= & {} r_{\infty }W_{1}^\infty ,\quad \beta _1=r_{\infty }^2 W_{11}^\infty ,\quad \gamma _1=\lambda _{2\infty }^2 W_{22}^\infty . \end{aligned}$$

In the limit $k\rightarrow 0$, the two roots of the above quadratic for $c^2$ are given by

$$\begin{aligned} 4c^2 = B \pm \sqrt{B^2-8 m \varOmega }=B \pm \sqrt{[2\gamma _1-m(\beta _1-\beta _0)]^2+8m (\alpha _1-\beta _0)^2}, \end{aligned}$$

(7)

where

$$\begin{aligned} B=2\gamma _1+m(\beta _1-\beta _0), \;\;\;\; \varOmega =\gamma _1 (\beta _1-\beta _0)-(\alpha _1-\beta _0)^2. \end{aligned}$$

If the membrane tube was stress-free prior to wave propagation ($r_{\infty }=1$, $\lambda _{2\infty }=1$), $\beta _0$ would be zero since it is equal to the radial pre-stress. To derive in this case the conditions that (7) defines four real values of the velocity c, we assume that the tube material is strongly elliptic, i. e.,

$$\begin{aligned} \beta _1>0, \quad \gamma _1>0,\quad \beta _1\gamma _1-\alpha _1^2>0. \end{aligned}$$

The last inequality ensures that for the stress-free membrane from (7) we have that four values of the speed c are real (two of them are minus the other two). These characteristic speeds correspond to four branches of long waves bifurcating from the quiescent state: two of them propagate to the right, and the other two propagate symmetrically to the left. As the internal pressure or axial stretch is increased, the four wave speeds remain real until $\varOmega =0$ at which two of the speeds vanish. It was shown in [21] that this is the condition for a localized bulging solution to bifurcate from the uniformly inflated state. Our current study is concerned with the parameter regime where $\varOmega$ remains positive. In our subsequent analysis, the two positive values of c defined by (7) are denoted by $c_1$ and $c_2$ ($>c_1$), respectively.

It can be easily shown [14] that the fluid Eq. (4) in this case can be integrated to yield

$$\begin{aligned} P=P_\infty +v_\mathrm{f0}\left( 1-\frac{r^4_\infty }{r^4}\right) ,\quad v_\mathrm{f}=c-\frac{ c \,r^2_\infty }{r^2}, \end{aligned}$$

(8)

where the constant $v_\mathrm{f0}$ is defined by

$$\begin{aligned} v_\mathrm{f0}=\frac{1}{2}b_\mathrm{f}c^2. \end{aligned}$$

It is also known [14, 17] that the equations in (3) together with (8) have two integrals. They can be obtained by integrating the first equation in (3), and $z'$ times the first equation in (3) added to $r'$ times the second equation in (3), respectively. They are given by

$$\begin{aligned}&\frac{W_2 z'}{\lambda _2} -\frac{1}{2} P^* r^2-\hat{c}^2 z'=C_1,\nonumber \\&\quad W-\lambda _2 W_2 +\frac{1}{2} \hat{c}^2 \lambda _2^2=C_2, \end{aligned}$$

(9)

where a prime denotes differentiation with respect to the traveling variable $\xi =Z-\hat{c}t$,

$$\begin{aligned} P^*=P_\infty +v_\mathrm{f0} \left( 1+\frac{\lambda ^4_{1 \infty }}{r^4} \right) , \end{aligned}$$

and the constants $C_1$ and $C_2$ can be determined by evaluating the corresponding left hands at the uniform state (5). These two equations are used in the next section to determine the values of $\lambda _1(0)$ and $\lambda _2(0)$ for each pair of $r_{ \infty }$ and $\lambda _{2 \infty }$ specified.

3 Characterization and computation of solitary waves

Differentiating Eqs. (9) with respect to $\xi$, we obtain a system of first-order ordinary differential equations [21, 22]:

$$\begin{aligned} \begin{aligned} \lambda _1'&= \lambda _2 \sin \phi , \\ \lambda _2'&= \frac{{W}_1-\lambda _2 {W}_{12}}{\bigl ({W}_{22}-\hat{c}^2\bigr )} \sin \phi ,\\ \phi '&= \frac{{W}_{1}}{{W}_{2}-\hat{c}^2\lambda _2} \cos \phi -\frac{P \lambda _1 \lambda _2}{{W}_{2}-\hat{c}^2\lambda _2}, \end{aligned} \end{aligned}$$

(10)

where the prime again denotes differentiation with respect to $\xi =Z-\hat{c} t$, $\phi$ is the angle between the meridian and the z-axis (so that $\tan \phi =dr/dz=r'/z'$). Without loss of generality, we may assume that the center of the symmetric localized traveling wave is located at $\xi =0$ so that $\phi (0)=0$. Then, if $\lambda _1(0)$ and $\lambda _2(0)$ are also known, the solitary wave solution can be determined by integrating the above system as an initial value problem.

For each set of fixed values of $r_\infty$, $\lambda _{2\infty }$, we may determine the bifurcation values of the speed c with the aid of (7), and as c is reduced or increased from each bifurcation value compute numerically a family of solitary wave solutions (recall that $c=\hat{c} \lambda _{2 \infty }$). Each family begins with $c=c_1$ or $c_2$ and ends with some finite value of c.

We recall that we non-dimensionalized the pressure using the scale $\mu H/R$, so if we use the corresponding typical values of the density for rubber $\rho =1190\,\text {kg}\,\text {m}^{-3}$ and of the shear modulus specified just below (2), and $R=5$ cm and $H=0.5$ cm for the radius and wall thickness, respectively, then the non-dimensionalized pressure $P=1$ corresponds to the dimensional pressure $p=0.414\cdot 10^5$ N m$^{-2}$.

We present our calculations for the Ogden strain–energy function (2) for a tube in the absence of axial stretch ($\lambda _{2\infty }=1$) and with different values of $r_{\infty }$ as shown in Table 1. The resulting values of the corresponding dimensionless and dimensional parameters are given in this table. According to [21], a standing solitary wave with $c=0$ can bifurcate sub-critically from the uniform state when $r_{\infty }=1.687$. The bulging amplitude grows as $r_{\infty }$ decreases and approaches its maximum at $r_{\infty }=1.127$. This means that such a standing solitary wave exists for each value of $r_{\infty }$ in the interval (1.127, 1.687) and will feature in our subsequent calculations.

Table 1 Velocities and pressures in the tube for the Ogden strain–energy function, $\lambda _{2\infty }=1$ and different radial strains $r_{\infty }$

Fig. 2

a Dependence of $\lambda _1(0)$ on c for $r_{\infty }=1.5$, $\lambda _{2\infty }=1$. The horizontal line $\lambda _1(0)=1.5$ corresponds to the trivial solution and B and D are the bifurcation points. Solitary waves exist only on branches AB (solitary waves of elevation) and FP (solitary waves of depression). b Neighborhood of the bifurcation point D [not visible in the scale of (a)]. Solitary waves exist only on branch DI. The point J lies on the branch DE in (a)

Fig. 3

Dependence of $\lambda _1(0)$ on c for $r_{\infty }=1.5$, $\lambda _{2\infty }=1$. B—the bifurcation point. Solitary waves exist only on branches AB and FP ($c_B=c_F$). Only non-solitary wave solutions can be found on segments BF and CP

Fig. 4

Phase portraits of typical solitary waves of depression on branch FP, solitary waves of elevation on branch AB and non-solitary wave solutions corresponding to the segment BF; 1 - $c=0$, 2 - $c=0.75\,c_1$, 3 - $c=c_1$, 4 - $c=1.05\,c_1$, 5 - $c=1.1\,c_1$, for $r_{\infty }=1.5$, $\lambda _{2\infty }=1$

Fig. 5

Plots of $\lambda _1-\lambda _1(\infty )$ (line 1) and $\lambda _2(0)-\lambda _2(\infty )$ (line 2) against $\xi$ corresponding to the solitary wave of elevation with $c=0$. The phase portrait of this solitary wave is shown by line 1 in Fig. 4

Fig. 6

Plots of $\lambda _1-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the non-solitary wave solution with $c=1.1\, c_1$. The phase portrait of this solution is shown by line 5 in Fig. 4

Since c enters in (9) through $c^2$, we consider only the case $c \ge 0$ when we compute the wave amplitude $\varDelta \equiv \lambda _1(0)-r_{\infty }$. The dependence of $\lambda _1(0)$ on c is shown in Fig. 2a for $r_{\infty }=1.5$. This value of $r_{\infty }$ lies in the above-mentioned interval, and so a solitary wave with $c=0$ exists. This solution corresponds to point A.

Points B ($c=c_1$) and D ($c=c_2$) are bifurcation points. Figure 2b gives a blow-up of a neighborhood of the point D, which is not visible in Fig. 2a. The point on the branch EG where $\varDelta =\lambda _1(0)-r_{\infty }=0$ is not a bifurcation point, because $\lambda _2(0)$ is non-trivial there, namely $\lambda _2(0)\ne \lambda _{2\infty }$. The point J belongs to the segment DE. However, solitary waves exist only on AB, FP and DI. Figure 2b shows that solitary waves on the branch DI have very small amplitudes.

Figure 3 repeats the bifurcation curve shown in Fig. 2a on a larger scale. In Fig. 4 we plot the phase portraits of some solutions, corresponding to a selection of points in Fig. 3.

Fig. 7

Plots of $\lambda _1-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the solitary wave of depression with $c=c_1$. The phase portrait of this solitary wave is shown by line 3 in Fig. 4

Fig. 8

Plots of $\lambda _1(0)-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the solution with $c=0.055$

Profiles of different solutions corresponding to the phase portraits in Fig. 4 are given in Figs. 5, 6, 7 and 8. In Fig. 5 we present plots of $\lambda _1(0)-\lambda _1(\infty )$ (line 1) and $\lambda _2(0)-\lambda _2(\infty )$ (line 2) against $\xi$ corresponding to the solitary wave of elevation with $c=0$. The phase portrait of this bulging solitary wave is shown by line 1 in Fig. 4. In Fig. 6 we plot the form of $\lambda _1(0)-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the non-solitary wave solution with $c=1.1\, c_1$ (branch BF in Figs. 2, 3). The phase portrait of this solution is shown by line 5 in Fig. 4. Figure 7 shows plots of $\lambda _1(0)-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the solitary wave of depression with $c=c_1$. The phase portrait of this solitary wave is shown by line 3 in Fig. 4. In Fig. 8 we show the profiles of non-solitary wave solutions belonging to the segment CP in Fig. 3.

Fig. 9

Solitary wave solutions for $\lambda _{1}(\xi )$ (left) and $\lambda _{2}(\xi )$ (right) corresponding to the branch DI when $r_{\infty }=1.5$, $\lambda _{2\infty }=1$ for which $c_2=1.615$. The curves correspond to (bottom-up) $c=1.61514, 1.61519, 1.61525, 1.61530$, respectively

For the branch AB in Fig. 2a, the amplitude $\varDelta$ increases to the maximal value when $c=0$; for the branch DI in Fig. 2b, the family of solitary waves ends at the point I where the solitary wave attains its maximum amplitude. Numerical evidence shows that when the critical point I is approached the solitary wave amplitude reaches its maximum value and the value of curvature of the $\lambda _2$-wave of depression at the maximal point of the wave tends to infinity (similar to the Stokes solitary wave of maximal amplitude on a water surface; see Fig. 9).

As $r_{\infty }$ is reduced from 1.5, the bifurcation curve in Fig. 3 evolves as follows. The point F moves toward the bifurcation point B, and the point P tends to the point C. At $r_{\infty }=\lambda _{1\infty }\approx 1.3$ the points F and B coalesce. The corresponding bifurcation curve is presented in Fig. 10. In this case, solitary waves of either elevation or depression exist in the vicinity of the point B. An example of their form is given in Fig. 11.

Fig. 10

Dependence of $\lambda _1(0)$ on c for $r_{\infty }\approx 1.3$, $\lambda _{2\infty }=1$. B – the bifurcation point. Solitary waves exist on branches AB and BP

Fig. 11

Solitary waves of elevation (upper curve) and depression (lower curve) at $c=0.2466$

When $r_{\infty }$ reaches the critical value $r_{\infty }^\mathrm{cr}\approx 1.127$ the bulging solitary waves with velocity $c=0$ are replaced by a kink wave for which the profile of $\lambda _1(\xi )$ near $\xi =0$ becomes very flat: both $\lambda _1'(0)$ and $\lambda _1''(0)$ are now zero. Figure 12 shows the family of standing solitary waves for different $r_{\infty }$.

Fig. 12

Profiles of standing solitary waves corresponding to $r_{\infty } =1.126753,\, 1.13,\,1.2,\,1.3,\,1.4,1.5$. Larger amplitudes of $\lambda _{1}(\xi )$ correspond to smaller values of $r_{\infty }$

As $r_{\infty }$ is reduced even further below $r_{\infty }^\mathrm{cr}$, the standing kink wave is replaced by a kink wave having a nonzero velocity. Figure 13 shows the form of one such running kink wave for $r_{\infty }=1.125$.

The bifurcation curve for $r_{\infty }=1.1$ is presented in Fig. 14. It can be seen that when $r_{\infty }$ passes through $\lambda _{10}\approx 1.3$ there now exists a point H above the bifurcation point B on the bifurcation curve such that no solitary waves exist on the segment BH.

Fig. 13

Profile of the running kink wave corresponding to $r_{\infty }=1.125$

Fig. 14

Dependence of $\lambda _1(0)$ on c for $r_{\infty }=1.1$, $\lambda _{2\infty }=1$. B—the bifurcation point. Solitary waves exist only on branches HK and BC ($c_B=c_H$) with point K corresponding to a kink wave solution. Only non-solitary wave solutions can be found on the segment BH

Fig. 15

The phase portraits of typical solitary waves of depression on the segment CB, solitary waves of elevation on the segment HK and non-solitary wave solutions corresponding to the segment BH. 1 - $c=0.1129$, 2 - $c=0.75\,c_1$, 3 - $c=c_1$, 4 - $c=1.01\,c_1$, 5 - $c=1.035\,c_1$, for $r_{\infty }=1.1$, $\lambda _{2\infty }=1$

The phase portraits of the solutions corresponding to a selection of points on the bifurcation curve in Fig. 14 are given in Fig. 15.

Fig. 16

Plots of $\lambda _1(0)-\lambda _1(\infty )$ and $\lambda _2(0)-\lambda _2(\infty )$ against $\xi$ corresponding to the solitary wave with $c=0.1129$. The phase portrait of this soliton is shown by line 1 in Fig. 15

Fig. 17

Plots of $\lambda _1(0)-\lambda _1(\infty )$ (upper curve) and $\lambda _2(0)-\lambda _2(\infty )$ (lower curve) against $\xi$ corresponding to the solitary wave of elevation with $c= c_1$. The phase portrait of this solitary wave is shown by line 3 in Fig. 15

Profiles of some of these solutions are presented in Figs. 16 and 17. As $r_{\infty }$ is further reduced toward 1 (corresponding to the stress-free initial state), the bifurcation curve evolves as follows. At $r_{\infty }\approx 1.07$, two new points L and F appear on the bifurcation curve such that solitary waves now exist for $c>c_1$ (corresponding to points on the segment BL), which was not the case before. As $r_{\infty }$ is reduced even further, the points L and F both move away from the bifurcation point B (Fig. 18).

Fig. 18

Dependence of $\lambda _1(0)$ on c for $r_{\infty }= 1.05$, $\lambda _{2\infty }=1$. Point B is the bifurcation point. Solitary waves exist on branches CF , BL and HK. Only non-solitary (periodic) wave solutions can be found on segments FB and LH

When the point L appears, point H starts moving to the left toward point K. With further reduction in $r_{\infty }$, first at $r_{\infty }\approx 1.04$, point H merges with point K. Then, at $\lambda _{1\infty }\approx 1$ from above, point F merges with point C. After that, only solitary waves on the segment BL remain.

We note that the parameter P actually denotes the dimensionless difference between internal and external pressures, so if the external pressure is the atmospheric one ($=10^5$ Pa), the internal pressure corresponding to $\lambda _{1\infty }=1.5$, for example, is $1.315\cdot 10^5$ Pa.

4 Spectral stability

We now proceed to the stability study of the solitary waves corresponding to the branch AB in Fig. 2a; the family on DI has very low amplitude and is not of interest for applications. Denote such fully nonlinear solitary wave solutions of (10) by $r=\bar{r}(\xi ), \; z=\bar{z}(\xi ), \; P=\bar{P}(\xi )$, $v_\mathrm{f}=\bar{v}_f(\xi )$, where $\xi =Z-\hat{c} t$ and $\hat{c}$ is now the velocity of the corresponding solitary wave and $c=\hat{c}\lambda _{2\infty }$. To study the stability of these solutions, we consider axisymmetric perturbations and write

$$\begin{aligned} r(\xi , t)= & {} \bar{r}(\xi )+\varPsi (\xi ) \mathrm{e}^{\eta t}, \;\; z(\xi ,t)=\bar{z}(\xi )+\varPhi (\xi ) \mathrm{e}^{\eta t},\\ P(\xi ,t)= & {} \bar{P}(\xi )+\varPi (\xi ) \mathrm{e}^{\eta t}, \;\;\;\; v_\mathrm{f}(\xi , t)=\bar{v}_{f}(\xi )+V(\xi ) \mathrm{e}^{\eta t}, \end{aligned}$$

where the mode functions $\varPsi (\xi ), \varPhi (\xi ), \varPi (\xi ), V(\xi )$ and the growth rate $\eta$ are to be determined.

On substituting these expressions into (3) and (4) and linearizing, we obtain

$$\begin{aligned}&\left[ \frac{1}{\bar{\lambda }_2} \bar{W}_2 \varPhi '+\frac{\bar{z}'}{\bar{\lambda }_2^3}\bigl (\bar{\lambda }_2\bar{W}_{22} -\bar{W}_2\bigr )\bigl (\bar{r}'\varPsi ' +\bar{z}'\varPhi '\bigr )+ \frac{\bar{z}'}{\bar{\lambda }_2} \bar{W}_{12} \varPsi \right] '\nonumber \\&\quad -\bar{P}(\bar{r}\varPsi '+\varPsi \bar{r}')- \bar{r}\bar{r}'\varPi = \eta ^2\varPhi -2\hat{c}\eta \varPhi '+\hat{c}^2\varPhi '',\nonumber \\&\left[ \frac{1}{\bar{\lambda }_2} \bar{W}_2 \varPsi '+\frac{\bar{r}'}{\bar{\lambda }_2^3}\bigl (\bar{\lambda }_2\bar{W}_{22} -\bar{W}_2\bigr )\bigl (\bar{r}'\varPsi ' +\bar{z}'\varPhi '\bigr )+\frac{\bar{r}'}{\bar{\lambda }_2} \bar{W}_{12} \varPsi \right] '\nonumber \\&\quad -\frac{1}{\bar{\lambda }_2} \bar{W}_{12} (\bar{r}'\varPsi '+\bar{z}'\varPhi ')-\varPsi \bar{W}_{11}+ \bar{P}(\bar{r}\varPhi '+\varPsi \bar{z}')+\bar{r}\bar{z}'\varPi =\eta ^2\varPsi -2\hat{c}\eta \varPsi '+\hat{c}^2\varPsi '',\nonumber \\&\qquad \qquad \eta (\bar{z}' \varPsi -\bar{r}' \varPhi )-c_1 \varPsi ' +\bar{v}_f \varPsi '+V \bar{r}'+\frac{1}{2} \bar{r} V'+\frac{1}{2} \varPsi \bar{v}'_f=0,\nonumber \\&\qquad \qquad \qquad b_\mathrm{f} \eta (\bar{z}'V -\bar{v}'_f \varPhi )-b_\mathrm{f} c_1 V' +b_\mathrm{f} (\bar{v}_f V'+ \bar{v}_f' V)+\varPi '=0, \end{aligned}$$

where a prime denotes differentiation with respect to $\xi$.

The equations for the perturbations can be written in the form

$$\begin{aligned} \mathbf{y}'=\mathcal{M}{} \mathbf{y}, \end{aligned}$$

(11)

where $\mathbf{y}=(\varPhi , \varPhi ', \varPsi , \varPsi ', \varPi , V)^T$ and $\mathcal{M}$ is a $6 \times 6$ matrix whose components are (numerically) known functions of $\xi$ and $\eta$. This system of equations is to be solved subject to the decay conditions $\mathbf{y} \rightarrow \mathbf{0}$ as $Z \rightarrow \pm \infty$. Denoting by $\mathcal{M}_\infty$ the limit of $\mathcal{M}$ as $\xi \rightarrow \pm \infty$ and substituting a trivial solution of the form $\mathbf{y}=\mathrm{e}^{\hat{k} \xi }{} \mathbf{r}$ into $\mathbf{y}'=\mathcal{M}_\infty \mathbf{y}$, we obtain the eigenvalue problem

$$\begin{aligned} (\mathcal{M}_\infty - \hat{k} I) \mathbf{r} =\mathbf{0}, \end{aligned}$$

(12)

where I is the $6 \times 6$ identity matrix. The equation det$(\mathcal{M}_\infty - \hat{k} I)=0$ would recover the dispersion relation (6) under the substitutions

$$\begin{aligned} \hat{k}=\mathrm{i}\, k \frac{\lambda _{2 \infty }}{r_\infty }, \;\;\;\; \eta =-\mathrm{i}\, \hat{c}\hat{k}=-\mathrm{i}\, k \frac{c}{r_\infty }. \end{aligned}$$

Recalling the properties of the dispersion relation (6) established in Sect. 3, we may immediately deduce that $\hat{k}$ can be imaginary only if $\eta$ is imaginary. This implies that if $\eta$ is confined to vary in the right half of the complex plane (further denoted as $\varOmega ^+$), $\hat{k}$ can never be purely imaginary, which means that for $\eta \in \varOmega ^+$ there are exactly three $\hat{k}$ in the right half of the complex plane and the same is true for the left half of the complex plane. Denote by $\hat{k}_1, \, \hat{k}_2, \, \hat{k}_3$ the three eigenvalues of (12) with negative real parts, and by $\mathbf{r}_1, \,\mathbf{r}_2, \,\mathbf{r}_3$ the associated right eigenvectors. There then exist solutions $\mathbf{y}_i(\xi )$ of (11) such that

$$\begin{aligned} \lim \limits _{\xi \rightarrow \infty } \mathrm{e}^{-\hat{k}_i \xi } \mathbf{y}_i(\xi )=\mathbf{r}_i, \;\;\;\; i=1, 2, 3. \end{aligned}$$

Alternatively, we may consider the exterior system [23]

$$\begin{aligned} \mathbf{y}^{\wedge '}=\mathcal{M}^\wedge \mathbf{y}^\wedge \end{aligned}$$

(13)

and its adjoint system

$$\begin{aligned} \mathbf{x}^{\wedge '}=-\mathbf{x}^\wedge (\mathcal{M}^\wedge ), \end{aligned}$$

(14)

where the components of the vector function $\mathbf{y}^\wedge$ consist of all the $3\times 3$ minors of the $6\times 3$ matrix $(\mathbf{y}_1, \mathbf{y}_2, \mathbf{y}_3)$, and the components of the $20\times 20$ matrix $\mathcal{M}^\wedge$ in terms of those of $\mathcal{M}$ have previously been derived in [16].

To construct the Evans function, we need to solve (13) and (14) subject to the conditions at both infinities

$$\begin{aligned} \lim \limits _{\xi \rightarrow \infty } \mathrm{e}^{-k^\wedge \xi } \mathbf{y}^\wedge (\xi )= & {} \mathbf{r}^\wedge ,\nonumber \\ \lim \limits _{\xi \rightarrow -\infty } \mathrm{e}^{k^\wedge \xi } \mathbf{x}^\wedge (\xi )= & {} \mathbf{l}^\wedge , \end{aligned}$$

(15)

where $k^\wedge =\hat{k}_1+ \hat{k}_2+\hat{k}_3$, $\mathbf{r}^\wedge$ and $\mathbf{l}^\wedge$ are right and left eigenvectors of $\mathcal{M}^\wedge _\infty$ associated with the eigenvalue $k^\wedge$, correspondingly. The Evans function is then defined by

$$\begin{aligned} D(\eta )=\mathbf{x}^\wedge (\eta ,\xi ) \cdot \mathbf{y}^\wedge (\eta ,\xi ), \end{aligned}$$

(16)

and it can be shown that the condition of existence of an unstable eigenvalue $\eta _0$ is equivalent to $D(\eta _0)=0$ [23].

With the aid of the software package Mathematica [24], we find that for each fixed value of $b_\mathrm{f}$ there exists at least one unstable eigenvalue on the real $\eta$-axis for c less than a threshold value $c_0=c_0(b_\mathrm{f},r_{\infty },\lambda _{2\infty })$. At the value $c_0$, this eigenvalue ceases to exist and there are no unstable eigenvalues on the real $\eta$-axis for $c>c_0$. Figure 19 illustrates these facts for the cases when $r_{\infty } = 1.5$, $\lambda _{2\infty } = 1$ and $\lambda _{2\infty } = 1.25$, and the material is modeled by the Ogden model.

Fig. 19

The dependence of the real eigenvalue $\eta$ on the velocity of the solitary wave c at which the function Evans function becomes zero. The upper and lower curves correspond to $\lambda _{2\infty }=1$ and $\lambda _{2\infty }=1.25$, respectively, and in both cases $r_{\infty }=1.5$; $R=5$ cm, $H=0.5$ cm; $\rho _\mathrm{f}=\rho =1000$ kg m$^{-3}$

The problem of finding the unstable discrete spectrum $\eta$ (eigenvalues located in the right half of the complex $\eta$-plane $\varOmega ^+$ ) out of the real axis is analogous to determining the zeroes of the Evans function $D(\eta )$ lying in $\varOmega ^+\setminus {\mathbb R}^+$. The number of zeroes of $D(\eta )$ can be computed with the help of the argument principle [25]

$$\begin{aligned} n^++\frac{1}{2}=\frac{1}{2\pi i}\int \limits _\gamma \frac{D'(z)}{D(z)}\,dz, \end{aligned}$$

(17)

where $n^+$ is the number of zeroes of $D(\cdot )$ inside a closed contour $\gamma \in \varOmega ^+$ traversed in the counterclockwise direction. Since $[D(\eta )]^*=D(\eta ^*)$ we choose the contour $\gamma$ as shown in Fig. 20a. The number of zeroes due to (17) can be expressed in terms of the change in the argument $\varphi$ (with the minus sign) of $D(\eta )$, $\eta (t)\in \gamma$ as the parameter $t\in (0,t_0)$ varies from the beginning to the end of the closed contour $\gamma$;

$$\begin{aligned} n^+=\frac{1}{2\pi }\bigl [\mathrm{arg}\,D(0)-\mathrm{arg}\,D(t_0)\bigr ]-\frac{1}{2}. \end{aligned}$$

The number of zeroes of $D(\eta )$ in $\varOmega ^+$ is determined by the number of rotations of the image of the contour $\gamma$ under the mapping $D(\cdot )$ around the origin. Therefore, we need to construct the function $D(\eta )$ on $\gamma$, i.e., numerically solve the ordinary linear equations (13), (14) under the conditions (15) with $\eta$ running in a considerably large segment AB of the positive imaginary axis, the closing circle BC and a segment CA of the real axis (Fig.20a).

Fig. 20

a The contour $\gamma$ on the complex $\eta$-plane; b the scaled image of $\gamma$ on $(\mathrm{Re}\,D(\eta )/K, \mathrm{Im}\,D(\eta )/K)$-plane, where $K=|D(\eta )|^{0.95} +2$

It was shown numerically that for $c>c_0$ the image of our contour $\gamma$ under the mapping $D(\cdot )$ does not wrap around the origin. In Fig. 20b the scaled image of the contour $\gamma$ is shown. Figure 21 gives the change of the argument of $D(\gamma )$ as the parameter t runs around the contour $\gamma$.

Fig. 21

The dependence of the rotation angle $\varphi$ of the vector connecting the current point on the contour $D(\gamma )$ and the origin from the length of the contour line between the point A and the current point on $\gamma$ (the essential parameter t)

Therefore, there are no unstable eigenvalues in $\varOmega ^+$ and solitary waves for $c>c_0$ are spectrally stable. This corresponds to the nonlinear orbital stability result of solitary waves of small amplitude [26], i. e., “beginners” of the respective family of solitary waves.

5 Conclusion

The present paper is devoted to the stability analysis of solitary waves in fluid-filled membrane tubes. The tube itself is modeled as a hyperelastic membrane, and the governing equations for the tube are written on the basis of the balance of forces acting on the surface element of the shell. The fluid is considered to be ideal and incompressible, and we treat the quasi-one-dimensional flow of it. The governing equations describe fluid flow in compliant (for example, rubber) pipes.

It is shown that provided the velocity of the solitary wave is greater than the threshold value $c_0>0$, it is spectrally stable. Previously it has been shown that standing solitary waves are unstable in form [16] and only a nonzero mean flow can stabilize it, but a standing wave would propagate as soon as a perturbation is applied, so we can speak only about the stability in form, which is typical for translationally invariant systems.

We consider the spectral stability with respect to axisymmetric perturbations of solitary waves with the help of the construction of the Evans function (16) and counting its zeros (coinciding with unstable eigenvalues of the linearization of governing equations about the solitary wave solution) in the right complex half-plane of the spectral parameter $\eta$, where it is analytic.

The bulging wave family bifurcating from the smallest positive root of (7) is considered. The amplitude of the solitary wave is greater the closer to zero its velocity c. We proceed as follows. First, we find that for velocities $c<c_0$, where $c_0$ is some small critical value, there exists at least one unstable mode with the associated eigenvalue $\eta ^+$ located at the positive real $\eta$-axis ${\mathbb R}^+_\eta$. It is possible that a smaller eigenvalue $\eta ^-\ll \eta ^+$ in ${\mathbb R}^+_\eta$ also exists, but it is so small that we cannot detect it with our numerical method. Therefore, the bulging wave for $c<c_0$ is exponentially unstable. When c increases through $c_0$, the unstable eigenvalue escapes from ${\mathbb R}^+_\eta$. This happens either as a result of merging the roots $\eta ^+$ and $\eta ^-$ for $c=c_0$, or as a result of $\eta ^+$ going to zero (Fig. 19). Second, for $c>c_0$ when there are no unstable eigenvalues on ${\mathbb R}^+_\eta$, we adopt the formulation of [25] and evaluate the Evans function on the contour (Fig. 20a) with sufficiently large radius. We cannot use the whole imaginary $\eta$-axis, because in our case the Evans function $D(\eta )$ is unbounded when $|\eta |\rightarrow \infty$. We then use the argument principle to count the number of the zeroes of the Evans function (coinciding with unstable eigenvalues) inside our contour (in the large domain of the right half of the complex $\eta$-plane $\varOmega ^+$ bounded by the imaginary $\eta$-axis) and found that there are no zeroes there. This concludes the proof that bulging solitary waves for $c>c_0$ are spectrally stable.

We have only considered the case of zero mean flow ($v_\mathrm{f\infty }=0$). It was found in [17] that a mean flow has a stabilizing effect on the standing solitary waves considered. We can expect that a mean flow will have a similarly stabilizing effect on the traveling solitary waves treated here.