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

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.


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][3][4][5][6][7][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][13][14][15].Other numerical methods such as differential transform method [15], semi-inverse method [16][17][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.

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]: 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 (1) can be written as [30] where  * =  2 is the normalized critical load and  * =  4  is the normalized restrained stiffness parameter.

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): Pin support (P): Free end (F): Guided end (G): 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][4][5][6][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 crosssection 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.

Numerical Discretization and Method of Solution
where  ()  are the weighting coefficients of the th-order derivative of the function () associated with points   .
Similarly, the weighting coefficients for the th-order derivative can be evaluated using the following equations: The weighting coefficients for second-and higher-order derivatives can also be evaluated using the following expression: These coefficients give the matrix  () which can be interpreted as follows: 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:  After calculating the matrix  (1) , weighting coefficient matrices for second-and higher-order derivatives can also be calculated using the following formulae: (1) ] [ (−1) ] = [ (−1) ] [ (1) ] . (14)

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]  (  )  = ∑ 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 firstorder derivative with respect to spatial coordinate  can be evaluated using the following equations: (1)   ,  = 1, 2, . . ., .
Similarly, the discretization of the higher-order derivatives of () with respect to space variable, , at any discrete point   can be expressed as The weighting coefficients for the th-order derivative at any discrete point   can be evaluated using the following equations: 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 In the present study, equally spaced sampling points are used.
Clamped support (C): Pin support (P): Free end (F): Guided end (G): 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.

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.

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  =  2 / and normalized restrained stiffness parameter  =  4 / for the columns with different end conditions.Stability criteria are as follows: C-P Column: cos  sin  −  sin  cos  = 0.
C-G Column: C-C Column: C-F Column: where 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,  =  2 /, are calculated for clamped-clamped beams with different normalized restrained stiffness parameters  =  4 / 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.
The critical buckling loads of the aforementioned columns have also been calculated by Atay and Cos ¸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.

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: 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,  =  2 / 0 , are calculated using LDQM for different geometric parameter values, , and normalized restrained stiffness parameters  =  4 / 0 .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.
= 2 case corresponds to the columns of parabolic varying width and constant height.For this case nondimensional critical buckling loads,  =  2 / 0 , are calculated using LDQM for different geometric parameter values, , and normalized restrained stiffness parameters  =  4 / 0 .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.
For the case of  = 3, column has constant width and linearly varying height.The numerical results for nondimensional critical buckling loads,  =  2 / 0 , are given in Table 7 for different geometric parameter values, , and normalized restrained stiffness parameters  =  4 / 0 .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.

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.

Power Law Gradient Assumption.
For power law gradient assumption, the variation of elastic modulus along the axial direction is taken as 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 ( 0 ) and 380 Gpa ( 1 ), respectively [37].Poisson's ratio is taken to be constant.The LDQM results for nondimensional critical buckling loads,  =  2 / 0 , are given in Table 8 for different material gradient index values, , and normalized restrained stiffness parameters  =  4 / 0 .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.

Exponential Law Gradient Assumption.
For exponential law gradient assumption, the variation of elastic modulus along the axial direction is taken to be where  is a material gradient index and takes the values of −1, −0.5, 0.5, and 1 in the present study.Negative values 2 ) for  = 0. 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,  =  2 / 0 , are given in Table 9 for different material gradient index values, , and normalized restrained stiffness parameters  =  4 / 0 .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.

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: Height of the column: Moment of inertia of the column: Elastic modulus of the column: where  0 , ℎ 0 ,  0 , and  0 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 () =  0 (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,  =  2 / 0  0 , 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.

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

Figure 3 :
Figure 3: Discretization of boundary and near boundary reference points for   = 7 and  = 21.
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.