# Inflation-induced aneurysm formation and evolution in graded cylindrical
tubes of arbitrary thickness

Yang Liu  tracy<EMAIL_ADDRESS>Liu Yang  Yu-Xin Xie<EMAIL_ADDRESS>Department
of Mechanics, School of Mechanical Engineering, Tianjin University, Tianjin
300354, China  Tianjin Key Laboratory of Modern Engineering Mechanics, Tianjin
300354, China

###### Abstract

We study the initiation and evolution of aneurysmal morphology in a
pressurized soft tube where the elastic modulus is non-uniform in the radial
direction. The primary deformation prior to instability is characterized
within the framework of nonlinear elasticity for a general material
constitution and a generic modulus gradient. To unravel the influence of
modulus gradient on aneurysm formation, we employ the incompressible Gent
model and select three representative modulus gradients, including a linear,
an exponential, and a sinusoidal function. In particular, the sinusoidal
distribution can be used to model actual artery structure. In addition, two
prototypical loading conditions are considered, namely, either the resultant
axial force or the axial length can be fixed. Based on an explicit bifurcation
condition in terms of the internal pressure and the resultant axial force for
aneurysm formation or localized bulging, an exhaustive theoretical analysis on
bulge initiation is carried out and the effect of geometric and material
parameters and modulus gradient on the critical stretch generating localized
bulging is revealed. It turns out that the modulus mismatch, as well as the
position of maximum modulus, can dramatically affect the onset of localized
bulging. Then we analytically elucidate the influence of modulus gradient on
bulge propagation according to Maxwell’s equal-area rule and conduct a finite
element analysis of bulge evolution. To perform post-bifurcation analysis, a
robust finite element model for graded Gent material is established in Abaqus
by UHYPER subroutine coding. Interestingly, it is found that a sinusoidally
distributed modulus has negligible influence on the critical stretch of bulge
initiation, the deformation process of bugle growth, and the maximum size of a
bulge. The current analysis not only can be used to qualitatively explain why
a healthy human artery evolves into a sandwich structure where the
intermediate layer is stiffest but also can provide useful insight into
localized instabilities such as necking and beading in graded structures.

###### keywords:

Graded tubes , Aneurysm formation , Bulge evolution , Nonlinear elasticity ,
Parametric study , Finite element analysis

##  1 Introduction

Inflation of a cylindrical rubber tube usually generates a prototypical
localization instability, known as localized bulging, and such an interesting
phenomenon was first recorded in detail and primarily analyzed by Mallock,
(1891). With applications to bulge initiation and propagation in pressure
vessels (Chater and Hutchinson,, 1984; Kamalarasa and Calladine,, 1988;
Kyriakides,, 2007) and aneurysm formation in human arteries (Rodríguez and
Merodio,, 2011; Sakalihasan et al.,, 2018), inflation of a cylindrical
hyperelastic tube offers a pertinent paradigm to study the mechanism of
localized instabilities (Kyriakides and Chang,, 1990, 1991; Hejazi et al.,,
2021). In addition, recent advances on inflated tubes and localized bulging
can shed light on bifurcation in arteries affected by Marfan’s syndrome
(Haughton and Merodio,, 2009; Merodio and Haughton,, 2010), bulge prevention
in energy harvesting devices (Bucchi and Hearn,, 2013), inflation-driven soft
robots (Jin et al.,, 2021), inflation of nematic elastomers (He et al.,,
2020), and bulge formation subject to additional effects of electric actuation
(Lu et al.,, 2015), swelling (Demirkoparan and Merodio,, 2017), magnetic field
(Reddy and Saxena,, 2018), and plasticity (Takla,, 2018).

In general, the evolution of localized bulging is mainly composed of three
typical periods, including bulge initiation, growth, and propagation (Wang et
al.,, 2019). Specifically, experimental evidence suggests that the diameter of
the bulge will reach a maximum in the propagation stage, and a two-phase
deformation can be observed (Kyriakides and Chang,, 1990). On the theoretical
side, the pioneering works on bifurcation analysis of inflated membrane tubes
and thick tubes were conducted by Haughton and Ogden, 1979 (a, b). Later, Fu
et al., (2008) found that a zero mode, which is usually viewed as another
uniform deformation, precisely generates localized bulging, and a bifurcation
condition was derived for a membrane tube based on the dynamical systems
theory. Since then, multiple fundamental problems were investigated, such as
solution stability (Pearce and Fu,, 2010; Fu and Xie,, 2010; Il’ichev and Fu,,
2012), influence of material constitution (Pearce,, 2012), imperfection
sensitivity (Fu and Xie,, 2012), and dynamic inflation (Rodríguez-Martínez et
al.,, 2015). Within the framework of nonlinear elasticity, an explicit
bifurcation condition for thick tubes was derived by Fu et al., (2016), and
the accuracy of membrane assumption was examined. This bifurcation condition
paves a convenient way to study localized bugling in rotating cylinders (Wang
et al.,, 2017), fiber-reinforced tubes (Wang and Fu,, 2018), and layered tubes
(Liu et al.,, 2019; Ye et al.,, 2019). Furthermore, a comprehensive
understanding of localized bulging can provide useful insight into localized
necking (Fu et al.,, 2018). We mention that the evolution of a bulge can only
be handled in the nonlinear regime. To perform post-bifurcation analysis for
inflated thick-walled cylinders, Alhayani et al., (2013) presented a numerical
procedure by virtue of nonlinear finite element models. Within the framework
of nonlinear elasticity, Ye et al., (2020) carried out a weakly nonlinear
analysis to deduce the amplitude equation for a bulge and discovered an
interesting parametric domain where necking occurs instead of bulging.
Recently, the analysis of localizations was further extended to situations
where surface tension is involved Wang, (2020); Emery and Fu, 2021a ; Emery
and Fu, 2021b ; Fu et al., (2021). Recent papers have shown a growing interest
in the community for such a paradigmatic localization problem.

In addition to the investigations mentioned above, some effort was dedicated
to establishing a reduced model applicable to inspect bulge formation and
propagation in cylindrical hyperelastic tubes. Lestringant and Audoly, (2018)
derived a one-dimensional model for the analysis of localized bulging in long
tubes making use of nonlinear membrane theory and asymptotic expansion. A
fundamental assumption is that localization can be viewed as a long-wavelength
periodic mode. Subsequently, Giudici and Biggins, (2020) proposed an energy
method in virtue of asymptotic expansions to track the nonlinear evolution of
an inhomogeneous solution and further to identify phase separation of typical
localization such as beading, bulging, and necking.

As mentioned earlier, localized bulging in inflated tubes can be viewed as an
analogue to aneurysm formation in arteries (Seshaiyer and Humphrey,, 2001; Fu
et al.,, 2012). Accordingly, theoretical or numerical analysis was performed
to supply mechanical insight into possible pathogenesis behind aneurysms, see
Alhayani et al., (2014); Fu and II’ichev, (2015); Varatharajan and DasGupta,
(2017, 2018) and the references therein. In particular, our previous analysis
reveals a fact that an aneurysm is more prone to appear as the innermost layer
of an artery stiffens (Ye et al.,, 2019). This qualitatively interprets why
the risk of aneurysm formation for a person becomes higher with aging (Kohn et
al.,, 2015).

Generally speaking, two loading conditions, i.e. either the resultant axial
force or the axial length is specified, are adopted in experiments. The former
was adopted in Kyriakides and Chang, (1990, 1991); Wang et al., (2019); Hejazi
et al., (2021) while the latter was used in Pamplona et al., (2006); Goncalves
et al., (2008); Guo et al., (2016); Wang et al., (2019); Hejazi et al.,
(2021). On the one hand, the resultant axial force can be fixed by suspending
a dead weight at one end while air comes into the tube from the other end. In
this scenario, the curve of pressure versus volume ratio possesses an
$N$-shape. Therefore, localized bulging is induced by a limit-point
instability (Alexander,, 1971; Kanner and Horgan,, 2007; Wang et al.,, 2018)
and the propagation pressure can be predicted by Maxwell’s equal-area rule
(Chater and Hutchinson,, 1984). On the other hand, the fixed axial length can
be acquired by imposing a pre-stretch on the tube and then fixing the total
length. In this situation, the pressure-volume ratio curve may be monotonic.
The effect of pre-stretch on pressurized thin tubes was examined by Mao et
al., (2014); Horný and Netušil, (2016). Furthermore, human arteries usually
suffer a fixed pre-stretch in vivo (Horný et al.,, 2014). In this study, we
also employ these two loading approaches.

It is pointed out that most existing literature on localized bulging is mainly
concerned with a homogeneous or piecewise homogeneous material. Yet
functionally graded materials (abbreviated by FGMs) with continuously varied
physical or mechanical properties have unfolded many unusual functions such as
high-temperature resistance and diminishing stress or strain concentration.
This new type of material first emerged in the middle of the 1980s by mixing
metal and ceramic phases in a manageable way. In doing so, it will acquire a
desirable material property that distributes spatially. Afterwards, the
mechanics of FGM has rapidly taken a center stage in solid mechanics (Ke et
al.,, 2008; Zhong et al.,, 2012; Jha et al.,, 2013). Especially, the
popularity of 3D printing techniques substantially facilitates the fabrication
of FGMs. Besides, human arteries chiefly consist of intima, media, and
adventitia (Holzapfel et al.,, 2000; Gasser et al.,, 2006), which can be
viewed as a sandwich structure and further a special graded structure. In most
practical applications, we intend to prevent localized bulging or aneurysm
formation, such as an Anaconda wave-energy extraction device (Bucchi and
Hearn,, 2013). Naturally, improving the stiffness of a homogeneous tube will
raise the critical pressure inducing localized bulging since pressure is
proportional to the elastic modulus. However, an extremely stiff structure
will lose the ability to deform subjected to a similar load and hence will
degrade certain structural or physiological performance, particularly for
arteries. Naturally, a compromised way is to partially stiffen. Then a well-
imposed question arises. Where is the optimal position to be strengthened, the
innermost surface, the outermost surface, or the middle part? In fact, the
question has been perfectly solved by nature. It is the intermediate layer
that is hardest, and this is also the optimization result of artery evolution
(Holzapfel et al.,, 2000; Chirinos,, 2012). But why? Previous studies have
shown the potential effect of modulus gradient form on the stress distribution
in functionally graded rubberlike cylinders and spheres (Batra,, 2011) or
pattern transition in growing graded tubular tissues (Liu et al.,, 2020). It
is therefore of fundamental significance to clarify the influence of material
inhomogeneity, different modulus gradients as well as the maximum modulus
mismatch on bulge initiation, growth, and propagation. This motivates the
current study. Although deformation or bifurcation analysis was carried by
Batra and Bahrami, (2009) for inflated and everted circular cylinders
comprised of a graded Mooney-Rivlin material and by Chen et al., (2017) for
pressurized graded cylindrical tubes, respectively, localization instability
was still not addressed. In this study, we aim at furnishing a thorough
analysis of localized bulging in graded hyperelastic tubes under the combined
action of internal pressure and axial stretching and unveiling the mechanism
behind structure optimization of arteries.

This paper is organized as follows. In Section 2, we characterize the primary
deformation generated by inflation for a general material constitution and a
generic modulus gradient and establish all required formulas used in the
bifurcation condition within the framework of nonlinear elasticity. Section 3
demonstrates the numerical strategy for identifying bulge initiation and
presents a detailed parametric study for the onset of localized bulging. A
finite element model is constructed and validated in Section 4. A nonlinear
analysis based on finite element analysis is performed in Section 5 and the
bulge propagation is investigated by Maxwell’s equal-area rule. Finally, some
conclusions are given in Section 6.

##  2 Primary deformation and bifurcation condition

Figure 1: (Color online) A graded tube with inner radius A, outer radius B and
length $L$ in its stress-free state.

As shown in Figure 1, a cylindrical tube with inner radius $A$ and outer
radius $B$ in its initial state is considered in this study. Since we focus on
a graded soft tube, it is assumed that the elastic modulus of the tube is no
longer uniformly distributed in the radial direction but is homogeneous in the
axial direction. When such a graded tube is subjected to an internal pressure
$P$ and an axial extension, a primary deformation is first generated, and the
inner and outer radii become $a$ and $b$, respectively. Then we define the
initial state as the reference configuration while the primary state the
current configuration. Assuming that the tube is composed of incompressible
hyperelastic material, we can use $W(\mathbf{F})$ to represent the strain-
energy function where $\mathbf{F}$ is the deformation gradient and is
subjected to a constraint $\operatorname{det}\mathbf{F}=1$, which is the so-
called incompressibility condition.

For the current problem, it is convenient to adopt cylindrical polar
coordinates in both the reference and current configurations. By determining
the centroid as the origin, the position coordinates of the same material
point in the reference and current states are denoted by $(R,\Theta,Z)$ and
$(r,\theta,z)$, respectively, and the tube occupies $A\leqslant R\leqslant B$
and $-L/2\leqslant Z\leqslant L/2$ in the reference state. Bearing in mind
that the primary deformation is axisymmetric, we can write the deformation
gradient as

$\mathbf{F}=\lambda_{1}\bm{e}_{r}\otimes\bm{e}_{r}+\lambda_{2}\bm{e}_{\theta}\otimes\bm{e}_{\theta}+\lambda_{3}\bm{e}_{z}\otimes\bm{e}_{z},$
(2.1)

where $\lambda_{i}$ $(i=1,2,3)$ is the principal stretch in $i$-direction and
the subscripts 1, 2 and 3 denote the $r$-, $\theta$\- and $z$-directions.
Furthermore, the common orthonormal basis
$\\{\bm{e}_{r},\bm{e}_{\theta},\bm{e}_{z}\\}$ is used in both states. In
particular, the first two principal stretches are equal to

$\lambda_{1}=\dfrac{\textrm{d}r}{\textrm{d}R},~{}~{}\lambda_{2}=\dfrac{r}{R},$
(2.2)

and the third one $\lambda_{3}$ is a constant since the primary deformation is
homogeneous in the $z$-direction.

The strain-energy function can be written as a function of these principal
stretches $W(\lambda_{1},\lambda_{2},\lambda_{3})$. Accordingly, the non-zero
components of Cauchy Stress tensor are given by

$\displaystyle\sigma_{ii}=\lambda_{i}W_{,i}-p,~{}~{}\textrm{no summation},$
(2.3)

where $p$ depicts the Lagrange multiplier enforcing the incompressibility
condition. Meanwhile, we have defined that a comma means derivative with
respect to the corresponding variable, for instance, $W_{,2}=\partial
W/\partial\lambda_{2}$.

It can be readily checked that the primary deformation is governed by the sole
equilibrium equation:

$\displaystyle\dfrac{\textrm{d}\sigma_{rr}}{\textrm{d}r}+\dfrac{\sigma_{rr}-\sigma_{\theta\theta}}{r}=0.$
(2.4)

In order to facilitate subsequent analysis, we define two parameters according
to the incompressibility condition

$\displaystyle\lambda=\lambda_{2}=\frac{r}{R},~{}\lambda_{z}=\lambda_{3}.$
(2.5)

Correspondingly, the first principal stretch is determined by
$\lambda_{1}=1/(\lambda\lambda_{z})$. In addition, the following relations can
be derived :

$\displaystyle
r^{2}=\lambda_{z}^{-1}(R^{2}-A^{2})+a^{2},~{}~{}\theta=\Theta,~{}~{}z=\lambda_{z}Z.$
(2.6)

Next, we replace $\lambda_{1}$ in $W$ and employ a reduced strain-energy
function $w$ which is determined by

$\displaystyle w(\lambda,\lambda_{z})\equiv
W(1/(\lambda\lambda_{z}),\lambda,\lambda_{z}).$ (2.7)

Subsequently, we find that

$\displaystyle
w_{,1}=W_{,2}-\frac{1}{\lambda^{2}\lambda_{z}}W_{,1},~{}~{}w_{,2}=W_{,3}-\frac{1}{\lambda\lambda_{z}^{2}}W_{,1}.$
(2.8)

Substituting the stress components in (2.3) into the governing equation (2.4)
and utilizing (2.8) yield

$\displaystyle\sigma_{rr}=\int_{r_{0}}^{r}\frac{\lambda
w_{,1}}{r}\textrm{d}r,$ (2.9)

where $r_{0}$ is a constant. Actually $w$ is dependent on $R$ for a graded
material. It can be deduced from $(\ref{eq2_5})_{1}$ and $(\ref{eq2_6})_{1}$
that $R=A\sqrt{\lambda_{z}\lambda_{a}^{2}-1}/\sqrt{\lambda_{z}\lambda^{2}-1}$.
Then we apply a variable exchange $\lambda=r/R$ and a new identity of the
Cauchy stress $\sigma_{rr}$ is given by

$\sigma_{rr}=\int_{\lambda_{0}}^{\lambda}\frac{w_{,1}}{1-\lambda^{2}\lambda_{z}}\textrm{d}\lambda,$
(2.10)

where $\lambda_{0}$ is another constant to be determined later and we have
replaced $R$ by
$A\sqrt{\lambda_{z}\lambda_{a}^{2}-1}/\sqrt{\lambda_{z}\lambda^{2}-1}$ such
that the integrand in (2.10) is only a function of $\lambda$ (see also
equation (2.15)).

Before we proceed further, it is appropriate to introduce two hoop stretches
at the inner surface and outer surface as follows

$\lambda_{a}=\frac{a}{A},~{}~{}\lambda_{b}=\frac{b}{B},$ (2.11)

and their relation can be obtained from (2.6)

$\lambda_{b}=\sqrt{\frac{B^{2}-A^{2}+A^{2}\lambda_{a}^{2}\lambda_{z}}{B^{2}\lambda_{z}}}.$
(2.12)

It is assumed that the outer surface of the tube is traction-free and the
inner surface suffers a pressure $P$. Therefore, we obtain

$\sigma_{rr}|_{\lambda=\lambda_{b}}=0,~{}~{}\sigma_{rr}|_{\lambda=\lambda_{a}}=-P.$
(2.13)

It then follows from the boundary condition $(\ref{eq2_13})_{1}$ that the
constant of integration $\lambda_{0}$ is identical to $\lambda_{b}$. In
addition, the other one $(\ref{eq2_13})_{2}$ furnishes

$P=\int_{\lambda_{b}}^{\lambda_{a}}\frac{w_{,1}}{1-\lambda^{2}\lambda_{z}}\textrm{d}\lambda.$
(2.14)

Actually, the pressure $P$ can be regarded as a function of $\lambda_{a}$ and
$\lambda_{z}$. Similarly, the resultant axial force $N$ in any cross-section
is also dependent on $(\lambda_{a},\lambda_{z})$ and holds the following form

$\displaystyle N=2\pi\int_{a}^{b}\sigma_{zz}r\textrm{d}r-P\pi a^{2}=\pi
A^{2}(\lambda_{a}^{2}\lambda_{z}-1)\int_{\lambda_{b}}^{\lambda_{a}}\frac{2\lambda_{z}w_{,2}-\lambda
w_{,1}}{(\lambda^{2}\lambda_{z}-1)^{2}}\lambda\textrm{d}\lambda.$ (2.15)

It should be mentioned that (2.14) and (2.15) are valid for both homogeneous
and graded materials. For graded tubes, the reduced strain-energy function in
these two expressions contains radius-dependent elastic modulus. Moreover,
(2.14) can be seen as a special case in Chen et al., (2017) where both the
inner and outer surfaces suffer pressures.

Currently, the primary deformation induced by inflation has been fully
characterized by the above two equations for a graded hyperelastic tube of
finite thickness. In practice, there are usually two different end conditions,
i.e. either the axial force $N$ or the axial length $\lambda_{z}L$ is fixed.
The former case can be attained by imposing an object of dead weight at one
end and leaving the other end free. Correspondingly, (2.14) and (2.15) provide
two algebraic equations for the two unknowns $\lambda_{a}$ and $\lambda_{z}$
and solving them for a given pressure and a given axial force yields the
solution that characterizes the primary deformation. This kind of loading type
has been used in the experimental investigations carried out by Kyriakides and
Chang, (1990, 1991); Guo et al., (2016); Wang et al., (2019); Hejazi et al.,
(2021). In the other loading condition, the tube is first uniformly stretched
and then the two ends are fastened. As a result, the axial length or the axial
stretch $\lambda_{z}$ is fixed. Furthermore, the axial force $N$ is not
controlled and may increase, decrease or remain constant during the inflation
process. We emphasize that this scenario corresponds to the $in$ $vivo$
situation for human arteries (Horný et al.,, 2014). Also, such an experimental
setup was employed by Pamplona et al., (2006); Goncalves et al., (2008); Wang
et al., (2019); Hejazi et al., (2021). Since $\lambda_{z}$ is prescribed, the
internal pressure $P$ is only a function of $\lambda_{a}$, and the primary
deformation is fully determined by equation (2.14).

In an earlier study by Chen et al., (2017), various bifurcation behaviors were
investigated based on the incremental theory for functionally graded tubes
subjected to internal pressure, external pressure, and end compression, and
periodic modes in the axial and hoop directions were analyzed. In this study,
we consider that the graded tube suffers a combined action of internal
pressure and an axial extension. Therefore, zero-mode may become preferred (Fu
et al.,, 2008), resulting in localized bulging in soft cylindrical tubes. On
the one hand, if either the resultant axial force $N$ or the prescribed axial
stretch $\lambda_{z}$ is greater enough, localized bulging can disappear (Liu
et al.,, 2019). On the other hand, if the axial stretch $\lambda_{z}$ is less
than a threshold value, global buckling can replace localized bulging (Lin et
al.,, 2020; Hejazi et al.,, 2021). In this paper, we intend to investigate
localized bulging (aneurysm formation) in graded tubes and further to
elucidate the influence of modulus gradient on bulge formation, growth, and
propagation. To this end, we merely concentrate on the situation that
localized bulging is preferred as a result of the primary bifurcation.

Bearing in mind that we have denoted the internal pressure $P$ and the
resultant axial force $N$ as functions of the stretches $\lambda_{a}$ and
$\lambda_{z}$, it is suitable to resort to an explicit bifurcation condition
that the Jacobian of $P$ and $F$ in terms of variables $\lambda_{a}$ and
$\lambda_{z}$ vanishes. This concise bifurcation condition was first proposed
by Fu et al., (2016) in the framework of finite elasticity. Later, its
applicability to localized bulging in bilayer tubes was verified numerically
using FE analysis in commercial software Abaqus (Liu et al.,, 2019; Ye et
al.,, 2019). Recently, Yu and Fu, (2021) offered analytical proof of this
bifurcation condition with two relatively weak assumptions. The first one is
that the strain-energy function is isotropic while the second one is that the
deformations before bulge initiation and in the propagation stage are both
axisymmetric. So, it also supports the validity of this bifurcation condition
when applied to graded tubes. Writing explicitly, we obtain the bifurcation
condition for localized bulging as follows

$J(P,N)=\frac{\partial P}{\partial\lambda_{a}}\frac{\partial
N}{\partial\lambda_{z}}-\frac{\partial P}{\partial\lambda_{z}}\frac{\partial
N}{\partial\lambda_{a}}=0.$ (2.16)

With the aid of (2.16), the first bifurcation points for fixed axial force and
fixed axial length can be identified, respectively, by the following equations

$\left\\{\begin{aligned} &\textrm{fixed axial force (loading type
I):}~{}~{}~{}J(P,N)=0,~{}~{}N(\lambda_{a},\lambda_{z})=N_{0},\\\
&\textrm{fixed axial length (loading type
II):}~{}~{}~{}J(P,N)=0,~{}~{}\lambda_{z}=\lambda_{z0},\end{aligned}\right.$
(2.17)

where $N_{0}$ and $\lambda_{z0}$ are given constants. For convenience, we
abbreviate loading type I (II) as LTI (LTII) from now on. Indeed, equation
$(\ref{eq2_17})_{1}$ furnishes two critical stretch values $\lambda_{a}^{c}$
and $\lambda_{z}^{c}$ for LTI where localized bulging initiates, and equation
$(\ref{eq2_17})_{2}$ determines one critical stretch value $\lambda_{a}^{c}$
for LTII. Subsequently, the critical pressures $P_{cr}$ for both loading types
can be calculated from equation (2.14).

Next, we suppose that the tube is composed of the incompressible Gent material
in all illustrative examples, and the strain-energy function is given by

$W=-\frac{\mu(R)}{2}J_{m}\textrm{ln}\bigg{(}1-\frac{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}-3}{J_{m}}\bigg{)},$
(2.18)

where $\mu(R)$ stands for the radius-dependent shear modulus for a graded tube
and $J_{m}$ is a material parameter representing the maximum extensibility. In
particular, we set $J_{m}=97.2$ which is typical for rubber (Gent,, 1996).
Note that this model can well capture the bulge initiation, bulge growth and
bulge propagation in finite element simulations (Liu et al.,, 2019; Ye et
al.,, 2019), and these three distinctive stages were commonly seen in
experiments (Wang et al.,, 2019).

As mentioned earlier, we are also concerned with how different modulus
gradients affect the onset of localized bulging and bulge evolution. We
emphasize that an accurate modulus distribution in a graded structure is
difficult to describe. When fabricated polydimethylsiloxane (PDMS) material is
exposed to UV radiation, Young’s modulus may have an exponential distribution,
resulting in a graded PDMS structure (Chen et al.,, 2018). Generally speaking,
the modulus distribution may be varied in different situations. A power law
was adopted in Batra and Bahrami, (2009) and a linear function was employed in
Chen et al., (2017). In our previous study for growth-induced surface
instabilities, it was found that different modulus distributions can alter
pattern transition (Liu et al.,, 2020). Moreover, it was shown in Ye et al.,
(2019) that a stiffer inner layer makes the structure less stable. If a
sandwich tube has a stiffer or softer core and the two faces share the same
shear modulus, there must exist a competition among these three layers and may
cause some interesting results. This intuitively motivates the investigation
on a non-monotonic modulus distribution, which can further be used to model
human arteries. As a result, we employ three modulus functions as follows

$\left\\{\begin{aligned} &\textrm{linear
distribution:}~{}\mu(R)=\mu_{B}\left[(\beta_{1}-1)\dfrac{R-B}{A-B}+1\right],\\\
&\textrm{exponential
distribution:}~{}\mu(R)=\mu_{B}\left[(\beta_{1}-1)\textrm{exp}\left(\zeta\dfrac{A-R}{B-A}\right)+1\right],\\\
&\textrm{sinusoidal
distribution:}~{}\mu(R)=\mu_{B}\left[(\beta_{2}-1)\sin\left(\dfrac{R-A}{B-A}\pi\right)+1\right],\end{aligned}\right.$
(2.19)

where $\beta_{1}=\mu(A)/\mu_{B}$ signifies the ratio of the shear modulus of
the inner surface to that of the outer surface and $\mu_{B}=\mu(B)$. The
$\beta>1$ means that the inner surface is stiffer than the outer surface. In
the sinusoidal distribution, we define $H=(B+A)/2$ and let
$\beta_{2}=\mu(H)/\mu_{B}$. Moreover, the parameter $\zeta$ in
$(\ref{eq2_19})_{2}$ refers to as the decay rate. In fact, letting $R=B$ in
$(\ref{eq2_19})_{1}$ yields
$\mu(B)=\mu_{B}\left((\beta_{1}-1)\textrm{exp}(-\zeta)+1\right)$. In order to
satisfies $\mu(B)=\mu_{B}$, we must have $\zeta=+\infty$. In our illustrative
examples, we specify $\zeta=30$ as an approximation such that exp$(-30)\approx
9.358\times 10^{-14}$. Meanwhile, the effect of $\zeta$ is out of the scope of
the current study. We then plot these three modulus distributions in Figure 2
to depict the details.

Figure 2: Three typical distributions of the shear modulus $\mu$ in the radial
direction. The parameters are given by $\beta_{1}=\beta_{2}=20$, $A/B=0.8$,
$\zeta=30$.

In this section, we have derived the expressions of the pressure $P$ and the
resultant axial force $N$ for the primary deformation of a graded soft tube.
Meanwhile, a theoretical framework for determining the bulge initiation is
established for a general material model and a general modulus gradient. In
the next section, a detailed parametric analysis will be carried out based on
the bifurcation condition (2.16).

##  3 Theoretical analysis of bulge initiation

We have indicated that this study only focuses on the situation where
localized mode is preferred. Accordingly, the prescribed axial force $N_{0}$
and the prescribed axial stretch $\lambda_{z0}$ in (2.17) should lie in a
proper domain. In the succeeding calculations, we specify $N_{0}=0$ and
$\lambda_{z0}=1.5$ and suppose that the analysis can also shed light on
qualitative effects of different parameters and modulus gradients for other
$N_{0}$ and $\lambda_{z0}$. To facilitate further analysis, we introduce some
dimensionless quantities as follows

$\displaystyle
A^{*}=\dfrac{A}{B},~{}~{}N^{*}=\dfrac{N}{\mu_{B}B^{2}},~{}~{}\mu^{*}=\dfrac{\mu}{\mu_{B}},~{}~{}P^{*}=\dfrac{P}{\mu_{B}},~{}~{}w^{*}=\dfrac{w}{\mu_{B}}.$
(3.1)

Correspondingly, the bifurcation condition for localized bulging is rewritten
as $J(P^{*},N^{*})=0$.

In practice, analytical expressions of the dimensionless pressure $P^{*}$ and
the dimensionless force $N^{*}$ can be derived for many homogeneous material
models such as the neo-Hookean model, Mooney-Rivlin model, Ogden model, and
Gent model. Accordingly, equation (2.17) contains two sets of nonlinear
algebraic equations, and solving them numerically can lead to the critical
stretch $\lambda_{a}^{c}$. However, when the shear modulus becomes position-
dependent in the Gent model, explicit formulas for $P^{*}$ and $N^{*}$ may be
difficult to acquire for a complex modulus distribution, for instance, an
exponential one. Therefore, we need to employ a numerical integration scheme.
To construct a general computation framework for identifying the first
bifurcation point, we will introduce our solution procedure first and then
illustrate the bifurcation results.

###  3.1 Solution strategy

Without loss of generality, we apply the exponential modulus to exhibit the
solution procedure, and all symbolic calculations are carried out in
$Mathematica$ (Mathematica,, 2019). This selection gives rise to that
exponential functions are involved in the integrands of (2.14) and (2.15). So
explicit primitive functions of these two integrations become difficult to
obtain. We shall resort to the built-in command $NIntegrate$ in $Mathematica$.
In the case of $N_{0}=0$, the two unknowns $\lambda_{a}$ and $\lambda_{z}$ can
be identified simultaneously according to $(\ref{eq2_17})_{1}$. In order to
generate effective numerical scheme, we define some necessary quantities as
follows

$\left\\{\begin{aligned}
&P_{1}(\lambda,\lambda_{a},\lambda_{z})=\dfrac{w^{*}_{,1}}{1-\lambda^{2}\lambda_{z}},\\\
&N_{1}(\lambda_{a},\lambda_{z})=\pi(A^{*})^{2}(\lambda^{2}_{a}\lambda_{z}-1),~{}~{}N_{2}(\lambda,\lambda_{a},\lambda_{z})=\lambda\dfrac{2\lambda_{z}w^{*}_{,2}-\lambda
w^{*}_{,1}}{(\lambda^{2}\lambda_{z}-1)^{2}}.\end{aligned}\right.$ (3.2)

By use of the above new notations, we rewrite $P^{*}$ and $N^{*}$ as

$\displaystyle
P^{*}=\int_{\lambda_{b}}^{\lambda_{a}}P_{1}(\lambda,\lambda_{a},\lambda_{z})\textrm{d}\lambda,~{}~{}N^{*}=N_{1}(\lambda_{a},\lambda_{z})\int_{\lambda_{b}}^{\lambda_{a}}N_{2}(\lambda,\lambda_{a},\lambda_{z})\textrm{d}\lambda.$
(3.3)

It is known that the bifurcation condition $J(P^{*},N^{*})=0$ actually
describes a curve in the $(\lambda_{a},\lambda_{z})$-plane. Similarly, the
loading condition $N^{*}=0$ provides another curve in the same plane. So, the
intersection point corresponds to the first bifurcation point where localized
bulging occurs. It is found from equation (2.12) that $\lambda_{b}$ is a
function of $\lambda_{a}$ and $\lambda_{z}$. According to the derivative of
integral with variable bounds, we arrive at

$\left\\{\begin{aligned} &\dfrac{\partial
P^{*}}{\partial\lambda_{a}}=\int_{\lambda_{b}}^{\lambda_{a}}\dfrac{\partial
P_{1}}{\partial\lambda_{a}}\textrm{d}\lambda+P_{1}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{a}}-P_{1}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{b}}\dfrac{\partial\lambda_{b}}{\partial\lambda_{a}},\\\
&\dfrac{\partial
P^{*}}{\partial\lambda_{z}}=\int_{\lambda_{b}}^{\lambda_{a}}\dfrac{\partial
P_{1}}{\partial\lambda_{z}}\textrm{d}\lambda-
P_{1}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{b}}\dfrac{\partial\lambda_{b}}{\partial\lambda_{z}},\\\
&\dfrac{\partial N^{*}}{\partial\lambda_{a}}=\dfrac{\partial
N_{1}}{\partial\lambda_{a}}\int_{\lambda_{b}}^{\lambda_{a}}N_{2}\textrm{d}\lambda+N_{1}\left(\int_{\lambda_{b}}^{\lambda_{a}}\dfrac{\partial
N_{2}}{\partial\lambda_{a}}\textrm{d}\lambda+N_{2}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{a}}-N_{2}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{b}}\dfrac{\partial\lambda_{b}}{\partial\lambda_{a}}\right),\\\
&\dfrac{\partial N^{*}}{\partial\lambda_{z}}=\dfrac{\partial
N_{1}}{\partial\lambda_{z}}\int_{\lambda_{b}}^{\lambda_{a}}N_{2}\textrm{d}\lambda+N_{1}\left(\int_{\lambda_{b}}^{\lambda_{a}}\dfrac{\partial
N_{2}}{\partial\lambda_{z}}\textrm{d}\lambda-
N_{2}(\lambda,\lambda_{a},\lambda_{z})\Big{|}_{\lambda=\lambda_{b}}\dfrac{\partial\lambda_{b}}{\partial\lambda_{z}}\right).\end{aligned}\right.$
(3.4)

Although all integrals in (3.4) can not be solved analytically, we emphasize
that it still formulates the theoretical foundation of our numerical scheme.
If all geometric and material parameters are specified, we could numerically
determine the implicit function of $\lambda_{z}(\lambda_{a})$ according to
$J(P^{*},N^{*})=0$.

Next, we briefly summarize the solution procedure by specifying $A^{*}=0.4$,
$\beta_{1}=30$, $\zeta=30$ and $J_{m}=97.2$. This parametric setting denotes a
thick grade tube where the inner surface is stiffer. We first give an initial
value $\lambda_{z}=1.001$ and seek a solution to $J(P^{*},N^{*})=0$ using
Newton’s iteration method. In doing so, a pair of solution
$\\{\lambda_{z}=1.001,\lambda_{a}=2.026\\}$ can be obtained. We then adopt a
loop algorithm starting from $\lambda_{z}=1.001$ with a small incremental
step. It is well-known that the classical Newton’s method for solving
nonlinear algebraic equations is sensitive to the initial guess. Especially
for complex problems, a proper initial guess can make the solution approach
more efficient. Under the assumption that the curve of $J(P^{*},N^{*})=0$ is
smooth, a nice initial guess can be found around the neighborhood of the
solution of the previous step. By doing so, the relation between $\lambda_{a}$
and $\lambda_{z}$ can be found. We emphasize that this solution procedure can
be applied to another nonlinear algebraic equation $N^{*}=0$, which offers the
second implicit function of $\lambda_{z}$ and $\lambda_{z}$. For simplicity,
extra computation details are omitted. In the case of fixed axial length, the
value $\lambda_{z}=1.5$ supplies a horizontal line in the
$(\lambda_{a},\lambda_{z})$-plane. We supply the graphical interpretation of
the solution strategy for the bulge initiation. Currently, determination of
the first bifurcation point is equivalent to identifying the intersection of
$N^{*}=0$ and $J(P^{*},N^{*})=0$ or the counterpart of $\lambda_{z}=1.5$ and
$J(P^{*},N^{*})=0$, as shown in Figure 3. For the prescribed parameters given
earlier, we find that $\lambda_{z}^{c}=1.182$ and $\lambda_{a}^{c}=1.932$ if
$N^{*}=0$ while $\lambda_{a}^{c}=1.852$ for $\lambda_{z}=1.5$.

So far, we have narrated the solution strategy for identifying the bulge
initiation for graded soft tubes of arbitrary thickness. Especially, it is a
numerical-scheme-based framework. Although we only employ the Gent model, this
framework also supplies a calculation platform for a general material model
and a generic modulus distribution. In the following analysis, we shall apply
the solution procedure to the three different modulus functions in (2.19) and
aim to elucidate the influence of modulus gradient as well as the position
where maximum modulus attains on the bulge initiation.

(a) The resultant axial force is fixed.

(b) The axial length is fixed.

Figure 3: Contour plots of the bifurcation condition and the corresponding
loading curves ($N^{*}=0$ or $\lambda_{z}=1.5$) when the exponential modulus
gradient $(\ref{eq2_19})_{2}$ is used, and the parameters are given by
$A^{*}=0.4$, $\beta_{1}=30$, $\zeta=30$ and $J_{m}=97.2$. The intersections
indicate the first bifurcation point resulting in localized bulging.

###  3.2 LTI: Fixed axial force

First of all, we consider the case of fixed axial force. The aim is to unravel
the effect of different geometrical and material parameters and various
modulus gradients on the bulge initiation. It is emphasized that there are at
most four free parameters in the graded Gent model, videlicet, the
dimensionless inner radius $A^{*}$, the ratio of shear modulus $\beta_{1}$ or
$\beta_{2}$, the $J_{m}$ and the decay rate $\zeta$ in $(\ref{eq2_19})_{2}$.
As introduced before, the last two parameters will be fixed by $J_{m}=97.2$
and $\zeta=30$ in all illustrative examples, and their influence is beyond the
scope of this study. As a result, there are two independent parameters $A^{*}$
and $\beta_{1}$ (or $\beta_{2}$). The former measures the thickness of the
graded tube while the latter displays the maximum modulus variation. In
addition to these parameters, another important factor is the modulus
distribution. Therefore, we shall assign a representative value of $\beta_{1}$
or $\beta_{2}$ and exhibit the results for different inner radii and repeat
this procedure when $A^{*}$ turns into the fixed parameter.

Applying the solution strategy outlined in the preceding subsection, the
dependences of the critical stretch $\lambda_{a}^{c}$ on the dimensionless
inner radius $A^{*}$ are shown in Figure 4. The modulus ratios are set to be
$5$ in Figure 4(a) and $30$ in Figure 4(b), respectively. It can be seen that
the critical stretch $\lambda_{a}^{c}$ always decreases monotonically with
increased $A^{*}$. Note that varying $A^{*}$ is equivalent to altering the
thickness of the tube. This implies the fact that localized bulging is
invariably easier to take place in a thinner tube, regardless of the modulus
distribution. The curves for the sinusoidal function are always higher,
leading to a greater critical stretch for a given $A^{*}$ and further a more
stable structure compared to the other two modulus distributions. We observe
that the curves for the sinusoidal function in both figures are almost
identical. It is therefore concluded that the modulus ratio $\beta_{2}$ has a
weak influence on the critical bifurcation stretch $\lambda_{a}^{c}$.
Furthermore, the dashed line is higher than the solid one in 4(a) while
becomes lower in 4(b). This indicates that a graded tube with the shear
modulus decaying exponentially from the inner surface is more stable than the
counterpart where the shear modulus decays linearly when $\beta_{1}=5$.
However, if $\beta_{1}=30$, a reverse inference can be obtained. Ultimately,
as $A^{*}$ goes to unity, the deviations among all curves almost vanish, so
modulus gradient can be ignored in extremely thin tubes.

(a) $\beta_{1}=\beta_{2}=5$.

(b) $\beta_{1}=\beta_{2}=30$.

Figure 4: The critical stretch $\lambda_{a}^{c}$ versus the dimensionless
inner radius $A^{*}$ when the three modulus functions in (2.19) are adopted.
The resultant axial force is given by $N^{*}=0$, and the other parameters are
given by $\zeta=30$ and $J_{m}=97,2$.

The above analysis primarily reveals the complexity of the influence of
modulus gradient on the bifurcation threshold. In the following study, the
modulus ratios $\beta_{1}$ and $\beta_{2}$ are taken into account to explore
the impact of the modulus ratio as well as the modulus distribution. Figure 5
displays the relations between $\lambda_{a}^{c}$ and the modulus ratio for
three modulus gradients by specifying $A^{*}=0.7$. The right figure highlights
the detail of the left one when the modulus ratio is around unity. Seen from
Figure 5, we find that the critical stretch approaches a limit value for all
modulus functions as $\beta_{1}$ or $\beta_{2}$ increases. If the shear
modulus has a sinusoidal distribution, the bulge initiation is nearly
independent of the modulus ratio $\beta_{2}$, especially when $\beta_{2}>1$.
Yet the critical stretch is always a decreasing function for the other two
modulus distributions. Specifically, if the shear modulus grows exponentially
from the inner surface ($\beta_{1}<1$), the bifurcation threshold only varies
in a very small range compared to the linearly increasing counterpart, and a
linear distribution always delays the occurrence of localized bulging.
Nevertheless, when $\beta_{1}>1$, the bifurcation threshold decreases fast and
starts to be impervious to the variation of $\beta_{1}$ as $\beta_{1}$
approaches practically 30 when the linear function $(\ref{eq2_19})_{1}$ is
selected. On the one hand, it is observed from figure 5 that these three
curves connect at $\beta_{1}=\beta_{2}=1$ where a graded tube is reduced to a
homogeneous one. On the other hand, the solid line and dashed line intersect
at another point and we define the corresponding horizontal coordinate as
$\beta_{c}$. It can be solved that $\beta_{c}=12.78$ for $A^{*}=0.7$.
Accordingly, the exponential distribution retards the bifurcation for
$1<\beta_{1}<\beta_{c}$ while diminishes the bulge initiation for
$\beta>\beta_{c}$. In other words, compared to a linearly decayed shear
modulus from the inner surface, an exponential counterpart gives rise to a
more stable structure as $\beta$ is lower than $\beta_{c}$ while resulting in
a less stable structure otherwise.

Figure 5: The critical stretch $\lambda_{a}^{c}$ versus the ratio of modulus
when the three modulus functions in (2.19) are adopted. (a) The modulus ratio
ranges from 0.01 to 101 and (b) the modulus ratio ranges from 0.01 to 1. The
right figure is a blowup of the right figure when the modulus ratio is small.
The resultant axial force is given by $N^{*}=0$, and the other parameters are
given by $A^{*}=0.7$, $\zeta=30$ and $J_{m}=97,2$.

To verify that the main characteristic in Figure 5 is irrelevant to the
special choice $A^{*}=0.7$, we further present the dependence of
$\lambda^{c}_{a}$ on the modulus ratio for three modulus functions in Figure 6
by considering a thicker tube where $A^{*}=0.5$. Likewise, the right part
exhibits a blowup when the modulus ratio is around unity. It can be seen that
Figure 6 differs from Figure 5 in the range of the vertical axis. In addition
to this small divergence, these two figures are qualitatively the same. So
similar conclusion can be drawn and we identify $\beta_{c}=11.22$ for
$A^{*}=0.5$.

In summary, if other conditions remain identical and only the thickness can be
varied, a higher thickness always delays bulge formation in a pressurized
graded tube. However, the effect of the modulus gradient is quite distinctive.
Although we only select three prototypical modulus functions as shown in
(2.19), the bifurcation threshold versus the modulus experiences various
tendencies. When the shear modulus either grows monotonically or declines
monotonically from the inner surface, the critical stretch $\lambda_{a}^{c}$
is constantly a decreasing function as $\beta_{1}$ increases and will attain a
constant value for large enough $\beta_{1}$. This implies that the bifurcation
threshold remains practically unaffected when the modulus mismatch between the
inner and outer surfaces is fairly high. Furthermore, these two curves
intersect at $\beta_{1}=1$ and $\beta_{1}=\beta_{c}>1$, leading to a
transition zone $1<\beta_{1}<\beta_{c}$ where a graded tube with shear modulus
distributing exponentially is more stable. Yet a non-monotonically sinusoidal
distribution of the shear modulus leads to another situation. The curve of
$\lambda_{a}^{c}$ versus $\beta_{c}$ is slowly increasing if $\beta_{2}<1$ and
becomes approximately a horizontal line once $\beta_{2}$ passes the value of
unity. Specifically, the critical stretch $\lambda_{a}^{c}$ is practically
identical to that of the homogeneous counterpart. As a result, this special
graded tube has similar deformation behavior to its homogeneous analogue. We
mention that a graded tube with the shear modulus distributed sinusoidally can
be viewed as an effectively continuous simulacrum of a sandwich tube. Thus, we
note that a sandwich tube with a hard internal core can resist the internal
pressure inducing aneurysm formation without influence on the critical
stretch. From the viewpoint of a structure optimization, such a composition
reduces material cost but increases the critical pressure giving rise to
aneurysm formation. Bearing in mind that human arteries are composed of a
sandwich structure where the core (intermediate layer) has the largest elastic
modulus (Holzapfel et al.,, 2000; Chirinos,, 2012), our theoretical prediction
agrees well with the natural evolution and optimization of human arteries.

Figure 6: The critical stretch $\lambda_{a}^{c}$ versus the ratio of modulus
when the three modulus functions in (2.19) are adopted. (a) The modulus ratio
ranges from 0.01 to 101 and (b) the modulus ratio ranges from 0.01 to 1. The
right figure is a blowup of the right figure when the modulus ratio is small.
The resultant axial force is given by $N^{*}=0$, and the other parameters are
given by $A^{*}=0.5$, $\zeta=30$ and $J_{m}=97,2$.

###  3.3 LTII: Fixed axial length

Secondly, we address the other commonly adopted loading condition where the
axial length is prescribed by stretching a rubber tube and then fixing two
ends. A noteworthy feature in LTII is that the curve of pressure against
stretch $\lambda_{a}$ or volume ratio $\lambda_{a}^{2}\lambda_{z}$ may lose
$N$-shape and further may be a monotonically increasing function.
Consequently, the limit-point theory fails to identify the bifurcation point,
and we resort to the bifurcation condition derived by Fu et al., (2016).

Previous studies of localized bulging in inflated bilayer tubes have indicated
that the effect of material and geometrical parameters on the bifurcation
threshold directing to localized bulging has no essential distinction (Liu et
al.,, 2019; Ye et al.,, 2019). Therefore, we anticipate that similar
bifurcation results as that in LTI will be observed here. To this end, we
specify the modulus ratio and plot the dependence of $\lambda_{a}^{c}$ on the
dimensionless inner radius $A^{*}$ in Figure 7 while Figure 8 displays
$\lambda_{a}^{c}$ against the modulus ratio $\beta_{1}$ or $\beta_{2}$. As
explained earlier, the axial stretch is given by $\lambda_{z}=1.5$ in all
illustrative examples. For comparison, Figure 7(a) shows the outcomes for
three modulus functions when $\beta_{1}=\beta_{2}=5$ while Figure 7(b)
exhibits the counterparts when $\beta_{1}=\beta_{2}=30$. As expected, the
curve for the sinusoidal function is the highest. In addition, the curve for a
linear function is lower than the one for exponential function as
$\beta_{1}=5$ but becomes higher when $\beta_{1}$ is transferred to 30. These
features are consistent with those in Figure 4, so we omit more detailed
descriptions.

Then we explore the influence of modulus ratio depicted in Figure 8. It can be
seen that the tendencies and shapes of all curves are consistent with those in
Figures 5 and 6. Similarly, we determine the second intersect of the solid
line and the dashed line as $\beta_{c}=13.01$.

(a) $\beta_{1}=\beta_{2}=5$.

(b) $\beta_{1}=\beta_{2}=30$.

Figure 7: The critical stretch $\lambda_{a}^{c}$ versus the dimensionless
inner radius $A^{*}$ when the three modulus functions in (2.19) are adopted.
The axial stretch is given by $\lambda_{z}=1.5$, and the other parameters are
given by $\zeta=30$ and $J_{m}=97,2$. Figure 8: The critical stretch
$\lambda_{a}^{c}$ versus the ratio of modulus when the three modulus functions
in (2.19) are adopted. (a) The modulus ratio ranges from 0.01 to 101 and (b)
the modulus ratio ranges from 0.01 to 1. The right figure is a blowup of the
right figure when the modulus ratio is small. The axial stretch is given by
$\lambda_{z}=1.5$, and the other parameters are given by $A^{*}=0.7$,
$\zeta=30$ and $J_{m}=97,2$.

At present, the influence of the tube thickness, the shear modulus ratio, and
the modulus distribution on bulge initiation is revealed by the use of the
bifurcation condition $J(P^{*},N^{*})=0$ and based on the solution procedure
established at the beginning of this section. We emphasize that in practical
applications not only the bulge formation but also the bulge evolution is of
great interest because material rapture often occurs in the growth or
propagation stage. However, theoretically characterizing the solution path for
bulge growth is extremely difficult. Only weakly nonlinear analysis can be
carried out for thick tubes in the framework of finite elasticity (Ye et al.,,
2020). Unfortunately, the information of weakly nonlinear behavior has limited
guidance in the bulge growth stage since the bifurcation nature of localized
bulging is subcritical, and experimentally observed amplitude is always
finite. To acquire more knowledge of bulge growth in inflated graded tubes, we
shall establish a finite element model in commercial software Abaqus in the
next section.

##  4 Finite element model

In this section, we pursue constructing a finite element (FE) model in Abaqus
to conduct a post-buckling analysis of localized bulging in a pressurized
graded tube. Although the incompressible Gent model has not been built in the
software, one could write UHYPER subroutine codes following the user guideline
(Abaqus,, 2013). In our previous studies (Liu et al.,, 2019; Ye et al.,,
2019), homogeneous Gent model was supplemented into Abaqus employing UHYPER
subroutine coding. Another important issue is to model modulus gradient. In
practice, we can discretize the tube into many sub-tubes in the thickness
direction and each sub-tube is composed of a homogeneous Gent material where
the shear modulus can be determined according to the given modulus function.
This is equivalent to approximating a continuous modulus function by a
staircase one. As the mesh size tends to zero, the deformation behavior of a
layered tube would perfectly coincide with that of the original graded one.
Based on this idea, pattern transition induced by volumetric growth in graded
structures was investigated using FE analysis based on this idea Liu et al.,
(2020). Moreover, another more convenient way is to make use of the linearly
temperature-dependent elastic modulus option in Abaqus (Diab and Kim,, 2014;
Chen et al.,, 2018). In this alternative approach, the temperature
distribution will exactly conform to the modulus distribution and the thermal
expansion coefficient must be set to be zero.

In principle, both approaches can be applied to simulate the deformation and
instability of a graded tube under the combined action of internal pressure
and axial extension in Abaqus. Notwithstanding, it should be pointed out that
the former method may need more grids to guarantee the quality of meshes and
may further require additional computation costs. So we employ the latter
method and write the UHYPER subroutine codes in Fortran for the graded Gent
model by creating a connection between the shear modulus and the temperature
field. Ultimately, a graded Gent model with arbitrary modulus function is well
established in Abaqus. In doing so, we can carry out FE simulations for the
localized bulging.

Before proceeding further, we shall check the validity of the FE model. For
that purpose, we calculate the bifurcation threshold using FE analysis and
compare the results with the corresponding theoretical predictions. To
eliminate the end effects of a grade tube with finite length, we construct all
FE models by specifying the length to overall diameter by 30, or equivalently
$L/B=60$. Furthermore, the built-in command “Buckle” can not be applied to
solve the eigenvalue problem when the critical mode is zero, and we utilize
the Module “Static, Riks” to perform a fully nonlinear analysis.

To trigger localized bulging, certain structures or physical imperfections
need to be involved in the FE model. In the case of fixed axial force, we
refer to the practical end conditions in an inflation experiment where both
the two ends are fastened and the radial movement is confined (Wang et al.,,
2019). So in the FE model, the radial displacements on the two ends are
restricted and this can be regarded as a geometrical imperfection. For fixed
axial length, the shear modulus at the center $Z=0$ is assigned to
$\mu(R)-0.001$ if $\beta_{1}$ (or $\beta_{2}$) is less than 30 or otherwise a
greater magnitude $\mu(R)-0.01$ shall be used. In addition, A quadratic 3D
hybrid element with reduced integration (C3D20RH) is employed in all
subsequent FE simulations.

For convenience, the FE model is validated by taking the linear modulus
function $(\ref{eq2_19})_{1}$ as an example, and all comparisons are
summarized in Figure 9. The theoretical solutions are denoted by solid lines
while the FE results are represented by red dots. The resultant axial force is
specified by $N^{*}=0$ in Figures 9(a) and 9(b) and the axial stretch is fixed
by $\lambda_{z}=1.5$ in Figures 9(c) and 9(d), respectively. It is found that
the critical stretch based on FEA agrees well with that from the theoretical
model. This implies the credibility and robustness of the established FE model
and the UHYPER subroutine codes for the graded Gent model.

(a) The modulus ratio is given by $\beta_{1}=30$.

(b) The inner radius is given by $A^{*}=0.7$.

(c) The modulus ratio is given by $\beta_{1}=30$.

(d) The inner radius is given by $A^{*}=0.7$.

Figure 9: (Color online) Comparisons of the critical stretch $\lambda_{a}^{c}$
between the theoretical predictions and the FE results when the shear modulus
varies linearly. The resultant axial force is fixed by $N^{*}=0$ in top
subfigures while the axial stretch is given by $\lambda_{z}=1.5$ in the bottom
subfigures. The red dots denote the FE solutions and the solid lines
correspond to the theoretical ones.

Now, a robust FE model for simulating localized bulging in pressurized graded
tubes has been formulated. In particular, the modulus gradient in the FE model
can be arbitrary. So the FE model paves a convenient way to implement a fully
nonlinear analysis of bulge evolution. In the next section, both bulge growth
and bulge propagation will be investigated for a graded tube.

##  5 Bulge growth and propagation

It is known that a complete process of localized bulging in pressurized tubes
contains three typical stages, consisting of bulge initiation, bulge growth,
and bulge propagation. When the internal pressure passes a critical value, a
bulge profile appears. With increased pressure, the diameter of the bulge
grows until it reaches a maximum. Then we define a new parameter
$\lambda_{a}^{m}$ that depicts the largest hoop stretch at the inner surface.
As a result, the $\lambda_{a}^{m}A^{*}$ indicates the maximum inner radius for
the bulge. Afterwards, the growing bulge ceases to expand in the radial
direction and only dilates in the axial direction at constant pressure, and
the instability enters into the propagation stage. We emphasize that an ideal
propagation stage only exists in an infinitely long tube, so the maximum
radius cannot be exactly reached in a tube of finite length.

In many actual applications, such as soft robots and flexible actuation, bulge
profile is harmful to the normal function of these devices, so the parametric
analysis for the onset of localized bulging can supply useful insight into
bulge suppression. If localized bulging has emerged, material rupture may
appear either in the growth stage or during the period of bulge propagation.
In practice, the rupture risk is possibly relevant to the size of a bulge or
the magnitude of circumferential stretch. For clinicians, the diameter of an
aneurysm can be used to evaluate the danger of aneurysm fracture (Humphrey and
Holzapfel,, 2012). It is therefore of great significance to estimate how large
a bulge can finally attain. As noted earlier, human arteries are composed of
three layers where the intermediate layer has the highest elastic modulus, and
we can approximately take a non-monotonic modulus distribution to
qualitatively elucidate the mechanism behind artery evolution. Hence, the
section aims to reveal the effect of modulus ratio and material gradient on
the maximum magnitude of the bulge and to provide further insight into how a
normal artery is constructed to resist aneurysm formation while keeping proper
elasticity.

Figure 10: (Color online) The dependence of pressure $P^{*}$ on volume ratio
$v$ for the primary deformation subjected to inflation of a graded tube where
the sinusoidal distribution $(\ref{eq2_19})_{3}$ is applied. The resultant
axial force is fixed by $N^{*}=0$ and the parameters are give by $A^{*}=0.7$
and $\beta_{2}=30$.

It is pointed out that the propagation stage is a two-phase deformation and is
composed of two uniform states joined by a transition zone. When the resultant
axial force is fixed, the curve of pressure versus $v$ or $\lambda_{a}$
presents an $N$-shape where the first stationary point corresponds to the
onset of localized bulging. The $v=\lambda_{a}^{2}\lambda_{z}$ denotes the
ratio of the volume enclosed in the undeformed tube to that in the inflated
state. In this scenario, the pressure under which a bulge starts to propagate
can be determined according to Maxwell’s equal-area rule (Chater and
Hutchinson,, 1984). To clearly illustrate this methodology, we consider a
graded tube with sinusoidally varying shear modulus subjected to inflation,
and the resultant axial force is fixed by $N^{*}=0$. Figure 10 sketches the
relation between the dimensionless pressure $P^{*}$ and the volume ratio $v$
when the dimensionless inner radius $A^{*}$ is set to be 0.7. The marked
pressure $P_{k}^{*}$ depicts the origination of bulge propagation. Maxwell’s
equal-area rule yields

$\displaystyle
P^{*}_{k}(v_{1}-v_{0})=\int_{v_{0}}^{v_{1}}P^{*}(v)\textrm{d}v,$ (5.1)

where $v_{0}$ corresponds to a state before localized bulging and $v_{1}$ the
counterpart after bulge has propagated two both ends. Referring to the
solution strategy outlined in Section 3.1, we can also program in
$Mathematica$ to identify the exact values of $P_{k}^{*}$, $v_{0}$ and
$v_{1}$, and the technicality is neglected for brevity. For instance, for the
parameters used in plotting Figure 10, we acquire that $P_{k}^{*}=3.318$,
$v_{0}=1.575$ and $v_{1}=192.87$. However, when the axial length is specified,
the pressure versus $v$ or $\lambda_{a}$ may become monotonic, so Maxwell’s
equal-area rule is no longer valid. Notwithstanding, it is naturally expected
that the principal conclusion for fixed axial force can offer qualitatively
guidance for other loading methods as usual.

Figure 11: Illustration of dependences of $\lambda_{a}^{m}$ on $\beta_{1}$
when the linear and exponential functions are employed or $\beta_{2}$ when the
sinusoidal function is applied. The resultant axial force is given by
$N^{*}=0$ and the dimensionless inner radius $A^{*}$ is prescribed by $0.7$.

In the subsequent analysis, we shall illuminate the influence of the modulus
ratio as well as the material gradient on the maximum circumferential stretch
$\lambda_{a}^{m}$ based on equation (5.1). Similar to the parametric analysis
for the onset of localized bulging, we also plot in Figure 11 the curves for
three modulus functions in (2.19) when the modulus ratio $\beta_{1}$ or
$\beta_{2}$ ranges from 1 to 101 and the dimensionless inner radius is given
by $A^{*}=0.7$. This implies that the shear modulus of the inner (or middle)
surface is equal to or larger than that of the outer surface. Moreover, it is
a wonder that the trend of each curve is identical to that in Figures 5 and 6.
Similarly, the sinusoidal distribution is highest among these curves,
indicating that a localized bugle can attain a greater radius in this case
compared to a linearly decayed or an exponentially decayed modulus
distribution if the peak value of shear modulus and the thickness remain the
same. Meanwhile, the curve of sinusoidal distribution increases slightly from
$\beta_{2}=1$ and soon becomes virtually flat as $\beta_{2}$ is around 10.
Consequently, such a special non-monotonic modulus gradient brings a
negligible change in the final deformed state. Note that these curves
intersect at $\beta_{1}=\beta_{2}=1$ where a homogeneous counterpart is
obtained. Furthermore, we only focus on the linear and exponential
distributions since the unique difference is the modulus gradient. It can be
seen that both curves are monotonically decreasing as $\beta_{1}>1$.
Especially, within the domain $1<\beta_{1}<\beta_{s}$ where $\beta_{s}=8.67$,
a pressurized graded tube with the shear modulus distributed exponentially can
finally undergo a larger bulge compared to the linear counterpart if the axial
force is specified. On the other hand, with increased $\beta_{1}$, the
dependence of the final magnitude of a bulge on the modulus ratio becomes
weaker and weaker until the modulus ratio has no influence.

Figure 12: (Color online) Graphical interpretation of the scaled amplitude
$\lambda_{a}(0)-\lambda_{\infty}$. The white horizontal line marks the axis of
the tube.

As of now, we have performed a theoretical analysis on the limit magnitude of
the bulge in a graded tube subject to combined fixed axial force and internal
pressure in the light of Maxwell’s equal area. An interesting discovery is
that a sinusoidal modulus function does practically not alter the final size
of a bulge compared to its homogeneous counterpart. Bearing in mind that a
similar situation takes place for the bulge initiation, we are left to
speculate what will happen in the bulge growth stage for a homogeneous tube
and a graded tube where the shear modulus distributes sinusoidally. To this
end, we apply the established FE model in the previous section and carry out a
fully nonlinear analysis. At first, we define the amplitude of localized
bulging. As illustrated in Figure 12, the hoop stretch is dependent not only
on the radial coordinate $R$ but also on the axial coordinate $Z$ after bulge
initiation. In a perfect system, namely, the tube is infinitely long, bulge
solution holds translation invariance. So we can assume that localized
bulging, associated with the greatest diameter of the bulge, always appears at
the center of the tube $Z=0$. Then the dimensionless radius at the inner
surface and $Z=0$ writes $\lambda_{a}(0)A^{*}$. Furthermore, we denote the
corresponding radius for the uniform inflation region by
$\lambda_{\infty}A^{*}$, where $\lambda_{\infty}$ stands for the hoop stretch
at the infinity. In a uniform inflation process, the $\lambda_{a}$ is
independent of $Z$ hence we get $\lambda_{a}=\lambda_{\infty}$. In an
illustrative FE example, we identify $\lambda_{\infty}$ in the non-bulging
region, as shown in Figure 12. The scaled amplitude of a bulge is given by
$\lambda_{a}(0)-\lambda_{\infty}$.

On specifying $A^{*}=0.7$, we show the bifurcation diagrams for a homogeneous
tube and a graded one in Figure 13 using FE analysis. The modulus ratio for
the graded tube is given by $\beta_{2}=30$, and the $\beta_{2}$ is reduced to
unity in a homogeneous tube. In particular, Figure 13(a) displays the results
when the resultant axial force is prescribed while Figure 13(b) plots the
counterparts when the axial length is fixed. From the bifurcation diagram, it
is observed that the amplitude of bulge is trivial in the uniform inflation
where $\lambda_{\infty}$ is reduced to the circumferential stretch in the
primary deformation $\lambda_{a}$. As the inflation continues, especially when
$\lambda_{\infty}$ passes the critical value $\lambda_{a}^{c}$, a non-trivial
solution resulting in localized bulging takes place and the scaled amplitude
$\lambda_{a}(0)-\lambda_{\infty}$ gets a positive value. Afterwards, the
amplitude increases with reduced $\lambda_{\infty}$. This implies the fact
that the existing bulge dilates in both the radial and the axial directions
while the non-bulging part shrinks. Remarkably, we again find an extremely
good coincidence for the deformation processes between the homogeneous tube
and the graded one no matter it is the resultant axial force $N^{*}$ or the
axial stretch $\lambda_{z}$ that is fixed. Combining the analysis for the
bulge initiation and the maximum size of a bulge, it is therefore summarized
that a graded tube with sinusoidally distributed shear modulus can undergo a
higher internal pressure and remain the deformation feature at the same time.
For instance, the deformation paths in Figure 13(a) are identical while the
associated critical pressures for a graded tube and a homogeneous one are
$P_{cr}^{*}=5.248$ and $P_{cr}^{*}=0.271$, respectively.

(a) The resultant axial force is fixed by $N^{*}=0$.

(b) The axial length is fixed by $\lambda_{z}=1.5$.

Figure 13: (Color online) Comparisons of the bifurcation diagrams between a
graded tube where the shear modulus function is given by $(\ref{eq2_19})_{3}$
and a homogeneous counterpart. The parameters are given by $A^{*}=0.7$ and
$\beta_{2}=30$. The black solid lines represent the results for the
homogeneous tube while the red dashed lines the counterparts for the graded
one. The profiles of these blue dots highlighted in the figure are displayed
in Figures 14 and 15, respectively.

We also highlight three pairs of
$\\{\lambda_{\infty},\lambda_{a}(0)-\lambda_{\infty}\\}$ by blue dots for the
graded tube in Figures 13(a) and 13(b), respectively. In particular, we select
a state before instability and two states after bulge initiation. The
corresponding coordinates
$\\{\lambda_{\infty},\lambda_{a}(0)-\lambda_{\infty}\\}$ read $\\{1.4,0\\}$,
$\\{1.267,4.004\\}$ and $\\{1.244,5.373\\}$ from the bottom to the top in
Figure 13(a) while in Figure 13(b) the counterparts are given by
$\\{1.4,0\\}$, $\\{1.209,3.253\\}$ and $\\{1.176,4.771\\}$. The deformed
profiles for the three highlighted points are exhibited in Figure 13(a). In
the right part, we further plot the cross-sections enclosed in the dashed
rectangular to offer a clear view of the stress distribution on the inner
surface. The aspect ratio of length to diameter is 30. As introduced in
Section 4, the radial displacement of two ends are restricted and this setup
serves as a geometric imperfection that can be used to trigger aneurysm
formation at a critical pressure. It can be seen that the end effect decays
very fast and has no influence on the deformation around the center where
localized bulging is expected to appear. As the stretch passes the critical
stretch $\lambda_{a}^{c}$, a bulge occurs at the center of the tube and starts
to expand in the radial and axial orientations. At the same time, the inner or
outer diameter at the non-bulging area (outside the dashed rectangular in
Figure 14) suffers a small but continuous drop, which explains why the curve
in the amplitude diagram turns back when $\lambda_{\infty}$ exceeds the
critical stretch $\lambda_{a}^{c}$. From the sectional view, it is also found
that the stress and strain are both concentrated on the center of the bulge,
which may cause potential rapture.

Figure 14: (Color online) Deformed profiles of the three marked points in
Figure 13(a) for the graded tube when the axial force is fixed by $N^{*}=0$.
(a) to (c) show the entire deformations at different periods while (d) to (e)
illustrate the sectional views. The values of $\lambda_{\infty}$ and
$\lambda_{a}(0)-\lambda_{\infty}$ are listed below the corresponding
subfigures. The parameters are given by $A^{*}=0.7$ and $\beta_{2}=30$.

Finally, we plot in Figure 15 the deformed configurations of all blue points
in Figure 13(b). The layout is similar to that in Figure 14. We emphasize that
an alternative physical defect that the shear modulus at the middle cross-
section $Z=0$ is reduced to $\mu(R)-0.01$ is employed in the case of fixed
axial length. In doing so, the center of the tube is slightly softer than
other positions. Then localized bulging occurs at the center if the internal
pressure attains a critical value. Additionally, the cross-sections in the
dashed rectangular are shown to depict evidently the stress distribution on
the inner surface.

Figure 15: (Color online) Deformed profiles of the three marked points in
Figure 13(b) for the graded tube when the axial stretch is fixed by
$\lambda_{z}=1.5$. (a) to (c) show the entire deformations at different
periods while (d) to (e) illustrate the sectional views. The values of
$\lambda_{\infty}$ and $\lambda_{a}(0)-\lambda_{\infty}$ are listed below the
corresponding subfigures. The parameters are given by $A^{*}=0.7$ and
$\beta_{2}=30$.

##  6 Concluding remarks

Localized bulging or aneurysm formation in inflated graded cylindrical tubes
of arbitrary thickness was investigated in detail within the framework of
finite elasticity. The primary deformation where a tube dilates uniformly was
analytically determined to formulate the expressions of the internal pressure
and the resultant axial force. Meanwhile, a concise bifurcation condition in
terms of the internal pressure and the resultant axial force has been applied
and a semi-analytical solution procedure was proposed to determine the onset
of localized bulging. Two typical loading conditions, defined by fixed axial
force and fixed axial length, were adopted throughout the paper. In
particular, the case of fixed axial force has potential applications in soft
pneumatic actuators (He et al.,, 2020; Li et al.,, 2022) while the other end
condition corresponds to the in vivo status of human arteries (Horný et al.,,
2014). For all illustrative examples, the incompressible Gent material was
engaged and three archetypal modulus gradients were used, containing a linear,
an exponential, and a sinusoidal function. Then a thorough theoretical
analysis was conducted for the critical stretch $\lambda_{a}^{c}$ at the inner
surface. Note that the structure is more stable with a higher
$\lambda_{a}^{c}$ in a volume control problem. In this sense, a thicker tube
is always less susceptible to suffering localized bulging regardless of a
concrete form of the material gradient. Remarkably, it turns out that the
critical stretch is sensitive to the specific modulus distribution applied. In
particular, a sinusoidal function almost does not influence $\lambda_{a}^{c}$,
no matter what values of the ratio of the shear modulus on the middle position
to that on the outer surface (denoted by $\beta_{2}$). Furthermore, a linear
or an exponential one will reduce $\lambda_{a}^{c}$ as the ratio of the shear
modulus on the inner surface to that on the outer surface (denoted by
$\beta_{1}$) increases. However, a larger enough $\beta_{1}$ ceases to affect
the critical stretch. It is emphasized that a qualitatively similar conclusion
can be obtained for both end conditions. The current analysis on bulge
initiation implies that not only the modulus mismatch but also the exact
position where maximum modulus attains is vital to the initiation of localized
bulging. Bearing in mind that an accurate measurement of elastic modulus of
human arteries is extremely difficult (Holzapfel et al.,, 2019) and we only
have a general understanding that the media (intermediate layer) of an artery
occupies a larger elastic modulus Holzapfel et al., (2000); Chirinos, (2012),
the current analysis concerning sinusoidally distributed shear modulus implies
that the deformation is almost insensitive to the practical value of the shear
modulus of the middle surface.

To track the evolution of a bulge, a finite element model incorporating
material inhomogeneity is established using Abaqus UHYPER subroutine coding.
This model can be used to perform nonlinear analysis for localized bulging in
inflated graded tubes for an arbitrary modulus gradient. The robustness of the
established FE model was validated by comparing the bulge initiations using
finite element analysis to the theoretical predictions when a linear modulus
gradient is applied. Ultimately, bulge propagation for the case of fixed axial
force was analytically studied according to Maxwell’s equal-area rule, and the
influence of material gradient, as well as the modulus ratio on the maximum
magnitude of a bulge, were clarified. Combing the analysis for the onset of
localized bulging, it is found that either the material gradient or the
modulus ratio has the same effect on the deformed profiles at the critical
state and at the propagation stage. Furthermore, we compared the bifurcation
diagrams between a graded tube where the shear modulus decays sinusoidally
from the middle surface to both lateral boundaries and its homogeneous
counterpart based on finite element simulations. Interestingly, these two
bifurcation diagrams are nearly identical for both loading types. Therefore,
this work provides an answer for one fundamental problem arising from
structure optimization, i.e. stiffening the intermediate surface enhances the
ability of a structure to resist internal pressure, and neither the critical
stretch nor the deformation path varies compared with the homogeneous
counterpart. This conclusion is perfectly consistent with the actual situation
of human arteries. It is expected that the current analysis would supply
useful insight into localization instabilities in graded structure as well as
into aneurysm prevention in functional soft devices.

## Acknowledgments

The work was supported by the National Natural Science Foundation of China
(Grant Nos 12072227, 12072225, and 12121002). The Abaqus simulations were
carried out on TianHe-1 (A) at the National Supercomputer Center in Tianjin,
China. We thank Prof. Yibin Fu at Keele University for valuable advice and
discussion.

## References

  * Abaqus, (2013) Abaqus, 2013. ABAQUS Analysis User¡¯s Manual, Version 6.13. Dassault Systèmes, Providence, RI, USA.
  * Alexander, (1971) Alexander, H., 1971. Tensile instability of initially spherical balloons. Int. J. Eng. Sci. 9, 151–160.
  * Alhayani et al., (2013) Alhayani, A.A., Giraldo, J.A., Rodríguez, J., Merodio, J., 2013. Computational modelling of bulging of inflated cylindrical shells applicable to aneurysm formation and propagation in arterial wall tissue. Finite Elem. Anal. Des. 73, 20–29.
  * Alhayani et al., (2014) Alhayani, A.A., Rodríguez, J., Merodio, J., 2014. Competition between radial expansion and axial propagation in bulging of inflated cylinders with application to aneurysms propagation in arterial wall tissue. Int. J. Eng. Sci. 85, 74–89.
  * Batra, (2011) Batra R.C., 2011. Material tailoring and universal relations for axisymmetric deformations of functionally graded rubberlike cylinders and spheres. Math. Mech. Solids 16(7), 729–738.
  * Batra and Bahrami, (2009) Batra, R. C., Bahrami, A., 2009. Inflation and eversion of functionally graded non-linear elastic incompressible circular cylinders. Int. J. Non-Linear Mech. 44, 311–323.
  * Bucchi and Hearn, (2013) Bucchi, A., Hearn., G.E., 2013. Delay or removal of aneurysm formation in the Anaconda wave energy extraction device. Renew. Energ. 55, 104–119.
  * Chater and Hutchinson, (1984) Chater, E., Hutchinson, J.W., 1984. On the propagation of bulges and buckles. ASME J. Appl. Mech. 51, 269–277.
  * Chen et al., (2017) Chen, W., Liu D., Kitipornchai, S., Yang, J., 2017. Bifurcation of pressurized functionally graded elastomeric hollow cylinders. Compos. Part B-Eng. 109, 259–276.
  * Chen et al., (2018) Chen, Z., Zhang, X., Song, J., 2018. Surface wrinkling of an elastic graded layer. Soft Matter 14, 8717–8723.
  * Chirinos, (2012) Chirinos, J. A., 2012. Arterial stiffness: Basic concepts and measurement techniques. J. of Cardiovasc. Trans. Res. 5, 243–255.
  * Demirkoparan and Merodio, (2017) Demirkoparan, H., Merodio, J., 2017. Bulging bifurcation of inflated circular cylinders of doubly fiber-reinforced hyperelastic material under axial loading and swelling. Math. Mech. Solids. 22, 666–682.
  * Diab and Kim, (2014) Diab, M., Kim, K.-S., 2014. Ruga-formation instabilities of a graded stiffness boundary layer in a neo-Hookean solid. Proc. R. Soc. A 470, 20140218.
  * (14) Emery, D., Fu, Y.B., 2021. Localised bifurcation in soft cylindrical tubes under axial stretching and surface tension. Int. J. Solids Struct. 219–220, 23–33.
  * (15) Emery, D., Fu, Y.B., 2021. Post-bifurcation behaviour of elasto-capillary necking and bulging in soft tubes. Proc. R. Soc. A 477, 20210311.
  * Fu et al., (2008) Fu, Y.B., Pearce, S.P., Liu, K.-K., 2008. Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation. Int. J. Non-linear Mech. 43, 697–706.
  * Fu and Xie, (2010) Fu, Y.B., Xie, Y.X., 2010. Stability of localized bulging in inflated membrane tubes under volume control. Int. J. Eng. Sci. 48, 1242–1252.
  * Fu et al., (2012) Fu, Y.B., Rogerson, G.A., Zhang, Y.T., 2012. Initiation of aneurysms as a mechanical bifurcation phenomenon. Int. J. Non-linear Mech. 47, 179–184.
  * Fu and Xie, (2012) Fu, Y.B., Xie, Y.X., 2012. Effects of imperfections on localized bulging in inflated membrane tubes. Phil. Trans. R. Soc. A. 370, 1896–1911.
  * Fu and II’ichev, (2015) Fu, Y.B., II’ichev, A.T. 2015. Localized standing waves in a hyperelastic membrane tube and their stabilization by a mean flow. Maths. Mech. Solids. 20,1198–1214.
  * Fu et al., (2016) Fu, Y.B., Liu, J.L., Francisco, G.S., 2016. Localized bulging in an inflated cylindrical tube of arbitrary thickness - the effect of bending stiffness. J. Mech. Phys. Solids. 90, 45–60.
  * Fu et al., (2018) Fu, Y.B., Xie, Y.X., Dorfmann, L., 2018. Localized necking of a dielectric membrane. Extreme Mech. Lett. 21, 44–48.
  * Fu et al., (2021) Fu, Y.B., Jin, L., Goriely, A., 2021. Necking, beading, and bulging in soft elastic cylinders. J. Mech. Phys. Solids. 147, 104250.
  * Gasser et al., (2006) Gasser, T.C., Ogden, R.W., Holzapfel, G.A., 2006. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface. 3, 15–35
  * Gent, (1996) Gent, A., 1996. A new constitutive relation for rubber. Rubber Chem. Technol. 69, 59–61.
  * Giudici and Biggins, (2020) Giudici, A., Biggins, J.S., 2020. Ballooning, bulging and necking: an exact solution for longitudinal phase separation in elastic systems near a critical point. Phys. Rev. E 102, 033007.
  * Goncalves et al., (2008) Goncalves, P.B., Pamplona, D.C., Lopes, S.R.X., 2008. Finite deformations of an initially stressed cylindrical shell under internal pressure. Int. J. Mech. Sci. 50, 92–103.
  * Guo et al., (2016) Guo, Z., Gattas, J., Wang S., Li, L., Albermani F., 2016. Experimental and numerical investigation of bulging behaviour of hyperelastic textured tubes. Int. J. Mech. Sci. 115–116, 665–675.
  * Haughton and Ogden, 1979 (a) Haughton, D.M., Ogden, R.W. 1979. Bifurcation of inflated circular cylinders of elastic material under axial loading-I. Membrane theory for thin-walled tubes. J. Mech. Phys. Solids. 27, 179–212.
  * Haughton and Ogden, 1979 (b) Haughton, D.M., Ogden, R.W. 1979. Bifurcation of inflated circular cylinders of elastic material under axial loading-II. exact theory for thick-walled tubes. J. Mech. Phys. Solids. 27, 489–512.
  * Haughton and Merodio, (2009) Haughton, D.M., Merodio, J., 2009. The elasticity of arterial tissue affected by Marfan’s syndrome. Mech. Res. Commun. 36, 659–668.
  * He et al., (2020) He, Q., Zheng, Y., Wang, Z., He, X., Cai, S., 2020. Anomalous inflation of a nematic balloon, J. Mech. Phys. Solids. 142, 104013.
  * Hejazi et al., (2021) Hejazi, M., Hsiang, Y., Srikantha Phani, A., 2021. Fate of a bulge in an inflated hyperelastic tube: theory and experiment. Proc. R. Soc. A 477, 20200837.
  * Holzapfel et al., (2000) Holzapfel, G.A., Gassera, T.C., Ogden, R.W., 2000. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elasticity. 61, 1–48.
  * Holzapfel et al., (2019) Holzapfel, G.A., Ogden, R.W., Sherifova, S., 2019. On fibre dispersion modelling of soft biological tissues: a review. Proc. R. Soc. A 475, 20180736.
  * Horný et al., (2014) Horný, L., Netuśil, M., Vonǎvková, T., 2014. Axial prestretch and circumferential distensibility in biomechanics of abdominal aorta, Biomech. Model. Mechanobiol. 13(4), 783–799.
  * Horný and Netušil, (2016) Horný, L., Netušil, M., 2016, How does axial prestretching change the mechanical response of nonlinearly elastic incompressible thin-walled tubes. Int. J. Mech. Sci. 106, 95–106.
  * Humphrey and Holzapfel, (2012) Humphrey, J.D., Holzapfel, G.A., 2012. Mechanics, mechanobiology, and modeling of human abdominal aorta and aneurysms. J. Biomech. 45 (5), 805–814.
  * Il’ichev and Fu, (2012) Il’ichev, A.T., Fu, Y.B., 2012. Stability of aneurysm solutions in a fluid-filled elastic membrane tube. Acta Mech. Sin. 28, 1209–1218.
  * Jha et al., (2013) Jha, D.K., Kant, T., Singh R.K., 2013. A critical review of recent research on functionally graded plates. Compos. Struct. 96, 833–849.
  * Jin et al., (2021) Jin, L., Forte, A.E., Bertoldi, K., 2021. Mechanical Valves for On-Board Flow Control of Inflatable Robots. Adv. Sci. 2021, 2101941.
  * Kamalarasa and Calladine, (1988) Kamalarasa, S., Calladine, C., 1988. Buckle propagation in submarine pipelines. Int. J. Mech. Sci. 30, 217–228.
  * Kanner and Horgan, (2007) Kanner, L.M., Horgan, C.O., 2007. Elastic instabilities for strain-stiffening rubber-like spherical and cylindrical thin shells under inflation. Int. J. Non-linear Mech. 42, 204–215.
  * Ke et al., (2008) Ke, L.L., Yang, J., Kitipornchai, S., Wang, Y.S., 2008. Electro-mechanical frictionless contact behavior of a functionally graded piezoelectric layered half-plane under a rigid punch. Int. J. Solids Struct. 45, 3313–3333.
  * Kohn et al., (2015) Kohn, J.C., Lampi, M.C., Reinhart-King, C.A., 2015. Age-related vascular stiffening: causes and consequences. Front. Genet. 6, 112.
  * Kyriakides, (2007) Kyriakides, S., 2007. Mechanics of offshore pipelines. Oxford, UK: Elsevier Science Ltd.
  * Kyriakides and Chang, (1990) Kyriakides, S., Chang, Y.-C., 1990. On the inflation of a long elastic tube in the presence of axial load. Int. J. Solids Struct. 26, 975–991.
  * Kyriakides and Chang, (1991) Kyriakides, S., Chang, Y.-C., 1991. The initiation and propagation of a localized instability in an inflated elastic tube. Int. J. Solids Struct. 27, 1085–1111.
  * Lestringant and Audoly, (2018) Lestringant, C., Audoly, B., 2018. A diffuse interface model for the analysis of propagating bulges in cylindrical balloons. Proc. R. Soc. A 474, 20180333.
  * Li et al., (2022) Li, K., Wang, Q., Xu, P., 2022. Inflation-induced torsion and bulging of a nematic elastomer balloon. Thin Wall. Struct. 170, 108621.
  * Lin et al., (2020) Lin, Z., Li, L., Ye, Y., 2020. Numerical simulation of localized bulging in an inflated hyperelastic tube with fixed ends. Int. J. Appl. Mech. 12, 20501185.
  * Liu et al., (2019) Liu, Y., Ye, Y., Althobaiti, A., Xie Y.-X., 2019. Prevention of localized bulging in an inflated bilayer tube. Int. J. Mech. Sci. 153–154, 359–368.
  * Liu et al., (2020) Liu, Y., Zhang, Z., Devillanova, G., Cai, Z., 2020. Surface instabilities in graded tubular tissues induced by volumetric growth. Int. J. Non-linear Mech. 127, 103612.
  * Lu et al., (2015) Lu, T., An, L., Li, J., Yuan, C., Wang, T.J., 2015. Electro-mechanical coupling bifurcation and bulging propagation in a cylindrical dielectric elastomer tube. J. Mech. Phys. Solids. 85, 160–175.
  * Mallock, (1891) Mallock, A., 1891. Note on the instability of India-rubber tubes and balloons when distended by fluid pressure. Proc. R. Soc. Lond. A. 49, 458–463.
  * Mao et al., (2014) Mao, G.Y., Li, T.F., Zou, Z.N., Qu, S.X., Shi, M.X., 2014. Prestretch effect on snap-through instability of short-length tubular elastomeric balloons under inflation. Int. J. Solids Struct. 51, 2109–2115.
  * Merodio and Haughton, (2010) Merodio, J., Haughton, D.M., 2010. Bifurcation of thick-walled cylindrical shells and the mechanical response of arterial tissue affected by Marfan’s syndrome. Mech. Res. Commun. 37, 1–6.
  * Pamplona et al., (2006) Pamplona, D.C., Goncalves, P.B., Lopes, S.R.X., 2006. Finite deformations of cylindrical membrane under internal pressure. Int. J. Mech. Sci. 48, 683–696.
  * Pearce and Fu, (2010) Pearce, S.P., Fu, Y.B., 2010. Characterisation and stability of localised bulging/necking in inflated membrane tubes. IMA J. Appl. Math. 75, 581–602.
  * Pearce, (2012) Pearce, S.P., 2012. Effect of strain-energy function and axial prestretch on the bulges, necks and kinks forming in elastic membrane tubes. Math. Mech. Solids.17, 860–875.
  * Reddy and Saxena, (2018) Reddy, N.H., Saxena, P., 2018. Localized bulging of rotating elastic cylinders and tubes. Instabilities in the axisymmetric magnetoelastic deformation of a cylindrical membrane. Int. J. Solids Struct. 136–137, 203–219.
  * Rodríguez and Merodio, (2011) Rodríguez, J., Merodio, J., 2011, A new derivation of the bifurcation conditions of inflated cylindrical membranes of elastic material under axial loading. Application to aneurysm formation. Mech. Res. Commun. 38, 203–210.
  * Rodríguez-Martínez et al., (2015) Rodríguez-Martínez, J.A., Fernández-Sáez, J., Zaera, R., 2015. The role of constitutive relation in the stability of hyper-elastic spherical membranes subjected to dynamic inflation. Int. J. Eng. Sci. 93, 31–45.
  * Sakalihasan et al., (2018) Sakalihasan, N., Michel, J.B., Katsargyris, A., Kuivaniemi, H., Defraigne, J.O., Nchimi, A., Powell, J.T., Yoshimura, K., Hultgren, R., 2018. Abdominal aortic aneurysms. Nat. Rev. Dis. Prim. 4, 1–22.
  * Seshaiyer and Humphrey, (2001) Seshaiyer, P., Humphrey, J.D., 2001. On the protective role of contact constraints in saccular aneurysms, J. Biomech. 34, 607–612.
  * Takla, (2018) Takla, M., 2018. Insight into elastic-plastic bifurcation of pressurized cylinders: Transition between bulging and necking; the line of catastrophic failure. Int. J. Mech. Sci. 148, 73–83.
  * Varatharajan and DasGupta, (2017) Varatharajan, N., DasGupta, A., 2017. Study of bifurcation in a pressurized hyperelastic membrane tube enclosed by a soft substrate. Int. J. Non-linear Mech. 95, 233–241.
  * Varatharajan and DasGupta, (2018) Varatharajan, N., DasGupta, A., 2018. Spectral stability of the bifurcation state of an arterial model with perivascular soft tissues. Mech. Res. Commun. 91, 7–12.
  * Wang, (2020) Wang, L., 2020. Axisymmetric instability of soft elastic tubes under axial load and surface tension. Int. J. Solids Struct. 191–192, 341–350.
  * Wang et al., (2017) Wang, J., Althobaiti, A, Fu Y., 2017. Localized bulging of rotating elastic cylinders and tubes. J. Mech. Mater. Struct.12(4), 545–561.
  * Wang and Fu, (2018) Wang, J., Fu, Y.B., 2018. Effect of double-fibre reinforcement on localized bulging of an inflated cylindrical tube of arbitrary thickness. J. Eng. Math, 109(1), 21–30.
  * Wang et al., (2019) Wang, S.B., Guo, Z.M., Zhou, L., Li, L.A., Fu, Y.B., 2019. An experimental study of localized bulging in inflated cylindrical tubes guided by newly emerged analytical results. J. Mech. Phys. Solids. 124, 536–554
  * Wang et al., (2018) Wang, T., Xu, F., Huo, Y., Potier-Ferry, M., 2018. Snap-through instabilities of pressurized balloons: pear-shaped bifurcation and localized bulging. Int. J. Non-linear Mech. 98, 137–144.
  * Mathematica, (2019) Wolfram Research Inc. 2019. Mathematica: version 12, Wolfram Research Inc. Champaign, IL.
  * Yu and Fu, (2021) Yu, X., Fu, Y., 2021. An analytic derivation of the bifurcation conditions for localization in hyperelastic tubes and sheets. arXiv:2109.04036v3.
  * Ye et al., (2019) Ye, Y., Liu, Y., Althobaiti, A., Xie, Y.-X., 2019. Localized bulging in an inflated bilayer tube of arbitrary thickness: Effects of the stiffness ratio and constitutive model. Int. J. Solids Struct. 176–177, 173–184.
  * Ye et al., (2020) Ye, Y., Liu, Y., Fu, Y, 2020. Weakly nonlinear analysis of localized bulging of an inflated hyperelastic tube of arbitrary wall thickness. J. Mech. Phys. Solids 135, 103804.
  * Zhong et al., (2012) Zhong, Z., Wu, L., Chen, W., 2012. Mechanics of functionally graded materials and structures. New York: Nova Science Publishers.