2. Problem formulation

In this section, we present the problem definition, along with the physical parameters relevant to the analysis and the small non-dimensional parameters. The examined liquid-filled chamber is illustrated in figure 1. Here, a spherical geometry is assumed (the validity of this assumption will be verified by numerical simulation in § 5.2). A spherical hyperelastic chamber with a stress-free radius of $r_{0}$ is connected to one or two rigid tubes with radius $a$ and length $\ell$. For simplicity, we assume identical tubes in the inlet and the outlet. Here, the flow field inside the chamber and tubes is considered incompressible, Newtonian, and with negligible inertial effects. The fluid's axial velocity inside the tube is $u_{z}$, and the volumetric flux rate is denoted by $q(t)$ (where $q_{in}(t)$ refers to the flow entering the body from the inlet tube and $q_{out}(t)$ refers to the flow moving from the body through the outlet tube). The relevant variables and parameters are: the time $t$; the axial coordinate and symmetry axis $z$; and the radial coordinate $s$ of the cylindrical system used to describe the tubes. Axisymmetry allows us to eliminate the azimuthal angle of the cylindrical system. Furthermore, the pressure and flow velocity fields of the entrapped fluid are $p(t,s,z)$ and $\boldsymbol {v}(t,s,z)$, while its constant density and dynamic viscosity are denoted by $\rho$ and $\mu$. The chamber's dynamics are approximated by a single degree of freedom, represented here by the chamber's instantaneous radius, denoted by $\eta (t)$. For spherical geometry, a coordinate system is chosen so that one of the coordinates remains constant on the boundary. Here, $\{r,\theta,\phi \}$ are the coordinates of a moving spherical system, located at the centre of the chamber, where $\theta$ is the polar angle, measured from the axis of symmetry to the radial coordinate $r$, and $\phi$ is the azimuthal angle, revolving around the axis of symmetry, $z$. The Cauchy stress tensor of the flow is denoted as $\boldsymbol {\sigma }_{\boldsymbol {f}}$. The stress-free shell's thickness is $w_{0}$, which is considered to be much smaller than the stress-free chamber's radius, namely $w_{0}\ll r_{0}$.

Figure 1. Illustration of the system under investigation, consisting of a hyperelastic chamber, as well as a fixed inlet tube and a massless outlet tube, both having identical radius and length. The bottom of the chamber is held fixed during the inflation, while its centre is allowed to move.

The following analysis utilizes three small parameters, including the ratio between the radius of the tubes and the radius of the stress-free chamber ($\epsilon_{a}$, denoted hereafter as the tube–chamber radii ratio),

(2.1)\begin{equation} \epsilon_{a} = \frac{a}{r_{0}}\ll 1, \end{equation}

the slenderness of the tubes ($\epsilon_{t}$, denoted hereafter as the tube slenderness),

(2.2)\begin{equation} \epsilon_{t} = \frac{a}{\ell}\ll 1 \end{equation}

and the last small parameter in the analysis is taken as the ratio between the viscous stresses and the overall pressure in the chamber ($\varepsilon$, denoted hereafter as the chamber viscous resistance parameter), defined by

(2.3)\begin{equation} \varepsilon = \frac{\mu v^{*}}{r_0 p^{*}}\ll 1, \end{equation}

where $v^{*}$ is the characteristic flow velocity in the chamber and $p^{*}$ is the characteristic pressure of the system. For the following analysis, we shall normalize the physical variables by considering the characteristic values of the problem as follows:

(2.4a–g)\begin{align} \boldsymbol{V}=\frac{\boldsymbol{v}}{v^{*}},\quad \widehat{\boldsymbol{\nabla}}=r_{0}\boldsymbol{\nabla},\quad \hat{\boldsymbol{\sigma}}_{f}=\frac{\boldsymbol{\sigma_{f}}}{p^{*}},\quad R=\frac{r}{r_{0}},\quad \lambda=\frac{\eta}{r_{0}},\quad P=\frac{p}{p^{*}},\quad T=\frac{t}{r_{0}/v^{*}}, \end{align}

where $\lambda (T)$ denotes the stretch of the chamber and $R$ is the normalized radial coordinate.

3. Constitutive model for a hyperelastic membrane

This section presents the constitutive law that governs the spherical shell dynamics. We consider a thin-shelled, spherical chamber made of incompressible hyperelastic isotropic material. Finite elasticity theory dictates a known form of the elastic strain energy density $\psi (\lambda )$, which depends only on the relative stretch $\lambda (T)$. Moreover, the elastic strain energy density satisfies $\psi (1)=0$. Different types of hyperelastic models differ in the type of material and the elastic strains experienced without failing. The most common models are neo-Hookean, Mooney–Rivlin, Ogden, Gent and Biological tissue (Ogden Reference Ogden1972).

The material is assumed to be incompressible, which leads to the relation between the pressurized and the stress-free states, given by $r^{2}w\approx r_{0}^{2}w_{0}$. Thanks to this relation, the chamber's instantaneous thickness is eliminated. To capture the chamber's bistability, we use the two-parameter Mooney–Rivlin model.

Under the above assumptions, and considering incompressibility, the normalized solid's Mooney–Rivlin strain energy function is given by (Ogden Reference Ogden1972; Beatty Reference Beatty1987)

(3.1)\begin{equation} \hat{\psi}(\lambda)=2\lambda^{2}+\frac{1}{\lambda^{4}}-3+\alpha\left(\lambda^{4}+\frac{2}{\lambda^{2}}-3\right), \end{equation}

where $\alpha =s_{2}/s_{1}$ is the ratio between two empirically determined constants, commonly denoted as the Mooney–Rivlin parameters.

In this study, the Mooney–Rivlin parameters are chosen as $s_{1}\approx 1.5$ MPa and $s_{2}\approx 0.15$ MPa (Beatty Reference Beatty1987; Treloar Reference Treloar1975). The normalized (3.1) was obtained by using the normalization $\psi ^{*}=s_{1}$ for the strain energy density function, and the magnitude of the parameter $\alpha$ is $O(10^{-1})$.

We first study the chamber's static behaviour, where the pressure (without flow) is dictated. In this case, both the stretch and the pressure are constant, denoted here by $\lambda _{SS}$ and $P_{SS}$, respectively. The behaviour mentioned above can be demonstrated by the overall effective potential energy of the system,

(3.2)\begin{equation} \mathcal{U}(\lambda_{SS};P_{SS})=\int_{\lambda_{SS}}\left(\frac{{\rm d}\hat{\psi}}{{\rm d}\xi}-\xi^{2}P_{SS}\right)\,\mathrm{d}\xi=\hat{\psi}(\lambda_{SS})-\frac{1}{3}\lambda_{SS}^{3}P_{SS}. \end{equation}

Based on (3.2), figure 2(a) shows a curve of the potential energy function where the constant pressure is $P_{SS}=2.75$. The Mooney–Rivlin relation (3.1), along with the steady version of the leading-order energy balance (3.2), formulated as $\partial \mathcal {U}/\partial \lambda _{SS}=0$ at constant pressure $P_{SS}$, yields a relation between stretch, $\lambda _{SS}$, and pressure, $P_{SS},$ in equilibrium condition

(3.3)\begin{equation} P_{SS}=\left.\left(\frac{1}{\lambda^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right)\right|_{\lambda_{SS}}=4 \left[\frac{1}{\lambda_{SS}}-\frac{1}{\lambda^{7}_{SS}}+\alpha\left(\lambda_{SS} -\frac{1}{\lambda^{5}_{SS}}\right)\right];\quad 0<\alpha\ll 1. \end{equation}

Figure 2. (a) The solid blue curve is a characteristic stretch–pressure curve of a single elastic chamber with $\alpha =0.1$. The solid orange curve is the effective potential energy function (3.2), corresponding to the constant pressure $P_{SS}=2.75$, illustrated by the dashed black line; green and red dots are the stable and unstable equilibrium radii. (b) Characteristic normalized stretch–pressure curve of a single elastic chamber according to (3.3) with $\alpha =0.1$. The extrema are marked, and the bistable region is marked in grey. Solid curves are stable branches and dashed curves are unstable ones. (c) The solid blue line is the exact solution of (3.3) calculated numerically. The dashed curves are the approximated solutions obtained in (A8)–(A9). (d) The evolution of the extremum values of pressure $P_{A},P_{B}$ and stretch $\lambda _{A},\lambda _{B}$ as a function of the small parameter $\alpha$. The solid curves are the exact values, and the dashed curves are the asymptotic approximations given in Appendix A ((A3), (A4) and (A7)).

This well known relation was extensively leveraged to describe the quasi-static inflation of spherical balloons (Treloar Reference Treloar1975; Beatty Reference Beatty1987; Ben-Haim et al. Reference Ben-Haim, Salem, Or and Gat2020) for spatially uniform pressures. As seen from figure 2(a), figure 2(a,b) showing the relation in (3.3) with $\alpha =0.1$, the uniform pressure of the chamber is not monotonic with respect to the radius. Therefore, the inverse relation, describing the chamber's radius as a function of the pressure, cannot be directly extracted. The curve $P_{SS}(\lambda _{SS})$ in figure 2(b) has two bifurcation points, described by a local maximum point at $(\lambda _{A},P_{A})$, and a local minimum point at $(\lambda _{B},P_{B})$. This figure shows a bifurcation, which occurs when the pressure enters or exits the range between the local extrema, $P_{B}< P_{SS}< P_{A}$, illustrated in grey. Asymptotic approximations for the bifurcation points of the equilibrium curve, $P_{A},P_{B},\lambda _{A}$ and $\lambda _{B}$, appear in Appendix A. The evolution of those extrema as a function of the small parameter $\alpha$ is presented in figure 2(d). Asymptotic approximations for the solution of the equilibrium equation (3.3), appear in Appendix A. Those approximations are plotted in figure 2(c) with dashed lines on the solid exact solution curves, represented by the inverse relation $\lambda _{SS}(P_{SS})$.

Analysing the effective potential energy $\mathcal {U}(\lambda _{SS};P_{ss})$ in (3.2), using the second derivative with respect to $\lambda _{SS}$, it can be proved that the right and left branches of $P_{SS}(\lambda _{SS})$ where $1<\lambda <\lambda _{A}$ or $\lambda _{B}<\lambda$ are stable equilibria and satisfy

(3.4)\begin{equation} \left.\frac{\partial^{2}\mathcal{U}}{\partial\lambda_{SS}^{2}}\right|_{P_{SS}}=\frac{{\rm d}P_{SS}}{{\rm d}\lambda_{SS}}>0. \end{equation}

Conversely, the intermediate branch $\lambda _{A}<\lambda _{SS}<\lambda _{B}$ is an unstable region satisfying $\partial ^{2}\mathcal {U}/\partial \lambda _{SS}^{2}<0$. This is precisely the bistability phenomenon.

4. Series solution of the flow field within an expanding chamber

In this section, the governing equations of the flow within the chamber will be formulated, as well as the problem's boundary conditions. An analytical series solution will then be presented, describing the velocity field and the flow's pressure distribution inside the spherical chamber.

4.1. Formulation and analysis of the governing equations

Under the assumptions discussed above, the momentum and continuity equations governing the fluid's behaviour expressed in the moving spherical frame,

(4.1a)$$\begin{gather} \rho\left(\frac{\partial{\boldsymbol{v}}}{\partial{t}}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}+\hat{\boldsymbol{z}}\frac{{\rm d}^{2}}{{\rm d}t^{2}}\sqrt{\eta^{2}-a^{2}}\right) ={-}\boldsymbol{\nabla} p+\mu\nabla^{2}\boldsymbol{v}-\rho g\hat{\boldsymbol{z}} , \end{gather}$$

(4.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} = 0 , \end{gather}$$

where the third term in the left-hand expression in the momentum equation (4.1a) describes the acceleration of the moving spherical frame, centred on the chamber's moving centre, relative to a stationary frame.

Assuming the flow in tubes is fully developed and axisymmetric, the volumetric flux is given by

(4.2)\begin{equation} q(t)={-}\frac{{\rm \pi} a^{4}}{8\mu} \frac{\partial p_{t}}{\partial{z}}, \end{equation}

where $\partial p_{t}/\partial z$ is the pressure gradient along the tube. Since the tube slenderness $\epsilon _{t}\ll 1$ we shall assume a constant pressure gradient. Normalization of (4.2) yields the characteristic flow rate as $q^{*}={\rm \pi} a^{2} u_z^{*}$ where $u_{z}^{*}=a^{2} p^{*}/\mu \ell$ is the characteristic axial component of the fluid velocity in the tube. An integral flow balance yields the relation between the characteristic velocity in the tube and the characteristic velocity of the flow within the chamber, as $v^{*}=\epsilon _{a}^{2} u_{z}^{*}$. Substituting the characteristic values into the chamber viscous resistance parameter (2.3), relates it to the other small parameters, as follows:

(4.3)\begin{equation} \varepsilon = \frac{a^{4}}{r_{0}^{3} \ell}=\epsilon^{3}_{a}\epsilon_{t}\ll 1. \end{equation}

From relation (4.3) it is clear that $\varepsilon$ is dependent merely on the geometry of the system, providing a simple relation between the hydrostatic and deviatoric stresses. Hence, an appropriate geometry can be defined in order to design an efficient and controllable system. We consider negligible gravity, i.e. $\rho gr_{0}/p^{*}\ll 1$ (where $g$ is the gravitational acceleration), and define a Reynolds number in the chamber as $Re=\rho v^{*}r_{0}/\mu \ll 1$. Since the Reynolds number is small, the flow's inertia may be neglected. Therefore, by utilizing the non-dimensional quantities specified in (2.4a–g), the fluid's motion (4.1) is governed by the Stokes equations for creeping flow with an implicit time variable,

(4.4a,b)\begin{equation} \widehat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{V}=0, \quad \widehat{\boldsymbol{\nabla}} P=\varepsilon \widehat{\nabla}^{2}\boldsymbol{V}+O(\varepsilon Re). \end{equation}

The validity of these equations is weakened at the vicinity of the connections to the tubes since in those regions, the characteristic velocity is approximately $u_{z}^{*}$ rather than $v^{*}$. At a distance $O(L)$ from the tube (where $a< L< r_{0}$), the velocity scale is $q^{*}/L^{2}$ and hence the Reynolds number is $\rho q^{*}/\mu L<\rho q^{*}/\mu a$. Thus, a sufficient condition for global neglect of inertia is

(4.5)\begin{equation} q^{*}\ll \frac{{\rm \pi}\mu a}{\rho}. \end{equation}

From the non-dimensional relation (4.4a,b), in the sphere, the pressure is spatially uniform at leading order and the viscous flow will generate small spatially varying corrections. Moreover, in order to get a better understanding of the chamber resistance small parameter's physical meaning, we may use the dynamical stress tensor in the fluid domain defined by the constitutive relation $\boldsymbol {\sigma _{f}}=-p\boldsymbol {I}+\mu [\boldsymbol {\nabla } \boldsymbol {v}+(\boldsymbol {\nabla } \boldsymbol {v})^{{\rm T}}]$ where $\boldsymbol {I}$ is the $3\times 3$ unit matrix. In the most general constitutive equation, $\boldsymbol {\sigma }_{f}$ consists of the linear and instantaneous dependence of the deviatoric stress, plus the hydrostatic stress, $-p\boldsymbol {I}$, stemming from the static pressure. Normalization of the total stress tensor yields

(4.6)\begin{equation} \boldsymbol{\hat{\sigma}_{f}}={-}P\boldsymbol{I}+\varepsilon[\boldsymbol{\widehat{\nabla}V}+\left(\boldsymbol{\widehat{\nabla} V}\right)^{{\rm T}}]. \end{equation}

As one can notice from (4.6) the velocity field is not included at the leading order. Mainly, the leading order of the problem is a case of fully developed uniform pressure without any velocities. Suppose the flow is dictated by controlled pressure or flux at the inlet, the velocity is generated, and additional small deviatoric stress is created, which quasi-statically leads the system to another hydrostatic state.

The axisymmetric Stokes equations (4.4a,b) can be solved in spherical polar coordinates using a series expansion (Happel & Brenner Reference Happel and Brenner2012). Consider the non-dimensional Stokes stream function $\varPsi (R,\theta,T)$. The flow velocity components $V_{R}$ and $V_{\theta }$ are related to the Stokes stream function $\varPsi (R,\theta,T)$ through

(4.7a,b)\begin{equation} V_{R}=\frac{1}{R^{2}\sin\theta} \frac{\partial\varPsi}{\partial\theta};\quad V_{\theta}={-}\frac{1}{R\sin\theta} \frac{\partial\varPsi}{\partial R}, \end{equation}

where $V_{R}$ and $V_{\theta }$ are the radial and tangential velocity components, respectively. By applying the curl operator to the momentum equation (4.4a,b) and using several simple algebraic manipulations, the Stokes equation can be reduced to a fourth-order biharmonic equation obtained in terms of the Stokes stream function as follows:

(4.8a)$$\begin{gather} E^{2}(E^{2}\varPsi) = 0 , \end{gather}$$

(4.8b)$$\begin{gather}\boldsymbol{\nabla} P ={-}\varepsilon\frac{\hat{\boldsymbol{\phi}}}{R\sin\theta}\times\boldsymbol{\nabla}(E^{2}\varPsi), \end{gather}$$

where in spherical coordinates,

(4.9)\begin{equation} E^{2}=\frac{\partial^{2}}{\partial{R^{2}}}+\frac{\sin\theta}{R^{2}} \frac{\partial}{\partial{\theta}} \left(\frac{1}{\sin\theta} \frac{\partial}{\partial{\theta}}\right). \end{equation}

Equation (4.8a) is solved by separation of variables, and (4.8b) is solved by integration with respect to the radial and tangential directions. However, for brevity we will not present the full calculation of the solution here. A solution for the stream function in spherical coordinates is of the form

(4.10)\begin{equation} \varPsi(R,\theta;T)=\mathcal{A}_{0}(T)+\sum_{n=2}^{\infty}[\mathcal{A}_{n}(T)R^{n}+\mathcal{C}_{n}(T)R^{n+2}]\mathcal{J}_{n}(\cos\theta) ,\end{equation}

where $\mathcal {A}_{n}(T)$ and $\mathcal {C}_{n}(T)$ are unknown functions, determined by the boundary conditions, and $\mathcal {J}_{n}(\xi )$ are the Gegenbauer functions of the first kind of order $n$ (and degree $-1/2$). Happel & Brenner (Reference Happel and Brenner2012) have exhaustively investigated the properties of these Gegenbauer functions in connection with the hydrodynamic application. For our present purposes, their properties can be deduced most readily from their relation with the corresponding Legendre functions of the first kind $\mathcal {P}_{n}(\xi )$ as

(4.11)\begin{equation} \mathcal{J}_{n}(\xi)=\frac{\mathcal{P}_{n-2}(\xi)-\mathcal{P}_{n}(\xi)}{2n-1}={-}\frac{1}{(n-1)!}\left(\frac{{\rm d}}{{\rm d}\xi}\right)^{n-2}\left(\frac{\xi^{2}-1}{2}\right)^{n-1};\quad n\geqslant 2. \end{equation}

In the degenerate cases $n=0,1$ we define $\mathcal {J}_{0}(\xi )=1$ and $\mathcal {J}_{1}(\xi )=-\xi$, respectively. Using the definition of the Stokes stream function (4.7a,b), we describe the series solution of the velocity field and pressure distribution as

(4.12a)$$\begin{gather} V_{R}(R,\theta;T) ={-}\sum_{n=2}^{\infty}\left[\mathcal{A}_{n}(T)R^{n-2}+\mathcal{C}_{n}(T)R^{n}\right]\mathcal{P}_{n-1}(\cos\theta), \end{gather}$$

(4.12b)$$\begin{gather}V_{\theta}(R,\theta;T) = \sum_{n=2}^{\infty}\left[n\mathcal{A}_{n}(T)R^{n-2}+(n+2)\mathcal{C}_{n}(T)R^{n}\right]\frac{\mathcal{J}_{n}(\cos\theta)}{\sin\theta} , \end{gather}$$

(4.12c)$$\begin{gather}P(R,\theta;T) = P_{0}(T)-\varepsilon\sum_{n=2}^{\infty}\left[\frac{2(2n+1)}{n-1}R^{n-1}\mathcal{C}_{n}(T)\right]\mathcal{P}_{n-1}(\cos\theta), \end{gather}$$

where $P_{0}(T)$ should be determined through the physical boundary conditions defined by pressure at the chamber's inlet and outlet. The unknown functions, $\mathcal {A}_{n}(T)$ and $\mathcal {C}_{n}(T)$, should be determined by requiring the kinematic boundary conditions of the flow to be satisfied. Note that the series solutions (4.12) are an exact solution of (4.8).

4.2. Formulation and discussion of the boundary conditions

Next, we are interested in finding the unknown functions, $\mathcal {A}_{n}(T)$ and $\mathcal {C}_{n}(T)$ by imposing two boundary conditions. The first boundary condition is obtained from the assumption that there is no penetration into the chamber's boundaries, and the second boundary condition is the no-slip condition. The Eulerian description of a material point located on the chamber's wall is given by $\boldsymbol {r} =\eta (t)\boldsymbol {\hat {r}}$, where $\boldsymbol {\hat {r}}=(\sin \theta \cos \phi,\sin \theta \sin \phi,\cos \theta )$ is the radial direction unit vector; thus, the material point's velocity is $\boldsymbol {\dot {\hat {r}}}=\dot {\eta }\boldsymbol {\hat {r}}+\eta \dot {\theta }\boldsymbol {\hat {\theta }}$, where $\boldsymbol {\hat {\theta }}=(\cos \theta \cos \phi,\cos \theta \sin \phi,-\sin \theta )$ is the polar direction unit vector. We recall that the fluid's behaviour is described in the moving spherical frame; therefore, the first term is a radial component, while the polar component is created due to the constraint of the rigid inlet/outlet tube. Assuming that the deformation is spherical, the following kinematic constraint must be satisfied:

(4.13)\begin{equation} \frac{\theta-\kappa\theta_i(t)}{\theta_f(t)-\theta}=\frac{\theta(0)-\kappa\theta_i(0)}{\theta_f(0)-\theta(0)}, \end{equation}

where $\kappa$ indicates whether there are both inlet and outlet tubes (when $\kappa =1$), or only an inlet tube (when $\kappa =0$). Moreover, $\theta _{i}(t)$ and $\theta _{f}(t)$ are angles corresponding to the connections between the tubes and the chamber (see figure 1). By simple geometric considerations, we get

(4.14a,b)\begin{equation} \theta_{i}(t)=\sin^{{-}1}\left(\frac{a}{\eta(t)}\right), \quad \theta_{f}(t)={\rm \pi}-\sin^{{-}1}\left(\frac{a}{\eta(t)}\right). \end{equation}

Considering the derivative of (4.13) with respect to time will lead to the relation between $\dot {\theta }$ and $\theta,\theta _i,\theta _f,\dot {\theta }_i$ and $\dot {\theta }_f$. Therefore, the material point's velocity can be rewritten as

(4.15)\begin{equation} \boldsymbol{\dot{\hat{r}}} =\boldsymbol{\hat{r}}\frac{{\rm d}\eta}{{\rm d}t}+\boldsymbol{\hat{\theta}}\frac{(\theta_{f}-\theta)\kappa\dot{\theta}_{i}+(\theta-\kappa\theta_{i})\dot{\theta}_{f}}{\theta_{f}-\kappa\theta_{i}}\eta(t). \end{equation}

The result obtained in (4.15) represents the change in the angle of a material point relative to the initial state. Consequently, the no-penetration and the no-slip conditions are defined in vector form as

(4.16) \begin{equation} \boldsymbol{v}\left(r=\eta(t),\theta;t\right) = \boldsymbol{\dot{\hat{r}}}+ \boldsymbol{\hat{z}}\left\{ \begin{array}{@{}ll} \kappa u_{z}^{(out)}(s_{out}=\eta\sin\theta) , & 0\leqslant\theta\leqslant\theta_{i}, \\ 0, & \theta_{i}<\theta<\theta_{f}, \\ u_{z}^{(in)}(s_{in}=\eta\sin\theta) , & \theta_{f}\leqslant\theta\leqslant{\rm \pi}. \end{array} \right. \end{equation}

Here, we assumed a spherical surface at the connection between the tubes and the chamber. The first vector in (4.16) is the spherical surface velocity that captures the inflation of the whole chamber, while the second component adds the fluid's velocity into or out from the tube. Note that the velocity of the chamber's centre (the velocity of the moving spherical frame, centred on the chamber's moving centre, relative to a stationary frame) should be subtracted from the fluid's velocity components that come into or out of the tube. However, the velocity of the moving spherical frame is negligible relative to the velocity of the flow in the tube. The wall's elasticity is expressed in the body's ability to increase the volume according to the material's constitutive laws (hyperelasticity). The asymptotic justification for neglecting the non-spherical deformation of the chamber appears in Appendix B.

A simple investigation of the boundary condition's derivative shows singularities at $\theta =\theta _{i}$ and $\theta =\theta _{f}$ where the known parabolic Hagen–Poiseuille relation is used to describe the flow inside the tube. Since the pressure distribution represented by the Legendre series (4.12) is the solution of the second-order differential equation, singularities lead to a divergent series. In order to avoid the singularities and to get a modified velocity profile that also takes into account the end effects, we assume a modified flow profile at the chamber's inlet (or outlet) as follows:

(4.17)\begin{equation} u_{z}^{(m)}(s,t;\epsilon_\omega) = u_{z}^{p}(s,t)+u_{z}^{c}(s,t;\epsilon_\omega). \end{equation}

Here, $u_{z}^{p}(s,t)$ is the known parabolic Hagen–Poiseuille profile

(4.18)\begin{equation} u_{z}^{p}(s,t) = \frac{2q(t)}{{\rm \pi} a^{2}}\left[1-\left(\frac{s}{a}\right)^{2}\right] \end{equation}

and $u_{z}^{c}(s,t;\epsilon _\omega )$ is an correction profile defined by

(4.19)\begin{align} u_{z}^{c}(s,t;\epsilon_\omega) = \left(\omega(s;\epsilon_\omega)-1\right)u_{z}^{p}(s,t)+\gamma\cdot\frac{q(t)}{{\rm \pi} a^{2}}\left[1-4\left(\frac{s}{a}\right)^{2}+3\left(\frac{s}{a}\right)^{4}\right]\omega(s;\epsilon_\omega), \end{align}

where $\gamma$ is a parameter determined by fitting the velocity profile obtained in a finite elements calculation. The function $\omega (s;\epsilon _\omega )$ is defined by $\omega (s;\epsilon _\omega )=-\tanh {(2/\epsilon _\omega )}+\tanh {((s/a+1)/\epsilon _\omega )}-\tanh {((s/a-1)/\epsilon _\omega )}$, and $0<\epsilon _\omega \ll 1$ is an arbitrary small parameter (denoted hereafter as the smoothing parameter). In fact, $\omega (s;\epsilon _\omega )$ is a weighted function making the velocity profile differentiable even at $\theta _i$ and $\theta _{f}$. Moreover, it can be shown that $\omega (s;\epsilon _\omega )=1+O(e^{-1/\epsilon _\omega })$ as $\epsilon _\omega \ll 1$ and $|s/a|\lesssim 1-O(\epsilon _\omega )$. The modified flow profile has four physical properties. The first is symmetry around the tube's radial coordinate, $s$. The second is the no-slip condition represented mathematically by zero velocity on the tube boundaries, $u_{z}^{(m)}(s=a,t)=0$. The third nulls the radial velocity gradient at the boundaries, $\partial u_{z}^{(m)}/\partial s=0$ where $s=a$. By neglecting $O(e^{-1/\epsilon _\omega })$ terms, the modified flow profile's fourth property is keeping the total flow equal to the Hagen–Poiseuille model's value (regardless of the choice of $\gamma$). Using non-dimensional parameters, the boundary conditions become

(4.20a)$$\begin{gather} V_{R}\left(R=\lambda\right)=\frac{{\rm d}\lambda}{{\rm d}T} + \frac{\cos\theta}{\lambda^{2}\tilde{\epsilon}^{2}} \left\{ \begin{array}{ll} \left.\kappa U_{z}^{(m)}\right|_{S=\tilde{\epsilon}^{{-}1}\sin\theta}, & 0\leqslant\theta\leqslant\sin^{{-}1}(\tilde{\epsilon}),\\ 0, & \sin^{{-}1}(\tilde{\epsilon})<\theta<{\rm \pi}-\sin^{{-}1}(\tilde{\epsilon}) \\ \left. U_{z}^{(m)}\right|_{S=\tilde{\epsilon}^{{-}1}\sin\theta}, & {\rm \pi}-\sin^{{-}1}(\tilde{\epsilon})\leqslant\theta\leqslant{\rm \pi}, \end{array} \right. \end{gather}$$

(4.20b)$$\begin{gather}V_{\theta}\left(R=\lambda\right)= \varGamma(\theta) -\frac{\sin\theta}{\lambda^{2}\tilde{\epsilon}^{2}} \left\{ \begin{array}{ll} \left.\kappa U_{z}^{(m)}\right|_{S=\tilde{\epsilon}^{{-}1}\sin\theta} , & 0\leqslant\theta\leqslant\sin^{{-}1}(\tilde{\epsilon}), \\ 0 , & \sin^{{-}1}(\tilde{\epsilon})<\theta<{\rm \pi}-\sin^{{-}1}(\tilde{\epsilon}), \\ \left. U_{z}^{(m)}\right|_{S=\tilde{\epsilon}^{{-}1}\sin\theta}, & {\rm \pi}-\sin^{{-}1}(\tilde{\epsilon})\leqslant\theta\leqslant{\rm \pi}, \end{array} \right. \end{gather}$$

where,

(4.21a,b)\begin{equation} \varGamma(\theta;\tilde{\epsilon})=\varTheta(\theta)\frac{\lambda\tilde{\epsilon}}{\sqrt{1-\tilde{\epsilon}^{2}}}\frac{{\rm d}\lambda}{{\rm d}T};\quad \varTheta(\theta)=\frac{(1+\kappa)\theta-\kappa(\theta_i+\theta_f)}{\theta_f-\kappa\theta_i}, \end{equation}

where $\boldsymbol {U_{z}^{(m)}=u_{z}^{(m)}}/u_{z}^{*}$ is the normalized modified axial velocity in tube, $S=s/a$ is the normalized cylindrical radial coordinate in tube and $\tilde {\epsilon }(T)=\epsilon _{a}/\lambda (T)\ll 1$.

The integral mass conservation equation is given by

(4.22)\begin{align} q_{in}-\kappa q_{out}=\frac{{\rm d}}{{\rm d}t}\left[\frac{(1+\kappa){\rm \pi} a^{2}}{3}\sqrt{\eta^{2}(t)-a^{2}}+\int_{\phi=0}^{2{\rm \pi}}\int_{\theta=\kappa\theta_{i}(t)}^{\theta_{f}(t)}\int_{r=0}^{\eta(t)}r^{2}\sin\theta\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\phi\right]. \end{align}

We neglect $O(\tilde {\epsilon }^{4})$ terms, related to the inlet and the outlet section; thus, (4.22) is normalized and simplified to

(4.23)\begin{equation} \frac{{\rm d}\lambda}{{\rm d}T}=\frac{1}{4\lambda^{2}}\left[Q_{in}(T)-\kappa Q_{out}(T)\right]. \end{equation}

The time-dependent unknown functions in (4.12), $\mathcal {A}_{n}(T)$ and $\mathcal {C}_{n}(T)$ are calculated by imposing the boundary conditions (4.20), where the latter are developed into a generalized Fourier series of Legendre or Gegenbauer polynomials,

(4.24a,b)\begin{equation} V_{R}\left(\lambda,\theta;T\right)=\sum_{n=1}^{\infty}\varLambda_{n}(T)\mathcal{P}_{n}(\cos\theta),\quad V_{\theta}\left(\lambda,\theta;T\right)=\frac{1}{\sin\theta}\sum_{n=2}^{\infty}\varphi_{n}(T)\mathcal{J}_{n}(\cos\theta). \end{equation}

The general Fourier coefficients are

(4.25)\begin{equation} \left. \begin{aligned} \varLambda_{n}(T) & = \frac{2n+1}{2}\int_{{-}1}^{1}{V_{R}\left(\lambda,\xi;T\right)\mathcal{P}_{n}(\xi)\,\mathrm{d}\xi}\equiv\frac{1}{\lambda^{2}\tilde{\epsilon}^{4}}\left[\tilde{\varLambda}_{n}^{(in)}Q_{in}+\tilde{\varLambda}_{n}^{(out)}Q_{out}\right],\\ \varphi_{n}(T) & = \frac{n(n-1)(2n-1)}{2}\int_{{-}1}^{1}{V_{\theta}\left(\lambda,\xi;T\right)\frac{\mathcal{J}_{n}(\xi)}{\sqrt{1-\xi^{2}}}\,\mathrm{d}\xi}\equiv\frac{1}{\lambda^{2}\tilde{\epsilon}^{4}}\left[\widetilde{\varphi}_{n}^{(in)}Q_{in}+\widetilde{\varphi}_{n}^{(out)}Q_{out}\right], \end{aligned} \right\} \end{equation}

where

(4.26) \begin{align} \left. \begin{aligned} \tilde{\varLambda}_{n}^{({\cdot})} & =\frac{2n+1}{2}\int_{I_{({\cdot})}}{ f(\xi;\tilde{\epsilon})\tilde{\omega}(\xi;\tilde{\epsilon})\mathcal{P}_{n}(\xi)\xi\,\mathrm{d}\xi},\\ \widetilde{\varphi}_{n}^{({\cdot})} & =\frac{n(n-1)(2n-1)}{2}\left[\frac{\varPhi^{({\cdot})}\lambda\tilde{\epsilon}^{5}}{4\sqrt{1-\tilde{\epsilon}^{2}}}\int_{{-}1}^{1}{\varTheta(\cos^{{-}1}(\xi))\frac{\mathcal{J}_{n}(\xi)}{\sqrt{1-\xi^{2}}}\,\mathrm{d}\xi}\right.\\ &\left.\quad-\int_{I_{({\cdot})}}{f(\xi;\tilde{\epsilon})\tilde{\omega}(\xi;\tilde{\epsilon}){J}_{n}(\xi)\,\mathrm{d}\xi}\right], \end{aligned} \right\} \end{align}

such that $f(\xi ;\tilde {\epsilon })=2(\xi ^{2}-1+\tilde {\epsilon }^{2})-\gamma (4-\tilde {\epsilon }^{2}-4\xi ^{2}+3\tilde {\epsilon }^{-2}(\xi ^{2}-1)^{2})$, the weighted function is $\tilde {\omega }(\xi ;\tilde {\epsilon })=\omega ({S=\tilde {\epsilon }^{-1}\sqrt {1-\xi ^{2}}})$,$\varPhi ^{(in)}=1$, $\varPhi ^{(out)}=-1$, and the integration's intervals are $I_{(in)}=[-1,-\sqrt {1-\tilde {\epsilon }^{2}}]$ and $I_{(out)}=[\sqrt {1-\tilde {\epsilon }^{2}},1]$. Note that the first integral expression of the $\widetilde {\varphi }_{n}$ coefficient in (4.26) is negligible with respect to the second integral expression; therefore, this term can be omitted for simplicity.

4.3. Formulation of the series solution

Substitution of (4.24a,b)–(4.26) into (4.12a) and (4.12b), where $R=\lambda$ yields two linear equations that define the unknown functions $\mathcal {A}_{n}(T)$ and $\mathcal {C}_{n}(T)$. This provides the solution for the velocity field and the solution is rewritten as

(4.27a)$$\begin{gather} V_{R} = \frac{1}{2}\sum_{n=2}^{\infty}\left[ \left((n+2)\chi^{{-}2}-n\right)\varLambda_{n-1}+\left(\chi^{2}-1\right)\varphi_{n}\right]\chi^{n}\mathcal{P}_{n-1}(\cos\theta), \end{gather}$$

(4.27b)$$\begin{gather}V_{\theta} = \frac{1}{2}\sum_{n=2}^{\infty}\left[ \left(\chi^{{-}2}-1\right)n(n+2)\varLambda_{n-1}+\left(n(1-\chi^{{-}2})+2\right)\varphi_{n}\right]\chi^{n}\frac{\mathcal{J}_{n}(\cos\theta)}{\sin\theta}, \end{gather}$$

where $\chi (R;T):=R/\lambda (T)\in [0,1]$. The general solution of the pressure distribution can also be rewritten as follows:

(4.28)\begin{equation} P(R,\theta;T)=P_{0}(T)-\frac{\varepsilon}{\lambda}\sum_{n=1}^{\infty}{\frac{(2n+3)\left((n+1)\varLambda_{n}+\varphi_{n+1}\right)}{n}\chi^{n}\mathcal{P}_{n}(\cos\theta)}. \end{equation}

This expression represents the pressure distribution inside the chamber, with unknown time-dependent functions. Now $P_{0}(T)$ may be determined according to the physical boundary conditions of the pressures acting on the flow – at the inlet and outlet. The stretch, $\lambda (T)$, (and therefore also $\varLambda _{n}(\lambda )$ and $\varphi _{n+1}(\lambda )$), may be determined by the wall's hyperelastic constitutive model. In the next sections, we present two analyses describing different physical cases, including the specific hyperelastic constitutive model we used.

6. The dynamic behaviour of two interconnected bistable chambers

In this section, we analyse the behaviour of a system consisting of two chambers, serially connected through slender tubes to a single inlet, whose flow rate is dictated and equal to $Q_{in}(T)\equiv Q(T)$. This system is instrumental for understanding the behaviour of interconnected bistable elements, and it sheds light on the capability to govern the constituent elements by employing a single input (Ben-Haim et al. Reference Ben-Haim, Salem, Or and Gat2020). In many works in which flow-controlled bistable elastic systems have been examined, the main assumption is uniform pressure within the elastic element (Treloar Reference Treloar1975; Dreyer et al. Reference Dreyer, Müller and Strehlow1982; Glozman et al. Reference Glozman, Hassidov, Senesh and Shoham2010; Ben-Haim et al. Reference Ben-Haim, Salem, Or and Gat2020). Here, we shall analyse the dynamics of such systems using the solution we developed in the previous sections, which considers the pressure distribution in the elastic element's inner space. This system's physics is described utilizing the analyses presented above, where we study identical tubes and chambers. The system under investigation combines two of the cases analysed in the previous section, as the flow from the inlet to the first chamber is dictated. We denote the stretch of the first and the second chambers as $\lambda _{1}(T),\lambda _{2}(T)$, respectively. From the integral mass rate balance,

(6.1)\begin{equation} Q(T)=4\lambda_{1}^{2}\frac{{\rm d}\lambda_{1}}{{\rm d}T}+4\lambda_{2}^{2}\frac{{\rm d}\lambda_{2}}{{\rm d}T}. \end{equation}

The stretch $\lambda _{2}(T)$ of the second chamber is governed by (5.12), where $P_{in}(T)$ is related to the pressure in the point in the tube relative to the connection with the (second) chamber. In our case $P_{in}^{(eff)}(T)=P(\lambda _{1},0;T)$ is the effective external pressure, where $P(R,\theta ;T)$ is given by (5.5). Recalling that the Fourier coefficients $\varLambda _{n}$ and $\varphi _{n}$ are linear combinations of the inlet the outlet flow rates (4.26), the effective external pressure is rewritten as

(6.2)\begin{equation} P_{in}^{(eff)}(T) = \left.\frac{1}{\lambda_{1}^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right|_{\lambda_{1}}+\frac{\varepsilon}{\lambda_{1}^{3}\tilde{\epsilon}_{1}^{4}}\left(4\lambda_{2}^{2}\frac{{\rm d}\lambda_{2}}{{\rm d}T}\varPi^{(out)}(\tilde{\epsilon}_{1})+Q(T)\varPi^{(in)}(\tilde{\epsilon}_{1}) \right), \end{equation}

where

(6.3)\begin{equation} \varPi^{({\cdot})}(\tilde{\epsilon}_{1}):=\sum_{n=1}^{\infty}\frac{(2n+3)\left((n+1)\tilde{\varLambda}_{n}^{({\cdot})}+\widetilde{\varphi}_{n+1}^{({\cdot})}\right)}{n}\left[-\frac{1}{2}\left(\mathbb{P}_{n}^{(in)}+\mathbb{P}_{n}^{(out)}\right)-1\right] \end{equation}

and $\tilde {\epsilon }_{1}(T)=\epsilon /\lambda _{1}(T)$. Utilizing the results achieved for the different values of $\tilde {\epsilon }_{1}$ yields the following approximations: $\varPi ^{(out)}(\tilde {\epsilon }_{1})=O(\tilde {\epsilon }_{1})$ and $\varPi ^{(in)}(\tilde {\epsilon }_{1})=O(\tilde {\epsilon }_{1}^{2})$. Using the governing equation of $\lambda _{2}(T)$ in (5.12) yields the second equation of motion,

(6.4)\begin{equation} \frac{{\rm d}\lambda_{2}}{{\rm d}T}=\left.\frac{1}{32\lambda_{2}^{2}}\left(\frac{1}{\lambda_{1}^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right|_{\lambda_{1}}-\left.\frac{1}{\lambda_{2}^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right|_{\lambda_{2}}\right)+O(\epsilon_t). \end{equation}

Equations (6.1) and (6.4) are a set of nonlinear coupled first-order differential equations that govern the evolution of the chamber's stretches $\lambda _{i}(T)$ under the single input $Q(T)$. In our previous work (Ben-Haim et al. Reference Ben-Haim, Salem, Or and Gat2020), we have presented an algorithm whose purpose is to bring the system from one equilibrium state to another by a single input. There, it is assumed that the process is quasi-static; thus, the chamber's pressure is uniform during the process. Thanks to the analysis presented in this work, it is possible to consider the pressure distribution and the chambers’ flow field during the dynamic process. The equilibrium state of this system is achieved when the flow rate is zero. In this case, the derivatives in time are equal to zero, and the differential equations degenerate into a single algebraic equation that defines the equilibrium curves. The equilibrium curves are defined by

(6.5)\begin{equation} \left.\left(\frac{1}{\lambda_{1}^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right)\right|_{\lambda_{1,SS}} = \left.\left(\frac{1}{\lambda_{2}^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right)\right|_{\lambda_{2,SS}}. \end{equation}

Equation (6.5) describes the equilibrium curves of the system presented in figure 10(a). Importantly, this nonlinear equation gives rise to two solutions. One solution is given by $\lambda _{1,SS}=\lambda _{2,SS}$, where the radii of both chambers are equal, whereas in the second one the radii are different, $\lambda _{1,SS}\neq \lambda _{2,SS}$, thanks to the bistability of the chambers. When the system is initially placed out of equilibrium (by setting zero input, or response to initial condition), the solution moves along the curve of constant total volume, $\lambda _{1}^{3}+\lambda _{2}^{3}=\text {const.}$, and converges toward stable equilibrium branches. In our previous work (Ben-Haim et al. Reference Ben-Haim, Salem, Or and Gat2020), we have presented an algorithm whose purpose is to bring the system from one equilibrium state to another using a single input. It assumed that the process is quasi-static; thus, the chamber's pressure is hydrostatic during the process. Thanks to the analysis presented in this work, it is possible to consider the pressure profile and the chambers’ flow field during the dynamic process. In figure 10, we used our algorithm and presented a scenario where the system undergoes an irreversible sequence of transitions between the chambers’ combined states while being controlled by a single input of flow rate $Q(T)$. The chosen input is piecewise constant, represented in figure 10(c). Figure 10(b) shows the system's trajectory in the $\{\lambda _{2},\lambda _{2}\}$ plane, overlaid on the equilibrium curves. The plots show how the system goes through the irreversible sequence of states. These state transitions are made possible by exploiting the following critical effect. When the state trajectory follows a stable branch and reaches a point where it becomes unstable, the trajectory rapidly ‘jumps’ and converges to a stable branch, moving very close to a cubic arc of constant total volume. Figure 11 shows the flow velocity and pressure distribution corresponding to an attractive instance in which the system passes through point A, presented in time and on the $\{\lambda _{2},\lambda _{1}\}$ plane, in figure 10. When the system reaches point A, the flow rate from the first chamber (the one connected to the inlet tube) towards the second chamber is spontaneously bigger than the inlet flow rate. Therefore, the first chamber is deflated.

Figure 10. (a) Equilibrium curves of the two-chamber system in the $\{\lambda _{1};\lambda _{2}\}$ plane. Black solid curves are stable branches and black dashed curves are unstable ones for $\alpha =0.1$. The evolution of the equilibrium curves by the $\alpha$ parameter is described by the grey curves. The blue line is the solution trajectories of numerical simulations of (6.1) and (6.4), overlaid on the branches of equilibrium curves. The red and green points describe the moments in which the flux is changed. (b) Time plots of chambers’ stretch $\lambda _{i}(T)$ obtained by numerical integration of the nonlinear dynamical system. (c) Time plot of inlet flow $Q$ for inflation and deflation in the case of two chambers.

Figure 11. Series solution of a system consisting of two chambers controlled by single inlet: (a) velocity field and streamlines; (b) pressure distribution. The flow velocity and pressure distribution corresponding to when the system passes through the point A shown in figure 10.

7. Concluding remarks

In this work, we analysed the dynamics of creeping flow in a bistable hyperelastic spherical chamber by calculating the velocity field and pressure distribution inside the chamber. The analytical results were compared with numerical simulations, which showed an excellent fit. From the normalization of the governing equations (4.1) and (4.2), we obtained the condition $\varepsilon =a^{4}/r_{0}^{3}\ell \ll 1$ which led to creeping flow inside the chambers. Moreover, we obtained the condition $q^{*}\ll \mu a/\rho$ for inertial effects to be negligible in the vicinity of the connection to the tubes.

In order to describe the coupled model of viscous–elastic dynamics, we first focused on the nonlinear constitutive elastic laws (the Mooney–Rivlin model). Next, using the analytical series solution, we studied the dynamic responses for two physical different cases. The first case was dictated volumetric flux. Based on the mechanical energy principle, we formulated the chamber's pressure distribution. We obtained the characteristic pressure in the elastic chambers as $p^{*}=w_{0}\psi ^{*}/r_{0}$ which depends on the hyperelastic model we used.

The pressure is spatially uniform at the leading order. Mainly, the leading order of the problem is a case of fully developed uniform pressure without any velocities. More spatial and temporal pressure and velocities are created in the dynamic case where the chamber is inflated or deflated. Based on the numeric simulation, we saw that close to the chamber's wall ($\chi \rightarrow 1$), the viscous flow will generate $O(\epsilon _{t})$ spatially varying corrections. In this analysis, we obtained a non-intuitive result – the pressure solution on the chamber's wall obtained its maximum value at an intermediate value of $\theta$ away from the poles $\theta =0$ and $\theta ={\rm \pi}$. This qualification may be critical in identifying the failure point of the chamber's wall.

In the second physical case, the flow was driven by an imposed inlet pressure. There we have simplified and reduced the time-varying dynamical equations to one compact equation depending on the elastic model, dictated pressure, chamber stretch and inlet tube slenderness (5.12). Although the finite element simulation showed that the chamber undergoes a slightly different deformation from an ideal sphere (making it slightly pear-shaped), the effect is small and the solution obtained in this work are an excellent approximation. However, in order to characterize the obtained geometric shape, more extensive analysis is required, which goes beyond the scope of this work.

In the last part of this work, we presented an investigation of a system consisting of two interconnected coupled chambers controlled by the flow at the chamber's inlet. This system demonstrates the bistability feature through which the chambers can display controlled transitions between different multistable states using a single input of controlled flow rate.

We have developed an analytical model that can serve as the basis for the physical understanding of acinar fluid mechanics in a rhythmically expanding spherical alveolus and its vicinity. These results allow the modelling of flows in applications such as the inflation of balloons in medical procedures and pulmonary drug delivery optimization to target specific lung regions. Moreover, the results might be leveraged to analyse the dynamics of particles inside spherical elastic chambers in low Reynolds flow as future work.