Compressibility-regulated elasto-capillary necking and bulging in soft solids

2.1. Three-dimensional finite-strain bulk–surface model

Consider a soft hyperelastic cylinder (with radius A and length 2 L in its reference configuration; A ≪ L ) subjected to axial loadings at its two ends and surface stress on its lateral surface. The cylindrical coordinate system with basis vectors { e θ , e z , e r } is used to characterize the axisymmetric deformation. Here, ( Θ , Z , R ) and ( θ , z , r ) are the coordinates of a material point in the reference and current configurations, respectively. The bulk deformation is characterized by the mapping from a material point X = Z e Z + R e R in the reference configuration to its current position x = z ( Z , R ) e z + r ( Z , R ) e r . For the lateral surface R = A , we denote the reference position as Y = X | R = A and its corresponding current position as y = x | R = A . The bulk deformation gradient is defined as F = ∂ x ∕ ∂ X . The surface deformation gradient F s is defined by the differentiation of the surface position y with respect to the reference surface coordinates, which can be expressed as the projection of F into the reference tangent space: F s = ∂ y ∂ Y = F | R = A P , where P = I − e R ⊗ e R is the projection tensor with I the 3D identity tensor. By the definitions, the deformation gradients of the bulk F and the surface F s are expressed as

F = ( r R 0 0 0 z Z z R 0 r Z r R ) , F s = ( r ( Z , A ) A 0 0 0 z Z ( Z , A ) 0 0 r Z ( Z , A ) 0 ) ,

(1)

respectively, where the quantities with subscripts in the matrices represent the partial derivatives with respect to the corresponding subscripts. It is clear that F s has a rank of 2, consistent with the two-dimensional (2D) nature of the surface. This matches the physical expectation that the surface is described by two coordinates ( Θ , Z ). The three principal stretches λ 1 , λ 2 , λ 3 for the bulk and the principal stretches λ 1 s , λ 2 s for the surface can be obtained by calculating the eigenvalues of F T F and F s T F s .

The description of surface deformation, decoupled from that of the underlying bulk solid, is realized within the framework of surface elasticity [26,50–54]. Assume that the deformations of the bulk and the surface are characterized by a strain energy function W ( λ 1 , λ 2 , λ 3 ) and a surface energy function Γ ( λ 1 s , λ 2 s ) , respectively. The total potential energy (scaled by 2 π ) of the cylinder in this 3D bulk–surface formation is taken as

E 3 d [ z , r ] = ∫ − L L ( ∫ 0 A ( W ( λ 1 , λ 2 , λ 3 ) − N z Z ) R dR + AΓ ( λ 1 s , λ 2 s ) ) dZ ,

(2)

where N is the resultant axial force per unit reference cross-sectional area, and the three terms involved represent the contributions of bulk energy, surface energy, and load potential, respectively.

2.2. Homogeneous deformation of a compressible cylinder

For the above elasto-capillary problem, the cylinder can undergo homogeneous deformation

z ( Z , R ) = λZ , r ( Z , R ) = aR ,

(3)

where λ and a are the constant axial stretch and radial stretch. For this deformation, the total potential energy (2) is reduced to 2 L Φ ( a , λ ) , where

Φ ( a , λ ) = 1 2 A 2 w ( a , λ ) + AΓ ( a , λ ) − 1 2 A 2 Nλ ,

(4)

with w ( a , λ ) : = W ( a , λ , a ) being the reduced strain energy.

When the potential energy (4) reaches its minimum, the two equilibrium equations can be established through the extremal conditions ∂ Φ ∕ ( ∂a ) = 0 and ∂ Φ ∕ ( ∂λ ) = 0 , which yields

w 1 ( a , λ ) + 2 A Γ 1 ( a , λ ) = 0 ,

(5)

w 2 ( a , λ ) + 2 A Γ 2 ( a , λ ) = N ,

(6)

with w i = ∂w ∕ ∂ λ i and Γ i = ∂Γ ∕ λ i s , i = 1 , 2 . With a loading parameter (axial force N or surface tension γ) prescribed, solving equations (5) and (6) for a and λ will fully determine the homogeneous deformation (3) of a compressible bulk material.

To systematically study the compressibility effect, three loading scenarios are investigated in later calculations: (I) varying N at fixed γ, (II) varying γ at fixed N, and (III) varying γ at fixed ends. Without loss of generality, we adopt the generalized neo-Hookean model for the bulk energy function

W = μ 2 ( I 1 − 3 − 2 ln ⁡ J ) + κ 2 ( 1 2 ( J 2 − 1 ) − ln ⁡ J ) ,

(7)

where μ is the shear modulus, κ = 2 μν ∕ ( 1 − 2 ν ) is the first Lamé parameter with ν ( ∈ [ − 1 , 1 ∕ 2 ] ) being the Poisson’s ratio of the bulk material, I 1 = λ 1 2 + λ 2 2 + λ 3 2 is the first principal invariant of the right Cauchy–Green tensor, and J = λ 1 λ 2 λ 3 is the determinant of the deformation gradient. The surface energy function is taken as

Γ = γ J s + μ s 2 ( I 1 s − 2 − 2 ln ⁡ J s ) + κ s 2 ( 1 2 ( J s 2 − 1 ) − ln ⁡ J s ) ,

(8)

where γ is the residual surface tension, μ s , κ s , I 1 s , and J s are the surface counterparts of μ, κ, I 1 , and J, respectively. More specifically, μ s is the surface shear modulus (or surface elasticity), κ s = 2 μ s ν s ∕ ( 1 − ν s ) represents the surface bulk modulus, ν s is the Poisson’s ratio of the surface, I 1 s = ( λ 1 s ) 2 + ( λ 2 s ) 2 , and J s = λ 1 s λ 2 s . Note that such a surface energy function with the three surface parameters (i.e., the residual surface tension γ, the surface elasticity μ s , and the surface Poisson’s ratio ν s ) yields strain-dependent surface stresses, which is the intrinsic feature of solid-like materials [22,23]. The physical admissibility of the surface parameters ( γ , μ s , ν s ) is governed by thermodynamic stability and the positive definiteness of the surface strain energy. The residual surface tension γ is typically positive, i.e., γ > 0 , representing the initial capillary-like contractile force. The surface shear modulus satisfies μ s ≥ 0 . The model recovers the classical liquid surface tension results as μ s → 0 . Furthermore, for the 2D surface hyperelasticity adopted here, the surface Poisson’s ratio falls within the range − 1 < ν s < 1 . Unlike the 3D bulk where ν = 1 ∕ 2 marks the volume-incompressibility limit ( κ → ∞ in equation (7)), the 2D surface only reaches the area-incompressibility limit as ν s → 1 ( κ s → ∞ in equation (8)). While the thermodynamic stability of the 3D bulk and 2D surface models permits the Poisson’s ratios to be negative, the present study focuses on the physically prevalent regime for soft solids where 0 < ν < 1 ∕ 2 and 0 < ν s < 1 .

In the following, we scale all lengths by A, all stresses by the shear modulus μ, the surface shear modulus μ s , and the surface tension γ by μA . As such, A and μ can be set equal to unity without loss of generality; we use the same symbols to denote the scaled quantities.

In the scenario of fixed γ and varying N, the homogeneous N-λ response can be derived by solving equation (6)

N = F ( λ ) : = ( ν ( 1 − a ( λ ) 4 λ 2 ) λ ( 2 ν − 1 ) − 2 μ s ν s ( a ( λ ) 2 λ 2 − 1 ) λ ( 2 ν s − 1 ) + 2 a ( λ ) γ + 2 ( λ 2 − 1 ) μ s λ + λ − 1 λ ) ,

(9)

wherein the radial stretch a ( λ ) is the implicit function of the radial stretch a derived from (5). The N-λ responses are shown in Figure 1(a), which exhibit up-down-up patterns.

OPEN IN VIEWER

Figure 1.

Homogeneous responses of the bulk–surface system under two distinct loading scenarios for various Poisson’s ratios ν. (a) Axial force N vs. stretch λ for fixed surface stresses. The Maxwell plateaus N = N M (dashed lines) are determined via the energy-equality criterion. (b) The corresponding evolution of radial stretch a ( λ ) for the curves in (a), illustrating the effect of compressibility. (c) Surface tension γ vs. stretch λ for fixed axial force. Note that the Maxwell plateaus γ = γ M (dashed lines) are derived from the total potential energy equality rather than the equal-area rule. (d) The corresponding evolution of radial stretch a ( λ ) for the curves in (c), highlighting the non-monotonic radial response under γ-controlled loading. The dashed curve a = λ − 1 ∕ 2 in (b) or (d) represents the theoretical limit for an incompressible bulk ( ν = 0 . 5 ).

In the scenario of fixed N and varying γ, the homogeneous γ-λ response can be determined analogously

γ = G ( λ ) : = 1 a ( λ ) ( ν ( 1 − a ( λ ) 4 λ 2 ) 2 λ ( 1 − 2 ν ) + μ s ν s ( a ( λ ) 2 λ 2 − 1 ) λ ( 2 ν s − 1 ) − ( λ 2 − 1 ) μ s λ − λ 2 + 1 2 λ + N 2 ) .

(10)

This relation is plotted in Figure 1(c), which can also result in up-down-up patterns if the parameters are set in a suitable range.

The up-down-up patterns observed in Figure 1(a) and (c), characterized by the presence of both a local maximum and a local minimum, signify the emergence of necking and bulging instabilities in the loading scenarios. This non-monotonicity is the mechanical precursor to localized deformation, as it enables the coexistence of two distinct phases at a constant Maxwell plateau. A comparative analysis of different Poisson’s ratios ν reveals that elasto-capillary responses are exceedingly sensitive to bulk compressibility. This sensitivity prompts a deeper investigation into how compressibility influences bifurcation and post-bifurcation behaviors under strain-dependent bulk–surface interactions.

It is worth noting that the radius response a ( λ ) plays a pivotal role in distinguishing the two scenarios. In scenario I (Figure 1(b)), a ( λ ) exhibits a monotonic decrease, where the radial deformation is passively enslaved to the axial stretch. However, in scenario II (Figure 1(d)), a ( λ ) displays a non-monotonic V-shaped response. This is a direct consequence of the radial equilibrium equation (5), which characterizes the coupling among a, λ and the surface tension γ. In scenario II, the surface tension γ varies with the deformation. This variation introduces a non-monotonic path in the radial equilibrium. The Maxwell state is determined by the equality of the total potential energy between two phases, i.e., Φ ( a ( λ 1 ) , λ 1 ) = Φ ( a ( λ 2 ) , λ 2 ) . In scenario II, the traditional Maxwell area rule becomes invalid for compressible materials. This failure stems from the loss of thermodynamic conjugacy between γ and λ. Specifically, because compressibility allows a to evolve independently, the geometric area under the γ-λ curve cannot capture the additional surface work and internal energy changes induced by volumetric deformation.

2.3. One-dimensional reduced model for localized deformation

As it is seen, the present bulk–surface problem may also admit a localized solution. In Zhu et al. [27], a variationally consistent asymptotic method [21,28–31] is employed to establish a 1D reduced model from the 3D bulk–surface model (2) by taking advantage of the small aspect ratio ϵ = A ∕ ( 2 L ) . From the 1D model, we can derive the asymptotically analytical solution for the localized deformation. Here, we present the main results. The derivation details can be referred to in [27].

The 1D reduced model is derived as

E 1D [ λ ] = ∫ − L L ( Φ ( λ ) + 1 2 B ( λ ) λ ′ ( Z ) 2 ) dZ + C ( λ ) λ ′ ( Z ) | − L L ,

(11)

where Φ ( λ ) ≜ Φ ( a ( λ ) , λ ) , noting that a ( λ ) is the radial stretch as an implicit function of the transverse equilibrium equation (5). The gradient moduli B ( λ ) and C ( λ ) are given by

B ( λ ) = A 4 a λ 2 ( ζ 2 − η 2 4 ζ + Γ 2 λA ) ,

(12)

C ( λ ) = − A 3 a λ Γ 2 η 4 ζ .

(13)

with

a λ = da ∕ dλ , ζ = λ w 2 + a Γ 1 ∕ A λ 2 − a 2 , η = a w 2 − λ Γ 1 ∕ A λ 2 − a 2 − Γ 12 + Γ 22 ∕ a λ A , Γ 12 = ∂ 2 Γ ∂ λ 1 s ∂ λ 2 s , Γ 22 = ∂ 2 Γ ( ∂ λ 2 s ) 2 .

(14)

The 1D energy functional reaches its extremum when the axial stretch λ ( Z ) satisfies

d Φ dλ − 1 2 B ′ ( λ ) λ ′ ( Z ) 2 − B ( λ ) λ ″ ( Z ) = 0 .

(15)

In line with the symmetry of the deformation about the assumed bifurcation initiation point Z = 0 , the initial value conditions for (15) are specified as follows [27]

λ ′ ( 0 ) = 0 , λ ( 0 ) = λ 0 .

(16)

The stretch λ 0 can be determined according to the bifurcation stage. By solving the coupled system (5), (15), and (16), we can obtain the corresponding localized solutions for a ( Z ) , λ ( Z ) .

It is pointed out that in [27], the elasto-capillary localized behaviors for incompressible bulk material (with ν = 1 ∕ 2 ) are investigated specifically. For incompressible material, the isochoric condition leads to a = λ − 1 2 and the transverse equilibrium equation (5) vanishes. In such a situation, solving the initial-value problem (15) and (16) is rather straightforward. In the current situation, the bulk material is considered as compressible ( 0 < ν < 1 ∕ 2 ). Coupling with the equation (5) complicates the problem. We will solve the coupled system for ( a , λ ) in the three loading scenarios, respectively, to systematically explore the effect of bulk compressibility on the critical bifurcation and post-bifurcation behaviors of the localized instability induced by the bulk–surface interactions. The 1D asymptotic results will be validated by 3D numerical simulations via the following Rayleigh–Ritz FE scheme.

2.4. Three-dimensional finite-element scheme for compressible bulk materials

In Zhu et al. [27], a Rayleigh–Ritz-based 3D FE scheme was proposed for the elasto-capillary problem of incompressible bulk materials, subjected to the constraint of the incompressibility condition J = 1 . Here, we unbound this incompressibility constraint and present a generalized FE scheme for compressible bulk materials.

In the 3D energy functional (2), the surface term ∫ − L L AΓ ( λ 1 s , λ 2 s ) dZ will induce ill-posedness at the end surface edges ( Z = ± L and R = A ) due to singular surface stress, posing significant challenges for 3D numerical simulation of localized deformations in the bulk–surface system. The ill-posedness can be resolved by incorporating the essential boundary condition z ( ± L , A ) = ± l into (2) using the Lagrange multiplier method, where l denotes half the current axial length of the cylinder. This boundary condition enforces homogeneous deformation far from the localized instability region, eliminating the unbalanced singular stress at the ends.

Given the symmetry of the deformations about Z = 0 , the 3D energy functional with the natural boundary condition replaced by the essential boundary condition is given by

E ~ [ z , r , q ] = ∫ 0 L ( ∫ 0 A WRdR + AΓ ) dZ − ∫ 0 A q ( R ) ( z ( L , R ) − l ) dR ,

(17)

where q ( R ) is the Lagrange multiplier representing the axial traction distribution (force per unit length along the circumference in the reference configuration) in the end section Z = L , satisfying the resultant axial force balance

2 π ∫ 0 A q ( R ) dR = N ⋅ π A 2 .

(18)

The above equation represents that the integral of the axial traction distribution q ( R ) over the cross-section equals the total applied axial force.

The computational domain [ 0 , L ] × [ 0 , A ] is meshed with n × m uniform four-node axisymmetric quadrilateral elements of side h = A ∕ m = L ∕ n . In each element, the fields z ( Z , R ) and a ( Z , R ) = r ( Z , R ) ∕ R are approximated using bilinear shape functions based on nodal values

z ( Z , R ) = z ij Q 1 ( x , y ) + z i + 1 , j Q 2 ( x , y ) + z i + 1 , j + 1 Q 3 ( x , y ) + z i , j + 1 Q 4 ( x , y ) , a ( Z , R ) = a ij Q 1 ( x , y ) + a i + 1 , j Q 2 ( x , y ) + a i + 1 , j + 1 Q 3 ( x , y ) + a i , j + 1 Q 4 ( x , y ) ,

(19)

while the Lagrange multiplier q ( R ) is estimated by linear interpolation

q ( R ) = q j + q j + 1 − q j h ( R − jh ) .

(20)

here, x = ( Z − ih ) ∕ h and y = ( R − jh ) ∕ h are local coordinates; z ij = z ( ih , jh ) , a ij = a ( ih , jh ) , and Q k ( x , y ) ( 1 ≤ k ≤ 4 ) are the shape functions

Q 1 ( x , y ) = ( 1 − x ) ( 1 − y ) , Q 2 ( x , y ) = x ( 1 − y ) , Q 3 ( x , y ) = xy , Q 4 ( x , y ) = ( 1 − x ) y .

(21)

The above interpolating functions are then substituted into the energy functional (17). The integration of each cell is further approximated by the four-point Gaussian product rule. In doing this, we can transform the energy functional (17) into a discrete form (as a function of z ij , a ij , q j )

E ~ [ z ij , a ij , q j ] , 0 ≤ i ≤ n , 0 ≤ j ≤ m .

(22)

When the discrete energy functional (22) reaches its minimum value, its partial derivatives with respect to z ij , a ij are equal to 0. Together with the restriction to q ( R ) in (18) and the boundary conditions z ( 0 , R ) = 0 and z ( L , R ) = l , we obtain 2 ( n + 1 ) ( m + 1 ) + 1 equations

z 0 , j = 0 , ∂ E ~ ∂ z i , j = 0 , z n , j = l , 1 ≤ i ≤ n − 1 , 0 ≤ j ≤ m ,

∂ E ~ ∂ a i , j = 0 , 0 ≤ i ≤ n , 0 ≤ j ≤ m , ∑ j = 0 n ∂ E ~ ∂ z n , j − 1 2 N A 2 = 0 .

(23)

In the above algebraic system, z ij , a ij , 0 ≤ i ≤ n , 0 ≤ j ≤ m are 2 ( n + 1 ) ( m + 1 ) unknowns. In loading scenario I, γ is fixed and N is set as the loading parameter, while in loading scenario II, N is fixed and γ is the loading parameter. In these two scenarios, the current cylinder length l is the last unknown in the system. In particular, in loading scenario III, l is fixed, and the strain softening may occur as the bifurcation proceeds (see section 5.3). Thus, N (or γ) is not monotonically varying. By contrast, the amplitude Λ continues to increase. We add a constraint into the above algebraic system

| a n , m − a 0 , m | = Λ ,

(24)

with the amplitude Λ designated as the loading parameter. In this formulation of 2 ( n + 1 ) ( m + 1 ) + 2 equations, N and γ are treated as two additional unknowns. Such a setup enables us to smoothly trace the complete post-bifurcation regime.

To find the solution for the above discrete system (23), a suitable initial guess is essential. An FE solution can be obtained with the use of a solution of the 1D model as an initial guess and set the previous FE solution as the initial guess for subsequent iterations. By doing so, the 3D numerical solution for the full post-bifurcation regime can be obtained recursively.

In what follows, we employ both the 1D asymptotic model and the 3D FE scheme to systematically investigate the effect of bulk compressibility (characterized by the Poisson’s ratio ν) on elasto-capillary instability across the three loading scenarios and elucidate its interaction with axial loading N, surface effects (characterized by the surface tension γ, the surface elasticity μ s , and the surface Poisson’s ratio ν s ), and geometric effect (characterized by the length-radius ratio L ∕ A ). We remark that in all numerical simulations, geometric dimensions are non-dimensionalized by the initial radius A (effectively setting A = 1 ). For a sufficiently large aspect ratio L ∕ A , the finite-length FE simulations can accurately capture the bifurcation behavior predicted by the 1D gradient model. This is attributed to the localized nature of the deformation, which decays exponentially and renders the boundary effects negligible at the center.

3.1. The compressibility effect on elasto-capillary necking threshold

The bifurcation condition can be determined by solving F ′ ( λ ) = 0 for the initiation stretch λ I and N I = F ( λ I ) . Alternatively, we can conduct a linear bifurcation (or weakly nonlinear) analysis by using the 1D reduced model (i.e., equations (5) and (15)) to deduce the bifurcation condition. In section 1, we conduct a weakly nonlinear analysis and derive the bifurcation condition (41). For the specific case of fixed surface tension, this condition is equivalent to requiring F ′ ( λ ) = 0 .

Figure 2 presents the continuous variation of necking initiation force N I with respect to the bulk Poisson’s ratio ν for different sets of surface parameters. It is found that N I is positively related to ν and the slope of the curve increases with the Poisson’s ratio for each set of surface parameters. This indicates that the ability of the bulk compressibility to affect the necking threshold rises when ν trends to higher level. Comparing the three subfigures, it can be observed that as Poisson’s ratio increases, the distance between the curves in Figure 2(a) is gradually widened, while the curves in Figure 2(b) and (c) tend to intersect with each other. This indicates that the surface tension has more dominant effect on N I if the cylinder is incompressible, while N I is more sensitive to surface elasticity and surface compressibility if the cylinder is compressible. Moreover, it is seen that the effect of surface elasticity μ s on the critical N I is more significant than the surface Poisson’s ratio ν s .

OPEN IN VIEWER

Figure 2.

The continuous variation of necking initiation force N I with respect to the bulk Poisson’s ratio ν with (a) γ = 8 , 9, 10, 11, μ s = 0 . 2 , and ν s = 0 . 2 . (b) γ = 9 , μ s = 0 , 0.1, 0.2, 0.3, and ν s = 0 . 2 . (c) γ = 9 , μ s = 0 . 2 , and ν s = 0 , 0.2, 0.4, 0.6.

Furthermore, we find that there exists a minimum value of γ, denoted by γ min , below which the response F ( λ ) is monotonically increasing and localization cannot occur. We denote the critical stretch at γ min by λ min . This minimum can be deduced by solving

F ′ ( λ ) = 0 , F ″ ( λ ) = 0 ,

(25)

for γ and λ. Note that the derivatives of the function F are related to the derivatives of the implicit function a ( λ ) , which should be derived from the implicit equation (5).

Figure 3 illustrates the continuous variation of γ min with respect to ν under the influence of μ s and ν s . It is seen that if the surface stress is strain-independent (when μ s = 0 ), γ min increases with ν. But unexpectedly, the response of γ min with respect to ν is not monotonically increasing when the surface elasticity is taken into account. γ min achieves a local minimum when the Poisson’s ratio is close to 0.45. This implies that the situations in the PR instabilities at constant and strain-dependent surface stress are rather different. On the contrary, the evolution of γ min with respect to ν remains non-monotonic regardless of the surface compressibility (see Figure 3(b)). We also find that γ min becomes increasingly responsive to alterations in μ s and ν s as the cylinder becomes more compressible.

OPEN IN VIEWER

Figure 3.

The continuous variation of γ min with respect to ν with (a) μ s = 0 , 0.1, 0.2, 0.3 and ν s = 0 . 2 , respectively. (b) μ s = 0 . 2 and ν s = 0 , 0.2, 0.4, 0.6, respectively.

3.2. The compressibility effect on elasto-capillary necking growth

For post-bifurcation, we first give the solution strategy of the 1D reduced model (5), (15), and (16) for tracking the necking amplitude responding to the applied axial force N with γ, μ s , ν s , and ν fixed. The stretch at the necking starting position λ 0 = λ ( 0 ) should be first determined. Indeed, the equilibrium equation (15) can be integrated once after multiplying the integration factor λ ′ ( Z ) , resulting in

Φ ( a , λ ) − 1 2 B ( λ ) λ ′ ( Z ) 2 = C 0 ,

(26)

where C 0 is a constant. By applying the implicit function (5) to the conservation law (26), λ 0 can be obtained by solving

Φ ( a 0 , λ 0 ) = Φ ( a ∞ , λ ∞ ) ,

(27)

where λ ∞ is the stretch at effectively infinity and can be determined by solving N = F ( λ ∞ ) and the implicit function (5) with N prescribed. By iteratively decreasing N from N I to N M , we can track the evolution of the necking amplitude a ∞ − a 0 in the full post-bifurcation regime. Furthermore, when the initial value λ 0 is specified, the solutions a ( Z ) and λ ( Z ) of the 1D ordinary differential system (5), (15), and (16) can be obtained numerically. Thus, the localization morphology of the cylinder at each post-bifurcation stage can be obtained.

By the above procedure, the N- λ 0 response and the necking amplitude varying with N ∕ N M in the post-bifurcation regime are shown in Figure 4. The corresponding 3D FE simulation results are also presented in Figure 4. The outcomes of the 1D model exhibit excellent agreements with the 3D FE results throughout the entire post-bifurcation regime. It can be seen that as the necking proceeds, the axial force N decreases slowly, while the stretch ratio λ 0 and the necking amplitude continue to increase.

OPEN IN VIEWER

Figure 4.

(a) Dependence of N on λ 0 = λ ( 0 ) . (b) Necking amplitudes against the ratio N ∕ N M , with γ = 8 , μ s = 0 . 2 , ν s = 0 . 2 , and ν = 0 . 4 , 0.3, 0.2, 0.1, respectively. The four open circles labeled 1–4 on the curve represent distinct necking stages, for which the corresponding localized solutions a ( Z ) are displayed in Figure 5.

To probe the necking evolution, we define the normalized necking extent as

ℓ = λ t − λ cr λ end − λ cr ,

(28)

where λ t , λ cr , and λ end represent the values of λ 0 at the current, the start, and the end of the necking stage. We take ℓ = 25 % , 50 % , 75 % , and 99 . 8 % (corresponding to the four open circles marked on the curve in Figure 4(a)), respectively, to characterize the necking evolution. The 1D asymptotic solutions for a ( λ ) and the corresponding 3D FE simulations for the four stages (with γ = 8 , μ s = 0 . 2 , ν s = 0 . 2 and ν = 0 . 1 ) are shown in Figure 5. Even for this highly compressible case ( ν = 0 . 1 ), the 1D model results are in near-perfect agreement with the FE simulation outcomes, validating the high accuracy of the 1D asymptotic reduction. It can be found that as necking evolves, the cylinder thins at the center and thickens at the ends, while the necking transition zone grows steeper. The deformed profiles clearly show a plateau near Z = 0 at the late stage of necking, which is a hallmark of localized necking instability in solid materials.

OPEN IN VIEWER

Figure 5.

The 1D asymptotic solutions (solid curves) and the corresponding 3D FE results (scattered dots) for the radial stretch a ( Z ) . Profiles are shown for four distinct necking stages, corresponding to the states marked in Figure 4(a), with γ = 8 , μ s = 0 . 2 , ν s = 0 . 2 , and a highly compressible bulk ν = 0 . 1 . The inset cylinder illustrates the final non-homogeneous necking morphology in the current configuration.

Comparing the curves for different ν in Figure 4, it is seen that compressibility will reduce the necking amplitude at the final stage (denoted as ( a ∞ − a 0 ) max ). We further plot the final amplitude against ν in Figure 6. As can be seen, there exists a particular value of ν (denoted by ν 0 ) allowing the final amplitude to reach its maximum value. For ν < ν 0 , the final amplitude increases with the Poisson’s ratio. Conversely, when ν > ν 0 , the amplitude experiences a sharp decline as the Poisson’s ratio further increases. This non-monotonic behavior originates from the competition between the elasto-capillary localization and the bulk volumetric constraints. Moderate compressibility optimally facilitates radial necking, while near-incompressibility suppresses it via the strict isochoric constraint.

OPEN IN VIEWER

Figure 6.

Dependence of the final necking amplitudes on ν for the case of γ = 8 , 9, 10, 11, respectively.

3.3. Compressibility-regulated guidelines for PR necking

From the above analysis, the bulk compressibility is demonstrated to play a pivotal role in regulating both the necking threshold and post-bifurcation evolution of PR instability in soft cylinders. The critical initiation force N I rises monotonically with ν. The minimum surface tension γ min required for necking shows a non-monotonic dependence on ν when surface elasticity is considered, distinct from the monotonic trend under strain-independent surface stress. The final necking amplitude exhibits a maximum at a critical ν 0 (around 0.45 in the current parameter setting). Overall, bulk compressibility is a primary controlling factor of the instability, and its effects are closely coupled with surface mechanical parameters.

For practical design and regulation of elasto-capillary necking, specific guidelines can be drawn from the above results. To increase the onset barrier for necking instability, raise the bulk Poisson’s ratio toward the nearly incompressible limit. To maximize the post-bifurcation necking amplitude, target the critical Poisson’s ratio ν 0 , and suppress excessive necking by using either nearly incompressible or highly compressible materials.

4.1. The compressibility effect on elasto-capillary bulging condition

Analogously to loading scenario I, there exists an N min in this loading scenario, under which bulging instability cannot occur. Keeping the other parameters unchanged, N min can be deduced by solving

G ′ ( λ ) = 0 , G ″ ( λ ) = 0 ,

(29)

wherein a is the implicit function a ( λ ) from the implicit equation (5).

To demonstrate the effect of bulk compressibility on N min , we fix μ s = ν s = 0 . 2 and derive the continuous N min -ν relationship, as shown in Figure 7(a). It is clear that N min increases monotonically with ν, indicating that a decrease in bulk compressibility (i.e., an increase in ν) elevates the threshold of the fixed axial force required for instability. Moreover, the initiation surface tension γ I is highly sensitive to both the Poisson’s ratio ν and the axial force N. We fix ν = 0 . 4995 (nearly incompressible), 0.25, and 0.05 and vary N continuously to obtain the γ I -N response curves (see Figure 7(b)). It can be observed that when ν ≤ 0 . 25 , the curves exhibit a near-vertical trajectory, indicating that γ I is extremely sensitive to changes in N in this regime: a slight increase in N leads to a sharp rise in γ I . In contrast, when the material approaches incompressibility ( ν → 0 . 5 , e.g., ν = 0 . 4995 ), γ I changes more gradually with increasing N, showing a weaker sensitivity to the axial force.

OPEN IN VIEWER

Figure 7.

Variations of (a) the critical N min with Poisson’s ratio ν, (b) the initiation surface tension γ I with the axial force N, for μ s = ν s = 0 . 2 . The red dashed line represents the γ I -N response when N = N min , revealing the intrinsic relationship between γ I and N across the entire range of Poisson’s ratios.

4.2. The compressibility effect on elasto-capillary bulging growth

In the post-bifurcation phase, the solution strategy for the amplitude is analogous to loading scenario I. We present the γ- λ 0 response and the bulging amplitude evolving with γ ∕ γ M for the cases of N = 2 . 3 , 2.35, 2.4, μ s = 0 . 2 , ν s = 0 . 2 , and ν = 0 . 05 , in Figure 8. It can be observed that the post-bifurcation response is also very sensitive to the fixed axial force N. The γ I ∕ γ M ratio and the bulging amplitude increase significantly with larger N. Notably, for applications where bulging is desired under a small axial force, compressible bulk materials are the preferred choice.

OPEN IN VIEWER

Figure 8.

(a) The γ- λ 0 response and (b) the dependence of bulging amplitudes on γ ∕ γ M , with N = 2 . 3 , 2.35, 2.4, μ s = 0 . 2 , ν s = 0 . 2 , and ν = 0 . 05 . The four open circles marked on the curve in (a) represent distinct bulging stages, for which the corresponding localized solutions a ( Z ) are displayed in Figure 9.

Furthermore, Figure 9 shows the solutions of the radial stretch a ( Z ) at four bulging extents ℓ = 25 % , 50 % , 75 % , and 99 . 8 % (corresponding to the four open circles marked in Figure 8). It is seen that the bulging front remains nearly stationary in the reference configuration during the early to intermediate post-bifurcation stages but spreads notably in the final phase of deformation, with a dramatic increase in the bulging amplitude. From the final bulging morphology of the cylinder in the current configuration (as shown by the inset in Figure 9), a significant discrepancy can be observed between the interface locations in the reference ( Z ≈ 10 ) and current ( z ≈ 4 ) configurations. This migration of the bulging front toward the center indicates a substantial axial contraction of the cylinder, which is driven by the localized large-amplitude radial expansion. In the entire post-bifurcation stage, the 1D model predictions agree excellently with the FE results, even for this highly compressible case ( ν = 0 . 05 ), which verifies the reliability of the model and allows efficient parametric studies to be conducted in place of computationally expensive 3D simulations.

OPEN IN VIEWER

Figure 9.

The 1D asymptotic solutions (solid curves) and the corresponding 3D FE results (scattered dots) for the radial stretch a ( Z ) . Profiles are shown for four distinct necking stages, corresponding to the states marked in Figure 8(a), with N = 2 . 35 , μ s = 0 . 2 , ν s = 0 . 2 , and a highly compressible bulk ν = 0 . 05 . The inset cylinder illustrates the final non-homogeneous bulging morphology in the current configuration.

Incidentally, we also investigate the relationship between bulk compressibility and the surface elasticity μ s . Here is a brief summary of our observations and conclusions. The bifurcation threshold γ I is inversely proportional to μ s . Increasing μ s promotes the occurrence while preventing the propagation of bulging instability. The critical γ I is sensitive to μ s when the Poisson’s ratio is sufficiently small. There is a μ s max above which bulging instability does not occur, and μ s max is inversely proportional to the Poisson’s ratio.

4.3. Compressibility-regulated guidelines for PR bulging

Under loading scenario II, the initiation and subsequent evolution of elasto-capillary bulging are governed by a synergistic interplay between the bulk compressibility ν and the axial preload N. The observed increase of N min ⁡ with ν reveals that bulk compressibility destabilizes the system. The critical surface tension γ I exhibits extreme sensitivity to both ν and N. The post-bifurcation analysis further reveals that the lack of volume conservation in compressible materials facilitates low-threshold triggering but simultaneously limits the bulging amplitude. In compressible materials, the PR bulging is rather sensitive to the fixed axial force N. In contrast, near-incompressible materials maximize the geometric non-linearity of the bulge due to the strict necessity of maintaining volume, thus achieving high-amplitude morphological responses.

These findings demonstrate that the choice of Poisson’s ratio ν facilitates a strategic trade-off between triggering sensitivity and deformation magnitude for regulating the PR bulging under loading scenario II. For applications requiring soft triggering or high sensitivity to external stimuli, materials with high compressibility ( ν ≪ 0 . 5 ) are optimal. This provides a wider operational window for triggering morphological transitions with minimal energetic cost, albeit at the expense of a limited bulging amplitude. For systems where large-scale geometric changes or significant tactile feedback is desired, near-incompressible materials ( ν → 0 . 5 ) are preferred. While these materials require a higher activation force (N), they are capable of achieving larger radial expansion amplitudes.

5.1. Solution strategy

The bifurcation condition ω ( λ c ) = 0 can be derived via a weakly nonlinear analysis by using the 1D model (see section 1). To trace the entire post-bifurcation behavior, the 1D model is discretized using finite difference methods and solved iteratively, with the analytical weakly nonlinear solution providing a primary initial guess.

Specifically, we divide the interval [ 0 , L ] into n equal parts to obtain [ Z 0 , Z 1 ] , [ Z 1 , Z 2 ] ,…, [ Z n − 1 , Z n ] , where Z j = jh , j = 1 , 2 , . . . , n with h = L ∕ n . For convenience, we define λ ( Z j ) = λ j . Thus, the fixed-length constraints (30) can be expressed as follows

1 2 λ 0 + ∑ j = 1 n − 1 λ j + 1 2 λ n − λ c L h = 0 .

(31)

The finite difference scheme for the 1D ordinary differential equation (15) is given by

A 2 μ 2 [ F ( λ j ) − N ] − 1 2 ∂B ∂λ ( λ j ) ( λ j + 1 − λ j − 1 2 h ) 2 − B ( λ j ) λ j + 1 − 2 λ j + λ j − 1 h 2 = 0 ,

(32)

where N = F ( λ ∞ ) and j = 0 , 1 , … , n .

Moreover, the finite difference discretization of the boundary conditions at the ends is also required. At Z = 0 , the symmetric deformation assumption λ ′ ( 0 ) = 0 yields

λ 1 − λ − 1 2 h = 0 ,

(33)

where λ − 1 = λ ( − h ) . At the end Z = L , the following soft asymptotic condition can be employed

λ ′ ( L ) + k ( λ ( L ) − λ ∞ ) = 0 ,

(34)

according to the exponential decaying solution (43) when the cylinder is effectively infinite long [31]. Thus, we have

λ n + 1 − λ n − 1 2 h + k ( λ n − λ ∞ ) = 0 .

(35)

Solving (33) and (35) for λ − 1 and λ n + 1 , respectively, and substituting them into the finite difference equation (15) for j = 0 , n yields the corresponding discrete boundary conditions

A 2 μ 2 [ F ( λ 0 ) − N ] − 2 B ( λ 0 ) λ 1 − λ 0 h 2 = 0 , A 2 μ 2 [ F ( λ n ) − N ] − 1 2 B ( λ n ) k 2 ( λ n − λ ∞ ) 2 − 2 B ( λ n ) λ n − 1 − λ n − hk ( λ n − λ ∞ ) h 2 = 0 .

(36)

Since strain softening may occur in this loading scenario, γ does not vary monotonically. Consequently, we treat λ 0 as the loading parameter rather than γ. The n + 2 equations (32), (36), and (31) constitute an algebraic system with n + 2 unknowns: λ 1 , λ 2 , … , λ n , λ ∞ , and γ. This discrete system is solved iteratively. First, we select a loading parameter λ 0 near the critical bifurcation point as the first loading step; an initial guess for the unknowns is then constructed using the analytical solution of the weakly nonlinear equation (see section 1). This allows us to find a numerical solution for the 1D model in the near-critical region. Following the computation of a near-critical numerical solution, subsequent steps use the preceding solution as an initial guess. This ensures the convergence of the finite difference system throughout the post-bifurcation regime.

5.2. The compressibility effect on bifurcation condition

Using the bifurcation condition (41) derived in section 1, we establish the relationship between the critical surface tension γ cr and the applied stretch λ c , as shown in Figure 10. A minimum critical surface tension, denoted γ m , occurs at a specific stretch λ m . For stretches below λ m , necking takes place, while bulging occurs above λ m . This trend aligns with earlier results for the incompressible case reported in [10] and [27]. Furthermore, the critical surface tension rises sharply as λ c moves away from λ m .

OPEN IN VIEWER

Figure 10.

For the case of μ s = ν s = 0 . 2 , plot of γ cr varying with λ c with ν = 0 . 05 , 0.15, 0.25, 0.35, 0.45, respectively. In the close-up view, the minimum point of each curve is marked by a round dot of the same color. Varying ν in the full range, the γ m − λ m response is marked by the black dashed curve. One can see that the black dashed curve intersects with the five round dots.

Mathematically, λ m and γ m can be obtained by solving

ω ( λ m ) = 0 , dω ( λ m ) dλ = ( ∂ω ∂a da dλ + ∂ω ∂λ ) | λ = λ m = 0 .

(37)

By varying the parameters ν, μ s , and ν s separately, we obtain the dependence of the necking-bulging transition stretch λ m on each parameter, as shown in Figure 11. It is observed that as the Poisson’s ratio ν increases, λ m grows nonlinearly from 0.5 to 1.2, whereas γ m first decreases and then increases. These trends indicate that the material’s compressibility modulates both the critical stretch required for instability and the associated surface energy needed to trigger either necking or bulging. On the contrary, λ m decreases gradually with increasing surface elasticity μ s and exhibits only a weak, nearly constant dependence on the surface Poisson’s ratio ν s . Hence, we can conclude that the competition between necking and bulging instabilities is most sensitive to variations in the bulk Poisson’s ratio ν, followed by the surface elasticity μ s , whereas the surface Poisson’s ratio ν s exerts only a marginal influence on the instability boundary.

OPEN IN VIEWER

Figure 11.

The dependence of the necking-bulging transition stretch λ m on (a) the Poisson’s ratio ν, (b) the surface elasticity μ s , and (c) the surface Poisson’s ratio ν s .

5.3. The compressibility effect on post-bifurcation behaviors

In the following, we explore the influence of bulk compressibility and geometric size on post-bifurcation behaviors, specifically necking when λ c < λ m and bulging when λ c > λ m . The Poisson’s ratio ν is taken as 0.25, 0.15, and 0.05, respectively; the length-radius ratio L ∕ A is set to 20, 40, and 100, respectively, with μ s = ν s = 0 . 2 held fixed. The corresponding critical stretches λ m for each ν are 0.6717, 0.6077, and 0.5684. Accordingly, we analyze the necking and bulging behaviors separately by taking λ c = 0 . 5 and λ c = 0 . 7 .

5.3.1. Post-bifurcation necking behavior

For λ 0 = 0 . 5 , the deformation proceeds as follows. Starting with γ = 0 , the cylinder is first stretched to λ 0 = 0 . 5 by a dimensionless axial force N of approximately 2.8295. As γ is gradually increased, the axial force N rises accordingly to maintain the prescribed stretch. The cylinder retains its uniform shape until γ reaches the critical value γ cr , beyond which necking localizes and develops.

Figure 12 plots γ ∕ γ cr against λ 0 , along with the evolution of the necking amplitude as γ ∕ γ cr increases. First, slender cylinders exhibit a more pronounced non-monotonicity in the γ ∕ γ cr − λ 0 response, indicating that significant strain softening occurs only when L ∕ A is sufficiently large. At advanced stages of necking, the curves for different L ∕ A nearly coincide, suggesting that the effect of L ∕ A diminishes once the bifurcation has propagated far enough. These trends align with earlier observations for incompressible materials [27]. Second, compressibility distinctly modulates both the size effect and the necking severity. Horizontal comparison of Figure 12 shows that necking in more compressible materials is less sensitive to L ∕ A . That means the size effect weakens as ν decreases. Moreover, at a fixed γ ∕ γ cr level, necking becomes more pronounced when the cylinder is more compressible (i.e., when ν is smaller). These results reveal that the length-radius ratio L ∕ A primarily influences strain-softening behavior, whereas bulk compressibility governs the size effect and the degree of necking.

OPEN IN VIEWER

Figure 12.

Post-bifurcation necking responses under fixed average stretch λ c = 0 . 5 for various aspect ratios L ∕ A and Poisson’s ratios ν. (a–c) γ ∕ γ cr - λ 0 response. (d–f) Evolution of the necking amplitude a L − a 0 vs. γ ∕ γ cr . The 1D model results (solid curves) are shown for L ∕ A = 20 , 40 , and 100. In the highly compressible case (c and f) where ν = 0 . 05 , FE results (diamonds) for L ∕ A = 40 are included, showing excellent agreement with the 1D model. The four open circles marked on the curve in (c) represent distinct necking stages, for which the corresponding localized profiles a ( Z ) are displayed in Figure 13. Other parameters are set as μ s = 0 . 2 and ν s = 0 . 2 .

The deformation morphologies at the four stages marked in Figure 12(c) are shown in Figure 13. A comparison between the 1D solution and the FE simulation in Figure 13(a) shows excellent agreement, confirming the validity of the 1D asymptotic results. As both ends are fixed, the necked region continuously thins, while the bulging part becomes increasingly prominent. In the current configuration, the necking front is observed to propagate continuously toward both ends of the cylinder as the deformation progresses.

OPEN IN VIEWER

Figure 13.

The 1D asymptotic solutions (solid curves) and the corresponding 3D FE results (scattered dots) for the radial stretch a ( Z ) . Profiles are shown for four distinct necking stages, corresponding to the states marked in Figure 12(c), with L ∕ A = 40 , μ s = 0 . 2 , ν s = 0 . 2 , and a highly compressible bulk ν = 0 . 05 . The inset cylinder illustrates the final non-homogeneous necking morphology in the current configuration.

5.3.2. Post-bifurcation bulging behavior

For λ c = 0 . 7 , the post-bifurcation bulging behavior is demonstrated in Figure 14. Initially, with surface tension set to zero, the cylinder is uniformly stretched to λ c = 0 . 7 under an axial force of approximately 1.3250 (given μ s = ν s = 0 . 2 ). Localized bulging initiates once γ reaches the critical value γ cr . As the bifurcation proceeds, the stretch λ 0 continues to decrease while the amplitude a 0 − a L steadily increases.

OPEN IN VIEWER

Figure 14.

Post-bifurcation bulging responses under fixed average stretch λ c = 0 . 7 for various aspect ratios L ∕ A and Poisson’s ratios ν. (a–c) γ ∕ γ cr - λ 0 response. (d–f) Evolution of the necking amplitude a 0 − a L vs. γ ∕ γ cr . The 1D model results (solid curves) are shown for L ∕ A = 20 , 40 , and 100. In the highly compressible case (c, f) where ν = 0 . 05 , FE results (diamonds) for L ∕ A = 40 are included, showing excellent agreement with the 1D model. The four open circles marked on the curve in (c) represent distinct bulging stages, for which the corresponding localized solutions a ( Z ) are displayed in Figure 15.

The influence of the length radius L ∕ A on the bulging response is similar to that observed in the previous necking situation. In contrast, the effect of Poisson’s ratio ν is opposite to that found for necking. A horizontal comparison of Figure 14 indicates that greater compressibility (smaller ν) not only enhances the sensitivity of the bifurcation outcome to variations in L ∕ A , but also reduces the necking degree.

In Figure 15, we present the bulging profiles a ( Z ) with L ∕ A = 40 , ν = 0 . 05 , together with corresponding FE simulation results (marked in Figure 14(c)). It is seen that the bulging region continues to thicken and propagate axially, while the remaining part of the cylinder undergoes further thinning to accommodate the global length restriction. This leads to an increasingly sharp transition between the expanded and constricted zones.

OPEN IN VIEWER

Figure 15.

The 1D asymptotic solutions (solid curves) and the corresponding 3D FE results (scattered dots) for the radial stretch a ( Z ) . Profiles are shown for four distinct bulging stages, corresponding to the states marked in Figure 14(c), with L ∕ A = 40 , μ s = 0 . 2 , ν s = 0 . 2 , and a highly compressible bulk ν = 0 . 05 . The inset cylinder illustrates the final non-homogeneous bulging morphology in the current configuration.

5.4. Compressibility-regulated guidelines for PR necking/bulging under fixed-end constraints

Under fixed-end conditions, the bulk compressibility (quantified by ν) and the length-radius ratio L ∕ A exert distinct influences on the elasto-capillary instability. The bulk compressibility primarily governs the onset threshold and localization intensity of necking or bulging (determined by the fixed average stretch λ c ), whereas L ∕ A dominates the evolution of strain-softening behavior. As ν decreases, necking becomes more pronounced while its sensitivity to geometric size diminishes; in contrast, during bulging, higher compressibility amplifies the size effect. These results reveal a powerful decoupling between geometric and material parameters, which can be exploited to independently tailor the onset mode, the strain-softening rate, and the instability intensity in compressible elastic cylinders.

Based on these findings, we propose the following guidelines: