Pulsatile flow in a thin-walled viscoelastic tube

Abstract

Low-inertia pulsatile flows in highly distensible viscoelastic vessels exist in many biological and engineering systems. However, many existing works focus on inertial pulsatile flows in vessels with small deformations. As such, here we study the dynamics of a viscoelastic tube at large deformation conveying low-Reynolds-number oscillatory flow using a fully coupled fluid–structure interaction computational model. We focus on a detailed study of the effect of wall (solid) viscosity and oscillation frequency on tube deformation, flow rate, phase shift and hysteresis, as well as the underlying flow physics. We find that the general behaviour is dominated by an elastic flow surge during inflation and a squeezing effect during deflation. When increasing the oscillation frequency, the maximum inlet flow rate increases and tube distention decreases, whereas increasing solid viscosity causes both to decrease. As the oscillation frequency approaches either $0$ (quasi-steady inflation cycle) or $\infty$ (steady flow), the behaviours of tubes with different solid viscosities converge. Our results suggest that deformation and flow rate are most affected in the intermediate range of solid viscosity and oscillation frequency. Phase shifts of deformation and flow rate with respect to the imposed pressure are analysed. We predict that the phase shifts vary throughout the oscillation; while the deformation always lags the imposed pressure, the flow rate may either lead or lag depending on the parameter values. As such, the flow rate shows hysteresis behaviour that traces either a clockwise or counterclockwise curve, or a mix of both, in the pressure–flow rate space. This directional change in hysteresis is fully characterised here in the appropriate parameter space. Furthermore, the hysteresis direction is shown to be predicted by the signs of the flow rate phase shifts at the crest and trough of the oscillation. A distinct change in the tube dynamics is also observed at high solid viscosity which leads to global or ‘whole-tube’ motion that is absent in purely elastic tubes.

1. Introduction

In both biological and engineering applications, many systems appear in which a distensible tube conveys pulsatile flow (Berger & Jou 2000; Grotberg & Jensen 2004). The most common example is the pulsatile flow in major arteries caused by the cardiac rhythms. Early theoretical works in this area broadly follow one of two approaches, the first being the use of a lumped-parameter model to predict the relationship between pressure and flow rate. The earliest such model was the Windkessel model, formulated by Frank (1899), consisting of a resistor and compliance element in parallel, representing the characteristic fluid resistance and arterial wall compliance, respectively. Several improvements over this two-element model have been made, with the three- and four-element models introducing the effects of wave propagation and inertia (Westerhof, Lankhaar & Westerhof 2009). Similar lumped parameter models may be derived by the addition or rearrangement of the elements, analogous to an electrical circuit (Zamir 2016a ; Soni et al. 2022). Such models are useful in describing the averaged pressure–flow rate relationship, but they are lacking in a detailed flow description and require constants to be fitted to experimental results.

The second common approach was to construct a quasi-one-dimensional (1-D) wave propagation model assuming small deformation, relating the fluid velocity to the pulse wave speed and pressure waveform. The large diameters and high fluid velocities in large arteries enabled early researchers to neglect fluid viscosity (Pedley 1980; Fung 1997; Fung 1997). The derivation of the wave speed in an elastic tube filled with inviscid fluid is often attributed to Moens (1878) and Korteweg (1878), though it can be traced back to Young (1808) (Fung 1997). Further works have also considered inviscid models (Parker & Jones 1990), while others include viscous friction based on an assumed velocity profile (San & Staples 2012; Roknujjaman et al. 2023). Roknujjaman et al. (2023) focused on wave propagation behaviour, finding that increased tube compliance leads to increased attenuation. The same behaviour was observed by Fancher & Katifori (2022), who applied a 1-D model to vessel networks to examine wave reflections. They also note that the amplitude of reflected waves increases with vessel stiffness, leading to large spikes in both pressure and flow rate. Other works have examined the effects of arterial taper (Anliker, Rockwell & Ogden 1971), which was found by Abdullateef, Mariscal-Harana & Khir (2021) to alter the pressure profile due to induced wave reflections.

A description of the flow field can be obtained by solving the governing equations directly. For this, the linearised Navier–Stokes equations coupled with the motion of the tube wall have been used, also allowing the consideration of fluid viscosity. The most well-known model for sinusoidal flow of viscous fluid in an elastic tube was given by Womersley (1955, 1957a ) and Morgan & Kiely (1954). Under the assumption of weakly deforming tubes, their solution provides a full axisymmetric description of the flow field and tube motion in response to an oscillatory pressure gradient, along with the wave speed. One interesting result of their work is that the predicted velocity leads the driving pressure in phase, even though they both oscillate with the same frequency. Furthermore, Womersley (1957b ) notes that this lead is highest ( $=\pi /4$ ) in the case of a rigid tube and reduces with increasing compliance.

Other theoretical models have been developed that are either built on those of Womersley (1955, 1957a ) and Morgan & Kiely (1954) (see Cox 1969) or use different approaches. Notable is the numerical work by Ling & Atabek (1972), who considered large deformation and nonlinear flow effects. Key assumptions of their model are that the pressure is constant in the radial direction and the pressure–radius relationship is known. These allow the wall and fluid motion to be decoupled. Further works have considered tubes of altered geometries, such as a smoothly expanding conduit considered by Pedrizzetti et al. (2002) using a perturbation approach with the zeroth-order solution being that of a rigid tube. Myers & Capper (2001) considered an elastic curved tube using a perturbation approach with the zeroth-order solution reducing to Womersley’s result.

This problem has also been explored numerically by fully coupled three-dimensional (3-D) fluid–structure interaction (FSI) simulations, many of which are inspired by arteries (Wang, Wood & Xu 2015; Kim, Park & Lim 2016; Mu et al. 2019). Chaniotis et al. (2010) studied pulsatile flows in various geometries, including curved tubes and bifurcations, seen in coronary arteries. More complex arterial geometries have also been considered, including a ventricle (Esmaily Moghadam et al. 2013) and the coronary arterial tree (Eslami et al. 2019). In particular, Eslami et al. (2019) found that the instantaneous wall shear stress differed between rigid and elastic geometries due to the propagation of flow and pressure waves. Further work by Bäumler et al. (2020) simulated pulsatile flow in an aortic dissection, exploring the effect of the elasticity of the dissection flap. They demonstrated that decreasing the stiffness greatly impacts the flow dynamics, showing up to a $20\,\%$ decrease in flow rate in the true lumen.

Physiological vessel walls are normally viscoelastic in nature which affects the transient dynamics. A key effect of the solid viscosity is a hysteresis behaviour of the imposed pressure and tube diameter (and, pressure and flow rate) which arises because of the different tube response during inflation and relaxation phases. The influence of solid viscosity was considered by several authors. Womersley (1957a ) refers to a modification made by Morgan & Kiely (1954) in their original work, replacing the elastic moduli and Poisson ratio with complex variants to represent the internal damping of the tube. It is shown that the solid viscosity tends to increase the damping of the travelling wave and decrease its velocity. A quasi-1-D model from Mal, Soni & Nayak (2024) also predicts pronounced damping of the pulse wave when the tube is treated as a modified Zener viscoelastic material. Similarly, the quasi-1-D model of Pontrelli (2002) predicts the solid viscosity attenuates high-frequency oscillations. Dragon & Grotberg (1991) followed a perturbation approach to model the same problem, although the focus of their work was not on the effect of viscoelasticity. They found that in a volume-cycled flow, the average flow rate is maximal at an intermediate value of oscillation frequency. A two-dimensional (2-D) model was formulated by Čanić et al. (2006a ,b) using perturbation and homogenisation theory, treating the tube as a weakly deforming Kelvin–Voigt viscoelastic material. Full FSI solutions for viscoelastic tubes from Wang et al. (2015) and Kim et al. (2016) showed the characteristic pressure–diameter hysteresis behaviour, with Wang et al. (2015) further noting a reduction in deformation when compared with the purely elastic case.

While the above studies were motivated by flow in large arteries, where inertia is dominant and deformation is not large, there are many examples of pulsatile or oscillatory flows in small vessels where the fluid viscous effects are dominant and deformation is large. Examples include microcirculatory vessels, initial lymphatics and alveolar ducts (Fung 1997; Denny & Schroter 1999; Fathy El-Amin 2016; Mallik, Mukherjee & Panchagnula 2020). Many recent studies have shown that blood vessels in the microcirculation exhibit pulsatile flow following cardiac rhythms (Zweifach 1974; Lee et al. 1994; Klassen et al. 1997; Kajiya et al. 2008; Rashid et al. 2012; Shih et al. 2015; Aby, Guevara-Torres & Schallek 2019; Bedggood & Metha 2019). These vessels are known to undergo passive dilation and relaxation in response to flow pulsation and active vasomotion for blood flow regulation. As discussed by Schmid-Schönbein et al. (1988; 1989), pulsation causes wave motion and a delayed flow response. However, these behaviours are not governed by fluid inertia, as is the case in large arteries, but by fluid viscosity and the microvessel’s distensibility. Vessel diameter changes are also significant; Shih et al. (2015) observed $40\,\%$ spontaneous change in mice cerebral microvessel diameter, while Yada et al. (1993) measured $20\,\%{-}40\,\%$ change in porcine cardiac microvessels. Lee & Schmid-Schönbein (1995) measured nearly $100\,\%$ change in skeletal capillary diameter under a transmural pressure change of $0$ – $100$ cm water. Toyota et al. (2005) and Kajiya et al. (2008) reported nearly $100\,\%$ change in cardiac microvessel diameter in a beating canine heart. For capillary vessels in the human retina, Neriyanuri et al. (2023) and Gu et al. (2018) reported the degree of flow pulsatility, defined as the fractional change of velocity over one cardiac cycle with respect to the mean velocity, whose average value is approximately $0.8$ . For lymphatic vessels, multiple pumping mechanisms exist to drive flow, creating a natural pulsatility (Moore Jr. & Bertram 2018). The smallest vessels in the lung (bronchioles and alveolar ducts) expand and contract in response to breathing (Denny & Schroter 1999; Mallik et al. 2020). In microfluidics, pulsatile flow is often used in deformable channels resulting in complex, transient flow dynamics (Iyer et al. 2015; Raj, Suthanthiraraj & Sen 2018; Dincau, Dressaire & Sauret 2020; Xia et al. 2021).

Microvessels also exhibit viscoelastic behaviour, which affects the transient dynamics of the coupled fluid–structure system (Fung 1981; Lee & Schmid-Schönbein 1990). Lee & Schmid-Schönbein (1990) explored these behaviours assuming a quasilinear standard solid model for the vessel wall material with a sinusoidally oscillating applied pressure. They assumed the fluid motion to be unidirectional and the flow rate to be constant for the whole tube and related to the local pressure gradient by Poiseuille’s formula. They found that the inlet flow rate exhibits hysteresis with respect to the applied pressure whose behaviour varied with oscillation frequency, primarily proceeding either clockwise or counterclockwise. The hysteresis direction was attributed to the balance between the solid viscosity and fluid viscosity, with the former biasing the hysteresis counterclockwise and the latter biasing it clockwise. Notably, at low frequencies when the time delays caused by both viscosities were of similar magnitude, the hysteresis curves crossed over themselves, forming a twisted figure-of-eight.

The above review suggests that, for pulsatile flows in viscoelastic microvessels undergoing large deformation, a fully coupled FSI model without any simplification of the governing equations of the fluid motion (e.g. unidirectionality, linearisation etc.) is lacking. As such, it is unknown how the consideration of large-deformation complete flow equations and a full fluid–structure coupling would alter the previous findings described above. Would the hysteresis transition from clockwise to counterclockwise as predicted by Lee & Schmid-Schönbein (1990) still exist? How would this behaviour change in the parameter space defined by solid viscosity and imposed oscillation frequency? How do tube deformation and flow rate evolve in the parameter space? Also lacking is a systematic study of the tube dynamics over the parameter space and its physical understanding. To fill this knowledge gap, in this work, we present a detailed numerical study of the viscoelastic response of a compliant tube subject to a pulsatile flow undergoing large deformation using a full FSI model and complete flow equations. We primarily focus on the role of solid viscosity and pulsation frequency, and limit to small but non-vanishing Reynolds number. The value of this study lies in the exploration of tubes with large deformations over a wide parameter space, which is lacking in the current literature. Our results suggest that the general tube and flow behaviour is dominated by elastic effects, and that deformation and flow rate are most affected in the intermediate range of solid viscosity and oscillation frequency. We further find that the phase shifts of both the deformation and flow rate with respect to the imposed pressure vary throughout the oscillation. Additionally, at high solid viscosity, there is a distinct change in the tube dynamics.

2. Methodology

Figure 1. (a) Problem set-up. (b) Standard linear solid viscoelastic model (Maxwell form). (c) Tube mesh close up. Parameters: $\beta =0.05$ , $De=50$ , $\alpha =0.11$ .

The problem set-up is given in figure 1(a). The undeformed tube is assumed to have a circular cross-section with radius and total length denoted by $R_0$ and $L$ . The vessel and fluid motion have two-way (i.e. strong) coupling; the vessel is deformed naturally (i.e. without any prescribed displacement) by the fluid motion and the flow, in turn, is affected by the movement of the vessel wall. The streamwise flow direction is $x$ , and $x=0$ and $L$ are the inlet and outlet, respectively. The radius of the deformed vessel varies both axially and in time as $R(x,t)$ . The vessel is immersed within a rectangular computational domain with the inlet and outlet coinciding with the left and right boundaries of the domain. The flow inside the vessel is driven by specifying an oscillatory pressure $P_{in}$ with frequency $f$ at the inlet, while the outlet pressure $P_{out}$ is held constant. An external pressure $P_{ext}$ is specified over the left and right boundaries of the computational domain outside the tube inlet and outlet. The fluids interior and exterior to the tube are assumed to have same density $\rho$ and viscosity $\mu$ . The tube is divided into three streamwise segments: two rigid segments at the entrance and exit, and a deforming part of length $L_{def}$ in between. The two ends of the deforming segment are thus ‘pinned’, i.e. they are constrained to the ends of the rigid segments and do not resist bending moments. The current model is capable of modelling the flow continuously through the transitions between the rigid and deforming segments. Therefore, no fluid boundary conditions are required at these locations.

2.1. Structural mechanics

The tube wall in the deforming section is materially isotropic and modelled as a standard linear viscoelastic solid with vanishing thickness. Then, the structural mechanics can be recast as a 2-D plane stress problem in the tangent plane of the tube surface. The undeformed and deformed states of the vessel surface are represented by coordinates of a material point as $\boldsymbol {X}$ and $\boldsymbol {x}(\boldsymbol {X},t)$ , respectively. The surface deformation gradient and Green strain tensors are then (Barthès-Biesel & Rallison 1981; Barthès-Biesel et al. 2002)

(2.1) \begin{equation} \hat {\mathbf {A}} = \left ( \mathbf {I} - \boldsymbol {nn} \right ) \cdot \frac {\partial \boldsymbol {x}} {\partial \boldsymbol {X}} \cdot \left ( \mathbf {I} - \boldsymbol {NN} \right ) \end{equation}

and

(2.2) \begin{equation} \hat {\mathbf {E}}=\frac {1}{2}({\hat {\mathbf {A}}}^{\mathbf {T}}\cdot \hat {\mathbf {A}}-(\mathbf {I}-\boldsymbol {NN})), \end{equation}

where $\boldsymbol {N}$ and $\boldsymbol {n}$ are the reference and deformed surface normal vectors, and $\mathbf {I}$ is the identity tensor. The 2-D Cauchy stress (i.e. surface traction tensor) is expressed as

(2.3) \begin{equation} \tilde {\boldsymbol {\unicode {x03C3}}}=\frac {1}{J_s}\hat {\mathbf {A}}\cdot \frac {\partial W_s}{\partial \hat {\mathbf {E}}}\cdot {\hat {\mathbf {A}}}^T, \end{equation}

where $J_s=\lambda _1\lambda _2$ is the surface area dilation and $W_s$ is the surface strain energy defined per unit area in the reference configuration (Pozrikidis 2003; Barthès-Biesel & Rallison 1981; Barthès-Biesel et al. 2002). The left Cauchy–Green surface deformation tensor, ${\hat {\mathbf {A}}}^T\cdot \hat {\mathbf {A}}$ , has two non-zero eigenvalues, $\lambda _1^2$ and $\lambda _2^2$ , where $\lambda _1$ and $\lambda _2$ are the stretch ratios along the principal axes on the surface. The corresponding principal strain components are $\frac {\lambda _i^2-1}{2}$ , $i=1,2$ . The principal traction components $\tilde {\sigma }_1$ and $\tilde {\sigma }_2$ can further be written as

(2.4) \begin{align} \tilde {\sigma }_1=\frac {1}{\lambda _2}\frac {\partial W_s}{\partial \lambda _1}, \qquad \tilde {\sigma }_2=\frac {1}{\lambda _1}\frac {\partial W_s}{\partial \lambda _2}. \end{align}

Due to material isotropy, the traction tensor can be expressed using the above components as

(2.5) \begin{equation} \tilde {\boldsymbol {\unicode {x03C3}}} = \tilde {\sigma }_1\,{\boldsymbol {e}}_1\otimes \boldsymbol {e}_1 +\tilde {\sigma }_2\,{\boldsymbol {e}}_2\otimes \boldsymbol {e}_2, \end{equation}

where $\boldsymbol {e}_1$ and $\boldsymbol {e}_2$ are the unit eigenvectors of ${\hat {\mathbf {A}}}^T\cdot \hat {\mathbf {A}}$ (Barthès-Biesel & Rallison 1981; Pozrikidis 2003).

The elastic response of the wall material is assumed to follow Hooke’s law, although the formulation and numerical methodology are applicable to hyperelastic models as well. The principal elastic tension $\tilde {\sigma }_1$ is then

(2.6) \begin{equation} \tilde {\sigma }_1=\frac {G_s}{1-\nu _s}[\lambda _1^2-1+\nu _s(\lambda _2^2-1)], \end{equation}

where $G_s$ is the shear modulus and $\nu _s$ is the Poisson ratio (Barthès-Biesel et al. 2002). The other component can be found by interchanging the indices.

The formulation for the viscoelastic stress is now briefly presented. The detailed derivation can be found elsewhere (Simo 1987; ABAQUS 2002; Yazdani & Bagchi 2013; Matteoli, Nicoud & Mendez 2021). The basic equation governing a linear isotropic viscoelastic material under small strain is given in terms of the Cauchy stress on the surface expressed by a hereditary integral

(2.7) \begin{equation} \tilde {\boldsymbol {\unicode {x03C3}}}(t) = \int _0^t 2G(t-s)\frac {{\rm d}\mathbf {e}}{{\rm d}s}{\rm d}s + \mathbf {I} \int _0^t K(t-s)\frac {{\rm d}\chi }{{\rm d}s}{\rm d}s, \end{equation}

where $\mathbf {e}$ and $\chi$ are the deviatoric and volumetric strains, and $G$ and $K$ are the small-strain shear and bulk relaxation moduli. For the most general linear viscoelastic model (the generalised Maxwell model), these relaxation moduli may be represented as a series of decaying exponentials (a Prony series) as (Fung 1981; Dill 2006)

(2.8) \begin{equation} G(t)=G_\infty +\sum _{i=1}^N G_i e^{-t/t_i^G}, \qquad K(t)=K_\infty +\sum _{i=1}^N K_i e^{-t/t_i^K}, \end{equation}

where $G_\infty$ and $K_\infty$ are the long-term relaxation moduli, $N$ is the number of terms in the series, and $t_i^G$ and $t_i^K$ are the relaxation times associated with each term. The instantaneous shear modulus is introduced as $G_0=G_\infty +\sum _{i=1}^N G_i$ . Neglecting the viscous effect of the volumetric strain, we also have $K(t)=K_\infty =K_0$ . Then, (2.7) becomes

(2.9) \begin{equation} \tilde {\boldsymbol {\unicode {x03C3}}}(t) = \tilde {\boldsymbol {\unicode {x03C3}}}_0^d(t) + \int _0^t 2\dot {G}(t-s)\mathbf {e}(t-s)\textrm{d}s + \tilde {\sigma }^v(t)\mathbf {I}, \end{equation}

where $\tilde {\boldsymbol {\unicode {x03C3}}}_0^d(t) = 2G_0 \mathbf {e}(t)$ is the instantaneous deviatoric stress, $\tilde {\sigma }^v$ is the volumetric stress and $\dot {G}={\rm d}G/{\rm d}t$ . A generalisation of the above to finite strain is written in terms of the Kirchhoff stress ( $\boldsymbol {\unicode {x03C4}}=\tilde {\boldsymbol {\unicode {x03C3}}} J_s$ ) on the surface as

(2.10) \begin{equation} \boldsymbol {\unicode {x03C4}}(t) = \boldsymbol {\unicode {x03C4}}_0^d(t) + \text {SYM} \left [ \int _0^t \frac {\dot {G}(s)}{G_0} \hat {\mathbf {A}}_t^{-1}(t-s) \cdot \boldsymbol {\unicode {x03C4}}_0^d(t-s) \cdot \hat {\mathbf {A}}_t(t-s) {\rm d}s \right ] +\tau ^v(t)\mathbf {I}, \end{equation}

where $\boldsymbol {\unicode {x03C4}}_0^d(t) = \tilde {\boldsymbol {\unicode {x03C3}}}_0^d(t)J_s$ , $\tilde {\boldsymbol {\unicode {x03C3}}}_0^d(t)$ is the instantaneous deviatoric part of $\tilde {\boldsymbol {\unicode {x03C3}}}$ in (2.5) and (2.6) with $G_s$ replaced by $G_0$ , $\tau ^v(t)$ is the corresponding volumetric stress, $\hat {\mathbf {A}}_t(t-s)=\partial \boldsymbol {x}(t-s)/\partial \boldsymbol {x}(t)$ is the deformation gradient of the configuration at time $(t-s)$ relative to time $t$ , and SYM indicates symmetry.

The simplest viscoelastic constitutive models are represented by a single linear spring and dashpot connected either in series (Maxwell model) or in parallel (Kelvin–Voigt model) (Dill 2006). These models are commonly used; however, they are somewhat limited as the Kelvin–Voigt model does not exhibit instantaneous elasticity (Dill 2006) and the Maxwell model does not show creep behaviour (Mal et al. 2024). As such, the viscoelastic constitutive model chosen for the wall material is the Maxwell form of the standard linear solid (SLS) model, also known as the Kelvin or Zener model (Fung 1981; Mal et al. 2024). This is a special case of the generalised Maxwell model (a parallel arrangement of a single spring and $N$ Maxwell elements) with a single Maxwell element, as shown in figure 1(b). The SLS model is capable of adequately describing the creep, relaxation and instantaneous elastic behaviour of the material with the minimum number of material constants. Both springs follow Hooke’s law. The single-spring branch represents the wall material shear elastic modulus $G_s$ , while the dashpot represents the solid viscosity $\mu _s$ . Thus, $N=1$ , $G_\infty =G_s$ and $G(t)=G_s+G_1e^{(-t/t_1^G)}$ , where the relaxation time $t_1^G=\mu _s/G_1$ . In the limit $G_1\rightarrow \infty$ , the SLS model becomes the Kelvin–Voigt model, which is represented by the spring $G_s$ and dashpot $\mu _s$ in parallel.

A finite element method is used to obtain $\boldsymbol {\unicode {x03C4}}$ (Shrivastava & Tang 1993). The surface of the vessel is discretised using triangular elements, the vertices (nodes) of which make up a Lagrangian framework (figure 1 c). For each element, a local coordinate system in the element plane is introduced. Then, $\hat {\mathbf {A}}$ , $\lambda _1$ , $\lambda _2$ and $\boldsymbol {\unicode {x03C4}}$ for each element are computed. The integral in (2.10) can be evaluated over a small time interval $\Delta t$ assuming $\boldsymbol {\unicode {x03C4}}_0^d(t)$ varies linearly.

For the fluid–structure coupling, it is useful to obtain the 2-D force at each Lagrangian node. The 2-D force $\boldsymbol {f}_i^{\varDelta }$ at the vertex $i$ of an element is obtained as

(2.11) \begin{equation} \boldsymbol {f}_i^{\varDelta }=\frac {\partial H_i}{\partial \hat {\boldsymbol {V}}}\cdot \boldsymbol {\unicode {x03C4}}\cdot {\hat {\mathbf {A}}}^{-T}, \end{equation}

where $\hat {\boldsymbol {V}}$ represents the local coordinate, $i=1,2,3$ are the three vertices of the element and $H_i$ is the shape function of vertex $i$ . The local force is then transferred to the global coordinate using the appropriate transformation matrix. Then, the resultant viscoelastic force density $\boldsymbol {f}_{ve}$ at any vertex is obtained by the vector resultant of the forces $\boldsymbol {f}_m$ contributed by $M$ surrounding elements which share that vertex:

(2.12) \begin{equation} \boldsymbol {f}_{ve}=\sum _{m\in M}{\boldsymbol {f}_m}. \end{equation}

2.2. Flow dynamics and fluid–structure coupling

The fluid is assumed to be incompressible, and the fluid motion is governed by the continuity and Navier–Stokes equations as

(2.13) \begin{equation} \nabla \cdot \boldsymbol {u}=0 \quad \text {and}\quad \rho \frac {D\boldsymbol {u}}{Dt} = -\nabla P + \mu \nabla ^2\boldsymbol {u}. \end{equation}

These equations are solved in the entire computation domain (figure 1 a), which is discretised using an Eulerian rectangular mesh of uniform size.

A hybrid of the sharp and diffuse interface immersed boundary methods is used for the fluid–structure coupling taking into account the advantages of each method. Details of this method are given by Krul & Bagchi (2024). For the axial direction ( $x$ ), no wall motion is allowed and the zero-velocity condition needs to be enforced on the $x$ -component of velocity. This is done using the sharp-interface ghost-node immersed boundary method (Mittal et al. 2008). The velocity constraint is imposed at the Eulerian nodes immediately neighbouring the vessel exterior, termed the ghost nodes (GN), as

(2.14) \begin{equation} u_{GN} = 2u_{w} - u_{IP}, \end{equation}

where $u$ represents the axial velocity component, $u_w$ is set to zero for no axial motion and $u_{IP}$ is the velocity of the image point which mirrors a GN over the vessel wall. Here, $u_{IP}$ can be obtained by interpolating the surrounding fluid velocity.

For the two normal directions ( $y$ and $z$ ), the tube is allowed to freely move and the continuous forcing immersed boundary method is used (Peskin 2002). In this approach, the force due to wall deformation is introduced in the Navier–Stokes equations as

(2.15) \begin{equation} \rho \frac {D\boldsymbol {u}}{Dt} = -\nabla P + \mu \nabla ^2\boldsymbol {u} + \boldsymbol {F}, \end{equation}

where

(2.16) \begin{equation} \boldsymbol {F} = \begin{bmatrix} 0 \\ {\displaystyle \int \limits _S{f_{ve_y}\delta (\boldsymbol {x} - \boldsymbol {x}^\prime ,t) {\rm d}\boldsymbol {x}^\prime } } \\ {\displaystyle \int \limits _S{f_{ve_z}\delta (\boldsymbol {x} - \boldsymbol {x}^\prime ,t) {\rm d}\boldsymbol {x}^\prime } } \end{bmatrix}, \end{equation}

where $f_{ve_y}$ and $f_{ve_z}$ are the $y$ and $z$ components of $\boldsymbol {f}_{ve}$ , $\boldsymbol {x}^\prime$ is any location on the vessel surface $S\in \mathbb {R}^3$ , $\boldsymbol {x}\in \mathbb {R}^3$ is the Eulerian variable, and $\delta$ is the 3-D Dirac delta function (Peskin 2002). A discrete form with a finite spread in $\mathbb {R}^3$ is used to numerically evaluate delta as (Peskin 2002)

(2.17) \begin{equation} \delta (\boldsymbol {x} - \boldsymbol {x}^\prime) \approx D (\boldsymbol {x} - \boldsymbol {x}^\prime) = \begin{cases} \frac {1}{\left ( 2\hslash \right )^3} \prod\limits_{i=1}^{3}{ \left [ 1 + \cos { \left \{ \frac {\pi }{\hslash } ( x_i - x_i^\prime ) \right \} } \right ] } & \text {if } | x_i - x_i^\prime | \le \hslash , \\ \hfil 0 & \text {otherwise}, \end{cases} \end{equation}

where $\hslash =2\varDelta _E$ and $\varDelta _E$ is the size of the unit Eulerian cell.

The vessel wall is advected in the $y$ and $z$ directions using the $v$ and $w$ components of the local fluid velocity as

(2.18) \begin{equation} \frac {{\rm d}\boldsymbol {x}^\prime }{{\rm d}t} = \boldsymbol {u}_w (\boldsymbol {x}^\prime , t), \end{equation}

where

(2.19) \begin{equation} \boldsymbol {u}_w = \begin{bmatrix} 0 \\ \displaystyle \int \limits _{\boldsymbol {x}}{v\delta (\boldsymbol {x} - \boldsymbol {x}^\prime ,t) {\rm d}\boldsymbol {x} } \\ \displaystyle \int \limits _{\boldsymbol {x}}{w\delta (\boldsymbol {x} - \boldsymbol {x}^\prime ,t) {\rm d}\boldsymbol {x}} \end{bmatrix}. \end{equation}

As the vessel wall deforms, some Eulerian mesh points exterior of the vessel may enter the vessel interior and vice versa. As a result, GNs and IPs must be updated at every step of the time marching for which a fast and parallelisable algorithm was developed (Krul & Bagchi 2024).

The flow solver is based on a four-step projection-based method for incompressible flows. A staggered-grid implementation is considered for the fluid velocity components and pressure. The diffusion terms are treated using the semi-implicit Crank–Nicolson scheme, and the nonlinear and force coupling terms are treated explicitly using the second-order Adams–Bashforth scheme:

(2.20) \begin{align} & \textit {Advection-diffusion}: && \rho \frac {\hat {\boldsymbol {u}}-\boldsymbol {u}^n}{\Delta t}=-\nabla P^n+\frac {\mu }{2}[\nabla ^2\hat {\boldsymbol {u}}\ +\nabla ^2\boldsymbol {u}^n]-\left [\frac {3}{2}\boldsymbol {N}^n-\frac {1}{2}\boldsymbol {N}^{n-1}\right ], \end{align}

(2.21) \begin{align} \notag &&& \text {where } \boldsymbol {N} = \rho \boldsymbol {u}\cdot \nabla \boldsymbol {u}-\boldsymbol {F}, \\ & \textit {Intermediate velocity}: && \rho \frac {\boldsymbol {u}^\ast -\hat {\boldsymbol {u}}}{\Delta t}=\nabla P^n, \end{align}

(2.22) \begin{align} & \textit {Poisson equation}: && \nabla ^2P^{n+1}=\rho \frac {\nabla \cdot \boldsymbol {u}^\ast }{\Delta t}, \end{align}

(2.23) \begin{align} & \textit {Pressure correction}: && \rho \frac {\boldsymbol {u}^{n+1}-\boldsymbol {u}^\ast }{\Delta t}=-\nabla P^{n+1}. \end{align}

The spatial derivatives in the advection–diffusion equation are treated using second-order discretisation. An alternating direction implicit (ADI) scheme is used to solve this equation.

2.3. Dimensionless parameters

The undeformed vessel radius $R_0$ is taken as the length scale and the centreline velocity $U_c$ of the Poiseuille flow in the undeformed vessel as the velocity scale. The pressure is scaled by $\rho U_c^2$ . Dimensionless variables are indicated by a $\ast$ . We introduce the ‘baseline’ pressure drop $\Delta P_0$ in the undeformed tube such that $U_c^\ast =1$ . The corresponding flow rate is denoted by $q_0$ . Also, we use a $\sim$ to denote normalised variables. Then,

(2.24) \begin{equation} \tilde P = \frac {P}{\Delta P_0}, \quad \tilde q = \frac {q}{q_0}, \end{equation}

where $P=P(x,t)$ and $q=q(x,t)$ are the pressure and flow rate in the deforming tube, respectively.

The oscillatory inlet pressure is prescribed as

(2.25) \begin{equation} \tilde P_{in} = \frac {P_{in}}{\Delta P_0} = \frac {1}{2}-\frac {1}{2}\cos (2\pi f^\ast t^\ast ), \end{equation}

where $f^\ast = f({R_0}/{U_c})$ is the dimensionless oscillation frequency, $T$ is the period and $t^\ast = t({U_c}/{R_0})$ . The outlet and external pressures are both kept at $0$ .

The relevant dimensionless parameters are

\begin{align*} \begin{array}{lll} &\displaystyle \text {Reynolds number} & Re = \frac {\rho U_c R_0}{\mu }\\ &\displaystyle \text {Womersley number} & \alpha = R_0\sqrt {\frac {2\pi f\rho }{\mu }}\\ &\displaystyle \text {Dimensionless wall elasticity (or deformability parameter)}\quad &\beta = \frac {\mu U_c}{G_s} \\ &\displaystyle \text {Deborah number} & De = \frac {\mu _s U_c}{G_s R_0} \\ \end{array} \end{align*}

The Reynolds number is the ratio of the inertial to viscous forces in the fluid, and the Womersley number represents the ratio of the fluid viscous and applied oscillation time scales. The deformability parameter is the ratio of the fluid viscous forces to the elastic forces and is a measure of tube compliance. The Deborah number is the ratio of the viscoelastic relaxation time $\tau _R={\mu _s}/{G_s}$ to the inertial time scale $\tau _I={R_0}/{U_c}$ . It is generally used for viscoelastic fluids, but in this case, it represents the strength of the wall viscosity, where $De=0$ corresponds to an elastic tube. The ranges of these parameters are shown in table 1 along with the tube geometry.

Table 1. Geometric and dimensionless parameter values.

2.4. Additional numerical details

The computational domain has dimensions $4\pi \times 1.48\pi \times 1.48\pi$ and consists of $280\times 104\times 104$ grid points. The tube surface is discretised by $55\,224$ triangular elements and the time step chosen is $\Delta t^\ast =5\times 10^{-4}$ .

A detailed convergence study of the mesh and domain size is presented by Krul & Bagchi (2024), which showed that the chosen mesh yields converged results. As an illustration of the robustness of the method, the tube surface mesh is shown in figure 1(c) for an extreme case for which high-wavenumber folds appear on the surface as the tube inflates to maximum. The folds are well resolved, and the mesh remains stable even after this and throughout multiple oscillation cycles.

Furthermore, the numerical results for a steady inflation of an elastic tube under a constant pressure were compared by Krul & Bagchi (2024) against the small deformation analytical models of Fung (1997), which are based on the law of Laplace, and Anand & Christov (2021), as well as a second-order analytical model accounting for large deformation. Excellent agreement between the numerical results and the second-order analytical model is observed for inflation as large as $75\,\%$ of the undeformed radius.

Additional validation for oscillatory pressure is given in § 3.1.

3. Results

The interest of this study is in the effect of wall viscosity and external oscillation at weak fluid inertia. As such, we restrict to $Re=0.1$ and vary $\beta$ , $De$ and $\alpha$ . We also set $G_1=100G_s$ in the following, except in § 3.3.6, where the effect of $G_1$ is presented.

Deformation of the tube is characterised by the circumferential strain defined as

(3.1) \begin{equation} \varepsilon (x,t) = \frac {R(x,t)-R_0}{R_0}, \end{equation}

which varies both in time and in the axial direction. The maximum strain at a given time instance $t$ , denoted as $\varepsilon _{max}(t)$ , is also used to characterise the tube deformation. The operators $\text {Max}\left \{\cdot \right \}$ , $\text {Min}\left \{\cdot \right \}$ and $\text {Avg}\left \{\cdot \right \}$ are the maximum, minimum and average values, respectively, with respect to time over one oscillation cycle. We further define an ‘amplitude’ as $\text {Amp}\left \{\cdot \right \} = \text {Max}\left \{\cdot \right \}-\text {Min}\left \{\cdot \right \}$ . To quantify the deformation and flow rate, we introduce the following notation:

\begin{align*} \begin{array}{@{}llcll} &\displaystyle \text{Maximum}\; \varepsilon \text{ over one period}\qquad &&&\Xi _{max} = \text {Max}\left \{\varepsilon _{max}\right \} \\ &\displaystyle \text {Amplitude of } \varepsilon \qquad &&&\Delta \Xi =\text {Amp}\left \{\varepsilon _{max}\right \} \\ &\displaystyle \text {Maximum inlet flow rate over one period}\qquad &&&Q_{max}=\text {Max}\{\tilde {q}_{in}\} \\ &\displaystyle \text {Amplitude of inlet flow rate}\qquad & &&\Delta Q=\text {Amp}\{\tilde {q}_{in}\} \end{array} \end{align*}

The model is first tested against small deformation theory in § 3.1, then the effects of varying $\alpha$ and $\beta$ are explored in § 3.2. The tube is purely elastic ( $De=0$ ) in both sections. This is followed by varying $De$ , $\alpha$ and $G_1$ in viscoelastic tubes in § 3.3.

3.1. Comparison with small deformation theory

Figure 2. Comparison against Womersley’s theory for an elastic tube. (a) Centreline velocity over time taken halfway through the tube. (b) Centreline velocity along the tube length at different time instances. The solid and dashed lines are the numerical and theoretical solutions, respectively. The different colours correspond to different time instances throughout the oscillation. Parameter: $\beta = 5\times 10^{-3}$ .

Figure 2 shows our simulation results compared against Womersley’s linearised analytical solution for an elastic tube (Womersley 1957a ; Zamir 2016a ). Here, we define $U=u_{cen}-\bar {u}_{cen}$ as the oscillatory component of the axial centreline velocity, where $\bar {u}_{cen}$ is the time-averaged centreline velocity. As per Womersley’s solution,

(3.2) \begin{align} U(x,t) = \text {Re}\left \{ \frac {B}{\rho c}\left [ 1-G\frac {1}{J_0(\Lambda )} \right ]e^{i\omega (t-x/c)} \right \}, \end{align}

where

(3.3) \begin{align} G = \frac {2+z(2\nu _s-1)}{z(2\nu _s-g)}, \quad g = \frac {2J_1(\Lambda )}{\Lambda J_0(\Lambda )}, \end{align}

(3.4) \begin{align} \Lambda = \left (\frac {i-1}{\sqrt {2}}\right )\alpha , \end{align}

$B$ is the amplitude of the pressure oscillation, $J_n$ are the Bessel functions of the first kind and Re $\{\cdot \}$ indicates the real part. We obtain $z$ by finding the roots of the equation

(3.5) \begin{equation} [ (g-1)(\nu _s^2-1) ]z^2 + \left [ \frac {\rho _w h}{\rho R_0}(g-1) + \left ( 2\nu _s-\frac {1}{2} \right )g-2 \right ]z + \frac {2\rho _w h}{\rho R_0} + g = 0, \end{equation}

from which the wave speed $c$ can be calculated as

(3.6) \begin{align} & c = \sqrt {\frac {2}{(1-\nu _s^2)z}}c_0, \end{align}

(3.7) \begin{align} & c_0^2 = \frac {(1+\nu _s)G_s}{\rho R_0}, \end{align}

where $\rho _w$ is the wall material density, $h$ is the wall thickness and $c_0$ is the Moens–Korteweg wave speed in an inviscid fluid. The larger of the two $z$ roots is used. We select $\beta = 5\times 10^{-3}$ in our simulations so that deformation remains small ( $\sim$ 5 % change in tube radius). Figure 2(a) shows the predicted oscillatory velocity waveform taken at the midpoint of the tube. A good agreement with the theory is evident. Figure 2(b) plots the centreline velocity along the tube axis at different time instances. These show decent agreement with some discrepancies that can be attributed to differences in the problem set-up. Our simulation uses a short tube with rigid inlet and outlet lengths, whereas Womersley’s solution assumes a long, fully elastic tube and linearised flow equations with the no-slip condition applied on the undeformed tube surface.

3.2. Elastic tube

Figure 3. General tube behaviour. (a)–(d) Axial velocity contours and streamlines over one period. (e)–(h) Pressure contours at the same time instances. Parameters: $De=0$ , $\alpha =0.22$ .

The time-periodic motion of an elastic tube under pulsatile flow is shown in figure 3 for one oscillation period. Note that there is an initial transience during the simulations, but all results are taken after the tube has settled into a time-periodic motion. During inflation, the increasing inlet pressure drives the fluid into the tube, creating a net flow and causing it to distend. Initially, the undeformed tube offers little resistance to the deformation and it therefore expands rapidly near the inlet. The flow rushes in to compensate, creating a flow surge at the inlet seen in figure 3(b). The surge causes the tube to over-inflate, creating a large restoring elastic force at its maximum inflation. As the external pressure decreases, this elastic force causes the tube to deflate rapidly, squeezing fluid out of the outlet (figure 3 c) and increasing the pressure inside the deforming section (figure 3 h). This pressure drives the fluid out of both ends of the tube once the applied pressure has decreased enough, causing a negative flow rate at the inlet (figure 3 d).

Note that the predicted flow is not unidirectional and pressure varies in the radial direction.

Figure 4. Elastic tube. (a) Circumferential strain over one period of oscillation. The inset indicates the time instances when the profiles were taken, with each point corresponding to the line of the same colour. The filled and empty circles correspond to the solid and dashed lines, respectively. Solid lines are during inflation; dashed lines are during deflation. Parameter: $\beta =0.09$ . (b),(c) Transient waveforms of (b) $\varepsilon _{max}$ and (c) inlet flow rate over one period for different $\beta$ . (d) $\Xi _{max}$ , $\Delta \Xi$ and maximum radius phase shifts (in radians) plotted against $\beta$ . (e) Maximum, amplitude and phase shifts (in radians) of the inlet flow rate versus $\beta$ . Parameter: $\alpha =0.5$ .

Figure 4(a) shows the $\varepsilon (x,t)$ profiles at several time instances. The maximum inflation occurs near the inlet of the tube experiencing higher pressure. The rate of inflation in this part is also more rapid than the rest of the tube, causing a relatively sharp peak. Downstream of this peak, the profile has concave curvature. This is not the case during deflation, where the profile has a more consistent curvature that is convex everywhere. The location of the maximum deformation gradually moves downstream during deflation, then rapidly jumps to the upstream end as the next inflation phase begins. Also, large deformation of the tube is evident as $\varepsilon _{max}$ is approximately $0.5$ .

In figure 4(b), $\varepsilon _{max}(t)$ is plotted for different $\beta$ along with the $\tilde {P}_{in}(t)$ waveform. The corresponding waveforms of $\tilde {q}_{in}(t)$ are shown in figure 4(c). As seen, deformation and flow rate increase with increasing $\beta$ . Figures 4(d) and 4(e) show $\Xi _{max}$ , $\Delta \Xi$ , $Q_{max}$ and $\Delta Q$ as functions of $\beta$ (left axis in the figures). The average deformation increases rapidly, while the amplitude of deformation saturates at relatively smaller $\beta$ (approximately $\beta =0.04$ ). This behaviour is seen because the inflation is controlled by the volume flow rate. At a higher average deformation, the volume required to inflate by the same amplitude is larger. Thus, the flow rate would need to increase significantly to accommodate that inflation within the oscillation period. The maximum and amplitude of the flow rate similarly increase with $\beta$ and tend to plateau for more compliant tubes, as seen in figure 4(e).

Figure 4(b) also indicates that there is a phase lag between the applied pressure and deformation. This is a consequence of the response time of the tube, as the tube and fluid flow take time to adjust to an equilibrium. This lag is represented by an apparent rightward shift of the $\varepsilon _{max}(t)$ waveform with respect to $\tilde {P}_{in}(t)$ . Unlike the deformation, however, the flow rate leads the applied pressure, as seen by a leftward shift in figure 4(c), even though the pressure is driving the flow. This is because the maximum flow corresponds to the maximum pressure gradient (Pedley 1980), which is not, in general, coincident with the maximum applied pressure.

To quantify the phase lag or lead, we write both the strain and flow rates as functions similar to the applied pressure, but with an added variable phase shift, $\phi (t^\ast )$ and $\psi (t^\ast )$ , as

(3.8) \begin{align} &\tilde q_{in} = \frac {\Delta Q}{2} \left [ 1-\cos \left (2\pi f^\ast t^\ast + \phi (t^\ast ) \right ) \right ] + \text {Min}\{\tilde q_{in}\}, \end{align}

(3.9) \begin{align} &\varepsilon _{max} = \frac {\Delta \Xi }{2} \left [ 1-\cos \left (2\pi f^\ast t^\ast + \psi (t^\ast ) \right ) \right ] + \text {Min}\{\varepsilon _{max}\}, \end{align}

where $-\pi \lt \phi (t^\ast ),\psi (t^\ast )\le \pi$ and are measured in radians. We do this by noting that both the strain and flow rate are periodic with the same frequency as the pressure. By this construction, the signs of $\phi (t^\ast )$ and $\psi (t^\ast )$ indicate whether the flow rate and deformation are leading $(\gt 0)$ or lagging behind $(\lt 0)$ the applied pressure at time $t^\ast$ . The time dependence of the phase shifts allows for deviations from simple harmonic motion, which makes (3.8) and (3.9) valid for any arbitrary periodic strain or flow rate waveform. As seen in figures 4(b) and 4(c), the phase shifts are varying in time. To concisely present the results, we select the phase shifts at the time instances corresponding to the maximum (crest) and minimum (trough) applied pressure and denote them as $\psi _c$ and $\psi _t$ for deformation, and $\phi _c$ and $\phi _t$ for flow rate. As seen in figure 4(d), both $\psi _c$ and $\psi _t$ are negative, and their magnitudes increase with increasing $\beta$ . Thus, the phase lag of deformation increases with more compliant vessels, given the higher deformability contributes to a longer response time. In contrast, $\phi _c$ and $\phi _t$ are positive, as seen in figure 4(e), and they decrease with increasing $\beta$ . As discussed by Womersley (1957b ), the maximum flow lead of $\pi /4$ occurs in the limiting case of a rigid tube, which is close to $\phi _c=0.554$ and $\phi _t=0.612$ for the lowest $\beta =0.005$ considered in the figure. The increasing deformability acts to reduce this lead.

It should be noted that the uniform leading of the flow rate for varying $\beta$ as predicted here (and by Womersley (1957b )) does not always occur, as the phase shift changes with the tube’s deformability, viscosity and oscillation frequency as discussed later.

Figure 5. Hysteresis in elastic tubes for different $\beta$ . (a) Hysteresis of the radius in an elastic tube with respect to the local centreline pressure, $x_R=3.36$ . (b) Hysteresis of the flow rate. The legend for each plot is in panel (b). Parameter: $\alpha =0.5$ .

The phase shift manifests a hysteresis in the time periodic motion. This is shown in figure 5(a) by plotting $\varepsilon (x_R,t)$ versus $P(x_R,t)$ , where $x_R$ is a fixed location near the maximum radius when the tube is fully inflated. The curves are traced out in the counterclockwise direction, as indicated by the arrows on the figure, indicating the deformation lags behind the applied pressure. The hysteresis curves get wider with increasing $\beta$ , indicating both the increasing amplitude of deformation and phase lag.

Note that it is generally expected that no hysteresis would occur in a purely elastic tube (Kim et al. 2016; Čanić et al. 2006a ). In the linearised analytical theory, the pressure varies only along the axial direction and the inflation at any $x$ depends on the local centreline pressure. Indeed, for the lowest $\beta =0.005$ considered, hysteresis is nearly absent. However, at large deformation, pressure varies radially, as was shown in figure 3, instead of being a function of axial distance alone. This effect increases with the deformability, creating a hysteresis even for elastic tubes.

Figure 5(b) presents the hysteresis in terms of the inlet flow rate and pressure. Here, the direction is clockwise because of the phase lead. Also, unlike the deformation, the flow rate hysteresis loop becomes thinner with increasing $\beta$ , due to the decreasing lead as noted earlier in figure 4(e).

Figure 6. Elastic tube. (a) $\Xi _{max}$ , $\Delta \Xi$ and deformation phase shifts plotted against $\alpha$ . (b) Maximum, amplitude and phase shifts of the inlet flow rate versus $\alpha$ . (c) Hysteresis of the radius in an elastic tube with respect to the local centreline pressure, $x_R=3.36$ . (d) Hysteresis of the flow rate. The legend for both panels (c) and (d) is in panel (d). Parameter: $\beta =0.05$ .

Next, we consider the effect of changing the Womersley number while keeping $\beta$ constant (figure 6). Here, $\Xi _{max}$ , $\Delta \Xi$ , $\psi _c$ and $\psi _t$ all decrease with increasing $\alpha$ , as plotted in figure 6(a). This is because the tube motion cannot keep up with the faster variation of the applied pressure. In the limit of $\alpha \rightarrow \infty$ , the tube would behave as if it were steadily inflated by the mean pressure ( $\tilde P_{in}=0.5$ ). Therefore, we should expect the amplitude of oscillation to decay to zero with increasing $\alpha$ , as observed. In contrast to the radius, the maximum volume flow rate at the inlet increases with $\alpha$ , as seen in figure 6(b). This is due to the flow surge during inflation, whose strength increases with $\alpha$ . The more rapidly increasing pressure means that a higher pressure is applied at lower deformation, and thus, low opposition from the elastic force. Further increasing $\alpha$ approaches the limit of the tube’s response time, causing both the maximum and amplitude of the inlet flow rate to plateau. Here, $\phi _c$ and $\phi _t$ initially increase at low $\alpha$ , after which they both decrease. In the quasi-steady limit, when both phase shifts are $0$ , a higher applied pressure corresponds to a higher inlet pressure gradient and, therefore, flow rate. This changes when increasing $\alpha$ as the tube’s rapid distention lowers the downstream pressure, moving the maximum pressure gradient earlier in the cycle. Increasing $\alpha$ increases the strength of the effect, further increasing the phase lead. However, increasing $\alpha$ also increases the response lag. This takes precedence at higher $\alpha$ and the phase shifts begin to decrease. This eventually causes a uniform lag when the phase shifts drop below zero. Additionally, both $\phi _c$ and $\phi _t$ trend towards each other, suggesting the phase shift $\phi (t^\ast )$ approaches a constant value for high $\alpha$ .

The deformation shows very weak hysteresis for different $\alpha$ , as seen in figure 6(c). In contrast, the flow rate hysteresis is noticeable in figure 6(d). The vertical tilt of the hysteresis curves increases with $\alpha$ as the strength of the flow surge increases. This also causes the curves to initially widen. Beyond a certain point, as $\alpha$ increases, the lag in the tube’s response causes the curves to start to narrow. In most cases, the hysteresis curves are either ellipses or narrow ovals; however, at very low $\alpha$ , the curve takes on a banana-like shape. This is because it is converging to the quasi-steady inflation curve for inlet flow rate versus applied pressure, for which the trend is nonlinear, as was shown in previous studies (Fung 1997; Anand & Christov 2021; Krul & Bagchi 2024).

3.3. Viscoelastic tube

In this section, we set $\beta =0.05$ throughout, and focus on the effect of $De$ and $\alpha$ . We first investigate the tube deformation and phase lag in § 3.3.1, then the flow rate and deformation hysteresis in § 3.3.2. The pressure propagation is explored in § 3.3.3, followed by the effect of oscillation frequency in § 3.3.4. Finally, we consider the interaction of viscoelasticity and oscillation frequency in § 3.3.5, and the effect of $G_1$ in § 3.3.6.

3.3.1. Deformation and phase lag

The qualitative behaviour of a purely elastic tube, as seen in figure 3, remains similar for a viscoelastic tube. However, a viscoelastic tube resists rapid motion, such as the sudden expansion during inflation and the rapid contraction during deflation. As such, the wall viscous forces will, in general, mitigate the effects of both the flow surge and the squeezing.

Figure 7. Viscoelastic tube. Time-varying tube profiles for (a) $De=5$ and (b) $De=25$ . The profiles are taken at the time instances indicated by the inset in panel (a). (c) $\Xi _{max}$ , $\Delta \Xi$ and maximum radius phase shifts plotted against $De$ . (d) Maximum, amplitude and phase shifts of the inlet flow rate versus $De$ . Parameter: $\alpha =0.16$ .

Figure 7 shows $\varepsilon (x,t)$ profiles for different $De$ . At higher $De$ , the profile downstream of the maximum deformation tends towards linear during both inflation and deflation. This is unlike the elastic case (figure 4 a), whose downstream curvature is concave during inflation and convex during deflation. The wall viscosity causes the motion to change from very localised to more global, where the entire tube tends to move simultaneously. This type of motion is the most favourable for the highly viscous wall, as it allows for the slowest changes in local deformation.

At higher $De$ , both the variation in radius and the maximum radius decrease, agreeing with the observations made by Wang et al. (2015). Both $\Xi _{max}$ and $\Delta \Xi$ are shown in figure 7(c) against $De$ , which shows this decaying trend. The strain amplitude trends towards zero, while the maximum strain approaches some non-zero constant value. This is consistent with the limit of $De\rightarrow \infty$ , when the tube would act as rigid with some constant deformed configuration.

Furthermore, at higher values of $De$ , such as in figure 7(b), folds appear on the vessel wall near the maximum radius during the inflation (also in figure 1 c). These are most pronounced in the middle of the inflation phase. This is when the tube deforms the fastest and the wall viscous forces are the highest. This behaviour is characteristic of the fluid-driven motion of a viscous membrane (Yazdani & Bagchi 2013).

One of the main consequences of increasing the wall viscosity is increasing the tube’s response time, causing the deformation to lag more behind the applied pressure (Čanić et al. 2006b ). The phase lag increases drastically at low $De$ before levelling out near $De\approx 50$ , as seen by the decreasing $\psi _c$ and $\psi _t$ in figure 7(c).

As $De$ increases, the emergence of folds on the tube wall affects the apparent phase shift at the crest and trough, leading to the sudden jump in $\psi _c$ and increase in $\psi _t$ .

Corresponding to the variation in tube wall motion, the inlet volume flow rate similarly varies with $De$ . The maximum flow rate and amplitude decrease with increasing $De$ , as shown in figure 7(d). The viscous resistance mitigates both the flow surge during inflation and the squeezing during deflation, lowering both the maximum and amplitude of the flow rate. The maximum and amplitude trend towards each other at high $De$ . This is because in the limit of infinite $De$ , the tube would behave as rigid, and $\tilde {q}_{in}$ would oscillate between zero and its maximum value, hence, the maximum and amplitude would be identical. Also plotted in figure 7(d) are $\phi _c$ and $\phi _t$ with both showing decaying trends. They are both positive for $De\le 15$ , indicating the flow rate leads the pressure. For $De\gt 25$ , $\phi _t$ remains positive but $\phi _c$ attains a small negative value. This indicates that, for a high enough $De$ , $\tilde {q}_{in}$ lags behind $\tilde {P}_{in}$ at or near the maximum inflation, as the high wall viscosity delays the tube motion.

3.3.2. Hysteresis

Figure 8. (a) Deformation and (b) flow rate hysteresis for varying $De$ . Inset shows zoomed-in cross-over behaviour. For both panels (a) and (b), the legend is in panel (a). Parameter: $\alpha =0.16$ .

Figure 8(a) shows the deformation hysteresis at varying $De$ , which is characteristic of the viscoelastic tube (Čanić et al. 2006b ). When increasing $De$ , the upward tilt of the curves decreases, signifying a decrease in the amplitude of oscillation, agreeing with the previously presented results. Also, the curves initially widen and then narrow. One can postulate that at infinite $De$ , the curve will become a horizontal line, as there is no tube motion.

Figure 8(b) shows the flow rate hysteresis for various $De$ . Unlike the deformation hysteresis, these curves proceed clockwise. When increasing $De$ , the tube motion is delayed and the curves become narrower. Their vertical tilt also decreases due to the mitigated flow surge, decreasing the amplitude of $\tilde {q}_{in}$ . Interestingly, at $De=25$ , the hysteresis curve creates a sharp point at the high-pressure end and further increasing the viscosity causes the curve to cross over itself. The inset in figure 8(b) shows a zoomed-in view of this cross-over with arrows indicating the direction in which the curve is traced out. Note that this transition corresponds to $\phi _c$ first dropping below zero.

Figure 9. Inlet flow rate hysteresis and phase shift. The flow rate versus the applied pressure is shown for (a) $De=5$ , $\alpha =0.08$ ; (b) $De=50$ , $\alpha =0.08$ ; (c) $De=1$ , $\alpha =1.06$ . (d–f) Phase shift over time for the same parameters as panels (a), (b) and (c), respectively. Note that $Q_{min}=\text {Min}\{\tilde {q}_{in}\}$ .

This directional change of hysteresis curves was also observed by Lee & Schmid-Schönbein (1990). As such, the hysteresis may proceed clockwise, counterclockwise or a combination in which the curve crosses over itself, creating a figure-of-eight. Figure 9 gives an example of each type of hysteresis as predicted by our study. The direction of hysteresis depends on the strength of the viscous (lagging) effects relative to the inherent flow lead, which depends on both $De$ and $\alpha$ . The figure-of-eight cross-over occurs when the influence of both the lagging and leading effects are similar in strength.

One notable feature is that the direction of hysteresis can be predicted using the phase shifts $\phi _c$ and $\phi _t$ . Figures 9(d)–9(f) show the variable phase shifts plotted alongside the waveforms of applied pressure and inlet flow rate. The inlet flow rate has been shifted and normalised to more clearly represent the phase shift. If both $\phi _c$ and $\phi _t$ are positive, the corresponding hysteresis curve proceeds clockwise and vice versa. The cross-over hysteresis in figure 9(b) corresponds to a negative $\phi _c$ and positive $\phi _t$ . The above can be summarised in the following compact form:

(3.10) \begin{align} & \phi _c \gt 0, \phi _t \gt 0 \implies \text {Clockwise hysteresis}, \end{align}

(3.11) \begin{align} & \phi _c \lt 0, \phi _t \gt 0 \implies \text {Cross-over hysteresis}, \end{align}

(3.12) \begin{align} & \phi _c \lt 0, \phi _t \lt 0 \implies \text {Counterclockwise hysteresis}. \end{align}

This holds true for the entire parameter space considered.

3.3.3. Pressure propagation

Figure 10. Axial velocity and pressure variation along the tube centreline. (a),(b) Centerline pressures for (a) $De=0$ and (b) $De=50$ . (c),(d) Corresponding centreline velocities. The dashed black lines show the upper and lower bounds of the dependent variable over the full oscillation cycle. The lines proceed in time in the order solid blue, red, green, dash-dotted blue, red and green, at times $t/T={3}/{12},{4}/{12},{5}/{12},{9}/{12},{10}/{12},{11}/{12}$ , respectively. Parameter: $\alpha =0.75$ .

Another feature that is altered by $De$ is the propagation of the pressure through the tube. To illustrate this effect, figure 10 shows the centreline pressure and velocity for two different tubes: one that is purely elastic and another that is highly viscous. The black dashed lines are the upper and lower bounds of the centreline velocity and pressure over the course of the cycle. Together, they create an envelope that encloses all possible values throughout the oscillation. The height of this envelope at any axial location is indicative of the strength of the velocity or pressure oscillation. In the elastic tube in figure 10(a), the rapid distention near the inlet damps the applied oscillations and prevents them from travelling further downstream. This leads to rapid inflation and deflation in only a localised region near the inlet of the tube. This behaviour is evident from the rapidly decreasing size of the pressure envelope near the inlet; the large applied pressure oscillations are completely damped halfway through the tube. In contrast, in the high viscosity tube, in figure 10(b), the pressure envelope reduces in size much more gradually, maintaining a non-zero height until the end of the tube is reached. This behaviour shows that the changing upstream pressure quickly propagates downstream and is felt over the entire length of the tube.

Another way to express this effect is that flow pulsatility propagates more effectively in tubes with higher viscosity. This applies to not only the pressure, but also the axial velocity. In figure 10(c), it is apparent that the purely elastic tube has large fluctuations in the flow near the inlet, but these quickly dissipate when moving downstream. The viscoelastic tube in figure 10(d) shows moderate oscillations near the inlet. There is some decay in the envelope, but, in contrast to the elastic tube, the oscillations persist along the entire length of the tube. Prior studies in stiffer tubes show that increased solid viscosity attenuates the pulse wave, especially at higher frequencies (Pontrelli 2002; Mal et al. 2024), suggesting that this increased pulsatility transmission appears only in highly deformable vessels. This may be due to the lowered wave speed and increased elastic damping associated with increased compliance (Zamir 2016b ; Roknujjaman et al. 2023).

3.3.4. Effect of oscillation frequency

Figure 11. (a) Time-varying tube profiles for $\alpha =0.5$ . The profiles advance in time in the order solid black, blue, red, green, dashed black, blue, red and green. These occur at time instances $t/T=0,{2}/{12},{4}/{12},{5}/{12},{6}/{12},{8}/{12},{10}/{12},{11}/{12}$ , respectively. $De=0.5$ . (b) Deformation hysteresis for varying $\alpha$ . $De=5$ .

The effect of varying $\alpha$ at a given $De$ is qualitatively similar to that seen earlier for $De=0$ in figure 6. Two notable behaviours exist when increasing $\alpha$ for the viscoelastic tube at a constant $De$ . First, increasing $\alpha$ causes the downstream portion of the tube to move less than the rest of the tube. This can be seen in figure 11(a) at $\alpha =0.5$ , where all of the profiles converge for $x^\ast \ge 10$ . The higher frequency oscillations are more readily damped by the elastic tube motion, as the pressure cannot build appreciably within each cycle. Note that this effect will change with the tube’s viscosity, as seen above in figure 10. Second, hysteresis of the viscoelastic tube at varying $\alpha$ , as shown in figure 11(b), is stronger than the elastic tube seen earlier in figure 6(c). Also, as $\alpha$ becomes large, the curve will become a horizontal line, as there is no tube motion.

3.3.5. Interaction of viscoelasticity and oscillation frequency

The interaction between the wall viscosity and enforced oscillation frequency creates different phenomena in the parameter space that are presented below.

Figure 12. $\Xi _{max}$ , $\Delta \Xi$ , $Q_{max}$ and $\Delta Q$ for different $De$ and $\alpha$ . The dashed black lines show the asymptotic limits for $\alpha \rightarrow 0$ and $\alpha \rightarrow \infty$ . For all panels, the legend is in panel (b).

In figures 12(a) and 12(b), we can see the combined effects of both $De$ and $\alpha$ on $\Xi _{max}$ and $\Delta \Xi$ . The curves are plotted against $\alpha$ and separated by $De$ . At both very low and very high values of $\alpha$ , they appear to approach constant values. Additionally, the curves for all different $De$ converge in these limits. In the limit as $\alpha \rightarrow 0$ , the motion is quasi-steady, implying that the motion is slow enough for the viscous resistance to be negligible. As such, we would expect that the tube motion is the same, regardless of $De$ . The maximum radius is that of the steady tube inflated by the maximum applied pressure, $\tilde {P}_{in}=1$ . Likewise, $\Delta \Xi$ approaches the same low $\alpha$ limit, given the quasi-steady cycle oscillates between maximum and zero inflation. Both $\Xi _{max}$ and $\Delta \Xi$ quickly decay with increasing $\alpha$ and the amplitude converges to zero. This decrease in motion means that the viscous effects become less relevant at higher $\alpha$ , hence why the different curves converge. Here, $\Xi _{max}$ converges to its value for the steady tube inflated by the mean pressure, $\tilde {P}_{in}=0.5$ . The relevant limits are marked in figures 12(a) and 12(b) as dashed black horizontal lines, showing the data match the convergence behaviour. From this, we can conclude that the wall viscosity has minimal effect on the maximum distention for both very low and very high $\alpha$ , but it has a large impact in the range $0.05\lesssim \alpha \lesssim 1$ .

Figures 12(c) and 12(d) show the combined effect of $De$ and $\alpha$ on $Q_{max}$ and $\Delta Q$ . For low $De$ , the volume flow rate initially increases with $\alpha$ due to the strengthening flow surge effect. Eventually, the oscillating pressure changes too quickly for the whole tube to respond and the full extent of the flow surge is no longer felt, causing $Q_{max}$ to plateau. The peak flow rate decreases for higher $De$ , as the flow surge is mitigated by the wall viscosity. Interestingly, when the wall viscosity is high enough, the flow rate initially decreases with $\alpha$ . In these cases, the tube viscosity is high enough to completely negate the elastic flow surge effect at lower $\alpha$ . The tube’s resistance to motion prevents it from fully deforming to experience the same $Q_{max}$ it had during quasi-steady inflation, thereby causing the decrease. Another important feature of figure 12(c) is that, at high $De$ , the inlet flow rate continues to increase over a large range of $\alpha$ . This counterintuitive behaviour is related to the pressure propagation increasing at higher viscosity (see figure 10). At higher values of viscosity, the instantaneous response approaches that of a rigid tube and the pressure builds throughout the tube before any changes in motion occur. When the viscous resistance subsides and the tube begins to move, the whole tube adjusts to the internal pressure simultaneously, as opposed to only the inlet section moving independently. This whole-tube motion causes a flow surge at the inlet. While the previous flow surge effect was attributed to purely elastic behaviour, this flow surge is due to the tube’s viscosity. Increasing $\alpha$ increases the strength of this effect, as a faster increasing pressure allows a higher pressure to build before the instantaneous viscous response gives way to motion. The bulk inflation of the tube therefore increases with $\alpha$ , causing the observed gradual increase in maximum inlet flow rate.

Here, $\Delta Q$ largely follows the same trends as $Q_{max}$ , as seen in figure 12(d). The main difference is that at higher $De$ , there is no large initial decrease with increasing $\alpha$ . This implies that, during deflation, the viscous wall forces do not counter the squeezing effect entirely, in contrast to the flow surge mitigation during inflation. This is because the flow surge is resisted by both the elastic and viscous wall forces, whereas the squeezing is driven by the elastic forces.

Figure 13. Average flow rate through tube over one period of oscillation. The orange box represents a quasi-steady inflation cycle ( $\alpha \rightarrow 0$ ). The legend is in figure 12(b).

In some applications, the average flow rate is of interest. Here, we introduce the parameter $q_\infty$ , which is the average flow rate as $\alpha \rightarrow \infty$ ; it is the flow rate through the deforming tube with the steady applied pressure of $\tilde {P}_{in}=0.5$ . In figure 13, the average flow rate normalised by $q_\infty$ is plotted against $\alpha$ for various values of $De$ . The normalised value obtained from our simulations tends towards $1$ as $\alpha$ increases for all values of $De$ . Additionally, the value at $\alpha =0$ is also shown, which is the average flow rate of the quasi-steady inflation cycle. It can also be seen that the curves for all $De$ approach this value as $\alpha$ becomes small. The figure further supports that the effect of viscoelasticity on the normalised flow rate is most pronounced at moderate values of oscillation frequency.

Figure 14. Direction of flow rate hysteresis due to the interaction of viscoelastic response time and enforced oscillation. The shapes indicate the direction that the hysteresis curve is drawn out in time, with circles, deltas and squares being clockwise, cross-over and counterclockwise, respectively.

The interaction of viscoelastic response time and enforced oscillation affects the direction of the flow rate hysteresis curves, as shown in figure 14, with circles as clockwise, squares as counterclockwise and deltas as cross-over. In most cases, the hysteresis is clockwise at low $\alpha$ . When $De\le 5$ , upon increasing $\alpha$ , the curves eventually begin to cross-over over a small range of $\alpha$ near $\alpha \approx 0.8$ . This crossing over behaviour serves as a transition between the hysteresis directions, such that further increasing $\alpha$ leads to counterclockwise hysteresis. Cross-over indicates that $\alpha$ is high enough that the viscoelastic response lag of the tube is on par with the inherent phase lead of the flow. It therefore occurs at lower $\alpha$ when $De$ increases, because the increased tube viscosity increases the response time. This pattern continues until $De=5$ , where the cross-over point suddenly shifts to a higher $\alpha$ . At this $De$ , the flow surge effect from viscous bulk tube motion begins to occur, biasing the hysteresis clockwise (leading). With increasing $De$ , the effect is strong enough that there is no cross-over hysteresis for $De\ge 10$ in the range of $\alpha$ tested. Interestingly, at high enough $De$ , the hysteresis curve crosses over itself at very low $\alpha$ , changing to clockwise as $\alpha$ increases. These cases correspond to the maximum flow rate initially decreasing as $\alpha$ increases. This follows from the strong viscous forces countering the flow surge effect during inflation, causing a net lag during the inflation stage.

3.3.6. Effect of $G_1$

Figure 15. (a) $\Xi _{max}$ , $\Delta \Xi$ and deformation phase shifts plotted against $G_1/G_s$ . (b) $Q_{max}$ , $\Delta Q$ and flow rate phase shifts versus $G_1/G_s$ . Parameters: $De=50$ , $\alpha =0.56$ .

So far, we have fixed $G_1/G_s = 100$ , where $G_1$ is the shear modulus of the spring-dashpot branch of the SLS model. Next, we consider the effect $G_1$ . In figure 15, we show the behaviour of the tube for different values of $G_1/G_s$ . Given the spring represented by $G_1$ is in series with the dashpot representing the solid viscosity, it is as if the value of $G_1$ determines the strength of the contribution of viscosity to the instantaneous motion. Even with non-zero viscosity, the construction of the standard linear solid model allows for a jump change in deformation due to the presence of a spring in both of the parallel elements. When $G_1$ is large, the jump change is not possible, as the dashpot prevents such motion. For this reason, the effect of changing $G_1$ is similar to that of increasing $De$ . As $G_1$ increases, $Q_{max}$ , $\Delta Q$ , $\Xi _{max}$ and $\Delta \Xi$ decrease and approach some constant value (figure 15). Both $\phi _c$ and $\phi _t$ first increase at low $G_1$ and then begin to decrease. Initially, the increased elasticity amplifies the flow surge and squeezing effects, increasing the flow rate’s lead of the applied pressure. However, beyond a certain stiffness, the strength of the viscosity increases. This increases the tube response time, therefore causing a lag that decreases the phase shifts. It can be seen that this point occurs near $G_1/G_s=50$ , which is roughly where $Q_{max}$ and $\Delta Q$ begin to converge. As noted earlier, the SLS model approaches the Kelvin–Voigt model as $G_1\rightarrow \infty$ . This convergence behaviour is implied in figure 15.

4. Conclusion

Using a fully coupled 3-D FSI model, we present a detailed numerical study of pulsatile flow in a thin-walled viscoelastic microvessel at large deformation, which was lacking in the literature. Our primary focus was to study the effect of solid viscosity and oscillation frequency at small but non-vanishing inertia on the tube deformation, flow rate, phase shift and hysteresis.

We found that the general behaviour is dominated by an elastic flow surge during inflation and squeezing effects during deflation. At lower $De$ , the flow surge strengthens with increasing oscillation frequency and weakens with increasing viscosity. As such, increasing the oscillation frequency causes the maximum inlet flow rate to increase, whereas increasing $De$ causes it to decrease. This changes when the solid viscosity is strong enough to fully counter the flow surge at low $\alpha$ . Two notable limits are approached at both low and high $\alpha$ for all values of $De$ . Here, $\alpha \rightarrow 0$ represents the limit of a quasi-steady inflation cycle. Conversely, as $\alpha \rightarrow \infty$ , the behaviour approaches the steady-state inflation of a tube by the mean pressure $\tilde P_{in}=0.5$ . The tube motion changes with $De$ and $\alpha$ in accordance with these limits. When varying the oscillation frequency, both the maximum and amplitude of the maximum radius properly approach their corresponding quasi-steady inflation (low $\alpha$ ) and steady flow (high $\alpha$ ) values for all $De$ . The average flow rate similarly approaches the proper values in these limits. In both of these limits, all values of $De$ converge. Additionally, as $De\rightarrow \infty$ , one can expect the tube to behave as rigid in some deformed configuration due to an infinite tube response time. Our results suggest that deformation and flow rate are most affected in the intermediate range of $De$ and $\alpha$ .

The phase lead of the flow rate over the pressure noted by Womersley (1957b ) has also been expanded upon here and our model predicts that the phase shift is not constant but varies throughout the oscillation. Additionally, the flow rate does not uniformly lead the pressure, but lags under certain conditions. In fact, this variable phase shift reliably predicts the different flow rate hysteresis behaviours discussed by Lee & Schmid-Schönbein (1990). The direction in which the hysteresis proceeds is represented by the signs of the phase shift at the crest and trough ( $\phi _c$ and $\phi _t$ ); two positive signs indicate clockwise, two negative signs indicate counterclockwise and differing signs indicate cross-over. This behaviour is highly dependent on both the tube’s viscosity and the oscillation frequency. At lower $De$ , the hysteresis direction transitions from clockwise at low $\alpha$ to counterclockwise at high $\alpha$ , with a small range of cross-over in between. In contrast, at moderate $De$ , the hysteresis remains clockwise for all $\alpha$ tested. At very high $De$ , the hysteresis crosses over at very low $\alpha$ , while it is clockwise elsewhere. This directional change in hysteresis is fully characterised here in the $Q_{max}$ , $\alpha$ and $De$ parameter space.

Additionally, a distinct change in the flow physics was predicted at high values of $De$ . At high $De$ , the maximum inlet flow rate increases over a large range of $\alpha$ after an initial decrease, with such behaviour strengthening at higher $De$ . This is due to a change in the tube dynamics wherein the instantaneous response is akin to that of a rigid tube, allowing the pressure to build appreciably throughout its length. This leads to global or ‘whole-tube’ motion, as opposed to the elastic case in which the region near the inlet inflates somewhat independently. Another interpretation is that the more elastic tube does not allow the effect of higher frequency oscillations to penetrate far into the tube, while the more viscous tube does.

The main limitation of the current model is that it does not allow solid motion in the axial direction. Therefore, problems with large axial deformation, such as extremely large, balloon-like deformations or axial buckling, cannot be achieved. However, such deformations are unlikely to appear in the biological and microfluidic applications discussed. In the future, this model can be readily applied to study tubes with nonlinear constitutive models, higher inertial flows and flows with suspended capsules.