Viscous flow fields in hyperelastic chambers

Abstract

Viscous flows in hyperelastic chambers are relevant to many biological phenomena such as inhalation into the lung's acinar region, and medical applications such as the inflation of a small chamber in minimally invasive procedures. In this work, we analytically study the viscous flow and elastic deformation created due to inflation of such spherical chambers from one or two inlets. Our investigation considers the shell's constitutive hyperelastic law coupled with the flow dynamics inside the chamber. For the case of a narrow tube filling a larger chamber, the pressure within the chamber involves a large spatially uniform part, and a small-order correction. We derive a closed-form expression for the inflation dynamics, accounting for the effect of elastic bistability. Interestingly, the obtained pressure distribution shows that the maximal pressure on the chamber's surface is greater than the pressure at the entrance to the chamber. The calculated series solution of the velocity and pressure fields during inflation is verified by using a fully coupled finite element scheme, resulting in excellent agreement. Our results allow the estimation of the chamber's viscous resistance at different pressures, thus enabling us to model the process of inflation and deflation.

1. Introduction

The inflation of elastic balloons has been extensively investigated in the past, mainly because the corresponding dynamics depend on both the flow and the balloon's material elasticity model. The inflation of a toy balloon or a spherical membrane was studied thoroughly by Beatty (1987). In his work, Beatty (1987) has presented an analysis of an incompressible, isotropic hyperelastic spherical pressurized membrane. According to his work, and similar results by Treloar (1975), the Mooney–Rivlin elasticity model successfully captures most of the overall physical effect. The majority of the research done so far considered hydrostatic uniform pressure distribution within the chamber and the determination of pressure as a constant parameter that uniformly affects the elastic walls (Treloar 1975; Needleman 1977; Beatty 1987; Mangan & Destrade 2015; Vandermarlière 2016; Hines, Petersen & Sitti 2017).

Balloons with controlled inflation are used in medical applications such as pleural pressure assessments (Milic-Emili et al. 1964) and enteroscopy (Yamamoto et al. 2001). A recent study by Manfredi et al. (2019) shows a promising biomedical application of a soft robot for a colonoscopy, which utilizes a double-balloon system for achieving inchworm-like crawling while bracing against the colonic walls. Haber et al. (2000) investigated alternating shear flow over a self-similar, rhythmically expanding hemispherical depression. Quasi-steady creeping flow in models of small airway units of the lung was investigated by Davidson & Fitz-Gerald (1972). Ilssar & Gat (2020) studied the inflation and deflation dynamics of a liquid-filled hyperelastic balloon, focusing on inviscid laminar flow. In those systems, the characteristic time it takes for the pressure to reach a constant uniform value in a chamber is assumed to be much shorter than the time it takes for the fluid to pass through the tubes (based on the viscous resistance). However, to assess the fluid and elastic shell's dynamics, a complete mathematical model describing the system's fluid–structure interaction at low Reynolds numbers is needed. The study of the fluid–structure interaction dynamics of low-Reynolds-number incompressible liquid flows and elastic structures may help introduce a new level of control in fluid–structure-based autonomous systems due to the presence of viscous force (Elbaz & Gat 2014, 2016).

In the soft-robotics field, recent studies show the propulsion of elastic structures embedded with internal cavities while controlling pressures or flow rates at the network's inlets (Fei & Gao 2014; Overvelde et al. 2015; Fei & Pang 2016; Gamus et al. 2017; Gorissen et al. 2019; Siefert et al. 2019; Ben-Haim et al. 2020; Salem et al. 2020). In the case of fluidic actuation, several works study variations of the well known ‘two-balloon system’, whereas others study networks of multiple connected chambers (Treloar 1975; Dreyer, Müller & Strehlow 1982; Glozman et al. 2010; Ben-Haim et al. 2020). As shown by these studies, for some given values of the pressure, multiple solutions for the volume are possible. Since the hyperelastic spherical membranes are multistable systems, it allows us to selectively inflate each balloon to one of its stable states by varying the input according to a particular carefully synthesized profile. Consequently, it can pave the way toward manufacturing soft robots that utilize minimal actuation to produce highly complex locomotion.

In this work, we examine the effect of elasticity on transient creeping flow in the bistable hyperelastic chambers. The chamber is assumed to be an ideal sphere. The Stokes equation governs the flow field, while Mooney–Rivlin constitutive laws model the elastic chamber. The fluidic pressure within the balloon is not uniform and cannot be directly determined from the known Mooney–Rivlin relation. The fluidic pressure distribution in the chamber is estimated by balancing the fluidic pressure with the total force on the elastic membrane. Since this force is obtained by integrating the pressure distribution (which depends on the angle $\theta$), it receives a different value than the value obtained in the hydrostatic models.

This work's structure is as follows. In § 2, the geometry, relevant parameters and physical assumptions are defined. In § 3, the hyperelastic Mooney–Rivlin constitutive law is presented. The strain energy function is analysed in order to present the bistable phenomenon. Section 4 presents closed-form solutions of the governing equations, describing the flow field within an expanding chamber. In § 5, two different physical cases are described. The case of dictated inlet mass rate coupled by the hyperelastic model is described in § 5.2, where numerical verification of the fully coupled model is presented. In § 5.3, we present the second case where the inlet pressure is dictated and the stretching of the hyperelastic sphere is governed by the flow dynamics. Section 6 examines the dynamic behaviour of two interconnected bistable chambers. Concluding remarks are presented in § 7.

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 1972).

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 1972; Beatty 1987)

(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 1987; Treloar 1975). 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 1975; Beatty 1987; Ben-Haim et al. 2020) 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 2012). 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 (2012) 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.

5. Results

In this section we obtain the full series solution under various different boundary conditions. First, we present an empirical estimation of $\gamma$ by a problem of flow from an injection tube into a half-space. Then, we present two analyses that describe different physical cases. In the first case, we present a chamber whose volumetric flow rate is dictated while the pressure is governed by the chamber's hyperelastic shell. In the second case, we present a chamber whose input pressure is dictated while the wall's hyperelastic law governs the time-varying chamber's radius.

5.1. Estimation of $\gamma$ – flow from an injection tube into a half-space

Assuming the radius of the tube is small relative to the chamber's radius ($\epsilon _{a}\ll 1$), the boundary's polar angle at the connection point is ${\rm \pi} -\theta _{f}=\sin ^{-1}{(\tilde {\epsilon })}= O(\tilde {\epsilon })$ for the inlet tube, and $\theta _{i}=\sin ^{-1}{(\tilde {\epsilon })}= O(\tilde {\epsilon })$ for the outlet tube. We focus on the flow field very close to the tube–chamber connection where $\chi \rightarrow 1$ and $\theta \rightarrow {\rm \pi}$ for the inlet tube or $\theta \rightarrow 0$ for the outlet tube. Since $\tilde {\epsilon }\ll 1$, it can be modelled as flow from an injection tube into a half-space.

Several works present numerical analyses for such problems (Weissberg 1962; Fitz-Gerald 1972; Tutty 1988; Sisavath, Jing & Zimmerman 2001), but there are no analytical solutions to the best of our knowledge. Thus, in order to examine the flow field in the connection region, we utilize finite element schemes devised in COMSOL Multiphysics. In these simulations, the entrapped fluid is modelled according to the Navier–Stokes equations, assuming that the flow is incompressible and isothermal. A comprehensive explanation of the numerical simulation will be described in the next part of the results section, while here we restrict ourselves to describing the basic geometry. We assumed a long tube (a hundred diameters in length) with unit radius. The boundary conditions are non-slip and non-penetration into the tube's edge.

The modified velocity profile we shall use was defined in (4.17); after setting the non-dimensional variables and neglecting $O(e^{-1/\epsilon _\omega })$ terms, the obtained normalized relation is

(5.1) \begin{equation} \frac{U_{z}^{(m)}}{Q(T)}=2(1-S^{2})+\gamma(1-4S^{2}+3S^{4}). \end{equation}

Hence, the parameter $\gamma$ should be found by fitting the profile (5.1) to the numerical simulation results, using the least squares method. The best fitted value is $\gamma =-0.1$. The velocity profile obtained in the numeric simulation and the modified velocity profile we fitted in (5.1) is shown in figure 3. All the parameters used in the simulations are elaborated in table 1. During the further analysis, we will use the same approximation (with $\gamma =-0.1$). This approximation will be valid as long as $\tilde {\epsilon }\rightarrow 0$.

Figure 3. Modified velocity profile (4.17) used as a boundary condition of the chamber solution. The grey dotted curve is the known parabolic Hagen–Poiseuille profile, obtained by setting $\gamma =0$ in the modified velocity profile (4.17). The blue line is the velocity profile obtained in the numeric computational fluid dynamics (CFD) simulation of injection tube into a half-space. The red dashed line is the modified velocity profile (5.1), with $\gamma =-0.1$ and $\epsilon _{\omega }=0.05$.

Table 1. Summary of physical parameters values used for plotting the analytical (series) solutions of the dynamic cases. The density and the kinematic viscosity are related to glycerol which is known as a viscous liquid. The geometric parameters are chosen as the typical value of hyperelastic small chambers used in Ben-Haim et al. (2020).

5.2. Case I – dictated inlet flux and hyperelastic wall model

In this case we dictate the volumetric rate into an elastic chamber. The integral mass conservation equation (4.22) is simplified to

(5.2)\begin{equation} \lambda(T)=\left[1+\frac{3}{4}\int_{0}^{T}\left(Q_{in}(\tau)-Q_{out}(\tau)\right)\,\mathrm{d}\tau\right]^{1/3}. \end{equation}

While the stretch function $\lambda (T)$ is known, as described in (5.2), the inlet pressure is not dictated and additional data regarding the pressure distribution is obtained from the hyperelastic constitutive relations. Here, it is worth emphasizing that any general constitutive elastic law can serve as a basis for subsequent development, even if it is not bistable or hyperelastic. Integrating the strain density function, $\psi (\lambda )$, over the volume of the chamber's spherical shell and keeping only leading-order terms yields the leading-order chamber's strain energy,

(5.3)\begin{equation} \iiint_{\mathbb{V}}{\psi}(\lambda)\,\mathrm{d}\mathbb{V}=4{\rm \pi} r_{0}^{2}w_{0}\psi^{*}\times\hat{\psi}(\lambda), \end{equation}

where $\mathbb {V}\approx w\times \mathbb {S}$ is the material (constant) volume of the thin shell, $\psi ^{*}$ is the characteristic value of $\psi$, and $\hat {\psi }$ is the normalized strain energy density function $\hat {\psi }(\lambda )=\psi (\lambda )/\psi ^{*}$. The work done by the surface traction acting between two states without body force is

(5.4)\begin{equation} \iint_{\mathbb{S}}\int_{\xi=r_{0}}^{\eta}\left(p(\xi,\theta)\,\mathrm{d}\mathbb{\boldsymbol{\mathbb{S}}}\right)\cdot\mathrm{d}\boldsymbol{\xi}=2{\rm \pi} r_{0}^{3}p^{*}\int_{\theta=\kappa\theta_{i}}^{\theta_{f}}\int_{\xi=1}^{\lambda}\xi^{2}P(\xi,\theta)\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\xi. \end{equation}

Substitution of the pressure's solution (4.28) into (5.4), and using the mechanical energy principle which states that the work done by the surface tractions acting between two equilibrium states without body force is balanced by the change in the total strain energy (Beatty 1987), allows us to fully define the pressure distribution in the chamber,

(5.5)\begin{align} P_{(I)}(R,\theta;T)&=\frac{1}{\lambda^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}-\frac{\varepsilon}{\lambda}\sum_{n=1}^{\infty} \frac{(2n+3)\left((n+1)\varLambda_{n}+\varphi_{n+1}\right)}{n}\nonumber\\ &\quad \times \left[\frac{1}{2}\left(\mathbb{P}_{n}^{(in)}+\kappa\mathbb{P}_{n}^{(out)}\right) +\chi^{n}\mathcal{P}_{n}(\cos\theta)\right], \end{align}

where $\mathbb {P}_{n}^{(\cdot )}$ is the zeroth moment of $\mathcal {P}_{n}(\xi )$ about an origin, defined as

(5.6a)$$\begin{gather} \mathbb{P}_{n}^{(in)}\equiv\int_{{-}1}^{-\sqrt{1-\tilde{\epsilon}^{2}}}{\mathcal{P}_{n}(\xi)\,\mathrm{d}\xi} = \frac{\mathcal{P}_{n+1}\left(-\sqrt{1-\tilde{\epsilon}^{2}}\right)-\mathcal{P}_{n-1}\left(-\sqrt{1-\tilde{\epsilon}^{2}}\right)}{2n+1}, \end{gather}$$

(5.6b)$$\begin{gather}\mathbb{P}_{n}^{(out)}\equiv\int_{\sqrt{1-\tilde{\epsilon}^{2}}}^{1}{\mathcal{P}_{n}(\xi)\,\mathrm{d}\xi} ={-}\frac{\mathcal{P}_{n+1}\left(\sqrt{1-\tilde{\epsilon}^{2}}\right)-\mathcal{P}_{n-1}\left(\sqrt{1-\tilde{\epsilon}^{2}}\right)}{2n+1}. \end{gather}$$

From order of magnitude analysis in (5.3) and (5.4), we obtain the characteristic pressure, which depends on the specific hyperelastic model, we use

(5.7)\begin{equation} p^{*}=\frac{w_{0}}{r_{0}}\psi^{*}. \end{equation}

According to the solution obtained in (5.5), the pressure distribution inside the chamber consists of two parts. The first part is the well known isotropic pressure obtained by Beatty (1987). This expression represents the isotropic pressure in the leading order, which is experienced by the chamber's elastic wall, assuming that the pressure is uniform and equal to $P_{S}(\lambda )=\lambda ^{-2}\times {\rm d}\psi /{\rm d}\lambda$. The second expression is the transient developing pressure profile, $\varepsilon P_{1}(R,\theta ;\tilde {\epsilon };T)$, which varies spatially and temporally.

5.2.1. Numerical verification of the fully coupled model

In order to validate the theoretical model, we have utilized commercially available software (COMSOL Multiphysics) in order to conduct finite element simulations considering the fully coupled dynamics of the system. In these simulations, the entrapped fluid is modelled according to the Navier–Stokes equations, assuming that the flow is incompressible and isothermal. Moreover, the shell is modelled according to the Mooney–Rivlin model, wherein, contrary to the theoretical model, it is not restricted to a spherical shape. All the parameters used in the simulations are elaborated in table 1. In all simulations, the geometry and boundary conditions are taken in correspondence with the problem statement in § 2. Further, the physics is described by the fluid–structure interaction module of COMSOL, referring to the chamber as a two-parameter Mooney–Rivlin hyperelastic solid. This is done while employing a moving mesh formulation to accommodate the changes of the fluid's domain. The coupling between the solid and the fluid is carried out by balancing the fluid's velocity and the time derivative of the solid's displacements and the normal components of their stress tensors at the interface between the solid and the fluid. The numerical schemes describe the deformation field of the solid utilizing second-order base functions, where those used to discretize the velocity field and pressure distribution of the fluid are cubic and quartic, respectively. Finally, the system's geometry is described as two-dimensional and axisymmetric to eliminate significant numerical errors and decrease the computational effort. Furthermore, the meshing of the geometry was enhanced until low sensitivity to further refinement was achieved. In the final meshing used in the simulations whose results are presented here, the solid is modelled by rectangular elements whose grid contains 200 tangential elements and six elements along with its thickness. Similarly, the tube is modelled by a regular grid having 30 radial and 103 axial elements. Finally, since the geometry of the fluid residing inside the chamber is of higher complexity, the latter is described by free triangular elements. The maximal size of these elements is restricted to 0.18 mm. On the boundaries of this region, where the sensitivity is higher, and the precision is of greater importance, the maximal size is further reduced to 0.1 mm.

While numerically simulating the dynamics of the fluid-filled chamber, utilizing the two constant parameters of the Mooney–Rivlin model used in the theoretical computations, we noticed discrepancies. By examining the numerically computed static pressure-stretch relation, it was apparent that these discrepancies might originate from errors in the hyperelastic module of COMSOL Multiphysics version 5.0, in which all the numerical simulations were carried out. Therefore, to compare the theoretical and numerically simulated dynamic responses of the fully coupled system, the hyperelastic model used in the numerical scheme had to be modified. To calibrate this model, we fitted the pressure-stretch relation of the numerical scheme at a static regime. Namely, we found the appropriate values of the two Mooney–Rivlin constant parameters, leading to the desired theoretical static pressure-stretch relation while inflating the chamber to different radii, corresponding to different stretches. The values of the pressure inside the chamber and its effective radii were taken a long time after the inflation was finished and after all dynamic effects had decayed, including the non-spherical modes. As a result, the values used in the simulations yield a relative error of less than $0.4\,\%$ in the pressure-stretch relation at the presented values.

Figure 4, which compares the theoretical and the numerically simulated pressure distribution, shows an excellent agreement, thus validating the analytical model and its underlying assumptions. In addition, a numerical investigation of the pressure distribution where $\chi \rightarrow 1$ (inner chamber's wall region) showed that the varies spatially and temporally correction is $O(\epsilon _{t})$. This result is consistent with the analytical result we obtained in the next section.

Figure 4. A numerical verification of the fully coupled model. The chamber's radius is 6 mm, and the volumetric flux is $q=200\,{\rm mm}^{3}\,{\rm s}^{-1}$. The dashed red line is obtained by the analytical (series) solution, which was calculated based on the first 100 terms in the series. The blue line is obtained in the CFD simulation. Excellent agreement is observed. (a) Contour curves of the pressure distribution obtained from the simulation (red curves) and the same contour curves obtained from the series solution (dashed blue curves). The black non-spherical (pear-shaped) boundary describes the elastic chamber wall obtained from the simulation. (b) The non-dimensional dynamic pressure along the symmetry axis $z$. (c) The non-dimensional dynamic pressure at the chamber wall as a function of the angle $\theta$. (d) The axial velocity of the flow relative to a cylindrical coordinate system. (e) The axial velocity component of the flow along the axis of symmetry.

5.2.2. Additional dynamic cases

Figure 5 presents the series solution of the flow velocity magnitude with streamlines and the pressure distribution, based on (4.27)–(5.5), and utilizing the parameters in table 1. The expected characteristic pressure is $p^{*}= O(1\,{\rm kPa})$, the characteristic flow velocity within the chamber is $v^{*}= O(1\,{\rm mm}\,{\rm s}^{-1})$ and $Re\approx 0.01$. The results are based on summation of $100$ terms in the series solution in (5.5). A remarkable and non-intuitive result was obtained from the pressure solution on the chamber's wall (see figure 4c). The maximum pressure obtained on the elastic wall is not obtained near the inlet tube. This result may be critical in identifying the failure point of the chamber's wall.

Figure 5. Series solution of single inlet chamber (case I): (a) velocity field and streamlines, where the colours describe the magnitude of the flow velocity, and the white dotted lines are the streamlines; (b) pressure distribution with white dash–dotted contour lines of constant pressure. The chamber's radius is 6 mm and the volumetric flux is $q=200\,{\rm mm}^{3}\,{\rm s}^{-1}$.

The theoretical solution also readily allows analysing chambers with two inlets, both with controlled volumetric flow rates. In the first case, whose typical flow velocity magnitude and pressure distributions are presented in figure 6, the chamber is inflated by equal flow rates from both inlets, whereas in the second case, whose behaviour is presented in figure 7, the flow rate in one of the inlets is reversed, meaning that the chamber is inflated and deflated simultaneously. In the first case, each streamline is directed to the chamber's wall; thus, the inflation rate is maximal. In the second case, even though the rate of the inlet and outlet are equal, and thus the volume are constant, a pressure distribution is developed that varies in space and not in time.

Figure 6. Series solution of a double inlet chamber (case I): (a) velocity field magnitude (colourmap) and streamlines (white dotted curves); (b) pressure distribution. The chamber's radius is 6 mm and the volumetric flux is $q=200\,{\rm mm}^{3}\,{\rm s}^{-1}$ in both inlets. In this case, each streamline is directed to the chamber's wall.

Figure 7. Series solution of a chamber with inlet and outlet (case I): (a) velocity field magnitude (colourmap) and streamlines (white dotted curves); (b) pressure distribution. The chamber's radius is 6 mm and the volumetric flux is $q=200\,{\rm mm}^{3}\,{\rm s}^{-1}$ in both inlet and outlet. Since the flow rate of the inlet and outlet are equal, a steady-state solution is obtained for the flow and pressure fields. After an initial transient, the flow and pressure fields reach a steady state.

5.3. Case II – dictated inlet pressure and hyperelastic wall model

In this case, we dictate only the inlet pressure, and solve for the fluidic pressure distribution, as well as the chamber's stretch, $\lambda (T)$.

We use integral mass conservation (4.22) in order to calculate the pressure at the location in which the tube connects with the chamber,

(5.8)\begin{equation} P_{C}(T)\approx P_{in}(T)-32\lambda^{2}\frac{{\rm d}\lambda}{{\rm d}T}. \end{equation}

Similarly to the previous case, we substitute $(R,\theta )=(\lambda,{\rm \pi} )$ in the pressure distribution solution (4.28) and equate it to $P_{C}(T)$ in order to solve for the function $P_{0}(T)$. The pressure distribution obtained by these means is

(5.9)\begin{align} &P_{(II)}(R,\theta;T)=P_{in}(T)-32\lambda^{2}\frac{{\rm d}\lambda}{{\rm d}T}\nonumber\\ &\quad +\frac{4\varepsilon}{\lambda\tilde{\epsilon}^{4}}\frac{{\rm d}\lambda}{{\rm d}T}\sum_{n=1}^{\infty}\frac{(2n+3)\left((n+1)\tilde{\varLambda}_{n}^{(in)} +\widetilde{\varphi}_{n+1}^{(in)}\right)}{n}\left[({-}1)^{n}-\chi^{n} \mathcal{P}_{n}(\cos\theta)\right]. \end{align}

Finally, using the mechanical energy principle (5.4), we derive a nonlinear ordinary differential equation governing the time-dependent stretch of the chamber,

(5.10)\begin{equation} \left[32\lambda^{2}-\frac{2\varepsilon}{\lambda\tilde{\epsilon}^{4}}\varUpsilon(\tilde{\epsilon})\right]\frac{{\rm d}\lambda}{{\rm d}T}=P_{in}(T)-\frac{1}{\lambda^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}, \end{equation}

where

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

The function $\varUpsilon (\tilde {\epsilon })$ was estimated numerically by summing the first $10^{5}$ terms in the series for several values of $\tilde {\epsilon }$ corresponding to $\tilde {\epsilon }\in [0.002,0.5]$.

By substituting the asymptotic approximation of $\varUpsilon$, into the differential equation that governs the chamber's stretch (5.10), we obtain its approximated explicit form given by

(5.12)\begin{equation} \frac{{\rm d}\lambda}{{\rm d} T}=\frac{1}{32\lambda^{2}}\left[P_{in}(T)-\frac{1}{\lambda^{2}}\frac{{\rm d}\hat{\psi}}{{\rm d}\lambda}\right]+O(\epsilon_{t}). \end{equation}

Next, we investigate the system's behaviour in (5.12), whose motion is governed by a controlled pressure inlet. In order to validate this model, we have compared its solution obtained utilizing the parameters in table 1 to finite element simulations carried out in COMSOL Multiphysics. The numerical scheme used here is similar to the one utilized in the previous section, but with the imposed flow rate replaced with an imposed piecewise constant pressure. Figure 8 compares the stretch of the chamber in time, as achieved theoretically from the asymptotic equation (5.12) and numerically from the simulation. Since the chamber is not restricted to be spherical in the numerical simulations, the stretch is taken to be its effective value given in Ilssar & Gat (2020) as $\eta _{eff}(t)=\sqrt {A_{S}(t)/4{\rm \pi} }$. Here, $A_{S}(t)$ is the body's surface area obtained in the non-spherical simulations, and $\eta _{eff}(t)$ is the effective radius of an ideal sphere having the same surface area as the non-spherical body. This figure shows an excellent agreement.

Figure 8. A numerical verification of the fully coupled model describing a chamber's dynamic responses to a pressure pulse imposed at a single inlet, as shown in the red line. The continuous blue curve represents the solution obtained in the numerical CFD simulation, and the dashed black line represents the solution obtained by the (5.12), developed in our analysis. An excellent match between the results can be observed.

In order to approximate the characteristic time constant of the system, enabling us to estimate the time it takes to reach a steady-state, we formulate a linear approximation. The linear system describing the system's dynamic response, close to an equilibrium point given by $(\lambda _{SS},P_{SS})$, when the pressure at the inlet is dictated to be $P_{ext}(\tau )$. The linear equation is solved analytically, leading to the following solution:

(5.13)\begin{equation} \lambda_{L}(\tau)=\lambda_{SS}+\left(\lambda(0)-\lambda_{SS}+\frac{1}{4\lambda_{SS}^{2}}\int_{0}^{\tau}{e^{\beta_{I}\tau'}\Delta P_{in}(\tau')}\,{\rm d}\tau'\right)e^{-\beta_{I}\tau}, \end{equation}

where $\Delta \lambda _{L}=\lambda _{L}(\tau )-\lambda _{SS}$ is a small stretch variation around its nominal value $\lambda _{L}$, $\Delta P_{in}(\tau )=P_{in}(\tau )-P_{SS}$ is the pressure variation from its nominal value, and

(5.14)\begin{equation} \beta_{I}=\frac{1}{4\lambda_{SS}^{2}}\frac{{\rm d}P_{SS}}{{\rm d}\lambda_{SS}}=\frac{P_{SS}}{2\lambda^{3}_{SS}}-\frac{3}{\lambda_{SS}^{4}}+\frac{9}{\lambda_{SS}^{10}}-\alpha\left(\frac{1}{\lambda_{SS}^{2}}-\frac{7}{\lambda_{SS}^{8}}\right). \end{equation}

From (5.14), it is clear that the solution branches of $P_{SS}$ in (3.3) are stable equilibria if and only if ${\rm d}P_{SS}/{\rm d}\lambda _{SS}>0$. The stability criterion obtained from the equation's linearization is identical to the stability criterion obtained from energetic considerations in (3.4). Moreover, from the linear solution, the relaxation time can be approximated as $T_{relax}=32/\beta _{I}$. Since the derivative ${\rm d}P_{SS}/{\rm d}\lambda _{SS}$ in the first branch (I) of the typical pressure-stretch curve, which is plotted in figure 2, is significantly higher than in the third branch (III), the dynamic response in the third region is much slower.

In order to achieve a better approximation, (5.12) is approximated by a quadratic Taylor expansion around the general equilibrium point. When the pressure at the inlet equals the steady-state pressure, $\Delta P_{in}(\tau )=0$, the solution of this equation under a small initial perturbation from equilibrium $\lambda _{Q}(T)$ is given by

(5.15)\begin{equation} \lambda_{Q}(\tau)=\lambda_{SS}+\frac{\beta_{I}}{-\beta_{II}+\left(\beta_{II}+\dfrac{\beta_{I}}{\lambda(0)-\lambda_{SS}}\right)e^{\beta_{I}\tau}}, \end{equation}

where

(5.16)\begin{equation} \beta_{II}={-}\frac{3P_{SS}}{4\lambda^{4}_{SS}}+\frac{6}{\lambda_{SS}^{5}}-\frac{45}{\lambda_{SS}^{11}}+\alpha\left(\frac{1}{\lambda_{SS}^{3}}-\frac{28}{\lambda_{SS}^{9}}\right). \end{equation}

In figure 9, the linear and the second-order approximations were displayed alongside the exact numerical solution of (5.12). The excellent agreement between the exact dynamic response and both approximations indicates that these two approximations are suitable for the prediction of the chamber's evolution.

Figure 9. Analytical Solutions: (a) the evolution of the chamber's stretch (solid blue curve), alongside its linear approximation (5.13) (dash–dotted green curve) and quadratic approximation (5.15) (dashed yellow curve), where the input pressure (dashed red curve) varies between its extremum values $P_{A}$ and $P_{B}$. (b) The equilibrium pressure-stretch curve (solid blue curve), alongside the lower (dashed red line) and higher (dashed green line) pressure values, corresponding to $P_{A}$ and $P_{B}$. The red dashed line shows the entry pressure value in the first and last stage ($P_{SS}=2.75$) and the green dashed line shows the pressure dictated in the middle stage ($P_{SS}=3.05$). The equilibrium points corresponding to the $P_{SS}=2.75$ pressure are marked with red marks where $\lambda _{SS}^{(I)}=1.25$ and $\lambda _{SS}^{(III)}=4.79$.

In § 3, we have shown that in the pressure range spanning between $P_{A}$ and $P_{B}$, there are three possible equilibrium radii for each constant pressure value. Multiple solutions can be exploited to switch from one equilibrium state to another under the same steady-state pressure while passing through the unstable, middle branch. An example of such a transition is shown in figure 9, where after increasing the input pressure and decreasing it back to its initial value, the chamber's radius does not return to its initial radius. Instead, it retains a larger radius, corresponding to the higher equilibrium state. Initially, the system converges to the equilibrium point $\lambda _{SS}^{I}$ corresponding to $P_{SS}$, in the first branch (I) which is closer to the initial conditions; then, the inlet pressure rises to a value higher than $P_{A}$ to move the chamber to another equilibrium point in the third branch (III). Finally, we decrease the pressure again to the same level as the initial step, $P_{SS}$, so that the system will converge to the second equilibrium point in the third branch (III). These results are consistent with the insights raised in our previous work (Ben-Haim et al. 2020). It is worth emphasizing at the end of this section that the general solutions obtained in relations (4.27), (5.5) and (5.9), are parametrically dependent on the dictated strain energy density function $\psi (\lambda )$, based on the chosen constitutive law. Here, we have chosen to present the results using the Mooney–Rivlin hyperelastic constitutive law to illustrate the phenomenon of bistability and to compare our results with some other works which have assumed uniform pressure. However, any other constitutive law for the elastic shell can be used in order to derive the solution for the pressure distribution and velocity field in the three cases illustrated in the section.

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. 2020). In many works in which flow-controlled bistable elastic systems have been examined, the main assumption is uniform pressure within the elastic element (Treloar 1975; Dreyer et al. 1982; Glozman et al. 2010; Ben-Haim et al. 2020). 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. 2020), 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. 2020), 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.