{'fraction_non_alphanumeric': 0.04163828808099402, 'fraction_numerical': 0.013327197422917625, 'mean_word_length': 4.544600938967136, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Viscous flow-fields in hyperelastic Chambers

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 bi-stability. 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 estimating the chamber's viscous resistance at different pressures, thus enabling us to model the process of inflation and deflation.


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 (Needleman 1977;Treloar 1975;Beatty 1987;Vandermarlière 2016;Hines et al. 2017;Mangan & Destrade 2015).
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 (Ben-Haim et al. 2020;Salem et al. 2020;Gamus et al. 2017;Siefert et al. 2019;Fei & Gao 2014;Fei & Pang 2016;Overvelde et al. 2015;Gorissen et al. 2019). 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 (Ben-Haim et al. 2020;Dreyer et al. 1982;Treloar 1975;Glozman et al. 2010). 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 multi-stable 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. 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 θ), 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 analyzed in order to present the bistable phenomena. 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 behavior of two interconnected bi-stable chambers. Concluding remarks are presented in § 7. Figure 1: Illustration of the system under investigation, consisting of a hyper-elastic 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 center is allowed to move.

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 liquidfilled 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 two rigid tubes with radius a and length ℓ. 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 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 v(t, s, z), while its constant density and dynamic viscosity are denoted ρ and µ. The chamber's dynamics are approximated by a single degree of freedom, represented here by the chamber's instantaneous radius, denoted η(t). For spherical geometry, a coordinate system is chosen so that one of the coordinates remains constant on the boundary. Here, {r, θ, φ} are the coordinates of a moving spherical system, located at the center of the chamber, where θ is the polar angle, measured from the axis of symmetry to the radial coordinate r, and φ is the azimuthal angle, revolving around the axis of symmetry, z. The Cauchy-stress tensor of the flow is denoted as σ 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 ≪ r 0 .
The following analysis utilizes three small parameters, including the ratio between the radius of the tubes and the radius of the stress-free chamber (ε a , denoted hereafter by tube-chamber radii ratio), ǫ a = a r 0 ≪ 1, (2.1) the slenderness of the tubes (ε t , denoted hereafter by tube slenderness), (2.2) and the last small parameter in the analysis is taken as the ratio between the viscous stresses and the overall pressure in the chamber (ε, denoted hereafter by chamber viscous resistance parameter ), defined by 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: where λ(T ) denotes the stretch of the chamber, and R is the normalized radial coordinate.

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 ψ(λ), which depends only on the relative stretch λ(T ). Moreover, the elastic strain energy density satisfies ψ(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 ≈ r 2 0 w 0 . Thanks to this relation, the chamber's instantaneous thickness is eliminated. To capture the chamber's bi-stability, 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), where α = 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 ≈ 1.5 MPa and s 2 ≈ 0.15 MPa (Beatty 1987;Treloar 1975). The normalized equation (3.1) was obtained by using the normalization ψ * = s 1 for the strain energy density function, and the magnitude of the parameter α is O(10 −1 ).
We first study the chamber's static behavior, where the pressure (without flow) is dictated. In this case, both the stretch and the pressure are constant, denoted here by λ SS and P SS , respectively. The behavior mentioned above can be demonstrated by the overall effective potential energy of the system, 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 ∂U/∂λ SS = 0 at constant pressure P SS , yields a relation between stretch, λ SS , and pressure, P SS , in equilibrium condition, This well-known relation was extensively leveraged to describe the quasi-static inflation of spherical balloons (Beatty 1987;Treloar 1975;Ben-Haim et al. 2020) for spatially uniform pressures. As seen from Figure 2(a), showing the relation in (3.3) with α = 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 (λ SS ) in Figure 2(b) has two bifurcation points, described by a local maximum point at (λ A , P A ), and a local minimum point at (λ 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 , λ A and λ B , appears in Appendix A. The evolution of those extrema as a function of the small parameter α is presented in Figure 2 Analyzing the strain energy function U(λ SS ; P ss ) in Equation (3.2), using the second derivative with respect to λ SS , it can be proven that the right and left branches of P SS (λ SS ) where 1 < λ < λ A or λ B < λ are stable equilibria and satisfy, Conversely, the intermediate branch λ A < λ SS < λ B is an unstable region satisfying ∂ 2 U/∂λ 2 SS < 0. This is precisely the bi-stability phenomenon.

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

Formulation and analysis of the governing equations
Under the assumptions discussed above, the momentum and continuity equations governing the fluid's behavior expressed in the moving spherical frame: where the third term in the left expression in the momentum equation (4.1a) describes the acceleration of the moving spherical frame, centered on the chamber's moving center, relative to a stationary frame. Assuming the flow in tubes is fully developed and axisymmetric, the volumetric flux is given by, where ∂p t /∂z is the pressure gradient along the tube. Since the tube slenderness ǫ t ≪ 1 we shall assume a constant pressure gradient. Normalization of (4.2) yields the characteristic flow rate as q * = πa 2 u * z where u * z = a 2 p * /µℓ 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 * = ǫ 2 a u * z . Substituting the characteristic values into the chamber viscous resistance parameter (2.3), relates it to the other small parameters, as follows From relation (4.3) it is clear that ε 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., ρgr 0 /p * ≪ 1 (where g is the gravitational acceleration), and define a Reynolds number in chamber as Re = ρv * r 0 /µ ≪ 1. Since the Reynolds number is small, the flow's inertia may be neglected. Therefore, by utilizing the non-dimensional quantities specified in (2.4) the fluid's motion (4.1) is governed by Stokes equations for creeping flow with an implicit time variable, 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 ρq * /µL < ρq * /µa. Thus, a sufficient condition for global neglect of inertia is, From the non-dimensional relation (4.4), 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 parameters' physical meaning, we may use the dynamical stress tensor in the fluid domain defined by the constitutive relation σ f = −pI + µ ∇v + (∇v) T where I is the 3 × 3 unit matrix. In the most general constitutive equation, σ f consisting of the linear and instantaneous dependence of the deviatoric stress, plus the hydrostatic stress, −pI, stemming from the static pressure. Normalization of the total stress tensor yields, As one can notice from equation (4.6) the velocity field is not included at the leadingorder. 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.4) can be solved in spherical polar coordinates using a series expansion (Happel & Brenner 2012). Consider the non-dimensional Stokes stream function Ψ (R, θ, T ). The flow velocity components V R and V θ are related to the Stokes stream function Ψ (R, θ, T ) through where V R and V θ are the radial and tangential velocity components, respectively. By applying the curl operator to the momentum equation (4.4) and using several simple algebraic manipulations, the Stokes equation can be reduced to a fourth-order bi-harmonic equation obtained in terms of the Stokes stream function as follows: where in spherical coordinates, (4.9) Equation (4.8a) is solved by separation of variables, and equation (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) where A n (T ) and C n (T ) are unknown functions, determined by the boundary conditions, and J n (ξ) 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 P n (ξ) as, In the degenerate cases n = 0, 1 we define J 0 (ξ) = 1 and J 1 (ξ) = −ξ, respectively. Using the definition of the Stokes stream function (4.7), we describe the series solution of the velocity field and pressure distribution as, 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, A n (T ) and 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, A n (T ) and 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 r(θ, t) = η(t)r, wherer = (sin θ cos φ, sin θ sin φ, cos θ) is the radial direction unit vector; thus, the material point's velocity isṙ =ηr + ηθθ, whereθ = (cos θ cos φ, cos θ sin φ, − sin θ) is the polar direction unit vector. We recall that the fluid's behavior 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, where κ indicates whether there are both inlet and outlet tubes (when κ = 1), or only an inlet tube (when κ = 0). Moreover, θ i (t) and θ f (t) are angles corresponding to the connections between the tubes and the chamber (see Figure 1). By simple geometric considerations, we get, . (4.14) Considering the derivative of equation (4.13) with respect to time will lead to the relation betweenθ and θ, θ i , θ f ,θ i andθ f . Therefore, the material point's velocity can be rewritten asṙ 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, Here, we assumed a spherical surface at the connection between the tubes and the chamber. The first vector in equation (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 center (the velocity of the moving spherical frame, centered on the chamber's moving center, relative to a stationary frame) should be subtracted by 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 (hyper-elasticity).
A simple investigation of the boundary condition's derivative shows singularities at θ = θ i and θ = θ 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 are leading 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: Here, u p z (s, t) is the known parabolic Hagen-Poiseuille profile (4.18) and u c z (s, t; ǫ ω ) is an correction profile defined by where γ is a parameter determined by fitting the velocity profile obtained in a finite elements calculation. The function ω(s; ǫ ω ) is defined by ω(s; ǫ ω ) = − tanh (2/ǫ ω ) + tanh (s/a + 1)/ǫ ω − tanh (s/a − 1)/ǫ ω , and 0 < ǫ ω ≪ 1 is an arbitrary small parameter (denoted hereafter by smoothing parameter ). In fact, ω(s; ǫ ω ) is a weighted function making the velocity profile differentiable even at θ i and θ f . Moreover, it can be shown that ω(s; ǫ ω ) = 1 + O(e −1/ǫω ) as ǫ ω ≪ 1 and |s/a| 1 − O(ǫ ω ). 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 where, The integral mass conservation equation is given by, We neglect O(ǫ 4 ) terms, related to the inlet and the outlet section; thus, equation (4.22) is normalized and simplified to (4.23) The time-dependent unknown functions in equation (4.12), A n (T ) and 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, The general Fourier coefficients are Note that the first integral expression of theφ n coefficient in (4.26) is negligible with respect to the second integral expression; therefore, this term can be omitted for simplicity.

Formulation of the series solution
Substitution of (4.24)-(4.26) into (4.12a-b) where R = λ yields two linear equations that define the unknown functions A n (T ) and C n (T ). This provides the solution for the velocity field and the solution is rewritten as, (n + 2)χ −2 − n Λ n−1 + χ 2 − 1 ϕ n χ n P n−1 (cos θ), (4.27a) The general solution of the pressure distribution can also be rewritten as follows, (2n + 3) (n + 1)Λ n + ϕ n+1 n χ n P n (cos θ). (4.28) This expression represents the pressure distribution inside the chamber, with unknown time-dependent functions. 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, λ(T ), (and therefore also Λ n (λ) and ϕ n+1 (λ)), 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.

Results
In this section we obtain the full series solution under various different boundary conditions. First, we present an empirical estimation of γ 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.

Estimation of γ -Flow from an injection tube into a half-space
Assuming the radius of the tube is small relative to the chamber's radius (ǫ a ≪ 1), the boundary's polar angle at the connection point is π − θ f = sin −1 (ǫ) = O(ǫ) for the inlet tube, and θ i = sin −1 (ǫ) = O(ǫ) for the outlet tube. We focus on the flow field very close to the tube-chamber connection where χ → 1 and θ → π for inlet tube or θ → 0 for outlet tube. Sinceǫ ≪ 1, it can be modeled as flow from an injection tube into a half-space.
Several works present numerical analyses for such problems (????), 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 modeled according to 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 (hundred diameters in length) with unit radius. The boundary condition are non-slip and non-penetration into the tube's edge.
The modified velocity profile we shall use was defined in equation (4.17); after setting the non-dimensional variables and neglecting O(e −1/ǫω ) terms, the obtained normalized relation is U Hence, the parameter γ should be found by fitting the profile (5.1) to the numerical simulation results, using the least-squares method. The best fitted value is γ = −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 γ = −0.1). This approximation will be valid as long asǫ → 0.

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 While the stretch function λ(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 hyper-elastic 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 bi-stable or hyperelastic. Integrating the strain density function, ψ(λ), over the volume   where V ≈ w·S is the material (constant) volume of the thin shell, ψ * is the characteristic value of ψ, andψ is the normalized strain energy density functionψ(λ) = ψ(λ)/ψ * . The work done by the surface traction acting between two states without body force is, 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), allow us to fully define the pressure distribution in the chamber, n is the zero th moment of P n (ξ) about an origin, defined as From order of magnitude analysis in (5.3) and (5.4), we obtain the characteristic pressure, which depends on the specific hyper-elastic model we use, 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 (λ) = λ −2 · dψ/dλ. The second expression is the transient developing pressure profile, εP 1 (R, θ;ǫ; T ), which varies spatially and temporally.

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 modeled according to Navier-Stokes equations, assuming that the flow is incompressible and isothermal. Moreover, the shell is modeled according to the Mooney-Rivlin model, wherein contrary to the theoretical model, it is not restricted to a spherical shape. All  Table 1. In all simulations, the geometry and boundary conditions are taken in correspondence with the problem statement in section § 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 modeled by rectangular elements whose grid contains 200 tangential elements and six elements along with its thickness. Similarly, the tube is modeled 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 Moony-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 Moony-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 period after the inflation was finished after all dynamic effects 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 pressuredistribution, shows an excellent agreement, thus validating the analytical model and its underlying assumptions. In addition, a numerical investigation of the pressure distribution where χ → 1 (inner chamber's wall region) displayed that the varies spatially and temporally correction is O(ǫ t ). This result is consistent with the analytical result we obtained in the next section. Figure 5 presents the series solution of the flow velocity magnitude with streamlines and the pressure distribution, based on equations (4.27) -(5.5), and utilizing the parameters in Table 1. The expected characteristic pressure is p * = O(1KPa), the characteristic flow velocity within the chamber is v * = O(1mm/s) and Re ≈ 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.

Additional dynamic cases
The theoretical solution also readily allows analyzing 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 behavior 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.

Case II -Dictated inlet pressure and hyper-elastic 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, λ(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, Similarly to the previous case, we substitute (R, θ) = (λ, π) in the pressure distribution solution (4.28c) 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 n (−1) n − χ n P n (cos θ) . (5.9) Finally, using the mechanical energy principle (5.4), we derive a nonlinear ordinary differential equation governing the time-dependent stretch of the chamber, The function Υ (ǫ) was estimated numerically by summing the first 10 5 terms in the series for several values ofǫ corresponding toǫ ∈ [0.002, 0.5]. By substituting the asymptotic approximation of Υ , into the differential equation that governs the chamber's stretch (5.10), we obtain its approximated explicit form given by Next, we investigate the system's behavior 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 η ef f (t) = A S (t)/4π. Here, A S (t) is the body's surface area obtained in the nonspherical simulations, and η ef f (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.
In order to approximate the characteristic time constant of the system, enabling us to estimate the period 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 (λ SS , P SS ), when the pressure at the inlet is dictated to be P ext (τ ). The linear equation is solved analytically, leading to the following solution: SSˆτ 0 e βI ·τ ′ △P in (τ ′ )dτ ′ e −βI ·τ . (5.13) where △λ L = λ L (τ ) − λ SS is a small stretch variation around its nominal value λ L , △P in (τ ) = P in (τ ) − P SS is the pressure variation from its nominal value, and (5.14) From (5.14), it is clear that the solution branches of P SS in (3.3) are stable equilibria if and only if dP SS /dλ 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/β I . Since the derivative dP SS /dλ 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, equation (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, △P in (τ ) = 0, the solution of this equation under small initial perturbation from equilibrium λ Q (T ) is given by: (5.16) In Figure (9), the linear and the second-order approximations were displayed alongside the exact numerical solution of equation (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.
In section 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 λ I SS 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). This 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 ψ(λ), based on the chosen constitutive law. Here, we have chosen to present the results using Mooney-Rivlin's hyperelastic constitutive law to illustrate the phenomenon of bi-stability 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.

The dynamic behavior of two interconnected bi-stable chambers
In this section, we analyze the behavior 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 ) ≡ Q(T ). This system is instrumental for understanding the behavior of interconnected bi-stable 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 bi-stable elastic systems have been examined, the main assumption is uniform pressure within the elastic element (Ben-Haim et al. 2020;Dreyer et al. 1982;Treloar 1975;Glozman et al. 2010). Here, we shall analyze 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 analyzed 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 λ 1 (T ), λ 2 (T ), respectively. From the integral mass rate The stretch λ 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 (ef f ) in (T ) = P (λ 1 , 0; T ) is the effective external pressure, where P (R, θ; T ) is given by (5.5). Recalling that the Fourier coefficients Λ n and ϕ n are linear combinations of the inlet the outlet flow rates (4.26), the effective external pressure is rewritten as andǫ 1 (T ) = ǫ/λ 1 (T ). Utilizing the results achieved for the different values ofǫ 1 yields the following approximations Π (out) (ǫ 1 ) = O(ǫ 1 ) and Π (in) (ǫ 1 ) = O(ǫ 2 1 ). Using the governing equation of λ 2 (T ) in (5.12) yields the second equation of motion, (6.4) Figure 11: Series solution of a system consisting of two chambers controlled by single inlet: Left -velocity field and stream-lines, Right -pressure distribution. The flow velocity and pressure distribution corresponding to when the system passes through the point A shown in Figure 10.
Equations (6.1) and (6.4) are a set of nonlinear coupled first-order differential equations that govern the evolution of the chambers' stretches λ 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) 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 λ 1,SS = λ 2,SS , where the radii of both chambers are equal, whereas in the second one the radii are different, λ 1,SS = λ 2,SS , thanks to the bi-stability 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, λ 3 1 + λ 3 2 = 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 quasistatic; 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 10(c). Figure 10(b) shows the system's trajectory in {λ 2 , λ 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 {λ 2 , λ 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.

Concluding Remarks
In this work, we analyzed the dynamics of creeping flow in a bi-stable hyperelastic spherical chamber by calculating the velocity field and pressure distribution inside the chamber. The analytical results were compared to numerical simulations, which showed an excellent fit. From the normalization of the governing equations (4.1) and (4.2), we obtained the condition ε = a 4 /r 3 0 ℓ ≪ 1 which led to creeping flow inside the chambers. Moreover, we obtained the condition q * ≪ µa/ρ 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 non-linear constitutive elastic laws (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 ψ * /t 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 (χ → 1), the viscous flow will generate O(ǫ 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 θ away from the poles θ = 0 and θ = π. 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 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 bi-stability feature through which the chambers can display controlled transitions between different multi-stable 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 modeling 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 analyze the dynamics of particles inside spherical elastic chambers in low-Reynolds flow as future work. equilibria exist only when 0 < α < 0.214. Those parameters are considered to lie within the physical domain of rubber-like materials, as well as for biological tissues. The next curve we examine is the solution of the equilibrium equation (3.3). Since the equation has no analytical solution, we shall find an asymptotic approximation by separating the solution into three main regions, as shown in Figure 2(c). In the first region, (I), the pressure in the chamber is lower than the maximum point, 0 < P SS < P A and 1 < λ SS < λ A . In the second region, (II), the pressure in the chamber is between the minimum and maximum points of the graph, P B < P SS < P A and λ A < λ SS < λ B . In the third region, (III), the pressure in the chamber is higher than the minimum point on the graph, P SS > P B and λ B < λ SS . In the first case (where P < P A ), the equilibrium point will be close to λ SS = 1. Ilssar & Gat (2020) present the following approximation to describe this branch by λ SS = 1 + δ 1 + δ 2 1 + O(δ 3 1 ), where δ 1 is formulated by, δ 1 (P SS ; α) = − 7P SS + 24 + 24α − 3 64αP SS + 192(α + 1) 2 − 21P 2 SS 8(7P SS − 33α − 21) .
(A 8) Next, the second stable equilibrium radius, which exists when the pressure is higher than the local minimum (P > P B ), is considered significantly larger than unity. Thus, in order to formulate an approximation for this equilibrium stretch, λ SS = δ −1 2 (α) ≫ 1 is substituted into (3.3). After some regular algebraic manipulation, a second-order algebraic equation in terms of λ SS is provided. The solution is given by δ 2 (P SS ; α) = P SS ± P 2 SS − 64α 8 .
The solution having a positive sign in front of the square root of the discriminant, is suitable for the solution of the third case in which P SS > P B , while the solution with the negative sign is suitable for solution of the unstable branch in which P B < P SS < P A . Those approximations are also plotted in Figure 2(c) with dash-lines on the solid curve representing the exact solution.
Appendix B. Asymptotic justification for neglecting the non-spherical deformation of 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, εP 1 (R, θ;ǫ; T ). This result was obtained analytically by (4.12) without the assumption of spherical deformation. Assuming that the chamber's radius is large relative to the tube's radius (ǫ a ≪ 1), and based on the coupled model of fluid-solid numerical simulation results (see Figure. 4), it is clear that the maximum pressure gradient is obtained in the inlet (or outlet) region and significantly decayed within the chamber's area. Since far away from the tube (χ → 1), the pressure is approximately hydrostatic, the real shape of the chamber will consist of an ideal sphere in addition to a small perturbation. Based on the numerical results we have obtained in § 5.2, the perturbation of the chamber's stretch is O(ǫ t ), so we shall assume a regular asymptotic approximation as follow, λ(θ; T ) = λ 0 (T ) + ǫ t λ 1 (θ; T ) + O(ǫ 2 t ).
(B 1) As we have already seen in section § 4, 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. According to the solution obtained by Beatty (1987), P S (λ 0 ) = λ −2 0 · dψ/dλ 0 is the well-known isotropic pressure which represents the hydro-static pressure in the leading order. On the other hand, from the CFD simulation, we saw that close to the chamber's wall (χ → 1), the viscous flow will generate O(ǫ t ) spatially varying corrections. Therefore, it is convenient to assume that pressure distribution can also be approximated as P (R, θ; T ) ≈ P S (λ 0 ) + ǫ tP (R, θ; T ) at χ → 1, (B 2) whereP ∼ O(1) and ǫ tP ∼ εP 1 close to the inner chamber's wall. In order to examine the effect of non-spherical deformation on the resulting pressure profile (without considering nonlinear hyperelastic analysis that is outside the scope of this work), we wish to use the regular asymptotic approximation (B 2) and set the radial perturbation (B 1), P (R = λ(θ; T ), θ; T ) ≈ P S (λ 0 ) + ǫ tP (R = λ(θ; T ), θ; T ) at χ → 1, (B 3) by taking the Taylor approximate, we get P (R = λ(θ; T ), θ; T ) = P S (λ 0 ) + ǫ tP + O(ǫ 2 t ) ∼ P S + εP 1 ; (B 4) hence, the O(ǫ t ) in the chamber's shape create O(ǫ 2 t ) pressure correction. In order to find the full fluid-structure interaction, the entire elastic equations must be solved coupled with the Stokes equation, but this analysis will be performed as future work.