Bulging initiation and propagation in fiber-reinforced swellable Mooney–Rivlin membranes

This article considers a thin-walled hollow cylinder, which is composed of a fibrous and swellable hyperelastic material. The fibers are arranged in two families and they are taken to be parallel within each fiber family. The two fiber families are also assumed to be mechanically equivalent and symmetrically disposed in the ground substance material. At each instant of the homogeneous swelling, the material is taken to be incompressible. This article studies the interplay of swelling, fiber orientation, and the mechanical properties of the constituents on the initiation as well as on the axial propagation of bulging.

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.

Geometric and mechanical conditions
Consider an incompressible elastic body in its undeformed reference configuration B r , in which the location of a material point is described by a position vector X. In the deformed configuration B, the location of the same material point is described by the position vector x, so that the mapping for this deformation is described by the deformation gradient tensor F = ∂x/∂X. The corresponding left and right Cauchy-Green deformation tensors are B = FF T and C = F T F. Now consider a cylindrical coordinate system that is defined by the three base unit vectors {E R , E , E Z } into the radial, circumferential, and axial directions, respectively. In the undeformed reference configuration B r , the location of a hollow cylinder can be described in terms of the cylindrical coordinates {R, , Z } as 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 X in the undeformed configuration then reads Prior to its bifurcation, the cylinder is subjected to an axial force N , an inflation pressure P, and a swelling field 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 where (r, θ = , z) are the cylindrical coordinates, {e r , e θ , e z } are the three base unit vectors of the deformed configuration, and is the length of the cylinder in the deformed configuration. The coordinates of the deformed configuration are where λ 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, λ θ = λ = r/R > 0, and the radial stretch then becomes Furthermore, in the cylindrical coordinate system the deformation gradient tensor and the left and right Cauchy-Green deformation tensors become respectively.

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 X of the cylinder contains both matrix and fiber material. In an incompressible and swellable solid, the Cauchy stress tensor can be expressed as where F is given in (6a), 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 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 C (and due to the symmetry, also of B) in Eq. (6b) have the forms 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). 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 M 1 = cos βe z + sin βe θ , M 2 = cos βe z − sin βe θ (10) describe the orientation of the two fiber families in the reference configuration using the angle β (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 M 1 and M 2 are where the invariants I 4 and I 5 are associated with the first fiber family (the one with orientation angle M 1 ) and the invariants I 6 and I 7 are associated with the second fiber family (the one with orientation angle 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 In order to simplify the notation, partial derivatives ofW shall be indicated by the subscripts

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 div σ = 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) 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] The radial normal stresses are negligible in the membrane approximation, hence we take σ rr = 0. Then the circumferential and axial components σ θθ and σ zz of the Cauchy stress tensor (7) become 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 in which δ is the increment. The bulging condition for a cylinder membrane of radius R and length L is derived in [7] and given by where The last term on the left-hand side of (16) vanishes for a cylinder of infinite length, L → ∞, which is the focus of this analysis.

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 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 μ 1 and μ 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) The material parameters k 1 and k 2 can be regarded as relative strength parameters for the fiber stiffness. Moreover k 1 andk 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 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 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 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 ), where I 5 is given in terms of the stretches in Eq. (11b).
Bulging in the limitk 1 → 0 andk 2 → 0: The exponential terms in the strain energy density function reduce to the standard forms in the limits In these limits, the contributions of the fibers to the bulging condition in (22) and (23) become and 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 In both figures, the blue curves show the results for μ 2 /μ 1 > 0, the red curves show the response for μ 2 /μ 1 < 0, and the black curve shows the neo-Hookean response (μ 2 = 0). b depicts relations between the circumferential and the axial stretches λ and λ z , respectively, which fulfill the bulging condition f I1 + f I2 = 0 for different ratios of μ 1 to μ 2 Different forms of (27) (for instance with μ 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.

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, This model has been applied in the modeling of the elastinous ground substance behavior of arterial tissue in works such as [34]. While μ 1 is usually taken to be positive, the parameter μ 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 λ and λ 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 μ 2 /μ 1 > 0, the red curves show the response for μ 2 /μ 1 < 0, and the black curve shows the neo-Hookean response (μ 2 = 0). Figure 2a shows different curves f /μ 2 1 vs. λ for an axial stretch of λ z = 1.1 and different ratios of μ 1 to μ 2 . These curves take larger values for f when μ 2 increases. The curve for μ 2 = −μ 1 does not take positive values for f /μ 2 1 . Figure 2b depicts relations between the circumferential and the axial stretches λ and λ z , which fulfill the bulging condition f I 1 + f I 2 = 0 for different ratios of μ 1 to μ 2 . Notice that for μ 2 = −μ 1 and for the depicted values for λ z and λ the bulging condition f = 0 is not fulfilled.

Fibers in neo-Hookean Matrix
In this section, we consider two fiber families that are embedded in a swellable neo-Hookean ground material (μ 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 .  (18)) Let us consider a neo-Hookean matrix (μ 2 = 0) with embedded fibers and k 1 = 0, k 2 = 0, for which the strain energy density (18) becomes

Fiber behavior in terms of I
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).
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 μ 2 1 ), for an axial stretch λ z = 1.1, β = π/4, and k 1 = μ 1 /100. Increasing values fork 1 or v lead to larger values of f . The solid lines represent results fork 1 → 0, which corresponds to the mechanical response of (27) for μ 2 = k 2 = 0. The dashed lines represent the results fork 1 = 0.4. Whenk 1 = 0.4 bulging may occur for v = 1, while ifk 1 = 0.4 and v = 1.2 bulging will not occur because f remains positive. Figure 3b shows different combinations for λ z and λ that fulfill the condition f I 1 + f I 4 = 0 for β = π/4 and k 1 = μ 1 /100. The different line styles represent the results for different values ofk 1 . This diagram shows that as k 1 → 0 the curves are decreasing so that for an axial stretch λ z we find a unique value of the azimuthal stretch λ on the depicted range. On the other hand, whenk 1 > 0 a value of the axial stretch may have multiple corresponding values of the azimuthal stretch λ associated with bulging.
The examples of Fig. 3 are restricted to a fixed value β = π/4 and the impact of the fiber orientation on the bulging for a fiber-reinforced neo-Hookean material (29) is studied in Fig. 4. Figure 4a shows the function (20) for f = f I 1 + f I 4 (normalized by μ 2 1 ) for an axial stretch λ z = 1.1, The different curves show that if β is sufficiently close to π/2, then bulging will not occur.
The polar diagram of Fig. 4b shows different combinations for β (circumferential orientation) and λ (radial values) that fulfill the condition f I 1 + f I 4 = 0 for k 1 = μ 1 /100,k 1 = 0.4, and λ z = 1 and the three swelling ratios v = 1, v = 1.1, and v = 1.2. This diagram shows that the smallest values of λ that fulfill the bulging condition f I 1 + f I 4 = 0 occur for β = 0 (β = ±π ), while the largest values of λ that fulfill this condition occur for β = ±π/2. The figure shows symmetries with respect to horizontal and vertical axes of the polar diagram.

Fiber behavior in terms of I
Let us now consider a neo-Hookean matrix with embedded fibers and k 1 = 0, k 2 = 0, for which the strain energy density function (18) becomes 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). 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 μ 2 1 ) for an axial stretch of λ z = 1.1, β = π/4, and k 1 = μ 1 /100. The solid lines represent results fork 2 → 0, and the dashed lines represent the results fork 2 = 0.001. Figure 5b shows different combinations for λ z and λ that fulfill the condition f I 1 + f I 5 = 0 for β = π/4, k 2 = μ 1 /100,k 2 → 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). Figure 6a shows the function (20) for f = f I 1 + f I 5 (normalized by the square of the shear modulus μ 1 ) for an axial stretch of λ z = 1.1, k 2 = μ 1 /100, v = 1,k 2 = 0, and for different values β.
The polar diagram of Fig. 6b shows different combinations for β (circumferential orientation) and λ (radial values) that fulfill the condition f I 1 + f I 5 = 0 for λ z = 1.00, k 2 = μ 1 /100,k 2 = 0.0, and two swelling ratios v = 1 and v = 1.2. This diagram shows the following for values π/2 ≥ β ≥ 0. As β increases from β = 0, the   values of λ associated with bulging increase until β reaches an angle close to π/6 from which as β increases the values of λ associated with bulging decrease.

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 μ 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 ≥ 0 and k 2 ≥ 0 that fulfill the bulging condition (20) for λ = 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 λ z = 1.1, three swelling values v = 1, v = 1.1, and v = 1.2 and different values of β. Fiber winding angle values β = π/6 and β = π/4 are used in Fig. 7a and b, respectively. Under these conditions there is only one solution. This is not the same for β = π/3. Under these circumstances, there are two positive roots fulfilling the bulging condition (20), which are shown in Fig. 7c   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.
In a last example, let us now turn to study a neo-Hookean matrix with embedded fibers and k 2 = 0, k 1 = μ 2 = 0, for which the strain energy density function (18) takes the form (30). Table 1c shows the pairs (λ θ1 , λ z1 ) and (λ θ2 , λ z2 ) for λ z1 = 1.1, and for the material parameters k 2 = 0.01μ 1 andk 2 = 0.01, and for different values of the swelling v that are determined from system of Eq. (31).

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 κ-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]. Table 1 Panel (a): Axial propagation of bulging under swelling for the Mooney-Rivlin material (28) The table shows the pairs (λ θ 1 , λ z1 ) and (λ θ 2 , λ z2 ) for λ z1 = 1.1 and different values for the swelling v that are determined from system of Eq. (31). The results are shown for μ 2 = 0 (neo-Hookean specialization) and for μ 2 = −μ 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 (λ θ 1 , λ z1 ) and (λ θ 2 , λ z2 ) for λ z1 = 1.1, for the material parameters k 1 = 0.05μ 1 andk 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 (λ θ 1 , λ z1 ) and (λ θ 2 , λ z2 ) for λ z1 = 1.1, for the material parameters k 2 = 0.01μ 1 andk 2 = 0.01, and for different values of the swelling v that are determined from system of Eq. ( The herein presented results may have implications in the initiation and propagation of aneurysms in the different layers of arterial tissue.