Nonlocal Analysis of the Flexural–Torsional Stability for FG Tapered Thin-Walled Beam-Columns

This paper addresses the flexural–torsional stability of functionally graded (FG) nonlocal thin-walled beam-columns with a tapered I-section. The material composition is assumed to vary continuously in the longitudinal direction based on a power-law distribution. Possible small-scale effects are included within the formulation according to the Eringen nonlocal elasticity assumptions. The stability equations of the problem and the associated boundary conditions are derived based on the Vlasov thin-walled beam theory and energy method, accounting for the coupled interaction between axial and bending forces. The coupled equilibrium equations are solved numerically by means of the differential quadrature method (DQM) to determine the flexural–torsional buckling loads associated to the selected structural system. A parametric study is performed to check for the influence of some meaningful input parameters, such as the power-law index, the nonlocal parameter, the axial load eccentricity, the mode number and the tapering ratio, on the flexural–torsional buckling load of tapered thin-walled FG nanobeam-columns, whose results could be used as valid benchmarks for further computational validations of similar nanosystems.


Introduction
Thin-walled beams with open cross-sections (e.g., channel, angle, I-and Tee-sections) carry an extensive variety of potential applications as structural components in various engineering fields (from civil to aeronautical engineering) since they offer high performances with a minimal weight. Moreover, thin-walled beams with varying cross-sections have been of great interest to designers and researches, especially in recent decades. The optimization of weight, the reduction in volume, and the improvement of both strength and stability represent some crucial reasons to increase their use as structural members. Due to the low torsion stiffness, a slender beam with a thin-walled cross-section subjected to an eccentric compressive axial force can buckle in the flexural-torsional mode. Thus, investigations about the stability of tapered thin-walled beams can be very complicated because of the coupled bending and torsional deformations involved, as well as the arbitrary variation in the geometrical properties along the longitudinal direction.
As far as advanced multi-phase composites are concerned, functionally graded materials (FGMs) represent a novel generation of composite materials, based on a smooth and gradual variation in the volume fraction of their constituent phases in any desired direction. Compared to traditional materials and laminated composites, FGMs possess some important advantages, primarily, multifunctionality, a high temperature-withstanding ability, the reduction or total removal of stress concentrations, together with the improved Nanomaterials 2021, 11,1936 2 of 27 strength and fracture toughness. Due to these favorable features, FGMs can represent ideal materials for the design of smart engineering systems and devices, which has motivated their recent extensive use in many engineering applications and modern industries, such as aerospace, automobile, optics, nuclear, electronic and turbine components.
With the recent development of nanotechnology, nanoscaled structural elements such as nanobeams and nanoplates are being widely used as key components in different modern engineering devices, including sensors, actuators, transistors, probes, and nanoelectromechanical systems (NEMS). This requires an appropriate study of the mechanical properties of similar structural systems, with even more complicated natures. The experimental tests demonstrate that classical continuum theories cannot be implemented for the exact analysis of nanostructures, as the size effect can play a significant role in their mechanical behavior. Thus, various higher-order size-dependent continuum theories, such as the modified couple stress theory [1], the surface energy theory [2] and nonlocal elasticity theory [3,4], have been expanded to model small-sized structures. Among these models, the nonlocal elasticity theory, as suggested by Eringen [3], has been widely used in the literature to investigate the stability, deformation and vibrational responses of nanostructural elements, assuming that the stress state at an arbitrary point in a body depends not only on the strain field at that point, but also on the strain fields at all points of the body. At the same time, FGMs have been increasingly applied in small-sized structures due to their superior mechanical properties. In such a context, over the past few years, several investigations have been performed to study the linear and nonlinear mechanical responses of nanosized structures made from homogenous or FGMs. Moreover, a large number of works can be found in the literature focusing on the elastic and/or inelastic static, vibration and instability behavior of beams with a thin-walled cross-section, due to their vast relevance in many engineering configurations. Among the most relevant works on the topic, Kitipornchai and Trahair [5] determined the flexural-torsional critical force of doubly-and/or singly-symmetric I-beams with a geometrical variation under non-uniform torsion. Wekezer [6,7] studied the stability of thin-walled beams with varying open sections based on shell theory strain tensors. Considering the influence of geometric nonlinearity, a finite element technique was suggested by Yang and Yau [8] to assess the buckling behavior of doubly symmetric tapered I-beams. Bradford and Cuk [9] adopted a novel finite element technique to determine the buckling limit state of web-tapered beams with a mono-symmetric I-section. In another study, web-tapered beams with a Tee-section were probed by Baker [10]. A finite element formulation was also applied by Rajasekaran [11,12] to approximate the linear stability resistance of tapered thin-walled beams. Similarly, a simple finite element solution was presented by Gupta et al. [13] and Ronagh et al. [14] to predict the lateral-torsional resistance of tapered I-beams. With the help of the total potential energy and Hamilton's principle, Chen [15] computed the vibrational properties of thin-walled beams with geometrical variation. An innovative finite element formulation was also proposed by Kim [16] to analyze the lateral-torsional buckling (LTB) and vibration behavior of beams with a tapered I-section under different boundary conditions. The shear deformation effect was also accounted within the formulation of Li [17] for the stability study of beams with a linearly variable cross-section under a compressive axial load. A nonlocal elasticity version was also suggested by Peddieson et al. [18] to elaborate a nonlocal Benoulli/Euler beam model. A semi-inverse approach was then employed by Elishakoff et al. [19] for the vibrational analysis of beams made of axially-inhomogeneous materials, whereas Refs. [20][21][22] represent some further useful contributions to the lateraltorsional stability study of thin-walled beams with doubly-and singly-symmetric I-sections under different boundary conditions. Taking into account small deformations and large displacements, Mohri et al. [23,24] analyzed the nonlinear flexural-torsional behavior of thin-walled beams with arbitrary cross-sections by employing the Galerkin method, while Samanta and Kumar [25] provided a shell finite element solution for the study of the distortional buckling resistance of beams with a singly-symmetric I-section under simply supports.
In the field of nonlocal differential elasticity methodology, Reddy [26] proposed some pioneering analytical solutions for the static, buckling and vibrational analyses of beams by considering different shear deformation theories. Some additional analytical outcomes for cantilever beams with linear tapered section were also presented by Challamel et al. [27]. Wang et al. [28] perused the flexural vibration problem of nano-and microbeams, following the assumptions of the nonlocal elasticity theory of Eringen in conjunction with the Timoshenko beam model. Many further works in the literature have successfully applied the Eringen nonlocal elasticity approach combined with different beam theories and numerical solution methods-see Refs. [29][30][31][32][33][34][35][36][37][38]. Among them, Pradhan and Sarkar [29] studied the deformation, instability and vibrational responses of an Euler-Bernoulli beam with variable geometrical and material properties. Aydogdu [30] derived a generalized nonlocal beam theory for the mechanical analysis of nanosize beams by means of the Eringen elasticity assumptions combined with different beam theories. In the same direction, Civalek and Akgöz [31] studied the free vibrational properties of microtubules, which problem was solved numerically based on the DQM. Danesh et al. [32] determined the equations of motion for the longitudinal vibration of nanorods with tapered cross-sections, and solved them via the DQM.Şimşek and Yurtcu [33] used the Timoshenko beam theory to survey the deformation and buckling capacity of nanobeams with varying materials. McCann et al. [34] studied the lateral buckling resistance of steel beam members under pure bending and with simply-supported ends, in presence of discrete elastic lateral restraints along their axial direction. Following the nonlocal continuum theoretical assumptions, a finite element formulation was suggested by Eltaher et al. [35,36] to assess the size effect on the mechanical response of nanobeams made of FG materials. An Euler-Bernoulli beam model was also proposed by Shahba et al. [37] to compute the critical axial forces and natural frequencies of tapered beams with axially non-homogeneous materials. Within the framework of large torsion, Benyamina et al. [38] developed a nonlinear formulation to analyze the lateral stability and buckling moment of tapered I-section beams under simply-simply supports. Among the different numerical strategies to handle similar problems, Nguyen et al. [39] proposed an approximate methodology to evaluate the critical moment of I-section beams in the presence of discrete torsional bracing. Attard and Kim [40] included the shear deformations to determine the lateral stability equations for isotropic beams with a thin-walled open section. Challamel and Wang [41] employed Bessel functions for an exact computation of the lateral-torsional buckling load of strip cantilever beam members subjected to an arbitrary loading distribution. A modified couple stress theory was differently combined with the first-order shear deformable beam model of Ke et al. [42] to describe the size effect on the dynamic stability of microbeams made of FGMs. A novel finite element solution was proposed by Borbon [43] to study the coupled vibrational responses of beams with non-symmetric thin-walled cross-sections, accounting for the possible influence of loading eccentricities, shear deformation and rotatory inertia, and a further approximate methodology was successfully introduced by Serna et al. [44] to study the elastic flexural buckling of non-uniform columns subjected to arbitrary axial forces. Akgoz and Civalek [45] surveyed the free vibrational problem of axially functionally graded (AFG) non-uniform microbeams based on a Euler-Bernoulli beam model and modified couple stress theory. In order to exhaustively assess the static and dynamic responses of beams made of FG piezoelectric materials, an improved three-noded beam element was formulated by Lezgy-Nazargah et al. [46], whereas a novel beam finite element was developed by Trahair [47] for the lateral stability analysis of cantilever tapered steel beams. Different examples of nonlocal models and numerical methods can be found in the literature for a large variety of coupled problems and engineering applications. In Refs. [48,49], the authors proposed a Timoshenko beam nonlocal model to assess the free vibrational response of magneto-electro-elastic nanobeams [48], also made of FGMs [49]. The von Kármán geometric nonlinearity was included within a first-order shear deformable beam model by Liu et al. [50] in a nonlocal elasticity context, to evaluate the buckling and post-buckling responses of nanobeams made of piezoelectric materials in thermo-electro-mechanical conditions. A third-order shear deformable beam theory was adopted by Nami et al. [51] for a thermal stability analysis of FG nanoplates. Among tapered member applications, a novel beam finite element was introduced by Mohri et al. [52], together with a large torsion assumption, to estimate the stability resistance of tapered thin-walled beams. A semianalytical procedure based on the Ritz technique was employed by Kuś [53] for analyzing the lateral stability of linearly tapered-web and/or flange doubly-symmetric I-beams. A finite element-based solution was proposed by Pandeya and Singhb [54] to survey the free vibrational behavior of a fixed-free nanobeam with a varying cross-section. According to the Eringen nonlocal theory and Euler-Bernoulli beam model, the nonlinear vibration of AFG nanobeams with a tapered section was exploited by Shafiei et al. [55], and a semianalytical finite strip procedure was implemented by Zhang et al. [56,57] for the study of the stability capacity of bars with an open and closed cross-section under an axial loading condition [56], accounting for the effect of lateral elastic braces on the overall stability response in Ref. [57]. Further studies on the nonlocal vibration, buckling, and post-buckling of size-dependent beams, rods and plates at different scales can be found in [58][59][60][61][62][63][64][65][66][67], both in an analytical and a numerical sense. More specifically, as far as thin-walled structures are concerned, novel efficient models and computational methods have been developed in the literature to treat even more complicated applications. Among the most recent works, a novel optimization methodology was proposed by Maalawi [68] to enhance the vibrational response of thin-walled box beams with varying material properties. An innovative finite element formulation was also suggested by Lezgy-Nazargah [69] based on the theory of a generalized layered global-local beam (GLGB), to carry out an elasto-plastic analysis of thin-walled beams with reduced computational effort. Nguyen et al. [70,71] derived an efficient finite element formulation to investigate the flexural-torsional stability and buckling response of FGM beams with a singly symmetric open section, in the framework of Vlasov's theory. Li et al. [72] applied the method of generalized differential quadrature to rigorously solve the bending, buckling and vibrational problems of AFG beams, accounting for nonlocal strain gradient theoretical assumptions. Moreover, Khaniki et al. [73][74][75][76][77][78][79][80] published several important contributions elated to the static, vibrational and buckling analysis of small-size beams with a constant or variable cross-section, made of homogenous and/or FGMs. A finite element approach was recently developed by Koutoati [81] to assess the static and free vibrations of multilayer composites and FG beams by means of different shear deformation beam theories. Following the first-order shear deformation theory, Glabisz et al. [82] formulated an innovate algorithm to analyze the stability and vibrational problem of nanobeams incorporating different end supports. Within a modified shear deformation theory context, in which it is not essential to use the shear correction factor, the stability and free vibration behavior of FG nanobeams were explored by Ebrahim et al. [83] using the Chebyshev-Ritz method. A double analytical and finite element solution has recently been proposed by Jrad et al. [84] to assess the triply coupled free vibrational responses of thin-walled beams under different boundary conditions. More recently, a third-order shear deformation theory was employed by Arefi and Civalek [85] to check for the static deformation of cylindrical nanoshells made from FG piezoelectric materials supported by a Pasternak elastic foundation.
Among the studies on tapered structures, Osmani and Meftah [86] studied the shear deformation effect on the buckling response of tapered I-shape beams under different loading conditions. An innovative methodology based on the classical energy approach was expounded by Chen et al. [87] for predicting the lateral buckling resistance of Ibeams with simple supports. Achref et al. [88] analytically assessed the higher-order instability loads of beams with thin-walled open cross-sections under different loading conditions by resorting to a classical finite element approach for comparative purposes. Different numerical approaches were applied in Refs. [89][90][91][92][93] for the linear stability and free vibrational study of homogenous and AFG tapered thin-walled beams with an open cross-section, subjected to different boundary conditions and arbitrary loading cases.
Based on the available literature, however, it seems that the flexural-torsional stability of AFG nanobeam-columns with tapered I-section has never been assessed. The current research is moving in this direction, and is aimed at probing the size-dependent buckling properties of AFG tapered nanobeams with a doubly-symmetric thin-walled cross-section, according to Vlasov assumptions. All the mechanical properties in the present work are graded in the longitudinal direction using the power function except, for the Poisson's ratio, wherein the small size effect is taken into account via the Eringen nonlocal elasticity theory. The nonlocal governing equations of the problem, together with the associated boundary conditions, are obtained by implementing the Vlosov model and the energy method, in order to account for the eccentricity effect of of a compressive axial loading from the centroid within the formulation. The DQM is here employed to solve the resulting stability equations in a strong form and to determine the flexural-torsional buckling load. Different numerical examples analyze the effects of several parameters, namely, the constituent volume fractions, tapering ratio, nonlocal parameter and mode number, on the flexural-torsional stability of AFG tapered nanobeams with an I-section subjected to simply supported boundary conditions. The work is organized as follows. After a preliminary description of the theoretical formulation (Section 2), we provide (in Section 3) the basic notions of the DQM, here applied as an efficient tool to solve the problem with reduced computational effort. In Section 4 we present the results from a large parametric investigation aimed at checking the sensitivity of the mechanical response to different input parameters, which is useful for design purposes. The main results and concluding remarks are discussed in Section 5.

Problem Definition
The following stability model represents an extension of the formulation proposed in Ref. [94] for non-prismatic thin-walled nanobeam-columns with an arbitrary distribution of the material properties in the axial direction, whose numerical outcomes could be useful for the development and design of thin-walled structures, such as scanning tunneling microscopes with nonuniform nanobeams at tunneling tips. Due to the rapid development of nanoscience, the stability of FG nanobeams with variable thin-walled cross sections represents one of their key design benefits, as here explored theoretically via nonconventional Eringen nonlocal elasticity, and numerically via the DQM.

Kinematics
Consider a straight tapered doubly symmetric I-beam made of non-homogeneous material, with variable properties along its longitudinal direction, as represented in Figure 1.
I-beams with simple supports. Achref et al. [88] analytically assessed the higher-order instability loads of beams with thin-walled open cross-sections under different loading conditions by resorting to a classical finite element approach for comparative purposes. Different numerical approaches were applied in Refs. [89][90][91][92][93] for the linear stability and free vibrational study of homogenous and AFG tapered thin-walled beams with an open cross-section, subjected to different boundary conditions and arbitrary loading cases.
Based on the available literature, however, it seems that the flexural-torsional stability of AFG nanobeam-columns with tapered I-section has never been assessed. The current research is moving in this direction, and is aimed at probing the size-dependent buckling properties of AFG tapered nanobeams with a doubly-symmetric thin-walled cross-section, according to Vlasov assumptions. All the mechanical properties in the present work are graded in the longitudinal direction using the power function except, for the Poisson's ratio, wherein the small size effect is taken into account via the Eringen nonlocal elasticity theory. The nonlocal governing equations of the problem, together with the associated boundary conditions, are obtained by implementing the Vlosov model and the energy method, in order to account for the eccentricity effect of of a compressive axial loading from the centroid within the formulation. The DQM is here employed to solve the resulting stability equations in a strong form and to determine the flexural-torsional buckling load. Different numerical examples analyze the effects of several parameters, namely, the constituent volume fractions, tapering ratio, nonlocal parameter and mode number, on the flexural-torsional stability of AFG tapered nanobeams with an I-section subjected to simply supported boundary conditions. The work is organized as follows. After a preliminary description of the theoretical formulation (Section 2), we provide (in Section 3) the basic notions of the DQM, here applied as an efficient tool to solve the problem with reduced computational effort. In Section 4 we present the results from a large parametric investigation aimed at checking the sensitivity of the mechanical response to different input parameters, which is useful for design purposes. The main results and concluding remarks are discussed in Section 5.

Problem Definition
The following stability model represents an extension of the formulation proposed in Ref. [94] for non-prismatic thin-walled nanobeam-columns with an arbitrary distribution of the material properties in the axial direction, whose numerical outcomes could be useful for the development and design of thin-walled structures, such as scanning tunneling microscopes with nonuniform nanobeams at tunneling tips. Due to the rapid development of nanoscience, the stability of FG nanobeams with variable thin-walled cross sections represents one of their key design benefits, as here explored theoretically via nonconventional Eringen nonlocal elasticity, and numerically via the DQM.

Kinematics
Consider a straight tapered doubly symmetric I-beam made of non-homogeneous material, with variable properties along its longitudinal direction, as represented in Figure 1.  The orthogonal right-hand Cartesian coordinate system (x, y, z) is adopted, wherein x denotes the longitudinal axis, and y and z are the first and second principal bending axes parallel to the flanges and web, respectively. The origin O of these axes is located at the centroid of the cross-section. In the current work, it is assumed that the height of the web and/or width of both flanges can vary linearly along the longitudinal direction (x-axis), while the thickness remains constant. In the case of doubly-symmetric thin-walled sections, the shear center coincides with the centroid. In this study, we consider only slender beams, such that shear deformations can be ignored in our formulation, together with the local and distortional deformations. Based on these assumptions and following the Vlasov model for non-uniform torsion [95], the displacement field for an arbitrary point on the beam can be expressed as In these equations, U is the axial displacement, V and W represent the lateral and vertical displacements (along the yand z-directions, respectively); u,v,w are the kinematic quantities defined at the reference surface; ω(y, z) stands for the warping function for the variable cross-section, which can be defined based on St. Venant torsional theory, and θ is the twisting angle. The Green strain tensor components in the large displacement include both the linear and the nonlinear strain parts, as follows where ε l ij denotes the linear part, and ε * ij refers to the quadratic nonlinear part. For thinwalled beams, the strain tensor components reduce to the following:

∂U ∂y
By using Equations (1)-(3) and considering a tapering geometry, the non-zero linear and nonlinear parts of the strain displacement field are defined as where r 2 = y 2 + z 2 . In this study, we consider a compressive axial load P acting at the end of the beam along the z-direction, together with an external bending moment acting around the major principal axis, M * y , while assuming a null bending moment M * z with respect to the z-axis. The most common cases of normal and shear stress associated with the external bending moment M * y and shear force V z are considered as where τ 0 xz is the mean value of the shear stress, σ 0 xx stands for the initial normal stress in the cross-section, and A and I y are the cross-sectional area and second moment of inertia around the y-axis, defined as follows:

Constitutive Relations
According to the Eringen nonlocal elasticity model [4], the stress at a point inside a body depends not only on the strain state at that point, but also on the strain states at all other points throughout the body. For homogenous and isotropic elastic solids, the nonlocal stress tensor σ at point x can be defined as where ε kl and C ijkl denote the linear strain components and the elastic stiffness coefficients, respectively. In addition, α(|x − x|, τ) is the nonlocal kernel function, |x − x| is the Euclidean distance, τ = e 0 a/l stands for the material parameter, where a is an internal characteristic length (e.g., lattice parameter, C-C bond length or granular distance) and l is an external characteristic length in the nanostructures (e.g., crack length, wavelength), and e 0 is a material constant, which is determined experimentally or in an approximate form by matching the dispersion curves of plane waves with those based on atomic lattice dynamics. It is possible to express the integral constitutive equation presented in Equation (7) in the following differential constitutive equation: where ∇ 2 is the Laplacian operator and µ = (e 0 a) 2 stands for the nonlocal parameter. For a nonlocal AFG I-beam, the nonlocal constitutive relations can be written as where E and G are the elastic and shear moduli, respectively, and σ xx , τ xy , and τ xz denote the Piola-Kirchhoff stress tensor components.

Equilibrium Equations
The principle of minimum total potential energy is applied to obtain the equilibrium equations together with the boundary conditions. For thin-walled beams, the total potential energy Π is expressed in its variational form by means of the elastic strain energy U l and the strain energy due to initial stress U 0 , Nanomaterials 2021, 11, 1936 8 of 27 Note that in a linear stability context, in the absence of an external force, the external work associated with the applied loads W e is equal to zero. At the same time, the variational form of the strain energy δU l is defined as where L and A stand for the element length and cross-sectional area, respectively, and δε l xx , δγ l xz and δγ l xy are the linear parts of the strain tensor in a variational form. By substituting Equation (4a-c) into Equation (11), the virtual elastic strain energy becomes By integration over the cross-sectional area, we get where N is the axial force, M y and M z denote the two bending moments, B ω is the bi-moment, and M sv is the St. Venant torsional moment. These stress resultants in Equation (13) are defined as Moreover, the variation in the strain energy due to the initial stresses can be stated as By introducing the first variation in the nonlinear strain-displacement relations, defined by Equation (4d-f), and the initial stresses (5a,b) in Equation (15), we get the following relation: At this stage, by integrating Equation (16) over the cross-section, the variation in the strain energy due to the initial stresses takes the following final form: or equivalently In Equation (18), I z is the second moment of inertia around the z-axis and r 0 is the polar radius gyration around the centroid, given by By introducing Equations (13) and (18) into Equation (10) and setting the coefficients of δu 0 , δv, δw, δθ, as to zero, we obtain the equilibrium equations under the following boundary conditions By substituting Equation (4a-c) into Equation (9) and the subsequent results into Equation (14), the stress resultants are obtained as In the previous expressions, J and I ω are the St. Venant torsion and warping constants, defined as This study is established in the context of small displacements and deformations. According to the linear stability, the nonlinear terms are also disregarded in the equilibrium equations. Based on these assumptions, the system of equilibrium equations for tapered I-beams under a nonlocal theory are finally derived by placing Equation (21) into Equation (19) EAu The related boundary conditions at the ends of the thin-walled nanobeam can be expressed as In the following section, a numerical solution procedure based on the DQM is applied to solve the governing equations for the flexural-torsional buckling of AFG nanobeams with varying I-sections, as has been successfully carried out in the literature for a large variety of problems [96][97][98][99][100][101][102][103].

Numerical Solution Method
Due to the varying cross-sectional mechanical properties, the resulting flexuraltorsional stability Equation (23a-d) for I-tapered nanobeams represent a system of threecoupled fourth-order differential equations with variable coefficients. Under these conditions, it is not possible to accurately estimate a general and straightforward closed-form solution. For such complicated problems, the DQM-based approach, as proposed for the first time by Bellman and Casti [96], is here employed as an efficient and easy tool to solve the coupled differential equations of the problem in a strong form. The basic concept of the proposed method relies on the possibility of discretizing the derivatives of a function with respect to a variable in differential equations at some fixed collocation points by means of a weighted linear summation of the function's values at its adjacent points. The governing equations, together with the associated boundary conditions, are thus transformed into a set of linear algebraic equations, which can be solved with the aid of a computational algorithm to derive an approximate solution for continuous differential equations. To this end, it is necessary to divide the computational region into a fixed number of grid points spanning the solution domain. The accuracy of this numerical approach depends on the number and types of selected sampling points, as also discussed in Refs. [97][98][99][100][101][102][103]. One of the best options for the sampling points in the stability and vibration analysis is the Chebyshev-Gauss-Lobatto points: where N is the total number of grid points in the longitudinal direction. According to DQM, the m th -order derivative of a function f (ξ) at a fixed grid point ξ i can be approximated as where f (ξ j ) refers to the functional value at grid points ξ j (i = 1, 2, . . . , N), and A (m) ij is the weighting coefficient for the m th -order derivative. The first-order derivative of the weighting coefficient A (1) ij is computed by the following algebraic formulation based on the Lagrangian interpolation polynomials, where The higher-order DQM weighting coefficients can be acquired from the first-order ones, as follows: In order to solve the stability equation by means of the differential quadrature approach, a dimensionless variable (ξ = x/L) is introduced. By the expansion of Equation (23), the governing equations of the problem take the following final discrete form: ij v j ) By rewriting the problem in matrix form, we get the following relation, (2) in which and δ jk is the Kronecker delta, defined as In Equation (31), the displacement vectors and the torsion angle vector are defined as The simple form of the final equation, Equation (31), can be stated as or [K] and [K G ] are 3N × 3N matrices, λ is the eigenvalues and {d} is the related eigenvectors. After the implementation of the boundary conditions, we compute the flexural-torsional buckling load from Equation (37), together with the associated vertical and lateral deflections and the twist angles of the AFG nanobeams.

Numerical Examples
In this section, we perform a parametric investigation to assess the sensitivity of the linear stability of AFG thin-walled nanobeam-columns (with a variable I-section and simply supported boundary conditions) to different material properties, as well as to different web and flange tapering parameters, mode numbers, nonlocal parameters, and axial load eccentricities. In what follows, we use the subscripts (•) 0 and (•) 1 to define the mechanical and geometrical properties of beams in their left (x = 0, ξ = 0) and right (x = L, ξ = 1) supports, respectively. The dimensionless buckling load parameter is determined as which accounts for simply supported, tapered beams with I-sections subjected to a compressive axial force. In this regard, it is presumed that the widths of both flanges, b 0 , and the web height, d 0 , of the I-section on the left side increase linearly up to b 1 = (1 + β)b 0 and d 1 = (1 + α)d 0 on the right side ( Figure 2). Thus, the flanges and web tapering ratios are defined as β = b 1 /b 0 − 1 and α = d 1 /d 0 − 1, respectively. Note that these two parameters (α, β) are non-negative variables and can change simultaneously or separately. At the same time, by equating (α, β) to zero, we revert to I-beams with a uniform cross-section. The geometrical schemes and dimensionless parameters are depicted in Figure 2. To perform the flexural-torsional buckling analysis, it is supposed that the compressive axial load is applied at three different positions: the top flange (TF) of the left side (i.e., for x = 0), the centroid, and the TF of the right side (i.e., for x = L).  For the same benchmark, the beams feature axially varying materials, ranging between pure ceramic on the side end and pure metal on the right side, according to a simple power-law function. More specifically, the ceramic phase is made of alumina ( ) 2 3 Al O with an equivalent Young's modulus 0 380 , whereas the metal phase For the same benchmark, the beams feature axially varying materials, ranging between pure ceramic on the side end and pure metal on the right side, according to a simple powerlaw function. More specifically, the ceramic phase is made of alumina (Al 2 O 3 ) with an equivalent Young's modulus E 0 = 380 GPa, whereas the metal phase is aluminum (Al) with an equivalent Young's modulus E 1 = 70 GPa, without considering the exact grain sizes and shapes of each material constituent. This means that the modulus of elasticity at an arbitrary coordinate is defined as where the power-law index m assumes a positive value, and is zero only in a pure metal member. We carry out a preliminary study aimed at defining the appropriate number of grid points within the domain to yield accurate results in terms of flexural-torsional buckling load. In the absence of further numerical nonlocal studies on the same thin-walled examples, the accuracy of our formulation is checked by comparing our results with predictions based on a classical finite element method, as performed via the commercial ANSYS code [104]. In detail, we evaluate the lowest values of the dimensionless buckling parameter (P nor ) for the same structure made of pure alumina with three different loading positions and different tapering ratios, α = β = 0 ÷ 1 by steps of 0.2 versus an increased number of grid points N. The main results are summarized in Table 1, where it seems that a total number of grid points N = 20 is sufficient to obtain the lowest normalized buckling load for different axial load positions and non-uniformity parameters. Based on results in Table 1, we can observe the good agreement between our mathematical DQM-based formulation and predictions made via the ANSYS code [104] for each selected loading case. After the validation phase of the model, we continue with a systematic study of the flexural-torsional buckling of AFG nanobeams with different input parameters, such as eccentric axial load, web and flange non-uniformity parameters, gradient index, mode number and nonlocality parameter. In order to assess the linear stability strength of AFG nanobeams with varying I-sections, we compute the lowest normalized flexural-torsional buckling loads (P nor ) of AFG tapered thin-walled nanobeams subjected to simply supported end conditions, as reported in Table 2, for different tapering ratios (α = β = 0, 0.3, 0.6, 0.9), material compositions (power-law exponent), nonlocal parameters (µ = 0 and 2), and three different loading positions. The compressive axial load can be applied on the TF of the left side (x = 0), at the centroid, and on the TF of the right side (x = L). In Figures 3-5, we represent the variation in P nor depending on Eringen's nonlocal parameters (ranging from 0 to 3) for thin-walled beams with homogenous materials or an FG beam with different gradient indexes m = 0.6, 1.3 and 2, while varying the tapering ratios from 0 to 0.9 and assuming three axial load positions, namely, on the TF for x = 0, on the centroid, and the TF for x = L. In this case, we consider a non-uniform beam with equal web height and flange width tapering ratios, α = β.    In Tables 3-5, we list the magnitude of the normalized flexural-torsional buckling parameter, P nor , for various combinations of web height and flange width tapering ratio, β and α, and nonlocal parameters (µ = 0, 1 and 3) with different non-homogenous indices (m = 0.6, 1.2 and 1.8). The contribution of a possible axial load eccentricity at the crosssection centroid on the buckling resistance is also taken into account. The normalized buckling parameters are respectively illustrated in Tables 3 and 4 for load positions on the TF at x = 0 and on the shear center, as well as in Table 5 for a load position on the TF at     In Tables 3-5, we list the magnitude of the normalized flexural-torsional buckling parameter, nor P , for various combinations of web height and flange width tapering ra-  Table 3. Effect of the power-law exponent and tapering parameter on the normalized flexural-torsional buckling load (P nor ) of simply supported thin-walled nanobeams with different nonlocal parameters (axial load applied on the TF at x = 0).   Table 4. Effect of the power-law exponent and tapering parameter on the normalized flexural-torsional buckling load (P nor ) of simply supported thin-walled nanobeams with different nonlocal parameters (axial load applied on the centroid).   Table 5. Effect of the power-law exponent and tapering parameter on the normalized flexural-torsional buckling load (P nor ) of simply supported thin-walled nanobeams with different nonlocal parameters (axial load applied on the TF at x = l).  Under the first assumed load position, Figure 6 illustrates the variation in the normalized buckling load with respect to the web tapering ratio α and the flange tapering parameter β for different nonlocal parameters and material compositions, i.e., for a pure ceramic and AFG with m = 1. The same analysis is repeated for an axial load located at x = L, whose results are plotted in Figure 7, where the reduction in the buckling load is observable with increasing values of µ and decreased tapering ratios of α and β. Moreover, in Figures 8-10 we plot the lowest buckling load P nor versus the tapering ratio for different values of gradient indexes, m, and nonlocal parameters (µ = 0, 1, 2 and 3), while assuming α = β, and three different positions of compressive load, as in the previous cases. Note that a beam subjected to a compressive axial force on the centroid can buckle in a pure torsional buckling mode. In this context, Table 6 addresses the influence of web and flange tapering ratios, material composition (power-law exponent), and nonlocal parameters (µ = 0, 0.5 and 1) on the normalized torsional buckling load of the tapered nanobeam. Based on the results for both local and nonlocal beams and all non-uniformity ratios, the stability strength improves as the non-homogeneity parameter increases. In other words, a higher flexural-torsional buckling capacity is obtained with an increased power index, m, due to the increased content of ceramic phase and increased gradient index, under a fixed ratio of α = β. As also visible in Figures 8-10, the rate of increase in the critical load with α, β is gradual for increased values of m. Based on a comparative evaluation of the results in Figures 8-10, with the same assumptions for α, β, m and µ, it seems that the load position has a significant effect on the stability strength of nanobeams with varying doubly symmetric I-sections, especially for higher values of the web tapering ratio. As also expected, the highest buckling capacity is obtained when the axial load is located exactly on the centroid, whereas the worst response corresponds to a compressive load applied on the TF owing to the compression of an initial bending moment resulting from a load eccentricity.  Both tables and figures clearly show that the non-uniformity parameter has a remarkable influence on the flexural-torsional buckling load. For each selected power-law exponent, nonlocal parameter and loading position, the stability values of prismatic beams with α = β = 0 and of double-tapered ones with α = β = 1 are the lowest and highest, respectively. By consulting Tables 3-5 and Figures 6 and 7, one can also note that the response is much more affected by the flange tapering parameter β than the web non-uniformity ratio α, since the lowest flexural-torsional buckling modes usually occur with the lowest axis moment of inertia.
For beam-columns subjected to an axial load on the centroid and on the TF for x = 0, it is found that the buckling parameter increases with increasing values of both α and β, due to the enhancement of the cross-section's geometrical properties along with an overall increase in flexural and torsional stiffness in the elastic member (see Tables 3 and 4). The sensitivity of the mechanical response is slightly different for I-beams with an axial load applied on the TF at x = L. As shown in Table 5 and Figure 7, the linear stability strength of beams with a constant flange width tends to reduce with increasing values of α. This is mainly related to the fact that the initial bending moment M * y due to axial load eccentricity is enhanced by increasing the web tapering ratio α. The effect of this phenomenon on the buckling resistance of web and flange tapered beams is quite negligible when all the section walls have the same non-uniformity ratio (i.e., for α = β).   (c) (d)   Table 6. Effect of the power-law exponent and tapering parameter on the normalized torsional buckling load (P nor ) of simply supported thin-walled nanobeams with different nonlocal parameters (axial load applied on the Centroid).  As was also expected, the nonlocal parameter shows a stiffness-softening effect and reduces the buckling strength for all the selected loading positions. Based on the plots in Figures 3-7, it seems that the effect of the Eringen's nonlocal parameter on the buckling response is more pronounced at higher tapering ratios and gradient indexes, especially for beams made of pure ceramic. For example, the normalized buckling load of prismatic members in Table 5 with m = 1.2 decreases by 36.5% when µ increases from 0 to 3. This can be explained by the fact that the flexural and torsional stiffness of simply supported tapered I-beams in a nonlocal theoretical context is inversely proportional to the Eringen's parameter. Usually, the introduction of a nonlocal effect increases the deflection response, or this increase is equivalent to the stiffness reduction in the structural member. Since the linear buckling resistance of beams is directly proportional to their stiffness, a meaningful decrease in the critical load is expected. In Figure 11, we represent the effect of the nonlocality parameter on the first four flexural-torsional buckling loads of nonlocal thinwalled beams with an I-section. In this way, we account for a compressive axial load applied at the TF for x = 0, while considering both an AFG prismatic beam with m = 1 and a homogeneous tapered beam with α = β = 0.5. Based on the plots in Figure 11, it is clear that the nonlocal parameter has more pronounced effects on higher flexural-torsional buckling modes when compared with the lowest ones. It can also be stated that it is necessary to rely on nonlocal theories for an accurate estimation of the flexural-torsional stability limit of nanosized thin-walled beams at higher buckling modes. In addition, it is clearly observable that the nonlocal parameter effect increases when the µ ranges between 0 and 0.9. Eringen's parameter. Usually, the introduction of a nonlocal effect increases the deflection response, or this increase is equivalent to the stiffness reduction in the structural member. Since the linear buckling resistance of beams is directly proportional to their stiffness, a meaningful decrease in the critical load is expected. In Figure 11, we represent the effect of the nonlocality parameter on the first four flexural-torsional buckling loads of nonlocal thin-walled beams with an I-section. In this way, we account for a compressive axial load applied at the TF for 0 x = , while considering both an AFG prismatic beam with 1 m = and a homogeneous tapered beam with 0.5 α β = = . Based on the plots in Figure 11, it is clear that the nonlocal parameter has more pronounced effects on higher flexural-torsional buckling modes when compared with the lowest ones. It can also be stated that it is necessary to rely on nonlocal theories for an accurate estimation of the flexural-torsional stability limit of nanosized thin-walled beams at higher buckling modes. In addition, it is clearly observable that the nonlocal parameter effect increases when the μ ranges between 0 and 0.9.

Conclusions
In this paper, we explore the flexural-torsional buckling of AFG nanobeams with a varying I-section by resorting to the Vlasov model and Eringen's nonlocal elasticity theory. The material properties vary in the axial direction of structural elements according to the power-law distribution of the material constituents. The principle of minimum potential energy is applied to determine the governing equilibrium equations and boundary conditions for AFG tapered thin-walled nanobeams subjected to eccentric axial loads. The governing equations of the problem are implemented and solved numerically by means of the DQM in order to determine the flexural-torsional buckling load. A broad systematic investigation checks for the influence of some important parameters, including the web and flange tapering ratios, the nonlocal parameter, the mode number, the axial load eccentricity and the non-homogeneity index, on the overall response of doubly symmetric tapered nanobeams subjected to simply supported boundary conditions. For all the selected loading positions, it is found that the flexural-torsional buckling capacity of nanobeams with a tapered I-section decreases as the nonlocal parameter increases, whereas the buckling load increases as the flange tapering ratio and ceramic phase,

Conclusions
In this paper, we explore the flexural-torsional buckling of AFG nanobeams with a varying I-section by resorting to the Vlasov model and Eringen's nonlocal elasticity theory. The material properties vary in the axial direction of structural elements according to the power-law distribution of the material constituents. The principle of minimum potential energy is applied to determine the governing equilibrium equations and boundary conditions for AFG tapered thin-walled nanobeams subjected to eccentric axial loads. The governing equations of the problem are implemented and solved numerically by means of the DQM in order to determine the flexural-torsional buckling load. A broad systematic investigation checks for the influence of some important parameters, including the web and flange tapering ratios, the nonlocal parameter, the mode number, the axial load eccentricity and the non-homogeneity index, on the overall response of doubly symmetric tapered nanobeams subjected to simply supported boundary conditions. For all the selected loading positions, it is found that the flexural-torsional buckling capacity of nanobeams with a tapered I-section decreases as the nonlocal parameter increases, whereas the buckling load increases as the flange tapering ratio and ceramic phase, Al 2 O 3 , increase. The effect of the flange tapering parameter β on the buckling capacity is more pronounced than that related to the web tapering ratio, α. As expected, the buckling capacity reaches its highest value in the absence of all possible eccentricity. In addition, the elastic buckling capacity of homogeneous double-tapered beams decreases as the nonlocal parameter increases, especially when compared to AFG prismatic I-beams. The small-scale effects become even more pronounced at higher buckling modes, such that they cannot be clearly disregarded when accurately defining the problem. In its present state, the formulation does not consider the grain sizes or the shapes of the alumina or aluminum components, but this will be considered in a more extended formulation that will include possible material anisotropies. A further extension of the proposed formulation will include the nonlinear effects on the coupled thermo-mechanical stability of tapered micro/nanosized systems, accounting for the possible presence of porosities and defects, together with different boundary and environmental conditions.