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

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,c1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0, c_1)$$\end{document} where c1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_1$$\end{document} 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 c0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} in the above interval such that solitary waves are spectrally stable if their speed is greater than c0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_0$$\end{document} and unstable otherwise.


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 → ∞ the ρ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 Kortewegde 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.
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.

Formulation of the problem
We model the tube as an incompressible, isotropic, hyperelastic, cylindrical membrane with constant undeformed radius R, thickness H and density ρ. 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] where a prime represents differentiation with respect to Z and h denotes the deformed thickness. The principal Cauchy stresses σ 1 , σ 2 , σ 3 in the deformed configuration for an incompressible material are given by is the strain-energy function,Ŵ i = ∂Ŵ /∂λ i andp is a Lagrange multiplier associated with the constraint of incompressibility. Utilizing the incompressibility constraint λ 1 λ 2 λ 3 = 1 and the membrane assumption of no stress through the thickness direction σ 3 = 0, we find and W 1 = ∂ W/∂λ 1 etc. For our numerical calculations, we shall use the Ogden material model [19] for whichŴ is given bŷ with α 1 = 1.3, α 2 = 5.0, α 3 = −2.0, μ 1 = 1.491, μ 2 = 0.003, μ 3 = −0.023, and μ = 4.225 kg cm −2 . The parameter μ is the ground-state shear modulus. We assume that the fluid filling the tube is inviscid and incompressible with constant density ρ f . It is initially at rest, and during wave propagation its axial velocity is denoted by v 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 , μ for the Cauchy stresses σ 1 , σ 2 , and the energy function W , μH/R for P, √ μ/ρ for v f , and R √ ρ/μ for the time. Using the same notation for the scaled dimensionless variables, we have [14,17] 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] where b f is defined by

The governing equations (3) and (4) admit the uniform solution
where W 1 = ∂ W/∂λ 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 ∞ to signify evaluation at λ 1 = r ∞ , 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 then the scaled wavenumber k and wave speed c =ĉλ 2∞ satisfy the dispersion relation [17], [20] where In the limit k → 0, the two roots of the above quadratic for c 2 are given by where If the membrane tube was stress-free prior to wave propagation (r ∞ = 1, λ 2∞ = 1), β 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., 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 Ω = 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 Ω 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 where the constant v f0 is defined by 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 where a prime denotes differentiation with respect to the traveling variable ξ = Z −ĉt, 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 λ 1 (0) and λ 2 (0) for each pair of r ∞ and λ 2∞ specified.

Characterization and computation of solitary waves
Differentiating Eqs. (9) with respect to ξ , we obtain a system of first-order ordinary differential equations [21,22]: where the prime again denotes differentiation with respect to ξ = Z −ĉt, φ is the angle between the meridian and the z-axis (so that tan φ = dr/dz = r /z ). Without loss of generality, we may assume that the center of the symmetric localized traveling wave is located at ξ = 0 so that φ(0) = 0. Then, if λ 1 (0) and λ 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 ∞ , λ 2∞ , 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 =ĉλ 2∞ ). 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 μH/R, so if we use the corresponding typical values of the density for rubber ρ = 1190 kg 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 · 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 (λ 2∞ = 1) and with different values of r ∞ 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 ∞ = 1.687. The bulging amplitude grows as r ∞ decreases and approaches its maximum at r ∞ = 1.127. This means that such a standing solitary wave exists for each value of r ∞ in the interval (1.127, 1.687) and will feature in our subsequent calculations.
Since c enters in (9) through c 2 , we consider only the case c ≥ 0 when we compute the wave amplitude Δ ≡ λ 1 (0) − r ∞ . The dependence of λ 1 (0) on c is shown in Fig. 2a for r ∞ = 1.5. This value of r ∞ 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 Δ = λ 1 (0) − r ∞ = 0 is  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.
For the branch AB in Fig. 2a, the amplitude Δ 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 λ 2 -wave of depression at the maximal point of  Fig. 4 the wave tends to infinity (similar to the Stokes solitary wave of maximal amplitude on a water surface; see Fig. 9).
As r ∞ 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 ∞ = λ 1∞ ≈ 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.
When r ∞ reaches the critical value r cr ∞ ≈ 1.127 the bulging solitary waves with velocity c = 0 are replaced by a kink wave for which the profile of λ 1 (ξ ) near ξ = 0 becomes very flat: both λ 1 (0) and λ 1 (0) are now zero. Figure 12 shows the family of standing solitary waves for different r ∞ .
As r ∞ is reduced even further below r 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 ∞ = 1.125.
The bifurcation curve for r ∞ = 1.1 is presented in Fig. 14. It can be seen that when r ∞ passes through λ 10 ≈ 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.
The phase portraits of the solutions corresponding to a selection of points on the bifurcation curve in Fig. 14 are given in Fig. 15.
Profiles of some of these solutions are presented in Figs. 16 and 17. As r ∞ is further reduced toward 1 (corresponding to the stress-free initial state), the bifurcation curve evolves as follows. At r ∞ ≈ 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 B L), which was not the case before. As r ∞ is reduced even further, the points L and F both move away from the bifurcation point B (Fig. 18).
When the point L appears, point H starts moving to the left toward point K . With further reduction in r ∞ , first at r ∞ ≈ 1.04, point H merges with point K . Then, at λ 1∞ ≈ 1 from above, point F merges with point C. After that, only solitary waves on the segment B L 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 λ 1∞ = 1.5, for example, is 1.315 · 10 5 Pa.

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 =r (ξ ), z =z(ξ ), P =P(ξ ), v f =v f (ξ ), where ξ = Z −ĉt andĉ is now the velocity of the corresponding solitary wave and c =ĉλ 2∞ . To study the stability of these solutions, we consider axisymmetric perturbations and write where the mode functions Ψ (ξ ), Φ(ξ ), Π (ξ ), V (ξ ) and the growth rate η are to be determined.
On substituting these expressions into (3) and (4) and linearizing, we obtain where a prime denotes differentiation with respect to ξ . The equations for the perturbations can be written in the form where y = (Φ, Φ , Ψ, Ψ , Π, V ) T and M is a 6 × 6 matrix whose components are (numerically) known functions of ξ and η. This system of equations is to be solved subject to the decay conditions y → 0 as Z → ±∞. Denoting by M ∞ the limit of M as ξ → ±∞ and substituting a trivial solution of the form y = ek ξ r into y = M ∞ y, we obtain the eigenvalue problem where I is the 6 × 6 identity matrix. The equation det(M ∞ −k I ) = 0 would recover the dispersion relation (6) under the substitutionsk Recalling the properties of the dispersion relation (6) established in Sect. 3, we may immediately deduce thatk can be imaginary only if η is imaginary. This implies that if η is confined to vary in the right half of the complex plane (further denoted as Ω + ),k can never be purely imaginary, which means that for η ∈ Ω + there are exactly threek in the right half of the complex plane and the same is true for the left half of the complex plane. Denote byk 1 ,k 2 ,k 3 the three eigenvalues of (12) with negative real parts, and by r 1 , r 2 , r 3 the associated right eigenvectors. There then exist solutions y i (ξ ) of (11) such that Alternatively, we may consider the exterior system [23] y ∧ = M ∧ y ∧ (13) and its adjoint system where the components of the vector function y ∧ consist of all the 3 × 3 minors of the 6 × 3 matrix (y 1 , y 2 , y 3 ), and the components of the 20 × 20 matrix M ∧ in terms of those of 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 where k ∧ =k 1 +k 2 +k 3 , r ∧ and l ∧ are right and left eigenvectors of M ∧ ∞ associated with the eigenvalue k ∧ , correspondingly. The Evans function is then defined by and it can be shown that the condition of existence of an unstable eigenvalue η 0 is equivalent to D(η 0 ) = 0 [23].
With the aid of the software package Mathematica [24], we find that for each fixed value of b f there exists at least one unstable eigenvalue on the real η-axis for c less than a threshold value c 0 = c 0 (b f , r ∞ , λ 2∞ ). At the value c 0 , this eigenvalue ceases to exist and there are no unstable eigenvalues on the real η-axis for c > c 0 . Figure 19 illustrates these facts for the cases when r ∞ = 1.5, λ 2∞ = 1 and λ 2∞ = 1.25, and the material is modeled by the Ogden model.
The problem of finding the unstable discrete spectrum η (eigenvalues located in the right half of the complex η-plane Ω + ) out of the real axis is analogous to determining the zeroes of the Evans function D(η) where n + is the number of zeroes of D(·) inside a closed contour γ ∈ Ω + traversed in the counterclockwise direction. Since [D(η)] * = D(η * ) we choose the contour γ as shown in Fig. 20a. The number of zeroes due to (17) can be expressed in terms of the change in the argument ϕ (with the minus sign) of D(η), η(t) ∈ γ as the parameter t ∈ (0, t 0 ) varies from the beginning to the end of the closed contour γ ; The number of zeroes of D(η) in Ω + is determined by the number of rotations of the image of the contour γ under the mapping D(·) around the origin. Therefore, we need to construct the function D(η) on γ , i.e., numerically solve the ordinary linear equations (13), (14) under the conditions (15) with η 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). It was shown numerically that for c > c 0 the image of our contour γ under the mapping D(·) does not wrap around the origin. In Fig. 20b the scaled image of the contour γ is shown. Figure 21 gives the change of the argument of D(γ ) as the parameter t runs around the contour γ .
Therefore, there are no unstable eigenvalues in Ω + 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.

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 η, 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 η + located at the positive real η-axis R + η . It is possible that a smaller eigenvalue η − η + in R + η 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 R + η . This happens either as a result of merging the roots η + and η − for c = c 0 , or as a result of η + going to zero (Fig. 19). Second, for c > c 0 when there are no unstable eigenvalues on R + η , 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 η-axis, because in our case the Evans function D(η) is unbounded when |η| → ∞. 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 η-plane Ω + bounded by the imaginary η-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 f∞ = 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.