# Buckling of Nonprismatic Column on Varying Elastic Foundation with Arbitrary Boundary Conditions

Authors: Ahmad A. Ghadban; Ahmed H. Al-Rahmani; Hayder A. Rasheed; Mohammed T. Albahttiti
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Buckling of nonprismatic single columns with arbitrary boundary conditions resting on a nonuniform elastic foundation may be considered as the most generalized treatment of the subject. The buckling differential equation for such columns is extremely difficult to solve analytically. Thus, the authors propose a numerical approach by discretizing the column into a finite number of segments. Each segment has constantsE (modulus of elasticity), I (moment of inertia), and β (subgrade stiffness). Next, an exact analytical solution is derived for each prismatic segment resting on uniform elastic foundation. These segments are then assembled in a matrix from which the critical buckling load is obtained. The derived formulation accounts for different end boundary conditions. Validation is performed by benchmarking the present results against analytical solutions found in the literature, showing excellent agreement. After validation, more examples are solved to illustrate the power and flexibility of the proposed method. Overall, the proposed method provides reasonable results, and the examples solved demonstrate the versatility of the developed approach and some of its many possible applications.


---

## Body

## 1. Introduction

By the beginning of the twentieth century, many researchers got into studying the specific cases of bucking using various methods such as continuous and lumped matrix analysis, Finite Element Method (FEM), and Boundary Element Method (BEM) [1]. Dinnik reported the exact solution of simply supported columns with monomial variation in stiffness and axial load in his paper published in 1932 [2]. Simply supported tapered columns were studied and the analytical solutions were presented by Gere and Carter in 1962 [3]. In 1970, Gallagher and Lee solved an axially loaded column with a variable flexural stiffness using Finite Element Method [2]. Additionally, Elishakoff and Bert in 1988 succeeded at improving Rayleigh’s method and obtaining an approximate solution for columns with variable stiffness [2]. A year later, Eisenberger presented the exact buckling solution for columns with variable cross-sections and variable axial load using any polynomial variation and using different boundary conditions [2]. In 1986, Ermopoulos published the solution for buckling of tapered bars under stepped axial loads [4]. Analysis on the buckling of uniform columns with different supports was performed by Iyengar in his book which was published in 1988 [5]. Elastic stability of elastically supported columns under the effect of distributed forces was studied by Lee and Kuo in 1991 [5]. In the same year, Arbabei and Li presented the solution for the buckling of nonprismatic elastic columns [6]. A more recent solution for the exact buckling of constrained stepped columns was presented by Yang and Park in 2003 [7]. Coşkun and Atay analyzed the critical buckling load for elastic columns of constant and variable cross-sections using variational iteration method in 2009 [6]. This variational iteration method produces an approximate solution for the presented problem. In 2012, Huang and Li published their solution for analytically determining the exact critical buckling loads of nonuniform columns [8]. Lee et al. published their calculations of natural frequencies and buckling loads for columns with intermediate multiple elastic springs in 2002 [9]. Atanackovic and Novakovic worked on Lagrange problem, in which they determined the optimal shape of an elastic column on elastic foundation using an area varying equation to represent the column shape [10]. Using Pontryagin’s maximum principle, Levy presented the optimal shape of simply supported columns on elastic foundations in 1990 [10]. Using Galerkin’s method, Lacarbonara solved for buckling and postbuckling of nonuniform nonlinearly elastic rods [11]. To the authors’ best knowledge, this is the first paper that provides a numerical treatment combining the buckling of nonprismatic columns resting on nonuniform elastic foundation with arbitrary end boundary conditions.


## 2. Formulation

### 2.1. The Buckling Equation

The buckling differential equation for the column shown in Figure1 can be derived using the bifurcation method resulting in the following equation:(1)EI∂4v∂x4+2∂EI∂x∂3v∂x3+∂2EI∂x2+P∂2v∂x2+βv=0.Figure 1
Nonprismatic columns with arbitrary boundary conditions resting on nonuniform elastic foundation.The boundary conditions (BCs) are as follows:(1) Shear
(2)
V
0
=
-
β
A
v
0
,
V
L
=
β
B
v
L
.(2) Moment
(3)
M
0
=
K
A
v
′
0
,
M
L
=
-
K
B
v
′
L
,where Ex is the column elastic modulus function, Ix is the column moment of inertia function, P is the column axial load, βx is the elastic foundation spring stiffness function, vx is the column lateral deflection function, x is the distance from the bottom of the column, βA is the bottom lateral discrete spring stiffness, βB is the top lateral discrete spring stiffness, KA is the bottom rotational discrete spring stiffness, KB is the top rotational discrete spring stiffness, L is the length of the column, Vx is the column shear force function, and Mx is the column moment function.Such a differential equation can prove to be extremely difficult to solve analytically for arbitraryE, I, and β functions. The authors are not aware of any analytical mathematical method that can solve this differential equation. Thus, the authors decided to solve the differential equation semianalytically by discretizing the column into a finite number of segments. Each segment has constants E, I, and β, as shown in Figure 2. E, I, and β values for each segment were taken to be the average of each function within that segment of the column. The averages were calculated using the general equation for averaging functions within an interval:(4)fave=∫abfxdxb-a.Figure 2
Idealized prismatic column segment resting on uniform elastic foundation.An exact solution was then obtained for each segment using the classical buckling equation of columns on an elastic foundation. The domain ofx for each segment was taken to be 0≤x≤Δx:(5)∂4vi∂x4+ki2∂2vi∂x2+β^i2vi=0,0≤x≤Δx,ki=PEiIi,β^i=βiEiIi,Δx=LN,where N is number of segments, Δx is segment length, and i is segment number.


### 2.2. Analytical Solution of the Segment Buckling Equation

The general solution to the segment differential equation with arbitrary boundary conditions is complicated and requires the consideration of four different cases. The solution to each case is shown below.Case 1 (βi=0).
The segment differential equation becomes(6)∂4vi∂x4+ki2∂2vi∂x2=0,0≤x≤Δx.
Letvi=eλix, which can be substituted into (6) resulting in(7)λi4+ki2λi2eλix=0→yieldsλi1,2=0,λi3,4=±jki,where(8)j=-1.
The deflection function is then obtained as follows:(9)vix=Aicos⁡kix+Bisin⁡kix+Cix+Di.Case 2 (ki4-4β^i2=0).
Letvi=eλix, which can be substituted into (5) resulting in(10)λi4+ki2λi2+β^i2eλix=0→yieldsλi1,2,3,4=±jki2.
The deflection function is then found to be(11)vix=Aicos⁡kix2+Bisin⁡kix2+Cixcos⁡kix2+Dixsin⁡kix2.Case 3 (ki4-4β^i2>0).
Letvi=eλix, which can be substituted into (5) resulting in(12)λi4+ki2λi2+β^i2eλix=0→yieldsλi2=-ki2±ki4-4β^i22→yieldsλi1,2=±jα1i,λi3,4=±jα2i,where(13)α1i=ki2+ki4-4β^i22,α2i=ki2-ki4-4β^i22.
The deflection function is then obtained to be as follows(14)vix=Aicos⁡α1ix+Bisin⁡α1ix+Cicos⁡α2ix+Disin⁡α2ix.Case 4 (ki4-4β^i2<0).
Letvi=eλix, which can be substituted into (5) resulting in(15)λi4+ki2λi2+β^i2eλix=0→yieldsλi2=-ki2±ki4-4β^i22=-ki22±j4β^i2-ki42.Letλi=ai+jbi implying λi2=ai2-bi2+j2aibi which results in the following two equations with two unknowns, aiandbi:(16)ai2-bi2=-ki22,(17)2aibi=±4β^i2-ki42.From (17) we get(18)ai=±4β^i2-ki44bi→yieldsai2=4β^i2-ki416bi2.Substituting into (16) we obtain(19)4β^i2-ki416bi2-bi2=-ki22→yields16bi4-8ki2bi2+ki4-4β^i2=0→yieldsbi2=ki2±2β^i4.We know thatki4-4β^i2<0→yields(ki2+2β^i)(ki2-2β^i)<0, but it is obvious that ki2+2β^i>0. Therefore, ki2-2β^i<0→yieldsbi2≠ki2-2β^i/4 since bi2>0. This leaves the following: bi2=ki2+2β^i/4→yieldsbi=±2β^i+ki2/2.Using (17) we can then obtain ai as follows:(20)ai=±2β^i-ki22.Therefore,λi is found to be(21)λi=±2β^i-ki22±j2β^i+ki22→yieldsλi1,2=γ1i±jgi,λi3,4=γ2i±jgi,where(22)gi=2β^i+ki22,γ1i=2β^i-ki22,γ2i=-2β^i-ki22.The deflection function is then obtained to be(23)vix=Aieγ1ixcosgix+Bieγ1ixsingix+Cieγ2ixcosgix+Dieγ2ixsingix,where Ai,Bi,Ci,andDi are the constant coefficients for each segment.After solving the segment differential equation, the following compatibility conditions were established at each node between the segments.(1) Displacement
(24)
v
i
0
=
v
i
-
1
Δ
x
,
v
i
Δ
x
=
v
i
+
1
0
.(2) Slope
(25)
v
i
′
0
=
v
i
-
1
′
Δ
x
,
v
i
′
Δ
x
=
v
i
+
1
′
0
.(3) Moment
(26)
E
i
I
i
v
i
′
′
0
=
E
i
-
1
I
i
-
1
v
i
-
1
′
′
Δ
x
,
E
i
I
i
v
i
′
′
Δ
x
=
E
i
+
1
I
i
+
1
v
i
+
1
′
′
0
.(4) Shear
(27)
E
i
I
i
v
i
′
′
′
0
+
P
v
i
′
0
=
E
i
-
1
I
i
-
1
v
i
-
1
′
′
′
Δ
x
+
P
v
i
-
1
′
Δ
x
,
E
i
I
i
v
i
′
′
′
Δ
x
+
P
v
i
′
Δ
x
=
E
i
+
1
I
i
+
1
v
i
+
1
′
′
′
0
+
P
v
i
+
1
′
0
.Finally, in order to obtain the coefficient matrix (P-matrix), the following boundary conditions were established.(1) Shear
(28)
E
1
I
1
v
1
′
′
′
0
+
P
v
1
′
0
=
-
β
A
v
1
0
,
E
N
I
N
v
N
′
′
′
Δ
x
+
P
v
N
′
Δ
x
=
β
B
v
N
Δ
x
.(2) Moment
(29)
E
1
I
1
v
1
′
′
0
=
K
A
v
1
′
0
,
E
N
I
N
v
N
′
′
Δ
x
=
-
K
B
v
N
′
Δ
x
.Solving the following system of equations would then render the buckling load:(30)PC=0,where P is the P-Matrix and C = A1  B1C1D1⋯AiBiCiDi⋯ANBNCNDNT.


### 2.3. TheP-Matrix

The structure of theP-matrix is as follows.(1) One Segment
(31)
4
×
4
B
C
M
a
t
r
i
x
.(2) Two Segments
(32)
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
4
×
4
F
i
r
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
4
×
4
L
a
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
.(3) Three Segments
(33)
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
4
×
4
F
i
r
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
8
×
4
I
n
t
e
r
m
e
d
i
a
t
e
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
4
×
4
L
a
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
.For the case of an assembly with more than three segments, the structure of theP-matrix would be the same as the three segments case but with additional intermediate segment compatibility submatrices that are shifted diagonally downwards. The buckling load would then be the load that produces a singular P-matrix.The submatrices inside theP-matrix are obtained by applying the compatibility and boundary conditions and are dependent on the case that governs each segment, as explained earlier. These submatrices were found to be as shown below.(a) Intermediate SegmentsCase 1 (βi=0).
(34)
1
0
0
1
0
k
i
1
0
-
P
0
0
0
0
0
P
0
-
cos
⁡
k
i
Δ
x
-
sin
⁡
k
i
Δ
x
-
Δ
x
-
1
k
i
sin
⁡
k
i
Δ
x
-
k
i
cos
⁡
k
i
Δ
x
-
1
0
P
cos
⁡
k
i
Δ
x
P
sin
⁡
k
i
Δ
x
0
0
0
0
-
P
0
.Case 2 (ki4-4β^i2=0).
(35)
1
0
0
0
0
k
i
2
1
0
-
P
2
0
0
2
E
i
I
i
k
i
0
k
i
P
2
2
-
P
2
0
-
cos
⁡
k
i
Δ
x
2
-
sin
⁡
k
i
Δ
x
2
-
Δ
x
cos
⁡
k
i
Δ
x
2
-
Δ
x
sin
⁡
k
i
Δ
x
2
k
i
2
sin
⁡
k
i
Δ
x
2
-
k
i
2
cos
⁡
k
i
Δ
x
2
Δ
x
k
i
2
sin
⁡
k
i
Δ
x
2
-
cos
⁡
k
i
Δ
x
2
-
Δ
x
k
i
2
cos
⁡
k
i
Δ
x
2
-
sin
⁡
k
i
Δ
x
2
P
2
cos
⁡
k
i
Δ
x
2
P
2
sin
⁡
k
i
Δ
x
2
Δ
x
P
2
cos
⁡
k
i
Δ
x
2
+
2
E
i
I
i
k
i
sin
⁡
k
i
Δ
x
2
Δ
x
P
2
sin
⁡
k
i
Δ
x
2
-
2
E
i
I
i
k
i
cos
⁡
k
i
Δ
x
2
k
i
P
2
2
sin
⁡
k
i
Δ
x
2
-
k
i
P
2
2
cos
⁡
k
i
Δ
x
2
Δ
x
k
i
P
2
2
sin
⁡
k
i
Δ
x
2
+
P
2
cos
⁡
k
i
Δ
x
2
-
Δ
x
k
i
P
2
2
cos
⁡
k
i
Δ
x
2
+
P
2
sin
⁡
k
i
Δ
x
2
.Case 3 (ki4-4β^i2>0).
(36)
1
0
1
0
0
α
1
i
0
α
2
i
-
E
i
I
i
α
1
i
2
0
-
E
i
I
i
α
2
i
2
0
0
P
α
1
i
-
E
i
I
i
α
1
i
3
0
P
α
2
i
-
E
i
I
i
α
2
i
3
-
cos
⁡
α
1
i
Δ
x
-
sin
⁡
α
1
i
Δ
x
-
cos
⁡
α
2
i
Δ
x
-
sin
⁡
α
2
i
Δ
x
α
1
i
sin
⁡
α
1
i
Δ
x
-
α
1
i
cos
⁡
α
1
i
Δ
x
α
2
i
sin
⁡
α
2
i
Δ
x
-
α
2
i
cos
⁡
α
2
i
Δ
x
E
i
I
i
α
1
i
2
cos
⁡
α
1
i
Δ
x
E
i
I
i
α
1
i
2
sin
⁡
α
1
i
Δ
x
E
i
I
i
α
2
i
2
cos
⁡
α
2
i
Δ
x
E
i
I
i
α
2
i
2
sin
⁡
α
2
i
Δ
x
P
α
1
i
-
E
i
I
i
α
1
i
3
sin
⁡
α
1
i
Δ
x
-
P
α
1
i
-
E
i
I
i
α
1
i
3
cos
⁡
α
1
i
Δ
x
P
α
2
i
-
E
i
I
i
α
2
i
3
sin
⁡
α
2
i
Δ
x
-
P
α
2
i
-
E
i
I
i
α
2
i
3
cos
⁡
α
2
i
Δ
x
.Case 4 (ki4-4β^i2<0).
(37)
1
0
1
0
γ
1
i
g
i
γ
2
i
g
i
-
P
2
2
E
i
I
i
γ
1
i
g
i
-
P
2
2
E
i
I
i
γ
2
i
g
i
-
E
i
I
i
β
^
i
γ
1
i
E
i
I
i
β
^
i
g
i
-
E
i
I
i
β
^
i
γ
2
i
E
i
I
i
β
^
i
g
i
-
e
γ
1
i
Δ
x
cos
⁡
g
i
Δ
x
-
e
γ
1
i
Δ
x
sin
⁡
g
i
Δ
x
-
e
γ
2
i
Δ
x
cos
⁡
g
i
Δ
x
-
e
γ
2
i
Δ
x
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
g
i
sin
⁡
g
i
Δ
x
-
γ
1
i
cos
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
-
g
i
cos
⁡
g
i
Δ
x
-
γ
1
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
g
i
sin
⁡
g
i
Δ
x
-
γ
2
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
-
g
i
cos
⁡
g
i
Δ
x
-
γ
2
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
P
2
cos
⁡
g
i
Δ
x
+
2
γ
1
i
E
i
I
i
g
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
P
2
sin
⁡
g
i
Δ
x
-
2
γ
1
i
E
i
I
i
g
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
P
2
cos
⁡
g
i
Δ
x
+
2
γ
2
i
E
i
I
i
g
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
P
2
sin
⁡
g
i
Δ
x
-
2
γ
2
i
E
i
I
i
g
i
cos
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
E
i
I
i
β
^
i
γ
1
i
cos
⁡
g
i
Δ
x
+
E
i
I
i
β
^
i
g
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
E
i
I
i
β
^
i
γ
1
i
sin
⁡
g
i
Δ
x
-
E
i
I
i
β
^
i
g
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
E
i
I
i
β
^
i
γ
2
i
cos
⁡
g
i
Δ
x
+
E
i
I
i
β
^
i
g
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
E
i
I
i
β
^
i
γ
2
i
sin
⁡
g
i
Δ
x
-
E
i
I
i
β
^
i
g
i
cos
⁡
g
i
Δ
x
.(b) First SegmentCase 1 (β1=0).
(38)
β
A
0
P
β
A
-
P
-
K
A
k
1
-
K
A
0
-
cos
⁡
k
1
Δ
x
-
sin
⁡
k
1
Δ
x
-
Δ
x
-
1
k
1
sin
⁡
k
1
Δ
x
-
k
1
cos
⁡
k
1
Δ
x
-
1
0
P
cos
⁡
k
1
Δ
x
P
sin
⁡
k
1
Δ
x
0
0
0
0
-
P
0
.Case 2 (k14-4β^12=0).
(39)
β
A
k
1
P
2
2
-
P
2
0
-
P
2
-
K
A
k
1
2
-
K
A
2
E
1
I
1
k
1
-
cos
⁡
k
1
Δ
x
2
-
sin
⁡
k
1
Δ
x
2
-
Δ
x
cos
⁡
k
1
Δ
x
2
-
Δ
x
sin
⁡
k
1
Δ
x
2
k
1
2
sin
⁡
k
1
Δ
x
2
-
k
1
2
cos
⁡
k
1
Δ
x
2
Δ
x
k
1
2
sin
⁡
k
1
Δ
x
2
-
cos
⁡
k
1
Δ
x
2
-
Δ
x
k
1
2
cos
⁡
k
1
Δ
x
2
-
sin
⁡
k
1
Δ
x
2
P
2
cos
⁡
k
1
Δ
x
2
P
2
sin
⁡
k
1
Δ
x
2
Δ
x
P
2
cos
⁡
k
1
Δ
x
2
+
2
E
1
I
1
k
1
sin
⁡
k
1
Δ
x
2
Δ
x
P
2
sin
⁡
k
1
Δ
x
2
-
2
E
1
I
1
k
1
cos
⁡
k
1
Δ
x
2
k
1
P
2
2
sin
⁡
k
1
Δ
x
2
-
k
1
P
2
2
cos
⁡
k
1
Δ
x
2
Δ
x
k
1
P
2
2
sin
⁡
k
1
Δ
x
2
+
P
2
cos
⁡
k
1
Δ
x
2
-
Δ
x
k
1
P
2
2
cos
⁡
k
1
Δ
x
2
+
P
2
sin
⁡
k
1
Δ
x
2
.Case 3 (k14-4β^12>0).
(40)
β
A
P
α
11
-
E
1
I
1
α
11
3
β
A
P
α
21
-
E
1
I
1
α
21
3
-
E
1
I
1
α
11
2
-
K
A
α
11
-
E
1
I
1
α
21
2
-
K
A
α
21
-
cos
⁡
α
11
Δ
x
-
sin
⁡
α
11
Δ
x
-
cos
⁡
α
21
Δ
x
-
sin
⁡
α
21
Δ
x
α
11
sin
⁡
α
11
Δ
x
-
α
11
cos
⁡
α
11
Δ
x
α
21
sin
⁡
α
21
Δ
x
-
α
21
cos
⁡
α
21
Δ
x
E
1
I
1
α
11
2
cos
⁡
α
11
Δ
x
E
1
I
1
α
11
2
sin
⁡
α
11
Δ
x
E
1
I
1
α
21
2
cos
⁡
α
21
Δ
x
E
1
I
1
α
21
2
sin
⁡
α
21
Δ
x
P
α
11
-
E
1
I
1
α
11
3
sin
⁡
α
11
Δ
x
-
P
α
11
-
E
1
I
1
α
11
3
cos
⁡
α
11
Δ
x
P
α
21
-
E
1
I
1
α
21
3
sin
⁡
α
21
Δ
x
-
P
α
21
-
E
1
I
1
α
21
3
cos
⁡
α
21
Δ
x
.Case 4 (k14-4β^12<0).
(41)
β
A
-
E
1
I
1
β
^
1
γ
11
E
1
I
1
β
^
1
g
1
β
A
-
E
1
I
1
β
^
1
γ
21
E
1
I
1
β
^
1
g
1
-
K
A
γ
11
-
P
2
-
K
A
-
2
E
1
I
1
γ
11
g
1
-
K
A
γ
21
-
P
2
-
K
A
-
2
E
1
I
1
γ
21
g
1
-
e
γ
11
Δ
x
cos
⁡
g
1
Δ
x
-
e
γ
11
Δ
x
sin
⁡
g
1
Δ
x
-
e
γ
21
Δ
x
cos
⁡
g
1
Δ
x
-
e
γ
21
Δ
x
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
g
1
sin
⁡
g
1
Δ
x
-
γ
11
cos
⁡
g
1
Δ
x
e
γ
11
Δ
x
-
g
1
cos
⁡
g
1
Δ
x
-
γ
11
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
g
1
sin
⁡
g
1
Δ
x
-
γ
21
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
-
g
1
cos
⁡
g
1
Δ
x
-
γ
21
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
P
2
cos
⁡
g
1
Δ
x
+
2
γ
11
E
1
I
1
g
1
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
P
2
sin
⁡
g
1
Δ
x
-
2
γ
11
E
1
I
1
g
1
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
P
2
cos
⁡
g
1
Δ
x
+
2
γ
21
E
1
I
1
g
1
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
P
2
sin
⁡
g
1
Δ
x
-
2
γ
21
E
1
I
1
g
1
cos
⁡
g
1
Δ
x
e
γ
11
Δ
x
E
1
I
1
β
^
1
γ
11
cos
⁡
g
1
Δ
x
+
E
1
I
1
β
^
1
g
1
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
E
1
I
1
β
^
1
γ
11
sin
⁡
g
1
Δ
x
-
E
1
I
1
β
^
1
g
1
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
E
1
I
1
β
^
1
γ
21
cos
⁡
g
1
Δ
x
+
E
1
I
1
β
^
1
g
1
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
E
1
I
1
β
^
1
γ
21
sin
⁡
g
1
Δ
x
-
E
1
I
1
β
^
1
g
1
cos
⁡
g
1
Δ
x
.(c) Last SegmentCase 1 (βN=0).
(42)
1
0
0
1
0
k
N
1
0
-
P
0
0
0
0
0
P
0
-
β
B
cos
⁡
k
N
Δ
x
-
β
B
sin
⁡
k
N
Δ
x
P
-
β
B
Δ
x
-
β
B
-
P
cos
⁡
k
N
Δ
x
-
K
B
k
N
sin
⁡
k
N
Δ
x
-
P
sin
⁡
k
N
Δ
x
+
K
B
k
N
cos
⁡
k
N
Δ
x
K
B
0
.Case 2 (kN4-4β^N2=0).
(43)
1
0
0
0
0
k
N
2
1
0
-
P
2
0
0
2
E
N
I
N
k
N
0
k
N
P
2
2
-
P
2
0
-
β
B
cos
⁡
k
N
Δ
x
2
-
k
N
P
2
2
sin
⁡
k
N
Δ
x
2
-
β
B
sin
⁡
k
N
Δ
x
2
+
k
N
P
2
2
cos
⁡
k
N
Δ
x
2
-
β
B
Δ
x
+
P
2
cos
⁡
k
N
Δ
x
2
-
Δ
x
k
N
P
2
2
sin
⁡
k
N
Δ
x
2
-
β
B
Δ
x
+
P
2
sin
⁡
k
N
Δ
x
2
+
Δ
x
k
N
P
2
2
cos
⁡
k
N
Δ
x
2
-
P
2
cos
⁡
k
N
Δ
x
2
-
k
N
2
K
B
sin
⁡
k
N
Δ
x
2
-
P
2
sin
⁡
k
N
Δ
x
2
+
k
N
2
K
B
cos
⁡
k
N
Δ
x
2
-
Δ
x
P
2
-
K
B
cos
⁡
k
N
Δ
x
2
-
2
E
N
I
N
k
N
+
K
B
Δ
x
k
N
2
sin
⁡
k
N
Δ
x
2
-
Δ
x
P
2
-
K
B
sin
⁡
k
N
Δ
x
2
+
2
E
N
I
N
k
N
+
K
B
Δ
x
k
N
2
cos
⁡
k
N
Δ
x
2
.Case 3 (kN4-4β^N2>0).
(44)
1
0
1
0
0
α
1
N
0
α
2
N
-
E
N
I
N
α
1
N
2
0
-
E
N
I
N
α
2
N
2
0
0
P
α
1
N
-
E
N
I
N
α
1
N
3
0
P
α
2
N
-
E
N
I
N
α
2
N
3
-
β
B
cos
⁡
α
1
N
Δ
x
-
P
α
1
N
-
E
N
I
N
α
1
N
3
sin
⁡
α
1
N
Δ
x
-
β
B
sin
⁡
α
1
N
Δ
x
+
P
α
1
N
-
E
N
I
N
α
1
N
3
cos
⁡
α
1
N
Δ
x
-
β
B
cos
⁡
α
2
N
Δ
x
-
P
α
2
N
-
E
N
I
N
α
2
N
3
sin
⁡
α
2
N
Δ
x
-
β
B
sin
⁡
α
2
N
Δ
x
+
P
α
2
N
-
E
N
I
N
α
2
N
3
cos
⁡
α
2
N
Δ
x
-
K
B
α
1
N
sin
⁡
α
1
N
Δ
x
-
E
N
I
N
α
1
N
2
cos
⁡
α
1
N
Δ
x
K
B
α
1
N
cos
⁡
α
1
N
Δ
x
-
E
N
I
N
α
1
N
2
sin
⁡
α
1
N
Δ
x
-
K
B
α
2
N
sin
⁡
α
2
N
Δ
x
-
E
N
I
N
α
2
N
2
cos
⁡
α
2
N
Δ
x
K
B
α
2
N
cos
⁡
α
2
N
Δ
x
-
E
N
I
N
α
2
N
2
sin
⁡
α
2
N
Δ
x
.Case 4 (kN4-4β^N2<0).
(45)
1
0
1
0
γ
1
N
g
N
γ
2
N
g
N
-
P
2
2
E
N
I
N
γ
1
N
g
N
-
P
2
2
E
N
I
N
γ
2
N
g
N
-
E
N
I
N
β
^
N
γ
1
N
E
N
I
N
β
^
N
g
N
-
E
N
I
N
β
^
N
γ
2
N
E
N
I
N
β
^
N
g
N
-
e
γ
1
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
1
N
cos
⁡
g
N
Δ
x
+
E
N
I
N
β
^
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
1
N
sin
⁡
g
N
Δ
x
-
E
N
I
N
β
^
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
2
N
cos
⁡
g
N
Δ
x
+
E
N
I
N
β
^
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
2
N
sin
⁡
g
N
Δ
x
-
E
N
I
N
β
^
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
P
2
-
K
B
γ
1
N
cos
⁡
g
N
Δ
x
+
K
B
+
2
E
N
I
N
γ
1
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
P
2
-
K
B
γ
1
N
sin
⁡
g
N
Δ
x
-
K
B
+
2
E
N
I
N
γ
1
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
P
2
-
K
B
γ
2
N
cos
⁡
g
N
Δ
x
+
K
B
+
2
E
N
I
N
γ
2
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
P
2
-
K
B
γ
2
N
sin
⁡
g
N
Δ
x
-
K
B
+
2
E
N
I
N
γ
2
N
g
N
cos
⁡
g
N
Δ
x
.


### 2.4. The Mode Shape

The mode shape can be obtained by substituting the buckling load into theP-matrix and setting one of the constant coefficients to be equal to 1. The original system looks like the following:(46)a11a12……………………a14Na21⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮P-Matrix⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮a4N1………………………a4N4NA1B1C1D1⋮AiBiCiDi⋮ANBNCNDN=0000⋮0000⋮0000.After lettingDN=1, the following system of equations is obtained. (47)a11a12…………a14N-1a21⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮Reduced  P-Matrix⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮a4N-11……………a4N-14N-1A1B1C1D1⋮AiBiCiDi⋮ANBNCN=-a14N-a24N⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮-a4N-14N.The constant coefficients can then be calculated by solving the system as follows:(48)CR=PR-1F,where PR is the reduced P-matrix, CR=A1B1C1D1⋯AiBiCiDi⋯ANBNCNT, and F=-a14N-a24N⋯⋯⋯-a4N-14NT.After substituting the constants back into their respective equations, one can obtain the deflection functions and plot the mode shape.


## 2.1. The Buckling Equation

The buckling differential equation for the column shown in Figure1 can be derived using the bifurcation method resulting in the following equation:(1)EI∂4v∂x4+2∂EI∂x∂3v∂x3+∂2EI∂x2+P∂2v∂x2+βv=0.Figure 1
Nonprismatic columns with arbitrary boundary conditions resting on nonuniform elastic foundation.The boundary conditions (BCs) are as follows:(1) Shear
(2)
V
0
=
-
β
A
v
0
,
V
L
=
β
B
v
L
.(2) Moment
(3)
M
0
=
K
A
v
′
0
,
M
L
=
-
K
B
v
′
L
,where Ex is the column elastic modulus function, Ix is the column moment of inertia function, P is the column axial load, βx is the elastic foundation spring stiffness function, vx is the column lateral deflection function, x is the distance from the bottom of the column, βA is the bottom lateral discrete spring stiffness, βB is the top lateral discrete spring stiffness, KA is the bottom rotational discrete spring stiffness, KB is the top rotational discrete spring stiffness, L is the length of the column, Vx is the column shear force function, and Mx is the column moment function.Such a differential equation can prove to be extremely difficult to solve analytically for arbitraryE, I, and β functions. The authors are not aware of any analytical mathematical method that can solve this differential equation. Thus, the authors decided to solve the differential equation semianalytically by discretizing the column into a finite number of segments. Each segment has constants E, I, and β, as shown in Figure 2. E, I, and β values for each segment were taken to be the average of each function within that segment of the column. The averages were calculated using the general equation for averaging functions within an interval:(4)fave=∫abfxdxb-a.Figure 2
Idealized prismatic column segment resting on uniform elastic foundation.An exact solution was then obtained for each segment using the classical buckling equation of columns on an elastic foundation. The domain ofx for each segment was taken to be 0≤x≤Δx:(5)∂4vi∂x4+ki2∂2vi∂x2+β^i2vi=0,0≤x≤Δx,ki=PEiIi,β^i=βiEiIi,Δx=LN,where N is number of segments, Δx is segment length, and i is segment number.


## 2.2. Analytical Solution of the Segment Buckling Equation

The general solution to the segment differential equation with arbitrary boundary conditions is complicated and requires the consideration of four different cases. The solution to each case is shown below.Case 1 (βi=0).
The segment differential equation becomes(6)∂4vi∂x4+ki2∂2vi∂x2=0,0≤x≤Δx.
Letvi=eλix, which can be substituted into (6) resulting in(7)λi4+ki2λi2eλix=0→yieldsλi1,2=0,λi3,4=±jki,where(8)j=-1.
The deflection function is then obtained as follows:(9)vix=Aicos⁡kix+Bisin⁡kix+Cix+Di.Case 2 (ki4-4β^i2=0).
Letvi=eλix, which can be substituted into (5) resulting in(10)λi4+ki2λi2+β^i2eλix=0→yieldsλi1,2,3,4=±jki2.
The deflection function is then found to be(11)vix=Aicos⁡kix2+Bisin⁡kix2+Cixcos⁡kix2+Dixsin⁡kix2.Case 3 (ki4-4β^i2>0).
Letvi=eλix, which can be substituted into (5) resulting in(12)λi4+ki2λi2+β^i2eλix=0→yieldsλi2=-ki2±ki4-4β^i22→yieldsλi1,2=±jα1i,λi3,4=±jα2i,where(13)α1i=ki2+ki4-4β^i22,α2i=ki2-ki4-4β^i22.
The deflection function is then obtained to be as follows(14)vix=Aicos⁡α1ix+Bisin⁡α1ix+Cicos⁡α2ix+Disin⁡α2ix.Case 4 (ki4-4β^i2<0).
Letvi=eλix, which can be substituted into (5) resulting in(15)λi4+ki2λi2+β^i2eλix=0→yieldsλi2=-ki2±ki4-4β^i22=-ki22±j4β^i2-ki42.Letλi=ai+jbi implying λi2=ai2-bi2+j2aibi which results in the following two equations with two unknowns, aiandbi:(16)ai2-bi2=-ki22,(17)2aibi=±4β^i2-ki42.From (17) we get(18)ai=±4β^i2-ki44bi→yieldsai2=4β^i2-ki416bi2.Substituting into (16) we obtain(19)4β^i2-ki416bi2-bi2=-ki22→yields16bi4-8ki2bi2+ki4-4β^i2=0→yieldsbi2=ki2±2β^i4.We know thatki4-4β^i2<0→yields(ki2+2β^i)(ki2-2β^i)<0, but it is obvious that ki2+2β^i>0. Therefore, ki2-2β^i<0→yieldsbi2≠ki2-2β^i/4 since bi2>0. This leaves the following: bi2=ki2+2β^i/4→yieldsbi=±2β^i+ki2/2.Using (17) we can then obtain ai as follows:(20)ai=±2β^i-ki22.Therefore,λi is found to be(21)λi=±2β^i-ki22±j2β^i+ki22→yieldsλi1,2=γ1i±jgi,λi3,4=γ2i±jgi,where(22)gi=2β^i+ki22,γ1i=2β^i-ki22,γ2i=-2β^i-ki22.The deflection function is then obtained to be(23)vix=Aieγ1ixcosgix+Bieγ1ixsingix+Cieγ2ixcosgix+Dieγ2ixsingix,where Ai,Bi,Ci,andDi are the constant coefficients for each segment.After solving the segment differential equation, the following compatibility conditions were established at each node between the segments.(1) Displacement
(24)
v
i
0
=
v
i
-
1
Δ
x
,
v
i
Δ
x
=
v
i
+
1
0
.(2) Slope
(25)
v
i
′
0
=
v
i
-
1
′
Δ
x
,
v
i
′
Δ
x
=
v
i
+
1
′
0
.(3) Moment
(26)
E
i
I
i
v
i
′
′
0
=
E
i
-
1
I
i
-
1
v
i
-
1
′
′
Δ
x
,
E
i
I
i
v
i
′
′
Δ
x
=
E
i
+
1
I
i
+
1
v
i
+
1
′
′
0
.(4) Shear
(27)
E
i
I
i
v
i
′
′
′
0
+
P
v
i
′
0
=
E
i
-
1
I
i
-
1
v
i
-
1
′
′
′
Δ
x
+
P
v
i
-
1
′
Δ
x
,
E
i
I
i
v
i
′
′
′
Δ
x
+
P
v
i
′
Δ
x
=
E
i
+
1
I
i
+
1
v
i
+
1
′
′
′
0
+
P
v
i
+
1
′
0
.Finally, in order to obtain the coefficient matrix (P-matrix), the following boundary conditions were established.(1) Shear
(28)
E
1
I
1
v
1
′
′
′
0
+
P
v
1
′
0
=
-
β
A
v
1
0
,
E
N
I
N
v
N
′
′
′
Δ
x
+
P
v
N
′
Δ
x
=
β
B
v
N
Δ
x
.(2) Moment
(29)
E
1
I
1
v
1
′
′
0
=
K
A
v
1
′
0
,
E
N
I
N
v
N
′
′
Δ
x
=
-
K
B
v
N
′
Δ
x
.Solving the following system of equations would then render the buckling load:(30)PC=0,where P is the P-Matrix and C = A1  B1C1D1⋯AiBiCiDi⋯ANBNCNDNT.


## 2.3. TheP-Matrix

The structure of theP-matrix is as follows.(1) One Segment
(31)
4
×
4
B
C
M
a
t
r
i
x
.(2) Two Segments
(32)
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
4
×
4
F
i
r
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
4
×
4
L
a
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
.(3) Three Segments
(33)
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
4
×
4
F
i
r
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
8
×
4
I
n
t
e
r
m
e
d
i
a
t
e
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
4
×
4
L
a
s
t
S
e
g
m
e
n
t
C
o
m
p
a
t
i
b
i
l
i
t
y
S
u
b
-
M
a
t
r
i
x
2
×
4
B
C
S
u
b
-
M
a
t
r
i
x
.For the case of an assembly with more than three segments, the structure of theP-matrix would be the same as the three segments case but with additional intermediate segment compatibility submatrices that are shifted diagonally downwards. The buckling load would then be the load that produces a singular P-matrix.The submatrices inside theP-matrix are obtained by applying the compatibility and boundary conditions and are dependent on the case that governs each segment, as explained earlier. These submatrices were found to be as shown below.(a) Intermediate SegmentsCase 1 (βi=0).
(34)
1
0
0
1
0
k
i
1
0
-
P
0
0
0
0
0
P
0
-
cos
⁡
k
i
Δ
x
-
sin
⁡
k
i
Δ
x
-
Δ
x
-
1
k
i
sin
⁡
k
i
Δ
x
-
k
i
cos
⁡
k
i
Δ
x
-
1
0
P
cos
⁡
k
i
Δ
x
P
sin
⁡
k
i
Δ
x
0
0
0
0
-
P
0
.Case 2 (ki4-4β^i2=0).
(35)
1
0
0
0
0
k
i
2
1
0
-
P
2
0
0
2
E
i
I
i
k
i
0
k
i
P
2
2
-
P
2
0
-
cos
⁡
k
i
Δ
x
2
-
sin
⁡
k
i
Δ
x
2
-
Δ
x
cos
⁡
k
i
Δ
x
2
-
Δ
x
sin
⁡
k
i
Δ
x
2
k
i
2
sin
⁡
k
i
Δ
x
2
-
k
i
2
cos
⁡
k
i
Δ
x
2
Δ
x
k
i
2
sin
⁡
k
i
Δ
x
2
-
cos
⁡
k
i
Δ
x
2
-
Δ
x
k
i
2
cos
⁡
k
i
Δ
x
2
-
sin
⁡
k
i
Δ
x
2
P
2
cos
⁡
k
i
Δ
x
2
P
2
sin
⁡
k
i
Δ
x
2
Δ
x
P
2
cos
⁡
k
i
Δ
x
2
+
2
E
i
I
i
k
i
sin
⁡
k
i
Δ
x
2
Δ
x
P
2
sin
⁡
k
i
Δ
x
2
-
2
E
i
I
i
k
i
cos
⁡
k
i
Δ
x
2
k
i
P
2
2
sin
⁡
k
i
Δ
x
2
-
k
i
P
2
2
cos
⁡
k
i
Δ
x
2
Δ
x
k
i
P
2
2
sin
⁡
k
i
Δ
x
2
+
P
2
cos
⁡
k
i
Δ
x
2
-
Δ
x
k
i
P
2
2
cos
⁡
k
i
Δ
x
2
+
P
2
sin
⁡
k
i
Δ
x
2
.Case 3 (ki4-4β^i2>0).
(36)
1
0
1
0
0
α
1
i
0
α
2
i
-
E
i
I
i
α
1
i
2
0
-
E
i
I
i
α
2
i
2
0
0
P
α
1
i
-
E
i
I
i
α
1
i
3
0
P
α
2
i
-
E
i
I
i
α
2
i
3
-
cos
⁡
α
1
i
Δ
x
-
sin
⁡
α
1
i
Δ
x
-
cos
⁡
α
2
i
Δ
x
-
sin
⁡
α
2
i
Δ
x
α
1
i
sin
⁡
α
1
i
Δ
x
-
α
1
i
cos
⁡
α
1
i
Δ
x
α
2
i
sin
⁡
α
2
i
Δ
x
-
α
2
i
cos
⁡
α
2
i
Δ
x
E
i
I
i
α
1
i
2
cos
⁡
α
1
i
Δ
x
E
i
I
i
α
1
i
2
sin
⁡
α
1
i
Δ
x
E
i
I
i
α
2
i
2
cos
⁡
α
2
i
Δ
x
E
i
I
i
α
2
i
2
sin
⁡
α
2
i
Δ
x
P
α
1
i
-
E
i
I
i
α
1
i
3
sin
⁡
α
1
i
Δ
x
-
P
α
1
i
-
E
i
I
i
α
1
i
3
cos
⁡
α
1
i
Δ
x
P
α
2
i
-
E
i
I
i
α
2
i
3
sin
⁡
α
2
i
Δ
x
-
P
α
2
i
-
E
i
I
i
α
2
i
3
cos
⁡
α
2
i
Δ
x
.Case 4 (ki4-4β^i2<0).
(37)
1
0
1
0
γ
1
i
g
i
γ
2
i
g
i
-
P
2
2
E
i
I
i
γ
1
i
g
i
-
P
2
2
E
i
I
i
γ
2
i
g
i
-
E
i
I
i
β
^
i
γ
1
i
E
i
I
i
β
^
i
g
i
-
E
i
I
i
β
^
i
γ
2
i
E
i
I
i
β
^
i
g
i
-
e
γ
1
i
Δ
x
cos
⁡
g
i
Δ
x
-
e
γ
1
i
Δ
x
sin
⁡
g
i
Δ
x
-
e
γ
2
i
Δ
x
cos
⁡
g
i
Δ
x
-
e
γ
2
i
Δ
x
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
g
i
sin
⁡
g
i
Δ
x
-
γ
1
i
cos
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
-
g
i
cos
⁡
g
i
Δ
x
-
γ
1
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
g
i
sin
⁡
g
i
Δ
x
-
γ
2
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
-
g
i
cos
⁡
g
i
Δ
x
-
γ
2
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
P
2
cos
⁡
g
i
Δ
x
+
2
γ
1
i
E
i
I
i
g
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
P
2
sin
⁡
g
i
Δ
x
-
2
γ
1
i
E
i
I
i
g
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
P
2
cos
⁡
g
i
Δ
x
+
2
γ
2
i
E
i
I
i
g
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
P
2
sin
⁡
g
i
Δ
x
-
2
γ
2
i
E
i
I
i
g
i
cos
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
E
i
I
i
β
^
i
γ
1
i
cos
⁡
g
i
Δ
x
+
E
i
I
i
β
^
i
g
i
sin
⁡
g
i
Δ
x
e
γ
1
i
Δ
x
E
i
I
i
β
^
i
γ
1
i
sin
⁡
g
i
Δ
x
-
E
i
I
i
β
^
i
g
i
cos
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
E
i
I
i
β
^
i
γ
2
i
cos
⁡
g
i
Δ
x
+
E
i
I
i
β
^
i
g
i
sin
⁡
g
i
Δ
x
e
γ
2
i
Δ
x
E
i
I
i
β
^
i
γ
2
i
sin
⁡
g
i
Δ
x
-
E
i
I
i
β
^
i
g
i
cos
⁡
g
i
Δ
x
.(b) First SegmentCase 1 (β1=0).
(38)
β
A
0
P
β
A
-
P
-
K
A
k
1
-
K
A
0
-
cos
⁡
k
1
Δ
x
-
sin
⁡
k
1
Δ
x
-
Δ
x
-
1
k
1
sin
⁡
k
1
Δ
x
-
k
1
cos
⁡
k
1
Δ
x
-
1
0
P
cos
⁡
k
1
Δ
x
P
sin
⁡
k
1
Δ
x
0
0
0
0
-
P
0
.Case 2 (k14-4β^12=0).
(39)
β
A
k
1
P
2
2
-
P
2
0
-
P
2
-
K
A
k
1
2
-
K
A
2
E
1
I
1
k
1
-
cos
⁡
k
1
Δ
x
2
-
sin
⁡
k
1
Δ
x
2
-
Δ
x
cos
⁡
k
1
Δ
x
2
-
Δ
x
sin
⁡
k
1
Δ
x
2
k
1
2
sin
⁡
k
1
Δ
x
2
-
k
1
2
cos
⁡
k
1
Δ
x
2
Δ
x
k
1
2
sin
⁡
k
1
Δ
x
2
-
cos
⁡
k
1
Δ
x
2
-
Δ
x
k
1
2
cos
⁡
k
1
Δ
x
2
-
sin
⁡
k
1
Δ
x
2
P
2
cos
⁡
k
1
Δ
x
2
P
2
sin
⁡
k
1
Δ
x
2
Δ
x
P
2
cos
⁡
k
1
Δ
x
2
+
2
E
1
I
1
k
1
sin
⁡
k
1
Δ
x
2
Δ
x
P
2
sin
⁡
k
1
Δ
x
2
-
2
E
1
I
1
k
1
cos
⁡
k
1
Δ
x
2
k
1
P
2
2
sin
⁡
k
1
Δ
x
2
-
k
1
P
2
2
cos
⁡
k
1
Δ
x
2
Δ
x
k
1
P
2
2
sin
⁡
k
1
Δ
x
2
+
P
2
cos
⁡
k
1
Δ
x
2
-
Δ
x
k
1
P
2
2
cos
⁡
k
1
Δ
x
2
+
P
2
sin
⁡
k
1
Δ
x
2
.Case 3 (k14-4β^12>0).
(40)
β
A
P
α
11
-
E
1
I
1
α
11
3
β
A
P
α
21
-
E
1
I
1
α
21
3
-
E
1
I
1
α
11
2
-
K
A
α
11
-
E
1
I
1
α
21
2
-
K
A
α
21
-
cos
⁡
α
11
Δ
x
-
sin
⁡
α
11
Δ
x
-
cos
⁡
α
21
Δ
x
-
sin
⁡
α
21
Δ
x
α
11
sin
⁡
α
11
Δ
x
-
α
11
cos
⁡
α
11
Δ
x
α
21
sin
⁡
α
21
Δ
x
-
α
21
cos
⁡
α
21
Δ
x
E
1
I
1
α
11
2
cos
⁡
α
11
Δ
x
E
1
I
1
α
11
2
sin
⁡
α
11
Δ
x
E
1
I
1
α
21
2
cos
⁡
α
21
Δ
x
E
1
I
1
α
21
2
sin
⁡
α
21
Δ
x
P
α
11
-
E
1
I
1
α
11
3
sin
⁡
α
11
Δ
x
-
P
α
11
-
E
1
I
1
α
11
3
cos
⁡
α
11
Δ
x
P
α
21
-
E
1
I
1
α
21
3
sin
⁡
α
21
Δ
x
-
P
α
21
-
E
1
I
1
α
21
3
cos
⁡
α
21
Δ
x
.Case 4 (k14-4β^12<0).
(41)
β
A
-
E
1
I
1
β
^
1
γ
11
E
1
I
1
β
^
1
g
1
β
A
-
E
1
I
1
β
^
1
γ
21
E
1
I
1
β
^
1
g
1
-
K
A
γ
11
-
P
2
-
K
A
-
2
E
1
I
1
γ
11
g
1
-
K
A
γ
21
-
P
2
-
K
A
-
2
E
1
I
1
γ
21
g
1
-
e
γ
11
Δ
x
cos
⁡
g
1
Δ
x
-
e
γ
11
Δ
x
sin
⁡
g
1
Δ
x
-
e
γ
21
Δ
x
cos
⁡
g
1
Δ
x
-
e
γ
21
Δ
x
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
g
1
sin
⁡
g
1
Δ
x
-
γ
11
cos
⁡
g
1
Δ
x
e
γ
11
Δ
x
-
g
1
cos
⁡
g
1
Δ
x
-
γ
11
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
g
1
sin
⁡
g
1
Δ
x
-
γ
21
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
-
g
1
cos
⁡
g
1
Δ
x
-
γ
21
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
P
2
cos
⁡
g
1
Δ
x
+
2
γ
11
E
1
I
1
g
1
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
P
2
sin
⁡
g
1
Δ
x
-
2
γ
11
E
1
I
1
g
1
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
P
2
cos
⁡
g
1
Δ
x
+
2
γ
21
E
1
I
1
g
1
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
P
2
sin
⁡
g
1
Δ
x
-
2
γ
21
E
1
I
1
g
1
cos
⁡
g
1
Δ
x
e
γ
11
Δ
x
E
1
I
1
β
^
1
γ
11
cos
⁡
g
1
Δ
x
+
E
1
I
1
β
^
1
g
1
sin
⁡
g
1
Δ
x
e
γ
11
Δ
x
E
1
I
1
β
^
1
γ
11
sin
⁡
g
1
Δ
x
-
E
1
I
1
β
^
1
g
1
cos
⁡
g
1
Δ
x
e
γ
21
Δ
x
E
1
I
1
β
^
1
γ
21
cos
⁡
g
1
Δ
x
+
E
1
I
1
β
^
1
g
1
sin
⁡
g
1
Δ
x
e
γ
21
Δ
x
E
1
I
1
β
^
1
γ
21
sin
⁡
g
1
Δ
x
-
E
1
I
1
β
^
1
g
1
cos
⁡
g
1
Δ
x
.(c) Last SegmentCase 1 (βN=0).
(42)
1
0
0
1
0
k
N
1
0
-
P
0
0
0
0
0
P
0
-
β
B
cos
⁡
k
N
Δ
x
-
β
B
sin
⁡
k
N
Δ
x
P
-
β
B
Δ
x
-
β
B
-
P
cos
⁡
k
N
Δ
x
-
K
B
k
N
sin
⁡
k
N
Δ
x
-
P
sin
⁡
k
N
Δ
x
+
K
B
k
N
cos
⁡
k
N
Δ
x
K
B
0
.Case 2 (kN4-4β^N2=0).
(43)
1
0
0
0
0
k
N
2
1
0
-
P
2
0
0
2
E
N
I
N
k
N
0
k
N
P
2
2
-
P
2
0
-
β
B
cos
⁡
k
N
Δ
x
2
-
k
N
P
2
2
sin
⁡
k
N
Δ
x
2
-
β
B
sin
⁡
k
N
Δ
x
2
+
k
N
P
2
2
cos
⁡
k
N
Δ
x
2
-
β
B
Δ
x
+
P
2
cos
⁡
k
N
Δ
x
2
-
Δ
x
k
N
P
2
2
sin
⁡
k
N
Δ
x
2
-
β
B
Δ
x
+
P
2
sin
⁡
k
N
Δ
x
2
+
Δ
x
k
N
P
2
2
cos
⁡
k
N
Δ
x
2
-
P
2
cos
⁡
k
N
Δ
x
2
-
k
N
2
K
B
sin
⁡
k
N
Δ
x
2
-
P
2
sin
⁡
k
N
Δ
x
2
+
k
N
2
K
B
cos
⁡
k
N
Δ
x
2
-
Δ
x
P
2
-
K
B
cos
⁡
k
N
Δ
x
2
-
2
E
N
I
N
k
N
+
K
B
Δ
x
k
N
2
sin
⁡
k
N
Δ
x
2
-
Δ
x
P
2
-
K
B
sin
⁡
k
N
Δ
x
2
+
2
E
N
I
N
k
N
+
K
B
Δ
x
k
N
2
cos
⁡
k
N
Δ
x
2
.Case 3 (kN4-4β^N2>0).
(44)
1
0
1
0
0
α
1
N
0
α
2
N
-
E
N
I
N
α
1
N
2
0
-
E
N
I
N
α
2
N
2
0
0
P
α
1
N
-
E
N
I
N
α
1
N
3
0
P
α
2
N
-
E
N
I
N
α
2
N
3
-
β
B
cos
⁡
α
1
N
Δ
x
-
P
α
1
N
-
E
N
I
N
α
1
N
3
sin
⁡
α
1
N
Δ
x
-
β
B
sin
⁡
α
1
N
Δ
x
+
P
α
1
N
-
E
N
I
N
α
1
N
3
cos
⁡
α
1
N
Δ
x
-
β
B
cos
⁡
α
2
N
Δ
x
-
P
α
2
N
-
E
N
I
N
α
2
N
3
sin
⁡
α
2
N
Δ
x
-
β
B
sin
⁡
α
2
N
Δ
x
+
P
α
2
N
-
E
N
I
N
α
2
N
3
cos
⁡
α
2
N
Δ
x
-
K
B
α
1
N
sin
⁡
α
1
N
Δ
x
-
E
N
I
N
α
1
N
2
cos
⁡
α
1
N
Δ
x
K
B
α
1
N
cos
⁡
α
1
N
Δ
x
-
E
N
I
N
α
1
N
2
sin
⁡
α
1
N
Δ
x
-
K
B
α
2
N
sin
⁡
α
2
N
Δ
x
-
E
N
I
N
α
2
N
2
cos
⁡
α
2
N
Δ
x
K
B
α
2
N
cos
⁡
α
2
N
Δ
x
-
E
N
I
N
α
2
N
2
sin
⁡
α
2
N
Δ
x
.Case 4 (kN4-4β^N2<0).
(45)
1
0
1
0
γ
1
N
g
N
γ
2
N
g
N
-
P
2
2
E
N
I
N
γ
1
N
g
N
-
P
2
2
E
N
I
N
γ
2
N
g
N
-
E
N
I
N
β
^
N
γ
1
N
E
N
I
N
β
^
N
g
N
-
E
N
I
N
β
^
N
γ
2
N
E
N
I
N
β
^
N
g
N
-
e
γ
1
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
1
N
cos
⁡
g
N
Δ
x
+
E
N
I
N
β
^
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
1
N
sin
⁡
g
N
Δ
x
-
E
N
I
N
β
^
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
2
N
cos
⁡
g
N
Δ
x
+
E
N
I
N
β
^
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
β
B
+
E
N
I
N
β
^
N
γ
2
N
sin
⁡
g
N
Δ
x
-
E
N
I
N
β
^
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
P
2
-
K
B
γ
1
N
cos
⁡
g
N
Δ
x
+
K
B
+
2
E
N
I
N
γ
1
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
1
N
Δ
x
P
2
-
K
B
γ
1
N
sin
⁡
g
N
Δ
x
-
K
B
+
2
E
N
I
N
γ
1
N
g
N
cos
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
P
2
-
K
B
γ
2
N
cos
⁡
g
N
Δ
x
+
K
B
+
2
E
N
I
N
γ
2
N
g
N
sin
⁡
g
N
Δ
x
-
e
γ
2
N
Δ
x
P
2
-
K
B
γ
2
N
sin
⁡
g
N
Δ
x
-
K
B
+
2
E
N
I
N
γ
2
N
g
N
cos
⁡
g
N
Δ
x
.


## 2.4. The Mode Shape

The mode shape can be obtained by substituting the buckling load into theP-matrix and setting one of the constant coefficients to be equal to 1. The original system looks like the following:(46)a11a12……………………a14Na21⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮P-Matrix⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮a4N1………………………a4N4NA1B1C1D1⋮AiBiCiDi⋮ANBNCNDN=0000⋮0000⋮0000.After lettingDN=1, the following system of equations is obtained. (47)a11a12…………a14N-1a21⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮Reduced  P-Matrix⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮a4N-11……………a4N-14N-1A1B1C1D1⋮AiBiCiDi⋮ANBNCN=-a14N-a24N⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮-a4N-14N.The constant coefficients can then be calculated by solving the system as follows:(48)CR=PR-1F,where PR is the reduced P-matrix, CR=A1B1C1D1⋯AiBiCiDi⋯ANBNCNT, and F=-a14N-a24N⋯⋯⋯-a4N-14NT.After substituting the constants back into their respective equations, one can obtain the deflection functions and plot the mode shape.


## 3. Sensitivity Analysis

As the proposed approach solves the differential equation for buckling numerically by discretizing the column into a finite number of segments, a sensitivity analysis was conducted to determine the optimum number of segments necessary to obtain accurate results. The sensitivity analysis was conducted for a simply supported column with both cross-section and soil stiffness varying to account for all approximations involved in the solution. The following functions were used:EI=x and  β=x2. The first 3 buckling loads and their associated mode shapes were determined. Table 1 shows the results for 10, 20, 50, 100, and 200 segments.Table 1
Sensitivity analysis results.
N
P
c
r
1
P
c
r
2
P
c
r
3
10
3.598021035
11.85057958
26.03460752
20
3.655748351
12.00379038
25.20788505
50
3.68272198
12.21549833
25.50181667
100
3.688362459
12.2799962
25.72594058
200
3.690129118
12.30323418
25.82633913The results, for the discretization segments above, were reasonably close; however, the results required a high segment count to stabilize. In the next section, it will be shown that this variation is not significant as the obtained results with a segment count of 20 provided excellent agreement with the results in the literature. Figure3 shows the first three mode shapes, which did not vary with the change in the number of segments.Figure 3
Mode shapes for column utilized in sensitivity analysis.


## 4. Validation of the Proposed Approach

In order to validate the proposed approach, a few cases were analyzed and compared with available solutions from the literature. For the first case, a column resting on elastic foundation is analyzed. A fixed-pinned column was evaluated with the following parameters:EI=1, L=1, and β=100 [10]. In this case, a single element was used in the analysis as no variation in any of the parameters occurs along the span. The critical load was determined to be 28.307, which matches the load obtained by Wang and his coauthors [12]. Figure 4 shows the first three mode shapes for this case and their associated buckling loads.Figure 4
Mode shapes for Case  1.The following cases are those for a nonprismatic column, as derived by Gere and Carter [13]. For Case  2, nonprismatic fixed-pinned column, as derived by Gere and Carter [13], was analyzed. The parameters were as follows: EI=1+3.5x/L3, L=1, and β=0, and the number of elements was taken to be 20. The buckling load obtained was 221.031, compared to 222.2 obtained by Gere’s approach. Figure 5 shows the mode shape obtained for this case.Figure 5
Mode shape for Case  2.In Case  3, the variable end conditions implementation was validated. This was accomplished by analyzing cases with different rotational spring values at the ends and comparing the results with the alignment charts. This was done for both the sway permitted and sway not permitted cases. The column analyzed had the following parameters:E=I=L=1. For the pinned-pinned reference case, the Euler buckling load (PE) was obtained to be 9.8696. The results obtained are shown in Table 2.Table 2
Results for Case  3.
k
a
k
b
Sway not permitted
Sway permitted
G
a
G
b
P
c
r
/
P
E (program)
P
c
r
/
P
E (charts)
% Error
G
a
G
b
P
c
r
/
P
E (program)
P
c
r
/
P
E (charts)
% error
40
20
0.05
0.1
3.471079
3.467779
0.095%
0.15
0.3
0.86597
0.865333
0.074%
10
10
0.2
0.2
2.853986
2.85336
0.022%
0.6
0.6
0.69959
0.700268
-
0.097
%
2
2
1
1
1.668095
1.669237
-
0.068
%
3
3
0.29998
0.300244
-
0.087
%
0.4
0.1
5
20
1.097091
1.096461
0.057%
15
60
0.04668
0.046709
-
0.057
%The results obtained showed excellent agreement with those obtained from the alignment charts. This validates the derivation of variable end supports and its implementation in the computer program.Finally, in Cases  4 and 5, the approach was validated for cases with nonprismatic sections [2] and variable soil stiffness [14], respectively. In Case  4, columns with different boundary conditions were analyzed as shown in Table 3. The parameters for these columns were as follows: EI=1+x/L, L=1, and β=0. The results obtained are also shown in Table 3.Table 3
Results for Case  4.
Boundary conditions
Critical load
Strong end
Weak end
Eisenberger [2]
Proposed approach
Free
Fixed
3.1177
3.118548250634
Fixed
Free
4.1242
4.1238
Hinged
Hinged
14.5112
14.5105
Hinged
Fixed
29.449
29.4571
Fixed
Hinged
29.4788
29.4734
Fixed
Fixed
57.394
57.4019For Case  5, a column resting on a variable elastic foundation was analyzed. The column had the following parameters:EI=1.5×105 N·m2, L=3 m, and β=4x-3x2+x3106N/m2. The critical load was obtained to be 148514.92 N, which is very close to what was obtained by Eisenberger and Clastornik [14].


## 5. Examples

### 5.1. Example  1: Column on Variable (Linear and Nonlinear) Elastic Foundation

After validating the proposed formulation in the previous section, the case of nonuniform soil pressure was investigated. A simply supported column was analyzed with 20 segments. The column’s parameters were as follows:EI=100 and L=10, and, similar to what Gere and Carter [13] did for the moment of inertia in nonprismatic columns, the soil stiffness (β) was taken to be as follows:(49)β=1+xLn,where n=1,2,or3.The soil pressure profiles are shown in Figure6.Figure 6
Soil pressure profiles forn=1,2,3.The first three buckling loads were determined for each case as shown in Table4. It is observed that the buckling loads experienced an increase with the increase in the exponent (n). This is reasonable because, as shown in Figure 6, the area under the soil pressure distribution increases as power n increases.Table 4
Critical buckling loads for nonuniform elastic foundation example.
n
P
c
r
1
P
c
r
2
P
c
r
3
1
25.0219
43.3211
90.5166
2
32.4221
45.9158
91.4651
3
40.7501
53.5283
93.133The mode shapes are shown in Figure7. The first mode shape shifts to the left as the exponent increases, which is expected as the relative soil stiffness is less in that area.Figure 7
Mode shapes for (a)n=1; (b) n=2; (c) n=3.
(a)
(b)
(c)


### 5.2. Example  2: Crack Propagation

In this example, a simply supported column was analyzed. As the maximum moment in simply supported structural elements develops at the midspan, it was assumed that crack initiation occurred at the center of the column. The crack was modeled by a change in the cross-section of column, making it nonprismatic, as shown in Figure8. This change was attained by reducing the moment of inertia (I). The propagation of the cracks in the column was modeled by increasing the cracked length (Lcr). The following constants were used for the columns in this example: E=1, I=1, L=1, and β=0. The matrix of the columns to be modeled in this example was generated by varying the crack length to column length ratio (Lcr/L) and moment of inertia of the cracked section to gross moment of inertia (Icr/I) ratio. The varying parameters from which the combinations were generated are shown in Table 5.Table 5
Cracking parameters for example 2.
L
c
r
/
L
0
0.1
0.2
0.3
0.4
0.5
I
c
r
/
I
1
0.95
0.6
0.4
—
—Figure 8
Column elevation for example 2.The analysis was performed, and the results for the columns are shown in Table6. For each moment of inertia ratio, the first column shows the buckling load obtained, while the second column shows the buckling load normalized against Euler’s buckling load. Euler’s buckling load is obtained for the uncracked prismatic column, which corresponds to a crack length ratio of 0 (row 1). The compression action in the column forces the crack to close; however, this causes an imperfection at that location which encourages the initiation of buckling from it.Table 6
Critical buckling loads and their normalized values for example 2.
I
c
r
/
I
0.95
0.6
0.4
L
c
r
/
L
P
c
r
P
c
r
/
P
e
P
c
r
P
c
r
/
P
e
P
c
r
P
c
r
/
P
e
0
9.869599
1
9.869599
1
9.869599
1
0.1
9.767448
0.98965
8.695862
0.881075
7.53605
0.763562
0.2
9.67209
0.979988
7.806693
0.790984
6.148034
0.622926
0.3
9.58769
0.971437
7.157659
0.725223
5.29744
0.536743
0.4
9.517094
0.964284
6.693711
0.678215
4.755724
0.481856
0.5
9.461763
0.958678
6.370312
0.645448
4.405123
0.446332In order to further analyze the obtained results, Figures9 and 10 were plotted. Figure 9 plots the normalized buckling loads against the cracked span ratios. A reduction in the buckling load as the cracked span ratio increased was observed. This is logical as the propagation of the cracks reduces the overall moment of inertia of the column and thus decreases the load required to buckle the column. Each series of data points was fitted with a second-degree polynomial function. These regression equations provided excellent correlation (R2=0.99).Figure 9
Normalized buckling load versus cracked span ratio for example 2.Figure 10
Normalized bucking load versus cracked moment of inertia ratio for example 2.Figure10 plots the normalized buckling loads against the cracked moment of inertia ratios. It was noted that the buckling load decreased as the cracked moment of inertia ratio decreased. This is rational because the smaller the ratio is, the larger the crack is and, thus, the smaller the overall moment of inertia becomes. An interesting trend observed was that the curve was bilinear for small cracked span ratios and approached the standard linear behavior as the span ratio increased.


### 5.3. Example  3: Column on Linearly Variable Elastic Foundation

In this example, a simply supported column was analyzed under different distributions of linear soil pressure. Details on the analyzed soil stiffness functions are shown in Table7 and Figure 11. These functions provide a constant area under the soil stiffness profile and only differ in the distribution of that area. For this example the length (L) was taken to be 10.Table 7
Soil stiffness functions for example 3.
Case #
β
x
Type
1
3
Uniform (rectangular)
2
2.25
+
0.15
x
Trapezoidal
3
1.5
+
0.3
x
Trapezoidal
4
0.75
+
0.45
x
Trapezoidal
5
0.6
x
TriangularFigure 11
Soil stiffness profile along column span for example 3.In order to facilitate the analysis, a few parameters will be defined. The slope of the soil stiffness profile(β′) is calculated by taking the derivative as follows:(50)β′=dβxdx.Another parameter is Kappa(κ), which relates the uniform soil stiffness (β1) to the stiffness of the column (EI) according to the following equation:(51)κ=β14EI4.For this example, columns withκ values that correspond to the following β/EI ratios were analyzed: 0.1, 1, 10, 100, and 200. The modulus of elasticity and moment of inertia were varied based on this ratio. The obtained buckling loads were normalized against Euler’s buckling load and are shown in Table 8.Table 8
Normalized buckling loads for example 3.
Case
1
2
3
4
5
β
′
0
0.15
0.3
0.45
0.6
κ
P
c
r
/
P
e
0.294283096
4.054726368
4.034826
3.965174
3.865672
3.746269
0.52331757
11.75376884
11.42211
10.78894
10.01508
9.140704
0.930604859
35.42211055
33.97487
31.13065
27.47236
22.9799
1.65487546
112.4170854
104.1457
92.38693
77.44221
57.70352
1.967989671
157.8341709
146.1256
128.598
106.3065
76.12563Additionally, the data points were plotted in Figures12 and 13. Figure 12 plots the normalized buckling load versus κ. It is observed that as the value of κ increased, the normalized buckling load increased. Also, at low κ values, which indicates that the EI is relatively high compared to β, the soil stiffness distribution did not have a significant effect on the buckling load. On the contrary, columns resting on relatively stiff soil experienced a significant variation in the buckling load based on the soil stiffness distribution (about 51%). Data in Table 8 show that as the slope (β′) increases, the buckling load decreases. Figure 13 plots the normalized buckling load versus β′. Lower κ results are represented by straight lines, which indicate again that the distribution of soil stiffness did not significantly affect the buckling load for these columns. As κ increased, the decrease in the critical load became more evident with the increase in β′.Figure 12
Normalized buckling load versusκ for example 3.Figure 13
Normalized buckling load versusβ′ for example 3.


## 5.1. Example  1: Column on Variable (Linear and Nonlinear) Elastic Foundation

After validating the proposed formulation in the previous section, the case of nonuniform soil pressure was investigated. A simply supported column was analyzed with 20 segments. The column’s parameters were as follows:EI=100 and L=10, and, similar to what Gere and Carter [13] did for the moment of inertia in nonprismatic columns, the soil stiffness (β) was taken to be as follows:(49)β=1+xLn,where n=1,2,or3.The soil pressure profiles are shown in Figure6.Figure 6
Soil pressure profiles forn=1,2,3.The first three buckling loads were determined for each case as shown in Table4. It is observed that the buckling loads experienced an increase with the increase in the exponent (n). This is reasonable because, as shown in Figure 6, the area under the soil pressure distribution increases as power n increases.Table 4
Critical buckling loads for nonuniform elastic foundation example.
n
P
c
r
1
P
c
r
2
P
c
r
3
1
25.0219
43.3211
90.5166
2
32.4221
45.9158
91.4651
3
40.7501
53.5283
93.133The mode shapes are shown in Figure7. The first mode shape shifts to the left as the exponent increases, which is expected as the relative soil stiffness is less in that area.Figure 7
Mode shapes for (a)n=1; (b) n=2; (c) n=3.
(a)
(b)
(c)


## 5.2. Example  2: Crack Propagation

In this example, a simply supported column was analyzed. As the maximum moment in simply supported structural elements develops at the midspan, it was assumed that crack initiation occurred at the center of the column. The crack was modeled by a change in the cross-section of column, making it nonprismatic, as shown in Figure8. This change was attained by reducing the moment of inertia (I). The propagation of the cracks in the column was modeled by increasing the cracked length (Lcr). The following constants were used for the columns in this example: E=1, I=1, L=1, and β=0. The matrix of the columns to be modeled in this example was generated by varying the crack length to column length ratio (Lcr/L) and moment of inertia of the cracked section to gross moment of inertia (Icr/I) ratio. The varying parameters from which the combinations were generated are shown in Table 5.Table 5
Cracking parameters for example 2.
L
c
r
/
L
0
0.1
0.2
0.3
0.4
0.5
I
c
r
/
I
1
0.95
0.6
0.4
—
—Figure 8
Column elevation for example 2.The analysis was performed, and the results for the columns are shown in Table6. For each moment of inertia ratio, the first column shows the buckling load obtained, while the second column shows the buckling load normalized against Euler’s buckling load. Euler’s buckling load is obtained for the uncracked prismatic column, which corresponds to a crack length ratio of 0 (row 1). The compression action in the column forces the crack to close; however, this causes an imperfection at that location which encourages the initiation of buckling from it.Table 6
Critical buckling loads and their normalized values for example 2.
I
c
r
/
I
0.95
0.6
0.4
L
c
r
/
L
P
c
r
P
c
r
/
P
e
P
c
r
P
c
r
/
P
e
P
c
r
P
c
r
/
P
e
0
9.869599
1
9.869599
1
9.869599
1
0.1
9.767448
0.98965
8.695862
0.881075
7.53605
0.763562
0.2
9.67209
0.979988
7.806693
0.790984
6.148034
0.622926
0.3
9.58769
0.971437
7.157659
0.725223
5.29744
0.536743
0.4
9.517094
0.964284
6.693711
0.678215
4.755724
0.481856
0.5
9.461763
0.958678
6.370312
0.645448
4.405123
0.446332In order to further analyze the obtained results, Figures9 and 10 were plotted. Figure 9 plots the normalized buckling loads against the cracked span ratios. A reduction in the buckling load as the cracked span ratio increased was observed. This is logical as the propagation of the cracks reduces the overall moment of inertia of the column and thus decreases the load required to buckle the column. Each series of data points was fitted with a second-degree polynomial function. These regression equations provided excellent correlation (R2=0.99).Figure 9
Normalized buckling load versus cracked span ratio for example 2.Figure 10
Normalized bucking load versus cracked moment of inertia ratio for example 2.Figure10 plots the normalized buckling loads against the cracked moment of inertia ratios. It was noted that the buckling load decreased as the cracked moment of inertia ratio decreased. This is rational because the smaller the ratio is, the larger the crack is and, thus, the smaller the overall moment of inertia becomes. An interesting trend observed was that the curve was bilinear for small cracked span ratios and approached the standard linear behavior as the span ratio increased.


## 5.3. Example  3: Column on Linearly Variable Elastic Foundation

In this example, a simply supported column was analyzed under different distributions of linear soil pressure. Details on the analyzed soil stiffness functions are shown in Table7 and Figure 11. These functions provide a constant area under the soil stiffness profile and only differ in the distribution of that area. For this example the length (L) was taken to be 10.Table 7
Soil stiffness functions for example 3.
Case #
β
x
Type
1
3
Uniform (rectangular)
2
2.25
+
0.15
x
Trapezoidal
3
1.5
+
0.3
x
Trapezoidal
4
0.75
+
0.45
x
Trapezoidal
5
0.6
x
TriangularFigure 11
Soil stiffness profile along column span for example 3.In order to facilitate the analysis, a few parameters will be defined. The slope of the soil stiffness profile(β′) is calculated by taking the derivative as follows:(50)β′=dβxdx.Another parameter is Kappa(κ), which relates the uniform soil stiffness (β1) to the stiffness of the column (EI) according to the following equation:(51)κ=β14EI4.For this example, columns withκ values that correspond to the following β/EI ratios were analyzed: 0.1, 1, 10, 100, and 200. The modulus of elasticity and moment of inertia were varied based on this ratio. The obtained buckling loads were normalized against Euler’s buckling load and are shown in Table 8.Table 8
Normalized buckling loads for example 3.
Case
1
2
3
4
5
β
′
0
0.15
0.3
0.45
0.6
κ
P
c
r
/
P
e
0.294283096
4.054726368
4.034826
3.965174
3.865672
3.746269
0.52331757
11.75376884
11.42211
10.78894
10.01508
9.140704
0.930604859
35.42211055
33.97487
31.13065
27.47236
22.9799
1.65487546
112.4170854
104.1457
92.38693
77.44221
57.70352
1.967989671
157.8341709
146.1256
128.598
106.3065
76.12563Additionally, the data points were plotted in Figures12 and 13. Figure 12 plots the normalized buckling load versus κ. It is observed that as the value of κ increased, the normalized buckling load increased. Also, at low κ values, which indicates that the EI is relatively high compared to β, the soil stiffness distribution did not have a significant effect on the buckling load. On the contrary, columns resting on relatively stiff soil experienced a significant variation in the buckling load based on the soil stiffness distribution (about 51%). Data in Table 8 show that as the slope (β′) increases, the buckling load decreases. Figure 13 plots the normalized buckling load versus β′. Lower κ results are represented by straight lines, which indicate again that the distribution of soil stiffness did not significantly affect the buckling load for these columns. As κ increased, the decrease in the critical load became more evident with the increase in β′.Figure 12
Normalized buckling load versusκ for example 3.Figure 13
Normalized buckling load versusβ′ for example 3.


## 6. Conclusion

In this paper, an approach to determine buckling loads and their associated mode shapes semianalytically for nonprismatic columns with variable boundary conditions resting on nonuniform elastic foundation was formulated. An exact solution for the segment buckling equation was derived under all possible cases that govern the differential equation. Coupling the solutions of the various segments resulted in a system of equations that can be used to numerically solve for the buckling load and plot the mode shape. Sensitivity analysis and benchmarking verification of the derived solution is performed next. Finally, three examples addressing the comparison of linear and nonlinear elastic foundation, crack propagation using nonprismatic column analysis, and investigating the effect of changing soil stiffness distributionβ on the buckling load are examined.



---

*Source: 5976098-2017-04-27.xml*