Page 1 of 20

International Journal for Numerical Methods in Engineering

Post-buckling Solutions of Hyper-elastic Beam by
Canonical Dual Finite Element Method
Kun Cai1, David Y. Gao2,3* and Qing H. Qin3
1

College of Water Resources and Architectural Engineering, Northwest A&F University, Yangling, 712100, China

2

School of Science, Information Technology and Engineering, University of Ballarat, Mt Helen, Victoria 3353, Australia

3

Research School of Engineering, Australia National University, Canberra ACT 0200, Australia

* Corresponding author’s email: d.gao@ballarat.edu.au (David Y. Gao).

Pe

Abstract: Post buckling problem of a large deformed beam is analyzed using canonical dual finite

er

element method (CD-FEM). The feature of this method is to choose correctly the canonical dual
stress so that the original non-convex potential energy functional is reformulated in a mixed
complementary energy form with both displacement and stress fields, and a pure complementary

Re

energy is explicitly formulated in finite dimensional space. Based on the canonical duality theory
and the associated triality theorem, a primal-dual algorithm is proposed, which can be used to find

vie

all possible solutions of this nonconvex post-buckling problem. Numerical results show that the
global maximum of the pure-complementary energy leads to a stable buckled configuration of the
beam. While the local extrema of the pure-complementary energy present unstable deformation

w

states, especially. We discovered that the unstable buckled state is very sensitive to the number of
total elements and the external loads. Theoretical results are verified through numerical examples

On

and some interesting phenomena in post-bifurcation of this large deformed beam are observed.

non-convex variational problem, global optimization

ly

Keywords: Canonical dual finite element method, post buckling, nonlinear beam model,

1. Beam model and motivation
Large deformed cantilever beam to be studied is shown in Fig. 1, which is subjected to a
distributed load q(x) and axial compressive force F. To solve this problem using FEM, the beam is
discretized with several beam elements, for the e-th element its two ends are marked as A and B. The
deflections of the two ends are wAe and wBe , respectively. The rotating angles of the two ends are

 Ae and  Be , respectively.

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 2 of 20

Fig. 1 Beam model (cantilever type)

Pe

In traditional nonlinear elastic beam theory, the stress in lateral direction (normal to the beam
axis) is neglected and the governing equations can be considered as the one-dimensional von

er

Karman model [1]:

 x, x  0

(1)

Re

EIw, xxxx   x w, xx  q( x)  0

The first equation implies a constant stress field  x which means that this von Karman beam

vie

model is actually a linear ordinary differential equation. However, if the thickness of a beam is quite
large, e.g., thickness is of 10% to 20% of beam span, the deformation in lateral direction couldn’t be
neglected. To solve this problem mathematically, a nonlinear beam model was proposed by Gao in

w

1996 [2] which is controlled by the follow nonlinear differential equation
EI  w, xxxx   Ew,2x w, xx  E w, xx  f ( x)  0, x   0 L 

(2)

On

where   3h(1  v 2 )  0 , I  2h3 3 ,   (1  )(1  v 2 ) F / E  0 and f ( x)  (1  v 2 )  q ( x) . The height
of the beam is 2h. E is the elastic modulus of material and v is the Poisson’s ratio. L is the length of

ly

the beam.

In this non-linear beam model, the following assumptions are considered: (a) cross-sections of
beam are initially uniform along the beam axis and have a symmetry axis about which bending
occurs; (b) cross sections remain perpendicular to the beam axis before and after deformation and
shear deformation is ignored (Kirchhoff-Love hypothesis); and (c) the beam is under moderately
large elastic deformation, i.e., w(x) ~ h/L and  ~ w,x (x) ([2-4]). Comparing with the
Euler–Bernoulli beam model, the present model ignores neither the stress nor the deformation of
cross section in the lateral direction. The axial displacement field of the present beam model is given
by the following equation [2]:

http://mc.manuscriptcentral.com/nme

Page 3 of 20

International Journal for Numerical Methods in Engineering

ux  

1

1   w,2x 
2
2h 1  

(3)

which implies that the axial deformation could be relatively larger. Therefore, this nonlinear beam
model can be used for studying both pre and post buckling analysis for a large class of real-world
problems in engineering and sciences ([5-9]).
The total potential energy 



p

p

 w U a  R  of the beam associated with Eq.(2) is given by

 w  0 
L

L
1
1
1

EI  w,2xx  E  w,4x  E  w,2x  dx   f ( x)  wdx
0
12
2
2


(4)

Pe

where U a is the admissible deformation space of beam in which certain necessary boundary
conditions are given. It is known from classic beam theory that the Euler buckling load can be

er

defined as

Re

Fcr  inf



wU a

L

0

EI  w,2xx dx



L

0

(5)

w,2x dx

If the end load F  Fcr , the beam is in an un-buckled state. The total potential energy  p is convex

vie

and the non-linear differential equation (2) has only one solution. However, if F  Fcr , the beam is
in a post-buckling state. In this case, the total potential energy 

p

is non-convex and Eq. (2) may

w

have at most three (strong) solutions [10] at each material point x   0 L  : two minimizers (one is
a global minimizer and the other is a local minimizer), corresponding to the two stable buckling

On

states, and one local maximizer, corresponding to an unstable buckling state. As we use numerical
methods to solve the following non-convex variational equation
p

 w,  w   0 
L

L
1
1
1

EI  w,2xx  E  w,4x  E  w,2x  dx   f ( x)  wdx = 0 ,
0
12
2
2


ly



(6)

we must encounter the non-uniqueness in a finite dimensional space. However, to find global
optimal solution of a non-convex problem is usually NP-hard due to the lack of a global optimality
condition [11].

It is well-known that for convex problems, the Hellinger–Reissner energy is a

saddle-point functional, which connects each primal (potential energy) variational problem with an
equivalent complementary dual problem. In contrast, for non-convex problems, the extremal
property of the generalized Hellinger–Reissner principle and the existence of a purely stress-based
complementary variational principle were two well-known debates existing in nonlinear elasticity

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 4 of 20

for over 40 years ([12-16]). The first problem was partially solved by Gao and Strang (1989) [17].
By introducing a so-called complementary gap function, they recovered a broken symmetry in
nonlinear governing equations of large deformation problems, and they proved that this gap function
provides a global optimality condition. The second open problem was solved by Gao ([3, 18]) and a
pure complementary energy variational principle was first proposed in both nonlinear beam theory [3]
and general nonlinear elasticity [19].

A general review on this history was given in [20].

In the work by Gao ([2, 3]), a canonical dual transformation was presented, i.e.,



E 2
w, x  E .
3

(7)

Pe

By using this canonical dual stress, the generalized total complementary energy of the beam

 : U a  S a  R can be expressed as
1

1

3



  w,      EI  w,2xx    w,2x 
  E   f ( x)  w  dx ,
0
2
4 E
2


er

L

2

(8)

Re

where U a , S a are admissible spaces of deflection w and the dual stress  , respectively.
The Gao-Strang complementary gap function [17] for this beam model is defined as

vie

L  EI


G ( w,  )   
w,2xx  w,2x dx .
0
2
 2


(9)

Clearly, if the beam is subjected an extensive axial load F, the stress field  should be positive all

w

over the domain and the gap function is convex in the displacement field w. In this case, the
generalized complementary energy   w, σ  is a saddle functional and the total potential energy

On

 p ( w ) is strictly convex, which leads to a unique global stable solution. If the axial load F is a
compressive force, the stress field  should be negative over [0, L]. However, as long as the axial

ly

stress is less than the Euler (pre-) buckling load, the gap function should be positive and the total
potential energy is still convex. Therefore, this positive gap function provides a global stability
criterion for general large deformation problems [17]. Seven years later it was discovered that the
negative gap function can be used to identify the biggest local extrema in post-buckling analysis [3].
It turns out that a so-called triality theory was proposed by Gao in 1997. Furthermore, a pure
complementary energy principle for finite elasticity theory was established in 1999 [18]. Since then,
the canonical duality theory was gradually developed [21]. This theory is composed mainly of 1)
canonical dual transformation, 2) a complementary-dual variational principle, and 3) the triality
theory. Detailed information on this theory and its extensive applications in nonconvex mechanics
and global optimization can be found in the monograph [21] as well as the review article [11].

http://mc.manuscriptcentral.com/nme

Page 5 of 20

International Journal for Numerical Methods in Engineering

Based on the Gao-Strang generalized complementary energy, the pure complementary energy of
this nonlinear beam can be obtained by

 d (σ )   ( w, σ )  w  ( w, σ )  0

(10)

which is defined on a statically admissible space S a . Then we have the following result.

Theorem 1: Complementary-dual principle
The complementary energy  d (σ ) is canonical dual to the total potential energy  p ( w ) in

 p (w) ,

 w, σ  U

 Sa is a critical point of   w, σ  , then w  U a is a critical point of

Pe

the sense that if

a

σ  S is a critical point of  d (σ ) and  p (w)   (w, σ )   d (σ ) .

er

In computational mechanics, it is well-known known that the traditional FEMs are based on the

Re

potential variational principle, which produces upper bound approaches to the related
boundary-valued problems. The dual finite element method is based on the complementary energy

vie

principle, which was original studied by Pian et al (1964, 1978, 2006) ([22-24]) and Belytschko et al
(1968, 1970, 1975) ([25-27]) mainly for infinitesimal deformation problems (although the buckling
analysis with finite prebuckling deformations had also been studied by Glaum, Belytschko and

w

Masur in 1975) [27]. The mixed/hybrid FEMs are based on the generalized Hellinger–Reissner
complementary energy principle, which have certain advantages for solving both elastic deformation

On

problems ([29, 30]) and large-scale structural plastic limit analysis [28].

Mathematically speaking, numerical discretization for nonconvex variational problems should

ly

lead to global optimization problems which could possess many local extrema. It is well-known in
computational science that traditional direct approaches for solving nonconvex minimization
problems in global optimization are fundamentally difficult or even impossible. Therefore, most of
nonconvex optimization problems are considered as NP-hard [11]. Unfortunately, this well-known
fact in computer science and global optimization is not fully recognized in computational mechanics.
It turns out that many local search finite element methods have been used for solving large
deformation problems.
The purpose of this article is to bridge the existing gap between computational mechanics and
global optimization by developing a canonical dual finite element method (CD-FEM) for large
deformation of beam model in Eq.(2). This method has been successfully applied for solving

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 6 of 20

nonconvex variational problems in phase transitions of solids governed by Landau-Ginzburg
equation [31] and post-buckling analysis of the nonlinear beam [32]. Due to the piece-wisely
constant stress interpolation and coarse mesh scheme, the CD-FEM proposed in [32] can be used
mainly for finding post-buckled solution of a beam under simple lateral loads. In this paper, our
purpose is to find all possible solutions in the post-buckling analysis of the nonlinear beam model
with complex lateral loads. Based on the canonical duality theory developed in Gao (2000b), we will
first show that by using independent cubic shape functions for deflection and linear interpolation for
dual stress field, the pure complementary energy function can be explicitly formulated in finite
dimensional space. By the triality theory proved recently [33], a canonical primal-dual algorithm is

Pe

then proposed which can be used to find all possible post-buckling solutions. Both stable and
unstable buckled states of the nonlinear beam are investigated using this CD-FEM and illustrated by

er

several examples.

2. Pure Complementary Energy and Triality Theory

Re

Suppose the domain  can be discretized into finite elements e , and each element has two
nodes. Each node has three unknown parameters, e.g., for node “A” of element “e” (see Fig. 1), they

vie

are deflection ( wAe ), rotating angular (  Ae ) and dual stress (  Ae ), for node “B”, they are wBe ,  Be
and  Be . The deformation field and the dual stress field are approximated, separately,
(11)

 eh ( x)  NT  σe

(12)

w

weh ( x)  N wT  we

On

where weT   wAe  Ae

wBe  Be  is the nodal displacement vector of the e-th element,

σeT   Ae  Be  is the nodal dual stress element. Their shape functions are as following

NT 

ly

1

2
 4 (1   ) (2   ) 


 Le (1   ) 2 (1   ) 
8

Nw  

 1 (1   ) 2 (2   ) 
4

L

 e (1   ) 2 (  1) 
8


1
(1   ) (1   )
2

(13)

(14)

The generalized total complementary energy given in Eq. (8) can be expressed in discretized

http://mc.manuscriptcentral.com/nme

Page 7 of 20

International Journal for Numerical Methods in Engineering

form as the following
m 1

1
e 1  2




 h we , σe      weT  Ge (σe )  we  σeT  K e  σe  λeT  σe  f eT  we  ce 
1
2

(15)

1
1
 w T  G (σ )  w  σ T  K  σ  λT  σ  f T  w  c
2
2

where m is the number of nodes in beam, w  R2m , σ  Rm are nodal deflection and dual stress
vectors, respectively; λ  Rm , f  R 2m , and c  R are defined by

3 Le
 3

λe   
 N ( x)  dx 
e 2
4





f ( x) N w  dx    f  x 



Pe

1

fe  

e

1

 

1
1


(16)

Le
L
  1  N w  e d
2

 2

(17)

er

3 2
 3 2 
ce   
 E  dx 
 ELe
e 4
4



(18)

Re

The compliance matrix K  R m  R m is composed by element matrices

vie

L
 3

Ke   
N  NT  dx  e
e 2 E
4 E



2 1
1 2 



(19)

and the Hessian matrix of the gap function G (σ )  R 2 m  R 2 m is obtained by assembling the

w

following matrix in each element
e

 EI  N   N    N   σ   N   N   dx
'' T
w

''
w

g13
g 23
g33

g14 
g 24 
g34 

g 44 

T

e

'
w

' T
w

(20)

ly

 g11 g12

g 22



 sym

'



On

Ge (σe )  

where

3Le
Le
3
3
3
3

 g11  2 t1  5 t2 ; g12  4 t1  20 (t2  t3 ); g13   2 t1  5 t2

3Le
Le
L2e
L2e

g

t

(
t

t
);
g

t

(2t2  t3 ); g 23   g12
 14
1
2
3
22
1
4
20
2
30


L2e
L2e
L2e
L2e
g

t

t
;
g

g
;
g


g
;
g

t

(2t2  t3 )
 24
1
2
33
11
34
14
44
1
4
60
2
30

and

http://mc.manuscriptcentral.com/nme

(21)

International Journal for Numerical Methods in Engineering

t1 

   Ae
   Ae
8EI
; t2  Be
; t3  Be
3
Le
2
2

Page 8 of 20

(22)

By the criticality condition

1
2




1
2

 Ξ h  w, σ     wT  G (σ )  w  σ T  K  σ  λT  σ  f T  w  c 
1

  G (σ )  w  f    w   w T  G, (σ )  w  K  σ  λ    σ  0
2


(23)

we obtain the following two equations:

Pe

G(σ) w  f  0

(24)

1 T
w  G,σ (σ )  w  K  σ  λ  0 ,
2

(25)

where G,σ stands for gradient of G with respect to  . Eq. (24) is actually a discretized

er

equilibrium equation of (2), while Eq. (25) is the inversed constitutive relation. Let Sah be a

Re

discretized feasible stress space such that G is invertible for any given σ  Sah . Then on the
discretized feasible deformation space U ah , the displacement vector w can be obtained by solving

vie

Eq. (26):

w  G 1 ( σ )  f

(26)

Substituting this into the generalized complementary energy Eq. (15), the discretized pure

w

complementary energy can be explicitly given by

1 T 1
1
f  G (σ )  f  σ T  K  σ  λT  σ  c
2
2

On

 d σ   

(27)

In order to identify both global and local extrema, we need the following subspaces:

Sah 

n

G (σ )


0 ,
0 ,

ly


 σ  R

Sah   σ  R n G (σ )

(28)
(29)

Where nσ is the dimension of discretized stress field, correspondingly nw is the dimension of
discretized deformation field. G(σ)

G(σ )

0 means that the matrix G is positive definite, while

0 stands for negative definite.

Theorem 2: Triality theory
Suppose

 w, σ 

is a critical point of Ξ  w , σ  .

(1) The critical point w  U ah is a global minimizer of  p ( w ) if and only if the critical point

http://mc.manuscriptcentral.com/nme

Page 9 of 20

International Journal for Numerical Methods in Engineering

σ  S ah  is a global maximizer of  d  σ  , i.e.,



p

 w   min 

p

 d σ    d σ 
 w   max
σS
h
a

wU ah

(30)

(2) The critical point w  U ah is a local maximizer of  p ( w ) if and only if σ  S ah  is a local
maximizer of  d  σ  .



p

 w   max 

p

 d σ    d σ 
 w   max
σS
h
a

wU ah

(31)

Pe

(3) If σ  S ah  and nw = nσ , then the critical point w  U ah is a local minimizer of  p ( w ) if and
only if σ is a local minimizer of  d  σ  . If σ  S ah  but nw > nσ (the case studied in this

er

paper), the vector w  G 1 (σ )  f is a saddle point of  p ( w ) on U ah , which is a local



p

Re

minimum only on a subspace of U ah such that the dimension of U s , nw_sub , equals nσ , i.e.

 w   min 
wU s

p

 d σ    d σ 
 w   min
σS

a

(32)

vie

This theorem plays an important role in non-convex mechanics and global optimization. It was

w

shown by Gao and Ogden that in nonconvex variational/boundary value problem of Ericksen’s
elastic bar, the global extremum solutions are usually nonsmooth which can’t be captured by using

On

any traditional Newton-type methods [34]. In the following sections, we will show that by this
triality theorem, both global and local extrema of the post-buckling beam can be obtained.

ly

3. Canonical Primal-Dual Algorithm and Flowchart

Based on the triality theory, a canonical primal-dual algorithm for finding all the post-buckling
configurations (both stable and unstable post-buckled states) can be proposed as the following.
Step 1: Initiate parameters: k=0, vectors ( σ (0) ) and matrices;
Step 2: calculate w ( k )  G 1 (σ ( k ) )  f (Eq.(26)) using FEM;
Step 3: Using tri-duality theorem to find stress fields ( σ ( k 1) ) corresponding three different
configurations:
(a) stable post-buckled configuration (Global Maximum)

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 10 of 20

max  ( w ( k ) , σ )
σ

T
1
subject to  w ( k )   G (σ )  w ( k )  0
2

(33)

(b) local unstable configuration (Local Maximum)

max  ( w ( k ) , σ )
σ

subject to

1 (k ) T
w   G (σ )  w ( k )  0

2

(34)

(c) unstable buckled configuration (Local Minimum)

min  ( w ( k ) , σ )
σ

Pe

subject to

Convergence check:
If σ ( k 1)  σ ( k )

σ ( k )  1.0 109

Else
k=k+1

vie

Re

go to step 5

If k< k*
go to step 5

w

Else
go to 2

End
Stop for post processing

ly

The maximum iteration k* equals 100 in this algorithm.

On

End

Step 5:

(35)

er

Step 4:

1 (k ) T
w   G (σ )  w ( k )  0

2

4. Numerical examples
In the following examples, the beams have the same span l =1.0m, height 2h = 0.1m. We assume
that the elastic modulus E = 1000Pa and the Poisson’s ratio  = 0.3.

4.1 Example 1-mesh dependence study
Our first example is a beam (L=1.0m) subjected to a concentrated lateral force with and
f ( x  0.5)  (1  v 2 )  q( x)  1.0 N (see Fig. 2). On the right end, the beam is subjected to a

compressive force F and   (1  )(1  v 2 ) F / E  0.003m 2 . To investigate the mesh dependency of

http://mc.manuscriptcentral.com/nme

Page 11 of 20

International Journal for Numerical Methods in Engineering

our results, the beam is discretized with 10, 20, 30, 40, 50 and 60 elements, respectively.

Fig. 2 Simply supported beam model

The curves with dark blue diamonds show the stable configurations (global maximum of

Pe

 d  σ  ) of the buckled beam with different number of elements. The maximum deflection at the
centre of the beam is near 0.05m for different mesh schemes, which implies the mesh-independency

er

of these stable configurations.

The curves with pink solid squares present the configurations (local maximum of  d  σ  ) of

Re

the beam. Clearly, all deflections are almost close to zero. It means that the beam is nearly in a pure
compressed deformation along the axial direction.

vie

The curves with black triangles demonstrate the unstable configurations (local minimum of

 d  σ  ) of the beam. These curves are smooth but very with the total elements used. Therefore, it

w

shows that these local unstable solutions are mesh-dependent. Especially, this unstable solution does
not appear till the total number of elements (NE) > 20 (see Fig 3). The reason for this

On

mesh-dependency is that this local unstable solution is a saddle point of the nonconvex total
potential energy which is sensitive to the discretization schemes of the original infinite dimensional
problem. This result shows that the triality theory plays an important role in post-buckling analysis.

ly

From the above analysis we can see that, as the total number of elements (NE) > 30, all the
solutions (especially the third type) converge to the certain values. Therefore, we will use NE = 40
elements with the same beam length for all the following examples.

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

0.06

Page 12 of 20

0.06
Global max

Global max

0.04

Local max

Beam deflection/m .

Beam deflection/m .

0.04

Local min

0.02
0.00
-0.02
-0.04

Local max
Local min

0.02
0.00
-0.02
-0.04

-0.06

-0.06
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

Beam Axis/m

(a) NE=10

Local max
Local min

0.02

Global max

0.00
-0.02

-0.06
0.0

0.2

0.4

0.6

Beam Axis/m

Re

-0.04

0.8

Local min

0.02
0.00
-0.02
-0.04
-0.06

1.0

0.0

0.2

0.06

0.4

0.6

0.8

1.0

Beam Axis/m

vie

(c) NE=30

(d) NE=40

0.06

Local max

Global max
Local max

0.04

Local min

0.02
0.00
-0.02
-0.04

Local min

0.02

On

Beam deflection/m .

0.04

w

Global max

Beam deflection/m .

Local max

0.04

er

Beam deflection/m .

0.04

1.0

0.06

Beam deflection/m .

Global max

0.8

(b) NE=20

Pe

0.06

0.6

Beam Axis/m

0.00

-0.02
-0.04
-0.06

0.0

0.2

0.4

0.6

0.8

Beam Axis/m

1.0

0.0

ly

-0.06

0.2

0.4

0.6

0.8

1.0

Beam Axis/m

(e) NE=50

(f) NE=60

Fig. 3 Mesh-dependency of solutions of the buckled beam

4.2 Example 2-Complicated lateral load
Let us to consider a clamped/clamped beam as shown in Fig. 4. The lateral load is assumed to be
piece wise uniform pressures (q(x)) such that f ( x)  (1  v 2 )  q( x)  0.1N / m . The left end is fixed,
while the deflection at the right end is specified as zero and there is a compressive force F. Two

http://mc.manuscriptcentral.com/nme

Page 13 of 20

International Journal for Numerical Methods in Engineering

cases are considered for the end load: (a)   (1  )(1  v 2 ) F / E  0.001m 2 ; (b)   0.05m2 .

Fig. 4 Clamped/simply supported beam model

0.0006

Pe

Global max
Local max
Local min

0.0002

er

Beam deflection/m

0.0004

0.0000

-0.0004

0.2

0.4
0.6
Beam Axis/m

vie

0.0

Re

-0.0002

0.8

1.0

Fig. 5 Three post-buckling configurations of beam as   0.01m2

w

0.0015

On

0.0005
0.0000

ly

Beam deflection/m

0.0010

-0.0005
-0.0010

Global max
Local max

-0.0015

Local min

-0.0020
0.0

0.2

0.4
0.6
Beam Axis/m

0.8

1.0

Fig. 6 Three post-buckling configurations of beam as   0.03m2

Since the lateral load is not uniformly distributed, this example can be used to study more
interesting phenomena of the post-buckling configurations subjected to a combination of the lateral
load (q(x)) and the end compression (F).

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 14 of 20

Fig. 5 and Fig. 6 present the three post-buckling configurations of the beam with   0.01m2
and   0.01m2 , respectively. We can see that the stable buckled beam (global max of  d  σ  ) is
rotationally symmetric with respect to the beam center, but the other two unstable solutions are not
symmetrical and are sensitive to the axial load, the bigger F, the larger axial deformation. Fig 6
shows the beam can have a deformation of over 20% of its length. This fact shows that the present
method can be used for large deformed post-buckling analysis.

4.3 Example 3-Different lateral loads

Pe

Now let us consider a simply supported beam as shown in Fig. 7. The right end of the beam is
subjected to a compressive force F such that   (1  )(1  v 2 ) F / E  0.01m 2 . Three types of lateral

er

loads are considered: (a) a concentrated force f ( x  0.5l )  (1  v 2 )  q  0.1N on the center of the
beam, (b) a uniformly distributed load f ( x)  0.1N / m , and

(c) piece-wise uniform load

(a)

w

vie

Re

f ( x)  0.1N / m .

ly

On

(b)

(c)

Fig. 7 Simply supported beam model

In this example, sensitivity of the post-buckling deformation with the lateral loads is
investigated. Fig. 8 and Fig 9 present stable and unstable post-buckling deformations of the beam
under concentrated and uniformed forces, respectively. Clearly, the stable buckled states are similar,

http://mc.manuscriptcentral.com/nme

Page 15 of 20

International Journal for Numerical Methods in Engineering

but not the unstable buckled states. While for the periodic lateral load, the stable and unstable
post-buckled configurations are similar in shapes but different in magnitudes (see Fig 10), which is
reasonable since the stable buckled state is a global minimal solution to the nonconvex problem, but
the unstable buckling state is only a stationary point of the total potential.

0.01

-0.01
-0.02
-0.03

Pe

Beam deflection/m .

0.00

-0.04

Global max

-0.05

er

Local min

-0.06

0.0

0.2

0.4
0.6
Beam Axis/m

0.8

1.0

Re

Fig. 8 Post-buckling configurations of the beam in case a)

vie

0.00

w

-0.01
-0.02

On

Beam deflection/m .

0.01

-0.03

Global max
Local min

-0.04
0.2

0.4
0.6
Beam Axis/m

ly

0.0

0.8

1.0

Fig. 9 Post-buckling configurations of the beam in case b)

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 16 of 20

0.0003

Beam deflection/m .

0.0002
0.0001
0.0000
-0.0001
-0.0002

Global max
Local min

-0.0003
0.0

0.2

0.4
0.6
Beam Axis/m

0.8

1.0

Pe

Fig. 10 Post-buckling configurations of the beam in case c)

er

4.4 Example 4-Coupling of lateral load and axial compression
Finally, we show some interesting phenomena of the unstable post-buckling state (i.e. the local

Re

minimum of  d  σ  ) for a simply supported beam (see Fig 11) with different concentrated loads:

compressive forces:

vie

f ( x  0.5l )  (1  v 2 )  q = 1.0N, 1.2N , and 1.45N , each one is combined with two different

(a)   (1  )(1  v 2 ) F / E  0.003m 2 and

(b)   0.008m2 . We discovered

that this unstable solution is very sensitive to the external loads as expected.

w
ly

On
Fig. 11 Simply supported beam model

From Fig. 12 we can see that for the given compressive load   0.003m2 , the stronger is the
concentrated load f , the smaller is the center deflection of the unstable buckled state. This case is
easy to understand as the strong concentrated load should push down the unstable buckled beam.
However, this phenomenon turns to the opposite way if the compressive load is increased to

  0.008m2 (see Fig 13). We believe that these phenomena need to have detailed study by both
numerical simulations and experiments. Carefully observation shows that all these unstable
post-buckled configurations are not symmetric to the center. For example, in Fig 12, we can see that

http://mc.manuscriptcentral.com/nme

Page 17 of 20

International Journal for Numerical Methods in Engineering

the left peak (x, y)=(0.25m, 0.0372m) of the black curve (f=1.2N) is lower than the right peak
(0.75m, 0.0466m). For the yellow curve (f=1.45 N), the left peak (0.2m, 0.048m) is higher than the
right peak (0.8m, 0.0463m).

These non symmetric deformations are due to the non symmetrical

boundary conditions.

0.050

0.030

Pe

Beam deflection/m .

0.040

0.020

f=1

0.010

f=1.2

er

f=1.45

0.000

0.0

0.2

0.4
0.6
Beam Axis/m

0.8

1.0

Re

Fig. 12 . Unstable post-buckling deflections of the beam as   0.003m2

vie

0.050

f=1.2

w

0.030
0.020

f=1.45

0.010
0.000
-0.010
0.0

0.2

0.4
0.6
Beam Axis/m

ly

On

Beam deflection/m .

0.040

f=1

0.8

1.0

Fig. 13 Unstable post-buckling deflections of the beam with   0.008m2

5. Conclusions
Based on the canonical duality theory, a mixed finite element method is proposed for solving a
post-buckling problem of a large deformed nonlinear beam. By using independent shape functions
for deflection and dual stress fields, a pure complementary energy function is explicitly formulated
in finite dimensional space.

Combining this pure complementary energy with the triality theory, a

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 18 of 20

canonical primal-dual algorithm is proposed which can be used to find not only the stable buckled
solution, i.e. the global minimal of the nonconvex total potential, but also the unbuckled state and the
unstable buckled state. Our numerical results show that the stable buckled state is indeed stable but
the unstable post-buckled solution is very sensitive to both the total number of meshes and the
external loads. Some interesting phenomena are discovered which deserve to have future study.

Acknowledgements
Financial supports by US Air Force Office of Scientific Research (Grant No. FA9550-10-1-0487),
the National Natural Science Foundation of China (Grant No. 50908190), Youth Talents Foundation

Pe

of Shaanxi Province (Grant No. 2011kjxx02), and the Human Resource Foundation of Northwest
A&F University (Grant No. QN2011125) are fully acknowledged.

er

References

Re

1. K. Washizu, 1968. Variational methods in elasticity and plasticity, Pergamon Press, New York.

vie

2. D.Y. Gao, 1996. Nonlinear elastic beam theory with application in contact problems and
variational approaches, Mechanics Research Communications 23(1):11–17.
3. D.Y. Gao, 1997. Dual extremum principles in finite deformation theory with applications

w

to post-buckling analysis of extended nonlinear beam theory, Applied Mechanics Reviews
50(11):S64–S71.

On

4. D.Y. Gao, 2000a. Finite deformation beam models and triality theory in dynamical
post-buckling analysis. Int. J. Non-Linear Mechanics 35:103–131.

beam, Electronic Journal of Differential Equations 27:1–12.

ly

5. F. M Bengue, M. Shillor, 2008. Regularity result for the problem of vibrations of a nonlinear

6. J. Machalova and H. Netuka, 2011. Bending of a nonlinear beam reposing on an unilateral
foundation, Appl. Comput. Mechanics 5:45–54.
7. K. T. Andrews, Y. Dumont, M. F. M’Bengue, J. Purcell and M. Shillor, 2012. Analysis and
simulations of a nonlinear elastic dynamic beam, Z. Angew. Math. Phys. 63:1005–1019.
8. K.L. Kuttler, J. Purcell, M. Shillor, 2012. Analysis and simulations of a contact problem for a
nonlinear dynamic beam with a crack, Quarterly Journal of Mechanics and Applied
Mathematics 65(1):1–25.
9. J. Ahn, K. L. Kuttler, and M. Shillor, 2012. Dynamic contact of two Gao beams, Electronic
Journal of Differential Equations 2012(194):1–42.

http://mc.manuscriptcentral.com/nme

Page 19 of 20

International Journal for Numerical Methods in Engineering

10. D.Y. Gao, R.W. Ogden, 2008a. Closed-form solutions, extremality and nonsmoothness criteria
in a large deformation elasticity problem, Zeitschrift fur and gewandte Mathematik und Physik
59:498–517.
11. D.Y. Gao, and H.D. Sherali, 2008. Canonical duality: Connection between nonconvex
mechanics and global optimization, in Advances in Appl. Mathematics and Global
Optimization, pp.249-316, Springer, 2009.
12. M. Levinson, 1965. The complementary energy theorem in finite elasticity, ASME
Transactions Journal of Applied Mechanics 32(4):826–828.
13. B. Fraeijs de Veubeke, 1972. A new variational principle for finite elastic displacements,

Pe

International Journal of Engineering Science 10:745–763.
14. W.T. Koiter, 1975. On the complementary energy theorem in non-linear elasticity theory,

er

Trends in Applications of Pure Mathematics to Mechanics, Pitman, London, pp. 207–232.
15. R.W. Ogden, 1975. A note on variational theorems in non-linear elastoplastics, Mathematical
proceedings of the Cambridge Philosophical Society 77:609–615.

Re

16. Z.H. Guo, 1980. The unified theory of variational principles in nonlinear elasticity, Archive of
Mechanics 32:577–596.

vie

17. D.Y. Gao, G. Strang, 1989. Geometric nonlinearity: potential energy, complementary energy
and the gap function, Quarterly Journal of Mechanics and Applied Mathematics
XLVII(3):487–504.

w

18. D.Y. Gao, 1999a. Pure complementary energy principle and triality theory in finite elasticity,
Mechanics Research Communications 26(1):31–37.

On

19. D.Y. Gao, 1999b. General analytic solutions and complementary variational principles for large
deformation nonsmooth mechanics, Meccanica 34:169–198.

ly

20. S.F. Li and A. Gupta, 2006. On dual configurational forces, Journal of Elasticity 84:13–31.
21. D.Y. Gao, 2000b. Duality principles in nonconvex systems: Theory, Methods and Applications,
Kluwer Academic Publishers, Dordrecht, Boston (Now Springer), 472pp.
22. T.H.H. Pian, 1964. Derivation of element stiffness matrices by assumed stress distributions,
AIAA Journal 2:1333–1336.
23. T.H.H. Pian, 1978. A historical note about hybrid elements, Int. J. Num. Meth. Engrg.
12:891–892.
24. T.H.H. Pian, C.C. Wu, 2006. Hybrid and incompatible finite element methods, Chapman and
Hall/CRC, Boca Raton.
25. P.G. Hodge and T. Belytschko, 1968. Numerical Methods for the Limit Analysis of Plates,
Journal of Applied Mechanics, 35:796–801.

http://mc.manuscriptcentral.com/nme

International Journal for Numerical Methods in Engineering

Page 20 of 20

26. T. Belytschko and P.G. Hodge, 1970. Plane Stress Limit Analysis by Finite Elements, Journal
of the Engineering Mechanics Division, ASCE 96:931–944.
27. L. Glaum, T. Belytschko and E.F. Masur, 1975. Buckling of Structures with Finite Prebuckling
Deformations - A Perturbation, Finite Element Analysis, Int. J. Solids and Struct.
11:1023–1033.
28. D.Y. Gao, 1988. Panpenalty finite element programming for limit analysis, Computers and
Structures 28:749–755.
29. QH Qin, 2004. Dual variational formulation for Trefftz finite element method of elastic
materials, Mech. Res. Commun. 31:321–330.

Pe

30. QH Qin and H. Wang. 2009. Matlab and C programming for Trefftz finite element methods,
Taylor & Francis, Boca Raton.

er

31. D.Y. Gao and H.F. Yu, 2008. Multi-scale modelling and canonical dual finite element method
in phase transitions of solids,

Int. J. Solids Struct. 45:3660–3673.

32. H.A.F.A. Santos, D.Y. Gao 2012. Canonical dual finite element method for solving

Re

post-buckling problems of a large deformation elastic beam, Int. J. Nonlinear Mechanics, 47 (2)
240-247.

vie

33. D.Y. Gao and C. Wu, 2012. On the Triality Theory for a Quartic Polynomial Optimization
Problem, J. Industrial and Management Optimization 8(1):229–242.
34. D.Y. Gao, R.W. Ogden, 2008b. Multiple solutions to non-convex variational problems with

w

implications for phase transitions and numerical computation, Quarterly Journal of Mechanics
and Applied Mathematics 61(4):497–522.

ly

On
http://mc.manuscriptcentral.com/nme

