Analysis of the Axial Vibration of Non-Uniform and Functionally Graded Rods via an Analytical-Based Numerical Approach

Abstract

In this study, an analytical-based numerical approach was proposed for the analysis of the free axial vibration of homogeneous and functionally graded rods with varying cross-sectional areas. The proposed approach is based on analytical approximation techniques, such as the Adomian decomposition method, variational iteration method, and homotopy perturbation method. However, the governing equations of the problems solved in this study were variable coefficient differential equations. These equations provide analytical solutions for strictly limited cases. Analytical approximation methods easily handle problems with uniform material properties and constant cross-sections, whereas with varying cross-sectional areas, the analytical integration process becomes a difficult task for the software. If the rod’s material is functionally graded with varying cross-sectional areas, the analytical integration process becomes a cumbersome task. The proposed approach eliminates all difficulties and requires computation within several seconds. The application of this method is straightforward, and the results obtained in this study are in excellent agreement with the solutions provided in the literature.

1. Introduction

The vibration of rods is of great interest and is still receiving considerable attention from researchers in structural, mechanical, and aeronautical engineering. The longitudinal vibration of rods with constant cross-sections has an exact well-known solution considering different boundary conditions. However, if the cross-sectional area varies with the length of the rod, the governing equation becomes a variable coefficient differential equation for which the exact solutions are very limited. In the case of functionally graded rods, additional functions are added to the governing equation, which makes obtaining a solution difficult.

Raman [1] used transformations to the equation of motion to transform it into a form similar to the one-dimensional Schroedinger equation in order to obtain the analytical solutions of several rods with variable cross-sections. This study examines a process that converts the Sturm–Liouville equation into a specific form of the equation of motion under study. By retracing the analysis, the shape of the cross-section of the rod can be identified for any second-order differential equation with a guaranteed solution existence.

Eisenberger [2] developed a technique called the exact element method, in which the dynamic axial stiffness of the rod is used to obtain exact solutions for the longitudinal vibration of variable cross-section rods. This technique is suitable for any polynomial variation in the cross-sectional area and mass distribution of a member. A comparison of the results of various examples is also available.

Abrate [3] transformed the equation of motion into a wave equation and determined the axial vibration frequencies using the Rayleigh–Ritz method. The free vibration of nonuniform beams with arbitrary boundary conditions and general shapes was examined. Simple equations are presented to determine the fundamental natural frequencies of beams with different end-support conditions.

Bapat [4] studied the vibration of rods composed of uniformly tapered sections with nonclassical boundary conditions. An efficient method was proposed to solve the vibration problem of rods consisting of N uniformly tapered sections with nonclassical boundary conditions. The proposed technique integrates the transfer matrix technique with the closed-form solution of the uniformly tapered rod, which results in a singular equation with only one unknown for every uniformly tapered section, requiring only one matrix multiplication for each section.

Kumar and Sujith [5] provided exact solutions for the longitudinal vibrations of non-uniform rods with specific cross-sectional area functions by applying a transformation to the governing equation. The study utilized special functions, including Bessel, Neumann, and trigonometric functions, to obtain solutions. Simple formulas were provided to predict the natural frequencies of non-uniform rods with varying end conditions, and the dependence on the taper was also discussed.

Li [6] presented an exact solution approach for the free longitudinal vibrations of one-step non-uniform rods with classical and non-classical boundary conditions using an appropriate functional transformation, where the distribution of mass is arbitrary, and the distribution of longitudinal stiffness is expressed as a functional relation with the mass distribution. The presented approach simplifies the governing differential equations for the free vibrations of one-step non-uniform rods, resulting in solvable differential equations for various functional relationships between the stiffness and mass.

Li [7] analyzed the longitudinal vibration of stepped non-uniform rods by describing the distribution of mass as arbitrary, and the distribution of longitudinal stiffness was expressed as a functional relation with respect to mass distribution. The differential equations governing the longitudinal free vibration of rods with variable cross-sections were reduced to Bessel’s equations or other analytically solvable differential equations by selecting appropriate expressions, such as power functions and exponential functions, for functional relations. Simple formulas were proposed to predict the longitudinal vibration frequencies and mode shapes of one-step rods with continuously varying cross sections. The transfer matrix technique and closed-form solutions of one-step non-uniform rods were integrated to generate a single-frequency equation for a multistep non-uniform rod with any number of steps.

Zeng and Bert [8] used the differential transformation method for vibration analyses of tapered bars with fixed-end conditions, assuming a linear variation for both cross-sectional are and mass.

Raj and Sujith [9] developed a family of closed-form solutions in terms of confluent hypergeometric functions for the longitudinal vibration of rods with variable cross-sectional areas, reducing the governing equation to confluent hypergeometric differential equations with a generic transformation. This study presents the eigenfrequencies of rods with certain area variations subjected to classical boundary conditions.

Elishakoff [10] provided exact solutions for the vibration and stability problems of non-uniform and inhomogeneous rods, beams, and plates.

Al-Kaisy et al. [11] studied the free vibration of a general non-uniform rod using the differential quadrature method to determine the non-dimensional natural frequency and normalized mode shapes of a non-uniform rod for free and clamped boundary conditions while accounting for the influence of varying cross-sectional area on vibration.

Provatidis [12] proposed a novel global collocation method for eigenvalue analyses of freely vibrated elastic structures where the proposed methodology was designed to handle various types of boundary conditions, including instances of two Dirichlet and one Dirichlet and one Neumann condition.

Calio and Elishakoff [13] developed a class of closed-form solutions for longitudinally vibrating inhomogeneous rods for a given distribution of material density that were clamped at one end and free at the other, yielding distributions of axial rigidity, which, together with a specific law of material density, satisfied the governing eigenvalue problem.

Arndt et al. [14] introduced an adaptive generalized finite element method for analyzing the free longitudinal vibrations of straight bars and trusses. The method involves enriching the standard finite element method space with functions that depend on the geometric and mechanical properties of the element.

Inaudi and Matusevich [15] devised a technique based on power series with domain partitioning, and it was presented in a matrix formulation, which effectively solved linear differential equations up to a desired degree of accuracy. This approach has been proposed as an alternative to other power series techniques employed in vibration analysis. This technique was presented in a study of the longitudinal vibration of a rod with a linearly varying cross-sectional area.

Guo and Yang [16] proposed a series solution for the vibration of arbitrary non-uniform rods with four types of profiles and variations in geometry or material properties and compared their results with the solutions obtained using the WKB method.

Yardimoglu and Aydin [17] used appropriate transformations to obtain exact solutions for the longitudinal natural vibration frequencies of rods with cross-sectional variations as the power of sinusoidal functions. The transformation reduces the governing equation to the associated Legendre equation, which is the frequency equation of a rod with a certain cross-sectional area variation and boundary conditions. The effects of variations in the cross-sectional area of the rods on the natural characteristics were also considered.

Shahba et al. [18] analyzed the longitudinal and transverse vibrations and stability of axially functionally graded beams using the finite element method.

Shahba and Rajasekaran [19] studied the free vibration and stability of axially functionally graded tapered Euler–Bernoulli beams, which were solved using the governing differential equations of motion and differential transform element method (DTEM) based on the differential transform method (DTM).

Shahba et al. [20] studied the structural analysis of axially functionally graded tapered beams from a mechanical point of view using the finite element method by introducing the concept of basic displacement functions.

Gan et al. [21] investigated longitudinal wave propagation in a rod with a variable cross-section using the transfer matrix method by establishing the equation of motion for the rod based on the elementary wave theory, the Love theory, and the Mindlin–Herrmann theory. Two types of rods with cross-sections varying in the exponential and polynomial forms were considered to illustrate the analytical predictions of the propagation characteristics of the longitudinal wave, and the results were compared with the results from the finite element analysis (FEA) method.

Hong et al. [22] presented a spectral element model for FGM axial bars from the governing equations of motion in the frequency domain using the variational method to analyze functionally graded material bars with respect to axial or longitudinal motions; the bars’ material properties vary in the radial direction according to the power law. The radial contraction was employed by adopting the Mindlin–Herrmann rod theory, and the model was verified by comparing it with finite element solutions.

Shokrollahi and Nejad [23] investigated the longitudinal free vibrations of non-uniform rods with nonlinear governing equations using discrete singular convolution, employing the regularized Shannon delta kernel.

Guo and Yang [24] proposed an iterative method that resulted in a series solution for the free and steady-state forced longitudinal vibrations of non-uniform rods. The convergence and linear independence of the proposed method were verified using convergence tests and the nonzero value of the corresponding Wronskian determinant.

Shali et al. [25] analyzed the axial vibration of non-uniform rods with different end conditions using the differential transform method.

Šalinić et al. [26] proposed a non-iterative computational technique to study the free vibrations of axially functionally graded tapered, stepped, and continuously segmented rods and beams with elastically restrained ends with attached masses. The proposed method was referred to as the symbolic–numeric method of initial parameters (SNMIPs) that stemmed from the modification of the iterative numerical method of initial parameters in differential forms known in the literature.

Celebi et al. [27] used the complementary function method in the spatial domain with a Laplace transform in the time domain for the forced vibration of cantilever rods having material properties and cross-sectional areas arbitrarily varying in the axial direction.

Pillutla et al. [28] studied the longitudinal vibrations of functionally graded rods with variable cross-sectional areas and material properties using the pseudospectral method.

Jedrysiak [29] considered the vibrations of microstructured periodic slender beams, and axial forces were considered to analyze the effect of the microstructure size of the beams on their vibrations, applying the general tolerance modelling method and standard modelling methods based on two various concepts—weakly slowly varying functions and slowly varying functions.

Jedrysiak [30] considered slender elastic nonperiodic beams with the axially functionally graded structure on the macro-level along the x-axis and a nonperiodic structure on the micro-level, applying the tolerance modelling method to derive the model equations of the general tolerance model and standard tolerance model and describing dynamics and stability for axially functionally graded beams with the microstructure.

The governing equation of axial rod vibrations is a differential equation with variable coefficients. Analytical solutions are available only for limited cases of cross-sectional area variations. Hence, alternative analysis techniques are employed for these situations, such as the weighted residual method, Ritz method, finite difference method, and finite element method. In the last two decades, some analytical approximation methods have gained much popularity in the solution of linear/nonlinear ordinary/partial differential equations and have been applied to various applied mechanics problems. Some of these analytical approximation methods include the differential transform method (DTM) [31], variational iteration method (VIM) [32], Adomian decomposition method (ADM) [33], homotopy perturbation method (HPM) [34], homotopy analysis method (HAM) [35], and optimal auxiliary functions method (OAFM) [36]. The application of these methods involves conducting analytical integration; however, they may be difficult to integrate, especially if expressions include singular and transcendental terms.

In this study, the difficulty in the analytical integration process is eliminated by dividing the problem domain into a number of subdomains such that in each domain, variable properties are assumed as constants that are computed at the center of the subdomain. This assumption accelerates the analytical approximation process. Specifically, VIM and ADM are used in view of the approach that is the subject of the present study.

2. Axial Vibration of Rods

Rod vibrations are a fundamental topic in textbooks on mechanical vibrations [35,36,37]. The rod shown in the following figure was considered.

In Figure 1, L is the length of the rod, x is the axial coordinate, w is the axial displacement, and P is the force acting on an infinitesimal element. The governing equation for the free axial vibration of the rod, as shown in Figure 1, is as follows:

$$\frac{\partial }{\partial \mathrm{x}}\left[\mathrm{E}\mathrm{A}(\mathrm{x})\frac{\partial \mathrm{w}(\mathrm{x},\mathrm{t})}{\partial \mathrm{x}}\right]=\rho \mathrm{A}(\mathrm{x})\frac{{\partial }^2\mathrm{w}(\mathrm{x},\mathrm{t})}{\partial {\mathrm{t}}^2}$$

(1)

Applying the separation of variables techniques with $\mathrm{w}(\mathrm{x},\mathrm{t})=\mathrm{u}(\mathrm{x})\mathrm{T}(\mathrm{t})$, the governing equation becomes the following:

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{x}}\left[\mathrm{E}\mathrm{A}(\mathrm{x})\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}\right]+{\omega }^2\rho \mathrm{A}(\mathrm{x})\mathrm{u}(\mathrm{x})=0$$

(2)

If the cross-sectional area is constant, that is, A(x) = A, Equation (2) takes the following form:

$$\frac{{\partial }^2\mathrm{u}(\mathrm{x})}{\partial {\mathrm{x}}^2}+{\left(\frac{\omega }{\mathrm{c}}\right)}^2\mathrm{u}(\mathrm{x})=0$$

(3)

where $\mathrm{c}=\sqrt{\mathrm{E}/\rho }$. Since Equation (3) is a differential equation with constant coefficients, it is much easier to solve than Equation (2).

If the bar is axially functionally graded, then the governing equation reads as follows:

$$\frac{\partial }{\partial \mathrm{x}}\left[\mathrm{E}(\mathrm{x})\mathrm{A}(\mathrm{x})\frac{\partial \mathrm{w}(\mathrm{x},\mathrm{t})}{\partial \mathrm{x}}\right]=\rho (\mathrm{x})\mathrm{A}(\mathrm{x})\frac{{\partial }^2\mathrm{w}(\mathrm{x},\mathrm{t})}{\partial {\mathrm{t}}^2}$$

(4)

Applying the separation of variables technique, Equation (4) becomes the following [19]:

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{x}}\left[\mathrm{E}(\mathrm{x})\mathrm{A}(\mathrm{x})\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}\right]+{\omega }^2\rho (\mathrm{x})\mathrm{A}(\mathrm{x})\mathrm{u}(\mathrm{x})=0$$

(5)

If the axial stiffness and density are constant, Equation (5) takes the form given in Equation (3). Equation (5) is now expanded to produce the following equation.

$$\frac{{\mathrm{d}}^2\mathrm{u}(\mathrm{x})}{\mathrm{d}{\mathrm{x}}^2}+\frac{1}{\mathrm{E}(\mathrm{x})\mathrm{A}(\mathrm{x})}\frac{\mathrm{d}}{\mathrm{d}\mathrm{x}}\left[\mathrm{E}(\mathrm{x})\mathrm{A}(\mathrm{x})\right]\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}+{\omega }^2\frac{\rho (\mathrm{x})}{\mathrm{E}(\mathrm{x})}\mathrm{u}(\mathrm{x})=0$$

(6)

For homogeneous rods, Equation (6) reduces into a simpler form.

$$\frac{{\mathrm{d}}^2\mathrm{u}(\mathrm{x})}{\mathrm{d}{\mathrm{x}}^2}+\frac{1}{\mathrm{A}(\mathrm{x})}\frac{\mathrm{d}}{\mathrm{d}\mathrm{x}}\left[\mathrm{A}(\mathrm{x})\right]\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}+{\left(\frac{\omega }{\mathrm{c}}\right)}^2\mathrm{u}(\mathrm{x})=0$$

(7)

The exact solutions of Equations (6) and (7) are available only for a limited number of cases. This study aimed to provide an analytical-based numerical approach for solving both equations for the computation of free axial vibration frequencies.

3. The Method

In this section, a brief explanation of the proposed technique is given.

3.1. ADM

ADM [33] is a powerful technique for solving linear/nonlinear ordinary/partial differential equations and has attracted significant attention in the applied sciences.

Consider an equation of the following form.

$$\mathrm{L}\mathrm{y}+\mathrm{R}\mathrm{y}+\mathrm{N}\mathrm{y}=\mathrm{f}(\mathrm{x})$$

(8)

where L is the linear operator of maximum order, N is the nonlinear operator, and R is the operator for the remaining terms. We assume that L is a second-order derivative: i.e., $\mathrm{L}={\mathrm{d}}^2/\mathrm{d}{\mathrm{x}}^2$. Then, the inverse operator of L is

$${\mathrm{L}}^{-1}\left(\cdot \right)={\displaystyle \underset{0}{\overset{\mathrm{x}}{\int }}{\displaystyle \underset{0}{\overset{\mathrm{t}}{\int }}\left(.\right)}}\mathrm{d}\tau \mathrm{d}\mathrm{t}$$

(9)

If all terms on the left-hand side of the equation in Equation (4) are taken to the right-hand side except for term Ly, then, by applying the inverse operator to both sides of the equation, the following relation is obtained:

$$\mathrm{y}(\mathrm{x})=\mathrm{y}(0)+\mathrm{x}{\mathrm{y}}^{\prime }(0)+\mathrm{g}(\mathrm{x})-{\mathrm{L}}^{-1}\mathrm{R}\mathrm{y}-{\mathrm{L}}^{-1}\mathrm{N}\mathrm{y}$$

(10)

where g(x) is obtained by integrating function f(x). The solution is defined as follows:

$$\mathrm{y}(\mathrm{x})={\displaystyle \sum _{\mathrm{n}=0}^{\infty }{\mathrm{y}}_{\mathrm{n}}(\mathrm{x})}$$

(11)

$$\mathrm{N}\mathrm{y}(\mathrm{x})={\displaystyle \sum _{\mathrm{n}=0}^{\infty }{\mathrm{A}}_{\mathrm{n}}(\mathrm{x})}$$

(12)

In Equation (8), A_n(x) is the nth Adomian polynomial and is defined as follows:

$${\mathrm{A}}_{\mathrm{k}}=\frac{1}{\mathrm{k}!}\frac{\partial }{\partial {\lambda }^{\mathrm{k}}}{\left[\mathrm{N}\left({\displaystyle \sum _{\mathrm{n}=0}^{\infty }{\mathrm{y}}_{\mathrm{n}}{\lambda }^{\mathrm{n}}}\right)\right]}_{\lambda =0}\mathrm{k}=0,1,2,\dots$$

(13)

Inserting Equations (11) and (12) into Equation (10), the following successive relations are obtained.

$${\mathrm{y}}_0(\mathrm{x})=\mathrm{y}(0)+\mathrm{x}{\mathrm{y}}^{\prime }(0)+\mathrm{g}(\mathrm{x})$$

(14)

$${\mathrm{y}}_{\mathrm{k}}(\mathrm{x})=-{\mathrm{L}}^{-1}\mathrm{R}{\mathrm{y}}_{\mathrm{k}-1}(\mathrm{x})-{\mathrm{L}}^{-1}{\mathrm{A}}_{\mathrm{k}-1}(\mathrm{x}) \mathrm{k}\ge 1$$

(15)

Finally, the solution can be defined in terms of an infinite series as follows:

$$\mathrm{y}(\mathrm{x})=\underset{\mathrm{n}\to \infty }{\lim }{\displaystyle \sum _{\mathrm{k}=0}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{k}}(\mathrm{x})}$$

(16)

After calculating the N terms, the solution to Equation (16) with these terms is called an Nth-order solution.

3.2. VIM

VIM is an analytical technique used to solve linear/nonlinear ordinary/partial differential equations [32]. This method rapidly converges to accurate results.

We consider the following differential equation:

$$\mathrm{L}\mathrm{y}+\mathrm{N}\mathrm{y}=\mathrm{f}(\mathrm{x})$$

(17)

where L is a linear operator, N is a nonlinear operator, and f(x) is a known non-homogeneous term. VIM proposes a correction function for the solution to Equation (17):

$${\mathrm{y}}_{\mathrm{n}}(\mathrm{x})={\mathrm{y}}_{\mathrm{n}-1}(\mathrm{x})+{\displaystyle \underset{0}{\overset{\mathrm{x}}{\int }}\lambda (\tau )}\left[\mathrm{L}{\mathrm{y}}_{\mathrm{n}-1}(\tau )+\mathrm{N}{\tilde{\mathrm{y}}}_{\mathrm{n}-1}(\tau )-\mathrm{f}(\tau )\right] \mathrm{n}\ge 1$$

(18)

where l is a generalized Lagrange multiplier based on variational theory. In addition, ${\mathrm{y}}_{\mathrm{n}}\left(\mathrm{x}\right)$ is the nth-order approximation for y(x), and ${\tilde{\mathrm{y}}}_{\mathrm{n}}\left(\mathrm{x}\right)$ is a limited variation of the nth-order approximation: that is, ${\delta \tilde{\mathrm{y}}}_{\mathrm{n}}\left(\tau \right)=0$. Considering the limited variation, if a variation is applied to Equation (18), the following relation is obtained.

$$\delta {\mathrm{y}}_{\mathrm{n}}(\mathrm{x})=\delta {\mathrm{y}}_{\mathrm{n}-1}(\mathrm{x})+\delta \left({\displaystyle \underset{0}{\overset{\mathrm{x}}{\int }}\lambda (\tau )}\mathrm{L}{\mathrm{y}}_{\mathrm{n}-1}(\tau )\right)$$

(19)

Expanding the parenthesis, the following equation can be written.

$$\delta {\mathrm{y}}_{\mathrm{n}}(\mathrm{x})=\delta {\mathrm{y}}_{\mathrm{n}-1}(\mathrm{x})+{\left[\lambda (\tau ){\displaystyle {\int }_0^{\tau }\mathrm{L}\delta {\mathrm{y}}_{\mathrm{n}-1}(\xi )\mathrm{d}\xi }\right]}_{\tau =0}^{\tau =\mathrm{x}}-{\displaystyle {\int }_0^{\mathrm{x}}{\lambda }^{\prime }(\tau )\left({\displaystyle {\int }_0^{\tau }\mathrm{L}\delta {\mathrm{y}}_{\mathrm{n}-1}(\xi )\mathrm{d}\xi }\right)}\mathrm{d}\tau$$

(20)

By applying stationary boundary conditions to Equation (20), the optimal value of the Lagrange multiplier is obtained. Once the Lagrange multiplier has been determined, the solution to Equation (17) can be calculated using the following successive approximations:

$${\mathrm{y}}_{\mathrm{n}}(\mathrm{x})={\mathrm{y}}_{\mathrm{n}-1}(\mathrm{x})+{\displaystyle \underset{0}{\overset{\mathrm{x}}{\int }}\lambda (\tau )}\left[\mathrm{L}{\mathrm{y}}_{\mathrm{n}-1}(\tau )+\mathrm{N}{\mathrm{y}}_{\mathrm{n}-1}(\tau )-\mathrm{f}(\tau )\right] \mathrm{n}\ge 1$$

(21)

The initial approximation y₀(x) is predefined and is chosen generally as the solution for the linear operator.

3.3. Subdomain-Based Numerical Solution Approach

The rod considered in this study was divided into subdomains, as shown in the figure below.

The governing equation in Equation (6) can be rearranged in the following form:

$$\frac{{\mathrm{d}}^2\mathrm{u}(\mathrm{x})}{\mathrm{d}{\mathrm{x}}^2}+\frac{{\mathrm{S}}^{\prime }(\mathrm{x})}{\mathrm{S}(\mathrm{x})}\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}+{\omega }^2\frac{\mathrm{m}(\mathrm{x})}{\mathrm{S}(\mathrm{x})}\mathrm{u}(\mathrm{x})=0$$

(22)

where $\mathrm{S}(\mathrm{x})=\mathrm{E}(\mathrm{x})\mathrm{A}(\mathrm{x})$ and $\mathrm{m}(\mathrm{x})=\rho (\mathrm{x})\mathrm{A}(\mathrm{x})$. The governing equation given in the ith domain is shown in Figure 2, which can be written as follows:

$$\frac{{\mathrm{d}}^2\mathrm{u}(\mathrm{x})}{\mathrm{d}{\mathrm{x}}^2}+\frac{{\mathrm{S}}^{\prime }({\tilde{\mathrm{x}}}_{\mathrm{i}})}{\mathrm{S}({\tilde{\mathrm{x}}}_{\mathrm{i}})}\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}+{\omega }^2\frac{\mathrm{m}({\tilde{\mathrm{x}}}_{\mathrm{i}})}{\mathrm{S}({\tilde{\mathrm{x}}}_{\mathrm{i}})}\mathrm{u}(\mathrm{x})=0$$

(23)

where ${\tilde{\mathrm{x}}}_{\mathrm{i}}$ is in the interval of the ith subdomain and may be performed in different ways. In this study, ${\tilde{\mathrm{x}}}_{\mathrm{i}}$ was assumed to be the axial coordinate of the center of the ith subdomain. Once ${\mathrm{S}}^{\prime }(\mathrm{x})/\mathrm{S}(\mathrm{x})$ and $\mathrm{m}(\mathrm{x})/\mathrm{S}(\mathrm{x})$ are calculated at ${\tilde{\mathrm{x}}}_{\mathrm{i}}$, they are assumed to be constant throughout the subdomain, which leads to a differential equation with constant coefficients. This type of arrangement simplifies the process described in Equations (15) and (21): for a sufficiently fine division of the problem domain, the central properties represent the sufficient properties of the related subdomain.

In addition to Equation (23), two continuity conditions should be satisfied between the two subdomains. We consider the subdomains shown in Figure 2. Between subdomain i and subdomain i − 1 at node i, the continuity conditions are as follows:

$${\mathrm{u}}_2^{\mathrm{i}-1}={\mathrm{u}}_1^{\mathrm{i}}$$

(24)

$${\left({\mathrm{u}}^{\prime }\right)}_2^{\mathrm{i}-1}={\left({\mathrm{u}}^{\prime }\right)}_1^{\mathrm{i}}$$

(25)

Finally, boundary conditions were applied for different end conditions. The boundary conditions for a fixed and free end are given as follows [35,36,37]:

$${\left.\mathrm{u}(\mathrm{x})\right|}_{\mathrm{at} \mathrm{fixed} \mathrm{end}}=0$$

(26)

$${\left.\frac{\mathrm{d}\mathrm{u}(\mathrm{x})}{\mathrm{d}\mathrm{x}}\right|}_{\mathrm{at} \mathrm{free} \mathrm{end}}=0$$

(27)

With respect to the initial approximation for the solution of Equation (23), each subdomain is described by the following equation, which satisfies the requirements for the initial approximation of VIM and ADM:

$${\mathrm{u}}_0^{\mathrm{i}}(\xi )={\mathrm{A}}_{\mathrm{i}}+{\mathrm{B}}_{\mathrm{i}}\xi$$

(28)

where ξ is the local coordinate. The Nth-order solution is sought within each subdomain according to Equations (15) and (21) based on ADM and VIM, respectively. A second-order solution for subdomain i via ADM with the initial approximation given in Equation (28) reads as follows:

$${\mathrm{u}}_{\mathrm{i}}^2(\xi )={\mathrm{A}}_{\mathrm{i}}\left(1-{\beta }_{\mathrm{i}}\frac{{\xi }^2}{2}+{\alpha }_{\mathrm{i}}{\beta }_{\mathrm{i}}\frac{{\xi }^3}{6}+{\beta }_{\mathrm{i}}^2\frac{{\xi }^4}{24}\right)+{\mathrm{B}}_{\mathrm{i}}\left(\xi -{\alpha }_{\mathrm{i}}\frac{{\xi }^2}{2}+{\alpha }_{\mathrm{i}}^2\frac{{\xi }^3}{6}-{\beta }_{\mathrm{i}}\frac{{\xi }^3}{6}+{\alpha }_{\mathrm{i}}{\beta }_{\mathrm{i}}\frac{{\xi }^4}{12}+{\beta }_{\mathrm{i}}^2\frac{{\xi }^5}{120}\right)$$

(29)

where ${\alpha }_{\mathrm{i}}={\mathrm{S}}^{\prime }({\tilde{\mathrm{x}}}_{\mathrm{i}})/\mathrm{S}({\tilde{\mathrm{x}}}_{\mathrm{i}})$ and ${\beta }_{\mathrm{i}}=\mathrm{m}({\tilde{\mathrm{x}}}_{\mathrm{i}})/\mathrm{S}({\tilde{\mathrm{x}}}_{\mathrm{i}})$. Then, the derivative of the second-order axial displacement with respect to the local coordinate would be as follows:

$$\frac{\mathrm{d}{\mathrm{u}}_{\mathrm{i}}^2(\xi )}{\mathrm{d}\xi }={\mathrm{A}}_{\mathrm{i}}\left(-{\beta }_{\mathrm{i}}\xi +{\alpha }_{\mathrm{i}}{\beta }_{\mathrm{i}}\frac{{\xi }^2}{2}+{\beta }_{\mathrm{i}}^2\frac{{\xi }^3}{6}\right)+{\mathrm{B}}_{\mathrm{i}}\left(1-{\alpha }_{\mathrm{i}}\xi +{\alpha }_{\mathrm{i}}^2\frac{{\xi }^2}{2}-{\beta }_{\mathrm{i}}\frac{{\xi }^2}{2}+{\alpha }_{\mathrm{i}}{\beta }_{\mathrm{i}}\frac{{\xi }^3}{3}+{\beta }_{\mathrm{i}}^2\frac{{\xi }^4}{24}\right)$$

(30)

Continuity conditions given in Equations (24) and (25) can be written with the second-order displacement in Equation (29) as follows:

$${\mathrm{A}}_{\mathrm{i}-1}\left(1-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{24}\right)+{\mathrm{B}}_{\mathrm{i}-1}\left(\xi -{\alpha }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{12}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^5}{120}\right)={\mathrm{A}}_{\mathrm{i}}$$

(31)

$${\mathrm{A}}_{\mathrm{i}-1}\left(-{\beta }_{\mathrm{i}-1}{\mathrm{L}}_{\mathrm{i}-1}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}\right)+{\mathrm{B}}_{\mathrm{i}-1}\left(1-{\alpha }_{\mathrm{i}-1}{\mathrm{L}}_{\mathrm{i}-1}+{\alpha }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{3}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{24}\right)={\mathrm{B}}_{\mathrm{i}}$$

(32)

where ${\mathrm{L}}_{\mathrm{i}-1}={\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}-1}$ is the length of the domain i − 1 (see Figure 2). Hence, second-order solutions in subdomains lead to the following relation in matrix form, which is obtained for node i.

$$\left[\begin{array}{rrrr}\hfill -\left(1-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{24}\right)& \hfill -\left(\xi -{\alpha }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{12}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^5}{120}\right)& \hfill 1& \hfill 0\\ \hfill -\left(-{\beta }_{\mathrm{i}-1}{\mathrm{L}}_{\mathrm{i}-1}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{6}\right)& \hfill -\left(1-{\alpha }_{\mathrm{i}-1}{\mathrm{L}}_{\mathrm{i}-1}+{\alpha }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}-{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^2}{2}+{\alpha }_{\mathrm{i}-1}{\beta }_{\mathrm{i}-1}\frac{{\mathrm{L}}_{\mathrm{i}-1}^3}{3}+{\beta }_{\mathrm{i}-1}^2\frac{{\mathrm{L}}_{\mathrm{i}-1}^4}{24}\right)& \hfill 0& \hfill 1\end{array}\right]\left\{\begin{array}{c}{\mathrm{A}}_{\mathrm{i}-1}\\ {\mathrm{B}}_{\mathrm{i}-1}\\ {\mathrm{A}}_{\mathrm{i}}\\ {\mathrm{B}}_{\mathrm{i}}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right\}$$

(33)

The fourth-order solution provides additional terms to the multiples of A_i and B_i given in the parentheses in Equation (29). In this study, second- and fourth-order solutions are used in each subdomain.

After applying all the required continuity and boundary conditions at each node, the following equation in the matrix form is obtained:

$$\left[\mathrm{K}(\omega )\right]\left\{\Lambda \right\}=\left\{0\right\}$$

(34)

In Equation (34), [K] is the coefficient matrix, which is a function of free vibration frequency w. {Λ} is a vector composed of coefficients A_i’s and B_i’s. A non-trivial solution to Equation (34) leads to the desired free vibration frequencies.

Below, the proposed approach is applied to a number of case studies available in the literature.

4. Numerical Results

In this section, the proposed technique is applied to available problems in the literature.

4.1. Case 1

Eisenberger [2] calculated the free axial vibration frequencies of a variable cross-section fixed–free rod with two different area functions. Both ADM and VIM are employed in the solution, and the results are compared below with different numbers of subdomains and solution orders within each subdomain. In the table, n represents the number of subdomains, and m represents the order of the solution within each subdomain.

From

Table 1. Comparison of results via ADM and VIM for $\mathrm{A}\left(\mathrm{x}\right)=2-\mathrm{x}$.

$\mathit{n}=10$				[2] First Five Frequencies
$\mathit{m}=2$		$\mathit{m}=4$
ADM	VIM	ADM	VIM
1.79466	1.79466	1.79410	1.79410	1.79401
4.80205	4.80205	4.80211	4.80211	4.80206
7.91027	7.91027	7.90899	7.90899	7.90896
11.05688	11.05688	11.03511	11.03511	11.03509
14.25171	14.25171	14.16801	14.16801	14.16799

and

Table 2. Comparison of results via ADM and VIM for $\mathrm{A}\left(\mathrm{x}\right)=3-4\mathrm{x}+2{\mathrm{x}}^2$.

$\mathit{n}=10$				[2] First Five Frequencies
$\mathit{m}=2$		$\mathit{m}=4$
ADM	VIM	ADM	VIM
1.97035	1.97035	1.97011	1.97011	1.97090
4.82192	4.82192	4.82058	4.82058	4.82076
7.91805	7.91805	7.91813	7.91813	7.91820
11.05315	11.05315	11.04139	11.04139	11.04144
14.23085	14.23085	14.17282	14.17282	14.17284

, it can be observed that ADM and VIM provide the same answers for each frequency with assumed values of n and m. Hence, from now on, the results of this study will be given as a single value for each of the cases considered. From Table 2, the second-order solution for the first frequency appears to be adequate, whereas for other frequencies, the fourth-order solution is required. The improvement of the solution with the number of subdomains is investigated below for both cases. Hence, the fourth-order solution was used in the solutions.

Table 3. Comparison of results with $\mathrm{A}\left(\mathrm{x}\right)=2-\mathrm{x}$ for a different number of subdomains.

This Study			[2] First Five Frequencies
$\mathit{n}=10$	$\mathit{n}=20$	$\mathit{n}=50$
1.79410	1.79403	1.79401	1.79401
4.80211	4.80207	4.80206	4.80206
7.90899	7.90897	7.90896	7.90896
11.03511	11.03510	11.03510	11.03509
14.16801	14.16799	14.16799	14.16799

and

Table 4. Comparison of results with $\mathrm{A}\left(\mathrm{x}\right)=3-4\mathrm{x}+2{\mathrm{x}}^2$ for a different number of subdomains.

$\mathit{n}=10$				[2] First Five Frequencies
$\mathit{m}=2$		$\mathit{m}=4$
ADM	VIM	ADM	VIM
1.97078	1.97088	1.97090	1.97090	1.97078
4.82071	4.82075	4.82076	4.82076	4.82071
7.91818	7.91820	7.91820	7.91820	7.91818
11.04143	11.04144	11.04144	11.04144	11.04143
14.17283	14.17284	14.17284	14.17284	14.17283

show that with a fine mesh, the results obtained using centroidal properties within each subdomain are in perfect agreement with the results in [2]. In conclusion, the proposed approach provides perfect results with the advantage of computation ease.

4.2. Case 2

Kumar and Sujith [5] considered specific cross-sectional area functions to determine the exact longitudinal vibration frequencies of non-uniform rods. However, numerical extraction from analytical expressions was computed incorrectly, and the results for various cases were corrected by other researchers [38,39]. The proposed approach was implemented with fourth-order solutions within each subdomain. The results are given in the tables below.

Table 5. Comparison of results for the fixed–fixed rod with $\mathrm{A}\left(\mathrm{x}\right)={\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)}^4$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	3.286029	3.286007	3.285998	3.474486	3.474339	3.474339
2	6.360710	6.360678	6.360671	6.480282	6.480034	6.480028
3	9.477253	9.477196	9.477180	9.561781	9.561370	9.561367
4	12.606002	12.605890	12.605802	12.671048	12.670362	12.670323
5	15.739881	15.739656	15.739648	15.792900	15.791752	15.791747
6	18.876445	18.876002	18.874533	18.921558	18.919655	18.919130

,

Table 6. Comparison of results for the fixed–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)}^4$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	0.824911	0.824969	0.824971	0.526525	0.526686	0.526694
2	4.600469	4.600454	4.600454	4.689503	4.689332	4.689329
3	7.789132	7.789097	7.789096	7.849000	7.848697	7.848695
4	10.949734	10.949659	10.949659	10.994129	10.993612	10.993610
5	14.101746	14.101592	14.101591	14.137111	14.136237	14.136235
6	17.250022	17.249709	17.249708	17.279715	17.278249	17.278247

and

Table 7. Comparison of results for the free–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)}^4$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	3.555887	3.555792	3.555788	4.041829	4.041340	4.041322
2	6.513149	6.513070	6.513068	6.857766	6.857191	6.857173
3	9.581360	9.581272	9.581270	9.830658	9.830000	9.829986
4	12.684744	12.684610	12.684608	12.877428	12.876563	12.876551
5	15.803120	15.802878	15.802877	15.959733	15.958463	15.958453
6	18.929253	18.928799	18.928798	19.061341	19.059368	19.059359

show very good agreement between the results of this study and those in [5,39]. Even for the largest mode with $n=20$, the maximum relative error is less than 0.00002. Hence, the proposed approach yields excellent results using the centroidal properties of each subdivision.

Different variations in the cross-sectional area were also considered in the same study [5].

In

Table 8. Comparison of results for fixed–fixed rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^2\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	2.978151	2.978187	2.978189	2.422827	2.422601	2.422727
2	6.203079	6.203097	6.203097	5.955421	5.956256	5.956376
3	9.371565	9.371576	9.371576	9.210152	9.210008	9.210127
4	12.526513	12.526518	12.526519	12.407305	12.406080	12.406195
5	15.676103	15.676100	15.676100	15.582357	15.580012	15.580119
6	18.823033	18.823011	18.823011	18.746557	18.743054	18.743152

,

Table 9. Comparison of results for fixed–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^2\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	1.517623	1.517637	1.517638	2.149426	2.148593	2.148560
2	4.702142	4.702145	4.702145	5.539852	5.535922	5.535762
3	7.848310	7.848311	7.848311	8.639936	8.633093	8.632812
4	10.991621	10.991621	10.991620	11.703833	11.695014	11.694640
5	14.134128	14.134123	14.134120	14.767991	14.758288	14.757860
6	17.276297	17.276824	17.276280	17.840719	17.831054	17.830600

and

Table 10. Comparison of results for free–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^2\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$, $\mathrm{b}=1$.

Mode	$\mathit{a}=1$			$\mathit{a}=2$
	This Study n = 20	This Study n = 100	[40]	This Study n = 20	This Study n = 100	[40]
1	3.309109	3.309071	3.309070	4.212406	4.209714	4.209604
2	6.375233	6.37509	6.375209	7.265756	7.260092	7.259860
3	9.487380	9.487363	9.487363	10.291902	10.283836	10.283498
4	12.613663	12.613649	12.613648	13.327822	13.318388	13.317980
5	15.745930	15.745914	15.745913	16.379162	16.369366	16.368917
6	18.881271	18.881240	18.881240	19.445177	19.435801	19.435335

, the results of this study agree very well with the analytical results [5,41]. The maximum relative error for the square sine function with $n=20$ is less than 0.00003. Consequently, the proposed approach worked very well in this case.

4.3. Case 3

Guo and Yang [16] obtained a series solution for a fixed–fixed rod using an exponential area function. The proposed method was applied to the case, and the first five frequencies for fixed–free and free–free rods were also computed for the researchers for further studies on the subject.

As shown in

Table 11. Comparison of results for the fixed–fixed rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\mathrm{e}}^{\mathrm{c}{\mathrm{x}}^2/{\mathrm{L}}^2}$, $\mathrm{L}=1$.

Mode	$\mathit{c}=0.5$		$\mathit{c}=1$		$\mathit{c}=2$
	This Study n = 100	Kummer Function [16]	This Study n = 100	Kummer Function [16]	This Study n = 100	Kummer Function [16]
1	3.231130281	3.231130281	3.339335867	3.339335867	3.603139793	3.603139793
2	6.329186675	6.329186675	6.387440255	6.387440254	6.539562685	6.539562675
3	9.455600357	9.455600357	9.494964059	9.494964058	9.598980448	9.598980427
4	12.58952819	12.58952819	12.61919022	12.61919022	12.69788674	12.69788671

, the results of this study are in perfect agreement with those in [9]. Guo et al. [9] only considered a fixed–fixed rod. As a contribution, we computed the first five frequencies of the fixed–free and free–free rods for the same area function using 100 subdomains in the solution region.

Table 12. First five frequencies of fixed–free and free–free rods with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\mathrm{e}}^{\mathrm{c}{\mathrm{x}}^2/{\mathrm{L}}^2}$, $\mathrm{L}=1$.

Mode	$\mathit{c}=0.5$		$\mathit{c}=1$		$\mathit{c}=2$
	Fixed–Free	Free–Free	Fixed–Free	Free–Free	Fixed–Free	Free–Free
1	1.414214	3.072491	1.263693	3.025089	0.985886	2.997101
2	4.667369	6.249688	4.638712	6.228916	4.631731	6.226225
3	7.827282	9.402573	7.810922	9.389054	7.809608	9.388313
4	10.976564	12.549750	10.965036	12.539695	10.964576	12.539391
5	14.122401	15.694679	14.113484	15.686666	14.113271	15.686512

presents the results below.

4.4. Case 4

Yardimoglu and Aydin [17] obtained exact solutions for the longitudinal natural vibration frequencies of rods with variable cross-sections. Yardimoglu considered the area function as a power of sine function, that is, ${\mathrm{s}\mathrm{i}\mathrm{n}}^n\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$, and calculated the frequencies of vibrations with different a, b, and n values for different end conditions. Nine different cases were considered [17]. Although the results of all cases in the paper agree very well with the proposed technique, only three of them are provided for comparison purposes. Below, the comparisons were chosen for the fixed values of $\mathrm{b}=1$ and $\mathrm{a}=\mathrm{\pi }-2\mathrm{b}$, and 100 subdomains were used in the analysis.

As shown in

Table 13. Comparison of results for the fixed–fixed rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^{\mathrm{n}}\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$.

Mode	$\mathit{n}=1$		$\mathit{n}=2$		$\mathit{n}=3$		$\mathit{n}=4$
	This Study	[17]	This Study	[17]	This Study	[17]	This Study	[17]
1	3.033656	3.033658	2.926834	2.926836	2.821242	2.821246	2.716998	2.717003
2	6.228474	6.228475	6.178606	6.178607	6.133694	6.133696	6.093837	6.093840
3	9.388171	9.388171	9.355383	9.355384	9.326455	9.326456	9.301423	9.301424
4	12.538877	12.538877	12.514409	12.514409	12.492984	12.492985	12.474620	12.474621
5	15.685953	15.685954	15.666425	15.666425	15.649387	15.649387	15.634848	15.634849
6	18.831208	18.831208	18.814954	18.814955	18.800802	18.800803	18.788757	18.788757

,

Table 14. Comparison of results for the fixed–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^{\mathrm{n}}\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$.

Mode	$\mathit{n}=1$		$\mathit{n}=2$		$\mathit{n}=3$		$\mathit{n}=4$
	This Study	[17]	This Study	[17]	This Study	[17]	This Study	[17]
1	1.568123	1.568123	1.560154	1.560155	1.547041	1.547041	1.529022	1.529023
2	4.715830	4.715830	4.726103	4.726102	4.743064	4.743064	4.766485	4.766484
3	7.856372	7.856372	7.863537	7.863537	7.875459	7.875459	7.892109	7.892108
4	10.997350	10.997350	11.002678	11.002678	11.011552	11.011552	11.023968	11.023967
5	14.138571	14.138571	14.142783	14.142783	14.149801	14.149801	14.159624	14.159624
6	17.279918	17.279918	17.283392	17.283392	17.289183	17.289183	17.297288	17.297288

and

Table 15. Comparison of results for the free–free rod with $\mathrm{A}\left(\mathrm{x}\right)={\mathrm{A}}_0{\sin }^{\mathrm{n}}\left(\mathrm{a}\mathrm{x}+\mathrm{b}\right)$.

Mode	$\mathit{n}=1$		$\mathit{n}=2$		$\mathit{n}=3$		$\mathit{n}=4$
	This Study	[17]	This Study	[17]	This Study	[17]	This Study	[17]
1	3.250524	3.250523	3.360331	3.360328	3.470895	3.470891	3.582098	3.582093
2	6.342611	6.342610	6.406607	6.406605	6.475018	6.475015	6.547679	6.547675
3	9.465160	9.465160	9.509270	9.509269	9.557058	9.557056	9.608469	9.608466
4	12.596871	12.596870	12.630356	12.630355	12.666807	12.666805	12.706199	12.706197
5	15.732444	15.732444	15.759386	15.759385	15.788776	15.788777	15.820607	15.820605
6	18.869994	18.869993	18.892515	18.892514	18.917112	18.917113	18.943782	18.943781

, the largest relative error was less than 2 × 10⁻⁶. The results obtained using the proposed technique agree very well with the analytical solutions [17]. Hence, the applied computational procedure was highly effective in this case.

4.5. Case 5

Shahba et al. [18], Shahba and Rajasekaran [19], and Pillutla et al. [28] studied the longitudinal free vibration of axially functionally graded rods using various computational approaches. They considered a tapered bar with a rectangular cross-section for which its breath taper ratio is c_b and height taper ratio is c_h, and its cross-section is given by $\mathrm{A}(\mathrm{x})=\left(1-{\mathrm{c}}_{\mathrm{b}}\overline{\mathrm{x}}\right)\left(1-{\mathrm{c}}_{\mathrm{h}}\overline{\mathrm{x}}\right)$, where $\overline{\mathrm{x}}=\mathrm{x}/\mathrm{L}$. The modulus of elasticity and mass density vary as $\mathrm{E}={\mathrm{E}}_0(1+\overline{\mathrm{x}})$ and $\rho ={\rho }_0\left(1+\overline{\mathrm{x}}+{\overline{\mathrm{x}}}^2\right)$. Numerical experiments showed that 20 subdomains were adequate for the first three frequencies. However, to determine the first six vibration frequencies, 100 subdomains were used in the calculations. In [18,19,28], the first three vibration frequencies of fixed–fixed and fixed–free rods were considered.

The efficiency of the proposed method can be observed in

Table 16. Comparison of results for fixed–fixed functionally graded rods.

ch	cb	${\mathit{\omega }}_1$			${\mathit{\omega }}_2$			${\mathit{\omega }}_3$
		[19]	[28]	This Study	[19]	[28]	This Study	[19]	[28]	This Study
0.0	0.0	2.8760	-	2.875963	5.7627	-	5.762720	8.6453	-	8.645279
	0.2	2.8631	-	2.863124	5.7562	-	5.756197	8.6409	-	8.640920
	0.4	2.8415	-	2.841461	5.7453	-	5.745320	8.6337	-	8.633663
	0.6	2.8023	-	2.802335	5.7251	-	5.725108	8.6201	-	8.620077
	0.8	2.7192	-	2.719200	5.6767	-	5.676634	8.5860	-	8.586022
0.2	0.2	2.8539	2.853926	2.853922	5.7515	5.751501	5.751493	8.6378	8.637841	8.637772
	0.4	2.8369	2.836936	2.836933	5.7430	5.742957	5.742966	8.6321	8.632147	8.632078
	0.6	2.8042	2.804223	2.804218	5.7260	5.726023	5.726014	8.6207	8.620747	8.620682
	0.8	2.7311	2.731119	2.731107	5.6829	5.682885	5.682900	8.5903	8.590337	8.590280
0.4	0.4	2.8260	2.825973	2.825970	5.7375	5.737498	5.737489	8.6284	8.628504	8.628435
	0.6	2.8016	2.801561	2.801557	5.7249	5.724929	5.724920	8.6200	8.620056	8.619992
	0.8	2.7415	2.741519	2.741508	5.6892	5.689170	5.689156	8.5947	8.594748	8.594692
0.6	0.6	2.7886	2.788644	2.788640	5.7188	5.718775	5.718765	8.6160	8.615999	8.615940
	0.8	2.7468	2.746813	2.746803	5.6942	5.694204	5.694192	8.5986	8.598631	8.598578
0.8	0.8	2.7340	2.734041	2.734030	5.6901	5.690056	5.690056	8.5965	8.596521	8.596473

and

Table 17. Comparison of results for fixed–free functionally graded rods.

ch	cb	${\mathit{\omega }}_1$			${\mathit{\omega }}_2$			${\mathit{\omega }}_3$
		[19]	[28]	This Study	[19]	[28]	This Study	[19]	[28]	This Study
0.0	0.0	1.1901	-	1.190121	4.2549	-	4.254908	7.1650	-	7.165005
	0.2	1.2503	-	1.250273	4.2683	-	4.268343	7.1728	-	7.172861
	0.4	1.3293	-	1.329260	4.2903	-	4.290291	7.1859	-	7.185878
	0.6	1.4400	-	1.440031	4.3333	-	4.333294	7.2123	-	7.212325
	0.8	1.6129	-	1.612889	4.4444	-	4.444389	7.2891	-	7.289126
0.2	0.2	1.3119	1.311936	1.311941	4.2841	4.284133	4.284152	7.1821	7.181803	7.182161
	0.4	1.3928	1.392789	1.392794	4.3090	4.309044	4.309063	7.1970	7.196715	7.197021
	0.6	1.5059	1.505925	1.505931	4.3559	4.355943	4.355962	7.2260	7.225773	7.225995
	0.8	1.6818	1.681840	1.681836	4.4726	4.472646	4.472622	7.3068	7.306703	7.306799
0.4	0.4	1.4759	1.475871	1.475877	4.3377	4.337722	4.337740	7.2143	7.214022	7.214275
	0.6	1.5917	1.591695	1.591701	4.3897	4.389690	4.389708	7.2466	7.246450	7.246626
	0.8	1.7706	1.770668	1.770663	4.5138	4.513827	4.513801	7.3329	7.332918	7.332981
0.6	0.6	1.7104	1.710443	1.710450	4.4487	4.448671	4.448690	7.2839	7.283791	7.283910
	0.8	1.8918	1.891820	1.891814	4.5832	4.583282	4.583249	7.3787	7.378728	7.378747
0.8	0.8	2.0723	2.072303	2.072286	4.7328	4.732920	4.732813	7.4888	7.488952	7.488819

. Due to the accurate results from the presented approach, additional information is presented in

Table 18. Additional non-dimensional frequencies for the functionally graded rod.

ch	cb	Fixed–Fixed Rod			Fixed–Free Rod
		${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$	${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$
0.0	0.0	11.527389	14.409395	17.291364	10.058501	12.946776	15.832722
	0.2	11.524117	14.406776	17.289181	10.064076	12.951100	15.836256
	0.4	11.518673	14.402421	17.285552	10.073349	12.958304	15.842147
	0.6	11.508451	14.394231	17.278720	10.092388	12.973163	15.854326
	0.8	11.482303	14.373052	17.260943	10.150175	13.019219	15.892511
0.2	0.2	11.521753	14.404884	17.287604	10.070687	12.956231	15.840450
	0.4	11.517481	14.401466	17.284756	10.081286	12.964471	15.847190
	0.6	11.508906	14.394595	17.279024	10.102171	12.980778	15.860559
	0.8	11.485529	14.375646	17.263112	10.162984	13.029246	15.900742
0.4	0.4	11.514750	14.399282	17.282936	10.093623	12.974071	15.855046
	0.6	11.508399	14.394194	17.278692	10.116990	12.992331	15.870023
	0.8	11.488913	14.378385	17.265409	10.182042	13.044192	15.913024
0.6	0.6	11.505375	14.391780	17.276683	10.144039	13.013508	15.887408
	0.8	11.492002	14.380928	17.267564	10.215716	13.070717	15.934868
0.8	0.8	11.490668	14.379970	17.266821	10.300149	13.138517	15.991278

and

Table 19. Non-dimensional frequencies for the free–free functionally graded rod.

ch	cb	${\mathit{\omega }}_1$	${\mathit{\omega }}_2$	${\mathit{\omega }}_3$	${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$
0.0	0.0	2.946681	5.794961	8.666278	11.543029	14.421883	17.301774
	0.2	2.934015	5.788683	8.662112	11.539910	14.419390	17.299698
	0.4	2.933673	5.788820	8.662252	11.540028	14.419489	17.299782
	0.6	2.962114	5.805445	8.673717	11.548734	14.426494	17.305639
	0.8	3.068501	5.880639	8.729109	11.592045	14.461882	17.335486
0.2	0.2	2.924899	5.784214	8.659156	11.537699	14.417624	17.298227
	0.4	2.929150	5.786684	8.660854	11.538987	14.418659	17.299092
	0.6	2.963729	5.806550	8.674503	11.549336	14.426982	17.306047
	0.8	3.078513	5.886890	8.733501	11.595402	14.464592	17.337756
0.4	0.4	2.939281	5.792205	8.664600	11.541814	14.420927	17.300985
	0.6	2.981602	5.816379	8.681183	11.554380	14.431030	17.309426
	0.8	3.106649	5.903695	8.745184	11.604301	14.471763	17.343756
0.6	0.6	3.033852	5.846773	8.702104	11.570253	14.443796	17.320096
	0.8	3.171214	5.944370	8.773819	11.626209	14.489453	17.358572
0.8	0.8	3.320973	6.058919	8.859947	11.693961	14.544927	17.405400

for the researchers for further studies on the subject.

In this study, additional computations are conducted for further research. These additional results are given in Table 18 and Table 19 below.

4.6. Case 6

Shahba et al. [20] and Pillutla et al. [28] studied the longitudinal free vibration of axially functionally graded rods with a cross-section given by $\mathrm{A}(\mathrm{x})=\left(1-\mathrm{c}\overline{\mathrm{x}}\right)$, where $\overline{\mathrm{x}}=\mathrm{x}/\mathrm{L}$. The modulus of elasticity and mass density varied according to $\mathrm{E}={\mathrm{E}}_0{\mathrm{e}}^{\overline{\mathrm{x}}}$ and $\rho ={\rho }_0{\mathrm{e}}^{\overline{\mathrm{x}}}$. As in the previous case, the numerical experiments showed that 20 subdomains were adequate for comparison. However, 100 subdomains were used to accurately determine the first six vibration frequencies. The results are summarized in the tables below.

As shown in

Table 20. Comparison of results for functionally graded rod.

c	Mode	Fixed–Fixed Rod			Fixed–Free Rod
		[20]	[28]	This Study n = 20	[20]	[28]	This Study n = 20
0.1	1	3.1757	3.172409	3.172409	1.2988	1.2985	1.298593
	2	6.3247	6.298648	6.298648	4.6478	4.637424	4.637424
	3	9.5228	9.435093	9.435095	7.8592	7.809505	7.809506
0.3	1	3.1514	3.148153	3.148153	1.3722	1.371958	1.371958
	2	6.3123	6.286365	6.286365	4.6656	4.655118	4.655118
	3	9.5144	9.426885	9.426885	7.8698	7.819942	7.819942
0.5	1	3.1120	3.108831	3.108818	1.4710	1.470676	1.470698
	2	6.2916	6.265918	6.265911	4.6983	4.687528	4.687538
	3	9.5005	9.413135	9.413131	7.8899	7.83957	7.839577
0.8	1	2.9780	2.975221	2.974925	1.7168	1.716251	1.716528
	2	6.2113	6.186568	6.186339	4.8486	4.836778	4.837135
	3	9.4427	9.357079	9.356909	7.9958	7.9433464	7.943773

, the results obtained from this study are in very good agreement with [28], even with 20 subdomains. However, for further research, additional computations were conducted, and the results shown in

Table 21. Non-dimensional frequencies for fixed–fixed rod.

c	${\mathit{\omega }}_1$	${\mathit{\omega }}_2$	${\mathit{\omega }}_3$	${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$
0.1	3.172409	6.298648	9.435093	12.574109	15.714154	18.854715
0.3	3.148154	6.286366	9.426885	12.567948	15.709223	18.850605
0.5	3.108831	6.265918	9.413135	12.557602	15.700935	18.843693
0.7	3.037652	6.225862	9.385513	12.536596	15.684013	18.829536
0.9	2.865805	6.107338	9.294870	12.463441	15.622847	18.777080

,

Table 22. Non-dimensional frequencies for fixed–free rod.

c	${\mathit{\omega }}_1$	${\mathit{\omega }}_2$	${\mathit{\omega }}_3$	${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$
0.1	1.298593	4.637424	7.809505	10.963902	14.112563	17.258642
0.3	1.371958	4.655118	7.819942	10.971323	14.118324	17.263351
0.5	1.470678	4.687529	7.839571	10.985386	14.129277	17.272319
0.7	1.614399	4.760055	7.886977	11.020261	14.156765	17.294968
0.9	1.853356	4.987153	8.076173	11.177288	14.289295	17.408893

and

Table 23. Non-dimensional frequencies for free–free rod.

C	${\mathit{\omega }}_1$	${\mathit{\omega }}_2$	${\mathit{\omega }}_3$	${\mathit{\omega }}_4$	${\mathit{\omega }}_5$	${\mathit{\omega }}_6$
0.1	3.174157	6.299530	9.435682	12.574551	15.714508	18.855010
0.3	3.168121	6.296532	9.433687	12.573055	15.713312	18.854014
0.5	3.182758	6.304872	9.439404	12.577387	15.716794	18.856923
0.7	3.246929	6.346320	9.468966	12.600161	15.735255	18.872421
0.9	3.458060	6.538647	9.630470	12.736129	15.851494	18.973406

were obtained according to different combinations of boundary conditions.

5. Conclusions

In this study, the axial free vibration of non-uniform, homogeneous, and functionally graded rods is considered. The governing equation of the problem is a variable coefficient differential equation, for which the analytical solutions are strictly limited. Analytical approximate solution methods are available for these types of problems. Although these methods are efficient, the analytical integration process becomes cumbersome or sometimes impossible for higher-order solutions to achieve better convergence. This study aims to propose a subdomain-based numerical solution approach to eliminate difficulties in the integration process while benefitting from the analytical formulation. Several case studies were considered, including polynomial, sinusoidal, and exponential variations of the cross-sectional area, for which the proposed technique led to excellent results. Functionally graded rods with variable cross-sectional areas and variable material properties were also considered. The proposed method was also very effective and yielded very good results. Numerical studies showed that 20 subdomains were sufficient for the first three vibration frequencies, whereas it was also sufficient for the first six vibration frequencies if the variation functions for the cross-sectional areas were not too complex. However, 100 subdomains were preferred for most comparative studies. Because the presented technique provided excellent results, additional computations were conducted for researchers who will aim to consider the problem in future research using different solution methods. In conclusion, the proposed approach is very effective in solving the vibration problems of homogeneous and functionally graded rods with variable cross-sectional areas.