# Buckling Analyses of Axially Functionally Graded Nonuniform Columns with Elastic Restraint Using a Localized Differential Quadrature Method

Authors: Yasin Yilmaz; Zekeriya Girgin; Savas Evran
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A localized differential quadrature method (LDQM) is introduced for buckling analysis of axially functionally graded nonuniform columns with elastic restraints. Weighting coefficients of differential quadrature discretization are obtained making use of neighboring points in forward and backward type schemes for the reference grids near the beginning and end boundaries of the physical domain, respectively, and central type scheme for the reference grids inside the physical domain. Boundary conditions are directly implemented into weighting coefficient matrices, and there is no need to use fictitious points near the boundaries. Compatibility equations are not required because the governing differential equation is discretized only once for each reference grid using neighboring points and variation of flexural rigidity is taken to be continuous in the axial direction. A large case of columns having different variations of cross-sectional profile and modulus of elasticity in the axial direction are considered. The results for nondimensional critical buckling loads are compared to the analytical and numerical results available in the literature. Some new results are also given. Comparison of the results shows the potential of the LDQM for solving such generalized eigenvalue problems governed by fourth-order variable coefficient differential equations with high accuracy and less computational effort.


---

## Body

## 1. Introduction

Beams and columns with variable flexural rigidities are commonly used in complex structures to achieve a better distribution of strength and weight and sometimes to satisfy architectural and functional requirements. The accurate prediction of dynamic behavior of beams and columns, particularly when the properties of the material and cross-section are variable, is of crucial importance in many areas of science and engineering such as civil, mechanical, biomedical, and aerospace engineering. Moreover, elastic restraints may also become an important issue for determination of critical buckling loads. The continuous change of the material properties can be achieved by gradually varying the volume fraction of the individual constituent materials from one point to another through any of the spatial coordinates. These types of materials are specifically called as functionally graded materials (FGMs). Continuity in FGM properties allows the elimination of interlayer delamination due to high stress and crack initiation and propagation in intermediate faces caused by high plastic deformation, frequently seen in laminated composites. The property variation of FGMs can be tailored to obtain the desired mechanical properties for different applications. With the advent of more advanced techniques, FGMs are now used in production of beam, plate, and shell structures that are emerging as promising structural elements in today’s industry (automotive and aircraft industry, space vehicles, machine elements, etc.) [1].There are many analytical and numerical studies for buckling analysis of beams and columns in the literature, but most of them are limited to the case of uniform case and a few studies take into account variation of both material properties and cross-section through the axial direction for columns with elastic restraints. Closed form solutions for these types of columns are only available for some special cases [2–8], and to obtain a general solution, numerical methods should be employed. As far as numerical methods are concerned, differential quadrature method (DQM) is claimed to be an efficient and powerful numerical technique with minimum computational effort and also an alternative to the finite difference and finite element methods. DQM, introduced firstly by Bellmann and his coworkers [9, 10], discretizes any derivative at a point in the solution domain by a weighted linear sum of function values along the direction of its respective coordinate. Since its advance, DQM has been successfully applied to a variety of problems in engineering science and a review can be found in [11]. The DQM and its modified versions are applied to buckling analyses of nonuniform and functionally graded columns in [12–15]. Other numerical methods such as differential transform method [15], semi-inverse method [16–18], functional perturbation method [19], variational iteration method [20, 21], homotopy perturbation method [22, 23], finite difference method [24], finite element method [25], integral-equation approach [26], dynamic stiffness method [27], Newton’s eigenvalue iteration method [28], and Fredholm integral technique [29] are also successfully applied to buckling analyses of nonuniform and axially functionally graded columns. Recently, Huang and Luo [30] proposed a method that uses power series and integral technique to obtain a polynomial characteristic equation in terms of critical buckling loads for axially inhomogeneous beams with elastic restraint. Then they solved this equation for smallest positive root using a commercial software.The key procedure in DQM lies in the determination of the weighting coefficients. Also, implementation of boundary conditions is another important issue. With this respect, many studies are conducted to overcome deficiencies of the DQM. Detailed information about the DQM and its development stages can be found in [31]. Even though, DQM has been regarded as an efficient numerical technique with high accuracy at the beginning, large number of grid numbers cannot be used, that is, more than 21, due to the stability problems. To overcome these deficiencies, more advanced version of DQM, localized differential quadrature method (LDQM) is proposed and applied to some engineering and physical problems such as 2D wave equation [32], 2D stream function formulation of Navier–Stokes equations [33], and mild slope equation [34]. The basic idea in the LDQ method is to apply DQ approximation to a small neighborhood of the grid point of interest rather than to the whole domain. The derivatives at each grid point are then approximated by a weighted sum of function values on its neighboring points, rather than on all of the grid points. Thus, a very accurate solution without losing stability can be obtained. The discussion and basic ideas beyond the LDQM are given in detail in [35].In this study, LQDM is applied to the buckling analyses of nonuniform axially functionally graded columns with elastic restraints for different boundary conditions. To the best of the authors’ knowledge, in the open literature, LDQM is not applied before to solutions of generalized eigenvalue problems governed by fourth-order variable coefficient differential equations. Fictitious points are not used for implementation of boundary conditions at the tip nodes and boundary conditions are directly implemented into weighting coefficient matrices. It is seen that problems encountered in original DQ method for implementation of clamped boundary conditions into weighting coefficient matrices are eliminated using LDQM. There is no restriction in the continuous variation of cross-section and elastic modulus in the axial direction. In the following at first, mathematical formulation of the problem is provided and an introduction is given for classical DQM and LDQM. Then numerical results for critical buckling loads are given for a large case of columns using LDQM.


## 2. Mathematical Formulation

### 2.1. Governing Differential Equation

A nonuniform axially functionally graded column of lengthL subjected to an axial compressive load P at the centroid axis with a continuous elastic restraint is considered (see Figure 1). The material properties and cross-section of the beam are assumed to vary continuously along the axial direction, x. According to Euler-Bernoulli beam theory, the effects of shear deformation and rotary inertia are neglected; thus the governing differential equation for buckling of an axially FG nonuniform beam on elastic foundation can be written as follows [8]:
(1)d2dx2[E(x)I(x)d2wdx2]+Pd2wdx2+kw=0,0<x<L,
where w(x) is the lateral displacement, E(x) and I(x) are axially varying modulus of elasticity and moment of inertia of the column, respectively, and k is the stiffness of the uniformly distributed lateral springs per unit length. Introducing nondimensional coordinate and deflection as
(2)ξ=xL,W(ξ)=w(x)L,
(1) can be written as [30]
(3)d2dξ2[E(ξ)I(ξ)d2Wdξ2]+α*d2Wdξ2+β*W=0,0<ξ<1,
where α*=PL2 is the normalized critical load and β*=kL4 is the normalized restrained stiffness parameter.Figure 1
Column with continuous elastic restraint.


### 2.2. Boundary Conditions

End supports of the beam directly affect the critical buckling load of the column. In the present study, four types of boundary conditions shown in Figure2 are considered [8]. The explicit expressions for the boundary conditions can be given in terms of W(ξ) as follows.Clamped support (C):(4a)W=0,dWdξ=0.Pin support (P):(4b)W=0,d2Wdξ2=0.Free end (F):(4c)d2Wdξ2=0,ddξ[E(ξ)I(ξ)d2Wdξ2]+α*dWdξ=0.Guided end (G):(4d)dWdξ=0,ddξ[E(ξ)I(ξ)d2Wdξ2]+α*dWdξ=0.Figure 2
Boundary conditions.The problem is to find the critical buckling load of the column using the fourth-order variable coefficient ordinary differential equation (3) subject to corresponding boundary conditions given in (4a), (4b), (4c), and (4d). A closed-form solution of the problem is not generally possible for arbitrarily varying coefficients except for some special cases [3–7] and it is of crucial importance to find the smallest root of the equations for general case in engineering applications. In Section 3, LDQM discretization of the governing differential equation and boundary conditions will be shown for general case of axial variation of material properties and cross-section and, in Section 4, various types of nonuniform axially FG beams will be considered and corresponding numerical solutions will be obtained for the critical buckling loads.


## 2.1. Governing Differential Equation

A nonuniform axially functionally graded column of lengthL subjected to an axial compressive load P at the centroid axis with a continuous elastic restraint is considered (see Figure 1). The material properties and cross-section of the beam are assumed to vary continuously along the axial direction, x. According to Euler-Bernoulli beam theory, the effects of shear deformation and rotary inertia are neglected; thus the governing differential equation for buckling of an axially FG nonuniform beam on elastic foundation can be written as follows [8]:
(1)d2dx2[E(x)I(x)d2wdx2]+Pd2wdx2+kw=0,0<x<L,
where w(x) is the lateral displacement, E(x) and I(x) are axially varying modulus of elasticity and moment of inertia of the column, respectively, and k is the stiffness of the uniformly distributed lateral springs per unit length. Introducing nondimensional coordinate and deflection as
(2)ξ=xL,W(ξ)=w(x)L,
(1) can be written as [30]
(3)d2dξ2[E(ξ)I(ξ)d2Wdξ2]+α*d2Wdξ2+β*W=0,0<ξ<1,
where α*=PL2 is the normalized critical load and β*=kL4 is the normalized restrained stiffness parameter.Figure 1
Column with continuous elastic restraint.


## 2.2. Boundary Conditions

End supports of the beam directly affect the critical buckling load of the column. In the present study, four types of boundary conditions shown in Figure2 are considered [8]. The explicit expressions for the boundary conditions can be given in terms of W(ξ) as follows.Clamped support (C):(4a)W=0,dWdξ=0.Pin support (P):(4b)W=0,d2Wdξ2=0.Free end (F):(4c)d2Wdξ2=0,ddξ[E(ξ)I(ξ)d2Wdξ2]+α*dWdξ=0.Guided end (G):(4d)dWdξ=0,ddξ[E(ξ)I(ξ)d2Wdξ2]+α*dWdξ=0.Figure 2
Boundary conditions.The problem is to find the critical buckling load of the column using the fourth-order variable coefficient ordinary differential equation (3) subject to corresponding boundary conditions given in (4a), (4b), (4c), and (4d). A closed-form solution of the problem is not generally possible for arbitrarily varying coefficients except for some special cases [3–7] and it is of crucial importance to find the smallest root of the equations for general case in engineering applications. In Section 3, LDQM discretization of the governing differential equation and boundary conditions will be shown for general case of axial variation of material properties and cross-section and, in Section 4, various types of nonuniform axially FG beams will be considered and corresponding numerical solutions will be obtained for the critical buckling loads.


## 3. Numerical Discretization and Method of Solution

### 3.1. Classical Differential Quadrature Method (DQM)

The method of differential quadrature (DQ) is based on an assumption that the derivative of a function with respect to a space variable at a given discrete point can be expressed as a weighted linear sum of the function values at all discrete points in the solution domain. Weighting coefficients of DQM can be calculated in several ways. In order to find simple algebraic expressions for the weighting coefficients without restricting the choice of grid points, the generalized DQM was developed by Shu and Richards [36]. Consider a function f(x) prescribed in a field domain a≤x≤b. Let f(xi) be the function values specified in a finite set of N discrete points xi(i=1,2,…,N) of the field domain. The rth-order derivative of f(x) at any discrete point xi can be written in DQ analog form as
(5)drf(xi)dxr=∑j=1Naij(r)f(xj),i=1,2,…,N,
where aij(r) are the weighting coefficients of the rth-order derivative of the function f(x) associated with points xi.In the generalized DQM, the test functions are assumed to be the Lagrange interpolation test functions such as(6)ℓj(x)=M(x)(x-xj)M(1)(xj),j=1,2,…,N,
where
(7)M(x)=∏m=1N(x-xm),M(1)(xj)=dM(xj)dx=∏m=1,m≠jN(xj-xm).The weighting coefficients for the first-order derivative can be obtained using Lagrange interpolation polynomial as follows:(8a)aij(1)=dℓj(xi)dx=M(1)(xi)(xi-xj)M(1)(xj),i,j=1,2,…,N,i≠j,(8b)aii(1)=dℓi(xi)dx=-∑j=1,j≠iNaij(1),i=1,2,…,N,and these coefficients give the [A(1)] matrix as
(9)[A(1)]=[a11(1)a12(1)⋯a1N(1)a21(1)a22(1)⋯a2N(1)⋮⋮⋱⋮aN1(1)aN2(1)⋯aNN(1)].
Similarly, the weighting coefficients for the rth-order derivative can be evaluated using the following equations:(10a)aij(r)=drℓj(xi)dxr=r(aii(r-1)aij(1)-aij(r-1)(xi-xj)),i,j=1,2,…,N,i≠j,r≥2,(10b)aii(r)=drℓi(xi)dxr=-∑j=1,j≠iNaij(r),i=1,2,…,N,r≥2.The weighting coefficients for second- and higher-order derivatives can also be evaluated using the following expression:
(11)aij(r)=∑k=1Naik(r-1)akj(1),i,j=1,2,…,N,r≥2.
These coefficients give the matrix A(r)  which can be interpreted as follows:
(12)[A(r)]=[a11(r)a12(r)⋯a1N(r)a21(r)a22(r)⋯a2N(r)⋮⋮⋱⋮aN1(r)aN2(r)⋯aNN(r)].Whenx is taken in the interval x∈[0,L] instead of x∈[0,1], the right-hand side of (14) must be divided by Lr which can be smaller or greater than 1. Then, the weighting coefficient matrix takes the following form:
(13)[A(r)]=(ddx)r=drdxr=1Lr[a11(r)a12(r)⋯a1N(r)a21(r)a22(r)⋯a2N(r)⋮⋮⋱⋮aN1(r)aN2(r)⋯aNN(r)].After calculating the matrixA(1), weighting coefficient matrices for second- and higher-order derivatives can also be calculated using the following formulae:
(14)[A(r)]=[A(1)][A(r-1)]=[A(r-1)][A(1)].


### 3.2. Localized Differential Quadrature Method (LDQM) and Sampling Points

In LDQM, firstly neighboring grids of any discrete pointxi in the computational domain should be determined according to the position of the point of interest and order of the first derivative approximation. For example if the function is discretized at the beginning boundary of the physical domain with respect to space variable x (at x=0) or nearby it and a sixth-order first derivative approximation is used that means 7 neighboring points are needed (Nx=7), the selection of the neighboring grid points should be forward type in the direction of the space variable, and if the function is discretized at the end boundary of the physical domain with respect to space variable x (at x=L) or nearby it, the selection of the neighboring grid points should be backward type. At the other interior reference points central type scheme is used. Figures 3 and 4 show the selection of neighboring points of a discrete point xi in the solution domain for the previously mentioned cases. Then, the discretization of the first-order derivative of a function f(x) with respect to space variable, x, at any discrete point xi can be approximated using a weighted linear combination of the function values at some of the neighboring reference points within the computational domain as [32]
(15)df(xi)dx=∑j∈Siaij(1)f(xj),i=1,2,…,N,
where Si represents the corresponding set of the neighboring nodes for the discrete grid point xi in the domain or at the boundaries, N is the total amount of grid points in the direction of x. Weighting coefficients of the first-order derivative with respect to spatial coordinate x can be evaluated using the following equations:Figure 3
Discretization of boundary and near boundary reference points forNx=7 and N=21.Figure 4
Discretization of an interior reference point forNx=7 and N=21.(16a)
a
i
j
(
1
)
=
∏
m
∈
S
i
,
m
≠
i
(
x
i
-
x
m
)
(
x
i
-
x
j
)
·
∏
m
∈
S
i
,
m
≠
j
(
x
j
-
x
m
)
,
i
=
1,2
,
…
,
N
,
j
∈
S
i
,
i
≠
j
,
(16b)
a
i
i
(
1
)
=
-
∑
j
∈
S
i
,
j
≠
i
a
i
j
(
1
)
,
i
=
1,2
,
…
,
N
.
Similarly, the discretization of the higher-order derivatives of f(x) with respect to space variable, x, at any discrete point xi can be expressed as
(17)drf(xi)dxr=∑j∈Siaij(r)f(xj),i=1,2,…,N,r≥2.
The weighting coefficients for the rth-order derivative at any discrete point xi can be evaluated using the following equations:(18a)aij(r)=r·(aijaii(r-1)-aij(r-1)(xi-xj)),i=1,2,…,N,j∈Si,i≠j,r≥2,(18b)aii(r)=-∑j∈Si,j≠iaij(r),i=1,2,…,N,r≥2.It is instructive to note that (18a) and (18b) are rewritten forms of (10a) and (10b) in the neighborhood of the reference point xi.A frequently used and convenient choice for the sampling points is that of the equally spaced sampling points which can be given in normalized coordinate of the space variablexi by
(19)ξi=i-1N-1,i=1,2,…,N.
In the present study, equally spaced sampling points are used.


### 3.3. LDQM Formulation of the Governing Differential Equation and Boundary Conditions

IntroducingD(ξ)=E(ξ)I(ξ) for flexural rigidity, (3) can be expanded into following form:
(20)D(ξ)d4Wdξ4+2dD(ξ)dξd3Wdξ3+d2D(ξ)dξ2d2Wdξ2+α*d2Wdξ2+β*W=0.
For the numerical computation equation (20) can be discretized using LDQM as follows:
(21)D(ξi)∑j∈Siaij(4)W(ξj)+2dD(ξ)dξ|ξ=ξi∑j∈Siaij(3)W(ξj)+d2D(ξ)dξ2|ξ=ξi∑j∈Siaij(2)W(ξj)+α*∑j∈Siaij(2)W(ξj)+β*W(ξi)=0,i=1,2,…,N.
Boundary conditions given in (4a), (4b), (4c), and (4d) can be discretized at ξ=0 using LDQM as follows.Clamped support (C):(22a)W(ξ1)=0,∑j∈Sia1j(1)W(xj)=0.Pin support (P):(22b)W(ξ1)=0,∑j∈Sia1j(2)W(xj)=0.Free end (F):(22c)∑j∈Sia1j(2)W(xj)=0,dD(ξ)dξ|ξ=0∑j∈Sia1j(2)W(xj)+D(ξ1)∑j∈Sia1j(3)W(xj)+α*∑j∈Sia1j(1)W(xj)=0.Guided end (G):(22d)∑j∈Sia1j(1)W(xj)=0,dD(ξ)dξ|ξ=0∑j∈Sia1j(2)W(xj)+D(ξ1)∑j∈Sia1j(3)W(xj)+α*∑j∈Sia1j(1)W(xj)=0.Boundary conditions at ξ=1 can also be discretized in a similar manner. Boundary conditions given in (22a), (22b), (22c), and (22d) can directly be substituted into weighting coefficient matrices as described by Shu [31]. It will be shown in the next section that the disadvantages of substituting boundary conditions into weighting coefficient matrices, encountered in classical DQM for clamped boundary condition, disappeared in LDQM. It is also worth mentioning that no compatibility equations are needed in LDQM because the governing differential equation is discretized only once for each reference grid using neighboring points and variation of flexural rigidity is continuous in the axial direction.Resulting set of algebraic equations can be put into matrix form to solve the generalized eigenvalue problem for critical buckling loads. Required matrix dimension is(N-2)×(N-2) for combination of clamped and simply supported end conditions since the function values at ξ=0 and ξ=1 are known. For combination of clamped and other boundary conditions (free and guided), the resulting matrix dimension is (N-1)×(N-1) since function value is only known at the clamped end.


## 3.1. Classical Differential Quadrature Method (DQM)

The method of differential quadrature (DQ) is based on an assumption that the derivative of a function with respect to a space variable at a given discrete point can be expressed as a weighted linear sum of the function values at all discrete points in the solution domain. Weighting coefficients of DQM can be calculated in several ways. In order to find simple algebraic expressions for the weighting coefficients without restricting the choice of grid points, the generalized DQM was developed by Shu and Richards [36]. Consider a function f(x) prescribed in a field domain a≤x≤b. Let f(xi) be the function values specified in a finite set of N discrete points xi(i=1,2,…,N) of the field domain. The rth-order derivative of f(x) at any discrete point xi can be written in DQ analog form as
(5)drf(xi)dxr=∑j=1Naij(r)f(xj),i=1,2,…,N,
where aij(r) are the weighting coefficients of the rth-order derivative of the function f(x) associated with points xi.In the generalized DQM, the test functions are assumed to be the Lagrange interpolation test functions such as(6)ℓj(x)=M(x)(x-xj)M(1)(xj),j=1,2,…,N,
where
(7)M(x)=∏m=1N(x-xm),M(1)(xj)=dM(xj)dx=∏m=1,m≠jN(xj-xm).The weighting coefficients for the first-order derivative can be obtained using Lagrange interpolation polynomial as follows:(8a)aij(1)=dℓj(xi)dx=M(1)(xi)(xi-xj)M(1)(xj),i,j=1,2,…,N,i≠j,(8b)aii(1)=dℓi(xi)dx=-∑j=1,j≠iNaij(1),i=1,2,…,N,and these coefficients give the [A(1)] matrix as
(9)[A(1)]=[a11(1)a12(1)⋯a1N(1)a21(1)a22(1)⋯a2N(1)⋮⋮⋱⋮aN1(1)aN2(1)⋯aNN(1)].
Similarly, the weighting coefficients for the rth-order derivative can be evaluated using the following equations:(10a)aij(r)=drℓj(xi)dxr=r(aii(r-1)aij(1)-aij(r-1)(xi-xj)),i,j=1,2,…,N,i≠j,r≥2,(10b)aii(r)=drℓi(xi)dxr=-∑j=1,j≠iNaij(r),i=1,2,…,N,r≥2.The weighting coefficients for second- and higher-order derivatives can also be evaluated using the following expression:
(11)aij(r)=∑k=1Naik(r-1)akj(1),i,j=1,2,…,N,r≥2.
These coefficients give the matrix A(r)  which can be interpreted as follows:
(12)[A(r)]=[a11(r)a12(r)⋯a1N(r)a21(r)a22(r)⋯a2N(r)⋮⋮⋱⋮aN1(r)aN2(r)⋯aNN(r)].Whenx is taken in the interval x∈[0,L] instead of x∈[0,1], the right-hand side of (14) must be divided by Lr which can be smaller or greater than 1. Then, the weighting coefficient matrix takes the following form:
(13)[A(r)]=(ddx)r=drdxr=1Lr[a11(r)a12(r)⋯a1N(r)a21(r)a22(r)⋯a2N(r)⋮⋮⋱⋮aN1(r)aN2(r)⋯aNN(r)].After calculating the matrixA(1), weighting coefficient matrices for second- and higher-order derivatives can also be calculated using the following formulae:
(14)[A(r)]=[A(1)][A(r-1)]=[A(r-1)][A(1)].


## 3.2. Localized Differential Quadrature Method (LDQM) and Sampling Points

In LDQM, firstly neighboring grids of any discrete pointxi in the computational domain should be determined according to the position of the point of interest and order of the first derivative approximation. For example if the function is discretized at the beginning boundary of the physical domain with respect to space variable x (at x=0) or nearby it and a sixth-order first derivative approximation is used that means 7 neighboring points are needed (Nx=7), the selection of the neighboring grid points should be forward type in the direction of the space variable, and if the function is discretized at the end boundary of the physical domain with respect to space variable x (at x=L) or nearby it, the selection of the neighboring grid points should be backward type. At the other interior reference points central type scheme is used. Figures 3 and 4 show the selection of neighboring points of a discrete point xi in the solution domain for the previously mentioned cases. Then, the discretization of the first-order derivative of a function f(x) with respect to space variable, x, at any discrete point xi can be approximated using a weighted linear combination of the function values at some of the neighboring reference points within the computational domain as [32]
(15)df(xi)dx=∑j∈Siaij(1)f(xj),i=1,2,…,N,
where Si represents the corresponding set of the neighboring nodes for the discrete grid point xi in the domain or at the boundaries, N is the total amount of grid points in the direction of x. Weighting coefficients of the first-order derivative with respect to spatial coordinate x can be evaluated using the following equations:Figure 3
Discretization of boundary and near boundary reference points forNx=7 and N=21.Figure 4
Discretization of an interior reference point forNx=7 and N=21.(16a)
a
i
j
(
1
)
=
∏
m
∈
S
i
,
m
≠
i
(
x
i
-
x
m
)
(
x
i
-
x
j
)
·
∏
m
∈
S
i
,
m
≠
j
(
x
j
-
x
m
)
,
i
=
1,2
,
…
,
N
,
j
∈
S
i
,
i
≠
j
,
(16b)
a
i
i
(
1
)
=
-
∑
j
∈
S
i
,
j
≠
i
a
i
j
(
1
)
,
i
=
1,2
,
…
,
N
.
Similarly, the discretization of the higher-order derivatives of f(x) with respect to space variable, x, at any discrete point xi can be expressed as
(17)drf(xi)dxr=∑j∈Siaij(r)f(xj),i=1,2,…,N,r≥2.
The weighting coefficients for the rth-order derivative at any discrete point xi can be evaluated using the following equations:(18a)aij(r)=r·(aijaii(r-1)-aij(r-1)(xi-xj)),i=1,2,…,N,j∈Si,i≠j,r≥2,(18b)aii(r)=-∑j∈Si,j≠iaij(r),i=1,2,…,N,r≥2.It is instructive to note that (18a) and (18b) are rewritten forms of (10a) and (10b) in the neighborhood of the reference point xi.A frequently used and convenient choice for the sampling points is that of the equally spaced sampling points which can be given in normalized coordinate of the space variablexi by
(19)ξi=i-1N-1,i=1,2,…,N.
In the present study, equally spaced sampling points are used.


## 3.3. LDQM Formulation of the Governing Differential Equation and Boundary Conditions

IntroducingD(ξ)=E(ξ)I(ξ) for flexural rigidity, (3) can be expanded into following form:
(20)D(ξ)d4Wdξ4+2dD(ξ)dξd3Wdξ3+d2D(ξ)dξ2d2Wdξ2+α*d2Wdξ2+β*W=0.
For the numerical computation equation (20) can be discretized using LDQM as follows:
(21)D(ξi)∑j∈Siaij(4)W(ξj)+2dD(ξ)dξ|ξ=ξi∑j∈Siaij(3)W(ξj)+d2D(ξ)dξ2|ξ=ξi∑j∈Siaij(2)W(ξj)+α*∑j∈Siaij(2)W(ξj)+β*W(ξi)=0,i=1,2,…,N.
Boundary conditions given in (4a), (4b), (4c), and (4d) can be discretized at ξ=0 using LDQM as follows.Clamped support (C):(22a)W(ξ1)=0,∑j∈Sia1j(1)W(xj)=0.Pin support (P):(22b)W(ξ1)=0,∑j∈Sia1j(2)W(xj)=0.Free end (F):(22c)∑j∈Sia1j(2)W(xj)=0,dD(ξ)dξ|ξ=0∑j∈Sia1j(2)W(xj)+D(ξ1)∑j∈Sia1j(3)W(xj)+α*∑j∈Sia1j(1)W(xj)=0.Guided end (G):(22d)∑j∈Sia1j(1)W(xj)=0,dD(ξ)dξ|ξ=0∑j∈Sia1j(2)W(xj)+D(ξ1)∑j∈Sia1j(3)W(xj)+α*∑j∈Sia1j(1)W(xj)=0.Boundary conditions at ξ=1 can also be discretized in a similar manner. Boundary conditions given in (22a), (22b), (22c), and (22d) can directly be substituted into weighting coefficient matrices as described by Shu [31]. It will be shown in the next section that the disadvantages of substituting boundary conditions into weighting coefficient matrices, encountered in classical DQM for clamped boundary condition, disappeared in LDQM. It is also worth mentioning that no compatibility equations are needed in LDQM because the governing differential equation is discretized only once for each reference grid using neighboring points and variation of flexural rigidity is continuous in the axial direction.Resulting set of algebraic equations can be put into matrix form to solve the generalized eigenvalue problem for critical buckling loads. Required matrix dimension is(N-2)×(N-2) for combination of clamped and simply supported end conditions since the function values at ξ=0 and ξ=1 are known. For combination of clamped and other boundary conditions (free and guided), the resulting matrix dimension is (N-1)×(N-1) since function value is only known at the clamped end.


## 4. Numerical Results and Discussions

In this section, LDQM is used to investigate the buckling behavior of axially FG nonuniform beams with a continuous elastic restraint. Firstly buckling behavior of homogeneous uniform columns with elastic restraints that has analytical solution is investigated to show the effectiveness of the method. Then the effects of variable cross-section and variable elastic modulus are investigated individually and both together.


### 4.1. Uniform Homogeneous Column with Elastic Restraints

For this case, (3) has an analytic solution and Wang et al. [8] have derived the stability criteria in nondimensional critical buckling load α=PL2/EI and normalized restrained stiffness parameter β=kL4/EI for the columns with different end conditions. Stability criteria are as follows:P-P Column:(23)sinT=0.C-P Column:(24)TcosTsinS-SsinTcosS=0.C-G Column:(25)TsinTcosS-ScosTsinS=0.C-C Column:(26)2ST[cosTcosS-1]+(T2+S2)sinTsinS=0.C-F Column:(27)[α(S2+T2)-2S2T2]cosTcosS-α(S2+T2)+(S4+T4)+ST[2α-(S2+T2)]sinTsinS=0,where
(28)S=α2-(α2)2-β,T=α2+(α2)2-β.
The stability criteria given in (23)–(27) are extremely nonlinear and finding the smallest root of the equation which is the critical buckling load is not easy for any assumed stiffness parameter. To show the effectiveness of the LDQM, nondimensional critical buckling loads, α=PL2/EI, are calculated for clamped-clamped beams with different normalized restrained stiffness parameters β=kL4/EI and compared to the exact results. In calculations the number of neighboring grids Nx is taken to be 11 and total number of grid points N is changed from 21 to 41. As can be seen from Table 1, critical buckling load converges to exact results as N increases from 11 to 41. In the following, neighboring grids Nx is taken to be 11 and total number of grid points N is taken to be 41 unless otherwise specified.Table 1
Non dimensional critical buckling loads (α=PL2/EI) for clamped-clamped uniform homogeneous columns with continuous elastic restraint (β=kL4/EI).
β
Exact [8]
DQMN=21
DQMN=31
DQMN=41
LDQMNx=11, N=21
LDQMNx=11, N=31
LDQMNx=11, N=41
0
39.47841760
39.47841735
39.46144027
—
39.47841219
39.47841756
39.47841760
50
43.26056590
43.26107726
38.55373074
—
43.26023292
43.26056213
43.26056583
100
47.00660086
47.00980042
37.89095672
—
47.00612534
47.00659571
47.00660072The critical buckling loads of the aforementioned columns have also been calculated by Atay and Coşkun [21] using a variation iteration approach and Huang and Luo [30] using power series and the integral technic. Nondimensional critical buckling loads for other boundary conditions are also obtained using LDQM and compared with the exact results [8], and with [21, 30]. The LDQM results given in Table 2 are in good agreement with the exact and other numerical results. It is also seen from Tables 1 and 2 that although the critical buckling load increases as the normalized restrained stiffness parameter increases for all of the boundary conditions, the column with C-G end condition is more sensitive to normalized restrained stiffness parameter change than the others.Table 2
Non dimensional critical buckling loads (α=PL2/EI) for uniform homogeneous columns with continuous elastic restraint (β=kL4/EI).
BC
β
=
0
β
=
50
β
=
100
Exact [8]
Reference [21]
Reference [30]
LDQM
Exact [8]
Reference [21]
Reference [30]
LDQM
Exact [8]
Reference [21]
Reference [30]
LDQM
C-P
20.1907
20.1908
20.190729
20.190729
24.2852
24.2855
24.285198
24.285198
28.3066
28.3080
28.306631
28.306631
P-P
9.8696
9.8696
9.869604
9.869604
14.9357
14.9357
14.935664
14.935664
20.0017
20.0017
20.001723
20.001723
C-G
9.8696
9.8696
9.869604
9.869604
23.5717
23.5722
23.571661
23.571659
32.6690
32.6482
32.668976
32.668976
C-F
2.4674
2.4674
2.467401
2.467401
8.8614
8.8614
8.861396
8.861396
11.9964
11.9964
11.996413
11.996413


### 4.2. Effect of Variable Cross-Section

Euler-Bernoulli columns with nonuniform cross-section along the axial direction are considered to see the effect of variable cross-section on critical buckling load. For this case elastic modulusE is taken to be constant and moment of inertia term changes according to a power law:
(29)I(ξ)=I0(1+aξ)p,a>-1,
where a is a geometric parameter and p=1, 2, and 3.First let us consider the casep=1 which corresponds to the columns of linearly varying width and constant height. Nondimensional critical buckling loads, α=PL2/EI0, are calculated using LDQM for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. The results are compared to available results in the literature and given in Tables 3 and 4. It is seen from Tables 3 and 4 that with a changes from negative to positive values, nondimensional critical buckling load increases since negative values of a correspond to decreasing width and positive values of a correspond to increasing width as ξ changes from 0 to 1. As normalized restrained stiffness parameter increases, critical buckling loads also increase for all of the boundary conditions as in the case of uniform columns.Table 3
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of linearly varying width and constant height (I(ξ)=I0(1+aξ)) for β=0.
BC
a
Reference [28]
Reference [3]
Reference [30]
LDQM
C-F
−0.5
—
—
—
2.06209226
0.5
—
—
—
2.81006133
1.0
3.12277257
3.11769623
3.11769623
3.11769628
F-C
−0.5
—
—
—
1.55884812
0.5
—
—
—
3.31117993
1.0
—
4.12418445
4.12418445
4.12418442
P-P
−0.5
—
—
—
7.25562477
0.5
—
—
—
12.25041131
1.0
14.50520092
14.51124954
14.51124954
14.51124954
C-P
−0.5
—
—
—
14.73942213
0.5
—
—
—
24.99361144
1.0
29.49596465
29.44896281
29.44896281
29.44896281
P-C
−0.5
—
—
—
14.72448141
0.5
—
—
—
24.99874044
1.0
—
29.47884426
29.47884426
29.47884426
C-C
−0.5
—
—
—
28.69697807
0.5
—
—
—
48.81145815
1.0
57.44524609
57.39395614
57.39395614
57.39395614
C-G
−0.5
—
—
—
7.16883539
0.5
—
—
—
12.19969805
1.0
—
—
—
14.33767077Table 4
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of linearly varying width and constant height (I(ξ)=I0(1+aξ)) with continuous elastic restraint (β=kL4/EI0).
BC
a
=
-
0.5
a
=
0.5
a
=
1
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
C-F
5.320851
7.249834
9.496180
7.002455
10.034833
13.926981
7.575957
10.969773
15.525985
F-C
5.484886
7.762992
10.365666
7.065179
9.787644
13.328385
7.822199
10.641702
14.499668
P-P
9.774589
12.290572
17.308835
14.778728
17.306502
22.359993
17.030858
19.549179
24.581143
C-P
16.663069
18.556531
22.240823
27.120454
29.236109
33.430687
31.622334
33.787588
38.091765
P-C
16.893794
19.045882
23.289167
26.982108
28.949869
32.834737
31.409689
33.326139
37.113062
C-C
30.583452
32.457228
36.164501
50.705246
52.591875
56.342983
59.283541
61.166905
64.914456
C-G
13.693532
18.852382
24.928160
19.887489
26.982550
38.296345
22.292684
29.808724
42.676370p
=
2 case corresponds to the columns of parabolic varying width and constant height. For this case nondimensional critical buckling loads, α=PL2/EI0, are calculated using LDQM for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. The results are compared to available results in the literature and given in Tables 5 and 6. Comparison of results given in Tables 3 and 4 with the results given in Tables 5 and 6 shows that critical buckling loads for p=2 case are higher than p=1 case for a>0 and critical buckling loads for p=2 case are lower than p=1 case for a<0 for all of the boundary conditions since decrease and increase in column width are influenced by a factor of 2 for the case of p=2. It is also seen from Tables 5 and 6 that although C-P and P-C columns have the same critical buckling loads for β=0, they have somewhat different critical buckling loads for β=25, 50, and 100.Table 5
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of parabolic varying width and constant height (I(ξ)=I0(1+aξ)2) for β=0.
BC
a
Reference [28]
Reference [3]
Reference [30]
LDQM
C-F
−0.5
—
—
—
1.68296637
0.5
—
—
—
3.17193009
1.0
3.85763006
3.83637692
3.83637692
3.83637679
F-C
−0.5
—
—
—
0.95909423
0.5
—
—
—
4.40655924
1.0
—
6.73186541
6.73186541
6.73186555
P-P
−0.5
—
—
—
5.19807211
0.5
—
—
—
15.07084413
1.0
20.79163330
20.79228846
20.79228846
20.79228846
C-P
−0.5
—
—
—
10.52729403
0.5
—
—
—
30.72425899
1.0
42.31345100
42.10917612
42.10917612
42.10917612
P-C
−0.5
—
—
—
10.52729403
0.5
—
—
—
30.72425899
1.0
—
42.10917612
42.10917612
42.10917612
C-C
−0.5
—
—
—
20.48084102
0.5
—
—
—
59.97123859
1.0
82.22844561
81.92336388
81.92336364
81.92336379
C-G
−0.5
—
—
—
5.19807212
0.5
—
—
—
15.07084413
1.0
—
—
—
20.79228847Table 6
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of parabolic varying width and constant height (I(ξ)=I0(1+aξ)2) with continuous elastic restraint (β=kL4/EI0).
BC
a
=
-
0.5
a
=
0.5
a
=
1
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
C-F
4.341228
5.715021
7.305691
7.704323
11.199778
15.954946
8.844363
13.022270
19.259820
F-C
4.814955
6.692851
8.861021
8.038782
10.842500
14.737209
10.168566
13.021924
17.364909
P-P
7.672218
10.127971
14.937938
17.585157
20.097766
25.116851
23.271880
25.748067
30.688871
C-P
12.293759
14.009160
17.265144
32.912173
35.092823
39.430646
44.369842
46.626819
51.129119
P-C
12.782280
15.020605
19.433707
32.627483
34.516476
38.249232
43.893645
45.666325
49.175037
C-C
22.316590
24.131698
27.694119
61.849231
63.721127
67.446145
83.766543
85.604838
89.266360
C-G
10.832136
14.608627
18.532817
23.114885
30.762149
44.091245
29.276062
37.551999
53.151712For the case ofp=3, column has constant width and linearly varying height. The numerical results for nondimensional critical buckling loads, α=PL2/EI0, are given in Table 7 for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. Looking at Table 7, it is seen that the increase in critical buckling loads for a>0 and the decrease in critical buckling loads for a<0 are more pronounced as compared to p=2 case because this time moment of inertia term is affected by a factor of 3. In Section 4.4, a column that has linearly varying width and height is also considered with varying elastic modulus.Table 7
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of constant width and linearly varying height (I(ξ)=I0(1+aξ)3) with continuous elastic restraint (β=kL4/EI0).
a
=
-
0.5
a
=
0.5
a
=
1
β
=
0
β
=
25
β
=
50
β
=
100
β
=
0
β
=
25
β
=
50
β
=
100
β
=
0
β
=
25
β
=
50
β
=
100
C-F
1.336426
3.406304
4.350010
5.487807
3.551010
8.394751
12.332922
18.020661
4.612119
10.082221
14.934526
22.864147
F-C
0.576515
4.175688
5.657646
7.558428
5.811363
9.281395
12.095892
16.252824
10.691414
13.774455
16.483421
20.954255
P-P
3.627812
6.022524
8.352143
12.576225
18.375617
20.866994
23.355361
28.321680
29.022499
31.437359
33.847167
38.650332
C-P
7.362241
8.951802
10.460281
13.205198
37.498111
39.736262
41.969764
46.421968
58.957688
61.272201
63.584952
68.205009
P-C
7.369711
9.677768
11.969157
16.483968
37.490417
39.306827
41.110477
44.677867
58.897925
60.520678
62.134215
65.333104
C-C
14.348488
16.093293
17.802510
21.098101
73.217187
75.065707
76.908867
80.578722
114.787912
116.546905
118.302355
121.800784
C-G
3.757964
8.394621
10.934803
13.536789
18.603820
26.955608
35.048357
49.998675
30.063711
38.940072
47.723673
64.917551


### 4.3. Effect of Material Nonhomogeneity in the Axial Direction

Euler-Bernoulli columns with material nonhomogeneity along the axial direction are considered to see the effect of variable elastic modulus on critical buckling loads. For this case moment of inertia term,I, is taken to be constant and elastic modulus varies according to a power law or an exponential low gradient assumption.


#### 4.3.1. Power Law Gradient Assumption

For power law gradient assumption, the variation of elastic modulus along the axial direction is taken as(30)E(ξ)=E0+(E1-E0)ξm,m>0,
where m is a material gradient index and takes the values of 0.25, 0.5, 1, 2, and 4 in the present study. Specifically, two constituent materials are taken to be Aluminum and Alumina. The elastic moduli of these constituents are 70 GPa (E0) and 380 Gpa (E1), respectively [37]. Poisson’s ratio is taken to be constant. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 8 for different material gradient index values, m, and normalized restrained stiffness parameters β=kL4/E0I. It is seen that although the constituent materials are not changed, the critical buckling loads change drastically as material gradient index value changes, that is, increase as m gets smaller values because as m goes to zero, the volume fraction of Alumina becomes dominant.Table 8
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for power law gradient assumption (E(ξ)=E0+(E1-E0)ξm).
BC
β
m
=
0.25
m
=
0.5
m
=
1
m
=
2
m
=
4
C-F
0
10.246957
7.735606
4.817137
3.371243
2.744415
50
19.568013
17.538824
14.864166
12.616815
10.701577
100
27.169221
25.136738
22.408392
18.876263
15.432611
P-P
0
45.748909
38.940772
28.774564
18.792973
13.011352
50
50.803306
43.960898
33.699710
23.651864
17.958863
100
55.857000
48.977666
38.610936
28.483003
22.886075
C-P
0
88.219496
73.141853
56.085968
40.610932
29.196866
50
92.644745
77.725784
60.691265
45.136468
33.591204
100
97.057453
82.298033
65.2873438
49.648452
37.960864
C-C
0
171.978825
141.782221
108.012834
79.801395
61.973263
50
194.878114
164.325089
111.724718
83.269291
65.289531
100
198.515086
167.927968
115.422272
86.713451
68.575996
C-G
0
43.162442
36.512699
26.886961
21.480866
17.319545
50
66.091837
59.354530
44.201739
38.490118
33.745313
100
79.619024
73.035781
60.849096
54.508426
48.551928


#### 4.3.2. Exponential Law Gradient Assumption

For exponential law gradient assumption, the variation of elastic modulus along the axial direction is taken to be(31)E(ξ)=E0eμξ,
where μ is a material gradient index and takes the values of −1, −0.5, 0.5, and 1 in the present study. Negative values of μ correspond to the case; the column is ceramic rich at ξ=0 and metal rich at ξ=1 whereas positive values of μ corresponds to the case; the column is metal rich at ξ=0 and ceramic rich at ξ=1. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 9 for different material gradient index values, μ, and normalized restrained stiffness parameters β=kL4/E0I. As can be seen from Table 9, critical buckling load increases for all of the boundary conditions with an increase in material gradient index as expected and it is also seen that restrained stiffness parameter has an important influence on critical buckling loads.Table 9
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for exponential law gradient assumption (E(ξ)=E0eμξ).
BC
β
μ
=
-
1
μ
=
-
0.5
μ
=
0.5
μ
=
1
C-F
0
1.782102
2.112127
2.844778
3.241181
50
6.272146
7.534741
10.200538
11.507320
100
8.112043
9.940415
14.224745
16.544646
P-P
0
5.826546
7.634493
12.587151
15.838195
50
10.824282
12.684208
17.638042
20.845819
100
15.772528
17.725908
22.685148
25.842687
C-P
0
11.988386
15.639912
25.782548
32.554915
50
15.640533
19.528086
30.050913
36.962919
100
19.113129
23.299501
34.275138
41.344887
C-C
0
23.490038
30.598353
50.448156
63.852545
50
27.172950
34.353512
54.212593
67.556554
100
30.785108
38.059830
57.948273
71.236233
C-G
0
5.972525
7.682708
12.666644
16.235007
50
16.438031
19.863476
27.642188
32.240813
100
21.152920
26.511574
39.345458
46.225436


### 4.4. Effects of Both Material Nonhomogeneity and Variable Cross-Section

To investigate the effects of varying elastic modulus and cross-section along the axial direction together, a column that has linearly varying width and height and linearly varying elastic modulus is considered.Width of the column:(32a)b(ξ)=b0(1-cbξ),0<cb<1.Height of the column:(32b)h(ξ)=h0(1-chξ),0<ch<1.Moment of inertia of the column:(32c)I(ξ)=I0(1-cbξ)(1-chξ)3.Elastic modulus of the column:(32d)E(ξ)=E0(1+ξ),where b0, h0, I0, and E0 are width, height, moment of inertia, and elastic modulus of the column at ξ=0, respectively, cb and ch are geometric parameters that correspond to width and height taper ratios, respectively.The case ofE(ξ)=E0(1+ξ) with linearly varying width and height is investigated by Shahba and Rajasekaran [15] using differential quadrature element method of lowest order (DQEL) for β=0 case, that is, there are no elastic restraints. In [15], 12 elements each consisting of 11 nodes were chosen; thus total number of reference grids were taken to be 121. In DQEL, four compatibility equations should be written at the connection points of the elements and required matrix dimension to solve the problem is 576×576. To be consistent with [15], in LDQM solution of the problem total number of reference grids and neighboring nodes are also taken to be 121 and 11, respectively, but now required matrix dimension to solve the problem is 121×121. The nondimensional critical buckling loads, α=PL2/E0I0, of the columns for C-F, C-C, P-P, and C-P end conditions are given in Tables 10, 11, 12, and 13, respectively. For β=0 case, the results are compared to results of [15] for available boundary conditions and found to be perfectly consistent. The results show the potential of LDQM for solution of generalized eigenvalue problems governed by fourth order varying coefficient ordinary differential equations with high accuracy and less computational effort.Table 10
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for C-F boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
3.1177
2.9497
2.7676
2.5652
2.3285
LDQM
3.117663
2.949668
2.767637
2.565141
2.328495
0.2
Reference  [15]
2.6225
2.4638
2.2915
2.0992
1.8725
LDQM
2.622455
2.463779
2.291515
2.099165
1.872518
0.4
Reference  [15]
2.1054
1.9585
1.7988
1.6200
1.4074
LDQM
2.105415
1.958547
1.798693
1.619887
1.407356
0.6
Reference  [15]
1.5522
1.4217
1.2798
1.1208
0.9309
LDQM
1.552205
1.421546
1.279989
1.120782
0.930901
0.8
Reference  [15]
0.9245
0.8217
0.7109
0.5883
0.4441
LDQM
0.924469
0.821641
0.710897
0.588251
0.444050
50
0
LDQM
10.969897
10.399534
9.730454
8.911626
7.828378
0.2
LDQM
9.235948
8.649698
7.968922
7.148865
6.086033
0.4
LDQM
7.195173
6.631250
5.989832
5.235064
4.284203
0.6
LDQM
4.857766
4.387832
3.866699
3.271051
2.544323
0.8
LDQM
2.344717
2.062718
1.758631
1.421514
1.024155
100
0
LDQM
15.525762
14.520910
13.367622
11.997451
10.254518
0.2
LDQM
12.557084
11.612707
10.545695
9.300364
7.742857
0.4
LDQM
9.374096
8.560820
7.656126
6.617654
5.337232
0.6
LDQM
6.093411
5.484473
4.816329
4.058816
3.136808
0.8
LDQM
2.854117
2.516107
2.150262
1.741823
1.254405Table 11
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for C-C boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
57.3940
51.7856
45.7356
38.9917
30.8922
LDQM
57.393956
51.785589
45.735646
38.991689
30.892202
0.2
Reference  [15]
41.9169
37.6023
32.9638
27.8171
21.6802
LDQM
41.916891
37.602328
32.963838
27.817083
21.680218
0.4
Reference  [15]
28.1794
25.0890
21.7813
18.1332
13.8242
LDQM
28.179398
25.089039
21.781290
18.133162
13.824196
0.6
Reference  [15]
16.3412
14.3958
12.3266
10.0639
7.4275
LDQM
16.341179
14.395830
12.326582
10.063919
7.427477
0.8
Reference  [15]
6.6801
5.7836
4.8399
3.8228
2.6649
LDQM
6.680051
5.783574
4.839860
3.822760
2.664885
50
0
LDQM
61.166910
55.582013
49.557537
42.839820
34.759070
0.2
LDQM
45.734227
41.416605
36.768958
31.599542
25.400286
0.4
LDQM
31.958849
28.829688
25.466443
21.730794
17.253425
0.6
LDQM
19.880964
17.847607
15.659665
13.222303
10.279687
0.8
LDQM
9.367432
8.320721
7.182753
5.897801
4.324708
100
0
LDQM
64.914455
59.351070
53.349176
46.652798
38.579025
0.2
LDQM
49.518379
45.193527
40.530226
35.326471
29.036478
0.4
LDQM
35.683768
32.505243
29.069896
25.216920
20.502165
0.6
LDQM
23.285078
21.132662
18.778906
16.088335
12.689894
0.8
LDQM
11.549700
10.288383
8.886242
7.271719
5.279100Table 12
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for P-P boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
14.5112
13.1398
11.6969
10.1451
8.3957
LDQM
14.511250
13.139786
11.696914
10.145084
8.395671
0.2
Reference  [15]
10.6860
9.5971
8.4543
7.2284
5.8498
LDQM
10.686041
9.597121
8.454289
7.228449
5.849840
0.4
Reference  [15]
7.2831
6.4715
5.6228
4.7164
3.7019
LDQM
7.283103
6.471476
5.622776
4.716370
3.701880
0.6
Reference  [15]
4.3287
3.7892
3.2283
2.6338
1.9748
LDQM
4.328693
3.789166
3.228315
2.633792
1.974761
0.8
Reference  [15]
1.8667
1.5950
1.3157
1.0239
0.7075
LDQM
1.866691
1.595022
1.315685
1.023937
0.707517
50
0
LDQM
19.549179
18.192929
16.760939
15.207995
13.417864
0.2
LDQM
15.752005
14.659979
13.503215
12.238938
10.750059
0.4
LDQM
12.298411
11.454703
10.551443
9.542443
8.289481
0.6
LDQM
9.058173
8.426077
7.721329
6.876642
5.694553
0.8
LDQM
5.246198
4.728938
4.126561
3.402176
2.479665
100
0
LDQM
24.581143
23.242950
21.824406
20.269884
18.419554
0.2
LDQM
20.817957
19.721682
18.544637
17.218652
15.521898
0.4
LDQM
17.286073
16.384686
15.372020
14.125521
12.234643
0.6
LDQM
13.406960
12.499393
11.357450
9.833780
7.688312
0.8
LDQM
6.909430
6.121287
5.250381
4.263196
3.070432Table 13
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (
β=kL4/E0I0) for C-P boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
29.448963
26.600568
23.560586
20.230563
16.374540
0.2
21.553772
19.335410
16.972868
14.391143
11.408619
0.4
14.537718
12.927659
11.218322
9.357427
7.216955
0.6
8.475484
7.447400
6.361291
5.186232
3.845576
0.8
3.499195
3.018030
2.514197
1.975492
1.371273
50
0
33.787588
30.872734
27.742199
24.278066
20.188428
0.2
25.672288
23.367048
20.888218
18.137705
14.868142
0.4
18.300248
16.579568
14.724334
12.655490
10.170102
0.6
11.626112
10.462193
9.200015
7.780077
6.049347
0.8
5.487561
4.864484
4.184411
3.415437
2.479706
100
0
38.091765
35.101502
31.866729
28.245336
23.874818
0.2
29.723039
27.314852
24.695929
21.738031
18.107927
0.4
21.920220
20.060135
18.018363
15.680300
12.748168
0.6
14.463191
13.118252
11.623244
9.887789
7.688601
0.8
6.910808
6.130247
5.269028
4.287546
3.090197


## 4.1. Uniform Homogeneous Column with Elastic Restraints

For this case, (3) has an analytic solution and Wang et al. [8] have derived the stability criteria in nondimensional critical buckling load α=PL2/EI and normalized restrained stiffness parameter β=kL4/EI for the columns with different end conditions. Stability criteria are as follows:P-P Column:(23)sinT=0.C-P Column:(24)TcosTsinS-SsinTcosS=0.C-G Column:(25)TsinTcosS-ScosTsinS=0.C-C Column:(26)2ST[cosTcosS-1]+(T2+S2)sinTsinS=0.C-F Column:(27)[α(S2+T2)-2S2T2]cosTcosS-α(S2+T2)+(S4+T4)+ST[2α-(S2+T2)]sinTsinS=0,where
(28)S=α2-(α2)2-β,T=α2+(α2)2-β.
The stability criteria given in (23)–(27) are extremely nonlinear and finding the smallest root of the equation which is the critical buckling load is not easy for any assumed stiffness parameter. To show the effectiveness of the LDQM, nondimensional critical buckling loads, α=PL2/EI, are calculated for clamped-clamped beams with different normalized restrained stiffness parameters β=kL4/EI and compared to the exact results. In calculations the number of neighboring grids Nx is taken to be 11 and total number of grid points N is changed from 21 to 41. As can be seen from Table 1, critical buckling load converges to exact results as N increases from 11 to 41. In the following, neighboring grids Nx is taken to be 11 and total number of grid points N is taken to be 41 unless otherwise specified.Table 1
Non dimensional critical buckling loads (α=PL2/EI) for clamped-clamped uniform homogeneous columns with continuous elastic restraint (β=kL4/EI).
β
Exact [8]
DQMN=21
DQMN=31
DQMN=41
LDQMNx=11, N=21
LDQMNx=11, N=31
LDQMNx=11, N=41
0
39.47841760
39.47841735
39.46144027
—
39.47841219
39.47841756
39.47841760
50
43.26056590
43.26107726
38.55373074
—
43.26023292
43.26056213
43.26056583
100
47.00660086
47.00980042
37.89095672
—
47.00612534
47.00659571
47.00660072The critical buckling loads of the aforementioned columns have also been calculated by Atay and Coşkun [21] using a variation iteration approach and Huang and Luo [30] using power series and the integral technic. Nondimensional critical buckling loads for other boundary conditions are also obtained using LDQM and compared with the exact results [8], and with [21, 30]. The LDQM results given in Table 2 are in good agreement with the exact and other numerical results. It is also seen from Tables 1 and 2 that although the critical buckling load increases as the normalized restrained stiffness parameter increases for all of the boundary conditions, the column with C-G end condition is more sensitive to normalized restrained stiffness parameter change than the others.Table 2
Non dimensional critical buckling loads (α=PL2/EI) for uniform homogeneous columns with continuous elastic restraint (β=kL4/EI).
BC
β
=
0
β
=
50
β
=
100
Exact [8]
Reference [21]
Reference [30]
LDQM
Exact [8]
Reference [21]
Reference [30]
LDQM
Exact [8]
Reference [21]
Reference [30]
LDQM
C-P
20.1907
20.1908
20.190729
20.190729
24.2852
24.2855
24.285198
24.285198
28.3066
28.3080
28.306631
28.306631
P-P
9.8696
9.8696
9.869604
9.869604
14.9357
14.9357
14.935664
14.935664
20.0017
20.0017
20.001723
20.001723
C-G
9.8696
9.8696
9.869604
9.869604
23.5717
23.5722
23.571661
23.571659
32.6690
32.6482
32.668976
32.668976
C-F
2.4674
2.4674
2.467401
2.467401
8.8614
8.8614
8.861396
8.861396
11.9964
11.9964
11.996413
11.996413


## 4.2. Effect of Variable Cross-Section

Euler-Bernoulli columns with nonuniform cross-section along the axial direction are considered to see the effect of variable cross-section on critical buckling load. For this case elastic modulusE is taken to be constant and moment of inertia term changes according to a power law:
(29)I(ξ)=I0(1+aξ)p,a>-1,
where a is a geometric parameter and p=1, 2, and 3.First let us consider the casep=1 which corresponds to the columns of linearly varying width and constant height. Nondimensional critical buckling loads, α=PL2/EI0, are calculated using LDQM for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. The results are compared to available results in the literature and given in Tables 3 and 4. It is seen from Tables 3 and 4 that with a changes from negative to positive values, nondimensional critical buckling load increases since negative values of a correspond to decreasing width and positive values of a correspond to increasing width as ξ changes from 0 to 1. As normalized restrained stiffness parameter increases, critical buckling loads also increase for all of the boundary conditions as in the case of uniform columns.Table 3
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of linearly varying width and constant height (I(ξ)=I0(1+aξ)) for β=0.
BC
a
Reference [28]
Reference [3]
Reference [30]
LDQM
C-F
−0.5
—
—
—
2.06209226
0.5
—
—
—
2.81006133
1.0
3.12277257
3.11769623
3.11769623
3.11769628
F-C
−0.5
—
—
—
1.55884812
0.5
—
—
—
3.31117993
1.0
—
4.12418445
4.12418445
4.12418442
P-P
−0.5
—
—
—
7.25562477
0.5
—
—
—
12.25041131
1.0
14.50520092
14.51124954
14.51124954
14.51124954
C-P
−0.5
—
—
—
14.73942213
0.5
—
—
—
24.99361144
1.0
29.49596465
29.44896281
29.44896281
29.44896281
P-C
−0.5
—
—
—
14.72448141
0.5
—
—
—
24.99874044
1.0
—
29.47884426
29.47884426
29.47884426
C-C
−0.5
—
—
—
28.69697807
0.5
—
—
—
48.81145815
1.0
57.44524609
57.39395614
57.39395614
57.39395614
C-G
−0.5
—
—
—
7.16883539
0.5
—
—
—
12.19969805
1.0
—
—
—
14.33767077Table 4
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of linearly varying width and constant height (I(ξ)=I0(1+aξ)) with continuous elastic restraint (β=kL4/EI0).
BC
a
=
-
0.5
a
=
0.5
a
=
1
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
C-F
5.320851
7.249834
9.496180
7.002455
10.034833
13.926981
7.575957
10.969773
15.525985
F-C
5.484886
7.762992
10.365666
7.065179
9.787644
13.328385
7.822199
10.641702
14.499668
P-P
9.774589
12.290572
17.308835
14.778728
17.306502
22.359993
17.030858
19.549179
24.581143
C-P
16.663069
18.556531
22.240823
27.120454
29.236109
33.430687
31.622334
33.787588
38.091765
P-C
16.893794
19.045882
23.289167
26.982108
28.949869
32.834737
31.409689
33.326139
37.113062
C-C
30.583452
32.457228
36.164501
50.705246
52.591875
56.342983
59.283541
61.166905
64.914456
C-G
13.693532
18.852382
24.928160
19.887489
26.982550
38.296345
22.292684
29.808724
42.676370p
=
2 case corresponds to the columns of parabolic varying width and constant height. For this case nondimensional critical buckling loads, α=PL2/EI0, are calculated using LDQM for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. The results are compared to available results in the literature and given in Tables 5 and 6. Comparison of results given in Tables 3 and 4 with the results given in Tables 5 and 6 shows that critical buckling loads for p=2 case are higher than p=1 case for a>0 and critical buckling loads for p=2 case are lower than p=1 case for a<0 for all of the boundary conditions since decrease and increase in column width are influenced by a factor of 2 for the case of p=2. It is also seen from Tables 5 and 6 that although C-P and P-C columns have the same critical buckling loads for β=0, they have somewhat different critical buckling loads for β=25, 50, and 100.Table 5
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of parabolic varying width and constant height (I(ξ)=I0(1+aξ)2) for β=0.
BC
a
Reference [28]
Reference [3]
Reference [30]
LDQM
C-F
−0.5
—
—
—
1.68296637
0.5
—
—
—
3.17193009
1.0
3.85763006
3.83637692
3.83637692
3.83637679
F-C
−0.5
—
—
—
0.95909423
0.5
—
—
—
4.40655924
1.0
—
6.73186541
6.73186541
6.73186555
P-P
−0.5
—
—
—
5.19807211
0.5
—
—
—
15.07084413
1.0
20.79163330
20.79228846
20.79228846
20.79228846
C-P
−0.5
—
—
—
10.52729403
0.5
—
—
—
30.72425899
1.0
42.31345100
42.10917612
42.10917612
42.10917612
P-C
−0.5
—
—
—
10.52729403
0.5
—
—
—
30.72425899
1.0
—
42.10917612
42.10917612
42.10917612
C-C
−0.5
—
—
—
20.48084102
0.5
—
—
—
59.97123859
1.0
82.22844561
81.92336388
81.92336364
81.92336379
C-G
−0.5
—
—
—
5.19807212
0.5
—
—
—
15.07084413
1.0
—
—
—
20.79228847Table 6
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of parabolic varying width and constant height (I(ξ)=I0(1+aξ)2) with continuous elastic restraint (β=kL4/EI0).
BC
a
=
-
0.5
a
=
0.5
a
=
1
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
β
=
25
β
=
50
β
=
100
C-F
4.341228
5.715021
7.305691
7.704323
11.199778
15.954946
8.844363
13.022270
19.259820
F-C
4.814955
6.692851
8.861021
8.038782
10.842500
14.737209
10.168566
13.021924
17.364909
P-P
7.672218
10.127971
14.937938
17.585157
20.097766
25.116851
23.271880
25.748067
30.688871
C-P
12.293759
14.009160
17.265144
32.912173
35.092823
39.430646
44.369842
46.626819
51.129119
P-C
12.782280
15.020605
19.433707
32.627483
34.516476
38.249232
43.893645
45.666325
49.175037
C-C
22.316590
24.131698
27.694119
61.849231
63.721127
67.446145
83.766543
85.604838
89.266360
C-G
10.832136
14.608627
18.532817
23.114885
30.762149
44.091245
29.276062
37.551999
53.151712For the case ofp=3, column has constant width and linearly varying height. The numerical results for nondimensional critical buckling loads, α=PL2/EI0, are given in Table 7 for different geometric parameter values, a, and normalized restrained stiffness parameters β=kL4/EI0. Looking at Table 7, it is seen that the increase in critical buckling loads for a>0 and the decrease in critical buckling loads for a<0 are more pronounced as compared to p=2 case because this time moment of inertia term is affected by a factor of 3. In Section 4.4, a column that has linearly varying width and height is also considered with varying elastic modulus.Table 7
Non dimensional critical buckling loads (α=PL2/EI0) for the columns of constant width and linearly varying height (I(ξ)=I0(1+aξ)3) with continuous elastic restraint (β=kL4/EI0).
a
=
-
0.5
a
=
0.5
a
=
1
β
=
0
β
=
25
β
=
50
β
=
100
β
=
0
β
=
25
β
=
50
β
=
100
β
=
0
β
=
25
β
=
50
β
=
100
C-F
1.336426
3.406304
4.350010
5.487807
3.551010
8.394751
12.332922
18.020661
4.612119
10.082221
14.934526
22.864147
F-C
0.576515
4.175688
5.657646
7.558428
5.811363
9.281395
12.095892
16.252824
10.691414
13.774455
16.483421
20.954255
P-P
3.627812
6.022524
8.352143
12.576225
18.375617
20.866994
23.355361
28.321680
29.022499
31.437359
33.847167
38.650332
C-P
7.362241
8.951802
10.460281
13.205198
37.498111
39.736262
41.969764
46.421968
58.957688
61.272201
63.584952
68.205009
P-C
7.369711
9.677768
11.969157
16.483968
37.490417
39.306827
41.110477
44.677867
58.897925
60.520678
62.134215
65.333104
C-C
14.348488
16.093293
17.802510
21.098101
73.217187
75.065707
76.908867
80.578722
114.787912
116.546905
118.302355
121.800784
C-G
3.757964
8.394621
10.934803
13.536789
18.603820
26.955608
35.048357
49.998675
30.063711
38.940072
47.723673
64.917551


## 4.3. Effect of Material Nonhomogeneity in the Axial Direction

Euler-Bernoulli columns with material nonhomogeneity along the axial direction are considered to see the effect of variable elastic modulus on critical buckling loads. For this case moment of inertia term,I, is taken to be constant and elastic modulus varies according to a power law or an exponential low gradient assumption.


### 4.3.1. Power Law Gradient Assumption

For power law gradient assumption, the variation of elastic modulus along the axial direction is taken as(30)E(ξ)=E0+(E1-E0)ξm,m>0,
where m is a material gradient index and takes the values of 0.25, 0.5, 1, 2, and 4 in the present study. Specifically, two constituent materials are taken to be Aluminum and Alumina. The elastic moduli of these constituents are 70 GPa (E0) and 380 Gpa (E1), respectively [37]. Poisson’s ratio is taken to be constant. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 8 for different material gradient index values, m, and normalized restrained stiffness parameters β=kL4/E0I. It is seen that although the constituent materials are not changed, the critical buckling loads change drastically as material gradient index value changes, that is, increase as m gets smaller values because as m goes to zero, the volume fraction of Alumina becomes dominant.Table 8
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for power law gradient assumption (E(ξ)=E0+(E1-E0)ξm).
BC
β
m
=
0.25
m
=
0.5
m
=
1
m
=
2
m
=
4
C-F
0
10.246957
7.735606
4.817137
3.371243
2.744415
50
19.568013
17.538824
14.864166
12.616815
10.701577
100
27.169221
25.136738
22.408392
18.876263
15.432611
P-P
0
45.748909
38.940772
28.774564
18.792973
13.011352
50
50.803306
43.960898
33.699710
23.651864
17.958863
100
55.857000
48.977666
38.610936
28.483003
22.886075
C-P
0
88.219496
73.141853
56.085968
40.610932
29.196866
50
92.644745
77.725784
60.691265
45.136468
33.591204
100
97.057453
82.298033
65.2873438
49.648452
37.960864
C-C
0
171.978825
141.782221
108.012834
79.801395
61.973263
50
194.878114
164.325089
111.724718
83.269291
65.289531
100
198.515086
167.927968
115.422272
86.713451
68.575996
C-G
0
43.162442
36.512699
26.886961
21.480866
17.319545
50
66.091837
59.354530
44.201739
38.490118
33.745313
100
79.619024
73.035781
60.849096
54.508426
48.551928


### 4.3.2. Exponential Law Gradient Assumption

For exponential law gradient assumption, the variation of elastic modulus along the axial direction is taken to be(31)E(ξ)=E0eμξ,
where μ is a material gradient index and takes the values of −1, −0.5, 0.5, and 1 in the present study. Negative values of μ correspond to the case; the column is ceramic rich at ξ=0 and metal rich at ξ=1 whereas positive values of μ corresponds to the case; the column is metal rich at ξ=0 and ceramic rich at ξ=1. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 9 for different material gradient index values, μ, and normalized restrained stiffness parameters β=kL4/E0I. As can be seen from Table 9, critical buckling load increases for all of the boundary conditions with an increase in material gradient index as expected and it is also seen that restrained stiffness parameter has an important influence on critical buckling loads.Table 9
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for exponential law gradient assumption (E(ξ)=E0eμξ).
BC
β
μ
=
-
1
μ
=
-
0.5
μ
=
0.5
μ
=
1
C-F
0
1.782102
2.112127
2.844778
3.241181
50
6.272146
7.534741
10.200538
11.507320
100
8.112043
9.940415
14.224745
16.544646
P-P
0
5.826546
7.634493
12.587151
15.838195
50
10.824282
12.684208
17.638042
20.845819
100
15.772528
17.725908
22.685148
25.842687
C-P
0
11.988386
15.639912
25.782548
32.554915
50
15.640533
19.528086
30.050913
36.962919
100
19.113129
23.299501
34.275138
41.344887
C-C
0
23.490038
30.598353
50.448156
63.852545
50
27.172950
34.353512
54.212593
67.556554
100
30.785108
38.059830
57.948273
71.236233
C-G
0
5.972525
7.682708
12.666644
16.235007
50
16.438031
19.863476
27.642188
32.240813
100
21.152920
26.511574
39.345458
46.225436


## 4.3.1. Power Law Gradient Assumption

For power law gradient assumption, the variation of elastic modulus along the axial direction is taken as(30)E(ξ)=E0+(E1-E0)ξm,m>0,
where m is a material gradient index and takes the values of 0.25, 0.5, 1, 2, and 4 in the present study. Specifically, two constituent materials are taken to be Aluminum and Alumina. The elastic moduli of these constituents are 70 GPa (E0) and 380 Gpa (E1), respectively [37]. Poisson’s ratio is taken to be constant. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 8 for different material gradient index values, m, and normalized restrained stiffness parameters β=kL4/E0I. It is seen that although the constituent materials are not changed, the critical buckling loads change drastically as material gradient index value changes, that is, increase as m gets smaller values because as m goes to zero, the volume fraction of Alumina becomes dominant.Table 8
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for power law gradient assumption (E(ξ)=E0+(E1-E0)ξm).
BC
β
m
=
0.25
m
=
0.5
m
=
1
m
=
2
m
=
4
C-F
0
10.246957
7.735606
4.817137
3.371243
2.744415
50
19.568013
17.538824
14.864166
12.616815
10.701577
100
27.169221
25.136738
22.408392
18.876263
15.432611
P-P
0
45.748909
38.940772
28.774564
18.792973
13.011352
50
50.803306
43.960898
33.699710
23.651864
17.958863
100
55.857000
48.977666
38.610936
28.483003
22.886075
C-P
0
88.219496
73.141853
56.085968
40.610932
29.196866
50
92.644745
77.725784
60.691265
45.136468
33.591204
100
97.057453
82.298033
65.2873438
49.648452
37.960864
C-C
0
171.978825
141.782221
108.012834
79.801395
61.973263
50
194.878114
164.325089
111.724718
83.269291
65.289531
100
198.515086
167.927968
115.422272
86.713451
68.575996
C-G
0
43.162442
36.512699
26.886961
21.480866
17.319545
50
66.091837
59.354530
44.201739
38.490118
33.745313
100
79.619024
73.035781
60.849096
54.508426
48.551928


## 4.3.2. Exponential Law Gradient Assumption

For exponential law gradient assumption, the variation of elastic modulus along the axial direction is taken to be(31)E(ξ)=E0eμξ,
where μ is a material gradient index and takes the values of −1, −0.5, 0.5, and 1 in the present study. Negative values of μ correspond to the case; the column is ceramic rich at ξ=0 and metal rich at ξ=1 whereas positive values of μ corresponds to the case; the column is metal rich at ξ=0 and ceramic rich at ξ=1. The LDQM results for nondimensional critical buckling loads, α=PL2/E0I, are given in Table 9 for different material gradient index values, μ, and normalized restrained stiffness parameters β=kL4/E0I. As can be seen from Table 9, critical buckling load increases for all of the boundary conditions with an increase in material gradient index as expected and it is also seen that restrained stiffness parameter has an important influence on critical buckling loads.Table 9
Non dimensional critical buckling loads (α=PL2/E0I) for the axially functionally graded uniform columns with continuous elastic restraint (β=kL4/E0I) for exponential law gradient assumption (E(ξ)=E0eμξ).
BC
β
μ
=
-
1
μ
=
-
0.5
μ
=
0.5
μ
=
1
C-F
0
1.782102
2.112127
2.844778
3.241181
50
6.272146
7.534741
10.200538
11.507320
100
8.112043
9.940415
14.224745
16.544646
P-P
0
5.826546
7.634493
12.587151
15.838195
50
10.824282
12.684208
17.638042
20.845819
100
15.772528
17.725908
22.685148
25.842687
C-P
0
11.988386
15.639912
25.782548
32.554915
50
15.640533
19.528086
30.050913
36.962919
100
19.113129
23.299501
34.275138
41.344887
C-C
0
23.490038
30.598353
50.448156
63.852545
50
27.172950
34.353512
54.212593
67.556554
100
30.785108
38.059830
57.948273
71.236233
C-G
0
5.972525
7.682708
12.666644
16.235007
50
16.438031
19.863476
27.642188
32.240813
100
21.152920
26.511574
39.345458
46.225436


## 4.4. Effects of Both Material Nonhomogeneity and Variable Cross-Section

To investigate the effects of varying elastic modulus and cross-section along the axial direction together, a column that has linearly varying width and height and linearly varying elastic modulus is considered.Width of the column:(32a)b(ξ)=b0(1-cbξ),0<cb<1.Height of the column:(32b)h(ξ)=h0(1-chξ),0<ch<1.Moment of inertia of the column:(32c)I(ξ)=I0(1-cbξ)(1-chξ)3.Elastic modulus of the column:(32d)E(ξ)=E0(1+ξ),where b0, h0, I0, and E0 are width, height, moment of inertia, and elastic modulus of the column at ξ=0, respectively, cb and ch are geometric parameters that correspond to width and height taper ratios, respectively.The case ofE(ξ)=E0(1+ξ) with linearly varying width and height is investigated by Shahba and Rajasekaran [15] using differential quadrature element method of lowest order (DQEL) for β=0 case, that is, there are no elastic restraints. In [15], 12 elements each consisting of 11 nodes were chosen; thus total number of reference grids were taken to be 121. In DQEL, four compatibility equations should be written at the connection points of the elements and required matrix dimension to solve the problem is 576×576. To be consistent with [15], in LDQM solution of the problem total number of reference grids and neighboring nodes are also taken to be 121 and 11, respectively, but now required matrix dimension to solve the problem is 121×121. The nondimensional critical buckling loads, α=PL2/E0I0, of the columns for C-F, C-C, P-P, and C-P end conditions are given in Tables 10, 11, 12, and 13, respectively. For β=0 case, the results are compared to results of [15] for available boundary conditions and found to be perfectly consistent. The results show the potential of LDQM for solution of generalized eigenvalue problems governed by fourth order varying coefficient ordinary differential equations with high accuracy and less computational effort.Table 10
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for C-F boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
3.1177
2.9497
2.7676
2.5652
2.3285
LDQM
3.117663
2.949668
2.767637
2.565141
2.328495
0.2
Reference  [15]
2.6225
2.4638
2.2915
2.0992
1.8725
LDQM
2.622455
2.463779
2.291515
2.099165
1.872518
0.4
Reference  [15]
2.1054
1.9585
1.7988
1.6200
1.4074
LDQM
2.105415
1.958547
1.798693
1.619887
1.407356
0.6
Reference  [15]
1.5522
1.4217
1.2798
1.1208
0.9309
LDQM
1.552205
1.421546
1.279989
1.120782
0.930901
0.8
Reference  [15]
0.9245
0.8217
0.7109
0.5883
0.4441
LDQM
0.924469
0.821641
0.710897
0.588251
0.444050
50
0
LDQM
10.969897
10.399534
9.730454
8.911626
7.828378
0.2
LDQM
9.235948
8.649698
7.968922
7.148865
6.086033
0.4
LDQM
7.195173
6.631250
5.989832
5.235064
4.284203
0.6
LDQM
4.857766
4.387832
3.866699
3.271051
2.544323
0.8
LDQM
2.344717
2.062718
1.758631
1.421514
1.024155
100
0
LDQM
15.525762
14.520910
13.367622
11.997451
10.254518
0.2
LDQM
12.557084
11.612707
10.545695
9.300364
7.742857
0.4
LDQM
9.374096
8.560820
7.656126
6.617654
5.337232
0.6
LDQM
6.093411
5.484473
4.816329
4.058816
3.136808
0.8
LDQM
2.854117
2.516107
2.150262
1.741823
1.254405Table 11
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for C-C boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
57.3940
51.7856
45.7356
38.9917
30.8922
LDQM
57.393956
51.785589
45.735646
38.991689
30.892202
0.2
Reference  [15]
41.9169
37.6023
32.9638
27.8171
21.6802
LDQM
41.916891
37.602328
32.963838
27.817083
21.680218
0.4
Reference  [15]
28.1794
25.0890
21.7813
18.1332
13.8242
LDQM
28.179398
25.089039
21.781290
18.133162
13.824196
0.6
Reference  [15]
16.3412
14.3958
12.3266
10.0639
7.4275
LDQM
16.341179
14.395830
12.326582
10.063919
7.427477
0.8
Reference  [15]
6.6801
5.7836
4.8399
3.8228
2.6649
LDQM
6.680051
5.783574
4.839860
3.822760
2.664885
50
0
LDQM
61.166910
55.582013
49.557537
42.839820
34.759070
0.2
LDQM
45.734227
41.416605
36.768958
31.599542
25.400286
0.4
LDQM
31.958849
28.829688
25.466443
21.730794
17.253425
0.6
LDQM
19.880964
17.847607
15.659665
13.222303
10.279687
0.8
LDQM
9.367432
8.320721
7.182753
5.897801
4.324708
100
0
LDQM
64.914455
59.351070
53.349176
46.652798
38.579025
0.2
LDQM
49.518379
45.193527
40.530226
35.326471
29.036478
0.4
LDQM
35.683768
32.505243
29.069896
25.216920
20.502165
0.6
LDQM
23.285078
21.132662
18.778906
16.088335
12.689894
0.8
LDQM
11.549700
10.288383
8.886242
7.271719
5.279100Table 12
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (β=kL4/E0I0) for P-P boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
Reference  [15]
14.5112
13.1398
11.6969
10.1451
8.3957
LDQM
14.511250
13.139786
11.696914
10.145084
8.395671
0.2
Reference  [15]
10.6860
9.5971
8.4543
7.2284
5.8498
LDQM
10.686041
9.597121
8.454289
7.228449
5.849840
0.4
Reference  [15]
7.2831
6.4715
5.6228
4.7164
3.7019
LDQM
7.283103
6.471476
5.622776
4.716370
3.701880
0.6
Reference  [15]
4.3287
3.7892
3.2283
2.6338
1.9748
LDQM
4.328693
3.789166
3.228315
2.633792
1.974761
0.8
Reference  [15]
1.8667
1.5950
1.3157
1.0239
0.7075
LDQM
1.866691
1.595022
1.315685
1.023937
0.707517
50
0
LDQM
19.549179
18.192929
16.760939
15.207995
13.417864
0.2
LDQM
15.752005
14.659979
13.503215
12.238938
10.750059
0.4
LDQM
12.298411
11.454703
10.551443
9.542443
8.289481
0.6
LDQM
9.058173
8.426077
7.721329
6.876642
5.694553
0.8
LDQM
5.246198
4.728938
4.126561
3.402176
2.479665
100
0
LDQM
24.581143
23.242950
21.824406
20.269884
18.419554
0.2
LDQM
20.817957
19.721682
18.544637
17.218652
15.521898
0.4
LDQM
17.286073
16.384686
15.372020
14.125521
12.234643
0.6
LDQM
13.406960
12.499393
11.357450
9.833780
7.688312
0.8
LDQM
6.909430
6.121287
5.250381
4.263196
3.070432Table 13
Non dimensional critical buckling loads (α=PL2/E0I0) for the axially functionally graded non-uniform columns with continuous elastic restraint (
β=kL4/E0I0) for C-P boundary condition.
β
c
h
c
b
0
0.2
0.4
0.6
0.8
0
0
29.448963
26.600568
23.560586
20.230563
16.374540
0.2
21.553772
19.335410
16.972868
14.391143
11.408619
0.4
14.537718
12.927659
11.218322
9.357427
7.216955
0.6
8.475484
7.447400
6.361291
5.186232
3.845576
0.8
3.499195
3.018030
2.514197
1.975492
1.371273
50
0
33.787588
30.872734
27.742199
24.278066
20.188428
0.2
25.672288
23.367048
20.888218
18.137705
14.868142
0.4
18.300248
16.579568
14.724334
12.655490
10.170102
0.6
11.626112
10.462193
9.200015
7.780077
6.049347
0.8
5.487561
4.864484
4.184411
3.415437
2.479706
100
0
38.091765
35.101502
31.866729
28.245336
23.874818
0.2
29.723039
27.314852
24.695929
21.738031
18.107927
0.4
21.920220
20.060135
18.018363
15.680300
12.748168
0.6
14.463191
13.118252
11.623244
9.887789
7.688601
0.8
6.910808
6.130247
5.269028
4.287546
3.090197


## 5. Conclusions

In this study an LDQM is proposed and applied to the critical buckling load analyses of axially functionally graded nonuniform columns with elastic restraint. The method can be applied for any type of nonhomogeneity in the axial direction either in cross-section or in material properties. No fictitious points are used and boundary conditions are directly substituted into weighting coefficient matrices. The matrix dimension in the solution procedure reduces drastically compared to other DQ methods since there is no need to write any compatibility equations. The introduced method can easily be extended to 2Dimensional problems.



---

*Source: 793062-2013-07-28.xml*