1 Introduction

Bulging of a pressurized cylinder has been studied in different works such as Fu et al. [1, 2] and Merodio and collaborators [3,4,5,6,7] in the context of volume-preserving material behavior. While Kyriakides and Chang [8] and Kyriakides [9] show that bulging can be the result of a limit load instability, Fu et al. [2] illustrate that localized bulging can occur even when there is no maximum pressure.

Formation of a localized bulge is a purely non-linear phenomenon of bifurcation. For isotropic elastomer tubes with one end closed, when air is forced into the tube from the other end, a necessary condition for bulging is that the pressure-change in volume response of a section of the constrained tube being inflated uniformly has an up–down behavior (see, e.g., [9, 10]). The bulge then grows radially unless the mechanical response recovers to a second stable branch, i.e., unless the pressure-change in volume response has an up–down–up behavior which is associated with axial propagation of bulging. The propagation of bulging into axial direction can be subdivided into two stages. In the first stage, the inner pressure remains fixed during the ensuing propagation of the bulging instability mode beyond the onset of bifurcation until a suitable configuration is obtained. Then in the second stage, the inflation pressure has to be increased for axial propagation of bulging to continue. The bulging propagation pressure in the first stage corresponds to the one obtained using the equal-area rule, i.e., the Maxwell (propagation) pressure. In the case of open-end hollow cylinders, axial propagation of bulging also involves the two periods captured for closed-end tubes (see Alhayani et al. [11]). Nevertheless, the so-called propagation pressure (first stage) is not captured with the Maxwell condition and furthermore, bulging can be obtained with no pressure maximum. This article continues the considerations from [12] and it focuses on these aspects in regard to swelling.

Due to different factors, such as injuries, inflammation, and hormonal changes, biological soft tissue may swell. Swelling (and shrinking) of soft tissue is accompanied with changes in the mechanical properties (see, e.g., Guo et al. [13] on swelling of arterial tissue). This change in the tissue volume has been taken into account in different mechanical modeling works, for example in the context of isotropic material behavior (see, e.g., Tsai et al. [14], and Pence & Tsai [15, 16] for homogeneous swelling and van der Sman [17] for inhomogeneous swelling), in the context of anisotropic material behavior due to fibrous constituents in a swellable matrix (see, e.g., Demirkoparan and Pence [18,19,20]). The effect of residual stresses and swelling in soft tissue modeling has been studied in different works such as Lanir [21, 22] and Sorrentino et al. [23]. The recent work by Topol et al. [24] models the interaction between soft tissue swelling of a ground substance material and strain stabilization of collagenous fiber to enzymatic degradation. This effect has been further investigated in the finite element study by Gou et al. [25] for fiber remodeling homeostasis in swelling cervical soft tissue. The recent works by Demirkoparan and Merodio [12, 26] study the swelling-induced bulging in hyperelastic fibrous tubes.

In the framework of the present work, we study bulging initiation and propagation in a tube open at both ends made of a swellable Mooney–Rivlin material. In the swellable ground substance material two families of fibers symmetrically disposed and with the same mechanical properties are embedded. Within a fiber family, the fibers are parallel to each other. The strain energy density function for the fibers is considered in both exponential and quadratic form. The use of exponential form for fibrous component has gained some popularity in the modeling of collagen fibers in arterial tissue. This is due to the fact that the fibers may be wavy when unloaded (see, e.g., the review by Holzapfel and Ogden [27]). This article studies the interplay of swelling, fiber orientation, and the mechanical properties of the constituents on bulging initiation and propagation, and it is organized as follows. In Sect. 2, the mathematical problem is described including the bulging condition for the materials at hand. In Sect. 3, the analysis of bulging initiation is carried out while Sect. 4 deals with axial propagation of bulging. Section 5 gives some final conclusions.

2 Geometric and mechanical conditions

Consider an incompressible elastic body in its undeformed reference configuration \(\mathcal {B}_r\), in which the location of a material point is described by a position vector \({\mathbf{X}}\). In the deformed configuration \(\mathcal {B}\), the location of the same material point is described by the position vector \({\mathbf{x}}\), so that the mapping for this deformation is described by the deformation gradient tensor \({\mathbf{F}}=\partial {\mathbf{x}}/\partial {\mathbf{X}}\). The corresponding left and right Cauchy–Green deformation tensors are \({\mathbf{B}}={\mathbf{F}}{\mathbf{F}}^\mathrm{{T}}\) and \({\mathbf{C}}={\mathbf{F}}^\mathrm{{T}}{\mathbf{F}}\).

Now consider a cylindrical coordinate system that is defined by the three base unit vectors \(\left\{ {\mathbf{E}}_R,{\mathbf{E}}_{\Theta },{\mathbf{E}}_{Z}\right\} \) into the radial, circumferential, and axial directions, respectively. In the undeformed reference configuration \(\mathcal {B}_r\), the location of a hollow cylinder can be described in terms of the cylindrical coordinates \(\left\{ R, \Theta , Z\right\} \) as

$$\begin{aligned} A\le R \le B, \qquad 0\le \Theta <2\pi ,\quad 0\le Z\le L, \end{aligned}$$
(1)

where A and B are the inner and outer radii, and L is the length of the cylinder, which may be finite or infinite. The position vector of a material point \({\mathbf{X}}\) in the undeformed configuration then reads

$$\begin{aligned} {\mathbf{X}}=R{\mathbf{E}}_R(\Theta )+Z{\mathbf{E}}_Z. \end{aligned}$$
(2)

Prior to its bifurcation, the cylinder is subjected to an axial force N, an inflation pressure P, and a swelling field

$$\begin{aligned} v=\text{ det }({\mathbf{F}}), \end{aligned}$$
(3)

which is the ratio of the volume of the swollen material to the volume of the unswollen material. Condition (3) ensures that at each instant of the swelling the material remains incompressible. Any location of the deformed cylinder is described by the position vector

$$\begin{aligned} {\mathbf{x}}=r {\mathbf{e}}_r(\theta )+z{\mathbf{e}}_z,\qquad 0\le \theta<2\pi ,\quad 0\le z <\ell , \end{aligned}$$
(4)

where \((r, \theta =\Theta , z)\) are the cylindrical coordinates, \(\left\{ {\mathbf{e}}_r,{\mathbf{e}}_{\theta },{\mathbf{e}}_{z}\right\} \) are the three base unit vectors of the deformed configuration, and \(\ell \) is the length of the cylinder in the deformed configuration. The coordinates of the deformed configuration are

$$\begin{aligned} r=r(R),\qquad \theta =\Theta ,\quad z=\lambda _zZ, \end{aligned}$$
(5)

where \(\lambda _z\) is the stretch ratio on the axial direction of the cylinder. The azimuthal stretch ratio corresponds to the ratio from the deformed to the undeformed radius, \(\lambda _{\theta }=\lambda =r/R>0\), and the radial stretch then becomes \(\lambda _r=v\lambda ^{-1}\lambda _z^{-1}=\partial r/\partial R\). From (3) it follows that \(v=(\partial r/\partial R)(r/R)\lambda _z\).

Furthermore, in the cylindrical coordinate system the deformation gradient tensor and the left and right Cauchy–Green deformation tensors become

$$\begin{aligned} {\mathbf{F}}&\;=\;&\dfrac{v}{\lambda \lambda _z}{\mathbf{e}}_r\otimes {\mathbf{E}}_R + \lambda {\mathbf{e}}_{\theta }\otimes {\mathbf{E}}_{\Theta } + \lambda _z{\mathbf{e}}_{z}\otimes {\mathbf{E}}_{Z}, \end{aligned}$$
(6a)
$$\begin{aligned} {\mathbf{B}}= & {} {\mathbf{C}}\;=\; \dfrac{v^2}{\lambda ^2\lambda _z^2}{\mathbf{e}}_r\otimes {\mathbf{e}}_r + \lambda ^2{\mathbf{e}}_{\theta }\otimes {\mathbf{e}}_{\theta } + \lambda _z^2{\mathbf{e}}_{z}\otimes {\mathbf{e}}_{z}, \end{aligned}$$
(6b)

respectively.

2.1 Material behavior

The material of the cylinder consists of a swellable isotropic matrix and reinforcing fibers. From the mechanical modeling perspective, the material is treated in a homogenized fashion so that any location \({\mathbf{X}}\) of the cylinder contains both matrix and fiber material. In an incompressible and swellable solid, the Cauchy stress tensor can be expressed as

$$\begin{aligned} \varvec{\sigma }=-q{\mathbf{I}}+\dfrac{1}{v}\dfrac{\partial W}{\partial {\mathbf{F}}}{\mathbf{F}}^\mathrm{{T}}, \end{aligned}$$
(7)

where \({\mathbf{F}}\) is given in (6a), \({\mathbf{I}}\) is the second-order identity tensor and q is a scalar that results from the constraint (3). The mechanical behavior of a material is described in terms of a strain energy density function

$$\begin{aligned} W=W(I_1,I_2, I_4, I_5, I_6, I_7, v), \end{aligned}$$
(8)

where \(I_1\),\(I_2\), \(I_4\), \(I_5\), \(I_6\), and \(I_7\) are invariants defined in what follows.

The three principal invariants of \({\mathbf{C}}\) (and due to the symmetry, also of \({\mathbf{B}}\)) in Eq. (6b) have the forms

$$\begin{aligned} I_1({\mathbf{C}})&\;=\;&\text{ tr }({\mathbf{C}})= \dfrac{v^2}{\lambda ^2\lambda _z^2} + \lambda ^2 + \lambda _z^2, \end{aligned}$$
(9a)
$$\begin{aligned} I_2({\mathbf{C}})&\;=\;&\frac{1}{2}\left[ \left[ \text{ tr }({\mathbf{C}})\right] ^2-\text{ tr }\left( {\mathbf{C}}^2\right) \right] =\dfrac{v^2}{\lambda _z^2}+\dfrac{v^2}{\lambda ^2}+\lambda ^2\lambda _z^2, \end{aligned}$$
(9b)
$$\begin{aligned} I_3({\mathbf{C}})&\;=\;&\text{ det }({\mathbf{C}})=v^2. \end{aligned}$$
(9c)

In modeling of compressible solids, the third invariant represents a change of material volume during the deformation, while in our modeling the material remains incompressible at all stages of the deformation. Hence the volume change in the material is solely due to swelling. This motivates the swelling v instead of \(I_3\) in the argument structure of the strain energy density function (8).

Fig. 1
figure 1

A cylindrical membrane. The material consists of a matrix and fibers. The fibers are arranged in two families of parallel fibers with equal mechanical properties, and these fiber families are in a symmetric and helical arrangement that guarantees the cylinder inflation described by Eq. (6) is an admissible deformation. The unit vectors \(\mathbf{{M}_1}\) and \(\mathbf{{M}_2}\) describe the orientation of the fibers from the two families in the reference configuration

Full size image

The fibers are arranged in two families of parallel fibers within each family with equal mechanical properties, and these fiber families are in a symmetric and helical arrangement that guarantees the cylinder inflation described by Eq. (6). The unit vectors

$$\begin{aligned} \mathbf{{M}_1}=\cos \beta {\mathbf{e}}_z+\sin \beta {\mathbf{e}}_{\theta },\qquad \mathbf{{M}_2}=\cos \beta {\mathbf{e}}_z-\sin \beta {\mathbf{e}}_{\theta } \end{aligned}$$
(10)

describe the orientation of the two fiber families in the reference configuration using the angle \(\beta \) (see Fig. 1). These fibers require further invariants to be introduced and used in the definition of the strain energy density function. In particular, the invariants associated with the orientations \(\mathbf{{M}_1}\) and \(\mathbf{{M}_2}\) are

$$\begin{aligned} I_4&\;=\;&{\mathbf{C}}:\mathbf{{M}_1}\otimes \mathbf{{M}_1} =\lambda _z^2\cos ^2\beta +\lambda ^2\sin ^2\beta , \end{aligned}$$
(11a)
$$\begin{aligned} I_5&\;=\;&{\mathbf{C}}^2:\mathbf{{M}_1}\otimes \mathbf{{M}_1} =\lambda _z^4\cos ^2\beta +\lambda ^4\sin ^2\beta , \end{aligned}$$
(11b)
$$\begin{aligned} I_6&\;=\;&{\mathbf{C}}:\mathbf{{M}_2}\otimes \mathbf{{M}_2} =\lambda _z^2\cos ^2\beta +\lambda ^2\sin ^2\beta , \end{aligned}$$
(11c)
$$\begin{aligned} I_7&\;=\;&{\mathbf{C}}^2:\mathbf{{M}_2}\otimes \mathbf{{M}_2} =\lambda _z^4\cos ^2\beta +\lambda ^4\sin ^2\beta , \end{aligned}$$
(11d)

where the invariants \(I_4\) and \(I_5\) are associated with the first fiber family (the one with orientation angle \(\mathbf{{M}_1}\)) and the invariants \(I_6\) and \(I_7\) are associated with the second fiber family (the one with orientation angle \(\mathbf{{M}_2}\)).

The invariants \(I_4\) and \(I_6\) can be interpreted as the squared magnitude of the fiber stretch ratios, while the invariants \(I_5\) and \(I_7\) are related to fiber shearing. Due to the symmetric arrangement of the fibers in the cylinder that results from (10) we obtain \(I_4 = I_6\) and \(I_5 = I_7\) (see Eq. (11)).

The expressions of the invariants in Eqs. (9) and (11) motivate the reformulation of the strain energy density function (8) in terms of the principal stretches as

$$\begin{aligned} W(I_1,I_2, I_4, I_5, I_6, I_7, v)={\bar{W}}(\lambda ,\lambda _z, v). \end{aligned}$$
(12)

In order to simplify the notation, partial derivatives of \({\bar{W}}\) shall be indicated by the subscripts as \({\bar{W}}_{\lambda _z}=\partial {\bar{W}}/\partial \lambda _z\), \({\bar{W}}_{\lambda _z\lambda _z}=\partial ^2 {\bar{W}}/\partial \lambda _z^2\), etc.

2.2 Bulging bifurcation in the membrane approximation

In the absence of body forces, as it is the case in the present work, the Cauchy stress tensor needs to satisfy \(\text{ div }\,\varvec{\sigma }={\mathbf{0}}\). Let us consider the cylinder in the membrane approximation in which it is taken a unique radius \(R=(A+B)/2\), the midsurface radius, and \(H=(B-A)\ll R\) is the thickness. In this membrane approximation, the relation between the inflation pressure P, the axial normal force N, and the deformation is given by [12]

$$\begin{aligned} P=\dfrac{H}{R}\dfrac{{\bar{W}}_{\lambda }}{\lambda \lambda _z},\qquad N={\bar{W}}_{\lambda _z}. \end{aligned}$$
(13)

The radial normal stresses are negligible in the membrane approximation, hence we take \(\sigma _{rr}=0\). Then the circumferential and axial components \(\sigma _{\theta \theta }\) and \(\sigma _{zz}\) of the Cauchy stress tensor (7) become

$$\begin{aligned} \sigma _{\theta \theta } =\lambda \dfrac{\partial {\bar{W}}}{\partial \lambda },\qquad \sigma _{zz} =\lambda _{z}\dfrac{\partial {\bar{W}}}{\partial \lambda _{z}}. \end{aligned}$$
(14)

In order to investigate the possible bulging bifurcation mode of the cylinder, we consider incremental displacements with respect to the deformed configuration (under equilibrium), of the form

$$\begin{aligned} \delta {\mathbf{u}}=\delta u_r(z){\mathbf{e}}_r+\delta u_z(z){\mathbf{e}}_z, \end{aligned}$$
(15)

in which \(\delta \) is the increment. The bulging condition for a cylinder membrane of radius R and length L is derived in [7] and given by

$$\begin{aligned} f(\bar{W},\lambda ,\lambda _z)+\left( \frac{2\pi R}{L}\right) ^2 \lambda ^2\lambda _z\bar{W}_{\lambda _z}\bar{W}_{\lambda _z\lambda _z}=0, \end{aligned}$$
(16)

where

$$\begin{aligned} f(\bar{W},\lambda ,\lambda _z)=\lambda _z^2\bar{W}_{\lambda _z\lambda _z}( \lambda ^2\bar{W}_{\lambda \lambda }- \lambda \bar{W}_{\lambda })-( \lambda \lambda _z\bar{W}_{\lambda \lambda _z}- \lambda \bar{W}_{\lambda })^2. \end{aligned}$$
(17)

The last term on the left-hand side of (16) vanishes for a cylinder of infinite length, \(L\rightarrow \infty \), which is the focus of this analysis.

3 Bulging for specific material models

Different hyperelastic models have been proposed that account for the mechanical behavior of fibrous biological soft tissue (see, e.g., Chagnon et al. [28]). Let us consider a strain energy density function in the form

$$\begin{aligned} \begin{array}{llllllllllll} W(I_1,I_2, I_4, I_5, v) &{}&{}=\displaystyle \frac{\mu _1}{2}\left( I_1-3v^{-2/3}\right) +\dfrac{\mu _2}{2}\left( I_2-3v^{-4/3}\right) \\ &{}&{}\quad +\, \dfrac{k_1}{{\bar{k}}_1} \left[ \exp \left( {\bar{k}}_1\left[ I_4-1\right] ^2\right) -1\right] +\dfrac{k_2}{{\bar{k}}_2} \left[ \exp \left( {\bar{k}}_2\left[ I_5-1\right] ^2\right) -1\right] . \end{array} \end{aligned}$$
(18)

The first two terms on the right-hand side of (18) define the isotropic behavior of the ground substance material, in which the fibers are embedded, and they describe a swellable modification of the classical Mooney–Rivlin material (see, e.g., [16]). The material parameters \(\mu _1\) and \(\mu _2\) are taken to be independent of the swelling, and we refer to works such as [24, 25] that discuss the changes of the material stiffness with the amount of swelling. The last two terms on the right-hand side of equation (18) account for the anisotropic material behavior due to the two fiber families. The expression given in (18) has already made use of the equal properties and the symmetric arrangement of the fibers that are associated with the first fiber family (and therefore with \(I_4\), \(I_5\)) and with the second fiber family (and therefore with \(I_6\), \(I_7\)), i.e., the following identities are already used in (18)

$$\begin{aligned} \dfrac{k_1}{{\bar{k}}_1} \left[ \exp \left( {\bar{k}}_1\left[ I_4-1\right] ^2\right) -1\right]= & {} \displaystyle \dfrac{k_1}{2{\bar{k}}_1} \sum _{i=4,6} \left[ \exp \left( {\bar{k}}_1\left[ I_i-1\right] ^2\right) -1\right] , \end{aligned}$$
(19a)
$$\begin{aligned} \dfrac{k_2}{{\bar{k}}_2} \left[ \exp \left( {\bar{k}}_2\left[ I_5-1\right] ^2\right) -1\right]= & {} \displaystyle \dfrac{k_2}{2{\bar{k}}_2} \sum _{i=5,7} \left[ \exp \left( {\bar{k}}_2\left[ I_i-1\right] ^2\right) -1\right] . \end{aligned}$$
(19b)

The material parameters \(k_1\) and \(k_2\) can be regarded as relative strength parameters for the fiber stiffness. Moreover \({\bar{k}}_1\) and \({\bar{k}}_2\) are also material parameters which may be regarded as fiber sensitivity parameters. The exponential terms account for fibers in an initially wavy form that have to be extended in order to fully bear loading. Specialization of the strain energy density function (18) has been applied in different works. For example, a non-swellable version of this material model has been applied in the work [29] for modeling of arteries.

For the material model (18) the bulging condition for an infinitely long cylinder (17) can be written as

$$\begin{aligned} f=f_{I_1}+f_{I_2}+f_{{I_4}}+f_{{I_5}}=0. \end{aligned}$$
(20)

The first two terms on the right-hand side of (20) account for the behavior of the matrix. In particular, the first term on the right-hand side of (20) results from the first term on the right-hand side of (18), and the second term of (20) results from the second term on the right-hand side of (18). In particular, for the material model (18), one can find that

$$\begin{aligned} f_{I_1}&\;=\;&\dfrac{\mu _1^2 (- \lambda _z^4 \lambda ^8 + 4 \lambda _z^4 \lambda ^2 v^2 + 6 \lambda _z^2 \lambda ^4 v^2 + 3 v^4)}{(\lambda _z^4 \lambda ^4)}, \end{aligned}$$
(21a)
$$\begin{aligned} f_{I_2}&\;=\;&\dfrac{\mu _2^2 (-\lambda _z^6 \lambda ^8 + 2 \lambda _z^4 \lambda ^4 v^2 - \lambda _z^2 v^4 + 12 \lambda ^2 v^4)}{\lambda _z^2 \lambda ^4}. \end{aligned}$$
(21b)

Moreover the third term in (20) accounts for the exponential response in (18) which is in turn associated with the invariant \(I_4\) (and \(I_6\)). One can also deduce that

$$\begin{aligned} \begin{array}{lllllllllll} f_{I_4}&{}=&{} \Bigg [\lambda _z^2\bigg \{\lambda ^2 \left[ 4 k_1 \sin ^2\beta \left[ I_4 - 1\right] + 8 k_1 \lambda ^2 \sin ^4\beta + 16 {\bar{k}}_1 k_1 \lambda ^2 \sin ^4\beta \left[ I_4 - 1\right] ^2\right] - 4 k_1 \lambda ^2 \sin ^2\beta \left[ I_4 - 1\right] \bigg \}\\ &{}&{}\times \left[ 4 k_1 \cos ^2\beta \left[ I_4 - 1\right] + 8 k_1 \lambda _z^2 \cos ^4\beta + 16 {\bar{k}}_1 k_1 \lambda _z^2 \cos ^4\beta \left[ I_4 - 1\right] ^2\right] \\ &{}&{} -\, \bigg \{\lambda _z \lambda \left[ 8 k_1 \lambda _z \lambda \cos ^2\beta \sin ^2\beta + 16 {\bar{k}}_1 k_1 \lambda _z \lambda \cos ^2\beta \sin ^2\beta \left[ I_4 - 1\right] ^2\right] \\ &{}&{} -\, 4 k_1 \lambda ^2 \sin ^2\beta \left[ I_4 - 1\right] \bigg \}^2\Bigg ] \exp \left( 2 {\bar{k}}_1 \left[ I_4 - 1\right] ^2\right) , \end{array} \end{aligned}$$
(22)

where \(I_4\) is given in terms of the stretches in Eq. (11a). The last term in (20) accounts for the exponential response in (18) that is associated with the invariant \(I_5\) (and \(I_7\)),

$$\begin{aligned} \begin{array}{lllllllllll} f_{I_5}&{}=&{}\Bigg [\lambda _z^2 \bigg \{\lambda ^2 \left[ 32 k_2 \lambda ^6 \sin ^4\beta + 24 k_2 \lambda ^2 \sin ^2\beta \left[ I_5 - 1\right] + 64 {\bar{k}}_2 k_2 \lambda ^6 \sin ^4\beta \left[ I_5 - 1\right] ^2\right] - 8 k_2 \lambda ^4 \sin ^2\beta \left[ I_5 - 1\right] \bigg \}\\ &{}&{}\times \left[ 32 k_2 \lambda _z^6 \cos ^4\beta + 24 k_2 \lambda _z^2 \cos ^2\beta \left[ I_5 - 1\right] + 64 {\bar{k}}_2 k_2 \lambda _z^6 \cos ^4\beta \left[ I_5 - 1\right] ^2\right] \\ &{}&{} -\, \bigg \{\lambda _z \lambda \left[ 32 k_2 \lambda _z^3 \lambda ^3 \cos ^2\beta \sin ^2\beta + 64 {\bar{k}}_2 k_2 \lambda _z^3 \lambda ^3 \cos ^2\beta \sin ^2\beta \left[ I_5 - 1\right] ^2\right] \\ &{}&{} -\, 8 k_2 \lambda ^4 \sin ^2\beta \left[ I_5 - 1\right] \bigg \}^2\Bigg ] \exp \left( 2 {\bar{k}}_2 \left[ I_5 - 1\right] ^2\right) , \end{array} \end{aligned}$$
(23)

where \(I_5\) is given in terms of the stretches in Eq. (11b).

Bulging in the limit \({\bar{k}}_1\rightarrow 0\) and \({\bar{k}}_2\rightarrow 0\): The exponential terms in the strain energy density function reduce to the standard forms in the limits

$$\begin{aligned} \lim _{{\bar{k}}_1\rightarrow 0} \dfrac{k_1}{{\bar{k}}_1} \left[ \exp \left( {\bar{k}}_1\left[ I_4-1\right] ^2\right) -1\right]&\; = \;&k_1 \left[ I_4-1\right] ^2, \end{aligned}$$
(24a)
$$\begin{aligned} \lim _{{\bar{k}}_2\rightarrow 0} \dfrac{k_2}{{\bar{k}}_2} \left[ \exp \left( {\bar{k}}_2\left[ I_5-1\right] ^2\right) -1\right]&\; = \;&k_2 \left[ I_5-1\right] ^2. \end{aligned}$$
(24b)

In these limits, the contributions of the fibers to the bulging condition in (22) and (23) become

$$\begin{aligned} \begin{array}{lllllll} f_{I_4}&{}=&{}\lambda _z^2 \left\{ \lambda ^2 \left[ 4 k_1 \sin ^2\beta \left[ I_4 - 1\right] + 8 k_1 \lambda ^2 \sin ^4\beta \right] - 4 k_1 \lambda ^2 \sin ^2\beta \left[ I_4 - 1\right] \right\} \\ &{}&{}\times \left[ 4 k_1 \cos ^2\beta \left[ I_4 - 1\right] + 8 k_1 \lambda _z^2 \cos ^4\beta \right] \\ &{}&{} - \left[ 4 k_1 \lambda ^2 \sin ^2\beta \left[ I_4 - 1\right] - 8 k_1 \lambda _z^2 \lambda ^2 \cos ^2\beta \sin ^2\beta \right] ^2, \end{array} \end{aligned}$$
(25)

and

$$\begin{aligned} \begin{array}{llllllllllllll} f_{I_5}&{}=&{}\lambda _z^2 \left\{ \lambda ^2 \left[ 32 k_2 \lambda ^6 \sin ^4\beta + 24 k_2 \lambda ^2 \sin ^2\beta \left[ I_5 - 1\right] \right] - 8 k_2 \lambda ^4 \sin ^2\beta \left[ I_5 - 1\right] \right\} \\ &{}&{}\times \left[ 32 k_2 \lambda _z^6 \cos ^4\beta + 24 k_2 \lambda _z^2 \cos ^2\beta \left[ I_5 - 1\right] \right] \\ &{}&{} - \left[ 8 k_2 \lambda ^4 \sin ^2\beta \left[ I_5 - 1\right] - 32 k_2 \lambda _z^4 \lambda ^4 \cos ^2\beta \sin ^2\beta \right] ^2. \end{array} \end{aligned}$$
(26)

The bulging condition function (20) with \(f_{I_1}\) and \(f_{I_2}\) in (21) and \(f_{I_4}\) and \(f_{I_5}\) in (25) and (26) corresponds to a strain energy density function in the form

$$\begin{aligned} \begin{array}{llllllllllll} W(I_1,I_2, I_4, I_5, v) =\displaystyle \frac{\mu _1}{2}\left( I_1-3v^{-2/3}\right) +\dfrac{\mu _2}{2}\left( I_2-3v^{-4/3}\right) +{k_1}\left[ I_4-1\right] ^2 +{k_2}\left[ I_5-1\right] ^2. \end{array} \end{aligned}$$
(27)

Different forms of (27) (for instance with \(\mu _2=k_2=0\), etc) have been used in the literature (see, e g., El Hamdaoui et al. [30, 31], Vinh et al. [32], and Goriely [33]). In what follows bulging for different models is analyzed.

3.1 Mooney–Rivlin behavior

Consider a swellable Mooney–Rivlin material having no fibers, which is described by the strain energy density function (18) with \(k_1=k_2=0\),

$$\begin{aligned} W(I_1,I_2,v)=\frac{\mu _1}{2}\left( I_1-3v^{-2/3}\right) +\dfrac{\mu _2}{2}\left( I_2-3v^{-4/3}\right) . \end{aligned}$$
(28)

This model has been applied in the modeling of the elastinous ground substance behavior of arterial tissue in works such as [34]. While \(\mu _1\) is usually taken to be positive, the parameter \(\mu _2\) can take positive or negative values in the mechanical modeling. In this case, the bulging condition (20) reduces to \(f=f_{I_1}+f_{I_2}=0\). As it is clear from (21a) and (21b), the resulting equation is a polynomial of order 8 in \(\lambda \) and \(\lambda _z\) (and order 4 polynomial in v). Therefore numerical methods are utilized to analyze the bulging condition.

Some results are shown in Fig. 2 for an unswollen material (\(v=1\)). In both panels, the blue curves show the results for \(\mu _2/\mu _1>0\), the red curves show the response for \(\mu _2/\mu _1<0\), and the black curve shows the neo-Hookean response (\(\mu _2=0\)).

Fig. 2
figure 2

Equation (20) \(f=f_{I_1}+f_{I_2}=0\) gives the bulging condition for a Mooney–Rivlin material in the absence of fibers. a shows different curves of values \(f/\mu _1^2\) vs. \(\lambda \) for an axial stretch of \(\lambda _z=1.1\) and different ratios of \(\mu _1\) to \(\mu _2\). In both figures, the blue curves show the results for \(\mu _2/\mu _1>0\), the red curves show the response for \(\mu _2/\mu _1<0\), and the black curve shows the neo-Hookean response (\(\mu _2=0\)). b depicts relations between the circumferential and the axial stretches \(\lambda \) and \(\lambda _z\), respectively, which fulfill the bulging condition \(f_{I_1}+f_{I_2}=0\) for different ratios of \(\mu _1\) to \(\mu _2\)

Full size image

Figure 2a shows different curves \(f/\mu _1^2\) vs. \(\lambda \) for an axial stretch of \(\lambda _z=1.1\) and different ratios of \(\mu _1\) to \(\mu _2\). These curves take larger values for f when \(\mu _2\) increases. The curve for \(\mu _2=-\mu _1\) does not take positive values for \(f/\mu _1^2\). Figure 2b depicts relations between the circumferential and the axial stretches \(\lambda \) and \(\lambda _z\), which fulfill the bulging condition \(f_{I_1}+f_{I_2}=0\) for different ratios of \(\mu _1\) to \(\mu _2\). Notice that for \(\mu _2=-\mu _1\) and for the depicted values for \(\lambda _z\) and \(\lambda \) the bulging condition \(f=0\) is not fulfilled.

3.2 Fibers in neo-Hookean Matrix

In this section, we consider two fiber families that are embedded in a swellable neo-Hookean ground material (\(\mu _2=0\)). Despite its simple form, the neo-Hookean model is often sufficient in describing the mechanical behavior of a ground material, in which the fibers are embedded. In the following parts, we study the bulging behavior of fibers with a mechanical behavior described by either \(I_4\) or \(I_5\).

3.2.1 Fiber behavior in terms of \(I_4\) (\(k_1\ne 0\), \(k_2=0\) in (18))

Let us consider a neo-Hookean matrix (\(\mu _2=0\)) with embedded fibers and \(k_1\ne 0\), \(k_2=0\), for which the strain energy density (18) becomes

$$\begin{aligned} \begin{array}{llllllllllll} W(I_1, I_4, v)= & {} \displaystyle \frac{\mu _1}{2}\left( I_1-3v^{-2/3}\right)+ & {} \dfrac{k_1}{{\bar{k}}_1} \left[ \exp \left( {\bar{k}}_1\left[ I_4-1\right] ^2\right) -1\right] . \end{array} \end{aligned}$$
(29)

The contribution of the fibers to the strain energy density function is described by the exponential term and the invariant \(I_4\). This model has gained some popularity in the modeling of mechanical behavior of arteries (see, e.g., [27, 35]). The exponential form is motivated by the observation that many biological materials have a common “exponential-shaped” stress–strain curve [36]. Some aspects of the bulging bifurcation for (29) have been studied in [26]. This article continues that investigation of the bulging bifurcation, which will then also serve as a reference for the other models analyzed in what follows in this article.

For the strain energy density function (29), (20) becomes \(f_{I_1}+f_{I_4}=0\), where \(f_{I_1}\) is given by Eq. (21a) and \(f_{I_4}\) is given by Eq. (22).

Fig. 3
figure 3

Bulging for a fiber-reinforced neo-Hookean material (29). a shows the function (20), which is \(f=f_{I_1}+f_{I_4}\) (normalized by the square of the shear modulus \(\mu _1\)), for an axial stretch \(\lambda _z=1.1\), \(\beta =\pi /4\), \(k_1=\mu _1/100\), and for three swelling ratios \(v=1\) (no swelling), \(v=1.1\), and \(v=1.2\). The solid lines represent results for \({\bar{k}_1}\rightarrow 0\) and the dashed lines represent results for \({\bar{k}_1}=0.4\). b shows different combinations for \(\lambda _z\) and \(\lambda \) that fulfill the condition \(f_{I_1}+f_{I_4}=0\) for \(\beta =\pi /4\), \(k_1=\mu _1/100\), and for three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\). The different line styles represent the results for different values of \({\bar{k}_1}\)

Full size image

Some conditions for bulging are studied in the two panels of Fig. 3, in which the results for the three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\) are shown using different line colors.

Figure 3a shows the function (20), which is \(f=f_{I_1}+f_{I_4}\) (normalized by \(\mu _1^2\)), for an axial stretch \(\lambda _z=1.1\), \(\beta =\pi /4\), and \(k_1=\mu _1/100\). Increasing values for \({\bar{k}}_1\) or v lead to larger values of f. The solid lines represent results for \({\bar{k}_1}\rightarrow 0\), which corresponds to the mechanical response of (27) for \(\mu _2=k_2=0\). The dashed lines represent the results for \({\bar{k}_1}=0.4\). When \({\bar{k}_1}=0.4\) bulging may occur for \(v=1\), while if \({\bar{k}_1}=0.4\) and \(v=1.2\) bulging will not occur because f remains positive.

Figure 3b shows different combinations for \(\lambda _z\) and \(\lambda \) that fulfill the condition \(f_{I_1}+f_{I_4}=0\) for \(\beta =\pi /4\) and \(k_1=\mu _1/100\). The different line styles represent the results for different values of \({\bar{k}_1}\). This diagram shows that as \({\bar{k}_1}\rightarrow 0\) the curves are decreasing so that for an axial stretch \(\lambda _z\) we find a unique value of the azimuthal stretch \(\lambda \) on the depicted range. On the other hand, when \({\bar{k}_1}>0\) a value of the axial stretch may have multiple corresponding values of the azimuthal stretch \(\lambda \) associated with bulging.

The examples of Fig. 3 are restricted to a fixed value \(\beta =\pi /4\) and the impact of the fiber orientation on the bulging for a fiber-reinforced neo-Hookean material (29) is studied in Fig. 4.

Fig. 4
figure 4

Impact of the fiber orientation on bulging for a fiber-reinforced neo-Hookean material (29). a shows the function (20) \(f=f_{I_1}+f_{I_4}\) (normalized by the square of the shear modulus \(\mu _1\)) for an axial stretch \(\lambda _z=1.1\), \(k_1=\mu _1/100\), \(v=1\), \({\bar{k}_1}=0.4\), and for different values for \(\beta \) on different ranges for \(f/\mu _1^2\). The polar diagram of b shows different combinations for \(\beta \) (circumferential orientation) and \(\lambda \) (radial values) that fulfill the condition \(f_{I_1}+f_{I_4}=0\) for \(k_1=\mu _1/100\), \({\bar{k}_1}=0.4\), and \(\lambda _z=1\) for the three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\)

Full size image

Figure 4a shows the function (20) for \(f=f_{I_1}+f_{I_4}\) (normalized by \(\mu _1^2\)) for an axial stretch \(\lambda _z=1.1\), \(k_1=\mu _1/100\), \(v=1\), \({\bar{k}_1}=0.4\), and for different values for \(\beta \) on different ranges for \(f/\mu _1^2\). In the case \(\beta =0\) the fibers are parallel to the the axial direction of the cylinder, \(\mathbf{{M}_1}=\mathbf{{M}_2}={\mathbf{e}}_z\), and in the case \(\beta =\pi /2\) the fibers are oriented in the circumferential direction of the cylinder, \(\mathbf{{M}_1}=-\mathbf{{M}_2}={\mathbf{e}}_{\theta }\). The different curves show that if \(\beta \) is sufficiently close to \(\pi /2\), then bulging will not occur.

The polar diagram of Fig. 4b shows different combinations for \(\beta \) (circumferential orientation) and \(\lambda \) (radial values) that fulfill the condition \(f_{I_1}+f_{I_4}=0\) for \(k_1=\mu _1/100\), \({\bar{k}_1}=0.4\), and \(\lambda _z=1\) and the three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\). This diagram shows that the smallest values of \(\lambda \) that fulfill the bulging condition \(f_{I_1}+f_{I_4}=0\) occur for \(\beta =0\) (\(\beta =\pm \pi \)), while the largest values of \(\lambda \) that fulfill this condition occur for \(\beta =\pm \pi /2\). The figure shows symmetries with respect to horizontal and vertical axes of the polar diagram.

3.2.2 Fiber behavior in terms of \(I_5\) (\(k_1=0\), \(k_2\ne 0\) in (18))

Let us now consider a neo-Hookean matrix with embedded fibers and \(k_1= 0\), \(k_2\ne 0\), for which the strain energy density function (18) becomes

$$\begin{aligned} \begin{array}{llllllllllll} W(I_1, I_5, v)= & {} \displaystyle \frac{\mu _1}{2}\left( I_1-3v^{-2/3}\right)+ & {} \dfrac{k_2}{{\bar{k}}_2} \left[ \exp \left( {\bar{k}}_2\left[ I_5-1\right] ^2\right) -1\right] . \end{array} \end{aligned}$$
(30)

In this model, the contribution of the fibers to the strain energy density function is described in exponential form by the invariant \(I_5\). For this strain energy density function the bulging condition (20) becomes \(f=f_{I_1}+f_{I_5}=0\), where \(f_{I_1}\) is given in Eq. (21a) and \(f_{I_5}\) is given in Eq. (22).

Fig. 5
figure 5

Bulging for a fiber-reinforced neo-Hookean material (30). a shows the function (20) for \(f=f_{I_1}+f_{I_5}\) (normalized by \(\mu _1^2\)) for an axial stretch of \(\lambda _z=1.1\), \(\beta =\pi /4\), \(k_2=\mu _1/100\), and for the three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\). The solid lines represent results for \({\bar{k}_2}\rightarrow 0\), and the dashed lines represent the results for \({\bar{k}_2}=0.001\)

Full size image

Figure 5 shows the bulging condition for three swelling ratios, which are highlighted by different line colors. In particular, Fig. 5a shows the function (20) for \(f=f_{I_1}+f_{I_5}\) (normalized by \(\mu _1^2\)) for an axial stretch of \(\lambda _z=1.1\), \(\beta =\pi /4\), and \(k_1=\mu _1/100\). The solid lines represent results for \({\bar{k}_2}\rightarrow 0\), and the dashed lines represent the results for \({\bar{k}_2}=0.001\).

Figure 5b shows different combinations for \(\lambda _z\) and \(\lambda \) that fulfill the condition \(f_{I_1}+f_{I_5}=0\) for \(\beta =\pi /4\), \(k_2=\mu _1/100\), \({\bar{k}_2}\rightarrow 0\), and for the three swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\).

Figure 6 depicts the impact of fiber orientation on bulging for the fiber-reinforced neo-Hookean material (30).

Fig. 6
figure 6

Impact of fiber orientation on bulging for a fiber-reinforced neo-Hookean material (30). a Shows the function (20) for \(f=f_{I_1}+f_{I_5}\) (normalized by the square of the shear modulus \(\mu _1\)) for an axial stretch of \(\lambda _z=1.1\), \(k_2=\mu _1/100\), \(v=1\), \({\bar{k}_2}=0\), and for different values for \(\beta \). The polar diagram of b shows different combinations for \(\beta \) (circumferential orientation) and \(\lambda \) (radial values) that fulfill the condition \(f_{I_1}+f_{I_5}=0\) for \(\lambda _z=1.00\), \(k_2=\mu _1/100\), \({\bar{k}_1}=0.0\), and for the two swelling ratios \(v=1\) and \(v=1.2\)

Full size image

Figure 6a shows the function (20) for \(f=f_{I_1}+f_{I_5}\) (normalized by the square of the shear modulus \(\mu _1\)) for an axial stretch of \(\lambda _z=1.1\), \(k_2=\mu _1/100\), \(v=1\), \({\bar{k}_2}=0\), and for different values \(\beta \).

The polar diagram of Fig. 6b shows different combinations for \(\beta \) (circumferential orientation) and \(\lambda \) (radial values) that fulfill the condition \(f_{I_1}+f_{I_5}=0\) for \(\lambda _z=1.00\), \(k_2=\mu _1/100\), \({\bar{k}_2}=0.0\), and two swelling ratios \(v=1\) and \(v=1.2\). This diagram shows the following for values \(\pi /2 \ge \beta \ge 0\). As \(\beta \) increases from \(\beta =0\), the values of \(\lambda \) associated with bulging increase until \(\beta \) reaches an angle close to \(\pi /6\) from which as \(\beta \) increases the values of \(\lambda \) associated with bulging decrease.

3.2.3 Fiber behavior in terms of both \(I_4\) and \(I_5\)

As a last example of this section, we now consider a swellable neo-Hookean matrix containing two families of symmetric fibers, for which the strain energy density is given by (27) with \(\mu _2=0\). In this case, the strain energy density function is formulated in terms of both \(I_4\) and \(I_5\).

Figure 7 shows combinations for \(k_1\ge 0\) and \(k_2\ge 0\) that fulfill the bulging condition (20) for \(\lambda =2\) with \(f_{I_1}\) given by (21a), \(f_{I_2}=0\), and \(f_{I_4}\) and \(f_{I_5}\) given by (25) and (26), respectively, with \(\lambda _z=1.1\), three swelling values \(v=1\), \(v=1.1\), and \(v=1.2\) and different values of \(\beta \). Fiber winding angle values \(\beta =\pi /6\) and \(\beta =\pi /4\) are used in Fig. 7a and b, respectively. Under these conditions there is only one solution. This is not the same for \(\beta =\pi /3\). Under these circumstances, there are two positive roots fulfilling the bulging condition(20), which are shown in Fig. 7c and d.

Fig. 7
figure 7

This figure shows combinations for \(k_1\ge 0\) and \(k_2\ge 0\) that fulfill the bulging condition (20) for \(\lambda =2\) with \(f_{I_1}\) given by (21a), \(f_{I_2}=0\), and \(f_{I_4}\) and \(f_{I_5}\) given by (25) and (26), respectively, where \(k_1\) and \(k_2\) are normalized by \(\mu _1\). These figures take \(\lambda _z=1.1\) and swelling ratios \(v=1\), \(v=1.1\), and \(v=1.2\) for the matrix. a is for \(\beta =\pi /6\) and b is or \(\beta =\pi /4\). There is a unique solution under these circumstances. On the other hand, c and d show solutions for \(\beta =\pi /3\)

Full size image

4 Axial propagation of bulging under swelling

Axial, quasi-static propagation of a bulging instability in thick-walled cylinders may be subdivided into two stages. In the first stage, the inner pressure remains constant during the ensuing propagation of the bulging instability mode beyond the onset of bifurcation until a suitable configuration is obtained. In the second stage, in a subsequent motion, for an ongoing axial propagation of the bulging the inflation pressure has to be increased. These configurations can be described by two uniform cylinders, one that is characterized by the pair of azimuthal and axial stretches \((\lambda _{\theta 1},\lambda _{z1})\) and other that is characterized by the pair \((\lambda _{\theta 2},\lambda _{z2})\), and by a transition zone between these cylinders (see Fig. 8). Furthermore, it shall be noticed that radial expansion of bulging is related to a decrease of pressure beyond the onset of bulging. These two pairs \((\lambda _{\theta 1},\lambda _{z 1})\) \((\lambda _{\theta 2},\lambda _{z 2})\) fulfill the conditions [12]

$$\begin{aligned}&\dfrac{\partial \hat{W}}{\partial \lambda _z} \Big |_{\lambda _z=\lambda _{z2}}^{\lambda _\theta =\lambda _{\theta 2}} - \dfrac{\lambda _{\theta 2}^{2}-\lambda _{\theta 1}^{2}}{2\lambda _{\theta 1}\lambda _{z1}} \,\,\dfrac{\partial \hat{W}}{\partial \lambda _\theta }\Big |_{\lambda _z=\lambda _{z1}}^{\lambda _\theta =\lambda _{\theta 1}}-\dfrac{\partial \hat{W}}{\partial \lambda _z}\Big |_{\lambda _z=\lambda _{z1}}^{\lambda _\theta =\lambda _{\theta 1}}=0, \end{aligned}$$
(31a)
$$\begin{aligned}&\hat{W}\Big |_{\lambda _z=\lambda _{z2}}^{\lambda _\theta =\lambda _{\theta 2}} -\lambda _{z2}\dfrac{\partial \hat{W}}{\partial \lambda _z}\Big |_{\lambda _z=\lambda _{z2}}^{\lambda _\theta =\lambda _{\theta 2}} - \hat{W}\Big |_{\lambda _z=\lambda _{z1}}^{\lambda _\theta =\lambda _{\theta 1}} + \lambda _{z1} \dfrac{\partial \hat{W}}{\partial \lambda _z}\Big |_{\lambda _z=\lambda _{z1}}^{\lambda _\theta =\lambda _{\theta 1}}=0, \end{aligned}$$
(31b)
$$\begin{aligned}&\lambda _{\theta 2}\lambda _{z2}\dfrac{\partial \hat{W}}{\partial \lambda _\theta } \Big |_{\lambda _z=\lambda _{z1}}^{\lambda _\theta =\lambda _{\theta 1}} -\lambda _{\theta 1}\lambda _{z1}\dfrac{\partial \hat{W}}{\partial \lambda _\theta } \Big |_{\lambda _z=\lambda _{z2}}^{\lambda _\theta =\lambda _{\theta 2}}=0, \end{aligned}$$
(31c)

that can be derived after some manipulations from either [2] or [11]. While the works [2, 11] do not include swelling, the system of Eq. (31) has been studied for an swellable fiber-reinforced matrix in [12], where the material has been taken to remain incompressible at all instants of swelling.

Fig. 8
figure 8

Cylindrical membrane under a constant inner pressure: bulge in equilibrium after bifurcation. Two uniform cylinders are shown, one that is characterized by the pair of azimuthal and axial stretches \((\lambda _{\theta 1},\lambda _{z 1})\) and other that is characterized by the pair \((\lambda _{\theta 2},\lambda _{z 2})\). These cylinder are connected by a transition zone between these cylinders

Full size image

In what follows, we study axial propagation of bulging by means of some numerical examples. Results are presented in the different parts of Table 1.

Table 1 Panel (a): Axial propagation of bulging under swelling for the Mooney–Rivlin material (28) The table shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\) and different values for the swelling v that are determined from system of Eq. (31). The results are shown for \(\mu _2=0\) (neo-Hookean specialization) and for \(\mu _2=-\mu _1/2\). Panel (b) studies axial propagation of bulging under swelling for a neo-Hookean matrix and fiber behavior in terms of \(I_4\) (29). It shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\), for the material parameters \(k_1=0.05\mu _1\) and \({\bar{k}}_1=0.1\), and for different values of the swelling v that are determined from system of Eq. (31). The results are consistent with those that have been obtained by Demirkoparan & Merodio [12] and Alhayani et al. [11]. Panel (c) studies axial propagation of bulging under swelling for a neo-Hookean matrix and fiber behavior in terms of \(I_5\) (30). It shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\), for the material parameters \(k_2=0.01\mu _1\) and \({\bar{k}}_2=0.01\), and for different values of the swelling v that are determined from system of Eq. (31)
Full size table

Consider a swellable Mooney–Rivlin material with a strain energy density given by Eq. (28). Table 1a shows results for axial propagation of bulging under swelling. This table shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\) and different values for the swelling v that are determined from system of Eq. (31). The results are shown for \(\mu _2=0\) (neo-Hookean specialization) and for \(\mu _2=-\mu _1/2\).

We now consider a neo-Hookean matrix with embedded fibers and \(k_1\ne 0\), \(k_2=\mu _2=0\), for which the strain energy density function (18) takes the form (29). Table 1b shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\), and for the material parameters \(k_1=0.05\mu _1\) and \({\bar{k}}_1=0.1\), and for different values of the swelling v that are determined from system of Eq. (31).

In a last example, let us now turn to study a neo-Hookean matrix with embedded fibers and \(k_2\ne 0\), \(k_1=\mu _2=0\), for which the strain energy density function (18) takes the form (30). Table 1c shows the pairs \((\lambda _{\theta 1},\lambda _{z 1})\) and \((\lambda _{\theta 2},\lambda _{z 2})\) for \(\lambda _{z 1}=1.1\), and for the material parameters \(k_2=0.01\mu _1\) and \({\bar{k}}_2=0.01\), and for different values of the swelling v that are determined from system of Eq. (31).

5 Concluding remarks

This article studies the effect of swelling, as well as other factors, on bulging initiation and propagation of pressurized membranes under axial loading. The material is taken to be hyperelastic and reinforced by two fiber families in helical arrangements.

In our present modeling, the fibers are taken to be parallel to each other within a fiber family. The article [12] studies bulging bifurcation for a fiber-reinforced cylinder, in which the fiber dispersion is quantified by the so-called \(\kappa \)-model [37]. Modeling of fiber dispersion has been beyond the scope of this article, and for different modeling techniques that account for fiber dispersion we refer to the review [38]. The mechanical properties of the constituents have been taken to be time-independent, and for example viscoelastic behavior has been neglected. In biological soft tissue, collagen fiber undergoes remodeling processes that are related to different biochemical and physical effects such as their stretch history [39]. The interplay of inflation behavior and their instabilities in combination with fiber remodeling has been studied in [40]. In this work, the hollow cylinder has been studied in the membrane treatment. The recent article studies the inflation of thick-walled cylinders with mechano-sensitive fiber remodeling [41].

The herein presented results may have implications in the initiation and propagation of aneurysms in the different layers of arterial tissue.