1 Introduction

The onset of self-excited oscillations in the cardiovascular and respiratory system of human beings is a common phenomenon which can happen under normal or pathological conditions [1]. The reason can be found in the mutual interaction between relatively flexible biological conduits and the physiological flows. This kind of oscillations is crucial in the generation mechanism of adventitious sounds in the lungs, like wheezing and crackling [2] and is relevant to the diagnostics of lung diseases like asthma or chronic obstructive pulmonary disease [3]. Other interesting occurrences of this phenomenon in the human body can be found in the oscillations of the retinal vein [4], in the voice production process [5], and in the arterial tree [6]. As a consequence, the current interest of the scientific community in understanding the physical mechanism of self-excited oscillation in the flexible conduits of the human body is wide and multidisciplinary, spanning over many different fields, from engineering [7] to medicine [8].

One simple yet effective model of an airway or blood vessel undergoing self-excited oscillations is a collapsible tube [9, 10]. Despite the simple geometry, collapsible tubes can reproduce the rich phenomenology of the biological conduits such as the buckling due to large intramural pressure [11] and the onset of self-excited oscillations [12]. The key parameter which drives the collapse of a flexible conduit is the intramural pressure [13] i.e. the difference between the internal and external pressure of the domain. If this difference become negative, the tube will start to collapse. From a more quantitative perspective, in absence of flow, the description of the collapse is given by the tube law [14, 15] i.e. that relates the intramural pressure with the area of the smallest cross section of the tube (one example is given in Fig. 1). As shown in Fig. 1, the collapse of a flexible tube is characterised by three different phases [16]. For low values of the intramural pressure, the collapse happens symmetrically and the area of the constriction changes slowly with respect to the pressure (pre-buckling phase). For higher (absolute) values of the intramural pressure, the cross section of the tube is characterised by a n-buckling configuration [17, 18], where n is the number of lobes in the cross section (an example with \(n=2\) is shown in Fig. 1). In this regime (post-buckling phase), the area of the central cross section is more sensitive with respect to variations of the intramural pressure. Finally, when the walls of the internal part of the domain (the lumen) touch each other, the dependence of the central area on the intramural pressure changes again (post-contact phase). The transitions between these phases happen at specific values of the intramural pressure called buckling critical pressure and contact critical pressure, respectively [2]. Due to its clinical relevance (for instance for the assessment of obstructive sleep apnoea syndrome [19, 20]), the dependence of the buckling critical pressure on the geometric and elastic parameters of an elastic tube has been extensively investigated in the literature [21]. Interestingly, an analytic expression was already derived by von Mises in 1914 [22]. However, since its derivation is strongly dependent on shell theory assumptions, von Mises’ equation overestimates the value of the buckling critical pressure for tube geometries relevant for biomedical applications [23]. For this reason, the buckling critical pressure has been deeply investigated [24] and it is still an open question [25]. In some more recent attempts, the value of the buckling critical pressure for tubes geometries relevant for clinical applications has been studied with respect to the thickness of the tube [26] and its length [27] by employing experiments and numerical simulations. Following the results obtained by Turzi [28] on circular rings, the authors of this work have employed Landau’s theory of phase-transitions to estimate the value of the buckling critical pressure in 3D domains and its dependence on the geometric parameters of the system [15].

Fig. 1
figure 1

An example of tube law. The blue circles highlight the transition to the post-buckling and post-contact configurations. THe snapshots of the cross sections show the state of the constriction in the different phases of the collapse. This work focuses on the analysis of post-contact phase. The area is normalised with respect to the area of the constriction at rest (colour figure online)

Full size image

The full post-buckling behaviour of a collapsible tube, including the contact, has a great interest for both clinicians and researchers interested in the modelling of human vessels [9]. Most of the literature focuses on the behaviour of the fluid flow in the post-buckling [29] and in the post-contact [30] phase. However, the value of the contact critical pressure and its dependence on the geometric parameters of the system has been investigated in a less intensive manner. In a seminal paper by Flaherty et al. [31] this analysis has been carried on for an elastic ring spanning the full post-buckling behaviour of the system until its complete closure. However, this work relies on strong geometrical assumptions (1D elastic ring) and the effects connected to the tube thickness, length, and axial pre-stretch are consequently neglected. The analysis of the contact critical pressure has nonetheless a great interest for the study of self-excited oscillations. Indeed, as shown in a recent experimental investigation by Gregory et al. [2], the onset of self-excited oscillations happens for values of the intramural pressure comparable with the contact critical pressure.

In this work, the dependence of the contact critical pressure of a collapsible tube on its geometric parameters is analysed. The method relies on the implementation of an experimentally validated 3D numerical model of a collapsible tube. Moreover, by employing the non-dimensional variables proposed in [2], a set of general equations for the contact critical pressure is derived. This work is limited to a solid-only analysis. However, the presented methodology can be employed without modifications in a fully-coupled fluid–structure interaction setting and will be object of future investigations.

2 Problem setup

In the next sections, an in-depth description of the problem setup in terms of its geometry and boundary conditions is given. The details of the numerical model implementation, its validation, and the computation of the tube law are discussed.

2.1 Geometry and boundary conditions

The system under investigation is a collapsible tube with a thick external wall and under the effect of an initial axial stretch. Axial pre-stretch is commonly observed in human vessels, with strain values up to 60% of the length at rest [32]. One advantage of the simple geometry of a collapsible tube is that it needs a relatively small number of parameters to be described. Given the value of the inner diameter of the tube D and its length at rest \(l_0\), the following three non-dimensional variables completely characterise the geometry of the system:

$$\begin{aligned} d = \frac{l_0}{D},\quad \gamma = \frac{h}{D},\quad l = \frac{L}{l_0}, \end{aligned}$$
(1)

where h indicates the thickness of the external wall of the tube, and L represents the length of the domain after the imposition of the axial pre-stretch. In the following, the triplet \((d,\gamma ,l)\) will be employed to indicate the specific conditions investigated within the geometrical parameters space of collapsible tubes. The range of values of these parameters has been chosen to focus the study on the human airways. In a seminal paper by Horsfield [33], the correspondence between a given airway in the lungs and a collapsible tube is described. Human airways show a large distribution of thickness values (\(\gamma \in [0.04,0.1]\)) and are relatively short (\(d\in [3,6]\)) [33]. Another important aspect is that they operates under the effect of a remarkably high axial pre-stretch (\(l\in [1.1, 1.8]\) [32]). Such large value is going to largely influence the effective stiffness of the system, and it is therefore interesting to analyse the sensitivity of the contact critical pressure with respect to the axial pre-stretch.

To ensure the well-posedness of the problem, the system is subject to the following boundary conditions (see also the schematics in Fig. 2). One short side of the tube is clamped. The other short side is initially axially stretched by a quantity which depends on the value of the parameter l and then clamped. The external walls of the system are subject to an inward isotropic pressure \(p_{\textrm{ext}}\) whose value increases linearly in time according to

$$\begin{aligned} p_{\textrm{ext}}(t) = \frac{p_{\textrm{max}}}{\tau }t, \end{aligned}$$
(2)

where \(p_{\textrm{max}}>0\) is the maximum value of the external pressure and \(t\in [0,\tau ]\). The values of \(p_{\textrm{max}}\) employed in the study vary between 4000 and 8000 Pa. These are the typical values present in the lungs under forced expiration, depending on the size of the patient [34], for instance in presence of asthma [35]. As discussed in Sect. 1, the collapse of the system is driven by the intramural pressure. In formulae, the intramural pressure is determined by the following relation

$$\begin{aligned} p_{\textrm{intr}} = p_{\textrm{int}} - p_{\textrm{ext}} \end{aligned}$$
(3)

where \(p_{\textrm{int}}\) is the pressure inside the lumen. Since in this study the fluid flow is neglected, the internal pressure vanishes and the intramural pressure is only determined by \(p_{\textrm{ext}}\) i.e. the external pressure.

Fig. 2
figure 2

Schematics of the boundary conditions imposed on the system. The black box represents the initial configuration of the tube, before the pre-stretching. The left side of the tube is clamped, whereas the right one is pre-stretched by l. The linearly growing isotropic pressure \(p_{\textrm{ext}}\) is represented by the red arrows. The colour represents the displacement in the axial direction (colour figure online)

Full size image

2.2 Numerical model and validation

Once the geometric parameters \((d,\gamma ,l)\) are specified, the numerical model of the corresponding collapsible tube is implemented in the software Star-CCM+ (Siemens, version 2210). The geometry is defined by means of the built-in CAD software. The cross section of the tube is designed as an ellipse whose major axis is r and minor axis is 0.99r. In this way, it is possible to control the direction of the buckling. The boundary conditions of the system are imposed as discussed in the previous section. A second-order unsteady solver is implemented to compute the displacement vector field of the collapsible tube at each time step. The time step is the same for all the simulations and is \(\Delta t=0.1\) s. To enhance the convergence of the solution, each time step is divided into 10 inner iterations. This number of inner iterations ensures an average residual error of the order of \(10^{-14}\) at the end of each time step. The domain is discretised by means of two radial layers, 64 angular layers and 50 longitudinal layers of hexahedral cells. A sensitivity study for the mesh has been performed in [15] on the same model, showing that such mesh setup ensures reasonable convergence. Moreover, previous validations of the flow solver have been carried out for constricted geometries mimicking flow through the vocal folds in terms of LES pressures as compared with pressure from experimental data on the system walls [36] or in the case of cardiovascular flows in the thoracic aorta [37].

The material law chosen to model the full post-buckling and post-contact behaviour of the system is a hyperelastic Neo-Hookean model. As shown below, such model is able to reproduce the experimental results. One advantage of Neo-Hookean elasticity is that its parameters are defined from the Poisson’s ratio and the Young’s modulus of the tube material and do not require any ad-hoc assumption [38]. The strain energy is defined as

$$\begin{aligned} U = \alpha (I_1-3)+\frac{\beta }{2}(D-1)^2, \end{aligned}$$
(4)

where D is the determinant of the deformation gradient, and \(I_1\) is the first invariant of the right Cauchy-Green deformation tensor.

$$\begin{aligned} \alpha =\frac{1}{4}\frac{E}{(\nu +1)},\quad \beta =\frac{E\nu }{(1+\nu )(1-2\nu )}, \end{aligned}$$
(5)

given E and \(\nu \) as the Young’s modulus and the Poisson’s ratio, respectively. In this study, the value of the elastic parameters are \(E=1\) MPa and \(\nu =0.49\), which are the same of the experimental setup used to validated the numerical model [39] (see below for further details). The deformation energy in Eq. (4) can be further improved by considering its regularisation with high order or micro-morphic models [40,41,42]. However, as discussed at the end of this subsection, the displacement computed by means of the model in Eq. (4) is able to replicate the experimental data. However, the comparison of the results in terms of different generalisation of elasticity theory [43,44,45] will be treated in a future work. Modelling the contact between different parts of the domain is an open challenge in computational mechanics [46, 47]. In this work, the contact between the internal walls of the tube is treated by employing a repulsive virtual plane (see Fig. 3). This virtual plane coincides with the major axis of the elliptic cross section of the tube (see Fig. 3). When the virtual plane and the internal walls of the tube are closer than a threshold distance (in this case \(5\times 10^{-6}\) m), the walls start to be affected by a repulsive force (\(|\vec {F}_{\textrm{rep}}|=10^{12}\) N) oriented with the same direction of the normal of the plane. By accurately tuning the threshold distance and the module of the repulsive force, it is possible to mimic the contact behaviour described mathematically in [31] and observed experimentally in [27]. These works described qualitatively the behaviour of a collapsible tube in its post-contact phase which is well reproduced by the numerical model. When the intramural pressure reaches the contact critical pressure value, a pair of points at the opposite sides of the cross section touches each other (see Fig. 3a). As the intramural pressure reaches larger negative values, the contact region becomes a straight-line segment (see Fig. 3b). For larger negative values of the intramural pressure, the length of the segment increases and tends to approximate the length of the perimeter of the cross section (see Fig. 3c).

Fig. 3
figure 3

The numerical model is able to reproduce the post-contact behaviour described in literature. The colour corresponds to the radial displacement. The black dashed line indicates the position of the repulsive plane employed to model the contact (colour figure online)

Full size image

The numerical model employed in this study is also able to reproduce quantitatively the full post-buckling and post-contact behaviour of a collapsible tube. The experimental validation of the model has been performed by using a set of publicly available experimental measurements [2, 39]. In the experimental routine, a collapsible tube with relevant values of the geometric parameters [15] has one side clamped and one axially stretched and then clamped. The intramural pressure is controlled by a syringe and monitored with a manometer. A 3D camera system is able to register the deformation of the collapsible tube and area of the central cross section of the tube is determined. The corresponding value of the intramural pressure is registered. This information is sufficient to reconstruct the experimental tube law of the collapsible tube from the pre-buckling configuration to the post-contact phase. The same geometry, initial and boundary conditions are implemented in the current numerical setup to replicate the experiments. The corresponding tube law is then compared with the experimental one (see next section for the details of the computation of the tube law). As shown in Fig. 4, the numerical results are able to reproduce the full experimental tube law.

Fig. 4
figure 4

Comparison [15] of the numerically computed tube laws with the corresponding ones measured experimentally [2, 39]. The error bars are computed as the absolute difference between the maximum and minimum value of the initial area, which has been measured 3 times

Full size image

2.3 Computing the tube law

The output of the simulation is the tube law of the system i.e. the relation between the intramural pressure and the corresponding central cross-sectional area. As discussed in the next section, the tube law will be employed to estimate the contact critical pressure. The central cross section (in the pre-stretched configuration) is of interest because is the most collapsed, being the furthest cross section from the constrained sides of the tube. At each time step, the value of the intramural pressure is registered as \(p^i_{\textrm{intr}} = -p^i_{\textrm{ext}}\), where the index \(i=1,\dots ,M\) represents the i-th time step. To compute the value of the corresponding area, the radial coordinates \((\theta ^i_j, r^i_j)\) of the central cross-sectional perimeter are registered at each time step, where the index j runs over the \(N=64\) mesh elements on the perimeter. It is then possible to compute the area \(A^i\) of the central cross section at the i-th time step by:

$$\begin{aligned} A^i = \frac{1}{2}\sum ^{N-1}_{j=1}\left| \theta ^i_{j+1}-\theta ^i_j\right| \left( r^i_j\right) ^2. \end{aligned}$$
(6)

The tube law is defined as the couple \(\left\{ p^i_{\textrm{intr}}\right\} _{i=1}^M\), \(\left\{ A^i_{\textrm{intr}}\right\} _{i=1}^M\).

3 Analysis of the contact critical pressure

At the transitions between the pre to the post-buckling phase and from the post-buckling to the post-contact phase (blue circles in Fig. 1), the tube law shows a clear change of regime connected to the different collapse mechanisms that the tube experiences when the intramural pressure reaches more negative values. In [15], the authors of this work have demonstrated that it is possible to treat the transition between the pre to post-buckling configurations as a second-order phase transition (see also the work by Turzi [28] on 1D rings). Moreover, the buckling critical pressure has been estimated by using Landau’s theory of phase transition. However, it does not appear possible to estimate the contact critical pressure by means of such methodology. The main reason is that there is no clear symmetry breaking in the system before and after the contact. Moreover, the nature of the contact is intrinsically discontinuous, i.e. the transition between the state of the system before and after the contact happens instantaneously and not continuously as in the case of a second-order phase transition. Consequently, the main assumptions of the methodology described in [15] are not valid for the treatment regarding the contact critical pressure. The discontinuous nature of this transition, however, can be exploited to address the problem. The main observation is that the first contact between the internal walls of the system corresponds to a sudden change in the slope of the tube law. It is possible to estimate the value of the contact critical pressure by looking for the corresponding discontinuity in the first derivative of the tube law computed with respect to the intramural pressure (see Fig. 5). It is important to take into account the effect of the repulsive virtual plane discussed in Sect. 2.2 (see Fig. 3) since it will interact with the internal walls of the tube before the actual contact happens. A careful analysis of this effect shows that the virtual plane induces two small additional discontinuities in the tube law (blue and green dots in Fig. 5). Consequently, the contact critical pressure is estimated as the pressure value at the end of the discontinuity in the derivative of the area with respect to the intramural pressure (i.e. the red dot in Fig. 5). The methodology to determine the contact critical pressure is the following:

  1. 1.

    For any triplet \((d,\gamma ,l)\) of geometric parameters, the corresponding numerical model is implemented according to Sect. 2.2.

  2. 2.

    The corresponding tube law, i.e. the two sets \(\{p_{\textrm{intr}}^i\}_{i=1}^M\), \(\{A^i\}_{i=1}^M\), is computed as discussed in Sect. 2.3.

  3. 3.

    The tube law is pre-processed. The pre-buckling region of the tube law is neglected. The derivative of the area of the central cross section with respect to the intramural pressure is numerically computed (an example is shown in Fig. 5) by using a Python algorithm employing the numpy.diff function [48].

  4. 4.

    The value of the contact critical pressure is estimated by identifying the discontinuity of the derivative of the tube law (the red dot Fig. 5). The contact critical pressure and the corresponding area are registered.

Fig. 5
figure 5

The main plot shows the derivative of the tube law with respect to the intramural pressure in a neighbourhood of the contact critical pressure. The smaller sub-figure shows the corresponding portion of the tube law. The blue and green dots highlight the spurious discontinuities due to the repulsive plane. The red dot indicates the value of the intramural pressure corresponding to the actual contact (colour figure online)

Full size image

As discussed in Sect. 4, the results obtained using this method are experimentally confirmed.

3.1 The role of the geometric parameters

By following the methodology described in the previous section, it is possible to study how the geometric parameters of the system affect the contact critical pressure. As first, the dependence of the contact critical pressure on the length-to-diameter ratio d is analysed. The values of the contact critical pressure \(p^{\textrm{cont}}_{\textrm{crit}}\) corresponding to \(d\in (3,3.5,4,4.5,5,5.5,6)\) are shown in Fig. 6a. The corresponding tube laws are plotted in Fig. 6b. The relation between \(p^{\textrm{cont}}_{\textrm{crit}}\) and d is well described by the following relation

$$\begin{aligned} -p^{\textrm{cont}}_{\textrm{crit}} = Ad^B+C \end{aligned}$$
(7)

where A, B, C are free parameters whose values are listed in Table 1. This model shows that for longer tubes the absolute value of the contact critical pressure decrease. This is compatible with the behaviour of the buckling critical pressure analysed in [15], whose estimated values are indicated as black boxes in Fig. 6b and in the following figures.

Fig. 6
figure 6

The left panel shows the values of the contact critical pressure estimated according to Sect. 3 and the fit model in Eq. (7). In the right panel, the tube laws associated with the different values of d are displayed. The black boxes are the values of the buckling critical pressures estimated in [15], while the black circles are the contact critical pressures

Full size image

Let now analyse the effect of the thickness-to-diameter ratio \(\gamma \) on the contact critical pressure. The values of the contact critical pressure \(p_{\textrm{crit}}^{\textrm{buckl}}\) corresponding to \(\gamma \in (0.05, 0.06, 0.07, 0.08, 0.09)\) are shown as circles in Fig. 7a. The corresponding tube laws are plotted in Fig. 7b. The relation between \(p_{\textrm{crit}}^{\textrm{buckl}}\) and \(\gamma \) is well described by the following function

$$\begin{aligned} -p_{\textrm{crit}}^{\textrm{buckl}} = A\gamma ^B + C \end{aligned}$$
(8)

where A, B, C are free parameters whose values are shown in Table 1.

Fig. 7
figure 7

The left panel shows the values of the contact critical pressure estimated according to Sect. 3 and the fit model in Eq. (8). In the right panel, the tube laws associated with the different values of \(\gamma \) are displayed. The black boxes are the values of the buckling critical pressures estimated in [15], while the black circles are the contact critical pressures

Full size image

Finally, the relation between the contact critical pressure \(p_{\textrm{crit}}^{\textrm{buckl}}\) and the pre-stretch ratio l is considered. The values of the contact critical pressures for \(l\in (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8)\) are shown in Fig. 8a. The corresponding tube laws are shown in Fig. 8b. The relation between the contact critical pressure and the pre-stretch ratio has been obtained via fitting procedure with the following model:

$$\begin{aligned} -p_{\textrm{crit}}^{\textrm{buckl}} = A\tanh (Bl) + C \end{aligned}$$
(9)

where A, B, C are free parameters whose values are listed in Table 1. Interestingly, the contact critical pressure tends asymptotically to a constant value for higher value of the pre-stretch ratio, as depicted in Fig 8a.

Fig. 8
figure 8

The left panel shows the values of the contact critical pressure estimated according to Sect. 3 and the fit model in Eq. (9). In the right panel, the tube laws associated with the different values of l are displayed. The black boxes are the values of the buckling critical pressures estimated in [15], while the black circles are the contact critical pressures

Full size image
Table 1 Values of the parameters obtained from the fit of the contact critical pressures with Eqs. (7–9)
Full size table

4 Non-dimensional equations

It is possible to employ the results of the analysis presented in the previous section to derive a system of general equations for the contact critical pressure and the corresponding cross-sectional area. These equations depend on the elastic and geometric parameters of any collapsible tube in a reasonable neighbourhood of the validated model. In [2], by using experimental data relevant for biomedical flows, Gregory et al. have demonstrated that by re-scaling the intramural pressure and the central cross-sectional area as shown in Eq. 10, tube laws corresponding to different geometric and elastic parameters would collapse in a single curve. The non-dimensional variables are defined as:

$$\begin{aligned} \hat{p}_{\textrm{intr}}=\frac{p_{\textrm{intr}}}{E/(1-\nu ^2)}l^{-2}\gamma ^{-1}d,\quad \hat{A}=\frac{A}{\pi r^2}l. \end{aligned}$$
(10)

By applying the definitions in Eq. 10 to the tube laws shown in Figs. 6b, 7 and 8b, the same behaviour is observed, as computational data tend to collapse, with a relatively small dispersion, around a single line (see Fig. 9).

Fig. 9
figure 9

(Left panel) All the tube laws analysed in this work. (Right panel) By scaling the intramural pressure and the cross-sectional area by means of Eq. (10), the tube laws tend to collapse on a single line [15]

Full size image

Let consider the set of all the contact critical pressures and the corresponding cross-sectional areas obtained by means of the method described in Sect. 3 spanning the geometric parameters \((d,\gamma ,l)\). By re-defining these values in terms of Eq. 10, the following two sets are derived

$$\begin{aligned} \left\{ \hat{p}_{\textrm{cont}}^k\right\} _{k=1}^R,\quad \left\{ \hat{A}_{\textrm{cont}}^k\right\} _{k=1}^R, \end{aligned}$$
(11)

which are defined as the sets of the all non-dimensional contact critical pressures and contact critical areas derived in this work. The index k runs on all the different (\(R=20\)) combinations of geometric parameters analysed in this study. By computing the corresponding average and the standard deviation and by combining the results with Eq. 10, the following general equations are derived:

$$\begin{aligned} \frac{p^{\textrm{crit}}_{\textrm{cont}}}{E/(1-\nu ^2)}d\left( l\right) ^{-2}\left( \gamma \right) ^{-1}=-0.104\pm 0.013;\qquad \frac{A^{\textrm{crit}}_{\textrm{cont}}}{\pi r^2}l = 0.309\pm 0.014. \end{aligned}$$
(12)

These values are compatible within the error with the ones estimated by Gregory [2] for the onset of self-exciting oscillations, here reported for the reader’s convenience:

$$\begin{aligned} \frac{\tilde{p}^{\textrm{crit}}_{\textrm{cont}}}{E/(1-\nu ^2)}d\left( l\right) ^{-2}\left( \gamma \right) ^{-1}=-0.12\pm 0.02;\qquad \frac{\tilde{A}^{\textrm{crit}}_{\textrm{cont}}}{\pi r^2}l = 0.27\pm 0.07. \end{aligned}$$
(13)

This observation experimentally confirms the methodology adopted in this work for the treatment of the contact of the internal walls of a collapsible tube and endorses the observation that the self-excited oscillations in collapsible tubes are triggered for values of the intramural pressure compatible with the contact critical pressure.

5 Conclusions

The dependence of the contact critical pressure of a collapsible tube on the geometric parameters of the system has been analysed in this work. The methodology is based on the post-processing of numerical data obtained via the implementation of the numerical model of a digital replica of the collapsible tube under given boundary conditions. The numerical model has been experimentally validated. A set of general non-dimensional equations for the contact critical pressure and the corresponding central cross-sectional area has been derived. The results are compatible with the experimental observation by Gregory et al. [2] that the value of the intramural pressure for the onset of self-excited oscillations in collapsible tubes is close to the contact critical pressure. The main advantage of the presented methodology is that it can be generalised to a two-way fluid–structure interaction analysis without further assumptions, which will be the aim of a following work. The main limitation of the presented modelling scheme is the absence of friction between the internal walls of the tube. Currently, the possibility to integrate the friction by means of a strategy similar to the one discussed in [49] in the present modelling scheme is under investigation. Another interesting perspective is to include in the modelling scheme the effects of damage [50] due to the repeated loading and large deformation.