## Abstract

In this paper, we study the corner-like formation in a slender nonlinearly elastic cylinder due to compression. Mathematically, this is a very challenging problem: one needs to study the bifurcation of the nonlinear field equations (complicated nonlinear partial differential equations) and show that there is a bifurcation leading to a solution with a corner-like profile. We also aim to obtain the asymptotic expression for this post-bifurcation solution. As far as we know, there is no available analytical method to obtain the post-bifurcation solution of nonlinear partial differential equation(s). Here we use a novel approach to tackle the present problem. Through a method of coupled series–asymptotic expansions, we manage to derive the normal form equation of the original nonlinear field equations, which can be written as a singular ordinary differential equation (ODE) system (the vector field has a singularity at one point). With welding end conditions, the problem is reduced to study the boundary-value problem of a singular ODE system. It seems that such a boundary-value problem of a singular ODE system was not formulated or studied before in the context of the present problem. With the help of phase planes, we manage to solve such a boundary-value problem. It turns out that there is a vertical singular line in phase planes, which causes the bifurcation to a corner-like solution. The expression for this post-bifurcation solution is also obtained. From the analytical results obtained, we reveal that the coupling effect of the material nonlinearity and geometrical size is the physical mechanism that causes the formation of a corner-like profile.

## 1. Introduction

Stabilities and instabilities are important topics in nonlinear elasticity and structures. Some often-appeared instability phenomena, such as the bending instability, barrelling instability, shear and necking instabilities, lateral buckling of a beam and buckling of a whirling rod, were considered by Antman (1995). Many authors have studied various stability and instability problems. For example, in Fu (2001), perturbation methods were used to carry out nonlinear stability analysis on some difficult problems such as the necking of an elastic plate under stretching. Here we study analytically the instability leading to a corner-like formation by investigating the uniaxial compression of a solid circular cylinder under welding end conditions.

Barrelling and buckling instabilities of compression of a solid circular cylinder or two-dimensional rectangular columns have been studied extensively, for example Simpson & Spector (1984*a*,*b*) and Davies (1989, 1991). The end condition used in these literatures was that the surface traction had no tangential component (which meant the end was greased or lubricated) and the boundary condition of the lateral surface was traction free. The analysis mentioned above was all linear analysis. Fu (2001) gave a list of contributions to linear stability analysis. Healey & Montes-Pizarro (2003) presented rigorous local and global bifurcation results for a compression of a three-dimensional nonlinear cylinder leading to a barrelling state by using the generalized degree designed in Healey & Simpson (1998) to overcome the specific difficulties of three-dimensional nonlinear elasticity. Their results are the first global bifurcation results for a problem from three-dimensional nonlinear elastostatics not governed by ordinary differential equations. With lubricated end conditions, they posed the problem in the space of solutions being twice Hölder continuous and proved the existence of non-trivial solutions bifurcating from trivial homogeneous solutions in such a space. With their global bifurcation results, it is possible to get the numerical post-bifurcation solutions by finite-element methods.

In this paper, we also study the uniaxial compression of a nonlinearly elastic cylinder with constant circular cross-section and focus on the corner-like instability. We consider the case that the ends of the cylinder are welded to rigid bodies. It seems that there is no study on stability and instability on a cylinder with such end conditions. The present work is partially motivated by an open problem mentioned in a review article by Beatty (1987), in which an experiment conducted by Willis (1948) was described. In the experiment, a sufficiently short and thick-walled cylindrical tube was compressed between lubricated rigid, parallel end plates (figure 1). Initially, both the inside and outside begin to ‘bulge outward to form convex, barrel-shaped surfaces at fairly small strain. But at approximately 15–20% compression, these convex surfaces become straight again near the platen ends while remaining convex in the major central portion of the sample. At a certain critical compressive load intensity, the inside wall reverses its curvature (Payne & Scott 1960) and collapses to form a cusp-shaped inside surface directed towards the outside boundary in the central plane, while the outside wall continues to bulge as usual (Beatty 1987). This phenomenon was called the Willis instability phenomenon. Beatty (1987) used ‘cusp’ to describe this instability. In fact, this cusp observed in the experiment can be regarded as a sharp crest with a very large second-order derivative. So, here we refer it as a ‘corner-like’ crest, that is, in this paper, we define a corner-like profile as one that is smooth but has a very large second-order derivative at one interior point. Actually, corner-like instabilities are widespread and can appear on many occasions such as compression of a block of sponge, compression of a soft drink can, compression of a cigarette filter, bending of a plastic tube, etc. There are many studies on various instabilities; however, it seems that there is no analytical study on the important corner-like instability. This could be because the corner-like instability is more difficult to study.1

The difficulty lies in that unlike most other instabilities the ‘post-buckling’ profile is nearly ‘non-smooth’ in a corner-like instability, in terms that the second-order derivative is very large at the crest. There are two important issues concerned with the corner-like instability. One is what is the physical mechanism that can cause the corner-like formation and the other is what kind of mathematics can be used to capture the sharp crest with a very large second-order derivative.

In this paper, we resolve these issues for an elastic solid cylinder subject to compression. The strain energy function of the cylinder is supposed to be given. The starting point is the exact three-dimensional field equations. Also, we suppose that the two ends are welded to rigid bodies, and as a result the boundary conditions are different from those considered by other works mentioned in the second paragraph. In addition, we aim to obtain the asymptotic solution for the post-bifurcation profile. To get the post-bifurcation solutions of nonlinear partial differential equations (PDEs) is a very difficult problem, and as far as we know, there is no analytical method available. Owing to the complexity of the governing equations, we use a formal approach to tackle the difficult problem of bifurcations of PDEs. Such a methodology has been introduced in Dai & Huo (2002), Dai & Fan (2004) and Dai & Cai (2006) to study other nonlinear problems. The basic idea is to use coupled series–asymptotic expansions to deduce the asymptotic normal form equation(s) for the nonlinear PDEs.

After proper scaling and a change of variables to the nonlinear system, a small variable and two small parameters are identified. A series expansion in this small variable reduces our problem. Further, we use the smallness of the parameter that characterizes the weak nonlinearity to deduce a nonlinear ordinary differential equation (ODE), which governs the leading-order axial strain, by using the regular perturbation method. Since all the other physical quantities can be expressed in terms of this axial strain, we can obtain the asymptotic solution for the three-dimensional displacement field.

It turns out that this nonlinear ODE, which is called the normal form equation of the original coupled nonlinear PDEs, can be written as a singular first-order system (the vector field is singular at one point). As a result, we have a boundary-value problem of a singular first-order ODE system at hand. While it was found that singular ODE systems of this type in an infinite domain can arise for travelling-wave solutions of certain nonlinear dispersive equations (Dai 1998, 2001; Dai & Huo 2000; Dai *et al.* 2004), it seems that the boundary-value problem of such a singular ODE system has not been formulated or studied before in the context of the present problem. It happens that the boundary conditions can also influence the bifurcations greatly.

We analyse this system with the help of phase planes. It turns out that there exists a vertical singular line in each phase plane and its relative position to equilibrium points plays an important role. The presence of this singular line, under a critical loading, can cause the appearance of a corner-like solution with a large second-order derivative. The asymptotic solution expressions are also deduced. From the analytical results obtained, we reveal that the coupling effect of the material nonlinearity and geometrical size is the mechanism that causes the formation of a corner-like crest.

The outline of this paper is as follows. In §2, we formulate the field equations and the boundary conditions for a hyperelastic material with a specific strain energy function together with jump conditions. In §§3 and 4, we use a novel method involving coupled series–asymptotic expansions to deduce the asymptotic normal form equation for the nonlinear field equations, which is a nonlinear ODE. In §5, by considering the energy and using the variational principle, the same ODE is obtained, which verifies our previous approach. In §6, we consider the bifurcation of the solutions under the welding end conditions. Finally, some conclusions are drawn at the end.

## 2. Field equations

We consider the axisymmetric deformations of a slender cylinder subject to an axial static force at two ends. The radius of the cylinder is *a* and the length is *l*. The geometry of this problem is described in figure 2. We take a cylindrical polar coordinate system and denote (*R*, *Θ*, *Z*) and (*r*, *θ*, *z*) as the coordinates of a material point in the reference and current configurations, respectively. The finite radial and axial displacements can be written as(2.1)We introduce the orthonormal bases associated with the cylindrical coordinates and denote these by *E*_{R}, *E*_{Θ}, *E*_{Z} and *e*_{r}, *e*_{θ}, *e*_{z} in the reference and current placements, respectively. The deformation gradient tensor is given in these orthonormal bases by(2.2)We consider the material nonlinearity up to the second order and suppose that the cylinder is composed of an isotropic compressible homogeneous material. For convenience, we consider the hyperelastic Murnaghan material, whose strain energy function *Φ* is(2.3)where is the Green strain tensor; *λ* and *μ* are known as the Lamé coefficients; and *ν*_{1}, *ν*_{2} and *ν*_{4} are other constitutive constants. Such a material has been studied in Dai & Fan (2004). The formula for the computation of the first Piola–Kirchhoff stress tensor containing terms up to the second-order material nonlinearity for an arbitrary strain energy function has been provided in Fu & Ogden (1999), which is given as(2.4)where =−, and are elastic moduli defined by(2.5)and is the identity tensor corresponding to a natural configuration. The expression of the component *S*_{zZ} corresponding to the strain energy function for the Murnaghan material is given as(2.6)and the expressions for other components are omitted for brevity. satisfies the following field equations:(2.7)(2.8)Substituting the components of into the above equations, we have two lengthy equations that are given by(2.9)(2.10)

In addition to the above two equations, we assume that the lateral surface of the cylinder is traction free. The stress components *S*_{rR} and *S*_{zR} should vanish on the lateral surface. Then we have the following boundary conditions:(2.11)

(2.12)

We consider the case that the two ends are welded to rigid plates with a compressive force applied at one end. So, we have(2.13)and(2.14)where *E* is Young's modulus so that *γ* is the dimensionless stress. (*Note*. In reality, usually one only knows the resultant, that is, is known.)

The focus is on the bifurcation of the nonlinear PDEs (2.9) and (2.10) under (2.11)–(2.14). We shall show that there is a bifurcation leading to a corner-like profile and deduce the asymptotic expression for the post-bifurcation solution.

Also, we note that for the solution obtained if its deformation gradient has a finite jump across a surface in the cylinder, this solution should satisfy the jump conditions that are associated with force balance and the continuity of the deformation, namely(2.15)must hold on the surface of discontinuity, where ^{±} and ^{±} denote its limiting values at a point on the surface of discontinuity; ** n** is its unit normal; and

**are all vectors tangent to the surface.**

*l*## 3. Transformed dimensionless equations

It is hard to deal with the complicated nonlinear PDEs (2.9) and (2.10) together with (2.11)–(2.14) directly. So, we will non-dimensionalize them to find small parameters and small variables. Then we can use the expansion method to proceed further. The approach adapted here is similar to those used in Dai & Huo (2002), Dai & Fan (2004) and Dai & Cai (2006).2 Here, we shall further consider the coupling effect of the material nonlinearity and geometrical size, and concentrate on the bifurcations of solutions in a boundary-value problem, which are not addressed in the previous studies.

First, we introduce a new set of dimensionless quantities through the following suitable scalings:(3.1)where *h* is a characteristic axial displacement that can be taken as the reduction of the distance between the ends in a uniform compression with a given *γ* (the specific value *h* does not matter and only its magnitude matters) and *ϵ* will be treated as a small parameter (i.e. we are considering a weak nonlinearity). Further, we introduce a very important change of variables as follows:(3.2)Then we can obtain dimensionless field equations and boundary conditions as follows:(3.3)(3.4)(3.5)(3.6)where , , , .

Then equations (3.3)–(3.6) compose a new system that is still very complicated and difficult to analyse directly. In the next section, by using series and asymptotic expansions, we derive the normal form equation of this complicated system in order to analyse the solution bifurcations.

## 4. Coupled series–asymptotic expansions and the normal form equation

In the §3, we find that does not appear explicitly or implicitly in the form of in (3.3)–(3.6). The importance of (3.2) is that the unknowns become functions of *s* and *x* and we only need to work with *s*; there is no need to consider the odd-power terms of . Here we consider a slender cylinder and the variable *s* is small. We can see that the whole problem depends on two small parameters *ϵ* and *δ*, one small variable *s* and the other variable *x*. Assume that *w* and *v* have the following Taylor expansions in the neighbourhood of *s*=0:(4.1)

(4.2)

A similar series expansion in Boström (2000),3 was used to derive a hierarchy of wave equations in linearly elastic rods. Here, we shall use both series and asymptotic expansions to tackle nonlinear problems. Substituting (4.1) and (4.2) into (3.3) and (3.4), the two left hands become series in the small variable *s*. All the coefficients should be zero, and, as a result, we have two sets of infinitely many equations. It is important to work with a closed system with finite equations, which can be achieved asymptotically by taking three equations from these two sets and further using the traction-free boundary conditions. For that purpose, first from the coefficient of *s*^{0} coming from (3.3), we obtain(4.3)Setting the coefficients of *s*^{0} and *s*^{1} coming from (3.4) to be zero, we obtain(4.4)(4.5)The above three equations contain five unknowns: *w*_{0}, *w*_{1}, *w*_{2}, *v*_{0} and *v*_{1}. To further obtain two equations, we substitute (4.1) and (4.2) into the boundary conditions (3.5) and (3.6) to get(4.6)(4.7)where *α*_{i} (*i*=1, 2, …, 19) are constants related to material constants (whose expressions are given in appendix A) and *O*(*ϵδ*^{2}, *δ*^{2}) include such terms as *v*_{0}, *v*_{1}, *v*_{2}, … and *w*_{0}, *w*_{1}, *w*_{2}, … and their derivatives with respect to *x* multiplied by *ϵδ*^{2}, *δ*^{2} or higher-order terms.

Now the governing equations (3.3) and (3.4) are changed into a one-dimensional system of differential equations (4.3)–(4.7) for the five unknowns *w*_{0}, *w*_{1}, *w*_{2}, *v*_{0} and *v*_{1} if we neglect the higher-order terms of *O*(*ϵδ*^{2}, *δ*^{2}). Here, we only look for solutions valid up to *O*(1), and it is permissible to neglect *O*(*ϵδ*^{2}, *δ*^{2}) terms. It should be pointed out that, as it is expected that this is a singular perturbation problem, one needs to keep *O*(1), *O*(*ϵ*), *O*(*δ*) and O(*ϵδ*) terms for obtaining *O*(1)-valid results. By using a regular perturbation method, we can express *w*_{1} by *w*_{0} and *v*_{0} from (4.4), then express *v*_{1} also by *w*_{0} and *v*_{0} from (4.3). Similarly, we substitute *w*_{1} and *v*_{1} into (4.5) and obtain the expression of *w*_{2} by the perturbation method. The results are given as(4.8)(4.9)(4.10)Substituting *w*_{1}, *v*_{1} and *w*_{2} into (4.6) and (4.7) and omitting the higher order terms yields the following two equations with only two unknowns *w*_{0} and *v*_{0}:(4.11)(4.12)By using the above two equations, we can express *v*_{0} in terms of *w*_{0} and its derivatives with respect to *x*. From (4.11), we have(4.13)From (4.12) and (4.11), we have(4.14)and(4.15)Substituting into (4.14) and (4.15), and then substituting (4.14) into *δ* terms of (4.13) and substituting (4.15) into *ϵ* terms of (4.13), we obtain(4.16)In these manipulations, the *O*(ϵ^{2}, *ϵδ*^{2}, *δ*^{2}) terms are omitted, which is permissible for obtaining *O*(1)-correct results. Now, from (4.12)–(4.16), we obtain(4.17)This is a single fourth-order nonlinear ordinary differential equation for the single unknown *w*_{0x}. It turns out that it can be integrated once, and as a result we obtain(4.18)where , and *A* is the integration constant.

It is important to identify the physical meaning of *A* since it is a physical problem and we would like to consider the bifurcations as physical parameters vary, not simply the existence of bifurcations. For that purpose, we consider the resultant force *T* acting on the cross-section that is planar and perpendicular to the cylinder axis in the reference configuration, which is given by(4.19)After tedious calculations, it is found that(4.20)We take the derivative with respect to *x* for (4.20). As (d*T*/d*x*)=0 for a static problem, we find that this yields exactly (4.12). Replacing *v*_{0} in (4.20) as above, we obtain(4.21)where is Young's modulus. Comparing (4.21) and (4.18), we have . Thus,(4.22)If we retain the original dimensional variable and let , then we have(4.23)where is the dimensionless engineering stress. We call (4.23) the normal form equation of the field equations (nonlinear PDEs) (2.9) and (2.10) together with the nonlinear boundary conditions (2.11) and (2.12) with a given axial resultant, since all terms needed to yield *O*(1)-valid results are kept in (4.23) and once *V* is found the three-dimensional displacement field (thus, strain and stress fields) can be found.

In the following section, we also use the energy approach to derive the same equation, which then verifies our aforementioned method.

## 5. The Euler–Lagrange equation

The normal form equation (4.23) is derived in the previous section based on the force balance equations (2.7) and (2.8). Now, we consider the energy and show that the same equation can be derived through the variational principle. The calculations are tremendously complicated in the present method (we use Mathematica for all the symbolic computations), and if the two derivations of the normal form equation (4.23) agree with each other, it also supports the correctness of our calculations. From (2.3), we can get the expression of the strain energy up to the third-order nonlinearity (which implies that the stress components are up to the second-order nonlinearity):(5.1)where *k*_{1}=(1/2)*λ*+(1/2)*ν*_{1}, *k*_{2}=(1/2)*λ*+*ν*_{1}+*ν*_{2}, *k*_{3}=(1/2)*λ*+*μ*+(1/2)*ν*_{1}+(1/4)*ν*_{4}, *k*_{4}=(1/2)*λ*+*μ*+*ν*_{1}+(1/3)*ν*_{2}+(1/3)*ν*_{4}, *k*_{5}=*μ*+*ν*_{1}+(1/2)*ν*_{4}.

The strain energy per unit length is given by(5.2)By using the expressions obtained in §§3 and 4, we can get the average strain energy over a cross-section in terms of *w*_{0} and *v*_{0}, which is given by(5.3)where *ξ*_{i}(*i*=1, …, 12) are constants related to the constitutive constants and their expressions can be found in appendix A.

By further using (4.11) and (4.12), we can reduce the above equation as(5.4)where(5.5)Without loss of generality, we take the length of the cylinder to be 1 and the total potential energy is then given by(5.6)Further by the variational principle, we have the following Euler–Lagrange equation:(5.7)which is(5.8)It can be seen that the above equation is just (4.23).

The expression of the total potential energy (5.6) itself is very important: when there are multiple solutions, the smallest energy criterion can be used to judge the preferred solution (configuration).

## 6. The boundary-value problem of a singular ODE system

Now, we conduct a detailed analysis on the normal form equation (4.23) together with welding end conditions. For that purpose, we rewrite it as a first-order system(6.1)

Since we only look for *O*(1)-correct results, it is sufficient to use *O*(1)-correct end boundary conditions. To the leading order, the welded ends imply that(6.2)From (4.16), also replacing the higher-order derivative terms by using (4.23), then letting *V*=*ϵw*_{0x}, we obtain(6.3)This equation has a unique negative root *Δ*. Then we have the following equivalent welding end conditions:(6.4)where(6.5)We note that *Δ* depends on *γ*.

We note that the vector field of the system (6.1) has a denominator term. If the material constant *D*_{2}>0, then in the compression case the vector field has a singular point at *V*=−(1/8*D*_{2}). Thus, this is a singular ODE system. So, (6.1) together with (6.4) forms a boundary-value problem of a singular system. The singular ODE system of this kind is rarely studied, and, as far as we know, was only analysed for travelling-wave solutions (in an infinite domain without boundary conditions or initial conditions) of certain nonlinear dispersive-wave equations (Dai 1998, 2001; Dai & Huo 2000; Dai *et al.* 2004).4 It appears that a boundary-value problem of a singular ODE system of the type has not been formulated or studied before in the context of the present problem. As we shall see the presence of this singularity can affect the bifurcation greatly. Also, from our derivations, it can be seen that this denominator term comes from the *δϵ*^{2} term in (4.18), which represents the coupling effect of the material nonlinearity and geometrical size. If we can show mathematically that this term leads to the bifurcation to a corner-like solution, then we have identified the mechanism of the corner-like instability: it is due to the interaction between the material nonlinearity and geometrical size, which causes the formation of a corner-like profile.

This singular ODE system contains a geometrical parameter *a*^{2}, two material constants *D*_{1} and *D*_{2} and a loading parameter *γ*, and also the boundary condition is related to material constants and *γ*. As a result, when the material constants are in different domains, as the engineering stress *γ* varies, there are a variety of bifurcation phenomena. Here, we concentrate on the case 4*D*_{2}>*D*_{1}>0, for which it is sufficient to illustrate the corner-like formation.

The equilibrium points of the system are given by(6.6)We consider that the loading *γ* (≤0) is in the interval (−1/4*D*_{2},0] such that (6.6)_{2} has two roots, say, denoted by *V*_{1} and *V*_{2} with *V*_{1}>*V*_{2}. In the *V*–*y* phase planes, there is also a vertical singular line *V*=−(1/8*D*_{2}), which, for 4*D*_{2}>*D*_{1}, is always located to the right of the equilibrium point (*V*_{2}, 0). It is easy to show that (*V*_{2}, 0) is always a saddle point. However, for the equilibrium point (*V*_{1},0), relative to the position of the vertical singular line, it can be a saddle point, an equilibrium point on the singular line and a centre point. As *γ* varies, there are three different phase planes, which are shown in figure 3, where *γ*_{c} will be defined later and *γ*_{p} is the value calculated from (6.6)_{2} for *V*=(1/8*D*_{2}), i.e.(6.7)

For the above phase planes, we have chosen *a*=0.25; *D*_{1}=2.4; *D*_{2}=0.65; the Poisson ratio *ν*=0.1. Then, correspondingly *γ*_{p}=−0.10355. We note that for a trajectory in a phase plane to be a solution of the boundary-value problem of (6.1) and (6.4) a necessary and sufficient condition is that it contacts the vertical line *V*=*Δ* twice and its *Z*-interval is exactly equal to 1. Next, we discuss the solution(s) in each phase plane separately.

### (a) *Case* (*a*) *γ*_{p}<*γ*<0

In this case, there is a unique solution denoted by trajectory 1 in the corresponding phase plane. Trajectory 1 crosses the *V*-axis at (*V*_{0}, 0). By referring to Gradshteyn & Ryzhik (1980), the solution expression is given by(6.8)where *Π* is the elliptic integral of the third kind and *F* is the elliptic integral of the second kind and(6.9)(6.10)(6.11)As *V*=*Δ* at *Z*=0, *V*_{0} is determined by(6.12)where

### (b) *Case* (*b*) *γ*=γ_{p}

In this case, there is a trajectory tangent to the singular line in the phase plane. There is a unique solution denoted by trajectory 1 in the phase plane. Trajectory 1 contacts the *V*-axis at (*V*_{0}, 0). The solution expression is given by(6.13)where *E*_{1} and *E*_{2} are the same as (6.9) and (6.10), respectively, and(6.14)As *V*=*Δ* at *Z*=0, *V*_{0} is determined by(6.15)where .

### (c) *Case* (*c*) *γ*_{c}<*γ*<*γ*_{p}

In the phase plane for this case and the following case, there is a trajectory crossing the singular line at points B and C and crossing the vertical line *V*=*Δ* at A and D.5

There is a unique solution represented by trajectory 1. It contacts the *V*-axis at (*V*_{0}, 0). This solution has the same expression as (6.13) and (6.14) and *V*_{0} is determined by (6.15).

### (d) *Case* (*d*) *γ*<*γ*_{c}

We define *γ*_{c} to be the critical stress value such that when *γ*=*γ*_{c} trajectory 2 in figure 3*c*, which starts from A, goes to B, jumps from B to C and finally arrives at D, is also a solution. In this case, there is another solution indexed by 1 in figure 3*c*. For trajectory 1, it contacts the *V*-axis at (*V*_{0}, 0), and its solution expression is the same as (6.13) and (6.14) and *V*_{0} is determined by (6.15). We find that for this solution at *Z*=0.5 the value of *V*_{ZZ} is relatively small. Indeed, for *γ*=*γ*_{c}, *V*_{ZZ}=1.66271 at *Z*=0.5 for the previously chosen parameters (*a*=0.25, *D*_{1}=2.4, *D*_{2}=0.65 and the Poisson ratio *ν*=0.1).

The solution expression corresponding to trajectory 2 is given by(6.16)where(6.17)(6.18)As *V*=*Δ* at *Z*=0, we have(6.19)which determines the value *γ*_{c}. For the previously chosen parameters, we find *γ*_{c}=0.103666.

Alternatively, (6.16) can be rewritten as(6.20)Obviously, this represents a non-smooth solution whose first-order derivative at *Z*=0.5 does not exist. This implies that the solution profile has a corner at *Z*=0.5. Since both the left and right derivatives exist, one can also view that at the corner the second-order derivative is infinite.

This non-smooth solution is a weak solution of the singular ODE system according to the definition given in Dai, Huang & Liu.6 However, for this solution, the jump conditions (2.15) cannot be satisfied exactly, thus we shall not take it as one solution of the field equations.7 Instead, we shall turn our attention on the smooth solution that has a very large second-derivative value (i.e. a corner-like solution profile) at one interior point.

### (e) *Case* (*e*) −(1/4*D*_{1})<*γ*<γ_{c}

When |*γ*| is slightly larger than |*γ*_{c}|, there are two solutions that can be indexed by 1 in the corresponding phase plane that is topologically the same as figure 3*c*. Both of them can be expressed by (6.13) and (6.14). Denote the crossing points with the *V*-axis of the two solution trajectories by and (*V*_{0}, 0), respectively, and the values of and *V*_{0} can be determined by (6.15). The difference between the two trajectories is that for one of them the point is very close to the singular line, which represents a sharp crest profile with a large second-order derivative *V*_{ZZ}.

To get the solution profile, we still use the previously chosen parameters. For , we find that *V*_{0}=−0.190148 and , which is very close to the singular line , and the corresponding total potential energy values for the two solutions are, respectively,(6.21)We can see that the trajectory 1 with has a smaller energy value and is thus a preferred solution. For this solution, we find that *V*_{ZZ}=291.11 at *Z*=0.5, indeed a large value. This solution curve corresponding to this trajectory is shown in figure 4*a*. Once *V* is found, the radial displacement *U*(*R*, *Z*) can be calculated by using (3.1)_{2}, (3.2)_{2}, (4.2), (4.9) and (4.16). The relation between the reference coordinate *Z* and current coordinate *z* is *z*=*Z*+*W*. While *W* can be calculated from (3.1)_{1}, (3.2)_{1}, (4.1), (4.8) and (4.10). Thus, a relation between *U* and the current coordinate *z* can be found for a given *R*. The surface profile of the cylinder in the current configuration is shown in figure 4*b*.

To summarize, from the results obtained in this section, we see that as the applied load increases, crossing *γ*_{c}, a bifurcation takes place: a solution with a flat profile around *Z*=0.5 (*V*_{Z}=0 and small *V*_{ZZ} value) jumps to another solution with a corner-like profile around *Z*=0.5 (in terms of a large *V*_{ZZ} value at *Z*=0.5). Obviously, this large second-order derivative value is caused by the presence of the vertical singular line. Thus, the coupling effect of the material nonlinearity and geometrical size is the mechanism that causes the corner-like formation.

## 7. Conclusions

It is a very difficult problem to study the bifurcation(s) of nonlinear partial differential equations and to obtain the post-bifurcation solution(s) analytically, and few results are available. Here, for the compression of a nonlinearly elastic cylinder, we manage to solve such a problem through a novel approach. By using coupled series–asymptotic expansions, we deduce the normal form equation of the original nonlinear field equations. It turns out that this normal form equation can be written as a singular ODE system. By a phase plane analysis, we study the bifurcation of this singular system under welding end conditions. To the authors' best knowledge, a boundary-value problem of a singular ODE system of this kind has not been formulated or studied before in the context of the present problem. Our results show that the vertical singular line in phase planes (which arises due to the singularity of the ODE system) plays an important role. In particular, its presence leads to the bifurcation to a corner-like solution with a very large second-order derivative value at one point. The expression of this post-bifurcation solution is also obtained. The mathematical results obtained also reveal that the coupling effect of the material nonlinearity and geometrical size is the physical mechanism that causes the corner-like formation in the cylinder. Some primary results (without details) with different boundary conditions were also announced in a short note (Dai & Wang 2007).

Finally, we point out that the method of the normal form equation adopted in this paper can be used to investigate other instability phenomena in thin/slender structures composed of nonlinearly elastic materials. Various problems will be studied by this approach in the future.

## Acknowledgements

The work described in this paper was fully supported by a grant from the Research Grants Council of HKSAR, China (project no. CityU 100906).

## Footnotes

↵In the work of Healey & Montes-Pizarro (2003), the post-bifurcation profile has corners at the two ends.

↵In Dai & Huo (2002), by the scaling and perturbation method, the authors derived asymptotically approximate model equations for nonlinear dispersive waves in incompressible elastic rods, which are consistent with the lateral boundary conditions. Then they compared the dispersion relations of the model equations with that obtained by using the Navier–Bernoulli hypothesis and the exact dispersion relation from the three-dimensional theory. It followed that the model equations can give a more accurate result. For solitary waves in the far field, they recovered the KdV equation. By some comparisons, it was found that the Navier–Bernoulli hypothesis is valid with small error.

↵In Dai & Fan (2004), by non-dimensionalizations of the field equations and series expansion of the axial and radial displacements and the regular perturbation method, two types of model equations are derived for weakly long waves in compressible elastic rods. By using the Navier–Bernoulli hypothesis with Love's relation, another two types of equations were obtained. After comparisons of the dispersion relation and the far-field equations for those four different model equations, it was found that the first two model equations were good approximations, while the last two may lead to small errors.

↵In Dai & Cai (2006), taking into account the radial displacement and traction-free boundary conditions, by a novel approach that combines a series expansion and an asymptotic expansion, the authors obtained a model equation, which was a nonlinear ODE governing the axial strain, from the three-dimensional field equations of a slender elastic cylinder. Some interesting analytical solutions for an infinite cylinder were found to describe the structure of a phase boundary, the nucleation process and the merge of two phase boundaries.

↵Boström assumed that the deformation is rotationally symmetric (i.e. axisymmetric) and used

*u*_{x}and*u*_{r}to denote the axial and radial displacement, respectively. The series expansions are*u*_{x}=*u*_{0}+*r*^{2}*u*_{2}+⋯,*u*_{r}=*ru*_{1}+*r*^{3}*u*_{3}+⋯, where*x*and*r*are axial and radial coordinates, respectively. From the rotationally symmetric three-dimensional equations of motion and boundary conditions on the lateral surface, a hierarchy of wave equations was derived.↵In Dai (1998), the author studied the travelling wave solutions of the equation

*u*_{t}−*u*_{xxt}+3*uu*_{x}=2*u*_{x}*u*_{xx}+*uu*_{xxx}, which can be written as a singular ODE, i.e.*u*_{ζ}^{2}=(*u*^{3}−*V**u*^{2}−2*d*_{1}*u*−2*d*_{2})/(*u*−*V*), where*d*_{1}and*d*_{2}are two integration constants and*ζ*=*x*−*Vt*. By phase plane analysis, there are three different types of phase planes due to the locations of the singular line. The author obtained 12 types of bounded travelling wave solutions including the important peakons and anti-peakons.↵In Dai & Huo (2000), the author studied solitary shock waves in compressible hyperelastic rods. By using a non-dimensionalization process and the reductive perturbation technique, the author derived a new type of nonlinear dispersive equation. A singular ODE was obtained when travelling waves were considered,

*u*_{ζ}^{2}=(−2*d*_{2}−2*d*_{1}*u*+*Vu*^{2}−*u*^{3})/(*u*−*V*) , where*d*_{1}and*d*_{2}are integration constants;*γ*is a parameter; and*ζ*=*x*−*Vt*. By phase plane analysis, there are many different phase planes when parameters vary. The author obtained periodic shock waves, periodic waves, solitary waves, etc. Significantly, solitary shock waves are found due to the singular lines in phase planes. For this kind of solution, its first-order derivatives are discontinuous at the wave peak.↵In Dai (2001), the author derived an approximate one-dimensional rod equation for a Mooney–Rivlin elastic rod. To study the travelling-wave solutions of the model equation, a singular ODE was found,

*λ*_{ζ}^{2}=(*λ*^{4}+*D*_{1}*λ*^{3}−2*D*_{2}*λ*^{2}−2*b*_{1}*λ*−*b*_{2})/(*λ*−*b*_{2}), where*λ*is the axial stretch;*ζ*=*z*−*Vt*,*b*_{1}and*b*_{2}are constants related to material; and*D*_{1}and*D*_{2}are constants of ratio of the integration constants to the parameters of material. Possible values of those constants were classified to get physically acceptable solutions. By phase plane analysis, the author obtained 10 types of travelling-wave solutions, including the important solitary cusp wave, which are localized with a discontinuity in the shear strain at the wave peak.↵In Dai

*et al*. (2004), the author mentioned three singular dynamics. One of those is , (3*u*^{2}−2*Vu*+2*d*−*γy*^{2})/2(*γu*−*V*) , where the dot (.) denotes the derivative respective to*ξ*=*x*−*Vt*and*d*is an integration constant. This is a model equation for weakly nonlinear long waves in a general compressible hyperelastic rod. By introducing the concept of weak solution for this kind of singular dynamics, the authors showed that compactons can arise in nonlinear elastic rods.↵For this trajectory, , and the numerator is also zero for

*V*=−(1/8*D*_{2}) so it can contact the singular line at B and C.↵In Dai

*et al.*(2004), the authors considered a general singular dynamics , where*x*∈*R*^{n−1}and*y*∈*R*, the functions ,*g*(*x*,*y*) and are smooth, that is, and dot (.) denotes the derivative respective to*t*. Suppose the zero point is a regular point of*r*(*x*,*y*), that is, for any (*x*,*y*)∈*R*^{n}, if*r*(*x*,*y*)=0, then rank . Then is an (*n*−1)-dimensional submanifold in*R*^{n}. Obviously, the vector field of the singular dynamic is not well defined on*M*due to approaching infinity. For solutions with discontinuities in the first-order derivative, which do not belong to the classical solutions, the authors introduced the following notion of weak solutions for the singular dynamics: if*x*(*t*) and*y*(*t*) are continuous or piecewise continuous in*R*and they satisfy and for every function , then (*x*(*t*),*y*(*t*)) is called a weak solution of the singular dynamics.↵From (2.15), the jump conditions here should be , and , where denotes the jump of the quantity. To show that the jump conditions for this non-smooth solution cannot be satisfied, we only need to show [

*W*_{R}]≠0. By using (6.16), it is easy to show that the leading order of [*W*_{R}] at*R*=*a*is approximated by 0.12 that is not zero.- Received November 6, 2007.
- Accepted February 13, 2008.

- © 2008 The Royal Society