The inflation of viscoelastic balloons and hollow viscera

For the first time, the problem of the inflation of a nonlinear viscoelastic thick-walled spherical shell is considered. Specifically, the wall has quasilinear viscoelastic constitutive behaviour, which is of fundamental importance in a wide range of applications, particularly in the context of biological systems such as hollow viscera, including the lungs and bladder. Experiments are performed to demonstrate the efficacy of the model in fitting relaxation tests associated with the volumetric inflation of murine bladders. While the associated nonlinear elastic problem of inflation of a balloon has been studied extensively, there is a paucity of studies considering the equivalent nonlinear viscoelastic case. We show that, in contrast to the elastic scenario, the peak pressure associated with the inflation of a neo-Hookean balloon is not independent of the shear modulus of the medium. Moreover, a novel numerical technique is described in order to solve the nonlinear Volterra integral equation in space and time originating from the fundamental problem of inflation and deflation of a thick-walled nonlinear viscoelastic shell under imposed pressure.

For the first time, the problem of the inflation of a nonlinear viscoelastic thick-walled spherical shell is considered. Specifically, the wall has quasilinear viscoelastic constitutive behaviour, which is of fundamental importance in a wide range of applications, particularly in the context of biological systems such as hollow viscera, including the lungs and bladder. Experiments are performed to demonstrate the efficacy of the model in fitting relaxation tests associated with the volumetric inflation of murine bladders. While the associated nonlinear elastic problem of inflation of a balloon has been studied extensively, there is a paucity of studies considering the equivalent nonlinear viscoelastic case. We show that, in contrast to the elastic scenario, the peak pressure associated with the inflation of a neo-Hookean balloon is not independent of the shear modulus of the medium. Moreover, a novel numerical technique is described in order to solve the 2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Understanding the mechanics of the classical problem of the inflation of a balloon requires the application of large deformation nonlinear elasticity theory, see, for example, Sec. 5 of [1], which refers back to the experiments of [2,3]. Of specific interest in the balloon inflation experiment is the resulting non-monotonicity of the pressure-stretch curve (the well-known and frequently experienced difficulty with the initial inflation of a balloon) and the associated critical stretch at which the pressure inside the balloon reaches an initial maximum followed by the subsequent more straightforward inflation after this initial maximum [4]. This initial maximal pressure condition gives rise to a non-trivial problem known in the literature as a limit-point instability [2,[5][6][7][8]. Furthermore, rather interestingly for balloons with neo-Hookean constitutive response, the critical stretch at which the pressure reaches the initial maximum is independent of the shear modulus of the medium, thereby exhibiting universal behaviour.
Understanding the balloon inflation and deflation problem has important applications. Recently, the US Food and Drug Administration (FDA), in order to treat obesity without the need for invasive surgery, approved [9,10] a silicone balloon device called a bioenteric intragastric balloon (BIB) (figure 1c). The BIB is delivered into the stomach via the mouth through a minimally invasive endoscopic procedure and inflated successively by fluid/air injection in order to take up space in the stomach and to induce early satiety. Recently, the FDA received several reports regarding some adverse events in patients with recently fitted liquid-filled BIBs of certain types, although it is not entirely clear what caused the complications in the patients [11][12][13]. One would suggest, however, that it is crucial to fully understand the inflation/deflation mechanism of BIBs [14]. This specific medical scenario suggests that a constitutive model devised for intragastric balloons should incorporate large deformations as well as viscous effects. The same considerations can be extended to atmospheric balloons (figure 1d), which are generally made of thin polymeric film and illustrate considerable viscoelastic behaviour [15,16]. Of potential relevance, therefore, is the study of the canonical viscoelastic balloon inflation problem.
The balloon problem is, of course, a very specific limiting case of the more general hollow thick-walled spherical annulus and this geometry arises as a simplification of many biological scenarios, e.g. in the context of hollow viscera, hence the dual study of the inflation of balloons and mammalian bladders in the paper from 1909 by Osborne & Sutherland [17]. In that paper, it was observed that there was a significant difference between these two scenarios and specifically it was noted that in the case of the bladder, in contrast to the balloon, no such initial local maximum pressure occurs; the curve is monotonic. Recent work by Mangan and Destrade has considered appropriate strain energy functions to model both thin-and thick-walled cases [18].
In biological applications, it should be noted that the materials in question are often strongly viscoelastic in the large deformation regime (indeed, viscoelastic hysteresis was noted by Osborne and Sutherland in their experiments). These include the inflation of organs such as the lungs, bladder, colon and also arteries and veins, noting that the latter three examples exhibit cylindrical rather than spherical symmetry of course. Of specific interest, in these cases, are the hysteretic pressure-volume curves that arise due to imposed volumetric changes or due to imposed pressures. Understanding the mechanical properties of living soft tissue is important especially in the context of diseased tissue, in order to prevent, for example, uncontrollable wall dilatation which can lead to aneurysm and to a surrounding wall rupture [19,20].
In the case of the bladder (figure 1b), inflation and deflation (or filling and voiding as they are often termed) are complex physiological processes driven by pressure differences inside and outside the bladder. Voiding is driven by the contraction of the detrusor muscle, which leads to an initially significant increase in the internal pressure, followed by a decrease upon the release of  urine from the bladder [21,22]. The spherical approximation during this process has been shown to be a fairly reasonable one [23,24]. In a related field, and similar geometry, it is well known that the pufferfish (figure 1a)u s e s self-inflation via water intake as a defence mechanism against predation. During inflation, thin rigid spines that are initially laterally located within the exterior dermis rotate orthogonally to the surface of the pufferfish, giving rise to an inflated spiky sphere [25,26]. This response means that they are much harder to swallow by would-be predators. Although, once again, the initial geometry is more complicated than a hollow sphere and has also been associated with an inhomogeneous response [27], this example serves to illustrate that such configurations are important in a broader biological context. It is also worth stressing here that although there have been numerous studies of the inflation mechanism and associated skin tissue response to loading, there appear to be no studies of deflation, during which viscoelastic effects will play an important role.
In a rather different limiting scenario, when the outer radius of the hollow sphere tends to infinity, the canonical hollow sphere problem becomes highly relevant to understanding the response of porous viscoelastic elastomers under hydrostatic pressure [28][29][30][31]. In closed-cell materials, voids are distributed throughout an elastomer, resulting in a compressible medium even when the host elastomer is incompressible. When the voids are in a dilute distribution, a first approximation is to neglect interaction between voids and determine the deformation of a single void under compression (see [28] such a medium, it will 'see' only the coupled effect of rigid body displacement and hydrostatic pressure. Essential to effective constitutive models of such porous viscoelastomers, therefore, is the ability to predict the response of a single void in a nonlinear viscoelastic medium under hydrostatic pressure. Such porous elastomers are of interest in a number of applications including protective materials, clothing and textiles and underwater acoustics, to name but a few. Recently, such materials and related structures have become the focus of the metamaterials community since they exhibit slow sound behaviour and associated strong resonances [32][33][34].
The importance of viscoelastic effects in large deformation elasticity has been stressed in a variety of recent studies, in particular those associated with soft biological tissues [35][36][37][38][39], polymers and rubbers [40,41]. Under very specific loading conditions, the viscoelastic nature of such materials can be neglected to a first approximation and the theory of nonlinear elasticity can be employed, where, under the hyperelasticity assumption, strain energy functions W can accommodate a variety of constitutive behaviours. In reality however, at finite-deformation rates, such approximations cannot be made, and a nonlinear theory of viscoelasticity is required.
Motivated by the plethora of application areas described above, here we study the canonical problem of the large radial viscoelastic deformation of a hollow, thick-walled sphere with initial inner and outer radii A and B, respectively, subjected to hydrostatic pressure both internally and externally by a quasi-static load. Inertia is therefore neglected. The limits A → B and B →∞ are therefore associated with the balloon and isolated void respectively. A fundamental analysis of such a problem in the nonlinear finite-deformation viscoelastic regime is currently lacking and therefore there is a paucity of models to understand the pressure-deformation curves associated with viscoelastic balloons, thick-walled viscoelastic shells, hollow viscera or isolated voids in porous viscoelastomers. One exception is the specific study of Wineman [42], who considered a rather special constitutive law associated with a spherical membrane. In particular, in a certain parameter regime, corresponding to α → 0 of that model, the medium becomes perfectly hyperelastic. Wineman studied the limit point instability associated with this material behaviour. A few years later, Calderer [43] considered the specific scenario of the radial elastic motion of a thick spherical shell under a constant pressure difference between the inner and outer surfaces. This is a particular case of the general (non-spherical) compressible medium problem studied by Ball [44], who showed that, for suitable initial conditions and pressures, no weak solution exists for all time, suggesting that the displacement becomes singular within a finite time. Later, Calderer considered the equivalent viscoelastic problem of a spherical shell [45,46], examining whether the presence of dissipation in the differential equation describing the deformation is able to prevent instabilities and ensure global existence of the solutions in time.
In order to account for material nonlinear viscoelasticity, here we assume a so-called quasilinear large deformation viscoelastic response [47]. Assuming spherical symmetry, two types of problems can then be studied: inflation/deflation due to either imposed strains (volumetric inflation) or imposed pressures, the former type being more straightforward to solve. Initially, then, here we consider two problems where volumetric inflation is assumed, meaning that the pressure difference and deformed radii can be determined straightforwardly from a single integral expression. First, we consider the fitting of the quasilinear viscoelastic (QLV) model to experimental data associated with volumetric inflation of murine bladders. Next, the canonical problem of the inflation of a balloon via volumetric control is described and compared, in particular, with the well-studied associated nonlinear elastic problem. To make progress in the more general framework of the inflation of a thick-walled viscoelastic spherical shell when pressure is imposed, the analysis is then confined to the case when the so-called instantaneous elastic stress is governed by a Mooney-Rivlin strain energy function, allowing study of, for example, a viscoelastic neo-Hookean response. More general constitutive models can be considered, but this case illustrates the modification to the nonlinear elastic result and is indicative of the fundamental differences between elasticity and viscoelasticity. More complex models, including, for example, anisotropy, can also be developed but, for the present purpose, they obscure the main message of this article, which is associated with the influence of viscoelasticity in nonlinear media.
The review paper by Wineman [48] has described a variety of possible viscoelastic constitutive models that incorporate finite deformation. A rather general model is that of Pipkin & Rogers [49], which has the advantage of allowing for strong nonlinearity and finite deformation. This model also incorporates coupling between relaxation and strain, where necessary, giving rise to models for strain-dependent relaxation. In order to simplify the approach and generate models that are more tractable, a simplified version of the Pipkin-Rogers constitutive model was suggested by Fung [50], and this approach is now known as QLV theory. This assumes that viscous relaxation rates are independent of the instantaneous local strain. Although this method has been criticized in recent years, a recent paper has shown that, although the model clearly cannot be valid for all materials over all deformation rates and strains, the criticisms put forward by previous authors were unfounded [47]. For example, in a number of publications, an incorrect QLV relation or stress measure was employed (the latter must be the second Piola-Kirchhoff stress to satisfy objectivity), or the incompressible limiting form was derived erroneously. QLV does, of course, have limitations; the fact that the relaxation functions are independent of strain is an important, possibly erroneous, assumption. However, the model appears to include enough detail to capture many of the essential elements of the physics while not being overly difficult to implement in the context of real-world applications.
Note that for an imposed deformation it is generally simple to derive the stress field since the stress is prescribed in an integral form with the deformation appearing in the integrand. On the other hand, when traction is imposed, the problem is more complicated since integral equations are generally required to be solved for the resulting deformation. Classical discretization methods used in the past [48,51,52] have been employed recently (with adapted and improved techniques) in order to solve such integral equations for homogeneous deformations [47,53,54]. These approaches permit wide utility of the models without the need for finite-element implementation. Nevertheless, inhomogeneous deformations, to our knowledge, have not been considered in this context, presumably because of the added complexity that arises from the presence of the spatial variation in the deformation. Even in the most simple cases, inhomogeneous deformations will lead to integral equations that are much more difficult to solve than those that arise for homogeneous deformations. However, every deformation that is controllable (i.e. a deformation which satisfies the balance equations of equilibrium with zero body force, supported by suitable surface tractions only) for homogeneous isotropic incompressible elastic materials, in the absence of body forces, is also a controllable deformation for a more general class of homogeneous isotropic incompressible materials known as simple materials (see [55][56][57]). We note that simple materials are defined as materials within which the stress at each point is determined by the histories of the deformation gradients at that point, and which therefore includes all materials with memory, i.e. viscoelastic materials.
Although inhomogeneous, the purely radial nonlinear elastic deformation associated with the inflation and deflation of an incompressible hollow sphere subjected to hydrostatic pressure is a universal solution belonging to 'Family 4' (see [58,59]). The problem is straightforward and its solution has long been known (see, for example, [28,58] and references therein). The equilibrium equation can be integrated exactly to yield a nonlinear equation for the internal deformed void radius a (or equivalently for the external radius b) in terms of the imposed pressure difference and the undeformed radius A (or the external undeformed radius B, respectively), in the form Here T rr and T θθ are the radial and tangential Cauchy stresses and we specify T rr =−p a on r = a and T rr =−p b on r = b (or when the outer hydrostatic pressure is applied in the far-field, i.e. B, b →∞,o nr =∞)( figure 2). Given a strain energy function, the integral on the right-hand side in (1.1) is determined in a straightforward manner and the deformed radius a is determined numerically [28]. The corresponding linear viscoelasticity problem is also straightforward, and since the associated constitutive law in that case may be inverted without difficulty, both imposed displacement and pressure conditions can be derived. This paper will focus therefore on the problem when the medium is QLV. The paper proceeds as follows. In §2, a summary of the equations governing the deformation of an incompressible isotropic QLV medium is presented for an arbitrary strain energy function. We then sequentially study a series of problems beginning with imposed volumetric strain and ending with imposed pressure. In §3, an assessment of the efficacy of the model is presented by fitting experimental data associated with the inflation of murine bladders when volumetric strains are imposed. The strongly related canonical thin-walled (balloon) inflation/deflation problem, when controlled by volumetric strain, is solved in §4, noting the additional effects that arise due to viscous effects over and above the classical nonlinear elastic balloon inflation problem. In §5, the general formulation is then restricted to the case of a neo-Hookean strain energy function in order to study the more difficult problem of inflation due to imposed pressure, with reference also being made to the more general Mooney-Rivlin strain energy function. The governing equation is a nonlinear Volterra integral equation and we describe a new root finding technique to determine the resulting radial stretch when pressure is prescribed. Some canonical problems associated with the inflation and deflation of a shell of finite thickness are subsequently solved using this procedure in §6. Concluding remarks are made in §7.

Governing equations
With reference to figure 2, consider a hollow sphere of inner and outer radii A and B, respectively, which is capable of large deformation and whose QLV constitutive behaviour will be specified shortly. The medium is assumed to be incompressible and isotropic. We consider the deformation of the medium when it is subjected to hydrostatic internal (p a (t)) and/or external (p b (t)) pressure in a time-dependent manner beginning at some reference time t 0 = 0. Similarly, we can consider the problem of a spherical void in an infinite host medium if we impose the external pressure on b, and take the limit B →∞(and hence b →∞).
We assume that the centre of the hollow sphere is located at the origin of a Cartesian coordinate system (X, Y, Z). The deformation of the void will be purely radial and therefore, working in spherical polar coordinates, we can represent this deformation in the form (

2.2)
Here I is the second-order identity tensor, x(s) = (x(s), y(s), z(s)) is the position vector of a generic particle P at time s ∈ (−∞, t]a n dX = (X, Y, Z) is its position at the reference time s = t 0 = 0. The quantity B ≡ FF T is the left Cauchy-Green strain tensor. Associated with this for an isotropic medium are the three principal strain invariants Since we assume that the motion begins at s = t 0 = 0, the polar components of the deformation gradient (2.2) associated with (2.1) are therefore given by For incompressible materials, the constraint of incompressibility Then Assuming next that an underlying elastic stress associated with the host can be given in terms of the derivative of an elastic potential W(I 1 , I 2 ), which is a so-called hyperelastic strain energy function, the constitutive law for the elastic Cauchy stress is where p is the Lagrange multiplier introduced by the incompressibility constraint and W i ≡ ∂W/∂I i . From the deformation gradient (2.4), the viscoelastic Cauchy stress for a QLV material can be written as [47] is associated with the deviatoric part of the Cauchy stress, C = F T F is the right Cauchy-Green strain tensor and D is the scalar relaxation function associated with the time-dependent response of the QLV material. Thus from (2.9) and using (2.4) the principal Cauchy stresses are and where from (2.10) and (2.7) it is straightforward to show that (2.14) We assume that the deformation acts on time scales such that inertia can be neglected and therefore the equation of motion reduces to div T = 0, (2.15) where body forces have also been neglected. Boundary conditions are The only non-trivial equation of (2.15) is the radial one, which when integrated with respect to r and imposing the boundary conditions (2.16) reduces to give (1.1). Writing the right-hand side in terms of the radial coordinate R in the reference configuration, this equation reduces to the form dR.
( 2.17) From (2.11) to (2.14) and (2.7) we have noting that the functions f 1 and f 2 are explicitly independent of the strain energy W. Two different types of problems then arise. The first, more straightforward type is to choose the strain energy function W, the shear modulus µ and relaxation function D(t) together with the initial internal and external radii A and B and impose the deformed radii a and b,i . e .we impose volumetric strain. We can then straightforwardly determine the pressure difference p b (t) − p a (t) as a function of time by determining the integral on the right-hand side of (2.17). The second type is a more complicated mathematical problem. In this case, our task is to determine α(t) from (2.17) when we impose W, µ, D, A and B as before but now where we impose pressures p a (t) and p b (t). From α(t) we can, of course, then determine a(t)a n db(t), the evolving internal and external boundaries of the hollow sphere. The case of imposed pressure is non-trivial because (2.18) is a nonlinear Volterra integral equation in space and time. In previous works, these have been solved by a numerical procedure that exploited the separable nature of the terms under the integral in terms of s and t [47,53]. Here, however, this separability does not arise and a new, modified approach is required. In §5, we will develop this approach for the commonly considered simple but illustrative and informative case when W is the neo-Hookean strain energy function [60]. Before this however we consider some problems when deformations are imposed and pressure differences are measured. We begin by illustrating that the QLV model can qualitatively model data obtained from experiments performed on the controlled filling of murine bladders. Furthermore, following this, the canonical problem of the inflation and deflation of a thin-walled viscoelastic balloon is studied, the elastic equivalent of which has been studied extensively in the past.

Quasilinear viscoelasticity model it to experimental investigations on urinary bladder illing (a) Animals
All experiments were performed using 12-16-week-old C57/BL6 male mice from Charles River (Margate, UK). Mice were acclimatized for 7 days in the laboratory animal husbandry unit under a 12 L:12 D cycle and had free access to water and food. Immediately before the onset of the protocol, mice were anaesthetized with isoflurane and humanely sacrificed by cervical dislocation according to UK Home Office legislation regulating Schedule 1 procedures (Scientific Procedure Act 1986). Ethics approval was obtained from the University of Sheffield Ethics Review Panel.

(b) Bladder preparation and isovolumetric experiments
Following euthanasia, the whole pelvic region including surrounding tissues was dissected from the animal and placed in a purpose-built recording chamber.

(c) Experimental protocol
Control distensions were always carried out at the start of the protocol. To do this, bladders were distended using isotonic saline (NaCl, 0.9%) at a rate of 100 µlmin −1 to a maximum pressure of 40 mmHg; at this point the infusion pump was stopped and rapid evacuation of the fluid occurred by opening the two-way catheter in the dome to the atmosphere. This was repeated several times at intervals of 10 min to assess the viability of the preparation and reproducibility of the intravesical pressure responses to distension. Once stable and reproducible responses were obtained (typically after three control distensions) the isovolumetric experiments were conducted.
To perform the isovolumetric experiments, bladders were filled to an intravesical pressure of 25 mmHg. At this point, the tap on the dome catheter was closed and the infusion pump was switched off. Bladders were left to equilibrate for 15-20 min at the intravesical volume, following which the bladder was emptied by opening the tap. Intravesical pressure data were continuously collected throughout the protocol and data points at 5 s intervals were analysed.

(d) Model it
We now fit the experimental data obtained using the procedure described above with a viscoelastic model using three methods: (i) QLV with an isotropic neo-Hookean strain energy function, (ii) QLV with an isotropic Fung strain energy function, and (iii) using linear viscoelasticity. As mentioned above, the bladders were filled at a constant rate of 100 µlmin −1 until a critical time t c , after which the volume was held constant. Therefore, an equation for the where the volume V is given in microlitres, the inner radius A is given in metres and the time t and critical time t c are given in seconds. The volume is related to the deformation parameter α(t) via the following equation: which is what we impose to drive the deformation for the first two methods. We use equations (2.17)-(2.19) to determine the pressure difference by substituting the following strain energy functions into the first equation of (2.19). The neo-Hookean strain energy function is given by where µ is the ground-state shear modulus of the material under consideration, which gives W 1 = µ/2, W 2 = 0. The isotropic Fung strain energy function is given by which gives W 1 = (a + bce c(I 1 −3) )/2, W 2 = 0. In order to reduce the number of fitting parameters for this model, we assume that a = 0, which means that, for consistency with linear elasticity, we must have bc = µ, and therefore W 1 = µe c( For the third method, we use linear viscoelasticity. The displacement vector u = u(r, t)e r for an incompressible linear viscoelastic material must satisfy div u = 0, and therefore we have The deformation of the inner wall is prescribed by V(t) using the following equation: The stress in an incompressible linear viscoelastic material is given by where p is a Lagrange multiplier associated with the incompressibility constraint and e is the linear strain tensor, which for our deformation is given by e = diag(du/dr, u/r, u/r). Using the above, we obtain which can be substituted into the equation below to obtain the pressure difference, For all three methods, we assume the non-dimensional relaxation function to be a one-term Prony series of the form where µ ∞ is the long-time shear modulus and τ is the relaxation time. We fit for the pressure difference p a (t) − p b (t) (note the sign change with respect to equations (2.17) and (3.9) due to the fact that we are considering a situation in which the internal pressure is greater than the external  pressure) using the three methods described above by using the parameters µ, µ ∞ , τ and c as fitting parameters. We assume that the bladder had an initial inner radius of 2 mm and outer radius of 2.3 mm and that the critical time that marked the cut-off of the filling phase was t c = 155 s. The results of the fitting process are plotted in figure 3. The predicted parameter values and the mean squared error associated with each fit are reported in the table below. The ability to more accurately capture the loading phase when using the Fung QLV model leads to a much smaller mean squared error than in the other two cases. Having illustrated that there is merit to the model, let us now consider two canonical problems in the sections to follow. The nonlinear elastic equivalents of each of the following problems are classical but the influence of viscoelastic effects has not yet been studied.

Inlation and delation of quasilinear viscoelastic balloons with imposed volumetric strain
We shall now consider a problem whose deformation is also volume controlled, but which is more canonical in its nature than the specific case of a bladder as described in the previous section. In particular, this problem has an equivalent well-studied analogy in nonlinear elasticity so that here we describe the additional effect of viscosity. This problem is the thin-walled limit of the hollow sphere. This is of interest in many scenarios, including biological membranes [61], and of course is strongly related to the bladder application considered in the previous section. The problem is also an extension of the classical elastic balloon inflation problem considered by numerous authors as described in the Introduction. The problem dictates that finite-deformation (visco)elasticity theory should be employed, but as discussed in the Introduction, in the purely elastic case with a neo-Hookean strain energy, the pressure-volume curve is non-monotonic [18] and further the peak pressure is independent of the shear modulus of the medium. The departure from neo-Hookean behaviour and the influence of the strain energy function has been studied extensively [ However, there is a lack of results on the subsequent behaviour when effects beyond elasticity are incorporated, e.g. viscoelasticity, with an exception associated with the study of Wineman [42,52], who considered a rather special viscoelastic constitutive law for the large deformation of a spherical membrane. Wineman [42] derived a necessary condition in order to guarantee the existence of a limit-point instability. However, in the limit α → 0 of that model, corresponding to a neo-Hookean material, the medium becomes perfectly hyperelastic, having a peak-pressure result that is consistent with the standard neo-Hookean universality result [4,18,62]. Let us consider the response in the context of a Mooney-Rivlin strain energy function where γ is a constant in the range −1/2 ≤ γ ≤ 1/2. The quasilinear constitutive law proposed here however remains viscoelastic regardless of the choice of the parameter γ and therefore we are able to study the influence of viscoelastic effects on a medium with a strain energy function that is associated with neo-Hookean behaviour (i.e. when γ = 1/2). As described in the Introduction, although this balloon inflation problem is canonical and therefore well studied in the elastic case, viscoelastic effects could potentially play an important role in, say, the context of the inflation and deflation of intragastric balloons, which represent one of the most frequently employed medical devices to treat obesity. The efficacy of intragastric balloons may be influenced by volume control [10][11][12][13] and creep, and these processes must be completely understood in order to better predict inflation/deflation behaviour. Related to this aim then, the canonical study of the inflation and deflation mechanisms of quasilinear viscoelastomeric balloons is presented here. By varying the elastic modulus µ ∞ and the relaxation time τ appearing in the relaxation function (3.10), we show the loss of neo-Hookean peak pressure universality, emphasizing the relevance of the pressure/volume control for materials that exhibit inelastic effects.
In order to present some specific results employing (2.17) to obtain the pressure difference, let us set the relaxation function to be the classical one-term Prony series as in (3.10). For convenience in this canonical problem, let us introduce the following dimensionless quantities: such that 1 <R <B and τ r is a reference relaxation time which we set as τ r = 1.0s. Let us fix ∆ = (B − 1) ≪ 1 in order to recover the theory of spherical membranes. To this extent from now on we set ∆ = 10 −6 , and we denote the critical stretch at which peak pressure occurs asã * . We consider inflation prescribed byã =ã(t) with p b = 0, as shown in figure 4a, with the resulting ∆p −ã curves presented in figure 4b for the cases of neo-Hookean (γ = 1/2) and Mooney-Rivlin (γ = 1/4, 1/8) strain energy functions and in the case whenμ ∞ = 0.5,τ = 0.1. Also note here that deflation curves can present a lower local maximum.
For the imposed inflation/deflation cycle depicted in figure 4a we now determine the pressure-stretch curves for balloons whose instantaneous elastic stress is derived from the neo-Hookean strain energy function when the non-dimensional relaxation timeτ varies withμ ∞ fixed (figure 5a) and alternatively when the non-dimensional long-time shear modulus parameterμ ∞ varies withτ fixed (figure 5c). This clearly illustrates the loss of universality of peak pressure once viscoelastic effects are taken into account, noting that for a purely elastic, neo-Hookean balloon this peak pressure occurs atã * =

Mathematical formulation and numerical scheme for inlation and delation due to imposed pressure
Let us now move on to the more complicated mathematical problem of imposing pressure and determining the deformed radii. Initially, we describe the general formulation of the problem in the context of the Mooney-Rivlin strain energy function (4.1) before focusing the numerical scheme on the case of a neo-Hookean material.

(a) Reduction of the governing equation
In the case when the strain energy function is of Mooney-Rivlin type (4.1), the expression for g in (2.19) 1 reduces to Using this in (2.18) and subsequently in (2.17) and interchanging the order of integration, (2.17) can be rewritten as where the functions J 1 (α(t)) and J 2 (α(s), α(t)) can be determined explicitly for the Mooney-Rivlin medium and are stated in appendix A. They take on a particularly simple form for the case of a neo-Hookean medium γ = 1/2, since then g(R, t) = µ in (5.1) reduces to a constant. Therefore here, in order to illustrate the method, we focus on the neo-Hookean strain energy function, where W does not depend on the second invariant I 2 . The procedure that we describe below can be modified to incorporate the more general Mooney-Rivlin form, or indeed other strain energy functions, but this specific, simple case is sufficient to illustrate an array of interesting behaviours.
In the case of a neo-Hookean medium then and and we stress here the dependence on α(s)a n dα(t), recalling that r 3 (R, s) = R 3 − α(s). It is straightforward to calculate the term J 1 (t) explicitly Recall that a limit of interest is B →∞, corresponding to a void in an unbounded medium and in this case (5.5) reduces to The function J 2 (α(s), α(t)) requires more care as we now describe. Start by defining the new variable and so du = (3/R)(1 − u)dR and Therefore, Next, letting β = α(s)/α(t) one can show that This is integrated to yield the explicit form Two important special cases of J 2 are when β = 0andβ = 1. In the former, i.e. when α(s) = 0, we have J 2 = 0, whereas in the latter, when s = t,wegetJ 2 (α(t), α(t)) ≡ J 1 (α(t)). Note furthermore that (5.11) is continuous at α(t) = 0 and it converges to and also J 2 (0, 0) = 0. For practical computations, when the strain is imposed the explicit expression (5.13) is very useful in order to ensure convergence of the scheme (if this is used even in the simplest strain imposed case). The above analysis can be extended in a straightforward fashion to the case of a Mooney-Rivlin medium, which yields the same expression (5.2) but now with alternative forms for J 1 and J 2 . Both of these are stated in appendix A. The function J 2 is evaluated explicitly but only in terms of a certain hypergeometric series.

(b) Numerical scheme for imposed pressure
Once the relaxation function D, material properties and relaxation times are given, for an imposed α, (5.2) predicts explicitly the applied pressure difference across the boundaries r = a, b required in order to maintain the given deformation. By contrast, for an imposed hydrostatic loading, (5.2) is a nonlinear (space dependent) Volterra integral equation, the solutions of which are non-trivial to determine. To solve the equation for a given time-dependent pressure difference, we discretize the integral and then pose the problem in terms of a root-finding scheme for each time step. To proceed, define the function and we seek values of α(t) such that R(t) = 0. At t = 0, clearly α(0) = 0. At the next time step t = ∆t ≪ 1, we have (p a (∆t) − p b (∆t)) + J 1 (α(∆t)) + ∆t 0 D ′ (∆t − s)J 2 (α(s), α(∆t)) ds. (5.15) Then anticipating what is to come, let R n = R(n∆t), α(n∆t) = α n and the imposed (scaled) pressure difference Given α 0 = 0, we need to determine α 1 such that R 1 = 0. Use the trapezium approximation to the integral to yield where we have introduced the notation D ′ n = D ′ (n∆t). More generally, keeping track of each α m , m < n we have Recalling that α 0 = 0, we use the result that J 2 = 0whenα(s) = 0 = α 0 and also employ (5.5) to find that Each successive α n is determined by ensuring that R n = 0. Note that this scheme can be employed regardless of the form of the functions J 1 and J 2 and so can be implemented in the case of a Mooney-Rivlin strain energy function, or other hyperelastic materials. Here, however, for the sake of illustration, we describe some specific implementations for the case of neo-Hookean media, using the expressions for J 1 and J 2 derived in the previous subsection above. 6. A parameter study on a hollow thick-walled quasilinear viscoelastic sphere under pressure Many physical models, including multiple examples of soft composite materials and biological soft tissues, particularly hollow viscera, take the approximate form of hollow thick-walled spheres. Let us first assess the accuracy and convergence of the scheme as we reduce the size of our time step ∆t. We shall work once again in terms of the non-dimensional parameters introduced in (4.2) and (4.3) and, unless otherwise specified, we takeμ ∞ = 0.5,τ = 1,B = 2. Figure 6 depicts the load and unload history as applied to the hollow sphere. After an initial load, then hold (creep) phase, followed by unloading, the sample is subjected to cyclic loading of the |sin 2t| type frequently employed to characterize dynamic viscoelastic behaviour. This is also useful for us in order to test the convergence of the time-step discretization ∆p. Predictions ofα are given in figure 6b according to the choice of the discretization ∆t. In the case when the stress gradient is not too large, even the discretization ∆t ∼ 0.50 appears to be in very good agreement with finer discretization time steps. Indeed, the prediction for ∆t ∼ 0.10 is in agreement with the ∆t ∼ 0.01 predictions almost everywhere, except the scenarios of very high-rate loading (which are not of interest in the quasi-static loading case), for example, here when |sin (2t)| approaches zero.
When quasi-static loading/unloading is considered, we can conclude that discretization with ∆t ∼ 0.10 is sufficient to guarantee the convergence of the numerical stress-strain predictions while keeping computational costs relatively low. Next, we analyse predictions forα when the relaxation time and the thickness of the hollow sphere are varied, according to the three specific load/unload pressure curves in figures 7a and 8a, 8c, and for which the discretization ∆t ∼ 0.1 will be considered.
In figure 7, we present results associated with the case when we subject the hollow viscoelastic thick-walled shell to a loading of the type given in figure 7a. Figure 7b illustrates the subsequent creep and recovery curves for fixedτ = 1 while varyingB,whereasfigure 7c shows predictions for a fixed thicknessB = 2 while varyingτ . In both cases, we takeμ ∞ = 0.5. Note that, for thick shells, curves converge relatively quickly to the limiting curve predicted for the spherical cavity case (B →∞), while figure 7c illustrates that larger relaxation timesτ imply slower relaxation rates and therefore, as expected, it takes more time for the material to fully recover. Finally, figure 7d illustrates the loading and unloading cycleα versus ∆p, indicating the typical hysteresis loops that arise.
In figure 8, we consider two further specific load/unload curves associated with cyclic pressure loading and unloading as depicted in figure 8a,c, where the latter is distinguished from the former by the rest periods between load cycles. In each load/unload figure we illustrate two different load amplitudes up to ∆p = 1( r e dc u r v e s )a n d∆p = 0.5 (black curves). The respective predictions forτ = 1, 2 are shown in figure 8b,d where red and black curves refer to the respective loading curves. In particular, we note that the rest periods allow the material to recover so that subsequentα cycles are very similar to the initial cycle. This is in contrast to the case where there is no rest period; in this case, the maximumα value on each subsequent cycle continues to increase due to viscoelastic creep.

Conclusion
The inflation and deflation of a nonlinear viscoelastic thick-walled spherical shell have been described in the context of incompressible QLV media. This canonical problem has a broad range of applications but is particularly important in the field of biological systems, including hollow viscera. The model described was shown to fit experimental data associated with the volumetric inflation of murine bladders very well for an appropriate strain energy function. Following this, the thin-walled shell limit, associated with a viscoelastic balloon, was considered, building on the significant body of work on the canonical, perfectly elastic balloon inflation problem. We concluded that, in contrast to the elastic balloon scenario, the peak pressure associated with inflation of a neo-Hookean viscoelastic balloon is not independent of the shear modulus of the medium.
In the scenario where pressure is imposed, a new formulation of the problem was required. The governing equation linking the inhomogeneous radial stretch to the imposed pressure difference is a nonlinear Volterra integral equation in the radial coordinate and time. We developed a novel numerical technique to solve this governing equation in the context of a viscoelastic Mooney-Rivlin material and applied it to the case of inflation/deflation of finite-thickness shells, for a range of pressure difference histories (across the shell wall). The great advantage of the proposed quasilinear constitutive model (and numerical solution scheme) is that it is extremely straightforward to implement, especially as it enables the required integration over the spatial variable to be accomplished analytically. Hence, it reduced what could be a time-consuming and complex multiple integration time-stepping exercise to a computationally highly efficient procedure, thus allowing for a comprehensive and wide-ranging parameter study associated with inhomogeneous deformations of nonlinear viscoelastic spherical shells.
Of future interest will be a study of the limit point instability for the viscoelastic balloon problem considered here, which would be an extension of the problem studied by Wineman [42].
Ethics. All experiments were performed using 12-16-week-old C57/BL6 male mice from Charles River (Margate, UK). Mice were acclimatized for 7 days in the laboratory animal husbandry unit under a 12 L:12 D cycle and had free access to water and food. Immediately before the onset of the protocol mice were anaesthetized with isoflurane and humanely sacrificed by cervical dislocation according to UK Home Office legislation regulating Schedule 1 procedures (Scientific Procedure Act 1986). Ethics approval was obtained from the University of Sheffield Ethics Review Panel.
Data accessibility. The experimental data are freely available to download and have the doi:10.17632/ 33n49wtnn2.1.