2.1. Fluid flow

We assume that the fluid flow is incompressible, axisymmetric and devoid of swirl. The flow velocity is denoted by $\boldsymbol {V}(R, Z, T) = \{V_R(R, Z, T), V_Z(R, Z, T)\}$ in cylindrical coordinates, and satisfies the incompressible Navier–Stokes equations

(2.1a)\begin{gather} \rho\left(\frac{\partial\boldsymbol{V}}{\partial T}+\boldsymbol{V}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{V}\right)={-}\boldsymbol{\nabla} P+\mu\,\nabla^2\boldsymbol{V}, \end{gather}

(2.1b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{V}=0 . \end{gather}

It is useful to cast the continuity equation in terms of the flux $Q(Z, T) = 2 {\rm \pi}\int _0^{R_t} V_Z R \,\mathrm {d} R$ by writing $\partial _Z Q + \partial _T ({\rm \pi} R_t^2) = 0$, which leads to

(2.2) \begin{equation} \frac{\partial^{}Q}{\partial Z} + 2 {\rm \pi}(R_0 + U)\,\frac{\partial^{}U}{\partial T} = 0. \end{equation}

The flow is instantaneously no-slip at the walls of the deforming tube ($\boldsymbol {V}= \partial _{T}{R_t} \boldsymbol {e}_R$ at $R = R_t$) and satisfies symmetry conditions at the centreline ($\partial _{R}{V_Z} = 0$, $V_R = 0$ at $R = 0$). The pressure oscillates as $P_i \cos (\omega T)$ at the inlet ($Z = 0$) and is zero at the outlet ($Z = L$).

2.2. Elastic response of the tube

To close the problem, it is necessary to relate the deformation of the tube to the stresses in the flow. For thin-walled tubes ($b \ll R_0$), resistance to deformation occurs primarily due to elastic stresses generated by circumferential stretching of the tube (Landau & Lifshitz Reference Landau and Lifshitz1986; Audoly & Pomeau Reference Audoly and Pomeau2018). Focusing on slender channels $(R_0 \ll L)$, and anticipating the analysis of the subsequent sections, we approximate the normal stress exerted by the fluid on the tube by the fluid pressure $P(Z,T)$. For small radial deformations $U \ll R_0$ and negligible axial deformations, these arguments lead to the pressure-displacement relation

(2.3)\begin{equation} U(Z,T) = \frac{R_0^2 (1 - \nu^2)}{E b}\,P(Z,T). \end{equation}

This approximation has been shown (Anand & Christov Reference Anand and Christov2020) to quantitatively capture 3-D DNS for the practically relevant case of nearly incompressible elastic solids ($\nu \approx 1/2$). More sophisticated membrane theories accounting for tension generated by viscous shear (Elbaz & Gat Reference Elbaz and Gat2014, Reference Elbaz and Gat2016; Rallabandi et al. Reference Rallabandi, Eggers, Herrada and Stone2021) have yielded results similar to (2.3) to leading order in the thickness of the tube walls. Models of the kind in (2.3), that relate local deformation to the local pressure (so-called Winkler foundation models), provide useful insight into flow–structure interactions in a wide range of systems (Dillard et al. Reference Dillard, Mukherjee, Karnal, Batra and Frechette2018). The bending rigidity and inertia of the thin-walled tubes are small and have been neglected in this analysis; their magnitudes are estimated in Appendix A.

2.3. Non-dimensionalization

We rescale the equations by introducing dimensionless spatial coordinates (lowercase) $z=Z/ L$ and $r=R/ R_0$, and a dimensionless time $t=\omega T$. A balance between pressure gradients and viscous stresses sets a characteristic axial velocity scale $P_{i} R_{0}^{2}/(\mu L)$. By continuity, the radial velocity scale is $P_{i}R_{0}^{3}/(\mu L^{2})$. We define a dimensionless pressure $p$, and dimensionless velocity components $v_z$ and $v_r$, according to

(2.4a–c)\begin{equation} p =\frac{P}{P_i}, \quad v_z = \frac{V_z}{P_i R_0^2/(\mu L)}\quad \mbox{and} \quad v_r=\frac{V_r}{P_i R_0^3/(\mu L^2)}. \end{equation}

The problem is governed by four independent dimensionless parameters,

(2.5a–d)\begin{equation} \varepsilon = \frac{R_0}{L}, \quad \mathcal{W}=\frac{ \rho R_0^2 \omega}{\mu}, \quad\varLambda = \frac{P_i R_0 (1 - \nu^2)}{E b} \quad \mbox{and} \quad \sigma = \frac{\mu L^2 \omega (1 - \nu^2)}{E b R_0}. \end{equation}

Here, $\varepsilon$ is the aspect ratio of the tube, and is small by construction. The Womersley number $\mathcal {W}$ is the ratio of a viscous diffusion time scale $T_d = \rho R_0^2/\mu$ to the driving time scale $\omega ^{-1}$. The elasticity of the tube walls enters in two distinct ways. The parameter $\varLambda$ characterizes the amplitude of radial deformation relative to the radius of the tube, and scales with the inlet pressure amplitude $P_i$. The linear elastic framework (2.3) implicitly assumes that $\varLambda$ is small, and we use this condition explicitly in the analysis of later sections. In addition to the deformation amplitude, the flow–structure interaction introduces a time scale $T_{EV} = \mu L^2 (1 - \nu^2)/(E b R_0)$ over which elastic deformations relax against the viscous resistance of the fluid (see e.g. Elbaz & Gat Reference Elbaz and Gat2014). The ratio of this relaxation time to the driving time defines the elasto-viscous parameter $\sigma$, which is analogous to the Deborah number in rheology. Similar quantities arise in other oscillatory fluid–elastic systems (Bickel Reference Bickel2007; Leroy & Charlaix Reference Leroy and Charlaix2011; Daddi-Moussa-Ider, Guckenberger & Gekle Reference Daddi-Moussa-Ider, Guckenberger and Gekle2016; Zhang et al. Reference Zhang, Arshad, Bertin, Almohamad, Raphael, Salez and Maali2022; Rallabandi Reference Rallabandi2024).

For typical parameters in arterial models and soft microfluidic systems, the aspect ratio $\varepsilon$ ranges from 0.001 to 0.1. The Womersley number $\mathcal {W}$ ranges from 1 to 100 at frequencies up to a few hundred Hz (Grotberg & Jensen Reference Grotberg and Jensen2004; Vishwanathan & Juarez Reference Vishwanathan and Juarez2020; Zhang & Rallabandi Reference Zhang and Rallabandi2023), while the small deformation amplitude $\varLambda$ falls within the range 0.01 to 0.1 (Kiran Raj et al. Reference Kiran Raj, DasGupta and Chakraborty2019; Pande et al. Reference Pande, Wang and Christov2023b). The elasto-viscous parameter $\sigma$ varies from 0.01 to 100, depending on the geometry of the tube ($\sigma$ is greater for longer tubes) and material properties; see Pande et al. (Reference Pande, Wang and Christov2023b). As we show in the following sections, the structure of the flow is governed primarily by $\mathcal {W}$ and $\sigma$, while $\varLambda$ sets the magnitude of the nonlinearity. We observe that $\mathcal {W}$ and $\sigma$ are similar in that they both represent the ratio of relaxation times to the driving time scale $\omega ^{-1}$. We will see later that $\mathcal {W}$ and $\sigma$ also control the spatial structure of the flow in similar ways.

3.1. Primary oscillatory flow

At $O(\varLambda ^0)$, advective inertia drops out of the momentum balance, yielding

(3.8)\begin{equation} \mathcal{W}\,\frac{\partial v_{1z}}{\partial t}={-}\frac{\partial p_1}{\partial z}+\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\, \frac{\partial v_{1z}}{\partial r}\right), \quad \text{with } v_{1z}|_{r = 1} = \left.\frac{\partial^{}v_{1z}}{\partial{r}}\right|_{r = 0} = 0. \end{equation}

The continuity equation (3.2) at $O(\varLambda ^0)$ is

(3.9) \begin{equation} \frac{\partial^{}q_1}{\partial{z}} + 2 {\rm \pi} \sigma\,\frac{\partial^{}p_1}{\partial{t}} = 0. \end{equation}

Motivated by boundary conditions (3.4) for pressure, we seek solutions in terms of complex phasors oscillating as ${\rm e}^{{\rm i} t}$, whose real parts solve the physical problem. Writing the pressure as

(3.10)\begin{equation} p_1(z, t) = \text{Re}[\mathcal{P}_1 (z)\, {\rm e}^{{\rm i} t}], \end{equation}

the oscillating axial velocity is found to be

(3.11)\begin{equation} v_{1z}= \frac{1}{\alpha^2}\,\frac{\partial \mathcal{P}_1}{\partial z} \left(\frac{I_0(r \alpha )}{ I_0(\alpha ) }- 1\right) {\rm e}^{{\rm i} t},\quad \text{with }\alpha=\sqrt{{\rm i}\, \mathcal{W}}, \end{equation}

where $I_n(\cdot )$ is a modified Bessel function of the first kind of order $n$. It will be understood that only the real parts of complex oscillating quantities are physically meaningful. Equation (3.11) is the well-known velocity profile identified by Womersley (Reference Womersley1955), which takes on a parabolic shape for small $\mathcal {W}$, and a plug-like profile with boundary layers of dimensional thickness $R_0 \mathcal {W}^{-1/2} = R_0 (\mu /(\rho \omega ))^{1/2}$ for large $\mathcal {W}$. From (3.6), the primary oscillatory flux is found to be

(3.12)\begin{equation} q_1=2{\rm \pi} \int^1_0 v_{1z} r \,\text{d}r={-}\frac{{\rm \pi}\,I_2(\alpha)}{\alpha^2\, I_0(\alpha)}\,\frac{\partial \mathcal{P}_1}{\partial z} \,\mathrm{e}^{\mathrm{i} t}. \end{equation}

Substituting (3.12) back into (3.9) leads to an equation governing the complex pressure:

(3.13)\begin{equation} \frac{\partial^2 \mathcal{P}_1}{\partial z^2}-k^2 \mathcal{P}_1=0,\quad \text{where } k=\sqrt{\mathrm{i}\sigma \, \frac{2\alpha^2\, I_0(\alpha) }{I_2(\alpha)}} \end{equation}

is a complex wavenumber. The above equation was also obtained by Ramachandra Rao (Reference Ramachandra Rao1983) in analysing linearized oscillatory flow in tapered elastic tubes, and a variant of it is given in Dragon & Grotberg (Reference Dragon and Grotberg1991). The primary pressure oscillates as $p_1= \cos t = \text {Re}[ {\rm e}^{{\rm i}t}]$ at the inlet $z =0$, and decreases to zero at the outlet ($\kern0.7pt p_1=0$ at $z=1$), yielding the boundary conditions $\mathcal {P}_1(0) = 1$, $\mathcal {P}_1(1)=0$. With these conditions, (3.13) admits the solution

(3.14)\begin{equation} \mathcal{P}_1(z) = \frac{\sinh\{k(1-z)\}}{\sinh{k}}. \end{equation}

Thus we see that the primary pressure in a deformable tube decays as a spatially oscillating exponential. The limit $\sigma \to 0$ recovers the rigid case, where the pressure decreases linearly along the tube.

The complex wavenumber $k$ scales as $\sigma ^{1/2}$ and depends on $\mathcal {W}$ through $\alpha$. Figure 2(a) shows the variation of $k/\sqrt {\sigma }$ with $\mathcal {W}$, indicating both real and imaginary parts. In the viscous limit ($\mathcal {W} \ll 1$), $k$ approaches a constant $4 \sqrt {\mathrm {i} \sigma }$, indicating comparable rates of spatial decay and oscillation. For $\mathcal {W} \gg 1$, the real part of $k$ approaches $\sqrt {\sigma }$, while the imaginary part diverges as ${2 \sigma \mathcal {W}}$. Asymptotic behaviours for both small and large $\mathcal {W}$ are indicated in figure 2(a). Figure 2(b) shows both real and imaginary parts of $\mathcal {P}_1$ for different $\sigma$ at a fixed $\mathcal {W}$. For short elastic relaxation times ($T_{EV} \ll \omega ^{-1}$, or small $\sigma$) the elastic solid responds instantaneously to the inlet pressure, and the pressure drops linearly in the tube, similar to a rigid channel. For longer relaxation times ($T_{EV} \gtrsim \omega ^{-1}$, or $\sigma \gtrsim 1$), local deformation of the solid lags the inlet pressure, leading to a non-monotonic pressure distribution in the channel with comparable real (in-phase relative to the inlet) and imaginary (out-of-phase relative to the inlet) contributions (see figure 2(b) for $\sigma = 1$). At large $\sigma$ (and small $\mathcal {W}$), the pressure distribution decays axially from the inlet over a length scale $\ell _{EV} = L \sigma ^{-1/2} \propto (E b R_0 /(\mu \omega ))^{1/2}$, with spatial oscillations. We observe that this is analogous to how the Womersley velocity profile exhibits boundary layers of thickness $R_0 \mathcal {W}^{-1/2}$ for large $\mathcal {W}$. Thus $\sigma$ controls the axial structure of the flow similarly to how $\mathcal {W}$ controls its radial structure. At larger values, $\mathcal {W}$ influences both radial and axial flow features; see (3.13), (3.14) and figure 2. Figure 2(c) depicts profiles of $p_1(z,t)$ at different times over an oscillation cycle for the particular combination $\sigma = \mathcal {W} = 1$. The pressure distribution travels along the channel in a wave-like pattern, gradually diminishing further away from the inlet of the tube. It is interesting to note that the elasticity of the tube strongly controls the oscillatory flow through a competition between relaxation and driving time scales, even though the amplitude of deformation has not yet entered at this order in the analysis.

Figure 2.

(a) Complex wavenumber $k$ as a function of $\mathcal {W}$, showing real and imaginary parts. Asymptotic expansions for large and small $\mathcal {W}$ are indicated by dashed curves. (b) Complex primary pressure $\mathcal {P}_1(z)$ at different $\sigma$ values, indicating real (solid) and imaginary (dashed) parts. The pressure decays linearly along the tube for small $\sigma$, whereas it decays exponentially for larger $\sigma$. (c) Pressure distribution at different times during an oscillation cycle for $\sigma =\mathcal {W}=1$, indicating wave-like behaviour.

3.2. Secondary steady flow

With insight into the oscillatory flow structure, we study the flow at $O(\varLambda )$, associated with finite deformation of the channel walls. We focus specifically on the rectified (time-averaged) part of this flow, which is associated with a net flux through the channel. We define the time average of a function over an oscillation cycle by $\langle\, f \rangle =(2 {\rm \pi})^{-1}\int _{\tau }^{\tau +2{\rm \pi} }f(t) \,\mathrm {d}t$. Averaging the momentum equation (3.1) and the no-slip condition (3.7), we find at $O(\varLambda )$ that

(3.15a)\begin{gather} \frac{\mathcal{W}}{\sigma} \langle \boldsymbol{v}_{1}\boldsymbol{\cdot} \boldsymbol{\nabla} v_{1z} \rangle={-}\frac{\partial \langle p_2 \rangle}{\partial z}+\frac{1}{r}\,\frac{\partial }{\partial r}\left(r\,\frac{\partial \langle v_{2z} \rangle}{\partial r}\right), \quad \text{with } \end{gather}

(3.15b)\begin{gather}\langle v_{2z} \rangle = v_{slip}(z) \stackrel{\text{def.}}{=} - \left.\left\langle p_1\, \frac{\partial^{}v_{1z}}{\partial{r}} \right\rangle\right|_{r =1} \quad \text{at } r = 1, \end{gather}

(3.15c)\begin{gather}\frac{\partial^{}}{\partial{r}} \langle v_{2z} \rangle = 0 \quad \text{at } r = 0. \end{gather}

Here, $v_{slip}(z)$ is an effective slip velocity satisfied by the secondary flow at the undeformed location of the tube walls, which arises as a consequence of the domain perturbation in (3.7). A useful rule for computing the time average of a product of oscillating complex quantities is $\langle \text {Re} [ A\, \mathrm {e}^{\mathrm {i} t}] \,\text {Re} [ B\, \mathrm {e}^{\mathrm {i} t}]\rangle = (1/2) \,\text {Re}[A^* B] = (1/2) \,\text {Re}[A B^*]$, where the asterisk denotes the complex conjugate (Longuet-Higgins Reference Longuet-Higgins1998). We apply this rule to the definition (3.15b) to calculate the effective slip velocity

(3.16)\begin{equation} v_{slip}(z) = {\rm Re} \left[\frac{k\, I_1(\alpha)}{2\alpha\,I_0(\alpha)\,|{\sinh{k}}|^2} \cosh\{k (1-z)\}\sinh\{k^*(1-z)\}\right].\end{equation}

Note that only the real parts of secondary (time-averaged) flow quantities written with complex variables hold physical significance. This will be implied from this point onwards unless otherwise indicated.

The problem (3.15) is linear in $\langle v_{2z} \rangle$ and can be decomposed into three parts: (i) the velocity generated from the advective ‘body force’, which we denote $v_{2z}^{adv}$; (ii) the velocity associated with the geometric nonlinearity due to deformation, which manifests as the effective slip at $r =1$; and (iii) the velocity due to the secondary pressure gradient. Solving (3.15) shows that the time-averaged axial velocity has the general structure

(3.17)\begin{equation} \langle v_{2z} \rangle=\frac{1}{4}\,\frac{\partial \langle p_2 \rangle}{\partial z}\,(r^2-1) + v_{slip} + v_{2z}^{adv}, \end{equation}

where the advective contribution $v_{2z}^{adv}$ satisfies

(3.18)\begin{align} \frac{\mathcal{W}}{\sigma}\,\langle \boldsymbol{v}_{1}\boldsymbol{\cdot} \boldsymbol{\nabla} v_{1z} \rangle=\frac{1}{r}\,\frac{\partial }{\partial r}\left(r\,\frac{\partial v_{2z}^{adv}}{\partial r}\right), \quad \text{with } v_{2z}^{adv}|_{r = 1} = 0 \quad \mbox{and} \quad \left.\frac{\partial^{}v_{2z}^{adv}}{\partial{r}} \right|_{r = 0} = 0. \end{align}

The system (3.18) lacks a simple analytical solution due to the complicated radial dependence (involving products of Bessel functions) of the advective body force, so we solve it numerically using finite difference techniques to find $v_{2z}^{adv}$.

To obtain the secondary pressure, we first insert (3.17) into (3.6), which yields a general expression for the mean flux:

(3.19)\begin{equation} \langle q_2 \rangle = 2{\rm \pi} \left(-\frac{1}{16}\,\frac{\langle p_2 \rangle}{\partial z}+\frac{1}{2}\,v_{slip} + \int_{0}^{1}v_{2z}^{adv}(r, z)\,r \,\mathrm{d}r\right). \end{equation}

We then average the continuity equation (3.2) and collect terms at $O(\varLambda )$ to find that the time-averaged flux satisfies $\partial _{z} {\langle q_2 \rangle }+ 2 {\rm \pi}\sigma \langle \partial _{t}{p_2} + p_1\,\partial _{t}{p_1} \rangle = 0$. Observing that $\langle \partial _{t}{p_2} \rangle = 0$ and $\langle p_1\,\partial _{t}{p_1} \rangle = \langle \partial _{t}(p_1^2/2) \rangle = 0$, we find, unsurprisingly, that $\langle q_2 \rangle$ is a constant. We use this result to integrate (3.19) along $z$, noting that $\langle p_2 \rangle$ must vanish exactly at both ends of the tube, to obtain

(3.20)\begin{equation} \langle q_2 \rangle = 2 {\rm \pi}\int_0^1 \left(\frac{1}{2}\,v_{slip}(z) + \int_{0}^{1}v_{2z}^{adv}(r, z)\,r \,\mathrm{d}r\right) \mathrm{d} z. \end{equation}

Rearranging (3.19), the time-averaged pressure is therefore

(3.21)\begin{equation} \langle p_2 \rangle(z) = 16 \int_0^z \left(\frac{1}{2}\, v_{slip}(x) + \int_{0}^{1}v_{2z}^{adv}(r, x)\,r\, \mathrm{d}r\right)\,{{\rm d}\kern0.7pt x} - \frac{8 z}{\rm \pi}\,\langle q_2 \rangle. \end{equation}

Using (3.16) and the numerical solution to (3.18), we implement the integrals in (3.20) and (3.21) numerically to obtain $\langle p_2 \rangle (z)$.

The axial distribution of the secondary pressure is shown in figure 3 for different $\sigma$ and $\mathcal {W}$. Symbols represent the numerical solutions discussed above, and solid curves represent an analytic theory that we detail in § 4. As $\sigma$ approaches zero for a fixed $\mathcal {W}$ (representing a nearly rigid channel with a short elasto-viscous relaxation time), the secondary pressure distribution becomes parabolic. However, as $\sigma$ increases (slower elasto-viscous relaxation), the secondary pressure distribution becomes asymmetric, with the peak pressure rising and occurring closer to the channel inlet. This symmetry breaking arises because the primary flow – which is responsible for driving the secondary flow – now decays exponentially from the inlet over the length scale $L/k$ (figure 2b). Figure 3(b) shows the dependence on $\mathcal {W}$ at a fixed value of $\sigma = 1$ (moderately flexible channels). As $\mathcal {W}$ increases, the peak of the secondary pressure increases up to $\mathcal {W} \approx 5$. For larger $\mathcal {W}$, the pressure begins to drop, deviating from the established pattern. Thus the pressure distribution depends sensitively on both $\sigma$ and $\mathcal {W}$. Advective inertial effects become noticeable for $\mathcal {W}$ as small as $0.5$, and dominate for larger $\mathcal {W}$ (see figure 6 in Appendix B). These features highlight the interplay between channel elasticity, fluid inertia and viscous resistance to flow, emphasizing the importance of both $\sigma$ and $\mathcal {W}$.

Figure 3.

Secondary flow pressure distribution for (a) fixed Womersley number $\mathcal {W}$ and (b) fixed elasto-viscous parameter $\sigma$. Symbols are results of numerical calculations, and curves represent the approximate analytic theory. Plot (a) shows the transition from parabolic distribution in quasi-rigid channels ($\sigma \to 0$) to asymmetric patterns in flexible channels ($\sigma$ increases). Plot (b) demonstrates increasing peak pressure with rising $\mathcal {W}$ in flexible channels for a range of $\mathcal {W}$, though the dependence becomes non-monotonic for $\mathcal {W} \gtrsim 5$.

As discussed previously, the flow varies axially over length scales set by $L$ and $L/k$, while radial gradients occur on scales $R_0$ and $R_0 \mathcal {W}^{-1/2}$ (where the latter quantity is the thickness of the oscillatory Stokes layer). The long-wave approximation used in the theory requires axial velocity gradients to be much smaller than radial ones, leading to the conditions $\varepsilon \ll 1$ and $\sigma \ll \varepsilon ^{-2}$. Furthermore, the theory requires small deformation amplitudes relevant to the relevant radial length scales of the flow, which yields the conditions $\varLambda \ll 1$ and $\mathcal {W} \ll \varLambda ^{-2}$. This condition on $\mathcal {W}$ also justifies neglecting the advective inertia of the secondary flow in (3.15) (Riley Reference Riley2001).