C Lestringant

B Audoly

A diffuse interface model for the analysis of propagating bulges in cylindrical balloons

With the aim to characterize the formation and propagation of bulges in cylindrical rubber balloons, we carry out an expansion of the non-linear axisymmetric membrane model assuming slow axial variations. We obtain a diffuse interface model similar to that introduced by van der Waals in the context of liquid-vapor phase transitions. This provides a quantitative basis to the well-known analogy between propagating bulges and phase transitions. The diffuse interface model is amenable to numerical as well as analytical solutions, including linear and non-linear bifurcation analyses. Comparisons to the original membrane model reveal that the diffuse interface model captures the bulging phenomenon very accurately, even for well-localized phase boundaries.

Introduction 1.Problem and background

Bulges in cylindrical rubber balloons are a classical example of localization in solid mechanics. When a balloon is inflated, an initial regime of uniform inflation is followed by the formation of a bulge: the bulge appears initially as a long-wavelength buckling mode that localizes rapidly and then grows locally, until it propagates and eventually invades the entire balloon [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]. As in other localization phenomena, the formation of a bulge reflects the non-convexity of the strain energy when restricted to homogeneous deformations [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF]: the onset of bulging occurs quickly after the maximum in the volume-pressure loading curve, and the propagation pressure can be predicted by Maxwell's equal-area rule [START_REF] Chater | On the propagation of bulges and buckles[END_REF].

Several other localization phenomena have been studied in solid mechanics, such as stress-induced phase transformations [START_REF] Ericksen | Equilirbium of bars[END_REF][START_REF] Bhattacharya | Microstructure of martensite[END_REF][START_REF] Fu | Characterization and stability of two-phase piecewise-homogeneous deformations[END_REF], the necking of bars [START_REF] Barenblatt | Neck propagation in polymers[END_REF][START_REF] Hutchinson | Bifurcation analysis of the onset of necking in an elastic/plastic cylinder under uniaxial tension[END_REF], shear bands in plastic materials [START_REF] Rice | The localization of plastic deformation[END_REF], kink bands in compressed fiber composites [START_REF] Wadee | Kink band instabilities in layered structures[END_REF][START_REF] Fu | Continuum-mechanical modelling of kink-band formation in fibre-reinforced composites[END_REF], as well as localized structures in thin elastic shells [START_REF] Power | Localization and propagation of instabilities in long shallow panels under external pressure[END_REF] and tape springs [START_REF] Seffen | Deployment dynamics of tape springs[END_REF][START_REF] Seffen | Folding and deployment of curved tape springs[END_REF].

These localization phenomena have been investigated based on two types of models, as discussed in [START_REF] Triantafyllidis | The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models[END_REF] for example. On the one hand, non-regularized models, also known as sharp interface models, make use of a classical strain energy functional depending solely on the strain: the onset of localization is associated with the loss of ellipticity of the equations of equilibrium at a critical value of the load [START_REF] Knowles | On the failure ot ellipticity and the emergence of discontinuous deformation gradients in plane finite elastostatics[END_REF]. Such models can typically predict the critical load, the formation of different phases and the orientation of the phase boundaries, but cannot predict their subsequent evolution, nor their number or distribution in space; they cannot resolve the displacement inside the localized region either. On the other hand, regularized models, also known as diffuse interface models, make use of a stored elastic energy functional depending on both the strain and the strain gradient: such models remedy the limitations of the non-regularized models, and in particular remain well posed beyond the onset of localization [START_REF] Triantafyllidis | A gradient approach to localization of deformation. I. Hyperelastic materials[END_REF].

Regularized models are often introduced heuristically, but can in some cases be justified mathematically. Such a justification has been done in the case of periodic elastic solids, such as elastic crystals [START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF][START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF], trusses made of elastic bars or beams [START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF], or elastic solids with a periodic micro-structure [START_REF] Triantafyllidis | The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models[END_REF][START_REF] Bacigalupo | Identification of higher-order continua equivalent to a Cauchy elastic composite[END_REF]. In these works on periodic solids, the ratio R/L 1 of microscopic cell size R to the macroscopic dimension L of the structure is used as an expansion parameter, and the homogenized properties of the periodic medium are derived through a systematic expansion in terms of the macroscopic strain and of its successive gradients.

The goal of this paper is to derive a one-dimensional, regularized model applicable to the analysis of axisymmetric bulges in cylindrical rubber balloons. It is part of a general effort to characterize localization phenomena occurring in slender structures, which have been much less studied than in periodic solids. In slender structures, regularized models can be derived by an asymptotic expansion as well, using now the aspect ratio R/L 1 as the expansion parameter where R is the typical transverse dimension of the structure and L its length. This approach has been carried out for the analysis of necking, and diffuse interface models have been derived asymptotically, first for a two-dimensional hyper-elastic strip by Mielke [START_REF] Mielke | Hamiltonian and Lagrangian flows on center manifolds, with application to elliptic variational problems[END_REF] and later for a general prismatic solid in three dimensions by Audoly & Hutchinson [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF]. These authors proposed an expansion method upon which we build ours. This approach has also been applied to soft elastic rods with surface tension which produce localization by the Plateau-Rayleigh instability [START_REF] Xuan | Plateau-Rayleigh instability in solids is a simple phase separation[END_REF].

Our asymptotic expansion starts from the axisymmetric membrane model, which has been used extensively to analyze bulges in cylindrical balloons [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF][START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]. Its outcome is a one-dimensional diffuse interface model, exactly similar to that introduced heuristically by van der Waals [START_REF] Van Der Waals | Thermodynamische Theorie der Kapillarität unter Voraussetzung stetiger Dichteänderung[END_REF] to analyze the liquid-vapor phase transitions at a mesoscopic level. The analogy between bulges in balloons and phase transitions has been known for a long time: Chater & Hutchinson [START_REF] Chater | On the propagation of bulges and buckles[END_REF] have adapted Maxwell's rule for the the coexistence of two phases to derive the pressure at which a bulge can propagate in a balloon, while Müller & Strehlow [START_REF] Müller | Rubber and Rubber Balloons[END_REF] have proposed a pedagogical introduction to the theory of phase transitions based on the mechanics of rubber balloons. Here, we push the analogy further, and show that the diffuse interface model can provide a quantitative description of bulges in balloons, not only accounting for the propagation pressure, but also for the domain boundary between the bulged and unbulged phases, as well as for its formation via a bifurcationborrowing from the theory of phase transitions, we will refer to this boundary as a 'diffuse interface'. The diffuse interface model is classical, tractable, and amenable to analytical bifurcation and post-bifurcation analysis, as we demonstrate. It is also simpler than the axisymmetric membrane model on which it is based.

There is a vast body of work on the bulging of cylindrical balloons, all of which have used the theory of axisymmetric membranes as a starting point. The stability and bifurcations from homogeneous solutions have been analyzed in [START_REF] Corneliussen | Finite Deformation of Elastic Membranes with Application to the Stability of an Inflated and Extended Tube[END_REF][START_REF] Shield | On the stability of finitely deformed elastic membranes, Part II: stability of inflated cylindrical and spherical membranes[END_REF][START_REF] Haughton | Bifurcation of inflated circular cylinders of elastic material under axial loading, part i: membrane theory for thin-walled tubes[END_REF]. Non-linear solutions comprising bulges have been derived in [START_REF] Yin | Non-uniform inflation of a cylindrical elastic membrane and direct determination of the strain energy function[END_REF]. The analysis of stability has been later extended to arbitrary incompressible hyperelastic materials, to various closure conditions at the ends of the tube, as well as to various type of loading controls based on either the internal pressure, the mass or the volume of the enclosed gas [START_REF] Chen | Stability and bifurcation of finite deformations of elastic cylindrical membranes, part I: Stability analysis[END_REF]. In a recent series of four papers, Fu et al. [START_REF] Fu | Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation[END_REF][START_REF] Fu | Stability of localized bulging in inflated membrane tubes under volume control[END_REF][START_REF] Pearce | Characterization and stability of localized bulging/necking in inflated membrane tubes[END_REF][START_REF] Fu | Effects of imperfections on localized bulging in inflated membrane tubes[END_REF] revisit the bifurcation problem, complement it with the weakly non-linear post-bifurcation analysis in the case of an infinite tube, and address imperfection sensitivity. Besides these theoretical studies, there has been a number of experimental and numerical papers on balloons. A compelling agreement between experiments and numerical simulations of the non-linear membrane model has been obtained by Kyriakides & Chang [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF][START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF], who provide detailed experimental and numerical results on the initiation, growth and propagation of bulges, highlighting the analogy with phase transitions. Given that the agreement between experiments and the non-linear membrane theory has already been covered thoroughly in this work, our focus here will be on comparing the diffuse interface model to the non-linear membrane model, using exactly the same material model as in Kyriakides' simulations and experiments.

Outline of the main results

Our work focusses on solutions to the non-linear axisymmetric membrane model that vary slowly in the axial direction, as happens typically at the onset of localization. A systematic expansion of the membrane energy is obtained in terms of the aspect-ratio parameter ε = R/L 1, where R is the initial radius of the balloon and L its initial length. The result reads

E[p, µ] = L 0 G 0 (p, µ(Z)) + 1 2 B 0 (p, µ(Z))µ 2 (Z) dZ (1)
where p denotes the (scaled) internal pressure, a control parameter in the experiments, Z is the axial coordinate and µ(Z) is a strain measure, see figure 1. Specifically, µ is the orthoradial stretch, defined as the ratio µ(Z) = r(Z) R of the current radius r to the initial radius. The potential G 0 appearing in the first term is a non-convex function of the stretch µ, much like in Ericksen's bar [START_REF] Ericksen | Equilirbium of bars[END_REF]: the non-regularized model for the balloon would correspond to the energy functional L 0 G 0 (p, µ(Z)) dZ. Values of p such that several minima of G 0 exist correspond to pressures for which different phases (associated with different values of the stretch µ) are in competition. The second term 1 2 B 0 µ 2 in the integrand is a correction of order ε 2 , that accounts for the energetic cost of inhomogeneity; in the theory of phase transitions, this is the term that would account for surface tension at an interface.

We provide simple and explicit formulas for both the potential G 0 characterizing homogeneous solutions, see §3, equations ( 9), (11) ( 13) and ( 14) in particular, and for the modulus B 0 of the regularizing term, see §4 and equation (20b).

The diffuse interface model is obtained by a systematic, formal expansion. It is asymptotically exact and does not rely on any unjustified kinematic assumptions: equation (1) approximates the energy of the original membrane model with an error of order ε 4 that is negligible compared to the smallest term retained, namely the gradient term of order ε 2 . By contrast, regularized models for slender structures have been proposed in earlier work starting from kinematic hypotheses, which appeared to be incorrect: see the treatment of necking in an elastic cylinder in [START_REF] Coleman | On the Rheology of Cold Drawing[END_REF] as well as the critical discussion in [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF].

Our derivation is based on a finite-strain membrane model. The non-linear features of the elastic constitutive law at finite strain are ultimately reflected in the diffuse interface model through the non-linear potential G 0 (p, µ) and through the dependence of the second gradient coefficient B 0 (p, µ) on the current strain µ. By contrast, an assumption of small strain has been used in previous work [START_REF] Dai | Phase transitions in a slender cylinder composed of an incompressible elastic material. I. Asymptotic model equation[END_REF][START_REF] Dai | An analytical study on the geometrical size effect on phase transitions in a slender compressible hyperelastic cylinder[END_REF] on the justifica-tion of a diffuse interface model to analyze phase transformations in an elastic cylinder: this assumption is questionable since the presence of coexisting phases involves finite variations of strain across the interface.

The outline of the paper is as follows. In §2, we introduce the non-linear membrane model. In §3 we analyze its homogeneous solutions, and derive an expression for the potential G 0 . Section 4 is the core of the paper, and establishes the diffuse interface model (1) by an asymptotic method. Section 5 derives solutions to the diffuse interface model using various methods, and compares them with the predictions of the original membrane model.

Non-linear membrane model

We consider a cylindrical membrane with uniform initial thickness H and radius R. We use the cylindrical coordinates (Z, θ) in reference configuration as Lagrangian variables. When subject to external load, the cylinder deforms into an axisymmetric membrane, see figure 1a. The cylindrical coordinates of a material point in actual configuration are written as (r(Z), θ), corresponding to a position x(Z, θ) = z(Z) e z + r(Z) e r (θ), where (e r (θ), e θ (θ), e z ) is the local cylindrical basis.

In the axisymmetric membrane theory, the deformation gradient is a 3 × 2 matrix which writes

F = µ e θ ⊗ e θ + (R µ e r + λ e z ) ⊗ e z
where we have defined an apparent axial stretch λ and the circumferential stretch µ as

λ(Z) = dz dZ (Z), µ(Z) = r(Z) R . (2)
The Green-Lagrange strain tensor E = 1 2 (F T • F -1) is a 2 × 2 diagonal matrix in the basis (e θ , e z ) tangent to the undeformed mid-surface: it will be represented compactly as a vector, whose entries are the diagonal components E θ and E z of the matrix,

E(λ, µ, µ ) = E θ E z = 1 2 µ 2 -1 λ 2 + (Rµ ) 2 -1 , (3)
where µ = dµ dZ is a stretch gradient, namely the axial gradient of circumferential stretch. A material model is now specified through a strain energy per unit volume w * (E θ , E z ). In previous work on axisymmetric membranes [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF][START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF][START_REF] Fu | Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation[END_REF], the 2d material model proposed by Ogden [START_REF] Ogden | Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubber-like solids[END_REF] for incompressible rubber has been used:

w * (E θ , E z ) = 3 i=1 S i α i αi θ + αi z + 1 θ z αi , with θ = √ 2E θ + 1 (circumfer. stretch), z = √ 2E z + 1 (axial stretch). (4)
We use this model as well for our numerical examples, with the same set of material parameters α i 's and S i 's as used in previous work, namely S 1 = 617 kPa, S 2 = 1.86kPa, S 3 = -9.79kPa, α 1 = 1.3, α 2 = 5.08 and α 3 = -2. All our results can be easily adapted to a different constitutive law. For this constitutive law, the initial shear modulus S ini can be obtained as S ini = in explicit form in §55.1, are natural: the applicable equilibrium condition will emerge automatically from the condition that the energy is stationary. As in the experiments of Kyriakides & Chang [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF], the membrane is subject to an interior pressure p * and to a stretching force F * applied along the axis, see figure 1. The total potential energy reads

E * memb = L 0 w * (E) 2 π R H dZ -π (R µ) 2 (λ dZ) p * -(λ dZ) F *
where 2 π R H dZ is the initial volume element (equal to the current volume element in the particular case of an incompressible material as discussed in this paper), π (R µ) 2 (λ dZ) = π r 2 dz is the current enclosed volume element, and λ dZ = dz is the current axial length element. We introduce a rescaled energy, denoted without an asterisk as

E memb = E * memb (2π R H) Sini : E memb [p, λ, µ] = L 0 w(E(λ, µ, µ )) -p e 2 λ µ 2 -F λ dZ. ( 5
)
The strain energy, the force and pressure have been rescaled as well, as w = w * Sini , F = F * 2πRHSini and p = p * Sini , respectively, and e = R H is an initial aspect ratio. In our numerical examples, we use the same value e = 55 16

as in [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]: even though this balloon is relatively thick prior to deformation, the non-linear membrane model has been checked to match the experimental results accurately in [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]. We also use the same value of the load

F = 1.149 (6)
as in these experiments. The parameter F will never be changed, and we do not keep track of how the various quantities depend on F ; the argument F will systematically be omitted in functions, as we did already in the left hand side of ( 5).

The functions λ(Z) and µ(Z) that make the energy (5) stationary yield the axisymmetric equilibra of the balloon. These solutions are obtained by a numerical method described in section 55.1, and are plotted as the double-struck grey curves in figure 3, where they are used as a reference.

Analysis of homogeneous solutions

Our general goal is to justify the diffuse interface model when λ(Z) and µ(Z) vary slowly as a function of Z. In this section, we start by considering the case where λ and µ do not depend on Z,

dλ dZ = 0, dµ dZ = 0. (7)
This corresponds to homogeneous solutions, i.e. to solutions with uniform inflation. These homogeneous solutions are well known, and are re-derived here for the sake of completeness. A catalog of such homogeneous solutions will be obtained, which plays a key role in the subsequent derivation of the diffuse interface model.

Kinematics of homogeneous solutions

For homogeneous solutions, the gradient term µ in (3) vanishes and the membrane strain reads

E 0 (λ, µ) = E θ 0 E z 0 = 1 2 µ 2 -1 λ 2 -1 . ( 8
)
All the quantities pertaining to homogeneous solutions are denoted using a subscript '0'. In the homogeneous case, the strain energy becomes

w 0 (λ, µ) = 1 S ini 3 i=1 S i α i λ αi + µ αi + 1 λ µ αi . ( 9
)
Of particular importance will be the second Piola-Kirchhoff membrane stress Σ 0 , defined as the gradient of the strain energy with respect to the strain:

Σ 0 (λ, µ) = Σ θ 0 Σ z 0 = ∂w ∂E θ ∂w ∂Ez E=E0 = 1 µ ∂w0 ∂µ (λ, µ) 1 λ ∂w0 ∂λ (λ, µ) . ( 10
)

Equilibrium of homogeneous solutions

In view of ( 5), the total potential energy of a homogeneous solution per unit reference length is

g 0 (p, λ, µ) = w 0 (λ, µ) -p e 2 λ µ 2 -F λ. (11)
Given the load parameters p (and F ) the equilibrium values of λ and µ are found by the condition of equilibrium in the axial and transverse directions,

∂g 0 ∂λ (p, λ, µ) = 0, ( 12a
)
∂g 0 ∂µ (p, λ, µ) = 0. ( 12b
)
We leave the load p left unspecified for the moment, and we view the axial equilibrium (12a) as an implicit equation for λ = λ 0 (p, µ) in terms of p and µ: by definition, λ 0 (p, µ) is the solution to the implicit equation

∂g 0 ∂λ (p, λ 0 (p, µ), µ) = 0. (13)
From now on, we will systematically eliminate λ = λ 0 (p, µ) in favor of the second unknown µ. Starting with the potential g 0 , we define a reduced potential G 0 as

G 0 (p, µ) = g 0 (p, λ 0 (p, µ), µ), (14)
as well as the stress n 0 dual to µ,

n 0 (p, µ) = - ∂G 0 ∂µ (p, µ). (15a)
This n 0 (p, µ) can be interpreted as an imbalance of hoop stress; it vanishes at equilibrium,

n 0 (p, µ) = 0. (15b) Indeed, we have n 0 = -∂G0 ∂µ = -dg0(p,λ0(p,µ),µ) dµ = -∂g0 ∂µ -∂λ0 ∂µ ∂g0
∂λ , where the both terms are zero by the equilibrium conditions (12b-13).

To summarize, we view λ as an internal variable slaved to the 'macroscopic' variable µ (the roles of λ and µ could be exchanged but the other way around would be more complicated as the mapping from λ to µ is not single-valued). A catalog of homogeneous solutions can be obtained by (i ) solving the axial equilibrium (12a) for λ = λ 0 (p, µ), (ii ) defining a reduced potential energy G 0 (p, µ) by ( 14), and (iii ) solving the equilibrium condition n 0 (p, µ) = 0 in the (p, µ) plane.

This program has been carried out and the results are shown in figure 2. The homogeneous stretch λ 0 (p, µ) and the potential G 0 (p, µ) are shown in parts a and b of the figure. In figure 2c, the pressure is plotted in terms of µ and is seen to increase, attain a local maximum p C = 0.1646, decrease, attain a local minimum p C = 0.1002, and finally increase again. The points of extremal pressure are where the onset of localization is expected to occur in a infinite medium (L = ∞) according to Considère's criterion [START_REF] Considère | Mémoire sur l'emploi du fer et de l'acier dans les constructions[END_REF]: we will refer to them as Considère points. For intermediate values of the pressure, p C < p < p C , the potential G 0 (p, µ) plotted in figure 2b has two minima and one maximum as a function of µ. The non-convexity of G 0 makes it possible for the bulged and unbulged domains to coexist, as recalled in the next section; the diffuse interface model derived later in §4 will be able to account for the boundary between these domains.

Maxwell's construction

In a first attempt to address inhomogeneities, we consider solutions made up of two phases, with respective properties (λ a = λ 0 (p, µ a ) , µ a ) and (λ b = λ 0 (p, µ b ) , µ b ). Discontinuities are allowed for the moment, their contribution to the energy being ignored: gradient term µ appearing in the membrane model are simply discarded. Let c denote the fraction of the phase 'a', and (1 -c) the fraction of the phase 'b', as measured after pulling everything back in the reference configuration.

Under these assumptions, the membrane energy (5) takes the form

E 0 (p, c, µ a , µ b ) = L (c G 0 (p, µ a ) + (1 -c) G 0 (p, µ b )).
Optimizing with respect to the µ i 's and to c successively, we find

n 0 (p, µ i ) = 0 (equilibrium in each phase) G 0 (p, µ b ) -G 0 (p, µ a ) = 0
(equilibrium between the two phases)

The second equation is analogous to the equality of chemical potential between two phases in equilibrium in the theory of liquid-vapor phase transformations. These equations can be solved for p and the µ i 's: in particular this selects a value of the pressure p = p M , known as Maxwell's pressure, where the two phases can coexist. The propagation pressure p M is a function of both the applied force F and of the constitutive model for the membrane, but this is implicit in our notation. For F = 1.149 and for the particular values of the constitutive parameters used here, we have obtained the Maxwell load as p M = 0.1087, see the red line joining the points labeled M a and M b in figure 2 Maxwell's equal-area rule for the propagation pressure can be rederived as follows. The quantity G 0 (p, µ b ) -G 0 (p, µ a ) appearing in ( 16) can be written as the integral of dG 0 along the curve corresponding to homogeneous solutions in the (p, µ) plane. Along this curve, ∂G0 ∂µ = -n 0 = 0 by (15a). Therefore, dG 0 = ∂G0 ∂p (p, µ) dp = dg0(p,λ0(p,µ),µ) dp dp = ∂g0 ∂p + ∂g0 ∂λ ∂λ0 ∂µ dp = ∂g0 ∂p dp after using [START_REF] Seffen | Deployment dynamics of tape springs[END_REF]. In view of [START_REF] Fu | Continuum-mechanical modelling of kink-band formation in fibre-reinforced composites[END_REF], this can be written as dG 0 = e 2 v 0 dp, where v 0 (p, µ) = λ 0 (p, µ)µ 2 = π r 2 dz π R 2 dZ denotes the ratio of the deformed to the undeformed volume of homogeneous solutions. Using [START_REF] Knowles | On the failure ot ellipticity and the emergence of discontinuous deformation gradients in plane finite elastostatics[END_REF], the variation of G 0 from one Maxwell point µ a to the other µ b is zero, and so µ b µa v 0 (p, µ) dp = 0. This equality implies the equality of the area of the shaded regions in figure 2d, which uses v 0 as the horizontal axis and p as the vertical axis.

Derivation of the diffuse interface model

We proceed to derive the diffuse interface model from the non-linear membrane theory. This reduction combines an assumption of scale separation, whereby the solution is assumed to vary on a length scale L much larger than the radius R, and the elimination of the unknown λ in favor of µ by means of the relation λ = λ 0 (p, µ).

Principle of the expansion

We assume scale separation and use the convention that the radius R is fixed and finite while L = R/ε goes to infinity: the solution is sought in terms of a scaled variable Z = ε Z through scaled functions λ and μ, where ε 1 is our expansion parameter,

λ ε (Z) = λ(ε Z), µ ε (Z) = μ(ε Z).
As a consequence of this scaling assumption, any derivative with respect to the slow axial variable Z entails a multiplication by the small parameter ε:

λ ε (Z) = λ(ε Z) = O(1), µ ε (Z) = μ(ε Z) = O(1), dλ ε dZ = d λ d Z (ε Z) = O(ε), dµ ε dZ = ε dμ d Z (ε Z) = O(ε).
For the sake of legibility, we drop the subscripts ε and remove any reference to the scaled functions λ and μ in the following: it will be sufficient for us to use the above order of magnitude estimates.

Derivation of the gradient effect by a formal expansion

The general expression of the membrane strain (3) can be split in two terms

E(λ, µ, µ ) = E 0 (λ, µ) + E 1 (µ ), (17a)
where the first one depends on the stretch and the second one on the stretch gradient,

E 0 (λ, µ) = 1 2 µ 2 -1 λ 2 -1 , E 1 (µ ) = 1 2 0 (Rµ ) 2 . ( 17b
)
In view of the results from the previous section, their orders of magnitude are

E 0 (λ, µ) = O(1) E 1 (µ ) = O(ε 2 ). ( 17c
)
In line with the fact that we use the finite elasticity theory, the strain E 0 is of order 1. E 1 being a small correction to E 0 , the strain energy density can be expanded as

w(E) = w(E 0 (λ, µ) + E 1 (µ )) = w(E 0 (λ, µ)) + ∂w ∂E (E 0 (λ, µ)) • E 1 (µ ) + O(|E 1 | 2 ) = w 0 (λ, µ) + Σ 0 (λ, µ) • E 1 (µ ) + O(ε 4 ) (18)
where we have used the definition of the membrane stress Σ 0 in [START_REF] Wadee | Kink band instabilities in layered structures[END_REF]. Inserting this into (5) yields the following approximation of the energy

E memb [p, µ] = L 0 g 0 (p, λ(Z), µ(Z)) + Σ 0 (λ(Z), µ(Z)) • E 1 (µ (Z)) dZ + O(L ε 4 ). ( 19
)
Note that the gradient of axial stretch λ does not appear in this expression.

Energy of the diffuse interface model

An important result, proved in appendix A, is that it is consistent, at this order of approximation, to replace the unknown λ(Z) with the axial stretch λ 0 (p, µ(Z)) of the homogeneous solution having the local value of µ(Z) as its circumferential stretch. This eliminates λ(Z) from the equations, and we obtain the diffuse interface model

E[p, µ] = L 0 G 0 (p, µ(Z)) dZ + 1 2 L 0 B 0 (p, µ(Z)) µ 2 (Z) dZ. ( 20a
)
where we have omitted terms of order L ε 4 and higher. The coefficient B 0 of the regularizing term has a simple expression which is found by identifying with [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF] and [START_REF] Wadee | Kink band instabilities in layered structures[END_REF] as

B 0 (p, µ(Z)) = R 2 1 λ ∂w 0 ∂λ (λ, µ(Z)) λ=λ0(p,µ(Z)) . ( 20b
)
This defines the regularizing term in terms of the energy w 0 (λ, µ) of homogeneous solutions, see [START_REF] Rice | The localization of plastic deformation[END_REF]. Even though this is implicit in our notation, both G 0 and B 0 depend on the force F . Equations (20a-20b) are our main result, and can be restated as follows. The energy E memb of the full non-linear membrane model can be approximated as the sum of (i ) the non-regularized energy G 0 dZ which depends on the stretch µ but not on its gradient, and is of order L, and (ii ) a much smaller correction

1 2
B 0 µ 2 dZ, of order ε 2 , that depends on the strain µ and as well as on its gradient µ = dµ dZ . These two terms provide an approximation of the full energy E of the non-linear membrane model which is accurate up to order L ε 4 .

Non-linear equilibrium of the diffuse interface model

The equilibrium equations are obtained from (20a) by the Euler-Lagrange method as

n 0 (p, µ(Z)) - 1 2 ∂B 0 ∂µ (p, µ(Z)) µ 2 (Z) + d dZ (B 0 (p, µ(Z)) µ (Z)) = 0. (21a)
In the absence of kinematic constraints, the variational method yields the natural conditions at the endpoints as well,

µ (0) = µ (L) = 0. ( 21b
)
Here, µ (L) = 0 is consistent with the symmetry condition at the center Z = 0 of the bulge. The equilibrium condition (21a) reduces to the condition (15b) applicable to homogeneous solutions, namely n 0 (p, µ) = 0, when the gradient effect is removed, by setting B 0 = 0.

Solution for a domain boundary in an infinite balloon

The existence of a first integral associated with the equilibrium (21a) has been noted by a number of authors such as Coleman & Newman [START_REF] Coleman | On the Rheology of Cold Drawing[END_REF]. It can be obtained by expanding the derivative in the last term in the right-hand side, and by multiplying the entire side by µ (Z); the result is d(-G0+B0 µ 2 ) dZ = 0. This shows that the following quantity is conserved:

-G 0 (p, µ(Z)) + B 0 (p, µ(Z)) µ 2 (Z) = C. ( 22
)
This equation can be used to solve for µ(Z) by quadrature. However, this method is impractical for numerical calculations as it involves evaluating integrals that are close to singular, even when the singular parts are taken care of analytically [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF]. This is why our numerical simulations in §55.1 use a direct integration method of the equilibrium (21a) rather than the quadrature method.

In the case of the boundary separating two domains in an infinite medium, however, the quadrature method is tractable. Then, the pressure matches Maxwell's pressure, p = p M , and µ(Z) tends to µ a and µ b for Z → ±∞, respectively. The value of C consistent with these asymptotic behaviors is the common value C = G 0 (p M , µ a ) = G 0 (p M , µ b ) of the potential, see [START_REF] Knowles | On the failure ot ellipticity and the emergence of discontinuous deformation gradients in plane finite elastostatics[END_REF]. The implicit equation [START_REF] Mielke | Hamiltonian and Lagrangian flows on center manifolds, with application to elliptic variational problems[END_REF] can then be plotted in the phase space (µ(Z), µ (Z)) using a contour plot method. We have checked that the resulting curve (not shown) falls on top of the dotted green curve labeled A in figure 3d, obtained by numerical integration of the equilibrium with a large but finite aspect ratio, L/R = 30: the analytical solution [START_REF] Mielke | Hamiltonian and Lagrangian flows on center manifolds, with application to elliptic variational problems[END_REF] in an infinite balloon provides an excellent approximation to a propagating interface in a finite balloon, as long as it is sufficiently remote from the endpoints. In the bifurcation diagram, the numerical solutions appears as a point A lying almost exactly on Maxwell's plateau, see figure 3a.

This analytical solution is also an excellent approximation to the domain boundary predicted by the original membrane model, as discussed below in §55.1.

Domain of validity of the diffuse interface model

The diffuse interface model has just been obtained from an expansion with respect to the small aspect ratio parameter ε = R/L. It is therefore valid for long balloons. Its domain of validity is in fact broader. Indeed, the only assumption which we have really used in the expansion is |R µ | 1, see (17b). In terms of the amplitude of variation of µ across the length, ∆µ ∼ L |µ |, the condition of validity of the expansion writes ε ∆µ 1. This shows that the diffuse interface model is applicable in two asymptotic limits: (i ) long balloons (small ε, finite ∆µ) and (ii ) short balloons close to bifurcation (finite ε, small ∆µ).

Comparison of the diffuse interface and membrane models

Using a formal expansion method, we have shown that the 2d non-linear axisymmetric membrane model ( §2) is asymptotically equivalent to the 1d diffuse interface model in [START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF]. This equivalence holds for 'slowly' varying solutions, i.e. when the axial gradients involve a length scale much larger than the tube radius, |dµ/dZ| 1/R. Here, we compare the predictions of the approximate diffuse interface model to those of the original membrane model. The goal is twofold. First, we verify our asymptotic expansion by checking consistency for slowly-varying solutions. Second, we push the diffuse interface model outside its domain of strict mathematical validity, by applying it to problems involving sharp boundaries and comparing to the predictions of the original membrane model.

Comparison of the full bifurcation diagrams

We start by comparing the bifurcation diagrams obtained with each one of the models for balloons of finite length L, see figure 3a. In this numerical example, both the membrane and the diffuse interval models use the constitutive law in (4), the standard set of material parameters listed below this equation, and the value of the pulling force F in (6). We limit our attention to solutions that are either homogeneous or comprise a single bulge centered at Z = 0: recall that the simulation domain (0, L) represents one half of a real balloon.

The equilibrium equations of the membrane model are obtained from the energy ( 5) by an Euler-Lagrange method, and are solved numerically. These equations of equilibrium and their numerical solution have already been documented in [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF], and we refer to this work for details; in our work, the solution branches were calculated using the path-following method from the AUTO-07p library [START_REF] Doedel | AUTO-07p: continuation and bifurcation software for ordinary differential equations[END_REF]. While [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF] used boundary conditions representing rigid plugs, we use instead the natural boundary conditions, namely the axial and radial equilibria

F = ∂w ∂ z z z -p e 2 θ
2 and µ = 0. These boundary conditions are relevant to the soft boundary device sketched in figure 1b, and are enforced at Z = L. At the center of bulge Z = 0, we impose the symmetry conditions µ = 0 and z = 0.

To solve the diffuse interface model ( 21) numerically, we first sample the functions n 0 (p, µ) and B 0 (p, µ) numerically. This tabulation is available as by-product of the analysis of homogeneous solutions from Section 3. Next, the solution branches are generated by solving the boundary value problem ( 21) using the pathfollowing library AUTO-07p. Alternatively, we tried solving this boundary value problem by the quadrature method described earlier, but it did not work well for the reason already explained.

The bifurcation diagrams are shown in figure 3a. The homogeneous solutions are plotted using the thick, dark blue curve: they are identical for both models, and are also identical to those derived earlier in figure 2d. Bifurcated solutions are shown as black double-struck curves (membrane model) and green dots (diffuse interface model) for different value of L = L/R. The bifurcation diagram uses the natural logarithm of the scaled volume v on the horizontal axis,

v = 1 L/R L 0 µ 2 λ 0 (p, µ) dZ R . ( 23
)
This is consistent with the definition of the scaled volume v 0 used in the analysis of homogeneous solutions.

For large values of the aspect-ratio L/R, the bifurcated branches display a plateau corresponding to Maxwell's pressure p M . The diffuse interface model appears to be highly accurate, as its bifurcation diagram is almost identical to that of the membrane model: in the figure, the green dots fall exactly onto the double-struck curves. Given that the diffuse interface model has been derived under an assumption of 'slow' axial variations, it could be expected that the models would agree near the bifurcation points (in the neighborhood of the dark blue curve) where the localization is mild, see §44.6. We did not anticipate the good agreement far from the bifurcation point, for configurations featuring relatively sharp interfaces such as that labeled A in the figure: for this solution, the largest value of the stretch gradient is 1.2, see figure 3d-even though this is not a small number, the diffuse interface model remains remarkably accurate.

Selected deformed configurations are plotted in figure 3c in real space: the predictions of both models are still indistinguishable, even inside the domain boundary. The predictions of the two models are not exactly identical, however: a small difference is visible when these solutions are represented in phase space, see figure 3d; in phase space, the subtle features of the interface are highlighted, while the uniform domains shrink to the points labeled 'a' and 'b' in the figure .  We do not provide a detailed analysis of the stability of the solutions of the diffuse interface model here, as it is already available from the literature (and it can be guessed from the bifurcation diagrams based on standard stability-exchange arguments). For a volume-controlled experiment for instance, the analysis of [START_REF] Carr | Structured phase transitions on a finite interval[END_REF] yields the stability diagram in figure 3b, which includes hysteresis loops. In particular, the solutions labelled A and A in figure 3a are unstable for this type of loading. Overall, the stability results based on the diffuse interface model appear to be exactly similar to those obtained from the membrane model [START_REF] Fu | Stability of localized bulging in inflated membrane tubes under volume control[END_REF].

To sum up, the diffuse interface model reproduces the entire bifurcation diagram of the original membrane model with good accuracy, even for well localized domain boundaries. In the following sections, we show that it is also well suited to linear and non-linear buckling analyses.

Onset of bulging: linear bifurcation analysis, finite length

We now compare the bifurcation load at the onset of bulging, as predicted by the diffuse interface model on the one hand, and by the membrane model on the other hand. The diffuse interface model yields a simple analytical prediction, that matches that of the membrane model exactly. The critical load of the diffuse interface model is derived by a classical linear bifurcation analysis as follows. Consider a perturbation to a homogeneous solution µ 0 , in the form µ(Z) = µ 0 + µ 1 (Z). Linearizing the equilibrium equation of the diffuse interface model in (21a) with respect to µ 1 , we obtain

∂n 0 ∂µ (p, µ 0 ) µ 1 (Z) + B 0 (p, µ 0 ) µ 1 (Z) = 0. ( 24
)
The boundary conditions are µ 1 (0) = 0 and µ 1 (L) = 0. The first critical mode µ 1 (Z) = cos π Z L corresponds to half a bulge in the simulation domain (0, L). When inserted into the above expression, this yields

∂n 0 ∂µ (p, µ 0 ) = B 0 (p, µ 0 ) R 2 π 2 L 2 (diffuse interface model). (25)
This equation must be solved together with the axial equilibrium condition for the unperturbed solution (15b), n 0 (p, µ 0 ) = 0. For any given value of the aspect-ratio L, the roots (p (L), µ (L) = µ 0 (p (L))) of these two equations define the critical parameters where the bifurcation occurs. The corresponding scaled volume can then be reconstructed as v (L) = v 0 (p (L), µ (L)). The dependence of the critical volume v on the aspect-ratio is shown by the green dots in figure 4a.

For comparison, we derive the bifurcation load predicted by the membrane model. The critical load for hard plugs has been obtained by Kyriakides & Chang [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]. Adapting their bifurcation analysis to soft plugs, we obtain the bifurcation condition as

λ 0 (p, µ 0 ) (p e λ 0 (p, µ 0 ) -w 0,µ 2 ) w 0,λ 2 + (w 0,λµ -p e µ 0 ) 2 w 0,λ w 0,λ 2 = π 2 L 2 (membrane model) (26)
where commas in subscripts denote partial derivatives of the homogeneous strain energy w 0 defined in [START_REF] Rice | The localization of plastic deformation[END_REF], all of which are evaluated at (λ, µ) = (λ 0 (p, µ 0 ), µ 0 ). Solving this equation together with the axial equilibrium yields p and µ = µ 0 (p ) as a function of L, as earlier.

The bifurcation loads of the diffuse interface model (green dots) and of the membrane model (solid dark gray curve) are compared in figure 4a. They agree exactly. This is not surprising as, close to bifurcation, the solutions of the membrane model depart from a uniform solution by an arbitrarily small perturbation, implying that the axial gradients are arbitrarily small: the assumptions underlying the diffuse interface model are satisfied close to the bifurcation point. The diffuse interface model captures exactly the retardation of the onset of buckling in balloons of finite length. Similarly, the critical load predicted by the one-dimensional diffuse interface model for the analysis of necking in solid cylinders has been found in [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF] to agree exactly with that based on the three-dimensional analysis of [START_REF] Hutchinson | Bifurcation analysis of the onset of necking in an elastic/plastic cylinder under uniaxial tension[END_REF].

For large values of the aspect-ratios L, the bifurcation equation ( 25) can be simplified by noticing that the left-hand side goes to zero, as the bifurcations takes place closer and closer to the Considère point (p C , µ C ) where the load is maximum, ∂n0 ∂µ (p C , µ C ) = 0. Expanding the left-hand side accordingly, one obtains from ( 25)

µ ≈ µ C + B 0 (p C , µ C ) R 2 ∂ 2 n0 ∂µ 2 (p C , µ C ) π 2 L 2 (diffuse interface model, large L). (27)
This yields the dotted curve in figure 4a, which is indeed consistent with the two other curves in the limit L → ∞.

Onset of bulging: weakly non-linear bifurcation analysis, finite length

Following the method of Lyapunov and Koiter [START_REF] Koiter | On the stability of elastic equilibrium[END_REF], an expansion of the bifurcated branch can be found by introducing a small arc-length parameter η and expanding p and µ as

µ = µ 0 (p) + η µ 1 (Z) + η 2 µ 2 (Z) + O(η 3 ) (28a) p = p + η 2 p 2 + O(η 3 ), (28b)
where µ 0 (p) denotes the branch of homogeneous solutions satisfying the equilibrium condition n 0 (p, µ 0 (p)) = 0, see (15b), p = p (L) is the critical pressure found by the linear stability analysis, µ 1 (Z) = cos πZ L is the linear bifurcation mode, and µ 2 (Z) and p 2 are higher-order corrections. The latter are now determined by inserting this expansion into the non-linear equilibrium [START_REF] Bacigalupo | Identification of higher-order continua equivalent to a Cauchy elastic composite[END_REF], and by solving it order by order in η.

It is actually preferable to work with the weak form of the equilibrium (principle of virtual work), which formally writes E ,µ (p, µ) • [μ] = 0, for any kinematically admissible virtual stretch μ. Here, E ,µ (p, µ) denotes the first variation of the total potential energy, defined as

E ,µ (p, µ) • [μ] = lim t→0 E(p, µ + t µ) -E(p, µ) t .
Higher-order variations of the energy are defined similarly. When the expansion ( 28) is inserted into the principle of virtual work, one obtains, at order η, the condition ∀μ, E ,µ 2 (p , µ 0 (p ))

• [µ 1 , μ] = 0. ( 29
)
We have recovered the bifurcation condition [START_REF] Xuan | Plateau-Rayleigh instability in solids is a simple phase separation[END_REF], which is automatically satisfied by the linear mode µ 1 (Z). At order η 2 the expansion yields

∀μ, E ,µ 2 (p , µ 0 (p )) • [µ 2 , μ] + 1 2 E ,µ 3 (p , µ 0 (p )) • [µ 1 , µ 1 , μ] = 0. ( 30
)
The first term in the left-hand side involves the tangent stiffness operator E ,µ 2 (p , µ 0 (p )) which is known to be singular by the bifurcation condition [START_REF] Shield | On the stability of finitely deformed elastic membranes, Part II: stability of inflated cylindrical and spherical membranes[END_REF]. Therefore a solvability condition must be verified before attempting to solve equation ( 30) for µ 2 (Z); it is derived by replacing μ with µ 1 and reads E ,µ 3 (p , µ 0 (p )) • [µ 1 , µ 1 , µ 1 ] = 0. In the left-hand side, the interactions of the three modes µ 1 produces harmonic waves with is similar to that derived Fu et al. [START_REF] Fu | Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation[END_REF] based on the full membrane model but its derivation is somewhat simpler.

In the limit L/R → ∞, the bifurcation takes place at the Considère point (µ C , p C ), where the pressure attains its maximum (the other bifurcation taking place at the minimum pressure p C , which can be treated similarly). Accordingly, the weakly bulged solution satisfies p ≈ p C and µ(Z) ≈ µ C , and the potential G 0 can be expanded as

G 0 (p, µ) = G C 0 + G C 0,p (p -p C ) + G C 0,µ (µ -µ C ) • • • + G C 0,p 2 2 (p -p C ) 2 + G C 0,pµ (p -p C ) (µ -µ C ) + G C 0,µ 2 2 (µ -µ C ) 2 + G C 0,µ 3 6 (µ -µ C ) 3 + • • •
The arguments appearing in subscript after a comma denote a partial derivative, while a superscript 'C' means that the function is evaluated at Considère's point (p C , µ C ). The values of all the coefficients G C 0 , G C 0,p , etc. are available from the analysis of homogeneous solutions ( §3).

In the right-hand side above, we can discard the terms that do not depend on µ, as well as the terms containing G C 0,µ which cancels by equation (15b), and that containing G C 0,µ 2 which vanishes at the maximum pressure p C . Accordingly, the energy (1) can be approximated as

E ≈ +∞ -∞ Cte(p) + G C 0,pµ (p -p C ) (µ -µ C ) + G C 0,µ 3 6 (µ -µ C ) 3 + 1 2 B C 0 dµ dZ 2 dZ. ( 32
)
The signs of the coefficients appearing in the integrand are important: for our particular constitutive law and using the results of Section 3, their numerical value is

G C 0,pµ = ∂ 2 G 0 ∂p ∂µ (p C , µ C ) = -9.366, G C 0,µ 3 = ∂ 3 G 0 ∂µ 3 (p C , µ C ) = -3.413, B C 0 = B 0 (p C , µ C ) = 0.8956.
A balance argument on the three last term in the integrand above suggests the change of variable

µ(Z) = µ C + (p C -p) 1/2 µ † µ(Z) (33a)
where

µ † = 2 |G C 0,pµ | |G C 0,µ 3 | 1/2 , Z = 2 |G C 0,pµ | |G C 0,µ 3 | (B C 0 ) 2 (p C -p) 1/4 Z. (33b)
In terms of the rescaled variables, the energy expansion [START_REF] Chen | Stability and bifurcation of finite deformations of elastic cylindrical membranes, part I: Stability analysis[END_REF] writes

E = 1 2 +∞ -∞ µ(Z) - µ 3 (Z) 3 + dµ dZ 2 dZ,
after dropping the term G 0 (p, µ C ) that is independent of µ, and rescaling the energy using a numerical constant. The weakly non-linear solutions are the stationary points µ(Z) of this energy functional. They can be analyzed based on the analogy with a mass moving in a potential U (µ) = -µ 2 + µ 3 6 , when Z is viewed as a time variable. A first type of solutions are those corresponding to the equilibria in the effective potential U (µ), namely µ(Z) = ±1: they yield an expansion of the branch of homogeneous solutions near the point of maximum pressure, as can be checked. A second type of solution corresponds to a soliton, i.e. to a non-constant but bounded solution; it can be derived by a quadrature method, using the conservation of the total mechanical energy of the mass in the effective potential. The result is

µ(Z) = -1 + 3 cosh 2 Z-Z0 2 . ( 34
)
This solution represents a weakly localized bulge centered about Z 0 . It is identical to that derived by Fu et al. [START_REF] Fu | Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation[END_REF] based on the full membrane model. The soliton solution [START_REF] Fu | Stability of localized bulging in inflated membrane tubes under volume control[END_REF] is plotted in figure 5b-c, and compared to the non-linear solution of the diffuse interface model: the weakly non-linear solution and the full non-linear solution agree asymptotically close to the bifurcation point, as expected. Unlike the weakly non-linear solution from the previous section, this solution captures the rapid localization of the bulge past the bifurcation point: from (33b), the width of the interface is Z ∼ (p -p C ) -1/4 near bifurcation.

To derive the weakly non-linear solution, we have retained some terms in the expansion of the energy (32), and omitted others such as

1 2 G C 0,µp 2 (µ-µ C ) (p-p C ) 2 or 1 2 B C 0,p (p-p C ) µ 2 .
This can be justified a posteriori, based on the scaling laws of the weakly non-linear solution: the scaling assumptions (µ -µ C ) ∼ (p C -p) 1/2 and Z ∼ (p C -p) -1/4 not only make the last three terms in ( 32) balanced (which is by design), they also make the other terms negligible, as can be checked. The expansion ( 32) is therefore consistent. By the argument given at the end of the previous section, the weakly non-linear solution just derived is unstable for long balloons (see also figure 3b) but stable for short balloons (see the bifurcated curve for L/R = 5.8 in figure 3b).

Final remarks and conclusion

We have proposed a diffuse interface model for the analysis of the formation, localization and propagation of bulges in cylindrical rubber balloons. The model has been derived from the non-linear membrane theory, and is asymptotically exact in the limit where the strain gradient dµ dZ is small compared to 1/R. Analytical and numerical solutions to the diffuse interface model have been obtained, showing good agreement with the predictions of the original non-linear membrane model, both at the onset of localization and for welllocalized solutions: in practice, the diffuse interface model remains accurate well beyond its domain of strict mathematical validity, i.e. even for relative large gradients such as those found at the boundary between a bulged and a non-bulged domain. Hopefully, our work will shed new light on the classical problem of bulging in elastic balloons, and will help highlighting its tight connection with the theory of phase transitions.

The model handles finite strain: only the strain gradients are assumed to be small. The elastic response of the material under finite strain gets reflected into the diffuse interface model through the non-linear potential G 0 (p, µ) and the non-linear strain gradient modulus B 0 (p, µ). Thanks to this feature, the predictions of the model are exact as far as linear bifurcation analyses are concerned, and remain accurate even well into the post-bifurcation regime. By contrast, expansion methods underlying bifurcation analyses typically assume that the solution is close to a homogeneous solution, and have a much narrower domain of applicability.

A simple expression for the coefficient B 0 (p, µ) of the strain gradient term has been established, see (20b). Remarkably, this coefficient is directly proportional to the pre-stress Σ 0 in the homogeneous solution, see [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF]. The same observation has been made concerning the strain-gradient model applicable to necking in a hyperelastic cylinder [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF]. In both cases, the contribution of the strain gradient to the elastic energy occurs through a 'geometric rigidity effect' -the vibration of a string under tension is another example of this geometric rigidity effect, whereby the pre-stress brings in an effective elastic stiffness.

The consistency of our results with those of Audoly & Hutchinson [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF] for a solid cylinder can be checked as follows. In (20a), we have derived the contribution of the strain gradient to the energy of a balloon using dimensionless quantities. In terms of the original (non-scaled) quantities, it can be rewritten as ∂λ , after restoring the initial area that had been scaled out. Identifying I = (2π R H) R 2 as the geometric moment of inertia of the annular cross-section in reference configuration, the contribution of the strain gradient to the energy can therefore be written as 1 2 L 0 I 1 λ ∂w * 0 ∂λ µ 2 (Z) dZ in an axisymmetric membrane. This is identical to the formula applicable to a cylinder as well, as can be checked by combining equations 2.28, 2.13 and 2.8 from [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF].

In the future, the following extensions of the model could be considered. While we have only been concerned with bifurcations in this paper, one could analyze the stability of the solutions based on the diffuse interface model as well: presumably, this would confirm the stability results obtained previously from the non-linear membrane model. Different boundary conditions than the natural boundary conditions could be used at the ends of the balloons: a hard plug can be enforced by the boundary condition µ = 1, while a rounded elastic cap can be prescribed through a non-cylindrical reference configuration R(Z); if, however, the prescribed initial profile R varies quickly with Z, there may be a conflict with our asymptotic procedure that assumes slow variations, and the diffuse interface model may not be able to account accurately for such boundary conditions. A coupling with an electrical field could also be introduced, as is relevant to the actuation of balloons made of dielectric elastomers [START_REF] Lu | Large conversion of energy in dielectric elastomers by electromechanical phase transition[END_REF][START_REF] Lu | Electro-mechanical coupling bifurcation and bulging propagation in a cylindrical dielectric elastomer tube[END_REF][START_REF] An | Experimental investigation of the electromechanical phase transition in a dielectric elastomer tube[END_REF]; we believe that a natural extension of our asymptotic expansion can be derived for dielectric balloons. In future work, we also hope to generalize our asymptotic reduction method and apply it to other localization phenomena occurring in slender structures, such as the formation of localized kinks in tape springs [START_REF] Seffen | Deployment dynamics of tape springs[END_REF][START_REF] Seffen | Folding and deployment of curved tape springs[END_REF].

Ethics statement. This work did not involve any ethical issue. Data accessibility statement. This work does not have any experimental data. Competing interests statement. We have no competing interests. Authors' contributions. Both authors have equally contributed to all aspects of the work. Funding. There has been no dedicated funding for this work. Acknowledgments. We would like to thank John Hutchinson for drawing our attention to the absence of a strain gradient model for the analysis of bulging in balloons, and for providing useful feedback on the manuscript.

A Elimination of λ(Z) from the diffuse interface model

This appendix provides details on the elimination of the unknown λ(Z) from the intermediate form [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF] of the energy, leading to the final form [START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF].

As the approximation of the energy [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF] contains no derivative of the axial stretch λ, optimizing with respect with the λ variable yields an algebraic (i.e. non-differential) problem for λ(Z) ∂g 0 ∂λ (p, λ(Z), µ(Z)) + ∂Σ 0 ∂λ (λ(Z), µ(Z))

• E 1 (µ (Z)) + O(ε 4 ) = 0. ( 35
)
The first term is of order ε 0 = 1, the following term is of order ε 2 by (17c), and we have omitted terms of order ε 4 and higher.

Let us solve [START_REF] Pearce | Characterization and stability of localized bulging/necking in inflated membrane tubes[END_REF] order by order for λ. At dominant order ε 0 = 1, λ(Z) is a solution of ∂g0 ∂λ (p, λ, µ(Z)) = 0, which we identify as the equilibrium condition (12a) applicable to homogeneous solutions. In [START_REF] Seffen | Deployment dynamics of tape springs[END_REF], the solution to this equation has be=en introduced as λ(Z) = λ 0 (p, µ(Z)), i.e. λ(Z) = λ 0 (p, µ(Z)) + O(ε 2 ).

To order ε 2 , the solution of equation ( 35) is found as λ(Z) = λ 0 (p, µ(Z)) + λ [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF] (Z) + O(ε 4 ), [START_REF] Fu | Effects of imperfections on localized bulging in inflated membrane tubes[END_REF] where λ [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF] = -∂Σ0 ∂λ • E 1 (µ ) / ∂ 2 g0 ∂λ 2 is a correction of order ε 2 , arising from the strain correction E 1 (µ ) = O(ε 2 ) (the expression of λ [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF] is given for the sake of completeness but is not needed in the following).

Inserting the expansion [START_REF] Fu | Effects of imperfections on localized bulging in inflated membrane tubes[END_REF] for λ(Z) into [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-d non-linear elasticity from periodic lattice models[END_REF] and expanding in series, we find E = G 0 (p, µ) + ∂g 0 ∂λ (p, λ 0 , µ) λ [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF] dZ + Σ 0 (λ 0 , µ) • E 1 (µ ) dZ + O(L ε 4

where we use λ 0 as a shorthand for λ 0 (p, µ(Z)). As the second term in the bracket vanishes by [START_REF] Seffen | Deployment dynamics of tape springs[END_REF], we are left with E = G 0 (p, µ) dZ + Σ 0 (λ 0 , µ) • E 1 (µ ) dZ + O(L ε 4 ).

Inserting the explicit expressions of Σ 0 and E 1 from ( 10) and (17b), one obtains the expressions announced in [START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF].

Figure 1 :

 1 Figure 1: Inflation of a cylindrical membrane: (a) reference configuration, (b) sketch of an equilibrium configuration when the membrane is subject to an axial force F and an internal pressure p. (c) Equivalent diffuse interface model derived in this paper: a bar having an order parameter µ(Z) undergoes a phase transformation.

Figure 2 :

 2 Figure 2: Analysis of the homogeneous solutions, with F = 1.149. (a) For any value of the pressure p, the axial equilibrium (13) yields an implicit curve in the (µ, λ) plane that defines the stretch λ 0 (p, µ) in terms of µ; the different curves correspond to p = 0.0885, 0.1002, 0.1087, 0.1187, 0.1285, 0.140, 0.1646, 0.200. (b) Reduced potential G 0 as a function of µ for the same set of values of p. Critical points are shown in red: Considère points C and C where the pressure is extremal (disks), and Maxwell construction (line and diamonds). (c) Homogeneous solutions in the (µ, p)-plane, as determined by solving the transverse equilibrium n 0 (p, µ) = 0. The dashed part of the curve between the Considère points is unstable, as ∂ 2 G 0 /∂µ 2 < 0. (d) Same set of homogeneous solutions, now represented in the (v 0 , p) plane where v 0 = µ 2 λ 0 (p, µ) is the ratio of the final to the initial volume. With this set of conjugate variables, Maxwell's rule applies and the two shaded regions have the same area.

Figure 3 :

 3 Figure 3: Comparison of the predictions of the non-linear membrane model and of the diffuse interface model, for F = 1.149. (a) Bifurcation diagrams for different values of the aspect-ratio L/R: homogeneous solutions (dark blue curve), bifurcated branches of the membrane model (double-struck black curves) and of the diffuse interface model (green dots); Considère's points (red dots) and Maxwell's plateau (dotted red line). Note that we use the logarithm of the scaled volume v on the horizontal axis, so that Maxwell's equal-area rule does not apply directly. (b) Sketch of the stability diagram for a volume-controlled experiment (not analyzed in this paper) and a particular value of L/R: stable (solid curves) vs. unstable (dotted curves) solutions; the hysteresis loop obtained for a quasi-static loading is shown by the arrows. (c) Comparison of the deformed configurations in physical space (z/R, µ = r/R) for L/R = 30, for the three configurations labeled A, A and A in (a), corresponding to (v, p) = (2.39, 0.146), (45., 0.109) and (77.43, 0.106) respectively. (d) Same solutions visualized in the phase space: a small discrepancy is visible in the center of the sharp interface (arrow). All plots except (b) are from numerical simulations.

Figure 4 :

 4 Figure 4: Linear and non-linear bifurcation analyses of the bulging instability based on the diffuse interface model, and comparison with the membrane model. (a) Critical volume v at the onset of bulging as a function of the aspect-ratio L = L/R. The diffuse interface model predicts the bifurcation load exactly (the green dots fall exactly onto the double-struck curve). Inset: bifurcation point shown in the bifurcation diagram as in figure 3a for a particular value L = L/R = 5.8. (b) Initial tangents to the bifurcated branches for different values of the aspect-ratio L = L/R, as predicted by a weakly non-linear expansion of the membrane model (green dotted lines); comparison with the non-linear branches predicted by the membrane model (double-struck curves).

Figure 5 :

 5 Figure 5: Comparison of the analytical weakly non-linear solutions (34) and the numerical non-linear solutions of the diffuse interface model, for an aspect-ratio L/R = 30. The simulation domain represents one half of a balloon having a single bulge in the center. (a) Position of a few selected configurations in the phase diagram; point A is the same configuration as in figure 3. (b) Comparison close to the bifurcation point. (c) Comparison for different bulging amplitudes in phase space: the weakly non-linear solution is asymptotically exact for small bulging amplitudes.

0 µ 2 (

 2 Z) dZ, where the non-scaled modulus is found from (20b) as B * 0 = (2π R H) R

i=1 α i S i . The domain 0 ≤ Z ≤ L represents one half of a balloon comprising a single bulge center at Z = 0, with symmetry conditions µ = λ = 0 enforced at Z = 0. At the other endpoint Z = L, we consider the ideal boundary conditions sketched in figure1b, whereby the terminal section of the balloon is resting and freely sliding on a planar 'plug'. These conditions would be difficult to achieve in experiments but they offer the advantage of being compatible with a uniform expansion of the membrane, which simplifies the analysis. By contrast, actual cylindrical balloons are typically closed up on their ends and cannot be inflated in a homogeneous manner due to end effects; these end effects could reproduced by employing different boundary conditions, but we prefer to ignore them. Note that Kyriakides & Chang[START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF][START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF] use a rigid plug condition on one end, µ(L) = 1, which is not realistic either. Our boundary conditions, sketched in figure1 and provided

wave-vector π/L and 3π/L, which all cancel out upon integration over the domain 0 ≤ Z ≤ L: the solvability condition is automatically verified (this is referred to as a symmetric system in bifurcation theory).

Next, equation [START_REF] Haughton | Bifurcation of inflated circular cylinders of elastic material under axial loading, part i: membrane theory for thin-walled tubes[END_REF] can be solved for µ 2 as µ 2 (Z) = µ 20 +µ 21 cos( π L Z)+µ 22 cos( 2 π L Z). The coefficient µ 21 remains unspecified at this order (in can be used to re-normalize the arc-length η), and the other coefficients are found as

In our notation, any quantity bearing a star is evaluated at the critical point (p , µ 0 (p )).

Next, the solvability condition at order η 3 yields the coefficient p 2 as

24 λ 0 -λ 0,p (n 0,µ ) 2 + n 0,µ λ 0,µ n 0,µ + λ 0 n 0,µ p -λ 0 n 0,µ 2 n 0,p .

The right-hand side is defined in terms of the properties of the homogeneous solution ( §3), and can evaluated numerically for any value of L: recall that all the quantities in the right-hand side are evaluated at (p (L), µ (L) = µ 0 (p (L)), where p (L) is the critical load as determined from the linear buckling analysis, see [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF].

Finally, an expansion of the volume factor v defined in ( 23) is obtained as follows. Observe that the integrand defining v is the function v 0 (p, µ) = µ 2 λ 0 (p, µ): inserting the expansions of p and µ from (28) into v 0 and averaging over Z, one derives an expansion of the volume factor as v(p, µ) = v(p , µ ) + v 2 η 2 + . . . , where the coefficient reads

The right-hand side depends on the properties of the homogeneous solutions, and can be evaluated numerically for any given value of the aspect-ratio L.

When the expansion of v is combined with that of p in (28b), we finally obtain the initial slope of the bifurcated branch as

where p 2 and v 2 have just been calculated. The corresponding tangents are plotted on figure 4(b) for various values of L. They agree very well with the bifurcated branches of the non-linear membrane model.

In figure 4b, it appears that the initial slope of the bifurcated branch changes sign by going through infinity for L/R ≈ 7.7, for both the membrane model and the diffuse interface model. By a classical stability exchange argument, this shows that the initial part of the bifurcated branch is unstable for long balloons, L/R > 7.7 but stable for short balloons, L/R < 7.7, according to both models. A similar fact has been noted in the elastic necking of bars [START_REF] Audoly | Analysis of necking based on a one-dimensional model[END_REF].

To sum up, the diffuse interface model is amenable to a weakly non-linear analysis which reproduces accurately the solutions of the original, fully non-linear membrane model.

Onset of localization: weakly non-linear analysis, infinite length

The domain of validity of the weakly non-linear expansion derived in the previous section is more and more limited when the aspect-ratio gets larger, L/R → ∞: in figure 4b, the domain where the tangent (dashed green line) yields a reasonable approximation to the bifurcated branch (black double-struck curve) shrinks when the aspect-ratio increases from L/R = 10 to 300. This is because for large values of L/R, the extended buckling mode localizes rapidly after bifurcation, a feature not captured by the analysis of the previous section. Here, we derive a different weakly non-linear solution of the diffuse interface model, assuming that the cylinder is infinitely long, L/R = ∞. This solution captures the quick localization of the bulges; it


Special issue on Foundations of Coordination Languages and Self-adaptive Systems

Preface: Special issue on Foundations of Coordination Languages and Self-adaptive Systems

Modern software systems are distributed, concurrent, mobile, and often involve the composition of heterogeneous components and of stand-alone services. Theoretical models, languages and tools for coordinating, composing and adapting services are hence required. They can indeed simplify the development of complex distributed service-based systems, enable functional correctness proofs and improve reusability and maintainability of such systems. This Special Issue (SI) is dedicated to Foundations of Coordination Languages and Self-adaptive Systems. This SI is a follow-up of the FOCLASA workshop that we organized at the 15th International Conference on Software Engineering and Formal Methods (SEFM 2017), held in Trento, Italy on September 4-8, 2017. The workshop received 16 full paper submissions. After a careful reviewing process, the international Program Committee decided to select seven papers for presentation during the workshop and inclusion in the FOCLASA'17 proceedings. From these papers, the six best papers were selected and invited for an extended version to this special issue. All the submissions went through a rigorous peer review process; four papers were finally accepted and are included in this special issue. These papers provide key insights on different aspects of coordination languages and selfadaptive systems.

The first paper, "Orchestrating Incomplete TOSCA Applications with Docker" by A. Brogi et al., presents first a TOSCA (Topology and Orchestration Specification for Cloud Applications) representation for multi-component applications, which can be used to specify the components forming an application, the dependencies occurring among them, and the software support that each component requires to effectively run. Second, the authors present TosKeriser, a tool that automatically completes TOSCA application specifications, by discovering and including Docker-based runtime environments providing the software support needed by the application components. This tool was applied successfully on two concrete case studies.

The second paper, "Towards Attack-Resistant Aggregate Computing Using Trust Mechanisms" by Casadei et al., focuses on security aspects in the context of Aggregate Computing. This paper has focused on attacks based on the diffusion of well-formed, factitious messages and the use of trust mechanisms was proposed to make algorithms resistant to them. After the sensitivity analysis of the proposed trust metrics, several experiments have been executed to verify the ability of an extended gradient algorithm adopting trust mechanisms to tolerate attacks issued by one or more fake nodes emitting random values capable of distorting the ideal gradient field. The investigation has also been complemented with an analysis of the effect, in terms of error and convergence, of some key influencing factors and parameters, for both plain trust and trust with recommendations.

The third paper, "Reasoning about Sensing Uncertainty and its Reduction in Decision-Making for Selfadaptation" by Cámara et al., describes an analysis technique based on model checking of stochastic multi-player games that enables one to quantify the benefits in adaptation performance of factoring sensing uncertainty explicitly into decision-making. The novelty of this approach resides in the definition of patterns to systematically encode environment sensing uncertainty in adaptation scenarios, and an analysis method to explore the solution space while factoring in such uncertainty. This ability shows potential for using these models in adaptive systems to make good decisions about when to model uncertainty explicitly, and when to take advantage of uncertainty reduction tactics.

The fourth paper, "Actor-Based Macroscopic Modeling and Simulation for Smart Urban Planning" by De Berardinis et al., focuses on Tangramob, an agent-based framework for simulating the evolution of urban traffic. The authors designed a simplified model architecture of Tangramob in Timed Rebeca, a timed extension of the actor-based modeling language Rebeca that supports timing features, and developed a tool-chain that generates runnable instances of this model starting from the same input files of Tangramob. Running Timed Rebeca models allows users to get an idea of how the mobility initiatives under study affect the traveling experience of commuters. Once a subset of these initiatives is identified according to user's goals and budget, it is possible to simulate them with Tangramob in order to get more detailed results.

Many people have contributed to the success of this special issue. Besides the authors of the papers, we would like to thank both the members of the Program Committee of the FOCLASA'17 workshop and the additional reviewers who kindly agreed to help us with the reviewing of the papers in this special issue. All carried out an excellent job during this demanding process: Pedro Alvarez, Farhad Arbab, Luis Barbosa, Simon Bliudze, Radu Calinescu, Javier Camara, Flavio De Paoli, Schahram Dustdar, Jean-Marie Jacquet, Sung-Shik Jongmans, Nima Kaviani, Narges Khakpour, Natallia Kokash, Michael Lienhardt, Alberto Lluch Lafuente, Sun Meng, Hernan C. Melgratti, Mohammad Mousavi, Ayoub Nouri, Marc Oriol, Pascal Poizat, Jose Proenca, Michael Sheng, Marjan Sirjani, Jacopo Soldani, Carolyn Talcott, Massimo Tivoli, Giuseppe Vizzari, Lina Ye, Gianluigi Zavattaro.

We hope that these papers will bring the reader novel insights and ideas.

The Science of Computer Programming (SCP) Guest Editors: Carlos Canal Universidad de Málaga, Spain Gwen Salaün Université Grenoble Alpes, France E-mail address: gwen.salaun@inria.fr