Free Vibration Analysis of Axially Functionally Graded (AFG) Cantilever Columns of a Regular Polygon Cross-Section with Constant Volume

: This paper deals with the transverse free vibration of axially functionally graded (AFG) cantilever columns under the inﬂuence of axial compressive load. The columns possessing a regular polygon in their cross-section are tapered and their material properties vary along the axis of the column. An emphasis is placed on the columns with constant volume for admissible geometries and materials. The governing di ﬀ erential equation of the problem is derived and solved using the direct integral approach in conjunction with the determinant search technique. The obtained results are in good agreement with those in the available literature and computed by ﬁnite element analysis. Numerical examples for the natural frequency and mode shape of the columns are presented to investigate the e ﬀ ects of parameters related to geometrical nonuniformity and material inhomogeneity.


Introduction
In a variety of structural engineering applications, columns are often built as one of the most important main components by which axial compressive forces, one of the main types of external loading, are supported [1]. As a result, over the past few decades, many researchers have devoted a lot of effort to improving the analysis of column structure systems. After the concept of functionally graded material (FGM) was established in 1984 by material scientists in Japan, recently, FGM is usually composited from ceramic and metallic materials, because these composites enhance the advantages of the materials, such as stronger mechanical performance, as well as better thermal resistance [2]. Therefore, FGM has been successfully applied for various engineering applications, such as aerospace, precision machinery, and biomedical structures. Due to the benefits of space utilization, esthetics, safety, optimization, and economy, tapered components are typically used in engineering practices [3]. In particular, taper elements behave differently from the uniform ones because of the variation of the cross-section along the axial coordinate yields effective stress distributions and a strong coupling between the stress resultants. This concept is important for optimizing the structural members and reducing the self-weight. Thus, by adopting tapered components, safe and economical designs are achieved.
In this respect, much research has been undertaken on the above-mentioned subject that deals with functionally graded cantilever columns. Generally, FGMs are divided into two types: laterally functionally graded material (LFGM) in which the mechanical properties, i.e., the Young's modulus and mass density, are composited laterally to the axial coordinate; and axially functionally graded of the material and geometrical properties of AFG cantilever columns on free vibration behaviors are extensively discussed. Section 5 summarizes this study and suggests areas for further study. Figure 1a is an AFG cantilever column with a length and a volume . From a geometrical point of view, the column is tapered, the cross-sectional shape is a regular polygon with an integer side number (≥ 3), as shown in Figure 1c, and the column volume is constant. The axial coordinate is zero at the clamped end, and the circumradius, area, and second moment of the regular polygonal cross-section are denoted by , , and , respectively. The Young's modulus and mass density are represented by and ,,and the flexural rigidity is denoted by (= ). At the clamped end ( = 0), , , and are denoted by , , and . At the free end ( = ), , , and are denoted by , ,AF and . The column is externally subjected to an axial compressive load less than the buckling load at the free end. When the column vibrates, the undeformed column axis depicted by the straight dashed line elastically deforms the mode shape depicted by the solid line in Figure 1a, defined in Cartesian coordinates ( , ). The cross-section of the deformed column is subjected to the dynamic shear force and bending moment , as well as the axial force . The column element shown in Figure 1b is loaded to the transverse inertia force and the rotatory inertia couple , since the column has mass. In this study, the free vibration is a harmonic motion in which each dynamic coordinate is proportional to sin( ). For example, , = sin( ) where (= ) is the transverse deflection, is the angular frequency in motion where the dynamic co0 (= 1,2,3, ⋯ ) is the mode number, and is the time.

Shown in
Using the equations of ∑ = 0 and ∑ ( , ) = 0 established from the free body diagram shown in Figure 1b When the column vibrates, the undeformed column axis depicted by the straight dashed line elastically deforms the mode shape depicted by the solid line in Figure 1a, defined in Cartesian coordinates (x, y). The cross-section of the deformed column is subjected to the dynamic shear force Q and bending moment M, as well as the axial force P. The column element shown in Figure 1b is loaded to the transverse inertia force F I and the rotatory inertia couple T, since the column has mass. In this study, the free vibration is a harmonic motion in which each dynamic coordinate is proportional to sin(ω i t). For example, y x,t = y x sin(ω i t) where y x (= y) is the transverse deflection, ω i is the angular frequency in motion where the dynamic co0 i(= 1, 2, 3, · · ·) is the mode number, and t is the time.
Using the equations of F y = 0 and M (x,y) = 0 established from the free body diagram shown in Figure 1b, the equations of the dynamic equilibrium are obtained as From Equation (2), the first derivative dM/dx is re-arranged as Combining the second derivative d 2 M/dx 2 obtained from Equation (3) and dQ/dx = F I in Equation (1) yields The bending moment M, transverse inertia force F I , and rotatory inertia couple T are expressed as [18,19] where the rotatory inertia index R is defined as From Equation (7), the first derivative dT/dx is obtained as Substituting Equations (6) and (9) into Equation (4) yields Equation (10), and from Equation (5), the second derivative d 2 M/dx 2 is obtained as Equation (11): Using Equations (10) and (11) and re-arranging against d 4 y/dx 4 yields the following equation, or Now consider the boundary conditions of the column ends (x = 0 and x = l). At the clamped end (x = 0), the deflection y and the rotation dy/dx are zero: At the free end (x = l), the bending moment M in Equation (5) and the shear force Q in Equation (3) together with Equations (5) and (7) are zero, that is In the equations presented, the variable functions of R f , A, and ρ are arbitrary, so if each respective function is given, the angular frequency can be determined. Now, to define variable functions mentioned above for the mathematical formulations. First, in order to define the function of Young's modulus E, the modular ratio m is introduced as There are various kinds of functions for Young's modulus E in the literature: linear [2,3,9,12,14,17], trigonometric [4], polynomial [5,[14][15][16], piecewise [6], exponential [7,10], and periodic [14] functions, etc. The linear function is selected in this study, and then the function of E at the coordinate x is expressed as [17] For the variable function of the mass density ρ, it is usual that ρ is equal to E [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17], and then the density ratio is the same of modular ratio m (see Figure 1c). Therefore, the mass density ratio m defined as a ratio of ρ f to ρ c and the linear function of ρ at the coordinate x can be written as For the variable function of the circumradius r, the taper ratio n is introduced as The function of r at the coordinate x is expressed as where F is an arbitrary function of x/l, F = F(x/l), but in terms of geometry, three kinds of taper functions F are selected in this study as follows.
for a linear taper, for a parabolic taper, and for a sinusoidal taper with n 1 = n − 1.
Using the function of r in Equation (20), variable functions of A and I for a k-sided regular polygonal cross-section at the coordinate x are obtained as where the constants of c 1 and c 2 for the regular polygon cross-section are: Using Equations (16) and (23), the variable function of flexural rigidity R f (= EI) at the coordinate x is established: Mathematics 2020, 8, 319 The volume V of the column can be determined as where the constant c 3 by the taper type is: for a linear taper, for a parabolic taper, and for a sinusoidal taper.
To facilitate numerical experiments, the following dimensionless system parameters are introduced: where (ξ, η) are the normalized Cartesian coordinates, λ is the volume ratio, p is the load parameter, b i is the buckling load parameter, and C i is the frequency parameter. Substituting Equations (18), (22), (23) and (25) into Equation (12) and using Equations (28)-(33) yields the fourth order ordinary dimensionless differential equation, or where for a linear taper, for a parabolic taper, and for a sinusoidal taper. Boundary conditions of Equations (13) and (14) in the dimensional form are transformed into the non-dimensional form using Equations (28)-(33), or For clamped end (ξ = 0), For free end (ξ = 1), The above fourth order ordinary differential Equation (34) with boundary conditions, Equations (36) and (37), governs the free vibration of AFG cantilever columns with a regular polygon cross-section and constant volume. In Equation (34), the taper type, side number k, modular ratio m, taper ratio n, volume ratio λ, and load parameter p are input parameters, while C i is the eigenvalue which is calculated with its mode shape (ξ i , η i ), using appropriate numerical methods.

Numerical Methods
Based on the above analyses, a FORTRAN computer program was written to compute the natural frequencies of the cantilever columns. The input column parameters were: (1) the geometrical properties of the circumradii r c and r f with the taper type, side number k(≥ 3), column length l, and column volume V; (2) the material properties, i.e., Young's moduli (E c , E f ) and mass densities ρ c , ρ f ; and (3) the axial compressive load P. These input parameters in the dimensional units can be shifted to the non-dimensional form, i.e., modular ratio m, taper ratio n, volume ratio λ, and load parameter p, as developed in the previous section. For integrating differential equations to calculate the mode shape (ξ i , η i ), the Runge-Kutta method [20], a direct integral method, was used, and for computing the frequency parameter C i , the determinant search method enhanced by the Regula-Falsi method [20] was used. This solution method of calculating eigenvalues, such as the natural frequencies of this study, from the initial value problem, is often used in the available literature [17]. For the sake of clarity, the numerical processes for solving the differential equation can be summarized as follows: (1) Define the column parameters of R, k, m, n, λ, and p with the taper type.
(2) Set a trial frequency C t in Equation (34) as a trial eigenvalue C i . The first trial C t is zero.
(3) Subject the initial boundary conditions of Equation (36) to the differential Equation (34) at ξ = 0 and assume two sets of unknown initial boundary conditions of (η , η ) at ξ = 0 as shown in Table 1, where ( ) is one derivative differential operator, etc. Table 1. Assumed initial conditions* at clamped end (ξ = 0) and obtained trials at free end (ξ = 1). Figures 1-4 in set no. [1] and [2] are assumed values and the subscript ' f ' stands for 'at the free end'.
(4) Integrate Equation (34) from ξ = 0 to ξ = 1 using the Runge-Kutta method. This result gives trial coordinates (η, η , η , η ) at the axial coordinate ξ. (5) Using the results of the two trial coordinates separately obtained by Set [1] and Set [2], the following linear combinations are established: [2] , η = η [1] + cη [2] , η = η [1] + cη [2] , η = η [1] + cη [2] (38) where c is a constant. If the trial C t assumed in Step (2) is a characteristic eigenvalue C i , the boundary conditions expressed in Equation (37) at the free end (ξ = 1) are satisfied for the linear combination equations, that is The above Equation (39) can be written in the matrix form: Since c 0, to satisfy the above equation, the following determinant D must be zero: The first convergence criterion is If this criterion is met, the trial C t is just a characteristic eigenvalue C i , and one should go to Step (7). (6) If not, increase the trial frequency C t = C t + ∆C t and perform Steps (2)-(5). During executions, note the sign of D 1 × D 2 , where D 1 is the determinant of the previous execution and D 2 is the determinant of the present execution. When the sign changes, the eigenvalue of C i lies between C t,1 and C t,2 , where C t,1 is the trial frequency corresponding to D 1 and C t,2 corresponds to D 2 . An advanced trial frequency C t,3 approaching closer to the eigenvalue C i can be obtained by the Regula-Falsi method, a solution method of the non-linear equation: The second convergence criterion is If the above criterion is met, the trial C t,3 is the eigenvalue C i for a given set of column parameters.
(8) In order to obtain higher frequencies, set a new C t = C i + ∆C t and repeat all above steps until the desired number of frequencies is computed. (9) Output the computed results of C i with the corresponding coordinates (η, η , η , η ) and stop calculation.

Numerical Experiments and Discussions
Prior to the numerical experiments, it is important to determine the suitable step size ∆ξ in the Runge-Kutta method for efficiently integrating efficiently differential Equation (34). To determine the appropriate ∆ξ, the convergence analysis was performed by the number of dividing column elements, 1/∆ξ, and its result is shown in Figure 2, where the input column parameters are presented. It has been observed that solutions of C i=1,2,3 with 1/∆ξ = 50 (∆ξ = 0.02) converge to those with 1/∆ξ = 200 (∆ξ = 0.005) within three significant numbers. In this study, all computations with 1/∆ξ = 100 (∆ξ = 0.01) in the parametric study were carried out on a PC without any difficulties.  In the available literature, the closed-form or numerical solutions to this problem are lacking, so that, for verification purpose, the selected results of this study are comparable to those of the generalpurpose software ADINA and those in the authors' previous work [17]. The predicted natural frequencies (= 2 ⁄ ) for = 4 and = ∞ with the linear taper are compared in Table 2. Here, AFGM is composited with pure aluminum (Al) at the clamped end and pure zirconia (ZrO2) at the free end, from which the Young's modulus and mass density can be defined along with using Equations (16) and (18). The column parameters are: = 0, = 1 m, = 0.0177 m 3 , = 0.5; = 70 GPa, = 2700 kg/m 3 (Al); = 140 GPa, = 5400 kg/m 3 (ZrO2); and = 0. From these parameters, natural frequencies = 107.8 Hz are obtained from predicted in this study. Results of this study, ADINA, and reference [17] in Table 2 are in good agreement within a 3.5% error. In these comparisons, the theory, including the numerical method developed herein, is verified.  In the available literature, the closed-form or numerical solutions to this problem are lacking, so that, for verification purpose, the selected results of this study are comparable to those of the general-purpose software ADINA and those in the authors' previous work [17]. The predicted natural frequencies f i (= ω i /2π) for k = 4 and k = ∞ with the linear taper are compared in Table 2. Here, AFGM is composited with pure aluminum (Al) at the clamped end and pure zirconia (ZrO 2 ) at the free end, from which the Young's modulus E and mass density ρ can be defined along with x using Equations (16) and (18). The column parameters are: R = 0, l = 1 m, V = 0.0177 m 3 , n = 0.5; E c = 70 GPa, ρ c = 2700 kg/m 3 (Al); E f = 140 GPa, ρ f = 5400 kg/m 3 (ZrO 2 ); and P = 0. From these parameters, natural frequencies f i = 107.8C i Hz are obtained from C i predicted in this study. Results of this study, ADINA, and reference [17] in Table 2 are in good agreement within a 3.5% error. In these comparisons, the theory, including the numerical method developed herein, is verified. Hereafter, the lowest three frequency parameters (i = 1, 2, 3) C i are computed for the numerical experiments. Also C i with p = 0 are computed in Tables 3-5 since this load case is the most practiced in practical engineering. First, selected analyses were conducted to determine the effects of rotatory inertia index R, side number k, and taper type on C i , and these results are listed in Tables 3-5. Note that the column parameters used in the parametric study are given in each table.    Table 3 shows the effect of the inertia index R on C i , where the volume ratio λ varies from λ = 0.005 to 0.05. From these results, the following findings are observed: (1) C i is always lower with rotatory inertia (R = 1) than without rotatory inertia (R = 0), as expected based on the free vibration analysis of structures made of conventional homogeneous materials [21]; (2) this frequency reduction is magnified by a higher mode number i and larger parameter λ; and (3) in practical column designs, the rotatory inertia couple reduces the frequency less than 0.6% for i = 1, less than 4.3% for i = 2 and less than 11.6% for i = 3, which cannot be negligible. Table 4 shows the effects of side number k on the frequency parameter C i , where the number k varies from k = 3 ∼ 8 and ∞. Hereafter, all numerical results included the rotatory inertia couple (R = 1). The frequency parameter C i with a smaller side number k is larger than C i with larger k. For an illustrative example, the fundamental frequency parameter C 1 of k = 3 (triangle) is 9.9% (1.552/1.412=1.099) larger than C 1 of k = ∞ even though the column volumes V are the same. When the k value increases, C i converges to C i of k = ∞. It is observed that C i for k = 8 (octagon) approaches C i for k = ∞ within 1.22%. From these results and others not shown, the effect of k is greatly reduced from i = 1 (critical mode) for the higher mode.
The effects of taper type on the frequency parameter C i , are shown in Table 5. The fundamental frequency parameter C 1 is larger in the order from the linear to sinusoidal to parabolic taper, while for the higher modes i = 2 and 3, C i are larger in the order from the parabolic to sinusoidal to linear taper. As an illustrative example, C 1 of the linear taper is 5.6% (1.445/1.369 = 1.056) larger than C 1 of the parabolic taper and 2.2% (1.445/1.414 = 1.022) larger than C 1 of the sinusoidal taper. In higher modes, the effect of taper type may be negligible.
The numerical results for the modular ratio m, taper ratio n, volume ratio λ, and load parameter p are presented in Figures 3-8, where the effect of rotatory inertia (R = 1) is included. Note that column parameters used in the parametric study are presented in each figure. Figure 3 presents the frequency curves of C i versus m. The frequency parameters C i decrease as the modular ratio m increases for all mode numbers i. The largest C i occurs at m = 0. Higher decreasing rates of C i are observed for the smaller values of m, particularly for m < 1. Larger values of m lead to smaller reduction rates of C i . The effects of taper type on the frequency parameter , are shown in Table 5. The fundamental frequency parameter 1 is larger in the order from the linear to sinusoidal to parabolic taper, while for the higher modes = 2 and 3, are larger in the order from the parabolic to sinusoidal to linear taper. As an illustrative example, 1 of the linear taper is 5.6% (1.445/1.369 = 1.056) larger than 1 of the parabolic taper and 2.2% (1.445/1.414 = 1.022) larger than 1 of the sinusoidal taper. In higher modes, the effect of taper type may be negligible.
The numerical results for the modular ratio , taper ratio , volume ratio , and load parameter are presented in Figures 3-8, where the effect of rotatory inertia ( = 1) is included. Note that column parameters used in the parametric study are presented in each figure. Figure 3 presents the frequency curves of versus . The frequency parameters decrease as the modular ratio increases for all mode numbers . The largest occurs at = 0. Higher decreasing rates of are observed for the smaller values of , particularly for < 1. Larger values of lead to smaller reduction rates of .    Figure 4 presents the frequency curves of C i versus n. The frequency parameters C i decrease as the modular ratio m increases for all mode numbers i. The characteristics of the frequency curves are similar to those of Figure 3. The largest C i occurs at n = 0 and the value of C i converges to a common value, i.e., C i = 0.   Figure 5 represents the three dimensional curved surface map of (m, n, C 1 ) with respect to the fundamental frequency parameter C 1 in the domain of 0 ≤ m ≤ 5 and 0 ≤ n ≤ 1 for a given set of column parameters of the linear taper, k = 4, λ = 0.03 and p = 0, which is the same as previous Figures 3 and 4. Figure 5 reflects the characteristics of both Figures 3 and 4. The C 1 value decreases as both m and n values increase, and therefore the maximum value of C 1,max = 6.462 occurs at the coordinates (m = 0,n = 0) and the smallest value of C 1,min = 0.626 at the coordinates (m = 5, n = 1), as shown in this figure. If the values m and n are infinitely extended, i.e., zero flexural rigidity R f = 0 at the clamped end (x = 0), the C 1 value converges the minimum value of C 1,min = 0. This means that the cantilever column with the above given set of column parameters has the fundamental frequency parameters C 1 in the range of 6.462 ≤ C 1 ≤ 0.   The frequency parameters decrease as the volume ratio increases. For the lower modes = 1 and 2, the effect of on is negligible, while for the higher mode = 3 not negligible. It is particularly noteworthy that, only in this figure, the angular frequency with the smaller is larger than with larger , since is proportional to √ (= √ 3 ⁄ ), i.e., = (√ √ ⁄ ⁄ ) (see Equation (33)).

Figure 5.
Surface map of (m, n, C 1 ) for linear taper, k = 4, λ = 0.03, and p = 0. Figure 6 presents the frequency curves of C i versus λ. The frequency parameters C i decrease as the volume ratio λ increases. For the lower modes i = 1 and 2, the effect of λ on C i is negligible, while for the higher mode i = 3 not negligible. It is particularly noteworthy that, only in this figure, the angular frequency ω i with the smaller C i is larger than ω i with larger C i , since ω i is proportional to √ λ = l 3 /V , i.e., ω i = ( √ λ E c /ρ c /l)C i (see Equation (33)).

Figure 5.
Surface map of ( , , 1 ) for linear taper, = 4, = 0.03, and = 0. Figure 6 presents the frequency curves of versus . The frequency parameters decrease as the volume ratio increases. For the lower modes = 1 and 2, the effect of on is negligible, while for the higher mode = 3 not negligible. It is particularly noteworthy that, only in this figure, the angular frequency with the smaller is larger than with larger , since is proportional to √ (= √ 3 ⁄ ), i.e., = (√ √ ⁄ ⁄ ) (see Equation (33)). Figure 6. Frequency parameter versus volume ratio curves. Figure 6. Frequency parameter C i versus volume ratio λ curves. Figure 7 shows the frequency curves of C i versus p. The frequency parameters C i decrease as the load parameter p increases. When C i decreases and reaches zero, the column buckles and then becomes static state, i.e., C i = 0. The corresponding p with C i = 0 is the buckling load parameter b i . For an example for the fundamental mode i = 1, the critical (i = 1) buckling load parameter b cr is b cr (= b 1 ) = 0.3562 is marked by in the horizontal p axis, i.e., C 1 = 0. Using this physical phenomenon, the buckling loads B with natural frequencies of zero can be determined [17]. It is particularly noteworthy that after buckling (p > b cr ), C i values are meaningless in practical engineering, since the column has already buckled.
Mathematics 2020, 8, 319 15 of 17 Figure 7 shows the frequency curves of versus . The frequency parameters decrease as the load parameter increases. When decreases and reaches zero, the column buckles and then becomes static state, i.e., = 0. The corresponding with = 0 is the buckling load parameter . For an example for the fundamental mode = 1, the critical ( = 1) buckling load parameter is (= 1 ) = 0.3562 is marked by ∎ in the horizontal axis, i.e., 1 = 0 . Using this physical phenomenon, the buckling loads with natural frequencies of zero can be determined [17]. It is particularly noteworthy that after buckling ( > ) , values are meaningless in practical engineering, since the column has already buckled. Now, consider the effects of column parameters on the vibration mode shapes. Figure 8 shows examples of the mode shapes for the given set of column parameters presented in this figure. In Figure 8a, three lowest ( = 1,2,3) mode shapes for the linear, parabolic, and sinusoidal tapers are shown. In Figure 8b, those for = 1, 2, 3 are shown. Three mode shapes depicted by solid, dashed, and dotted curves, respectively, in Figure 8a,b, are much different from each other and therefore, the Now, consider the effects of column parameters on the vibration mode shapes. Figure 8 shows examples of the mode shapes for the given set of column parameters presented in this figure. In Figure 8a, three lowest (i = 1, 2, 3) mode shapes for the linear, parabolic, and sinusoidal tapers are shown. In Figure 8b, those for m = 1, 2, 3 are shown. Three mode shapes depicted by solid, dashed, and dotted curves, respectively, in Figure 8a,b, are much different from each other and therefore, the effects of the taper type and modular ratio n on the mode shapes are significant. From these mode shapes, the positions of nodes and maximum amplitudes of the free vibrations are understood. Now, consider the effects of column parameters on the vibration mode shapes. Figure 8 shows examples of the mode shapes for the given set of column parameters presented in this figure. In Figure 8a, three lowest ( = 1,2,3) mode shapes for the linear, parabolic, and sinusoidal tapers are shown. In Figure 8b, those for = 1, 2, 3 are shown. Three mode shapes depicted by solid, dashed, and dotted curves, respectively, in Figure 8a,b, are much different from each other and therefore, the effects of the taper type and modular ratio on the mode shapes are significant. From these mode shapes, the positions of nodes and maximum amplitudes of the free vibrations are understood.

Concluding Remarks
This paper presents a unique numerical approach for analyzing the free vibration of AFG cantilever columns subjected to an axial compressive force. The column is tapered, the cross-sectional shape is a regular polygon, and the volume of the column is constant. A mathematical model of governing differential equations for such columns is formulated based on the dynamic equilibrium of the free body diagram of the column element subjected to the rotatory inertia couple as well as the transverse inertia force. For numerical experiments, the linear functions of Young's modulus and mass density are selected and three kinds of taper functions of the column are chosen: the linear, parabolic, and sinusoidal taper. The governing equation is numerically integrated by the direct integral method for computing the mode shape, and the determinant search method enhanced by the Regula-Falsi method is used to compute the natural frequencies. For verification purposes, the predicted natural frequencies are compared with those obtained by the general-purpose software ADINA and reference [17]. Effects of the taper type, side number, modular ratio, taper ratio, and volume ratio on the natural frequencies and mode shapes are discussed. For further study, free vibration analysis of AFG columns supported by various end conditions combined with hinged and clamped ends, not considered in this study, is required.