Abstract

This paper discusses the evaluation of instabilities on the quasi-static equilibrium path of fluid-loaded pre-stretched cylindrical membranes and the switching to a secondary branch at a bifurcation point. The membrane is represented by only the in-plane stress components, for which an incompressible, isotropic hyperelastic material model is assumed. The free inflation problem yields governing equations and boundary conditions, which are discretized by finite differences and solved by a Newton–Raphson method. An incremental arclength-cubic extrapolation method is used to find generalized equilibrium paths, with different parametrizations. Limit points and bifurcation points are observed on the equilibrium path when fluid level is seen as the controlling parameter. An eigen-mode injection method is employed to switch to a secondary equilibrium branch at the bifurcation point. A limit point with respect to fluid level is observed for a partially filled membrane when the aspect ratio (length/radius) is high, whereas for smaller aspect ratios, the limit point with respect to fluid level is observed at over-filling. Pre-stretch is observed to have a stiffening effect in the pre-limit zone and a softening effect in the post-limit zone.

1. Introduction

Fluid-containing closed membranes find wide applications in several engineering and medical fields, such as liquid storage tanks and liquid filled tube tyres of tractors (liquid ballasting). Fluid-loaded membranes as barriers for wave attenuation [1] and to control the flow of fluid in streams by inflatable dams [2] are advanced applications. In the field of medical technology, the study of fluid-filled membranes is of paramount importance for proper understanding of thin-walled capsules, tissues, cells and liposomes [3–5].

Alongside with free inflation [6–10], constrained inflation of membranes exhibits several aspects of interesting behaviour due to material and geometric nonlinearities, constraining properties and contact conditions. The study of constrained inflated membranes is of practical importance in contexts such as thermoforming, blow moulding and balloon angioplasty. These phenomena are studied for different geometries of air-inflated membranes [11–19].

The analysis of large deformations of fluid-filled membranes is reported in [20–23]. More recently, researchers [24] have presented the static and dynamic analysis of an initially stressed fluid-filled cylindrical membrane. They performed an experimental investigation, which was compared with numerical simulations to show the behaviour of cylindrical membranes under the action of internal hydrostatic pressure. Fluid-filled inextensible cylindrical membranes under internal fluid pressure and external weight are studied in [25]. Hsieh et al. [26] studied theoretical and experimental vibrations of an inextensible cylindrical membrane inflated by a liquid. Recently, a static analysis of a cylindrical membrane under the influence of internal and external fluid pressures in the context of inflatable cylindrical rubber dams is proposed in [27].

Although fluid-filled cylindrical membranes have been frequently considered in past research, the effects from pre-stretches, aspect ratios and fluid density on fluid-filled cylindrical membranes have not been fully described. The softening/stiffening phenomena observed in pre-stretched air-inflated membranes [9,18] need to be investigated for finite inflation of cylindrical membranes under the action of hydrostatic pressure. A previous investigation [24] does not consider static over-pressure created after the membrane is completely filled. Even if bifurcation of membranes have been shown experimentally [24], theoretical analysis failed to show limit point and bifurcation instabilities. The geometric, material and force nonlinearities complicate the problems further and response becomes sensitive to parameters like material constants, pre-stretch, aspect ratio and fluid density. These issues motivate us to present a study on instability phenomena associated with fluid-loaded cylindrical membranes.

The stability investigation of inflated membranes is concentrated on three main topics: limit point instability, bifurcation instability and wrinkling. The latter, corresponding to more or less localized regions of compressive stresses, often occur for inflated membrane structures, and can be seen as related to geometric effects. Important contributions on this topic are given in [28–31]. Limit points and bifurcation points are the two major types of critical points in equilibrium paths. Detailed investigations of limit point instabilities, with some experimental results, are found in [6–9,32,33]. Alexander [6] showed an occurrence of a non-spherical mode after limit point pressure experimentally, which is essentially a secondary equilibrium solution. Recently, Eriksson & Nordmark [34] have shown a non-unique response in terms of instabilities of Mooney–Rivlin membranes for specific boundary conditions. At a bifurcation point, several methods have been proposed for calculation of the secondary branches [35–38]. An eigen-mode injection method for calculation of secondary branches proposed by Wagner & Wriggers [36] is used in this paper, but the evaluation of a general equilibrium path [39] is also discussed. The injection method has been used for inflation of two connected rubber balloons [40]. Further results on stability and bifurcation of membranes can be found in [41–46]. Recently, an important study on bifurcation of shells is presented in [47], cf. also references therein.

In this work, finite inflation of pre-stretched cylindrical membranes under the action of internal hydrostatic pressure is studied. The thin homogeneous, isotropic, hyperelastic membrane is modelled as a neo-Hookean solid affected by only in-plane stresses. Only quasi-static equilibrium solutions are considered, where dynamic and thermal effects are neglected. There are two ways to control fluid-loaded cylindrical membranes by the fluid level or by the fluid volume. We have chosen fluid level as controlling parameter for our study as this is well defined in the numerical computations, but corresponding volumes are easily calculated from obtained equilibrium solutions, allowing comparison to experiments with fluid volume as controlling parameters [24].

The equilibrium equations are obtained by a variational formulation for the free inflation problem and discretized by finite differences. Axisymmetric solutions are obtained by a one-dimensional finite difference approach, whereas asymmetric solutions are obtained by a two-dimensional finite difference approach. The resulting nonlinear algebraic equations are solved by a Newton–Raphson method. An incremental arclength-cubic extrapolation method is used to find generalized equilibrium paths. The simulations based on the efficient model proposed here are for a few instances verified by more demanding three-dimensional finite-element-based solutions, the model described in [8].

2. Problem formulation

(a) Kinematics of deformation

Consider an initial stress-free homogeneous and isotropic cylindrical membrane defined by a radius A, length L/δ and constant initial thickness H with rigid discs of radius A attached to the ends as shown in figure 1. The membrane is then subjected to uniform axial traction to pre-stretch it to a new length L, where δ is an axial pre-stretch parameter. The radius A and length L are seen as fixed for all cylindrical membranes considered here, but pre-stretch is variable.

Schematic of a fluid-loaded cylindrical membrane during inflation. (a) Unstretched membrane and (b) pre-stretching and inflation of membrane.

The undeformed configuration of the membrane is represented by the coordinates R, Θ and Z in the radial, circumferential and axial directions, respectively. The membrane is then inflated with a fluid of density ρ up to fluid level zw. Assuming general asymmetric configuration, the radial, circumferential and axial coordinates of a material point on the membrane after inflation may be represented by, respectively, r(Z,Θ), θ(Z,Θ) and z(Z,Θ). As for a general configuration of an inflated membrane, the principal stretch directions are no longer constant and not aligned with the coordinate axes, it is convenient to represent the strain energy for the membrane in terms of the invariants of the Cauchy–Green deformation tensor B.

The three-dimensional metric tensor for the undeformed cylindrical membrane is
G=1000A20001,2.1
while the three-dimensional metric tensor for the deformed cylindrical membrane is
g=r,Z2+r2θ,Z2+z,Z2r,Zr,Θ+r2θ,Zθ,Θ+z,Zz,Θ0r,Zr,Θ+r2θ,Zθ,Θ+z,Zz,Θr,Θ2+r2θ,Θ2+z,Θ2000λ32,2.2
where an index following a comma denotes differentiation. The stretch λ3=h/H is the principal stretch in the thickness direction with h and H the deformed and undeformed thicknesses of the membrane.

The Cauchy–Green deformation tensor for inflated membrane is given by
B=G−1g2.3
with strain invariants given by
I1=Bii=Giigii,I2=BiiBjj−BjiBij2=I12−GijGpqgipgjq2andI3=∣B∣=∣g∣∣G∣2.4
with the summation convention adopted. The principal stretches are then obtained as the square roots of the eigenvalues of B.

3. The variational formulation

(a) Material strain energy

The strain energy density function (per unit undeformed volume) of an isotropic, homogeneous, hyperelastic material may be expressed as E^=E^(I1,I2,I3). In this work, we consider the neo-Hookean constitutive model for which the strain energy density function is given by E^=(μ/2)(I1−3), where μ is the initial shear modulus. The total elastic potential can then be written as
WM=∫0L/δ∫02πHμ2(I1−3)AdΘdZ,
referring to initial material volume and where the incompressibility assumption is related to the metric tensor by the relations,
I3=1,∣G∣=∣g∣,λ32=G11G22G33(g11g22−g21g12).3.1

(b) Potential energy of the fluid

The membrane is loaded with a fluid of density ρ up to the level zw. The potential energy of the fluid may then be written as −∫PdV, where P is a hydrostatic pressure. The volume element dV in the spatial coordinate system is dV =r2/2 dθ dz and can be written in material coordinates as dV =r2/2 (θ,Θz,Z−θ,Zz,Θ) dΘ dZ.

The potential energy of the fluid for a general configuration can be written as
WP=−∫0L/δ∫02πρg(zw−z)H(zw−z)r22(θ,Θz,Z−θ,Zz,Θ)dΘdZ,3.2
where the Heaviside step function H is used to limit the potential energy of fluid to only the fluid-filled region with z≤zw and g is the gravitational constant, here chosen as g=9.81 m s−2.

(c) Equations of equilibrium

The total potential energy of the system is given by
Π=WM+WP=∫0L/δ∫02πLdΘdZ,3.3
where L is the Lagrangian which is given in terms of three field variables r(Z,Θ), θ(Z,Θ) and z(Z,Θ), and also the fluid level zw. Thus, the Euler–Lagrange equations of equilibrium can be written formally as
∂L∂r−∂∂Z∂L∂r,Z−∂∂Θ∂L∂r,Θ=0,∂L∂θ−∂∂Z∂L∂θ,Z−∂∂Θ∂L∂θ,Θ=0and∂L∂z−∂∂Z∂L∂z,Z−∂∂Θ∂L∂z,Θ=0.3.4
These are three coupled second-order nonlinear partial differential equations. The corresponding boundary conditions for the problem are obtained as
r(0,Θ)=A,θ(0,Θ)=Θ,z(0,Θ)=0,rLδ,Θ=A,θLδ,Θ=ΘandzLδ,Θ=L.3.5

4. Numerical formulation for axisymmetric solution

On the primary equilibrium path, the membrane deforms axisymmetrically, and the circumferential field variable becomes θ=Θ. The potential energy Π in equation (3.3) is thereby stated in terms of field variables r and z, which are functions of the independent axial coordinate Z.

The boundary value problem described by equations (3.4) and (3.5) is discretized by one-dimensional finite differences. The independent coordinate Z is discretized with m intervals in [0,L/δ] along the axial direction, which gives m−1 internal nodal points. The first and second derivatives of r and z are approximated from the discrete values as
r,Z=ri+1−ri−12ΔZ,z,Z=zi+1−zi−12ΔZ4.1
and
(r,Z),Z=ri+1−2ri+ri−1ΔZ2,(z,Z),Z=zi+1−2zi+zi−1ΔZ2,4.2
where ΔZ=(L/δ)/m. The governing equations give 2(m−1) algebraic equations after discretization. Equations (3.4) become a set of algebraic equations in the radial and axial directions, respectively,
F1i(ri−1,ri,ri+1,zi−1,zi,zi+1)=0foralli∈[2,m]andF2i(ri−1,ri,ri+1,zi−1,zi,zi+1)=0foralli∈[2,m]4.3
with boundary conditions,
r1=A,rm+1=A,z1=0andzm+1=L.4.4
The nonlinear set of algebraic equations (4.3) is directly solved by the Newton–Raphson method. For inflation of hyperelastic membranes, limit points are encountered when the fluid level zw is the controlling parameter. This makes convergence difficult at limit points with respect to zw, where the loading characteristics change. A suitably close initial guess is required for the Newton–Raphson method to converge. In order to find this and an automatic tracking of loading characteristics, we have used a combination of an incremental arc-length method with a cubic extrapolation technique, [48]. Let X0, X1, X2 and X3 be known solution vectors fulfilling equilibrium, where Xk={(zw)k,(r2)k,…,(rm)k,(z2)k,…,(zm)k} and k is the index of the equilibrium state. An initial guess vector X4 for the next solution can be found from a cubic parametric fitting
X4=∑k=03∏j=0,j≠k3y4−yjyk−yjXk,4.5
where y is a defined arc length measuring function, for which the four previous arc length are y0=0,y1=A1,y2=y1+A2,y3=y2+A3, and the desired new arc measure is y4=y3+ΔA. We choose ΔA=A3 and Ai given by the Euclidean norm of differences as
Ak=|Xk−Xk−1|,4.6
representing a secant arc-length in the parametric space. With this initial guess vector X4, the solution for fixed-fluid level (zw)4 is obtained by the Newton–Raphson method even if the fluid-level derivative changes sign in the new step. The above algorithm is repeated for each new equilibrium state and is found to be sufficiently robust to cross the limit points automatically, even if more elaborate parametrization of the path can be used, [48]. It is noted that the above method will be able to pass limit points by a finite step, but is not able to isolate and identify any critical solutions on the path. The accurate treatment of these is described below.

5. Eigenvalue problems and symmetry

(a) Primary branch

The cylindrical membrane has C∞v symmetry, and also the primary equilibrium solution has this symmetry. The symmetry ensures that the system remains unchanged under rotations about the Z-axis and reflections about planes passing through the Z-axis. We can show that it is possible to decompose the eigenspaces around primary solutions into one- and two-dimensional subspaces, thus we can see that some eigenvalues, that correspond to a one-dimensional subspace are single, whereas eigenvalues corresponding to two-dimensional subspace are double. The basis for the eigenspaces consists of

— Single eigenvector withC∞vsymmetry. A zero eigenvalue typically corresponds to a limit point and the solution still remains on the same branch.

— Single eigenvector withC∞symmetry. The C∞ symmetry ensures that the solution remains unchanged under rotations but not under reflections. A zero eigenvalue typically corresponds to a pitchfork bifurcation, and indicates some kind of twisting instability which is not observed for our system.

— Double eigenvectors withCnvsymmetry. The Cnv symmetry ensures that the solution remains unchanged under rotations by multiples of 2π/n, and certain reflections by plane angles whose differences are multiples of 2π/n. A (double) zero eigenvalue will typically produce a bifurcation point, with secondary branches being a family of solutions differing by rotations.

(b) Secondary branch

On secondary equilibrium branches, the solutions are no longer axisymmetric. They, however, maintain symmetry of type C∞ or Cnv. The previous experimental investigations [24] shows a bifurcation of the membrane where it maintains C1v symmetry, so we are also looking for the same by taking a symmetry plane passing through Θ=0. We can decompose the eigenspace of secondary solution into one-dimensional subspaces, thus we can see that eigenvalues corresponding to one-dimensional subspaces are single. The basis for the eigenspace consists of

— Single eigenvector with C1v symmetry, where perturbation field variables u, v and w are even, odd and even functions of Θ, respectively. It gives a single zero eigenvalue characteristic for a limit point. We have restricted ourselves for this kind of eigenvectors, which maintains C1v symmetry, by imposing special boundary conditions.

— Single eigenvector which breaks C1v symmetry, where perturbation field variables u, v and w are odd, even and odd functions of Θ, respectively. It always gives a zero eigenvalue corresponding to the infinitesimal rotated solution, and a different zero eigenvalue corresponding to pitchfork bifurcation. We have not considered these types of eigenvectors for our analysis.

6. Stability of primary branch and critical points

As mentioned in the Introduction, the equilibrium paths of fluid-filled membrane show critical points where the behaviour of the membrane changes. The two types of critical points identified on an equilibrium path are limit points, where the slope of the loading curve will change, and bifurcation points, where new equilibrium paths emerge. The axisymmetric equilibrium path is computed by using a numerical formulation described in §4. The primary axisymmetric equilibrium path is then complemented by a secondary path emerging from a bifurcation point. To find these critical points on the primary equilibrium path, we statically perturbed the equilibrium solutions.

(a) Detection of limit point and stability of primary branch

In order to find a limit point on the primary path, we have to perturb the inflated shape, but only axisymmetrically. The solution on the primary path is perturbed with field variables u(Z), w(Z). Let the inflated shape of a spatially perturbed configuration about equilibrium be given by
r(Z)=r0(Z)+ϵu(Z)+O(ϵ2)andz(Z)=z0(Z)+ϵw(Z)+O(ϵ2),6.1
where ϵ is perturbation magnitude, while r0(Z) and z0(Z) are deformed field variables taken from an equilibrium solution. The chosen perturbation field is of C∞v symmetry.

The Lagrangian L in equation (3.3) is expanded in terms of field variables up to quadratic order. The perturbed Lagrangian L^ is the term with coefficient of ϵ2 from the expanded Lagrangian L. The linearized Euler–Lagrangian equations are then obtained in terms of field variables u and w and written as
∂L^∂u−∂∂Z∂L^∂u,Z=0and∂L^∂w−∂∂Z∂L^∂w,Z=0.6.2
The differential equations in u and w are then discretized by one-dimensional finite differences and with introduced boundary conditions u(0)=u(L/δ)=w(0)=w(L/δ)=0. The set of linear algebraic equations for this perturbation is, with a similar discretization as above, written as
[K]{u2,u3,u4,…,w2,w3,w4,…,}T={0}T,6.3
where the subscripts represent nodal numbers and [K] is the perturbed structural stiffness matrix. The eigenvalues of this matrix will decide a limit point and the stability of the primary branch with respect to axisymmetric perturbations.

The conclusions are as follows:

— if all eigenvalues are positive: the equilibrium state is stable for assumed perturbations;

— if one eigenvalue is zero: the equilibrium state is most probably a limit point with respect to fluid level; and

— if at least one eigenvalue is negative: the equilibrium state is unstable for assumed perturbations.

(b) Detection of bifurcation point and branch switching

The perturbation field for the detection of limit point and subsequent determination of stability of the primary branch is shown in §6a. Here, we perturb the equilibrium shape with new field variables u(Z), v(Z) and w(Z), with the help of sinusoidal functions (sinnθ,cosnθ) along the circumferential direction, which can give the bifurcation points on the primary branch. As the primary branch consists of axisymmetric solutions, we have to perturb the inflated shape asymmetrically in order to find the bifurcation point.
r(Z,Θ)=r0(Z)+ϵ1u(Z)cosnΘ+ϵ2u(Z)sinnΘ+O(ϵ12)+O(ϵ22),θ(Z,Θ)=Θ+ϵ1v(Z)sinnΘ−ϵ2v(Z)cosnΘ+O(ϵ12)+O(ϵ22)andz(Z,Θ)=z0(Z)+ϵ1w(Z)cosnΘ+ϵ2w(Z)sinnΘ+O(ϵ12)+O(ϵ22),6.4
where ϵ is a perturbation magnitude, while r0(Z) and z0(Z) are again deformed field variables taken from the primary equilibrium solution. The chosen perturbation field is a two-dimensional subspace of the primary equilibrium solutions and produces double eigenvalues with eigenvectors that have Cnv symmetry. The pair of eigenspace basis vectors {u(Z)cosnΘ,v(Z)sinnΘ,w(Z)cosnΘ} and {u(Z)sinnΘ,−v(Z)cosnΘ,w(Z)sinnΘ} give the same eigenvalue, which has reflection symmetry about a vertical plane passing through Θ=0 and Θ=π/2n. The linear combination of these eigenvectors can give another eigenvector which is also a solution and having a different reflection symmetry plane, but all are equivalent solutions. It is thus sufficient to study one eigenvector by selecting a reflection symmetry plane passing through Θ=0 to π. Here, we are only looking for bifurcation which maintains C1v symmetry. So we modify the above perturbation field, and the inflated shape of an asymmetrically perturbed configuration about equilibrium with n=1 can be written as
r(Z,Θ)=r0(Z)+ϵu(Z)cosΘ+O(ϵ2),θ(Z,Θ)=Θ+ϵv(Z)sinΘ+O(ϵ2)andz(Z,Θ)=z0(Z)+ϵw(Z)cosΘ+O(ϵ2).6.5
The chosen perturbation field has the lowest possible angle dependency. The modified perturbation field gives only one eigenvalue zero at a bifurcation point. From the experimental observations [24], the chosen perturbation field is sufficient to predict a bifurcation point accurately.

The Lagrangian L in equation (3.3) is expanded in terms of field variables up to quadratic order. The perturbed Lagrangian L^ is the term with coefficient of ϵ2 from the expanded Lagrangian L and integrated over Θ. The linearized Euler–Lagrangian equations are obtained in terms of the field variables u, v and w and written as
∂L^∂u−∂∂Z∂L^∂u,Z=0,∂L^∂v−∂∂Z∂L^∂v,Z=0and∂L^∂w−∂∂Z∂L^∂w,Z=0.6.6
The differential equations are again discretized by one-dimensional finite difference expressions, introducing boundary conditions u(0)=u(L/δ)=v(0)=v(L/δ)=w(0)=w(L/δ)=0. The set of linear algebraic equations is written as
[K]{u2,u3,…,v2,v3,…,w2,w3,…,}T={0}T,6.7
where [K] is the perturbed structural stiffness matrix. The eigenvalues of this matrix will indicate a bifurcation point if at least one eigenvalue is zero. If at least one eigenvalue is negative, the equilibrium state is unstable to additional deflections within the assumed perturbation space.

In case of a bifurcation point, a secondary branch will emerge from the primary branch at this point, on which the inflated shape is no longer axisymmetric. In order to obtain equilibrium solutions on this secondary branch, a general asymmetric kinematic description is needed together with a numerical solution procedure described below. To initiate the secondary branch, one needs an asymmetric initial guess to the direction of the secondary branch. For this, we employ an eigen-mode injection method proposed in [36] and used by researchers [37,49]. The eigenvector corresponding to a single zero eigenvalue at the bifurcation point indicates the direction of a secondary branch. The new perturbed field variables V¯={u¯,v¯,w¯}T are given by
V¯=βψ∥ψ∥,6.8
where ψ is the eigenvector of the perturbed structural stiffness matrix [K] corresponding to the zero eigenvalue and β a scaling parameter, chosen to give a robust asymmetric initial guess. The asymmetric initial guess in the direction of the secondary branch for three field variables r, θ and z are
r(Z,Θ)=r0(z)+u¯(Z)cosΘ,θ(Z,Θ)=Θ+v¯(Z)sinΘandz(Z,Θ)=z0(Z)+w¯(Z)cosΘ.6.9

7. Numerical formulation for asymmetric solution

On the secondary branch, the membrane deforms asymmetrically. The field variables are r, θ and z, which are functions of the independent axial and circumferential material coordinates Z and Θ. Two-dimensional finite differences are used to discretize the governing equations in the axial and circumferential directions. The resulting set of highly nonlinear and coupled algebraic equations are solved by a Newton–Raphson method. An incremental arclength-cubic extrapolation method is employed to find the generalized equilibrium paths, as described above.

The imposed boundary conditions ensure that r, θ and z remain even, odd and even functions of Θ, respectively, which maintains C1v symmetry of secondary equilibrium solutions. There are an infinite number of bifurcated solutions possible along the circumferential direction. Here, we fixed our secondary branch by considering a symmetry plane passing through Θ=0 to π. The boundary-value problem described by equations (3.4) and (3.5) was discretized by finite differences. The axial independent coordinate Z was discretized with m intervals along the axial direction and the circumferential independent coordinate Θ with n intervals along the circumferential direction ranging from Θ=0 to Θ=π. The number of internal node points thereby were (m−1)(n+1) in addition to 2(n+1) boundary node points. The derivatives of r, θ and z were approximated from the discrete values as, for example,
r,Z=ri+1,j−ri−1,j2ΔZ,(r,Θ),Θ=ri,j+1−2ri,j+ri,j−1ΔΘ2and(r,Z),Θ=ri+1,j+1+ri−1,j−1−ri−1,j+1−ri+1,j−14ΔZΔΘ,7.1
where ΔZ=(L/δ)/m and ΔΘ=π/n. and subscripts i and j represent nodal numbers from i=1 at Z=0 in the axial direction and j=1 at Θ=0 in the circumferential direction. Equations (3.4) become a set of nonlinear coupled algebraic equations in the radial, circumferential and axial directions, respectively,
F1i(ri−1,j−1,…,θi−1,j−1,…,zi−1,j−1,…,zm+1,n+2,zw)=0foralli∈[2,m],j∈[1,n+1],F2i(ri−1,j−1,…,θi−1,j−1,…,zi−1,j−1,…,zm+1,n+2,zw)=0foralli∈[2,m],j∈[1,n+1]andF3i(ri−1,j−1,…,θi−1,j−1,…,zi−1,j−1,…,zm+1,n+2,zw)=0foralli∈[2,m],j∈[1,n+1]7.2
with boundary conditions
r1,j=A,rm+1,j=A,θ1,j=(j−1)ΔΘandθm+1,j=(j−1)ΔΘ,z1,j=0,zm+1,j=L.7.3
Continuity and periodicity conditions also need to be assigned for field variables, according to
ri,0=ri,2,ri,n+2=ri,n,Θi,0=−Θi,2andΘi,n+2=2π−Θi,n,zi,0=zi,2,zi,n+2=zi,n.7.4
After modification of equation (7.2) through the introduction of equations (7.3) and (7.4), the set of nonlinear algebraic equations was directly solved by a Newton–Raphson method, as above.

Although fluid level is used as controlling parameter for the primary path, this cannot be used as algorithmic controlling parameter here, as initial bifurcating solutions on secondary path have almost the same fluid level. Therefore, the radial displacement at (Z=L/(2δ),Θ=0) was used as controlling parameter. One global equation was added, making the fluid level a depending variable of the applied radial displacement rc. The augmenting equation was written
Fg(rc,r(m+2)/2,1)=(r(m+2)/2,1−rc)=0,7.5
where subscript (m+2)/2 of r represents the considered component. An incremental arc-cubic extrapolation method was used to extrapolate new initial guesses and values of the controlling parameter rc as described above. For this case Xk={(r(m+2)/2,1)k,(r)k,…,(θ)k,…,(z)k,…,(zw)k}. The prediction for the first step on the secondary path was given by an eigen-mode injection method as described by equation (6.9) with the increment in the controlling parameter. This choice of a radial displacement at the midsection as the controlling step parameter has been found fully reliable, even if it is not the fastest growing displacement component during buckling which would be found further down in the cylindrical membrane.

8. Stability of secondary branch

The primary branch was obtained from the method presented in §4 and its stability was checked with the help of the perturbation technique described in §6. The asymmetric configuration of the membrane on a secondary branch was obtained from the method presented in §7. This section describes a perturbation method to ascertain stability of the secondary branch.

Here, we perturb the asymmetric equilibrium shape with new field variables u(Z,Θ), v(Z,Θ) and w(Z,Θ), where they are even, odd and even functions of Θ, respectively. The inflated shape of a perturbed configuration about an asymmetric equilibrium is thereby given by
r(Z,Θ)=r0(Z,Θ)+ϵu(Z,Θ)+O(ϵ2),θ(Z,Θ)=θ0(Z,Θ)+ϵv(Z,Θ)+O(ϵ2)andz(Z,Θ)=z0(Z,Θ)+ϵw(Z,Θ)+O(ϵ2),8.1
where ϵ is a perturbation magnitude, while r0(Z,Θ), θ0(Z,Θ) and z0(Z,Θ) are deformed field variables taken from the secondary equilibrium solution. They are even, odd and even functions of Θ, respectively. The chosen perturbation field is a one-dimensional subspace of secondary equilibrium solutions. Without loss of generality, we will assume the reflection plane passing through Θ=0 to π. The Lagrangian L in equation (3.3) is expanded in terms of field variables up to quadratic order. The perturbed Lagrangian L^ is the term with coefficient of ϵ2 from the expanded Lagrangian L. The linearized Euler–Lagrangian equations are obtained in terms of the field variables u, v and w and written as
∂L^∂u−∂∂Z∂L^∂u,Z−∂∂Θ∂L^∂u,Θ=0,∂L^∂v−∂∂Z∂L^∂v,Z−∂∂Θ∂L^∂v,Θ=0and∂L^∂w−∂∂Z∂L^∂w,Z−∂∂Θ∂L^∂w,Θ=0.8.2
The differential equations were again discretized by two-dimensional finite differences, as in §7. The boundary conditions were u1,j=um+1,j=v1,j=vm+1,j=w1,j=wm+1,j=0 and continuity-periodicity conditions are ui,0=ui,2, ui,n+2=ui,n, vi,0=−vi,2, vi,n+2=−vi,n, wi,0=wi,2 and wi,n+2=wi,n. The set of linear algebraic equations were written as
[K]{u2,u3,…,v2,v3,…,w2,w3,…,}T={0}T,8.3
where [K] is the perturbed structural stiffness matrix. The eigenvalues of this matrix decide the C1v stability of the secondary branch.

The conclusions are

— if all eigenvalues are positive: the asymmetric equilibrium state is stable;

— if one eigenvalue is zero: the asymmetric equilibrium state is most probably a limit point with respect to fluid level; and

— if at least one eigenvalue is negative: the asymmetric equilibrium state is unstable.

9. Detection of wrinkling

For an isotropic material, the principal stress directions are the same as the principal stretch directions, so the principal stress resultants T1 and T2 for neo-Hookean material can be written as [9],
T1=μHλ1λ2−1λ13λ23andT2=μHλ2λ1−1λ13λ23,9.1
where subscripts 1,2,3 of principal stretches and stresses represents first, second and third principal directions, respectively. For a general configuration of inflated membrane, the first and second principal directions are not aligned to meridional and circumferential directions, but the third principal direction is aligned to the thickness direction. A membrane is in an impending wrinkling state when either T1=0 or T2=0. These conditions in kinematic constraints are written as λ1=1/λ2 or λ2=1/λ1. The results obtained in the present study were monitored for the wrinkling but no wrinkling occurred for the studied geometries.

10. Results and discussion

The general parameters chosen for this study are A=0.02 m, H=0.08 mm and shear modulus μ=0.4225 MPa. For better comparison and understanding of characteristics, all the results are presented in non-dimensional form, with non-dimensional parameters directly written on axis labels. The values of pre-stretch are taken as δ=1.5, δ=2 and δ=3. The final length of the membrane (aspect ratio) is taken as L/A=8 or L/A=24. The densities of fluid are taken as ρ=500 kg m−3, ρ=1000 kg m−3 (water) and ρ=2000 kg m−3. The parameters chosen for the study do not necessarily refer to any practical problem, but are merely chosen as representative values.

(a) Stability of equilibrium branches and critical points

In structural mechanics, a response curve (load–deflection curve) is a graphical representation of an equilibrium path, representing the overall quasi-static behaviour of a structure. For this problem, the fluid level zw is taken as the representative load parameter and the radial displacement rmid at node ((m+2)/2,1) is taken as the representative displacement quantity. The inherently nonlinear response to pressurization can give instabilities in the form of critical points with respect to the loading parameter. Response curves for fluid-loaded neo-Hookean membranes with different aspect ratios are shown in figure 2. It is clear that the tangent to the equilibrium path becomes parallel to the deflection axis at the limit points. At a bifurcation point, the tangent to the equilibrium path is not unique. In the present problem, a bifurcation allows transition from symmetric to asymmetric configuration. The zero eigenvalue of the perturbed stiffness matrix [K] obtained from §6a indicates a limit point, while one of [K] obtained from §6b indicates a bifurcation point.

The primary and secondary equilibrium paths are shown in figure 2, along with critical points with respect to the fluid level. The continuous dark line shows the primary branch, while the dashed line shows the secondary branch emerging from a bifurcation point. The secondary branches shown in figures 2 and 3 are only for positive increments in rmid, but the secondary branch continues in the opposite direction with negative increments in rmid (not shown in the figures), in symmetric solutions but the values for Θ=0 and Θ=π swapped. For the pre-stretch ratios included in figure 2, the structure always shows a limit point before the first bifurcation point on the path, but this is not true for higher or even very low pre-stretch values as observed from simulations done in the three-dimensional finite-element model described in [8].

It is clear from the perturbation analysis that the primary branch before the limit point is stable when fluid level is seen as the controlling parameter. At the limit point, one eigenvalue passes through zero and system goes into the unstable post-limit behaviour. The instability of the secondary branch is ascertained according to §8, as one eigenvalue of the perturbed stiffness matrix remains negative. which indicates that the secondary branch is unstable. Based on the signs of eigenvalues, the bifurcation point is super-critical with respect to fluid level.

The inflation of fluid-loaded cylindrical membranes can be performed by controlling, for example, either fluid level or fluid volume, with level most convenient for simulations and volume for experiments. The stability conclusions are highly dependent on the controlling parameter [8]. To check stability of equilibrium paths and to draw the relevance of present results with experimentally observed results [24], the variation of fluid volume inside membrane (V/A3) with fluid level zw/A is shown in figure 3. No limit point is observed on the primary branch for neo-Hookean membrane when fluid volume is seen as controlling parameter. A bifurcation point is observed on the primary branch from which a secondary branch emerges. The primary branch is stable with respect to fluid volume up to the bifurcation point, as observed experimentally [24]. Our results thereby emphasize the observation noted in the previous experimental literature, [24]. It also provides a method for investigation of critical points and tracking of secondary branches from bifurcation points. It is noted that the individual eigenvalue must be traced along the equilibrium path, not just the stiffness matrix determinant, to allow the stability conclusions.

(b) Stretch-induced softening and stiffening

Figure 3 shows that, at a given fluid level, the volume of fluid decreases as the pre-stretch is increased in the pre-limit zone (stiffening before limit-point) and volume increases with pre-stretch in the post-limit zone (softening past the limit-point). This stretch-induced softening in a post-limit zone is counterintuitive, as we normally expect a membrane to become stiffer with pre-stretch. However, [18] has shown that the pre-stretch induces a softening effect for an air-inflated cylindrical membrane in the pre-limit and post-limit zones due to a decrease in circumferential stress from pre-stretching. From the above results, it is evident that pre-stretch induces a stiffening effect in the pre-limit zone with hydrostatic pressure. The reason is that, as the axial pre-stretch softens the cylindrical membrane in the radial direction, gas has a tendency to inflate the membrane more in the radial direction aiming at a spherical shape, giving a softening effect. With hydrostatic pressure, the membrane deforms more axially (which is the stiffer direction) than in the radial direction, which causes a stiffening effect in the pre-limit zone with pre-stretch.

(c) Effect of fluid densities

The response curves of the pre-stretched cylindrical membrane with different fluid densities, as representative of the parameter variations, are shown in figure 4. As fluid density is halved, the fluid levels at the limit and bifurcation points increase. With increasing fluid density, the limit and bifurcation points occur earlier. From figures 2 and 4, it is clear that the response is highly dependent on pre-stretch, aspect ratio and fluid density. Figure 4 shows the shifting of the limit point fluid level with decrease in fluid density. With the same pre-stretch δ=2, the fluid levels at the limit and bifurcation points increase for a cylindrical membrane with the lighter fluid.

(d) Principal stress resultants

The pre-stretched cylindrical membrane is very sensitive to hydrostatic pressure, which may give instability in the form of wrinkling, dependent on pre-stretch, boundary conditions and loading parameters. It is necessary to evaluate the principal stress resultants of the deformed membranes.

The variations of principal stress resultants T1/μH and T2/μH along the axial coordinate for axisymmetric and asymmetric configurations are shown in figures 5 and 6, respectively. Figure 5 shows that the circumferential stress resultant increases as fluid level increases, while the meridional principal stress resultant decreases. The reason behind this is that the circumferential stretch increases faster than the meridional stretch with increase in fluid level, and that this causes reduction in the meridional stress resultant equation (9.1). The meridional stress resultant T1/μH is always greater than zero for the deformed membrane. The circumferential stress resultant T2/μH is always greater than zero for a cylindrical membrane for aspect ratio L/A=8. However, from figure 5b, the circumferential stress resultants near the center and the upper end remain constant and low for axisymmetric deformation, but they never become negative, indicating that no wrinkling occurs for the studied geometries. The variations of principal stress resultants for the asymmetric configuration are different, as stresses are no longer constant along the circumferential direction (figure 6). The circumferential stress resultant is constant and low at the centre and the upper end. It is observed that with addition of fluid, meridional stretch λ1, circumferential stretch λ2, and circumferential stress resultant T2/μH increase, while the meridional stress resultant T1/μH decreases in the pre-stretched cylindrical membrane. The possibility of wrinkling thereby diminishes for the computed pressure range and chosen study parameters.

Variation of principal stress resultants Ti/μH with axial coordinate Z/A for fluid-filled pre-stretched cylindrical membrane. (a) L/A=8, ρgA2/μH=0.1160 and (b) L/A=24, ρgA2/μH=0.1160. The instruction to read this figure is, the farthest curve from the fluid-level nomenclature is related to the first value of fluid level. (Online version in colour.)

(e) Membrane profiles

On the primary branch, the membrane deforms axisymmetrically, while it breaks the symmetry on the secondary branch. Three-dimensional profiles of membranes on the primary and secondary branches are shown in figures 7 and 8. The profile of the membrane on the primary branch before and after limit point are shown in figure 7 at a specific fluid level. Before the limit point, the membrane tends to deform uniformly, figure 7a. After the limit point as fluid level recedes, the deformation of the membrane near the upper end reduces and the membrane deformation is more focused to the bottom, figure 7b. After the limit point, localization effects are predominant. The post-bifurcation shape of the membrane at zw/A=7.336 for a cylindrical membrane with aspect ratio L/A=8 and pre-stretch δ=2 is shown in figure 8a. This clearly shows that the membrane bulges downwards after bifurcation and, as no contact condition is introduced at Z=0, even goes below the lower rigid ring. The increased radial displacement, as discussed above, occurs primarily at Θ=0. The post-bifurcation shape of the membrane at fluid level zw/A=7.779 with aspect ratio L/A=24 is shown in figure 8b, with the same observations as above. This asymmetric equilibrium configuration is unstable with fluid level seen as controlling parameter.

11. Concluding remarks

An instability investigation on fluid-containing cylindrical membranes has been presented. The effects of aspect ratio, pre-stretch and fluid density on the deformation field of the membrane are shown. The important conclusions reported in this paper are as follows:

— thin hyperelastic membranes can show several forms of instabilities when subjected to an internal hydrostatic pressure of variable intensity;

— it is necessary to qualify the occurrence of critical points with respect to the chosen controlling parameter, as the stability conclusions are highly dependent on the precise control. The numerical investigation in this paper is carried out by considering fluid level as primary controlling parameter. Limit points and bifurcation points with respect to fluid level are detected;

— the perturbation technique with fluid level as controlling parameter shows that a pre-limit primary branch is stable and post-limit primary branch is unstable. The secondary branch starting from a bifurcation point is non-axisymmetric and unstable. Singularity of the perturbed stiffness matrix indicates the occurrence of a critical point on the equilibrium branch;

— occurrences of limit and bifurcation points depend on pre-stretch and fluid density in free inflation. Variations of fluid density shift the limit and bifurcation points;

— pre-stretch produces a stiffening effect in the pre-limit zone and a softening effect in the post-limit zone for a neo-Hookean material model;

— as principal stress resultants do not become negative, wrinkling does not occur for a fluid-loaded pre-stretched cylindrical membrane in the considered parameters and pressure range;

— with shorter final length L/A (smaller aspect ratio), the cylindrical membrane shows limit points with respect to fluid level when the membrane is filled above its height, while for larger aspect ratios the limit point occurs for partial filling at small to moderate pre-stretch ratios;

— on the secondary equilibrium branch for a neo-Hookean membrane, the localized deformation is predominant, which requires dense and non-uniform meshing for accurate results; and

— although not shown here, additional simulations have shown that a hardening parameter in a constitutive material model, e.g. a non-zero C2 coefficient in a Mooney–Rivlin model, will quantitatively and qualitatively change the response to loading.

It is believed that the mentioned conclusions are at least qualitatively relevant for studying engineering applications of fluid-loaded cylindrical membranes as well as for studying various geometries found in living organisms. The instability investigation on a cylindrical membrane with different boundary conditions and control parameters will be a future direction of studies.

Funding statement

We greatly acknowledge funding from KTH, Stockholm.

Author contributions

A.P. built up the numerical model in Mathematica and performed simulations with feedback from A.N. and A.E. A.N. verified some results in Comsol Multiphysics and A.E. in an independent FEM model in Matlab. All authors have contributed in drafting manuscript.