Buckling and Post-Buckling of Graphene Tubes

Mohamed B. Elgindi∗, Dongming Wei†, Yu Liu‡, and Hailan Xu§

September 21, 2014

Abstract

Buckling and post-buckling of a long, thin, inextensible tube made of graphene material

acted upon by a uniform normal pressure are considered. A system of differential equations

governing the equilibrium states of such a tube is derived. Perturbation solutions are provided

for the cases of pressure values close to the critical buckling pressures. A Matlab numerical

solver based on Newton’s and shooting methods are developed. The buckling and post-

buckling shapes of the deformed tube subject to various levels of pressure are presented.

Keywords: Graphene, Buckling and Post-Buckling, Deformation of cylindrical

tubes, Deformation of circular tubes, Differential system of equations

1 Introduction

The creation of graphene made significant contributions in many industries, such as in the elec-

tronic devices, the automobile, the aircraft, and the ship industries. Some developments focus on

replacing some critical components in existing electronic devices by graphene made components

in order to make products with higher integrity. The design of new devices with structural com-

ponents made of graphene requires understanding of the strength capabilities of such components

both for small (buckling) and for large (post-buckling) deformations due to small and large forces,

respectively. The purposes of this paper are to determine the buckling loads and the buckling

∗M. B. Elgindi are with Department of Mathematics, Faculty of Science, Texas A&M University at Qatar, Doha,Qatar, Email: Mohamed.elgindi@qatar.tamu.edu.

†D. Wei is with Department of Mathematics, School of Science and Technology, Nazarbayev University,53 Kabanbay Batyr Ave., Astana, 010000, Republic of Kazakhstan, phone: +7 7052972342, email: dong-ming.wei@nu.edu.kz, and he is on leave from University of New Orleans, New Orleans, LA 70148, USA.

‡Y. Liu is with the Department of Electrical Engineering, University of New Orleans, New Orleans, LA 70148,US., Email: lyu2@uno.edu, Phone: +1-504-261-1126. Corresponding author.

§H. Xu is an Engineer of Chongqing Vehicle Test & Research Institute, Chongqing 401122, China, Email:cjxuhailan@cmhk.com.

1

and the post-buckling shapes for a circular, long, inextensible tube made of graphene material

subject to a normal uniform external pressure. This study should be useful in the design of micro

as well as macro mechanical devices subject to mechanical loading conditions. Buckling and post-

buckling of circular tubes made of linear (elastic) material subject to uniform and non-uniform

external pressure are examined analytically, and numerically in many previous works, for exam-

ple, see [1, 2, 3, 4] and the references therein. The corrugation, collapse and stability of carbon

nano-tubes are examined by [5, 6, 7]. Several mathematical formulations for the equations that

describe the equilibrium states of such tubes are given in these references.

Linear (elastic) materials are characterized by the relationship

σ = Eϵ (1)

where σ, E and ϵ are the stress, Young’s modulus and the strain, respectively. Among other

things, for linear material, the critical buckling loads are known to be given by [1, 2, 3]

pc =EI

R3(N2 − 1), (2)

where N ≥ 2 represents the number of the axes of symmetry of the non-circular shapes. An

elastic tube deforms when external stress is applied and recovers to its original shape after the

stress is removed. This deformation does not alter the molecular structure before and after stress,

see [8, 9].

A strain hardening or work hardening is a process to strengthen the material through the

plastic deformation. Plastic deformation is associated with dislocation within the microstructure

of the material, which may result in additional strengthening of the material (see, for example,

[9, 10]). Stress σ and strain ϵ are nonlinearly related for such material. Hollomon’s stress-strain

law is one of the commonly used mathematical description of work hardening isotropic material

σ = K |ϵ|n−1ϵ, 0 < n ≤ 1 (3)

where K is the Bulk modulus, n is the strain hardening exponent. In [8], the authors have studied

the buckling and post-buckling of the Hollomon tubes. It is found that the critical buckling loads

2

are given by

pc =nKInR2+n

(N2 − 1), (4)

where In =∫A|y|n+1dydz is the generalized area moment of inertia, and N ≥ 2 is as in (2).

Post-buckling shapes are also presented in [8].

The purpose of this paper is to extend the methodology of [8] to the buckling analysis of a

graphene tube. The stress and strain equation for graphene has been established as the following

quadratic law [11, 12, 13, 14, 15]:

σ = Eϵ+D |ϵ| ϵ, (5)

where E is Young’s modulus, D = − E2

4σmax, and σmax is the material’s maximal shear stress. As in

[8], due to symmetry, we may consider the deformation of a typical cross section of the tube, that

is, the deformation of a circular, thin, inextensible ring subject to an external uniform pressure

acting normally in the plane of the tube. As will be shown in Section 3, for the graphene tube,

the critical buckling loads are given by

pc =EI1R3

(1− EI22I1Rσmax

)(N2 − 1), (6)

where In =∫A|y|n+1dydz, n = 1, 2, is the generalized area moment of inertia, and N ≥ 2 is as in

(2) and (4). This critical loads formula is a novel extension of the other two found previously in

[1, 2, 3, 8]. Post-buckling shapes are also presented in this paper. We find that numerical results

obtained by using ABAQUS compares favorably with our critical loads given by (6), see Figure 5.

The layout of the paper is as follows. In Section 2 we present the equilibrium equations for a

graphene tube under uniform external pressure as a system of differential equations. In Section 3

we present a perturbation analysis valid for small deformation, prove formula (6) and obtain an

approximate solution. In Section 4 we present numerical simulations for the post-buckling behavior

of a graphene tube at various loading conditions and compare our results to those obtained by

ABAQUS. In Section 5 we give some concluding remarks.

3

2 Mathematical Formulations

The deformation of a circular ring made of graphene material that follows the quadratic law (5),

under a uniform external pressure p acting normally to the surface of the ring (Figure 1), is

considered in this work. The ring is assumed thin and inextensible. The ring remains circular for

Figure 1: The left figure shows an illustration of a thin graphene ring under uniform and normalexternal pressure. The right figure shows forces and moments applied on an infinitesimal grapheneelement of the ring.

small p > 0. As p increases beyond a critical pressure pc, i.e., the first buckling load, the circular

state is no longer stable and the ring deforms into a buckled state. Based on the methods in

[1, 2, 8], by balancing the forces and moments applied on an infinitesimal element of length ds

(see Figure 1), the following differential equations are obtained

dx

ds= cos θ,

dy

ds= sin θ, (7)

and

d2M

ds2+ k

∫kdM

dsds− p = 0, (8)

where (x(s), y(s)) are the coordinates of a point with arc-length s to the origin, M and k = dθds

are the moment and the curvature, respectively.

The moment of linear material is proportional to the curvature. This relationship, however,

becomes more complicated for non-linear material. Let z-axis be in the outer-normal direction

4

extended from the neutral axis of the ring, then the strain ϵ is

ϵ = zk, (9)

with the bending moment M given by

M =

∫A

σzdA. (10)

where A is the area of the deformed ring. Combining (9), (5) and (10) yields

M =

∫A

(Ezk +D|zk|zk))zdA. (11)

The nth moment of inertia In (n = 1, 2) is given below (see [4]):

In =

∫A

zn+1dA,

which together with (11) lead to

M = EI1k +DI2|k|k. (12)

Equation (12) can be viewed as a generalized Euler-Bernoulli bending moment (moment of inertia)

for a graphene tube. See [16, 17] for some similar equations for buckling of column and beams.

Combination of Eq. (7), (8) and (12) leads to the following equilibrium equation:

d2

ds2(EI1k +DI2 |k| k) + k

∫kd

ds(EI1k +DI2 |k| k) ds− p = 0, (13)

together with Eq. (7). Equations (13) and (7), through normalization, can be rewritten using

non-dimensional variables and parameters as

d2

ds2(k − α|k|k) + 1

2k3 − 2α

3|k|3k − ck − λ = 0, (14)

and

dx

ds= cos θ,

dy

ds= sin θ, (15)

5

where x, y and s are normalized by the tube radius R, and α = |D|I2EI1R

, λ = PR3

EI1and c is an

integration constant to be determined. For N ≥ 2 axes of symmetry the associated boundary

conditions are:

k′(0) = 0, k′(2π

N

)= 0,

∫ 2πN

0

k(s)ds =2π

N. (16)

The constant c is determined along with the solutions of (14) and (15) subject to the conditions

in (16) and x(0) = y(0) = 0.

3 Perturbation of Circular Solution

For values of λ less that the first buckling load λc, the ring remains circular, and as λ exceeds λc,

it buckles into non-circular states. For λ close to λc, without loss of generality, we may assume

that k ≥ 0, and hence Eq. (14) simplifies to

d2

ds2(k − αk2) +

1

2k3 − 2α

3k4 − ck − λ = 0. (17)

To this end, we seek asymptotic expansions in powers of 0 < δ ≪ 1 for k, λ and c in the forms:

k = k0 + δk1 + δ2k2 + δ3k3 + · · · ,

λ = λ0 + δλ1 + δ2λ2 + δ3λ3 + · · · ,

c = c0 + δc1 + δ2c2 + δ3c3 + · · · , (18)

satisfying the boundary conditions

k′i(0) = 0, k′i

(2π

N

)= 0,

∫ 2πN

0

ki(s)ds =2π

N, i = 0, 1, 2, 3, · · · . (19)

Note that k0(s) ≡ 1, c0 = 12 − 2α

3 − λ is a solution of the boundary value problem (14), (15) and

(16) for every λ > 0, which we name as the basic solution. The bifurcation points and bifurcation

branches off this basic solution for each N ≥ 2 are determined in this section. Note that the

coefficients of δ in (18) are functions of s and are independent of δ. By substituting (18) into (17)

and equating the coefficients of the powers of δ (up to the fourth-order) to zero, we obtain the

6

following system of equations

d2

ds2(k0 − αk20) +

1

2k30 −

2α

3k40 − c0k0 − λ0 = 0, (20)

d2

ds2(k1 − 2αk0k1) +

3

2k1k

20 −

8α

3k1k

30 − c1k0 − c0k1 − λ1 = 0, (21)

d2

ds2(k2 − 2αk0k2 − αk21

)+

3

2k2k

20 +

3

2k21k0

− 8α

3k2k

30 − 4αk21k

20 − c2k0 − c1k1 − c0k2 − λ2 = 0, (22)

d2

ds2(k3 − 2αk1k2 − 2αk0k3) +

3

2k3k

20 + 3k0k1k2 +

1

2k31

− 8α

3k3k

30 −

8α

3k0k

31 − 8αk2k1k

20 − c3k0 − c2k1 − c1k2 − c0k3 − λ3 = 0, (23)

Since k0 ≡ 1, (20) gives:

c0 =1

2− 2α

3− λ0, (24)

where λ0 is the critical buckling load to be determined from the solvability condition of (21).

Using (24), Eq. (21) becomes:

d2

ds2k1 + β2k1 =

c1 + λ1

1− 2α, (25)

where

β2 =

(3

2(1− 2α)− 8α

3(1− 2α)− c0

1− 2α

), (26)

which has non-trivial solution:

k1 = cos(βs), (27)

provided that:

λ0 = (1− 2α)(N2 − 1), c1 = λ1 = 0. (28)

7

Equation (28) gives the critical loads presented by formula (6) of a graphene material tubes

(through normalization) for each N ≥ 2. With this, (22) gives:

d2

ds2k2 + β2k2 =

c2 + λ2 − 32k

21 + 4αk21 + α d2

ds2 k21

1− 2α

=c2 + λ2 + 2α− 3

4 +(2α− 2αβ2 − 3

4

)cos(2βs)

1− 2α(29)

whose solution is

k2 =2α− 2αβ2 − 3

4

−3(1− 2α)β2(cos(2βs)− cos(βs)) , (30)

provided that c1 = 0, λ1 = 0 and hence

λ2 =3

4− 2α− c2. (31)

Similarly, (23) gives:

d2

ds2k3 + β2k3 =

c3 + λ3 + c2k1 − (3− 8α)k1k2 − ( 12 − 8α3 )k31 + 2α d2

ds2 k1k2

(1− 2α), (32)

which, by Fredholm alternative theory, has a solution if

c2 =2α− 2αβ2 − 3

4

−3(1− 2α)β2

(αβ2 +

3

2− 4α

)+

(3

8− 2α

)(33)

and

λ3 =3

4− 2α−

2α− 2αβ2 − 34

−3(1− 2α)β2

(αβ2 +

3

2− 4α

)−

(3

8− 2α

)(34)

Therefore, we obtain the following perturbation expansions for k(s), c and λ:

k =1 + δ cos(βs) + δ22α− 2αβ2 − 3

4

−3(1− 2α)β2(cos(2βs)− cos(βs)) (35)

c =1

2− 2α

3− (1− 2α)(β2 − 1) + δ2

[2α− 2αβ2 − 3

4

−3(1− 2α)β2

(αβ2 +

3

2− 4α

)+

(3

8− 2α

)](36)

λ =(1− 2α)(β2 − 1) + δ2[3

4− 2α−

2α− 2αβ2 − 34

−3(1− 2α)β2

(αβ2 +

3

2− 4α

)−(3

8− 2α

)](37)

8

Observe that if β = N and α = 0, then Eq. (35)–(37) give the same expansions obtained earlier

in [1, 2, 3] for linear material.

4 Numerical Simulations

We have shown that the critical pressures are given by Eq. (28) for N ≥ 2. k(s) remains positive

for small deformation, however, k(s) may become negative for large deformation and we need to

keep the absolute value sign in (14). This makes this equation different from a usual ordinary

differential equation. Therefore, we need to use special techniques to handle this situation. We

can rewrite (14) as:

d2

ds2(k − α |k| k) + 1

2k3 − 2α

3|k|3 k − ck − λ = 0, (38)

and let

v = ϕ(k) = k − α |k| k. (39)

Then

k = ϕ−1(v) =

1−

√1−4αv2α , v ≥ 0

−1+√1+4αv

2α , v < 0(40)

We write Eq. (8) and (39) as the following system:

y′1

y′2

y′3

y′4

y′5

=

k

y3

−k3

2 + 2α3 |k|3 k + ck + λ

cos(y1)

sin(y1)

, (41)

where k is given in Eq. (40), y1 = θ, y2 = v, y3 = dv/ds, y4 = x and y5 = y with the boundary

conditions

y1(0) = y3(0) = y4(0) = y5(0) = 0, y1(2π

N) =

2π

N, y3(

2π

N) = 0. (42)

9

−1.5 −1 −0.5 0 0.5 1 1.50

0.5

1

1.5

2

2.5

x−axis

y−axis

λ=5.1

λ=4.5

λ=3.5

λ=3.0

−1 −0.5 0 0.5 10

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

x−axis

y−axis

λ=10

λ=20.541

λ=15

λ=8

Figure 2: Post buckling deformation for N = 2 (left) and N = 3 (right), with α = 0.01.

−1.5 −1 −0.5 0 0.5 1 1.50

0.5

1

1.5

2

2.5

x−axis

y−axis

λ=3

λ=3.5

λ=4

λ=4.4036

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

x−axis

y−axis

λ=15.148

λ=12

λ=8

λ=9

Figure 3: Post buckling deformation for N = 2 (left) and N = 3 (right), with α = 0.05.

We use the initial value solver ode15s in Matlab together with Newton Raphson’s method to

approximate the non-circular solutions of (41)–(42). We use the perturbation solutions of Section

3 to generate initial guesses for λ close to λc. For N = 2 and 3, we present the deformations

for λ beyond the corresponding critical pressures in Figures 2 and 3 for α = 0.01 and α = 0.05,

respectively. In practice, deformations corresponding to N = 2 are the most stable configuration,

especially for long tubes and where deviations are small. The relationships between the buckling

loads and the tube radius are given in Figure 4, where our results compare favorably with the

results by ABAQUS, which uses the Rik-arc length method.

10

0 0.2 0.4 0.6 0.8 1 1.2 1.4

x 10−6

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

R (m)

(log1

0) C

ritic

al B

uckl

ing

Load

(N

/m)

N=2N=3N=2 ABAQUSN=3 ABAQUS

Figure 4: Buckling load vs. ring radius R. The radius of the ring’s wall is chosen to be r = 10nm.Note that the buckling load is scaled by log10, and close agreement has been reached between oursolutions and the results obtained from ABAQUS.

5 Conclusions

A system of differential equations that govern the equilibrium states of an inextensible graphene

tube under uniform and normal external pressure is derived. Perturbation and numerical solutions

are obtained to understand buckling and post-buckling behavior for such graphene tube. The

model and the results obtained extend the previous results of [1, 2, 3, 8] (most importantly

with respect to the plastic deformation) for the (graphene) constitutive quadratic law (5). In

particular, the novel analytic formula for the critical buckling loads (6) is provided. To authors’

best knowledge, our work is the first attempt to study the buckling and the post-buckling of a

graphene material tube. This attempt is based only on the mathematical formulation and analysis.

Experimental work is highly desired and is one of our goals for future work. Although experimental

data are the fundamental evidence to support our formulation and results, in the absence of these

required data at this time, the work in this paper may be considered as an initial development

in analysing structural behavior and failure of graphene tubes due to an applied uniform normal

external pressure.

Acknowledgement

This research was partially supported by the Research Office at Texas A&M University at Qatar.

11

References

[1] M. B. Elgindi, D. H. Y. Yen, and C. Y. Wang, “Deformation of a thin-walled cylindrical tube

submerged in a liquid,” Journal of Fluids and Structures, vol. 6, no. 3, pp. 353–370, 1992.

[2] J. E. Flaherty, J. B. Keller, and S. I. Rubinow, “Post buckling behavior of elastic tubes and

rings with opposit sides in contact,” SIAM Journal on Applied Mathematics, vol. 23, no. 4,

pp. 446–455, 1972.

[3] I. Tadjbakhsh and F. Odeh, “Equilibrium states of elastic rings,” Journal of Mathematical

Analysis and Applications, vol. 18, no. 1, pp. 59–74, 1967.

[4] F. Bloom and D. Coffin, Handbook of Thin Plate Buckling and Postbuckling. New York, NY:

Chapman and Hall, 2000.

[5] J. Xiao, B. Liu, Y. Huang, J. Zuo, K. C. Hwang, and M. F. Yu, “Collapse and stability of

single- and multi-wall carbon nanotubes,” Nanotechnology, vol. 18, 2007.

[6] M. Wang, X. Zhang, Q. S. Zheng, and Y. Liu, “A membrane-spring model for carbon nan-

otubes with Van DerWaals interaction between non-bonded atoms,” Nanotechnology, vol. 18,

2007.

[7] H. Shima and M. Sato, “Multiple radial corrugations in multiwalled carbon nanotubes under

pressure,” Nanotechnology, vol. 19, 2008.

[8] M. B. Elgindi, D. Wei, Y. Liu, K. Kamran, and H. Xu, “Buckling and deformation of Hol-

lomon’s power-law tubes,” Thin-Walled Structures, vol. 74, pp. 213–221, 2014.

[9] D. C. William and G. R. David, Fundamentals of Materials Science and Engineering: An

Integrated Approach. John Wiley and Sons, 2011.

[10] F. Campbell, Elements of Metallurgy and Engineering Alloys. ASM International, 2008.

[11] Q. Lu and M. A. B. R. Huang, “Elastic bending modulus of monolayer graphene,” Journal

of Physics D: Applied Physics, vol. 42, no. 10, 2009.

[12] E. Cadelano, P. L. Palla, S. Giordano, and L. Colombo, “Nonlinear elasticity of monolayer

graphene,” Physical Review Letters, vol. 102, 2009.

[13] E. W. Malina, “Mechanical behavior of atomically thin graphene sheets using atomic force

microscopy nanoindentation,” Master’s thesis, University of Vermont, 2011.

12

[14] E. Cadelano, “Graphene under strain-a combined continuum-atomistic approach,” Ph.D. dis-

sertation, Universita degli studi di Cagliari, 2011.

[15] I. A. Ovid’ko, “Mechanical properties of graphene,” Reviews on Advanced Materials Science,

vol. 34, pp. 1–11, 2013.

[16] D. Wei and Y. Liu, “Analytic and finite element solutions of the power-law Euler-Bernoulli

beams,” Finite element in analysis and design, vol. 52, pp. 31–40, 2012.

[17] D. Wei, A. Sarria, and M. B. Elgindi, “Critical buckling loads of the perfect Hollomon’s

power-law columns,” Mechanics Research Communications, vol. 47, pp. 69–76, January 2010.

13