Research Article  Open Access
Buckling of Nonprismatic Column on Varying Elastic Foundation with Arbitrary Boundary Conditions
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 constants (modulus of elasticity), (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.
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 crosssections 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 crosssections 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 Figure 1 can be derived using the bifurcation method resulting in the following equation:
The boundary conditions (BCs) are as follows:
(1) Shear
(2) Momentwhere is the column elastic modulus function, is the column moment of inertia function, is the column axial load, is the elastic foundation spring stiffness function, is the column lateral deflection function, is the distance from the bottom of the column, is the bottom lateral discrete spring stiffness, is the top lateral discrete spring stiffness, is the bottom rotational discrete spring stiffness, is the top rotational discrete spring stiffness, is the length of the column, is the column shear force function, and is the column moment function.
Such a differential equation can prove to be extremely difficult to solve analytically for arbitrary , , 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 , , and , as shown in Figure 2. , , 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:
An exact solution was then obtained for each segment using the classical buckling equation of columns on an elastic foundation. The domain of for each segment was taken to be :where is number of segments, is segment length, and 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 (). The segment differential equation becomesLet , which can be substituted into (6) resulting inwhereThe deflection function is then obtained as follows:
Case 2 (). Let , which can be substituted into (5) resulting inThe deflection function is then found to be
Case 3 (). Let , which can be substituted into (5) resulting inwhereThe deflection function is then obtained to be as follows
Case 4 (). Let , which can be substituted into (5) resulting in
Let implying which results in the following two equations with two unknowns, :
From (17) we get
Substituting into (16) we obtain
We know that , but it is obvious that . Therefore, since . This leaves the following: .
Using (17) we can then obtain as follows:
Therefore, is found to bewhere
The deflection function is then obtained to bewhere 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
(2) Slope
(3) Moment
(4) Shear
Finally, in order to obtain the coefficient matrix (Pmatrix), the following boundary conditions were established.
(1) Shear
(2) Moment
Solving the following system of equations would then render the buckling load:where is the Matrix and = .
2.3. The Matrix
The structure of the matrix is as follows.
(1) One Segment
(2) Two Segments
(3) Three Segments
For the case of an assembly with more than three segments, the structure of the 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 matrix.
The submatrices inside the 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 Segments
Case 1 ().
Case 2 ().
Case 3 ().
Case 4 ().
(b) First Segment
Case 1 ().
Case 2 ().
Case 3 ().
Case 4 ().
(c) Last Segment
Case 1 ().
Case 2 ().
Case 3 ().
Case 4 ().
2.4. The Mode Shape
The mode shape can be obtained by substituting the buckling load into the matrix and setting one of the constant coefficients to be equal to 1. The original system looks like the following:
After letting , the following system of equations is obtained.
The constant coefficients can then be calculated by solving the system as follows:where is the reduced matrix, , and .
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 crosssection and soil stiffness varying to account for all approximations involved in the solution. The following functions were used: and . 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.

The 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. Figure 3 shows the first three mode shapes, which did not vary with the change in the number of segments.
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 fixedpinned column was evaluated with the following parameters: , , and [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.
The following cases are those for a nonprismatic column, as derived by Gere and Carter [13]. For Case 2, nonprismatic fixedpinned column, as derived by Gere and Carter [13], was analyzed. The parameters were as follows: , , and , 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.
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: . For the pinnedpinned reference case, the Euler buckling load was obtained to be 9.8696. The results obtained are shown in Table 2.

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: , , and . The results obtained are also shown in Table 3.

For Case 5, a column resting on a variable elastic foundation was analyzed. The column had the following parameters: N·m^{2}, m, and . 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: and , 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:where
The soil pressure profiles are shown in Figure 6.
The first three buckling loads were determined for each case as shown in Table 4. It is observed that the buckling loads experienced an increase with the increase in the exponent . This is reasonable because, as shown in Figure 6, the area under the soil pressure distribution increases as power increases.

The mode shapes are shown in Figure 7. 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.
(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 crosssection of column, making it nonprismatic, as shown in Figure 8. This change was attained by reducing the moment of inertia . The propagation of the cracks in the column was modeled by increasing the cracked length . The following constants were used for the columns in this example: , , , and . The matrix of the columns to be modeled in this example was generated by varying the crack length to column length ratio and moment of inertia of the cracked section to gross moment of inertia ratio. The varying parameters from which the combinations were generated are shown in Table 5.

The analysis was performed, and the results for the columns are shown in Table