Buckling load analysis of cracked curved beams using differential 

quadrature element method 
 

M. Zare1, A. Asnafi2*1 
1Department of Mechanical Engineering, Shahid Chamran University, Ahvaz, 61355, 

Iran. 
2Hydro-Aeronautical Research Center, Shiraz University, Shiraz, 71348-51154, Iran. 

 

 

ABSTRACT  

 

This paper investigated the buckling load of a cracked curved beam subjected to external 

excitations considering the effects of shear deformations and geometric nonlinearity due 

to large deformations. The governing nonlinear equations of motion were derived. The 

analysis in stationary case was developed for each half of the beam and then the 

differential quadrature element method (DQEM) was used to discretize and solve the 

problem. The resulting nonlinear system of equation was analyzed using continuity 

conditions between the beam segments and an arc length strategy. To verify the validity 

of the proposed method, the beam was modeled using the finite element method. The 

agreement between the results showed the accuracy of the proposed method in prediction 

the buckling load of the beam.  Finally, the effect of crack parameters (depth and location) 

on the buckling load was investigated. As the results showed, the crack depth and buckling 

load related conversely. Furthermore, the closer the crack to the midpoint, the less load 

was required to make the beam undergo buckling.  

 

 

KEYWORDS: Differential quadrature element method (DQEM), Curved beam vibration, 

Buckling load, Arc length method, Non-linear analysis. 

 

1.0  INTRODUCTION 

 

Curved beams have been applied in a variety of industries to be used as mechanical 

devices, building arches and so on. Obviously, over time, such structures can be damaged 

and one of the most prevalent damages is crack. Moreover, the stability of structures is of 

high significance; so, it is a great idea to have a deeper insight into the stability analysis 

of cracked structural elements.  

Much research has been done about the stability of beams as well as their damage 

phenomena. Bradford et al. (Bradford, 2002) studied the in-plane elastic stability of 

arches under a central concentrated load analytically. The stability of a cracked beam 

exited by a follower load was done by Wang (Wang, 2004). In that research, the critical 

load was obtained on the basis of the variation of resonant frequencies. In-Soo et.al (Son, 

2007) studied the natural frequencies of a cracked beam exposed to a follower force. In 

their study the buckling load of the beam was also calculated. The static stability analysis 

of a uniform column with multiple cracks has been done by Caddemi et.al (Caddemi, 

                                                           
*Corresponding Author. Email: asnafi@shirazu.ac.ir 

mailto:email@address.ac.ir


 

 

2013). In that research the buckling modes as well as the corresponding buckling loads 

were presented. In another study done by Nikpour (Nikpour, 1990), the buckling 

phenomenon of a beam with an edged-notched composite was analyzed. In the mentioned 

study, the local compliance due to the crack was considered as a function of the crack-tip 

stress intensity factors and the elastic properties of material. The buckling and post 

buckling analysis of curved beams under distributed, concentrated and thermal loads were 

done by Eslami (Eslami, 2018).  In that study, different boundary conditions and loads 

were studied and then the pre and post buckling loads of rings were obtained and 

discussed.  Attard (Attard, 1986) presented two new finite element formulations for 

obtaining the lateral buckling load of elastic beams under static loads. Torabi et al. 

(Torabi, 2014) studied the free vibration of a Timoshenko beam with multiple cracks 

using differential quadrature element method (DQEM). In that study, they revealed that 

how crack parameters influenced the natural frequencies. 

In this study the buckling load of a general cracked curved beam with a radial 

concentrated force at middle point of the beam is investigated. The problem is solved 

considering the static analysis based on the differential quadrature element and arc length 

methods and finally, a formula is proposed for the buckling load with respect the radial 

displacement of the mid-point of the beam. Also, the effect of crack parameters, depth 

and location, on the buckling load is studied.  

 

 

2.0  EQUATION OF MOTION OF THE CRACKED CURVED BEAM 

 

The equations of motion for a curved beam's post-buckled state, taking into account the 

effects of shear deformation and rotary inertia, as well as, the extension of the neutral 

axis, can be written as (Nikpour, 1990): 
1

( )
Q

N AU
R S S




 
   
 

 (1) 

1
( )

N
Q AW

R S S




 
  
 

 (2) 

( 1)
M Q N

N Q A
S kAG EA

 


   


 (3) 

(( 1) cos( ) ( )sin( ) 1)
N U W W U

EA S R S R
 

 
     

 
 (4) 

( ( 1)sin( ) ( ) cos( ))
Q U W W U

kAG S R S R
 

 
     

 
 (5) 

M EI
S





 (6) 

where dot means the derivative with respect to time. As shown in Fig. 1, which represents 

an element of a curved beam, W, U and φ denote the radial and tangential displacements 

and the angle of rotation. M, N and Q represent the bending moment, normal and shear 

distributed forces respectively. Moreover, A, I, γ, G, E and k are the structural properties 

of the beam, which denote area and moment of inertia of the cross section, mass density 

per unit volume of the beam material, shear and Young's modulus, and shear factor of the 

cross section, respectively. In this figure β is the angular location of the crack.  



 

 

 
Figure 1. Load and displacement components of a curved beam element (Karaagac, 2011) 

 

In order to include the effect of crack, a rotational spring with stiffness (K0) is assumed 

as (Cerri, 2008): 

3 3

0

1 1
, ( ) , ( ) 2

( ) / 12 12

D
D

cD

EI EI
K EI E b h h EI E bh

EI EI h



   


 (7) 

The equations of motion are solved by DQE method which is a new numerical method 

for rapidly solving linear and nonlinear differential equations (Appendix A). 

At equilibrium state, the terms containing time evolutions in Eqs. (1-6) are eliminated. 

Then, by applying DQ discretization to the equations of motion at an interior node mi, in 

an element i, the following equilibrium equations are obtained: 

 
( ) ( )1

( ) ( ) 0
i i

i

i i

e ei m m
em

Q
N

R S S

 
  

 
 (8) 

( ) ( )1
( ) ( ) 0

i i

i

i i

e ei m m
em

N
Q

R S S

 
  

 
 (9) 

( ) ( ) ( )
( ) ( ) ( 1) 0

i i i

i i

i i i

e e ei im m m
e em m

M Q N
N Q

S KAG EA


   


 (10) 

( ) ( ) ( ) ( ) ( )
(( 1) cos(( ) ) ( )sin(( ) ) 1) 0

i i i i i

i i

i i i i i

e e e e ei im m m m m
e em m

N U W W U

EA S R S R
 

 
      

 
 (11) 

( ) ( ) ( ) ( ) ( )
( ( 1)sin(( ) ) ( ) cos(( ) )) 0

i i i i i

i i

i i i i i

e e e e ei im m m m m
e em m

Q U W W U

kAG S R S R
 

 
      

 
 (12) 

( )
( ) 0

i

i

i

ei i m
em

M EI
S


 


 (13) 

Now, the continuity conditions are applied to each interface of the discretized segments 

of the beam. The continuity conditions make some relations between the radial and 

tangential displacements, the angular rotation, the normal and shear forces and the 

bending moment of adjacent elements.  

The radial and tangential displacements and the angular rotation continuity conditions 

at the inter-element boundary of two adjacent elements i and i + 1, except for the crown 

of the beam and the cracked section, are expressed as: 
1 1 1

1 1 1
, ,i i i

i i i i i i

N N N
W W U U  

  
    (14) 

In addition, the normal and shear forces and the bending moment continuity conditions 

at the inter-element boundary of two adjacent elements i and i + 1 can be expressed, 

respectively, as: 



 

 

1 1 1 1
1 1 1 11 1 1 1

1 1

(( 1) cos( ) ( )sin( ) 1)

(( 1) cos( ) ( )sin( ) 1)

i i i i

i i

i i i i

i i i iN N N N

N N

i i i i
i i i i

E A

E A

u w w u

S R S R

u w w u

S R S R

 

 
   

   

 
     

 

 
    

 

 (15) 

1 1 1 1
1 1 1 1 11 1 1 1

1 1

( ( 1)sin( ) ( ) cos( ))

( ( 1)sin( ) ( )

 

cos( ))

i i i i

i i

i i i i

i i i i iN N N N

N N

i i i i
i i i i i

u w w u

S R S R

u w w u

k A G

S R S
k A

R
G

 

 
   

    

 
     

 

 
    

 

 
(16) 

1
1 1 1

1

 
i

i

i i
i i i iN

i i

N

E I E I
S S

 


 
 


   

(17) 

The continuity conditions at the crown of the arch in the equilibrium state are: 
1 1 1

1 1 0 1 0

1 1

1 1, ,
2

 , , ,i i i

i i

i i i i i i

N N N

i ii i

N N
u u m

i
S

U U W w W w

S S S

 

 

  

   
  

   

   

 (18) 

where, w0 is the radial displacement of the middle point of the beam. Also, the continuity 

conditions at the cracked section is: 
1 1 1

1 1 0 1
, , ( )i i i

i i i i i i

cN N N
W W U U K M 

  
     (19) 

Where Mc  is the bending moment at the cracked section. Finally, the boundary 

conditions for a beam clamped at both ends are: 
1 1 1

1 1 1
 0    0     0 W U    (20) 

 0 0 0i i i
m m m

N N N
W U     

 (21) 

 

To compute the buckling load, the radial displacement of the crown of the beam (w0) is 

used as the input of the arc length strategy to obtain the concentrated load in each half of 

the arch. In this manner, by investigation the normal force N, shear force Q, and the 

angular rotation  at the arch crown, the concentrated load in each half of the beam is 
calculated as (Zhu, 2014): 

 

Now, by applying the continuity conditions at the crown point of the beam, the value of 

buckling force is obtained as: 

 

 

3.0     RESULTS AND DISCUSSION  

 

Without any loss of generality, using Eqs. (22 & 23), the magnitude of the concentrated 

load versus the radial displacement of the middle point (w0), for a clamped-clamped beam 

with the properties of Table.1 was obtained and plotted in Fig. 2. It is to be noted that to 

ensure about the accuracy of the proposed method, a finite element simulation was also 

done and the results obtained throughout above methods were compared.   

1 1 1 1

1 2 1 1 1 1
 sin( ) cos( ), sin( ) cos( )

i i i i i i i i

N N N N
F N Q F N Q   

   
     (22) 

1 2
,

2
c

m
P F F i    (22) 



 

 

 
Table 1. Mechanical properties of the curved beam 

Property Notation Value 

Radius of the beam axis R 83 (cm) 

Opening angle of the beam θ 40 (deg) 

Height of the cross section h 0.5 (cm) 

Base of the cross section b 2 (cm) 

Young’s Modulus E 11 (GPa) 

Poisson’s Ratio v 0.3 

Density γ 7800 (kg/m3) 

 

A good agreement between the results obtained by our model and those obtained 

through FE modeling (ANSYS) was seen (see Fig. 2) which confirmed the accuracy of 

the proposed method.  

 

 
 

Figure 2. Variation of the concentrated load versus the crown radial displacement 

 

 

Note also that in the FE simulation, the 3-node element BEAM189 was used to mesh and 

analyze the nonlinear buckling behavior of the beam. In the analysis, 100 elements were 

used to achieve more accurate results. See the finite element model along with the 

structured meshes in Fig.3.  

 

0

100

200

300

0 0.02 0.04 0.06 0.08 0.1

P
(N

)

w0(m)

FEM DQEM



 

 

 
(a) 

 
 

 

 

 

 

(b) 

 
Figure 3. a: The finite element model of the beam. b: Assumed element for meshing the 

model 

 

3.1     The effect of crack depth on the buckling load 

 

In Fig. 4, the changing of the buckling load for the beam introduced in Table. 1 by a 

variation in the crack depth was shown. In this figure the relative location of the crack 

(the proportion of the angular location of the crack to the opening angle of the beam) is 

0.17. Note that in this figure the relative depth is defined as the crack depth to the height 

of the cross section. 

 

 
 

Figure 4. Changing of the buckling load versus the midpoint displacement for different crack 

depths 

 

 

0

100

200

300

0 0.02 0.04 0.06 0.08

P
(N

)

w0 (m)

relative depth=.1

relative depth=.3

relative depth=.5



 

 

 

As the figure shows, increasing the crack depth make a decrease in the buckling 

load which is due to the decrease in the flexural stiffness of the beam.   

 

3.2     The effect of crack location on the buckling load 

 

In Fig.5, it is shown how the buckling load changes with the variation of the 

crack location for the aforementioned beam with relative crack depth of 0.5. It 

represents that as the crack location gets closer to the middle point of the beam, 

the buckling load decreases. This is mainly because the bending moment 

increases by moving toward the midpoint of the beam. Also, the crack effect is 

related to the bending moment as Eq.(19) shows; so, the buckling load 

decreases. 

 

 
Figure 5. Changing of the buckling load versus the midpoint displacement for different crack 

locations 

 

 

4.0     CONCLUSION  

 

The differential quadrature element method along with an arc length strategy were used 

to obtain the buckling load of a cracked curved beam in this article. Using the equation 

of motion in stationary case, and applying the continuity conditions between adjacent 

segments of the beam, a formula for buckling load was proposed. The proposed method 

was firstly validated by a finite element simulation. After that, the effects of the crack 

depth and location on the buckling load were studied. It was shown that the buckling load 

became smaller as the depth of the crack increased or the crack location approached to 

the beam crown.  

 

 

5.0     CONFLICT OF INTEREST 

0

100

200

300

0 0.02 0.04 0.06 0.08 0.1

P
(N

)

w0 (m)

relative location=0.25

relative location=.083

relative location=0.42



 

 

 

We declare that we have no conflict of interest. 

 

 

 

6.0    REFERENCES 
 

Attard, M. M. (1986). Lateral buckling analysis of beams by the FEM. Computers & 

structures, 23(2), 217-231. 

Bradford, M. A. (2002). In-plane elastic stability of arches under a central concentrated 

load. Journal of engineering mechanics, 128(7), 710-719. 

Caddemi, S. C. (2013). The influence of multiple cracks on tensile and compressive 

buckling of shear deformable beams. International Journal of Solids and 

Structures, 50(20-21), 3166-3183. 

Cerri, M. N. (2008). Vibration and damage detection in undamaged and cracked circular 

arches: experimental and analytical results. Journal of sound and vibration, 

314(1-2), 83-94. 

Eslami, M. R. (2018). Buckling and Post-buckling of Curved Beams and Rings. In 

Buckling and Postbuckling of Beams, Plates, and Shells (pp. 11-188). Springer. 

Karaagac, C. O. (2011). Crack effects on the in-plane static and dynamic stabilities of a 

curved beam with an edge crack. Journal of Sound and Vibration, 330(8), 1718-

1736. 

Nikpour, K. (1990). Buckling of cracked composite columns. International Journal of 

Solids and Structures, 26(12), 1371-1386. 

Son, I. S. (2007). Stability Analysis of Cracked Cantilever Beam with Tip Mass and 

Follower Force. Transactions of the Korean Society for Noise and Vibration 

Engineering, 17(7), 605-610. 

Torabi, K. A. (2014). A DQEM for transverse vibration analysis of multiple cracked non-

uniform Timoshenko beams with general boundary conditions. Computers & 

Mathematics with Applications, 6(3), 527-541. 

Wang, Q. (2004). A comprehensive stability analysis of a cracked beam subjected to 

follower compression. International Journal of Solids and Structures, 41(18-19), 

4875-4888. 

Zhu, J. A. (2014). In-plane nonlinear buckling of circular arches including shear 

deformations. Archive of Applied Mechanics, 84(12), 1841-1860. 

 

APPENDIX A. DIFFERENTIAL QUADRATURE ELEMENT METHOD 

 



 

 

The DQEM is a new numerical method for rapidly solving linear and nonlinear 

differential equations. The DQEM is based on the DQ method, an approximate method 

for expressing partial derivatives of a function at a point located in the domain of the 

function, as the weighted linear sum of the values of the variable function at all the defined 

precision points in the derivation direction. Eq. (A.1) is the mathematical representation 

of the DQ expansion: 

where f is the function, N is the number of precision points, xi is the precision associated 

with the i-th point of the function domain, and  represents the weighting coefficients used 

for finding the first derivative of the function at the i-th precision point of the function 

domain represented as xi above. In the DQEM, the studied structure is divided into several 

elements. Then, the continuity conditions are applied on the inter-element boundary of 

two adjacent elements and the boundary conditions of the beam as well as the governing 

equations on each element, using the differential quadrature method. According to Eq. 

(A.1), two important factors in the DQ method include: calculation of the DQ weighting 

coefficients, and selecting the precision points. The Lagrangian functions were used to 

compute the weighted coefficients, and the Gauss–Lobatto Chebyshev polynomial was 

used for selecting the precision points. 

(1)

 

1

|            , 1, ,3, , 2
i

N

x x ij j

j

df
C f i j N

dx




    (A.1)