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

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.

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 length 𝐿 subjected to an axial compressive load 𝑃 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, 𝑥. 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]: 𝑑²𝑑𝑥²𝐸(𝑥)𝐼(𝑥)𝑑²𝑤𝑑𝑥²𝑑+𝑃²𝑤𝑑𝑥²+𝑘𝑤=0,0<𝑥<𝐿, (1) where 𝑤(𝑥) is the lateral displacement, 𝐸(𝑥) and 𝐼(𝑥) are axially varying modulus of elasticity and moment of inertia of the column, respectively, and 𝑘 is the stiffness of the uniformly distributed lateral springs per unit length. Introducing nondimensional coordinate and deflection as 𝑥𝜉=𝐿(𝜉)=𝑤(𝑥),𝑊𝐿, (2) (1) can be written as [30] 𝑑²𝑑𝜉²𝐸(𝜉)𝐼(𝜉)𝑑²𝑊𝑑𝜉²+𝛼^∗𝑑²𝑊𝑑𝜉²+𝛽^∗𝑊=0,0<𝜉<1, (3) where 𝛼^∗=𝑃𝐿² is the normalized critical load and 𝛽^∗=𝑘𝐿⁴ 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 Figure 2 are considered [8]. The explicit expressions for the boundary conditions can be given in terms of 𝑊(𝜉) as follows. Clamped support (C):

 𝑊=0,𝑑𝑊𝑑𝜉=0. (4a) Pin support (P): 𝑑𝑊=0,²𝑊𝑑𝜉²=0. (4b) Free end (F): 𝑑²𝑊𝑑𝜉²𝑑=0,𝐸(𝜉)𝐼(𝜉)𝑑𝑑𝜉²𝑊𝑑𝜉²+𝛼^∗𝑑𝑊𝑑𝜉=0. (4c) Guided end (G): 𝑑𝑊𝑑𝑑𝜉=0,𝐸(𝜉)𝐼(𝜉)𝑑𝑑𝜉²𝑊𝑑𝜉²+𝛼^∗𝑑𝑊𝑑𝜉=0. (4d)

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 𝑓(𝑥) prescribed in a field domain 𝑎≤𝑥≤𝑏. Let 𝑓(𝑥_𝑖) be the function values specified in a finite set of 𝑁 discrete points 𝑥_𝑖(𝑖=1,2,…,𝑁) of the field domain. The 𝑟th-order derivative of 𝑓(𝑥) at any discrete point 𝑥_𝑖 can be written in DQ analog form as 𝑑^𝑟𝑓𝑥_𝑖𝑑𝑥^𝑟=^𝑁_𝑗=1𝑎^(𝑟)_𝑖𝑗𝑓𝑥_𝑗,𝑖=1,2,…,𝑁, (5) where 𝑎^(𝑟)_𝑖𝑗 are the weighting coefficients of the 𝑟th-order derivative of the function 𝑓(𝑥) associated with points 𝑥_𝑖.

In the generalized DQM, the test functions are assumed to be the Lagrange interpolation test functions such as ℓ_𝑗(𝑥)=𝑀(𝑥)𝑥−𝑥_𝑗𝑀⁽¹⁾𝑥_𝑗,𝑗=1,2,…,𝑁, (6) where 𝑀(𝑥)=^𝑁_𝑚=1𝑥−𝑥_𝑚,𝑀⁽¹⁾𝑥_𝑗=𝑑𝑀𝑥_𝑗=𝑑𝑥^𝑁_{𝑚=1,𝑚≠𝑗}𝑥_𝑗−𝑥_𝑚. (7)

The weighting coefficients for the first-order derivative can be obtained using Lagrange interpolation polynomial as follows:𝑎⁽¹⁾_𝑖𝑗=𝑑ℓ_𝑗𝑥_𝑖=𝑀𝑑𝑥⁽¹⁾𝑥_𝑖𝑥_𝑖−𝑥_𝑗𝑀⁽¹⁾𝑥_𝑗,𝑎𝑖,𝑗=1,2,…,𝑁,𝑖≠𝑗,_(1)𝑖𝑖=𝑑ℓ_𝑖𝑥_𝑖𝑑𝑥=−^𝑁_{𝑗=1,𝑗≠𝑖}𝑎⁽¹⁾_𝑖𝑗,𝑖=1,2,…,𝑁, (8a)(8b)and these coefficients give the [𝐴⁽¹⁾] matrix as 𝐴⁽¹⁾=𝑎⁽¹⁾₁₁𝑎⁽¹⁾₁₂⋯𝑎⁽¹⁾_1𝑁𝑎⁽¹⁾₂₁𝑎⁽¹⁾₂₂⋯𝑎⁽¹⁾_2𝑁𝑎⋮⋮⋱⋮_(1)𝑁1𝑎⁽¹⁾_𝑁2⋯𝑎⁽¹⁾_𝑁𝑁. (9) Similarly, the weighting coefficients for the 𝑟th-order derivative can be evaluated using the following equations:𝑎^(𝑟)_𝑖𝑗=𝑑^𝑟ℓ_𝑗𝑥_𝑖𝑑𝑥^𝑟=𝑟𝑎^(𝑟−1)_𝑖𝑖𝑎⁽¹⁾_𝑖𝑗−𝑎^(𝑟−1)_𝑖𝑗𝑥_𝑖−𝑥_𝑗,𝑎𝑖,𝑗=1,2,…,𝑁,𝑖≠𝑗,𝑟≥2,_{(𝑟)𝑖𝑖}=𝑑^𝑟ℓ_𝑖𝑥_𝑖𝑑𝑥^𝑟=−^𝑁_{𝑗=1,𝑗≠𝑖}𝑎^(𝑟)_𝑖𝑗,𝑖=1,2,…,𝑁,𝑟≥2. (10a)(10b)The weighting coefficients for second- and higher-order derivatives can also be evaluated using the following expression: 𝑎^(𝑟)_𝑖𝑗=^𝑁_𝑘=1𝑎^(𝑟−1)_𝑖𝑘𝑎⁽¹⁾_𝑘𝑗,𝑖,𝑗=1,2,…,𝑁,𝑟≥2. (11) These coefficients give the matrix 𝐴^(𝑟) which can be interpreted as follows: 𝐴^(𝑟)=𝑎^(𝑟)₁₁𝑎^(𝑟)₁₂⋯𝑎^(𝑟)_1𝑁𝑎^(𝑟)₂₁𝑎^(𝑟)₂₂⋯𝑎^(𝑟)_2𝑁𝑎⋮⋮⋱⋮_(𝑟)𝑁1𝑎^(𝑟)_𝑁2⋯𝑎^(𝑟)_𝑁𝑁. (12)

When 𝑥 is taken in the interval 𝑥∈[0,𝐿] instead of 𝑥∈[0,1], the right-hand side of (14) must be divided by 𝐿^𝑟 which can be smaller or greater than 1. Then, the weighting coefficient matrix takes the following form: 𝐴^(𝑟)=𝑑𝑑𝑥^𝑟=𝑑^𝑟𝑑𝑥^𝑟=1𝐿^𝑟𝑎^(𝑟)₁₁𝑎^(𝑟)₁₂⋯𝑎^(𝑟)_1𝑁𝑎^(𝑟)₂₁𝑎^(𝑟)₂₂⋯𝑎^(𝑟)_2𝑁𝑎⋮⋮⋱⋮_(𝑟)𝑁1𝑎^(𝑟)_𝑁2⋯𝑎^(𝑟)_𝑁𝑁. (13)

After calculating the matrix 𝐴⁽¹⁾, weighting coefficient matrices for second- and higher-order derivatives can also be calculated using the following formulae: 𝐴^(𝑟)=𝐴⁽¹⁾𝐴^(𝑟−1)=𝐴^(𝑟−1)𝐴⁽¹⁾. (14)

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

In LDQM, firstly neighboring grids of any discrete point 𝑥_𝑖 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 𝑥 (at 𝑥=0) or nearby it and a sixth-order first derivative approximation is used that means 7 neighboring points are needed (𝑁_𝑥=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 𝑥 (at 𝑥=𝐿) 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 𝑥_𝑖 in the solution domain for the previously mentioned cases. Then, the discretization of the first-order derivative of a function 𝑓(𝑥) with respect to space variable, 𝑥, at any discrete point 𝑥_𝑖 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] 𝑑𝑓𝑥_𝑖𝑑𝑥=_{𝑗∈𝑆𝑖}𝑎⁽¹⁾_𝑖𝑗𝑓𝑥_𝑗,𝑖=1,2,…,𝑁, (15) where 𝑆_𝑖 represents the corresponding set of the neighboring nodes for the discrete grid point 𝑥_𝑖 in the domain or at the boundaries, 𝑁 is the total amount of grid points in the direction of 𝑥. Weighting coefficients of the first-order derivative with respect to spatial coordinate 𝑥 can be evaluated using the following equations:

Figure 3 

Discretization of boundary and near boundary reference points for 𝑁𝑥=7 and 𝑁=21.

Figure 4 

Discretization of an interior reference point for 𝑁𝑥=7 and 𝑁=21.

𝑎⁽¹⁾_𝑖𝑗=∏_{𝑚∈𝑆𝑖,𝑚≠𝑖}𝑥_𝑖−𝑥_𝑚𝑥_𝑖−𝑥_𝑗⋅∏_{𝑚∈𝑆𝑖,𝑚≠𝑗}𝑥_𝑗−𝑥_𝑚,𝑖=1,2,…,𝑁,𝑗∈𝑆_𝑖𝑎,𝑖≠𝑗,_(1)𝑖𝑖=−_{𝑗∈𝑆𝑖,𝑗≠𝑖}𝑎⁽¹⁾_𝑖𝑗,𝑖=1,2,…,𝑁. (16a)(16b) Similarly, the discretization of the higher-order derivatives of 𝑓(𝑥) with respect to space variable, 𝑥, at any discrete point 𝑥_𝑖 can be expressed as 𝑑^𝑟𝑓𝑥_𝑖𝑑𝑥^𝑟=_{𝑗∈𝑆𝑖}𝑎^(𝑟)_𝑖𝑗𝑓𝑥_𝑗,𝑖=1,2,…,𝑁,𝑟≥2. (17) The weighting coefficients for the 𝑟th-order derivative at any discrete point 𝑥_𝑖 can be evaluated using the following equations:𝑎^(𝑟)_𝑖𝑗=𝑟⋅𝑎_𝑖𝑗𝑎^(𝑟−1)_𝑖𝑖−𝑎^(𝑟−1)_𝑖𝑗𝑥_𝑖−𝑥_𝑗,𝑖=1,2,…,𝑁,𝑗∈𝑆_𝑖𝑎,𝑖≠𝑗,𝑟≥2,_{(𝑟)𝑖𝑖}=−_{𝑗∈𝑆𝑖,𝑗≠𝑖}𝑎^(𝑟)_𝑖𝑗,𝑖=1,2,…,𝑁,𝑟≥2. (18a)(18b)It is instructive to note that (18a) and (18b) are rewritten forms of (10a) and (10b) in the neighborhood of the reference point 𝑥_𝑖.

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 variable 𝑥_𝑖 by 𝜉_𝑖=𝑖−1𝑁−1,𝑖=1,2,…,𝑁. (19) In the present study, equally spaced sampling points are used.

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

Introducing 𝐷(𝜉)=𝐸(𝜉)𝐼(𝜉) for flexural rigidity, (3) can be expanded into following form: 𝐷(𝜉)𝑑⁴𝑊𝑑𝜉⁴(𝜉)+2𝑑𝐷𝑑𝑑𝜉³𝑊𝑑𝜉³+𝑑²𝐷(𝜉)𝑑𝜉²𝑑²𝑊𝑑𝜉²+𝛼^∗𝑑²𝑊𝑑𝜉²+𝛽^∗𝑊=0. (20) For the numerical computation equation (20) can be discretized using LDQM as follows: 𝐷𝜉_{𝑖𝑗∈𝑆𝑖}𝑎⁽⁴⁾_𝑖𝑗𝑊𝜉_𝑗(𝜉)+2𝑑𝐷𝑑𝜉_{𝜉=𝜉𝑖𝑗∈𝑆𝑖}𝑎⁽³⁾_𝑖𝑗𝑊𝜉_𝑗+𝑑²𝐷(𝜉)𝑑𝜉²_{𝜉=𝜉𝑖𝑗∈𝑆𝑖}𝑎⁽²⁾_𝑖𝑗𝑊𝜉_𝑗+𝛼^∗_{𝑗∈𝑆𝑖}𝑎⁽²⁾_𝑖𝑗𝑊𝜉_𝑗+𝛽^∗𝑊𝜉_𝑖=0,𝑖=1,2,…,𝑁. (21) Boundary conditions given in (4a), (4b), (4c), and (4d) can be discretized at 𝜉=0 using LDQM as follows. Clamped support (C):𝑊𝜉₁=0,_{𝑗∈𝑆𝑖}𝑎⁽¹⁾_1𝑗𝑊𝑥_𝑗=0. (22a) Pin support (P): 𝑊𝜉₁=0,_{𝑗∈𝑆𝑖}𝑎⁽²⁾_1𝑗𝑊𝑥_𝑗=0. (22b) Free end (F): _{𝑗∈𝑆𝑖}𝑎⁽²⁾_1𝑗𝑊𝑥_𝑗(𝜉)=0,𝑑𝐷𝑑𝜉_{𝜉=0𝑗∈𝑆𝑖}𝑎⁽²⁾_1𝑗𝑊𝑥_𝑗+𝐷𝜉_{1𝑗∈𝑆𝑖}𝑎⁽³⁾_1𝑗𝑊𝑥_𝑗+𝛼^∗_{𝑗∈𝑆𝑖}𝑎⁽¹⁾_1𝑗𝑊𝑥_𝑗=0. (22c) Guided end (G): _{𝑗∈𝑆𝑖}𝑎⁽¹⁾_1𝑗𝑊𝑥_𝑗(𝜉)=0,𝑑𝐷𝑑𝜉_{𝜉=0𝑗∈𝑆𝑖}𝑎⁽²⁾_1𝑗𝑊𝑥_𝑗+𝐷𝜉_{1𝑗∈𝑆𝑖}𝑎⁽³⁾_1𝑗𝑊𝑥_𝑗+𝛼^∗_{𝑗∈𝑆𝑖}𝑎⁽¹⁾_1𝑗𝑊𝑥_𝑗=0. (22d)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 (𝑁−2)×(𝑁−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 (𝑁−1)×(𝑁−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 𝛼=𝑃𝐿²/𝐸𝐼 and normalized restrained stiffness parameter 𝛽=𝑘𝐿⁴/𝐸𝐼 for the columns with different end conditions. Stability criteria are as follows: P-P Column: sin𝑇=0. (23) C-P Column: 𝑇cos𝑇sin𝑆−𝑆sin𝑇cos𝑆=0. (24) C-G Column: 𝑇sin𝑇cos𝑆−𝑆cos𝑇sin𝑆=0. (25) C-C Column: 2𝑆𝑇[cos𝑇cos𝑆−1]+𝑇²+𝑆²sin𝑇sin𝑆=0. (26) C-F Column: 𝛼𝑆²+𝑇²−2𝑆²𝑇²cos𝑇cos𝑆−𝛼𝑆²+𝑇²+𝑆⁴+𝑇⁴+𝑆𝑇2𝛼−𝑆²+𝑇²sin𝑇sin𝑆=0, (27)where 𝑆=𝛼2−𝛼2²−𝛽,𝑇=𝛼2+𝛼2²−𝛽. (28) 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, 𝛼=𝑃𝐿²/𝐸𝐼, are calculated for clamped-clamped beams with different normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸𝐼 and compared to the exact results. In calculations the number of neighboring grids 𝑁_𝑥 is taken to be 11 and total number of grid points 𝑁 is changed from 21 to 41. As can be seen from Table 1, critical buckling load converges to exact results as 𝑁 increases from 11 to 41. In the following, neighboring grids 𝑁_𝑥 is taken to be 11 and total number of grid points 𝑁 is taken to be 41 unless otherwise specified.

Table 1 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸𝐼) for clamped-clamped uniform homogeneous columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸𝐼).

The 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 (𝛼=𝑃𝐿²/𝐸𝐼) for uniform homogeneous columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸𝐼).

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 modulus 𝐸 is taken to be constant and moment of inertia term changes according to a power law: 𝐼(𝜉)=𝐼₀(1+𝑎𝜉)^𝑝,𝑎>−1, (29) where 𝑎 is a geometric parameter and 𝑝=1, 2, and 3.

First let us consider the case 𝑝=1 which corresponds to the columns of linearly varying width and constant height. Nondimensional critical buckling loads, 𝛼=𝑃𝐿²/𝐸𝐼₀, are calculated using LDQM for different geometric parameter values, 𝑎, and normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸𝐼₀. 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 𝑎 changes from negative to positive values, nondimensional critical buckling load increases since negative values of 𝑎 correspond to decreasing width and positive values of 𝑎 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 (𝛼=𝑃𝐿²/𝐸𝐼₀) for the columns of linearly varying width and constant height (𝐼(𝜉)=𝐼₀(1+𝑎𝜉)) for 𝛽=0.

Table 4 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸𝐼₀) for the columns of linearly varying width and constant height (𝐼(𝜉)=𝐼₀(1+𝑎𝜉)) with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸𝐼₀).

𝑝=2 case corresponds to the columns of parabolic varying width and constant height. For this case nondimensional critical buckling loads, 𝛼=𝑃𝐿²/𝐸𝐼₀, are calculated using LDQM for different geometric parameter values, 𝑎, and normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸𝐼₀. 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 𝑝=2 case are higher than 𝑝=1 case for 𝑎>0 and critical buckling loads for 𝑝=2 case are lower than 𝑝=1 case for 𝑎<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 𝑝=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 (𝛼=𝑃𝐿²/𝐸𝐼₀) for the columns of parabolic varying width and constant height (𝐼(𝜉)=𝐼₀(1+𝑎𝜉)²) for 𝛽=0.

Table 6 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸𝐼₀) for the columns of parabolic varying width and constant height (𝐼(𝜉)=𝐼₀(1+𝑎𝜉)²) with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸𝐼₀).

For the case of 𝑝=3, column has constant width and linearly varying height. The numerical results for nondimensional critical buckling loads, 𝛼=𝑃𝐿²/𝐸𝐼₀, are given in Table 7 for different geometric parameter values, 𝑎, and normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸𝐼₀. Looking at Table 7, it is seen that the increase in critical buckling loads for 𝑎>0 and the decrease in critical buckling loads for 𝑎<0 are more pronounced as compared to 𝑝=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 (𝛼=𝑃𝐿²/𝐸𝐼₀) for the columns of constant width and linearly varying height (𝐼(𝜉)=𝐼₀(1+𝑎𝜉)³) with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸𝐼₀).

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, 𝐼, 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 𝐸(𝜉)=𝐸₀+𝐸₁−𝐸₀𝜉^𝑚,𝑚>0, (30) where 𝑚 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 (𝐸₀) and 380 Gpa (𝐸₁), respectively [37]. Poisson’s ratio is taken to be constant. The LDQM results for nondimensional critical buckling loads, 𝛼=𝑃𝐿²/𝐸₀𝐼, are given in Table 8 for different material gradient index values, 𝑚, and normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸₀𝐼. 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 𝑚 gets smaller values because as 𝑚 goes to zero, the volume fraction of Alumina becomes dominant.

Table 8 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸₀𝐼) for the axially functionally graded uniform columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸₀𝐼) for power law gradient assumption (𝐸(𝜉)=𝐸₀+(𝐸₁−𝐸₀)𝜉^𝑚).

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) 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, 𝛼=𝑃𝐿²/𝐸₀𝐼, are given in Table 9 for different material gradient index values, 𝜇, and normalized restrained stiffness parameters 𝛽=𝑘𝐿⁴/𝐸₀𝐼. 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 (𝛼=𝑃𝐿²/𝐸₀𝐼) for the axially functionally graded uniform columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸₀𝐼) for exponential law gradient assumption (𝐸(𝜉)=𝐸₀𝑒^𝜇𝜉).

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:𝑏(𝜉)=𝑏₀1−𝑐_𝑏𝜉,0<𝑐_𝑏<1. (32a) Height of the column: ℎ(𝜉)=ℎ₀1−𝑐_ℎ𝜉,0<𝑐_ℎ<1. (32b) Moment of inertia of the column: 𝐼(𝜉)=𝐼₀1−𝑐_𝑏𝜉1−𝑐_ℎ𝜉³. (32c) Elastic modulus of the column: 𝐸(𝜉)=𝐸₀(1+𝜉), (32d)where 𝑏₀, ℎ₀, 𝐼₀, and 𝐸₀ are width, height, moment of inertia, and elastic modulus of the column at 𝜉=0, respectively, 𝑐_𝑏 and 𝑐_ℎ are geometric parameters that correspond to width and height taper ratios, respectively.

The case of 𝐸(𝜉)=𝐸₀(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, 𝛼=𝑃𝐿²/𝐸₀𝐼₀, 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 (𝛼=𝑃𝐿²/𝐸₀𝐼₀) for the axially functionally graded non-uniform columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸₀𝐼₀) for C-F boundary condition.

Table 11 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸₀𝐼₀) for the axially functionally graded non-uniform columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸₀𝐼₀) for C-C boundary condition.

Table 12 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸₀𝐼₀) for the axially functionally graded non-uniform columns with continuous elastic restraint (𝛽=𝑘𝐿⁴/𝐸₀𝐼₀) for P-P boundary condition.

Table 13 

Non dimensional critical buckling loads (𝛼=𝑃𝐿²/𝐸₀𝐼₀) for the axially functionally graded non-uniform columns with continuous elastic restraint ( 𝛽=𝑘𝐿⁴/𝐸₀𝐼₀) for C-P boundary condition.

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.

Acknowledgment

The authors would like to thank Pamukkale University Scientific Research Council for supporting this study under project contract no. 2011BSP013.