Swirling fluid flow in flexible, expandable elastic tubes: variational approach, reductions and integrability

TL;DR

This paper develops a variational theory for fluid-structure interaction in elastic tubes, incorporating vorticity effects, and derives simplified models revealing stability thresholds and connections to classical equations.

Contribution

It introduces a novel variational framework that accounts for vorticity in elastic tube fluid dynamics, reducing complex equations to a single compact form and linking to well-known integrable equations.

Findings

Constant vorticity leads to a simplified single-equation model.

Solutions exhibit instability when vorticity exceeds a threshold.

Stable regimes relate to Boussinesq, KdV, and Monge-Ampère equations.

Abstract

Many engineering and physiological applications deal with situations when a fluid is moving in flexible tubes with elastic walls. In the real-life applications like blood flow, there is often an additional complexity of vorticity being present in the fluid. We present a theory for the dynamics of interaction of fluids and structures. The equations are derived using the variational principle, with the incompressibility constraint of the fluid giving rise to a pressure-like term. In order to connect this work with the previous literature, we consider the case of inextensible and unshearable tube with a straight centerline. In the absence of vorticity, our model reduces to previous models considered in the literature, yielding the equations of conservation of fluid momentum, wall momentum and the fluid volume. We show that even when the vorticity is present, but is kept at a constant value, the case of an inextensible, unshearable and straight tube with elastics walls carrying a fluid allows an alternative formulation, reducing to a single compact equation for the back-to-labels map instead of three conservation equations. That single equation shows interesting instability in solutions when the vorticity exceeds a certain threshold. Furthermore, the equation in stable regime can be reduced to Boussinesq-type, KdV and Monge-Amp\`ere equations equations in several appropriate limits, namely, the first two in the limit of long time and length scales and the third one in the additional limit of the small cross-sectional area. For the unstable regime, we numerical solutions demonstrate the spontaneous appearance of large oscillations in the cross-sectional area.

Equations

$$\mathsf{L}=\mathsf{L}(\Lambda ,\dot{\Lambda },\mathbf{r},\dot{\mathbf{r}}):TG\to \mathbb{R},$$

$$\delta {\int }_0^T\mathsf{L}(\Lambda ,\dot{\Lambda },\mathbf{r},\dot{\mathbf{r}})\text{mbox}dt=0,$$

$$\mathbf{\gamma }={\Lambda }^{-1}\dot{\mathbf{r}},$$

$$\mathbf{\Gamma }={\Lambda }^{-1}{\mathbf{r}}^{\prime },$$

$$\delta \mathbf{\omega }=\frac{\partial \mathbf{\Sigma }}{\partial t}+\mathbf{\omega }\times \mathbf{\Sigma },\delta \mathbf{\gamma }=\frac{\partial \mathbf{\Psi }}{\partial t}+\mathbf{\gamma }\times \mathbf{\Sigma }+\mathbf{\omega }\times \mathbf{\Psi },$$

$$\delta \mathbf{\Omega }=\frac{\partial \mathbf{\Sigma }}{\partial s}+\mathbf{\Omega }\times \mathbf{\Sigma },\delta \mathbf{\Gamma }=\frac{\partial \mathbf{\Psi }}{\partial s}+\mathbf{\Gamma }\times \mathbf{\Sigma }+\mathbf{\Omega }\times \mathbf{\Psi },$$

$$\delta {\int }_0^T\ell (\mathbf{\omega },\mathbf{\gamma },\mathbf{\Omega },\mathbf{\Gamma })\text{mbox}dt=0,$$

$$\left\{\begin{array}[]{l}\displaystyle\left(\partial_{t}+\boldsymbol{\omega}\times\right)\frac{\delta\ell}{\delta\boldsymbol{\omega}}+\boldsymbol{\gamma}\times\frac{\delta\ell}{\delta\boldsymbol{\gamma}}+\left(\partial_{s}+\boldsymbol{\Omega}\times\right)\frac{\delta\ell}{\delta\boldsymbol{\Omega}}+\boldsymbol{\Gamma}\times\frac{\delta\ell}{\delta\boldsymbol{\Gamma}}=0\\ \displaystyle\left(\partial_{t}+\boldsymbol{\omega}\times\right)\frac{\delta\ell}{\delta\boldsymbol{\gamma}}+\left(\partial_{s}+\boldsymbol{\Omega}\times\right)\frac{\delta\ell}{\delta\boldsymbol{\Gamma}}=0,\end{array}\right.$$

$${\frac{\delta \ell }{\delta \mathbf{\Omega }}|}_{s=0,L}=0,{\frac{\delta \ell }{\delta \mathbf{\Gamma }}|}_{s=0,L}=0.$$

$${\partial }_t\mathbf{\Omega }=\mathbf{\omega }\times \mathbf{\Omega }+{\partial }_s\mathbf{\omega }\text{and}{\partial }_t\mathbf{\Gamma }+\mathbf{\omega }\times \mathbf{\Gamma }={\partial }_s\mathbf{\gamma }+\mathbf{\Omega }\times \mathbf{\gamma }.$$

$$s=\phi (t,s_0),s\in [0,L].$$

$$\mathcal{Q}:=\mathcal{F}$$

$$\times \left\{\phi :{\phi }^{-1}[0,L]\to [0,L]\mid \text{φ diffeomorphism}\right\}.$$

$$K_{\mathrm{rod}}=\frac{1}{2}{\int }_0^L\left(\alpha |\mathbf{\gamma }{|}^2+a{\dot{R}}^2+\mathbb{I}(R)\mathbf{\omega }\cdot \mathbf{\omega }\right)|\mathbf{\Gamma }|\text{mbox}ds,$$

$$\frac{d}{dt}\mathbf{r}(t,\phi (t,s_0))$$

$$={\partial }_t\mathbf{r}(t,\phi (t,s_0))+{\partial }_s\mathbf{r}(t,\phi (t,s_0))u(t,\phi (t,s_0)),$$

$$u(t,s)=\left({\partial }_t\phi \circ {\phi }^{-1}\right)(t,s),s\in [0,L].$$

$$K_{\mathrm{fluid},\mathrm{transl}}=\frac{1}{2}{\int }_{{\phi }^{-1}(0,t)}^{{\phi }^{-1}(L,t)}\rho Q_0(s_0){\left|\frac{d}{dt}\mathbf{r}(t,\phi (t,s_0))\right|}^2\text{mbox}d{s}_0,$$

$$\frac{D\lambda }{Dt}=({\partial }_t+u{\partial }_s)\lambda =0,\lambda :={\int }_{C(t)}\mathbf{v}\cdot \text{mbox}d\mathbf{l}=0,$$

$$({\partial }_t+u{\partial }_s)\lambda =0\iff \lambda ={\lambda }_0\circ {\phi }^{-1}(s,t).$$

$${\partial }_t\xi +u{\partial }_s\xi -\xi {\partial }_{s}u=0.$$

$$\Phi =\frac{2\pi \rho {\int }_0^{R}r{v}_{\mathrm{rot}}^2\text{mbox}dr}{\lambda (s,t{)}^2}|\mathbf{\Gamma }|,\Phi =\Phi (\mathbf{\Omega },\mathbf{\Gamma },R,R^{\prime },\dots ).$$

$$\Phi =\frac{1}{{\lambda }^2}\int 2\pi \rho {\int }_{ϵ}^{R}r{v}_{\mathrm{rot}}^2\text{mbox}dr=\frac{1}{2\pi }\log \frac{R}{ϵ}$$

$$K_{\mathrm{fluid},\mathrm{rot}}=\frac{\pi }{2}\int \rho \Phi {\lambda }^2\text{mbox}ds.$$

$$K_{\mathrm{fluid}}=\frac{1}{2}{\int }_0^L\rho Q{\left|\mathbf{\gamma }+\mathbf{\Gamma }u\right|}^2+\rho \Phi {\lambda }^2\text{mbox}ds,$$

$$E_{\rm rod}=\frac{1}{2}\int_{0}^{L}\Big{(}\mathbb{J}\boldsymbol{\Omega}\!\cdot\!\boldsymbol{\Omega}+\lambda(R)|\boldsymbol{\Gamma}-\boldsymbol{\chi}|^{2}+2F(R,R^{\prime},R^{\prime\prime})\Big{)}|\boldsymbol{\Gamma}|\mbox{d}s\,,$$

$$\rho Q=(\rho Q_0\circ {\phi }^{-1}){\partial }_s{\phi }^{-1}.$$

$$Q=Q_0\circ {\phi }^{-1}{\partial }_s{\phi }^{-1}\Rightarrow {\partial }_{t}Q+{\partial }_s(Qu)=0.$$

$$\mathsf{L}=\mathsf{L}\big{(}\Lambda,\dot{\Lambda},\mathbf{r},\dot{\mathbf{r}},\varphi,\dot{\varphi},R,\dot{R}\big{)}:T\mathcal{Q}\rightarrow\mathbb{R},\quad\mathsf{L}=K_{\rm rod}+K_{\rm fluid}-E_{\rm rod}-E_{\rm int},$$

$$\ell (\mathbf{\omega },\mathbf{\gamma },\mathbf{\Omega },\mathbf{\Gamma },u,Q,\lambda ,R,\dot{R})$$

$$\displaystyle=\int_{0}^{L}\Big{[}\Big{(}\frac{1}{2}\alpha|\boldsymbol{\gamma}|^{2}+\frac{1}{2}\mathbb{I}(R)\boldsymbol{\omega}\!\cdot\!\boldsymbol{\omega}+\frac{1}{2}a\dot{R}^{2}-F(R,R^{\prime})-\frac{1}{2}\mathbb{J}\boldsymbol{\Omega}\!\cdot\!\boldsymbol{\Omega}$$

$$\displaystyle\qquad\qquad-\frac{1}{2}S(R)|\boldsymbol{\Gamma}-\boldsymbol{\chi}|^{2}\Big{)}|\boldsymbol{\Gamma}|+\frac{1}{2}\rho Q\left|\boldsymbol{\gamma}+\boldsymbol{\Gamma}u\right|^{2}+\frac{1}{2}\rho\Phi\lambda^{2}s-p_{\rm ext}Q\Big{]}\mbox{d}s$$

$$\displaystyle=:\int_{0}^{L}\Big{[}\ell_{0}(\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\Gamma},u,Q,R,\dot{R},R^{\prime},\lambda)-p_{\rm ext}Q\Big{]}\mbox{d}s\,.$$

$$\delta {\int }_0^T\mathsf{L}(\Lambda ,\dot{\Lambda },\mathbf{r},\dot{\mathbf{r}},\phi ,\dot{\phi },\lambda ,R,\dot{R})\text{mbox}dt=0,$$

Full text

Keywords: Variational methods, Fluid-Structure Interaction, Collapsible tubes carrying fluid, Spiral blood flow, Boussinesq equation, Monge-Ampère equation

Declarations of interest: none.

1. 1 Introduction

2. 2 Exact geometric theory for flexible tubes conveying incompressible fluid

3. 2.1 Background for geometrically exact rod theory

4. 2.2 Definition of the configuration space for the tube with expandable walls conveying fluid

5. 2.3 Definition of the Lagrangian

6. 2.4 Variational principle and equations of motion

7. 3 Tubes with a straight centerline: dynamics and reduction to a single equation

8. 3.1 Equations of motion for inextensible and unshearable tubes

9. 3.2 The case of uniform circulation, energy, and the new constant of motion

10. 3.3 Equations of motion for the back-to-labels map

11. 3.4 The case of constant vorticity

12. 3.5 Energy conservation and additional constant of motion

13. 3.6 Linear stability and instability

14. 3.7 Numerical simulations: stable case

15. 4 Asymptotic analysis of equation (3.35)

16. 4.1 Small cross section (collapse) and reduction to Monge-Ampère equation

17. 4.2 Approximation by Boussinesq-type equations

18. 4.3 Reduction to Korteweg-deVries (KdV) equation

19. 5 Conclusions

20. 6 Acknowledgements

1 Introduction

Models of flexible tubes conveying fluid have been under investigation for many decades, because of their importance to industrial and biomedical applications. Broadly speaking, the research can be put in two different categories. The first category of papers focused on the instability of centerline of the tubes due to the complex interaction of the fluid inside the tube and the tube’s elasticity, most often treated as an elastic rod. For such systems, an instability appears when the flow rate through the tube exceeds a certain critical value, sometimes referred to as the garden hose instability. While this phenomenon has been known for a very long time, the quantitative research in the field started around 1950 [1]. We believe that Benjamin [2, 3] was the first to formulate a quantitative theory for the 2D dynamics of the initially straight tubes by considering a linked chain of tubes conveying fluids and using an augmented Hamilton principle of critical action that takes into account the momentum of the jet leaving the tube. A continuum equation for the linear disturbances was then derived as the limit of the discrete system. Interestingly, Benjamin’s approach is the inverse to the derivation of variational integrators for continuum equations [4, 5, 6] This linearized equation for the initially straight tubes was further considered by Gregory and Païdoussis [7].

These initial developments formed the basis for further stability analysis of this problem for finite, initially straight tubes [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. The linear stability theory has shown a reasonable agreement with experimentally observed onset of the instability [18, 10, 19, 20, 21]. Nonlinear deflection models were also considered in [22, 14, 23, 24], and the compressible (acoustic) effects in the flowing fluid in [25]. Alternatively, a more detailed 3D theory of motion was developed in [26] and extended in [27], based on a modification of the Cosserat rod treatment for the description of elastic dynamics of the tube, while keeping the cross-section of the tube constant and orthogonal to the centerline. In particular, [27] analyzes several non-straight configurations, such as tube hanging under the influence of gravity, both from the point of view of linear stability and nonlinear behavior. Unfortunately, this Cosserat-based theory could not easily incorporate the effects of the cross-sectional changes in the dynamics. Some authors have treated the instability from the point of view of the follower force approach, which treats the system as an elastic beam, ignoring the fluid motion, with a force that is always tangent to the end of the tube. Such a force models the effect of the jet leaving the nozzle [28], showing interesting dynamics of solutions. However, once the length of the tube becomes large, the validity of the follower force approach has been questioned, see [29] for a lively and thorough discussion. For the history of the development of this problem in the Soviet/Russian literature, we refer the reader to the monograph [30] (still only available in Russian). To briefly touch upon the developments in Russian literature that have been published in parallel with their western counterparts, and perhaps less known in the west, we refer the reader to the selection of papers, with most of them now available in English [31, 32, 33, 34, 35, 36, 37, 38, 39, 40].

It is also worth noting the developments in the theory of initially curved pipes conveying fluid, which has been considered in in some detail in the literature. The equations of motion for such theory were initially derived using the balance of elastic forces from tube’s deformation and fluid forces acting on the tube when the fluid is moving along a curved line in space. In the western literature, we shall mention the earlier work [41], followed with more detailed studies [42, 43, 44] which developed the theory suited for both extensible and inextensible tubes and discussed the finite-element method realization of the problem. We shall also mention [45, 46] deriving a variational approach for the planar motions of initially circular tubes, although the effect of curved fluid motion was still introduced as extra forces through the Lagrange-d’Alembert principle. In the Soviet/Russian literature, [36] developed the rod-based theory of oscillations and [39] considered an improved treatment of forces acting on the tubes. Most of the work has been geared towards the understanding of the planar cases with in-plane vibrations as the simplest and most practically relevant situations (still, however, leading to quite complex formulas).

In spite of substantial amount of literature in previous work in the area, many questions remain difficult to answer using the traditional approach. Most importantly, it is very difficult (and perhaps impossible) to extend the previous theory to accurately take into account the changes in the cross-sectional area of the tube, also called the collapsible tube case. In many previous works, the effects of cross-sectional changes have been considered through the quasi-static approximation: if $A(s,t)$ is the local cross-section area, and $u(s,t)$ is the local velocity of the fluid, with $s$ being the coordinate along the tube and $t$ the time, then the quasi-static assumption states that $uA=$ const, see, e.g., [22, 47]. Unfortunately, this simple law is not correct in general and should only be used for steady flows. This problem has been solved in the fully variational derivation [48, 49], where a geometrically exact setting for dealing with a variable cross-section depending on tube’s deformation was developed and studied, showing the important effects of the cross-sectional changes on both linear and nonlinear dynamics. The nonlinear theory was derived from a variational principle in a rigorous geometric setting and for general Lagrangians. It can incorporate general boundary conditions and arbitrary deviations from equilibrium in the three-dimensional space. From a mathematical point of view, the Lagrangian description of these systems involves both left-invariant (elastic) and right-invariant (fluid) quantities. This method has since been applied to the linear stability of initially curved (helical) tubes in [50], where it was shown that the geometric method of [48, 49] is guaranteed to yield an equation with constant coefficients for the linear stability analysis for any helical initial state. This theory further allowed consistent variational approximations of the solutions, both from the point of view of deriving simplified reduced models and developing structure preserving numerical schemes [51], echoing the original development of the equations by Benjamin [2, 3] for constant cross-section.

The second category of papers concerns with the flow of fluid through tubes which can dynamically change its cross-section. In these works, the cross-section, or its radius, in the case of circular cross-sections, serves as an additional dynamic variable. The equations of motion then involve the balance of fluid momentum, wall momentum and conservation law for the fluid. The theory for such tubes is indispensable in biomedical applications, such as arterial flows, and the key initial progress in theoretical understanding is based on balancing the conservation laws in a collapsible tube, see [52]. Further work in this field explored applications to arterial [53, 54, 55, 56, 57, 58, 59] and lung flows [60, 61, 62]. Analytical studies for such flows are usually limited to cases when the centerline of the tube is straight. In addition, for analytical progress, one needs to further assume that material particles on the wall can only move normally to the centerline during the dynamics (see, however, [63] for non-trivial longitudinal displacement). We should also mention the instability through the neck formation and self-sustaining flow pulsations suggested by [64] and experimentally measured in [65], the treatment of forward and backward running waves in arteries [66, 67] and the treatment of poroelastic walls [68]. While substantial progress in the analysis of the flow has been achieved so far, it was difficult to describe analytically the general dynamics of 3D deformations of the tube involving e.g. the combination of shear, transversal deformation, and extensions. A lively discussion of the dichotomy between the needs of numerical solutions of Navier-Stokes equations in realistic geometries, and the importance of theoretical developments based on simplified models can be found in [69]. For more informations about the application of collapsible tubes to biological flows, we refer the reader to the reviews [70, 71, 72].

A fully variational theory of motion of tubes with expandable walls in three dimensions was derived in [73], both for compressible and incompressible flow motion inside the tube. The major difference with the earlier models described by the variational approach [48, 49] was the ability to treat compressible flows, and incorporate the coupled dynamics of the walls, the tube and the fluid. This model can treat arbitrary motion of the tubes in three dimensions, using arbitrary elasticities and mass distribution of the walls. For compressible flows, the model is also capable to incorporate the motion of shock waves in a moving tubes with expandable walls. It was shown that all simplified analytical models of flows with expandable walls can be obtained as particular cases of the more generalized model, with the introduction of appropriate friction terms. As far as we are aware, there is no way to achieve such generality of motion with non-variational models. Indeed, one can guess the right terms in simplified geometries using force balance, but trying to balance forces in more complex problems, such as deforming centerline and expandable walls, is fraught with possibilities for errors. In contrast, variational methods provide a systematic way of deriving the equations and their subsequent analysis, such as Hamiltonian formulation and the existence of Poisson brackets. We thus believe that the the problems in fluid-structure interactions can benefit greatly from the application of modern variational methods and geometric mechanics.

Most of the analytical aspects of previous works discussed above have treated the flow of fluid in a simplified manner, treating one characteristic streamwise velocity for each point of the centerline. However, the fluid flow regimes with vorticity are frequently encountered around stents in arteries, have been observed both experimentally and numerically. They are also known as ’spiral flows’, or ’spiral laminar flows’, and are believed to play a crucial role in arterial dynamics. The study of these flows with applications to blood flow dynamics started with the pioneering observation by Stonebridge & Drophy [74] and continued by many authors [75, 76, 77, 78, 79, 80]. Numerical studies of vorticity transport in complex geometries relevant to arteries can be found in e.g. [81], but as far as we are aware, there has been no theoretical studies including the wall dynamics, fluid dynamics and vorticity, especially from the variational point of view. The goal of this paper is precisely to fill this theoretical gap. More precisely, we shall derive a fully three-dimensional flow of mean vorticity transport based on Kelvin’s theorem of circulation in fluid mechanics. We will mostly focus on the variational derivation, and thus most of our considerations will be dedicated to the inviscid models. The friction terms essential for blood flow applications can be added later as they have to be incorporated in the variational principle using the Lagrange-d’Alembert’s method of external forces. Our model will be valid for arbitrary three-dimensional deformations of the centerline and arbitrary elasticity laws satisfied by the wall. As we show, the simplified mathematical treatment for the flow can produce interesting mechanical effects and instability caused by the interaction of wall and fluid. As far as practical applications of the theory presented here, we view it as a relatively complete model of industrial flows in flexible and expandable tubes with high vorticity, and a first step towards understanding the complete structure of inertial terms in arterial models, to which terms containing the internal viscous friction and wall shear can later be added.

A major role in this manuscript will be played by the back-to-label Lagrangian map for the fluid part, which maps the location of fluid particle at a given point on the tube at the time $t$ to its original location at $t=0$. The back-to-labels map has an important role in the development of continuum mechanics. Of particular value to our studies here is the Clebsch representation of back-to-labels map which was used to derive variational and Hamiltonian structure of continuum mechanics. One of the pioneering examples of this approach was developed in [82], where Poisson brackets and Clebsch representation for fluid mechanics, magnetohydrodynamics (MHD), and elasticity theory were derived. Nonlinear stability theory of fluid flows using back-to-labels maps was investigated in [83, 84], and applications to quasi-geostrophic flows was studied in [85]. A multi-symplectic formulation of fluid mechanics for ideal fluids using back-to-labels maps was further derived in [86]. In these papers, variational derivatives with respect to back-to-labels map was crucial in obtaining the laws of motion and proving the Hamiltonian structure of the flow. Because of the features of 3D hydrodynamics, the expression of the incompressibility condition in terms of back-to-labels map is rather awkward, and thus the works cited here used variations with constrained Lagrangian and needed the pressure as the corresponding Lagrange multiplier. Our problem has the advantage that one can explicitly connect the cross-sectional area and the spatial derivatives of back-to-labels map, and substitute into the Lagrangian directly, leading to the unconstrained Euler-Lagrange equations.

Structure of the paper and main results

The paper is structured as follows. In Section 2, we derive the full three-dimensional equations of motion for flexible tubes carrying incompressible fluid which has streamwise velocity and and vorticity. We have tried to keep the derivation self-contained and pedagogical, and build up our theory from the theory of exact geometric rods with no fluid, to the theory of fluid-conveying tubes with expandable walls, introducing the swirling motion of the fluid in the process. The fluid volume incompressibility condition leads to the appearance of a pressure-like term in the equations of motion. In Section 3, we consider the case of the inextensible and unshearable tube, and a static straight centerline when the walls can only move normal to the centerline, which is the case often considered in the literature for theoretical studies. The equations of motion then reduce to a generalization of the equations familiar from the previous literature, with four equations of motion describing the conservation of the fluid and wall momenta and the fluid’s mass, and the transport of vorticity (or, more precisely, velocity circulation), with the corrections due to the vorticity terms. These equations explicitly contain the fluid pressure as the Lagrange multiplier for incompressibility. We show the conservation of energy for these equations and also derive a new constant of motion linear in momenta in (3.9), that seems to have not been observed in the previous works on the subject. By using an alternative derivation using the Lagrangian back-to-labels map, we show how to reduce this equation to a single nonlinear scalar equation of motion in Section 3.3. In the linearized form, the equation shows interesting stability properties, in particular, instability for a sufficiently large velocity in the tube. We perform the numerical studies of the solution in both linearly stable and unstable regimes. For the linearly stable regime, in Section 4 we utilize the slow-time and large-wavelength approximation to first develop a reduction to a Monge-Ampére’s equation, then a Boussinesq-like equation and finally the KdV equation. We compare the numerics of the approximate models with the full equations and discuss the relevance of KdV approximations for the left- and right-running waves.

2 Exact geometric theory for flexible tubes conveying incompressible fluid

In this section, we first quickly review the Lagrangian variational formulation for geometrically exact rods without fluid motion as the foundation of the theory. Then, we extend this variational formulation to incorporate the motion of the incompressible fluid inside the tube and the motion of the wall, and finish by introducing the vorticity. To achieve this goal we need to first identify the configuration manifold of the system (which is infinite-dimensional), as well as the convective and spatial variables for the tube and the fluid. Our model is then obtained by an application of the Hamilton principle, reformulated in convective variables for the tube and in spatial variables for the fluid. For rigorous justification of the variational approach to fluid-structure interactions see [49, 73].

2.1 Background for geometrically exact rod theory

Here we briefly review the theory of geometrically exact rods following the approach developed, on the Hamiltonian side, in [87] and subsequently in the Lagrangian framework more appropriate to this article in [88, 89]. A more comprehensive introduction is contained in [49, 50] to which we refer the reader for the details. The purely elastic (i.e., rods carrying no fluid) geometrically exact theory is equivalent to the Cosserat’s rods [90], the equivalence of these approaches was shown in [87, 89].

The configuration of the rod deforming in the ambient space $\mathbb{R}^{3}$ is defined by specifying the position of its line of centroids by means of a map $\mathbf{r}(t,s)\in\mathbb{R}^{3}$, and by giving the orientation of the cross-section at that point. Here $t$ is the time and $s\in[0,L]$ is a parameter along the strand that does not need to be arclength. The orientation of the cross-section is given by a moving basis $\{\mathbf{e}_{i}(t,s)\mid i=1,2,3\}$ attached to the cross section relative to a fixed frame $\{\mathbf{E}_{i}\mid i=1,2,3\}$. The moving basis is described by an orthogonal transformation $\Lambda(t,s)\in SO(3)$ such that $\mathbf{e}_{i}(t,s)=\Lambda(t,s)\mathbf{E}_{i}$. We interpret the maps $\Lambda(t,s)$ and $\mathbf{r}(t,s)$ as a curve $t\mapsto(\Lambda(t),\mathbf{r}(t))\in G$ in the infinite dimensional Lie group $G=\mathcal{F}([0,L],SO(3)\times\mathbb{R}^{3})$ of $SO(3)\times\mathbb{R}^{3}$-valued smooth maps defined on $[0,L]$ 111The notation $\mathcal{F}([0,L],V)$ refers to the set of functions defined on the interval $s\in[0,L]$ taking the values in the set $V$.. The Lie group $G$ is the configuration manifold for the geometrically exact rod.

Following Hamilton’s principle, given a Lagrangian function

[TABLE]

defined on the tangent bundle $TG$ of the configuration Lie group $G$ defined above, the equations of motion are the Euler-Lagrange equations obtained by the critical action principle

[TABLE]

for arbitrary variations $\delta\Lambda$ and $\delta\mathbf{r}$ vanishing at $t=0,T$. It turns out that the Lagrangian of geometrically exact rods can be exclusively expressed in terms of the convective variables

[TABLE]

see [87], where $\boldsymbol{\gamma}(t),\boldsymbol{\omega}(t)\in\mathcal{F}([0,L],\mathbb{R}^{3})$ are the linear and angular convective velocities and $\boldsymbol{\Gamma}(t),\boldsymbol{\Omega}(t)\in\mathcal{F}([0,L],\mathbb{R}^{3})$ are the linear and angular convective strains. This gives rise to a Lagrangian $\ell=\ell(\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\Gamma}):\mathcal{F}([0,L],\mathbb{R}^{3})^{4}\rightarrow\mathbb{R}$ written exclusively in terms of convective variables. Here, we treat all four variables as infinite-dimensional functions of $t$, which are themselves $\mathbb{R}^{3}$-valued functions of $s$ for any fixed $t$, see [89] for the detailed mathematical exposition of the method. For the moment, we leave the Lagrangian function unspecified, we will give its explicit expression later in §2.3 for the case of fluid-conveying tubes.

The equations of motion in convective description are obtained by writing the critical action principle (2.1) in terms of the Lagrangian $\ell$. This is accomplished by computing the constrained variations of $\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\Gamma}$ induced by the free variations $\delta\Lambda,\delta\mathbf{r}$ via the definitions (2.2). We find

[TABLE]

where $\boldsymbol{\Sigma}(t,s)=\Lambda(t,s)^{-1}\delta\Lambda(t,s)\in\mathbb{R}^{3}$ and $\boldsymbol{\Psi}(t,s)=\Lambda(t,s)^{-1}\delta\mathbf{r}(t,s)\in\mathbb{R}^{3}$ are arbitrary functions vanishing at $t=0,T$. Hamilton’s principle (2.1) induces the variational principle

[TABLE]

with respect to the constrained variations (2.3), (2.4), which in turn yields the reduced Euler-Lagrange equations

[TABLE]

together with the boundary conditions

[TABLE]

If one of the extremity (say $s=0$) of the rod is kept fixed, i.e., $\mathbf{r}(t,0)=\mathbf{r}_{0}$, $\Lambda(t,0)=\Lambda_{0}$ for all $t$, then only the boundary condition at $s=L$ arises above.

From their definition (2.2), the convective variables verify the compatibility conditions

[TABLE]

Remark 2.1** (Lagrangian reduction by symmetry)**

The process of passing from the Lagrangian (or material) representation in terms of $(\Lambda,\dot{\Lambda},\mathbf{r},\dot{\mathbf{r}})$ with variational principle (2.1) to the convective representation in terms of $(\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\Gamma})$ with constrained variational principle (2.5) can be understood via a Lagrangian reduction process by symmetries. It has been carried out in [89] and is based on the affine Euler-Poincaré reduction theory of [91].

Remark 2.2** (On the functional form of the Lagrangian)**

Note that in all our theoretical considerations we will keep the Lagrangian in the general form, as we are interested in the symmetry-reduction approach to the fully three dimensional problem, rather than in the derivation of the equations of motion in a particular reduced setting, e.g., restricted to two dimensions, straight line, etc. We believe that such an approach based on Lagrangian mechanics yields the simplest possible treatment of the elastic, three-dimensional deformation of the tube.

2.2 Definition of the configuration space for the tube with expandable walls conveying fluid

We now incorporate the motion of the fluid inside the tube, the vorticity and the motion of the wall of the tube by extending the geometrically exact framework. Recall that the geometrically exact rod (without fluid) consists of a left invariant system and, therefore, can be written in terms of convective variables. On the other hand, the fluid is a right invariant system, naturally written in terms of spatial variables. The coupling of these two systems therefore yields the interesting combined fluid-structure interaction involving both convective and spatial variables but whose left and right invariances are broken by the coupling constraint. Without the vorticity, our derivation in this chapter is equivalent to [73], see also [92] for a more detailed and pedagogical exposition.

In addition to the rod variables $(\Lambda,\mathbf{r})\in\mathcal{F}([0,L],SO(3)\times\mathbb{R}^{3})$ considered above, the configuration manifold for the fluid-conveying tube also contains the Lagrangian description of the fluid. It is easier to start by defining the back-to-label map, which is an embedding $\psi:[0,L]\rightarrow\mathbb{R}$, assigning to a current fluid label particle $s\in[0,L]$ located at $\mathbf{r}(s)$ in the tube, its Lagrangian label $s_{0}\in\mathbb{R}$. Its inverse $\varphi:=\psi^{-1}$ gives the current configuration of the fluid in the tube. A time dependent curve of such maps thus describes the fluid motion in the tube, i.e.,

[TABLE]

We now include the motion of the wall of the tube, as a reaction to the fluid motion and pressure. In order to incorporate this effect in the simplest possible case of a tube with a circular cross-section, let us consider the tube radius $R(t,s)$ to be a free variable. In this case, the Lagrangian depends on $R$, as well as on its time and space derivatives $\dot{R}$ and $R^{\prime}$, respectively. Let us assume that $R$ can lie on an interval $I_{R}$, for example, $I_{R}=\mathbb{R}_{+}$ (the set of positive numbers). As we will see below, the presence of the swirl will add an additional scalar function to the variables of the Lagrangian, namely, the circulation of fluid velocity measured at a given cross-section. Then, the configuration manifold for the fluid-conveying tube is given by the infinite dimensional manifold

[TABLE]

Note that the domain of definition of the fluid motion $s=\varphi(t,s_{0})$ is time dependent, i.e., we have $\varphi(t):[a(t),b(t)]\rightarrow[0,L]$, for $\varphi(t,a(t))=0$ and $\varphi(t,b(t))=L$. The time dependent interval $[a(t),b(t)]$ contains the labels of all the fluid particles that are present in the tube at time $t$.

2.3 Definition of the Lagrangian

Let us now turn our attention to the derivation of the Lagrangian of the fluid-conveying geometrically exact tube with expandable walls. The Lagrangian will be computed as the sum of kinetic energy of rod and fluid and negative of the potential energy of the rod’s deformation.

Kinetic energy of the rod.

The kinetic energy of the elastic rod is the function $K_{\rm rod}$ given by

[TABLE]

where $\alpha$ is the linear density of the tube and $\mathbb{I}(R)$ is the local moment of inertia of the tube. The term $\frac{1}{2}a\dot{R}^{2}$ describes the kinetic energy of the radial motion of the tube.

We now derive the total kinetic energy of the fluid. In material representation, the total velocity of the fluid particle with label $s_{0}$ is given by

[TABLE]

where the Eulerian velocity is defined by

[TABLE]

Therefore, the kinetic energy of the translational motion of the fluid reads

[TABLE]

where $Q_{0}(s_{0})\mbox{d}s_{0}$ is the initial infinitesimal volume of fluid.

Circulation of velocity and Kelvin’s theorem

In order to go beyond the theory derived in [73], we introduce the swirling motion in the Lagrangian. We will introduce the simplest model of swirling. We assume that the shape of the tube at cross-section is approximately circular, and there is an axisymmetric swirling motion in addition to the streamline velocity $u(s,t)$. The trajectories of the fluid particles are locally helical, although they could be quite complex globally, i.e. extended for long $t$.

Take $C_{\alpha}$ to be a contour that circles a given cross-section close to the boundary at a Lagrangian point $\alpha$. Assuming that the deformations of the centerline remain small, the cross-section roughly circular, and the viscosity can be neglected. the cross-section remains normal to the centerline. Under this assumption, $C(t)$ will always be close to the boundary of a cross-section at all other times $t$ at the point $s=\varphi(\alpha,t)$. For an inviscid incompressible fluid considered here, Kelvin’s theorem states that the circulation $\lambda(s,t)$ is conserved along the flow, in other words, for a three-dimensional fluid velocity $\boldsymbol{v}$,

[TABLE]

where $D/Dt$ is the full derivative of the integral quantity. A more general version of Kelvin’s theorem uses the fluid momentum instead of [93]. In what follows we shall use $\lambda$ as an additional dynamic quantity in the Lagrangian to model the swirling motion of the flow. We thus obtain from (2.12) the following approximate conservation law for the swirling motion.

[TABLE]

Remark 2.3** (On the description of swirling motion through vorticity)**

One could be tempted to introduce a ’typical’ vorticity in the streamwise direction $\xi(s,t)$ instead of circulation $\lambda(s,t)$ introduced above. Then, the vorticity advection written in the local coordinates will state

[TABLE]

If the velocity field is assumed to be smooth in every cross-section, then we can write $\lambda\simeq\xi Q$. Differentiating $\lambda$ defined that way, we obtain exactly (2.13).

While this calculation is formally correct, it is difficult to justify it physically, since the swirling motion in inviscid fluids may have a singularity in the velocities at the centerline of the tube. Thus, the concept of a ’typical’ vorticity for a cross-section is ill-defined even for the simplest problem of a perfectly cylindrical straight tube with the flow independent of $s$. In that case, for every cross-section, the problem reduces to the problem of the motion of a point vortex in a circle inviscid fluid, and a point vortex at the center with a given circulation value is a solution [94, 95]. However, Kelvin’s theorem is still valid, even when the vorticity is infinite at a set of points, as long as the circulation of the vorticity is finite. Thus, we believe that the circulation $\lambda(s,t)$ defined by (2.12) with the evolution equation given by (2.13) is superior characterization of swirl as compared to the notion of a ’typical’ vorticity. Note that physical cut-offs must be introduced to utilize the notion of the kinetic energy, as, technically speaking, the kinetic energy of the fluid driven by a point vortex is infinite. **

The rotational component of kinetic energy can be then approximated by integrating the rotational component of velocity over a cross-sectional area, and then over the whole tube. Note that dimensionally, $\lambda^{2}\simeq R^{2}v_{\rm rot}^{2}$, so $\rho\lambda^{2}$ already has the *dimensions * of the kinetic energy of rotation per unit length of the tube. The exact value of the kinetic energy depends on the velocity profile of the radial velocity as a function of $r$, namely, $v_{r}(r,s,t)$. In general, this velocity profile is a complex function of the tube’s shape and the local vorticity profile of the fluid. We thus define the profile shape function for every cross-section $C_{s}$

[TABLE]

In general, the shape function $\Phi$ is dimensionless with values of order 1. It describes the deviation of the velocity profile from its given equilibrium value during the dynamics. The profile shape function $\Phi$ can only be obtained through either numerical simulations or experiments, and is thus a modeling component of the theory. Our equations of motion will be derived by the variational principle and will be valid for any $\Phi$. While the exact nature of $\Phi$ and its dependence on the flow is to be investigated in future work, although some conclusions can be drawn from simple considerations. For example, for a point vortex on the plane with circulation $\lambda$, the rotational velocity is $v_{\rm rot}=\frac{\lambda}{2\pi r}$, so the kinetic energy of the fluid’s rotation in a circle of radius $R$ is infinite. However, if we introduce a smoothing cut-off of radius $\epsilon$ around the vortex attributed to viscosity, we arrive to

[TABLE]

In the numerical simulations in this paper we will use (2.16) for a fixed value of $\epsilon$. Then,

[TABLE]

Using (2.10) together with the change of variables $s=\varphi(t,s_{0})$, we can rewrite $K_{\rm fluid}$ as

[TABLE]

where $\rho$ is the mass density of the fluid per unit volume, in units Mass/Length3, and $Q(t,s)$ is the area of the tube’s cross section, in units Length2.

Elastic energy.

The potential energy due to elastic deformation is a function of $\boldsymbol{\Omega}$, $\boldsymbol{\Gamma}$ and $R$. While the equations will be derived for an arbitrary potential energy, we shall assume the simplest possible quadratic expression for the calculations, namely,

[TABLE]

where $\boldsymbol{\chi}\in\mathbb{R}^{3}$ is a fixed vector denoting the axis of the tube in the reference configuration, $\mathbb{J}$ is a symmetric positive definite $3\times 3$ matrix, which may depend on $R$, $R^{\prime}$ and $R^{\prime\prime}$, and $\lambda(R)$ is the stretching rigidity of the tube. The stretching term proportional to $\lambda(R)$ can take the more general form $\mathbb{K}(\boldsymbol{\Gamma}-\boldsymbol{\chi})\cdot(\boldsymbol{\Gamma}-\boldsymbol{\chi})$, where $\mathbb{K}$ is a $3\times 3$ tensor. The part of this expression for the elastic energy containing the first two terms in (2.18) is commonly used for a Cosserat elastic rod, but more general functions of deformations $\boldsymbol{\Gamma}$ are possible, see [50] for a more detailed discussion of possible forms of the potential energy. A particular case is a quadratic function of $\boldsymbol{\Gamma}$ leading to a linear dependence between stresses and strains. We have also introduced the elastic energy of wall $F(R,R^{\prime},R^{\prime\prime})$ which can be explicitly computed for simple elastic tubes. In general $F$ depends on higher derivatives, such as $R^{\prime\prime}$. We shall derive the equations of motion in their general form for arbitrary $F$, but for particular examples used in simulations we will use $F(R,R^{\prime})$ for simplicity.

Mass conservation.

Before we write the final expression for the Lagrangian, let us discuss the question of the mass conservation since it will be used as a constraint to derive the equations of motion, with the appropriate Lagrange multiplier playing the role of fluid’s pressure. We shall assume that the fluid fills the tube completely, and the fluid velocity at each given cross-section is aligned with the axis of the tube. Since we are assuming a one-dimensional approximation for the fluid motion inside the tube, the mass density per unit length $\rho Q=\rho A|\boldsymbol{\Gamma}|$, has to verify

[TABLE]

Since $\rho=$const, and the fluid is incompressible, we deduce the conservation law

[TABLE]

The origin of this conservation law is illustrated on the left panel of Figure 2.1. A Lagrangian point $a$ is mapped to $s=\varphi(a,t)$, and the Lagrangian point $a^{\prime}$ is mapped to $s^{\prime}=\varphi(a^{\prime},t)$. Inversely, a point on a centerline $s$ will be mapped to the initial Lagrangian point $a=\varphi^{-1}(s,t)$. Taking $a^{\prime}=a+\mbox{d}a$, we see that the infinitesimal physical volume is described at the physical space by $Q\mbox{d}s$, and the corresponding volume at $t=0$ by $Q_{0}(a)\partial_{s}\varphi^{-1}(s,t)$, with $a=\varphi^{-1}(s,t)$, giving exactly (2.20).

The key new feature of this Section is in the introduction of the swirling motion in the kinetic energy part of the Lagrangian. The particle trajectories for such a swirling motion are approximately helical, as illustrated by a sketch on the right panel of Figure 2.1.

Lagrangian.

From all the expressions given above, we obtain the Lagrangian of the fluid-conveying given by

[TABLE]

and defined on the tangent bundle $T\mathcal{Q}$ of the configuration space $\mathcal{Q}$, see (2.9). Note that all the arguments of $\mathsf{L}$ are functions of $s$, so we don’t need to include the spatial derivatives of $(\Lambda,\mathbf{r},R)$ explicitly as variables in $\mathsf{L}$. These spatial derivatives appear explicitly in the expression of the integrand of the reduced Lagrangian (2.22) below. Assuming there is a uniform external pressure $p_{\rm ext}$ acting on the tube, the Lagrangian expressed in terms of the variables $\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\Gamma},u,Q,\lambda$ reads

[TABLE]

where $S(R)$ is the function describing the potential energy of tube’s stretching along its centerline as a function of $R$. Our derivation of the equations of motion (2.26) below will be valid for the Lagrangian $\ell$ being an arbitrary function of its variables. Of course, the numerical simulations will require specifying the exact functional form of the Lagrangian $\ell$ for every particular case.

2.4 Variational principle and equations of motion

The equations of motion are obtained from the Hamilton principle applied to the Lagrangian (2.21), namely

[TABLE]

for arbitrary variations $\delta\Lambda,\delta\mathbf{r},\delta\varphi,\delta R$ vanishing at $t=0,T$. In terms of the symmetry reduced Lagrangian $\ell$, this variational principle becomes

[TABLE]

for variations (2.3), (2.4), and $\delta u$, $\delta\xi$ and $\delta\lambda$ computed as

[TABLE]

where $\eta=\delta\varphi\circ\varphi^{-1}$. Note that $\eta(t,s)$ is an arbitrary function vanishing at $t=0,T$. Note that the Lagrangian function in (2.24) will contain the incompressibility constraint (2.20) with the Lagrange multiplier $p$, which we drop from the notation in (2.24) for simplicity. A lengthy computation (see [73] for details) yields the system

[TABLE]

where the symbols ${\delta\ell}/{\delta\boldsymbol{\omega}},{\delta\ell}/{\delta\boldsymbol{\Gamma}},...$ denote the functional derivatives of $\ell$ relative to the $L^{2}$ pairing, and we introduced the notations

[TABLE]

Note that the first equation arises from the terms proportional to $\boldsymbol{\Sigma}$ in the variation of the action functional and thus describes the conservation of angular momentum. The second equation arises from the terms proportional to $\boldsymbol{\psi}$ and describes the conservation of linear momentum. We remind the reader that the quantities $\boldsymbol{\Sigma}$ and $\boldsymbol{\psi}$ are defined by $\boldsymbol{\Sigma}=(\Lambda^{T}\delta\Lambda)^{\vee}$ and $\boldsymbol{\psi}=\Lambda^{T}\delta\mathbf{r}$, with $(\boldsymbol{\Sigma}(s,t),\boldsymbol{\psi}(s,t))$ for fixed $s$ and $t$ taking values in the Lie algebra $\mathfrak{se}(3)$.

The third equation is obtained by collecting the terms proportional to $\eta$ and describes the conservation of fluid momentum. The fourth equation comes from collecting the terms proportional to $\delta R$ and describes the elastic deformation of the walls due to the pressure. Finally, the last four equations arise from the four definitions $\boldsymbol{\Omega}=(\Lambda^{-1}\Lambda_{s})^{\vee}$, $\boldsymbol{\Gamma}=\Lambda^{-1}\mathbf{r}^{\prime}$, $Q=(Q_{0}\circ\varphi^{-1})\partial_{s}\varphi^{-1}$, and $\lambda=\lambda_{0}\circ\varphi^{-1}$.

We will derive the equations for general $Q=Q(\boldsymbol{\Omega},\boldsymbol{\Gamma},R)$, keeping in mind the fully three-dimensional dynamics, but will only use it for the tubes satisfying $Q=\pi R^{2}|\boldsymbol{\Gamma}|$ which is a good approximation for mostly straight tubes. For the choice

[TABLE]

the functional derivatives are computed as

[TABLE]

where $\partial\ell_{0}/\partial\boldsymbol{\Omega},\partial\ell_{0}/\partial\boldsymbol{\Gamma}$, … denote the ordinary partial derivatives of $\ell_{0}$, whose explicit form can be directly computed from the expression of $\ell_{0}$ in (2.22).

Theorem 2.4

For the Lagrangian $\ell$ in (2.22), the variational principle (2.24) with constrained variations (2.3), (2.4), (2.25) yields the equations of motion

[TABLE]

together with appropriate boundary conditions enforcing vanishing of the variations of the boundary terms.

Proof. The system is obtained by replacing the expression of the functional derivatives (2.27) in the system (2.26).

Note that the term $\frac{\delta\ell_{0}}{\delta R}$ contains both the terms coming from the elastic energy and the derivatives of the shape function $\Phi(R)$. Thus, the vorticity contributes to the pressure term and momentum balance of the wall, as expected.

Remark 2.5** (On the formal correspondence between the entropy and circulation)**

One can notice the formal correspondence between (2.26) and the corresponding equations for the incompressible fluid with thermal energy derived in [73]. Mathematically, the appearance of equations is equivalent when the circulation $\lambda$ is associated with entropy $S$ in [73], which is due to the same way the circulation and entropy are advected. That association is purely mathematical: there is no physical correspondence between these quantities, and they enter the Lagrangian in a very different way. The circulation enters the kinetic energy part of the Lagrangian, and entropy enters the internal energy part. Thus, in spite of the apparent formal similarity between these equations, the use of any particular physically relevant Lagrangian leads to very different equations for the fluids containing swirling vs thermal energy. **

Equations (2.28) have to be solved as a system of nonlinear partial differential equations, since all equations are coupled through the functional form of the Lagrangian (2.22). This can be seen, for example, from computing the derivative

[TABLE]

which appears in the the balance of fluid momentum, i.e., the third equation in (2.28). Thus, the fluid momentum involves the radius, stretch $\boldsymbol{\Gamma}$, and linear velocity $\boldsymbol{\gamma}$, as well as the fluid’s velocity.

If we define the rescaled fluid momentum as

[TABLE]

the third equation in (2.28) can be (deceptively) simply written as

[TABLE]

which is reminiscent of the equations of motion for one-dimensional fluid with the additional circulation terms describing the evolution of vorticity.

For $Q(R,\boldsymbol{\Gamma})=A(R)|\boldsymbol{\Gamma}|$, the equations (2.28) can be further simplified since

[TABLE]

Notice also that for $A(R)=\pi R^{2}$, we have $\mbox{d}A/\mbox{d}R=2\pi R$, the circumference of a circle.

3 Tubes with a straight centerline: dynamics and reduction to a single equation

In this section we shall specify our model to the case of an inextensible and unshearable expandable tube. We then focus on straight expandable tubes with no rotational motion and compare our model with previous works. This case has been studied extensively in the literature, mostly in the context of blood flow involving incompressible fluid. We also show that these simplified models arise from a variational principle using the back-to-labels map. The explicit use of back-to-labels map allows to derive a single equation of motion without the need for the pressure-like Lagrange multiplier term $p$.

3.1 Equations of motion for inextensible and unshearable tubes

The dynamics of an inextensible and unshearable, but expandable, tube conveying compressible fluid is obtained by imposing the constraint $\boldsymbol{\Gamma}(t,s)=\boldsymbol{\chi}$, for all $t,s$. If we denote by $\mathbf{z}$ the Lagrange multiplier associated to this constraint and add the term $\int_{0}^{T}\int_{0}^{L}\mathbf{z}\cdot(\boldsymbol{\Gamma}-\boldsymbol{\chi})\mbox{d}s\,\mbox{d}t$ to the action functional in our variational principle (2.24), the first two equations in (2.28) will change to

[TABLE]

The physical meaning of $\mathbf{z}$ is the reaction force enforcing the inextensibility constraint.

Particular simple solutions of this system can be obtained by assuming the axis of the tube being straight and no rotational motion. In that case, the inextensibility constraint leads to $\mathbf{z}=z\boldsymbol{\chi}$, the direction along the axis. For such particular solutions, the angular momentum equation is satisfied identically, and the linear momentum equation reduces to a one dimensional equation for the reaction force $z$. Therefore, we only need to compute the fluid momentum equation and the Euler-Lagrange equations for $R$. From the third and fourth equations in (2.28), we get for $A=A(R)=\pi R^{2}$:

[TABLE]

together with the conservation of mass and vorticity

[TABLE]

This gives a system of four equations for the four variables $u$, $R$, $\lambda$, and $p$.

It is often more convenient to use the cross-sectional area $A$ instead of the radius $R$, and express the Lagrangian in terms of $A$ and its derivatives, as well as $u$ and $\lambda$. To explicitly denote the change of variables, we shall denote that Lagrangian as $\ell_{A}$, and the Lagrangian function as $\ell_{0,A}$, so $\ell_{A}=\int\ell_{0,A}\mbox{d}s$. The equations of motion are then:

[TABLE]

Without the vorticity, (3.3) represent the familiar model of tube with elastic walls incorporating the wall inertia [56, 57]. For example, one can consider the following form of the Lagrangian function $\ell_{0,A}:$

[TABLE]

where $K_{A}(A,A_{t})$ is the kinetic energy of the wall motion expressed in terms of $A$ and its time derivatives, $F_{A}(A,A_{s})$ is the elastic energy function expressed in terms of $A$ and its spatial derivatives, and $\Phi_{A}$ is the vorticity shape function $\Phi(R)$ expressed in terms of $A$. If the dependence of the Lagrangian function $\ell_{0,A}$ on $A_{t}$ and $A_{s}$ is neglected, the second equation of (3.3) yields the so-called wall pressure law giving an explicit connection between pressure and area, see [72] and references therein. This wall pressure law $p(A)$ can then be substituted in the first equation of (3.3), and, in conjunction with the conservation law (the first equation of (3.2)) allows to obtain a closed set of two coupled equations for $A$ and $u$. We shall not drop the inertia terms in (3.3) in what follows and pursue the more general case of a general Lagrangian dependence on $A$ and its time and $s$-derivatives.

3.2 The case of uniform circulation, energy, and the new constant of motion

Let us now analyze the equations of motion (3.3) in the case when $\lambda=\lambda_{*}=$const. Then, $\lambda=\lambda_{*}$ everywhere in the flow. Physically, such system can be realized by producing a constant vortex circulation at tube’s entry using a mechanical device, e.g. an idealized propeller. In that case, we can define $\overline{\ell}_{0,A}=\ell_{0,A}$ with $\lambda$ set to be equal to the value of $\lambda_{*}$. Equations of motion for the general Lagrangian function $\overline{\ell}_{0,A}(A,A_{t},A_{s},u)$ are then given as

[TABLE]

One would expect that there is an energy-like quantity $E$ which is conserved on solutions of (3.5). Indeed, we can prove the following

Lemma 3.1** (Energy conservation)**

The quantity

[TABLE]

is conserved on solutions of (3.5). **

Proof. This result is proven by computing the time derivative of $E$ and integrating by parts, using the equations of motion (3.5), assuming that the boundary terms vanish when integrating by parts, and utilizing the constraint, i.e., the last equation of (3.5). This proof closely follows that of Lemma 3.4 below, and it is relatively straightforward, so we do not present is here.

Remark 3.2** (On the physical meaning of $E$)**

One can notice that the quantity $E$ defined by (3.33) is, in general, not equal to the kinetic plus potential energy (also sometimes called the ’physical’ energy). For example, for the Lagrangian defined by (3.4),

[TABLE]

where as the kinetic plus potential energy will be given by the expression

[TABLE]

The energy defined by (3.7) is conserved, whereas the quantity defined by (3.8) is not. Changing time derivatives to the corresponding momenta in (3.7) will give the Hamiltonian of the system. The non-conservation of ’physical’ energy, i.e., the kinetic plus potential energy, is typical for systems with rotation, see [96]. We believe it is difficult to guess the rather non-intuitive form (3.7) without appealing to the variational principle used in the derivation of (3.31). **

What is rather surprising, however, is that one can find there is another constant of motion. We prove the following

Theorem 3.3** (Additional constant of motion)**

The following quantity is conserved on the solutions of (3.5)

[TABLE]

Proof. Differentiating the quantity defined by (3.9) with respect to time, we obtain:

[TABLE]

In the second line of (3.10) we have used equations (3.5) for time derivatives of momenta, the third line is obtained by the integration by parts and neglecting the boundary terms, and noticing that the expression under the integrand in the last line is (minus) the full derivative of the quantity $\overline{\ell}_{0,A}+(p-p_{\rm ext})A$ with respect to $s$.

It is hard to readily assign a physical relevance to the quantity (3.9), beyond it being a linear combination of the wall and fluid momenta. The expression is valid for an arbitrary Lagrangian functions, and, apart from the kinetic energy coming from the rotation and being proportional to the parameter $\lambda_{*}^{2}$, which can be set to [math], is quite general to describe a majority of particular choices of kinetic and potential energy for models in the literature, see e.g. [69, 72] and references within. It was thus rather surprising to us that in spite of a large body of literature on the topic, we have failed to find any mention of the conserved quantity (3.9). It is probably due to the fact that the physical meaning of the constant of motion (3.9) is not apparent, and one can hardly guess it from the equations of motion written as the balance of mass and momenta conservation using considerations from the general physical principles.

We now turn our attention to the alternative description of motion utilizing the back-to-labels map. As it turns out, this description is highly beneficial from the mathematical point of view, as it allows to combine several conservation laws into a single equation.

3.3 Equations of motion for the back-to-labels map

Let us now show how to reduce the the fluid and wall momenta equations, and the conservation laws, to a single differential equation using the back-to-labels map. The procedure of this method can be explained as follows. Let us denote, for shortness, $\psi(s,t)=\varphi^{-1}(s,t)$ to be the back-to-labels map, giving the Lagrangian label of a fluid particle at the point $s$ at time $t$. For a given initial cross-section $A_{0}(a)$, the pressure term is entering the equations as the Lagrange multiplier for the constraint $A(s,t)=A_{0}\circ\psi\partial_{s}\psi=A_{0}(\psi)\partial_{s}\psi$. Instead of enforcing this holonomic constraint using the Lagrange multiplier $p$, in this Section, we will use $\psi(s,t)$ as the free variable, express all the other variables, such as the cross-section, in terms of $\psi$ in the Lagrangian, and take the variations with respect to $\psi$. It will allow us to compute the equations of motion as the explicit Euler-Lagrange equations in terms of variables $\psi$, and not the Euler-Poincaré equations as we used in (2.28).

We can define the function $\Psi(\alpha)$ as the antiderivative of $A_{0}(\alpha)$ and rewrite the conservation law as

[TABLE]

In general, $\Psi(\alpha)$ is a given function given by the initial profile of the cross-section. Since $A_{0}(\alpha)>0$, $\Psi(\alpha)$ is a monotonically increasing function of $\alpha$. We then use (3.11) to rewrite the conservation law

[TABLE]

Setting the integration constant $C(t)=0$ in the above formula, for example, due to the boundary conditions, we find the expression of $u$ in terms of the back-to-labels map $\psi$ as

[TABLE]

The Lagrangian $\ell_{A}$ and the corresponding action $S_{A}$ for an inextensible, unshearable tube with the straight centerline is written as

[TABLE]

for example, see the explicit version of $\ell_{A}$ given in (3.4). We have denoted, as before, the Lagrangian itself as $\ell_{A}$, which involves the integral over the $s$-domain. The function that is being integrated, i.e., the integrand of the Lagrangian, is denoted with a subscript [math], i.e., $\ell_{0,A}$.

Since, according to (2.20), $A=A_{0}(\psi)\partial_{s}\psi$, using (3.13) and the constraint $A=A_{0}(\psi)\partial_{s}\psi$, we can express the Lagrangian for the tube (3.14) as a function of $\psi$ only, and derive the equations of motion as a single Euler-Lagrange equation obtained by variations of the action with respect to $\psi$. Let us illustrate this method in more details, including the vorticity in the tube.

In order to explicitly separate the effect of fluid’s motion, let us define $b=\psi(s,t)-s$, since for the absence of fluid’s motion $\psi(s,t)=s$. Then, from (3.13), the fluid velocity is expressed as

[TABLE]

In what follows we consider, for simplicity, the initial cross-section to be constant and the tube being initially uniform, i.e, $A_{0}(\alpha)=A_{0}=$const, corresponding to $\Psi(\alpha)=\alpha$. In principle, we can extend the theory for an arbitrary initial cross-section $A_{0}(\alpha)>0$, although this will introduce unnecessary algebraic complexity into equations.

Instead of the full cross-section $A$, it is more convenient to work with the relative deviation of the cross-section from the equilibrium, as this approach leads to the most algebraically simple equations. More precisely, let us define $a$ to be the relative deviation from the equilibrium value according to $A=A_{0}(1+a)$. Then, the constraint $A=A_{0}\circ\varphi^{-1}\partial_{s}\varphi^{-1}$ becomes

[TABLE]

In general, we will have $\ell_{a}=\int\ell_{0,a}(a,a_{t},a_{s},u,\lambda)\mbox{d}s$. Just as in (3.14), we have dropped the dependence of the Lagrangian on $a_{ss}$, the inclusion of this dependence is trivial and only leads to unnecessary algebraic complexity. We re-define the Lagrangian using (3.16), by substituting $a=b_{s}$, $a_{t}=b_{ts}$, $a_{s}=b_{ss}$ into $\ell_{a}$, and further employing (3.15) to express the velocity in terms of derivatives of $b$:

[TABLE]

We then take variations of the corresponding action defined as $S_{b}=\int\ell_{b}\mbox{d}s\mbox{d}t$ We notice that the variations of the velocity in terms of $\delta b$ are given as

[TABLE]

and hence

[TABLE]

Starting with the Lagrangian $\ell_{a}=\int\ell_{0,a}(a,a_{t},a_{s},u,\lambda)\mbox{d}s$, and remembering that $a=b_{s}$, we compute the Euler-Lagrange equations from the critical action principle as

[TABLE]

Hence, since $\delta b$ is arbitrary, the Euler-Lagrange equation corresponding to the variations with respect to $b$ gives

[TABLE]

We can simplify (3.21) further and bring it to a more compact form. Notice that due to the conservation law for volume (3.2) for $A=A_{0}(1+a)=A_{0}(1+b_{s})$, we have

[TABLE]

see also (3.12) for more general derivation of a non-constant initial cross-section $A_{0}$. In the expression above, we will set the integration function $C(t)=0$ assuming appropriate boundary conditions. This is the simplest choice we shall use in the remainder of the paper. Multiplying equation (3.21) by $(1+b_{s})$ and using (3.21), we can write the complete system as

[TABLE]

We can now make approximations in the Lagrangian directly, and choosing any physically relevant Lagrangian function $\ell_{0,a}$ in (3.23) will produce a consistent system with the Lagrangian structure. To be concrete, let us take the potential energy of the tube’s wall to be $\frac{1}{2}(Pa_{s}^{2}+Qa^{2}+F(a))$, where $F(a)$ is a smooth function describing the higher-order terms in the elastic potential, for example, $F(a)=Q_{2}a^{4}$. Thus, $F(a)=o(a^{2})$ as $a\rightarrow 0$. The first two terms of this potential energy are in fact the elastic part of the potential energy used in the generalized string model, see e.g. [56, 57], which is well-accepted in the literature. The last term $F(a)$, corresponding to the nonlinear elasticity, is one way to stabilize the dynamics in the unstable regime, as we shall see below. For further calculations, we take the Lagrangian to be

[TABLE]

This Lagrangian is, in essence, a particular version of (3.4) reformulated in terms of relative deviation of the area $a$ instead of the area $A$ itself. The first term in (3.24) is the kinetic energy of the wall deformation; the second and third terms are the potential energy of the wall deformation. The fourth term is the kinetic energy of the fluid $\frac{1}{2}\rho Au^{2}$ and the last term is the kinetic energy of fluid rotation, with $\Phi(a)$ being a dimensionless quantity (shape function) of order 1. Technically, $K=K(a)$, $P=P(a)$ and $Q=Q(a)$ are functions of $a$, however, in the approximation (3.24) we have put them to be constants for simplicity since we assumed $a$ to be sufficiently small during the evolution. In that case, the equations (3.23) become

[TABLE]

Notice that we do not have a pressure appearing in (3.25) since the constraint (3.16) is enforced explicitly by substitution in the Lagrangian. The fluid momentum equation, i.e., the first equation of (3.25), now becomes a single equation for the quantity $b$, which is coupled with the evolution equation for $\lambda$. After some algebra, we obtain

[TABLE]

System (3.26) defines two coupled PDEs for the unknowns $b$ and $\lambda$. It is interesting to compare this system with the one obtained by the Euler-Poincaré method (3.3), which explicitly involves pressure, and has to be solved as a system of four equations for $(u,\lambda,A,p)$. Notice also that (3.3) does not define an evolution equation for pressure $p$, which has to be found in order to enforce the incompressibility. In our opinion, equation (3.26) is advantageous for both theoretical and numerical studies, since there is no need to compute the pressure term.

3.4 The case of constant vorticity

Let us now consider the case of Section 3.2, where the circulation generated by the source at the origin of the tube is held constant by a mechanical device, so $\lambda=\lambda_{0}$. The general Lagrangian with $\lambda=\lambda_{*}=$const in (3.21) gives

[TABLE]

where, as before, we have denoted $\overline{\ell}_{0,a}=\ell_{0,a}$ to be the Lagrangian function $\ell_{0,a}$ evaluated at $\lambda=\lambda_{0}$. Similarly, for a particular choice of the Lagrangian $\overline{\ell}_{0,a}$ given by (3.24), the last equation of (3.26) describing the evolution of $\lambda$ is satisfied identically, and the equation reduces to a single PDE for the variable $b(s,t)$. We shall mostly analyze a particular form (3.26) in what follows; however, the conservation of energy we will derive will be valid for an arbitrary functions $\overline{\ell}_{0,a}$.

If we choose the relevant time and length scales $T$ and $L$ respectively, and defined the non-dimensionalized variables according to

[TABLE]

We will choose the length and time scales so that $K_{*}=1$ and $P_{*}=1$ in (3.28), although other choices of time and length scales are also possible. Combining the coefficients containing the background vorticity $\lambda_{*}$ in a single coefficient $\zeta$ defined as

[TABLE]

For example, for the cross-section approximation of the point vortex (2.16), we obtain

[TABLE]

Keeping the notations $b$, $t$ and $s$ for the non-dimensionalized $b$, time and $s$, and performing some algebra, we obtain

[TABLE]

Equation (3.31) gives, as far as we know, a new result expressing the evolution of the collapsible tube with compliant walls carrying vorticity in terms of back-to-labels mab $b(s,t)$. While it is difficult to observe the variable $b$ in an experiment, its spatial derivative $a=b_{s}$ can readily be observed. The rest of the paper is dedicated to the analysis of (3.31), both in general case, and in the asymptotic approximations of long wavelengths and large times.

We start with the analysis of the full equation (3.31) and its general counterpart (3.21) first.

3.5 Energy conservation and additional constant of motion

Energy conservation

Let us define $\ell_{0,b}$ to be $\overline{\ell}_{0,a}$, i.e., the Lagrangian function $\overline{\ell}_{0,a}$ with the substitution $a=b_{s}$, $a_{t}=b_{st}$, $a_{s}=b_{ss}$. Since (3.27) is derived as an Euler-Lagrange equation for the Lagrangian

[TABLE]

it conserves the energy-like function as stated in the following

Lemma 3.4** (Energy conservation)**

For appropriate periodic boundary conditions, the equation (3.31) conserves the quantity

[TABLE]

Proof. We prove this theorem by computing the evolution of energy density $e$ defined by (3.33).

[TABLE]

The expression in the square brackets is exactly the Euler-Lagrange equation for the Lagrangian (3.32), so it vanishes identically. Thus, we have $\partial_{t}e+\partial_{s}J_{e}=0$ where the energy flux $J_{e}$ is defined by the equation in the parenthesis on the last line of (3.34). Thus, $\int e_{t}\mbox{d}s$ vanishes for periodic boundary conditions of the conditions where $J_{e}=0$ on the boundaries, for example, $s\in(-\infty,\infty)$ with all variables tending to [math] as $s\rightarrow\pm\infty$.

In the case when at the undisturbed state the fluid is moving with a constant velocity $U$, it is is convenient to write (3.31) in the form $b(s,t)=-Ut+B(s,t)$ so $B(s,t)=0$ if there are no deviations of the flow from the steady flow $u=U$. Then, (3.31) becomes

[TABLE]

In addition to the energy conservation, one can prove the appropriate equivalent of (3.9) reformulated in terms of the back-to-labels map $B$. Namely, we have the following

Lemma 3.5** (Additional constant of motion)**

For periodic boundary conditions, or boundary conditions where $B\rightarrow 0$ when $s\rightarrow\pm\infty$, equation (3.35) conserves the following quantity:

[TABLE]

Proof. This can be seen by re-writing (3.35) in the form $\partial_{t}I+\partial_{s}J$ for some quantity $J$ depending on $B$, $B_{t}$, $B_{s}$, $B_{ss}$ etc.. The result is valid for arbitrary functions $\zeta(a)$ and $f(a)$, and arbitrary velocity $U$.

3.6 Linear stability and instability

Let us now turn our attention to the linear stability analysis of equation (3.31). Since we assume that $f(b_{s})=o(b_{s})$, we can write the linearized version of that equation, by assuming small amplitude of $b$ with moderate values for $t$- and $s$-derivatives. For convenience, we also define $\zeta_{0}:=\zeta(0)$ to be the linearized value of $\zeta$ at $b=0$. Then, (3.31) gives

[TABLE]

Note that on the left of that equation, we obtain exactly the Helmholtz operator acting on $b$. On the right of that equation, we have a Helmholtz operator for $\zeta_{0}>0$ and anti-Helmholtz operator for $\zeta_{0}<0$. We will consider cases of $\zeta_{0}$ being possibly positive or negative, keeping in mind that for physical reasons, $\zeta_{0}>0$ due to the expression (3.30). While the systems values $\zeta_{0}<0$ may exist, their physical origin still remains uncertain, and thus all our numerical solutions will be computed for $\zeta_{0}>0$ (3.30) unless explicitly stated.

To compute the dispersion relation for no external flow, assume $b=b_{0}e^{i(ks-\omega t)}$ and substitute into (3.37). We obtain:

[TABLE]

For large $k$, the waves propagate with the constant speed (phase velocity) of $\pm 1$. If $\zeta\geq 0$, all wavelengths are stable. If $\zeta_{0}<0$, then there is a band of instability $0<k<\sqrt{-\zeta}$ where the waves in that band grow exponentially. Long wavelengths $k\rightarrow 0$ propagate with the phase and group velocities $\sqrt{\zeta_{0}}$ and short wavelengths $k\rightarrow\infty$ propagate with the speed $1$.

Let us now compute the stability for the case when the base flow has the circulation $\lambda_{*}$ and velocity $U$ in the dimensionless variables. The linearization of the terms in the equation (3.35) proceeds as follows:

[TABLE]

Then, the linearization of (3.31) up to order $\epsilon$ is given by

[TABLE]

Substitution of $B=e^{ikx-i\omega t}$ gives the dispersion relation connecting $\omega$ and $k$

[TABLE]

We can thus solve for $\omega=\omega(k)$ as

[TABLE]

We want to find the areas of stability and instability, i.e., the values of parameters $U$ and $\zeta_{0}$ such that ${\rm Im}\omega=0$ (stability) or ${\rm Im}\omega>0$ (instability). Because of the $\pm$ sign in (3.42), it is sufficient to find the values of parameters when $\omega$ becomes complex for instability. The square root (discriminant) in (3.42) is non-negative if and only if the quadratic polynomial $f(z)=-U^{2}z+(1+z)(\zeta_{0}+z)$ is non-negative for all $z=k^{2}\geq 0$. Therefore, the system is unstable if and only if

[TABLE]

Since the minimum of the function on the right-hand side of (3.43) is achieved at $k^{2}=\sqrt{\zeta_{0}}$ (we remind that $\zeta_{0}$ is assumed to be positive), then the condition of instability is given as

[TABLE]

For values of $U$ satisfying (3.44), there is a band of unstable wavelengths of $k$ satisfying (3.43). When $U$ is slightly above $U_{*}$, the instability band of $k$ is centered around $k_{*}=\zeta_{0}$. Thus, at the bifurcation point $U=U_{*}$, the typical unstable wavelength is $k\sim\sqrt{\zeta_{0}}$.

We shall also note that for for the long wavelength $k\rightarrow 0$, equation (3.41) gives $\omega\sim v_{\pm}k$, where the wave speed $v_{\pm}=-U\pm\sqrt{\zeta_{0}}$. The long wavelength are thus always stable for $\zeta_{0}>0$. In the short wave limit $k\rightarrow\infty$, equation (3.41) gives $\omega=\sim\pm k$, so the short waves propagate with speeds of $\pm 1$.

The results of linear stability analysis presented above are summarized on Figure 3.1. We take $\zeta_{0}=1$ corresponding to the critical velocity $U_{*}=2$. The top and the bottom panels of this figure show, respectively, the real part of $\omega(k)$ computed from (3.42). On the left panels, we show the results with $U=1$, which is the stable case. For the stable case, ${\rm Im}\omega(k)=0$ for all $k$. The right panels show the results for $U=3$ which is the unstable case. There is a band of instability centered around $k\sim\zeta_{0}=1$ as expected from the discussion after (3.43).

3.7 Numerical simulations: stable case

The results of simulations of (3.35) with periodic boundary conditions for $a(s,t)=B_{s}(s,t)$ are presented on Figure 3.2. We have taken $Q_{*}=1$, and $\zeta(a)$ given by (3.30). For the stable case, we have also set the nonlinear terms to [math], i.e., $f(a)=0$. The parameters of the simulations on the Figure are $U=1$ and $\lambda_{*}=1$, with $\zeta_{0}\simeq 1.1592$, corresponding to the stable case. We present $a=B_{s}$ since $a(s,t)$ is the physical observable, being the deviation of the relative area from its equilibrium value.

On the left panel we present $a=B_{s}$ inferred from the solution $B(s,t)$ of the full equation (3.35) starting with a localized disturbance $B(s,0)=e^{-0.1s^{2}}$, $B_{t}(s,0)=0$. On the right panel, we present $a=B_{s}$ inferred from the solution for the Boussinesq-like approximation of the full equation (4.13), derived below in Section 4.2 with the same parameter values and initial conditions as on the left panel.

To check the accuracy of our calculations, on Figure 3.3, we also show the energy defined by (3.7) (left panel) and the conserved quantity $I$ defined by (3.36) (right panel). The energy defined in (3.7) is preserved to the relative accuracy of $10^{-4}$. For comparison, on the right panel we also present the terms $B_{t}$ and $B_{ss}B_{ts}$, (green dotted and blue dashed lines) which are of order one. The constant of motion defined by (3.36), shown as a solid red line, is preserved to expected precision, about $10^{-4}$.

If $|U|>1+\sqrt{(}\zeta_{0})$, the solution becomes unstable, and if there is no potential energy to limit the growth of the solution, it quickly reaches the singular (collapsed) state $a=-1$. In reality, any physical system will have a potential energy preventing the singularity. One method is to introduce higher derivatives in the elasticity function $F$, and another way is to consider the potential energy preventing both indefinite growth of $a$ and collapse to $a=-1$. We chose the same parameters as in (3.2) except

[TABLE]

The initial conditions at $t=0$ are taken to be a random perturbation for $B(s,0)$ with the amplitude $10^{-3}$ and $B_{t}(s,0)=0$. The results of the simulations for $a=B_{s}$ are demonstrated on Figure 3.4. The solution grows exponentially from the unstable equilibrium state $a=0$, saturating at an amplitude $|a|\sim 0.5$ and yielding complex non-steady dynamics.

4 Asymptotic analysis of equation (3.35)

4.1 Small cross section (collapse) and reduction to Monge-Ampère equation

One of the most interesting regimes of evolution occurs when the cross-section almost collapses to 0, in other words $1+B_{s}>0$ and $|1+B_{s}|\ll 1$. This regime is important for physiological flow applications and is studied quite extensively in the literature, see, for example, recent review [97] and the references therein. We show that in this regime, our equation (3.35) is approximated by the celebrated Monge-Ampère equation. We proceed as follows. For $1+B_{s}$ being small and positive, and $B_{ss}$, $B_{sstt}$, $B_{ssss}$ remaining finite, the asymptotic equation following from (3.35) is

[TABLE]

We now use the substitution

[TABLE]

With this substitution, (4.1) integrates once with respect to $s$. Dropping the arbitrary function of time appearing under this integration, we obtain the following equation for $K(s,t)$:

[TABLE]

which is the (real) homogeneous Monge-Ampère equation with $(s,t)$ being independent variables. For the approximation (4.1) to be valid, we assume $K$ to be order $1$, and since $B_{s}+1$ is small, we need that $K_{ss}$ should be small. We can thus assume that $K=K(\epsilon S,\epsilon\tau)$ where $\epsilon\ll 1$ is a small parameter and the variables $S,\tau$ are of order 1. Then the equation (4.3) transforms to

[TABLE]

i.e. satisfies the same Monge-Ampère equation. Thus, any solution $K(S,\tau)$ of (4.4) will generate a solution $K(s,t)$ describing the collapse of the tube. For physical reasons, we need $K_{SS}>0$ for all solutions since $1+B_{s}=K_{ss}>0$. For that solution, the neglected terms in (3.35) are order $\epsilon^{3}$ or smaller in order than the terms kept in (4.1), so the approximation is consistent.

As is known from the literature, equation (4.4) possesses a rich structure of solutions and is integrable [98, 99, 100]. For example, for any smooth function $g(x)$ with $g^{\prime\prime}(x)>0$, $K(S,\tau)=g(S-c\tau)$ is a solution of (4.4) for any constant $c$. Thus, the collapse of the tube is strongly dependent on the initial conditions and is unlikely to be described by a universal law.

4.2 Approximation by Boussinesq-type equations

In this section, we are interested in a consistent approximation of (3.35) for small amplitude and long wavelength approximation which is a common approximation for hydrodynamic system with moving boundaries. Perhaps the most celebrated examples are the Korteveg-de Vries (KdV) and Boussinesq equations. We start with the reduction of our system to the Boussinesq-type equation. Such reduced models have been used in the analysis of pulse propagation in arteries, see [101, 102].

In order to keep the description self-contained, let us briefly introduce the integrable Boussinesq equation (BE),

[TABLE]

Note that the notations $u,\psi,\lambda$ in this section are not related to variables from the previous sections. We hope no confusion arises from the clash of notations, since these are the variables used commonly in the field of integrable equations.

Equation (4.5) arises in many various physical applications, most notably in describing propagation of the long waves in shallow water [103, 104]. The BE is integrable by the inverse scattering method, with the Lax pair representation that originates from Zakharov [105]:

[TABLE]

Here $\psi$ is the eigenfunction of the spectral problem and $\lambda$ is the spectral parameter. The solutions are explicitly described in terms of solitons which are computed by Hirota [106]. It is more convenient to rescale the variables in (4.5)

[TABLE]

where $\alpha$ and $\beta$ are some constants. The BE equation (4.5) then becomes

[TABLE]

Let us now turn our attention to the equation (3.35) for $\zeta=$const and no nonlinear terms in the potential, i.e $f(a)=0$. We get the following nonlinear equation, written here for reference:

[TABLE]

This equation will form the basis of derivation of approximate limiting cases in this chapter.

Let us consider a small amplitude and large wavelength propagation regime, where $B\rightarrow\varepsilon B,$ $\partial_{t}\rightarrow\varepsilon\partial_{t},$ $\partial_{x}\rightarrow\varepsilon\partial_{x}$ where $\varepsilon$ is a small parameter that defines the scale for the corresponding quantities. By keeping terms up to and including $\varepsilon^{2}$ in (4.9) we obtain

[TABLE]

The leading order terms in this equation are

[TABLE]

giving an operator equality to the leading order

[TABLE]

when the operators are acting on the solution $B$ of (4.10) . The plus sign is for the right - running waves, the minus sign is for the left running waves. Then, $c_{+}$ is, to the leading order, the propagation speed of the right-running waves and $c_{-}$ is the speed of the left-running waves.

Next, we use the expression for $\partial_{t}$ from (4.12) to make an approximation of the last term in (4.10)

[TABLE]

Then, this approximation of (4.10) gives an equation simultaneously valid for both left- and right-running waves:

[TABLE]

This equation provides a universal, Boussinesq-type approximation of equation (3.35). We present a numerical solution of this equation on the right panel of the Figure 3.2 to compare with the solution of full equation (3.35) with the same initial conditions. We observe that there is a very good quantitative agreement between these solutions for this propagation regime. Note that (4.13) is substantially algebraically simpler than (3.35). We do not yet know whether the equation (4.13) is integrable, but from the algebraic structure of this equation, we believe that it is much more amenable to analytical studies than the full (3.35).

While it is not known whether (4.13) is integrable, a further approximation can be made to reduce it to an integrable Boussinesq case. We can make the following approximation using (4.12):

[TABLE]

This approximation, substituted in (4.13) gives the following equation:

[TABLE]

Remembering that $a:=B_{s}$ is the relative deviation of the cross-secitonal area, we can write (4.15) as

[TABLE]

In order to write (4.16) in the canonical form (4.5), we make the following change of variables:

[TABLE]

In the new variables, (4.16) transforms to

[TABLE]

which is of the form

[TABLE]

The $+$ sign in front of $\beta^{2}$ term corresponds to the classical Boussinesq equation (4.8), and the $-$ sign corresponds to the alternative sign Boussinesq equation obtained from (4.8) by the substitution $x\rightarrow ix$ and $t\rightarrow it$, which is also integrable. The constants $\alpha$ and $\beta$ from (4.20) are given by

[TABLE]

Note that including $\alpha=0$ makes the equation (4.8) a linear equation. While this reduction is interesting, to describe both left- and right-running waves at the same time, we believe that equation (4.13) is more universal since it does not contain the $\pm$ signs.

4.3 Reduction to Korteweg-deVries (KdV) equation

It is also interesting to investigate the appropriateness of the use of Korteweg-deVries (KdV) approximation for this problem. While the universal Boussinesq-like equation (4.13) is valid for both left and right-travelling waves, its integrable counterpart (4.19) is written separately for the left- and right-running waves. However, each version of equation (4.19) itself has both the left- and the right-running waves. Thus, considering both left- and right version of (4.19) generates too much spurious behavior compared to the full system. In order to obtain an adequate integrable model for those waves, we reduce (3.35) to the KdV equation for waves running in each direction as follows. The KdV approximation for the fluid flow with elastic walls has been considered in the literature before, especially in the context of solitary wave propagation and variable properties of the arteries [107, 108] and reflection of the waves from the branching points [109]. Thus, we believe it is important to consider the appropriateness of KdV approximation in our case as well.

We are going to proceed as follows. Suppose that $B$ satisfies a KdV equation of the form

[TABLE]

Here $p,q,c$ are yet undetermined constants. A subindex $\pm$ will denote right/left running wave, and we shall enforce $c=c_{\pm}=U\pm\sqrt{\zeta}$ as defined by (4.12). We will drop the $\pm$ sign in the derivation below for brevity. From (4.22) we necessarily have $B_{t}=-cB_{s}+\mathcal{O}(\varepsilon^{2})$. Differentiating (4.22) with respect to $t$ gives

[TABLE]

Next, differentiating (4.22) with respect to $s$ gives

[TABLE]

From (4.23), (4.24) we arrive to

[TABLE]

Finally, adding equations in (4.25) we obtain

[TABLE]

which matches (4.15). Indeed, we take $c^{2}-2Uc=\zeta-U^{2}$, or $c=c_{\pm}=U\pm\sqrt{\zeta}$, as expected for the speeds of the left- and right-running waves (4.12). For the parameters $p$ and $q$, we obtain

[TABLE]

The equation for $B$ is then

[TABLE]

where we have implicitly assumed that $c$ has the label of $c_{\pm}$ for the left ($-$) and right ($+$) running waves. The KdV equation for $a=B_{s}$ is

[TABLE]

In order to compare the accuracy of the equation (4.28), we perform the simulation with the same initial conditions and parameters as in Figure 3.2, separating the solutions into the left- and right-running waves, i.e., taking the plus and minus signs in the corresponding equation (4.28). The results for $a=B_{s}$ are presented on Figure 4.1. Note that the speed of propagation is computed quite well as compared to the left panel of Figure 3.2, i.e., solution of the full equation (3.35). However, the details of the evolution do not quite match, especially because of the contribution of short wavelength ripples in the KdV equation, which is supposed to be the long-wavelength approximation only. Still, we believe that it is useful to consider the KdV for the propagation of the wave pattern in either direction.

5 Conclusions

In this paper, we have derived several novel results for the description of fluid with vorticity flowing in tubes with expandable walls. As our approach is variational in nature, it is applicable to the flow of essentially inviscid fluids. The friction terms can be introduced later in the model by using the Lagrange-d’Alembert’s principle of external forces. We have intentionally avoided discussion of the friction terms here, as the nature of fluid friction in blood is very complex and is certainly beyond the scope of this paper. That being said, in order to apply our theory successfully to the blood flow, the friction between the fluid and the walls will need to be introduced. This is an important feature of the practical applicability of the theory developed here and should be considered later.

On a more theoretical side, in this manuscript we have concentrated on the purely Lagrangian description of the dynamics. An equivalent Hamiltonian formulation can be derived by constructing an appropriate Poisson bracket. Since the description through the back-to-labels map comes as the Euler-Lagrange equation, the Poisson bracket for the back-to-label description will be canonical. We have not touched upon the Hamiltonian description here to keep the paper concise; the Hamiltonian description will be considered in the follow-up paper.

Another interesting direction of study, having both a theoretical and a practical value, is a more detailed incorporation of the dynamics of the media outside the tube. In our paper, we have modelled this media as simply exerting a constant pressure, but having no dynamics of its own. This may be a very good approximation for tubes conveying fluid when the outside media is an ideal gas kept at constant pressure. However, in reality, the dynamics of tubes submerged in a resisting media will cause a non-trivial response from the media beyond constant pressure, and thus affect the dynamics in a different ways. This response could be treated in the variational framework if, for example, the outside media is modelled either as an elastic body, or an ideal fluid. This interesting and challenging direction of study will be also considered in the future work.

6 Acknowledgements

We acknowledge inspiration from fruitful and productive discussions with D. D. Holm, F. Gay-Balmaz and T. S. Ratiu. The research of VP was partially supported by the NSERC Discovery Grant program and the University of Alberta. This work has also been made possible by the awarding of a James M Flaherty Visiting Professorship from the Ireland Canada University Foundation, with the assistance of the Government of Canada/avec l’appui du gouvernement du Canada.